Nullspace extraction using mumps

script.py

from mpi4py import MPI

try:
    from petsc4py import PETSc

    import dolfinx

    if not dolfinx.has_petsc:
        print("This demo requires DOLFINx to be compiled with PETSc enabled.")
        exit(0)
except ModuleNotFoundError:
    print("This demo requires petsc4py.")
    exit(0)

# +
import numpy as np

import ufl
from dolfinx import la
from dolfinx.fem import (
    Expression,
    Function,
    FunctionSpace,
    dirichletbc,
    form,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, GhostMode, create_box, locate_entities_boundary

dtype = PETSc.ScalarType

msh = create_box(
    MPI.COMM_WORLD,
    [np.array([0.0, 0.0, 0.0]), np.array([2.0, 1.0, 1.0])],
    [16, 16, 16],
    CellType.tetrahedron,
    ghost_mode=GhostMode.shared_facet,
)

# Create a centripetal source term $f = \rho \omega^2 [x_0, \, x_1]$:

ω, ρ = 300.0, 10.0
x = ufl.SpatialCoordinate(msh)
f = ufl.as_vector((ρ * ω**2 * x[0], ρ * ω**2 * x[1], 0.0))

# +
E = 1.0e9
ν = 0.3
μ = E / (2.0 * (1.0 + ν))
λ = E * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))


def σ(v):
    """Return an expression for the stress σ given a displacement field"""
    return 2.0 * μ * ufl.sym(ufl.grad(v)) + λ * ufl.tr(ufl.sym(ufl.grad(v))) * ufl.Identity(len(v))



V = functionspace(msh, ("Lagrange", 1, (msh.geometry.dim,)))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = form(ufl.inner(σ(u), ufl.grad(v)) * ufl.dx)
L = form(ufl.inner(f, v) * ufl.dx)

facets = locate_entities_boundary(
    msh, dim=2, marker=lambda x: np.isclose(x[0], 0.0) | np.isclose(x[1], 1.0)
)
bc = dirichletbc(
    np.zeros(3, dtype=dtype), locate_dofs_topological(V, entity_dim=2, entities=facets), V=V
)

bcs = []
A = assemble_matrix(a, bcs=bcs)
A.assemble()

b = assemble_vector(L)
apply_lifting(b, [a], bcs=bcs)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)

for bc in bcs:
    bc.set(b.array_w)

ns = build_nullspace(V)
opts = PETSc.Options()

opts["mat_mumps_icntl_24"] = 1 # Zero pivot detection
opts["mat_mumps_icntl_25"] = 3 # Which nullspace to extract

opts["mat_mumps_icntl_4"] = 2 # Verbosity
opts["mat_mumps_icntl_2"] = 1 # std out blaaah
opts["mat_mumps_cntl_3"] = 1e-6 # Threshold factor
opts["ksp_type"] = "preonly"

opts["pc_type"] = "lu"
opts["pc_factor_mat_solver_type"] = "mumps"
solver = PETSc.KSP().create(msh.comm)
solver.setFromOptions()

solver.setOperators(A)
uh = Function(V)

solver.setMonitor(lambda _, its, rnorm: print(f"Iteration: {its}, rel. residual: {rnorm}"))
print(uh.x.array)
solver.solve(b, uh.x.petsc_vec)
solver.view()
print(uh.x.array)

with dolfinx.io.XDMFFile(msh.comm, "out_elasticity/nullspace.xdmf", "w") as file:
    file.write_mesh(msh)
    file.write_function(uh)