/* 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>
InverseMatrix<Matrix,Preconditioner>::
InverseMatrix (const Matrix &m,
const Preconditioner &preconditioner)
- :
- matrix (&m),
- preconditioner (preconditioner)
+ :
+ matrix (&m),
+ preconditioner (preconditioner)
{}
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()));
}
}
BlockSchurPreconditioner (
const TrilinosWrappers::BlockSparseMatrix &S,
const InverseMatrix<TrilinosWrappers::SparseMatrix,
- PreconditionerMp> &Mpinv,
+ PreconditionerMp> &Mpinv,
const PreconditionerA &Apreconditioner);
void vmult (TrilinosWrappers::BlockVector &dst,
BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S,
const InverseMatrix<TrilinosWrappers::SparseMatrix,
- PreconditionerMp> &Mpinv,
+ PreconditionerMp> &Mpinv,
const PreconditionerA &Apreconditioner)
:
darcy_matrix (&S),
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));
}
+ // @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
{
};
+ // @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>
{
}
+ // @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>
}
+ // @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)
{
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)
:
{}
+ // @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 ()
{
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 ()
}
}
+
+ // @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 ()
}
+
+ // @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 ()
{
}
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;
}
}
+
+ // @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 ()
{
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 ()
{
}
+
+ // @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 ()
{
}
}
+
+
+ // @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>::
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,
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>::
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);
}
}
+ // @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 ()
{
{
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(),
}
+
+
+ // @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
}
}
+
+
+ // @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>::
// 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>::
(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();
}
}
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;
}
+
+ // @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
{
}
+
+ // @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 ()
}
+ // @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 ()
{
{
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);
}
projects / dealii-svn.git / commitdiff