Stokes Solver

We want to solve the following Stokes block system.

\[\begin{split}\begin{bmatrix} K & G \\ G^T & C \end{bmatrix} \begin{bmatrix} u\\ p \end{bmatrix} = \begin{bmatrix} f\\ h \end{bmatrix}.\end{split}\]

If we apply Gaussian elimination to the above as a 2x2 block matrix system we can write this as:

\[\begin{split}\begin{bmatrix} K & G\\ 0 & S \end{bmatrix} \begin{bmatrix} u\\ p \end{bmatrix} = \begin{bmatrix} f\\ \hat{h} \end{bmatrix},\end{split}\]

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).

solver.set_inner_method("mumps")
solver.options.scr.ksp_type="cg"
solver.set_penalty(1.0e7)
solver.solve()
solver.print_stats()


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



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()


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