]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add Chih-Che's comments.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Nov 2010 02:54:00 +0000 (02:54 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Nov 2010 02:54:00 +0000 (02:54 +0000)
git-svn-id: https://svn.dealii.org/trunk@22846 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-43/step-43.cc

index de380e660af2d28b9a33768033c081b22b60950b..1e1f91297837b14d9de6ba1ea476a3febc10b687 100644 (file)
 /*    further information on this license.                        */
 
 
+                                // @sect3{Include files}
+
+                                // The first step, as always, is to include
+                                // the functionality of these well-known
+                                // deal.II library files and some C++ header
+                                // files.
+                                // 
+                                // In this program, we use a tensor-valued
+                                // coefficient. Since it may have a spatial
+                                // dependence, we consider it a tensor-valued
+                                // function. The following include file
+                                // provides the TensorFunction class that
+                                // offers such functionality:
+                                // 
+                                // Then we need to include some header files
+                                // that provide vector, matrix, and
+                                // preconditioner classes that implement
+                                // interfaces to the respective Trilinos
+                                // classes, which has been used in
+                                // step-31. In particular, we will need
+                                // interfaces to the matrix and vector
+                                // classes based on Trilinos as well as
+                                // Trilinos preconditioners:
+                                // 
+                                // At the end of this top-matter, we import
+                                // all deal.II names into the global
+                                // namespace:
 #include <base/quadrature_lib.h>
 #include <base/logstream.h>
 #include <base/utilities.h>
 
 using namespace dealii;
 
+
+                                // @sect3{The InverseMatrix class template}
+
+                                // This part is exactly the same as that used in step-31.
+
+                                // @sect3{Schur complement preconditioner}
+
+                                // This part for the Schur complement
+                                // preconditioner is almost the same as that
+                                // used in step-31. The only difference is
+                                // that the original variable name
+                                // stokes_matrix is replaced by another name
+                                // darcy_matrix to satisfy our problem.
 namespace LinearSolvers
 {
   template <class Matrix, class Preconditioner>
@@ -78,9 +118,9 @@ namespace LinearSolvers
   InverseMatrix<Matrix,Preconditioner>::
   InverseMatrix (const Matrix &m,
                  const Preconditioner &preconditioner)
-                 :
-                 matrix (&m),
-                 preconditioner (preconditioner)
+                 :
+                 matrix (&m),
+                 preconditioner (preconditioner)
   {}
 
 
@@ -99,11 +139,11 @@ namespace LinearSolvers
 
     try
       {
-       cg.solve (*matrix, dst, src, preconditioner);
+       cg.solve (*matrix, dst, src, preconditioner);
       }
     catch (std::exception &e)
       {
-       Assert (false, ExcMessage(e.what()));
+       Assert (false, ExcMessage(e.what()));
       }
   }
 
@@ -114,7 +154,7 @@ namespace LinearSolvers
       BlockSchurPreconditioner (
         const TrilinosWrappers::BlockSparseMatrix     &S,
         const InverseMatrix<TrilinosWrappers::SparseMatrix,
-                            PreconditionerMp>         &Mpinv,
+       PreconditionerMp>         &Mpinv,
         const PreconditionerA                         &Apreconditioner);
 
       void vmult (TrilinosWrappers::BlockVector       &dst,
@@ -135,7 +175,7 @@ namespace LinearSolvers
   BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
   BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix  &S,
                            const InverseMatrix<TrilinosWrappers::SparseMatrix,
-                                               PreconditionerMp>      &Mpinv,
+                          PreconditionerMp>      &Mpinv,
                            const PreconditionerA                      &Apreconditioner)
                   :
                   darcy_matrix            (&S),
@@ -147,8 +187,8 @@ namespace LinearSolvers
 
   template <class PreconditionerA, class PreconditionerMp>
   void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
-     TrilinosWrappers::BlockVector       &dst,
-     const TrilinosWrappers::BlockVector &src) const
+    TrilinosWrappers::BlockVector       &dst,
+    const TrilinosWrappers::BlockVector &src) const
   {
     a_preconditioner.vmult (dst.block(0), src.block(0));
     darcy_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
@@ -158,6 +198,66 @@ namespace LinearSolvers
 }
 
 
+                                // @sect3{The TwoPhaseFlowProblem class}
+
+                                // The definition of the class that defines
+                                // the top-level logic of solving the
+                                // time-dependent advection-dominated
+                                // two-phase flow problem (or
+                                // Buckley-Leverett problem
+                                // [Buckley 1942]) is mainly based on
+                                // three tutorial programs (step-21, step-31,
+                                // step-33). The main difference is that,
+                                // since adaptive operator splitting is
+                                // considered, we need a bool-type variable
+                                // solve_pressure_velocity_part to tell us
+                                // when we need to solve the pressure and
+                                // velocity part, need another bool-type
+                                // variable
+                                // previous_solve_pressure_velocity_part to
+                                // determine if we have to cumulate
+                                // micro-time steps that we need them to do
+                                // extrapolation for the total velocity, and
+                                // some solution vectors
+                                // (e.g. nth_darcy_solution_after_solving_pressure_part
+                                // and
+                                // n_minus_oneth_darcy_solution_after_solving_pressure_part)
+                                // to store some solutions in previous time
+                                // steps after the solution of the pressure
+                                // and velocity part.
+                                // 
+                                // The member functions within this class
+                                // have been named so properly so that
+                                // readers can easily understand what they
+                                // are doing.
+                                // 
+                                // Like step-31, this tutorial uses two
+                                // DoFHandler objects for the darcy system
+                                // (presure and velocity) and
+                                // saturation. This is because we want it to
+                                // run faster, which reasons have been
+                                // described in step-31.
+                                // 
+                                // There is yet another important thing:
+                                // unlike step-31. this step uses one more
+                                // ConstraintMatrix object called
+                                // darcy_preconditioner_constraints. This
+                                // constraint object only for assembling the
+                                // matrix for darcy preconditioner includes
+                                // hanging node constrants as well as
+                                // Dirichlet boundary value
+                                // constraints. Without this constraint
+                                // object for the preconditioner, we cannot
+                                // get the convergence results when we solve
+                                // darcy linear system.
+                                // 
+                                // The last one variable indicates whether
+                                // the matrix needs to be rebuilt the next
+                                // time the corresponding build functions are
+                                // called. This allows us to move the
+                                // corresponding if into the function and
+                                // thereby keeping our main run() function
+                                // clean and easy to read.
 template <int dim>
 class TwoPhaseFlowProblem
 {
@@ -263,6 +363,11 @@ class TwoPhaseFlowProblem
 };
 
 
+                                // @sect3{Pressure right hand side, Pressure boundary values and saturation initial value classes}
+
+                                // This part is directly taken from step-21
+                                // so there is no need to repeat the same
+                                // descriptions.
 template <int dim>
 class PressureRightHandSide : public Function<dim>
 {
@@ -362,6 +467,17 @@ SaturationInitialValues<dim>::vector_value (const Point<dim> &p,
 }
 
 
+                                // @sect3{Permeability models}
+
+                                // In this tutorial, we still use two
+                                // permeability models previous used in
+                                // step-21 so we refrain from excessive
+                                // comments about them. But we want to note
+                                // that if ones use the Random Medium model,
+                                // they can change one parameter called the
+                                // number of high-permeability regions/points
+                                // to increase the amount of permeability in
+                                // the computational domain.
 namespace SingleCurvingCrack
 {
   template <int dim>
@@ -479,6 +595,13 @@ namespace RandomMedium
 }
 
 
+                                // @sect3{Physical quantities}
+
+                                // The implementations of all the physical
+                                // quantities such as total mobility
+                                // $\lambda_t$ and fractional flow of water
+                                // $F$ are taken from step-21 so again we
+                                // don't have do any comment about them.
 double mobility_inverse (const double S,
                          const double viscosity)
 {
@@ -506,6 +629,30 @@ double get_fractional_flow_derivative (const double S,
   return numerator / denomerator;
 }
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem}
+
+                                // The constructor of this class is an
+                                // extension of the constructor in step-21
+                                // and step-31. We need to add the various
+                                // variables that concern the saturation. As
+                                // discussed in the introduction, we are
+                                // going to use $Q_2 \times Q_1$
+                                // (Taylor-Hood) elements again for the darcy
+                                // system, which element combination fulfills
+                                // the Ladyzhenskaya-Babuska-Brezzi (LBB)
+                                // conditions
+                                // [Brezzi and Fortin 1991, Chen 2005], and $Q_1$
+                                // elements for the saturation. However, by
+                                // using variables that store the polynomial
+                                // degree of the darcy and temperature finite
+                                // elements, it is easy to consistently
+                                // modify the degree of the elements as well
+                                // as all quadrature formulas used on them
+                                // downstream. Moreover, we initialize the
+                                // time stepping, variables related to
+                                // operator splitting as well as the option
+                                // for matrix assembly and preconditioning:
 template <int dim>
 TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem (const unsigned int degree)
                 :
@@ -534,6 +681,63 @@ TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem (const unsigned int degree)
 {}
 
 
+                                // @sect3{TwoPhaseFlowProblem<dim>::setup_dofs}
+
+                                // This is the function that sets up the
+                                // DoFHandler objects we have here (one for
+                                // the darcy part and one for the saturation
+                                // part) as well as set to the right sizes
+                                // the various objects required for the
+                                // linear algebra in this program. Its basic
+                                // operations are similar to what authors in
+                                // step-31 did.
+                                // 
+                                // The body of the function first enumerates
+                                // all degrees of freedom for the darcy and
+                                // saturation systems. For the darcy part,
+                                // degrees of freedom are then sorted to
+                                // ensure that velocities precede pressure
+                                // DoFs so that we can partition the darcy
+                                // matrix into a $2 \times 2$ matrix. Like
+                                // step-31, the present step does not perform
+                                // any additional DoF renumbering.
+                                // 
+                                // Then, we need to incorporate hanging node
+                                // constraints and Dirichlet boundary value
+                                // constraints into
+                                // darcy_preconditioner_constraints. However,
+                                // this constraints are only set to the
+                                // pressure component since the Schur
+                                // complement preconditioner that corresponds
+                                // to the porous media flow operator in
+                                // non-mixed form, $-\nabla \cdot [\mathbf K
+                                // \lambda_t(S)]\nabla$. Therefore, we use a
+                                // component_mask that filters out the
+                                // velocity component, so that the
+                                // condensation is performed on pressure
+                                // degrees of freedom only.
+                                // 
+                                // After having done so, we count the number
+                                // of degrees of freedom in the various
+                                // blocks:
+                                // 
+                                // The next step is to create the sparsity
+                                // pattern for the darcy and saturation
+                                // system matrices as well as the
+                                // preconditioner matrix from which we build
+                                // the darcy preconditioner. As in step-31,
+                                // we choose to create the pattern not as in
+                                // the first few tutorial programs, but by
+                                // using the blocked version of
+                                // CompressedSimpleSparsityPattern. The
+                                // reason for doing this is mainly memory,
+                                // that is, the SparsityPattern class would
+                                // consume too much memory when used in three
+                                // spatial dimensions as we intend to do for
+                                // this program. So, for this, we follow the
+                                // same way as step-31 did and we don't have
+                                // to repeat descriptions again for the rest
+                                // of the member function.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::setup_dofs ()
 {
@@ -686,6 +890,53 @@ void TwoPhaseFlowProblem<dim>::setup_dofs ()
   saturation_rhs.reinit (n_s);
 }
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner}
+
+                                // This function assembles the matrix we use
+                                // for preconditioning the darcy system. What
+                                // we need are a vector matrix weighted by
+                                // $\left(\mathbf{K} \lambda_t\right)^{-1}$
+                                // on the velocity components and a mass
+                                // matrix weighted by $\left(\mathbf{K}
+                                // \lambda_t\right)$ on the pressure
+                                // component. We start by generating a
+                                // quadrature object of appropriate order,
+                                // the FEValues object that can give values
+                                // and gradients at the quadrature points
+                                // (together with quadrature weights). Next
+                                // we create data structures for the cell
+                                // matrix and the relation between local and
+                                // global DoFs. The vectors phi_u and
+                                // grad_phi_p are going to hold the values of
+                                // the basis functions in order to faster
+                                // build up the local matrices, as was
+                                // already done in step-22. Before we start
+                                // the loop over all active cells, we have to
+                                // specify which components are pressure and
+                                // which are velocity.
+                                // 
+                                // The creation of the local matrix is rather
+                                // simple. There are only a term weighted by
+                                // $\left(\mathbf{K} \lambda_t\right)^{-1}$
+                                // (on the velocity) and a mass matrix
+                                // weighted by $\left(\mathbf{K}
+                                // \lambda_t\right)$ to be generated, so the
+                                // creation of the local matrix is done in
+                                // two lines. Once the local matrix is ready
+                                // (loop over rows and columns in the local
+                                // matrix on each quadrature point), we get
+                                // the local DoF indices and write the local
+                                // information into the global matrix. We do
+                                // this by directly applying the constraints
+                                // (i.e. darcy_preconditioner_constraints)
+                                // from hanging nodes locally and Dirichlet
+                                // boundary conditions with zero values. By
+                                // doing so, we don't have to do that
+                                // afterwards, and we don't also write into
+                                // entries of the matrix that will actually
+                                // be set to zero again later when
+                                // eliminating constraints.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner ()
@@ -775,6 +1026,32 @@ TwoPhaseFlowProblem<dim>::assemble_darcy_preconditioner ()
     }
 }
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::build_darcy_preconditioner}
+
+                                // This function generates the inner
+                                // preconditioners that are going to be used
+                                // for the Schur complement block
+                                // preconditioner. The preconditioners need
+                                // to be regenerated at every saturation time
+                                // step since they contain the independent
+                                // variables saturation $S$ with time.
+                                // 
+                                // Next, we set up the preconditioner for the
+                                // velocity-velocity matrix
+                                // $\mathbf{M}^{\mathbf{u}}$ and the Schur
+                                // complement $\mathbf{S}$. As explained in
+                                // the introduction, we are going to use an
+                                // IC preconditioner based on a vector matrix
+                                // (which is spectrally close to the darcy
+                                // matrix $\mathbf{M}^{\mathbf{u}}$) and
+                                // another based on a Laplace vector matrix
+                                // (which is spectrally close to the
+                                // non-mixed pressure matrix
+                                // $\mathbf{S}$). Usually, the
+                                // TrilinosWrappers::PreconditionIC class can
+                                // be seen as a good black-box preconditioner
+                                // which does not need any special knowledge.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::build_darcy_preconditioner ()
@@ -791,6 +1068,128 @@ TwoPhaseFlowProblem<dim>::build_darcy_preconditioner ()
 
 }
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::assemble_darcy_system}
+
+                                // This is the function that assembles the
+                                // linear system for the darcy system.
+                                // 
+                                // Regarding the technical details of
+                                // implementation, the procedures are similar
+                                // to those in step-22 and step-31 we reset
+                                // matrix and vector, create a quadrature
+                                // formula on the cells, and then create the
+                                // respective FEValues object. For the update
+                                // flags, we require basis function
+                                // derivatives only in case of a full
+                                // assembly, since they are not needed for
+                                // the right hand side; as always, choosing
+                                // the minimal set of flags depending on what
+                                // is currently needed makes the call to
+                                // FEValues::reinit further down in the
+                                // program more efficient.
+                                // 
+                                // There is one thing that needs to be
+                                // commented Â¡V since we have a separate
+                                // finite element and DoFHandler for the
+                                // saturation, we need to generate a second
+                                // FEValues object for the proper evaluation
+                                // of the saturation solution. This isn't too
+                                // complicated to realize here: just use the
+                                // saturation structures and set an update
+                                // flag for the basis function values which
+                                // we need for evaluation of the saturation
+                                // solution. The only important part to
+                                // remember here is that the same quadrature
+                                // formula is used for both FEValues objects
+                                // to ensure that we get matching information
+                                // when we loop over the quadrature points of
+                                // the two objects.
+                                // 
+                                // The declarations proceed with some
+                                // shortcuts for array sizes, the creation of
+                                // the local matrix, right hand side as well
+                                // as the vector for the indices of the local
+                                // dofs compared to the global system.
+                                // 
+                                // Note that in its present form, the
+                                // function uses the permeability implemented
+                                // in the RandomMedium::KInverse
+                                // class. Switching to the single curved
+                                // crack permeability function is as simple
+                                // as just changing the namespace name.
+                                // 
+                                // Here's the an important step: we have to
+                                // get the values of the saturation function
+                                // of the previous time step at the
+                                // quadrature points. To this end, we can use
+                                // the FEValues::get_function_values
+                                // (previously already used in step-9,
+                                // step-14 and step-15), a function that
+                                // takes a solution vector and returns a list
+                                // of function values at the quadrature
+                                // points of the present cell. In fact, it
+                                // returns the complete vector-valued
+                                // solution at each quadrature point,
+                                // i.e. not only the saturation but also the
+                                // velocities and pressure:
+                                // 
+                                // Next we need a vector that will contain
+                                // the values of the saturation solution at
+                                // the previous time level at the quadrature
+                                // points to assemble the source term in the
+                                // right hand side of the momentum
+                                // equation. Let's call this vector
+                                // old_saturation_values.
+                                // 
+                                // The set of vectors we create next hold the
+                                // evaluations of the basis functions as well
+                                // as their gradients and symmetrized
+                                // gradients that will be used for creating
+                                // the matrices. Putting these into their own
+                                // arrays rather than asking the FEValues
+                                // object for this information each time it
+                                // is needed is an optimization to accelerate
+                                // the assembly process, see step-22 for
+                                // details.
+                                // 
+                                // The last two declarations are used to
+                                // extract the individual blocks (velocity,
+                                // pressure, saturation) from the total FE
+                                // system.
+                                // 
+                                // Now start the loop over all cells in the
+                                // problem. We are working on two different
+                                // DoFHandlers for this assembly routine, so
+                                // we must have two different cell iterators
+                                // for the two objects in use. This might
+                                // seem a bit peculiar, since both the darcy
+                                // system and the saturation system use the
+                                // same grid, but that's the only way to keep
+                                // degrees of freedom in sync. The first
+                                // statements within the loop are again all
+                                // very familiar, doing the update of the
+                                // finite element data as specified by the
+                                // update flags, zeroing out the local arrays
+                                // and getting the values of the old solution
+                                // at the quadrature points. Then we are
+                                // ready to loop over the quadrature points
+                                // on the cell.
+                                // 
+                                // Once this is done, we start the loop over
+                                // the rows and columns of the local matrix
+                                // and feed the matrix with the relevant
+                                // products.
+                                // 
+                                // The last step in the loop over all cells
+                                // is to enter the local contributions into
+                                // the global matrix and vector structures to
+                                // the positions specified in
+                                // local_dof_indices. Again, we let the
+                                // ConstraintMatrix class do the insertion of
+                                // the cell matrix elements to the global
+                                // matrix, which already condenses the
+                                // hanging node constraints.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::assemble_darcy_system ()
 {
@@ -881,9 +1280,9 @@ void TwoPhaseFlowProblem<dim>::assemble_darcy_system ()
                 }
 
               local_rhs(i) += (-phi_p[i] * pressure_rhs_values[q])*
-                               darcy_fe_values.JxW(q);
+                             darcy_fe_values.JxW(q);
             }
-         }
+       }
 
       for (unsigned int face_no=0;
            face_no<GeometryInfo<dim>::faces_per_cell;
@@ -924,6 +1323,23 @@ void TwoPhaseFlowProblem<dim>::assemble_darcy_system ()
     }
 }
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::assemble_saturation_system}
+
+                                // This function is to assemble the linear
+                                // system for the saturation transport
+                                // equation. It includes two member
+                                // functions: assemble_saturation_matrix ()
+                                // and assemble_saturation_rhs (). The former
+                                // function that assembles the saturation
+                                // left hand side needs to be changed only
+                                // when grids have been changed since the
+                                // matrix is filled only with basis
+                                // functions. However, the latter that
+                                // assembles the right hand side must be
+                                // changed at every saturation time step
+                                // since it depends on an unknown variable
+                                // saturation.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::assemble_saturation_system ()
 {
@@ -937,6 +1353,22 @@ void TwoPhaseFlowProblem<dim>::assemble_saturation_system ()
   assemble_saturation_rhs ();
 }
 
+
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::assemble_saturation_matrix}
+
+                                // This function is easily understood since
+                                // it only forms a simple mass matrix for the
+                                // left hand side of the saturation linear
+                                // system by basis functions phi_i_s and
+                                // phi_j_s only. Finally, as usual, we enter
+                                // the local contribution into the global
+                                // matrix by specifying the position in
+                                // local_dof_indices. This is done by letting
+                                // the ConstraintMatrix class do the
+                                // insertion of the cell matrix elements to
+                                // the global matrix, which already condenses
+                                // the hanging node constraints.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::assemble_saturation_matrix ()
 {
@@ -983,6 +1415,57 @@ void TwoPhaseFlowProblem<dim>::assemble_saturation_matrix ()
 }
 
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::assemble_saturation_rhs}
+
+                                // This function is to assemble the right
+                                // hand side of the saturation transport
+                                // equation. Before assembling it, we have to
+                                // call two FEValues objects for the darcy
+                                // and saturation systems respectively and,
+                                // even more, two FEFaceValues objects for
+                                // the both systems because we have a
+                                // boundary integral term in the weak form of
+                                // saturation equation. For the FEFaceValues
+                                // object of the saturation system, we also
+                                // enter the normal vectors with an update
+                                // flag update_normal_vectors.
+                                // 
+                                // Next, before looping over all the cells,
+                                // we have to compute some parameters
+                                // (e.g. global_u_infty, global_S_variasion,
+                                // and global_Omega_diameter) that the
+                                // artificial viscosity $\nu$ needs, which
+                                // desriptions have been appearing in
+                                // step-31.
+                                // 
+                                // Next, we start to loop over all the
+                                // saturation and darcy cells to put the
+                                // local contributions into the global
+                                // vector. In this loop, in order to simplify
+                                // the implementation in this function, we
+                                // generate two more functions: one is
+                                // assemble_saturation_rhs_cell_term and the
+                                // other is
+                                // assemble_saturation_rhs_boundary_term,
+                                // which is contained in an inner boudary
+                                // loop. The former is to assemble the
+                                // integral cell term with neccessary
+                                // arguments and the latter is to assemble
+                                // the integral global boundary $\Omega$
+                                // terms. It should be noted that we achieve
+                                // the insertion of the cell or boundary
+                                // vector elements to the global vector in
+                                // the two functions rather than in this
+                                // present function by giving these two
+                                // functions with a common argument
+                                // local_dof_indices, and two arguments
+                                // saturation_fe_values darcy_fe_values for
+                                // assemble_saturation_rhs_cell_term and
+                                // another two arguments
+                                // saturation_fe_face_values
+                                // darcy_fe_face_values for
+                                // assemble_saturation_rhs_boundary_term.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs ()
 {
@@ -1046,6 +1529,19 @@ void TwoPhaseFlowProblem<dim>::assemble_saturation_rhs ()
     }
 }
 
+
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_cell_term}
+
+                                // In this function, we actually compute
+                                // every artificial viscosity for every
+                                // element. Then, with the artificial value,
+                                // we can finish assembling the saturation
+                                // right hand side cell integral
+                                // terms. Finally, we can pass the local
+                                // contributions on to the global vector with
+                                // the position specified in
+                                // local_dof_indices.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::
@@ -1107,8 +1603,8 @@ assemble_saturation_rhs_cell_term (const FEValues<dim>             &saturation_f
                          old_grad_saturation_solution_values[q] * grad_phi_i_s
                          +
                          old_s * phi_i_s)
-                         *
-                         saturation_fe_values.JxW(q);
+                       *
+                       saturation_fe_values.JxW(q);
       }
 
   saturation_constraints.distribute_local_to_global (local_rhs,
@@ -1116,6 +1612,16 @@ assemble_saturation_rhs_cell_term (const FEValues<dim>             &saturation_f
                                                      saturation_rhs);
 }
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::assemble_saturation_rhs_boundary_term}
+
+                                // In this function, we have to give
+                                // upwinding in the global boundary faces,
+                                // i.e. we impose the Dirichlet boundary
+                                // conditions only on inflow parts of global
+                                // boundary, which has been described in
+                                // step-21 so we refrain from giving more
+                                // descriptions about that.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::
@@ -1155,11 +1661,11 @@ assemble_saturation_rhs_boundary_term (const FEFaceValues<dim>             &satu
         local_rhs(i) -= time_step *
                         normal_flux *
                         f_saturation((is_outflow_q_point == true
-                        ?
-                        old_saturation_solution_values_face[q]
-                        :
-                        neighbor_saturation[q]),
-                        viscosity) *
+                                     ?
+                                     old_saturation_solution_values_face[q]
+                                     :
+                                     neighbor_saturation[q]),
+                                    viscosity) *
                         saturation_fe_face_values.shape_value (i,q) *
                         saturation_fe_face_values.JxW(q);
     }
@@ -1169,6 +1675,56 @@ assemble_saturation_rhs_boundary_term (const FEFaceValues<dim>             &satu
 }
 
 
+                                // @sect3{TwoPhaseFlowProblem<dim>::solve}
+
+                                // This function is to implement the operator
+                                // splitting algorithm. At the beginning of
+                                // the implementation, we decide whther to
+                                // solve the pressure-velocity part by
+                                // running an a posteriori criterion, which
+                                // will be described in the following
+                                // function. If we get the bool variable true
+                                // from that function, we will solve the
+                                // pressure-velocity part for updated
+                                // velocity. Then, we use GMRES with the
+                                // Schur complement preconditioner to solve
+                                // this linear system, as is described in the
+                                // Introduction. After solving the velocity
+                                // and pressure, we need to keep the
+                                // solutions for linear extrapolations in the
+                                // future. It is noted that we always solve
+                                // the pressure-velocity part in the first
+                                // three micro time steps to ensure accuracy
+                                // at the beginning of computation, and to
+                                // provide starting data to linearly
+                                // extrapolate previously computed velocities
+                                // to the current time step.
+                                // 
+                                // On the other hand, if we get a false
+                                // variable from the criterion, we will
+                                // directly use linear extrapolation to
+                                // compute the updated velocity for the
+                                // solution of saturation later.
+                                // 
+                                // Next, like step-21, this program need to
+                                // compute the present time step.
+                                // 
+                                // Next, we need to use two bool variables
+                                // solve_pressure_velocity_part and
+                                // previous_solve_pressure_velocity_part to
+                                // decide whether we stop or continue
+                                // cumulating the micro time steps for linear
+                                // extropolations in the next iteration. With
+                                // the reason, we need one variable
+                                // cumulative_nth_time_step for keeping the
+                                // present aggregated micro time steps and
+                                // anther one n_minus_oneth_time_step for
+                                // retaining the previous micro time steps.
+                                // 
+                                // Finally, we start to calculate the
+                                // saturation part with the use of the
+                                // incomplete Cholesky decomposition for
+                                // preconditioning.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::solve ()
 {
@@ -1183,11 +1739,11 @@ void TwoPhaseFlowProblem<dim>::solve ()
 
       {
         const LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
-                                           TrilinosWrappers::PreconditionIC>
+         TrilinosWrappers::PreconditionIC>
           mp_inverse (darcy_preconditioner_matrix.block(1,1), *Mp_preconditioner);
 
         const LinearSolvers::BlockSchurPreconditioner<TrilinosWrappers::PreconditionIC,
-                                                      TrilinosWrappers::PreconditionIC>
+         TrilinosWrappers::PreconditionIC>
           preconditioner (darcy_matrix, mp_inverse, *Amg_preconditioner);
 
         SolverControl solver_control (darcy_matrix.m(),
@@ -1290,6 +1846,33 @@ void TwoPhaseFlowProblem<dim>::solve ()
 
 }
 
+
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::determine_whether_to_solve_pressure_velocity_part}
+
+                                // This function is to implement an a
+                                // posteriori criterion in
+                                // \eqref{eq:recompute-criterion} for
+                                // adaptive operator splitting. As mentioned
+                                // in step-31, we use two FEValues objects
+                                // initialized with two cell iterators that
+                                // we walk in parallel through the two
+                                // DoFHandler objects associated with the
+                                // same Triangulation object; for these two
+                                // FEValues objects, we use of course the
+                                // same quadrature objects so that we can
+                                // iterate over the same set of quadrature
+                                // points, but each FEValues object will get
+                                // update flags only according to what it
+                                // actually needs to compute.
+                                // 
+                                // In addition to this, if someone doesn't
+                                // want to perform their simulation with
+                                // operator splitting, they can lower the
+                                // criterion value (default value is $5.0$)
+                                // down to zero ad therefore numerical
+                                // algorithm becomes the original IMPES
+                                // method.
 template <int dim>
 bool
 TwoPhaseFlowProblem<dim>::determine_whether_to_solve_pressure_velocity_part () const
@@ -1358,6 +1941,16 @@ TwoPhaseFlowProblem<dim>::determine_whether_to_solve_pressure_velocity_part () c
     }
 }
 
+
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::compute_refinement_indicators}
+
+                                // This function is to to compute the
+                                // refinement indicator in
+                                // \eqref{eq:refinement_indicator} for each
+                                // cell and its implementation is similar to
+                                // that contained in step-33. There is no
+                                // need to repeat descriptions about it.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::
@@ -1389,6 +1982,20 @@ compute_refinement_indicators (Vector<double> &refinement_indicators) const
 //  std::cout << "max_refinement_indicator =" << max_refinement_indicator << std::endl;
 }
 
+
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::refine_grid}
+
+                                // This function is to decide if every cell
+                                // is refined or coarsened with computed
+                                // refinement indicators in the previous
+                                // function and do the interpolations of the
+                                // solution vectors. The main difference from
+                                // the previous time-dependent tutorials is
+                                // that there is no need to do the solution
+                                // interpolations if we don't have any cell
+                                // that is refined or coarsend, saving some
+                                // additional computing time.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::
@@ -1411,9 +2018,9 @@ refine_grid (const Vector<double> &refinement_indicators)
            (std::fabs(refinement_indicators(cell_no)) > saturation_value))
          cell->set_refine_flag();
        else
-           if ((cell->level() > double(n_refinement_steps)) &&
-               (std::fabs(refinement_indicators(cell_no)) < 0.75 * saturation_value))
-             cell->set_coarsen_flag();
+         if ((cell->level() > double(n_refinement_steps)) &&
+             (std::fabs(refinement_indicators(cell_no)) < 0.75 * saturation_value))
+           cell->set_coarsen_flag();
       }
   }
 
@@ -1429,7 +2036,7 @@ refine_grid (const Vector<double> &refinement_indicators)
 
     for (; cell!=endc; ++cell)
       if (cell->refine_flag_set())
-           ++number_of_cells_refine;
+       ++number_of_cells_refine;
       else
        if (cell->coarsen_flag_set())
          ++number_of_cells_coarsen;
@@ -1519,6 +2126,17 @@ refine_grid (const Vector<double> &refinement_indicators)
 }
 
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::output_results}
+
+                                // This function to process the output
+                                // data. We only store the results when we
+                                // actually solve the pressure and velocity
+                                // part at the present time step. The rest of
+                                // the implementation is similar to that
+                                // output function in step-31, which
+                                // implementations has been explained in that
+                                // tutorial.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::output_results ()  const
 {
@@ -1620,6 +2238,12 @@ void TwoPhaseFlowProblem<dim>::output_results ()  const
 }
 
 
+
+                                // @sect3{TwoPhaseFlowProblem<dim>::THE_REMAINING_FUNCTIONS}
+
+                                // The remaining functions that have been
+                                // used in step-31 so we don't have to
+                                // describe their implementations.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::project_back_saturation ()
@@ -1803,6 +2427,18 @@ compute_viscosity (const std::vector<double>          &old_saturation,
 }
 
 
+                                // @sect3{TwoPhaseFlowProblem<dim>::run}
+
+                                // In this function, we follow the structure
+                                // of the same function partly in step-21 and
+                                // partly in step-31 so again there is no
+                                // need to repeat it. However, since we
+                                // consider the simulation with grid
+                                // adaptivity, we need to compute a
+                                // saturation predictor, which implementation
+                                // was first used in step-33, for the
+                                // function that computes the refinement
+                                // indicators.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::run ()
 {
@@ -1852,7 +2488,7 @@ void TwoPhaseFlowProblem<dim>::run ()
         {
           predictor_saturation_solution = saturation_solution;
           predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution);
-          Vector<double> refinement_indicators (triangulation.n_active_cells());
+         Vector<double> refinement_indicators (triangulation.n_active_cells());
           compute_refinement_indicators(refinement_indicators);
           refine_grid(refinement_indicators);
         }

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.