diff --git a/docs/conf.py b/docs/conf.py index b6aa721..d1a37b7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -92,6 +92,7 @@ html_theme = "pydata_sphinx_theme" html_theme_options = { + "show_nav_level": 2, "logo": { "text": f"ebtelplusplus {version}", }, @@ -119,7 +120,6 @@ html_sidebars = { "bibliography*": [], "development*": [], - "comparison*": [], } # Render inheritance diagrams in SVG graphviz_output_format = "svg" diff --git a/docs/development.rst b/docs/development.rst index 0417046..bf0ce70 100644 --- a/docs/development.rst +++ b/docs/development.rst @@ -16,7 +16,7 @@ Next, clone your fork, git clone https://github.com//ebtelPlusPlus.git -Because ``ebtelplusplus`` implements the actual simulation code in C++ and thus is not a pure Python package, installation of this package requires first compiling the C++ and building the needed binaries. +Because ``ebtelplusplus`` implements the actual simulation code in C++ and thus is not a pure Python package, installation of this package requires first compiling the C++ code and building the needed binaries. To do this, you will first need to have the `Boost `__ package installed (at least v1.53). There are various ways to install this package, including: @@ -29,6 +29,12 @@ There are various ways to install this package, including: * - `conda-forge `__ - any - ``conda install -c conda-forge libboost-devel`` + * - APT + - Linux (Debian-based) + - ``apt install libboost-all-dev`` + * - YUM + - Linux (RHEL) + - ``yum install boost-devel`` * - `Homebrew `__ - macOS - ``brew install boost`` diff --git a/docs/index.rst b/docs/index.rst index b9e0326..fcc70a7 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -6,23 +6,63 @@ ebtelplusplus :hidden: generated/gallery/index - comparison + topic_guides/index development reference bibliography -``ebtelplusplus`` (or ``ebtel++``) is a two-fluid version of the original enthalpy-based thermal evolution of loops (EBTEL) model implemented in C++ and wrapped in Python. -The EBTEL model, originally developed by :cite:t:`klimchuk_highly_2008`, efficiently computes spatially-averaged, time-dependent plasma parameters ( e.g. temperature, pressure, density) of dynamically-heated coronal loops. -It is often desirable to compute solutions for a large number of coronal loops, but the spatial and temporal scales needed to solve the full 1D-hydrodynamic equations lead to long compute times for even 1D hydrodynamic codes. -EBTEL computes quick and accurate solutions for spatially-averaged quantities, allowing efficient insight into how these monolithic structures evolve. -EBTEL also calculates the differential emission measure (DEM) for both the transition region and the corona. -Details regarding this formulation can be found in :cite:t:`klimchuk_highly_2008`. +``ebtelplusplus`` (or ``ebtel++``) is an implementation of the enthalpy-based thermal evolution of loops (EBTEL) model implemented in C++ and wrapped in Python. + +.. grid:: 1 2 2 2 + :gutter: 2 + + .. grid-item-card:: Examples + :link: generated/gallery + :text-align: center + + :material-outlined:`palette;8em;sd-text-secondary` + + Examples of how to set up a simulation + + .. grid-item-card:: Topic Guides + :link: ebtelplusplus-topic-guide + :link-type: ref + :text-align: center + + :material-outlined:`school;8em;sd-text-secondary` + + Detailed explanations about the EBTEL model + + .. grid-item-card:: Contributing + :link: ebtelplusplus-development + :link-type: ref + :text-align: center + + :material-outlined:`code;8em;sd-text-secondary` + + Instructions for how to contribute to `ebtelplusplus` + + .. grid-item-card:: Reference + :link: ebtelplusplus-reference + :link-type: ref + :text-align: center + + :material-outlined:`menu_book;8em;sd-text-secondary` + + Technical description of the inputs, outputs, and behavior of each component of ebteplusplus .. note:: If you are looking for the original IDL implementation, the repository for that code can be found `here `__. +The EBTEL model, originally developed by :cite:t:`klimchuk_highly_2008`, efficiently computes spatially-averaged, time-dependent plasma parameters ( e.g. temperature, pressure, density) of dynamically-heated coronal loops. +It is often desirable to compute solutions for a large number of coronal loops, but the spatial and temporal scales needed to solve the full 1D-hydrodynamic equations lead to long compute times for even 1D hydrodynamic codes. +EBTEL computes quick and accurate solutions for spatially-averaged quantities, allowing efficient insight into how these monolithic structures evolve. +EBTEL also calculates the differential emission measure (DEM) for both the transition region and the corona. +Details regarding this formulation can be found in :cite:t:`klimchuk_highly_2008`. +For a more detailed overview of the model, see :ref:`ebtelplusplus-topic-guide-derivation`. + Installation ------------ @@ -40,38 +80,18 @@ Once you've successfully installed ``ebtelplusplus``, see the Example Gallery fo Citation -------- -If you use ``ebtelplusplus`` in any published work, please include citations to the following papers: +If you use ``ebtelplusplus`` in any published work, we ask that you please include citations to the following papers: * :cite:t:`klimchuk_highly_2008` * :cite:t:`cargill_enthalpy-based_2012` * :cite:t:`cargill_enthalpy-based_2012-1` * :cite:t:`barnes_inference_2016` +If you make use of the variable abundance feature in the radiative losses, please also cite :cite:t:`reep_modeling_2024`. Additionally, please add the following line within your methods, conclusion or acknowledgements sections: *This research used version X.Y.Z (software citation) of the ebtelplusplus open source software package.* The software citation should be the specific `Zenodo DOI`_ for the version used within your work. -The citation for most current version on Zenodo is, - -.. code:: bibtex - - @software{Barnes2024, - author = {Barnes, Will and - Klimchuk, James and - Cargill, Peter and - Bradshaw, Stephen and - Reep, Jeffrey and - Schonfeld, Samuel and - Collazzo, Rowan}, - title = {rice-solar-physics/ebtelPlusPlus: v0.2}, - month = jul, - year = 2024, - publisher = {Zenodo}, - version = {v0.2}, - doi = {10.5281/zenodo.12675386}, - url = {https://doi.org/10.5281/zenodo.12675386} - } - .. _Zenodo DOI: https://zenodo.org/records/12675374 diff --git a/docs/reference.rst b/docs/reference.rst index dc46966..49123c0 100644 --- a/docs/reference.rst +++ b/docs/reference.rst @@ -3,6 +3,16 @@ API Reference ============= +ebtelplusplus +------------- + .. automodapi:: ebtelplusplus + :no-heading: + :no-main-docstr: + :no-inheritance-diagram: + +ebtelplusplus.models +-------------------- .. automodapi:: ebtelplusplus.models + :no-heading: diff --git a/docs/references.bib b/docs/references.bib index 53c6aa8..e7fd4f8 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -106,3 +106,17 @@ @article{reep_modeling_2024 annotation = {ADS Bibcode: 2024ApJ...970L..41R} } +@article{reep_geometric_2022, + title = {Geometric {{Assumptions}} in {{Hydrodynamic Modeling}} of {{Coronal}} and {{Flaring Loops}}}, + author = {Reep, Jeffrey W. and {Ugarte-Urra}, Ignacio and Warren, Harry P. and Barnes, Will T.}, + year = {2022}, + month = jul, + journal = {The Astrophysical Journal}, + volume = {933}, + pages = {106}, + issn = {0004-637X}, + doi = {10.3847/1538-4357/ac7398}, + urldate = {2023-08-29}, + annotation = {ADS Bibcode: 2022ApJ...933..106R} +} + diff --git a/docs/comparison.rst b/docs/topic_guides/comparison.rst similarity index 99% rename from docs/comparison.rst rename to docs/topic_guides/comparison.rst index 1409c2d..a5230b2 100644 --- a/docs/comparison.rst +++ b/docs/topic_guides/comparison.rst @@ -1,4 +1,4 @@ -.. _ebtelplusplus-comparison: +.. _ebtelplusplus-topic-guide-comparison: Why ebtelplusplus? ====================== diff --git a/docs/topic_guides/derivation.rst b/docs/topic_guides/derivation.rst new file mode 100644 index 0000000..597f364 --- /dev/null +++ b/docs/topic_guides/derivation.rst @@ -0,0 +1,272 @@ +.. _ebtelplusplus-topic-guide-derivation: + +Deriving the EBTEL Equations +============================ + +This page will briefly summarize the derivation of the version of the EBTEL +model used in `ebtelplusplus` which includes cross-sectional area expansion +and two-fluid effects. +For a more detailed explanation of the EBTEL model, including the physical +motivation behind the assumptions, see :cite:t:`klimchuk_highly_2008`. +Additionally, :cite:t:`barnes_inference_2016` provides a detailed discussion +of the two-fluid implementation while :cite:t:`cargill_static_2022` discusses +the inclusion of cross-sectional area expansion. + +The 1D Hydrodynamic Equations +----------------------------- + +Deriving the EBTEL model involves taking spatial integrals, +over both the transition region (TR) and corona, of the one-dimensional (1D) +hydrodynamic equations that govern the evolution of the dynamically-heated +solar atmosphere confined along a magnetic field line. +This yields a set of so-called "zero-dimensional" (0D) equations which +describe the evolution in time of the coronally-averaged temperature, density, +and pressure. +We begin by listing the relevant 1D hydrodynamic equations. +The 1D density, electron energy, and ion energy equations are given by, + +.. math:: + :label: density-1d + + \frac{\partial n}{\partial t} + \frac{1}{A}\frac{\partial}{\partial s}(Anv) = 0, + +.. math:: + :label: electron-energy-1d + + \frac{\partial E_e}{\partial t} + \frac{1}{A}\frac{\partial}{\partial s}(A(E_e + p_e)v) = &-\frac{1}{A}\frac{\partial}{\partial s}(AF_e) + v\frac{\partial p_e}{\partial s} - n^2\Lambda \\ + &+ \frac{k_Bn}{\gamma-1}\nu_{ie}(T_i - T_e) + Q_e, + +.. math:: + :label: ion-energy-1d + + \frac{\partial E_i}{\partial t} + \frac{1}{A}\frac{\partial}{\partial s}(A(E_i + p_i)v) = &-\frac{1}{A}\frac{\partial}{\partial s}(AF_i) - v\frac{\partial p_e}{\partial s} + m_i n v g_{\parallel} + \frac{1}{A}\left(\frac{4A\eta v}{3}\frac{\partial v}{\partial s}\right) \\ + &+ \frac{k_Bn}{\gamma-1}\nu_{ie}(T_e - T_i) + Q_i, + +where :math:`s` is the field-aligned spatial coordinate, :math:`n` is the density, +:math:`v` is the velocity, :math:`A` is the cross-sectional area, :math:`F_e` and +:math:`F_i` are the electron and ion thermal conductive fluxes, respectively, +:math:`T_e` and :math:`T_i` are the electron and ion temperatures respectively, +:math:`\nu_{ie}` is the electron-ion collision frequency, +:math:`\Lambda` describes the energy lost to radiation, and :math:`Q_e,Q_i` describe +the *ad hoc* heating applied to the electrons and ions, respectively.These equations +are closed by set of equations of state for the electron and ion energies, :math:`E_e,E_i`, +and pressures, :math:`p_e,p_i`, + +.. math:: + :label: eqs-of-state + + E_e &= \frac{p_e}{\gamma - 1}, \\ + E_i &= \frac{p_i}{\gamma - 1} + \frac{1}{2}m_inv^2, \\ + p_e &= nk_BT_e, \\ + p_i &= nk_BT_i. + +For additional details about these 1D hydrodynamic equations, see :cite:t:`reep_geometric_2022`. + +Using the above equations of state coupled with the assumptions that (1) all flows are +subsonic and (2) gravity is negligible for loops of length :math:`<150` Mm +:cite:p:`{see Section 2 of}klimchuk_highly_2008`, we can simplify :eq:`electron-energy-1d` +and :eq:`ion-energy-1d` to, + +.. math:: + :label: electron-energy-1d-simple + + \frac{A}{\gamma-1}\frac{\partial p_e}{\partial t} + \frac{\gamma}{\gamma-1}\frac{\partial}{\partial s}(Ap_ev) = &-\frac{\partial}{\partial s}(AF_e) + Av\frac{\partial p_e}{\partial s} - An^2\Lambda \\ + &+ Ak_Bn\nu_{ie}(T_i - T_e) + AQ_e, + +.. math:: + :label: ion-energy-1d-simple + + \frac{A}{\gamma-1}\frac{\partial p_i}{\partial t} + \frac{\gamma}{\gamma-1}\frac{\partial}{\partial s}(Ap_iv) = -\frac{\partial}{\partial s}(AF_i) - Av\frac{\partial p_e}{\partial s} + Ak_Bn\nu_{ie}(T_e - T_i) + AQ_i, + +We can now apply the combined methodology of both :cite:t:`barnes_inference_2016` and +:cite:t:`cargill_static_2022` to Eqs. :eq:`density-1d`, :eq:`electron-energy-1d-simple`, +and :eq:`ion-energy-1d-simple` to derive the EBTEL equations including both two-fluid +effects and cross-sectional expansion. + +The EBTEL Electron Pressure Equation +------------------------------------ + +To derive the EBTEL electron pressure equation, we begin by taking a spatial integral +over Eq. :eq:`electron-energy-1d-simple` from the base of the corona to the apex of a +semi-circular loop that is symmetric about the apex, + +.. math:: + :label: e-energy-coronal-integral + + \frac{A_cL_c}{\gamma-1}\frac{d p_{e,c}}{dt} - \frac{\gamma}{\gamma-1}(Ap_ev)_0 = (AF_e)_0 + A_cL_cQ_{e,c} + A_c\psi_c - A_cR_c, + +where :math:`c` denotes an average taken over the coronal portion of the loop, :math:`0` +denotes evaluation at the TR-corona interface, :math:`L_c` is the coronal portion of the +loop half-length, and, + +.. math:: + :label: psi-corona + + \psi_c = \frac{1}{A_c}\int_c\mathrm{d}s\,Av\frac{\partial p_e}{\partial s} + \frac{1}{A_c}\int_c\mathrm{d}s\,\frac{Ak_Bn\nu_{ie}}{\gamma-1}(T_i - T_e), + +.. math:: + :label: losses-corona + + R_c = \frac{1}{A_c}\int_c\mathrm{d}s\,An^2\Lambda. + +Note that the coronal integral :math:`\int_c\mathrm{d}s=\int_{s=L_{TR}}^{s=L}\mathrm{d}s`, +where :math:`L_{TR}` is the length of the TR and :math:`L=L_{TR}+L_c` is the total loop +half-length from the bottom of the TR to the apex of the loop. Additionally, because the +loop is assumed symmetric about the apex and isolated from the lower atmosphere, the +velocity and heat flux terms vanish at those locations. + +Similarly, we can integrate Eq. :eq:`electron-energy-1d-simple` over the TR, + +.. math:: + :label: e-energy-tr-integral + + \frac{A_{TR}L_{TR}}{\gamma-1}\frac{d p_{e,TR}}{dt} + \frac{\gamma}{\gamma-1}(Ap_ev)_0 = &-(AF_e)_0 + A_{TR}L_{TR}Q_{e,TR} \\ + &+ A_{TR}\psi_{TR} - A_{TR}R_{TR}, + +where :math:`TR` denotes an average taken over the TR portion of the loop, :math:`\psi_{TR}` +and :math:`R_{TR}` have the same form as Eqs. :eq:`psi-corona` and :eq:`losses-corona`, respectively, +and the TR integral :math:`\int_{TR}\mathrm{d}s=\int_{s=0}^{s=L_{TR}}`. + +Following the approach of :cite:t:`cargill_static_2022`, we add Eqs. +:eq:`e-energy-coronal-integral` and :eq:`e-energy-tr-integral` and let :math:`p_e=p_{e,c}=p_{e,TR}` +and :math:`Q_e=Q_{e,c}=Q_{e,TR}` to get the EBTEL electron pressure equation, + +.. math:: + :label: ebtel-electron-pressure + + \boxed{\frac{1}{\gamma-1}\frac{dp_e}{dt} = Q_e + \frac{\psi_c}{L_*}\left(1+\frac{A_{TR}\psi_{TR}}{A_c\psi_c}\right) - \frac{R_c}{L_*}\left(1+c_1\frac{A_{TR}}{A_c}\right)}, + +where :math:`L_* = L_c + (A_{TR}/A_c)L_{TR}` and :math:`c_1=R_{TR}/R_c`. Eq. :eq:`ebtel-electron-pressure` describes the time-evolution of the spatially-averaged electron pressure. + +The EBTEL Ion Pressure Equation +------------------------------------ + +To derive the EBTEL ion pressure equation, we apply the same procedure as above to Eq. +:eq:`ion-energy-1d-simple`. The spatial integral of Eq. :eq:`ion-energy-1d-simple` over +the coronal portion of the loop is, + +.. math:: + :label: i-energy-c-integral + + \frac{A_cL_c}{\gamma-1}\frac{d p_{i,c}}{dt} - \frac{\gamma}{\gamma-1}(Ap_iv)_0 = (AF_i)_0 + A_cL_cQ_{i,c} - A_c\psi_c, + +and for the TR portion of the loop is, + +.. math:: + :label: i-energy-tr-integral + + \frac{A_{TR}L_{TR}}{\gamma-1}\frac{d p_{i,TR}}{dt} + \frac{\gamma}{\gamma-1}(Ap_iv)_0 = -(AF_i)_0 + A_{TR}L_{TR}Q_{i,TR} - A_{TR}\psi_{TR}. + +As above, we can add Eqs. :eq:`i-energy-c-integral` and :eq:`i-energy-tr-integral` together +and let :math:`p_i=p_{i,c}=p_{i,TR}` and :math:`Q_i=Q_{i,c}=Q_{i,TR}` to obtain the EBTEL +ion pressure equation, + +.. math:: + :label: ebtel-ion-pressure + + \boxed{\frac{1}{\gamma-1}\frac{dp_i}{dt} = Q_i - \frac{\psi_c}{L_*}\left(1 + \frac{A_{TR}\psi_{TR}}{A_c\psi_c}\right)}. + +Eq. :eq:`ebtel-ion-pressure` describes the time-evolution of the spatially-averaged ion pressure. + +The EBTEL Density Equation +------------------------------------ + +Lastly, we derive the EBTEL density equation. +We begin by taking a spatial integral of Eq. :eq:`density-1d` over the coronal portion of the loop, + +.. math:: + + A_cL_c\frac{dn}{dt} - (Anv)_0 = 0 + +and using the equation of state for :math:`p_e` from Eq. :eq:`eqs-of-state`, + +.. math:: + :label: density-coronal-integral + + A_cL_c\frac{dn}{dt} = \frac{(Ap_ev)_0}{k_BT_{e,0}}. + +The quantity :math:`(Ap_ev)_0` is the area-weighted electron enthalpy flux at the TR-corona interface. +We can derive an expression for this term by substituting Eq. :eq:`ebtel-electron-pressure` into +Eq. :eq:`e-energy-tr-integral` and doing a lot of algebra, + +.. math:: + :label: electron-enthalpy-flux + + \frac{\gamma}{\gamma-1}(Ap_ev)_0 = -\frac{A_{TR}L_cR_c}{L_*}\left(c_1 - \frac{L_{TR}}{L_c}\right) + \frac{A_{TR}L_c\psi_c}{L_*}\left(\frac{\psi_{TR}}{\psi_c} - \frac{L_{TR}}{L_c}\right) - (AF_e)_0. + +Similarly, we can derive an expression for the area-weighted ion enthalpy flux by substituting +Eq. :eq:`ebtel-ion-pressure` into Eq. :eq:`i-energy-tr-integral`, + +.. math:: + :label: ion-enthalpy-flux + + \frac{\gamma}{\gamma-1}(Ap_iv)_0 = -\frac{A_{TR}L_c\psi_c}{L_*}\left(\frac{\psi_{TR}}{\psi_c} - \frac{L_{TR}}{L_c}\right) - (AF_i)_0. + +Adding Eqs. :eq:`electron-enthalpy-flux` and :eq:`ion-enthalpy-flux`, + +.. math:: + :label: total-enthalpy-flux + + \frac{\gamma}{\gamma-1}(Ap_ev)_0 + \frac{\gamma}{\gamma-1}(Ap_iv)_0 = -\frac{A_{TR}L_cR_c}{L_*}\left(c_1 - \frac{L_{TR}}{L_c}\right) - (AF_e)_0 - (AF_i)_0. + +Additionally, we define :math:`\xi\equiv T_e/T_i` and again use the electron and ion equations +of state from Eq. :eq:`eqs-of-state` to find, + +.. math:: + :label: temperature-ratio + + \xi = \frac{T_e}{T_i} = \frac{T_{e,0}}{T_{i,0}} = \frac{A_0n_0k_BT_{e,0}v_0}{A_0n_0k_BT_{i,0}v_0} = \frac{(Ap_ev)_0}{(Ap_iv)_0}. + +Substituting this into Eq. :eq:`total-enthalpy-flux` gives, + +.. math:: + :label: electron-enthalpy-flux-simple + + (Ap_ev)_0 = -\frac{\xi(\gamma - 1)}{\gamma(\xi + 1)}\left(\frac{A_{TR}L_cR_c}{L_*}\left(c_1 - \frac{L_{TR}}{L_c}\right) + (AF_e)_0 - (AF_i)_0\right). + +Finally, substituting Eq. :eq:`electron-enthalpy-flux-simple` into Eq. :eq:`density-coronal-integral` +yields the EBTEL density equation, + +.. math:: + :label: ebtel-density + + \boxed{\frac{dn}{dt} = -\frac{(\gamma - 1)\xi c_2}{(\xi + 1)\gamma c_3 k_B L_c T_e}\left(\frac{A_{TR}L_c}{A_cL_*}R_c\left(c_1 - \frac{L_{TR}}{L_c}\right) + \frac{A_0}{A_c}(F_{e,0} + F_{i,0})\right)}, + +where :math:`c_2,c_3` are constants relating the base temperature to the spatially-averaged coronal +temperature. Eq. :eq:`ebtel-density` describes the time-evolution of the spatial averaged density. +In summary, Eqs. :eq:`ebtel-electron-pressure`, :eq:`ebtel-ion-pressure`, and :eq:`ebtel-density` +comprise the two-fluid equations including cross-sectional area expansion. + +Limiting Behavior +----------------- + +Below, we briefly describe how Eqs. :eq:`ebtel-electron-pressure`, :eq:`ebtel-ion-pressure`, and +:eq:`ebtel-density` reduce to the other versions of the EBTEL model. + +Constant Cross-section +++++++++++++++++++++++ + +Under the assumption of constant cross-section, :math:`A_c=A_{TR}=A_0` and :math:`L_{TR}\ll L_c` +such that :math:`L_*\approx L_c = L`. As such, the EBTEL equations simplify to, + +.. math:: + + \frac{dp_e}{dt} &= (\gamma-1)Q_e + \frac{\gamma-1}{L}(\psi_c+\psi_{TR} - R_c(1+c_1)), \\ + \frac{dp_i}{dt} &= (\gamma-1)Q_i - \frac{\gamma-1}{L}(\psi_c + \psi_{TR}), \\ + \frac{dn}{dt} &= -\frac{(\gamma - 1)\xi c_2}{(\xi + 1)\gamma c_3 k_B L_c T_e}(c_1R_c + F_{e,0} + F_{i,0}). + +Note that these expressions are equivalent to the two-fluid EBTEL equations as given in :cite:t:`barnes_inference_2016`. + +Single-fluid +++++++++++++ + +Under the single-fluid assumption, :math:`T_e=T_i` at all times. +Using this and adding Eqs. :eq:`ebtel-electron-pressure` and :eq:`ebtel-ion-pressure`, + +.. math:: + + \frac{dp}{dt} &= (\gamma - 1)\left(Q - \frac{R_c}{L_*}\left(1 + \frac{A_{TR}}{A_c}c_1\right)\right), \\ + \frac{dn}{dt} &= -\frac{(\gamma - 1)c_2}{2\gamma c_3 k_B L_c T}\left(\frac{A_{TR}L_c}{A_cL_*}R_c\left(c_1 - \frac{L_{TR}}{L_c}\right) + \frac{A_0}{A_c}F_0\right). + +Note that these expressions are equivalent to the expanding cross-section EBTEL equations given +in :cite:t:`cargill_static_2022`. diff --git a/docs/topic_guides/index.rst b/docs/topic_guides/index.rst new file mode 100644 index 0000000..348dafd --- /dev/null +++ b/docs/topic_guides/index.rst @@ -0,0 +1,10 @@ +.. _ebtelplusplus-topic-guide: + +Topic Guides +============ + +.. toctree:: + :maxdepth: 1 + + comparison + derivation diff --git a/ebtelplusplus/__init__.py b/ebtelplusplus/__init__.py index 3b2b232..7437c69 100644 --- a/ebtelplusplus/__init__.py +++ b/ebtelplusplus/__init__.py @@ -1,11 +1,11 @@ """ ebtelplusplus: Zero-dimensional hydrodynamics of coronal loops """ -from ebtelplusplus.high_level import build_configuration, run +from ebtelplusplus.high_level import EbtelResult, run try: from ebtelplusplus._version import __version__ except ImportError: __version__ = "unknown" -__all__ = ["run", "build_configuration", "__version__"] +__all__ = ["run", "EbtelResult", "__version__"] diff --git a/ebtelplusplus/high_level.py b/ebtelplusplus/high_level.py index eb5ca23..5612e8e 100644 --- a/ebtelplusplus/high_level.py +++ b/ebtelplusplus/high_level.py @@ -2,45 +2,29 @@ High level interface to running an ebtel++ simulation """ import astropy.units as u +import dataclasses +import typing from astropy.utils.data import get_pkg_data_path -from collections import namedtuple import ebtelplusplus._low_level -from ebtelplusplus.models import DemModel, PhysicsModel, SolverModel +from ebtelplusplus.models import DemModel, HeatingModel, PhysicsModel, SolverModel -__all__ = ["run", "build_configuration"] - - -_UNITS_MAPPING = { - 'time': 's', - 'electron_temperature': 'K', - 'ion_temperature': 'K', - 'density': 'cm-3', - 'electron_pressure': 'dyne cm-2', - 'ion_pressure': 'dyne cm-2', - 'velocity': 'cm s-1', - 'heat': 'erg cm-3 s-1', - 'dem_temperature': 'K', - 'dem_tr': 'cm-5 K-1', - 'dem_corona': 'cm-5 K-1', - 'electron_thermal_conduction': 'erg cm-2 s-1', - 'ion_thermal_conduction': 'erg cm-2 s-1', - 'radiative_loss': 'erg cm-3 s-1', - 'tr_corona_radiative_loss_ratio': '', -} - - -EbtelResult = namedtuple('EbtelResult', - _UNITS_MAPPING.keys(), - defaults=[None,]*len(_UNITS_MAPPING)) +__all__ = ["run", "EbtelResult"] @u.quantity_input -def run(total_time: u.s, loop_length: u.cm, heating, physics=None, solver=None, dem=None): +def run(total_time: u.s, loop_length: u.cm, heating=None, physics=None, solver=None, dem=None): """ - Run an ebtel++ simulation + Run an ebtelplusplus simulation + + To run an ebtelplusplus simulation, at a minimum, you must specify a total simulation time + and a loop half length. Additional input parameters are configured through high-level objects + as listed below. Most users will likely want to modify the ``heating`` input by constructing + a `~ebtelplusplus.models.HeatingModel` that specifies the amount of energy injected into the + simulation as a function of time. Additional inputs are documented in the docstring of each + respective input model. Parameters ---------- @@ -48,35 +32,42 @@ def run(total_time: u.s, loop_length: u.cm, heating, physics=None, solver=None, Total duration of the simulation loop_length: `~astropy.units.Quantity` Loop half length - heating_model: `~ebtelplusplus.models.HeatingModel` + heating: `~ebtelplusplus.models.HeatingModel`, optional Configuration of heating model - physics_model: `~ebtelplusplus.models.PhysicsModel`, optional + physics: `~ebtelplusplus.models.PhysicsModel`, optional Configuration parameters related to the physics of the simulation - solver_model: `~ebtelplusplus.models.SolverModel`, optional + solver: `~ebtelplusplus.models.SolverModel`, optional Configuration parameters related to the numerical solver - dem_model: `~ebtelplusplus.models.DemModel`, optional + dem: `~ebtelplusplus.models.DemModel`, optional Configuration parameters related to the DEM calculation Returns ------- - results: `EbtelResult` - Data structure holding the results of the ebtel++ calculation + results: `~ebtelplusplus.EbtelResult` + Data structure holding the results of the ebtelplusplus calculation """ + if heating is None: + heating = HeatingModel() if solver is None: solver = SolverModel() if physics is None: physics = PhysicsModel() if dem is None: dem = DemModel() - config = build_configuration(total_time, loop_length, heating, physics, solver, dem) + config = _build_configuration(total_time, loop_length, heating, physics, solver, dem) results = ebtelplusplus._low_level.run(config) - return EbtelResult(**{k: u.Quantity(v, _UNITS_MAPPING[k]) for k,v in results.items()}) + return EbtelResult(results, config) @u.quantity_input -def build_configuration(total_time:u.s, loop_length:u.cm, heating, physics, solver, dem): +def _build_configuration(total_time:u.s, loop_length:u.cm, heating, physics, solver, dem): """ - Helper function for building a dictionary of ebtel++ configuration options + Helper function for building a dictionary of ebtelplusplus configuration options + + .. note:: + It is not necessary to use this function if you are running a simulation using + `~ebtelplusplus.run`. However, it may be useful to use this function directly + if you want to store the resulting inputs a dictionary. Parameters ---------- @@ -102,3 +93,95 @@ def build_configuration(total_time:u.s, loop_length:u.cm, heating, physics, solv **dem.to_dict(), 'heating': heating.to_dict(), } + + +@dataclasses.dataclass +class EbtelResult: + """ + Result of an ebtelplusplus simulation + + .. note:: + + This class is not meant to be instantiated directly. + Rather, it is meant to be returned by `~ebtelplusplus.run`. + + Parameters + ---------- + time: `~astropy.units.Quantity` + Simulation time at each time step. Has shape ``(N,)`` + electron_temperature: `~astropy.units.Quantity` + The temperature of the electrons at each time step + ion_temperature: `~astropy.units.Quantity` + The temperature of the ions at each time step + density: `~astropy.units.Quantity` + The density of the electrons and ions at each time step + electron_pressure: `~astropy.units.Quantity` + The pressure of the electrons at each time step + ion_pressure: `~astropy.units.Quantity` + The pressure of the ions at each time step + total_pressure: `~astropy.units.Quantity` + The total pressure at each time step. This is the sum of + the electron and ion pressures and is equivalent to the + pressure returned for single-fluid models. + velocity: `~astropy.units.Quantity` + The velocity at each timestep. This parameter should be + used cautiously as the EBTEL does not model the actual + transport of material through the spatial extent of the + loop. + heat: `~astropy.units.Quantity` + The energy released into the loop at each time step. + This is the amount of heating as specified in the + input `~ebtelplusplus.models.HeatingModel`. + electron_thermal_conduction: `~astropy.units.Quantity` + Electron thermal conductive flux at the TR-corona interface + at each time step. + ion_thermal_conduction: `~astropy.units.Quantity` + Ion thermal conductive flux at the TR-corona interface + at each time step. + radiative_loss: `~astropy.units.Quantity` + Energy lost due to radiation at each time step + tr_corona_radiative_loss_ratio: `~astropy.units.Quantity` + The ratio of the average radiative losses in the TR and the + corona. In the EBTEL model, this is often referred to as :math:`c_1`. + dem_temperature: `~astropy.units.Quantity` + The temperature grid on which the differential emission measure + distribution (DEM) is computed. Has shape ``(M,)``. These can be + interpreted as the centers of the temperature bins. + dem_tr: `~astropy.units.Quantity` + The TR DEM. Has shape ``(N,M)``. For more details of how this is + calculated, see Section 3 and the Appendix of :cite:t:`klimchuk_highly_2008`. + dem_corona: `~astropy.units.Quantity` + The coronal DEM. Has shape ``(N,M)``. + inputs: `dict` + All model inputs used to run the simulation. + """ + inputs: dict + time: u.Quantity[u.s] + electron_temperature: u.Quantity[u.K] + ion_temperature: u.Quantity[u.K] + density: u.Quantity[u.cm**(-3)] + electron_pressure: u.Quantity[u.dyne*u.cm**(-2)] + ion_pressure: u.Quantity[u.dyne*u.cm**(-2)] + velocity: u.Quantity[u.cm/u.s] + heat: u.Quantity[u.erg*u.cm**(-3)*u.s**(-1)] + electron_thermal_conduction: u.Quantity[u.erg*u.cm**(-2)*u.s**(-1)] + ion_thermal_conduction: u.Quantity[u.erg*u.cm**(-2)*u.s**(-1)] + radiative_loss: u.Quantity[u.erg*u.cm**(-3)*u.s**(-1)] + tr_corona_radiative_loss_ratio: u.Quantity[u.dimensionless_unscaled] + dem_temperature: u.Quantity[u.K] = None + dem_tr: u.Quantity[u.K**(-1)*u.cm**(-5)] = None + dem_corona: u.Quantity[u.K**(-1)*u.cm**(-5)] = None + # derived fields + total_pressure: u.Quantity[u.dyne*u.cm**(-2)] = dataclasses.field(init=None) + + def __init__(self, results, inputs): + # Use type hinting above to assign units to outputs; this is so + # we only have to specify the units of the outputs in one place + type_hints = type(self).__annotations__ + for k, v in results.items(): + _, unit = typing.get_args(type_hints[k]) + setattr(self, k, u.Quantity(v, unit)) + self.inputs = inputs + + def __post_init__(self): + self.total_pressure = self.electron_pressure + self.ion_pressure diff --git a/ebtelplusplus/models.py b/ebtelplusplus/models.py index c510534..ad2f313 100644 --- a/ebtelplusplus/models.py +++ b/ebtelplusplus/models.py @@ -1,5 +1,8 @@ """ -Models for configuring ebtel++ runs +Classes for configuring inputs to `ebtelplusplus.run`. +All model inputs have default values such that each model only need be instantiated +using the values of the parameters a user wants to change. +The default values of each input are listed below. """ import astropy.units as u import dataclasses @@ -18,7 +21,7 @@ @dataclasses.dataclass class PhysicsModel: """ - ebtel++ input parameters related to the physics of the simulation + ebtelplusplus input parameters related to the physics of the simulation Parameters ---------- @@ -87,7 +90,7 @@ def to_dict(self): @dataclasses.dataclass class SolverModel: """ - ebtel++ input parameters related to the numerical solver + ebtelplusplus input parameters related to the numerical solver Parameters ---------- @@ -121,11 +124,11 @@ def to_dict(self): @dataclasses.dataclass class DemModel: """ - ebtel++ input parameters related to differential emission measure (DEM) calculation. + ebtelplusplus input parameters related to differential emission measure (DEM) calculation. - Optionally, ebtel++ can can also calculate the differential emission - measure (DEM) in both the transition region and the corona. See sections 2. - 2 and 3 of :cite:t:`klimchuk_highly_2008` for the details of this + Optionally, ebtelplusplus can can also calculate the differential emission + measure (DEM) in both the transition region and the corona. See sections 2.2 + and 3 of :cite:t:`klimchuk_highly_2008` for the details of this calculation. Note that this will result in much longer computation times. Parameters @@ -164,7 +167,7 @@ class HeatingEvent: Each event has a linear rise phase, a constant phase, and a linear decay phase. Using this format, it is easy to specify either - symmetric or asymmetric events of many different shapes + symmetric or asymmetric events of many different shapes. Parameters ---------- @@ -255,9 +258,9 @@ def __init__(self, rise_start, duration, rate): @dataclasses.dataclass class HeatingModel: """ - ebtel++ input parameters for the time-dependent heating + ebtelplusplus input parameters for time-dependent heating - The ebtel++ time-dependent heating model is parameterized by a + The ebtelplusplus time-dependent heating model is parameterized by a series of discrete events (`HeatingEvent`) combined with a constant background heating rate. Furthermore, this energy can be injected into either the electrons or the ions or some admixture of the two. @@ -274,8 +277,8 @@ class HeatingModel: List of `HeatingEvent` objects that parameterize the energy injected into the loop by a series of discrete heating events. """ - background: u.Quantity[u.erg/(u.cm**3*u.s)] - partition: float = 1.0 + background: u.Quantity[u.erg/(u.cm**3*u.s)] = 1e-6 * u.Unit('erg cm-3 s-1') + partition: float = 0.5 events: list[HeatingEvent] = None def __post_init__(self): diff --git a/ebtelplusplus/tests/test_compare_idl.py b/ebtelplusplus/tests/test_compare_idl.py index 1c8701f..0dce7d3 100644 --- a/ebtelplusplus/tests/test_compare_idl.py +++ b/ebtelplusplus/tests/test_compare_idl.py @@ -7,7 +7,6 @@ import ebtelplusplus from ebtelplusplus.models import ( - DemModel, HeatingModel, PhysicsModel, SolverModel, @@ -57,16 +56,11 @@ def test_compare_idl_single_event(physics_model, heating_model, physics=physics_model, solver=solver_model) - config_dict = ebtelplusplus.build_configuration(total_time, - loop_length, - heating_model, - physics_model, - solver_model, - DemModel()) + r_idl = read_idl_test_data( 'idl_single_event.txt', ebtel_idl_path, - config_dict + r_cpp.inputs, ) if plot_idl_comparisons: plot_comparison(r_cpp, r_idl) @@ -98,13 +92,7 @@ def test_compare_idl_area_expansion( heating_model, physics=physics_model, solver=solver_model) - config_dict = ebtelplusplus.build_configuration(total_time, - loop_length, - heating_model, - physics_model, - solver_model, - DemModel()) - r_idl = read_idl_test_data(f'idl_area_expansion_{A_c=}_{A_0=}_{A_tr=}.txt', ebtel_idl_path, config_dict) + r_idl = read_idl_test_data(f'idl_area_expansion_{A_c=}_{A_0=}_{A_tr=}.txt', ebtel_idl_path, r_cpp.inputs) if plot_idl_comparisons: plot_comparison(r_cpp, r_idl) assert u.allclose(r_cpp.electron_temperature, r_idl['temperature'], rtol=RTOL) diff --git a/ebtelplusplus/tests/test_solver.py b/ebtelplusplus/tests/test_solver.py index 227fdd1..853b35a 100644 --- a/ebtelplusplus/tests/test_solver.py +++ b/ebtelplusplus/tests/test_solver.py @@ -99,6 +99,7 @@ def test_insufficient_heating(bad_heating, static_solver,): def test_NaNs_in_solver(use_adaptive_solver): solver = SolverModel(use_adaptive_solver=use_adaptive_solver) heating = HeatingModel( + partition=1.0, background=1e-6*u.Unit('erg cm-3 s-1'), events=[TriangularHeatingEvent(0.0*u.s, 200*u.s, -10*u.Unit('erg cm-3 s-1'))] ) diff --git a/pyproject.toml b/pyproject.toml index 205d172..f7f14bc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -135,6 +135,8 @@ lines_between_types = 1 [tool.codespell] skip = "*.fts,*.fits,venv,*.pro,*.bib,*.asdf,*.json" +ignore-words-list = "ebtel,hydrad" + [tool.ruff] target-version = "py310"