Hi Colin,

On 24/04/16 00:28, Colin Cotter wrote:
> Hi Firedrakers,
>
>   I'm revisiting the solver approach in the vertical slice model, to
> try to incorporate it better into dcore. A vertical slice model is
> obtained from a 3D model, but you assume that all the fields are
> independent of y, whilst retaining a fully 3D velocity vector.
>
> There are two ways to achieve this:
> 1) Build a mesh that is only one element wide and periodic in the
> y-direction.
> 2) Split the velocity space into a part which is tangential to the
> slice plane (u,w), and a part which is normal to the slice plane (v).
> In compatible FEM this means making a 2D H(div) space for the
> tangential part, and a scalar DG space for v.
>
> So far we have been doing 2) in a way that requires rather specialised
> code that is different from what we use in 3D, which makes the code
> very bloated and hard to test.
>
> First of all, is 1) possible in Firedrake? If not, how difficult would
> it be to make it possible? This would be my preferred solution.

So you would like a mesh, I'll draw quads:

x---o
|   |
|   |
y---p

Where the vertical edges and the vertices are periodic (i.e. x and o
are topologically the same and similarly y and p).

I think you might be able to do this if you talk directly to the
dmplex API.  First you make the topological cells you want, then you
build a coordinate section (a functionspace) on that DM and set up
your coordinates appropriately.  We'd then have to do a little bit of
plumbing to make sure that your could tell firedrake that the
coordinate field petsc provides is already discontinuous.

> If 1) is out, then I'd like to improve our approach to 2) so that it
> is less intrusive in our code. I'd really like to make a sum of spaces
> so that we get a 3D vector with the in plane components coming from an
> H(div) space and the out-of-plane components coming from a DG space.

This would imply that different basis functions in your element pull
back differently.  That is Right Out at the moment.  If you would
*only* like to treat these parts of the space as (u, v, w) when
solving you could do the following.  Make a 3-variable mixed space:

W = XZ * Y * P

In your solver options do:

'pc_fieldsplit_0_fields': '0, 1',
'pc_fieldsplit_1_fields': '2',
# normal fieldsplit options for 0 and 1 where 0 now refers to the
# XZY space and 1 to the P space.

You have to turn off nested matrices for this to work.

> Then, I could call this whole thing my velocity space and e.g. solve
> mixed velocity pressure problems using it. Is *this* possible? If not,
> how difficult would it be to make it possible?

Basis functions with different pullback is difficult I think, and
requires serious re-engineering all through the form compiler
pipeline.  The solver trick I just showed should "work right now".

Cheers,

Lawrence