Re: [firedrake] Nonlinear multi-physics example with fully-coupled Newton solve

22 Sep 2015

      Dear Lawrence,

Just as an explanation: I'm trying to mimic a simple finite volume scheme (two-point flux approximation)

So in my nonlinear form I want to integrate the fluxes

v_\alpha = -kr\alpha/\mu\alpha*K*(\nabla p\alpha - \rho\alpha*g)

over the faces of my elements (which is not yet in my code)

Henrik
...
-----Ursprüngliche Nachricht-----
...
Von: firedrake-bounces@imperial.ac.uk [mailto:firedrake-
...
bounces@imperial.ac.uk] Im Auftrag von Lawrence Mitchell
...
Gesendet: 22 September 2015 11:58
...
An: firedrake@imperial.ac.uk
...
Betreff: Re: [firedrake] Nonlinear multi-physics example with fully-
...
coupled Newton solve
...

...
Dear Henrik,
...

...
I will try and address the bigger comments later (I am in a course
...
today), but a brief note on question (3) is below.
...

...
...
On 22 Sep 2015, at 09:02, Buesing, Henrik <HBuesing@eonerc.rwth-<mailto:HBuesing@eonerc.rwth-aachen.de>
...
aachen.de<mailto:HBuesing@eonerc.rwth-aachen.de>> wrote:
...
...
...
...
Dear Lawrence,
...
...
...
...
I guess, the answer is both. I tried to setup an example for two-phase
...
flow in porous media (see twophase.py). I have attached the system I
...
want to setup (see System.pdf). The primary variables are pw and Sn.
...
...
...
...
pw ... water pressure
...
...
Sn ... gas saturation
...
...
phi ... porosity
...
...
rho_\alpha ... density of phase alpha
...
...
kr\alpha ... relative mobility of phase alpha, e.g. Brooks-Corey or
...
...
just linear as in twophase.py mu\alpha ... viscosity of phase alpha pc
...
...
= 0 ... no capillary pressure
...
...
\vec{g}=(0,0,-g)=(0,0,0) ... no gravity K ... permeability
...
...
...
...
My boundary conditions are Dirichlet for the right boundary and gas
...
inflow at the bottom middle.
...
...
...
...
1) Regarding formulation of the problem:
...
...
...
...
a) I would like to upwind krw, and krn, i.e. always take the upstream
...
value. In the case of no capillarity and no gravity, I only look at
...
pressure gradients, e.g.
...
...
take krw_ij = krw_j if pw_j -pw_i >=0
...
...
= krw_i if pw_j -pw_i < 0
...
...
I had a look at the steady-state advection example, but I didn't
...
manage to transfer it to my problem.
...
...
...
...
b) I would like to include gravity. I tried g =
...
interpolate(Expression("(0.0,0.0,-grav)"),DG2), but got a "Can't add
...
expressions with different shapes, when doingdot(grad(pw)-
...
rhow*g,grad(v)).
...
...
...
...
...
...
2) Regarding solvers:
...
...
...
...
a)  The linear solver always converges within 1 iteration to a
...
residual of 0.0 (bcgs+hypre/boomeramg). That's quite weird. If I use a
...
direct method (preonly+lu) I get NaNs. So, I guess, there is still
...
something wrong.
...
...
b)  When the linear solver takes no step, than the nonlinear solver is
...
happy with that and just converges. This leads to the point, where
...
nothing happens at all.
...
...
...
...
3) Regarding output
...
...
a) Can I deactivate the pyop2 INFO output?
...

...
Set:
...

...
parameters["pyop2_options"]["log_level"] = "WARNING"
...

...

...
Cheers,
...

...
Lawrence