From: Daniel Arndt Date: Tue, 10 Apr 2018 22:41:31 +0000 (+0200) Subject: Store the pressure mass matrix separately in step-22 X-Git-Tag: v9.0.0-rc1~193^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ea4b87cb937d929c1428d3c3d0042933e3d25297;p=dealii.git Store the pressure mass matrix separately in step-22 --- diff --git a/examples/step-22/doc/intro.dox b/examples/step-22/doc/intro.dox index 32def7a482..c23e62d78e 100644 --- a/examples/step-22/doc/intro.dox +++ b/examples/step-22/doc/intro.dox @@ -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. diff --git a/examples/step-22/doc/results.dox b/examples/step-22/doc/results.dox index 52b4cfb5bb..122fb898f2 100644 --- a/examples/step-22/doc/results.dox +++ b/examples/step-22/doc/results.dox @@ -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 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 pmass_preconditioner; pmass_preconditioner.initialize (pressure_mass_matrix, SparseILU::AdditionalData()); diff --git a/examples/step-22/step-22.cc b/examples/step-22/step-22.cc index 56a3dc7d95..33fddfecaf 100644 --- a/examples/step-22/step-22.cc +++ b/examples/step-22/step-22.cc @@ -103,12 +103,14 @@ namespace Step22 // @sect3{The StokesProblem 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 hanging_node_constraints into - // constraints. + // are nearl same as used there. The only difference is that we have an + // additional member pressure_mass_matrixthat 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 + // hanging_node_constraints into constraints. template class StokesProblem { @@ -134,6 +136,8 @@ namespace Step22 BlockSparsityPattern sparsity_pattern; BlockSparseMatrix system_matrix; + SparseMatrix pressure_mass_matrix; + BlockVector solution; BlockVector 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 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 // phi_p[i] * phi_p[j] 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., // block(0,0) 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 preconditioner; - preconditioner.initialize (system_matrix.block(1,1), + preconditioner.initialize (pressure_mass_matrix, SparseILU::AdditionalData()); InverseMatrix,SparseILU > - 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