Use mass scalars while computing EOS quantities in NSCBC (#610)

### Description
We currently do not pass mass scalars as an input to the
`ComputeSoundSpeed` function in NSCBC.

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [ ] I have added a link to any related issues see (see above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [ ] I have added tests for any new physics that this PR adds to the
code.
- [ ] I have tested this PR on my local computer and all tests pass.
- [ ] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR.

---------

Co-authored-by: Piyush Sharda <psharda@RSAA-043527.local>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

src/NSCBC_inflow.hpp

@@ -39,18 +39,26 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto dQ_dx_inflow_x1_lower(quokka::valarray<
 	const Real w = Q[3];
 	const Real P = Q[4];
 	const Real Eint_aux = Q[5];
+
 	amrex::GpuArray<Real, HydroSystem<problem_t>::nscalars_> s{};
 	for (int i = 0; i < HydroSystem<problem_t>::nscalars_; ++i) {
 		s[i] = Q[6 + i];
 	}
 
-	const Real T = quokka::EOS<problem_t>::ComputeTgasFromEint(rho, quokka::EOS<problem_t>::ComputeEintFromPres(rho, P));
-	const Real Eint_aux_t = quokka::EOS<problem_t>::ComputeEintFromTgas(rho, T_t);
+	static constexpr int nmscalars_ = Physics_Traits<problem_t>::numMassScalars;
+	amrex::GpuArray<Real, nmscalars_> massScalars;
+	for (int n = 0; n < nmscalars_; ++n) {
+		// indexing here follows primVars
+		massScalars[n] = Q[HydroSystem<problem_t>::primScalar0_index + n];
+	}
+
+	const Real T = quokka::EOS<problem_t>::ComputeTgasFromEint(rho, quokka::EOS<problem_t>::ComputeEintFromPres(rho, P, massScalars), massScalars);
+	const Real Eint_aux_t = quokka::EOS<problem_t>::ComputeEintFromTgas(rho, T_t, massScalars);
 
 	const Real du_dx = dQ_dx_data[1];
 	const Real dP_dx = dQ_dx_data[4];
 
-	const Real c = quokka::EOS<problem_t>::ComputeSoundSpeed(rho, P);
+	const Real c = quokka::EOS<problem_t>::ComputeSoundSpeed(rho, P, massScalars);
 	const amrex::Real M = std::clamp(std::sqrt(u * u + v * v + w * w) / c, 0., 1.);
 
 	const Real eta_2 = 2.;
@@ -168,16 +176,23 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setInflowX1LowerLowOrder(const amrex::I
 	// compute primitive vars from valid region
 	quokka::valarray<amrex::Real, N> const Q_i = HydroSystem<problem_t>::ComputePrimVars(consVar, ilo, j, k);
 
+	static constexpr int nmscalars_ = Physics_Traits<problem_t>::numMassScalars;
+	amrex::GpuArray<Real, nmscalars_> massScalars;
+	for (int n = 0; n < nmscalars_; ++n) {
+		// indexing here follows primVars
+		massScalars[n] = Q_i[HydroSystem<problem_t>::primScalar0_index + n];
+	}
+
 	// compute centered ghost values
 	const Real rho = Q_i[0];
-	const Real Eint = quokka::EOS<problem_t>::ComputeEintFromTgas(rho, T_t);
+	const Real Eint = quokka::EOS<problem_t>::ComputeEintFromTgas(rho, T_t, massScalars);
 	quokka::valarray<amrex::Real, N> Q_im1{};
 	Q_im1[0] = rho; // extrapolate density
 	Q_im1[1] = u_t; // prescribe velocity
 	Q_im1[2] = v_t;
 	Q_im1[3] = w_t;
-	Q_im1[4] = quokka::EOS<problem_t>::ComputePressure(rho, Eint); // prescribe temperature
-	Q_im1[5] = Eint;					       // prescribe temperature
+	Q_im1[4] = quokka::EOS<problem_t>::ComputePressure(rho, Eint, massScalars); // prescribe temperature
+	Q_im1[5] = Eint;							    // prescribe temperature
 	for (int i = 0; i < HydroSystem<problem_t>::nscalars_; ++i) {
 		Q_im1[6 + i] = s_t[i]; // prescribe passive scalars
 	}

src/NSCBC_outflow.hpp

@@ -35,6 +35,13 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto dQ_dx_outflow(quokka::valarray<amrex::R
 	const amrex::Real w = Q[3];
 	const amrex::Real P = Q[4];
 
+	static constexpr int nmscalars_ = Physics_Traits<problem_t>::numMassScalars;
+	amrex::GpuArray<Real, nmscalars_> massScalars;
+	for (int n = 0; n < nmscalars_; ++n) {
+		// indexing here follows primVars
+		massScalars[n] = Q[HydroSystem<problem_t>::primScalar0_index + n];
+	}
+
 	// normal derivatives
 	const amrex::Real drho_dx = dQ_dx_data[0];
 	const amrex::Real du_dx = dQ_dx_data[1];
@@ -53,7 +60,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto dQ_dx_outflow(quokka::valarray<amrex::R
 	const amrex::Real dw_dz = dQ_dz_data[3];
 	const amrex::Real dP_dz = dQ_dz_data[4];
 
-	const amrex::Real c = quokka::EOS<problem_t>::ComputeSoundSpeed(rho, P);
+	const amrex::Real c = quokka::EOS<problem_t>::ComputeSoundSpeed(rho, P, massScalars);
 	const amrex::Real M = std::clamp(std::sqrt(u * u + v * v + w * w) / c, 0., 1.);
 	const amrex::Real beta = M;
 	const amrex::Real K = 0.25 * c * (1 - M * M) / L_x; // must be non-zero for well-posedness