`NavierStokesBase::getForce` (HIT): getForceVerbose min/max loops dereference device memory on the host and read past `Scal`

[WeiqunZhang]

Location: Tutorials/HIT/NS_getForce.cpp:55, Tutorials/HIT/NS_getForce.cpp:106, Tutorials/HIT/NS_getForce.cpp:743
Severity: Low — diagnostic-only path (ns.getForceVerbose=1), but it segfaults GPU builds and reads out of bounds
Category: GPU | Memory-UB
Based on commit bb697bf5 (line numbers refer to that tree).

Problem

PR #204 rewrote the verbose block of Source/NS_getForce.cpp:94-107,182-192 to use
State.min<RunOn::Gpu>(n) / max bounded by nComp(). The HIT override, which replaces that
file in the HIT build, kept the pre-fix hand-rolled loops for Vel, Scal and force:

// Tutorials/HIT/NS_getForce.cpp:55-56,151-157
   const Real* VelDataPtr  = Vel.dataPtr();
   const Real* ScalDataPtr = Scal.dataPtr(scalScomp);
   ...
               for (int n=0; n<NUM_SCALARS; n++) {
                  int cell = ((n*kx+k)*jx+j)*ix+i;
                  Real s = ScalDataPtr[cell];

MacProj::mac_sync_compute (Source/MacProj.cpp:595) passes rhoMF[Smfi], a one-component
FAB (FillPatchIterator rho_fpi(...,Density,1), MacProj.cpp:535), as Scal with
scalScomp=0, so the loop reads NUM_SCALARS-1 components past the FAB. In a GPU build the
FAB data live in device memory and the host dereference faults.

Impact

Tutorials/HIT with ns.getForceVerbose=1: segfault on the first getForce call with
USE_CUDA/HIP; on CPU an out-of-bounds read on every mac-sync call once amr.max_level>0.

Suggested fix

Mirror the base implementation: drop the raw pointers and the hand-rolled loops, and use
FArrayBox::min/max<RunOn::Gpu> bounded by each FAB's own component count.

--- a/Tutorials/HIT/NS_getForce.cpp
+++ b/Tutorials/HIT/NS_getForce.cpp
@@ -53,7 +53,4 @@
 {
 
-   const Real* VelDataPtr  = Vel.dataPtr();
-   const Real* ScalDataPtr = Scal.dataPtr(scalScomp);
-
    const Real  grav     = gravity;
    const int*  f_lo     = force.loVect();
@@ -104,67 +101,13 @@
 #endif
 
-      Vector<Real> velmin(AMREX_SPACEDIM), velmax(AMREX_SPACEDIM);
-      Vector<Real> scalmin(NUM_SCALARS), scalmax(NUM_SCALARS);
-      for (int n=0; n<AMREX_SPACEDIM; n++) {
-          velmin[n]= 1.e234;
-          velmax[n]=-1.e234;
-      }
-      int ix = v_hi[0]-v_lo[0]+1;
-      int jx = v_hi[1]-v_lo[1]+1;
-#if (AMREX_SPACEDIM == 3)
-      int kx = v_hi[2]-v_lo[2]+1;
-      for (int k=0; k<kx; k++) {
-#endif
-         for (int j=0; j<jx; j++) {
-            for (int i=0; i<ix; i++) {
-               for (int n=0; n<AMREX_SPACEDIM; n++) {
-#if (AMREX_SPACEDIM == 3)
-                  int cell = ((n*kx+k)*jx+j)*ix+i;
-#else
-                  int cell = (n*jx+j)*ix+i;
-#endif
-                  Real v = VelDataPtr[cell];
-                  if (v<velmin[n]) velmin[n] = v;
-                  if (v>velmax[n]) velmax[n] = v;
-               }
-            }
-         }
-#if (AMREX_SPACEDIM == 3)
-      }
-#endif
-      for (int n=0; n<AMREX_SPACEDIM; n++)
+      for (int n=0; n<Vel.nComp(); n++)
          amrex::Print() << "Vel  " << n << " min/max "
-                        << velmin[n] << " / " << velmax[n] << '\n';
+                        << Vel.min<RunOn::Gpu>(n) << " / "
+                        << Vel.max<RunOn::Gpu>(n) << '\n';
 
-      for (int n=0; n<NUM_SCALARS; n++) {
-         scalmin[n]= 1.e234;
-         scalmax[n]=-1.e234;
-      }
-      ix = s_hi[0]-s_lo[0]+1;
-      jx = s_hi[1]-s_lo[1]+1;
-#if (AMREX_SPACEDIM == 3)
-      kx = s_hi[2]-s_lo[2]+1;
-      for (int k=0; k<kx; k++) {
-#endif
-         for (int j=0; j<jx; j++) {
-            for (int i=0; i<ix; i++) {
-               for (int n=0; n<NUM_SCALARS; n++) {
-#if (AMREX_SPACEDIM == 3)
-                  int cell = ((n*kx+k)*jx+j)*ix+i;
-#else
-                  int cell = (n*jx+j)*ix+i;
-#endif
-                  Real s = ScalDataPtr[cell];
-                  if (s<scalmin[n]) scalmin[n] = s;
-                  if (s>scalmax[n]) scalmax[n] = s;
-               }
-            }
-         }
-#if (AMREX_SPACEDIM == 3)
-      }
-#endif
-      for (int n=0; n<NUM_SCALARS; n++)
-         amrex::Print() << "Scal " << n << " min/max " << scalmin[n]
-                        << " / " << scalmax[n] << '\n';
+      for (int n=scalScomp; n<Scal.nComp(); n++)
+         amrex::Print() << "Scal " << n << " min/max "
+                        << Scal.min<RunOn::Gpu>(n) << " / "
+                        << Scal.max<RunOn::Gpu>(n) << '\n';
    } //end if(getForceVerbose)
 
@@ -741,38 +684,12 @@
 
    if (ParallelDescriptor::IOProcessor() && getForceVerbose) {
-      Vector<Real> forcemin(ncomp);
-      Vector<Real> forcemax(ncomp);
+      // For scalar-only calls force holds ncomp components starting at 0.
       for (int n=0; n<ncomp; n++) {
-         forcemin[n]= 1.e234;
-         forcemax[n]=-1.e234;
+         const int fcomp = ( scomp<AMREX_SPACEDIM ) ? scomp+n : n;
+         amrex::Print() << "Force " << n+scomp << " min/max "
+                        << force.min<RunOn::Gpu>(fcomp) << " / "
+                        << force.max<RunOn::Gpu>(fcomp) << '\n';
       }
 
-      int ix = f_hi[0]-f_lo[0]+1;
-      int jx = f_hi[1]-f_lo[1]+1;
-#if (AMREX_SPACEDIM == 3)
-      int kx = f_hi[2]-f_lo[2]+1;
-      for (int k=0; k<kx; k++) {
-#endif
-         for (int j=0; j<jx; j++) {
-            for (int i=0; i<ix; i++) {
-               for (int n=0; n<ncomp; n++) {
-#if (AMREX_SPACEDIM == 3)
-                  int cell = ((n*kx+k)*jx+j)*ix+i;
-#else
-                  int cell = (n*jx+j)*ix+i;
-#endif
-                  Real f = force.dataPtr()[cell];
-                  if (f<forcemin[n]) forcemin[n] = f;
-                  if (f>forcemax[n]) forcemax[n] = f;
-               }
-            }
-         }
-#if (AMREX_SPACEDIM == 3)
-      }
-#endif
-      for (int n=0; n<ncomp; n++)
-         amrex::Print() << "Force " << n+scomp << " min/max " << forcemin[n]
-                        << " / " << forcemax[n] << '\n';
-
       amrex::Print() << "NavierStokesBase::getForce(): Leaving..."
                      << '\n' << "---" << '\n';