Enforce checks and limits on mass scalars (#378)

* return invalid state if mass scalars have negative partial densities

* enforce limits on mass scalars

* use if constexpr

---------

Co-authored-by: Piyush Sharda <[email protected]>

src/hydro_system.hpp

@@ -30,6 +30,9 @@
 #include "radiation_system.hpp"
 #include "valarray.hpp"
 
+// Microphysics headers
+#include "extern_parameters.H"
+
 // this struct is specialized by the user application code
 //
 template <typename problem_t> struct HydroSystem_Traits {
@@ -427,6 +430,18 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem<problem_t>::isStateValid(am
 	const amrex::Real rho = cons(i, j, k, density_index);
 	bool isDensityPositive = (rho > 0.);
 
+	bool isMassScalarPositive = true;
+	if constexpr (nmscalars_ > 0) {
+		amrex::GpuArray<Real, nmscalars_> massScalars_ = RadSystem<problem_t>::ComputeMassScalars(cons, i, j, k);
+
+		for (int idx = 0; idx < nmscalars_; ++idx) {
+			if (massScalars_[idx] < 0.0) {
+				isMassScalarPositive = false;
+				break;  // Exit the loop early if any element is not positive
+			}
+		}
+
+	}
 	// when the dual energy method is used, we *cannot* reset on pressure
 	// failures. on the other hand, we don't need to -- the auxiliary internal
 	// energy is used instead!
@@ -441,7 +456,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem<problem_t>::isStateValid(am
 #endif
 	// return (isDensityPositive && isPressurePositive);
 
-	return isDensityPositive;
+	return isDensityPositive && isMassScalarPositive;
 }
 
 template <typename problem_t>
@@ -726,6 +741,17 @@ void HydroSystem<problem_t>::EnforceLimits(amrex::Real const densityFloor, amrex
 			state[bx](i, j, k, x3Momentum_index) *= rescale_factor;
 		}
 
+		// Enforcing Limits on mass scalars
+		if (nmscalars_ > 0) {
+
+		    for (int idx = 0; idx < nmscalars_; ++idx) {
+		        if (state[bx](i, j, k, scalar0_index + idx) < 0.0) {
+		            state[bx](i, j, k, scalar0_index + idx) = small_x * rho;
+		        }
+		    }
+
+		}
+
 		// Enforcing Limits on temperature estimated from Etot and Ekin
 		if (!HydroSystem<problem_t>::is_eos_isothermal()) {
 			// recompute gas energy (to prevent P < 0)