Hi everyone
I’ve run into some quite peculiar behaviour of the LinearSolver with LU factorisation and otherwise default settings.
I’m trying to solve the Stokes problem
-\Delta \boldsymbol{u} + \nabla p = \boldsymbol{0}\\
-\operatorname{div} \boldsymbol{u} = 0
with boundary conditions
\boldsymbol{u} = \boldsymbol{u}_D \qquad \text{on } \Gamma_D\\
-\frac{\partial \boldsymbol{u}}{\partial \boldsymbol{n}} + p \boldsymbol{n} = \boldsymbol{0} \qquad \text{on } \Gamma_N.
Depending on the mesh or the choice of elements, the LinearSolver often computes inf values, even though both the continuous and the discrete problems should be well-posed. If the code is run in parallel, or if "pc_factor_mat_solver_type" is set to "mumps" or "umfpack", the solution is correct. I cannot find any systematics behind this:
Unit square
"pc_factor_mat_solver_type" |
P2-P1 elements | P2-P0 elements |
|---|---|---|
| default |  |
|
"mumps" |
 |
|
"umfpack" |
|
|
DFG 2-D benchmark
"pc_factor_mat_solver_type" |
P2-P1 elements | P2-P0 elements |
|---|---|---|
| default | |
|
"mumps" |
|
|
"umfpack" |
|
|
MWE on unit square
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import TrialFunctions, TestFunctions, div, grad, inner, dx
domain = mesh.create_unit_square(
comm=MPI.COMM_WORLD,
nx=10,
ny=10,
cell_type=mesh.CellType.triangle,
diagonal=mesh.DiagonalType.crossed
)
# Function spaces
V_el = element("P", domain.topology.cell_name(), 2, shape=(2, ))
Q_el = element("DG", domain.topology.cell_name(), 0)
# Q_el = element("P", domain.topology.cell_name(), 1)
VQ = fem.functionspace(domain, mixed_element([V_el, Q_el]))
V, _ = VQ.sub(0).collapse()
# Dirichlet BCs for left, top and bottom walls
uD = fem.Function(V)
uD.interpolate(lambda x: (x[1]*(1. - x[1]), 0.*x[0]))
dirichlet_facets = mesh.locate_entities_boundary(
msh=domain,
dim=1,
marker=lambda x: np.isclose(x[0], 0.) | np.isclose(x[1], 1.) | np.isclose(x[1], 0.)
)
dirichlet_dofs = fem.locate_dofs_topological((VQ.sub(0), V), 1, dirichlet_facets)
bc = fem.dirichletbc(uD, dirichlet_dofs, VQ.sub(0))
# Weak form
(u, p) = TrialFunctions(VQ)
(v, q) = TestFunctions(VQ)
a = inner(grad(u), grad(v)) * dx - p * div(v) * dx - q * div(u) * dx
L = inner(fem.Constant(domain, (0., 0.)), v) * dx
lu_solver = {"ksp_type": "preonly", "pc_type": "lu"} # gives inf
# lu_solver["pc_factor_mat_solver_type"] = "mumps" # works
# lu_solver["pc_factor_mat_solver_type"] = "umfpack" # works
problem = LinearProblem(a, L, bcs=[bc], petsc_options=lu_solver)
uph = problem.solve()
uh, ph = uph.sub(0).collapse(), uph.sub(1).collapse()
u_error2_form = fem.form(inner(uh - uD, uh - uD) * dx + inner(grad(uh - uD), grad(uh - uD)) * dx)
u_error2_subdomain = fem.assemble_scalar(u_error2_form)
u_error2 = MPI.COMM_WORLD.reduce(u_error2_subdomain, op=MPI.SUM)
u_min = MPI.COMM_WORLD.reduce(np.min(uh.x.array), op=MPI.MIN)
u_max = MPI.COMM_WORLD.reduce(np.max(uh.x.array), op=MPI.MAX)
if MPI.COMM_WORLD.rank == 0:
print("Velocity error:", np.sqrt(u_error2))
print("Velocity min:", u_min)
print("Velocity max:", u_max)
MWE on DFG 2-D benchmark domain
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import gmsh
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import TrialFunctions, TestFunctions, div, grad, inner, dx
# Meshing
gmsh.initialize()
channel = gmsh.model.occ.addRectangle(x=0, y=0, z=0, dx=2.2, dy=.41)
obstacle = gmsh.model.occ.addDisk(xc=.2, yc=.2, zc=0, rx=.05, ry=.05)
gmsh.model.occ.cut([(2, channel)], [(2, obstacle)])
gmsh.model.occ.synchronize()
domain_entities = [entity[1] for entity in gmsh.model.getEntities(dim=2)]
gmsh.model.addPhysicalGroup(2, domain_entities, tag=1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", .05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", .05)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(dim=2)
domain, cell_tags, facet_tags = io.gmshio.model_to_mesh(
model=gmsh.model,
comm=MPI.COMM_WORLD,
rank=0,
gdim=2
)
gmsh.finalize()
# Function spaces
V_el = element("P", domain.topology.cell_name(), 2, shape=(2, ))
Q_el = element("DG", domain.topology.cell_name(), 0)
# Q_el = element("P", domain.topology.cell_name(), 1)
VQ = fem.functionspace(domain, mixed_element([V_el, Q_el]))
V, _ = VQ.sub(0).collapse()
# Dirichlet BCs for obstacle, left, top and bottom walls
uD = fem.Function(V)
uD.interpolate(lambda x: (x[1]*(.41 - x[1]), 0.*x[0]))
dirichlet_facets = mesh.locate_entities_boundary(
msh=domain,
dim=1,
marker=lambda x: ~np.isclose(x[0], 2.2) | np.isclose(x[1], 0.) | np.isclose(x[1], .41)
)
dirichlet_dofs = fem.locate_dofs_topological((VQ.sub(0), V), 1, dirichlet_facets)
bc = fem.dirichletbc(uD, dirichlet_dofs, VQ.sub(0))
# Weak form
(u, p) = TrialFunctions(VQ)
(v, q) = TestFunctions(VQ)
a = inner(grad(u), grad(v)) * dx - p * div(v) * dx - q * div(u) * dx
L = inner(fem.Constant(domain, (0., 0.)), v) * dx
lu_solver = {"ksp_type": "preonly", "pc_type": "lu"} # gives inf with P2-P0 elements
# lu_solver["pc_factor_mat_solver_type"] = "mumps" # works
# lu_solver["pc_factor_mat_solver_type"] = "umfpack" # works
problem = LinearProblem(a, L, bcs=[bc], petsc_options=lu_solver)
uph = problem.solve()
uh, ph = uph.sub(0).collapse(), uph.sub(1).collapse()
u_error2_form = fem.form(inner(uh - uD, uh - uD) * dx + inner(grad(uh - uD), grad(uh - uD)) * dx)
u_error2_subdomain = fem.assemble_scalar(u_error2_form)
u_error2 = MPI.COMM_WORLD.reduce(u_error2_subdomain, op=MPI.SUM)
u_min = MPI.COMM_WORLD.reduce(np.min(uh.x.array), op=MPI.MIN)
u_max = MPI.COMM_WORLD.reduce(np.max(uh.x.array), op=MPI.MAX)
if MPI.COMM_WORLD.rank == 0:
print("Velocity error:", np.sqrt(u_error2))
print("Velocity min:", u_min)
print("Velocity max:", u_max)
Using conda with FEniCSx version 0.9.0 and PETSc version 3.22.3
Any ideas on what’s going on here?