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