From b03f5f785e90502e5039c728ba63c9c1a6d9eb21 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 24 Nov 2010 02:54:00 +0000 Subject: [PATCH] Add Chih-Che's comments. git-svn-id: https://svn.dealii.org/trunk@22846 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-43/step-43.cc | 686 +++++++++++++++++++++++++++- 1 file changed, 661 insertions(+), 25 deletions(-) diff --git a/deal.II/examples/step-43/step-43.cc b/deal.II/examples/step-43/step-43.cc index de380e660a..1e1f912978 100644 --- a/deal.II/examples/step-43/step-43.cc +++ b/deal.II/examples/step-43/step-43.cc @@ -12,6 +12,33 @@ /* 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 #include #include @@ -54,6 +81,19 @@ 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 @@ -78,9 +118,9 @@ namespace LinearSolvers InverseMatrix:: 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 &Mpinv, + PreconditionerMp> &Mpinv, const PreconditionerA &Apreconditioner); void vmult (TrilinosWrappers::BlockVector &dst, @@ -135,7 +175,7 @@ namespace LinearSolvers BlockSchurPreconditioner:: BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, const InverseMatrix &Mpinv, + PreconditionerMp> &Mpinv, const PreconditionerA &Apreconditioner) : darcy_matrix (&S), @@ -147,8 +187,8 @@ namespace LinearSolvers template void BlockSchurPreconditioner::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 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 class PressureRightHandSide : public Function { @@ -362,6 +467,17 @@ SaturationInitialValues::vector_value (const Point &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 @@ -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::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 TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) : @@ -534,6 +681,63 @@ TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) {} + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::setup_dofs () { @@ -686,6 +890,53 @@ void TwoPhaseFlowProblem::setup_dofs () saturation_rhs.reinit (n_s); } + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::assemble_darcy_preconditioner () @@ -775,6 +1026,32 @@ TwoPhaseFlowProblem::assemble_darcy_preconditioner () } } + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::build_darcy_preconditioner () @@ -791,6 +1068,128 @@ TwoPhaseFlowProblem::build_darcy_preconditioner () } + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::assemble_darcy_system () { @@ -881,9 +1280,9 @@ void TwoPhaseFlowProblem::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::faces_per_cell; @@ -924,6 +1323,23 @@ void TwoPhaseFlowProblem::assemble_darcy_system () } } + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::assemble_saturation_system () { @@ -937,6 +1353,22 @@ void TwoPhaseFlowProblem::assemble_saturation_system () assemble_saturation_rhs (); } + + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::assemble_saturation_matrix () { @@ -983,6 +1415,57 @@ void TwoPhaseFlowProblem::assemble_saturation_matrix () } + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::assemble_saturation_rhs () { @@ -1046,6 +1529,19 @@ void TwoPhaseFlowProblem::assemble_saturation_rhs () } } + + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem:: @@ -1107,8 +1603,8 @@ assemble_saturation_rhs_cell_term (const FEValues &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 &saturation_f saturation_rhs); } + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem:: @@ -1155,11 +1661,11 @@ assemble_saturation_rhs_boundary_term (const FEFaceValues &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 &satu } + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::solve () { @@ -1183,11 +1739,11 @@ void TwoPhaseFlowProblem::solve () { const LinearSolvers::InverseMatrix + TrilinosWrappers::PreconditionIC> mp_inverse (darcy_preconditioner_matrix.block(1,1), *Mp_preconditioner); const LinearSolvers::BlockSchurPreconditioner + TrilinosWrappers::PreconditionIC> preconditioner (darcy_matrix, mp_inverse, *Amg_preconditioner); SolverControl solver_control (darcy_matrix.m(), @@ -1290,6 +1846,33 @@ void TwoPhaseFlowProblem::solve () } + + + // @sect3{TwoPhaseFlowProblem::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 bool TwoPhaseFlowProblem::determine_whether_to_solve_pressure_velocity_part () const @@ -1358,6 +1941,16 @@ TwoPhaseFlowProblem::determine_whether_to_solve_pressure_velocity_part () c } } + + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem:: @@ -1389,6 +1982,20 @@ compute_refinement_indicators (Vector &refinement_indicators) const // std::cout << "max_refinement_indicator =" << max_refinement_indicator << std::endl; } + + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem:: @@ -1411,9 +2018,9 @@ refine_grid (const Vector &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 &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 &refinement_indicators) } + + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::output_results () const { @@ -1620,6 +2238,12 @@ void TwoPhaseFlowProblem::output_results () const } + + // @sect3{TwoPhaseFlowProblem::THE_REMAINING_FUNCTIONS} + + // The remaining functions that have been + // used in step-31 so we don't have to + // describe their implementations. template void TwoPhaseFlowProblem::project_back_saturation () @@ -1803,6 +2427,18 @@ compute_viscosity (const std::vector &old_saturation, } + // @sect3{TwoPhaseFlowProblem::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 void TwoPhaseFlowProblem::run () { @@ -1852,7 +2488,7 @@ void TwoPhaseFlowProblem::run () { predictor_saturation_solution = saturation_solution; predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution); - Vector refinement_indicators (triangulation.n_active_cells()); + Vector refinement_indicators (triangulation.n_active_cells()); compute_refinement_indicators(refinement_indicators); refine_grid(refinement_indicators); } -- 2.39.5