Applying Dirichlet BCs to a matrix-vector product RHS

General
1

Dear all,

I am a little confused with the apply_lifting function and the apparent need of updating the ghosts with reverse-scatter after calling it.
I’m aware that usually the process of applying Dirichlet BCs is this:

# ... define a, L, bc
A = assemble_matrix(a, bcs=[bc])
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], bcs=[[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])

However, if the RHS vector b is not the result of assemble_vector but is computed as a preassembled matrix-vector product, the “add-reverse scatter” does not apply.
Is it possible to use apply_lifting at all in this case? If not, how to apply BCs instead?

I suspect that the usual application of BCs to the LHS matrix A during its assembly is not an option and that the rows of A have to be manipulated using PETSc directly. If that’s true, I’d appreciate pointers to code examples :slight_smile:

For illustration consider the MWE of a simple Poisson problem where the RHS is constructed by multiplying a mass matrix M with some vector (the line M.mult(b, u.vector)):

import dolfinx
import dolfinx.fem.petsc
import numpy as np
import ufl
from dolfinx import fem
from mpi4py import MPI
from petsc4py import PETSc

N = 32
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)
V = fem.FunctionSpace(mesh, ("CG", 1))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

a = fem.form(ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx)
mass = fem.form(ufl.inner(u, v) * ufl.dx)


def walls(x):
    return np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1))


wall_dofs = fem.locate_dofs_geometrical(V, walls)
bc = fem.dirichletbc(fem.Constant(mesh, PETSc.ScalarType(1.0)), wall_dofs, V)

A = fem.petsc.assemble_matrix(a, bcs=[bc])
A.assemble()

M = fem.petsc.assemble_matrix(mass)
M.assemble()

# create RHS vector b = M*u
u = fem.Function(V)
u.interpolate(lambda x: np.full((x.shape[1],), 3.0))
u.x.scatter_forward()

b = u.vector.duplicate()

M.mult(u.vector, b)

fem.petsc.apply_lifting(b, [a], [[bc]])
# b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(b, [bc])
b.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

ksp = PETSc.KSP().create()
ksp.setType("preonly")
pc = ksp.getPC()
pc.setType("lu")

ksp.setOperators(A)
ksp.solve(b, u.vector)
u.x.scatter_forward()

print(f"rank {MPI.COMM_WORLD.rank}:  {u.x.norm() = }")

The results differ when this is run on 1 or more procs.

Thank you!

2

I would simply assemble the lifting contribution into a separate vector, use scatter-reverse and then add the two vectors together.

1 Like