From: bangerth Date: Tue, 7 Feb 2012 15:32:48 +0000 (+0000) Subject: More reshuffling and documenting. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=28d95ad984acae1dee87b76270d782ae141e9acd;p=dealii-svn.git More reshuffling and documenting. git-svn-id: https://svn.dealii.org/trunk@25006 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-43/step-43.cc b/deal.II/examples/step-43/step-43.cc index 13ffd2320a..d7d897d87f 100644 --- a/deal.II/examples/step-43/step-43.cc +++ b/deal.II/examples/step-43/step-43.cc @@ -77,6 +77,295 @@ namespace Step43 using namespace dealii; + // @sect3{Pressure right hand side, pressure boundary values and saturation initial value classes} + + // The following part is taken + // directly from step-21 so there is + // no need to repeat the + // descriptions found there. + template + class PressureRightHandSide : public Function + { + public: + PressureRightHandSide () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + + + template + double + PressureRightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } + + + template + class PressureBoundaryValues : public Function + { + public: + PressureBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + + template + double + PressureBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const + { + return 1-p[0]; + } + + + template + class SaturationBoundaryValues : public Function + { + public: + SaturationBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + }; + + + + template + double + SaturationBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const + { + if (p[0] == 0) + return 1; + else + return 0; + } + + + template + class SaturationInitialValues : public Function + { + public: + SaturationInitialValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + }; + + + template + double + SaturationInitialValues::value (const Point &/*p*/, + const unsigned int /*component*/) const + { + return 0; + } + + + template + void + SaturationInitialValues::vector_value (const Point &p, + Vector &values) const + { + for (unsigned int c=0; cn_components; ++c) + values(c) = SaturationInitialValues::value (p,c); + } + + + // @sect3{Permeability models} + + // In this tutorial, we still use + // the two permeability models + // previously used in step-21 so we + // again refrain from commenting in + // detail about them. + namespace SingleCurvingCrack + { + template + class KInverse : public TensorFunction<2,dim> + { + public: + KInverse () + : + TensorFunction<2,dim> () + {} + + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + }; + + + template + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p + class KInverse : public TensorFunction<2,dim> + { + public: + KInverse () + : + TensorFunction<2,dim> () + {} + + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + + private: + static std::vector > centers; + + static std::vector > get_centers (); + }; + + + + template + std::vector > + KInverse::centers = KInverse::get_centers(); + + + template + std::vector > + KInverse::get_centers () + { + const unsigned int N = (dim == 2 ? + 40 : + (dim == 3 ? + 100 : + throw ExcNotImplemented())); + + std::vector > centers_list (N); + for (unsigned int i=0; i(rand())/RAND_MAX; + + return centers_list; + } + + + + template + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p= 0) && (S<=1), + ExcMessage ("Saturation is outside its physically valid range.")); + + return S*S / ( S * S + viscosity * (1-S) * (1-S)); + } + + + double fractional_flow_derivative (const double S, + const double viscosity) + { + Assert ((S >= 0) && (S<=1), + ExcMessage ("Saturation is outside its physically valid range.")); + + const double temp = ( S * S + viscosity * (1-S) * (1-S) ); + + const double numerator = 2.0 * S * temp + - + S * S * + ( 2.0 * S - 2.0 * viscosity * (1-S) ); + const double denominator = std::pow(temp, 2.0); + + const double F_prime = numerator / denominator; + + Assert (F_prime >= 0, ExcInternalError()); + + return F_prime; + } + + // @sect3{Helper classes for solvers and preconditioners} // In this first part we define a @@ -245,10 +534,10 @@ namespace Step43 // pressure variable. We need this // because we are building a // Laplace matrix for the pressure - // *as an approximation of the - // Schur complement) which is only - // definite if boundary conditions - // are applied. + // as an approximation of the Schur + // complement) which is only + // positive definite if boundary + // conditions are applied. // // The collection of member // functions and variables thus @@ -281,7 +570,8 @@ namespace Step43 void solve (); void compute_refinement_indicators (const TrilinosWrappers::Vector &predicted_saturation_solution, Vector &refinement_indicators) const; - void refine_mesh (const unsigned int max_grid_level, + void refine_mesh (const unsigned int min_grid_level, + const unsigned int max_grid_level, const Vector &indicator); void output_results () const; @@ -364,300 +654,49 @@ namespace Step43 double AOS_threshold; std_cxx1x::shared_ptr Amg_preconditioner; - std_cxx1x::shared_ptr Mp_preconditioner; + std_cxx1x::shared_ptr Mp_preconditioner; bool rebuild_saturation_matrix; - }; - - - // @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 - { - public: - PressureRightHandSide () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - }; - - - - template - double - PressureRightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const - { - return 0; - } - - - template - class PressureBoundaryValues : public Function - { - public: - PressureBoundaryValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - }; - - - template - double - PressureBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const - { - return 1-p[0]; - } - - - template - class SaturationBoundaryValues : public Function - { - public: - SaturationBoundaryValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - }; - - - - template - double - SaturationBoundaryValues::value (const Point &p, - const unsigned int /*component*/) const - { - if (p[0] == 0) - return 1; - else - return 0; - } - - - template - class SaturationInitialValues : public Function - { - public: - SaturationInitialValues () : Function(1) {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void vector_value (const Point &p, - Vector &value) const; + // At the very end we declare a + // variable that denotes the + // material model. Compared to + // step-21, we do this here as + // a member variable since we + // will want to use it in a + // variety of places and so + // having a central place where + // such a variable is declared + // will make it simpler to + // replace one class by another + // (e.g. replace + // RandomMedium::KInverse by + // SingleCurvingCrack::KInverse). + const RandomMedium::KInverse k_inverse; }; - template - double - SaturationInitialValues::value (const Point &/*p*/, - const unsigned int /*component*/) const - { - return 0; - } - - - template - void - SaturationInitialValues::vector_value (const Point &p, - Vector &values) const - { - for (unsigned int c=0; cn_components; ++c) - values(c) = SaturationInitialValues::value (p,c); - } - - - // @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 - class KInverse : public TensorFunction<2,dim> - { - public: - KInverse () - : - TensorFunction<2,dim> () - {} - - virtual void value_list (const std::vector > &points, - std::vector > &values) const; - }; - - - template - void - KInverse::value_list (const std::vector > &points, - std::vector > &values) const - { - Assert (points.size() == values.size(), - ExcDimensionMismatch (points.size(), values.size())); - - for (unsigned int p=0; p - class KInverse : public TensorFunction<2,dim> - { - public: - KInverse () - : - TensorFunction<2,dim> () - {} - - virtual void value_list (const std::vector > &points, - std::vector > &values) const; - - private: - static std::vector > centers; - - static std::vector > get_centers (); - }; - - - - template - std::vector > - KInverse::centers = KInverse::get_centers(); - - - template - std::vector > - KInverse::get_centers () - { - const unsigned int N = (dim == 2 ? - 40 : - (dim == 3 ? - 100 : - throw ExcNotImplemented())); - - std::vector > centers_list (N); - for (unsigned int i=0; i(rand())/RAND_MAX; - - return centers_list; - } - - - - template - void - KInverse::value_list (const std::vector > &points, - std::vector > &values) const - { - Assert (points.size() == values.size(), - ExcDimensionMismatch (points.size(), values.size())); - - for (unsigned int p=0; p::TwoPhaseFlowProblem} // The constructor of this class is an - // extension of the constructor in step-21 + // extension of the constructors 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 + // (Taylor-Hood) elements again for the Darcy + // system, an element combination that 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 + // 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 + // time stepping variables related to // operator splitting as well as the option // for matrix assembly and preconditioning: template @@ -692,58 +731,67 @@ namespace Step43 // This is the function that sets up the // DoFHandler objects we have here (one for - // the darcy part and one for the saturation + // 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 + // operations are similar to what // step-31 did. // // The body of the function first enumerates - // all degrees of freedom for the darcy and - // saturation systems. For the darcy part, + // 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. + // DoFs so that we can partition the Darcy + // matrix into a $2 \times 2$ matrix. // - // Then, we need to incorporate hanging node - // constraints and Dirichlet boundary value + // 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: + // darcy_preconditioner_constraints. + // The boundary condition + // constraints are only set on 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$, acts only + // on the pressure + // variable. Therefore, we use a + // component_mask that filters out + // the velocity component, so that + // the condensation is performed on + // pressure degrees of freedom + // only. // - // 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 + // After having done so, we count + // the number of degrees of freedom + // in the various blocks. This + // information is then used 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 + // 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 () @@ -897,11 +945,20 @@ namespace Step43 } - // @sect3{TwoPhaseFlowProblem::assemble_darcy_preconditioner} + // @sect3{Assembling matrices and preconditioners} + + // The next few functions are + // devoted to setting up the + // various system and + // preconditioner matrices and + // right hand sides that we have to + // deal with in this program. + + // @sect4{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 + // for preconditioning the Darcy system. What + // we need are a vector mass matrix weighted by // $\left(\mathbf{K} \lambda_t\right)^{-1}$ // on the velocity components and a mass // matrix weighted by $\left(\mathbf{K} @@ -922,27 +979,41 @@ namespace Step43 // 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 + // 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 + // \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. + // that takes care of hanging node + // and zero Dirichlet boundary + // condition constraints. By doing + // so, we don't have to do that + // afterwards, and we later don't + // have to use + // ConstraintMatrix::condense and + // MatrixTools::apply_boundary_values, + // both functions that would need + // to modify matrix and vector + // entries and so are difficult to + // write for the Trilinos classes + // where we don't immediately have + // access to individual memory + // locations. template void TwoPhaseFlowProblem::assemble_darcy_preconditioner () @@ -963,9 +1034,6 @@ namespace Step43 const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); - const RandomMedium::KInverse k_inverse; -// const SingleCurvingCrack::KInverse k_inverse; - std::vector > k_inverse_values (n_q_points); Tensor<2,dim> k_value; @@ -1033,7 +1101,7 @@ namespace Step43 } - // @sect3{TwoPhaseFlowProblem::build_darcy_preconditioner} + // @sect4{TwoPhaseFlowProblem::build_darcy_preconditioner} // This function generates the inner // preconditioners that are going to be used @@ -1049,7 +1117,7 @@ namespace Step43 // 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 + // (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 @@ -1075,10 +1143,10 @@ namespace Step43 } - // @sect3{TwoPhaseFlowProblem::assemble_darcy_system} + // @sect4{TwoPhaseFlowProblem::assemble_darcy_system} // This is the function that assembles the - // linear system for the darcy system. + // linear system for the Darcy system. // // Regarding the technical details of // implementation, the procedures are similar @@ -1169,7 +1237,7 @@ namespace Step43 // 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 + // 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 @@ -1228,8 +1296,6 @@ namespace Step43 const PressureRightHandSide pressure_right_hand_side; const PressureBoundaryValues pressure_boundary_values; - const RandomMedium::KInverse k_inverse; -// const SingleCurvingCrack::KInverse k_inverse; std::vector pressure_rhs_values (n_q_points); std::vector boundary_values (n_face_q_points); @@ -1330,7 +1396,7 @@ namespace Step43 } - // @sect3{TwoPhaseFlowProblem::assemble_saturation_system} + // @sect4{TwoPhaseFlowProblem::assemble_saturation_system} // This function is to assemble the linear // system for the saturation transport @@ -1361,7 +1427,7 @@ namespace Step43 - // @sect3{TwoPhaseFlowProblem::assemble_saturation_matrix} + // @sect4{TwoPhaseFlowProblem::assemble_saturation_matrix} // This function is easily understood since // it only forms a simple mass matrix for the @@ -1422,12 +1488,12 @@ namespace Step43 - // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs} + // @sect4{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 + // 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 @@ -1446,7 +1512,7 @@ namespace Step43 // step-31. // // Next, we start to loop over all the - // saturation and darcy cells to put 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 @@ -1535,7 +1601,7 @@ namespace Step43 - // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs_cell_term} + // @sect4{TwoPhaseFlowProblem::assemble_saturation_rhs_cell_term} // In this function, we actually compute // every artificial viscosity for every @@ -1598,7 +1664,7 @@ namespace Step43 const Tensor<1,dim> grad_phi_i_s = saturation_fe_values.shape_grad (i, q); local_rhs(i) += (time_step * - f_saturation(old_s,viscosity) * + fractional_flow(old_s,viscosity) * present_u * grad_phi_i_s - @@ -1617,7 +1683,7 @@ namespace Step43 } - // @sect3{TwoPhaseFlowProblem::assemble_saturation_rhs_boundary_term} + // @sect4{TwoPhaseFlowProblem::assemble_saturation_rhs_boundary_term} // In this function, we have to give // upwinding in the global boundary faces, @@ -1664,12 +1730,12 @@ namespace Step43 for (unsigned int i=0; i old_saturation_after_solving_pressure (n_q_points); std::vector present_saturation (n_q_points); - const RandomMedium::KInverse k_inverse; -// const SingleCurvingCrack::KInverse k_inverse; - std::vector > k_inverse_values (n_q_points); double max_global_aop_indicator = 0.0; @@ -2017,9 +2080,11 @@ namespace Step43 template void TwoPhaseFlowProblem:: - refine_mesh (const unsigned int max_grid_level, + refine_mesh (const unsigned int min_grid_level, + const unsigned int max_grid_level, const Vector &refinement_indicators) { + //TODO: use a useful refinement criterion, in much the same way as we do in step-31 { typename DoFHandler::active_cell_iterator cell = saturation_dof_handler.begin_active(), @@ -2034,7 +2099,8 @@ namespace Step43 (std::fabs(refinement_indicators(cell_no)) > saturation_refinement_threshold)) cell->set_refine_flag(); else - if (std::fabs(refinement_indicators(cell_no)) < 0.75 * saturation_refinement_threshold) + if ((static_cast(cell->level()) > min_grid_level) && + (std::fabs(refinement_indicators(cell_no)) < 0.75 * saturation_refinement_threshold)) cell->set_coarsen_flag(); } } @@ -2307,10 +2373,10 @@ namespace Step43 for (unsigned int i=0; i refinement_indicators (triangulation.n_active_cells()); compute_refinement_indicators(predicted_saturation_solution, refinement_indicators); - refine_mesh (initial_refinement + n_pre_refinement_steps, + refine_mesh (initial_refinement, + initial_refinement + n_pre_refinement_steps, refinement_indicators); } @@ -2564,7 +2632,15 @@ namespace Step43 -int main () + // @sect3{The main function} + // + // The main function looks almost the + // same as in all other programs. In + // particular, it is essentially the + // same as in step-31 where we also + // explain the need to initialize the + // MPI subsystem. +int main (int argc, char *argv[]) { try { @@ -2573,6 +2649,8 @@ int main () deallog.depth_console (0); + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv); + TwoPhaseFlowProblem<3> two_phase_flow_problem(1); two_phase_flow_problem.run (); }