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()