Thank you for the iTutorials! They are very helpful to me.
I think there may be an issue with the Navier Stokes tutorial, though. When I copied the code into Jupyter and ran it, I did not see any velocity or pressure output in Netgen. I think the iTutorial code breaks because there is something wrong with the inverse âinv_stokesâ â perhaps the matrix generated from âstokesâ is not invertible as defined. When I added print(u.vec) after u.vec.data += inv_stokes * resI saw many entries with value -nan.
After comparing the code to the code in py_tutorials/navierstokes.py, I noticed one important difference: In input block [5], the definition of the stokes function passed to the integrator is like this stokes = nu*grad(ux)*grad(vx)+nu*grad(uy)*grad(vy) - div_u*q - div_v*p When I changed it to the corresponding line in navierstokes.py stokes = nu*grad(ux)*grad(vx)+nu*grad(uy)*grad(vy) + div_u*q + div_v*p - 1e-10*p*qthe example worked correctly.
But the code in ânavierstokes.pyâ has the opposite signs for the divergence terms and an additional p*q term so it doesnât seem to match the derivation presented in the iTutorial. Iâm not sure what to make of thisâŚ?
Also, it may be helpful to others to change the iTutorial code from Draw(velocity,mesh,"u") to Draw(Norm(velocity),mesh,"velocity", sd=3) as in ânavierstokes.pyâ and have norm(velocity) display by default (instead of pressure) for consistency.
Edit: After making the code change above, I can see the pressure distribution, but it seems that the values may have the wrong sign. My intuition is that the pressure should be higher upstream of the obstruction. I confirmed this at the âfeatflowâ site.
1.) About the regularisation: Indeed the -1e-10 pq is important if you solve with the internal sparsecholesky or with umfpack as the matrix is singular otherwise! If you use the pardiso solver some magic is happening there so you will still get a solutionâŚ
2.) About the mismatch of the sign: When you start with the Stokes equations ( -\nu \Delta u + \nabla p = f) and do integration by parts you end up with a term like - \int_\Omega \div(v) * p. So actually the version with the minus in front of the divergence is the proper one. The positive sign here is just to get a nice looking positive matrix, so only for inner peace
Hi Thanks for explaining all this! I confirmed that by just adding the uv term, the matrix is non-singular and the solution (for both velocity and pressure) appears to be correct. I need to think about why subtracting a small multiple of uv has the desired effect on the matrix and its physical meaning in this equationâŚ
As far as the signs of the divergence terms, it seems that the price for inner peace is that the signs of the pressure values are more or less reversed from what they should be!
I hope you have added â1e-8 * p * qâ and not uv? If you have added only a uv term then you still end up with a singular matrix?
The problem is, that the lower right block (so the part of the pressure space) of the matrix is a big zero block, and thus you canât invert your matrix when using the sparsecholesky solver. When you use the umpfack solver, which would be able to handle such a matrix, you still have the problem that you might have not fixed your pressure, thus you have no unique solution (= singular matrix)!
Think of a simple stokes problem with hom. bc. Then your pressure would be in L^2_0 (mean value = 0) to guarantee a unique pressure. As we have not taken care of that when we defined the pressure space we have to guarantee a unique solution in a different way. Thatâs where the regularisation comes in⌠Note that it would be also possible to define a dirichlet BC for your pressure.
There i have added another lagrange parameter to get a pressure in L^2_0. When you use something like that donât forget that we have an outflow boundary when we talk about NavierStokes, so you might have to ad some boundary termsâŚ
