Unexpected behaviour on mass matrix inversion for Stokes flow

Main Forum Old-Forum - do not post here
1

Greetings to the NGSolve community,

I am currently using ngsolve for some numerical experiments in my PhD thesis. So far, I am quite fond of its possibilities and the nice python interface. Thanks to the programmers for this!
I constantly try to enlargen my knowledge about FEM and numerics in general, but am still at its very beginnings, so I hope this question is not too stupid.

Debugging my code I boiled some problems I had down to the following snippet:

from ngsolve import *
from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
mesh = Mesh( geo.GenerateMesh(maxh=0.05))

V = VectorH1(mesh, order=3, dirichlet="wall|outlet|inlet")
Q = H1(mesh, order=2)
X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()
    
a = BilinearForm(X) 
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)*q+div(v)*p)
a.Assemble() 
    
gfu = GridFunction(X)
gfu.components[0].Set((1,1)) 
gfu.components[1].Set(1)
fdofs = X.FreeDofs()
for i in range(len(gfu.vec)):
    if fdofs[i] == False:
        gfu.vec[i] = 0
print("gfu.vec ", min(gfu.vec), max(gfu.vec))
    
t1 = gfu.vec.CreateVector()
t2 = GridFunction(X)
t1.data = a.mat * gfu.vec
t2.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="umfpack") * t1
print("t2 ", min(t2.vec), max(t2.vec))

As you can see, I am considering the Stokes problem with no-slip boundary conditions everywhere
and check how good the stiffness matrix is inverted. As a rather basic
test I compare the minimum and maximum value of t2.vec and gfu.vec and expect them to be the same up to some (small) deviations. This works perfectly if I leave out the inlet or outlet boundary in the declaration of the Dirichlet boundary. Considering full no-slip boundary, however, the inverse is done quite poorly.
I am confused by this result and up to my knowledge I am not considering an ill-posed problem,
so in theory this should work. Therefore, it seems like I am using ngsolve wrong, but cannot figure out what is missing.
I have to admin that I am not using the latest version of ngsolve. Mine is
NGSolve-6.2.1902-78-g34df28fd if this clarifies anything.

I will be quite grateful for any hint.
Best regards,

Philipp

2

Hi Philipp,

using no-slip condition on the whole boundary leads to the problem that the pressure p is not unique any more as you can add a constant c and p+c still fulfils the equation (Integrate div(v)*p by parts).

Therefore, to make the pressure unique again there are two possibilities:

  1. add a small “regularization” -1e-10pq to your bilinear form
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)*q+div(v)*p -1e-10*p*q)
  1. use a Lagrange multiplier to enforce that p has zero mean value, see attached file

With both strategies your problem with the “poor” inverse should be solved.

Best
Michael

https://ngsolve.org/media/kunena/attachments/889/stokes_mass.py

Attachment: stokes_mass.py

3

Hi Michael,

wow, that was fast. Thank you so much! Of course, now I remember this little caveat,
but it just didn’t came to me. I tried your solution and it works perfectly.
Just one more question: What exactly is the difference between NumberSpace and, e.g., H1? I tried the Lagrange multiplier approach with the latter and it works also.

Thanks again.
Philipp