- Thank you received: 0
Unexpected behaviour on mass matrix inversion for Stokes flow
4 years 6 months ago - 4 years 6 months ago #2670
by philipp
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:
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
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:
Code:
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
Last edit: 4 years 6 months ago by philipp.
- mneunteufel
-
- Offline
- Premium Member
-
Less
More
- Thank you received: 59
4 years 6 months ago #2671
by mneunteufel
Replied by mneunteufel on topic Unexpected behaviour on mass matrix inversion for Stokes flow
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-10*p*q to your bilinear form
2) 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
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-10*p*q to your bilinear form
Code:
a += SymbolicBFI(InnerProduct(grad(u),grad(v))+div(u)*q+div(v)*p -1e-10*p*q)
2) 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
Attachments:
The following user(s) said Thank You: philipp
4 years 6 months ago #2673
by philipp
Replied by philipp on topic Unexpected behaviour on mass matrix inversion for Stokes flow
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
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
Time to create page: 0.104 seconds