Hello,
I wonder whether it is possible to use apply_lifting in conjunction with a block system to enforce Dirichlet boundary conditions. I illustrated my issue in the following code example:
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.mesh import create_unit_interval
from dolfinx.fem import functionspace, Function, form, Constant, apply_lifting, set_bc
import dolfinx.fem.petsc as petsc
from basix.ufl import element
from ufl import Measure, TestFunction, inner, grad, derivative
msh = create_unit_interval(MPI.COMM_WORLD, 10)
dx = Measure("dx", domain=msh, metadata={"quadrature_degree":1})
cg1_element = element("Lagrange", msh.basix_cell(), 1)
V = functionspace(msh, cg1_element)
W = functionspace(msh, cg1_element)
v = Function(V)
v_t = Function(V)
v_test = TestFunction(V)
f = Function(V)
v_exact = Function(V)
w = Function(W)
w_t = Function(W)
w_test = TestFunction(W)
g = Function(W)
w_exact = Function(W)
# The following two functions are from
# https://gitlab.com/rafinex-external-rifle/fenicsx-pctools/
def functions_to_vec(u: list[Function], x):
"""Copies functions into block vector."""
if x.getType() == "nest":
for i, subvec in enumerate(x.getNestSubVecs()):
u[i].x.petsc_vec.copy(subvec)
subvec.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
else:
offset = 0
for i in range(len(u)):
size_local = u[i].x.petsc_vec.getLocalSize()
with x.localForm() as loc:
loc.array[offset : offset + size_local] = u[i].vector.array_r
offset += size_local
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
def vec_to_functions(x, u: list[Function]):
"""Copies block vector into functions."""
if x.getType() == "nest":
for i, subvec in enumerate(x.getNestSubVecs()):
subvec.copy(u[i].x.petsc_vec)
u[i].vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
else:
offset = 0
for i in range(len(u)):
size_local = u[i].x.petsc_vec.getLocalSize()
u[i].x.petsc_vec.array[:] = x.array_r[offset : offset + size_local]
offset += size_local
u[i].x.petsc_vec.ghostUpdate(
addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD
)
class PETScTSSolver:
def __init__(self, ts_problem, solution_vars, bcs):
self.solution_vars = solution_vars
self.v = solution_vars[0]
self.w = solution_vars[1]
self.func_space_v = self.v.ufl_function_space()
self.func_space_w = self.w.ufl_function_space()
self.comm = self.func_space_v.mesh.comm
self.v_t = Function(self.func_space_v)
self.w_t = Function(self.func_space_w)
self.sigma_t = Constant(self.func_space_v.mesh, 0.0)
self.f_v = ts_problem.F_v(self.v, self.v_t, self.w, self.w_t)
self.f_w = ts_problem.F_w(self.v, self.v_t, self.w, self.w_t)
self.f_vv = self.sigma_t * ts_problem.dFvdvt(self.v, self.v_t, self.w, self.w_t) + ts_problem.dFvdv(self.v, self.v_t, self.w, self.w_t)
self.f_vw = self.sigma_t * ts_problem.dFvdwt(self.v, self.v_t, self.w, self.w_t) + ts_problem.dFvdw(self.v, self.v_t, self.w, self.w_t)
self.f_wv = self.sigma_t * ts_problem.dFwdvt(self.v, self.v_t, self.w, self.w_t) + ts_problem.dFwdv(self.v, self.v_t, self.w, self.w_t)
self.f_ww = self.sigma_t * ts_problem.dFwdwt(self.v, self.v_t, self.w, self.w_t) + ts_problem.dFwdw(self.v, self.v_t, self.w, self.w_t)
self.bcs = bcs
self.ts = PETSc.TS().create(self.comm)
def evalIFunction(self, ts, t, x, xdot, F):
with F.localForm() as f_local:
f_local.set(0.0)
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
vec_to_functions(x, self.solution_vars)
vec_to_functions(xdot, [self.v_t, self.w_t])
petsc.assemble_vector_block(F, L=[form(self.f_v), form(self.f_w)], a=[[form(self.f_vv), form(self.f_vw)],[form(self.f_wv), form(self.f_ww)]])
apply_lifting(F, a=[[form(self.f_vv), form(self.f_vw)],[form(self.f_wv), form(self.f_ww)]], bcs=[self.bcs], x0=x, alpha=-1.0)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(F, bcs=self.bcs, x0=x, scale=-1.0)
F.assemble()
return True
def solve(self, t0, dt_init, max_t, max_steps):
self.ts.setIFunction(self.evalIFunction)
x0 = petsc.create_vector_block(form([self.f_v, self.f_w]))
self.ts.setSolution(x0)
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(x0)
return True
class WeakFormPDE():
def F_v(self, v, v_t, w, w_t):
self.v = v
self.v_t = v_t
self.w = w
self.w_t = w_t
F_v = self.v_t * v_test * dx + inner(grad(v_test), grad(self.v)) * dx + self.w * v_test * dx
return F_v
def F_w(self, v, v_t, w, w_t):
self.v = v
self.v_t = v_t
self.w = w
self.w_t = w_t
F_w = self.w_t * w_test * dx + inner(grad(w_test), grad(self.w)) * dx + self.v * w_test * dx
return F_w
def dFvdv(self, v, v_t, w, w_t):
return derivative(self.F_v(v, v_t, w, w_t), v)
def dFvdvt(self, v, v_t, w, w_t):
return derivative(self.F_v(v, v_t, w, w_t), v_t)
def dFvdw(self, v, v_t, w, w_t):
return derivative(self.F_v(v, v_t, w, w_t), w)
def dFvdwt(self, v, v_t, w, w_t):
return derivative(self.F_v(v, v_t, w, w_t), w_t)
def dFwdv(self, v, v_t, w, w_t):
return derivative(self.F_w(v, v_t, w, w_t), v)
def dFwdvt(self, v, v_t, w, w_t):
return derivative(self.F_w(v, v_t, w, w_t), v_t)
def dFwdw(self, v, v_t, w, w_t):
return derivative(self.F_w(v, v_t, w, w_t), w)
def dFwdwt(self, v, v_t, w, w_t):
return derivative(self.F_w(v, v_t, w, w_t), w_t)
Dirichlet_bcs = []
petsc_ts_solver = PETScTSSolver(WeakFormPDE(), [v, w], Dirichlet_bcs)
petsc_ts_solver.solve(t0=0, dt_init=1e-3, max_t=1, max_steps=1e3)
While executing apply_lifting, the above code causes an exception, which I do not understand. When removing the two lines
apply_lifting(F, a=[[form(self.f_vv), form(self.f_vw)],[form(self.f_wv), form(self.f_ww)]], bcs=[self.bcs], x0=x, alpha=-1.0)
set_bc(F, bcs=self.bcs, x0=x, scale=-1.0)
from the code it works. However, then also the lifting procedure is not applied. Any ideas what needs t be changed here?
Thank you for your help.