This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision | ||
problems:parstokes [2013/01/08 01:36] gag [Approximate solution of systems with A and M] |
problems:parstokes [2019/04/13 13:12] (current) tzwinger [Introduction] |
||
---|---|---|---|
Line 1: | Line 1: | ||
====== ParStokes: Iterative Solver for Large-Scale Glaciological Applications ====== | ====== ParStokes: Iterative Solver for Large-Scale Glaciological Applications ====== | ||
+ | |||
+ | ==== Alternatives ==== | ||
+ | Since March 2019, a new [[problems: | ||
==== Introduction ==== | ==== Introduction ==== | ||
- | The linear | + | The linear |
- | | A G | |v| |f| | + | | A B^T | |v| |f| |
- | |-D C | |p| = |g| (1) | + | | B |
where the v-component of the solution vector represents the velocities, | 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, | + | and the p-component the pressure. The matrix |
- | C a stabilization term, G the gradient and D=-G^T | + | '' |
While direct solvers are a robust way of solving (1), their computational cost in 3D | 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 | 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 | for large 3D applications. Iterative methods are therefore needed. A new module called | ||
- | `ParStokes' | + | **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 [[http:// | ||
Line 23: | Line 27: | ||
method (called the 'outer iteration' | method (called the 'outer iteration' | ||
followed in the ParStokes module is to use GCR (the generalized conjugate residual method) | followed in the ParStokes module is to use GCR (the generalized conjugate residual method) | ||
- | as outer iteration, and a preconditioner of the block diagonal | + | as outer iteration, and a preconditioner of the block upper triangular |
- | | + | | A B^T | |
- | | 0 M |. (2) | + | P = | 0 M |. (2) |
Here M denotes the mass matrix, scaled by the element-wise fluidity of the current | Here M denotes the mass matrix, scaled by the element-wise fluidity of the current | ||
Line 45: | Line 49: | ||
A sif-file for the ParStokes module should have (at least) three Solver sections, one | 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 | + | for the coupled problem, one for the velocity system with '' |
- | system with M. An example could look like this: | + | system with '' |
Solver 1 | Solver 1 | ||
Equation = " | Equation = " | ||
- | Procedure = "DummyRoutine" "DummyRoutine" | + | |
- | Variable = -dofs 3 " | + | |
- | Variable Output = False | + | Variable = -dofs 3 " |
- | Exec Solver = " | + | |
- | + | ||
- | Element = "p:1 b:4" | + | |
- | Bubbles in Global System = False | + | |
- | + | ||
- | Linear System Symmetric = True | + | |
- | Linear System Scaling = True | + | |
- | Linear System Row Equilibration = Logical True | + | |
- | Linear System Solver = Iterative | + | |
- | Linear System Max Iterations = 250 | + | |
- | Linear System Preconditioning = ILU0 | + | |
- | Linear System Convergence Tolerance = 1.0e-6 | + | |
- | Linear System Abort Not Converged = False | + | |
- | Skip Compute Nonlinear Change = Logical True | + | |
- | Back Rotate N-T Solution = Logical False | + | |
- | Linear System Timing = True | + | |
End | End | ||
Solver 2 | Solver 2 | ||
Equation = " | Equation = " | ||
- | Procedure = "DummyRoutine" "DummyRoutine" | + | |
- | Variable = -dofs 1 " | + | |
- | Variable Output = False | + | Variable = -dofs 1 " |
- | Exec Solver = " | + | |
- | + | ||
- | Element = "p:1 b:4" | + | |
- | Bubbles in Global System = False | + | |
- | + | ||
- | 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 | + | |
- | Skip Compute Nonlinear Change = Logical True | + | |
- | Back Rotate N-T Solution = Logical False | + | |
- | Linear System Timing = True | + | |
End | End | ||
Line 98: | Line 70: | ||
Procedure = " | Procedure = " | ||
Element = "p:1 b: | Element = "p:1 b: | ||
- | Bubbles in Global System = False | ||
Variable = FlowVar | Variable = FlowVar | ||
Variable Dofs = 4 | Variable Dofs = 4 | ||
Line 108: | Line 79: | ||
!----------------------------------------------- | !----------------------------------------------- | ||
Block Preconditioning = Logical True | Block Preconditioning = Logical True | ||
- | Linear System Scaling = Logical True | ||
- | Linear System Row Equilibration = Logical True | ||
- | Linear System Solver = " | ||
Linear System GCR Restart = Integer 200 | Linear System GCR Restart = Integer 200 | ||
Linear System Max Iterations = 200 | Linear System Max Iterations = 200 | ||
Line 118: | Line 86: | ||
Nonlinear System Newton After Tolerance = 1.0e-3 | Nonlinear System Newton After Tolerance = 1.0e-3 | ||
End | End | ||
+ | | ||
+ | == Boundary Conditions == | ||
- | == Remarks == | + | 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: | ||
- | a) The entry ' | + | ! bedrock: |
- | so that the ParStokes module finds them. The function DummyRoutine is an empty subroutine | + | Boundary Condition 2 |
- | that you have to supply, it may look like this: | + | Name = " |
- | + | Target Boundaries(1) = 5 | |
- | < | + | |
- | | + | |
+ | | ||
+ | V 1 = Real 0.0 | ||
+ | V 2 = Real 0.0 | ||
+ | V 3 = Real 0.0 | ||
+ | | ||
| | ||
- | USE DefUtils | + | In the case of sliding conditions, the keyword "Slip Coefficient" |
- | USE SolverUtils | + | Surface Traction {1,2,3} = Real |
- | USE ElementUtils | + | (the same as the keyword '' |
- | IMPLICIT NONE | + | |
- | TYPE(Solver_t) :: Solver | + | (same as '' |
- | TYPE(Model_t) :: Model | + | |
- | REAL(KIND=dp) :: dt | + | |
- | LOGICAL :: TransientSimulation | + | |
- | + | ||
- | END SUBROUTINE DummyRoutine | + | |
- | </ | + | |
- | and is compiled as usual, by | + | == Remarks == |
- | + | a) All solver routines presented above are included in the basic Elmer(/Ice) distribution | |
- | elmerf90 -o DummyRoutine.so DummyRoutine.f90 | + | |
b) We propose to use linear elements for all variables here, augmented | b) We propose to use linear elements for all variables here, augmented | ||
- | with fourth order bubble functions to achieve stability. While the high | + | with four bubble functions |
order stabilization makes the assembly somewhat costly, it is important for | order stabilization makes the assembly somewhat costly, it is important for | ||
- | achieving a stable solution. The corresponding | + | achieving a stable solution. The corresponding |
Element = "p:1 b: | Element = "p:1 b: | ||
- | Bubbles in Global System = False | + | |
+ | The setting (see item f of this list) | ||
+ | | ||
+ | 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 [[http:// | ||
- | c) the scaling of the linear systems is essential in order to get | + | c) The scaling of the linear systems is essential in order to get |
appropriate stopping criteria for the nested iterations. For the | appropriate stopping criteria for the nested iterations. For the | ||
outer iteration (Solver 3 above), one should generally use | outer iteration (Solver 3 above), one should generally use | ||
Line 159: | Line 135: | ||
Linear System Row Equilibration = Logical True | Linear System Row Equilibration = Logical True | ||
- | For the systems with A and M it is possible to preserve symmetry by setting | + | d) For the systems with '' |
- | Linear System Scaling = Logical True | + | Linear System Scaling = Logical True |
- | Linear System Row Equilibration = Logical False | + | Linear System Row Equilibration = Logical False |
In that case CG can be used instead of BiCGStab. | In that case CG can be used instead of BiCGStab. | ||
Line 175: | Line 151: | ||
symmetry causes the problems. | symmetry causes the problems. | ||
- | d) It is possible to replace A in the preconditioner by the block diagonal | + | e) It is possible to replace A in the preconditioner by the block diagonal |
matrix | matrix | ||
Line 195: | Line 171: | ||
outer iteration decreases. | outer iteration decreases. | ||
- | e) The variables in the solver sections for the preconditioner | + | f) ParStokes by default does only solve the Stokes |
- | identified as being the same as in the coupled system (Solver 3), 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 | + | |
- | ! bedrock: | + | Bubbles in Global System |
- | + | ||
- | Boundary Condition 2 | + | to avoid the by default implemented elimination of bubbles from the global system. In that case one also should take care that the '' |
- | Name = " | + | |
- | | + | |
- | | + | |
- | 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 | + | |
+ | g) It should be possible to restart a ParStokes simulation | ||
==== Approximate solution of systems with A and M ==== | ==== Approximate solution of systems with A and M ==== | ||
+ | |||
+ | 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 | The solution of systems with M is relatively easy in practice, we recommend using a CG | ||
Line 226: | Line 193: | ||
Variable Output = False | Variable Output = False | ||
Exec Solver = " | Exec Solver = " | ||
- | Element = "p:1 b:4" | ||
- | Bubbles in Global System = False | ||
Linear System Symmetric = True | Linear System Symmetric = True | ||
Linear System Scaling = True | Linear System Scaling = True | ||
Line 236: | Line 201: | ||
Linear System Preconditioning = Diagonal | Linear System Preconditioning = Diagonal | ||
Linear System Residual Output = 10 | Linear System Residual Output = 10 | ||
- | Skip Compute Nonlinear Change = Logical True | ||
- | Back Rotate N-T Solution = Logical False | ||
Linear System Timing = True | Linear System Timing = True | ||
End | End | ||
Line 257: | Line 220: | ||
this variant, set in the according solver section | this variant, set in the according solver section | ||
- | Linear System Use HYPRE = Logical True | + | Linear System Use HYPRE = Logical True |
The parameters | The parameters |