IncfloVelFill HIGH J: the Burggraf (probtype 16) lid profile 16x^2(1-x)^2 is assigned to the normal velocity v instead of the tangential velocity u since the PR #122 refactor; benchmark.burggraf runs the wrong problem

[WeiqunZhang]

Severity: high · Category: physics-numerics
Locations: src/prob/prob_bc.H:187-192, src/prob/prob_bc.H:203-213, src/prob/prob_init_fluid.cpp:1149-1154, src/derive/incflo_error.cpp:130-135, test_no_eb_2d/benchmark.burggraf, test_no_eb_3d/benchmark.burggraf
Based on commit 46de3367 (line numbers refer to that tree).

The defect

The regularised lid-driven cavity of Burggraf has exact solution
u = 8 f(x) g'(y), v = -8 f'(x) g(y) with f = x^4-2x^3+x^2, g = y^4-y^2 (this is what
init_burggraf, prob_init_fluid.cpp:1149-1154, and DiffFromExact, incflo_error.cpp:130-135, use). On the
lid y = 1 this gives the tangential velocity u = 16 f(x) = 16(x^4-2x^3+x^2) and normal velocity v = 0.

Until 60881d3 (PR #122, 2024-08-29) the fill functor did exactly that:

if (16 == probtype) {
    if (orig_comp+nc == 0) { // x velocity
        vel(i,j,k,dcomp+nc) = 16.0 * (x*x*x*x - 2.0 * x*x*x + x*x);
    }
}

The refactor in #122 rewrote every face block around a norm_vel variable ("This may modify the normal
velocity for specific problems") and turned the Burggraf case into (prob_bc.H:187-192)

if (16 == probtype) {
    amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
    norm_vel = amrex::Real(16) * (x*x*x*x - amrex::Real(2) * x*x*x + x*x);
}

which is then applied by the generic branch (prob_bc.H:205-213) as
if (orig_comp+nc == dir /* dir = 1 */) vel = norm_vel; else vel = bcv_vel[yhi][orig_comp+nc];.
So the y-hi ghost cells now hold v = 16 x^2 (1-x)^2 (a normal velocity through the lid, pointing out of
the domain) and u = bcv_vel[yhi][0] = 1 (the constant from yhi.velocity = 1. 0. in the inputs).

Why it matters

benchmark.burggraf (2D and 3D no-EB) no longer solves the Burggraf problem: the lid is a uniform-speed
lid (u = 1) with a prescribed normal velocity profile through it. The ghost values feed the Godunov
edge states, the tensor diffusion solve (Dirichlet, get_diffuse_velocity_bc for no_slip_wall) and the
nodal projection (whose BC on a wall is Neumann for phi, so the spurious normal inflow is not even
removed). test_convergence_burggraf compares against the true Burggraf solution and cannot converge.
This is a regression introduced by #122, not a pre-existing modelling choice.

How to reach it

test_no_eb_2d/benchmark.burggraf, test_no_eb_3d/benchmark.burggraf (probtype 16, yhi.type = nsw,
yhi.velocity = 1. 0.), DIM=2 or 3, no EB, CPU or GPU; test_no_eb_*/test_convergence_burggraf.

Suggested fix

Restore the pre-#122 behaviour: leave norm_vel alone for probtype 16 and set the x component of the
high-y ghost cells to the lid profile; the normal component keeps the wall value (0, since
init_bcs zeroes the normal component of a no_slip_wall velocity).

Diff checked with git apply --check (audit/notes/T4-scratch/011.diff):

--- a/src/prob/prob_bc.H
+++ b/src/prob/prob_bc.H
@@ -184,13 +184,6 @@
                 int dir = 1;
                 amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][dir];
 
-                // This may modify the normal velocity for specific problems
-                if (16 == probtype)
-                {
-                    amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
-                    norm_vel = amrex::Real(16) * (x*x*x*x - amrex::Real(2) * x*x*x + x*x);
-                }
-
 #if (AMREX_SPACEDIM == 3)
                 if (1102 == probtype)
                 {
@@ -207,6 +200,11 @@
                     {
                         if ( orig_comp+nc == dir ) {
                             vel(i,j,k,dcomp+nc) = norm_vel;
+                        } else if (16 == probtype && orig_comp+nc == 0) {
+                            // Burggraf lid: tangential velocity u = 16 x^2 (1-x)^2 on the
+                            // high-y wall; the normal velocity keeps the wall value.
+                            amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
+                            vel(i,j,k,dcomp+nc) = amrex::Real(16) * (x*x*x*x - amrex::Real(2) * x*x*x + x*x);
                         } else {
                             vel(i,j,k,dcomp+nc) =  bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][orig_comp+nc];
                         }

[asalmgren]

closed by #265