Hello!
I am trying to solve the mixed poisson equation based on this demo with periodic boundary conditions in y-direction on the unit square with dolfinx_mpc.
tl;dr:
Non-periodic code works with (RT, DG) and (CG, CG) spaces for (\sigma, u). Periodic code produces non-periodic result with (CG, CG) spaces and does not work at all with (RT, DG) spaces.
Question: How to make it work with periodicity and (RT, DG) spaces?
I made this MWE for a non-periodic case to show that in principle my code works:
File non-periodic.py:
from mpi4py import MPI
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, io, mesh
from ufl import Measure, TestFunctions, TrialFunctions, div, inner
from dolfinx.fem.petsc import LinearProblem
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"
}
msh = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
# RT and DG spaces are used. But this example also works with CG for both elements.
k = 2
VE = element("RT", msh.basix_cell(), k, shape=(msh.geometry.dim,), dtype=default_real_type)
SE = element("DG", msh.basix_cell(), k - 1, dtype=default_real_type)
ME = mixed_element([VE, SE])
V = fem.functionspace(msh, ME)
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
f = 1.0 # Source
dx = Measure("dx", msh)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx
problem = LinearProblem(a, L, petsc_options=petsc_options)
w_h = problem.solve()
sigma_h, u_h = w_h.split()
VIS = fem.functionspace(msh, ("CG", 1))
u_vis = fem.Function(VIS)
u_vis.interpolate(u_h)
with io.XDMFFile(msh.comm, "results/u.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(u_vis)
This code has RT elements for the \sigma-space and DG for the u-space. However chosing CG elements for both spaces also works. The result looks as follows:
Now I tried to extend the code to include periodicity, but the result looks exactly the same as in the above screenshot, i.e. the periodic BC is not applied:
File periodic.py:
from mpi4py import MPI
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, io, mesh
from ufl import Measure, TestFunctions, TrialFunctions, div, inner
import numpy as np
from dolfinx_mpc import LinearProblem, MultiPointConstraint
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"
}
def periodic_boundary(x):
return np.isclose(x[1], 1.0)
def periodic_relation(x):
out_x = np.zeros_like(x)
out_x[0] = x[0]
out_x[1] = 1.0 - x[1]
return out_x
msh = mesh.create_unit_square(MPI.COMM_WORLD, 64, 64)
# Setting CG for both elements leads to a non periodic solution as if periodicity
# was never set. Changing to RT and DG elements does not work at all because an
# error is thrown.
k = 2
VE = element("CG", msh.basix_cell(), k, shape=(msh.geometry.dim,), dtype=default_real_type)
SE = element("CG", msh.basix_cell(), k - 1, dtype=default_real_type)
ME = mixed_element([VE, SE])
V = fem.functionspace(msh, ME)
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
f = 1.0
dx = Measure("dx", msh)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx
# Periodicity in y-direction:
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V.sub(0), periodic_boundary, periodic_relation, bcs = []) # sigma
mpc.create_periodic_constraint_geometrical(V.sub(1), periodic_boundary, periodic_relation, bcs = []) # u
mpc.finalize()
problem = LinearProblem(a, L, mpc, bcs = [], petsc_options=petsc_options)
w_h = problem.solve()
sigma_h, u_h = w_h.split()
VIS = fem.functionspace(msh, ("CG", 1))
u_vis = fem.Function(VIS)
u_vis.interpolate(u_h)
with io.XDMFFile(msh.comm, "results/u.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(u_vis)
For the periodic case I chose CG elements for both spaces because using RT and DG causes an error. You can find the error message below.
Error message:
Traceback (most recent call last):
File "/mnt/d/sim/mixed_poisson/question/periodic.py", line 53, in <module>
mpc.create_periodic_constraint_geometrical(V.sub(0), periodic_boundary, periodic_relation, bcs = []) # sigma
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/pbi/miniconda3/envs/fenicsx/lib/python3.13/site-packages/dolfinx_mpc/multipointconstraint.py", line 290, in create_periodic_constraint_geometrical
mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical(
V._cpp_object, indicator, relation, bcs, scale, True
)
RuntimeError: Cannot evaluate dof coordinates - this element does not have pointwise evaluation.
This question is related to one of my other questions: