Compute thermodynamic derivatives and implement dens+pres mode in pri…

…mordial chem EOS (AMReX-Astro#1302) Add the eos_input_rp mode for primordial chem EOS. Also, compute the sound speed and stores it in state.cs if pressure is defined.
maxpkatz · Aug 9, 2023 · 9cb1d87 · 9cb1d87
1 parent cec7336
commit 9cb1d87
Show file tree

Hide file tree

Showing 3 changed files with 23 additions and 5 deletions.
diff --git a/EOS/gamma_law/actual_eos.H b/EOS/gamma_law/actual_eos.H
@@ -198,7 +198,7 @@ void actual_eos (I input, T& state)
 
     // Compute the pressure simply from the ideal gas law, and the
     // specific internal energy using the gamma-law EOS relation.
-    Real pressure = state.rho * state.T * (C::k_B / (state.mu * m_nucleon));
+    Real pressure = state.rho * state.T * C::k_B / (state.mu * m_nucleon);
     Real energy = pressure / (eos_gamma - 1.0) * rhoinv;
     if constexpr (has_pressure<T>::value) {
         state.p = pressure;
@@ -265,6 +265,7 @@ void actual_eos (I input, T& state)
 
             // sound speed
             state.cs = std::sqrt(eos_gamma * state.p * rhoinv);
+            state.G = 0.5 * (1.0 + eos_gamma);
         }
     }
 

diff --git a/EOS/primordial_chem/actual_eos.H b/EOS/primordial_chem/actual_eos.H
@@ -214,11 +214,16 @@ void actual_eos (I input, T& state)
 
     case eos_input_rp:
         // dens, pressure, and xmass are inputs
-#ifndef AMREX_USE_GPU
-        amrex::Error("eos_input_rp is not supported");
-#endif
 
-          break;
+        if constexpr (has_pressure<T>::value) {
+            dens = state.rho;
+
+            // stop the integration if the pressure < 0
+            AMREX_ASSERT(state.p > 0.);
+            eint = state.p * sum_gammasinv / dens;
+            temp = eint / (sum_gammasinv * gasconstant * sum_Abarinv);
+        }
+        break;
 
     case eos_input_ps:
         // pressure, entropy, and xmass are inputs
@@ -269,11 +274,21 @@ void actual_eos (I input, T& state)
     if constexpr (has_pressure<T>::value) {
         Real pressure = state.rho * eint / sum_gammasinv;
         state.p = pressure;
+
+        state.dpdT = pressure / temp;
+        state.dpdr = pressure / dens;
+        state.cs = std::sqrt((1.0 + 1.0/sum_gammasinv) * state.p /state.rho);
+        state.G = 0.5 * (1.0 + (1.0 + 1.0/sum_gammasinv));
     }
 
     Real dedT = sum_gammasinv * sum_Abarinv * gasconstant;
+    Real dedr = 0.0_rt;
     if constexpr (has_energy<T>::value) {
         state.dedT = dedT;
+        state.dedr = dedr;
+        if constexpr (has_pressure<T>::value) {
+            state.dpde = state.dpdT / state.dedT;
+        }
     }
 
 }

diff --git a/interfaces/eos_type.H b/interfaces/eos_type.H
@@ -27,6 +27,7 @@ struct eos_t:eos_base_t {
     amrex::Real dsdr{};
     amrex::Real dpde{};
     amrex::Real dpdr_e{};
+    amrex::Real G{}; // fundamental derivative (Thompson 1971), we use equation 2.24 of Menikoff & Plohr 1989
 
     amrex::Real cv{};
     amrex::Real cp{};
@@ -145,6 +146,7 @@ struct chem_eos_t:eos_base_t {
     amrex::Real dpdr_e{};
     amrex::Real dedT{};
     amrex::Real dedr{};
+    amrex::Real G{};
 
     amrex::Real mu{};
     amrex::Real cs{};