From a75af546bcf8182117783cd797a5a360fcaca275 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 30 Oct 2006 16:09:31 +0000 Subject: [PATCH] Finish git-svn-id: https://svn.dealii.org/trunk@14126 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-21/doc/results.dox | 171 +++++++++++++++++++++++ 1 file changed, 171 insertions(+) diff --git a/deal.II/examples/step-21/doc/results.dox b/deal.II/examples/step-21/doc/results.dox index 91a95c86b2..3e46d87170 100644 --- a/deal.II/examples/step-21/doc/results.dox +++ b/deal.II/examples/step-21/doc/results.dox @@ -1,8 +1,179 @@

Results

+If we run the program, we get the following kind of output: +
+
+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.
+
+...
+
+
+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 solution 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 +SingleCurvingCrack::KInverse 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 +n_refinement_steps=7, 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 +RandomMedium::KInverse, 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 +n_refinement_steps=5, 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 +
+
+      TwoPhaseFlowProblem<2> two_phase_flow_problem(0);
+
+
+in the main function to +
+
+      TwoPhaseFlowProblem<3> two_phase_flow_problem(0);
+
+
+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. + + + +

Possibilities for extensions

+ + +

Solvers

+ +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. + + +

Adaptivity

+ +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. -- 2.39.5