From bfd95042549d77c173cc3978adef2460ac728f0b Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 11 Apr 2018 18:32:50 +0200 Subject: [PATCH] Adapt Wolfgang's approach --- examples/step-22/doc/intro.dox | 27 +----- examples/step-22/doc/results.dox | 2 + examples/step-22/step-22.cc | 150 +++++++++++++++++++------------ 3 files changed, 95 insertions(+), 84 deletions(-) diff --git a/examples/step-22/doc/intro.dox b/examples/step-22/doc/intro.dox index c23e62d78e..d1d28f9909 100644 --- a/examples/step-22/doc/intro.dox +++ b/examples/step-22/doc/intro.dox @@ -473,32 +473,7 @@ 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 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. - \end{array}\right) -@f} -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. +pressure variables and we will store it in a separate object. diff --git a/examples/step-22/doc/results.dox b/examples/step-22/doc/results.dox index 122fb898f2..a1f44fec28 100644 --- a/examples/step-22/doc/results.dox +++ b/examples/step-22/doc/results.dox @@ -537,6 +537,8 @@ With this said, we turn to the realization of the solver call with GMRES with $k=100$ temporary vectors: @code + const SparseMatrix &pressure_mass_matrix + = preconditioner_matrix.block(1,1); 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 33fddfecaf..cb9b1e1055 100644 --- a/examples/step-22/step-22.cc +++ b/examples/step-22/step-22.cc @@ -103,9 +103,10 @@ namespace Step22 // @sect3{The StokesProblem class template} // This is an adaptation of step-20, so the main class and the data types - // 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. + // are nearly the same as used there. The only difference is that we have an + // additional member preconditioner_matrix, that is used for + // preconditioning the Schur complement, and a corresponding sparsity pattern + // preconditioner_sparsity_pattern. // 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 @@ -136,7 +137,8 @@ namespace Step22 BlockSparsityPattern sparsity_pattern; BlockSparseMatrix system_matrix; - SparseMatrix pressure_mass_matrix; + BlockSparsityPattern preconditioner_sparsity_pattern; + BlockSparseMatrix preconditioner_matrix; BlockVector solution; BlockVector system_rhs; @@ -446,7 +448,7 @@ namespace Step22 { A_preconditioner.reset (); system_matrix.clear (); - pressure_mass_matrix.clear (); + preconditioner_matrix.clear (); dof_handler.distribute_dofs (fe); DoFRenumbering::Cuthill_McKee (dof_handler); @@ -510,27 +512,30 @@ namespace Step22 << " (" << n_u << '+' << n_p << ')' << std::endl; - // The next task is to allocate a sparsity pattern for the system matrix - // we will create. We could do this in the same way as in step-20, - // i.e. directly build an object of type SparsityPattern through - // DoFTools::make_sparsity_pattern. However, there is a major reason not - // to do so: In 3D, the function DoFTools::max_couplings_between_dofs - // yields a conservative but rather large number for the coupling between - // the individual dofs, so that the memory initially provided for the - // creation of the sparsity pattern of the matrix is far too much -- so - // much actually that the initial sparsity pattern won't even fit into the - // physical memory of most systems already for moderately-sized 3D - // problems, see also the discussion in step-18. Instead, we first build - // a temporary object that uses a different data structure that doesn't - // require allocating more memory than necessary but isn't suitable for - // use as a basis of SparseMatrix or BlockSparseMatrix objects; in a - // second step we then copy this object into an object of - // BlockSparsityPattern. This is entirely analogous to what we already did - // in step-11 and step-18. + // The next task is to allocate a sparsity pattern for the system matrix we + // will create and one for the preconditioner matrix. We could do this in + // the same way as in step-20, i.e. directly build an object of type + // SparsityPattern through DoFTools::make_sparsity_pattern. However, there + // is a major reason not to do so: + // In 3D, the function DoFTools::max_couplings_between_dofs yields a + // conservative but rather large number for the coupling between the + // individual dofs, so that the memory initially provided for the creation + // of the sparsity pattern of the matrix is far too much -- so much actually + // that the initial sparsity pattern won't even fit into the physical memory + // of most systems already for moderately-sized 3D problems, see also the + // discussion in step-18. Instead, we first build temporary objects that use + // a different data structure that doesn't require allocating more memory + // than necessary but isn't suitable for use as a basis of SparseMatrix or + // BlockSparseMatrix objects; in a second step we then copy these objects + // into objects of type BlockSparsityPattern. This is entirely analogous to + // what we already did in step-11 and step-18. In particular, we make use of + // the fact that we will never write into the $(1,1)$ block of the system + // matrix and that this is the only block to be filled for the + // preconditioner matrix. // - // All this is done inside a new scope, which - // means that the memory of dsp will be released once the - // information has been copied to sparsity_pattern. + // All this is done inside new scopes, which means that the memory of + // dsp will be released once the information has been copied to + // sparsity_pattern. { BlockDynamicSparsityPattern dsp (2,2); @@ -541,15 +546,50 @@ namespace Step22 dsp.collect_sizes(); - DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false); + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); + + for (unsigned int c=0; c preconditioner_coupling (dim+1, dim+1); + + for (unsigned int c=0; c void StokesProblem::assemble_system () { system_matrix=0; system_rhs=0; - pressure_mass_matrix = 0.; + preconditioner_matrix = 0; QGauss quadrature_formula(degree+2); @@ -589,6 +629,7 @@ namespace Step22 const unsigned int n_q_points = quadrature_formula.size(); FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix local_preconditioner_matrix (dofs_per_cell, dofs_per_cell); Vector local_rhs (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); @@ -636,6 +677,7 @@ namespace Step22 { fe_values.reinit (cell); local_matrix = 0; + local_preconditioner_matrix = 0; local_rhs = 0; right_hand_side.vector_value_list(fe_values.get_quadrature_points(), @@ -656,9 +698,11 @@ namespace Step22 { local_matrix(i,j) += (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) - div_phi_u[i] * phi_p[j] - - phi_p[i] * div_phi_u[j] - + phi_p[i] * phi_p[j]) + - phi_p[i] * div_phi_u[j]) * fe_values.JxW(q); + + local_preconditioner_matrix(i,j) += (phi_p[i] * phi_p[j]) + * fe_values.JxW(q); } // For the right-hand side we use the fact that the shape @@ -681,16 +725,7 @@ namespace Step22 } } - // Note that in the above computation of the local matrix contribution - // we added the term phi_p[i] * phi_p[j] , yielding a - // pressure mass matrix in the $(1,1)$ block of the matrix as - // 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). This block will - // finally end up in the pressure_mass_matrix object. - // - // Note also that operator* is overloaded for symmetric tensors, + // Note that operator* is overloaded for symmetric tensors, // yielding the scalar product between the two tensors in the first // line of the local matrix contribution. @@ -698,28 +733,27 @@ namespace Step22 // simultaneously use the ConstraintMatrix object to apply Dirichlet // boundary conditions and eliminate hanging node constraints, as we // discussed in the introduction), we have to be careful about one - // thing, though. We have only built half of the local matrix - // because of symmetry, but we're going to save the full system matrix - // in order to use the standard functions for solution. This is done + // thing, though. We have only built half of the local matrices + // because of symmetry, but we're going to save the full matrices + // in order to use the standard functions for solving. This is done // by flipping the indices in case we are pointing into the empty part - // of the local matrix. + // of the local matrices. for (unsigned int i=0; iget_dof_indices (local_dof_indices); constraints.distribute_local_to_global (local_matrix, local_rhs, local_dof_indices, system_matrix, system_rhs); + constraints.distribute_local_to_global (local_preconditioner_matrix, + local_dof_indices, + preconditioner_matrix); } - // 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 @@ -798,11 +832,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 (pressure_mass_matrix, + preconditioner.initialize (preconditioner_matrix.block(1,1), SparseILU::AdditionalData()); InverseMatrix,SparseILU > - m_inverse (pressure_mass_matrix, preconditioner); + m_inverse (preconditioner_matrix.block(1,1), 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 -- 2.39.5