From 20c4236d511562de102391863658b3454cc49cc1 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 10 Jan 2024 11:15:53 -0500 Subject: [PATCH 1/2] start of a simple 2-T gamma-law EOS --- EOS/gamma_law_2T/Make.package | 1 + EOS/gamma_law_2T/README.md | 26 +++ EOS/gamma_law_2T/_parameters | 3 + EOS/gamma_law_2T/actual_eos.H | 320 ++++++++++++++++++++++++++ networks/general_null/gammalaw_2T.net | 10 + 5 files changed, 360 insertions(+) create mode 100644 EOS/gamma_law_2T/Make.package create mode 100644 EOS/gamma_law_2T/README.md create mode 100644 EOS/gamma_law_2T/_parameters create mode 100644 EOS/gamma_law_2T/actual_eos.H create mode 100644 networks/general_null/gammalaw_2T.net diff --git a/EOS/gamma_law_2T/Make.package b/EOS/gamma_law_2T/Make.package new file mode 100644 index 0000000000..885d445eb8 --- /dev/null +++ b/EOS/gamma_law_2T/Make.package @@ -0,0 +1 @@ +CEXE_headers += actual_eos.H diff --git a/EOS/gamma_law_2T/README.md b/EOS/gamma_law_2T/README.md new file mode 100644 index 0000000000..c422cbee5f --- /dev/null +++ b/EOS/gamma_law_2T/README.md @@ -0,0 +1,26 @@ +# gamma_law_2T + +This is an equation of state for a 2-temperature fluid. We assume +that there is a composition that has neutrals and ions, with +effectively the same mass. We call these collectively the "heavies". +The electrons are assumed to have zero mass but still contribute to +the thermodynamics. + +We assume that there is a single gamma for all components of the EOS. + +At the moment, we assume only a single neutral (this has Z = 0) and a +single ion (this has Z > 0). But this can be generalized as needed. + +This also requires 2 pieces of auxiliary data: + + * the specific internal energy of all the "heavies" + * the specific internal energy of the electrons. + +It is suggested that you use the `general_null` network with the +`gammalaw_2T.net` inputs file. + +There is really no single temperature, so inputs like `eos_input_re` +that want to return a unique temperature will instead get an effective +average temperature, such that calling the EOS with this average +temperature will give the same input e. + diff --git a/EOS/gamma_law_2T/_parameters b/EOS/gamma_law_2T/_parameters new file mode 100644 index 0000000000..fd20f05eb6 --- /dev/null +++ b/EOS/gamma_law_2T/_parameters @@ -0,0 +1,3 @@ +@namespace: eos + +eos_gamma real 5.d0/3.d0 diff --git a/EOS/gamma_law_2T/actual_eos.H b/EOS/gamma_law_2T/actual_eos.H new file mode 100644 index 0000000000..68d7345563 --- /dev/null +++ b/EOS/gamma_law_2T/actual_eos.H @@ -0,0 +1,320 @@ +#ifndef ACTUAL_EOS_H +#define ACTUAL_EOS_H + +#include +#include +#include +#include + +using namespace eos_rp; + +// This is a constant gamma equation of state, using an ideal gas. +// +// The gas may either be completely ionized or completely neutral. +// +// The ratio of specific heats (gamma) is allowed to vary. NOTE: the +// expression for entropy is only valid for an ideal MONATOMIC gas +// (gamma = 5/3). + +const std::string eos_name = "gamma_law_2T"; + +inline +void actual_eos_init() { + + // make sure we have the needed aux data + + // ensure that there are only 2 species, one with Z = 0 + + // ensure that the species have the same atomic mass + + // constant ratio of specific heats + if (eos_gamma <= 0.0) { + amrex::Error("gamma_const cannot be < 0"); + } + +} + + + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +bool is_input_valid (I input) +{ + static_assert(std::is_same::value, "input must be an eos_input_t"); + + bool valid = true; + + if (input == eos_input_th) { + valid = false; + } else if (input == eos_input_tp) { + valid = false; + } else if (input == eos_input_ps) { + valid = false; + } else if (input == eos_input_ph) { + valid = false; + } else if (input == eos_input_th) { + valid = false; + } + + return valid; +} + + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void actual_eos (I input, T& state) +{ + static_assert(std::is_same::value, "input must be an eos_input_t"); + + // for the moment, we will assume just ions and neutrals + static_assert(NumSpec == 2); + + // get the index of the ions + short i_ion{-1}; + short i_neut{-1}; + if (zion[0] > 0) { + i_ion = 0; + i_neut = 1; + } else { + i_ion = 1; + i_neut = 0; + } + + // Get the mass of a nucleon from m_u. + const Real m_nucleon = C::m_u; + + // The mean molecular weight here is constructed to allow us to + // define an average temperature. This is not otherwise used in the + // thermodynamics. + + if constexpr (has_xn::value) { + state.mu = (zion[i_ion] + 1) * state.xn[i_ion] / aion[i_ion] + + state.xn[i_neut] / aion[i_neut]; + state.mu = 1.0 / state.mu + } + + // we can always get the heavy and electron temperatures, since the + // aux data is their internal energies. Note: this assumes that + // aion[i_ion] == aion[i_neut] + + Real T_h = aion[i_ion] * m_nucleon * (eos_gamma - 1.0_rt) / + ((state.xn[i_ion] + state.xn[i_neut]) * C::k_B) * state.aux[AuxZero::e_h]; + + Real T_e = aion[i_ion] * m_nucleon * (eos_gamma - 1.0_rt) / + (state.xn[i_ion] * zion[i_ion] * C::k_B) * state.aux[AuxZero::e_e]; + + // For all EOS input modes EXCEPT eos_input_rt, first compute dens + // and temp as needed from the inputs. + + switch (input) + { + + case eos_input_rt: + + // dens, temp and xmass are inputs + + // We don't need to do anything here + break; + + case eos_input_rh: + + // dens, enthalpy, and xmass are inputs + + // Solve for the temperature: + // h = e + p/rho = (p/rho)*[1 + 1/(gamma-1)] = (p/rho)*gamma/(gamma-1) + + if constexpr (has_enthalpy::value) { + state.T = (state.h * state.mu * m_nucleon / C::k_B)*(eos_gamma - 1.0)/eos_gamma; + } + + break; + + case eos_input_tp: + + // temp, pres, and xmass are inputs + +#ifndef AMREX_USE_GPU + amrex::Error("EOS: eos_input_tp has not been implemented gamma_law_2T EOS."); +#endif + + break; + + case eos_input_rp: + + // dens, pres, and xmass are inputs + + // Solve for the temperature: + // p = rho k T / (mu m_nucleon) + + if constexpr (has_pressure::value) { + state.T = state.p * state.mu * m_nucleon / (C::k_B * state.rho); + } + + break; + + case eos_input_re: + + // dens, energy, and xmass are inputs + + // Solve for the temperature + // e = k T / [(mu m_nucleon)*(gamma-1)] + + if constexpr (has_energy::value) { + state.T = state.e * state.mu * m_nucleon * (eos_gamma - 1.0) / C::k_B; + } + + break; + + case eos_input_ps: + + // pressure, entropy, and xmass are inputs + +#ifndef AMREX_USE_GPU + amrex::Error("EOS: eos_input_ps has not been implemented gamma_law_2T EOS."); +#endif + + break; + + case eos_input_ph: + + // pressure, enthalpy and xmass are inputs + +#ifndef AMREX_USE_GPU + amrex::Error("EOS: eos_input_ph has not been implemented gamma_law_2T EOS."); +#endif + + break; + + case eos_input_th: + + // temperature, enthalpy and xmass are inputs + + // This system is underconstrained. + + #ifndef AMREX_USE_GPU + amrex::Error("EOS: eos_input_th is not a valid input for the gamma law EOS."); + #endif + + break; + + default: + + #ifndef AMREX_USE_GPU + amrex::Error("EOS: invalid input."); + #endif + + break; + + } + + // Now we have the density and temperature (and mass fractions / + // mu), regardless of the inputs. + + Real Tinv = 1.0 / state.T; + Real rhoinv = 1.0 / state.rho; + + // 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 energy = pressure / (eos_gamma - 1.0) * rhoinv; + if constexpr (has_pressure::value) { + state.p = pressure; + } + if constexpr (has_energy::value) { + state.e = energy; + } + + // enthalpy is h = e + p/rho + if constexpr (has_enthalpy::value) { + state.h = energy + pressure * rhoinv; + } + + // entropy (per gram) of an ideal monoatomic gas (the Sackur-Tetrode equation) + // NOTE: this expression is only valid for gamma = 5/3. + if constexpr (has_entropy::value) { + const Real fac = 1.0 / std::pow(2.0 * M_PI * C::hbar * C::hbar, 1.5); + + state.s = (C::k_B / (state.mu * m_nucleon)) * + (2.5 + std::log((std::pow(state.mu * m_nucleon, 2.5) * rhoinv) * + std::pow(C::k_B * state.T, 1.5) * fac)); + } + + // Compute the thermodynamic derivatives and specific heats + if constexpr (has_pressure::value) { + state.dpdT = state.p * Tinv; + state.dpdr = state.p * rhoinv; + } + if constexpr (has_energy::value) { + state.dedT = state.e * Tinv; + state.dedr = 0.0; + } + if constexpr (has_entropy::value) { + state.dsdT = 1.5 * (C::k_B / (state.mu * m_nucleon)) * Tinv; + state.dsdr = - (C::k_B / (state.mu * m_nucleon)) * rhoinv; + } + if constexpr (has_enthalpy::value) { + state.dhdT = state.dedT + state.dpdT * rhoinv; + state.dhdr = 0.0; + } + + if constexpr (has_xne_xnp::value) { + state.xne = 0.0; + state.xnp = 0.0; + } + if constexpr (has_eta::value) { + state.eta = 0.0; + } + if constexpr (has_pele_ppos::value) { + state.pele = 0.0; + state.ppos = 0.0; + } + + if constexpr (has_energy::value) { + state.cv = state.dedT; + + if constexpr (has_pressure::value) { + state.cp = eos_gamma * state.cv; + + state.gam1 = eos_gamma; + + state.dpdr_e = state.dpdr - state.dpdT * state.dedr * (1.0 / state.dedT); + state.dpde = state.dpdT * (1.0 / state.dedT); + + // sound speed + state.cs = std::sqrt(eos_gamma * state.p * rhoinv); + if constexpr (has_G::value) { + state.G = 0.5 * (1.0 + eos_gamma); + } + } + } + + if constexpr (has_dpdA::value) { + state.dpdA = - state.p * (1.0 / state.abar); + } + if constexpr (has_dedA::value) { + state.dedA = - state.e * (1.0 / state.abar); + } + + if (eos_assume_neutral) { + if constexpr (has_dpdZ::value) { + state.dpdZ = 0.0; + } + if constexpr (has_dedZ::value) { + state.dedZ = 0.0; + } + } else { + if constexpr (has_dpdZ::value) { + state.dpdZ = state.p * (1.0 / (1.0 + state.zbar)); + } + if constexpr (has_dedZ::value) { + state.dedZ = state.e * (1.0/(1.0 + state.zbar)); + } + } +} + +inline +void actual_eos_finalize() { + +} + +#endif diff --git a/networks/general_null/gammalaw_2T.net b/networks/general_null/gammalaw_2T.net new file mode 100644 index 0000000000..fbac39c699 --- /dev/null +++ b/networks/general_null/gammalaw_2T.net @@ -0,0 +1,10 @@ +# this is the null description for a network containing H only +# + +# name short-name aion zion + neutrals N 1.0 0.0 + ions I 1.0 1.0 + +# auxiliary variables +__aux_e_h +__aux_e_e From 7829c8ac0200f821166de5c48e0f9f1ccdf380db Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 19 Jan 2024 11:22:33 -0500 Subject: [PATCH 2/2] more updates --- EOS/gamma_law_2T/actual_eos.H | 36 +++++++++-------------------------- 1 file changed, 9 insertions(+), 27 deletions(-) diff --git a/EOS/gamma_law_2T/actual_eos.H b/EOS/gamma_law_2T/actual_eos.H index 68d7345563..37f0c86137 100644 --- a/EOS/gamma_law_2T/actual_eos.H +++ b/EOS/gamma_law_2T/actual_eos.H @@ -35,7 +35,6 @@ void actual_eos_init() { } - template AMREX_GPU_HOST_DEVICE AMREX_INLINE bool is_input_valid (I input) @@ -46,6 +45,8 @@ bool is_input_valid (I input) if (input == eos_input_th) { valid = false; + } else if (input == eos_input_rh) { + valid = false; } else if (input == eos_input_tp) { valid = false; } else if (input == eos_input_ps) { @@ -111,7 +112,7 @@ void actual_eos (I input, T& state) case eos_input_rt: - // dens, temp and xmass are inputs + // dens, xmass are inputs, and we have T_h and T_e from above // We don't need to do anything here break; @@ -120,12 +121,9 @@ void actual_eos (I input, T& state) // dens, enthalpy, and xmass are inputs - // Solve for the temperature: - // h = e + p/rho = (p/rho)*[1 + 1/(gamma-1)] = (p/rho)*gamma/(gamma-1) - - if constexpr (has_enthalpy::value) { - state.T = (state.h * state.mu * m_nucleon / C::k_B)*(eos_gamma - 1.0)/eos_gamma; - } +#ifndef AMREX_USE_GPU + amrex::Error("EOS: eos_input_tp has not been implemented gamma_law_2T EOS."); +#endif break; @@ -142,26 +140,14 @@ void actual_eos (I input, T& state) case eos_input_rp: // dens, pres, and xmass are inputs - - // Solve for the temperature: - // p = rho k T / (mu m_nucleon) - - if constexpr (has_pressure::value) { - state.T = state.p * state.mu * m_nucleon / (C::k_B * state.rho); - } + // we already have the 2 temperatures from above break; case eos_input_re: // dens, energy, and xmass are inputs - - // Solve for the temperature - // e = k T / [(mu m_nucleon)*(gamma-1)] - - if constexpr (has_energy::value) { - state.T = state.e * state.mu * m_nucleon * (eos_gamma - 1.0) / C::k_B; - } + // we already have the 2 temperatures from above break; @@ -207,11 +193,7 @@ void actual_eos (I input, T& state) } - // Now we have the density and temperature (and mass fractions / - // mu), regardless of the inputs. - - Real Tinv = 1.0 / state.T; - Real rhoinv = 1.0 / state.rho; + // Now we have the density and temperature (and mass fractions) // Compute the pressure simply from the ideal gas law, and the // specific internal energy using the gamma-law EOS relation.