This tutorial is part of multiphenicsx.

Tutorial 04: computation of the inf-sup constant for a Stokes problem discretization

In this tutorial we compare the computation of the inf-sup constant of a Stokes problem using standard assembly with mixed function spaces and block assembly.

import typing

import basix.ufl
import dolfinx.cpp
import dolfinx.fem
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import slepc4py.SLEPc
import ufl
import viskex

import multiphenicsx.fem
import multiphenicsx.fem.petsc

Mesh

mesh = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 32, 32)

# Create connectivities required by the rest of the code
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

def wall(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
    """Determine the position of the wall."""
    return np.logical_or(  # type: ignore[no-any-return]
        x[1] < 0 + np.finfo(float).eps, x[1] > 1 - np.finfo(float).eps)


boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, wall)

viskex.dolfinx.plot_mesh(mesh)

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Function spaces

V_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim, ))
Q_element = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)

Auxiliary function for eigenvector normalization

def normalize(u1: dolfinx.fem.Function, u2: dolfinx.fem.Function, p: dolfinx.fem.Function) -> None:
    """Normalize an eigenvector."""
    scaling_operations: list[tuple[  # type: ignore[no-any-unimported]
        dolfinx.fem.Function, typing.Callable[[dolfinx.fem.Function], ufl.Form],
        typing.Callable[[petsc4py.PETSc.ScalarType], petsc4py.PETSc.ScalarType]
    ]] = [
        # Scale functions with a W^{1,1} (for velocity) or L^1 (for pressure) norm to take away
        # possible sign differences.
        (u1, lambda u: (u.dx(0) + u.dx(1)) * ufl.dx, lambda x: x),
        (u2, lambda u: (u.dx(0) + u.dx(1)) * ufl.dx, lambda x: x),
        (p, lambda p: p * ufl.dx, lambda x: x),
        # Normalize functions with a H^1 (for velocity) or L^2 (for pressure) norm.
        (u1, lambda u: ufl.inner(ufl.grad(u), ufl.grad(u)) * ufl.dx, lambda x: np.sqrt(x)),
        (u2, lambda u: ufl.inner(ufl.grad(u), ufl.grad(u)) * ufl.dx, lambda x: np.sqrt(x)),
        (p, lambda p: ufl.inner(p, p) * ufl.dx, lambda x: np.sqrt(x))
    ]
    for (function, bilinear_form, postprocess) in scaling_operations:
        scalar = postprocess(mesh.comm.allreduce(
            dolfinx.fem.assemble_scalar(dolfinx.fem.form(bilinear_form(function))), op=mpi4py.MPI.SUM))
        function.x.petsc_vec.scale(1. / scalar)
        function.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)

Standard formulation using a mixed function space

def run_monolithic() -> tuple[np.float64, dolfinx.fem.Function, dolfinx.fem.Function, dolfinx.fem.Function]:
    """Run standard formulation using a mixed function space."""
    # Function spaces
    W_element = basix.ufl.mixed_element([V_element, Q_element])
    W = dolfinx.fem.functionspace(mesh, W_element)

    # Test and trial functions: monolithic
    vq = ufl.TestFunction(W)
    (v, q) = ufl.split(vq)
    up = ufl.TrialFunction(W)
    (u, p) = ufl.split(up)

    # Variational forms
    lhs = (ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - ufl.inner(p, ufl.div(v)) * ufl.dx
           - ufl.inner(ufl.div(u), q) * ufl.dx)
    rhs = - ufl.inner(p, q) * ufl.dx

    # Define restriction for DOFs associated to homogenous Dirichlet boundary conditions
    dofs_W = np.arange(0, W.dofmap.index_map.size_local + W.dofmap.index_map.num_ghosts)
    W0 = W.sub(0)
    V, _ = W0.collapse()
    bdofs_V = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundary_facets)[0]
    restriction = multiphenicsx.fem.DofMapRestriction(W.dofmap, np.setdiff1d(dofs_W, bdofs_V))

    # Assemble lhs and rhs matrices
    A = multiphenicsx.fem.petsc.assemble_matrix(
        dolfinx.fem.form(lhs), restriction=(restriction, restriction))
    A.assemble()
    B = multiphenicsx.fem.petsc.assemble_matrix(
        dolfinx.fem.form(rhs), restriction=(restriction, restriction))
    B.assemble()

    # Solve
    eps = slepc4py.SLEPc.EPS().create(mesh.comm)
    eps.setOperators(A, B)
    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GNHEP)
    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.TARGET_REAL)
    eps.setTarget(1.e-5)
    eps.getST().setType(slepc4py.SLEPc.ST.Type.SINVERT)
    eps.getST().getKSP().setType("preonly")
    eps.getST().getKSP().getPC().setType("lu")
    eps.getST().getKSP().getPC().setFactorSolverType("mumps")
    eps.solve()
    assert eps.getConverged() >= 1

    # Extract leading eigenvalue and eigenvector
    vr = dolfinx.cpp.fem.petsc.create_vector_block([(restriction.index_map, restriction.index_map_bs)])
    vi = dolfinx.cpp.fem.petsc.create_vector_block([(restriction.index_map, restriction.index_map_bs)])
    eigv = eps.getEigenpair(0, vr, vi)
    r, i = eigv.real, eigv.imag
    assert abs(i) < 1.e-10
    assert r > 0., "r = " + str(r) + " is not positive"
    print("Inf-sup constant (monolithic): ", np.sqrt(r))

    # Transform eigenvector into eigenfunction
    r_fun = dolfinx.fem.Function(W)
    vr.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
    with r_fun.x.petsc_vec.localForm() as r_fun_local, \
            multiphenicsx.fem.petsc.VecSubVectorWrapper(vr, W.dofmap, restriction) as vr_wrapper:
        r_fun_local[:] = vr_wrapper
    u_fun = r_fun.sub(0).collapse()
    (u_fun_1, u_fun_2) = (u_fun.sub(0).collapse(), u_fun.sub(1).collapse())
    p_fun = r_fun.sub(1).collapse()
    normalize(u_fun_1, u_fun_2, p_fun)

    eps.destroy()
    return (r, u_fun_1, u_fun_2, p_fun)

(eig_m, u_fun_1_m, u_fun_2_m, p_fun_m) = run_monolithic()

Inf-sup constant (monolithic):  0.6058505264323524

viskex.dolfinx.plot_scalar_field(u_fun_1_m, "u1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(u_fun_2_m, "u2")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(p_fun_m, "p")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Block formulation using two independent function spaces

def run_block() -> tuple[np.float64, dolfinx.fem.Function, dolfinx.fem.Function, dolfinx.fem.Function]:
    """Run block formulation using two independent function spaces."""
    # Function spaces
    V = dolfinx.fem.functionspace(mesh, V_element)
    Q = dolfinx.fem.functionspace(mesh, Q_element)

    # Test and trial functions
    (v, q) = (ufl.TestFunction(V), ufl.TestFunction(Q))
    (u, p) = (ufl.TrialFunction(V), ufl.TrialFunction(Q))

    # Variational forms
    lhs = [[ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, - ufl.inner(p, ufl.div(v)) * ufl.dx],
           [- ufl.inner(ufl.div(u), q) * ufl.dx, None]]
    rhs = [[None, None],
           [None, - ufl.inner(p, q) * ufl.dx]]
    rhs[0][0] = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0)) * ufl.inner(u, v) * ufl.dx

    # Define restriction for DOFs associated to homogenous Dirichlet boundary conditions
    dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)
    bdofs_V = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)
    dofs_Q = np.arange(0, Q.dofmap.index_map.size_local + Q.dofmap.index_map.num_ghosts)
    restriction_V = multiphenicsx.fem.DofMapRestriction(V.dofmap, np.setdiff1d(dofs_V, bdofs_V))
    restriction_Q = multiphenicsx.fem.DofMapRestriction(Q.dofmap, dofs_Q)
    restriction = [restriction_V, restriction_Q]

    # Assemble lhs and rhs matrices
    A = multiphenicsx.fem.petsc.assemble_matrix_block(
        dolfinx.fem.form(lhs), bcs=[], restriction=(restriction, restriction))
    A.assemble()
    B = multiphenicsx.fem.petsc.assemble_matrix_block(
        dolfinx.fem.form(rhs), bcs=[], restriction=(restriction, restriction))
    B.assemble()

    # Solve
    eps = slepc4py.SLEPc.EPS().create(mesh.comm)
    eps.setOperators(A, B)
    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GNHEP)
    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.TARGET_REAL)
    eps.setTarget(1.e-5)
    eps.getST().setType(slepc4py.SLEPc.ST.Type.SINVERT)
    eps.getST().getKSP().setType("preonly")
    eps.getST().getKSP().getPC().setType("lu")
    eps.getST().getKSP().getPC().setFactorSolverType("mumps")
    eps.solve()
    assert eps.getConverged() >= 1

    # Extract leading eigenvalue and eigenvector
    vr = dolfinx.cpp.fem.petsc.create_vector_block(
        [(restriction_.index_map, restriction_.index_map_bs) for restriction_ in restriction])
    vi = dolfinx.cpp.fem.petsc.create_vector_block(
        [(restriction_.index_map, restriction_.index_map_bs) for restriction_ in restriction])
    eigv = eps.getEigenpair(0, vr, vi)
    r, i = eigv.real, eigv.imag
    assert abs(i) < 1.e-10
    assert r > 0., "r = " + str(r) + " is not positive"
    print("Inf-sup constant (block): ", np.sqrt(r))

    # Transform eigenvector into eigenfunction
    (u_fun, p_fun) = (dolfinx.fem.Function(V), dolfinx.fem.Function(Q))
    vr.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
    with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(vr, [V.dofmap, Q.dofmap], restriction) as vr_wrapper:
        for vr_wrapper_local, component in zip(vr_wrapper, (u_fun, p_fun)):
            with component.x.petsc_vec.localForm() as component_local:
                component_local[:] = vr_wrapper_local
    (u_fun_1, u_fun_2) = (u_fun.sub(0).collapse(), u_fun.sub(1).collapse())
    normalize(u_fun_1, u_fun_2, p_fun)

    eps.destroy()
    return (r, u_fun_1, u_fun_2, p_fun)

(eig_b, u_fun_1_b, u_fun_2_b, p_fun_b) = run_block()

Inf-sup constant (block):  0.6058505264323903

viskex.dolfinx.plot_scalar_field(u_fun_1_b, "u1")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.
MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(u_fun_2_b, "u2")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

viskex.dolfinx.plot_scalar_field(p_fun_b, "p")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.

MESA: error: ZINK: failed to choose pdev
glx: failed to create drisw screen

Error computation between standard and block formulations

def run_error(
    eig_m: np.float64, eig_b: np.float64, u_fun_1_m: dolfinx.fem.Function, u_fun_1_b: dolfinx.fem.Function,
    u_fun_2_m: dolfinx.fem.Function, u_fun_2_b: dolfinx.fem.Function, p_fun_m: dolfinx.fem.Function,
    p_fun_b: dolfinx.fem.Function
) -> None:
    """Compute errors between the mixed and block cases."""
    err_inf_sup = np.abs(np.sqrt(eig_b) - np.sqrt(eig_m)) / np.sqrt(eig_m)
    print("Relative error for inf-sup constant equal to", err_inf_sup)
    assert np.isclose(err_inf_sup, 0., atol=1.e-8)
    eigenvector_operations: list[tuple[  # type: ignore[no-any-unimported]
        dolfinx.fem.Function, dolfinx.fem.Function, typing.Callable[[dolfinx.fem.Function], ufl.Form], str
    ]] = [
        (u_fun_1_b, u_fun_1_m, lambda u: ufl.inner(ufl.grad(u), ufl.grad(u)) * ufl.dx, "velocity 1"),
        (u_fun_2_b, u_fun_2_m, lambda u: ufl.inner(ufl.grad(u), ufl.grad(u)) * ufl.dx, "velocity 2"),
        (p_fun_b, p_fun_m, lambda p: ufl.inner(p, p) * ufl.dx, "pressure")
    ]
    for (fun_b, fun_m, squared_norm_form, component_name) in eigenvector_operations:
        err_fun = np.sqrt(mesh.comm.allreduce(
            dolfinx.fem.assemble_scalar(dolfinx.fem.form(squared_norm_form(fun_b - fun_m))), op=mpi4py.MPI.SUM))
        print("Relative error for ", component_name, "component of eigenvector equal to", err_fun)
        assert np.isclose(err_fun, 0., atol=1.e-6)

run_error(eig_m, eig_b, u_fun_1_m, u_fun_1_b, u_fun_2_m, u_fun_2_b, p_fun_m, p_fun_b)

Relative error for inf-sup constant equal to 6.267160922640535e-14

Relative error for  velocity 1 component of eigenvector equal to 4.166542203869248e-11
Relative error for  velocity 2 component of eigenvector equal to 9.307500898476484e-11

Relative error for  pressure component of eigenvector equal to 2.810011638283906e-11