From a595c36d9aa8d47461123f7958fd5a0c6baf65d0 Mon Sep 17 00:00:00 2001 From: Ryan Grove Date: Tue, 26 Apr 2016 22:22:52 -0400 Subject: [PATCH] Addressed most of Wolfgang's comments with the exception of "const" --- examples/step-55/CMakeLists.txt | 2 +- examples/step-55/doc/intro.dox | 174 ++++++++++++++++++++++--------- examples/step-55/doc/results.dox | 8 +- examples/step-55/doc/tooltip | 2 +- examples/step-55/step-55.cc | 114 +++++++++++--------- 5 files changed, 196 insertions(+), 104 deletions(-) diff --git a/examples/step-55/CMakeLists.txt b/examples/step-55/CMakeLists.txt index 5a32dc83cd..eef6addfe9 100644 --- a/examples/step-55/CMakeLists.txt +++ b/examples/step-55/CMakeLists.txt @@ -23,7 +23,7 @@ SET(TARGET_SRC CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) -FIND_PACKAGE(deal.II 8.3 QUIET +FIND_PACKAGE(deal.II 8.4 QUIET HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} ) IF(NOT ${deal.II_FOUND}) diff --git a/examples/step-55/doc/intro.dox b/examples/step-55/doc/intro.dox index f4923c9c77..600d636203 100644 --- a/examples/step-55/doc/intro.dox +++ b/examples/step-55/doc/intro.dox @@ -1,54 +1,104 @@ - +

Introduction

Stokes Problem

-The purpose of this tutorial is to create an efficient linear solver for the Stokes equation and compare it to alternative approaches. Using FGMRES with geometric multigrid as a precondtioner for the velocity block, we see that the linear solvers used in Step-22 cannot keep up since multigrid is the only way to get $O(n)$ solve time. Using the Timer class, we collect some statistics to compare setup times, solve times, and number of iterations. We also compute errors to make sure that what we have implemented is correct. - -Let $u \in H_0^1 = \{ u \in H^1(\Omega), u|_{\partial \Omega} = 0 \}$ and $p \in L_*^2 = \{ p_f \in L^2(\Omega), \int_\Omega p_f = 0 \}$. The Stokes equations read as follows in non-dimensionalized form: - -@f{eqnarray*} - - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p & =& f \\ - - \nabla \cdot u &=& 0 -@f} - -Note that we are using the deformation tensor instead of $\Delta u$ (a detailed desription of the difference between the two can be found in Step-22, but in summary, the deformation tensor is more physical as well as more expensive). - -

Linear solver and preconditioning issues

-The weak form of the discrete equations naturally leads to the following linear system for the nodal values of the velocity and pressure fields: +The purpose of this tutorial is to create an efficient linear solver +for the Stokes equation and compare it to alternative +approaches. Using FGMRES with geometric multigrid as a precondtioner +for the velocity block, we see that the linear solvers used in Step-22 +cannot keep up since multigrid is the only way to get $O(n)$ solve +time. Using the Timer class, we collect some statistics to compare +setup times, solve times, and number of iterations. We also compute +errors to make sure that what we have implemented is correct. + +Let $u \in H_0^1 = \{ u \in H^1(\Omega), u|_{\partial \Omega} = 0 \}$ +and $p \in L_*^2 = \{ p_f \in L^2(\Omega), \int_\Omega p_f = 0 +\}$. The Stokes equations read as follows in non-dimensionalized form: + +@f{eqnarray*} - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + +(\nabla \textbf{u})^T\right] + \nabla p & =& f \\ - \nabla \cdot u &=& +0 @f} + +Note that we are using the deformation tensor instead of $\Delta u$ (a +detailed desription of the difference between the two can be found in +Step-22, but in summary, the deformation tensor is more physical as +well as more expensive). + +

Linear solver and preconditioning issues

The weak form of +the discrete equations naturally leads to the following linear system +for the nodal values of the velocity and pressure fields: @f{eqnarray*} -\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) \left(\begin{array}{c} U \\ P \end{array}\right) = \left(\begin{array}{c} F \\ 0 \end{array}\right), -@f} +\left(\begin{array}{cc} A & B^T \\ B & 0 +\end{array}\right) \left(\begin{array}{c} U \\ P \end{array}\right) = +\left(\begin{array}{c} F \\ 0 \end{array}\right), @f} + +Our goal is to compare several solver approaches. In contrast to the +way in which step-22 solves the Stokes equation, we instead attack the +block system at once using a direct solver or FMGRES with an efficient +preconditioner. The idea is as follows: if we find a block +preconditioner $P$ such that the matrix @f{eqnarray*} +\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) P^{-1} @f} + +is simple, then an iterative solver with that preconditioner will +converge in a few iterations. Notice that we are doing right +preconditioning for this. Using the Schur complement $S=BA^{-1}B^T$, +we find that -Our goal is to compare several solver approaches. In contrast to the way in which step-22 solves the Stokes equation, we instead attack the block system at once using a direct solver or FMGRES with an efficient preconditioner. The idea is as follows: if we find a block preconditioner $P$ such that the matrix @f{eqnarray*} -\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) P^{-1} -@f} - -is simple, then an iterative solver with that preconditioner will converge in a few iterations. Notice that we are doing right preconditioning for this. Using the Schur complement $S=BA^{-1}B^T$, we find that - -@f{eqnarray*} - P^{-1} = \left(\begin{array}{cc} \hat{A} & B^T \\ 0 & \hat{S} \end{array}\right)^{-1} -@f} +P^{-1} = \left(\begin{array}{cc} \hat{A} & B^T \\ 0 & + \hat{S} \end{array}\right)^{-1} @f} is a good choice. It is important to note that -@f{eqnarray*} - P = \left(\begin{array}{cc} A^{-1} & 0 \\ 0 & I \end{array}\right) \left(\begin{array}{cc} I & B^T \\ 0 & -I \end{array}\right) \left(\begin{array}{cc} I & 0 \\ 0 & S^{-1} \end{array}\right). -@f} - -Since $P$ is aimed to be a preconditioner only, we shall use approximations to the inverse of the Schur complement $S$ and the matrix $A$. Therefore, in the above equations, $-M_p=\hat{S} \approx S$, where $M_p$ is the pressure mass matrix and is solved by using CG + ILU, and $\hat{A}$ is an approximation of $A$ obtained by one of multiple methods: CG + ILU, just using ILU, CG + GMG (Geometric Multigrid as described in Step-16), or just performing a few V-cycles of GMG. The inclusion of CG is more expensive, in general. - -As a comparison, instead of FGMRES, we also use the direct solver UMFPACK to compare our results to. If you want to use UMFPACK as a solver, it is important to note that since you have a singular system (since the integral of mean pressure being equal to zero not implemented), we set the first pressure node equal to zero since the direct solver can not handle the singular system like the other methods could. +@f{eqnarray*} +P = +\left(\begin{array}{cc} A^{-1} & 0 \\ 0 & I \end{array}\right) +\left(\begin{array}{cc} I & B^T \\ 0 & -I \end{array}\right) +\left(\begin{array}{cc} I & 0 \\ 0 & S^{-1} \end{array}\right). @f} + +Since $P$ is aimed to be a preconditioner only, we shall use +approximations to the inverse of the Schur complement $S$ and the +matrix $A$. Therefore, in the above equations, $-M_p=\hat{S} \approx +S$, where $M_p$ is the pressure mass matrix and is solved by using CG ++ ILU, and $\hat{A}$ is an approximation of $A$ obtained by one of +multiple methods: CG + ILU, just using ILU, CG + GMG (Geometric +Multigrid as described in Step-16), or just performing a few V-cycles +of GMG. The inclusion of CG is more expensive, in general. + +As a comparison, instead of FGMRES, we also use the direct solver +UMFPACK to compare our results to. If you want to use UMFPACK as a +solver, it is important to note that since you have a singular system +(since the integral of mean pressure being equal to zero not +implemented), we set the first pressure node equal to zero since the +direct solver can not handle the singular system like the other +methods could.

Reference Solution

-The domain, right hand side, and boundary conditions we implemented were chosen for their simplicity and the fact that they made it possible for us to compute errors using a reference solution. We apply Dirichlet boundary condtions for the whole velocity on the whole boundary of the domain Ω=[0,1]×[0,1]×[0,1]. To enforce the boundary conditions we can just use our reference solution that we will now define. - -Let $u=(u_1,u_2,u_3)=(2\sin (\pi x), - \pi y \cos (\pi x),- \pi z \cos (\pi x))$ and $p = \sin (\pi x)\cos (\pi y)\sin (\pi z)$. - -If you look up in the deal.ii manual what is needed to create a class inherited from Function@, you will find not only a value function, but vector_value, value_list, etc. Different things you use in your code will require one of these particular functions. This can be confusing at first, but luckily the only thing you actually need to implement is value. The other ones have default implementations inside deal.ii and will be called on their own as long as you implement value correctly. - -Notice that our reference solution fulfills $\nabla \cdot u = 0$. In addition, the pressure is chosen to have a mean value of zero. For the Method of Manufactured Solutions of Step-7, we need to find $\bf f$ such that: +The domain, right hand side, and boundary conditions we implemented +were chosen for their simplicity and the fact that they made it +possible for us to compute errors using a reference solution. We +apply Dirichlet boundary condtions for the whole velocity on the whole +boundary of the domain Ω=[0,1]×[0,1]×[0,1]. To enforce the boundary +conditions we can just use our reference solution that we will now +define. + +Let $u=(u_1,u_2,u_3)=(2\sin (\pi x), - \pi y \cos (\pi x),- \pi z \cos +(\pi x))$ and $p = \sin (\pi x)\cos (\pi y)\sin (\pi z)$. + +If you look up in the deal.ii manual what is needed to create a class +inherited from Function@, you will find not only a +value function, but vector_value, value_list, etc. Different things +you use in your code will require one of these particular +functions. This can be confusing at first, but luckily the only thing +you actually need to implement is value. The other ones have default +implementations inside deal.ii and will be called on their own as long +as you implement value correctly. + +Notice that our reference solution fulfills $\nabla \cdot u = 0$. In +addition, the pressure is chosen to have a mean value of zero. For +the Method of Manufactured Solutions of Step-7, we need to find $\bf +f$ such that: @f{align*} {\bf f} = - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p. @@ -56,16 +106,36 @@ Notice that our reference solution fulfills $\nabla \cdot u = 0$. In addition, t Using the reference solution above, we obtain: -@f{eqnarray*} -{\bf f} &=& (2 \pi^2 \sin (\pi x),- \pi^3 y \cos(\pi x),- \pi^3 z \cos(\pi x))\\ -& & + (\pi \cos(\pi x) \cos(\pi y) \sin(\pi z) ,- \pi \sin(\pi y) \sin(\pi x) \sin(\pi z), \pi \cos(\pi z) \sin(\pi x) \cos(\pi y)) -@f} - -

Computing Errors

-Because we do not enforce the mean pressure to be zero for our numerical solution in the linear system, we need to postprocess the solution after solving. To do this we use the compute_mean_value function to compute the mean value of the pressure to subtract it from the pressure. - -

DoF Handlers

-Geometric multigrid needs to know about the finite element system for the velocity. Since this is now part of the entire system, it is no longer easy to access. The reason for this is that there is currently no way in deal.ii to ask, "May I have just part of a DoF handler?" So in order to answer this request for our needs, we have to create a new DoF handler for just the velocites and assure that it has the same ordering as the DoF Handler for the entire system so that you can copy over one to the other. - -

Differences from Step-22

-The main difference between Step-55 and Step-22 is that we use block solvers instead of the Schur Complement approach used in step-22. Details of this approach can be found under the Block Schur complement preconditioner subsection of the Possible Extensions section of Step-22. For the preconditioner of the velocity block, we borrow a class from ASPECT called BlockSchurPreconditioner that has the option to solve for the inverse of $A$ or just apply one preconditioner sweep for it instead, which provides us with an expensive and cheap approach, respectively. \ No newline at end of file +@f{eqnarray*} +{\bf f} &=& (2 \pi^2 \sin (\pi x),- \pi^3 y \cos(\pi +x),- \pi^3 z \cos(\pi x))\\ & & + (\pi \cos(\pi x) \cos(\pi y) +\sin(\pi z) ,- \pi \sin(\pi y) \sin(\pi x) \sin(\pi z), \pi \cos(\pi +z) \sin(\pi x) \cos(\pi y)) @f} + +

Computing Errors

+Because we do not enforce the mean +pressure to be zero for our numerical solution in the linear system, +we need to postprocess the solution after solving. To do this we use +the compute_mean_value function to compute the mean value +of the pressure to subtract it from the pressure. + +

DoF Handlers

+Geometric multigrid needs to know about the +finite element system for the velocity. Since this is now part of the +entire system, it is no longer easy to access. The reason for this is +that there is currently no way in deal.ii to ask, "May I have just +part of a DoF handler?" So in order to answer this request for our +needs, we have to create a new DoF handler for just the velocites and +assure that it has the same ordering as the DoF Handler for the entire +system so that you can copy over one to the other. + +

Differences from Step-22

+The main difference between +Step-55 and Step-22 is that we use block solvers instead of the Schur +Complement approach used in step-22. Details of this approach can be +found under the Block Schur complement preconditioner subsection of +the Possible Extensions section of Step-22. For the preconditioner of +the velocity block, we borrow a class from ASPECT called +BlockSchurPreconditioner that has the option to solve for the inverse +of $A$ or just apply one preconditioner sweep for it instead, which +provides us with an expensive and cheap approach, respectively. diff --git a/examples/step-55/doc/results.dox b/examples/step-55/doc/results.dox index af80c1a67a..cbe2bc2e17 100644 --- a/examples/step-55/doc/results.dox +++ b/examples/step-55/doc/results.dox @@ -150,9 +150,13 @@ As can be seen from the table, -1. UMFPACK uses large amounts of memory, especially in 3d. Also, UMFPACK timings do not scale with problem size. +1. UMFPACK uses large amounts of memory, especially in 3d. Also, +UMFPACK timings do not scale with problem size. -2. The number of iterations for $A$ increase for ILU with refinement leading to worse then linear scaling in solve time. In contrast, the number of inner iterations for $A$ stay constant with GMG leading to nearly perfect scaling in solve time. +2. The number of iterations for $A$ increase for ILU with refinement +leading to worse then linear scaling in solve time. In contrast, the +number of inner iterations for $A$ stay constant with GMG leading to +nearly perfect scaling in solve time. 3. GMG needs slightly more memory than ILU. diff --git a/examples/step-55/doc/tooltip b/examples/step-55/doc/tooltip index 2c54b81f63..d33c3abeef 100644 --- a/examples/step-55/doc/tooltip +++ b/examples/step-55/doc/tooltip @@ -1 +1 @@ -Precondtioners: Geometric multigrid vs ILU +Geometric multigrid preconditioners for the Stokes problem diff --git a/examples/step-55/step-55.cc b/examples/step-55/step-55.cc index e2c021e839..f66d09562c 100644 --- a/examples/step-55/step-55.cc +++ b/examples/step-55/step-55.cc @@ -61,7 +61,7 @@ #include #include -// We need to include this class for all of the timings between ILU and Multigrid +// We need to include the followign file for all of the timings between ILU and Multigrid #include // This includes the files necessary for us to use geometric Multigrid @@ -80,7 +80,7 @@ namespace Step55 using namespace dealii; // In order to make it easy to switch between the different solvers that are being - // used in Step-55, an enum was created that can be passed as an argument to the + // used in Step-55, we declare an enum that can be passed as an argument to the // constructor of the main class. struct SolverType { @@ -95,13 +95,16 @@ namespace Step55 // we decided to separate the implementations for 2d and 3d using // template specialization. We do this to make it easier for us to debug // as well as its aesthetic value. + // + // Please note that the first dim components are the velocity components + // and the last is the pressure. template class Solution : public Function { public: Solution () : Function(dim+1) {} virtual double value (const Point &p, - const unsigned int component) const; + const unsigned int component = 0) const; virtual Tensor<1,dim> gradient (const Point &p, const unsigned int component = 0) const; }; @@ -111,9 +114,11 @@ namespace Step55 Solution<2>::value (const Point<2> &p, const unsigned int component) const { + Assert (component <= dim, ExcIndexRange(component,0,dim+1)); + using numbers::PI; - double x = p(0); - double y = p(1); + const double x = p(0); + const double y = p(1); if (component == 0) return sin (PI * x); @@ -122,7 +127,6 @@ namespace Step55 if (component == 2) return sin (PI * x) * cos (PI * y); - Assert (false, ExcMessage ("Component out of range in Solution")); return 0; } @@ -131,10 +135,12 @@ namespace Step55 Solution<3>::value (const Point<3> &p, const unsigned int component) const { + Assert (component <= dim, ExcIndexRange(component,0,dim+1)); + using numbers::PI; - double x = p(0); - double y = p(1); - double z = p(2); + const double x = p(0); + const double y = p(1); + const double z = p(2); if (component == 0) return 2.0 * sin (PI * x); @@ -145,7 +151,6 @@ namespace Step55 if (component == 3) return sin (PI * x) * cos (PI * y) * sin (PI * z); - Assert (false, ExcMessage ("Component out of range in Solution")); return 0; } @@ -155,6 +160,8 @@ namespace Step55 Solution<2>::gradient (const Point<2> &p, const unsigned int component) const { + Assert (component <= dim, ExcIndexRange(component,0,dim+1)); + using numbers::PI; double x = p(0); double y = p(1); @@ -174,8 +181,7 @@ namespace Step55 return_value[0] = PI * cos (PI * x) * cos (PI * y); return_value[1] = - PI * sin (PI * x) * sin(PI * y); } - else - Assert (false, ExcMessage ("Component out of range in Solution")); + return return_value; } @@ -184,6 +190,8 @@ namespace Step55 Solution<3>::gradient (const Point<3> &p, const unsigned int component) const { + Assert (component <= dim, ExcIndexRange(component,0,dim+1)); + using numbers::PI; double x = p(0); double y = p(1); @@ -213,8 +221,7 @@ namespace Step55 return_value[1] = - PI * sin (PI * x) * sin(PI * y) * sin (PI * z); return_value[2] = PI * sin (PI * x) * cos (PI * y) * cos (PI * z); } - else - Assert (false, ExcMessage ("Component out of range in Solution")); + return return_value; } @@ -234,6 +241,8 @@ namespace Step55 RightHandSide<2>::value (const Point<2> &p, const unsigned int component) const { + Assert (component <= dim, ExcIndexRange(component,0,dim+1)); + using numbers::PI; double x = p(0); double y = p(1); @@ -244,7 +253,6 @@ namespace Step55 if (component == 2) return 0; - Assert (false, ExcMessage ("Component out of range in RightHandSide")); return 0; } @@ -254,6 +262,8 @@ namespace Step55 RightHandSide<3>::value (const Point<3> &p, const unsigned int component) const { + Assert (component <= dim, ExcIndexRange(component,0,dim+1)); + using numbers::PI; double x = p(0); double y = p(1); @@ -267,7 +277,6 @@ namespace Step55 if (component == 3) return 0; - Assert (false, ExcMessage ("Component out of range in RightHandSide")); return 0; } @@ -432,30 +441,30 @@ namespace Step55 void compute_errors (); void output_results (const unsigned int refinement_cycle) const; - const unsigned int degree; - SolverType::type solver_type; + const unsigned int degree; + SolverType::type solver_type; - Triangulation triangulation; - FESystem fe; - FESystem velocity_fe; - DoFHandler dof_handler; - DoFHandler velocity_dof_handler; + Triangulation triangulation; + FESystem velocity_fe; + FESystem fe; + DoFHandler dof_handler; + DoFHandler velocity_dof_handler; - ConstraintMatrix constraints; + ConstraintMatrix constraints; - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix system_matrix; - SparseMatrix pressure_mass_matrix; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + SparseMatrix pressure_mass_matrix; - BlockVector solution; - BlockVector system_rhs; + BlockVector solution; + BlockVector system_rhs; MGLevelObject mg_sparsity_patterns; - MGLevelObject > mg_matrices; - MGLevelObject > mg_interface_matrices; + MGLevelObject > mg_matrices; + MGLevelObject > mg_interface_matrices; MGConstrainedDoFs mg_constrained_dofs; - TimerOutput computing_timer; + TimerOutput computing_timer; }; @@ -466,9 +475,9 @@ namespace Step55 degree (degree), solver_type (solver_type), triangulation (Triangulation::maximum_smoothing), - fe (FE_Q(degree+1), dim, // Finite element for whole system - FE_Q(degree), 1), velocity_fe (FE_Q(degree+1), dim), // Finite element for velocity-only + fe (velocity_fe, 1, // Finite element for whole system + FE_Q (degree), 1), dof_handler (triangulation), velocity_dof_handler (triangulation), computing_timer (std::cout, TimerOutput::summary, @@ -479,7 +488,7 @@ namespace Step55 // @sect4{StokesProblem::setup_dofs} -// This function sets up things differently based on if you want to use ILU or GMG as a preconditioner. +// This function sets up things based on if you want to use ILU or GMG as a preconditioner. template void StokesProblem::setup_dofs () { @@ -491,7 +500,10 @@ namespace Step55 // We don't need the multigrid dofs for whole problem finite element dof_handler.distribute_dofs(fe); - // This first creates and array (0,0,1) which means that it first does everything with index 0 and then 1 + // In the following code block, we first create an array of length dim+1 + // that is initialized to all zeros; we then set the pressure vector + // component to 1. This allows us to keep our velocities together + // and separate from the pressure. std::vector block_component (dim+1,0); block_component[dim] = 1; @@ -514,7 +526,9 @@ namespace Step55 // Distribute only the dofs for velocity finite element velocity_dof_handler.distribute_dofs(velocity_fe); - // Multigrid only needs the dofs for velocity + // Multigrid only needs the dofs for velocity. This does not clear the + // mg_interface_matrices object. Instead, it actually clears the level + // objects it stores. velocity_dof_handler.distribute_mg_dofs(velocity_fe); typename FunctionMap::type boundary_condition_function_map; @@ -613,7 +627,9 @@ namespace Step55 system_matrix=0; system_rhs=0; - double mass_factor = (solver_type == SolverType::UMFPACK) ? 0.0 : 1.0; + // The following bollean is used to signify when you want to assemble the mass matrix + // inside the (1,1) block, which is the case when you are not using UMFPACK + bool assemble_pressure_mass_matrix = (solver_type == SolverType::UMFPACK) ? false : true; QGauss quadrature_formula(degree+2); @@ -671,7 +687,7 @@ namespace Step55 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] - + mass_factor * phi_p[i] * phi_p[j]) + + (assemble_pressure_mass_matrix ? phi_p[i] * phi_p[j] : 0)) * fe_values.JxW(q); } @@ -819,9 +835,11 @@ namespace Step55 // @sect4{StokesProblem::solve} -// This function sets up things differently based on if you want to use ILU or GMG as a preconditioner. Both methods share -// the same solver (GMRES) but require a different preconditioner to be assembled. Here we time not only the entire solve -// function, but we separately time the set-up of the preconditioner as well as the GMRES solve. +// This function sets up things differently based on if you want to use ILU +// or GMG as a preconditioner. Both methods share the same solver (GMRES) +// but require a different preconditioner to be assembled. Here we time not +// only the entire solve function, but we separately time the set-up of the +// preconditioner as well as the GMRES solve. template void StokesProblem::solve () { @@ -870,14 +888,14 @@ namespace Step55 SparseILU A_preconditioner; A_preconditioner.initialize (system_matrix.block(0,0)); - SparseILU pmass_preconditioner; - pmass_preconditioner.initialize (pressure_mass_matrix); + SparseILU S_preconditioner; + S_preconditioner.initialize (pressure_mass_matrix); const BlockSchurPreconditioner, SparseILU > preconditioner (system_matrix, pressure_mass_matrix, A_preconditioner, - pmass_preconditioner, + S_preconditioner, use_expensive); computing_timer.leave_subsection(); @@ -936,8 +954,8 @@ namespace Step55 PreconditionMG, MGTransferPrebuilt > > A_Multigrid(velocity_dof_handler, mg, mg_transfer); - SparseILU pmass_preconditioner; - pmass_preconditioner.initialize (pressure_mass_matrix, + SparseILU S_preconditioner; + S_preconditioner.initialize (pressure_mass_matrix, SparseILU::AdditionalData()); const BlockSchurPreconditioner< @@ -946,7 +964,7 @@ namespace Step55 preconditioner (system_matrix, pressure_mass_matrix, A_Multigrid, - pmass_preconditioner, + S_preconditioner, use_expensive); computing_timer.leave_subsection(); -- 2.39.5