Source code for firedrake.preconditioners.massinv

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

__all__ = ("MassInvPC", )

[docs] 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"``. """
[docs] 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
[docs] def update(self, pc): pass
[docs] def apply(self, pc, X, Y): self.ksp.solve(X, Y)
# Mass matrix is symmetric applyTranspose = apply
[docs] def view(self, pc, viewer=None): super(MassInvPC, self).view(pc, viewer) viewer.printfASCII("KSP solver for M^-1\n") self.ksp.view(viewer)
[docs] def destroy(self, pc): if hasattr(self, "ksp"): self.ksp.destroy()