From 8b31f17de93afea56d77cd78cef0de71fd46d3cc Mon Sep 17 00:00:00 2001 From: David Wells Date: Fri, 10 May 2019 23:17:28 -0400 Subject: [PATCH] Modernize step-21. 1. Ranged for-loops 2. Improve discussion relative to step-20 3. Use DGQ, not DQ, for standard discontinuous elements 4. inline function definitions --- examples/step-21/doc/intro.dox | 8 +- examples/step-21/doc/results.dox | 2 +- examples/step-21/step-21.cc | 327 ++++++++++++------------------- 3 files changed, 125 insertions(+), 212 deletions(-) diff --git a/examples/step-21/doc/intro.dox b/examples/step-21/doc/intro.dox index 90b63a47ca..4c1f2c8d96 100644 --- a/examples/step-21/doc/intro.dox +++ b/examples/step-21/doc/intro.dox @@ -238,7 +238,7 @@ cell term to get an equation as follows: In each time step, we then apply the mixed finite method of @ref step_20 "step-20" to the velocity and pressure. To be well-posed, we choose Raviart-Thomas spaces $RT_{k}$ for $\mathbf{u}$ and discontinuous elements of -class $DQ_{k}$ for $p$. For the saturation, we will also choose $DQ_{k}$ +class $DGQ_{k}$ for $p$. For the saturation, we will also choose $DGQ_{k}$ spaces. Since we have discontinuous spaces, we have to think about how to evaluate @@ -267,8 +267,8 @@ fluxes can also be found in step-12.

Linear solvers

The linear solvers used in this program are a straightforward extension of the -ones used in step-20. Essentially, we simply have to extend -everything from +ones used in step-20 (but without LinearOperator). Essentially, we simply have +to extend everything from two to three solution components. If we use the discrete spaces mentioned above and put shape functions into the bilinear forms, we arrive at the following linear system to be solved for time step $n+1$: @@ -294,7 +294,7 @@ B & 0 & 0\\ @f] where the individual matrices and vectors are defined as follows using shape functions $\mathbf v_i$ (of type Raviart Thomas $RT_k$) for -velocities and $\phi_i$ (of type $DG_k$) for both pressures and saturations: +velocities and $\phi_i$ (of type $DGQ_k$) for both pressures and saturations: @f{eqnarray*} M^u(S^n)_{ij} &=& \left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{v}_i,\mathbf diff --git a/examples/step-21/doc/results.dox b/examples/step-21/doc/results.dox index cf2fb3232d..66345340db 100644 --- a/examples/step-21/doc/results.dox +++ b/examples/step-21/doc/results.dox @@ -42,7 +42,7 @@ the pressure once is significantly reduced. The final observation concerns the number of iterations needed to solve for the saturation, i.e. one. This shouldn't surprise us too much: the matrix we have to solve with is the mass matrix. However, this is the mass matrix for -the $DQ_0$ element of piecewise constants where no element couples with the +the $DGQ_0$ element of piecewise constants where no element couples with the degrees of freedom on neighboring cells. The matrix is therefore a diagonal one, and it is clear that we should be able to invert this matrix in a single CG iteration. diff --git a/examples/step-21/step-21.cc b/examples/step-21/step-21.cc index 4b97da08e5..897de3cadf 100644 --- a/examples/step-21/step-21.cc +++ b/examples/step-21/step-21.cc @@ -144,21 +144,15 @@ namespace Step21 : Function(1) {} - virtual double value(const Point & p, - const unsigned int component = 0) const override; + virtual double value(const Point & /*p*/, + const unsigned int /*component*/ = 0) const override + { + return 0; + } }; - template - double - PressureRightHandSide::value(const Point & /*p*/, - const unsigned int /*component*/) const - { - return 0; - } - - // @sect4{Pressure boundary values} // The next are pressure boundary values. As mentioned in the introduction, @@ -171,19 +165,14 @@ namespace Step21 : Function(1) {} - virtual double value(const Point & p, - const unsigned int component = 0) const override; + virtual double value(const Point &p, + const unsigned int /*component*/ = 0) const override + { + return 1 - p[0]; + } }; - template - double - PressureBoundaryValues::value(const Point &p, - const unsigned int /*component*/) const - { - return 1 - p[0]; - } - // @sect4{Saturation boundary values} @@ -200,25 +189,18 @@ namespace Step21 : Function(1) {} - virtual double value(const Point & p, - const unsigned int component = 0) const override; + virtual double value(const Point &p, + const unsigned int /*component*/ = 0) const override + { + if (p[0] == 0) + return 1; + else + return 0; + } }; - template - double - SaturationBoundaryValues::value(const Point &p, - const unsigned int /*component*/) const - { - if (p[0] == 0) - return 1; - else - return 0; - } - - - // @sect4{Initial data} // Finally, we need initial data. In reality, we only need initial data for @@ -241,29 +223,19 @@ namespace Step21 {} virtual double value(const Point & p, - const unsigned int component = 0) const override; + const unsigned int component = 0) const override + { + return Functions::ZeroFunction(dim + 2).value(p, component); + } virtual void vector_value(const Point &p, - Vector & value) const override; + Vector & values) const override + { + Functions::ZeroFunction(dim + 2).vector_value(p, values); + } }; - template - double InitialValues::value(const Point & p, - const unsigned int component) const - { - return Functions::ZeroFunction(dim + 2).value(p, component); - } - - - template - void InitialValues::vector_value(const Point &p, - Vector & values) const - { - Functions::ZeroFunction(dim + 2).vector_value(p, values); - } - - // @sect3{The inverse permeability tensor} @@ -289,34 +261,30 @@ namespace Step21 : 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())); + virtual void + value_list(const std::vector> &points, + std::vector> & values) const override + { + Assert(points.size() == values.size(), + ExcDimensionMismatch(points.size(), values.size())); - for (unsigned int p = 0; p < points.size(); ++p) - { - values[p].clear(); + for (unsigned int p = 0; p < points.size(); ++p) + { + values[p].clear(); - const double distance_to_flowline = - std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0])); + const double distance_to_flowline = + std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0])); - const double permeability = - std::max(std::exp(-(distance_to_flowline * distance_to_flowline) / - (0.1 * 0.1)), - 0.01); + const double permeability = + std::max(std::exp(-(distance_to_flowline * distance_to_flowline) / + (0.1 * 0.1)), + 0.01); - for (unsigned int d = 0; d < dim; ++d) - values[p][d][d] = 1. / permeability; - } - } + for (unsigned int d = 0; d < dim; ++d) + values[p][d][d] = 1. / permeability; + } + } + }; } // namespace SingleCurvingCrack @@ -364,60 +332,50 @@ namespace Step21 virtual void value_list(const std::vector> &points, - std::vector> & values) const override; - - private: - static std::vector> centers; + std::vector> & values) const override + { + Assert(points.size() == values.size(), + ExcDimensionMismatch(points.size(), values.size())); - static std::vector> get_centers(); - }; + for (unsigned int p = 0; p < points.size(); ++p) + { + values[p].clear(); + double permeability = 0; + for (unsigned int i = 0; i < centers.size(); ++i) + permeability += std::exp(-(points[p] - centers[i]).norm_square() / + (0.05 * 0.05)); + const double normalized_permeability = + std::min(std::max(permeability, 0.01), 4.); - template - std::vector> - KInverse::centers = KInverse::get_centers(); + for (unsigned int d = 0; d < dim; ++d) + values[p][d][d] = 1. / normalized_permeability; + } + } + private: + static std::vector> centers; - template - std::vector> KInverse::get_centers() - { - const unsigned int N = - (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented())); + static std::vector> 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 < N; ++i) - for (unsigned int d = 0; d < dim; ++d) - centers_list[i][d] = static_cast(rand()) / RAND_MAX; + std::vector> centers_list(N); + for (unsigned int i = 0; i < N; ++i) + for (unsigned int d = 0; d < dim; ++d) + centers_list[i][d] = static_cast(rand()) / RAND_MAX; - return centers_list; - } + 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 < points.size(); ++p) - { - values[p].clear(); - - double permeability = 0; - for (unsigned int i = 0; i < centers.size(); ++i) - permeability += - std::exp(-(points[p] - centers[i]).norm_square() / (0.05 * 0.05)); - - const double normalized_permeability = - std::min(std::max(permeability, 0.01), 4.); - - for (unsigned int d = 0; d < dim; ++d) - values[p][d][d] = 1. / normalized_permeability; - } - } + std::vector> + KInverse::centers = KInverse::get_centers(); } // namespace RandomMedium @@ -458,44 +416,44 @@ namespace Step21 class InverseMatrix : public Subscriptor { public: - InverseMatrix(const MatrixType &m); + InverseMatrix(const MatrixType &m) + : matrix(&m) + {} - void vmult(Vector &dst, const Vector &src) const; + void vmult(Vector &dst, const Vector &src) const + { + SolverControl solver_control(std::max(src.size(), 200), + 1e-8 * src.l2_norm()); + SolverCG<> cg(solver_control); + + dst = 0; + + cg.solve(*matrix, dst, src, PreconditionIdentity()); + } private: const SmartPointer matrix; }; - template - InverseMatrix::InverseMatrix(const MatrixType &m) - : matrix(&m) - {} - - - - template - void InverseMatrix::vmult(Vector & dst, - const Vector &src) const - { - SolverControl solver_control(std::max(src.size(), 200), - 1e-8 * src.l2_norm()); - SolverCG<> cg(solver_control); - - dst = 0; - - cg.solve(*matrix, dst, src, PreconditionIdentity()); - } - - class SchurComplement : public Subscriptor { public: SchurComplement(const BlockSparseMatrix & A, - const InverseMatrix> &Minv); + const InverseMatrix> &Minv) + : system_matrix(&A) + , m_inverse(&Minv) + , tmp1(A.block(0, 0).m()) + , tmp2(A.block(0, 0).m()) + {} - void vmult(Vector &dst, const Vector &src) const; + void vmult(Vector &dst, const Vector &src) const + { + system_matrix->block(0, 1).vmult(tmp1, src); + m_inverse->vmult(tmp2, tmp1); + system_matrix->block(1, 0).vmult(dst, tmp2); + } private: const SmartPointer> system_matrix; @@ -506,32 +464,21 @@ namespace Step21 - SchurComplement::SchurComplement( - const BlockSparseMatrix & A, - const InverseMatrix> &Minv) - : system_matrix(&A) - , m_inverse(&Minv) - , tmp1(A.block(0, 0).m()) - , tmp2(A.block(0, 0).m()) - {} - - - void SchurComplement::vmult(Vector & dst, - const Vector &src) const - { - system_matrix->block(0, 1).vmult(tmp1, src); - m_inverse->vmult(tmp2, tmp1); - system_matrix->block(1, 0).vmult(dst, tmp2); - } - - - class ApproximateSchurComplement : public Subscriptor { public: - ApproximateSchurComplement(const BlockSparseMatrix &A); + ApproximateSchurComplement(const BlockSparseMatrix &A) + : system_matrix(&A) + , tmp1(A.block(0, 0).m()) + , tmp2(A.block(0, 0).m()) + {} - void vmult(Vector &dst, const Vector &src) const; + void vmult(Vector &dst, const Vector &src) const + { + system_matrix->block(0, 1).vmult(tmp1, src); + system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1); + system_matrix->block(1, 0).vmult(dst, tmp2); + } private: const SmartPointer> system_matrix; @@ -540,23 +487,6 @@ namespace Step21 }; - ApproximateSchurComplement::ApproximateSchurComplement( - const BlockSparseMatrix &A) - : system_matrix(&A) - , tmp1(A.block(0, 0).m()) - , tmp2(A.block(0, 0).m()) - {} - - - void ApproximateSchurComplement::vmult(Vector & dst, - const Vector &src) const - { - system_matrix->block(0, 1).vmult(tmp1, src); - system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1); - system_matrix->block(1, 0).vmult(dst, tmp2); - } - - // @sect3{TwoPhaseFlowProblem class implementation} @@ -720,10 +650,7 @@ namespace Step21 const FEValuesExtractors::Scalar pressure(dim); const FEValuesExtractors::Scalar saturation(dim + 1); - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); local_matrix = 0; @@ -878,10 +805,7 @@ namespace Step21 const FEValuesExtractors::Scalar saturation(dim + 1); - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { local_rhs = 0; fe_values.reinit(cell); @@ -936,8 +860,7 @@ namespace Step21 fe_face_values.get_quadrature_points(), neighbor_saturation); else { - const typename DoFHandler::active_cell_iterator neighbor = - cell->neighbor(face_no); + const auto neighbor = cell->neighbor(face_no); const unsigned int neighbor_face = cell->neighbor_of_neighbor(face_no); @@ -1093,18 +1016,11 @@ namespace Step21 switch (dim) { case 2: - solution_names.emplace_back("u"); - solution_names.emplace_back("v"); - solution_names.emplace_back("p"); - solution_names.emplace_back("S"); + solution_names = {"u", "v", "p", "S"}; break; case 3: - solution_names.emplace_back("u"); - solution_names.emplace_back("v"); - solution_names.emplace_back("w"); - solution_names.emplace_back("p"); - solution_names.emplace_back("S"); + solution_names = {"u", "v", "w", "p", "S"}; break; default: @@ -1169,10 +1085,7 @@ namespace Step21 Vector(dim + 2)); double max_velocity = 0; - typename DoFHandler::active_cell_iterator cell = - dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell != endc; ++cell) + for (const auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); fe_values.get_function_values(solution, solution_values); -- 2.39.5