Here is a MWE that reproduces the issue:
# In this minimal working example, we try to solve the PDE system
# \partial_t u -\nabla**2 u + v - f(x,t) = 0
# \partial_t v -\nabla**2 v + u - g(x,t) = 0
# with Dirichlet boundary conditions.
# Method of manufactured solutions:
# v(x,t) = -(x-0.5)**2 + 0.25 * exp(-t)
# w(x,t) = -(x-0.5)**2 - 0.25 * exp(-2*t)
# yields
# f(x,t) = -0.25*exp(-t) + 2 + (x-0.5)**2 - 0.25*exp(-2*t)
# g(x,t) = -0.5*exp(-2*t) - 2 - (x-0.5)**2 + 0.25*exp(-t)
from ufl import Measure, TestFunctions, split, derivative, inner, grad
from dolfinx import mesh
from dolfinx.mesh import create_unit_interval
from basix.ufl import element, mixed_element
from dolfinx.fem import functionspace, Function, Constant, form, locate_dofs_topological, dirichletbc, apply_lifting, set_bc
import dolfinx.fem.petsc as petsc
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
class v_exact_fun():
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = -(x[0]-0.5)**2 + 0.25*np.exp(-self.t)
return values
class w_exact_fun():
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = (x[0]-0.5)**2 - 0.25*np.exp(-2*self.t)
return values
class f_fun():
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = -0.25*np.exp(-self.t) + 2 + (x[0]-0.5)**2 - 0.25*np.exp(-2*self.t)
return values
class g_fun():
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 0.5*np.exp(-2*self.t) - 2 - (x[0]-0.5)**2 + 0.25*np.exp(-self.t)
return values
class PETScTSSolver:
# F(t, u, u') = 0
def __init__(self, ts_problem, u, bcs):
self.u = u
self.func_space = u.ufl_function_space()
self.comm = self.func_space.mesh.comm
self.u_t = Function(self.func_space)
self.sigma_t = Constant(self.func_space.mesh, 0.0)
self.f = ts_problem.F(self.u, self.u_t) # Residual (linear form)
self.f_j = self.sigma_t * ts_problem.dFdut(self.u, self.u_t) + ts_problem.dFdu(self.u, self.u_t) # Jacobian
self.functional_vector_petsc = petsc.create_vector(form(self.f))
self.Jacobian_matrix_petsc = petsc.create_matrix(form(self.f_j))
self.bcs = bcs
self.ts = PETSc.TS().create(self.comm)
def evalIFunction(self, ts, t, x, xdot, F):
bc_v_l_fun.interpolate(v_exact_fun(t))
f.interpolate(f_fun(t))
g.interpolate(g_fun(t))
x.copy(self.u.x.petsc_vec)
xdot.copy(self.u_t.x.petsc_vec)
with F.localForm() as f_local:
f_local.set(0.0)
petsc.assemble_vector(F, form(self.f))
apply_lifting(F, a=[form(self.f_j)], bcs=[self.bcs], x0=[self.u.x.petsc_vec], alpha=-1.0)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(F, bcs=self.bcs, x0=self.u.x.petsc_vec, scale=-1.0)
F.assemble()
return True
def evalIJacobian(self, ts, t, x, xdot, sigma_t, J, P):
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.x.petsc_vec)
self.u.x.petsc_vec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
xdot.copy(self.u_t.x.petsc_vec)
self.sigma_t.value = sigma_t
J.zeroEntries()
petsc.assemble_matrix(J, form(self.f_j), self.bcs)
J.assemble()
return True
def solve(self, t0, dt_init, max_t, max_steps):
y = Function(self.func_space)
y.sub(0).interpolate(v_exact_fun(0))
y.sub(1).interpolate(w_exact_fun(0))
self.ts.setIFunction(self.evalIFunction, f=self.functional_vector_petsc)
self.ts.setIJacobian(self.evalIJacobian, J=self.Jacobian_matrix_petsc)
self.ts.setSolution(y.x.petsc_vec)
self.ts.setTime(t0)
self.ts.setTimeStep(dt_init)
self.ts.setMaxTime(max_t)
self.ts.setType(PETSc.TS.Type.BDF)
self.ts.setTolerances(atol=1e-5, rtol=1e-5)
self.ts.setMaxSteps(max_steps)
self.ts.solve(y.x.petsc_vec)
y.x.petsc_vec.copy(self.u.x.petsc_vec)
return True
class WeakFormPDE():
def F(self, u, u_t):
self.v, self.w = split(u)
self.v_t, self.w_t = split(u_t)
F_v = self.v_t * v_test * dx + inner(grad(v_test), grad(self.v)) * dx + self.w * v_test * dx - f * v_test * dx
F_w = self.w_t * w_test * dx + inner(grad(w_test), grad(self.w)) * dx + self.v * w_test * dx - g * w_test * dx
return F_v + F_w
def dFdu(self, u, u_t):
return derivative(self.F(u, u_t), u)
def dFdut(self, u, u_t):
return derivative(self.F(u, u_t), u_t)
msh = create_unit_interval(MPI.COMM_WORLD, 9)
dx = Measure("dx", domain=msh, metadata={"quadrature_degree":1})
cg1_elem = element("Lagrange", msh.basix_cell(), 1)
mix_elem = mixed_element(2 * [cg1_elem])
mix_elem_space = functionspace(msh, mix_elem)
v_space, v_dofmap = mix_elem_space.sub(0).collapse()
w_space, w_dofmap = mix_elem_space.sub(1).collapse()
v_test, w_test = TestFunctions(mix_elem_space)
u = Function(mix_elem_space)
u_t = Function(mix_elem_space)
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_connectivity(fdim, tdim)
boundary_conditions = [
(lambda x: np.isclose(x[0], 0), 'bdry_l'),
(lambda x: np.isclose(x[0], 1), 'bdry_r'),
]
dof_boundaries_v = {}
dof_boundaries_w = {}
for marker, dof_name in boundary_conditions:
facets_bdry = mesh.locate_entities_boundary(msh, dim=fdim, marker=marker)
dof_boundaries_v[dof_name] = locate_dofs_topological(V=mix_elem_space.sub(0), entity_dim=fdim, entities=facets_bdry)
dof_boundaries_w[dof_name] = locate_dofs_topological(V=mix_elem_space.sub(1), entity_dim=fdim, entities=facets_bdry)
bc_v_l_fun = Function(v_space)
bc_v_l = dirichletbc(value=bc_v_l_fun, dofs=dof_boundaries_v['bdry_l'])
Dirichlet_bcs = [bc_v_l]
#Dirichlet_boundary_values = {'bdry_l': 0.0}
#Dirichlet_bcs = []
#for dof_key, value in Dirichlet_boundary_values.items():
# Dirichlet_bcs.append(dirichletbc(value=value, dofs=dof_boundaries_v[dof_key], V=mix_elem_space.sub(0)))
f = Function(v_space)
g = Function(w_space)
v_exact = Function(functionspace(msh, cg1_elem))
w_exact = Function(functionspace(msh, cg1_elem))
PETScTSSolver_instance = PETScTSSolver(WeakFormPDE(), u, Dirichlet_bcs)
PETScTSSolver_instance.solve(t0=0, dt_init=1e-10, max_t=1, max_steps=1e5)
As stated in my first post, when using zero Dirichlet boundary conditions via
Dirichlet_boundary_values = {'bdry_l': 0.0}
Dirichlet_bcs = []
for dof_key, value in Dirichlet_boundary_values.items():
Dirichlet_bcs.append(dirichletbc(value=value, dofs=dof_boundaries_v[dof_key], V=mix_elem_space.sub(0)))
and disabling bc_v_l_fun.interpolate(v_exact_fun(t)) within evalIFunction instead of using
bc_v_l_fun = Function(v_space)
bc_v_l = dirichletbc(value=bc_v_l_fun, dofs=dof_boundaries_v['bdry_l'])
Dirichlet_bcs = [bc_v_l]
it works.
Thank you for your help.