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})
-<a name="Intro"></a>
+<a name="Intro"></a>
<h1>Introduction</h1>
<h3> Stokes Problem </h3>
-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).
-
-<h3> Linear solver and preconditioning issues </h3>
-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).
+
+<h3> Linear solver and preconditioning issues </h3> 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.
<h3> Reference Solution </h3>
-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 <code>Function@<dim@></code>, 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 <code>Function@<dim@></code>, 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.
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}
-
-<h3> Computing Errors </h3>
-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 <code>compute_mean_value</code> function to compute the mean value of the pressure to subtract it from the pressure.
-
-<h3> DoF Handlers </h3>
-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.
-
-<h3> Differences from Step-22 </h3>
-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}
+
+<h3> Computing Errors </h3>
+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 <code>compute_mean_value</code> function to compute the mean value
+of the pressure to subtract it from the pressure.
+
+<h3> DoF Handlers </h3>
+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.
+
+<h3> Differences from Step-22 </h3>
+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.
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.
-Precondtioners: Geometric multigrid vs ILU
+Geometric multigrid preconditioners for the Stokes problem
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/grid/grid_out.h>
-// 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 <deal.II/base/timer.h>
// This includes the files necessary for us to use geometric Multigrid
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
{
// 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 <int dim>
class Solution : public Function<dim>
{
public:
Solution () : Function<dim>(dim+1) {}
virtual double value (const Point<dim> &p,
- const unsigned int component) const;
+ const unsigned int component = 0) const;
virtual Tensor<1,dim> gradient (const Point<dim> &p,
const unsigned int component = 0) const;
};
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);
if (component == 2)
return sin (PI * x) * cos (PI * y);
- Assert (false, ExcMessage ("Component out of range in Solution"));
return 0;
}
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);
if (component == 3)
return sin (PI * x) * cos (PI * y) * sin (PI * z);
- Assert (false, ExcMessage ("Component out of range in Solution"));
return 0;
}
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);
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;
}
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);
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;
}
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);
if (component == 2)
return 0;
- Assert (false, ExcMessage ("Component out of range in RightHandSide"));
return 0;
}
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);
if (component == 3)
return 0;
- Assert (false, ExcMessage ("Component out of range in RightHandSide"));
return 0;
}
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<dim> triangulation;
- FESystem<dim> fe;
- FESystem<dim> velocity_fe;
- DoFHandler<dim> dof_handler;
- DoFHandler<dim> velocity_dof_handler;
+ Triangulation<dim> triangulation;
+ FESystem<dim> velocity_fe;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+ DoFHandler<dim> velocity_dof_handler;
- ConstraintMatrix constraints;
+ ConstraintMatrix constraints;
- BlockSparsityPattern sparsity_pattern;
- BlockSparseMatrix<double> system_matrix;
- SparseMatrix<double> pressure_mass_matrix;
+ BlockSparsityPattern sparsity_pattern;
+ BlockSparseMatrix<double> system_matrix;
+ SparseMatrix<double> pressure_mass_matrix;
- BlockVector<double> solution;
- BlockVector<double> system_rhs;
+ BlockVector<double> solution;
+ BlockVector<double> system_rhs;
MGLevelObject<SparsityPattern> mg_sparsity_patterns;
- MGLevelObject<SparseMatrix<double> > mg_matrices;
- MGLevelObject<SparseMatrix<double> > mg_interface_matrices;
+ MGLevelObject<SparseMatrix<double> > mg_matrices;
+ MGLevelObject<SparseMatrix<double> > mg_interface_matrices;
MGConstrainedDoFs mg_constrained_dofs;
- TimerOutput computing_timer;
+ TimerOutput computing_timer;
};
degree (degree),
solver_type (solver_type),
triangulation (Triangulation<dim>::maximum_smoothing),
- fe (FE_Q<dim>(degree+1), dim, // Finite element for whole system
- FE_Q<dim>(degree), 1),
velocity_fe (FE_Q<dim>(degree+1), dim), // Finite element for velocity-only
+ fe (velocity_fe, 1, // Finite element for whole system
+ FE_Q<dim> (degree), 1),
dof_handler (triangulation),
velocity_dof_handler (triangulation),
computing_timer (std::cout, TimerOutput::summary,
// @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 <int dim>
void StokesProblem<dim>::setup_dofs ()
{
// 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<unsigned int> block_component (dim+1,0);
block_component[dim] = 1;
// 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<dim>::type boundary_condition_function_map;
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<dim> quadrature_formula(degree+2);
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);
}
// @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 <int dim>
void StokesProblem<dim>::solve ()
{
SparseILU<double> A_preconditioner;
A_preconditioner.initialize (system_matrix.block(0,0));
- SparseILU<double> pmass_preconditioner;
- pmass_preconditioner.initialize (pressure_mass_matrix);
+ SparseILU<double> S_preconditioner;
+ S_preconditioner.initialize (pressure_mass_matrix);
const BlockSchurPreconditioner<SparseILU<double>, SparseILU<double> >
preconditioner (system_matrix,
pressure_mass_matrix,
A_preconditioner,
- pmass_preconditioner,
+ S_preconditioner,
use_expensive);
computing_timer.leave_subsection();
PreconditionMG<dim, Vector<double>, MGTransferPrebuilt<Vector<double> > >
A_Multigrid(velocity_dof_handler, mg, mg_transfer);
- SparseILU<double> pmass_preconditioner;
- pmass_preconditioner.initialize (pressure_mass_matrix,
+ SparseILU<double> S_preconditioner;
+ S_preconditioner.initialize (pressure_mass_matrix,
SparseILU<double>::AdditionalData());
const BlockSchurPreconditioner<
preconditioner (system_matrix,
pressure_mass_matrix,
A_Multigrid,
- pmass_preconditioner,
+ S_preconditioner,
use_expensive);
computing_timer.leave_subsection();
projects / dealii.git / commitdiff