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..37f0c86137 --- /dev/null +++ b/EOS/gamma_law_2T/actual_eos.H @@ -0,0 +1,302 @@ +#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_rh) { + 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, xmass are inputs, and we have T_h and T_e from above + + // We don't need to do anything here + break; + + case eos_input_rh: + + // dens, enthalpy, 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_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 + // we already have the 2 temperatures from above + + break; + + case eos_input_re: + + // dens, energy, and xmass are inputs + // we already have the 2 temperatures from above + + 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) + + // 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