Skip to content

Commit

Permalink
more EOS work
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale committed Jan 13, 2025
1 parent 9846e10 commit b3a5a6c
Show file tree
Hide file tree
Showing 3 changed files with 247 additions and 223 deletions.
258 changes: 35 additions & 223 deletions Docs/source/eos.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
===========
Expand All @@ -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

Expand Down Expand Up @@ -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 <typename T>
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
Expand All @@ -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.
Loading

0 comments on commit b3a5a6c

Please sign in to comment.