Since March 2019, a new vectorized Stokes solver which further is capable of utilizing the library version of the block preconditioning is available. We recommend using this one.
The linear system that needs to be solved for ice dynamics in Elmer/Ice are of the form
| A B^T | |v| |f| | B C | |p| = |g| (1)
where the v-component of the solution vector represents the velocities,
and the p-component the pressure. The matrix A
is an elasticity-like operator,
C
a stabilization term, with B
the discrete version of the operator -div.
While direct solvers are a robust way of solving (1), their computational cost in 3D grows with the square of the number of unknowns, and their memory usage becomes excessive for large 3D applications. Iterative methods are therefore needed. A new module called ParStokes is therefore available in Elmer now, which is dedicated to this purpose. Detailed description of its functionality can be found in chapter 24 of the Elmer Models Manual.
The general strategy for the iterative solution of (1) is to use a Krylov-subspace method (called the 'outer iteration'), combined with a preconditioner P. The approach followed in the ParStokes module is to use GCR (the generalized conjugate residual method) as outer iteration, and a preconditioner of the block upper triangular form
| A B^T | P = | 0 M |. (2)
Here M denotes the mass matrix, scaled by the element-wise fluidity of the current approximation of the solution. It can be shown that, under certain regularity assumptions, this preconditioner leads to a convergence rate of the outer iteration which is independent of the mesh-size h. However, the application of the inverse of P still requires the exact solution of linear systems with A and M. These operations would still be too expensive with direct solvers for large 3D problems, so further approximations are made. Note that the inverse of P has to be applied in each outer iteration, so the approximate solution of systems with A and M should be very cheap compared to the exact solution of a system (1).
The ParStokes module is a dedicated solver for ice applications and is not as flexible as the NavierStokes module. Certain choices such as the outer iterative solver and the use of high-order bubble functions for stabilization are hard-coded.
A sif-file for the ParStokes module should have (at least) three Solver sections, one
for the coupled problem, one for the velocity system with A
, and one for the pressure
system with M
. An example could look like this:
Solver 1 Equation = "Velocity Preconditioning" Exec Solver = "Before Simulation" Procedure = "VelocityPrecond" "VelocityPrecond" Variable = -dofs 3 "V" ! that would be default End
Solver 2 Equation = "Pressure Preconditioning" Exec Solver = "Before Simulation" Procedure = "PressurePrecond" "PressurePrecond" Variable = -dofs 1 "p" ! that would be default End
Solver 3 Equation = "Stokes" Procedure = "ParStokes" "StokesSolver" Element = "p:1 b:4" Variable = FlowVar Variable Dofs = 4 Convective = Logical False Block Diagonal A = Logical True Use Velocity Laplacian = Logical True !----------------------------------------------- ! Keywords related to the block preconditioning !----------------------------------------------- Block Preconditioning = Logical True Linear System GCR Restart = Integer 200 Linear System Max Iterations = 200 Linear System Convergence Tolerance = 1.0e-6 Nonlinear System Max Iterations = 100 Nonlinear System Convergence Tolerance = 1.0e-5 Nonlinear System Newton After Tolerance = 1.0e-3 End
The variables in the solver sections for the preconditioner (Solver 1 and 2) are not identified as being the same as in the coupled system (Solver 3 in example above), so boundary conditions have to be assigned to all the components twice. For instance, to specify a no-slip condition at the bedrock the according section in the sif-file might look like this:
! bedrock: Boundary Condition 2 Name = "bedrock" Target Boundaries(1) = 5 FlowVar 1 = Real 0.0 FlowVar 2 = Real 0.0 FlowVar 3 = Real 0.0 V 1 = Real 0.0 V 2 = Real 0.0 V 3 = Real 0.0 End
In the case of sliding conditions, the keyword “Slip Coefficient” will be recognised by the preconditioner as well as the solver. Stress conditions are imposed using the keywords
Surface Traction {1,2,3} = Real
(the same as the keyword Pressure {1,2,3}
in Navier-Stokes) for component-wise setting, or
Normal Surface Traction = Real
(same as External Pressure
in Navier-Stokes) for imposing a normal force (e.g., sea pressure).
a) All solver routines presented above are included in the basic Elmer(/Ice) distribution
b) We propose to use linear elements for all variables here, augmented with four bubble functions per velocity component to achieve stability. While the high order stabilization makes the assembly somewhat costly, it is important for achieving a stable solution. The corresponding line in all three solver sections are
Element = "p:1 b:4"
The setting (see item f of this list)
Bubbles in Global System = False
is not needed, as it is the default in any of the solvers.
P2-P1/Q2-Q1 approximation offers an alternative for using p-bubbles (see section 24.4 of Elmer Models manual).
c) The scaling of the linear systems is essential in order to get appropriate stopping criteria for the nested iterations. For the outer iteration (Solver 3 above), one should generally use
Linear System Scaling = Logical True Linear System Row Equilibration = Logical True
d) For the systems with A
and M
it is possible to preserve symmetry by setting (both are default settings)
Linear System Scaling = Logical True Linear System Row Equilibration = Logical False
In that case CG can be used instead of BiCGStab.
note: while the matrices A and M are in principle symmetric, the symmetry may be destroyed by boundary conditions. The flag
Linear System Symmetric = True
will yield a symmetric system matrix for the most common boundary conditions, but in case CG fails to converge it may be worthwhile to try BiCGStab instead to see if a loss of symmetry causes the problems.
e) It is possible to replace A in the preconditioner by the block diagonal matrix
| A_{11} 0 0 | A_b= | 0 A_{22} 0 | | 0 0 A_{33} |
where couplings between the velocity components are neglected. This is selected by setting
Block Diagonal A = Logical True
To make the systems with A even cheaper, one can in addition set
Use Velocity Laplacian = Logical True
in which case the diagonal blocks $A_{ii}$ are replaced by scaled Laplace operators. The price you pay for these approximations is an increase in the number of outer iterations. The gain is that the systems are easier to solve, so the number of inner iterations in each outer iteration decreases.
f) ParStokes by default does only solve the Stokes (and not Navier-Stokes) equation. If one wishes to include acceleration terms, setting the keyword “Convective = Logical True” will do so. Also, 2 Picard iterations are always performed, independently from what is specified in the sif-file. The line “Reset Newton…” is not necessary. In this case, it would be good to include
Bubbles in Global System = True
to avoid the by default implemented elimination of bubbles from the global system. In that case one also should take care that the Element
keyword contains the same bubble degree for all three solvers.
g) It should be possible to restart a ParStokes simulation
In the example above, no numerical setting for the solution of the blocks M and A are given. The defualt settings (see source code) should suffice in most cases. Here, we just present some tweaks for these parts:
The solution of systems with M is relatively easy in practice, we recommend using a CG iterative solver with diagonal preconditioning or ILU0, and solving the system to relatively high accuracy:
Solver 2 Equation = "Pressure Preconditioning" Procedure = "DummyRoutine" "DummyRoutine" Variable = -dofs 1 "P" Variable Output = False Exec Solver = "before simulation" Linear System Symmetric = True Linear System Scaling = True Linear System Solver = iterative Linear System Iterative Method = CG Linear System Max Iterations = 1000 Linear System Convergence Tolerance = 1.0e-6 Linear System Preconditioning = Diagonal Linear System Residual Output = 10 Linear System Timing = True End
The systems with A are considerably harder to solve, but the accuracy to which they are solved is not so critical for the overall convergence of the outer iteration. The matrix A is related to a linear elasticity problem, and for such systems the standard ILU preconditioning in Elmer tends to lack robustness in parallel. This is related to the large null space (rigid body modes) of the matrices on each processor partition. For moderately sized problems (say, 16 partitions), ILU0 or ILU1 are likely to give good results. For large-scale applications we advocate using one of the parallel linear algebra packages interfaced by Elmer, HYPRE or Trilinos. Below we suggest two alternative solvers for parallel computations.
A preconditioner for A that does not suffer from the singularity of the partition-wise matrices is the parallel ILU solver Euclid, available through the HYPRE library. To use this variant, set in the according solver section
Linear System Use HYPRE = Logical True
The parameters
Linear System Iterative Method Linear System Preconditioning Linear System Convergence Tolerance Linear System Max Iterations
are used to determine the HYPRE strategy (BiCGStab and CG are implemented as solvers, and ILU0, ILU1, ILU2 etc should be used as preconditioner) This solver is likely to be slower for small numbers of partitions (np) but more robust wrt. np and more scalable (up to a few hundred processors).
In order to achieve good parallel performance for large-scale applications we propose a rather coarse approximation of the 'inverse of A' operation by a single multigrid cycle. Elmer should be compiled with support for Trilinos. Then set the parameters
Linear System Use Trilinos = Logical True Trilinos Parameter File = String input.xml
The parameter file input.xml contains settings for the smoothed aggregation solver ML. We are currently experimenting with this solver and will publish suggestions on their choice for achieving robustness later on. Some simple examples of using Trilinos with Elmer (with xml input files) can be found in the source code under
elmerfem/fem/examples/TrilinosSolvers