`do_enforce_minimum_density` can divide by zero when rescaling passives

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

Start in Source/hydro/advection_util.cpp at do_enforce_minimum_density, especially the passive-field rescaling around lines 667, 689, and 691. Check the zero and non-positive-density paths, including GPU availability, and verify that floor enforcement cannot produce non-finite passive-field values.