`do_enforce_minimum_density` can divide by zero when rescaling passives

[WeiqunZhang]

Summary

When URHO < small_dens, do_enforce_minimum_density rescales passively-advected conserved fields by small_dens / URHO. If URHO == 0 (or is non-positive),
this causes division by zero or non-finite scaling before the state reset is completed.

Location

• Source/hydro/advection_util.cpp:689
• Source/hydro/advection_util.cpp:691
• Source/hydro/advection_util.cpp:667

Problem Details

Current code path:

if (state_arr(i,j,k,URHO) < small_dens) {
    ...
    for (int ipassive = 0; ipassive < npassive; ipassive++) {
        state_arr(i,j,k,n) *= (small_dens / state_arr(i,j,k,URHO));
    }
    ...
}

A zero-density check exists only inside a verbosity-gated CPU block; it can be skipped (e.g., verbose_warnings == 0) and is unavailable on GPU. The divide
therefore remains reachable.

Impact

• Potential Inf/NaN in passively advected fields during floor enforcement.
• Corrupt composition/auxiliary state fed into EOS (xn, aux) and subsequent updates.
• Hard-to-diagnose instability because the routine is intended to recover from bad states.

Suggested Patch

diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp
--- a/Source/hydro/advection_util.cpp
+++ b/Source/hydro/advection_util.cpp
@@
-      for (int ipassive = 0; ipassive < npassive; ipassive++) {
-          const int n = upassmap(ipassive);
-          state_arr(i,j,k,n) *= (small_dens / state_arr(i,j,k,URHO));
-      }
+      const Real rho_old = state_arr(i,j,k,URHO);
+      if (rho_old > 0.0_rt) {
+          const Real rescale = small_dens / rho_old;
+          for (int ipassive = 0; ipassive < npassive; ipassive++) {
+              const int n = upassmap(ipassive);
+              state_arr(i,j,k,n) *= rescale;
+          }
+      } else {
+          // Non-positive density: avoid generating inf/nan in passives.
+          for (int ipassive = 0; ipassive < npassive; ipassive++) {
+              const int n = upassmap(ipassive);
+              state_arr(i,j,k,n) = 0.0_rt;
+          }
+      }

Prepared by Codex

[zingale]

the issue here is that we don't like to call amrex::Error in a GPU kernel. We could do a reduction on an integer flag and abort on the host if we detect that any zone was bad, or use an atomic, like we do here:
https://github.com/AMReX-Astro/Microphysics/blob/6ff6460dd33b68ed64907213a11845eeb9bb9fb4/unit_test/test_react/main.cpp#L214

But we could also just ignore this -- the code will crash shortly after anyway.