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.
// @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
{
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
+ SparseMatrix<double> pressure_mass_matrix;
+
BlockVector<double> solution;
BlockVector<double> system_rhs;
// 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,
{
A_preconditioner.reset ();
system_matrix.clear ();
+ pressure_mass_matrix.clear ();
dof_handler.distribute_dofs (fe);
DoFRenumbering::Cuthill_McKee (dof_handler);
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);
{
system_matrix=0;
system_rhs=0;
+ pressure_mass_matrix = 0.;
QGauss<dim> quadrature_formula(degree+2);
- 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
// 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
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
// 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
// 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
projects / dealii.git / commitdiff