Hi, I am trying to understand how to run my simulations with periodic boundary conditions in parallel. While I can run this example in parallel, I don’t understand which part I am missing in my simple MWE here. I believe I am missing a step for the MPC setup. It runs fine normally, but in parallel, I get the following error (part of a longer error):
<>/lib/python3.11/site-packages/dolfinx/la/petsc.py", line 189, in assign
dolfinx.la.petsc.assign(u.x.array, x)
_x.array_w[:] = x0
~~~~~~~~~~^^^
_x.array_w[:] = x0
ValueError: could not broadcast input array from shape (8691,) into shape (9861,)
File "<>/lib/python3.11/functools.py", line 909, in wrapper
~~~~~~~~~~^^^
ValueError: could not broadcast input array from shape (8613,) into shape (9792,)
_x.array_w[:] = x0
~~~~~~~~~~^^^
ValueError: could not broadcast input array from shape (8655,) into shape (9129,)
return dispatch(args[0].__class__)(*args, **kw)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "<>/lib/python3.11/site-packages/dolfinx/la/petsc.py", line 189, in assign
_x.array_w[:] = x0
~~~~~~~~~~^^^
ValueError: could not broadcast input array from shape (8625,) into shape (9519,)
MWE setup:
from mpi4py import MPI
import numpy as np
from dolfinx import fem, mesh, io, default_scalar_type
import dolfinx_mpc
import ufl
comm = MPI.COMM_WORLD
# --- Mesh and function space ---
domain = mesh.create_unit_cube(comm, 20, 20, 20, ghost_mode=mesh.GhostMode.shared_facet )
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
V = fem.functionspace(domain, ("Lagrange", 1, (3,)))
# --- Neo-Hookean strain energy ---
E, nu = 100.0, 0.3
mu, lmbda = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu))
# We solve for the fluctuation (which is trivial in this MWE)
u = fem.Function(V, name="fluctuation")
F_macro = fem.Constant(domain, np.eye(3, dtype=default_scalar_type))
F = ufl.variable(F_macro + ufl.grad(u))
J = ufl.det(F)
psi = (mu / 2) * (ufl.inner(F, F) - 3) - mu * ufl.ln(J) + (lmbda / 2) * ufl.ln(J)**2
P = ufl.diff(psi, F)
# --- Weak form ---
v = ufl.TestFunction(V)
residual = ufl.inner(P, ufl.grad(v)) * ufl.dx
# --- Periodic BCs: fix origin corner, tie +faces to -faces ---
bcs = [fem.dirichletbc(
np.zeros(3, dtype=default_scalar_type),
fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0) & np.isclose(x[1], 0) & np.isclose(x[2], 0)),
V
)]
def periodic_relation(x):
out = x.copy()
out[0, np.isclose(x[0], 1)] = 0.0
out[1, np.isclose(x[1], 1)] = 0.0
out[2, np.isclose(x[2], 1)] = 0.0
return out
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(
V,
lambda x: np.isclose(x[0], 1) | np.isclose(x[1], 1) | np.isclose(x[2], 1),
periodic_relation, bcs
)
mpc.finalize()
# --- Solver ---
problem = dolfinx_mpc.NonlinearProblem(
residual, u, bcs=bcs, mpc=mpc,
J=ufl.derivative(residual, u, ufl.TrialFunction(V)),
petsc_options={
"ksp_type": "gmres",
"pc_type": "hypre",
"snes_type": "newtonls",
"snes_rtol": 1e-8,
"snes_atol": 1e-8,
}
)
# --- Homogenized stress forms ---
P_forms = [fem.form(P[i, j] * ufl.dx) for i in range(3) for j in range(3)]
# --- Uniaxial compression ---
for step, eps in enumerate([0.02, 0.05, 0.10, 0.15, 0.20]):
F_macro.value[:] = np.diag([1.0, 1.0 - eps, 1.0]).astype(default_scalar_type)
_, converged, _ = problem.solve()
u.x.scatter_forward()
P_avg = np.array([
comm.allreduce(fem.assemble_scalar(f), MPI.SUM) for f in P_forms
]).reshape(3, 3) # vol = 1.0 for unit cube
if comm.rank == 0:
print(f"eps={eps:.2f}: P_22={P_avg[1,1]:6.2f}, converged={converged > 0}")