ComputeStrain_S/T high-order z-derivatives miss small-denominator protection

[asalmgren]

Summary

In stretched/terrain strain kernels, high-order z-derivative coefficients repeatedly use:

• idz0 = 1.0 / dz0
• c3 = 2.0 / (f - f2) with f = dz1/dz0 + 2

without any explicit lower bound checks on dz0 or (f - f2).

Why this may be a robustness issue

• In highly distorted or pathological vertical meshes, very small dz0 or near-singular coefficient denominators can produce Inf/Na
  N and destabilize diffusion.
• Even if typical production meshes are safe, explicit guards improve failure behavior and diagnostics.

Example locations

• Source/Diffusion/ERF_ComputeStrain_S.cpp:323-328,349-354,375-380,401-406
• Source/Diffusion/ERF_ComputeStrain_T.cpp:469-473,505-509,533-537,570-574,599-603,633-637

Suggested hardening patch (representative)

diff --git a/Source/Diffusion/ERF_ComputeStrain_S.cpp b/Source/Diffusion/ERF_ComputeStrain_S.cpp
@@
-            Real idz0 = 1.0 / dz0;
-            Real f    = (dz1 / dz0) + 2.0;
-            Real f2   = f*f;
-            Real c3   = 2.0 / (f - f2);
-            Real c2   = -f2*c3;
-            Real c1   = -(1.0-f2)*c3;
-
-            Real du_dz = (c1 * u(i,j,k-1) + c2 * u(i,j,k) + c3 * u(i,j,k+1))*idz0;
+            constexpr Real tiny = 1.0e-12;
+            Real du_dz;
+            if (amrex::Math::abs(dz0) <= tiny) {
+                du_dz = (u(i,j,k) - u(i,j,k-1)) / amrex::max(dz_ptr[k], tiny);
+            } else {
+                Real idz0 = 1.0 / dz0;
+                Real f    = (dz1 / dz0) + 2.0;
+                Real denom = f - f*f;
+                if (amrex::Math::abs(denom) <= tiny) {
+                    du_dz = (u(i,j,k) - u(i,j,k-1)) * idz0;
+                } else {
+                    Real c3 = 2.0 / denom;
+                    Real c2 = -(f*f) * c3;
+                    Real c1 = -(1.0 - f*f) * c3;
+                    du_dz = (c1*u(i,j,k-1) + c2*u(i,j,k) + c3*u(i,j,k+1)) * idz0;
+                }
+            }