How to apply Dirichlet boundary conditions using a matrix free method

Question

In a matrix-free method, there is no step like A = assemble(a) hence no DirichletBC.apply(A). The PETScKrylovSolver takes no boundary condition argument. How to apply a Dirichlet boundary condition in this case?

(example from http://fenicsproject.org/qa/7541/ufl-unsupported-operand-type-s-for-form-and-form)

a = inner(grad(u), grad(v))* dx
b = f * v * ds(1)
p = a + u * v * dx

A = assemble(a)
B = assemble(b)
P = assemble(p)

class CompositeOperator(object):
    def __init__(self, A, B):
        self.A_mat = as_backend_type(A).mat()
        self.B_vec = as_backend_type(B).vec()
    def mult(self, mat, x, y):
        self.A_mat.mult(x, y)
        a = self.B_vec.dot(x)
        y.axpy(a, self.B_vec)

    def as_petscmatrix(self):
        from petsc4py import PETSc
        mat = PETSc.Mat().createPython(self.A_mat.getSizes(), 
                                       comm = mpi_comm_world())

        mat.setPythonContext(self)
        return PETScMatrix(mat)


A_B = CompositeOperator(A, B).as_petscmatrix()

solver = PETScKrylovSolver("cg", "hypre_amg")
solver.set_operators(A_B, P)

Comments

Hi, you could consider the get_boundary_values of the DirichletBC class and use the map(boundary dofs and boundary values) to modify the matvec/mult methods of your operators.

Answer 1

When you apply DirichletBC to a matrix, the diagonal entries for the relevant dofs are set to 1, and all off-diagonal entries in the same rows are set to zero. For your operator you get the same effect by applying DirichletBC to the matrix A and modifying the operator to take two vectors, with homogenized BC applied to one of them.

Note that this approach breaks symmetry.

a = inner(grad(u), grad(v))* dx
b = f * v * ds(1)
p = a + u * v * dx

A = assemble(a)
B = assemble(b)
P = assemble(p)

class CompositeOperator(object):
    def __init__(self, A, B, C):
        self.A_mat = as_backend_type(A).mat()
        self.B_vec = as_backend_type(B).vec()
        self.C_vec = as_backend_type(C).vec()

    def mult(self, mat, x, y):
        self.A_mat.mult(x, y)
        a = self.C_vec.dot(x)
        y.axpy(a, self.B_vec)

    def as_petscmatrix(self):
        from petsc4py import PETSc
        mat = PETSc.Mat().createPython(self.A_mat.getSizes(), 
                                       comm = mpi_comm_world())

        mat.setPythonContext(self)
        return PETScMatrix(mat)

# Apply boundary conditions
bc = DirichletBC(V, 1., "x[1]<DOLFIN_EPS")
bc.apply(A)
bc.apply(P)
C = B.copy()
bc.homogenize()
bc.apply(B)


A_B = CompositeOperator(A, B, C).as_petscmatrix()

solver = PETScKrylovSolver("cg", "hypre_amg")
solver.set_operators(A_B, P)

You have to apply the DirichletBC to the right hand side before homogenizing it.

Comments

That seems like a good solution. Two problems still though:

1) In my example, I have several matrices A in the CompositeOperator object, so I should apply the boundary condition to one of them, and somehow make the corresponding rows in the others zero (including the diagonal element). Another option would be to multiply the Dirichlet values by the amount of such matrices A. Neither seems like an elegant solution though?

2) I added some lines to make the sample code executable. I also set the Krylov method and preconditioner to default since 'hypre_amg' is not installed. When running the code below, I get:

*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 73.
*** Where:   This error was encountered inside /home/juanlu/miniconda/conda-bld/work/dolfin-1.4.0/dolfin/la/PETScKrylovSolver.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  

Granted, the example does not relate to anything physical, but why it doesn't work? Am I passing the wrong types to the solver? Or is it possible that there is an error in the build for anaconda: the folder /home/juanlu/... is not mine but probably belongs to the maintainer of fenics on anaconda?

Thanks

from dolfin import *
mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))* dx
f = Expression('x[0]*x[1]')
b = f * v * dx
p = a + u * v * dx

A = assemble(a)
B = assemble(b)
P = assemble(p)

class CompositeOperator(object):
    def __init__(self, A, B, C):
        self.A_mat = as_backend_type(A).mat()
        self.B_vec = as_backend_type(B).vec()
        self.C_vec = as_backend_type(C).vec()

    def mult(self, mat, x, y):
        self.A_mat.mult(x, y)
        a = self.C_vec.dot(x)
        y.axpy(a, self.B_vec)

    def as_petscmatrix(self):
        from petsc4py import PETSc
        mat = PETSc.Mat().createPython(self.A_mat.getSizes(),
                                       comm = mpi_comm_world())

        mat.setPythonContext(self)
        return PETScMatrix(mat)

# Apply boundary conditions
bc = DirichletBC(V, 1., "x[1]<DOLFIN_EPS")
bc.apply(A)
bc.apply(P)
C = B.copy()
bc.homogenize()
bc.apply(B)


A_B = CompositeOperator(A, B, C).as_petscmatrix()

solver = PETScKrylovSolver('default', 'default')
solver.set_operators(A_B, P)

solution = Function(V)
solver.solve(solution.vector(), C)

---

1) You can zero out the relevant rows of a matrix with

bc.zero(A)

Alternatively, you could assemble all the integrals into one matrix.

I don't know what you would consider to be a more elegant solution? It breaks symmetry, but usually conjugate gradients still works fine with this kind of asymmetry.

2) The code runs fine for me (dolfin master, petsc 3.5.4, petsc4py 3.5.1).

---

1) Mathematically I have no problems with the solution (my problem is not symmetric to begin with anyway.) What I meant was not elegant in that it requires a user to abandon yet another abstraction level and manually edit rows corresponding to a DirichletBC object. Clearly, it was naive of me to believe the developers had not already provided a method to do just that :-) So thanks!

2) It also works on my computer at home (1.6.0 from FEniCS PPA) so that must be a separate problem.