]> https://gitweb.dealii.org/ - dealii.git/commitdiff
examples/step-20: Actually use LinearOperator in the example step
authorMatthias Maier <tamiko@43-1.org>
Fri, 8 Feb 2019 19:51:36 +0000 (13:51 -0600)
committerMatthias Maier <tamiko@43-1.org>
Tue, 12 Feb 2019 00:46:55 +0000 (18:46 -0600)
examples/step-20/step-20.cc

index 8710392acae2a60350925c3e248d2c63d08560e2..20ddcb621d42806a67a750ed640087a02cfe2ebc 100644 (file)
@@ -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
  */
 
 
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/function.h>
+
 #include <deal.II/lac/block_vector.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/block_sparse_matrix.h>
 #include <deal.II/lac/solver_cg.h>
 #include <deal.II/lac/precondition.h>
 
+// The only two new header files that deserve some attention are those for
+// the LinearOperator and PackagedOperation classes:
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/packaged_operation.h>
+
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/tria_accessor.h>
@@ -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 <code>InverseMatrix</code> 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
-  // <code>Subscriptor</code> and, as mentioned above, it stores a pointer to
-  // the underlying matrix with a <code>SmartPointer</code> object. This class
-  // also appears in step-21 and a more advanced version of it appears in
-  // step-22.
-  template <class MatrixType>
-  class InverseMatrix : public Subscriptor
-  {
-  public:
-    InverseMatrix(const MatrixType &m);
-
-    void vmult(Vector<double> &dst, const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const MatrixType> matrix;
-  };
-
-
-  template <class MatrixType>
-  InverseMatrix<MatrixType>::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 <class MatrixType>
-  void InverseMatrix<MatrixType>::vmult(Vector<double> &      dst,
-                                        const Vector<double> &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<unsigned int>(src.size(), 200),
-                                 1e-8 * src.l2_norm());
-    SolverCG<>    cg(solver_control);
-
-    dst = 0;
-
-    cg.solve(*matrix, dst, src, PreconditionIdentity());
-  }
-
-
-  // @sect4{The <code>SchurComplement</code> class}
+  // @sect4{MixedLaplace::solve}
 
-  // The next class is the Schur complement class. Its rationale has also been
-  // discussed in length in the introduction. Like <code>InverseMatrix</code>,
-  // this class is derived from Subscriptor and stores SmartPointer&nbsp;s
-  // pointing to the system matrix and <code>InverseMatrix</code> wrapper.
-  //
-  // The <code>vmult</code> 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 <code>vmult</code> function required these
-  // vectors, not the <code>vmult</code> function itself), we allocate them
-  // directly, rather than going through the <code>VectorMemory</code>
-  // 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 <code>mutable</code>.
-  //
-  // The rest of the (short) implementation of this class is straightforward
-  // if you know the order of matrix-vector multiplications performed by the
-  // <code>vmult</code> 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 <int dim>
+  void MixedLaplaceProblem<dim>::solve()
   {
-  public:
-    SchurComplement(const BlockSparseMatrix<double> &          A,
-                    const InverseMatrix<SparseMatrix<double>> &Minv);
-
-    void vmult(Vector<double> &dst, const Vector<double> &src) const;
-
-  private:
-    const SmartPointer<const BlockSparseMatrix<double>>           system_matrix;
-    const SmartPointer<const InverseMatrix<SparseMatrix<double>>> m_inverse;
-
-    mutable Vector<double> 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<double> &          A,
-    const InverseMatrix<SparseMatrix<double>> &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 <code>op_M_inv</code> operator:
 
-  void SchurComplement::vmult(Vector<double> &      dst,
-                              const Vector<double> &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 <code>ApproximateSchurComplement</code> 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<double> &A);
+    const auto op_M_inv = inverse_operator(op_M, solver_M, preconditioner_M);
 
-    void vmult(Vector<double> &dst, const Vector<double> &src) const;
+    // This puts us in the position to be able to declare the Schur
+    // complement <code>op_S</code> and the approximate Schur complement
+    // <code>op_aS</code>:
 
-  private:
-    const SmartPointer<const BlockSparseMatrix<double>> 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<double> tmp1, tmp2;
-  };
+    // We now create a preconditioner out of <code>op_aS</code> 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<double> &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<double> &      dst,
-                                         const Vector<double> &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 <int dim>
-  void MixedLaplaceProblem<dim>::solve()
-  {
-    InverseMatrix<SparseMatrix<double>> inverse_mass(system_matrix.block(0, 0));
-    Vector<double>                      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<double>  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<ApproximateSchurComplement> 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);
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.