Dear community,
I’m solving a 2D Stokes problem with a Brinkman-type penalization using FEniCSx and dolfinx_mpc.
The domain is the unit square [0,1]×[0,1] with fully periodic boundary conditions applied via multipoint constraints mpc. Velocity and pressure are both discretized with P1 elements, with a stabilization. The system is assembled using nested matrices and solved with a field-split preconditioner, following some topics posted on the forum previously (for example).
While periodicity works as expected on the edges, I’m encountering issues at the domain corners: the DOFs there seem not to be properly constrained, resulting in discontinuities or mismatches. I’ve implemented a third periodic constraint to handle the corners, but the issue persists. Any ideas on what might be going wrong or how to correctly handle corner DOFs with dolfinx_mpc?
Here I attached the full code, thank you in advance.
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from dolfinx import fem, mesh
from dolfinx.fem import Function, functionspace, Constant, dirichletbc, form, assemble_scalar
from dolfinx.mesh import create_rectangle, CellType, DiagonalType, locate_entities_boundary, meshtags
from dolfinx.io import XDMFFile
from dolfinx_mpc import (MultiPointConstraint,
assemble_matrix, assemble_vector,
create_matrix_nest, create_vector_nest,
assemble_matrix_nest, assemble_vector_nest)
# --- Periodic MPC utility functions ---
def define_periodic_space_bcs(domain, Lx, Ly, f_space, dim, value):
toll = 1e-6
point_dof = fem.locate_dofs_geometrical(
f_space, lambda x: np.isclose(x[0], 0.0, atol=toll) & np.isclose(x[1], 0.0, atol=toll)
)
if dim == 1:
return [dirichletbc(Constant(domain, PETSc.ScalarType(value)), point_dof, f_space)]
else:
return [dirichletbc(np.zeros(dim, dtype=PETSc.ScalarType), point_dof, f_space)]
def define_periodic_space(domain, Lx, Ly, f_space, bcs):
toll = 1e-6
fdim = domain.topology.dim - 1
mpc = MultiPointConstraint(f_space)
def is_slave_x(x): return np.isclose(x[0], Lx, atol=toll)
def is_master_x(x): return np.isclose(x[0], 0.0, atol=toll)
def is_slave_y(x): return np.isclose(x[1], Ly, atol=toll)
def is_master_y(x): return np.isclose(x[1], 0.0, atol=toll)
def map_slave_x(x):
out = x.copy()
out[0] -= Lx
idx_corner = np.logical_or(is_slave_y(x), is_master_y(x))
out[:, idx_corner] = np.nan
return out
def map_slave_y(x):
out = x.copy()
out[1] -= Ly
idx_corner = np.logical_or(is_slave_x(x), is_master_x(x))
out[:, idx_corner] = np.nan
return out
def map_corner(x):
out = x.copy()
idx = np.logical_and(is_slave_x(x), is_slave_y(x))
out[:, ~idx] = np.nan
out[0, idx] -= Lx
out[1, idx] -= Ly
return out
tag_x, tag_y, tag_corner = 2, 3, 4
facets_x = locate_entities_boundary(domain, fdim, is_slave_x)
facets_y = locate_entities_boundary(domain, fdim, is_slave_y)
mt_x = meshtags(domain, fdim, facets_x, np.full(len(facets_x), tag_x, dtype=np.int32))
mt_y = meshtags(domain, fdim, facets_y, np.full(len(facets_y), tag_y, dtype=np.int32))
mpc.create_periodic_constraint_topological(f_space, mt_x, tag_x, map_slave_x, bcs)
mpc.create_periodic_constraint_topological(f_space, mt_y, tag_y, map_slave_y, bcs)
facets_corner = locate_entities_boundary(domain, fdim,
lambda x: np.isclose(x[0], Lx, atol=toll) & np.isclose(x[1], Ly, atol=toll)
)
mt_corner = meshtags(domain, fdim, facets_corner,
np.full(len(facets_corner), tag_corner, dtype=np.int32))
mpc.create_periodic_constraint_topological(f_space, mt_corner, tag_corner, map_corner, bcs)
mpc.finalize()
return mpc
# --- Main execution ---
if __name__ == "__main__":
# Create mesh
Lx = Ly = 1.0
n = 80
domain = create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [Lx, Ly]], [n, n], CellType.triangle)
# Define function spaces
V = functionspace(domain, ("Lagrange", 1, (2,)))
Q = functionspace(domain, ("Lagrange", 1))
# Periodic MPCs
bcs_V = define_periodic_space_bcs(domain, Lx, Ly, V, 2, 0.0)
bcs_Q = define_periodic_space_bcs(domain, Lx, Ly, Q, 1, 0.0)
mpc_V = define_periodic_space(domain, Lx, Ly, V, bcs_V)
mpc_Q = define_periodic_space(domain, Lx, Ly, Q, bcs_Q)
# Trial and test functions
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
p, q = ufl.TrialFunction(Q), ufl.TestFunction(Q)
# Density and material interpolation
rho = Function(functionspace(domain, ("DG", 0)))
rho.interpolate(lambda x: np.where((x[0] - 0.5)**2 + (x[1] - 0.5)**2 < 0.2 ** 2, 1.0, 0.0))
z_min, z_max = 0.0, 1e6
qRAMP = 0.01
zeta = z_max + (z_min - z_max) * (1 - rho) * (1 + qRAMP) / (1 - rho + qRAMP)
dx = ufl.dx(domain)
alpha_stab = 1.0 / 6.0
# Bilinear forms
a_uu = (ufl.inner(ufl.grad(u), ufl.grad(v)) + zeta * ufl.inner(u, v)) * dx
a_up = -ufl.inner(ufl.grad(p), v) * dx
a_pu = -ufl.inner(ufl.grad(q), u) * dx
a_pp = -alpha_stab * ufl.CellDiameter(domain)**2 * ufl.inner(ufl.grad(p), ufl.grad(q)) * dx
a = [[form(a_uu), form(a_up)], [form(a_pu), form(a_pp)]]
# RHS for unit forcing in x-direction
Lx_u = ufl.inner(ufl.as_vector([1.0, 0.0]), v) * dx
Lp = Constant(domain, PETSc.ScalarType(0.0)) * q * dx
L = [form(Lx_u), form(Lp)]
# Assemble system matrix
A = create_matrix_nest(a, [mpc_V, mpc_Q])
A00 = A.getNestSubMatrix(0, 0)
A01 = A.getNestSubMatrix(0, 1)
A10 = A.getNestSubMatrix(1, 0)
A11 = A.getNestSubMatrix(1, 1)
assemble_matrix(a[0][0], (mpc_V, mpc_V), bcs=[], A=A00)
assemble_matrix(a[0][1], (mpc_V, mpc_Q), bcs=[], A=A01)
assemble_matrix(a[1][0], (mpc_Q, mpc_V), bcs=bcs_Q, A=A10)
assemble_matrix(a[1][1], (mpc_Q, mpc_Q), bcs=bcs_Q, A=A11)
for block in [A00, A01, A10, A11]:
block.assemble()
A.assemble()
# Assemble RHS
b = create_vector_nest(L, [mpc_V, mpc_Q])
assemble_vector_nest(b, L, [mpc_V, mpc_Q])
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# Create preconditioner (block-diagonal)
P = PETSc.Mat().createNest([[A00, None], [None, A11]])
P.assemble()
# Solve system
x = b.copy()
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A, P)
ksp.setType("fgmres")
pc = ksp.getPC()
pc.setType("fieldsplit")
pc.setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
nested_IS = P.getNestISs()
pc.setFieldSplitIS(("u", nested_IS[0][0]), ("p", nested_IS[0][1]))
ksp_u, ksp_p = pc.getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("hypre")
ksp_u.getPC().setHYPREType("boomeramg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("hypre")
ksp_p.getPC().setHYPREType("boomeramg")
ksp.setTolerances(rtol=1e-8, atol=1e-8)
ksp.setFromOptions()
ksp.solve(b, x)
for x_sub in x.getNestSubVecs():
x_sub.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
# Extract velocity and pressure
u_sol = Function(mpc_V.function_space)
u_sol.x.petsc_vec.setArray(x.getNestSubVecs()[0].array)
u_sol.x.scatter_forward()
mpc_V.backsubstitution(u_sol)
p_sol = Function(mpc_Q.function_space)
p_sol.x.petsc_vec.setArray(x.getNestSubVecs()[1].array)
p_sol.x.scatter_forward()
mpc_Q.backsubstitution(p_sol)
# Output
with XDMFFile(domain.comm, "velocity.xdmf", "w") as xdmf_u:
xdmf_u.write_mesh(domain)
xdmf_u.write_function(u_sol)
with XDMFFile(domain.comm, "pressure.xdmf", "w") as xdmf_p:
xdmf_p.write_mesh(domain)
xdmf_p.write_function(p_sol)