From: Matthias Maier Date: Fri, 8 Feb 2019 19:51:36 +0000 (-0600) Subject: examples/step-20: Actually use LinearOperator in the example step X-Git-Tag: v9.1.0-rc1~347^2~8 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f21656e3fd265a36fd53616888bffca7601600e5;p=dealii.git examples/step-20: Actually use LinearOperator in the example step --- diff --git a/examples/step-20/step-20.cc b/examples/step-20/step-20.cc index 8710392aca..20ddcb621d 100644 --- a/examples/step-20/step-20.cc +++ b/examples/step-20/step-20.cc @@ -14,7 +14,8 @@ * --------------------------------------------------------------------- * - * Author: Wolfgang Bangerth, Texas A&M University, 2005, 2006 + * Authors: Wolfgang Bangerth, Texas A&M University, 2005, 2006; + * (port to LinearOperator:) Matthias Maier, 2019 */ @@ -27,12 +28,18 @@ #include #include #include + #include #include #include #include #include +// The only two new header files that deserve some attention are those for +// the LinearOperator and PackagedOperation classes: +#include +#include + #include #include #include @@ -559,206 +566,91 @@ namespace Step20 // @sect3{Linear solvers and preconditioners} - // The linear solvers and preconditioners we use in this example have been - // discussed in significant detail already in the introduction. We will - // therefore not discuss the rationale for these classes here any more, but - // rather only comment on implementational aspects. - - // @sect4{The InverseMatrix class template} - - // There are a few places in this program where we will need either the - // action of the inverse of the mass matrix or the action of the inverse of - // the approximate Schur complement. Rather than explicitly calling - // SolverCG::solve every time that we need to solve such a system, we will - // wrap the action of either inverse in a simple class. The only things we - // would like to note are that this class is derived from - // Subscriptor and, as mentioned above, it stores a pointer to - // the underlying matrix with a SmartPointer object. This class - // also appears in step-21 and a more advanced version of it appears in - // step-22. - template - class InverseMatrix : public Subscriptor - { - public: - InverseMatrix(const MatrixType &m); - - void vmult(Vector &dst, const Vector &src) const; - - private: - const SmartPointer matrix; - }; - - - template - InverseMatrix::InverseMatrix(const MatrixType &m) - : matrix(&m) - {} - + // The linear solvers and preconditioners we use in this example have + // been discussed in significant detail already in the introduction. We + // will therefore not discuss the rationale for our approach here any + // more, but rather only comment on some remaining implementational + // aspects. - template - void InverseMatrix::vmult(Vector & dst, - const Vector &src) const - { - // To make the control flow simpler, we recreate both the ReductionControl - // and SolverCG objects every time this is called. This is not the most - // efficient choice because SolverCG instances allocate memory whenever - // they are created; this is just a tutorial so such inefficiencies are - // acceptable for the sake of exposition. - 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()); - } - - - // @sect4{The SchurComplement class} + // @sect4{MixedLaplace::solve} - // The next class is the Schur complement class. Its rationale has also been - // discussed in length in the introduction. Like InverseMatrix, - // this class is derived from Subscriptor and stores SmartPointer s - // pointing to the system matrix and InverseMatrix wrapper. - // - // The vmult function requires two temporary vectors that we do - // not want to re-allocate and free every time we call this function. Since - // here, we have full control over the use of these vectors (unlike above, - // where a class called by the vmult function required these - // vectors, not the vmult function itself), we allocate them - // directly, rather than going through the VectorMemory - // mechanism. However, again, these member variables do not carry any state - // between successive calls to the member functions of this class (i.e., we - // never care what values they were set to the last time a member function - // was called), we mark these vectors as mutable. - // - // The rest of the (short) implementation of this class is straightforward - // if you know the order of matrix-vector multiplications performed by the - // vmult function: - class SchurComplement : public Subscriptor + // As already outlined in the introduction, the solve function consists + // essentially of two steps. First, we have to form the first equation + // involving the Schur complement and solve for the pressure (component 1 + // of the solution). Then, we can reconstruct the velocities from the + // second equation (component 0 of the solution). + template + void MixedLaplaceProblem::solve() { - public: - SchurComplement(const BlockSparseMatrix & A, - const InverseMatrix> &Minv); - - void vmult(Vector &dst, const Vector &src) const; - - private: - const SmartPointer> system_matrix; - const SmartPointer>> m_inverse; - - mutable Vector tmp1, tmp2; - }; + // As a first step we declare references to all block components of the + // matrix, the right hand side and the solution vector that we will + // need. + const auto &M = system_matrix.block(0, 0); + const auto &B = system_matrix.block(0, 1); + const auto &F = system_rhs.block(0); + const auto &G = system_rhs.block(1); - 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()) - {} + auto &U = solution.block(0); + auto &P = solution.block(1); + // Then, we will create corresponding LinearOperator objects and create + // the op_M_inv operator: - 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); - } + const auto op_M = linear_operator(M); + const auto op_B = linear_operator(B); + ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10); + SolverCG<> solver_M(reduction_control_M); + PreconditionJacobi<> preconditioner_M; - // @sect4{The ApproximateSchurComplement class} + preconditioner_M.initialize(M); - // The third component of our solver and preconditioner system is the class - // that approximates the Schur complement with the method described in the - // introduction. We will use this class to build a preconditioner for our - // system matrix. - class ApproximateSchurComplement : public Subscriptor - { - public: - ApproximateSchurComplement(const BlockSparseMatrix &A); + const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M); - void vmult(Vector &dst, const Vector &src) const; + // This puts us in the position to be able to declare the Schur + // complement op_S and the approximate Schur complement + // op_aS: - private: - const SmartPointer> system_matrix; + const auto op_S = transpose_operator(op_B) * op_M_inv * op_B; + const auto op_aS = + transpose_operator(op_B) * linear_operator(preconditioner_M) * op_B; - mutable Vector tmp1, tmp2; - }; + // We now create a preconditioner out of op_aS that + // applies a few number of CG iterations (until a very modest relative + // reduction of $10^{-16}$ is reached): + ReductionControl reduction_control_aS(2000, 1.e-18, 1.0e-6); + SolverCG<> solver_aS(reduction_control_aS); + PreconditionIdentity preconditioner_aS; + const auto preconditioner_S = + inverse_operator(op_aS, solver_aS, preconditioner_aS); - ApproximateSchurComplement::ApproximateSchurComplement( - const BlockSparseMatrix &A) - : system_matrix(&A) - , tmp1(A.block(0, 0).m()) - , tmp2(A.block(0, 0).m()) - {} + // Now on to the first equation. The right hand side of it is + // $B^TM^{-1}F-G$, which is what we compute in the first few lines. We + // then solve the first equation with a CG solver and the + // preconditioner we just declared. + const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G; - 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); - } + SolverControl solver_control_S(2000, 1.e-12); + SolverCG<> solver_S(solver_control_S); - // @sect4{MixedLaplace::solve} + const auto op_S_inv = inverse_operator(op_S, solver_S, preconditioner_S); - // After all these preparations, we can finally write the function that - // actually solves the linear problem. We will go through the two parts it - // has that each solve one of the two equations, the first one for the - // pressure (component 1 of the solution), then the velocities (component 0 - // of the solution). - template - void MixedLaplaceProblem::solve() - { - InverseMatrix> inverse_mass(system_matrix.block(0, 0)); - Vector tmp(solution.block(0).size()); + P = op_S_inv * schur_rhs; - // Now on to the first equation. The right hand side of it is - // $B^TM^{-1}F-G$, which is what we compute in the first few lines: - { - SchurComplement schur_complement(system_matrix, inverse_mass); - Vector schur_rhs(solution.block(1).size()); - inverse_mass.vmult(tmp, system_rhs.block(0)); - system_matrix.block(1, 0).vmult(schur_rhs, tmp); - schur_rhs -= system_rhs.block(1); - - // Now that we have the right hand side we can go ahead and solve for the - // pressure, using our approximation of the inverse as a preconditioner: - SolverControl solver_control(solution.block(1).size(), - 1e-12 * schur_rhs.l2_norm()); - SolverCG<> cg(solver_control); - - ApproximateSchurComplement approximate_schur(system_matrix); - InverseMatrix approximate_inverse( - approximate_schur); - cg.solve(schur_complement, - solution.block(1), - schur_rhs, - approximate_inverse); - - std::cout << solver_control.last_step() - << " CG Schur complement iterations to obtain convergence." - << std::endl; - } + std::cout << solver_control_S.last_step() + << " CG Schur complement iterations to obtain convergence." + << std::endl; // After we have the pressure, we can compute the velocity. The equation // reads $MU=-BP+F$, and we solve it by first computing the right hand // side, and then multiplying it with the object that represents the // inverse of the mass matrix: - { - system_matrix.block(0, 1).vmult(tmp, solution.block(1)); - tmp *= -1; - tmp += system_rhs.block(0); - inverse_mass.vmult(solution.block(0), tmp); - } + U = op_M_inv * (F - op_B * P); }