[firedrake] Parallel issues and projection issues

[William Booker]

Dear Thomas,

The full problem I'm considering has additional terms, including adding another unknown for the pressure in the system.
Would we able to extend the forward elimination and backwards reconstruction to 3 variables?

I'll write up the full problem and send the compressible stratified code along with,

Thanks

Wil




________________________________
From: firedrake-bounces at imperial.ac.uk <firedrake-bounces at imperial.ac.uk> on behalf of Gibson, Thomas <t.gibson15 at imperial.ac.uk>
Sent: 01 February 2017 08:57:36
To: firedrake
Subject: Re: [firedrake] Parallel issues and projection issues


What Lawrence is suggesting looks accurate. Using his definitions of P, Q, R, and S, then the Schur compliment of the block S in that matrix [P Q; R S] gives you a system for u and rho, with LHS operator (P - QS^-1R). The spaces are discontinuous and hence the inverse can be computed cell-local. This is precisely what the Slate package is for in Firedrake.

It's not a problem that some of the matrices are mixed; the problem is that P - QS^-1R is still mixed. For Slate to solve this, we need PyOP2 to write into mixed dats (which it currently can't do). However, if you write out P - QS^-1R in terms of A, B, C, D, E, F and G (as defined by Lawrence) and eliminate, say rho, then you would have a system in terms of one unknown (u). And unless I'm missing something, you can just use the computed u to recover rho. Both of these steps (forward elimination and backwards reconstruction) can be done in Slate.

I would have to see the full finite element problem before I claim Slate can do this. Further edits may be required for handling the `th(f) = theta*dot(f, n)('+') + (1 - theta)*(dot(f, n)('-'))` term. If you're interested in trying this, then we can have a more detailed discussion.

Best regards,
- Thomas
________________________________
From: firedrake-bounces at imperial.ac.uk <firedrake-bounces at imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell at imperial.ac.uk>
Sent: 31 January 2017 18:20:46
To: firedrake
Subject: Re: [firedrake] Parallel issues and projection issues



On 31/01/17 17:44, Lawrence Mitchell wrote:
>
>
> On 31/01/17 12:49, William Booker wrote:
>> Dear David,
>>
>> The choice of theta is meant to be arbitrary in the scheme, and we
>> have had thoughts of looking into an optimal theta for implementations.
>>
>> So I would prefer to keep the asymmetry in the scheme.
>
> So the point is that if theta = 0.5 is not optimal then it will be
> *mesh dependent* (because if theta != 0.5 then the best choice is not
> invariant under renumbering of the mesh).
>
>> I'll implement the conditional form in my DIV definition, and get back
>> to you about the results.
>>
>>
>> Just a quick query, would the facet normal n not have to be assigned
>> to - or + in that expression?
>
> Yes, you want to take an (arbitrarily) restricted side of the expression:
>
> so use:
>
> dot(f, n)('+')
>
>
>> Also, is there a simpler way to implement the projections that I have ?
>> Going forward with a system that effectively doubles the amount of
>> variables I have is a bit cumbersome.
>
> So the physical variables are effectively diagnosed from the
> variational derivatives?
>
> But then the evolution equation says how the physical variable evolves
> according to the derivatives?
>
> Presumably then the elimination of one set, or the other, of these
> variables requires that you invert some operator.  Is this operator
> cell-local?

Having taken a closer look, I think I see what you have:

given trial functions

u, rho, du, drho

and test functions

v, sigma, dv, dsigma

you have a matrix:

[A 0 0 B
 0 C D 0
 E 0 F 0
 0 G 0 H]

If I read the code right where:

A = <u, v>
C = <rho, sigma>
F = <du, dv>
H = <drho, dsigma>

and writing

th(f) = theta*dot(f, n)('+') + (1 - theta)*(dot(f, n)('-'))

B = <grad drho, v> + << jump(drho), th(u) >>
D = <du, grad sigma> + << th(du), jump(sigma) >>

and

E = <u, dv>
G = <rho, dsigma>

But you'd like to solve a system only for u and rho.

So writing

P = [A 0
     0 C]

Q = [0 B
     D 0]

R = [E 0
     0 G]

S = [F 0
     0 H]

You have:

[P Q
 R S]

Block eliminating for the latter row gives (may have got this the
wrong way round)

P - Q S^{-1} R

Which is a system you can solve, and only contains the u, rho
variables, but includes the discrete projections as part of Q and R.

Fortunately, since you're all DG and S is just a mass matrix, S^{-1}
is cell local.

This kind of stuff is exactly what Thomas Gibson has been working on
recently for hybridisation methods.  But it's possible that it will
work for you.

I will let him chime in on whether everything could be done, noting
that P, Q, S and R are all mixed matrices.

I don't see a way other than this, because this elimination is only
weak, not strong (the relationships between u and du, and rho and drho
do not hold quad-point-wise).

Perhaps this is what Colin had thought?

Lawrence

-------------- next part --------------
HTML attachment scrubbed and removed