[firedrake] High-contrast Stokes problem

[Lawrence Mitchell]

Hi Loic,
> On 6 Feb 2017, at 10:31, Loic Gouarin <loic.gouarin at math.u-psud.fr> wrote:
> 
> Hi,
> 
> I would like to solve a Stokes problem with a viscosity which is not constant and could be high in some parts of the domain (as in sinker problem). The solvers that I want to use are in this article: 
> 
> http://dx.doi.org.ezproxy.math.cnrs.fr/10.1016/j.cam.2013.10.016
> 
> I would like to use a P0 interpolation for the viscosity function. 
> 
> Could you tell me how can I set this function into my variational formulation ?


You can make a P0 space into which you interpolate your viscosity, and then use that in the form:

P0 = FunctionSpace(mesh, "DG", 0)

nu = Function(P0)

nu.interpolate(whatever)

Now use nu as normal:

a = nu*inner(grad(u), grad(v))*dx ...


> The other question is to construct the S matrix when you use the Schur complement as a preconditioner. I have to set the pressure mass matrix for S and scale it with the inverse of the viscosity. How can I do that ?

You can provide a separate form that is used to construct the preconditioning operator.  Say something like:

aP = nu*inner(grad(u), grad(v))*dx + (1/nu)*p*q*dx

solve(a == L, w, Jp=aP)

And pass solver parameters:

-pc_type fieldsplit
-pc_fieldsplit_type schur
# We didn't provide off diagonal blocks in the preconditioning
# matrix
-pc_use_amat
# Use a11 block from aP to provide preconditioning operator for S.
-pc_fieldsplit_schur_precondition a11

If you have really high contrast (and/or many sinkers), you will probably find that the viscosity-weighted pressure mass matrix gives mesh independent convergence, but quite high iteration counts.  In those case I think you want one of the "BFBt"-like preconditioners (e.g. https://arxiv.org/abs/1607.03936)

Cheers,

Lawrence