Thanks Lawrence for your answer. I works perfectly and we can't do more simple !!! I have a last question. I don't want to use a Q2-Q1 elements for the Stokes problem but a 4Q1-Q1 elements which means that linear elements for velocity live on a mesh that is once more refined globally than the elements for pressure. When I do that: N = 64 meshu = UnitSquareMesh(2*N, 2*N, quadrilateral=True) meshp = UnitSquareMesh(N, N, quadrilateral=True) V = VectorFunctionSpace(meshu, "DQ", 1) W = FunctionSpace(meshp, "DQ", 1) Z = V * W I have an error: ValueError: All function spaces must be defined on the same mesh! Thanks, Loic Le 06/02/2017 à 11:59, Lawrence Mitchell a écrit :

Hi Loic,

...

On 6 Feb 2017, at 10:31, Loic Gouarin <loic.gouarin@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 _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Tel: 01 69 15 60 14 http://www.math.u-psud.fr/~gouarin https://github.com/gouarin

Thanks Lawrence. Loic Le 06/02/2017 à 14:19, Lawrence Mitchell a écrit :

Hi Loic,

...

On 6 Feb 2017, at 12:58, Loic Gouarin <loic.gouarin@math.u-psud.fr> wrote:

Thanks Lawrence for your answer.

I works perfectly and we can't do more simple !!!

I have a last question. I don't want to use a Q2-Q1 elements for the Stokes problem but a 4Q1-Q1 elements which means that linear elements for velocity live on a mesh that is once more refined globally than the elements for pressure.

Right now we don't support these kind of "isoQ1" elements. I think it would be good to do so (since, amongst other things, they're good for preconditioning high order elliptic problems). The problem you encounter is that right now there's no support for building problems like this where some of the function spaces come from different meshes.

We have intentions to support this, but it's quite a lot of work to get right.

Cheers,

Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Tel: 01 69 15 60 14 http://www.math.u-psud.fr/~gouarin https://github.com/gouarin