[firedrake] Parallel issues and projection issues

[Gibson, Thomas]

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