From 83d8d298eede36995b59292e973519bb10e1557d Mon Sep 17 00:00:00 2001 From: bangerth Date: Sat, 28 Oct 2006 19:26:18 +0000 Subject: [PATCH] Finish documenting git-svn-id: https://svn.dealii.org/trunk@14119 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-21/step-21.cc | 1071 ++++++++++++++------------- 1 file changed, 566 insertions(+), 505 deletions(-) diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc index 7fd03abb90..981b069468 100644 --- a/deal.II/examples/step-21/step-21.cc +++ b/deal.II/examples/step-21/step-21.cc @@ -11,18 +11,18 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ - // This program is an adaptation of step-20 - // and includes some technique of DG method - // from step-12. A good part of the program - // is therefore very similar to step-20 and - // we will not comment again on these - // parts. Only the new stuff will be - // discussed in more detail. - - // @sect3{Include files} - - // All of these include files have been used - // before: + // This program is an adaptation of step-20 + // and includes some technique of DG method + // from step-12. A good part of the program + // is therefore very similar to step-20 and + // we will not comment again on these + // parts. Only the new stuff will be + // discussed in more detail. + + // @sect3{Include files} + + // All of these include files have been used + // before: #include #include #include @@ -57,50 +57,50 @@ #include #include - // 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: + // 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: #include - // The last step is as in all - // previous programs: + // The last step is as in all + // previous programs: using namespace dealii; // @sect3{The TwoPhaseFlowProblem class} - // This is the main class of the program. It - // is close to the one of step-20, but with a - // few additional functions: - // - //
    - //
  • assemble_rhs_S assembles the - // right hand side of the saturation - // equation. As explained in the - // introduction, this can't be integrated - // into assemble_rhs since it depends - // on the velocity that is computed in the - // first part of the time step. - // - //
  • get_maximal_velocity does as its - // name suggests. This function is used in - // the computation of the time step size. - // - //
  • project_back_saturation resets - // all saturation degrees of freedom with - // values less than zero to zero, and all - // those with saturations greater than one - // to one. - //
- // - // The rest of the class should be pretty - // much obvious. The viscosity variable - // stores the viscosity $\mu$ that enters - // several of the formulas in the nonlinear - // equations. + // This is the main class of the program. It + // is close to the one of step-20, but with a + // few additional functions: + // + //
    + //
  • assemble_rhs_S assembles the + // right hand side of the saturation + // equation. As explained in the + // introduction, this can't be integrated + // into assemble_rhs since it depends + // on the velocity that is computed in the + // first part of the time step. + // + //
  • get_maximal_velocity does as its + // name suggests. This function is used in + // the computation of the time step size. + // + //
  • project_back_saturation resets + // all saturation degrees of freedom with + // values less than zero to zero, and all + // those with saturations greater than one + // to one. + //
+ // + // The rest of the class should be pretty + // much obvious. The viscosity variable + // stores the viscosity $\mu$ that enters + // several of the formulas in the nonlinear + // equations. template class TwoPhaseFlowProblem { @@ -138,14 +138,14 @@ class TwoPhaseFlowProblem }; - // @sect3{Equation data} + // @sect3{Equation data} - // @sect4{Pressure right hand side} - // At present, the right hand side of the - // pressure equation is simply the zero - // function. However, the rest of the program - // is fully equipped to deal with anything - // else, if this is desired: + // @sect4{Pressure right hand side} + // At present, the right hand side of the + // pressure equation is simply the zero + // function. However, the rest of the program + // is fully equipped to deal with anything + // else, if this is desired: template class PressureRightHandSide : public Function { @@ -153,7 +153,7 @@ class PressureRightHandSide : public Function PressureRightHandSide () : Function(1) {}; virtual double value (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; }; @@ -161,16 +161,16 @@ class PressureRightHandSide : public Function template double PressureRightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { return 0; } - // @sect4{Pressure boundary values} - // The next are pressure boundary values. As - // mentioned in the introduction, we choose a - // linear pressure field: + // @sect4{Pressure boundary values} + // The next are pressure boundary values. As + // mentioned in the introduction, we choose a + // linear pressure field: template class PressureBoundaryValues : public Function { @@ -178,29 +178,29 @@ class PressureBoundaryValues : public Function PressureBoundaryValues () : Function(1) {}; virtual double value (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; }; template double PressureBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { return 1-p[0]; } - // @sect4{Saturation boundary values} + // @sect4{Saturation boundary values} - // Then we also need boundary values on the - // inflow portions of the boundary. The - // question whether something is an inflow - // part is decided when assembling the right - // hand side, we only have to provide a - // functional description of the boundary - // values. This is as explained in the - // introduction: + // Then we also need boundary values on the + // inflow portions of the boundary. The + // question whether something is an inflow + // part is decided when assembling the right + // hand side, we only have to provide a + // functional description of the boundary + // values. This is as explained in the + // introduction: template class SaturationBoundaryValues : public Function { @@ -208,7 +208,7 @@ class SaturationBoundaryValues : public Function SaturationBoundaryValues () : Function(1) {}; virtual double value (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; }; @@ -216,7 +216,7 @@ class SaturationBoundaryValues : public Function template double SaturationBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const + const unsigned int /*component*/) const { if (p[0] == 0) return 1; @@ -226,25 +226,25 @@ SaturationBoundaryValues::value (const Point &p, - // @sect4{Initial data} - - // Finally, we need initial data. In reality, - // we only need initial data for the - // saturation, but we are lazy, so we will - // later, before the first time step, simply - // interpolate the entire solution for the - // previous time step from a function that - // contains all vector components. - // - // We therefore simply create a function that - // returns zero in all components. We do that - // by simply forward every function to the - // ZeroFunction class. Why not use that right - // away in the places of this program where - // we presently use the InitialValues - // class? Because this way it is simpler to - // later go back and choose a different - // function for initial values. + // @sect4{Initial data} + + // Finally, we need initial data. In reality, + // we only need initial data for the + // saturation, but we are lazy, so we will + // later, before the first time step, simply + // interpolate the entire solution for the + // previous time step from a function that + // contains all vector components. + // + // We therefore simply create a function that + // returns zero in all components. We do that + // by simply forward every function to the + // ZeroFunction class. Why not use that right + // away in the places of this program where + // we presently use the InitialValues + // class? Because this way it is simpler to + // later go back and choose a different + // function for initial values. template class InitialValues : public Function { @@ -252,10 +252,10 @@ class InitialValues : public Function InitialValues () : Function(dim+2) {}; virtual double value (const Point &p, - const unsigned int component = 0) const; + const unsigned int component = 0) const; virtual void vector_value (const Point &p, - Vector &value) const; + Vector &value) const; }; @@ -263,7 +263,7 @@ class InitialValues : public Function template double InitialValues::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { return ZeroFunction(dim+2).value (p, component); } @@ -272,7 +272,7 @@ InitialValues::value (const Point &p, template void InitialValues::vector_value (const Point &p, - Vector &values) const + Vector &values) const { ZeroFunction(dim+2).vector_value (p, values); } @@ -280,23 +280,23 @@ InitialValues::vector_value (const Point &p, - // @sect3{The inverse permeability tensor} + // @sect3{The inverse permeability tensor} - // As announced in the introduction, we - // implement two different permeability - // tensor fields. Each of them we put into a - // namespace of its own, so that it will be - // easy later to replace use of one by the - // other in the code. + // As announced in the introduction, we + // implement two different permeability + // tensor fields. Each of them we put into a + // namespace of its own, so that it will be + // easy later to replace use of one by the + // other in the code. - // @sect4{Single curving crack permeability} + // @sect4{Single curving crack permeability} - // The first function for the permeability - // was the one that models a single curving - // crack. It was already used at the end of - // step-20, and its functional form is given - // in the introduction of the present - // tutorial program: + // The first function for the permeability + // was the one that models a single curving + // crack. It was already used at the end of + // step-20, and its functional form is given + // in the introduction of the present + // tutorial program: namespace SingleCurvingCrack { template @@ -304,88 +304,88 @@ namespace SingleCurvingCrack { public: virtual void value_list (const std::vector > &points, - std::vector > &values) const; + std::vector > &values) const; }; template void KInverse::value_list (const std::vector > &points, - std::vector > &values) const + std::vector > &values) const { Assert (points.size() == values.size(), - ExcDimensionMismatch (points.size(), values.size())); + ExcDimensionMismatch (points.size(), values.size())); for (unsigned int p=0; pget_centers that computes the list of - // center points when called. - // - // Note that this class works just fine in - // both 2d and 3d, with the only difference - // being that we use more points in 3d: by - // experimenting we find that we need more - // exponentials in 3d than in 2d (we have - // more ground to cover, after all, if we - // want to keep the distance between centers - // roughly equal), so we choose 40 in 2d and - // 100 in 3d. For any other dimension, the - // function does presently not know what to - // do so simply throws an exception - // indicating exactly this. + // @sect4{Random medium permeability} + + // This function does as announced in the + // introduction, i.e. it creates an overlay + // of exponentials at random places. There is + // one thing worth considering for this + // class. The issue centers around the + // problem that the class creates the centers + // of the exponentials using a random + // function. If we therefore created the + // centers each time we create an object of + // the present type, we would get a different + // list of centers each time. That's not what + // we expect from classes of this type: they + // should reliably represent the same + // function. + // + // The solution to this problem is to make + // the list of centers a static member + // variable of this class, i.e. there exists + // exactly one such variable for the entire + // program, rather than for each object of + // this type. That's exactly what we are + // going to do. + // + // The next problem, however, is that we need + // a way to initialize this variable. Since + // this variable is initialized at the + // beginning of the program, we can't use a + // regular member function for that since + // there may not be an object of this type + // around at the time. The C++ standard + // therefore says that only non-member and + // static member functions can be used to + // initialize a static variable. We use the + // latter possibility by defining a function + // get_centers that computes the list of + // center points when called. + // + // Note that this class works just fine in + // both 2d and 3d, with the only difference + // being that we use more points in 3d: by + // experimenting we find that we need more + // exponentials in 3d than in 2d (we have + // more ground to cover, after all, if we + // want to keep the distance between centers + // roughly equal), so we choose 40 in 2d and + // 100 in 3d. For any other dimension, the + // function does presently not know what to + // do so simply throws an exception + // indicating exactly this. namespace RandomMedium { template @@ -393,7 +393,7 @@ namespace RandomMedium { public: virtual void value_list (const std::vector > &points, - std::vector > &values) const; + std::vector > &values) const; private: static std::vector > centers; @@ -413,15 +413,15 @@ namespace RandomMedium KInverse::get_centers () { const unsigned int N = (dim == 2 ? - 40 : - (dim == 3 ? - 100 : - throw ExcNotImplemented())); + 40 : + (dim == 3 ? + 100 : + throw ExcNotImplemented())); std::vector > centers_list (N); for (unsigned int i=0; i(rand())/RAND_MAX; + centers_list[i][d] = static_cast(rand())/RAND_MAX; return centers_list; } @@ -431,46 +431,46 @@ namespace RandomMedium template void KInverse::value_list (const std::vector > &points, - std::vector > &values) const + std::vector > &values) const { Assert (points.size() == values.size(), - ExcDimensionMismatch (points.size(), values.size())); + ExcDimensionMismatch (points.size(), values.size())); for (unsigned int p=0; p Tensor<1,dim> extract_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) + const unsigned int i, + const unsigned int q) { Tensor<1,dim> tmp; @@ -506,8 +506,8 @@ extract_u (const FEValuesBase &fe_values, template double extract_div_u (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) + const unsigned int i, + const unsigned int q) { double divergence = 0; for (unsigned int d=0; d &fe_values, template Tensor<1,dim> extract_grad_s (const FEValuesBase &fe_values, - const unsigned int i, - const unsigned int q) + const unsigned int i, + const unsigned int q) { Tensor<1,dim> tmp; for (unsigned int d=0; d &fe_values, // @sect3{Linear solvers and preconditioners} - // The linear solvers we use are also - // completely analogous to the ones used in - // step-20. The following classes are - // therefore copied verbatim from there. + // The linear solvers we use are also + // completely analogous to the ones used in + // step-20. The following classes are + // therefore copied verbatim from there. template class InverseMatrix : public Subscriptor { @@ -617,7 +617,7 @@ class SchurComplement : public Subscriptor SchurComplement:: SchurComplement (const BlockSparseMatrix &A, - const InverseMatrix > &Minv) + const InverseMatrix > &Minv) : system_matrix (&A), m_inverse (&Minv), @@ -674,30 +674,30 @@ void ApproximateSchurComplement::vmult (Vector &dst, // @sect3{TwoPhaseFlowProblem class implementation} - // Here now the implementation of the main - // class. Much of it is actually copied from - // step-20, so we won't comment on it in much - // detail. You should try to get familiar - // with that program first, then most of what - // is happening here should be mostly clear. + // Here now the implementation of the main + // class. Much of it is actually copied from + // step-20, so we won't comment on it in much + // detail. You should try to get familiar + // with that program first, then most of what + // is happening here should be mostly clear. // @sect4{TwoPhaseFlowProblem::TwoPhaseFlowProblem} // First for the constructor. We use $RT_k - // \times DG_k \times DG_k$ spaces. The time + // \times DQ_k \times DQ_k$ spaces. The time // step is set to zero initially, but will be // computed before it is needed first, as // described in a subsection of the // introduction. template TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) - : - degree (degree), + : + degree (degree), fe (FE_RaviartThomas(degree), 1, FE_DGQ(degree), 1, - FE_DGQ(degree), 1), - dof_handler (triangulation), - n_refinement_steps (5), - time_step (0), + FE_DGQ(degree), 1), + dof_handler (triangulation), + n_refinement_steps (5), + time_step (0), viscosity (0.2) {} @@ -724,16 +724,16 @@ void TwoPhaseFlowProblem::make_grid_and_dofs () DoFTools::count_dofs_per_component (dof_handler, dofs_per_component); const unsigned int n_u = dofs_per_component[0], n_p = dofs_per_component[dim], - n_s = dofs_per_component[dim+1]; + n_s = dofs_per_component[dim+1]; std::cout << "Number of active cells: " - << triangulation.n_active_cells() - << std::endl + << triangulation.n_active_cells() + << std::endl << "Number of degrees of freedom: " - << dof_handler.n_dofs() + << dof_handler.n_dofs() << " (" << n_u << '+' << n_p << '+'<< n_s <<')' - << std::endl - << std::endl; + << std::endl + << std::endl; const unsigned int n_couplings = dof_handler.max_couplings_between_dofs(); @@ -805,11 +805,11 @@ void TwoPhaseFlowProblem::assemble_system () QGauss face_quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | + update_values | update_gradients | update_q_points | update_JxW_values); FEFaceValues fe_face_values (fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_q_points | update_JxW_values); + update_values | update_normal_vectors | + update_q_points | update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; @@ -842,67 +842,67 @@ void TwoPhaseFlowProblem::assemble_system () local_matrix = 0; local_rhs = 0; - // Here's the first significant - // difference: 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: + // Here's the first significant + // difference: 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: fe_values.get_function_values (old_solution, old_solution_values); - // Then we also have to get the values - // of the pressure right hand side and - // of the inverse permeability tensor - // at the quadrature points: + // Then we also have to get the values + // of the pressure right hand side and + // of the inverse permeability tensor + // at the quadrature points: pressure_right_hand_side.value_list (fe_values.get_quadrature_points(), - pressure_rhs_values); + pressure_rhs_values); k_inverse.value_list (fe_values.get_quadrature_points(), k_inverse_values); - // With all this, we can now loop over - // all the quadrature points and shape - // functions on this cell and assemble - // those parts of the matrix and right - // hand side that we deal with in this - // function. The individual terms in - // the contributions should be - // self-explanatory given the explicit - // form of the bilinear form stated in - // the introduction: + // With all this, we can now loop over + // all the quadrature points and shape + // functions on this cell and assemble + // those parts of the matrix and right + // hand side that we deal with in this + // function. The individual terms in + // the contributions should be + // self-explanatory given the explicit + // form of the bilinear form stated in + // the introduction: for (unsigned int q=0; q phi_i_u = extract_u (fe_values, i, q); - const double div_phi_i_u = extract_div_u (fe_values, i, q); + const double div_phi_i_u = extract_div_u (fe_values, i, q); const double phi_i_p = extract_p (fe_values, i, q); - const double phi_i_s = extract_s (fe_values, i, q); - const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q); + const double phi_i_s = extract_s (fe_values, i, q); + const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q); for (unsigned int j=0; j phi_j_u = extract_u (fe_values, j, q); - const double div_phi_j_u = extract_div_u (fe_values, j, q); + const double div_phi_j_u = extract_div_u (fe_values, j, q); const double phi_j_p = extract_p (fe_values, j, q); const double phi_j_s = extract_s (fe_values, j, q); - + local_matrix(i,j) += (phi_i_u * k_inverse_values[q] * - mobility_inverse(old_s,viscosity) * phi_j_u + mobility_inverse(old_s,viscosity) * phi_j_u - div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u - + phi_i_s * phi_j_s) + + phi_i_s * phi_j_s) * fe_values.JxW(q); } @@ -911,32 +911,32 @@ void TwoPhaseFlowProblem::assemble_system () } - // Next, we also have to deal with the - // pressure boundary values. This, - // again is as in step-20: + // Next, we also have to deal with the + // pressure boundary values. This, + // again is as in step-20: for (unsigned int face_no=0; - face_no::faces_per_cell; - ++face_no) - if (cell->at_boundary(face_no)) - { - fe_face_values.reinit (cell, face_no); - - pressure_boundary_values - .value_list (fe_face_values.get_quadrature_points(), - boundary_values); - - for (unsigned int q=0; q - phi_i_u = extract_u (fe_face_values, i, q); - - local_rhs(i) += -(phi_i_u * - fe_face_values.normal_vector(q) * - boundary_values[q] * - fe_face_values.JxW(q)); - } - } + face_no::faces_per_cell; + ++face_no) + if (cell->at_boundary(face_no)) + { + fe_face_values.reinit (cell, face_no); + + pressure_boundary_values + .value_list (fe_face_values.get_quadrature_points(), + boundary_values); + + for (unsigned int q=0; q + phi_i_u = extract_u (fe_face_values, i, q); + + local_rhs(i) += -(phi_i_u * + fe_face_values.normal_vector(q) * + boundary_values[q] * + fe_face_values.JxW(q)); + } + } // The final step in the loop // over all cells is to @@ -947,10 +947,10 @@ void TwoPhaseFlowProblem::assemble_system () for (unsigned int i=0; i::assemble_system () } - // So much for assembly of matrix and right - // hand side. Note that we do not have to - // interpolate and apply boundary values - // since they have all been taken care of in - // the weak form already. + // So much for assembly of matrix and right + // hand side. Note that we do not have to + // interpolate and apply boundary values + // since they have all been taken care of in + // the weak form already. - // @sect4{TwoPhaseFlowProblem::assemble_rhs_S} + // @sect4{TwoPhaseFlowProblem::assemble_rhs_S} - // As explained in the introduction, we can - // only evaluate the right hand side of the - // saturation equation once the velocity has - // been computed. We therefore have this - // separate function to this end. + // As explained in the introduction, we can + // only evaluate the right hand side of the + // saturation equation once the velocity has + // been computed. We therefore have this + // separate function to this end. template void TwoPhaseFlowProblem::assemble_rhs_S () { QGauss quadrature_formula(degree+2); QGauss face_quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_q_points | update_JxW_values); + update_values | update_gradients | + update_q_points | update_JxW_values); FEFaceValues fe_face_values (fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_q_points | update_JxW_values); + update_values | update_normal_vectors | + update_q_points | update_JxW_values); FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, - update_values); + update_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.n_quadrature_points; @@ -1014,107 +1014,107 @@ void TwoPhaseFlowProblem::assemble_rhs_S () fe_values.get_function_values (old_solution, old_solution_values); fe_values.get_function_values (solution, present_solution_values); - // First for the cell terms. These are, - // following the formulas in the - // introduction, $(S^n,\sigma)-(F(S^n) - // \mathbf{v}^{n+1},\nabla sigma)$, - // where $\sigma$ is the saturation - // component of the test function: + // First for the cell terms. These are, + // following the formulas in the + // introduction, $(S^n,\sigma)-(F(S^n) + // \mathbf{v}^{n+1},\nabla sigma)$, + // where $\sigma$ is the saturation + // component of the test function: for (unsigned int q=0; q present_u; - for (unsigned int d=0; d grad_phi_i_s = extract_grad_s(fe_values, i, q); - - local_rhs(i) += (time_step * - f_saturation(old_s,viscosity) * - present_u * - grad_phi_i_s - + - old_s * phi_i_s) - * - fe_values.JxW(q); - } - - // Secondly, we have to deal with the - // flux parts on the face - // boundaries. This was a bit more - // involved because we first have to - // determine which are the influx and - // outflux parts of the cell - // boundary. If we have an influx - // boundary, we need to evaluate the - // saturation on the other side of the - // face (or the boundary values, if we - // are at the boundary of the domain). - // - // All this is a bit tricky, but has - // been explained in some detail - // already in step-9. Take a look there - // how this is supposed to work! + for (unsigned int i=0; i present_u; + for (unsigned int d=0; d grad_phi_i_s = extract_grad_s(fe_values, i, q); + + local_rhs(i) += (time_step * + f_saturation(old_s,viscosity) * + present_u * + grad_phi_i_s + + + old_s * phi_i_s) + * + fe_values.JxW(q); + } + + // Secondly, we have to deal with the + // flux parts on the face + // boundaries. This was a bit more + // involved because we first have to + // determine which are the influx and + // outflux parts of the cell + // boundary. If we have an influx + // boundary, we need to evaluate the + // saturation on the other side of the + // face (or the boundary values, if we + // are at the boundary of the domain). + // + // All this is a bit tricky, but has + // been explained in some detail + // already in step-9. Take a look there + // how this is supposed to work! for (unsigned int face_no=0; face_no::faces_per_cell; - ++face_no) - { - fe_face_values.reinit (cell, face_no); - - fe_face_values.get_function_values (old_solution, old_solution_values_face); - fe_face_values.get_function_values (solution, present_solution_values_face); - - if (cell->at_boundary(face_no)) - saturation_boundary_values - .value_list (fe_face_values.get_quadrature_points(), - neighbor_saturation); - else - { - const typename DoFHandler::active_cell_iterator - neighbor = cell->neighbor(face_no); - const unsigned int - neighbor_face = cell->neighbor_of_neighbor(face_no); - - fe_face_values_neighbor.reinit (neighbor, neighbor_face); - - fe_face_values_neighbor - .get_function_values (old_solution, - old_solution_values_face_neighbor); - - for (unsigned int q=0; qat_boundary(face_no)) + saturation_boundary_values + .value_list (fe_face_values.get_quadrature_points(), + neighbor_saturation); + else + { + const typename DoFHandler::active_cell_iterator + neighbor = cell->neighbor(face_no); + const unsigned int + neighbor_face = cell->neighbor_of_neighbor(face_no); + + fe_face_values_neighbor.reinit (neighbor, neighbor_face); + + fe_face_values_neighbor + .get_function_values (old_solution, + old_solution_values_face_neighbor); + + for (unsigned int q=0; q present_u_face; - for (unsigned int d=0; d= 0); - - for (unsigned int i=0; i present_u_face; + for (unsigned int d=0; d= 0); + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); for (unsigned int i=0; i::solve () Vector tmp2 (solution.block(2).size()); - // First the pressure, using the pressure - // Schur complement of the first two - // equations: + // First the pressure, using the pressure + // Schur complement of the first two + // equations: { m_inverse.vmult (tmp, system_rhs.block(0)); system_matrix.block(1,0).vmult (schur_rhs, tmp); @@ -1157,14 +1157,14 @@ void TwoPhaseFlowProblem::solve () SolverControl solver_control (system_matrix.block(0,0).m(), - 1e-12*schur_rhs.l2_norm()); + 1e-12*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); cg.solve (schur_complement, solution.block(1), schur_rhs, preconditioner); std::cout << " " - << solver_control.last_step() + << solver_control.last_step() << " CG Schur complement iterations for pressure." << std::endl; } @@ -1178,43 +1178,43 @@ void TwoPhaseFlowProblem::solve () m_inverse.vmult (solution.block(0), tmp); } - // Finally, we have to take care of the - // saturation equation. The first business - // we have here is to determine the time - // step using the formula in the - // introduction. Knowing the shape of our - // domain and that we created the mesh by - // regular subdivision of cells, we can - // compute the diameter of each of our - // cells quite easily (in fact we use the - // linear extensions in coordinate - // directions of the cells, not the - // diameter). The maximal velocity we - // compute using a helper function defined - // below: + // Finally, we have to take care of the + // saturation equation. The first business + // we have here is to determine the time + // step using the formula in the + // introduction. Knowing the shape of our + // domain and that we created the mesh by + // regular subdivision of cells, we can + // compute the diameter of each of our + // cells quite easily (in fact we use the + // linear extensions in coordinate + // directions of the cells, not the + // diameter). The maximal velocity we + // compute using a helper function defined + // below: time_step = std::pow(0.5, double(n_refinement_steps)) / - get_maximal_velocity(); + get_maximal_velocity(); - // The next step is to assemble the right - // hand side, and then to pass everything - // on for solution. At the end, we project - // back saturations onto the physically - // reasonable range: + // The next step is to assemble the right + // hand side, and then to pass everything + // on for solution. At the end, we project + // back saturations onto the physically + // reasonable range: assemble_rhs_S (); { SolverControl solver_control (system_matrix.block(2,2).m(), - 1e-8*system_rhs.block(2).l2_norm()); + 1e-8*system_rhs.block(2).l2_norm()); SolverCG<> cg (solver_control); cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), - PreconditionIdentity()); - + PreconditionIdentity()); + project_back_saturation (); - + std::cout << " " - << solver_control.last_step() + << solver_control.last_step() << " CG iterations for saturation." - << std::endl; + << std::endl; } @@ -1224,10 +1224,16 @@ void TwoPhaseFlowProblem::solve () // @sect4{TwoPhaseFlowProblem::output_results} - // There is nothing surprising here: + // There is nothing surprising here. Since + // the program will do a lot of time steps, + // we create an output file only every fifth + // time step. template void TwoPhaseFlowProblem::output_results () const -{ +{ + if (timestep_number % 5 != 0) + return; + std::vector solution_names; switch (dim) { @@ -1235,7 +1241,7 @@ void TwoPhaseFlowProblem::output_results () const solution_names.push_back ("u"); solution_names.push_back ("v"); solution_names.push_back ("p"); - solution_names.push_back ("S"); + solution_names.push_back ("S"); break; case 3: @@ -1243,7 +1249,7 @@ void TwoPhaseFlowProblem::output_results () const solution_names.push_back ("v"); solution_names.push_back ("w"); solution_names.push_back ("p"); - solution_names.push_back ("S"); + solution_names.push_back ("S"); break; default: @@ -1266,27 +1272,53 @@ void TwoPhaseFlowProblem::output_results () const + // @sect4{TwoPhaseFlowProblem::project_back_saturation} + + // In this function, we simply run over all + // saturation degrees of freedom and make + // sure that if they should have left the + // physically reasonable range, that they be + // reset to the interval $[0,1]$. To do this, + // we only have to loop over all saturation + // components of the solution vector; these + // are stored in the block 2 (block 0 are the + // velocities, block 1 are the pressures). + // + // It may be instructive to note that this + // function almost never triggers when the + // time step is chosen as mentioned in the + // introduction. However, if we choose the + // timestep only slightly larger, we get + // plenty of values outside the proper + // range. Strictly speaking, the function is + // therefore unnecessary if we choose the + // time step small enough. In a sense, the + // function is therefore only a safety device + // to avoid situations where our entire + // solution becomes unphysical because + // individual degrees of freedom have become + // unphysical a few time steps earlier. template void TwoPhaseFlowProblem::project_back_saturation () { for (unsigned int i=0; i 1) - { - std::cout << "xxx " << solution.block(2)(i) << std::endl; - solution.block(2)(i) = 1; - } + solution.block(2)(i) = 1; } + // @sect4{TwoPhaseFlowProblem::get_maximal_velocity} + // The following function is used in + // determining the maximal allowable time + // step. What it does is to loop over all + // quadrature points in the domain and find + // what the maximal magnitude of the velocity + // is. template double TwoPhaseFlowProblem::get_maximal_velocity () const @@ -1296,9 +1328,9 @@ TwoPhaseFlowProblem::get_maximal_velocity () const = quadrature_formula.n_quadrature_points; FEValues fe_values (fe, quadrature_formula, - update_values); + update_values); std::vector > solution_values(n_q_points, - Vector(dim+2)); + Vector(dim+2)); double max_velocity = 0; typename DoFHandler::active_cell_iterator @@ -1310,14 +1342,14 @@ TwoPhaseFlowProblem::get_maximal_velocity () const fe_values.get_function_values (solution, solution_values); for (unsigned int q=0; q velocity; - for (unsigned int i=0; i velocity; + for (unsigned int i=0; i::get_maximal_velocity () const // @sect4{TwoPhaseFlowProblem::run} - // This is the final function of our - // main class. It's only job is to - // call the other functions in their order: + // This is the final function of our main + // class. Its brevity speaks for + // itself. There are only two points worth + // noting: First, the function projects the + // initial values onto the finite element + // space at the beginning; the + // VectorTools::project function doing this + // requires an argument indicating the + // hanging node constraints. We have none in + // this program (we compute on a uniformly + // refined mesh), but the function requires + // the argument anyway, of course. So we have + // to create a constraint object. In its + // original state, constraint objects are + // unsorted, and have to be sorted (using the + // ConstraintMatrix::close function) before + // they can be used. This is what we do here, + // and which is why we can't simply call the + // VectorTools::project function with an + // anonymous temporary object + // ConstraintMatrix() as the + // second argument. + // + // The second point worth mentioning is that + // we only compute the length of the present + // time step in the middle of solving the + // linear system corresponding to each time + // step. We can therefore output the present + // end time of a time step only at the end of + // the time step. template void TwoPhaseFlowProblem::run () { make_grid_and_dofs(); - ConstraintMatrix constraints; - constraints.close(); - - std::list production_rate; - std::list production_time; - - Vector tmp (old_solution.size()); - VectorTools::project (dof_handler, constraints, QGauss(degree+2),InitialValues(),tmp); - std::copy (tmp.begin(), tmp.end(), old_solution.begin()); + { + ConstraintMatrix constraints; + constraints.close(); + + VectorTools::project (dof_handler, + constraints, + QGauss(degree+2), + InitialValues(), + old_solution); + } timestep_number = 1; double time = 0; @@ -1350,7 +1410,7 @@ void TwoPhaseFlowProblem::run () do { std::cout << "Timestep " << timestep_number - << std::endl; + << std::endl; assemble_system (); @@ -1361,9 +1421,9 @@ void TwoPhaseFlowProblem::run () time += time_step; ++timestep_number; std::cout << " Now at t=" << time - << ", dt=" << time_step << '.' - << std::endl - << std::endl; + << ", dt=" << time_step << '.' + << std::endl + << std::endl; } while (time <= 250); } @@ -1371,12 +1431,13 @@ void TwoPhaseFlowProblem::run () // @sect3{The main function} - // In the main function, we pass the degree - // of the finite element space to the - // constructor of the TwoPhaseFlowProblem. - // Here, we use zero-th degree elements, - // i.e. $RT_0\times DQ_0 \times DQ_0$. The - // rest is as in all the other programs. + // That's it. In the main function, we pass + // the degree of the finite element space to + // the constructor of the TwoPhaseFlowProblem + // object. Here, we use zero-th degree + // elements, i.e. $RT_0\times DQ_0 \times + // DQ_0$. The rest is as in all the other + // programs. int main () { try @@ -1389,25 +1450,25 @@ int main () catch (std::exception &exc) { std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; + << "----------------------------------------------------" + << std::endl; std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; + << "----------------------------------------------------" + << std::endl; std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; return 1; } -- 2.39.5