Source code for firedrake.preconditioners.massinv

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

__all__ = ("MassInvPC", )




class MassInvPC(PCBase):

    needs_python_pmat = True

    """A matrix free operator that inverts the mass matrix in the provided space.

    Internally this creates a PETSc KSP object that can be controlled
    by options using the extra options prefix ``Mp_``.

    For Stokes problems, to be spectrally equivalent to the Schur
    complement, the mass matrix should be weighted by the viscosity.
    This can be provided (defaulting to constant viscosity) by
    providing a field defining the viscosity in the application
    context, keyed on ``"mu"``.
    """


    def initialize(self, pc):
        from firedrake import TrialFunction, TestFunction, dx, assemble, inner, parameters
        prefix = pc.getOptionsPrefix()
        options_prefix = prefix + "Mp_"
        # 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("MassInvPC only makes sense if test and trial space are the same")

        V = test.function_space()

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

        u = TrialFunction(V)
        v = TestFunction(V)
        # Handle vector and tensor-valued spaces.

        # 1/mu goes into the inner product in case it varies spatially.
        a = inner(1/mu * u, v)*dx

        opts = PETSc.Options()
        mat_type = opts.getString(options_prefix + "mat_type",
                                  parameters["default_matrix_type"])

        A = assemble(a, form_compiler_parameters=context.fc_params,
                     mat_type=mat_type, options_prefix=options_prefix)

        Pmat = A.petscmat
        Pmat.setNullSpace(P.getNullSpace())
        tnullsp = P.getTransposeNullSpace()
        if tnullsp.handle != 0:
            Pmat.setTransposeNullSpace(tnullsp)

        ksp = PETSc.KSP().create(comm=pc.comm)
        ksp.incrementTabLevel(1, parent=pc)
        ksp.setOperators(Pmat)
        ksp.setOptionsPrefix(options_prefix)
        ksp.setFromOptions()
        ksp.setUp()
        self.ksp = ksp




    def update(self, pc):
        pass




    def apply(self, pc, X, Y):
        self.ksp.solve(X, Y)


    # Mass matrix is symmetric
    applyTranspose = apply



    def view(self, pc, viewer=None):
        super(MassInvPC, self).view(pc, viewer)
        viewer.printfASCII("KSP solver for M^-1\n")
        self.ksp.view(viewer)




    def destroy(self, pc):
        if hasattr(self, "ksp"):
            self.ksp.destroy()