Hello everyone, I’m trying to solve the 1D advection equation, given by:
I’d like to add periodic boundary conditions, and from my understanding dolfinx_mpc is the recommended tool to tackle this. However, I’m running into an error when I use the dolfinx_mpc function create_periodic_constraint_topological.
Specifically, the error I receive is:
TypeError: create_periodic_constraint_topological(): incompatible function arguments. The following argument types are supported:
1. (arg0: dolfinx::fem::FunctionSpace, arg1: dolfinx::mesh::MeshTags<int>, arg2: int, arg3: Callable[[numpy.ndarray[numpy.float64]], numpy.ndarray[numpy.float64]], arg4: List[dolfinx::fem::DirichletBC<double>], arg5: float, arg6: bool) -> dolfinx_mpc.cpp.mpc.mpc_data
Invoked with: <dolfinx.cpp.fem.FunctionSpace object at 0x7fd54c4b38f0>, <dolfinx.mesh.MeshTagsMetaClass object at 0x7fd54c35d6c0>, 2, <function periodic_relation at 0x7fd54c3c9510>, [], 1, False
Below is my MWE:
#from this import d
from types import CellType
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem import FunctionSpace
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import (CellType, create_unit_cube, locate_entities_boundary,meshtags)
from dolfinx import geometry
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint
from dolfinx.common import Timer, TimingType, list_timings
from ufl import ds, dx, grad, inner, SpatialCoordinate
from mpi4py import MPI
from petsc4py import PETSc
from matplotlib import pyplot as plt
from collections import OrderedDict
import os
start = 0
stop = 1
num_intervals = 100
delta_x = (start - stop) / num_intervals
t = 0.0
T = 2
num_timesteps = 100
ts = np.linspace(0, 1, num_timesteps+1)
dt = (T - t) / num_timesteps
adv = 0.4
u0 = lambda x: np.sin(x[0])
# u0 = lambda x: x[0] ** 2
degree = 1
domain = mesh.create_unit_interval(MPI.COMM_WORLD, num_intervals)
V = FunctionSpace(domain, ("CG", degree))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
ut = fem.Function(V)
ut.interpolate(u0)
a = (u + dt * adv * ufl.grad(u)[0])*v*ufl.dx
L = ut * v * ufl.dx
uh = fem.Function(V)
uh.interpolate(u0)
bilinear = fem.form(a)
linear = fem.form(L)
bcs = []
def periodic_boundary(x):
return np.array(np.isclose(x[0], 1))
def periodic_relation(x):
out_x = np.zeros(x.shape)
out_x[0] = 1 - x[0]
out_x[1] = x[1]
out_x[2] = x[2]
return out_x
facets = locate_entities_boundary(domain, domain.topology.dim - 1, periodic_boundary)
arg_sort = np.argsort(facets)
mt = meshtags(domain, domain.topology.dim - 1, facets[arg_sort], np.full(len(facets), 2, dtype=np.int64))
with Timer("~PERIODIC: Initialize MPC"):
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V, mt, 2, periodic_relation, bcs, 1.0)
mpc.finalize()
A = fem.petsc.assemble_matrix(bilinear, bcs=bcs)
A.assemble()
b = fem.petsc.create_vector(linear)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
for i in range(num_timesteps):
t += dt
t = round(t, 4)
# Update L reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
fem.petsc.assemble_vector(b, linear)
# set bcs, solve
fem.petsc.set_bc(b, bcs)
solver.solve(b, uh.vector)
uh.x.scatter_forward()
# Update solution at previous time step (ut)
ut.x.array[:] = uh.x.array
The error message seems to suggest that the parameters I pass to create_periodic_constraint_topological function don’t match what is expected, however I’m not able to determine where the mismatch is.
Also, I attempted to locally run the PBC demo at https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic3d_topological.py, and I run into the same exact error without any changes to the demo. Could this be some sort of versioning issue? I’m running dolfinx v0.6.0 and dolfinx_mpc v0.6.1 (both installed via conda) and Python 3.10.12 .
Does anyone have any ideas of what could be going wrong here? I would really appreciate anyone’s help on this matter! Thank you in advance.