Tests: abort on NaN solutions in linear solver tests (#5671)

The tests judge their answers with a max-norm, but that reduction drops
NaNs: std::max(a,b) is (a < b) ? b : a, so a NaN leaves the accumulator
untouched. An all-NaN solution reports a max-norm error of 0, and since
ctest reads only the exit code, a failed solve could pass the suite.

Check the solution with MultiFab::contains_nan after each solve and
abort if it trips. ngrow=0, since ghost cells may hold signaling NaNs.

ABecLaplacian_F is left alone; the Fortran interface has no binding.

Author: WeiqunZhang

Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp

@@ -1293,21 +1293,33 @@ namespace {
             int const jm = lowerNeighbor(j,1);
             int const jp = upperNeighbor(j,1);
             Real const xc = xa(i,j,k,n);
+
+            Real hpz = Real(1.0);
+            Real hmz = Real(1.0);
+            if constexpr (UseEB) {
+                hpz = (this->edgecent[1](i,j,k) == Real(1.0))
+                    ? Real(1.0) : Real(1.0)+Real(2.0)*this->edgecent[1](i,j,k);
+                hmz = (this->edgecent[1](i,j-1,k) == Real(1.0))
+                    ? Real(1.0) : Real(1.0)-Real(2.0)*this->edgecent[1](i,j-1,k);
+            }
+
             Real out;
-            Real scale = Real(1.0);
+            Real scale;
 
             if (r == Real(0.0)) {
                 Real const sigp = this->radialEdgeSigma(i,j,k,r);
                 if constexpr (UseEB) {
                     if (this->levset(i+1,j,k) >= Real(0.0)) {
                         Real const hp = (this->edgecent[0](i,j,k) == Real(1.0))
                             ? Real(1.0) : Real(1.0)+Real(2.0)*this->edgecent[0](i,j,k);
-                        out = -Real(4.0)*sigp*xc/(dr*dr*hp*hp);
-                        scale = hp;
+                        scale = amrex::min(hp,hmz,hpz);
+                        out = -Real(4.0)*sigp*mlebndfdlap_scaled_h2inv(scale,hp)*xc/(dr*dr);
                     } else {
-                        out = Real(4.0)*sigp*(xa(ip,j,k,n)-xc)/(dr*dr);
+                        scale = amrex::min(Real(1.0),hmz,hpz);
+                        out = scale*Real(4.0)*sigp*(xa(ip,j,k,n)-xc)/(dr*dr);
                     }
                 } else {
+                    scale = Real(1.0);
                     out = Real(4.0)*sigp*(xa(ip,j,k,n)-xc)/(dr*dr);
                 }
             } else {
@@ -1319,52 +1331,44 @@ namespace {
                     hm = (this->edgecent[0](i-1,j,k) == Real(1.0))
                         ? Real(1.0) : Real(1.0)-Real(2.0)*this->edgecent[0](i-1,j,k);
                 }
+                scale = amrex::min(hm,hp,hmz,hpz);
 
                 Real const sigp = this->radialEdgeSigma(i,j,k,r);
                 Real const sigm = this->radialEdgeSigma(i-1,j,k,r);
                 Real tmp;
                 if constexpr (UseEB) {
                     tmp = (this->levset(i+1,j,k) < Real(0.0))
-                        ? sigp*(xa(ip,j,k,n)-xc)*(r+Real(0.5)*dr)
-                        : -sigp*xc/hp*(r+Real(0.5)*hp*dr);
+                        ? scale*sigp*(xa(ip,j,k,n)-xc)*(r+Real(0.5)*dr)
+                        : -sigp*mlebndfdlap_scaled_hinv(scale,hp)*xc*(r+Real(0.5)*hp*dr);
                     tmp += (this->levset(i-1,j,k) < Real(0.0))
-                        ? sigm*(xa(im,j,k,n)-xc)*(r-Real(0.5)*dr)
-                        : -sigm*xc/hm*(r-Real(0.5)*hm*dr);
+                        ? scale*sigm*(xa(im,j,k,n)-xc)*(r-Real(0.5)*dr)
+                        : -sigm*mlebndfdlap_scaled_hinv(scale,hm)*xc*(r-Real(0.5)*hm*dr);
                 } else {
                     tmp = sigp*(xa(ip,j,k,n)-xc)*(r+Real(0.5)*dr)
                         + sigm*(xa(im,j,k,n)-xc)*(r-Real(0.5)*dr);
                 }
                 out = tmp*Real(2.0)/((hp+hm)*r*dr*dr);
-                scale = amrex::min(hm,hp);
-            }
-
-            Real hp = Real(1.0);
-            Real hm = Real(1.0);
-            if constexpr (UseEB) {
-                hp = (this->edgecent[1](i,j,k) == Real(1.0))
-                    ? Real(1.0) : Real(1.0)+Real(2.0)*this->edgecent[1](i,j,k);
-                hm = (this->edgecent[1](i,j-1,k) == Real(1.0))
-                    ? Real(1.0) : Real(1.0)-Real(2.0)*this->edgecent[1](i,j-1,k);
             }
 
             Real const sigp = this->axialEdgeSigma(i,j,k,r);
             Real const sigm = this->axialEdgeSigma(i,j-1,k,r);
             Real tmp;
             if constexpr (UseEB) {
                 tmp = (this->levset(i,j+1,k) < Real(0.0))
-                    ? sigp*(xa(i,jp,k,n)-xc) : -sigp*xc/hp;
+                    ? scale*sigp*(xa(i,jp,k,n)-xc)
+                    : -sigp*mlebndfdlap_scaled_hinv(scale,hpz)*xc;
                 tmp += (this->levset(i,j-1,k) < Real(0.0))
-                    ? sigm*(xa(i,jm,k,n)-xc) : -sigm*xc/hm;
+                    ? scale*sigm*(xa(i,jm,k,n)-xc)
+                    : -sigm*mlebndfdlap_scaled_hinv(scale,hmz)*xc;
             } else {
                 tmp = sigp*(xa(i,jp,k,n)-xc) + sigm*(xa(i,jm,k,n)-xc);
             }
-            out += tmp*Real(2.0)/((hp+hm)*dz*dz);
-            scale = amrex::min(scale,hm,hp);
+            out += tmp*Real(2.0)/((hpz+hmz)*dz*dz);
 
             if (r != Real(0.0)) {
-                out -= alpha*xc/(r*r);
+                out -= scale*alpha*xc/(r*r);
             }
-            return out*scale;
+            return out;
         }
 
         Real dr;
@@ -1480,54 +1484,63 @@ namespace {
             }
 
             Real const xc = xa(i,j,k,n);
-            int const im = lowerNeighbor(i,0);
-            int const ip = upperNeighbor(i,0);
             Real const hpx = (edgecent[0](i,j,k) == Real(1.0))
                 ? Real(1.0) : Real(1.0)+Real(2.0)*edgecent[0](i,j,k);
             Real const hmx = (edgecent[0](i-1,j,k) == Real(1.0))
                 ? Real(1.0) : Real(1.0)-Real(2.0)*edgecent[0](i-1,j,k);
+            Real const hpy = (edgecent[1](i,j,k) == Real(1.0))
+                ? Real(1.0) : Real(1.0)+Real(2.0)*edgecent[1](i,j,k);
+            Real const hmy = (edgecent[1](i,j-1,k) == Real(1.0))
+                ? Real(1.0) : Real(1.0)-Real(2.0)*edgecent[1](i,j-1,k);
+#if (AMREX_SPACEDIM == 3)
+            Real const hpz = (edgecent[2](i,j,k) == Real(1.0))
+                ? Real(1.0) : Real(1.0)+Real(2.0)*edgecent[2](i,j,k);
+            Real const hmz = (edgecent[2](i,j,k-1) == Real(1.0))
+                ? Real(1.0) : Real(1.0)-Real(2.0)*edgecent[2](i,j,k-1);
+            Real const scale = amrex::min(hmx,hpx,hmy,hpy,hmz,hpz);
+#else
+            Real const scale = amrex::min(hmx,hpx,hmy,hpy);
+#endif
+
+            int const im = lowerNeighbor(i,0);
+            int const ip = upperNeighbor(i,0);
             Real const sigxp = this->edgeSigmaX(i,j,k);
             Real const sigxm = this->edgeSigmaX(i-1,j,k);
             Real tmp = (levset(i+1,j,k) < Real(0.0))
-                ? sigxp*(xa(ip,j,k,n)-xc) : -sigxp*xc/hpx;
+                ? sigxp*scale*(xa(ip,j,k,n)-xc)
+                : -sigxp*mlebndfdlap_scaled_hinv(scale,hpx)*xc;
             tmp += (levset(i-1,j,k) < Real(0.0))
-                ? sigxm*(xa(im,j,k,n)-xc) : -sigxm*xc/hmx;
+                ? sigxm*scale*(xa(im,j,k,n)-xc)
+                : -sigxm*mlebndfdlap_scaled_hinv(scale,hmx)*xc;
             Real y = beta[0]*tmp*Real(2.0)/(hpx+hmx);
-            Real scale = amrex::min(hmx,hpx);
 
             int const jm = lowerNeighbor(j,1);
             int const jp = upperNeighbor(j,1);
-            Real const hpy = (edgecent[1](i,j,k) == Real(1.0))
-                ? Real(1.0) : Real(1.0)+Real(2.0)*edgecent[1](i,j,k);
-            Real const hmy = (edgecent[1](i,j-1,k) == Real(1.0))
-                ? Real(1.0) : Real(1.0)-Real(2.0)*edgecent[1](i,j-1,k);
             Real const sigyp = this->edgeSigmaY(i,j,k);
             Real const sigym = this->edgeSigmaY(i,j-1,k);
             tmp = (levset(i,j+1,k) < Real(0.0))
-                ? sigyp*(xa(i,jp,k,n)-xc) : -sigyp*xc/hpy;
+                ? sigyp*scale*(xa(i,jp,k,n)-xc)
+                : -sigyp*mlebndfdlap_scaled_hinv(scale,hpy)*xc;
             tmp += (levset(i,j-1,k) < Real(0.0))
-                ? sigym*(xa(i,jm,k,n)-xc) : -sigym*xc/hmy;
+                ? sigym*scale*(xa(i,jm,k,n)-xc)
+                : -sigym*mlebndfdlap_scaled_hinv(scale,hmy)*xc;
             y += beta[1]*tmp*Real(2.0)/(hpy+hmy);
-            scale = amrex::min(scale,hmy,hpy);
 
 #if (AMREX_SPACEDIM == 3)
             int const km = lowerNeighbor(k,2);
             int const kp = upperNeighbor(k,2);
-            Real const hpz = (edgecent[2](i,j,k) == Real(1.0))
-                ? Real(1.0) : Real(1.0)+Real(2.0)*edgecent[2](i,j,k);
-            Real const hmz = (edgecent[2](i,j,k-1) == Real(1.0))
-                ? Real(1.0) : Real(1.0)-Real(2.0)*edgecent[2](i,j,k-1);
             Real const sigzp = this->edgeSigmaZ(i,j,k);
             Real const sigzm = this->edgeSigmaZ(i,j,k-1);
             tmp = (levset(i,j,k+1) < Real(0.0))
-                ? sigzp*(xa(i,j,kp,n)-xc) : -sigzp*xc/hpz;
+                ? sigzp*scale*(xa(i,j,kp,n)-xc)
+                : -sigzp*mlebndfdlap_scaled_hinv(scale,hpz)*xc;
             tmp += (levset(i,j,k-1) < Real(0.0))
-                ? sigzm*(xa(i,j,km,n)-xc) : -sigzm*xc/hmz;
+                ? sigzm*scale*(xa(i,j,km,n)-xc)
+                : -sigzm*mlebndfdlap_scaled_hinv(scale,hmz)*xc;
             y += beta[2]*tmp*Real(2.0)/(hpz+hmz);
-            scale = amrex::min(scale,hmz,hpz);
 #endif
 
-            return y*scale;
+            return y;
         }
 
         Array4<Real const> levset;

Tests/LinearSolvers/ABecLap_SP/MyTest.cpp

@@ -28,6 +28,15 @@ MyTest::solve ()
     } else {
         amrex::Abort("Unknown prob_type");
     }
+
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    for (int ilev = 0; ilev < int(solution.size()); ++ilev) {
+        if (solution[ilev].contains_nan(0, solution[ilev].nComp(), 0)) {
+            amrex::Abort("MyTest::solve: solution contains NaN on level "
+                         + std::to_string(ilev));
+        }
+    }
 }
 
 void

Tests/LinearSolvers/ABecLaplacian_C/MyTest.cpp

@@ -26,6 +26,14 @@ MyTest::solve ()
 #ifdef AMREX_USE_HYPRE
     if (use_mlhypre) {
         solveMLHypre();
+        // A failed solve often returns NaNs.  Check for them explicitly,
+        // because the max-norm checks used by these tests silently drop NaNs.
+        for (int ilev = 0; ilev < int(solution.size()); ++ilev) {
+            if (solution[ilev].contains_nan(0, solution[ilev].nComp(), 0)) {
+                amrex::Abort("MyTest::solve: solution contains NaN on level "
+                             + std::to_string(ilev));
+            }
+        }
         return;
     }
 #endif
@@ -45,6 +53,15 @@ MyTest::solve ()
     } else {
         amrex::Abort("Unknown prob_type");
     }
+
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    for (int ilev = 0; ilev < int(solution.size()); ++ilev) {
+        if (solution[ilev].contains_nan(0, solution[ilev].nComp(), 0)) {
+            amrex::Abort("MyTest::solve: solution contains NaN on level "
+                         + std::to_string(ilev));
+        }
+    }
 }
 
 void

Tests/LinearSolvers/CellEB/MyTest.cpp

@@ -117,6 +117,15 @@ MyTest::solve ()
         mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), tol_rel, tol_abs);
     }
 
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    for (int ilev = 0; ilev < int(phi.size()); ++ilev) {
+        if (phi[ilev].contains_nan(0, phi[ilev].nComp(), 0)) {
+            amrex::Abort("MyTest::solve: solution contains NaN on level "
+                         + std::to_string(ilev));
+        }
+    }
+
     if (verbose) {
         Vector<MultiFab> res(max_level+1);
         for (int ilev = 0; ilev <= max_level; ++ilev) {

Tests/LinearSolvers/CellEB2/MyTest.cpp

@@ -164,6 +164,15 @@ MyTest::solve ()
         }
     }
 
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    for (int ilev = 0; ilev < int(phi.size()); ++ilev) {
+        if (phi[ilev].contains_nan(0, phi[ilev].nComp(), 0)) {
+            amrex::Abort("MyTest::solve: solution contains NaN on level "
+                         + std::to_string(ilev));
+        }
+    }
+
     if (verbose > 0) {
         for (int ilev = 0; ilev <= max_level; ++ilev)
         {

Tests/LinearSolvers/CellOverset/MyTest.cpp

@@ -77,6 +77,12 @@ MyTest::solve ()
 
     // In region with overset mask = 0, phi has valid solution and rhs is zero.
     Real mlmg_err = mlmg.solve({&phi}, {&rhs}, 1.e-11, 0.0);
+
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    if (phi.contains_nan(0, phi.nComp(), 0)) {
+        amrex::Abort("MyTest::solve: solution contains NaN");
+    }
 }
 
 void

Tests/LinearSolvers/CurlCurl/MyTest.cpp

@@ -108,6 +108,16 @@ MyTest::solve ()
         mlmg.solve({&solution}, {&rhs}, tol_rel, tol_abs);
     }
 
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.  This must
+    // happen before the loop below, which overwrites solution with the error.
+    for (int idim = 0; idim < 3; ++idim) {
+        if (solution[idim].contains_nan(0, solution[idim].nComp(), 0)) {
+            amrex::Abort("MyTest::solve: solution contains NaN in direction "
+                         + std::to_string(idim));
+        }
+    }
+
     amrex::Print() << "  Number of cells: " << n_cell << '\n';
     for (int idim = 0; idim < 3; ++idim) {
         MultiFab::Subtract(solution[idim], exact[idim], 0, 0, 1, 0);

Tests/LinearSolvers/EBABecLapOverset/MyTest.cpp

@@ -148,6 +148,15 @@ MyTest::solve (bool use_overset_mask)
         mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), tol_rel, tol_abs);
     }
 
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    for (int ilev = 0; ilev < int(phi.size()); ++ilev) {
+        if (phi[ilev].contains_nan(0, phi[ilev].nComp(), 0)) {
+            amrex::Abort("MyTest::solve: solution contains NaN on level "
+                         + std::to_string(ilev));
+        }
+    }
+
     if (verbose) {
         Vector<MultiFab> res(max_level+1);
         for (int ilev = 0; ilev <= max_level; ++ilev) {

Tests/LinearSolvers/EBConvergenceTest/MyTest.cpp

@@ -71,6 +71,15 @@ MyTest::solve ()
     const Real tol_abs = 0.0;
     mlmg.solve(amrex::GetVecOfPtrs(phi), amrex::GetVecOfConstPtrs(rhs), tol_rel, tol_abs);
 
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    for (int ilev = 0; ilev < int(phi.size()); ++ilev) {
+        if (phi[ilev].contains_nan(0, phi[ilev].nComp(), 0)) {
+            amrex::Abort("MyTest::solve: solution contains NaN on level "
+                         + std::to_string(ilev));
+        }
+    }
+
     for (int ilev = 0; ilev <= max_level; ++ilev) {
         amrex::VisMF::Write(phi[0], "phi-"+std::to_string(ilev));
     }

Tests/LinearSolvers/EBTensor/MyTest.cpp

@@ -122,6 +122,12 @@ MyTest::solve ()
 //    solution.setVal(0.0);
     mlmg.solve({&solution}, {&rhs}, tol_rel, tol_abs);
 
+    // A failed solve often returns NaNs.  Check for them explicitly, because
+    // the max-norm checks used by these tests silently drop NaNs.
+    if (solution.contains_nan(0, solution.nComp(), 0)) {
+        amrex::Abort("MyTest::solve: solution contains NaN");
+    }
+
     MultiFab error(grids, dmap, 1, 0, MFInfo(), *factory);
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
         amrex::Print() << "\n";