mixed laplacian + static condensation

[Guosheng Fu]

Hey,

I have a question on static condensation for mixed laplacian.
(My goal is to implement the hybrid hdiv-Stokes element of Cockburn/Sayas for Brinkman's equation)

Here I am solving a mixed-laplacian equation with neumann bdry.
I used the "number" fe space to hack the average-pressure zero condition suggested by C. Lehrenfeld.

I would like to assemble the system with static condensation with the flag
eliminate_internal = True
in the bilinear form.
My condensed global system would have three type of dofs:
velocity (only on edges) + pressure (1 per element) + "number"

I don't know how to make this happen for either the Hdiv space or the L2 space...
I checked the source code a bit, and found a todo list in l2hofespace.cpp... :<

Best,
Guosheng

Here is the code without static condensation, correct order of convergence is observed.

PS: how about adding Raviart-Thomas Hdiv element (and also the Hcurl version...) on tet/trig into the source code?? :>

from ngsolve import *
from netgen.geom2d import unit_square
from numpy import pi
from time import sleep

mesh = Mesh(unit_square.GenerateMesh(maxh=1.6))
mesh.Refine()
mesh.Refine()

order = 2

exact_u = CoefficientFunction((-2*pi*cos(2*pi*x)*sin(2*pi*y),
                   -2*pi*sin(2*pi*x)*cos(2*pi*y)))

exact_p = CoefficientFunction(sin(2*pi*x)*sin(2*pi*y))
source = 8*pi*pi*exact_p

# BDM for velocity
V2 = HDiv(mesh, order=order, dirichlet=[1,2,3,4])
# Pk-1 for pressure
if (order < 1): #RT0-P0
  V3 = L2(mesh, order= 0)
else : #BDMk-Pk-1
  V3 = L2(mesh, order= order-1)

# FIXME heck a constant
V0 = FESpace("number",mesh)

V = FESpace ([V2,V3, V0])
u, p, p0 = V.TrialFunction()
v, q, q0 = V.TestFunction()

n = specialcf.normal(mesh.dim)

a = BilinearForm (V,flags = { "eliminate_internal" : False })
#a = BilinearForm (V,flags = { "eliminate_internal" : True })
a += SymbolicBFI (u*v - p*div(v)+ div(u)*q ) # mixed-laplace
a += SymbolicBFI (p0*q - p * q0) # averate-zero pressure

f = LinearForm(V)
f += SymbolicLFI(source*q)

gfu = GridFunction(V)

store = []
def SolveBVP():
    V.Update()
    print("ndof:", V.ndof)
    gfu.Update()

    gfu.components[0].Set(exact_u,BND) # interpolate bdry

    a.Assemble()
    f.Assemble()

    # direct implementation
    inv_a = a.mat.Inverse(V.FreeDofs())
    res = f.vec.CreateVector()
    res.data = f.vec - a.mat*gfu.vec
    gfu.vec.data += inv_a * res

    err_u =  sqrt (Integrate ( (gfu.components[0]-exact_u)*(gfu.components[0]-exact_u), mesh))
    err_p =  sqrt (Integrate ( (gfu.components[1]-exact_p)*(gfu.components[1]-exact_p), mesh))
    print (err_u, err_p)
    store.append ( (V.ndof, mesh.ne, err_u,err_p) )
    #Draw (CoefficientFunction(gfu.components[0:2]), mesh, "u-p", 
      #autoscale=False, sd=2, min=0,max=1)
    #input("finish solveBVP")

SolveBVP()

while V.ndof < 50000:  
    mesh.Refine()
    SolveBVP()

i = 1
while i < len(store) :
    rate_u = -log(store[i-1][2]/store[i][2])/log(store[i-1][1]/store[i][1])*2
    rate_p = -log(store[i-1][3]/store[i][3])/log(store[i-1][1]/store[i][1])*2
    print("%2.2f, %2.2e | %2.2f, %2.2e" % \
      (rate_u, store[i-1][2], rate_p, store[i-1][3]))
    i =  i+1
    sleep(0.5)

[schruste]

Dear Guosheng,

You need to add the piecewise constant pressure to the degrees of freedom that are not remove via static condensation. In NGSolve this is done by marking the "coupling type" of a dof. The Hdiv-dofs that are on facets are (at least) INTERFACE_DOFS, interior bubbles are LOCAL_DOFS. As L2-dofs are always LOCAL_DOFS they are eliminated in your problem and you get garbage. You can set the couplingtype of a dof from python. In your case the following would do the trick:

for el in V.components[1].Elements():
 V.SetCouplingType(V.components[0].ndof+el.dofs[0],COUPLING_TYPE.INTERFACE_DOF)
  V.FreeDofs()[V.components[0].ndof+el.dofs[0]] = True
  V.FreeDofs(True)[V.components[0].ndof+el.dofs[0]] = True

This does a loop over the L2-space-Elements and marks the first dof on an element (this is the constant!) as an INTERFACE_DOF.

Next, you have to consider the "free dofs" (V.FreeDofs()). This is a bitarray which marks dofs as "free" if they are not bound to essential conditions. This is important for the linear solvers. Further, in the case that you use static condensation, we have "external free dofs" (you access those with V.FreeDofs(True)). These bitarrays are usually set during the Update-Phase of the FESpace. As we change the structure of the dofs after the Update-Phase we also have to redo some of the Update-work, i.e. we correct the FreeDofs-Flag manually. Note that V.Update() overwrites the couplingtypes s.t. we have to do the previous steps after each V.Update().

To do the correct thing with linear systems you may simply use

    BVP(bf=a,lf=f,gf=gfu,pre=c,maxsteps=1000).Do()

if you have set up a corresponding preconditioner c.

Best,
Christoph

[Guosheng Fu]

Dear Christoph,

Thanks for the quick response.
Ah..ha
Here are two follow up questions:
1, how to choose between iterative solvers in BVP? I guess there is a flag for cg/gmres...couldn't find it

2, another bug solving the linear system manually, I plot the pressure for order = 2, the values near boundary is wrong... don't know why.... i suspect the line on boundary conditions is wrong...
res.data = f.vec - a.mat*gfu.vec

Best,
Guosheng

[schruste]

Dear Guosheng,

Yeah, I forgot to solve the interior bubble equations.... and I confused the order of boundary value homogenization and static condensation. Here is code which should work:

    inv_a = a.mat.Inverse(V.FreeDofs(static_condensation), inverse="pardiso")
    res = f.vec.CreateVector()
    update = f.vec.CreateVector()

    if static_condensation:
        f.vec.data += a.harmonic_extension_trans * f.vec

    #homogenize w.r.t. boundary conditions
    res.data = f.vec - a.mat*gfu.vec
    #add correction (boundary conditions)
    gfu.vec.data += inv_a * res

    if static_condensation:
        gfu.vec.data += a.harmonic_extension * gfu.vec
        gfu.vec.data += a.inner_solve * f.vec

Good luck,
Christoph

[Guosheng Fu]

Good! It works.
the inner solve do not need boundary data...

I guess CG+direct pardiso preconditioner would work for me for a while... :)

[schruste]

I forgot: On how to choose the solvers: So far the python-interface only provides CG.
The only other iterative solver exported to the python-interface so far is the QMRSolver which you could use (or you write your own one in python).
To use a different the QMRSolver you should use the aforementioned approach and replace inv_a with inv_a = QMRSolver(....)

Best,
Christoph