Navier Stokes iTutorial

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 after I 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 When I changed it to the corresponding line in navierstokes.py the 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 to 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.

Hi,

1.) About the regularisation: Indeed the -1e-10 *p*q 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

I hope that helped.

Hi Thanks for explaining all this! I confirmed that by just adding the u*v 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 u*v 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 u*v? If you have added only a u*v 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.

You can have a look at the files here: ngsolve.org/user-meeting2017/resources/lederer_navierstokes.zip

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...

Yes, sorry for the typo -- I meant to write p*q term not u*v. Thank you for the additional insight!