init_ConvectedVortex: meanFlowDir=+/-2 swaps u/v; meanFlowDir=0 leaves velocity unset

[WeiqunZhang]

Severity: medium/low · Category: correctness, physics-numerics · Fix order: 12 of 21 — fix this 12th.

_{Filenames are numbered in reverse fix order: 001 = fix last, 021 = fix first. This file is 010.}

Locations: Source/prob/prob_init.cpp:662, Source/prob/prob_init.cpp:650

Based on commit 9bf664bf (line numbers refer to that tree).

Both defects live in the single switch(IC.meanFlowDir) of init_ConvectedVortex, and both are transcription losses from the Fortran CoVo kernel this routine replaced.

The defect

Source/prob/prob_init.cpp:662 — init_ConvectedVortex cases meanFlowDir=2/-2 assign vel_x=v_vort and vel_y=u_vort, swapping the vortex components while keeping the deltax/deltay formulas, which yields a non-solenoidal, essentially irrotational field (div = Fpsi(dy^2-dx^2)/r^4 != 0) instead of a vortex.

Source/prob/prob_init.cpp:650 — The switch(IC.meanFlowDir) in init_ConvectedVortex has no case 0 and no default, and meanFlowDir defaults to 0 (prob_init.H:41), so if prob.meanFlowDir is not set the velocity is never assigned and the vortex is silently dropped.

Why it matters

F023: probtype=8 with prob.meanFlowDir=2 (mean flow in y): initial velocity has O(forcevort) divergence and quadrupole vorticity; the initial nodal projection immediately deforms the field, so the isentropic-vortex benchmark (analytic solution = translated IC) fails; cases 1 and 3 use the correct (u_vort, v_vort) assignment, showing 2/-2 are copy-paste errors.

F064: probtype=8 without prob.meanFlowDir in the inputs (it is only pp.query'd): velocity stays at the S_new.setVal(0.0) value everywhere, so the run silently starts from rest with no vortex and no warning, unlike init_DoubleShearLayer/init_TaylorGreen which Abort on missing parameters.

Suggested fix

Restore the pre-port kernel: the deleted Fortran CoVo (Exec/eb_CoVo_2d/PROB_2D.F90, commit 673e01c, lines 580-598) had case 2: vel(i,j,1)=u_vort; vel(i,j,2)=meanFlowMag+v_vort plus the sign-flipped -2 twin, so the C++ port simply transcribed those two branches with u_vort/v_vort swapped. Only the mean-flow term should change with direction; the Gaussian vortex is rotationally symmetric.

That Fortran also defaulted meanFlowDir=1 and validated it in probinit, aborting on anything outside {+/-1,+/-2,+/-3}. Decide the contract: if 0 is to mean "vortex at rest", make it an explicit case 0 and still validate on the host, because a bare default: silently accepts typos like meanFlowDir=4. Put that check before the ParallelFor, beside the DoubleShearLayer/TaylorGreen Aborts (lines 360, 523) — amrex::Abort cannot run inside the device lambda. Both changes touch the same switch and should land together; a Tutorials y-direction input would keep case +/-2 covered.

For Source/prob/prob_init.cpp:662 (F023):

--- a/Source/prob/prob_init.cpp
+++ b/Source/prob/prob_init.cpp
@@ -659,13 +659,13 @@
                       vel(i,j,k,2) = w_vort);
          break;
       case 2 :
-         AMREX_D_TERM(vel(i,j,k,0) = v_vort;,
-                      vel(i,j,k,1) = IC.meanFlowMag + u_vort;,
+         AMREX_D_TERM(vel(i,j,k,0) = u_vort;,
+                      vel(i,j,k,1) = IC.meanFlowMag + v_vort;,
                       vel(i,j,k,2) = w_vort);
          break;
       case -2 :
-         AMREX_D_TERM(vel(i,j,k,0) = v_vort;,
-                      vel(i,j,k,1) = -IC.meanFlowMag + u_vort;,
+         AMREX_D_TERM(vel(i,j,k,0) = u_vort;,
+                      vel(i,j,k,1) = -IC.meanFlowMag + v_vort;,
                       vel(i,j,k,2) = w_vort);
          break;
       case 3 :

_{Un-swaps the vortex components in cases 2/-2 so each velocity component gets its own vortex term, matching the correct twins (cases 1/-1/3/-3): only the mean-flow term moves to the y-component. Restores a divergence-free vortex plus uniform y translation. No API dependency.}

For Source/prob/prob_init.cpp:650 (F064):

--- a/Source/prob/prob_init.cpp
+++ b/Source/prob/prob_init.cpp
@@ -628,6 +628,9 @@
 {
   const auto domlo = amrex::lbound(domain);
 
+  if ( IC.meanFlowDir < -3 || IC.meanFlowDir > 3 )
+    amrex::Abort("\n    init_ConvectedVortex: prob.meanFlowDir must be 0 (no mean flow) or +/-1, +/-2, +/-3\n    in the inputs file.");
+
   amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
   {
     AMREX_D_TERM(Real x = problo[0] + (i - domlo.x + 0.5)*dx[0];,
@@ -648,6 +651,11 @@
     // Fill Velocity
     //
     switch(IC.meanFlowDir) {
+      case 0  :
+         AMREX_D_TERM(vel(i,j,k,0) = u_vort;,
+                      vel(i,j,k,1) = v_vort;,
+                      vel(i,j,k,2) = w_vort);
+         break;
       case 1  :
          AMREX_D_TERM(vel(i,j,k,0) = IC.meanFlowMag + u_vort;,
                       vel(i,j,k,1) = v_vort;,

_{Makes the default meanFlowDir=0 mean "pure vortex, no mean flow" (consistent with meanFlowMag default 0) and adds a host-side Abort for out-of-range values, in the style of init_DoubleShearLayer's direction check. Maintainer must decide whether 0 should instead be an error demanding an explicit prob.meanFlowDir; if so, change the test to include 0 and drop case 0.}

_{Diff(s) are against 9bf664bf, written from the current source and verified only with git apply --check — never compiled, never run, never applied to the tree. Treat them as precise intent, not tested patches.}

Verification evidence

F023 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Lines 662/667: case 2: vel(i,j,k,0) = v_vort;, vel(i,j,k,1) = IC.meanFlowMag + u_vort; (same swap in case -2), while cases 1/-1/3/-3 use (u_vort, v_vort). With u_vort=-Fdy/r^2psi, v_vort=Fdx/r^2psi, the swapped field has div = Fpsi(dy^2-dx^2)/r^4 != 0 and quadrupole radial flow, not a vortex. Not a rotation (deltas unrotated, no sign flip). Reachable via prob.meanFlowDir=2 (pp.query, line 72); tutorials use 1 so it is latent but supported. git log -S shows only the original commit e98bbe9 touched this; no fix landed.

F064 — confirmed (one verifier lens)

Lens 1 (refutation attempt): prob_init.H:41 int meanFlowDir = 0;; prob_init.cpp:72 only pp.query("meanFlowDir", IC.meanFlowDir); switch at 650-681 has cases 1,-1,2,-2,3,-3, no case 0, no default. Velocity's only prior write is S_new.setVal(0.0) (line 83), so with meanFlowDir unset velocity stays zero while density/tracers are filled — silent start from rest. Comparison holds: init_DoubleShearLayer aborts on missing prob.direction (line 360) and init_TaylorGreen on missing prob.velocity_factor (line 523).

---

_{Based on commit 9bf664bf, which is also the tree the audit verified against. From an automated audit of Source/, Tutorials/ and Util/. Audit finding ids: F023, F064. Reviewer unit(s): EB+Force+Prob. Nothing here was compiled or run — the failure scenarios are code reasoning, so the reaching configuration above is the cheapest way to confirm or refute it.}