Question: applying of BC's in PCD

[ErikM23]

Dear all,

I am new in firedrake and I would like to ask you about applying of boundary conditions in PCD preconditioner pcd.py

I have got a 3 questions:

1. How to transfer the boundary conditions from my program to pcd1.py (here, pcd1.py should be modified version of the origin pcd.py)

2. How to apply BC's (specially, how to apply BC's to PETCs vectors "a, b" in the function "apply(pc, x, y)" - when I use something like " self.bcs.apply(a)" it doesnt work).

3. Is it possible to apply DirichletBC to matrix and get '0' on the diagonal instead of '1' ?

Thank you in advance.

from firedrake.preconditioners.base import PCBase
from firedrake.petsc import PETSc

__all__ = ("PCDPC", )

class PCDPC(PCBase):

    needs_python_pmat = True

    r"""A Pressure-Convection-Diffusion preconditioner for Navier-Stokes.

    This preconditioner approximates the inverse of the pressure schur
    complement for the Navier-Stokes equations by.

    .. math::

       S^{-1} \sim K^{-1} F_p M^{-1}

    Where :math:`K = \\nabla^2`, :math:`M = \mathbb{I}` and
    :math:`F_p = 1/\mathrm{Re} \\nabla^2 + u\cdot\\nabla`.

    The inverse of :math:`K` is approximated by a KSP which can be
    controlled using the options prefix ``pcd_Kp_``.

    The inverse of :math:`M` is similarly approximated by a KSP which
    can be controlled using the options prefix ``pcd_Mp_``.

    :math:`F_p` requires both the Reynolds number and the current
    velocity.  You must provide these with options using the glbaol
    option ``Re`` for the Reynolds number and the prefixed option
    ``pcd_velocity_space`` which should be the index into the full
    space that gives velocity field.

    .. note::

       Currently, the boundary conditions applied to the PCD operator
       are correct for characteristic velocity boundary conditions,
       but sub-optimal for in and outflow boundaries.
    """
    def initialize(self, pc):
        from firedrake import TrialFunction, TestFunction, dx, \
            assemble, inner, grad, split, Constant, parameters
        from firedrake.assemble import allocate_matrix, create_assembly_callable
        if pc.getType() != "python":
            raise ValueError("Expecting PC type python")
        prefix = pc.getOptionsPrefix() + "pcd_"

        # we assume P has things stuffed inside of it
        _, P = pc.getOperators()
        context = P.getPythonContext()

        test, trial = context.a.arguments()
        if test.function_space() != trial.function_space():
            raise ValueError("Pressure space test and trial space differ")

        Q = test.function_space()

        p = TrialFunction(Q)
        q = TestFunction(Q)

        # !!!!!!!!!!!!!!!!!!!! get BC's from the main program:
        # !!!!!!!!!!!!!!!!!!!! self.bcs = context.appctx["bcsPCD"]  ??? Is it correct ???

        mass = p*q*dx

        # Regularisation to avoid having to think about nullspaces.
        stiffness = inner(grad(p), grad(q))*dx + Constant(1e-6)*p*q*dx

        opts = PETSc.Options()
        # we're inverting Mp and Kp, so default them to assembled.
        # Fp only needs its action, so default it to mat-free.
        # These can of course be overridden.
        # only Fp is referred to in update, so that's the only
        # one we stash.
        default = parameters["default_matrix_type"]
        Mp_mat_type = opts.getString(prefix+"Mp_mat_type", default)
        Kp_mat_type = opts.getString(prefix+"Kp_mat_type", default)
        self.Fp_mat_type = opts.getString(prefix+"Fp_mat_type", "matfree")

        Mp = assemble(mass, form_compiler_parameters=context.fc_params,
                      mat_type=Mp_mat_type,
                      options_prefix=prefix + "Mp_")

        # !!!!!!!!!!!!!!!!!!!!!  applying the BC's to 'Kp' : is it enough to add 'bcs = self.bcs' ???
        Kp = assemble(stiffness, form_compiler_parameters=context.fc_params,
                      mat_type=Kp_mat_type,
                      options_prefix=prefix + "Kp_")

        # FIXME: Should we transfer nullspaces over.  I think not.

        Mksp = PETSc.KSP().create(comm=pc.comm)
        Mksp.incrementTabLevel(1, parent=pc)
        Mksp.setOptionsPrefix(prefix + "Mp_")
        Mksp.setOperators(Mp.petscmat)
        Mksp.setUp()
        Mksp.setFromOptions()
        self.Mksp = Mksp

        Kksp = PETSc.KSP().create(comm=pc.comm)
        Kksp.incrementTabLevel(1, parent=pc)
        Kksp.setOptionsPrefix(prefix + "Kp_")
        Kksp.setOperators(Kp.petscmat)
        Kksp.setUp()
        Kksp.setFromOptions()
        self.Kksp = Kksp

        state = context.appctx["state"]

        Re = context.appctx.get("Re", 1.0)

        velid = context.appctx["velocity_space"]

        u0 = split(state)[velid]
        fp = 1.0/Re * inner(grad(p), grad(q))*dx + inner(u0, grad(p))*q*dx

        self.Re = Re
        self.Fp = allocate_matrix(fp, form_compiler_parameters=context.fc_params,
                                  mat_type=self.Fp_mat_type,
                                  options_prefix=prefix + "Fp_")
        self._assemble_Fp = create_assembly_callable(fp, tensor=self.Fp,
                                                     form_compiler_parameters=context.fc_params,
                                                     mat_type=self.Fp_mat_type)
        self._assemble_Fp()
        Fpmat = self.Fp.petscmat
        self.workspace = [Fpmat.createVecLeft() for i in (0, 1)]

    def update(self, pc):
        self._assemble_Fp()

    def apply(self, pc, x, y):
        a, b = self.workspace
        self.Mksp.solve(x, a)
        self.Fp.petscmat.mult(a, b)

        # !!!!!!!!!!!!!!!!!! How to apply bcs to vector 'b'  or 'a' ???
        # !!!!!!!!!!!!!!!!!! 'a' and 'b' don't belong to function_space

        self.Kksp.solve(b, y)

    def applyTranspose(self, pc, x, y):
        a, b = self.workspace
        self.Kksp.solveTranspose(x, b)
        self.Fp.petscmat.multTranspose(b, a)
        self.Mksp.solveTranspose(y, a)

[wence-]

Hi!

Incidentally, we'd be very happy to accept modifications to the current PCD implementation that better treat boundary conditions. I think the best unified explanation and treatment of how they should be implemented is in chapter 3 of Jan Blechta's thesis

How to transfer the boundary conditions from my program to pcd1.py (here, pcd1.py should be modified version of the origin pcd.py)

Your approach of sending them through context.appctx is reasonable.

How to apply BC's (specially, how to apply BC's to PETCs vectors "a, b" in the function "apply(pc, x, y)" - when I use something like " self.bcs.apply(a)" it doesnt work).

You need to turn the PETSc Vec into a function, which necessitates some copying. If you have the space V around then:

tmp = Function(V)
with tmp.dat.vec_wo as v:
    x.copy(v)
# Now tmp contains the value from `x`
bc.apply(tmp)
with tmp.dat.vec_ro as v:
   v.copy(x)
# Now x has the bcs applied.

Is it possible to apply DirichletBC to matrix and get '0' on the diagonal instead of '1' ?

This is not presently supported, but it would be easy. Right now we always set the diagonal to 1 in assemble. We need a way of communicating what the value should be through the interface. Probably the best way is for DirichletBC to gain an optional keyword argument that is used to provide the value that is put on the diagonal.

[ErikM23]

Dear Lawrence,

Thank you for the answer and hints.

Maybe one more question: I am not sure where exactly in assemble I can set the diagonal to 0 instead of 1.

[wence-]

Maybe one more question: I am not sure where exactly in assemble I can set the diagonal to 0 instead of 1.

For example here: https://github.com/firedrakeproject/firedrake/blob/master/firedrake/assemble.py#L447

You would provide the extra keyword argument diag_val=some_value