From fe202ce4e57db62d048d587e170d4d538af2564b Mon Sep 17 00:00:00 2001
From: psharda <piyush.suneet@gmail.com>
Date: Wed, 2 Aug 2023 00:05:03 +0200
Subject: [PATCH] implement full equation for G

---
 src/EOS.hpp  | 15 +++++++++++----
 src/HLLC.hpp |  8 ++++----
 2 files changed, 15 insertions(+), 8 deletions(-)

diff --git a/src/EOS.hpp b/src/EOS.hpp
index e7f43d8d7..c2e6a012b 100644
--- a/src/EOS.hpp
+++ b/src/EOS.hpp
@@ -241,8 +241,11 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeOtherDeriva
 	amrex::Real deint_dP = NAN;
 	// compute derivative of density w/r/t pressure, given density and pressure
 	amrex::Real dRho_dP = NAN;
-	// compute derivative of pressure w/r/t density, given density and pressure (needed for the fundamental derivative G)
+	// 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;
 
 #ifdef PRIMORDIAL_CHEM
 	eos_t chemstate;
@@ -264,7 +267,9 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeOtherDeriva
 	deint_dRho = chemstate.dedr;
 	deint_dP = 1.0 / chemstate.dpde;
 	dRho_dP = 1.0 / (chemstate.dpdr * C::k_B / boltzmann_constant_);
-	G_gamma = (chemstate.cs * chemstate.cs) * rho / P;
+	dP_dRho_s = chemstate.cs * chemstate.cs;
+	G_gamma = dP_dRho_s * rho / P;
+	d2P_dRho2_s = chemstate.d2pdr2_s;
 
 #else
 	if constexpr (gamma_ != 1.0) {
@@ -276,10 +281,12 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeOtherDeriva
 		deint_dRho = estate.dedr;
 		deint_dP = 1.0 / estate.dpde;
 		dRho_dP = 1.0 / (estate.dpdr * C::k_B / boltzmann_constant_);
-		G_gamma = (estate.cs * estate.cs) * rho / P;
+		dP_dRho_s = estate.cs * estate.cs;
+		G_gamma = dP_dRho_s * rho / P;
+		d2P_dRho2_s = estate.d2pdr2_s;
 	}
 #endif
-	return std::make_tuple(deint_dRho, deint_dP, dRho_dP, G_gamma);
+	return std::make_tuple(deint_dRho, deint_dP, dRho_dP, dP_dRho_s, G_gamma, d2P_dRho2_s);
 }
 
 template <typename problem_t>
diff --git a/src/HLLC.hpp b/src/HLLC.hpp
index 72eba4e4b..cbe5830c7 100644
--- 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, G_gamma_L] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sL.rho, sL.P, sL.massScalar);
-		auto [dedr_R, dedp_R, drdp_R, G_gamma_R] = quokka::EOS<problem_t>::ComputeOtherDerivatives(sR.rho, sR.P, sR.massScalar);
+		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);
 
 		// equation A.5a of Kershaw+1998
 		// need specific internal energy here
@@ -57,8 +57,8 @@ 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.); // 'fundamental derivative' for ideal gases
-		const double G_R = 0.5 * (G_gamma_R + 1.);
+		const double G_L = 0.5 * (G_gamma_L + 1.0 + (std::pow(sL.rho, 2)/(sL.P*G_gamma_L))*d2pdr2_s_L + (sL.rho/(sL.P*G_gamma_L))*dpdr_s_L - (std::pow(sL.rho/sL.P, 2)/G_gamma_L)*std::pow(dpdr_s_L, 2)); // 'fundamental derivative' for ideal gases
+		const double G_R = 0.5 * (G_gamma_R + 1.0 + (std::pow(sR.rho, 2)/(sR.P*G_gamma_R))*d2pdr2_s_R + (sR.rho/(sR.P*G_gamma_R))*dpdr_s_R - (std::pow(sR.rho/sR.P, 2)/G_gamma_R)*std::pow(dpdr_s_R, 2));
 
 		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.);