dolfin team mailing list archive

[David Ham <question143195@xxxxxxxxxxxxxxxxxxxxx>]

New question #143195 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/143195


Suppose I'm solving the linear shallow water equations with discontinuous velocity and continuous height:

V = VectorFunctionSpace(mesh, 'DG', 1, dim=2)# Velocity space
 
H = FunctionSpace(mesh, 'CG', 2)             # Height space

W=V*H                                        # Mixed space of both.

(v, q)=TestFunctions(W)
(u, h)=TrialFunctions(W)

M=inner(v,q)*dx+inner(u,h)*dx

# Divergence term.
Ct=-inner(u,grad(q))*dx+inner(avg(u),jump(q,n))*dS
# Pressure gradient operator
C=-g*Depth*adjoint(Ct)                            

# Coriolis term
F=params["f"]*inner(v,as_vector([-u[1],u[0]]))*dx  

If I then apply backward Euler timesteping, I get a system whose block structure is:

/   M_u-F  C    \    /u_new \         / rhs_u\
|                       |  |             |    =   |          |
\   Ct         M_h/  \ h_new /         \ rhs_h/

This system is (a) big  and (b) non-symmetric. However, the standard trick at this stage is to use static condensation to turn this into a discrete wave equation which is much smaller and is SPD:

(-Ct(M_u-F)^(-1)C + M_h) h_new = Ct(M_u-F)^(-1)rhs_u + rhs_h

The reason this is affordable is that M_u-F contains no derivatives and therefore no facet integrals. This means that the inversion commutes with assembly so one can directly assemble the inverse matrix by first inverting the local tensor  and then conducting global assembly. (M_u-F)^(-1) has the same sparsity as M_u-F.

So, the question is, how do I get Dolfin to assemble (M_u-F)^(-1) rather than M_u-F?

(The next question will be, can I multiply sparse matrices as above in Dolfin, but one thing at a time!)

Many thanks,

David

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.