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
+<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*}
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}
+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
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} \widetilde{A} & B^T \\ 0 &
+ \widetilde{S} \end{array}\right)^{-1} @f}
-is a good choice. It is important to note that
+is a good choice, where $\widetilde{A}$ is an approximation of $A$, $\widetilde{S}$ is an approximation of $S$, and
@f{eqnarray*}
P =
\left(\begin{array}{cc} A^{-1} & 0 \\ 0 & I \end{array}\right)
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
+matrix $A$. Therefore, in the above equations, $-M_p=\widetilde{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
++ ILU, and $\widetilde{A^{-1}}$ is an approximation of $A^{-1}$ obtained by one of
+multiple methods: CG with ILU as preconditioner, just using ILU, CG with GMG (Geometric
+Multigrid as described in step-16) as a precondtioner, 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.
+solver, it is important to set the first pressure node equal to zero
+to avoid the system being singular (recall that the Stokes equation
+itself only determines the pressure up to a constant when using only
+Dirichlet boundary conditions). If we do not do this, then the direct
+solver will produce an error message whereas the iterative solvers
+quietly solve it anyway.
<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
+The domain, right hand side, and boundary conditions we implement
+are chosen for their simplicity and the fact that they make 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
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
+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
+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
+the Method of Manufactured Solutions of step-7, we need to find $\bf
f$ such that:
@f{align*}
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
+the VectorTools::compute_mean_value() function to compute the mean value
of the pressure to subtract it from the pressure.
-<h3> DoF Handlers </h3>
+<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
+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>
+<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
+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 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
projects / dealii.git / commitdiff