There are several issues in this code.
- Trying to use an MPC to map to a DirichletBC:
If the solution is periodic, there should be no DirichletBC, i.e.
def dirichlet_boundary_periodic(
x: NDArray[Union[np.float32, np.float64]],
) -> NDArray[np.bool_]:
return np.full_like(x[0], False, dtype=bool)
You need to be careful with the edge mappings. I attach a code that produces a sensible (and consistent) formulation:
from __future__ import annotations
from typing import Union
from mpi4py import MPI
import dolfinx.fem as fem
import numpy as np
from dolfinx import default_scalar_type, plot
from dolfinx.common import Timer
from dolfinx.fem import Function, assemble_scalar, form
from dolfinx.mesh import (
CellType,
create_unit_cube,
locate_entities_boundary,
meshtags,
exterior_facet_indices,
)
from numpy.typing import NDArray
from ufl import (
SpatialCoordinate,
TestFunction,
TrialFunction,
derivative,
dx,
grad,
inner,
cos,
sin,
pi,
div,
)
from dolfinx_mpc import NonlinearProblem, LinearProblem, MultiPointConstraint
def benchmark_periodic3D(celltype: CellType, N: int, periodic: bool):
# Create mesh and finite element
mesh = create_unit_cube(MPI.COMM_WORLD, N, N, N, celltype)
tdim = mesh.topology.dim
V = fem.functionspace(mesh, ("Lagrange", 2))
tol = float(5e1 * np.finfo(default_scalar_type).resolution)
print(
f"\nTolerance: {tol}, Celltype: {celltype}, Nº cells: {mesh.topology.index_map(mesh.topology.dim).size_local}"
)
def dirichlet_boundary_total(
x: NDArray[Union[np.float32, np.float64]],
) -> NDArray[np.bool_]:
return np.logical_or(
np.logical_or(
np.logical_or(
np.isclose(x[0], 0.0, atol=tol), np.isclose(x[0], 1.0, atol=tol)
),
np.logical_or(
np.isclose(x[1], 0.0, atol=tol), np.isclose(x[1], 1.0, atol=tol)
),
),
np.logical_or(
np.isclose(x[2], 0.0, atol=tol), np.isclose(x[2], 1.0, atol=tol)
),
)
def dirichlet_boundary_periodic(
x: NDArray[Union[np.float32, np.float64]],
) -> NDArray[np.bool_]:
return np.full_like(x[0], False, dtype=bool)
# Create Dirichlet boundary condition
def u_ex(x):
return (
np.cos(2 * np.pi * x[0])
* np.cos(4 * np.pi * x[1])
* np.cos(2 * np.pi * x[2])
)
u_L = Function(V)
u_L.interpolate(u_ex)
dirichletboundary = (
dirichlet_boundary_periodic if periodic else dirichlet_boundary_total
)
geometrical_dofs = fem.locate_dofs_geometrical(V, dirichletboundary)
bc = fem.dirichletbc(u_L, geometrical_dofs) # , V)
bcs = [bc]
print(f"Nº Dirichlet dofs: {len(geometrical_dofs)}")
def PeriodicBoundary(axis):
"""
Full surface minus dofs constrained by BCs
"""
def complicated_mpc(x):
condition = np.isclose(x[axis], 1.0, atol=tol)
if axis == 0:
return condition
elif axis == 1:
return np.logical_and(
condition, np.logical_and(x[0] > tol, x[0] < 1.0 - tol)
)
elif axis == 2:
return np.logical_and(
condition,
np.logical_and(
np.logical_and(x[0] > tol, x[0] < 1.0 - tol),
np.logical_and(x[1] > tol, x[1] < 1.0 - tol),
),
)
else:
raise RuntimeError("Axis should be 0, 1, or 2")
return complicated_mpc
def create_periodic_map(axis):
def periodic_relation(x):
out_x = x.copy()
out_x[axis] = 0
return out_x
return periodic_relation
with Timer("~~Periodic: Compute mpc condition"):
mpc = MultiPointConstraint(V)
for axis in range(3):
mpc.create_periodic_constraint_geometrical(
V, PeriodicBoundary(axis), create_periodic_map(axis), bcs, tol=tol
)
mpc.finalize()
# Define variational problem
u = Function(mpc.function_space, name="Solution_u")
v = TestFunction(V)
x = SpatialCoordinate(mesh)
u_soln = cos(2 * pi * x[0]) * cos(4 * pi * x[1]) * cos(2 * pi * x[2])
f = -div(grad(u_soln))
F = inner(grad(u), grad(v)) * dx - inner(f, v) * dx
u_ = TrialFunction(V)
a = inner(grad(u_), grad(v)) * dx
L = inner(f, v) * dx
# Sanity check that the MPC class has some constraints to impose
num_slaves_global = mesh.comm.allreduce(len(mpc.slaves), op=MPI.SUM)
num_masters_global = mesh.comm.allreduce(len(mpc.masters.array), op=MPI.SUM)
assert (not periodic) | ((num_slaves_global > 0) & (num_masters_global > 0))
print(
f"Periodicity: {periodic}, Nº slaves: {num_slaves_global}, Nº masters: {num_masters_global}"
)
tol = np.finfo(u.x.array.dtype).resolution
petsc_options_nonlinear = {
"snes_type": "ksponly",
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"snes_atol": tol,
"snes_rtol": tol,
"ksp_error_if_not_converged": True,
"snes_error_if_not_converged": True,
"snes_monitor": None,
}
petsc_options_linear = {
"ksp_type": "preonly",
"mat_mumps_icntl_24": 1,
"mat_mumps_icntl_25": 0,
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"ksp_error_if_not_converged": True,
}
problem = LinearProblem(a, L, mpc, bcs=bcs, petsc_options=petsc_options_linear)
u_h = problem.solve()
u_avg_analytical = mesh.comm.allreduce(
assemble_scalar(form(u_soln * dx)), op=MPI.SUM
)
u_avg = mesh.comm.allreduce(assemble_scalar(form(u_h * dx)), op=MPI.SUM)
u_h.x.array[:] -= u_avg - u_avg_analytical
import dolfinx
with dolfinx.io.VTXWriter(
mesh.comm, f"linear_mpc_{periodic=}_{N}_3D.bp", [u_h]
) as writer:
writer.write(0.0)
l2_error_local = assemble_scalar(form((u_h - u_soln) ** 2 * dx))
l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM)
print(f" - Error L2 local: {l2_error_local}")
print(f" - Error L2 global: {l2_error_global}")
print("NONLINEAR:")
J = derivative(F, u, TrialFunction(V))
problem = NonlinearProblem(
F, u, mpc=mpc, bcs=bcs, J=J, petsc_options=petsc_options_nonlinear
)
_, converged, num_iterations = problem.solve()
print(f" - Converged: {converged}, iterations: {num_iterations}")
l2_error_local = assemble_scalar(form((u - u_soln) ** 2 * dx))
l2_error_global = mesh.comm.allreduce(l2_error_local, op=MPI.SUM)
print(f" - Error L2 local: {l2_error_local}")
print(f" - Error L2 global: {l2_error_global}")
with dolfinx.io.VTXWriter(
mesh.comm, f"nonlinear_mpc_{periodic=}_{N}_3D.bp", [u]
) as writer:
writer.write(0.0)
if __name__ == "__main__":
celltypes = [CellType.tetrahedron, CellType.hexahedron]
Ne = [5, 10, 20]
for celltype, N in zip(celltypes, Ne):
for periodic in [True, False]:
benchmark_periodic3D(celltype, N, periodic)
print("\n")