[petsc-users] Stokes-Brinkmann equation preconditioner

[Lawrence Mitchell]

> On 4 Oct 2019, at 10:46, Matthew Knepley via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> On Thu, Oct 3, 2019 at 6:34 PM Salazar De Troya, Miguel via petsc-users <petsc-users at mcs.anl.gov> wrote:
> I am trying to solve the Stokes equation with the Brinkman term to simulate a solid material. My intention is to implement the preconditioner in this paper:https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426 (section 2.6) 
> 
> 
> Link does not work for me.

Try https://onlinelibrary.wiley.com/doi/epdf/10.1002/fld.426 (mail clients are awful)

>  
> where they solve for the velocity and substitute that expression in the pressure equation. They end up solving a system of the type B K^-1 B^T, i.e. the Schur complement of the problem. For this system of equations, they argue that the preconditioner in page 11 is perfect for a given constant Brinkman penalty term.
> 
>  
> 
> Because I am solving for velocity and pressure without doing any substitution, I thought I could use a PC fieldsplit type Schur (full factorization)
> 
> 
> Yes, this will form the exact factorization and a matrix-free form of the Schur complement.
>  
> and provide the preconditioner in the paper to solve the Schur complement.
> 
> 
> Yes, you can provide a user-defined PC for the Schur complement.
>  
> My question is, should I provide this preconditioner through PCFieldSplitSetSchurPre or through fieldsplit_1_pc_type (probably through the Firedrake interface as inhttps://www.firedrakeproject.org/demos/stokes.py.html) ?
> 
> 
> The name PCFieldSplitSetSchurPre seems to be very misleading. You do not use it to provide a _preconditioner_. You use it to determine
> the _preconditioning matrix_ from which the actual preconditioner is built. The preconditioner itself is defined using -fieldpslit_1_pc_type.
> Since I do not know what the preconditioner looks like, I cannot say what preconditioner matrix you would want. Since Firedrake can construct
> any operator for you, you might not care about the matrix we pass to you.

The action of the PC is defined by:

M^{-1} q = \phi_q + \mu q

\mu is some number, and \phi_q solves the BVP

-div \alpha^{-1} grad \phi_q = q

grad \phi_q \dot n = 0
\int \phi_q = 0

So your PC for the fieldsplit_1 needs to:

1. Discretise and solve this BVP to compute \phi_q
2. Add \mu q.

I.e. the action of PCApply is:

y <- \mu x + Laplace^{-1} x

(Aside, surely a primal-dual error has been made in the analysis here)

You could, I think, do this by providing the discretisation of this laplacian and using an additive PCComposite (although I don't know what PC you would use to just scale the input, it's easy to write one though).

Cheers,

Lawrence