Hi Justin,

> On 30 Jan 2019, at 20:47, Justin Chang <jychang48@gmail.com> wrote:
>
> Hi guys,
>
> (This might actually be a PETSc question, but I know there are some PETSc folks that also peruse through this mailing list)
>
> I'm looking again through the Preconditioning saddle-point systems demo, one example shown was providing the scour complement approximation for the mixed poisson problem. I see that the bilinear form for the discontinuous laplacian (plus the flux mass term) is set as the preconditioning Mat for KSPSetOperators() and that these PETSc options are utilized:
>
> parameters = {
>     "ksp_type": "gmres",
>     "ksp_rtol": 1e-8,
>     "pc_type": "fieldsplit",
>     "pc_fieldsplit_type": "schur",
>     "pc_fieldsplit_schur_fact_type": "full",
>     "pc_fieldsplit_off_diag_use_amat": True,
>     "fieldsplit_0_ksp_type": "preonly",
>     "fieldsplit_0_pc_type": "ilu",
>     "fieldsplit_1_ksp_type": "preonly",
>     "fieldsplit_1_pc_type": "hypre"
> }
>
> This raises a couple questions from myself:
>
> 1) Why is "pc_fieldsplit_off_diag_use_amat": True necessary? In the aP term, I see a dot(sigma, tau)*dx term in the bilinear, so wouldn't that actually be the amat?

The provided preconditioning matrix only has entries in the diagonal blocks (the 1-1, and 2-2 blocks). But for "full" factorisation, you need the off-diagonal blocks. So off_diag_use_amat says to get those from A, rather than P.

> 2) I don't see -pc_fieldsplit_schur_precondition set anywhere meaning the default is the a11 matrix. In the original bilinear form, a11 would be zero but because you passed in aP, a11 now pick up the laplacian term?

Yes, exactly. A^{-1} and S^{-1} are constructed from pmat (in petsc parlance).

> 3) Can this easily have been done instead using "pc_fieldsplit_schur_precondition": "user" and the user would provide a python context that builds only the schur complement approximation and applies the action?

Yes, in fact, I should update this demo (or perhaps you would like to?) to use the now preferred way of doing things, which is to make a Python PC.

You have multiple options here. If the right thing to do is to just assemble some matrix and invert that (as in this example), then you write (in up to the minute Firedrake) a little python class that inherits from AuxiliaryOperatorPC:

class DGLaplacian(AuxiliaryOperatorPC):
    def form(self, pc, test, trial):
        # make aP here
        return (aP, list_of_any_boundary_conditions_or_None)

Then no need to provide the aP and instead you just say:

parameters = {
  ...
  "fieldsplit_1": {
     "pc_type": "python",
     # This can name any class in a module Python can import
     "pc_python_type": "__main__.DGLaplacian",
     "aux_pc_type": "whatever_you_want"
  }
}

The AuxiliaryOperatorPC assembles the operator you provide in the form method and makes a PC with that operator, using "aux_" as a prefix.

You can see an example here: https://firedrakeproject.org/demos/geometric_multigrid.py.html

If you need to do more complicated things than just provide a form (for example control how you apply the PC) then we have a setup to do this too. You instead inherit from PCBase and need to provide a class:

class MyPC(PCBase):
   def initialize(self, pc):
      # Called for initial setup

   def update(self, pc):
      # Called when the operator changes

   def apply(self, pc, x, y):
      # Compute y <- pc(x)

   def applyTranspose(self, pc, x, y):
      # Compute y <- pc^T(x)
      # You normally don't need this, so can just "raise NotImplementedError"
      # If there's no better idea


There is some discussion in the manual here: https://firedrakeproject.org/matrix-free.html, or with more words in the SISC paper: http://arxiv.org/abs/1706.01346

To see things in action, look at some of the examples in:

https://github.com/firedrakeproject/firedrake/tree/master/firedrake/preconditioners

Or some fancier things in:

https://github.com/firedrakeproject/firedrake/blob/master/firedrake/slate/static_condensation/scpc.py

Hope that's helpful!

Lawrence