- Thank you received: 0
Forum Message
We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.
Notice
The forum is in read only mode.
Unexpected behaviour on mass matrix inversion for Stokes flow
3 years 11 months ago - 3 years 11 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: 3 years 11 months ago by philipp.
- mneunteufel
- Offline
- Premium Member
-
Less
More
- Thank you received: 59
3 years 11 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
3 years 11 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.172 seconds
