club derivative functions into one to increase speed

quokka-astro · Jul 31, 2023 · e638244 · e638244
1 parent 240a65d
commit e638244
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 96 deletions.
diff --git a/src/EOS.hpp b/src/EOS.hpp
@@ -9,6 +9,7 @@
 
 #include <cmath>
 #include <optional>
+#include <tuple>
 
 #include "AMReX_Array.H"
 #include "AMReX_GpuQualifiers.H"
@@ -56,16 +57,7 @@ template <typename problem_t> class EOS
 	    -> amrex::Real;
 
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
-	ComputeeintDensDerivative(amrex::Real rho, amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {})
-	    -> amrex::Real;
-
-	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
-	ComputeeintPresDerivative(amrex::Real rho, amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {})
-	    -> amrex::Real;
-
-	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
-	ComputeDensPresDerivative(amrex::Real rho, amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {})
-	    -> amrex::Real;
+	ComputeOtherDerivatives(amrex::Real rho, amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {});
 
 	[[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto
 	ComputePressure(amrex::Real rho, amrex::Real Eint, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars = {}) -> amrex::Real;
@@ -241,88 +233,13 @@ EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re
 
 template <typename problem_t>
 AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto
-EOS<problem_t>::ComputeeintDensDerivative(const amrex::Real rho, const amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
-    -> amrex::Real
+EOS<problem_t>::ComputeOtherDerivatives(const amrex::Real rho, const amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
 {
 	// compute derivative of specific internal energy w/r/t density, given density and pressure
 	amrex::Real deint_dRho = NAN;
-
-#ifdef PRIMORDIAL_CHEM
-	eos_t chemstate;
-	chemstate.rho = rho;
-	chemstate.p = P;
-	// initialize array of number densities
-	for (int ii = 0; ii < NumSpec; ++ii) {
-		chemstate.xn[ii] = -1.0;
-	}
-
-	if (massScalars) {
-		const auto &massArray = *massScalars;
-		for (int nn = 0; nn < nmscalars_; ++nn) {
-			chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho)
-		}
-	}
-
-	eos(eos_input_rp, chemstate);
-	deint_dRho = chemstate.dedr;
-#else
-	if constexpr (gamma_ != 1.0) {
-		chem_eos_t estate;
-		estate.rho = rho;
-		estate.p = P;
-		estate.mu = mean_molecular_weight_ / C::m_u;
-		eos(eos_input_rp, estate);
-		deint_dRho = estate.dedr;
-	}
-#endif
-	return deint_dRho;
-}
-
-template <typename problem_t>
-AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto
-EOS<problem_t>::ComputeeintPresDerivative(const amrex::Real rho, const amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
-    -> amrex::Real
-{
 	// compute derivative of specific internal energy w/r/t density, given density and pressure
 	amrex::Real deint_dP = NAN;
-
-#ifdef PRIMORDIAL_CHEM
-	eos_t chemstate;
-	chemstate.rho = rho;
-	chemstate.p = P;
-	// initialize array of number densities
-	for (int ii = 0; ii < NumSpec; ++ii) {
-		chemstate.xn[ii] = -1.0;
-	}
-
-	if (massScalars) {
-		const auto &massArray = *massScalars;
-		for (int nn = 0; nn < nmscalars_; ++nn) {
-			chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho)
-		}
-	}
-
-	eos(eos_input_rp, chemstate);
-	deint_dP = 1.0 / chemstate.dpde;
-#else
-	if constexpr (gamma_ != 1.0) {
-		chem_eos_t estate;
-		estate.rho = rho;
-		estate.p = P;
-		estate.mu = mean_molecular_weight_ / C::m_u;
-		eos(eos_input_rp, estate);
-		deint_dP = 1.0 / estate.dpde;
-	}
-#endif
-	return deint_dP;
-}
-
-template <typename problem_t>
-AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto
-EOS<problem_t>::ComputeDensPresDerivative(const amrex::Real rho, const amrex::Real P, const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)
-    -> amrex::Real
-{
-	// compute derivative of specific internal energy w/r/t density, given density and pressure
+	// compute derivative of density w/r/t pressure, given density and pressure
 	amrex::Real dRho_dP = NAN;
 
 #ifdef PRIMORDIAL_CHEM
@@ -342,20 +259,26 @@ EOS<problem_t>::ComputeDensPresDerivative(const amrex::Real rho, const amrex::Re
 	}
 
 	eos(eos_input_rp, chemstate);
+	deint_dRho = chemstate.dedr;
+	deint_dP = 1.0 / chemstate.dpde;
 	dRho_dP = 1.0 / (chemstate.dpdr * C::k_B / boltzmann_constant_);
+
 #else
 	if constexpr (gamma_ != 1.0) {
 		chem_eos_t estate;
 		estate.rho = rho;
 		estate.p = P;
 		estate.mu = mean_molecular_weight_ / C::m_u;
 		eos(eos_input_rp, estate);
+		deint_dRho = estate.dedr;
+		deint_dP = 1.0 / estate.dpde;
 		dRho_dP = 1.0 / (estate.dpdr * C::k_B / boltzmann_constant_);
 	}
 #endif
-	return dRho_dP;
+	return std::make_tuple(deint_dRho, deint_dP, dRho_dP);
 }
 
+
 template <typename problem_t>
 AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputePressure(amrex::Real rho, amrex::Real Eint,
 									      const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars)

diff --git a/src/HLLC.hpp b/src/HLLC.hpp
@@ -44,19 +44,13 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLC(quokka::HydroState<N_scalars, N_ms
 
 	if (gamma != 1.0) {
 
-		dedr_L = quokka::EOS<problem_t>::ComputeeintDensDerivative(sL.rho, sL.P, sL.massScalar);
-		dedr_R = quokka::EOS<problem_t>::ComputeeintDensDerivative(sR.rho, sR.P, sR.massScalar);
+		auto [dedr_L, dedp_L, drdp_L] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sL.rho, sL.P, sL.massScalar);
+		auto [dedr_R, dedp_R, drdp_R] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sR.rho, sR.P, sR.massScalar);
 
 		// equation A.5a of Kershaw+1998
 		// need specific internal energy here
 		const double C_tilde_rho = 0.5 * ((sL.Eint / sL.rho) + (sR.Eint / sR.rho) + sL.rho * dedr_L + sR.rho * dedr_R);
 
-		dedp_L = quokka::EOS<problem_t>::ComputeeintPresDerivative(sL.rho, sL.P, sL.massScalar);
-		dedp_R = quokka::EOS<problem_t>::ComputeeintPresDerivative(sR.rho, sR.P, sR.massScalar);
-
-		drdp_L = quokka::EOS<problem_t>::ComputeDensPresDerivative(sL.rho, sL.P, sL.massScalar);
-		drdp_R = quokka::EOS<problem_t>::ComputeDensPresDerivative(sR.rho, sR.P, sR.massScalar);
-
 		// equation A.5b of Kershaw+1998
 		const double C_tilde_P = 0.5 * ((sL.Eint / sL.rho) * drdp_L + (sR.Eint / sR.rho) * drdp_R + sL.rho * dedp_L + sR.rho * dedp_R);