]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Get rid of IterativeInverse in step-20.
authorDavid Wells <wellsd2@rpi.edu>
Sat, 20 Feb 2016 23:41:57 +0000 (18:41 -0500)
committerDavid Wells <wellsd2@rpi.edu>
Sun, 21 Feb 2016 14:04:45 +0000 (09:04 -0500)
IterativeInverse was marked as deprecated in April 2015 in favor of
LinearOperator. This puts this step in a strange place; we should not
use deprecated features in tutorials, but at the same time we cannot use
its replacement, which requires C++11 support. To resolve this, this
step just provides its own replacement for a special case of
IterativeInverse.

examples/step-20/doc/intro.dox
examples/step-20/step-20.cc

index 63ca74a0f728a09d3251dfdfb3c6f32c20da34ce..9823e7097ff9b8f95301149ef3cd564e1ff83157 100644 (file)
@@ -434,26 +434,25 @@ function <code>vmult()</code> that is all that the SolverCG class
 requires from an object representing a matrix. We can make our life a
 bit easier by also introducing an object that represents $M^{-1}$ and
 that has its own <code>vmult()</code> function that, if called, solves
-the linear system with $M$; in fact, such a class already exists in
-deal.II: this is accomplished by using the class
-IterativeInverse. Using it, the class that implements the Schur
-only needs to offer the <code>vmult()</code>
+the linear system with $M$. Using this (which we will implement as the
+<code>InverseMatrix</code> class in the body of the program), the class
+that implements the Schur only needs to offer the <code>vmult()</code>
 function to perform a matrix-vector multiplication, using the algorithm
 above. Here are again the relevant parts of the code:
 
 @code
-class SchurComplement
+class SchurComplement : public Subscriptor
 {
   public:
-    SchurComplement (const BlockSparseMatrix<double> &A,
-                     const IterativeInverse<Vector<double> > &Minv);
+    SchurComplement (const BlockSparseMatrix<double>            &A,
+                     const InverseMatrix<SparseMatrix<double> > &inverse_mass);
 
     void vmult (Vector<double>       &dst,
                 const Vector<double> &src) const;
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const IterativeInverse<Vector<double> > > m_inverse;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double> > > inverse_mass;
 
     mutable Vector<double> tmp1, tmp2;
 };
@@ -463,7 +462,7 @@ void SchurComplement::vmult (Vector<double>       &dst,
                              const Vector<double> &src) const
 {
   system_matrix->block(0,1).vmult (tmp1, src); // multiply with the top right block: B
-  m_inverse->vmult (tmp2, tmp1);               // multiply with M^-1
+  inverse_mass->vmult (tmp2, tmp1);            // multiply with M^-1
   system_matrix->block(1,0).vmult (dst, tmp2); // multiply with the bottom left block: B^T
 }
 @endcode
@@ -490,37 +489,31 @@ matrix and the mass matrix, respectively:
 template <int dim>
 void MixedLaplaceProblem<dim>::solve ()
 {
-  PreconditionIdentity identity;
-  IterativeInverse<Vector<double> > m_inverse;
-  m_inverse.initialize(system_matrix.block(0,0), identity);
-  m_inverse.solver.select("cg");
-  static ReductionControl inner_control(1000, 0., 1.e-13);
-  m_inverse.solver.set_control(inner_control);
-
+  InverseMatrix<SparseMatrix<double> > inverse_mass (system_matrix.block(0,0));
   Vector<double> tmp (solution.block(0).size());
 
   {
+    SchurComplement schur_complement (system_matrix, inverse_mass);
     Vector<double> schur_rhs (solution.block(1).size());
-
-    m_inverse.vmult (tmp, system_rhs.block(0));
+    inverse_mass.vmult (tmp, system_rhs.block(0));
     system_matrix.block(1,0).vmult (schur_rhs, tmp);
     schur_rhs -= system_rhs.block(1);
 
-    SolverControl solver_control (system_matrix.block(0,0).m(),
-                                  1e-6*schur_rhs.l2_norm());
-    SolverCG<>    cg (solver_control);
+    SolverControl solver_control (solution.block(1).size(),
+                                  1e-12*schur_rhs.l2_norm());
+    SolverCG<> cg (solver_control);
 
-    cg.solve (SchurComplement(system_matrix, m_inverse),
-              solution.block(1),
-              schur_rhs,
-              PreconditionIdentity());
+    PreconditionIdentity preconditioner;
+    cg.solve (schur_complement, solution.block(1), schur_rhs,
+              preconditioner);
   }
+
   {
     system_matrix.block(0,1).vmult (tmp, solution.block(1));
     tmp *= -1;
     tmp += system_rhs.block(0);
 
-    m_inverse.vmult (solution.block(0), tmp);
+    inverse_mass.vmult (solution.block(0), tmp);
   }
 }
 @endcode
@@ -568,7 +561,7 @@ this looks almost as expensive as solving with $S$ right away. However, note
 that in the inner iteration, we do not have to calculate $M^{-1}$, but only
 the inverse of its diagonal, which is cheap.
 
-The next step is to define a class that represents the approximate Schur
+The next step is to define a class that represents this approximate Schur
 complement. This should look very much like the Schur complement class itself,
 except that it doesn't need the object representing $M^{-1}$ any more
 since we can compute the inverse of the diagonal of $M$ on the fly:
@@ -576,26 +569,28 @@ since we can compute the inverse of the diagonal of $M$ on the fly:
 @code
 class ApproximateSchurComplement : public Subscriptor
 {
-  public:
-    ApproximateSchurComplement (const BlockSparseMatrix<double> &A);
+public:
+  ApproximateSchurComplement (const BlockSparseMatrix<double> &A);
 
-    void vmult (Vector<double>       &dst,
-                const Vector<double> &src) const;
+  void vmult (Vector<double>       &dst,
+              const Vector<double> &src) const;
 
-  private:
-    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+private:
+  const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
 
-    mutable Vector<double> tmp1, tmp2;
+  mutable Vector<double> tmp1, tmp2;
 };
 
 
-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);
-}
+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);
+  }
 @endcode
 
 Note how the <code>vmult</code> function differs in simply doing one Jacobi sweep
@@ -606,18 +601,14 @@ the diagonal of $M$; in other words, the operation $({\textrm{diag}\
 }M)^{-1}x$ on a vector $x$ is exactly what the function
 SparseMatrix::precondition_Jacobi above does.)
 
-With all this, we already have the preconditioner: it should be the inverse of
-the approximate Schur complement, i.e. we need code like this:
+With all this, we nearly already have the preconditioner: it should be the
+inverse of the approximate Schur complement. We implement this with
+the <code>InverseMatrix</code> class:
 
 @code
-    ApproximateSchurComplement
-      approximate_schur_complement (system_matrix);
-
-    IterativeInverse<Vector<double> >
-      preconditioner;
-    preconditioner.initialize(approximate_schur_complement, identity);
-    preconditioner.solver.select("cg");
-    preconditioner.solver.set_control(inner_control);
+  ApproximateSchurComplement approximate_schur (system_matrix);
+  InverseMatrix<ApproximateSchurComplement> approximate_inverse
+  (approximate_schur);
 @endcode
 
 That's all!
@@ -626,30 +617,21 @@ Taken together, the first block of our <code>solve()</code> function will then
 look like this:
 
 @code
-    Vector<double> schur_rhs (solution.block(1).size());
-
-    m_inverse.vmult (tmp, system_rhs.block(0));
-    system_matrix.block(1,0).vmult (schur_rhs, tmp);
-    schur_rhs -= system_rhs.block(1);
-
-    SchurComplement
-      schur_complement (system_matrix, m_inverse);
-
-    ApproximateSchurComplement
-      approximate_schur_complement (system_matrix);
-
-    IterativeInverse<Vector<double> >
-      preconditioner;
-    preconditioner.initialize(approximate_schur_complement, identity);
-    preconditioner.solver.select("cg");
-    preconditioner.solver.set_control(inner_control);
-
-    SolverControl solver_control (solution.block(1).size(),
-                                  1e-12*schur_rhs.l2_norm());
-    SolverCG<> cg (solver_control);
-
-    cg.solve (schur_complement, solution.block(1), schur_rhs,
-              preconditioner);
+      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);
+
+      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);
 @endcode
 
 Note how we pass the so-defined preconditioner to the solver working on the
index 40260643a8cd2b722435fd9f502fa3f5401b287b..90fa8cc12710866f904a6b9198f51f13f55d1b89 100644 (file)
@@ -556,15 +556,62 @@ namespace Step20
   // 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.
+  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)
+  {}
+
+
+  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(src.size(), static_cast<std::size_t> (200)),
+                                  1e-8*src.l2_norm());
+    SolverCG<>    cg (solver_control);
+
+    dst = 0;
+
+    cg.solve (*matrix, dst, src, PreconditionIdentity());
+  }
+
 
   // @sect4{The <code>SchurComplement</code> class template}
 
   // The next class is the Schur complement class. Its rationale has also been
-  // discussed in length in the introduction. The only things we would like to
-  // note is that the class, too, is derived from the <code>Subscriptor</code>
-  // class and that as mentioned above it stores pointers to the entire block
-  // matrix and the inverse of the mass matrix block using
-  // <code>SmartPointer</code> objects.
+  // 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
@@ -583,25 +630,26 @@ namespace Step20
   class SchurComplement : public Subscriptor
   {
   public:
-    SchurComplement (const BlockSparseMatrix<double> &A,
-                     const IterativeInverse<Vector<double> > &Minv);
+    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 IterativeInverse<Vector<double> > > m_inverse;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double> > > m_inverse;
 
     mutable Vector<double> tmp1, tmp2;
   };
 
 
-  SchurComplement::SchurComplement (const BlockSparseMatrix<double> &A,
-                                    const IterativeInverse<Vector<double> > &Minv)
+  SchurComplement
+  ::SchurComplement (const BlockSparseMatrix<double>            &A,
+                     const InverseMatrix<SparseMatrix<double> > &Minv)
     :
     system_matrix (&A),
-    m_inverse (&Minv),
+    inverse_mass (&Minv),
     tmp1 (A.block(0,0).m()),
     tmp2 (A.block(0,0).m())
   {}
@@ -619,29 +667,9 @@ namespace Step20
   // @sect4{The <code>ApproximateSchurComplement</code> class template}
 
   // The third component of our solver and preconditioner system is the class
-  // that approximates the Schur complement so we can form a an InverseIterate
-  // object that approximates the inverse of the Schur complement. It follows
-  // the same pattern as the Schur complement class, with the only exception
-  // that we do not multiply with the inverse mass matrix in
-  // <code>vmult</code>, but rather just do a single Jacobi
-  // step. Consequently, the class also does not have to store a pointer to an
-  // inverse mass matrix object.
-  //
-  // We will later use this class as a template argument to the
-  // IterativeInverse class which will in turn want to use it as a
-  // template argument for the PointerMatrix class. The latter class
-  // has a function that requires us to also write a function that
-  // provides the product with the transpose of the matrix this object
-  // represents. As a consequence, in the code below, we also
-  // implement a <tt>Tvmult</tt> function here that represents the
-  // product of the transpose matrix with a vector. It is easy to see
-  // how this needs to be implemented here: since the matrix is
-  // symmetric, we can as well call <code>vmult</code> wherever the
-  // product with the transpose matrix is required. (Note, however,
-  // that even though we implement this function here, there will in
-  // fact not be any need for it as long as we use SolverCG as the
-  // solver since that solver does not ever call the function that
-  // provides this operation.)
+  // 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:
@@ -649,8 +677,6 @@ namespace Step20
 
     void vmult (Vector<double>       &dst,
                 const Vector<double> &src) const;
-    void Tvmult (Vector<double>       &dst,
-                 const Vector<double> &src) const;
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
@@ -659,85 +685,58 @@ namespace Step20
   };
 
 
-  ApproximateSchurComplement::ApproximateSchurComplement (const BlockSparseMatrix<double> &A)
-    :
+
+  ApproximateSchurComplement::ApproximateSchurComplement
+  (const BlockSparseMatrix<double> &A) :
     system_matrix (&A),
     tmp1 (A.block(0,0).m()),
     tmp2 (A.block(0,0).m())
   {}
 
 
-  void ApproximateSchurComplement::vmult (Vector<double>       &dst,
-                                          const Vector<double> &src) const
+  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);
   }
 
-
-  void ApproximateSchurComplement::Tvmult (Vector<double>       &dst,
-                                           const Vector<double> &src) const
-  {
-    vmult (dst, src);
-  }
-
-
-
   // @sect4{MixedLaplace::solve}
 
   // 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). Both parts need an object representing the inverse mass
-  // matrix and an auxiliary vector, and we therefore declare these objects at
-  // the beginning of this function.
+  // of the solution).
   template <int dim>
   void MixedLaplaceProblem<dim>::solve ()
   {
-    ReductionControl inner_control(1000, 0., 1.e-13);
-    PreconditionIdentity identity;
-    IterativeInverse<Vector<double> > m_inverse;
-    m_inverse.initialize(system_matrix.block(0,0), identity);
-    m_inverse.solver.select("cg");
-    m_inverse.solver.set_control(inner_control);
-
+    InverseMatrix<SparseMatrix<double> > inverse_mass (system_matrix.block(0,0));
     Vector<double> tmp (solution.block(0).size());
 
     // 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 declare the
-    // objects representing the Schur complement, its approximation, and the
-    // inverse of the approximation. Finally, we declare a solver object and
-    // hand off all these matrices and vectors to it to compute block 1 (the
-    // pressure) of the solution:
+    // 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());
-
-      m_inverse.vmult (tmp, system_rhs.block(0));
+      inverse_mass.vmult (tmp, system_rhs.block(0));
       system_matrix.block(1,0).vmult (schur_rhs, tmp);
       schur_rhs -= system_rhs.block(1);
 
-
-      SchurComplement
-      schur_complement (system_matrix, m_inverse);
-
-      ApproximateSchurComplement
-      approximate_schur_complement (system_matrix);
-
-      IterativeInverse<Vector<double> >
-      preconditioner;
-      preconditioner.initialize(approximate_schur_complement, identity);
-      preconditioner.solver.select("cg");
-      preconditioner.solver.set_control(inner_control);
-
-
+      // 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);
+      SolverCG<> cg (solver_control);
 
+      ApproximateSchurComplement approximate_schur (system_matrix);
+      InverseMatrix<ApproximateSchurComplement> approximate_inverse
+      (approximate_schur);
       cg.solve (schur_complement, solution.block(1), schur_rhs,
-                preconditioner);
+                approximate_inverse);
 
       std::cout << solver_control.last_step()
                 << " CG Schur complement iterations to obtain convergence."
@@ -753,7 +752,7 @@ namespace Step20
       tmp *= -1;
       tmp += system_rhs.block(0);
 
-      m_inverse.vmult (solution.block(0), tmp);
+      inverse_mass.vmult (solution.block(0), tmp);
     }
   }
 

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.