define G directly in microp

quokka-astro · Aug 2, 2023 · 4458bb2 · 4458bb2
1 parent 1735c03
commit 4458bb2
Show file tree

Hide file tree

Showing 3 changed files with 8 additions and 14 deletions.
diff --git a/extern/Microphysics b/extern/Microphysics
diff --git a/src/EOS.hpp b/src/EOS.hpp
@@ -243,9 +243,8 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeOtherDeriva
 	amrex::Real dRho_dP = NAN;
 	// compute derivative of pressure w/r/t density at constant entropy, given density and pressure (needed for the fundamental derivative G)
 	amrex::Real dP_dRho_s = NAN;
-	amrex::Real G_gamma = NAN;
-	// compute second derivative of pressure w/r/t density at constant entropy, given density and pressure (needed for the fundamental derivative G)
-	amrex::Real d2P_dRho2_s = NAN;
+	// fundamental derivative
+	amrex::Real G = NAN;
 
 #ifdef PRIMORDIAL_CHEM
 	eos_t chemstate;
@@ -268,8 +267,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeOtherDeriva
 	deint_dP = 1.0 / chemstate.dpde;
 	dRho_dP = 1.0 / (chemstate.dpdr * C::k_B / boltzmann_constant_);
 	dP_dRho_s = chemstate.cs * chemstate.cs;
-	G_gamma = dP_dRho_s * rho / P;
-	d2P_dRho2_s = chemstate.d2pdr2_s;
+	G = chemstate.G;
 
 #else
 	if constexpr (gamma_ != 1.0) {
@@ -282,11 +280,10 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeOtherDeriva
 		deint_dP = 1.0 / estate.dpde;
 		dRho_dP = 1.0 / (estate.dpdr * C::k_B / boltzmann_constant_);
 		dP_dRho_s = estate.cs * estate.cs;
-		G_gamma = dP_dRho_s * rho / P;
-		d2P_dRho2_s = estate.d2pdr2_s;
+		G = estate.G;
 	}
 #endif
-	return std::make_tuple(deint_dRho, deint_dP, dRho_dP, dP_dRho_s, G_gamma, d2P_dRho2_s);
+	return std::make_tuple(deint_dRho, deint_dP, dRho_dP, dP_dRho_s, G);
 }
 
 template <typename problem_t>

diff --git a/src/HLLC.hpp b/src/HLLC.hpp
@@ -44,8 +44,8 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLC(quokka::HydroState<N_scalars, N_ms
 	double S_R = NAN;
 	if (gamma != 1.0) {
 
-		auto [dedr_L, dedp_L, drdp_L, dpdr_s_L, G_gamma_L, d2pdr2_s_L] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sL.rho, sL.P, sL.massScalar);
-		auto [dedr_R, dedp_R, drdp_R, dpdr_s_R, G_gamma_R, d2pdr2_s_R] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sR.rho, sR.P, sR.massScalar);
+		auto [dedr_L, dedp_L, drdp_L, dpdr_s_L, G_L] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sL.rho, sL.P, sL.massScalar);
+		auto [dedr_R, dedp_R, drdp_R, dpdr_s_R, G_R] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sR.rho, sR.P, sR.massScalar);
 
 		// equation A.5a of Kershaw+1998
 		// need specific internal energy here
@@ -57,9 +57,6 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLC(quokka::HydroState<N_scalars, N_ms
 		// equation 4.12 of Kershaw+1998
 		cs_tilde = std::sqrt((1.0 / C_tilde_P) * (H_tilde - 0.5 * vsq_tilde - C_tilde_rho));
 
-		const double G_L = 0.5 * (G_gamma_L + 1.0 + 1.0 + (sL.rho/dpdr_s_L)*d2pdr2_s_L - (sL.rho*dpdr_s_L)/sL.P); // 'fundamental derivative' for ideal gases
-		const double G_R = 0.5 * (G_gamma_R + 1.0 + 1.0 + (sR.rho/dpdr_s_R)*d2pdr2_s_R - (sR.rho*dpdr_s_R)/sR.P);
-
 		const double s_NL = 0.5 * G_L * std::max(dU, 0.); // second-order wavespeed correction
 		const double s_NR = 0.5 * G_R * std::max(dU, 0.);