start of a simple 2-T gamma-law EOS

EOS/gamma_law_2T/Make.package

@@ -0,0 +1 @@
+CEXE_headers += actual_eos.H

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.
+

EOS/gamma_law_2T/_parameters

@@ -0,0 +1,3 @@
+@namespace: eos
+
+eos_gamma            real     5.d0/3.d0

EOS/gamma_law_2T/actual_eos.H

@@ -0,0 +1,302 @@
+#ifndef ACTUAL_EOS_H
+#define ACTUAL_EOS_H
+
+#include <string>
+#include <extern_parameters.H>
+#include <fundamental_constants.H>
+#include <cmath>
+
+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 <typename I>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+bool is_input_valid (I input)
+{
+  static_assert(std::is_same<I, eos_input_t>::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 <typename I, typename T>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void actual_eos (I input, T& state)
+{
+    static_assert(std::is_same<I, eos_input_t>::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<T>::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<T>::value) {
+        state.p = pressure;
+    }
+    if constexpr (has_energy<T>::value) {
+        state.e = energy;
+    }
+
+    // enthalpy is h = e + p/rho
+    if constexpr (has_enthalpy<T>::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<T>::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<T>::value) {
+        state.dpdT = state.p * Tinv;
+        state.dpdr = state.p * rhoinv;
+    }
+    if constexpr (has_energy<T>::value) {
+        state.dedT = state.e * Tinv;
+        state.dedr = 0.0;
+    }
+    if constexpr (has_entropy<T>::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<T>::value) {
+        state.dhdT = state.dedT + state.dpdT * rhoinv;
+        state.dhdr = 0.0;
+    }
+
+    if constexpr (has_xne_xnp<T>::value) {
+        state.xne = 0.0;
+        state.xnp = 0.0;
+    }
+    if constexpr (has_eta<T>::value) {
+        state.eta = 0.0;
+    }
+    if constexpr (has_pele_ppos<T>::value) {
+        state.pele = 0.0;
+        state.ppos = 0.0;
+    }
+
+    if constexpr (has_energy<T>::value) {
+        state.cv = state.dedT;
+
+        if constexpr (has_pressure<T>::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<T>::value) {
+                state.G = 0.5 * (1.0 + eos_gamma);
+            }
+        }
+    }
+
+    if constexpr (has_dpdA<T>::value) {
+        state.dpdA = - state.p * (1.0 / state.abar);
+    }
+    if constexpr (has_dedA<T>::value) {
+        state.dedA = - state.e * (1.0 / state.abar);
+    }
+
+    if (eos_assume_neutral) {
+        if constexpr (has_dpdZ<T>::value) {
+            state.dpdZ = 0.0;
+        }
+        if constexpr (has_dedZ<T>::value) {
+            state.dedZ = 0.0;
+        }
+    } else {
+        if constexpr (has_dpdZ<T>::value) {
+            state.dpdZ = state.p * (1.0 / (1.0 + state.zbar));
+        }
+        if constexpr (has_dedZ<T>::value) {
+            state.dedZ = state.e * (1.0/(1.0 + state.zbar));
+        }
+    }
+}
+
+inline
+void actual_eos_finalize() {
+
+}
+
+#endif

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