Spherical benchmarks: mismatch between analytical and computed pressure values

[gthyagi]

Hi @knepley,

I am running spherical benchmarks from the following papers:

1. Paper1
2. Paper2

I successfully reproduced the benchmark results from the first paper, but I'm encountering difficulties with the second. The numerical pressure values deviate significantly from the analytical ones. The primary difference between the benchmarks in these two papers lies in the initial conditions, specifically the magnitude of the velocity on the outer surface and the magnitude of the density difference.

For both benchmarks, I am using identical solver settings.

Analytical (left) and numerical (right) pressure solution:

Analytical (left) and numerical (right) velocity solution:

Here are the scripts I am using:

Paper 1: Benchmark script
Paper 2: Benchmark script

[knepley]

Okay, the first thing to do is see if we are converging. If we converge, but the solution is wrong, then something is wrong in the residual definition. If we are not converging, something is wrong in the solver. So please -snes_monitor -snes_converged_reason -ksp_monitor_true_residual -ksp_converged_reason to the run and post the output.

[gthyagi]

stokes.petsc_options["ksp_monitor_true_residual"] = None
stokes.petsc_options["snes_monitor"] = None

0 SNES Function norm 1.601023304539e+08
    Residual norms for stokes_ solve.
    0 KSP Residual norm 1.601023304539e+08
    Residual norms for stokes_ solve.
    0 KSP unpreconditioned resid norm 1.601023304539e+08 true resid norm 1.601023304539e+08 ||r(i)||/||b|| 1.000000000000e+00
    1 KSP Residual norm 1.370586276036e-07
    1 KSP unpreconditioned resid norm 1.370586276036e-07 true resid norm 1.355832398677e-07 ||r(i)||/||b|| 8.468536309454e-16
  1 SNES Function norm 1.370752142493e-07

stokes.snes.getConvergedReason()
stokes.snes.ksp.getConvergedReason()

3
2

Note: @knepley, only the pressure values do not match. The velocity values look fine.

[knepley]

Okay, the SNES residual is small, so there is a definition problem. Pressure has an arbitrary zero. How are you setting the 0 of pressure, and how does your analytic solution do it?

[gthyagi]

This is what paper says:

"The analytical velocity solution is prescribed on both internal (r = R1) and external (r = R2) boundaries. Given these boundary conditions, the pressure is determined up to a constant. Both codes allow for a surface normalisation (the average pressure along the surface has a prescribed value - in this case zero), and a volume normalisation (the average pressure over the whole volume has a prescribed zero value). Since I have shown here above that both are identically nul for the pressure field, the choice of pressure normalisation does not matter. Also, the boundary conditions preclude the presence of a pure rotational mode of numerical origin (Zhong et al., 2008)."

Currently, I am not setting the 0 of the pressure. I am not sure how to set it. At the moment, Underworld only supports setting boundary conditions on velocity field.

In the case of the Annulus, I generated the density field (i.e., from analytical equation) for the Stokes system and applied velocity conditions at the inner and outer radii. This approach resulted in numerical results consistent with the analytical ones, falling within the expected error range. However, when applying the same method to the spherical case, the outcomes appear anomalous.

[lmoresi]

@gthyagi, You have a very accurate velocity solution, so the first thing we need to understand is whether there is some null space in the pressure that does not appear in the numerical residual but does cause you problems. The pressure null-space is any pressure field that does not excite a flow. As the benchmark paper points out, the surface average values are both indeterminate and they remove them. This implies subtracting a hydrostatic pressure gradient which is also in the null space (hence the name).

The way we set things up, the pressure null space does not cause numerical convergence problems, but you do need to deal with it after the fact in the same way as the analytic solutions do. I am a little bit puzzled by what I am seeing here.

[gthyagi]

Actually the pressure issue is causing problem in L2-norm error convergence in velocity.

We need to deal with the issue before solving Stokes system.

[lmoresi]

Oh, my bad: I misunderstood the issue. I thought the velocity solution was well behaved.

[gthyagi]

@lmoresi Here is the benchmark script. Please go through it to help diagnose the pressure null space issue.
https://github.com/underworldcode/underworld3-documentation/blob/main/Jupyterbook/Notebooks/Examples-Spherical-Stokes/Ex_Stokes_Spherical_Benchmark_Thieulot.py.