Stokes Solver¶
We want to solve the following Stokes block system.
If we apply Gaussian elimination to the above as a 2x2 block matrix system we can write this as:
where \(S=G^{T}K^{-1}G-C\) is the Schur complement and \(\hat{h}=G^{T}K^{-1}f -h\).
This system is now solved first for the pressure using the Schur complement matrix, \(S\). Then a backsolve for the velocity gives the complete solution.
Note that wherever \(K^{-1}\) appears, the inverse is never explicitly calculated but is achieved via a PETSc solve method. While solving for the pressure, there are necessarily solves using \(K\) inside of the matrix \(S\). We often refer to these as ‘inner’ solves.
Basic usage of the Stokes solver class involves being able to easily set up the inner solves in a few different ways (Setting up the pressure solve is more advanced).
To illustrate some basic usage let’s set up a simple problem to solve.
import underworld as uw from underworld import function as fn res=128 mesh = uw.mesh.FeMesh_Cartesian(elementRes=(res,res)) velocityField = mesh.add_variable(2) pressureField = mesh.subMesh.add_variable(1) velocityField.data[:] = (0.,0.) pressureField.data[:] = 0. # We are going to make use of one of the existing analytic solutions so that we may easily # obtain functions for a viscosity profile and forcing terms. # Exact solution solCx with defaults sol = fn.analytic.SolCx() stokesSystem = uw.systems.Stokes(velocityField,pressureField,sol.fn_viscosity,sol.fn_bodyforce,conditions=sol.get_bcs(velocityField)) # Now we create a solver. solver=uw.systems.Solver(stokesSystem) # The Stokes solver will use multigrid as a preconditioner along with # PETSc's 'fgmres' Krylov method by default for the 'inner' solve. solver.solve() # Now let's set up the inner solve to do a direct solve. # Note that the `lu` direct solver will not work in parallel. solver.set_inner_method("lu") solver.solve()
Let’s run underworld’s help function on the solver.configure function.
help(solver.set_inner_method)
Help on method set_inner_method in module underworld.systems._bsscr: set_inner_method(self, solve_type='mg') method of underworld.systems._bsscr.StokesSolver instance Configure velocity/inner solver (A11 PETSc prefix). Available options below. Note that associated solver software (for example mumps) must be installed. - mg : Geometric multigrid (default). - nomg : Disables multigrid. - lu : LU direct solver (serial only). - mumps : MUMPS parallel direct solver. - superludist : SuperLU parallel direct solver. - superlu : SuperLU direct solver (serial only).
We can see all the of the options that are configured in the solver
using the list()
functions for each component of the solver. These
are the most useful ones.
print("System Level Options:") solver.options.main.list() print("") print("Schur Complement Solve Options:") solver.options.scr.list() # Specifics for the print("") print("Inner (velocity) Solve Options:") solver.options.A11.list() # Specifics for the inner (velocity) solve print("") print("Multigrid (where enabled) Options:") solver.options.mg.list() # Options relevant to multigrid (if chosen)
System Level Options:
('remove_constant_pressure_null_space', False)
('ksp_k2_type', 'NULL')
('change_backsolve', False)
('penalty', 0.0)
('pc_type', 'none')
('force_correction', True)
('k_scale_only', True)
('Q22_pc_type', 'uw')
('change_A11rhspresolve', False)
('ksp_type', 'bsscr')
('rescale_equations', False)
('restore_K', False)
Schur Complement Solve Options:
('ksp_type', 'fgmres')
('ksp_rtol', 1e-05)
Inner (velocity) Solve Options:
('pc_type', 'lu')
('_mg_active', False)
('ksp_type', 'preonly')
Multigrid (where enabled) Options:
('active', False)
('levels', 8)
Further information about options is available via the help()
Python
function:
help(solver.options.A11)
The options generally follow PETSc
naming conventions.
A useful trick is to be able to imitate the classic “penalty method” which can be very efficient with modest-sized (2D) problems.
In the penalty method, we add a term to the weak form of the Stokes equation which penalises \(\lambda | \nabla \cdot \mathbf{u}| > 0\) and where \(\lambda\) is a sufficiently large constant to enforce the constraint. Typically \(10^7\) is considered sufficient.
The problem with this method is that the condition number of the system is severely compromised by adding the penalty term and standard iterative methods do not work well.
Our solvers have been configured with the penalty term and, for sufficiently robust choices of the inner solver, this can help solve problems faster (by improving pressure convergence).
An indestructible solver like lu
or mumps
(Mumps is a direct
solve that will work in parallel) can use very large penalties. Hence we
can recreate the penalty method as follows (though it still follows the
pattern of the Schur complement solver, while the classical method takes
some shortcuts).
try: solver.set_inner_method("mumps") solver.options.scr.ksp_type="cg" solver.set_penalty(1.0e7) solver.solve() solver.print_stats() except RuntimeError: # If the above fails, "mumps" probably isn't installed pass
[1;35m
Pressure iterations: 4
Velocity iterations: 1 (presolve)
Velocity iterations: 4 (pressure solve)
Velocity iterations: 1 (backsolve)
Velocity iterations: 6 (total solve)
SCR RHS setup time: 4.8789e-01
SCR RHS solve time: 5.6970e-03
Pressure setup time: 1.0583e-03
Pressure solve time: 2.6979e-02
Velocity setup time: 8.5831e-06 (backsolve)
Velocity solve time: 5.4877e-03 (backsolve)
Total solve time : 5.7443e-01
Velocity solution min/max: 0.0000e+00/0.0000e+00
Pressure solution min/max: 0.0000e+00/0.0000e+00
[00m
Now let’s go back to using multigrid. We can use a penalty here too, but the gigantic numbers won’t work.
solver.set_inner_method("mg") solver.set_penalty(1.0) solver.solve() solver.print_stats()
[1;35m
Pressure iterations: 17
Velocity iterations: 6 (presolve)
Velocity iterations: 88 (pressure solve)
Velocity iterations: 5 (backsolve)
Velocity iterations: 99 (total solve)
SCR RHS setup time: 2.8383e-02
SCR RHS solve time: 8.8551e-02
Pressure setup time: 1.5616e-03
Pressure solve time: 8.6988e-01
Velocity setup time: 9.5367e-06 (backsolve)
Velocity solve time: 3.8740e-02 (backsolve)
Total solve time : 1.1070e+00
Velocity solution min/max: 0.0000e+00/0.0000e+00
Pressure solution min/max: 0.0000e+00/0.0000e+00
[00m