77
88#include " eos_type.H"
99#include < AMReX.H>
10- #include < iostream>
1110#include < network.H>
1211#include < fundamental_constants.H>
1312#include < extern_parameters.H>
1413#include < cmath>
1514#include < actual_eos_data.H>
16- #include < ostream>
1715
1816const std::string eos_name = " multigamma"
1917
@@ -68,27 +66,20 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE
6866void actual_eos (I input, T& state)
6967{
7068 static_assert (std::is_same_v<I, eos_input_t >, " input must be either eos_input_rt or eos_input_re"
71- const amrex::Real gasconstant = C::n_A * C::k_B;
72- const amrex::Real protonmass = C::m_p;
73- // Special gamma factors
69+ constexpr amrex::Real protonmass = C::m_p;
70+
7471 amrex::Real sum_Abarinv = 0 .0_rt;
75- amrex::Real sum_gammasinv = 0 .0_rt;
72+ amrex::Real sum_ni_fi_over_2 = 0 .0_rt;
7673 amrex::Real rhotot = 0 .0_rt;
7774
7875 for (int n = 0 ; n < NumSpec; ++n) {
79- rhotot += state.xn [n]*spmasses[n];
80- // std::cout << "actual_eos()--> xn[" << n << "]: " << state.xn[n] << ", spmasses[" << n << "]: " << spmasses[n] << std::endl;
81- }
82- // std::cout << "Density in EOS: " << rhotot << std::endl;
83-
84- for (int n = 0 ; n < NumSpec; ++n) {
76+ rhotot += state.xn [n] * spmasses[n];
77+ sum_ni_fi_over_2 += state.xn [n] / (gammas[n] - 1.0 );
8578 sum_Abarinv += state.xn [n];
86- sum_gammasinv += (state.xn [n]*protonmass/rhotot) * (1.0 /(gammas[n]-1.0 ));
8779 }
88- // std::cout << "actual_eos()--> rhotot: " << rhotot << ", state.xn[0]: " << state.xn[0] << ", state.xn[1]: " << state.xn[1] << ", state.xn[2]: " << state.xn[2] << std::endl ;
80+ sum_ni_fi_over_2 /= sum_Abarinv ;
8981 sum_Abarinv *= protonmass/rhotot;
9082 state.mu = 1.0 / sum_Abarinv;
91- sum_gammasinv /= sum_Abarinv;
9283
9384 // -------------------------------------------------------------------------
9485 // For all EOS input modes EXCEPT eos_input_rt, first compute dens
@@ -98,7 +89,6 @@ void actual_eos (I input, T& state)
9889 amrex::Real temp = NAN ;
9990 amrex::Real dens = NAN ;
10091 amrex::Real eint = NAN ;
101- // std::cout << "Eos was called" << std::endl;
10292 switch (input) {
10393
10494 case eos_input_rt:
@@ -107,8 +97,7 @@ void actual_eos (I input, T& state)
10797 // We don't need to do anything here
10898 temp = state.T ;
10999 dens = state.rho ;
110- eint = sum_gammasinv * sum_Abarinv * gasconstant * state.T ;
111- // std::cout << "actual_eos() --> eos_input_rt: temp: " << temp << ", dens: " << dens << ", eint: " << eint << ", sum_gammasinv: " << sum_gammasinv << ", sum_Abarinv: " << sum_Abarinv << std::endl;
100+ eint = sum_ni_fi_over_2 * C::k_B * temp / rhotot;
112101 break ;
113102
114103 case eos_input_re:
@@ -121,9 +110,8 @@ void actual_eos (I input, T& state)
121110
122111 // stop the integration if the internal energy < 0
123112 AMREX_ASSERT (state.e > 0 .);
124- temp = state.e /( sum_gammasinv * gasconstant * sum_Abarinv );
113+ temp = state.e * rhotot / (sum_ni_fi_over_2 * C::k_B );
125114 eint = state.e ;
126- // std::cout << "actual_eos() --> eos_input_re: temp: " << temp << ", dens: " << dens << ", eint: " << eint << ", sum_gammasinv: " << sum_gammasinv << ", sum_Abarinv: " << sum_Abarinv << std::endl;
127115 }
128116
129117 break ;
@@ -151,8 +139,8 @@ void actual_eos (I input, T& state)
151139
152140 // stop the integration if the pressure < 0
153141 AMREX_ASSERT (state.p > 0 .);
154- eint = state.p * sum_gammasinv / dens;
155- temp = eint / (sum_gammasinv * gasconstant * sum_Abarinv );
142+ eint = state.p * sum_ni_fi_over_2 / dens;
143+ temp = eint * rhotot / (sum_ni_fi_over_2 * C::k_B );
156144 }
157145 break ;
158146
@@ -203,18 +191,18 @@ void actual_eos (I input, T& state)
203191 }
204192
205193 if constexpr (has_pressure<T>::value) {
206- amrex::Real pressure = state.rho * eint / sum_gammasinv ;
194+ amrex::Real pressure = state.rho * eint / sum_ni_fi_over_2 ;
207195 state.p = pressure;
208196
209197 state.dpdT = pressure / temp;
210198 state.dpdr = pressure / dens;
211- state.cs = std::sqrt ((1.0 + 1.0 /sum_gammasinv ) * state.p /state.rho );
199+ state.cs = std::sqrt ((1.0 + 1.0 /sum_ni_fi_over_2 ) * state.p /state.rho );
212200 if constexpr (has_G<T>::value) {
213- state.G = 0.5 * (1.0 + (1.0 + 1.0 /sum_gammasinv ));
201+ state.G = 0.5 * (1.0 + (1.0 + 1.0 /sum_ni_fi_over_2 ));
214202 }
215203 }
216204
217- amrex::Real dedT = sum_gammasinv * sum_Abarinv * gasconstant ;
205+ amrex::Real dedT = sum_ni_fi_over_2 * C::k_B / rhotot ;
218206 amrex::Real dedr = 0 .0_rt;
219207 if constexpr (has_energy<T>::value) {
220208 state.dedT = dedT;
0 commit comments