store mass scalars in HydroState

quokka-astro · Jul 31, 2023 · 3d1d8c6 · 3d1d8c6
1 parent 5157159
commit 3d1d8c6
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 2 deletions.
diff --git a/src/HydroState.hpp b/src/HydroState.hpp
@@ -15,8 +15,9 @@ template <int N> struct HydroState {
 	double E;		      // total energy density
 	double Eint;		      // internal energy density
 	std::array<double, N> scalar; // passive scalars
+	std::array<double, N> massScalar; // mass scalars
 };
 
 } // namespace quokka
 
-#endif // HYDROSTATE_HPP_
+#endif // HYDROSTATE_HPP_
diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp
@@ -937,12 +937,17 @@ void HydroSystem<problem_t>::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu
 		sR.E = E_R;
 		sR.Eint = Eint_R;
 
-		// The remaining components are passive scalars, so just copy them from
+		// The remaining components are mass scalars and passive scalars, so just copy them from
 		// x1LeftState and x1RightState into the (left, right) state vectors U_L and
 		// U_R
 		for (int n = 0; n < nscalars_; ++n) {
 			sL.scalar[n] = x1LeftState(i, j, k, scalar0_index + n);
 			sR.scalar[n] = x1RightState(i, j, k, scalar0_index + n);
+			// also store mass scalars separately
+			if (n < nmscalars_) {
+				sL.massScalar[n] = x1LeftState(i, j, k, scalar0_index + n);
+				sR.massScalar[n] = x1RightState(i, j, k, scalar0_index + n);
+			}
 		}
 
 		// difference in normal velocity along normal axis