From b3a5a6cec58de1b84f7836dbebafa339176e1e73 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 13 Jan 2025 10:57:58 -0500 Subject: [PATCH] more EOS work --- Docs/source/eos.rst | 258 ++++------------------------ Docs/source/eos_implementations.rst | 211 +++++++++++++++++++++++ Docs/source/index.rst | 1 + 3 files changed, 247 insertions(+), 223 deletions(-) create mode 100644 Docs/source/eos_implementations.rst diff --git a/Docs/source/eos.rst b/Docs/source/eos.rst index 9d8c69229..0c8778505 100644 --- a/Docs/source/eos.rst +++ b/Docs/source/eos.rst @@ -63,222 +63,17 @@ The *eos_type* passed in is one of .. tip:: The EOS implementations make heavy use of templating to - ``compile-out'' the thermodynamic quantities that are not needed + "compile-out" the thermodynamic quantities that are not needed (depending on the input type). This can greatly increase performance. As such, you should pick the smallest EOS structure (``eos_re_t``, ``eos_rep_t``, ...) that contains the thermodynamic information you need. -Available Equations of State -============================ - -.. index:: eos_t - -The following equations of state are available in Microphysics. - -.. note:: - - Except where noted, each of these EOSs will provide the full - thermodynamic data (including all derivatives) in the ``eos_t`` - type. - - -gamma_law ---------- - -``gamma_law`` represents a gamma law gas, with -equation of state: - -.. math:: p = (\gamma - 1) \rho e. - -:math:`\gamma` is specified by the runtime parameter ``eos.eos_gamma``. For -an ideal gas, this represents the ratio of specific heats. The gas is -assumed to be ideal, with the pressure given by - -.. math:: p = \frac{\rho k T}{\mu m_u} - -where :math:`k` is Boltzmann’s constant and :math:`\mu` is the mean molecular -weight, calculated from the composition, :math:`X_k`. This EOS assumes -the gas is either completely neutral (``eos.assume_neutral = 1``), -giving: - -.. math:: \mu^{-1} = \sum_k \frac{X_k}{A_k} - -or completely ionized (``eos.assume_neutral = 0``), giving: - -.. math:: \mu^{-1} = \sum_k \left ( 1 + Z_k \right ) \frac{X_k}{A_k} - -The entropy comes from the Sackur-Tetrode equation. Because of the -complex way that composition enters into the entropy, the entropy -formulation here is only correct for a :math:`\gamma = 5/3` gas. - - -polytrope ---------- - -``polytrope`` represents a polytropic fluid, with equation of -state: - -.. math:: p = K \rho^\gamma. - -The gas is also assumed to obey the above gamma law relation -connecting the pressure and internal energy. Therefore :math:`\rho` is the -only independent variable; there is no temperature dependence. The -user either selects from a set of predefined options reflecting -physical polytropes (e.g. a non-relativistic, fully degenerate -electron gas) or inputs their own values for :math:`K` and :math:`\gamma` -via ``eos.polytrope_K`` and ``eos.polytrope_gamma``. - -The runtime parameter ``eos.polytrope_type`` selects the pre-defined -polytropic relations. The options are: - -- ``eos.polytrope_type = 1``: sets :math:`\gamma = 5/3` and - - .. math:: K = \left ( \frac{3}{\pi} \right)^{2/3} \frac{h^2}{20 m_e m_p^{5/3}} \frac{1}{\mu_e^{5/3}} - - where :math:`mu_e` is the mean molecular weight per electron, specified via ``eos.polytrope_mu_e`` - - This is the form appropriate for a non-relativistic - fully-degenerate electron gas. - -- ``eos.polytrope_type = 2``: sets :math:`\gamma = 4/3` and - - .. math:: K = \left ( \frac{3}{\pi} \right)^{1/3} \frac{hc}{8 m_p^{4/3}} \frac{1}{\mu_e^{4/3}} - - This is the form appropriate for a relativistic fully-degenerate - electron gas. - -ztwd ----- - -``ztwd`` is the zero-temperature degenerate electron equation -of state of Chandrasekhar (1935), which is designed to describe -white dward material. The pressure satisfies the equation: - -.. math:: p(x) = A \left( x(2x^2-3)(x^2 + 1)^{1/2} + 3\, \text{sinh}^{-1}(x) \right), - -with :math:`A = \pi m_e^4 c^5 / (3 h^3)`. Here :math:`x` is a dimensionless -measure of the Fermi momentum, with :math:`\rho = B x^3` and :math:`B = 8\pi \mu_e -m_p m_e^3 c^3 / (3h^3)`, where :math:`\mu_e` is the mean molecular weight -per electron and :math:`h` is the Planck constant. - -The enthalpy was worked out by Hachisu (1986): - -.. math:: h(x) = \frac{8A}{B}\left(x^2 + 1\right)^{1/2}. - -(note the unfortunate notation here, but this :math:`h` is specific -enthalpy). The specific internal energy satisfies the standard -relationship to the specific enthalpy: - -.. math:: e = h - p / \rho. - -Since the pressure-density relationship does not admit a closed-form -solution for the density in terms of the pressure, if we call the EOS -with pressure as a primary input then we do Newton-Raphson iteration -to find the density that matches this pressure. - -multigamma ----------- - -``multigamma`` is an ideal gas equation of state where each -species can have a different value of :math:`\gamma`. This mainly affects -how the internal energy is constructed as each species, represented -with a mass fraction :math:`X_k` will have its contribution to the total -specific internal energy take the form of :math:`e = p/\rho/(\gamma_k - 1)`. -The main thermodynamic quantities take the form: - -.. math:: - - \begin{aligned} - p &= \frac{\rho k T}{m_u} \sum_k \frac{X_k}{A_k} \\ - e &= \frac{k T}{m_u} \sum_k \frac{1}{\gamma_k - 1} \frac{X_k}{A_k} \\ - h &= \frac{k T}{m_u} \sum_k \frac{\gamma_k}{\gamma_k - 1} \frac{X_k}{A_k}\end{aligned} - -We recognize that the usual astrophysical :math:`\bar{A}^{-1} = \sum_k -X_k/A_k`, but now we have two other sums that involve different -:math:`\gamma_k` weightings. - -The specific heats are constructed as usual, - -.. math:: - - \begin{aligned} - c_v &= \left . \frac{\partial e}{\partial T} \right |_\rho = - \frac{k}{m_u} \sum_k \frac{1}{\gamma_k - 1} \frac{X_k}{A_k} \\ - c_p &= \left . \frac{\partial h}{\partial T} \right |_p = - \frac{k}{m_u} \sum_k \frac{\gamma_k}{\gamma_k - 1} \frac{X_k}{A_k}\end{aligned} - -and it can be seen that the specific gas constant, :math:`R \equiv c_p - -c_v` is independent of the :math:`\gamma_i`, and is simply :math:`R = -k/m_u\bar{A}` giving the usual relation that :math:`p = R\rho T`. -Furthermore, we can show - -.. math:: - - \Gamma_1 \equiv \left . \frac{\partial \log p}{\partial \log \rho} \right |_s = - \left ( \sum_k \frac{\gamma_k}{\gamma_k - 1} \frac{X_k}{A_k} \right ) \bigg / - \left ( \sum_k \frac{1}{\gamma_k - 1} \frac{X_k}{A_k} \right ) = - \frac{c_p}{c_v} \equiv \gamma_\mathrm{effective} - -and :math:`p = \rho e (\gamma_\mathrm{effective} - 1)`. - -This equation of state takes several runtime parameters that can set -the :math:`\gamma_i` for a specific species. The parameters are: - -- ``eos.eos_gamma_default``: the default :math:`\gamma` to apply for all - species - -- ``eos.species_X_name`` and ``eos.species_X_gamma``: set the - :math:`\gamma_i` for the species whose name is given as - ``eos.species_X_name`` to the value provided by ``eos.species_X_gamma``. - Here, ``X`` can be one of the letters: ``a``, ``b``, or - ``c``, allowing us to specify custom :math:`\gamma_i` for up to three - different species. - -helmholtz ---------- - -``helmholtz`` contains a general, publicly available stellar -equation of state based on the Helmholtz free energy, with -contributions from ions, radiation, and electron degeneracy, as -described in :cite:`timmes:1999`, :cite:`timmes:2000`, :cite:`flash`. - -.. note:: - - Our implementation of the ``helmholtz`` EOS has been modified - extensively from the original Fortran source. It has been - made threadsafe and makes heavy use of C++ templating to optimize - the evaluation of thermodynamic quantities. - -The ``helmholtz`` EOS has the ability to perform a Newton-Raphson -iteration so that if we call the EOS with, e.g., density and energy, -and iterate over temperature until we find the temperature -that matches this density–energy combination. If we cannot find an -appropriate temperature, we will reset it to ``small_temp``, which -needs to be set in the equation of state wrapper module in the code -calling this. - -.. index:: eos.use_eos_coulomb, eos.eos_input_is_constant, eos.eos_ttol, eos.eos_dtol, eos.prad_limiter_rho_c, eos.prad_limiter_delta_rho - -The following runtime parameters affect the EOS: - -* ``eos.use_eos_coulomb`` - -* ``eos.eos_input_is_constant`` : when inverting the EOS for find the - density and/or temperature that match the inputs, there is a choice - of whether to update the inputs to match the final density / - temperature, respecting thermodynamic consistency. If - ``eos_input_is_constant=1`` is set (the default), then we leave the - input thermodynamic quantities unchanged, respecting energy - conservation. - -* ``eos.eos_ttol``, ``eos.eos_dtol`` - -* ``eos.prad_limiter_rho_c``, ``eos.prad_limiter_delta_rho`` +.. tip:: + You can also pass a ``burn_t`` struct into the EOS, although this + will give you access to a much smaller range of data. -We thank Frank Timmes for permitting us to modify his code and -publicly release it in this repository. Composition =========== @@ -295,6 +90,11 @@ Auxiliary Composition .. index:: USE_AUX_THERMO +.. note:: + + The main use-case for the auxiliary composition is when using a reaction + network together with the tabulated NSE state. + With ``USE_AUX_THERMO=TRUE``, we interpret the composition from the auxiliary variables. For ``eos_t eos_state``, the auxiliary variables are @@ -351,14 +151,19 @@ extends the non-"extra" variants with these additional fields. The composition derivatives can be used via the ``composition_derivatives()`` function in ``eos_composition.H`` -to compute :math:`\partial p/\partial X_k |_{\rho, T, X_j}`, :math:`\partial e/\partial X_k |_{\rho, T, X_j}`, and :math:`\partial h/\partial X_k |_{\rho, T, X_j}`. +to compute :math:`\partial p/\partial X_k |_{\rho, T, X_j}`, :math:`\partial e/\partial X_k |_{\rho, T, X_j}`, and :math:`\partial h/\partial X_k |_{\rho, T, X_j}`. This has the interface: + +.. code:: c++ + + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + eos_xderivs_t composition_derivatives (const T& state) + Initialization and Cutoff Values ================================ -Input Validation -================ The EOS will make sure that the inputs are within an acceptable range, (e.g., ``small_temp`` :math:`< T <` ``maxT``). If they are not, then it @@ -368,19 +173,26 @@ If you are calling the EOS with ``eos_input_re``, and if :math:`e < 10^{-200}`, then it calls the EOS with ``eos_input_rt`` with T = max ( ``small_temp``, T ). -User’s are encourage to do their own validation of inputs before calling -the EOS. +.. note:: + + User’s are encourage to do their own validation of inputs before calling + the EOS. EOS Structure ============= -Each EOS should have two main routines by which it interfaces to the -rest of CASTRO. At the beginning of the simulation, -``actual_eos_init`` will perform any initialization steps and save -EOS variables (mainly ``smallt``, the temperature floor, and -``smalld``, the density floor). These variables are stored in the -main EOS module of the code calling these routines. This would be the -appropriate time for, say, loading an interpolation table into memory. +Each EOS should have two main routines through which it interfaces to the +rest of Microphysics. + +* ``actual_eos_init()`` : At the beginning of the simulation, + ``actual_eos_init`` will perform any initialization steps and save + EOS variables (mainly ``smallt``, the temperature floor, and + ``smalld``, the density floor). These variables are stored in the + main EOS module of the code calling these routines. + + This is also where an EOS with tables would read in the tables + and initialize the memory they are stored in. -The main evaluation routine is called ``actual_eos``. It should -accept an eos_input and an ``eos_t`` state; see Section :ref:`data_structures`. +* ``actual_eos()`` : this is the main evaluation routine. It should + accept an ``eos_input_t`` specifying the thermodynamic inputs and a + struct (like ``eos_t``) that stores the thermodynamic quantities. diff --git a/Docs/source/eos_implementations.rst b/Docs/source/eos_implementations.rst new file mode 100644 index 000000000..ef99156b9 --- /dev/null +++ b/Docs/source/eos_implementations.rst @@ -0,0 +1,211 @@ +**************************** +Available Equations of State +**************************** + +.. index:: eos_t + +The following equations of state are available in Microphysics. + +.. note:: + + Except where noted, each of these EOSs will provide the full + thermodynamic data (including all derivatives) in the ``eos_t`` + type. + + +gamma_law +========= + +``gamma_law`` represents a gamma law gas, with +equation of state: + +.. math:: p = (\gamma - 1) \rho e. + +:math:`\gamma` is specified by the runtime parameter ``eos.eos_gamma``. For +an ideal gas, this represents the ratio of specific heats. The gas is +assumed to be ideal, with the pressure given by + +.. math:: p = \frac{\rho k T}{\mu m_u} + +where :math:`k` is Boltzmann’s constant and :math:`\mu` is the mean molecular +weight, calculated from the composition, :math:`X_k`. This EOS assumes +the gas is either completely neutral (``eos.assume_neutral = 1``), +giving: + +.. math:: \mu^{-1} = \sum_k \frac{X_k}{A_k} + +or completely ionized (``eos.assume_neutral = 0``), giving: + +.. math:: \mu^{-1} = \sum_k \left ( 1 + Z_k \right ) \frac{X_k}{A_k} + +The entropy comes from the Sackur-Tetrode equation. Because of the +complex way that composition enters into the entropy, the entropy +formulation here is only correct for a :math:`\gamma = 5/3` gas. + + +polytrope +========= + +``polytrope`` represents a polytropic fluid, with equation of +state: + +.. math:: p = K \rho^\gamma. + +The gas is also assumed to obey the above gamma law relation +connecting the pressure and internal energy. Therefore :math:`\rho` is the +only independent variable; there is no temperature dependence. The +user either selects from a set of predefined options reflecting +physical polytropes (e.g. a non-relativistic, fully degenerate +electron gas) or inputs their own values for :math:`K` and :math:`\gamma` +via ``eos.polytrope_K`` and ``eos.polytrope_gamma``. + +The runtime parameter ``eos.polytrope_type`` selects the pre-defined +polytropic relations. The options are: + +- ``eos.polytrope_type = 1``: sets :math:`\gamma = 5/3` and + + .. math:: K = \left ( \frac{3}{\pi} \right)^{2/3} \frac{h^2}{20 m_e m_p^{5/3}} \frac{1}{\mu_e^{5/3}} + + where :math:`mu_e` is the mean molecular weight per electron, specified via ``eos.polytrope_mu_e`` + + This is the form appropriate for a non-relativistic + fully-degenerate electron gas. + +- ``eos.polytrope_type = 2``: sets :math:`\gamma = 4/3` and + + .. math:: K = \left ( \frac{3}{\pi} \right)^{1/3} \frac{hc}{8 m_p^{4/3}} \frac{1}{\mu_e^{4/3}} + + This is the form appropriate for a relativistic fully-degenerate + electron gas. + +ztwd +==== + +``ztwd`` is the zero-temperature degenerate electron equation +of state of Chandrasekhar (1935), which is designed to describe +white dward material. The pressure satisfies the equation: + +.. math:: p(x) = A \left( x(2x^2-3)(x^2 + 1)^{1/2} + 3\, \text{sinh}^{-1}(x) \right), + +with :math:`A = \pi m_e^4 c^5 / (3 h^3)`. Here :math:`x` is a dimensionless +measure of the Fermi momentum, with :math:`\rho = B x^3` and :math:`B = 8\pi \mu_e +m_p m_e^3 c^3 / (3h^3)`, where :math:`\mu_e` is the mean molecular weight +per electron and :math:`h` is the Planck constant. + +The enthalpy was worked out by Hachisu (1986): + +.. math:: h(x) = \frac{8A}{B}\left(x^2 + 1\right)^{1/2}. + +(note the unfortunate notation here, but this :math:`h` is specific +enthalpy). The specific internal energy satisfies the standard +relationship to the specific enthalpy: + +.. math:: e = h - p / \rho. + +Since the pressure-density relationship does not admit a closed-form +solution for the density in terms of the pressure, if we call the EOS +with pressure as a primary input then we do Newton-Raphson iteration +to find the density that matches this pressure. + +multigamma +========== + +``multigamma`` is an ideal gas equation of state where each +species can have a different value of :math:`\gamma`. This mainly affects +how the internal energy is constructed as each species, represented +with a mass fraction :math:`X_k` will have its contribution to the total +specific internal energy take the form of :math:`e = p/\rho/(\gamma_k - 1)`. +The main thermodynamic quantities take the form: + +.. math:: + + \begin{aligned} + p &= \frac{\rho k T}{m_u} \sum_k \frac{X_k}{A_k} \\ + e &= \frac{k T}{m_u} \sum_k \frac{1}{\gamma_k - 1} \frac{X_k}{A_k} \\ + h &= \frac{k T}{m_u} \sum_k \frac{\gamma_k}{\gamma_k - 1} \frac{X_k}{A_k}\end{aligned} + +We recognize that the usual astrophysical :math:`\bar{A}^{-1} = \sum_k +X_k/A_k`, but now we have two other sums that involve different +:math:`\gamma_k` weightings. + +The specific heats are constructed as usual, + +.. math:: + + \begin{aligned} + c_v &= \left . \frac{\partial e}{\partial T} \right |_\rho = + \frac{k}{m_u} \sum_k \frac{1}{\gamma_k - 1} \frac{X_k}{A_k} \\ + c_p &= \left . \frac{\partial h}{\partial T} \right |_p = + \frac{k}{m_u} \sum_k \frac{\gamma_k}{\gamma_k - 1} \frac{X_k}{A_k}\end{aligned} + +and it can be seen that the specific gas constant, :math:`R \equiv c_p - +c_v` is independent of the :math:`\gamma_i`, and is simply :math:`R = +k/m_u\bar{A}` giving the usual relation that :math:`p = R\rho T`. +Furthermore, we can show + +.. math:: + + \Gamma_1 \equiv \left . \frac{\partial \log p}{\partial \log \rho} \right |_s = + \left ( \sum_k \frac{\gamma_k}{\gamma_k - 1} \frac{X_k}{A_k} \right ) \bigg / + \left ( \sum_k \frac{1}{\gamma_k - 1} \frac{X_k}{A_k} \right ) = + \frac{c_p}{c_v} \equiv \gamma_\mathrm{effective} + +and :math:`p = \rho e (\gamma_\mathrm{effective} - 1)`. + +This equation of state takes several runtime parameters that can set +the :math:`\gamma_i` for a specific species. The parameters are: + +- ``eos.eos_gamma_default``: the default :math:`\gamma` to apply for all + species + +- ``eos.species_X_name`` and ``eos.species_X_gamma``: set the + :math:`\gamma_i` for the species whose name is given as + ``eos.species_X_name`` to the value provided by ``eos.species_X_gamma``. + Here, ``X`` can be one of the letters: ``a``, ``b``, or + ``c``, allowing us to specify custom :math:`\gamma_i` for up to three + different species. + +helmholtz +========= + +``helmholtz`` contains a general, publicly available stellar +equation of state based on the Helmholtz free energy, with +contributions from ions, radiation, and electron degeneracy, as +described in :cite:`timmes:1999`, :cite:`timmes:2000`, :cite:`flash`. + +.. note:: + + Our implementation of the ``helmholtz`` EOS has been modified + extensively from the original Fortran source. It has been + made threadsafe and makes heavy use of C++ templating to optimize + the evaluation of thermodynamic quantities. + +The ``helmholtz`` EOS has the ability to perform a Newton-Raphson +iteration so that if we call the EOS with, e.g., density and energy, +and iterate over temperature until we find the temperature +that matches this density–energy combination. If we cannot find an +appropriate temperature, we will reset it to ``small_temp``, which +needs to be set in the equation of state wrapper module in the code +calling this. + +.. index:: eos.use_eos_coulomb, eos.eos_input_is_constant, eos.eos_ttol, eos.eos_dtol, eos.prad_limiter_rho_c, eos.prad_limiter_delta_rho + +The following runtime parameters affect the EOS: + +* ``eos.use_eos_coulomb`` + +* ``eos.eos_input_is_constant`` : when inverting the EOS for find the + density and/or temperature that match the inputs, there is a choice + of whether to update the inputs to match the final density / + temperature, respecting thermodynamic consistency. If + ``eos_input_is_constant=1`` is set (the default), then we leave the + input thermodynamic quantities unchanged, respecting energy + conservation. + +* ``eos.eos_ttol``, ``eos.eos_dtol`` + +* ``eos.prad_limiter_rho_c``, ``eos.prad_limiter_delta_rho`` + + +We thank Frank Timmes for permitting us to modify his code and +publicly release it in this repository. diff --git a/Docs/source/index.rst b/Docs/source/index.rst index b4ca4e7f1..869b7f7ee 100644 --- a/Docs/source/index.rst +++ b/Docs/source/index.rst @@ -49,6 +49,7 @@ of state routines. :hidden: eos + eos_implementations transport .. toctree::