]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Store the pressure mass matrix separately in step-22
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Tue, 10 Apr 2018 22:41:31 +0000 (00:41 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Tue, 10 Apr 2018 22:46:14 +0000 (00:46 +0200)
examples/step-22/doc/intro.dox
examples/step-22/doc/results.dox
examples/step-22/step-22.cc

index 32def7a4821a0d81d44e2e90d2e003409d6d43c9..c23e62d78e40c7fe97ae35a83090a1fe9d534214 100644 (file)
@@ -473,33 +473,32 @@ program to see that indeed the number of CG iterations does not change as we
 refine the mesh).
 
 So all we need in addition to what we already have is the mass matrix on the
-pressure variables. We could do that by building this matrix on the
-side in a separate data structure. However, it is worth remembering
-that although we build the system matrix
+pressure variables and we will store it in a separate object. However, it is
+worth remembering that we never actually do matrix-vector products or any other
+operations that considers all blocks with the entire system matrix
 @f{eqnarray*}
   \left(\begin{array}{cc}
-    A & B^T \\ B & 0
+    A & B^T \\ B & 0.
   \end{array}\right)
 @f}
-as one object (of type BlockSparseMatrix), we never actually do
-matrix-vector products with this matrix, or any other operations that
-consider the entire matrix. Rather, we only build it in this form for
-convenience (because it reflects the structure of the FESystem finite
-element and associated DoFHandler object) but later only operate on
-the $(0,0),(0,1)$, and $(1,0)$ blocks of this matrix. In other words,
-our algorithm so far entirely ignores the $(1,1)$ (pressure-pressure)
-block as it is empty anyway.
-
-Now, as mentioned, we need a pressure mass matrix to precondition the
-Schur complement and that conveniently the pressure-pressure block of
-the matrix we build anyway is currently empty and ignored. So what we
-will do is to assemble the needed mass matrix in this space; this does
-change the global system matrix but since our algorithm never operates
-on the global matrix and instead only considers individual blocks,
-this fact does not affect what we actually compute. Later, when
-solving, we then precondition the Schur complement with $M_p^{-1}$ by
-doing a few CG iterations on the well-conditioned pressure mass matrix
-$M_p$ stored in the $(1,1)$ block.
+Rather, we only build it in this form for convenience (because it reflects the
+structure of the FESystem finite element and associated DoFHandler object) but
+later only operate on the $(0,0),(0,1)$, and $(1,0)$ blocks of this matrix. In
+other words, our algorithm so far entirely ignores the $(1,1)$
+(pressure-pressure) block as it is empty anyway.
+
+Now, as mentioned, we need a pressure mass matrix to precondition the Schur
+complement and that conveniently the pressure-pressure block of the matrix we
+build anyway is currently empty and ignored. So what we will do is to assemble
+the needed mass matrix in this space; this does change the global system matrix
+but since our algorithm never operates on the global matrix and instead only
+considers individual blocks, this fact does not affect what we actually compute.
+Nevertheless, we decided for clafity to copy this block after assembling in a
+separate object and clear the (1,1) block again. This way the system matrix
+could also be used in its entirety, e.g. in a direct solver. Later, when
+solving, we then precondition the Schur complement with $M_p^{-1}$ by doing
+a few CG iterations on the well-conditioned pressure mass matrix $M_p$
+stored separately.
 
 
 
index 52b4cfb5bb81695bb994237b5e8c3197ac1250e3..122fb898f2c1fb364e9f57b2bab72656def53121 100644 (file)
@@ -537,11 +537,6 @@ With this said, we turn to the realization of the solver call with GMRES with
 $k=100$ temporary vectors:
 
 @code
-      SparseMatrix<double> pressure_mass_matrix;
-      pressure_mass_matrix.reinit(sparsity_pattern.block(1,1));
-      pressure_mass_matrix.copy_from(system_matrix.block(1,1));
-      system_matrix.block(1,1) = 0;
-
       SparseILU<double> pmass_preconditioner;
       pmass_preconditioner.initialize (pressure_mass_matrix,
         SparseILU<double>::AdditionalData());
index 56a3dc7d9581c0adf6eb67d231c3583cf762dd5b..33fddfecaf5addcdb6fd4a7cd95c70a48d8b03cb 100644 (file)
@@ -103,12 +103,14 @@ namespace Step22
   // @sect3{The <code>StokesProblem</code> class template}
 
   // This is an adaptation of step-20, so the main class and the data types
-  // are the same as used there. In this example we also use adaptive grid
-  // refinement, which is handled in analogy to step-6. According to the
-  // discussion in the introduction, we are also going to use the
-  // ConstraintMatrix for implementing Dirichlet boundary conditions. Hence,
-  // we change the name <code>hanging_node_constraints</code> into
-  // <code>constraints</code>.
+  // are nearl same as used there. The only difference is that we have an
+  // additional member <code>pressure_mass_matrix</code>that is used for
+  // preconditioning the Schur complement.
+  // In this example we also use adaptive grid refinement, which is handled
+  // in analogy to step-6. According to the discussion in the introduction,
+  // we are also going to use the ConstraintMatrix for implementing Dirichlet
+  // boundary conditions. Hence, we change the name
+  // <code>hanging_node_constraints</code> into <code>constraints</code>.
   template <int dim>
   class StokesProblem
   {
@@ -134,6 +136,8 @@ namespace Step22
     BlockSparsityPattern      sparsity_pattern;
     BlockSparseMatrix<double> system_matrix;
 
+    SparseMatrix<double>      pressure_mass_matrix;
+
     BlockVector<double> solution;
     BlockVector<double> system_rhs;
 
@@ -401,8 +405,8 @@ namespace Step22
   // releases the pointer to the preconditioner object (if the shared pointer
   // pointed at anything at all at this point) since it will definitely not be
   // needed any more after this point and will have to be re-computed after
-  // assembling the matrix, and unties the sparse matrix from its sparsity
-  // pattern object.
+  // assembling the matrix, and unties the sparse matrices from their sparsity
+  // pattern objects.
   //
   // We then proceed with distributing degrees of freedom and renumbering
   // them: In order to make the ILU preconditioner (in 3D) work efficiently,
@@ -442,6 +446,7 @@ namespace Step22
   {
     A_preconditioner.reset ();
     system_matrix.clear ();
+    pressure_mass_matrix.clear ();
 
     dof_handler.distribute_dofs (fe);
     DoFRenumbering::Cuthill_McKee (dof_handler);
@@ -540,9 +545,11 @@ namespace Step22
       sparsity_pattern.copy_from (dsp);
     }
 
-    // Finally, the system matrix, solution and right hand side are created
-    // from the block structure as in step-20:
+    // Finally, the system matrix, the pressure mass matrix, the solution and
+    // the right hand side vector are created from the block structure
+    // similar to the approach in step-20:
     system_matrix.reinit (sparsity_pattern);
+    pressure_mass_matrix.reinit (sparsity_pattern.block(1,1));
 
     solution.reinit (2);
     solution.block(0).reinit (n_u);
@@ -567,6 +574,7 @@ namespace Step22
   {
     system_matrix=0;
     system_rhs=0;
+    pressure_mass_matrix = 0.;
 
     QGauss<dim>   quadrature_formula(degree+2);
 
@@ -651,7 +659,6 @@ namespace Step22
                                           - phi_p[i] * div_phi_u[j]
                                           + phi_p[i] * phi_p[j])
                                          * fe_values.JxW(q);
-
                   }
 
                 // For the right-hand side we use the fact that the shape
@@ -680,7 +687,8 @@ namespace Step22
         // discussed in the introduction. That this term only ends up in the
         // $(1,1)$ block stems from the fact that both of the factors in
         // <code>phi_p[i] * phi_p[j]</code> are only non-zero when all the
-        // other terms vanish (and the other way around).
+        // other terms vanish (and the other way around). This block will
+        // finally end up in the pressure_mass_matrix object.
         //
         // Note also that operator* is overloaded for symmetric tensors,
         // yielding the scalar product between the two tensors in the first
@@ -705,6 +713,13 @@ namespace Step22
                                                 system_matrix, system_rhs);
       }
 
+    // As dicussed previously, we built a pressure mass matrix in the $(1,1)$
+    // block of the system matrix. Now, correct for this by copying the block
+    // to the pressure_mass_matrix object and clear it in the system matrix
+    // afterwards.
+    pressure_mass_matrix.copy_from(system_matrix.block(1,1));
+    system_matrix.block(1,1) = 0.;
+
     // Before we're going to solve this linear system, we generate a
     // preconditioner for the velocity-velocity matrix, i.e.,
     // <code>block(0,0)</code> in the system matrix. As mentioned above, this
@@ -762,8 +777,7 @@ namespace Step22
 
       // Now to the preconditioner to the Schur complement. As explained in
       // the introduction, the preconditioning is done by a mass matrix in the
-      // pressure variable.  It is stored in the $(1,1)$ block of the system
-      // matrix (that is not used anywhere else but in preconditioning).
+      // pressure variable.
       //
       // Actually, the solver needs to have the preconditioner in the form
       // $P^{-1}$, so we need to create an inverse operation. Once again, we
@@ -784,11 +798,11 @@ namespace Step22
       // 1.2. It needs about twice the number of iterations, but the costs for
       // its generation are almost negligible.
       SparseILU<double> preconditioner;
-      preconditioner.initialize (system_matrix.block(1,1),
+      preconditioner.initialize (pressure_mass_matrix,
                                  SparseILU<double>::AdditionalData());
 
       InverseMatrix<SparseMatrix<double>,SparseILU<double> >
-      m_inverse (system_matrix.block(1,1), preconditioner);
+      m_inverse (pressure_mass_matrix, preconditioner);
 
       // With the Schur complement and an efficient preconditioner at hand, we
       // can solve the respective equation for the pressure (i.e. block 0 in

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.