Example 18: Stokes problem; regularization term

[ullmannsven]

Hi,

I have a question regarding the regularization chosen in the Stokes problem (Example 18) https://github.com/kinnala/scikit-fem/blob/10.0.2/docs/examples/ex18.py , where you use 1e-6 * C, where C = asm(mass, basis['p']) in the definition of the system matrix K in line 70. I was wondering if there is any physical intuition for the value 1e-6, and is there any physical intuition why exactly the mass matrix is chosen there?

Best and thanks in advance, Sven

[kinnala]

I believe this corresponds to using the penalty method to impose $\int p^2 ~dx = 0$.

As you might be aware, the pressure in the Stokes problem with full Dirichlet BC is unique up to a constant, so the constant has to be fixed in one way or another to avoid singular matrix.

I did not make the example myself, so I cannot comment if there was something behind the value 1e-6. My guess is that it is just a suitable small number to not affect the solution too much. As you might also know, the penalty method is "inconsistent" in the sense that it will always slightly change the solution, so smaller number you use there, the better chances you have not affecting the velocity field.

[kinnala]

Some other techniques I am aware of:

1. Using additional Lagrange multiplier to impose $\int p~dx = 0$
2. Fixing value of $p$ at some node
3. Using some iterative method which does not care about a single zero eigenvalue of the matrix.

[gdmcbain]

I did not make the example myself,

Hello. I made thie example, so I should defend it; I'm sorry it's not as well documented as it could be. Apparently we discussed it on Gitter six years ago, but I don't know whether that's recorded. The regularization or penalization is very standard, so I don't have an immediate reference. It's probably the simplest and most common way to discretize the (Navier–)Stokes equations with finite elements. I'm guessing a good reference would be something like

• Ern, A., Guermond, J.-L. (2004). Theory and Practice of Finite Elements. Springer

but I don't have it to hand.

As to the choice of the value of the penalization constant, the "folklore" that was passed down to me when I was learning FEM was that basically this kind of trick halved ones precision. It would set a level of noisy near-zero divergence that was acceptable and then other physical significance variations in the code would be resolved relative to that noise rather than to the noise-level set by say using double-precision floating-point arithmetic. Given that, the idea was to choose the penalization near the square-root of the machine-epsilon so that the top half of the significant digits were retained for the solution while the rest were given over to whatever was being penalized. This was just told to me like that, in words, without formal mathematical justification; I grumbled and accepted it and moved on until I found methods which allowed me to drop it; see below. I expect Ern & Guermond do have a proper analysis, or

• Elman, H. C., Silvester, D. J.,, Wathen, A. J. (2014). Finite Elements and Fast Iterative Solvers : with Applications in Incompressible Fluid Dynamics, 2nd edn. Oxford University Press

(again, not to hand, sorry).

As @kinnala says:

Using some iterative method which does not care about a single zero eigenvalue of the matrix.

and this in, in my opinion, the way to go. This is demonstrated in Example 32: block diagonally preconditioned Stokes solver

scikit-fem/docs/examples/ex32.py

Line 118 in a785c12

[B, None]], 'csr') # no pressure block required for minres

as it says in the comment there. The idea of using minres

scikit-fem/docs/examples/ex32.py

Line 140 in a785c12

solver=solver_iter_krylov(minres, verbose=True, M=M)),

was taken directly from Elman et al.—very clever, very elegant, works a treat, no more penalization, don't look back.

[ullmannsven]

Thanks a lot for the fast and helpful replies! I will check out example 32, I think this is exactly what I need! Thanks