Stokes HDG + average-zero pressure condition

[Guosheng Fu]

Hello again,

Here is my second question:
I am solving the stokes equation with smooth solution:
-laplace u + grad p = f
div u = 0 (*)
u|_bdry= g_bdry
ave-p = 0

C. Lehrenfeld's Hdiv-space is used. The implementation is done without static condensation.
Again, I am using fespace "number" to take care of average pressure zero.
The convergence rate doesn't make sense this time:
order = 1 looks fine,
order = 2 gives me the following rates:
0.96, 9.59e-03 | 0.60, 2.83e-01
3.02, 1.18e-03 | 2.16, 6.32e-02
0.04, 1.15e-03 | 0.75, 3.77e-02 <-----------------this doesn't make any sense
6.01, 1.79e-05 | 3.46, 3.42e-03
3.01, 2.22e-06 | 2.04, 8.29e-04

order =3 gives me:
4.48, 7.07e-04 | 3.48, 4.26e-02
4.12, 4.05e-05 | 3.28, 4.37e-03
-5.41, 1.72e-03 | -3.97, 6.84e-02 <-----------------this doesn't make any sense
5.72, 3.27e-05 | 3.47, 6.17e-03 <-----------------this doesn't make any sense

Then, I did another experiment: Relaxed average-zero pressure condition by changing the divergence-free equation (*) with
div(u) + gamma p = gamma * exact_p
I commented out the average-zero pressure line, and adding back pressure-mass matrix (the lines ended with#FIXME wrong pde)
This time, I got neat convergence rates.
e.g. order = 3 gives me:
4.39, 4.44e-04 | 2.13, 2.40e-02
4.04, 2.71e-05 | 3.35, 2.35e-03
4.03, 1.66e-06 | 2.99, 2.95e-04
4.01, 1.03e-07 | 3.00, 3.69e-05

Finally, I modifed the boundary condition so that pressure is determined:
I changed normal velocity bdry to be the following bdry
(grad u.n).n - p = something
I commented out average-zero and gamma terms, and changed Hdiv space without bdry and added back correct bdry source term on LFI (the lines ended with #FIXME weird bdry)
This time, nice convergence rates again:
e.g. order = 3 gives me:
4.48, 6.05e-04 | 3.04, 3.31e-02
4.04, 3.67e-05 | 3.18, 3.65e-03
4.01, 2.27e-06 | 3.00, 4.56e-04
4.00, 1.42e-07 | 3.00, 5.70e-05

The above tests leads me to believe there is an internal bug for the "number" fespace...
But, I am very confused why the average-zero "number" fespace trick works for the laplace case, but not for the stokes case...

Best,
Guosheng

Here is the sample 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()

order = 2
gamma = 1

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

exact_du1 = CoefficientFunction((-pi*sin(pi*x)*sin(pi*y),
                   pi*cos(pi*x)*cos(pi*y))) # grad(u0)
exact_du2 = CoefficientFunction((-pi*cos(pi*x)*cos(pi*y),
                   pi*sin(pi*x)*sin(pi*y)))   # grad(u1)

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

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

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

# tangential velocity
V4 = FESpace ("vectorfacet", mesh, order=order, dirichlet=[1,2,3,4])

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

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

gradu = CoefficientFunction ( (grad(u),), dims=(2,2) )
gradv = CoefficientFunction ( (grad(v),), dims=(2,2) )

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

def tang(vec):
    return vec - (vec*n)*n

a = BilinearForm (V,flags = { "eliminate_internal" : False })
#a = BilinearForm (V,flags = { "eliminate_internal" : True })
a += SymbolicBFI (InnerProduct(gradu, gradv) - p*div(v)+ div(u)*q )
a += SymbolicBFI ( InnerProduct ( gradu.trans * n,  tang(vhat-v) ), element_boundary=True )
a += SymbolicBFI ( InnerProduct ( gradv.trans * n,  tang(uhat-u) ), element_boundary=True )
a += SymbolicBFI ( 5*(order+1)*(order+1)/h * InnerProduct ( tang(vhat-v),  tang(uhat-u) ), element_boundary=True )
#a += SymbolicBFI (p0*q - p * q0) # averate-zero pressure
#a += SymbolicBFI (gamma*p*q + p0 * q0) # FIXME wrong pde
a += SymbolicBFI (p0 * q0) # FIXME weird bdry

f = LinearForm(V)
f += SymbolicLFI(source*v)
#f += SymbolicLFI(gamma*exact_p*q) # FIXME wrong pde
f += SymbolicLFI((((exact_du1 * n)*n[0]+(exact_du2 * n)*n[1]
          -exact_p)*v*n), BND) # FIXME weird bdry

gfu = GridFunction(V)

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

    gfu.components[0].Set((exact_u*n)*n,BND) # interpolate normal bdry
    gfu.components[2].Set(tang(exact_u),BND) # interpolate tangential 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=-1,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(2)
    rate_p = log(store[i-1][3]/store[i][3])/log(2)
    print("%2.2f, %2.2e | %2.2f, %2.2e" % \
      (rate_u, store[i][2], rate_p, store[i][3]))
    i =  i+1
    sleep(0.5)

[Philip Lederer]

Hey Guosheng!

i tried your example and observed the same problems with ndof=8443. When you draw the pressure you can see that something goes wrong there. At the moment i don't know exactly that's happening, but i don't think it's the numer FESpace because I tried the same with Taylor Hood and everthing worked fine.

[Anonymous]

Hey Guosheng!

Joachim tried the same on his mac and had no problem. We realised that he uses the umfpack solver instead of the pardiso so we think that this leads to the observed problems. You can either regularise by adding

a += SymbolicBFI (1e-8*p0 * q0)

to your Bilinearform or change the signs in your bilinearform to get a symmetric system:

:::python
a = BilinearForm (V,flags = { "eliminate_internal" : False }, symmetric=True")
a += SymbolicBFI (InnerProduct(gradu, gradv) + pdiv(v)+ div(u)q )
a += SymbolicBFI ( InnerProduct ( gradu.trans * n, tang(vhat-v) ), element_boundary=True )
a += SymbolicBFI ( InnerProduct ( gradv.trans * n, tang(uhat-u) ), element_boundary=True )
a += SymbolicBFI ( 5(order+1)(order+1)/h * InnerProduct ( tang(vhat-v), tang(uhat-u) ), element_boundary=True )
a += SymbolicBFI (p0*q + p * q0)

Somehow the problems (at least on my computer) dissaper, maybe also because of the the flag "symmetric=True", as now the symmetric Pardiso is called. Don't forget to change the sign of the pressure where you calculate the pressure error,... it took me a while to realize that ;)

Another possibility is to use a CG solver:

c = Preconditioner(a, type="direct", flags={ "inverse": "pardiso"} )
c.Update()
bvp = BVP(bf=a, lf=f, gf=sol, pre=c, maxsteps=20)
bvp.Do()

...and finally, of course you can also install umfpack :D

all the best,
Philip

[Guosheng Fu]

Hi Philip,

Oh... it's a solver issue...thanks!!
I was trying to install umfpack...mumps...no good luck....I give up
Should have done that in Brenel with you guys... :(

All right, I will symmetrize the system wheneven possible when using pardiso... I use to trust direct solver 100%, not anymore

Here is a funny test:
I used CG solver with pardiso preconditioner

The same problem is tested as my original post, I solve either use static condensation or not.
Both approaches give me same error up to order 3, I noticed the direct approach need a few iterations to converge for order=3 while static condensed approach just need 1 solve.
However, for order = 4 or 5, the direct approach fail to converge, no issue for static condensed approach

It seems pardiso preconditioner + CG solver is a good and safe solver to tell us whether it was doing the "right" thing... :)

Best,
Guosheng

[Philip Lederer]

At least it works now...
For UMFPACK download the package from here:
http://faculty.cse.tamu.edu/davis/suitesparse.html

then do a make (maybe you have to install some additional packages).
After that add
-DUSE_UMFPACK=ON
and
-DCMAKE_PREFIX_PATH= "...path to the suite sparse folder where you used make"
to your cmake flags for ngsolve and compile!

now umfpack should be available like that

inv_a = a.mat.Inverse(V.FreeDofs(), inverse="umfpack")
res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_a * res

(jsut tested it again, indeed, umfpack works!)

[Guosheng Fu]

Good!
I will try once again then.... :)

Since we are on the solver discussion.
What other (cool) solver other than your inhouse ones is, or is planed, to hook into ngsolve?
Is there any plan for addting Trilinos or Petsc in the near future?

[Guosheng Fu]

Philip,

Just an update, I installed umpack as you suggested.
It works, is more stable, but is really slow and memory consuming comparing with pardiso... did you observe similar behavior of the solver?

Best,
Guosheng

[Philip Lederer]

...note to myself: log in before you answer ;)

[schruste]

side remark:
If you update NGSolve now, you will observe that the "error barrier" is decreased a bit so that you obtain your optimal order eoc for one additional level of refinement ;).

[Guosheng Fu]

Why is that the case?
You improved pardiso??

[schruste]

No, something simpler. A technicality on how to compute the gradient of Hdiv-functions sufficiently accurate.