]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finish
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 Oct 2006 16:09:31 +0000 (16:09 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 Oct 2006 16:09:31 +0000 (16:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@14126 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-21/doc/results.dox

index 91a95c86b25af9e1dbf66faf0ef2e81f5d76e22e..3e46d87170c369c05fb7588f3b97ef28962ad026 100644 (file)
@@ -1,8 +1,179 @@
 <a name="Results"></a> <h1>Results</h1>
 
+If we run the program, we get the following kind of output:
+<pre>
+<code>
+Number of active cells: 1024
+Number of degrees of freedom: 4160 (2112+1024+1024)
+
+Timestep 1
+   22 CG Schur complement iterations for pressure.
+   1 CG iterations for saturation.
+   Now at t=0.0326742, dt=0.0326742.
+
+Timestep 2
+   17 CG Schur complement iterations for pressure.
+   1 CG iterations for saturation.
+   Now at t=0.0653816, dt=0.0327074.
+
+Timestep 3
+   17 CG Schur complement iterations for pressure.
+   1 CG iterations for saturation.
+   Now at t=0.0980651, dt=0.0326836.
+
+...
+</code>
+</pre>
+As we can see, the time step is pretty much constant right from the start,
+which indicates that the velocities in the domain are not strongly dependent
+on changes in saturation, although they certainly are through the factor
+$\lambda(S)$ in the pressure equation.
+
+Our second observation is that the number of CG iterations needed to solve the
+pressure Schur complement equation drops from 22 to 17 between the first and
+the second time step (in fact, it remains around 17 for the rest of the
+computations). The reason is actually simple: Before we solve for the pressure
+during a time step, we don't reset the <code>solution</code> variable to
+zero. The pressure (and the other variables) therefore have the previous time
+step's values at the time we get into the CG solver. Since the velocities and
+pressures don't change very much as computations progress, the previous time
+step's pressure is actually a good initial guess for this time step's
+pressure. Consequently, the number of iterations we need once we have computed
+the pressure once is significantly reduced.
+
+The final observation concerns the number of iterations needed to solve for
+the saturation, i.e. one. This shouldn't surprise us too much: the matrix we
+have to solve with is the mass matrix. However, this is the mass matrix for
+the $DQ_0$ element of piecewise constants where no element couples with the
+degrees of freedom on neighboring cells. The matrix is therefore a diagonal
+one, and it is clear that we should be able to invert this matrix in a single
+CG iteration.
+
+
+With all this, here are a few movies that show how the saturation progresses
+over time. First, this is for the single crack model, as implemented in the
+<code>SingleCurvingCrack::KInverse</code> class:
 
 @image html step-21.centerline.gif
 
+As can be seen, the water rich fluid snakes its way mostly along the
+high-permeability zone in the middle of the domain, whereas the rest of the
+domain is mostly impermeable. This and the next movie are generated using
+<code>n_refinement_steps=7</code>, leading to a $128\times 128$ mesh with some
+16,000 cells and about 66,000 unknowns in total.
+
+
+The second movie shows the saturation for the random medium model of class
+<code>RandomMedium::KInverse</code>, where we have randomly distributed
+centers of high permeability and fluid hops from one of these zones to
+the next:
+
 @image html step-21.random2d.gif
 
+
+Finally, here is the same situation in three space dimensions, on a mesh with
+<code>n_refinement_steps=5</code>, which produces a mesh of some 32,000 cells
+and 167,000 degrees of freedom:
+
 @image html step-21.random3d.gif
+
+To repeat these computations, all you have to do is to change the line
+<pre>
+<code>
+      TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
+</code>
+</pre>
+in the main function to
+<pre>
+<code>
+      TwoPhaseFlowProblem<3> two_phase_flow_problem(0);
+</code>
+</pre>
+The visualization uses a cloud technique, where the saturation is indicated by
+colored but transparent clouds for each cell. This way, one can also see
+somewhat what happens deep inside the domain. A different way of visualizing
+would have been to show isosurfaces of the saturation evolving over
+time. There are techniques to plot isosurfaces transparently, so that one can
+see several of them at the same time like the layers of an onion.
+
+So why don't we show such isosurfaces? The problem lies in the way isosurfaces
+are computed: they require that the field to be visualized is continuous, so
+that the isosurfaces can be generated by following contours at least across a
+single cell. However, our saturation field is piecewise constant and
+discontinuous. If we wanted to plot an isosurface for a saturation $S=0.5$,
+chances would be that there is no single point in the domain where that
+saturation is actually attained. If we had to define isosurfaces in that
+context at all, we would have to take the interfaces between cells, where one
+of the two adjacent cells has a saturation greater than and the other cell a
+saturation less than 0.5. However, it appears that most visualization programs
+are not equipped to do this kind of transformation.
+
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+
+<h4>Solvers</h4>
+
+At present, the program is not particularly fast: the 2d random medium
+computation took about a day for the 1,000 or so time steps. The corresponding
+3d computation took almost two days for 800 time steps. The reason why it
+isn't faster than this is twofold. First, we rebuild the entire matrix in
+every time step, although some parts such as the $B$, $B^T$, and $M^S$ blocks
+never change. 
+
+Second, we could do a lot better with the solver and
+preconditioners. Presently, we solve the Schur complement $B^TM^u(S)^{-1}B$
+with a CG method, using $[B^T (\textrm{diag}(M^u(S)))^{-1} B]^{-1}$ as a
+preconditioner. Applying this preconditioner is expensive, since it involves
+solving a linear system each time. This may have been appropriate for @ref
+step_20 "step-20", where we have to solve the entire problem only
+once. However, here we have to solve it hundreds of times, and in such cases
+it is worth considering a preconditioner that is more expensive to set up the
+first time, but cheaper to apply later on.
+
+One possibility would be to realize that the matrix we use as preconditioner,
+$B^T (\textrm{diag}(M^u(S)))^{-1} B$ is still sparse, and symmetric on top of
+that. If one looks at the flow field evolve over time, we also see that while
+$S$ changes significantly over time, the pressure hardly does and consequently
+$B^T (\textrm{diag}(M^u(S)))^{-1} B \approx B^T (\textrm{diag}(M^u(S^0)))^{-1}
+B$. In other words, the matrix for the first time step should be a good
+preconditioner also for all later time steps.  With a bit of
+back-and-forthing, it isn't hard to actually get a representation of it as a
+SparseMatrix object. We could then hand it off to the SparseMIC class to form
+a sparse incomplete Cholesky decomposition. To form this decomposition is
+expensive, but we have to do it only once in the first time step, and can then
+use it as a cheap preconditioner in the future. We could do better even by
+using the SparseDirectUMFPACK class that produces not only an incomplete, but
+a complete decomposition of the matrix, which should yield an even better
+preconditioner. 
+
+Finally, why use the approximation $B^T (\textrm{diag}(M^u(S)))^{-1} B$ to
+precondition $B^T M^u(S)^{-1} B$? The latter matrix, after all, is the mixed
+form of the Laplace operator on the pressure space, for which we use linear
+elements. We could therefore build a separate matrix $A^p$ on the side that
+directly corresponds to the non-mixed formulation of the Laplacian, for
+example using the bilinear form $(\mathbf{K}\lambda(S^n) \nabla
+\varphi_i,\nabla\varphi_j)$. We could then form an incomplete or complete
+decomposition of this non-mixed matrix and use it as a preconditioner of the
+mixed form.
+
+Using such techniques, it can reasonably be expected that the solution process
+will be faster by at least an order of magnitude.
+
+
+<h4>Adaptivity</h4>
+
+Adaptivity would also clearly help. Looking at the movies, one clearly sees
+that most of the action is confined to a relatively small part of the domain
+(this particularly obvious for the saturation, but also holds for the
+velocities and pressures). Adaptivity can therefore be expected to keep the
+necessary number of degrees of freedom low, or alternatively increase the
+accuracy.
+
+On the other hand, adaptivity for time dependent problems is not a trivial
+thing: we would have to change the mesh every few time steps, and we would
+have to transport our present solution to the next mesh every time we change
+it (something that the SolutionTransfer class can help with). These are not
+insurmountable obstacles, but they do require some additional coding and more
+than we felt comfortable was worth packing into this tutorial program.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.