diff --git a/.codespell-ignore-words b/.codespell-ignore-words index 7f2a4307b2..d83562f9aa 100644 --- a/.codespell-ignore-words +++ b/.codespell-ignore-words @@ -1,3 +1,4 @@ +aprox blocs bloc inout @@ -15,3 +16,6 @@ crate vie rock4 rock2 +hep +ine + diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 0e8b2aea54..779583c71b 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -41,8 +41,8 @@ jobs: make realclean make NETWORK_DIR=ignition_reaclib/URCA-simple USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 4 - - name: compile test_nse_net (ase) + - name: compile nse_net_cell (ase) run: | export PATH=/usr/local/nvidia/bin:/usr/local/cuda/bin:${PATH} - cd unit_test/test_nse_net + cd unit_test/nse_net_cell make USE_CUDA=TRUE COMP=gnu USE_MPI=FALSE -j 4 diff --git a/.github/workflows/jac_cell.yml b/.github/workflows/jac_cell.yml index 91a05a46cd..5be0992002 100644 --- a/.github/workflows/jac_cell.yml +++ b/.github/workflows/jac_cell.yml @@ -34,7 +34,7 @@ jobs: - name: Run jac_cell run: | cd unit_test/jac_cell - ./main3d.gnu.ex inputs_aprox13 > test.out + ./main3d.gnu.ex inputs > test.out - name: Compare to stored output run: | diff --git a/.github/workflows/nse_net.yml b/.github/workflows/nse_net.yml index fc1124fe61..9680177f4e 100644 --- a/.github/workflows/nse_net.yml +++ b/.github/workflows/nse_net.yml @@ -26,42 +26,42 @@ jobs: sudo apt-get update -y -qq sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0 - - name: Compile, test_nse_net (NSE_NET, ase) + - name: Compile, nse_net_cell (NSE_NET, ase) run: | - cd unit_test/test_nse_net + cd unit_test/nse_net_cell make realclean make -j 4 - - name: Run test_nse_net (NSE_NET, ase) + - name: Run nse_net_cell (NSE_NET, ase) run: | - cd unit_test/test_nse_net + cd unit_test/nse_net_cell ./main3d.gnu.ex inputs_ase amrex.fpe_trap_{invalid,zero,overflow}=1 > test.out - name: Print backtrace - if: ${{ failure() && hashFiles('unit_test/test_nse_net/Backtrace.0') != '' }} - run: cat unit_test/test_nse_net/Backtrace.0 + if: ${{ failure() && hashFiles('unit_test/nse_net_cell/Backtrace.0') != '' }} + run: cat unit_test/nse_net_cell/Backtrace.0 - name: Compare to stored output (NSE_NET, ase) run: | - cd unit_test/test_nse_net + cd unit_test/nse_net_cell diff -I "^Initializing AMReX" -I "^AMReX" -I "^reading in reaclib rates" test.out ci-benchmarks/nse_net_unit_test.out - - name: Compile, test_nse_net/make_table (NSE_NET, ase, make_table) + - name: Compile, nse_net_cell/make_table (NSE_NET, ase, make_table) run: | - cd unit_test/test_nse_net/make_table + cd unit_test/nse_net_cell/make_table make realclean make -j 4 - - name: Run, test_nse_net/make_table (NSE_NET, ase, make_table) + - name: Run, nse_net_cell/make_table (NSE_NET, ase, make_table) run: | - cd unit_test/test_nse_net/make_table + cd unit_test/nse_net_cell/make_table ./main3d.gnu.ex amrex.fpe_trap_{invalid,zero,overflow}=1 > test.out - name: Print backtrace - if: ${{ failure() && hashFiles('unit_test/test_nse_net/make_table/Backtrace.0') != '' }} - run: cat unit_test/test_nse_net/make_table/Backtrace.0 + if: ${{ failure() && hashFiles('unit_test/nse_net_cell/make_table/Backtrace.0') != '' }} + run: cat unit_test/nse_net_cell/make_table/Backtrace.0 - name: Compare to stored output (NSE_NET, ase, make_table) run: | - cd unit_test/test_nse_net/make_table + cd unit_test/nse_net_cell/make_table diff -I "^Initializing AMReX" -I "^AMReX" -I "^reading in reaclib rates" test.out ci-benchmarks/nse_net_make_table_unit_test.out diff --git a/.github/workflows/test_partition_functions.yml b/.github/workflows/test_partition_functions.yml index 2752920610..8099a19721 100644 --- a/.github/workflows/test_partition_functions.yml +++ b/.github/workflows/test_partition_functions.yml @@ -28,17 +28,17 @@ jobs: - name: Compile run: | - cd unit_test/test_part_func + cd unit_test/part_func_cell make clean make -j 4 - - name: Run test_part_func + - name: Run part_func_cell run: | - cd unit_test/test_part_func + cd unit_test/part_func_cell ./main3d.gnu.ex > test.out - name: Compare to stored output run: | - cd unit_test/test_part_func + cd unit_test/part_func_cell diff -I "^Initializing AMReX" -I "^AMReX" -I "^reading in reaclib rates" test.out ci-benchmarks/part_func.out diff --git a/.gitignore b/.gitignore index 1cbd05a575..250e331460 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ *.ind Docs/*/*.pdf Docs/*/*.ps +Docs/source/runtime_parameters.rst MaestroUsersGuide.out MaestroUsersGuide.old @@ -79,7 +80,7 @@ extern.f90 build/ doxy_files Docs/source/runtime_parameters.rst - +Docs/source/changelog.md # C++ parameters extern_parameters.H diff --git a/CHANGES.md b/CHANGES.md index 07e961ead5..df0a3906f5 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,4 +1,55 @@ -# 25.01 +# Changelog + +## 25.02 + + * documentation updates (#1678, #1690, #1692, #1695, #1700, #1702, + #1703, #1707, #1709, #1712, #1713, #1714, #1715, #1716, #1723, + #1726, #1732, #1733) including renaming the docs directory to + `Docs` (#1689) and the addition of a CITATION.md (#1731) + + * codespell fixes (#1724) + + * rename `test_nse_net` -> `nse_net_cell`; `test_part_func` -> + `part_func_cell` (#1729) + + * remove old testing scripts (#1727) + + * update `test_react` `README.md` (#1722) + + * implement Debye-Huckel screening and allow it to be used as a test + for whether screening is needed (#1688) + + * remove the `use_raw_inputs` option from the EOS (#1721) + + * remove some old Fortran references (#1718, #1719) + + * switch from `std::clamp` to `amrex::Clamp` due to HIP compiler + issues (#1711) + + * reorganize the He nets -- they are now all under + `networks/he-burn` and share a common python setup (#1687, #1710) + also update these nets with pynucastro (#1685). This also + add a new network with 31. + + Renamed nets are: + * `subch_base` -> `he-burn/he-burn-18a` + * `subch_simple` -> `he-burn/he-burn-22a` + * `He-C-Fe-group` -> `he-burn/he-burn-36a` + * `CNO_He_burn` -> `he-burn/cno-he-burn-33a` + + * fix the Chapman-Jouguet detonation utility (#1699) + + * update the mailmap (#1706) + + * switch `std::pow(x, 1./3.)` to `std::cbrt(x)` (#1705) + + * remove `do_acc` option (#1708) + + * make the breakout EOS check if 1/mu is defined (#1694) + + * remove old Doxygen (#1697) and skynet (#1698) scripts from `util/` + +## 25.01 * update HIP/CUDA dependences to include sparse libraries (#1686) @@ -6,7 +57,7 @@ * update the integration and NSE docs (#1682) -# 24.12 +## 24.12 * documentation improvements (#1661, #1667, #1670) @@ -18,7 +69,7 @@ * `burn_cell` can now initialize all mass fractions to be equal (#1665) -# 24.10 +## 24.10 * metal chemistry updates (#1648) with ices (#1650) and cosmic rays (#1651) @@ -26,7 +77,7 @@ * doc updates (#1652) -# 24.09 +## 24.09 * Improvements to the primordial chemistry network and the addition of a new version that includes metals and dust (#1642, #1644) @@ -35,16 +86,16 @@ * documentation improvements (#1637) - * outputting the burn_t now prints the mass fractions / number densities + * outputting the `burn_t` now prints the mass fractions / number densities in scientific notation (#1643) * improvements to the looping and zeroing of the Jacobian in the integrators (#1636, #1640) -# 24.08 +## 24.08 * autodiff is now used with the templated reaction networks (#1614) - + some autodiff clean-ups and derivative fixes (#1604, #1612, + and some autodiff clean-ups and derivative fixes (#1604, #1612, #1613, #1616, #1619, #1633) * we can now output warnings from GPUs if you compile with @@ -88,11 +139,11 @@ * a fast log and fast pow approximation was added (#1591) - * the primordial_chem network now uses the fast math routines (#1605) + * the `primordial_chem network` now uses the fast math routines (#1605) * fix potential Inf in constexpr linear algebra (#1603) -# 24.07 +## 24.07 * added an autodiff library and converted all of the screening functions to use autodiff for the thermodynamic derivatives @@ -106,7 +157,7 @@ * fix return code for PrimordialChem unit test (#1590) - * NSE optimizations (including chabrier1998 screening) #1585 + * NSE optimizations (including `chabrier1998` screening) (#1585) * remove "using namespace amrex" from most headers (#1584) @@ -118,7 +169,7 @@ * retry tolerances now default to use the same values as the first attempt, unless they are explicitly set in an inputs file (#1573) -# 24.06 +## 24.06 * added the ability to access the runtime parameters via a struct. This will eventually be used to remove the dependency on globals @@ -143,7 +194,7 @@ * added an `eos_rh_t` EOS type (#1539) -# 24.05 +## 24.05 * Runtime parameters can now be type `bool` (#1536) @@ -155,7 +206,7 @@ * Update the pynucastro networks to cache derived rate partition functions (#1529) -# 24.04 +## 24.04 * A new `test_screening_templated` unit test was added -- this works with any of the templated networks. (#1525) @@ -169,7 +220,7 @@ * A `reinterpret_cast` in `rhs.H` was removed (#1435) -# 24.03 +## 24.03 * pivoting in the linear algebra routines can now be disabled (#1454) @@ -192,7 +243,7 @@ * added a zone-by-zone retry capability to the burner (#969) -# 24.02 +## 24.02 * Lots of general code cleaning from coverity and clang-tidy (#1450, #1452, #1453, #1460, #1459, #1457, #1458) @@ -200,9 +251,9 @@ * Fixed a bug in the VODE pivoting when a cached Jacobian is used (#1456) - * Added reverse rates to CNO_extras (#1445) + * Added reverse rates to `CNO_extras` (#1445) - * Sync networks up with pynucastro to get constexpr mion/bion + * Sync networks up with pynucastro to get `constexpr` mion/bion (#1437) * NSE+SDC improvements (#1431) @@ -210,12 +261,12 @@ * Start of moving the runtime parameters from globals to structs (#1422) -# 24.01 +## 24.01 * The quantum corrections for the Chabrier screening are now optional (#1428) - * We've replaced std::pow() with amrex::Math::powi for integer + * We've replaced `std::pow()` with `amrex::Math::powi` for integer powers for better GPU performance (#1432) * `in_nse` now works with an `eos_t` for tabular NSE (#1424) @@ -239,23 +290,23 @@ * constant T evolution was fixed (#1408) - * An unsafe reinterpret_cast was removed from linear algebra (#1412) + * An unsafe `reinterpret_cast` was removed from linear algebra (#1412) * Methods for computing T derivatives of an NSE table quantity were added (#1407) -# 23.12 +## 23.12 * The SDC+NSE update now includes plasma neutrino losses (#1357, #1400) * The default tabular NSE interpolation is now cubic (#1399) - * Self-consistent NSE now requires chabrier1998 or null screening + * Self-consistent NSE now requires `chabrier1998` or `null` screening (#1398) - * A new network, subch_base, was added that further simplifies - subch_simple (#1393) + * A new network, `subch_base`, was added that further simplifies + `subch_simple` (#1393) * A slightly larger network for Urca was added (#1365) @@ -265,9 +316,9 @@ * A bug was fixed in the neutrino cooling that was introduced in an optimization last release (#1380) -# 23.11 +## 23.11 - * The sneut5 neutrino cooling term was cleaned up (#1371, #1372, + * The `sneut5` neutrino cooling term was cleaned up (#1371, #1372, #1373, #1374, #1375, #1377, #1378, #1379) * The number of predictor-corrector iterations for the SDC+NSE algorithm @@ -276,16 +327,16 @@ * The Urca network now includes a more accurate rate for neutron decay and electon-capture onto a proton. (#1359) - * The He-C-Fe-group network now includes the positron parts of the + * The `He-C-Fe-group` network now includes the positron parts of the weak reaction rates (#1360) - * A check was added to ensure that the helm_table.dat is valid on + * A check was added to ensure that the `helm_table.dat` is valid on reading (#1355) -# 23.10 +## 23.10 * The simplified-SDC and true-SDC code paths for integration - have been merged (#1338, #1340, #1341). + have been merged (#1338, #1340, #1341). * All pynucastro networks have been updated with the latest version of pynucastro (2.1.0) (#1342) @@ -295,22 +346,22 @@ * `NUM_EXTRA_SPECIES` was removed (#1321) -# 23.09 +## 23.09 - * The file NETWORK_PROPERTIES has been removed from each network, + * The file `NETWORK_PROPERTIES` has been removed from each network, as the legacy screening method is no longer used. (#1310) - * The rprox network was updated and the Jacobian was fixed (#1300) + * The `rprox` network was updated and the Jacobian was fixed (#1300) - * The primordial_chem EOS now can take density and pressure as + * The `primordial_chem` EOS now can take density and pressure as inputs (#1302) -# 23.07 +## 23.07 - * The preprocessor variable EXTRA_THERMO has been removed. + * The preprocessor variable `EXTRA_THERMO` has been removed. Use cases that depend on dpdA/dpdZ or dedA/dedZ should use - eos_extra_t, which is a container that holds all of the - entities in eos_t as well as these derivatives wrt A and Z. (#1229) + `eos_extra_t`, which is a container that holds all of the + entities in `eos_t` as well as these derivatives wrt A and Z. (#1229) * added the ability to scale the energy we integrate by the initial energy in the ODE integration (#1224) @@ -318,25 +369,25 @@ * added an implementation of the Gershgorin circle theorem for estimating the spectral radius of our ODE system (#1222) - * removed SDC_EVOLVE_ENTHALPY -- this was not being used (#1204) + * removed `SDC_EVOLVE_ENTHALPY` -- this was not being used (#1204) -# 23.06 +## 23.06 * Added a new Runge-Kutta-Chebyshev integrator (#1191) * Lots of clean-up to the primordial chem network (#1180, #1181 #1198) - * Namespaces for the runtime parameters are now required in C++ ( + * Namespaces for the runtime parameters are now required in C++ (#1056) * The SDC+NSE update for tabular NSE was fixed -- we were previously computing the energy release incorrectly (#1092) -# 23.05 +## 23.05 - * The abort_on_failure runtime parameter has been removed (#1174) + * The `abort_on_failure` runtime parameter has been removed (#1174) -# 23.04 +## 23.04 * added preliminary CMake support (#1151, #1164, #1166) @@ -344,15 +395,15 @@ * clang-tidy code clean-ups(#1141, #1153, #1156) - * burn_t now stores whether we entered NSE (#1144, #1147) + * `burn_t` now stores whether we entered NSE (#1144, #1147) - * burn_t now store chemical potentials for NSE (#1149) + * `burn_t` now store chemical potentials for NSE (#1149) * some NSE solver updates to make it easier to enter NSE (#1138, #1139) * a new CNO + rp network for XRBs (#1145) -# 23.03 +## 23.03 * updated all of the pynucastro networks to the latest pynucastro version (#1134, #1136) @@ -365,30 +416,30 @@ * fixed some bugs in the NSE solver and made the hybrid Powell solver more robust (#1122) -# 23.02 +## 23.02 * `T_fixed` now works with NSE (#1098, #1111) * `USE_NSE` was changed to `USE_NSE_TABLE` (#1108) -# 23.01 +## 23.01 - * a new test, burn_cell_primordial_chem/, works with the primordial + * a new test, `burn_cell_primordial_chem`, works with the primordial chemistry (#1064) - * burn_cell and burn_cell_sdc now work with aux data with NSE + * `burn_cell` and `burn_cell_sdc` now work with `aux` data with NSE (#1084, #1094) * a screening implementation from Chabrier & Potekhin (1998) was added (#1085) - * test_react can now output the burn state that took the longest to + * `test_react` can now output the burn state that took the longest to evaluate (#967) * an initial implementation of adaptive nuclear statistic equilibrium was added (#983) -# 22.12 +## 22.12 * A first order backward Euler integrator was added that works with both Strang and simplified-SDC integration (#504, #1041, #1042, #1045) @@ -402,15 +453,15 @@ * An issue with reading the helmholtz table on GPUs was fixed (#1020) -# 22.11 +## 22.11 * use of the auxiliary state to define composition is now enabled - via USE_AUX_THERMO and the preprocessor variable AUX_THERMO + via `USE_AUX_THERMO` and the preprocessor variable `AUX_THERMO` (#1003) -# 22.10 +## 22.10 - * A null screening routine was added to disable screening for any network at + * A `null` screening routine was added to disable screening for any network at compile time. (#992) * An option to disable the clipping of species in the VODE integration @@ -426,36 +477,36 @@ * screening in the approximate rates in pynucastro nets was fixed (#978) -# 22.09 +## 22.09 * An NSE solver was added (#963) - * A new network, subch_simple, was added that further simplifies + * A new network, `subch_simple`, was added that further simplifies He/C burning (#964) -# 22.08 +## 22.08 - * The subch network was renamed `subch_approx` and the subch2 + * The `subch` network was renamed `subch_approx` and the `subch2` network was renamed `subch_full` (#947) -# 22.07 +## 22.07 * Two new screening formulations have been added for reaction rates, based on Chugunov, DeWitt, and Yakovlev 2007 and Chugunov and DeWitt 2009. These can be used with any network by setting - SCREEN_METHOD at compile time.(#887) + `SCREEN_METHOD` at compile time. (#887) -# 22.06 +## 22.06 - * The subch2 network now has runtime parameters allowing for + * The `subch2` network now has runtime parameters allowing for some key rates to be disabled (#921). -# 22.05 +## 22.05 - * The subch2 network was improved by adding some missing C+C, C+O, + * The `subch2` network was improved by adding some missing C+C, C+O, and O+O rates. (#915) -# 22.04 +## 22.04 * aprox networks now use a templated C++ righthand side formulation that builds the ODE system at compile time. (#802) @@ -463,7 +514,7 @@ * pynucastro networks were regenerated to take advantage of recent optimizations (#901) -# 22.02 +## 22.02 * The Microphysics repo was moved to the AMReX-Astro github organization: https://github.com/amrex-astro/Microphysics @@ -475,18 +526,18 @@ * Fortran support has been removed from the runtime parameter scripts (#869) -# 22.01 +## 22.01 * we added back in support for the "Nonaka plot". This outputs the state in the RHS routine for a single zone during the reaction network integration. (#830) - * we removed the xrb_simple network. This was never used in any + * we removed the `xrb_simple` network. This was never used in any science calculations (#827) * the simplified-SDC step rejection logic in VODE was improved (#818) -# 21.12 +## 21.12 * all of the pynucastro networks were regenerated with the latest pynucastro and converted to C++. Performance was also improved @@ -494,32 +545,32 @@ * a bug was fixed in the VODE step rejection logic (#815) - * we added USE_MICROPHYSICS_DEBUG that defines MICROPHYSICS_DEBUG to + * we added `USE_MICROPHYSICS_DEBUG` that defines `MICROPHYSICS_DEBUG` to turn on more verbosity to help with debugging (#817) -# 21.11 +## 21.11 - * burn_cell was not correctly doing substepping in some cases. + * `burn_cell` was not correctly doing substepping in some cases. This has been fixed (#784) - * With Intel compilers, logical runtime parameters in Fortran - were not being correctly cast to int (#789) + * With Intel compilers, `logical` runtime parameters in Fortran + were not being correctly cast to `int` (#789) * Simplified-SDC now works with Fortran nets (#786) -# 21.09 +## 21.09 - * Added a new nova network (nova2) with pp and (hot-)CNO and some + * Added a new nova network (`nova2`) with pp and (hot-)CNO and some breakout reactions (#751) - * Some fixes to the NSE bailout in aprox19 (#739, #753, #755) and + * Some fixes to the NSE bailout in `aprox19` (#739, #753, #755) and the relaxation check on the NSE criteria (#754) - * Added a new unit test for single-zone SDC (burn_cell_sdc) (#744) + * Added a new unit test for single-zone SDC (`burn_cell_sdc`) (#744) -# 21.08 +## 21.08 - * test_react can now be run with a uniform composition to test GPU + * `test_react` can now be run with a uniform composition to test GPU performance (#734) * the numerical Jacobian now uses a more robust one-sided difference @@ -531,46 +582,46 @@ * for the NSE bailout, we can now relax the conditions needed to enter NSE after a failed burn (#702) -# 21.07 +## 21.07 - * The C++ networks now implement abort_on_failure functionality + * The C++ networks now implement `abort_on_failure` functionality (#697) -# 21.06 +## 21.06 * The ability to use a system BLAS library was removed (#675) * An equation of state for hypervelocity impacts was added (#645) -# 21.05 +## 21.05 * For aprox19 + NSE, we now "bail out" of the integration immediately if the state enters NSE, and then do the rest of the update through the NSE table. (#658) - * The old gamma_law EOS was removed and gamma_law_general was - renamed gamma_law. The old gamma_law EOS have a very reduced + * The old `gamma_law` EOS was removed and `gamma_law_general` was + renamed `gamma_law`. The old `gamma_law` EOS have a very reduced subset of thermodynamic quantities that it computed, for efficiency purposes. This is no longer needed now that we have - templated the EOSes and have different eos_t data types (#653). + templated the EOSes and have different `eos_t` data types (#653). * Integration for simplified-SDC was interpreting rtol incorrectly. This has been fixed (#643) - * Screening for the 3-alpha reaction in the subch, subch2, and nova + * Screening for the 3-alpha reaction in the `subch`, `subch2`, and `nova` networks was fixed (#627, #634, #635) -# 21.04 +## 21.04 * We added a new mechanism to recover a failed burn when the state - tries to enter NSE during the evolution, when using the aprox19 + + tries to enter NSE during the evolution, when using the `aprox19` + NSE network. Now it will capture the failure and redo the burn if it satisfies the NSE criteria (#628) * We updated the VODE logic for rejecting a step to consider mass fractions for both simplified-SDC and true SDC burns (#619) -# 21.03 +## 21.03 * We now integrate internal energy (e) directly instead of integrating temperature (T) for the thermodynamic evolution. T is @@ -580,7 +631,7 @@ * simplified-SDC can be used with the NSE table in aprox19 now (#423, #497) -# 21.02 +## 21.02 * Fortran support for the VODE integrator has been removed (#538) @@ -592,58 +643,58 @@ * Fortran support for simplified-SDC in the VODE integrator has been removed. (#492) -# 21.01 +## 21.01 * Microphysics now requires C++17 (gcc >= 7, CUDA >= 11). (#485) * The BS integrator was removed. This was Fortran only, doesn't support SDC integration, and not well used. (#488) -# 20.12 +## 20.12 - * The default absolute tolerance for species (atol_spec) has been - increased to 1.e-8 (from 1.e-12). (#422) + * The default absolute tolerance for species (`atol_spec`) has been + increased to `1.e-8` (from `1.e-12`). (#422) * An interface has been added for C++ integrators to call the RHS from a network that only has a Fortran implementation. This allows - the use of USE_CXX_REACTIONS = TRUE for any network (however, CUDA + the use of `USE_CXX_REACTIONS = TRUE` for any network (however, CUDA is not currently supported for this case). (#419) -# 20.11 +## 20.11 - * The aprox19 + NSE network was ported to C++ (#362) + * The `aprox19` + NSE network was ported to C++ (#362) * The simplified-SDC code path was ported to C++ (#389) -# 20.10 +## 20.10 * An option to use NSE instead of integrating the reaction - network has been added to the aprox19 network. (#332) + network has been added to the `aprox19` network. (#332) * The BS integrator no longer supports simplified-SDC (#393) - * The triple_alpha_plus_cago network switch to using binding + * The `triple_alpha_plus_cago` network switch to using binding energies in MeV, consistent with the aprox nets (#354) -# 20.09 +## 20.09 - * Unit tests now write a job_info file (#383) + * Unit tests now write a `job_info` file (#383) - * A new single-zone EOS test routine was created as unit_test/eos_cell + * A new single-zone EOS test routine was created as `unit_test/eos_cell` (#382) - * The gamma_law eos (not gamma_law_general) now fills the sound + * The `gamma_law` eos (not `gamma_law_general`) now fills the sound speed, entropy, and derivatives for more inputs (#374) - * The rprox network now has screening (#377) + * The `rprox` network now has screening (#377) - * The NETWORK_PROPERTIES file was split to put the number of - auxiliary species into its own file, NAUX_NETWORK. This allows + * The `NETWORK_PROPERTIES` file was split to put the number of + auxiliary species into its own file, `NAUX_NETWORK`. This allows us to put if-logic into the file to choose the number of - auxiliary quantities based on make setting (like USE_NSE). + auxiliary quantities based on make setting (like `USE_NSE`). (#370) -# 20.08 +## 20.08 * Several of the unit tests had separate C++ and Fortran implementations. These have been unified (#343, #344, #345) @@ -654,68 +705,68 @@ change by more than a factor of 2, or an abundance < 0 or > 1, as well as timesteps where the temperature ends up negative. (#350) -# 20.07 +## 20.07 - * The "master" branch has been renamed "main" (#333) + * The `master` branch has been renamed `main` (#333) - * NETWORK_PROPERTIES now includes the number of Aux quantities (#330) + * `NETWORK_PROPERTIES` now includes the number of Aux quantities (#330) -# 20.06 +## 20.06 - * For integration with simplified SDC, we now interpret atol_spec + * For integration with simplified SDC, we now interpret `atol_spec` as an absolute tolerance on X alone instead of (rho X) (#311) - * burn_cell can now use the C++ burner if compiled with - USE_CXX_REACTIONS=TRUE and run with do_cxx = 1. (#313) + * `burn_cell` can now use the C++ burner if compiled with + `USE_CXX_REACTIONS=TRUE` and run with `do_cxx = 1`. (#313) - * The original burn_cell (which used the F90 BoxLib build system) - is removed and replaced with burn_cell_C (which uses the newer + * The original `burn_cell` (which used the F90 BoxLib build system) + is removed and replaced with `burn_cell_C` (which uses the newer build system). (#316) * The analytic Jacobian with simplified SDC now is written in terms of the conserved fluid state and works for a wide range of problems (#228) -# 20.05 +## 20.05 - * We now have an option for using sparse storage for aprox13 in C++ + * We now have an option for using sparse storage for `aprox13` in C++ (#307) - * iso7 and aprox13 are now available as a C++ network (#303, 305) + * `iso7` and `aprox13` are now available as a C++ network (#303, #305) - * species names are available as an enum in network_properties.H (#304) + * species names are available as an `enum` in `network_properties.H` (#304) - * The screening on O16+O16 in iso7 was fixed (#302) + * The screening on O16+O16 in `iso7` was fixed (#302) * The VODE integrator is now available in C++ (#299) -# 20.04 +## 20.04 - * The wion network property was removed (#294) + * The `wion` network property was removed (#294) * There are new unit tests for the screening and aprox rates modules (both C++ and Fortran interfaces). - * The screening routines were ported to C++ (#290) and the screenz - routine was removed in favor of screen5 (#293) + * The screening routines were ported to C++ (#290) and the `screenz` + routine was removed in favor of `screen5` (#293) - * a new method, is_input_valid, was added to all EOSes (both C++ + * a new method, `is_input_valid`, was added to all EOSes (both C++ and Fortran interfaces) that can be used to query whether an EOS - supports a particular input mode (e.g. eos_input_rp). (#291) + supports a particular input mode (e.g. `eos_input_rp`). (#291) - * The aprox rates used with iso7, aprox13, aprox19, and aprox21 + * The aprox rates used with `iso7`, `aprox13`, `aprox19`, and `aprox21` have been converted to C++ (#288) * We've rewritten the VODE integrator to remove all "go to" - statements (#275, 276, 278, 280, 281, 282, 283, 284, 285, 286, - 287) + statements (#275, #276, #278, #280, #281, #282, #283, #284, #285, + #286, #287) - * We removed the ability to have nspec_evolve < nspec. This + * We removed the ability to have `nspec_evolve` < `nspec`. This feature was not widely used and greatly complicated the code - paths.(#279) + paths. (#279) -# 20.03 +## 20.03 * The nuclei information for both Fortran and C++ is now automatically generated from a network inputs file at compile @@ -723,37 +774,37 @@ constant (#253, 258). * The license for StarKiller Microphyscs was made explicit and - a license.txt file was added (#267) + a `license.txt` file was added (#267) * A framework for pure C++ EOSes has been created and a pure C++ - unit test, test_eos_C, is available to test these. (#246) The - following EOSes have been ported to C++: ztwd (#268), multigamma - (#265), polytrope (#264), gamma_law (#263), helmholtz (#262), - gamma_law_general (#246), rad_power_law (#269), breakout (#270) + unit test, `test_eos_C`, is available to test these. (#246) The + following EOSes have been ported to C++: `ztwd` (#268), `multigamma` + (#265), `polytrope` (#264), `gamma_law` (#263), `helmholtz` (#262), + `gamma_law_general` (#246), `rad_power_law` (#269), `breakout` (#270) - * The GPackage.mak files that were a remnant of the old + * The `GPackage.mak` files that were a remnant of the old BoxLib F90 build system have been removed. They were not maintained. (#212). * All of the interface files have been collected together - in the interfaces/ dir. (#240) + in the `interfaces/` dir. (#240) - * The network C++ headers have been renamed network_properties.H + * The network C++ headers have been renamed `network_properties.H` and the nuclei information (aion and zion) have been added. (#244) -# 20.02 +## 20.02 - * Added a C++ header file, actual_network.H, that defines the + * Added a C++ header file, `actual_network.H`, that defines the network size. This is the start of making making the microphysics routines available in C++. * regenerated the pynucastro networks with the latest weak rate formulations from pynucastro. -# 20.01 +## 20.01 - * The burn_t type no longer includes ydot or jac -- this allows + * The `burn_t` type no longer includes `ydot` or `jac` -- this allows us to optimize the memory access on GPUs (#220) * The radiation pressure contribution to the Helmholtz EOS has @@ -764,11 +815,11 @@ * The original VODE integrator was removed and the Fortran 90 version VODE90 was renamed to VODE. (#221) - * The test_react unit tests no longer require a composition inputs - file (xin*). They now create the composition profile at runtime. + * The `test_react` unit tests no longer require a composition inputs + file (`xin*`). They now create the composition profile at runtime. (#211) -# 19.12 +## 19.12 * Simplified SDC integration now uses the same retry strategy as the default (non-SDC) integration. (#215) @@ -778,7 +829,7 @@ switch to the BS integrator if loosening the tolerances does not allow the burn to complete. (#201) - * The parameter ode_max_steps was made consistent in VODE and + * The parameter `ode_max_steps` was made consistent in VODE and VODE90; in some places it was being ignored. (#214) * The helmholtz EOS was restructured, splitting the different @@ -786,60 +837,60 @@ accesses. (#200) * The derivatives with respect to mass fraction (dpdX, dedX, dhdX) - were removed from eos_t and are now available through a new type, - eos_xderivs_t and the composition_derivatives() routine. (#207) + were removed from `eos_t` and are now available through a new type, + `eos_xderivs_t` and the `composition_derivatives()` routine. (#207) - * A bug in the screening of the C12+C12 and O16+O16 rates in iso7 - was fixed. + * A bug in the screening of the O16+O16 rate in iso7 was + fixed. (#204) - * The test_eos unit test now outputs all of the variables in the - eos_t type. + * The `test_eos` unit test now outputs all of the variables in the + `eos_t` type. -# 19.11 +## 19.11 * VODE90 now works with the simplified SDC time step algorithms, and the preprocessor option for this SDC was changed to - SIMPLIFIED_SDC (#194) + `SIMPLIFIED_SDC` (#194) - * rprox now works on GPUs + * `rprox` now works on GPUs -# 19.10 +## 19.10 - * The iso7 network was ported to GPUs (#172) + * The `iso7` network was ported to GPUs (#172) * VODE90 now better agrees with VODE (#192) - * When building with multiple integrators, the contents of the rpar - modules could clash. This has been fixedc (#136) + * When building with multiple integrators, the contents of the `rpar` + modules could clash. This has been fixed. (#136) * A module for making the "Nonaka plot" tracking the evolution of a quantity during the burn was added, and is enabled with - USE_NONAKA_PLOT=TRUE + `USE_NONAKA_PLOT=TRUE` (#190) -# 19.08 +## 19.08 - * A new network, subch2, was added that combines aprox13 and the - subch networks. (#184) + * A new network, `subch2`, was added that combines `aprox13` and the + `subch` networks. (#184) -# 19.05 +## 19.05 - * The aprox21 network was missing the analytic Jacobian term for + * The `aprox21` network was missing the analytic Jacobian term for the derivative of He4 with respect to Ni56. This is fixed. (#175) * The numerical Jacobian module, used by the BS and VBDF integrators had some wrong scalings. These have now been fixed (#179, #180) -# 19.01 +## 19.01 * the docs are now automatically build from the sphinx source using travis on github. -# 18.12 +## 18.12 * simplify conductivity interface to match the eos interface by moving the conductivity into the eos type -# 18.11 +## 18.11 * new python bindings have been added @@ -850,114 +901,113 @@ was being applied incorrectly. -# 18.10 +## 18.10 - * test_eos and test_react now both work on GPUs (using the AMReX + * `test_eos` and `test_react` now both work on GPUs (using the AMReX `gpu` branch) * the intermediate blending of the weak and strong screening regimes was wrong, and has been fixed. We've also synced some parameters - up to agree with those in MESA and Kepler. (#149, 150) + up to agree with those in MESA and Kepler. (#149, #150) - * eos_input_is_constant is now set to true for the helmholtz EOS. + * `eos_input_is_constant` is now set to `true` for the helmholtz EOS. This mean that the EOS inputs will not be modified after the EOS call. This is good for conserving energy in a hydro code, but the tradeoff is a small (to root finding tolerance) inconsistency in - the thermodynamic state. + the thermodynamic state. (#154) -# 18.09 +## 18.09 - * The Helmholtz parameters ttol and dtol (controlling the error - for the Newton iteration when in a mode other than eos_input_rt) - are now runtime parameters in the extern namelist as eos_ttol - and eos_dtol. + * The Helmholtz parameters `ttol` and `dtol` (controlling the error + for the Newton iteration when in a mode other than `eos_input_rt`) + are now runtime parameters in the extern namelist as `eos_ttol` + and `eos_dtol`. -# 18.08 +## 18.08 - * the unit tests (test_react, test_sdc, and test_eos) have been + * the unit tests (`test_react`, `test_sdc`, and `test_eos`) have been ported from the Fortran to C++ build system in AMReX. This will allow us to test the GPU framework in AMReX. -# 18.07 +## 18.07 * added CUDA support to the VODE90 integrator, the helmeos, and the - networks aprox13, aprox19, aprox21, ignition_simple, - C-burn-simple, URCA-simple. + networks `aprox13`, `aprox19`, `aprox21`, `ignition_simple`, + `C-burn-simple`, `URCA-simple`. * Ported the unit test frameworks to FBoxLib -# 18.05 +## 18.05 * lots of documentation updates - * some fixes to the numerical Jacobian involving X vs. Y + * some fixes to the numerical Jacobian involving X vs. Y (#100) - * a new subCh network for He burning was added. + * a new `subCh` network for He burning was added. * implemented the new c12(a,g)o16 nuclear reaction rate and its - corresponding inverse from the work of Deboer et al. 2017 - (https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.89.035007). - To use the new rate, user must set `use_c12ag_deboer17` to `true`. - This rate is only usable in the `aprox13`, `aprox19`, `aprox21`, - and `iso7` reaction rate networks. Closes issue #44. + corresponding inverse from the work of Deboer et al. 2017 (Rev Mod + Phys 89, 035007, 2017). To use the new rate, user must set + `use_c12ag_deboer17` to `true`. This rate is only usable in the + `aprox13`, `aprox19`, `aprox21`, and `iso7` reaction rate + networks. (#56) - * a routine util/cj_detonation was added to compute the + * a routine `util/cj_detonation` was added to compute the Chapman-Jouguet detonation velocity for any of the networks * the burn retry strategy now sticks with the current integrator and - uses looser tolerances before switching to a different integrator. + uses looser tolerances before switching to a different integrator. (#96) -# 18.04 +## 18.04 * pynucastro (https://github.com/pynucastro/pynucastro) can now generate reaction networks compatible with StarKiller. See the - subch network. + `subch` network. -# 17.11 +## 17.11 * a new option to boost the reaction rates has been added - to the integrators (PR #64) + to the integrators (#64) * we now disable some composition derivatives in the EOS by default, for performance and memory reasons. They can be re-enabled by defining the preprocessor variable - EXTRA_THERMO (PR #59) + `EXTRA_THERMO` (#59) -# 17.10 +## 17.10 * the compositional derivatives are no longer available by default from the EOS. To get these, set the preprocessor - variable EXTRA_THERMO. This change was done for performance + variable `EXTRA_THERMO`. This change was done for performance reasons. - * the aprox19 and aprox21 networks no longer use a numerical + * the `aprox19` and `aprox21` networks no longer use a numerical Jacobian by default, as this was found to result in some - bad numerical issues in VODE (PR #49) + bad numerical issues in VODE (#49) - * the maximum temperature for reactions, MAX_TEMP, is now + * the maximum temperature for reactions, `MAX_TEMP`, is now an adjustable input parameter rather than being hardcoded - at 1.d11. + at `1.d11`. * the Helmholtz EOS table is now read by the IO processor and - broadcast to other processors (PR #53) + broadcast to other processors (#53) * the VODE integrator now does some additional checks on the - state to ensure consistency (PR #47) + state to ensure consistency (#47) -# 17.09 +## 17.09 * a new rety mechanism was implemented that allows a different - integrator to be used if the primary integrator fails + integrator to be used if the primary integrator fails (#39) * the electron Ni56 electron capture rates and energy losses were updated from Mazurek (1973) to LMP (2000). Thanks to - Carl Fields for this contribution. Pull request #40 + Carl Fields for this contribution. (#40) +## 17.08 -# 17.08 - - * fix to aprox21 from Aron Michel (HITS) that fills in missing + * fix to `aprox21` from Aron Michel (HITS) that fills in missing reactions * updated the helmholtz EOS to use the latest table from Frank @@ -967,38 +1017,34 @@ * add stellar conductivities from Frank Timmes - -# 17.06 +## 17.06 * a new Fortran 90 port of VODE has been added * the unit tests now require AMReX instead of BoxLib to build - -# 17.01 +## 17.01 * we've removed the option to integrate molar fractions and instead the ODE system always operates on mass fractions (the networks return derivatives of molar fractions and they are automatically converted). +## 16.12 -# 16.12 - - * a new unit test, test_sdc, was created to test the SDC interface + * a new unit test, `test_sdc`, was created to test the SDC interface to the networks - * we now rely on the network module to provide aion_inv (1/aion) + * we now rely on the network module to provide `aion_inv` (1/aion) * the VODE integrator now supports SDC integration - -# 16.09 +## 16.09 * num_rate_groups is now a property of the individual networks * a new integration method, Rosenbrock, was added to the BS - option (set ode_method) + option (set `ode_method`) * the number of RHS and Jac evaluations is now passed out of the burner through the burn_t type for diagnostic and @@ -1007,21 +1053,20 @@ * support for spectral deferred correction coupling of the burner and hydro was added to the BS integrator +## 16.08 -# 16.08 - - * Microphysics/eos/ has been renamed Microphysics/EOS/ to better + * `Microphysics/eos/` has been renamed `Microphysics/EOS/` to better conform to the conventions used in Castro and Maestro * the User's Guide has been extensively updated * OpenMP and OpenACC directives have been added to the unit tests - * the BS integrator's type, bs_t, has now contains a burn_t - internally, simplifying the conversion from bs_t to call the - actual_rhs/jac + * the BS integrator's type, `bs_t`, has now contains a `burn_t` + internally, simplifying the conversion from `bs_t` to call the + `actual_rhs/jac` - * the rates() component of burn_t was removed. We no longer + * the `rates()` component of `burn_t` was removed. We no longer rely on rate caching * we now store the simulation time at the start of the burn as @@ -1029,38 +1074,38 @@ time * the species derivatives (dh/dX and de/dX) and enthalpy were - removed from the burn_t + removed from the `burn_t` * a new option to integrate of X instead of Y was added - (integrate_molar_fraction = F) + (`integrate_molar_fraction = F`) - * integration of networks with nspec_evolve < nspec were fixed + * integration of networks with `nspec_evolve` < `nspec` were fixed to now apply the algrebic constraint relating mass fractions - through a new update_unevolved_species() function + through a new `update_unevolved_species()` function - * the electron capture rate on Ni56 used by aprox19 and aprox21 was + * the electron capture rate on Ni56 used by `aprox19` and `aprox21` was fixed * the BS integrator can now use the initial timestep estimation - algorithm that VODE uses 9use_timestep_estimator = T) + algorithm that VODE uses (`use_timestep_estimator = T`) * a centered difference numerical Jacobian option was added -# 16.07 +## 16.07 - * we now use MICROPHYSICS_HOME instead of MICROPHYSICS_DIR as the - environment variable to point to the Microphysics/ directory. + * we now use `MICROPHYSICS_HOME` instead of `MICROPHYSICS_DIR` as the + environment variable to point to the `Microphysics/` directory. - * there are now two standalone unit tests, test_react and test_eos + * there are now two standalone unit tests, `test_react` and `test_eos` that don't need Maestro or Castro to compile. * new burn mode that limits numerically unstable burning. - * UsersGuide/ was renamed to Docs/ to be consistent with the other + * `UsersGuide/` was renamed to `Docs/` to be consistent with the other BoxLib codes * the energy equation now uses an offset to help with the BS ODE integration convergence - * the runtime parameter small_x now is owned by the network + * the runtime parameter `small_x` now is owned by the network diff --git a/CITATION.md b/CITATION.md new file mode 100644 index 0000000000..93941cbdc4 --- /dev/null +++ b/CITATION.md @@ -0,0 +1,2 @@ +If you use AMReX-Astro Microphysics, please cite the Zenodo DOI. +The latest version is always here: https://doi.org/10.5281/zenodo.2620544 diff --git a/Docs/Makefile b/Docs/Makefile index cc4d477864..bf411458d3 100644 --- a/Docs/Makefile +++ b/Docs/Makefile @@ -21,6 +21,8 @@ clean: html: ./rp.py > source/runtime_parameters.rst + curl -L -H 'Accept: application/x-bibtex' https://zenodo.org/api/records/2620544 > source/zenodo.bibtex.txt + python parse_changelog.py ../CHANGES.md > source/changelog.md doxygen Doxyfile breathe-apidoc --o source doxy_files/xml @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/Docs/parse_changelog.py b/Docs/parse_changelog.py new file mode 100644 index 0000000000..6305a01ec3 --- /dev/null +++ b/Docs/parse_changelog.py @@ -0,0 +1,25 @@ +import argparse +import re + +PR_URL_BASE = r"https://github.com/AMReX-Astro/Microphysics/pull/" + +pr = re.compile(r"(\#)(\d+)") + + +def doit(clfile): + + with open(clfile) as cl: + for line in cl: + new_line = re.sub(pr, rf"[\g<0>]({PR_URL_BASE}\g<2>)", line) + print(new_line, end="") + + +if __name__ == "__main__": + + parser = argparse.ArgumentParser() + parser.add_argument("changelog", type=str, nargs=1, + help="ChangeLog file") + + args = parser.parse_args() + + doit(args.changelog[0]) diff --git a/Docs/source/CNO_He_burn.png b/Docs/source/CNO_He_burn.png deleted file mode 100644 index 24c454e17c..0000000000 Binary files a/Docs/source/CNO_He_burn.png and /dev/null differ diff --git a/Docs/source/burn_cell.rst b/Docs/source/burn_cell.rst new file mode 100644 index 0000000000..7ccfdde5e8 --- /dev/null +++ b/Docs/source/burn_cell.rst @@ -0,0 +1,177 @@ +.. _sec:burn_cell: + +************* +``burn_cell`` +************* + +.. index:: burn_cell + +``burn_cell`` is a simple one-zone burn that will evolve a state with +a network for a specified amount of time. This can be used to +understand the timescales involved in a reaction sequence or to +determine the needed ODE tolerances. This is designed to work +with the Strang-split integration wrappers. The system that is evolved +has the form: + +.. math:: + + \begin{align*} + \frac{dX_k}{dt} &= \dot{\omega}_k(\rho, X_k, T) \\ + \frac{de}{dt} &= \epsilon(\rho, X_k, T) + \end{align*} + +with density held constant and the temperature found via the equation of state, +$T = T(\rho, X_k, e)$. + + +.. note:: + + Since the energy evolves due to the heat release (or loss) + from reactions, the temperature will change during the burn + (unless ``integrator.call_eos_in_rhs=0`` is set). + + +Getting Started +=============== + +The ``burn_cell`` code is located in +``Microphysics/unit_test/burn_cell``. An inputs file which sets the +default parameters for your choice of network is needed to run the +test. There are a number of inputs files in the unit test directory +already with a name list ``inputs_network``, where ``network`` +is the network you wish to use for your testing. These can be +used as a starting point for any explorations. + + +Setting the thermodynamics +-------------------------- + +The parameters that affect the thermodynamics are: + +* ``unit_test.density`` : the initial density + +* ``unit_test.temperature`` : the initial temperature + +* ``unit_test.small_temp`` : the low temperature cutoff used in the equation of state + +* ``unit_test.small_dens`` : the low density cutoff used in the equation of state + +The composition can be set either by setting each mass fraction explicitly via the +parameters, ``unit_test.X1``, ``unit_test.X2``, ..., +e.g.: + +:: + + unit_test.X1 = 0.5 + unit_test.X2 = 0.2 + unit_test.X3 = 0.2 + unit_test.X4 = 0.1 + +where parameters up to ``X35`` are available. If the values don't sum to ``1`` +initially, then the test will do a normalization. This normalization can be +disabled by setting: + +:: + + unit_test.skip_initial_normalization = 1 + + +Alternately, the composition can be set automatically by initializing all +of the mass fractions equally (to $1/N$, where $N$ is the number of species), +by setting: + +:: + + unit_test.uniform_xn = 1 + + +Controlling time +---------------- + +The test will run unit a time ``unit_test.tmax``, outputting the state +at regular intervals. The parameters controlling the output are: + +* ``unit_test.tmax`` : the end point of integration. + +* ``unit_test.tfirst`` : the first time interval to output. + +* ``unit_test.nsteps`` : the number of steps to divide the integration into, + logarithmically-spaced. + +If there is only a single step, ``unit_test.nsteps = 1``, then we integrate +from $[0, \mathrm{tmax}]$. + +If there are multiple steps, then the first output will be at a time +$\mathrm{tmax} / \mathrm{nsteps}$, and the steps will be +logarithmically-spaced afterwards. + + +Integration parameters +---------------------- + +The tolerances, choice of Jacobian, and other integration parameters +can be set via the usual Microphysics runtime parameters, e.g. +``integrator.atol_spec``. + + +Building and Running the Code +============================= + +The code can be built simply as: + +.. prompt:: bash + + make + +and the network and integrator can be changed using the normal +Microphysics build system parameters, e.g., + +.. prompt:: bash + + make NETWORK_DIR=aprox19 INTEGRATOR_DIR=rkc + +The build process will automatically create links in the build +directory to the EOS table and any reaction rate tables needed by your +choice of network. + + +.. important:: + + You need to do a ``make clean`` before rebuilding with a different + network or integrator. + + +To run the code, in the ``burn_cell`` directory run:: + + ./main3d.gnu.ex inputs + +where ``inputs`` is the name of your inputs file. + +Working with Output +=================== + +.. note:: + + For this part, we'll assume that the default ``aprox13`` and + ``VODE`` options were used for the network and integrator, and the + test was run with ``inputs.aprox13``. + +As the code runs, it will output to ``stdout`` details of the initial +and final state and the number of integration steps taken (along with whether +the burn was successful). The full history of the thermodynamic state will also be output to a file, +``state_over_time.txt``, with each line corresponding to one of the +``nsteps`` requested in the time integration. + +The script ``plot_burn_cell.py`` can be used to visualize the evolution: + +.. prompt:: bash + + python plot_burn_cell.py state_over_time.txt + +This will generate the following figure: + +.. figure:: state.png + :alt: An example of a plot output by the burn_cell unit test. + +Only the most abundant species are plotted. The number of species to plot and the +limits of $X$ can be set via runtime parameters (see ``python plot_burn_cell.py -h``). diff --git a/Docs/source/changes.rst b/Docs/source/changes.rst new file mode 100644 index 0000000000..6b5fae42c2 --- /dev/null +++ b/Docs/source/changes.rst @@ -0,0 +1,2 @@ +.. mdinclude:: ./changelog.md + diff --git a/Docs/source/citing.rst b/Docs/source/citing.rst new file mode 100644 index 0000000000..de4c2cc602 --- /dev/null +++ b/Docs/source/citing.rst @@ -0,0 +1,9 @@ +Citing Microphysics +=================== + +.. mdinclude:: ../../CITATION.md + +The bibtex for the latest version is included below (updated automatically): + +.. literalinclude:: ./zenodo.bibtex.txt + :language: bibtex diff --git a/Docs/source/conf.py b/Docs/source/conf.py index aeb914d8bc..e92103c066 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -1,5 +1,4 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- # # Microphysics documentation build configuration file, created by # sphinx-quickstart on Tue Oct 23 11:59:54 2018. @@ -23,7 +22,6 @@ import subprocess import sys -import breathe sys.path.insert(0, os.path.abspath('.')) @@ -31,8 +29,7 @@ def get_version(): prog = shlex.split("git describe --tags --abbrev=0") p0 = subprocess.Popen(prog, stdout=subprocess.PIPE, stderr=subprocess.PIPE) - stdout0, stderr0 = p0.communicate() - rc = p0.returncode + stdout0, _ = p0.communicate() stdout = stdout0.decode('utf-8') return stdout @@ -47,21 +44,22 @@ def get_version(): # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = ['sphinx.ext.autodoc', - 'sphinx.ext.mathjax', - 'sphinx_math_dollar', - 'sphinx.ext.viewcode', - 'sphinxcontrib.bibtex', - 'nbsphinx', - 'numpydoc', - 'IPython.sphinxext.ipython_console_highlighting', - 'sphinx.ext.githubpages', - 'sphinx_copybutton', - 'sphinx-prompt', - 'sphinx_rtd_theme', - 'breathe'] + 'sphinx.ext.mathjax', + 'sphinx_math_dollar', + 'sphinx.ext.viewcode', + 'sphinxcontrib.bibtex', + 'nbsphinx', + 'numpydoc', + 'IPython.sphinxext.ipython_console_highlighting', + 'sphinx.ext.githubpages', + 'sphinx_copybutton', + 'sphinx-prompt', + 'sphinx_mdinclude', + 'sphinx_rtd_theme', + 'breathe'] breathe_projects = { - "microphysics":"../doxy_files/xml", + "microphysics": "../doxy_files/xml", } breathe_default_project = "microphysics" @@ -109,7 +107,7 @@ def get_version(): # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This patterns also effect to html_static_path and html_extra_path -exclude_patterns = [] +exclude_patterns = ["changelog.md"] # The name of the Pygments (syntax highlighting) style to use. pygments_style = 'sphinx' @@ -128,7 +126,7 @@ def get_version(): mathjax3_config["tex"]["macros"] = {} -with open('mathsymbols.tex', 'r') as f: +with open('mathsymbols.tex') as f: for line in f: macros = re.findall(r'\\newcommand{\\(.*?)}(\[(\d)\])?{(.+)}', line) for macro in macros: diff --git a/Docs/source/data_structures.rst b/Docs/source/data_structures.rst index 8a4c4ccc3c..1e130d4b7d 100644 --- a/Docs/source/data_structures.rst +++ b/Docs/source/data_structures.rst @@ -43,7 +43,7 @@ network, and they will come in through the ``network_properties.H`` header file. There is a lot more information that can be saved here, such as the partial derivatives of the thermodynamic state variables with respect to each other. To see a complete list, examine the ``eos_type.H`` -file: ``Castro/Microphysics/interfaces/eos_type.H``. +file. Networks ======== @@ -111,10 +111,30 @@ to access the different components of the state: * ``neqs`` : the total number of variables we are integrating. - It is assumed that the first ``nspec`` are the species. + It is assumed that the first ``nspec`` are the species. * ``net_ienuc`` : the index of the specific internal energy in the solution vector +For :ref:`sdc-evolution`, it also defines integer indices for the +``burn_t y[]`` array: + +* ``SFS``: the first species + +* ``SEINT`` the energy + +and then a number of components that are not evolved: + +* ``SRHO`` density + +* ``SMX``, ``SMY``, ``SMZ`` : the momenta + +* ``SEDEN`` : the total energy density + +* ``SFX`` : the first auxiliary quantity + +with the total number of state variables ``SVAR`` and the number of evolved +variables ``SVAR_EVOLVE``. + Integrators =========== diff --git a/Docs/source/eos.rst b/Docs/source/eos.rst index dee9ee809f..0c8778505c 100644 --- a/Docs/source/eos.rst +++ b/Docs/source/eos.rst @@ -2,206 +2,18 @@ Equations of State ****************** -In this chapter on equations of state, we list the EOS routines -available for your use, and then we describe the correct structure of -an EOS module in case you want to build your own. +The general interface to the equation of state is: -Available Equations of State -============================ - -.. index:: eos_t - -The following equations of state are available in Microphysics. -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} +.. code:: c++ -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 + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void eos (const I input, T& state, bool use_raw_inputs = false) -.. math:: +where ``input`` specifies which thermodynamic quantities are taken as +the input and ``state`` is a C++ struct that holds all of the +thermodynamic information. - \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`. - -We have modified this EOS a bit to fit within the context of our -codes. The vectorization is explicitly thread-safe for use with OpenMP -and OpenACC. In addition, we have added the ability to perform a -Newton-Raphson iteration so that if we call the EOS with density and -energy (say), then we will 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. However, there is a choice of whether -to update the energy to match this temperature, respecting -thermodynamic consistency, or to leave the energy alone, respecting -energy conservation. This is controlled through the -``eos.eos_input_is_constant`` parameter in your inputs file. - -We thank Frank Timmes for permitting us to modify his code and -publicly release it in this repository. - -stellarcollapse ---------------- - -``stellarcollapse`` is the equation of state module provided -on http://stellarcollapse.org. It is designed -to be used for core-collapse supernovae and is compatible with a -large number of equations of state that are designed to describe -matter near nuclear density. You will need to download an -appropriate interpolation table from that site to use this. Interface and Modes =================== @@ -248,16 +60,28 @@ The *eos_type* passed in is one of * ``chem_eos_t`` : adds some quantities needed for the primordial chemistry EOS and explicitly does not include the mass fractions. -In general, you should use the type that has the smallest set of -information needed, since we optimize out needless quantities at -compile type (via C++ templating) for ``eos_re_t`` and ``eos_rep_t``. +.. tip:: -.. note:: + The EOS implementations make heavy use of templating to + "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. + +.. 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. + + +Composition +=========== - All of these modes require composition as an input. Usually this is - via the set of mass fractions, ``eos_t.xn[]``, but if ``USE_AUX_THERMO`` - is set to ``TRUE``, then we instead use the auxiliary quantities - stored in ``eos_t.aux[]``. +All input modes for ``eos()`` require a composition. Usually this is +via the set of mass fractions, ``eos_t xn[]``, but if ``USE_AUX_THERMO`` +is set to ``TRUE``, then we instead use the auxiliary quantities +stored in ``eos_t.aux[]``. .. _aux_eos_comp: @@ -266,16 +90,22 @@ 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. -The auxiliary variables are +For ``eos_t eos_state``, the auxiliary variables are + -* ``eos_state.aux[iye]`` : electron fraction, defined as +* ``eos_state.aux[AuxZero::iye]`` : electron fraction, defined as .. math:: Y_e = \sum_k \frac{X_k Z_k}{A_k} -* ``eos_state.aux[iabar]`` : the average mass of the nuclei, :math:`\bar{A}`, defined as: +* ``eos_state.aux[AuxZero::iabar]`` : the average mass of the nuclei, :math:`\bar{A}`, defined as: .. math:: @@ -283,7 +113,7 @@ The auxiliary variables are In many stellar evolutions texts, this would be written as :math:`\mu_I`. -* ``eos_state.aux[ibea]`` : the binding energy per nucleon (units of +* ``eos_state.aux[AuxZero::ibea]`` : the binding energy per nucleon (units of MeV), defined as .. math:: @@ -292,11 +122,17 @@ The auxiliary variables are where :math:`B_k` is the binding energy of nucleus :math:`k` -Given a composition of mass fractions, the function -``set_aux_comp_from_X(state_t& state)`` will initialize these -auxiliary quantities. +Given a composition of mass fractions, the function: + +.. code:: c++ + + template + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void set_aux_comp_from_X(state_t& state) + +will initialize the auxiliary data. -The equation of state also needs :math:`\bar{Z}` which is easily computed as +Many equations of state also need :math:`\bar{Z}` which is easily computed as .. math:: @@ -315,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 @@ -332,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_cell.rst b/Docs/source/eos_cell.rst new file mode 100644 index 0000000000..7399884eaa --- /dev/null +++ b/Docs/source/eos_cell.rst @@ -0,0 +1,76 @@ +.. _sec:eos_cell: + +************ +``eos_cell`` +************ + +.. index:: eos_cell + +``eos_cell`` simply calls the equation of state on an input density, temperature, +and composition, and then outputs the full thermodynamic state. This is mainly +used to understand the thermodynamics one might encounter in a simulation +when using a particular EOS. + +Getting Started +=============== + +The ``eos_cell`` code is located in +``Microphysics/unit_test/eos_cell``. An inputs file which sets the +default parameters for your thermodynamic state is needed to run the +test. + +Setting the thermodynamics +-------------------------- + +The parameters that affect the thermodynamics are: + +* ``unit_test.density`` : the initial density + +* ``unit_test.temperature`` : the initial temperature + +* ``unit_test.small_temp`` : the low temperature cutoff used in the equation of state + +* ``unit_test.small_dens`` : the low density cutoff used in the equation of state + +The composition can be set in the same way as in ``burn_cell``, either +by setting each mass fraction explicitly via the parameters, +``unit_test.X1``, ``unit_test.X2``, ..., or forcing them to be all +equal via ``unit_test.uniform_xn=1``. + + +Building and Running the Code +============================= + +The code can be built simply as: + +.. prompt:: bash + + make + +.. note:: + + Even though there are no reactions, a network is still required, + and can be set via the ``NETWORK_DIR`` build variable. By default, + the ``aprox13`` network is used. + + The network choice serves only to set the composition, and a + ``general_null`` network may also be used. + +The build process will automatically create links in the build +directory to any required EOS table. + +To run the code, in the ``eos_cell`` directory run:: + + ./main3d.gnu.ex inputs_eos + +where ``inputs_eos`` is the provided inputs file. You may edit the +thermodynamic state in that file prior to running. + + +Output +====== + +All output is directed to ``stdout`` and simply lists the entries in the +full ``eos_t`` datatype. Example output is shown below: + +.. literalinclude:: ../../unit_test/eos_cell/ci-benchmarks/eos_helmholtz.out diff --git a/Docs/source/eos_implementations.rst b/Docs/source/eos_implementations.rst new file mode 100644 index 0000000000..f8ea10538e --- /dev/null +++ b/Docs/source/eos_implementations.rst @@ -0,0 +1,294 @@ +**************************** +Available Equations of State +**************************** + +.. index:: eos_t, EOS_DIR + +The following equations of state are available in Microphysics. + +.. note:: + + The EOS is chosen at compile-time via the ``EOS_DIR`` make + variable. + +.. note:: + + Except where noted, each of these EOSs will provide the full + thermodynamic data (including all derivatives) in the ``eos_t`` + type. + + +``breakout`` +============ + +The ``breakout`` EOS is essentially the same as ``gamma_law``, but it gets +its composition information from the auxiliary data. In particular, +it expects an auxiliary quantity named ``invmu`` which is the inverse +of the mean molecular weight: + +.. math:: + + \frac{1}{\mu} = \sum_k \frac{Z_k X_k}{A_k} + +The ``general_null`` network provides this when used with the ``breakout.net`` +network inputs. + +``gamma_law`` +============= + +.. index:: eos.eos_gamma, eos.assume_neutral + +``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. + + +``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`` : do we include Coulomb corrections? This + is enabled by default. Coulomb corrections can cause problems in + some regimes, because the implementation in ``helmholtz`` doesn't + have the correct asymptotic behavior and can lead to negative + pressures or energies. + +* ``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`` : these are the tolerances + for temperature and density used by the Newton solver when + inverting the EOS. + +* ``eos.prad_limiter_rho_c``, ``eos.prad_limiter_delta_rho`` : by + default, radiation pressure is included in the optically-thick, LTE + limit (with $p_\gamma = (1/3)a T^4$). At low densities, this can + cause issues, leading to an artificially high soundspeed dominated + by radiation when, in fact, we should be optically thin. These + parameters allow us turn off the radiation component smoothly, + starting at a density ``eos.prad_limiter_rho_c`` and transitioning + via a $\tanh$ profile to zero over a scale + ``eos.prad_limiter_delta_rho``. + +We thank Frank Timmes for permitting us to modify his code and +publicly release it in this repository. + +``metal_chem`` +============== + +This is a multi-gamma equation of state for metal ISM chemistry. + +``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: + +.. index:: eos.eos_gamma_default + +- ``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. + + + +``polytrope`` +============= + +.. index:: eos.polytrope_K, eos.polytrope_gamma, eos.polytrope_type, eos.polytrope_mu_e + +``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. + + +``primordial_chem`` +=================== + +This is a version of the multi-gamma equation of state that models primordial chemistry. + +``rad_power_law`` +================= + +This is an artificial equation of state for radiation transport test problems. It uses +a parameterization of the specific heat at constant volume: + +.. math:: + + c_v = A \rho^m T^{-n} + +and energy: + +.. math:: + + e = \frac{A}{1 - n} \rho^m T^{1-n} + +where the runtime parameters provide the constants: + +* ``eos.eos_const_c_v`` $= A$ + +* ``eos.eos_c_v_exp_m`` $= m$ + +* ``eos.eos_c_v_exp_n`` $= n$ + + +``tillotson`` +============= + +This is an equation of state for hypervelocity impacts based on :cite:`tillotson:1962`. + + +``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. diff --git a/Docs/source/getting_started.rst b/Docs/source/getting_started.rst index 02aab8f7a8..ab6482767a 100644 --- a/Docs/source/getting_started.rst +++ b/Docs/source/getting_started.rst @@ -2,6 +2,29 @@ Getting Started *************** +Requirements +============ + +Microphysics requires + +* A C++17 or later compilers +* AMReX (https://github.com/amrex-codes/amrex) +* python (>= 3.10) +* GNU make + +optional dependencies are: + +* CUDA (>= 11) +* ROCm (>= 6.3.1 --- earlier versions have register allocation bugs) + +Microphysics is meant to be compiled into an application code as part +of its build process, with the network, EOS, integrator, and more +picked at compile time. As such, there is not a single library that +can be built and linked against. + +Below we describe how to use Microphysics in a "standalone" fashion, +using the provided unit tests, and as part of an application code. + Standalone ========== @@ -35,27 +58,42 @@ one-zone burn. In ``Microphysics/`` do: cd unit_test/burn_cell make -This will create an executable called ``main3d.gnu.ex``. Then you can run it as: +This will create an executable called ``main3d.gnu.ex``. +By default, the test is built with the 13-isotope ``aprox13`` network, +``helmholtz`` EOS, and VODE integrator. + + +Then you can run it as: .. prompt:: bash - ./main3d.gnu.ex inputs_aprox21 + ./main3d.gnu.ex inputs_aprox13 -By default, the test is built with the 21-isotope ``aprox21`` network. -Here ``inputs_aprox21`` is the inputs file that sets options. +Here ``inputs_aprox13`` is the inputs file that sets options. +This will output information about the starting and final state to the +terminal and produce a file ``state_over_time.txt`` that contains the +thermodynamic history at different points in time. +.. note:: + + See the :ref:`sec:burn_cell` documentation for more details on this + unit test and how to visualize the output. Running with AMReX Application Code =================================== .. index:: MICROPHYSICS_HOME -Getting started with Microphysics using either CASTRO or MAESTROeX is +Getting started with Microphysics using either `CASTRO +`_ or `MAESTROeX +`_ is straightforward. Because the modules here are already in a format that the AMReX codes understand, you only need to provide to the code calling these routines their location on your system. The code will do -the rest. To do so, define the ``MICROPHYSICS_HOME`` environment +the rest. + +First we need to define the ``MICROPHYSICS_HOME`` environment variable, either at a command line or (if you use the bash shell) through your ``~/.bashrc``, e.g.: @@ -70,6 +108,4 @@ rely on the Microphysics ``Make.Microphysics_extern`` makefile stub source to the build. All of the interfaces that these codes use are found in ``Microphysics/interfaces/``. -Other codes can use Microphysics in the same fashion. Unit tests in -``Microphysics/unit_test/`` provide some examples of using the -interfaces. +Other AMReX-based codes can use Microphysics in the same fashion. diff --git a/Docs/source/index.rst b/Docs/source/index.rst index d94f7b1b8b..a366bd78c9 100644 --- a/Docs/source/index.rst +++ b/Docs/source/index.rst @@ -14,11 +14,14 @@ The original design was to support codes based on the `AMReX `_ adaptive mesh refinement library :cite:`amrex_joss`, including `CASTRO `_ :cite:`castro_I`, `MAESTROeX -`_ :cite:`maestroex`, and Quokka :cite:`quokka`. These all have a +`_ :cite:`maestroex`, and, later, Quokka :cite:`quokka`. These all have a consistent interface and the separate Microphysics repository allows them to share the same equation of state, reaction networks, and more. -Later, Microphysics was adopted by the `Quokka `_ -simulation code. + +Microphysics is written in C++ as a (mostly) header-only library, making +extensive use of templating and C++ `constexpr` features to be performant +on both CPU and GPU architectures. It is compatible with both NVIDIA CUDA +and AMD HIP/ROCm. While there are a number of unit tests that exercise the functionality, Microphysics is primarily intended to be used along with another simulation @@ -26,11 +29,13 @@ code. At the moment, the interfaces and build stubs are compatible with the AMReX codes and use the AMReX build system. -A number of the routines contained here we authored by other people. -We bundle them here with permission, usually changing the interfaces -to be compatible with our standardized interface. We in particular -thank Frank Timmes for numerous reaction networks and his equation -of state routines. +.. note:: + + A number of the routines contained here we authored by other people. + We bundle them here with permission, usually changing the interfaces + to be compatible with our standardized interface. We in particular + thank Frank Timmes for numerous reaction networks and his equation + of state routines. .. toctree:: @@ -50,6 +55,7 @@ of state routines. :hidden: eos + eos_implementations transport .. toctree:: @@ -94,6 +100,8 @@ of state routines. :caption: References :hidden: + citing + changes zreferences .. toctree:: diff --git a/Docs/source/jac_cell.rst b/Docs/source/jac_cell.rst new file mode 100644 index 0000000000..51db309555 --- /dev/null +++ b/Docs/source/jac_cell.rst @@ -0,0 +1,98 @@ +.. _sec:jac_cell: + +************ +``jac_cell`` +************ + +.. index:: jac_cell, numerical_jacobian.H + +``jac_cell`` is used to test the accuracy of an analytic Jacobian +provided by a network, but comparing to a finite-difference +approximation. Given a thermodynamic state, ``jac_cell`` evaluates +the analytic Jacobian as well as a numerical approximation to it, and +then outputs (to ``stdout``) the Jacobian elements side-by-side for +each method of computing the Jacobian. + +.. note:: + + Some integrators, like VODE, have their own numerical Jacobian + routines. This test uses the implementation in + ``integration/util/numerical_jacobian.H``, which closely follows the ideas + of the VODE numerical approximation, as described in :cite:`lsode`. + + +Getting Started +=============== + +The ``jac_cell`` code is located in +``Microphysics/unit_test/jac_cell``. An inputs file which sets the +default parameters for your thermodynamic state is needed to run the +test. + +Setting the thermodynamics +-------------------------- + +The parameters that affect the thermodynamics are: + +* ``unit_test.density`` : the initial density + +* ``unit_test.temperature`` : the initial temperature + +* ``unit_test.small_temp`` : the low temperature cutoff used in the equation of state + +* ``unit_test.small_dens`` : the low density cutoff used in the equation of state + +While the mass fractions can be set individually (using +``unit_test.X1``, ``unit_test.X2``, ...), it is recommended to use +``unit_test.uniform_xn=1`` to initialize all the mass fractions to be +equal. + + +Building and Running the Code +============================= + +The code can be built simply as: + +.. prompt:: bash + + make + +.. note:: + + By default, this will build with the ``aprox13`` network, which + uses the C++ templating method (:ref:`sec:templated_rhs`) of + building the analytic Jacobian at compile-time. The network can be + changed via the build parameter ``NETWORK_DIR``. + +.. important:: + + The screening routines provided by Microphysics do not return the + derivative of the screening factor with respect to composition. This + means that the analytic Jacobian provided by the network will ignore + this contribution, but the numerical Jacobian will implicitly include + it. + + To compare just the network Jacobian elements, it is suggested that + you build with ``SCREEN_METHOD=null``. + + +The build process will automatically create links in the build +directory to any required EOS table. + +To run the code, in the ``jac_cell`` directory run:: + + ./main3d.gnu.ex inputs + +where ``inputs`` is the provided inputs file. You may edit the +thermodynamic state in that file prior to running. + + +Output +====== + +All output is directed to ``stdout``. Example output is shown below +(this output includes screening, highlighting the difference that the +screening composition derivatives make): + + +.. literalinclude:: ../../unit_test/jac_cell/ci-benchmarks/jac_cell_aprox13.out diff --git a/Docs/source/networks.rst b/Docs/source/networks.rst index 908a6ee200..4411495b70 100644 --- a/Docs/source/networks.rst +++ b/Docs/source/networks.rst @@ -2,40 +2,35 @@ Available Reaction Networks *************************** +A network defines the composition, which is needed by the equation +of state and transport coefficient routines. Even if there are no +reactions taking place, a network still needs to be defined, so +Microphysics knows the properties of the fluid. -``iso7``, ``aprox13``, ``aprox19``, and ``aprox21`` -=================================================== +.. index:: NETWORK_DIR -These are alpha-chains (with some other nuclei) from Frank Timmes. -These networks share common rates (from ``Microphysics/rates``), -plasma neutrino loses (from ``Microphysics/neutrinos``), and -electron screening (from ``Microphysics/screening``). +.. note:: -Energy generation. ------------------- + The network is set at compile-time via the ``NETWORK_DIR`` + make variable. -These networks store the total binding energy of the nucleus in MeV as -``bion(:)``. They then compute the mass of each nucleus in grams as: +.. tip:: -.. math:: M_k = (A_k - Z_k) m_n + Z_k (m_p + m_e) - B_k + If reactions can be ignored, then the ``general_null`` network can + be used --- this simply defines a composition with no reactions. -where :math:`m_n`, :math:`m_p`, and :math:`m_e` are the neutron, proton, and electron -masses, :math:`A_k` and :math:`Z_k` are the atomic weight and number, and :math:`B_k` -is the binding energy of the nucleus (converted to grams). :math:`M_k` -is stored as ``mion(:)`` in the network. - -The energy release per gram is converted from the rates as: - -.. math:: \epsilon = -N_A c^2 \sum_k \frac{dY_k}{dt} M_k - \epsilon_\nu - -where :math:`N_A` is Avogadro’s number (to convert this to “per gram”) -and :math:`\edotnu` is the neutrino loss term (see :ref:`neutrino_loss`). +.. note:: + Many of the networks here are generated using `pynucastro + `_ using the ``AmrexAstroCxxNetwork`` + class. ``general_null`` ================ -``general_null`` is a bare interface for a nuclear reaction network -- +.. index:: general_null + +``general_null`` is a bare interface for a nuclear reaction network --- no reactions are enabled. The data in the network is defined at compile type by specifying an inputs file. For example, @@ -75,36 +70,268 @@ The name of the inputs file by one of two make variables: GENERAL_NET_INPUTS := /path/to/file/triple_alpha_plus_o.net +.. index:: network_properties.H + At compile time, the "`.net`" file is parsed and a network header ``network_properties.H`` is written using the python script ``write_network.py``. The make rule for this is contained in -``Make.package``. +``Microphysics/networks/Make.package``. + + +``iso7``, ``aprox13``, ``aprox19``, and ``aprox21`` +=================================================== + +These are alpha-chains (with some other nuclei) based on the `original +Fortran networks from Frank Timmes +`_. These +networks share common rates from ``Microphysics/rates`` and are +implemented using the templated C++ network infrastructure. + +These networks approximate a lot of the links, in particular, +combining $(\alpha, p)(p, \gamma)$ and $(\alpha, \gamma)$ into a +single effective rate. + +The available networks are: + +* ``iso7`` : contains $\isotm{He}{4}$, $\isotm{C}{12}$, + $\isotm{O}{16}$, $\isotm{Ne}{20}$, $\isotm{Mg}{24}$, $\isotm{Si}{28}$, + $\isotm{Ni}{56}$ and is based on :cite:`iso7`. + +* ``aprox13`` : adds $\isotm{S}{32}$, $\isotm{Ar}{36}$, $\isotm{Ca}{40}$, $\isotm{Ti}{44}$, $\isotm{Cr}{48}$, $\isotm{Fe}{52}$ + +* ``aprox19`` : adds $\isotm{H}{1}$, $\isotm{He}{3}$, $\isotm{N}{14}$, $\isotm{Fe}{54}$, + $\mathrm{p}$, $\mathrm{n}$. Here, $\mathrm{p}$ participates only in the photodisintegration rates at high mass number, and is distinct from $\isotm{H}{1}$. + +* ``aprox21`` : adds $\isotm{Cr}{56}$, $\isotm{Fe}{56}$. This is designed to reach + a lower $Y_e$ than the other networks, for use in massive star simulations. Note + that the link to $\isotm{Cr}{56}$ is greatly approximated. + + +These networks store the total binding energy of the nucleus in MeV as +``bion(:)``. They then compute the mass of each nucleus in grams as: + +.. math:: M_k = (A_k - Z_k) m_n + Z_k (m_p + m_e) - B_k + +where :math:`m_n`, :math:`m_p`, and :math:`m_e` are the neutron, proton, and electron +masses, :math:`A_k` and :math:`Z_k` are the atomic weight and number, and :math:`B_k` +is the binding energy of the nucleus (converted to grams). :math:`M_k` +is stored as ``mion(:)`` in the network. + +The energy release per gram is converted from the rates as: + +.. math:: \epsilon = -N_A c^2 \sum_k \frac{dY_k}{dt} M_k - \epsilon_\nu + +where :math:`N_A` is Avogadro’s number (to convert this to “per gram”) +and :math:`\epsilon_\nu` is the neutrino loss term (see :ref:`neutrino_loss`). + + ``CNO_extras`` ============== -This network replicates the popular [MESA "cno_extras" -network](https://docs.mesastar.org/en/latest/net/nets.html) which is +This network replicates the popular `MESA "cno_extras" +network `_ which is meant to study hot-CNO burning and the start of the breakout from CNO -burning. - -We add ${}^{56}\mathrm{Fe}$ as an inert nucleus to allow this to be -used for X-ray burst simulations. +burning. This network is managed by pynucastro. .. figure:: cno_extras_hide_alpha.png :align: center +.. note:: + + We add ${}^{56}\mathrm{Fe}$ as an inert nucleus to allow this to be + used for X-ray burst simulations (not shown in the network diagram + above). + + +nova networks +============= + +The ``nova`` and ``nova2`` networks both are intended for modeling classical novae. + + +* ``nova`` focuses just on CNO/hot-CNO: + + .. figure:: nova.png + :align: center + +* ``nova2`` expands ``nova`` by adding the pp-chain nuclei: + + .. figure:: nova2.png + :align: center + + +He-burning networks +=================== + +This is a collection of networks meant to model He burning. The are inspired by the +"aprox"-family of networks, but contain more nuclei/rates, and are managed by +pynucastro. + +One feature of these networks is that they include a bypass rate for +:math:`\isotm{C}{12}(\alpha, \gamma)\isotm{O}{16}` discussed in +:cite:`ShenBildsten`. This is appropriate for explosive He burning. +That paper discusses the sequences: + +* :math:`\isotm{C}{14}(\alpha, \gamma)\isotm{O}{18}(\alpha, + \gamma)\isotm{Ne}{22}` at high temperatures (T > 1 GK). We don't + consider this. + +* :math:`\isotm{N}{14}(\alpha, \gamma)\isotm{F}{18}(\alpha, + p)\isotm{Ne}{21}` is the one they consider important, since it produces + protons that are then available for :math:`\isotm{C}{12}(p, + \gamma)\isotm{N}{13}(\alpha, p)\isotm{O}{16}`. + +This leaves :math:`\isotm{Ne}{21}` as an endpoint, which we connect to +the other nuclei by including :math:`\isotm{Na}{22}`. + +For the :math:`\isotm{C}{12} + \isotm{C}{12}`, :math:`\isotm{C}{12} + +\isotm{O}{16}`, and :math:`\isotm{O}{16} + \isotm{O}{16}` rates, we +also need to include: + +* :math:`\isotm{C}{12}(\isotm{C}{12},n)\isotm{Mg}{23}(n,\gamma)\isotm{Mg}{24}` + +* :math:`\isotm{O}{16}(\isotm{O}{16}, n)\isotm{S}{31}(n, \gamma)\isotm{S}{32}` + +* :math:`\isotm{O}{16}(\isotm{C}{12}, n)\isotm{Si}{27}(n, \gamma)\isotm{Si}{28}` + +Since the neutron captures on those +intermediate nuclei are so fast, we leave those out and take the +forward rate to just be the first rate. We do not include reverse +rates for these processes. + +These networks also combine some of the +:math:`A(\alpha,p)X(p,\gamma)B` links with :math:`A(\alpha,\gamma)B`, +allowing us to drop the intermediate nucleus :math:`X`. Some will +approximate $A(n,\gamma)X(n,\gamma)B$ into an effective +$A(nn,\gamma)B$ rate (double-neutron capture). + +The networks are named with a descriptive name, the number of nuclei, +and the letter ``a`` if they approximate $(\alpha, p)(p,\gamma)$, +the letter ``n`` if they approximate double-neutron capture, and the +letter ``p`` if they split the protons into two groups (one for +photo-disintegration). + + +``he-burn-18a`` +--------------- + +.. note:: + + This network was previously called ``subch_base``. + +This is the simplest network and is similar to ``aprox13``, but includes +a better description of $\isotm{C}{12}$ and $\isotm{O}{16}$ burning, as +well as the bypass rate for $\isotm{C}{12}(\alpha,\gamma)\isotm{O}{16}$. + +It has the following features / simplifications: + +* $\isotm{Cl}{35}$, $\isotm{K}{39}$, $\isotm{Sc}{43}$, + $\isotm{V}{47}$, $\isotm{Mn}{51}$, and $\isotm{Co}{55}$ are approximated + out of the $(\alpha, p)(p, \gamma)$ links. + +* The nuclei :math:`\isotm{N}{14}`, :math:`\isotm{F}{18}`, + :math:`\isotm{Ne}{21}`, and :math:`\isotm{Na}{22}` are not included. + This means that we do not capture the :math:`\isotm{N}{14}(\alpha, + \gamma)\isotm{F}{18}(\alpha, p)\isotm{Ne}{21}` rate sequence. + +* The reverse rates of :math:`\isotm{C}{12}+\isotm{C}{12}`, + :math:`\isotm{C}{12}+\isotm{O}{16}`, :math:`\isotm{O}{16}+\isotm{O}{16}` are + neglected since they're not present in the original aprox13 network + +* The :math:`\isotm{C}{12}+\isotm{Ne}{20}` rate is removed + +* The :math:`(\alpha, \gamma)` links between :math:`\isotm{Na}{23}`, + :math:`\isotm{Al}{27}` and between :math:`\isotm{Al}{27}` and + :math:`\isotm{P}{31}` are removed, since they're not in the + original aprox13 network. + +The network appears as: + +.. figure:: ../../networks/he-burn/he-burn-18a/he-burn-18a.png + :align: center + +The nuclei in gray are those that have been approximated about, but the links +are effectively accounted for in the approximate rates. + +There are 2 runtime parameters that can be used +to disable rates: + +* ``network.disable_p_c12__n13`` : if set to ``1``, then the rate + :math:`\isotm{C}{12}(p,\gamma)\isotm{N}{13}` and its inverse are + disabled. + +* ``network.disable_he4_n13__p_o16`` : if set to ``1``, then the rate + :math:`\isotm{N}{13}(\alpha,p)\isotm{O}{16}` and its inverse are + disabled. + +Together, these parameters allow us to turn off the sequence +:math:`\isotm{C}{12}(p,\gamma)\isotm{N}{13}(\alpha, p)\isotm{O}{16}` that +acts as a bypass for :math:`\isotm{C}{12}(\alpha, \gamma)\isotm{O}{16}`. + +``he-burn-22a`` +--------------- + +.. note:: + + This network was previously called ``subch_simple``. + + +This builds on ``he-burn-18a`` by including the +:math:`\isotm{N}{14}(\alpha, \gamma)\isotm{F}{18}(\alpha, +p)\isotm{Ne}{21}` rate sequence, which allows an enhancement to the +:math:`\isotm{C}{12}(p, \gamma)\isotm{N}{13}(\alpha, p)\isotm{O}{16}` +rate due to the additional proton release. + +.. figure:: ../../networks/he-burn/he-burn-22a/he-burn-22a.png + :align: center + +.. warning:: Due to inclusion of the rate sequence, + ${}^{14}\mathrm{N}(\alpha, \gamma){}^{18}\mathrm{F}(\alpha, + p){}^{21}\mathrm{Ne}$, there is an artificial end-point at + ${}^{22}\mathrm{Na}$. + +Like ``he-burn-18a``, there are 2 runtime parameters that can disable +the rates for the $\isotm{C}{12}(p,\gamma)\isotm{N}{13}(\alpha, +p)\isotm{O}{16}$ sequence. + +``he-burn-31anp`` +----------------- + +This builds on ``he-burn-22a`` by adding some iron-peak nuclei. It no longer +approximates out $\isotm{Mn}{51}$ or $\isotm{Co}{55}$, and includes approximations +to double-neutron capture. Finally, it splits the protons into two groups, +with those participating in reactions with mass numbers > 48 treated as a separate +group (for photo-disintegration reactions). + +The iron group here resembles ``aprox21``, but has the addition of stable $\isotm{Ni}{58}$ +and doesn't include the approximation to $\isotm{Cr}{56}$. + +.. figure:: ../../networks/he-burn/he-burn-31anp/he-burn-31anp.png + :align: center + + +``he-burn-36a`` +--------------- + +This has the most complete iron-group, with nuclei up to $\isotm{Zn}{60}$ and no approximations +to the neutron captures. This network can be quite slow. + +.. figure:: ../../networks/he-burn/he-burn-36a/he-burn-36a.png + :align: center + + ``CNO_He_burn`` -=============== +--------------- This network is meant to study explosive H and He burning. It combines the ``CNO_extras`` network (with the exception of the inert ${}^{56}\mathrm{Fe}$ -with the ``subch_simple`` network. This allows it to capture hot-CNO and +with the ``he-burn-22a`` network. This allows it to capture hot-CNO and He burning. -.. figure:: CNO_He_burn.png +.. figure:: ../../networks/he-burn/cno-he-burn-33a/cno-he-burn-33a.png :align: center ``ECSN`` @@ -116,9 +343,15 @@ It includes various weak rates that are important to this process. .. figure:: ECSN.png :align: center +C-ignition networks +=================== + +There are a number of networks that have been developed for exploring +carbon burning in near-Chandrasekhar mass which dwarfs. + ``ignition_chamulak`` -===================== +--------------------- This network was introduced in our paper on convection in white dwarfs as a model of Type Ia supernovae :cite:`wdconvect`. It models @@ -127,11 +360,6 @@ and captures the effects of a much larger network by setting the ash state and energetics to the values suggested in :cite:`chamulak:2008`. -.. _energy-generation.-1: - -Energy generation. ------------------- - The binding energy, :math:`q`, in this network is interpolated based on the density. It is stored as the binding energy (ergs/g) *per nucleon*, with a sign convention that @@ -142,7 +370,7 @@ binding energies are negative. The energy generation rate is then: (this is positive since both :math:`q` and :math:`dY/dt` are negative) ``ignition_reaclib`` -==================== +-------------------- This contains several networks designed to model C burning in WDs. They include: @@ -158,7 +386,7 @@ This contains several networks designed to model C burning in WDs. They include ``ignition_simple`` -=================== +------------------- This is the original network used in our white dwarf convection studies :cite:`lowMach4`. It includes a single-step @@ -178,22 +406,6 @@ of (Graboske 1973) for weak screening and the work of (Alastuey 1978 and Itoh 1979) for strong screening. -nova networks -============= - -The ``nova`` and ``nova2`` networks both are intended for modeling classical novae. - - -* ``nova`` focuses just on CNO/hot-CNO: - - .. figure:: nova.png - :align: center - -* ``nova2`` expands ``nova`` by adding the pp-chain nuclei: - - .. figure:: nova2.png - :align: center - ``powerlaw`` ============ @@ -266,113 +478,3 @@ This network was used for the X-ray burst studies in This is a 2 reaction network for helium burning, capturing the :math:`3`-:math:`\alpha` reaction and :math:`\isotm{C}{12}(\alpha,\gamma)\isotm{O}{16}`. Additionally, :math:`^{56}\mathrm{Fe}` is included as an inert species. - - -subch networks -============== - -The subch networks recreate an ``aprox13`` -alpha-chain + including a bypass rate for :math:`\isotm{C}{12}(\alpha, -\gamma)\isotm{O}{16}` discussed in :cite:`ShenBildsten`. This is appropriate -for explosive He burning. - -:cite:`ShenBildsten` discuss the sequences: - -* :math:`\isotm{C}{14}(\alpha, \gamma)\isotm{O}{18}(\alpha, - \gamma)\isotm{Ne}{22}` at high temperatures (T > 1 GK). We don't - consider this. - -* :math:`\isotm{N}{14}(\alpha, \gamma)\isotm{F}{18}(\alpha, - p)\isotm{Ne}{21}` is the one they consider important, since it produces - protons that are then available for :math:`\isotm{C}{12}(p, - \gamma)\isotm{N}{13}(\alpha, p)\isotm{O}{16}`. - -This leaves :math:`\isotm{Ne}{21}` as an endpoint, which we connect to -the other nuclei by including :math:`\isotm{Na}{22}`. - -For the :math:`\isotm{C}{12} + \isotm{C}{12}`, :math:`\isotm{C}{12} + -\isotm{O}{16}`, and :math:`\isotm{O}{16} + \isotm{O}{16}` rates, we -also need to include: - -* :math:`\isotm{C}{12}(\isotm{C}{12},n)\isotm{Mg}{23}(n,\gamma)\isotm{Mg}{24}` - -* :math:`\isotm{O}{16}(\isotm{O}{16}, n)\isotm{S}{31}(n, \gamma)\isotm{S}{32}` - -* :math:`\isotm{O}{16}(\isotm{C}{12}, n)\isotm{Si}{27}(n, \gamma)\isotm{Si}{28}` - -Since the neutron captures on those -intermediate nuclei are so fast, we leave those out and take the -forward rate to just be the first rate. We do not include reverse -rates for these processes. - - -``subch_simple`` ----------------- - -``subch_simple`` uses the ideas above but approximates some -of the rates by -combining some of the :math:`A(\alpha,p)X(p,\gamma)B` links with -:math:`A(\alpha,\gamma)B`, allowing us to drop the intermediate -nucleus :math:`X`. We do this for :math:`\isotm{Cl}{35}`, -:math:`\isotm{K}{39}`, :math:`\isotm{Sc}{43}`, :math:`\isotm{V}{47}`, -:math:`\isotm{Mn}{51}`, and :math:`\isotm{Co}{55}`. - -Further simplifications include: - -* The reverse rates of :math:`\isotm{C}{12}+\isotm{C}{12}`, - :math:`\isotm{C}{12}+\isotm{O}{16}`, :math:`\isotm{O}{16}+\isotm{O}{16}` are - neglected since they're not present in the original aprox13 network - -* The :math:`\isotm{C}{12}+\isotm{Ne}{20}` rate is removed - -* The :math:`(\alpha, \gamma)` links between :math:`\isotm{Na}{23}`, - :math:`\isotm{Al}{27}` and between :math:`\isotm{Al}{27}` and - :math:`\isotm{P}{31}` are removed, since they're not in the - original aprox13 network. - -The network appears as: - -.. figure:: subch_simple.png - :align: center - -The nuclei in gray are those that have been approximated about, but the links -are effectively accounted for in the approximate rates. - -.. warning:: Due to inclusion of the rate sequence, - ${}^{14}\mathrm{N}(\alpha, \gamma){}^{18}\mathrm{F}(\alpha, - \mathrm{p}){}^{21}\mathrm{Ne}$, there is an artificial end-point at - ${}^{22}\mathrm{Na}$. - -``subch_base`` --------------- - -``subch_base`` is the simplest subch network. It is created to reconcile the -artificial end-point at :math:`\isotm{Na}{22}`. This is done by excluding -:math:`\isotm{N}{14}`, :math:`\isotm{F}{18}`, :math:`\isotm{Ne}{21}`, -and :math:`\isotm{Na}{22}`. These nuclei were added to include -:math:`\isotm{N}{14}(\alpha, \gamma)\isotm{F}{18}(\alpha, p)\isotm{Ne}{21}` -rate sequence, which allows an enhancement to the -:math:`\isotm{C}{12}(p, \gamma)\isotm{N}{13}(\alpha, p)\isotm{O}{16}` -rate due to the additional proton release. However, we find the effect is not -extremely significant. - -.. figure:: subch_base.png - :align: center - -disabling rates ---------------- - -For all subch networks, there are 2 runtime parameters that can be used -to disable rates: - -* ``network.disable_p_c12__n13`` : if set to ``1``, then the rate - :math:`\isotm{C}{12}(p,\gamma)\isotm{N}{13}` and its inverse are - disabled. - -* ``network.disable_he4_n13__p_o16`` : if set to ``1``, then the rate - :math:`\isotm{N}{13}(\alpha,p)\isotm{O}{16}` and its inverse are - disabled. - -Together, these parameters allow us to turn off the sequence -:math:`\isotm{C}{12}(p,\gamma)\isotm{N}{13}(\alpha, p)\isotm{O}{16}` that -acts as a bypass for :math:`\isotm{C}{12}(\alpha, \gamma)\isotm{O}{16}`. diff --git a/Docs/source/nse.rst b/Docs/source/nse.rst index 74d55f0ae8..38c5403ab8 100644 --- a/Docs/source/nse.rst +++ b/Docs/source/nse.rst @@ -262,6 +262,11 @@ reaction network. It solves for the chemical potentials of the proton and neutron and from there gets the abundances of each of the nuclei under the assumption of NSE, following the procedure outlined in :cite:`Calder_2007`. +.. important:: + + Self-consistent NSE does not support the templated C++ networks + (like ``aprox13``). You should use a pynucastro-generated network. + The solve is done using a port of the hybrid Powell method from MINPACK (we ported the solver to templated C++). diff --git a/Docs/source/ode_integrators.rst b/Docs/source/ode_integrators.rst index 94ed58f72e..6da6d8d6a3 100644 --- a/Docs/source/ode_integrators.rst +++ b/Docs/source/ode_integrators.rst @@ -15,9 +15,14 @@ provide a routine to convert from the integrator’s internal representation to the ``burn_t`` type required by the ``actual_rhs`` and ``actual_jac`` routine. -The name of the integrator can be selected at compile time using -the ``INTEGRATOR_DIR`` variable in the makefile. Presently, -the allowed options are: +.. index:: INTEGRATOR_DIR + +.. note:: + + The integrator is chosen at compile-time using + the ``INTEGRATOR_DIR`` variable in the makefile. + +Presently, allowed integrators are: * ``BackwardEuler``: an implicit first-order accurate backward-Euler method. An error estimate is done by taking 2 half steps and diff --git a/Docs/source/one_zone_tests.rst b/Docs/source/one_zone_tests.rst index 43a1e9ea48..8a34952ef4 100644 --- a/Docs/source/one_zone_tests.rst +++ b/Docs/source/one_zone_tests.rst @@ -3,180 +3,13 @@ One Zone Tests ************** There are several tests that let you call the EOS or reaction network -on a single zone to inspect the output directly. +on a single zone to inspect the output directly. Some of these +have analysis scripts, which we describe in the next sections. +.. toctree:: + :maxdepth: 1 + :hidden: -``burn_cell`` -============= - -.. index:: burn_cell - -``burn_cell`` is a simple one-zone burn that will evolve a state with -a network for a specified amount of time. This can be used to -understand the timescales involved in a reaction sequence or to -determine the needed ODE tolerances. This is designed to work -with the Strang-split integration wrappers. The system that is evolved -has the form: - -.. math:: - - \begin{align*} - \frac{dX_k}{dt} &= \dot{\omega}_k(\rho, X_k, T) \\ - \frac{de}{dt} &= \epsilon(\rho, X_k, T) - \end{align*} - -with density held constant and the temperature found via the equation of state, -$T = T(\rho, X_k, e)$. - - -.. note:: - - Since the energy evolves due to the heat release (or loss) - from reactions, the temperature will change during the burn - (unless ``integrator.call_eos_in_rhs=0`` is set). - - -Getting Started ---------------- - -The ``burn_cell`` code is located in -``Microphysics/unit_test/burn_cell``. An inputs file which sets the -default parameters for your choice of network is needed to run the -test. There are a number of inputs files in the unit test directory -already with a name list ``inputs_network``, where ``network`` -is the network you wish to use for your testing. These can be -used as a starting point for any explorations. - - -Setting the thermodynamics -^^^^^^^^^^^^^^^^^^^^^^^^^^ - -The parameters that affect the thermodynamics are: - -* ``unit_test.density`` : the initial density - -* ``unit_test.temperature`` : the initial temperature - -* ``unit_test.small_temp`` : the low temperature cutoff used in the equation of state - -* ``unit_test.small_dens`` : the low density cutoff used in the equation of state - -The composition can be set either by setting each mass fraction explicitly via the -parameters, ``unit_test.X1``, ``unit_test.X2``, ..., -e.g.: - -:: - - unit_test.X1 = 0.5 - unit_test.X2 = 0.2 - unit_test.X3 = 0.2 - unit_test.X4 = 0.1 - -where parameters up to ``X35`` are available. If the values don't sum to ``1`` -initially, then the test will do a normalization. This normalization can be -disabled by setting: - -:: - - unit_test.skip_initial_normalization = 1 - - -Alternately, the composition can be set automatically by initializing all -of the mass fractions equally (to $1/N$, where $N$ is the number of species), -by setting: - -:: - - unit_test.init_species_all_equal = 1 - - -Controlling time -^^^^^^^^^^^^^^^^ - -The test will run unit a time ``unit_test.tmax``, outputting the state -at regular intervals. The parameters controlling the output are: - -* ``unit_test.tmax`` : the end point of integration. - -* ``unit_test.tfirst`` : the first time interval to output. - -* ``unit_test.nsteps`` : the number of steps to divide the integration into, - logarithmically-spaced. - -If there is only a single step, ``unit_test.nsteps = 1``, then we integrate -from $[0, \mathrm{tmax}]$. - -If there are multiple steps, then the first output will be at a time -$\mathrm{tmax} / \mathrm{nsteps}$, and the steps will be -logarithmically-spaced afterwards. - - -Integration parameters -^^^^^^^^^^^^^^^^^^^^^^ - -The tolerances, choice of Jacobian, and other integration parameters -can be set via the usual Microphysics runtime parameters, e.g. -``integrator.atol_spec``. - - -Building and Running the Code ------------------------------ - -The code can be built simply as: - -.. prompt:: bash - - make - -and the network and integrator can be changed using the normal -Microphysics build system parameters, e.g., - -.. prompt:: bash - - make NETWORK_DIR=aprox19 INTEGRATOR_DIR=rkc - -The build process will automatically create links in the build -directory to the EOS table and any reaction rate tables needed by your -choice of network. - - -.. important:: - - You need to do a ``make clean`` before rebuilding with a different - network or integrator. - - -To run the code, enter the burn_cell directory and run:: - - ./main3d.gnu.ex inputs - -where ``inputs`` is the name of your inputs file. - -Working with Output -------------------- - -.. note:: - - For this part, we'll assume that the default ``aprox13`` and - ``VODE`` options were used for the network and integrator, and the - test was run with ``inputs.aprox13``. - -As the code runs, it will output to ``stdout`` details of the initial -and final state and the number of integration steps taken (along with whether -the burn was successful). The full history of the thermodynamic state will also be output to a file, -``state_over_time.txt``, with each line corresponding to one of the -``nsteps`` requested in the time integration. - -The script ``plot_burn_cell.py`` can be used to visualize the evolution: - -.. prompt:: bash - - python plot_burn_cell.py state_over_time.txt - -This will generate the following figure: - -.. figure:: state.png - :alt: An example of a plot output by the burn_cell unit test. - -Only the most abundant species are plotted. The number of species to plot and the -limits of $X$ can be set via runtime parameters (see ``python plot_burn_cell.py -h``). + burn_cell.rst + eos_cell.rst + jac_cell.rst diff --git a/Docs/source/refs.bib b/Docs/source/refs.bib index 0137c25cd4..c8c3159a66 100644 --- a/Docs/source/refs.bib +++ b/Docs/source/refs.bib @@ -721,6 +721,20 @@ @ARTICLE{itoh:1996 adsnote = {Provided by the SAO/NASA Astrophysics Data System} } +@ARTICLE{iso7, + author = {{Timmes}, F.~X. and {Hoffman}, R.~D. and {Woosley}, S.~E.}, + title = "{An Inexpensive Nuclear Energy Generation Network for Stellar Hydrodynamics}", + journal = {\apjs}, + keywords = {Hydrodynamics, Methods: Numerical, Nuclear Reactions, Nucleosynthesis, Abundances, Stars: General}, + year = 2000, + month = jul, + volume = {129}, + number = {1}, + pages = {377-398}, + doi = {10.1086/313407}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2000ApJS..129..377T}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} @article{amrex_joss, doi = {10.21105/joss.01370}, @@ -772,3 +786,19 @@ @ARTICLE{quokka adsnote = {Provided by the SAO/NASA Astrophysics Data System} } +@ARTICLE{tillotson:1962, + author = {{Tillotson}, J. H.}, + title = "{Metallic Equations of State for Hypervelocity Impact}", + year = {1962}, + journal = {Internal Report, General Atomics}, + url = {https://apps.dtic.mil/sti/citations/AD0486711} +} + +@article{lsode, + title = {Description and {Use} of {LSODE}, the {Livermore} {Solver} for {Ordinary} {Differential} {Equations}}, + language = {en}, + year = {1993}, + author = {{Radhakrishnan}, Krishnan and {Hindmarsh}, Alan C.}, + journal = {Lawrence Livermore National Laboratory Report UCRL-ID-113855}, + pages = {124} +} \ No newline at end of file diff --git a/Docs/source/rp_intro.rst b/Docs/source/rp_intro.rst index 027de8f709..70f9034f77 100644 --- a/Docs/source/rp_intro.rst +++ b/Docs/source/rp_intro.rst @@ -6,7 +6,7 @@ The behavior of the network and EOS are controlled by many runtime parameters. These parameters are defined in plain-text files ``_parameters`` located in the different directories that hold the microphysics code. At compile time, a script in the AMReX build -system, findparams.py, locates all of the ``_parameters`` files that +system, ``findparams.py``, locates all of the ``_parameters`` files that are needed for the given choice of network, integrator, and EOS, and assembles all of the runtime parameters into a set of header files (using the ``write_probin.py`` script). @@ -32,4 +32,3 @@ noted in separate tables. .. toctree:: runtime_parameters - diff --git a/Docs/source/screening.rst b/Docs/source/screening.rst index dba3ec7b6a..7d9cbd25ce 100644 --- a/Docs/source/screening.rst +++ b/Docs/source/screening.rst @@ -94,6 +94,22 @@ The options are: weak screening regime, :math:`\Gamma < 0.1`, and strong screening regime, :math:`1 \lesssim \Gamma \lesssim 160`. +.. index:: screening.enable_debye_huckel_skip, screening.debye_huckel_skip_threshold + +* ``debye_huckel`` : + + This is just the Debye-Hückel weak-screening limit from + :cite:`chugunov:2009`. + + While it can be used on its own (by building with + ``SCREEN_METHOD=debye_huckel``, it is really meant to be used as a + test to determine whether a more extensive screening approximation + should be used. By setting ``screening.enable_debye_huckel_skip``, + we first compute this weak-screening approximation and then, if it + is larger than ``screening.debye_huckel_skip_threshold``, the full + screening factor is computed (using the method specified via + ``SCREEN_METHOD``). + * ``null`` : This disables screening by always returning 1 for the screening diff --git a/Docs/source/sdc.rst b/Docs/source/sdc.rst index 1d9b52e29d..2eed54cceb 100644 --- a/Docs/source/sdc.rst +++ b/Docs/source/sdc.rst @@ -1,3 +1,5 @@ +.. _sdc-evolution: + ***************************** Spectral Deferred Corrections ***************************** diff --git a/Docs/source/subch_base.png b/Docs/source/subch_base.png deleted file mode 100644 index ea6ab93d7d..0000000000 Binary files a/Docs/source/subch_base.png and /dev/null differ diff --git a/Docs/source/subch_simple.png b/Docs/source/subch_simple.png deleted file mode 100644 index 8ecc229a30..0000000000 Binary files a/Docs/source/subch_simple.png and /dev/null differ diff --git a/Docs/source/templated_networks.rst b/Docs/source/templated_networks.rst index c6bb55cdae..50f419e545 100644 --- a/Docs/source/templated_networks.rst +++ b/Docs/source/templated_networks.rst @@ -1,3 +1,5 @@ +.. _sec:templated_rhs: + ********************************* Templated Network Righthand Sides ********************************* @@ -384,16 +386,11 @@ Linear Algebra The VODE integrator needs routines to do LU factorization and back substitution. We build off of the linpack ``dgefa`` and ``dgesl`` -routines, but because we know at compile time which Jacobian terms are -non-zero, we are able to use ``constexpr`` for-loops to only do the -calculations on non-zero elements. This greatly reduces the amount of work -in the linear algebra. - -Note: +routines, using ``constexpr`` for-loops that know the network size +at compile time. -* Currently we are still storing a dense Jacobian -- we just skip computation - on the elements that are 0. +.. note:: -* These routines do not perform pivoting. This does not seem to be an - issue for the types of matrices we solve with reactions (since they are - all of the form :math:`I - \tau J`, where :math:`tau` is the timestep). + These routines do not perform pivoting. This does not seem to be an + issue for the types of matrices we solve with reactions (since they are + all of the form :math:`I - \tau J`, where :math:`tau` is the timestep). diff --git a/Docs/source/unit_tests.rst b/Docs/source/unit_tests.rst index 5a36b2c4af..6c01f664c7 100644 --- a/Docs/source/unit_tests.rst +++ b/Docs/source/unit_tests.rst @@ -115,12 +115,12 @@ by options in the input file. One-zone tests ============== -.. index:: burn_cell, burn_cell_primordial_chem, burn_cell_sdc, eos_cell, jac_cell, nse_table_cell, test_ase, test_part_func +.. index:: burn_cell, burn_cell_primordial_chem, burn_cell_sdc, eos_cell, jac_cell, nse_table_cell, nse_net_cell, part_func_cell * ``burn_cell`` : given a $\rho$, $T$, and $X_k$, integrate a reaction network through a specified time - and output the new state. + and output the new state. See :ref:`sec:burn_cell` for more information. * ``burn_cell_primordial_chem`` : @@ -133,27 +133,27 @@ One-zone tests * ``eos_cell`` : given a $\rho$, $T$, and $X_k$, call the equation of state and print out - the thermodynamic information. + the thermodynamic information. See :ref:`sec:eos_cell` for more information. * ``jac_cell`` : for a single thermodynamic state, compute the analytic Jacobian (using the functions provided by the network) and a numerical finite-difference approximation to the Jacobian and output them, - element-by-element, to the display. + element-by-element, to the display. See :ref:`sec:jac_cell` for more information. * ``nse_table_cell`` : given a $\rho$, $T$, and $Y_e$, evaluate the NSE state via table interpolation and print it out. -* ``test_ase`` : +* ``nse_net_cell`` : for the self-consistent NSE, take a $\rho$, $T$, and $Y_e$, and solve for the NSE state. Then check the NSE condition to see if we are actually satisfying the NSE criteria for the network. -* ``test_part_func`` +* ``part_func_cell`` exercise the partition function interpolation for a few select nuclei. diff --git a/EOS/README.md b/EOS/README.md index 6f8af294d4..fb8dadeeeb 100644 --- a/EOS/README.md +++ b/EOS/README.md @@ -1,20 +1,31 @@ The following equations of state are provided: -* `gamma_law`: a general constant-gamma ideal gas. You can - set whether or not the gas is ionized via a runtime parameter. +* `breakout`: this is essentially the same as `gamma_law`, but it sets + the mean molecular weight based on ``eos_t aux[]`` instead of the + mass fractions. -* `helmholtz`: a general stellar equation of state consisting of - ions, degenerate/relativitic electrons, and radiation. This is Frank +* `gamma_law` : a general constant-gamma ideal gas. You can set + whether or not the gas is ionized via a runtime parameter. + +* `helmholtz`: a general stellar equation of state consisting of ions, + degenerate/relativitic electrons, and radiation. This is Frank Timmes's EOS with a modified interface and thread-safety. +* `metal_chem`: a multi-gamma EOS for metal ISM chemistry. + * `multigamma`: a gamma-law EOS where each species can have its own gamma (ratio of specific heats). * `polytrope`: a simple polytropic EOS: `P = K rho**gamma` -* `stellarcollapse`: the nuclear equation of state from - stellarcollapse.org +* `primordial_chem`: a version of multigamma for primordial ISM + chemistry + +* `rad_power_law`: an artificial EOS for radiation transport test + problems. + +* `tillotson`: an EOS for hypevelocity impacts based on Tillotson + (1962) * `ztwd`: a zero-temperature white dwarf equation of state -* `primordial_chem`: a version of multigamma for primordial ISM chemistry diff --git a/integration/utils/numerical_jacobian.H b/integration/utils/numerical_jacobian.H index 9822380410..d396c69bcd 100644 --- a/integration/utils/numerical_jacobian.H +++ b/integration/utils/numerical_jacobian.H @@ -104,7 +104,7 @@ void numerical_jac(BurnT& state, const jac_info_t& jac_info, JacNetArray2D& jac) w = integrator_rp::rtol_spec * std::abs(yj) + integrator_rp::atol_spec; - // the incremement we use in the derivative is defined in the LSODE paper, Eq. 3.35 + // the increment we use in the derivative is defined in the LSODE paper, Eq. 3.35 amrex::Real dy = amrex::max(std::sqrt(U) * std::abs(yj), r0 * w); diff --git a/interfaces/eos.H b/interfaces/eos.H index c5ecfd7cf4..3abccc2c77 100644 --- a/interfaces/eos.H +++ b/interfaces/eos.H @@ -394,14 +394,13 @@ void check_inputs (const I input, T& state) template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void eos (const I input, T& state, bool use_raw_inputs = false) +void eos (const I input, T& state) { static_assert(std::is_same_v, "input must be an eos_input_t"); // Input arguments bool has_been_reset = false; - bool use_composition_routine = true; // Local variables @@ -411,15 +410,9 @@ void eos (const I input, T& state, bool use_raw_inputs = false) } #endif - if (use_raw_inputs) { - use_composition_routine = false; - } - if constexpr (has_xn::value) { - if (use_composition_routine) { - // Get abar, zbar, etc. - composition(state); - } + // Get abar, zbar, etc. + composition(state); } // Force the inputs to be valid. diff --git a/networks/he-burn/cno-he-burn-33a/cno-he-burn-33a.png b/networks/he-burn/cno-he-burn-33a/cno-he-burn-33a.png new file mode 100644 index 0000000000..e0d2e6c75b Binary files /dev/null and b/networks/he-burn/cno-he-burn-33a/cno-he-burn-33a.png differ diff --git a/networks/he-burn/cno-he-burn-33a/cno_he_burn_33a.py b/networks/he-burn/cno-he-burn-33a/cno_he_burn_33a.py index 3e09f97928..bd2c6150be 100644 --- a/networks/he-burn/cno-he-burn-33a/cno_he_burn_33a.py +++ b/networks/he-burn/cno-he-burn-33a/cno_he_burn_33a.py @@ -39,14 +39,15 @@ def doit(): rho = 1.e6 T = 1.e9 - net.plot(rho, T, comp, outfile="cno-he-burn-33a.png", - rotated=True, hide_xalpha=True, curved_edges=True, - size=(1500, 450), - node_size=500, node_font_size=11, node_color="#337dff", node_shape="s", - Z_range=(1, 29)) + fig = net.plot(rho, T, comp, + rotated=True, hide_xalpha=True, curved_edges=True, + size=(1500, 450), + node_size=500, node_font_size=11, node_color="#337dff", node_shape="s", + Z_range=(1, 29)) net.write_network() + fig.savefig("cno-he-burn-33a.png", bbox_inches="tight") if __name__ == "__main__": doit() diff --git a/networks/he-burn/he-burn-31anp/he-burn-31anp.png b/networks/he-burn/he-burn-31anp/he-burn-31anp.png index c869945c08..25d6c93aee 100644 Binary files a/networks/he-burn/he-burn-31anp/he-burn-31anp.png and b/networks/he-burn/he-burn-31anp/he-burn-31anp.png differ diff --git a/networks/he-burn/he-burn-31anp/he_burn_31anp.py b/networks/he-burn/he-burn-31anp/he_burn_31anp.py index 6774ab48a3..2c80930974 100644 --- a/networks/he-burn/he-burn-31anp/he_burn_31anp.py +++ b/networks/he-burn/he-burn-31anp/he_burn_31anp.py @@ -46,7 +46,7 @@ def doit(): size=(1800, 900), node_size=500, node_shape="s", node_color="#337dff", node_font_size=10) - fig.savefig("he-burn-31anp.png") + fig.savefig("he-burn-31anp.png", bbox_inches="tight") net.write_network() diff --git a/networks/he-burn/he-burn-36a/he-burn-36a.png b/networks/he-burn/he-burn-36a/he-burn-36a.png index 47578149b8..f08c0bb8c6 100644 Binary files a/networks/he-burn/he-burn-36a/he-burn-36a.png and b/networks/he-burn/he-burn-36a/he-burn-36a.png differ diff --git a/networks/he-burn/he-burn-36a/he_burn_36a.py b/networks/he-burn/he-burn-36a/he_burn_36a.py index 04493ecea4..8b22c854e6 100644 --- a/networks/he-burn/he-burn-36a/he_burn_36a.py +++ b/networks/he-burn/he-burn-36a/he_burn_36a.py @@ -40,7 +40,7 @@ def doit(): size=(1800, 900), node_size=500, node_shape="s", node_color="#337dff", node_font_size=10) - fig.savefig("he-burn-36a.png") + fig.savefig("he-burn-36a.png", bbox_inches="tight") net.write_network() diff --git a/requirements.txt b/requirements.txt index ac3ebbf470..376a90128a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,3 +11,4 @@ numpy>=1.13.3 sphinx-prompt sphinx-math-dollar sphinx-copybutton +sphinx-mdinclude diff --git a/screening/_parameters b/screening/_parameters index ca5a6b79d1..672c0a9a41 100644 --- a/screening/_parameters +++ b/screening/_parameters @@ -1,3 +1,5 @@ @namespace: screening enable_chabrier1998_quantum_corr bool 0 +enable_debye_huckel_skip bool 0 +debye_huckel_skip_threshold real 1.01e0 diff --git a/screening/screen.H b/screening/screen.H index ffeacc374f..54999aaddf 100644 --- a/screening/screen.H +++ b/screening/screen.H @@ -7,6 +7,7 @@ #define SCREEN_METHOD_chugunov2007 2 #define SCREEN_METHOD_chugunov2009 3 #define SCREEN_METHOD_chabrier1998 4 +#define SCREEN_METHOD_debye_huckel 5 #include #include @@ -152,6 +153,33 @@ fill_plasma_state(plasma_state_t& state, const number_t& temp, state.gamma_e_fac = gamma_e_constants * std::cbrt(state.n_e); } +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +number_t debye_huckel (const plasma_state_t& state, + const scrn::screen_factors_t& scn_fac) +{ + // Calculates the Debye-Huckel enhancement factor for weak Coloumb coupling, + // as listed in the Appendix of Chugunov and DeWitt 2009, PhRvC, 80, 014611 + + // input: + // state = plasma state (T, rho, abar, zbar, etc.) + // scn_fac = screening factors for A and Z + + amrex::Real z1z2 = scn_fac.z1 * scn_fac.z2; + amrex::Real z_ratio = state.z2bar / state.zbar; + + // Gamma_e from eq. 6 + number_t Gamma_e = state.gamma_e_fac / state.temp; + + // eq. A1 + number_t h_DH = z1z2 * admath::sqrt(3 * amrex::Math::powi<3>(Gamma_e) * z_ratio); + + // machine limit the output + constexpr amrex::Real h_max = 300.e0_rt; + number_t h = admath::min(h_DH, h_max); + return admath::exp(h); +} + #if SCREEN_METHOD == SCREEN_METHOD_screen5 template AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -620,11 +648,20 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE number_t actual_screen(const plasma_state_t& state, const scrn::screen_factors_t& scn_fac) { - number_t scor; + number_t scor = 1.0_rt; #if SCREEN_METHOD == SCREEN_METHOD_null // null screening amrex::ignore_unused(state, scn_fac); - scor = 1.0_rt; + return scor; +#endif + if (screening_rp::enable_debye_huckel_skip) { + scor = debye_huckel(state, scn_fac); + if (scor <= screening_rp::debye_huckel_skip_threshold) { + return scor; + } + } +#if SCREEN_METHOD == SCREEN_METHOD_debye_huckel + scor = debye_huckel(state, scn_fac); #elif SCREEN_METHOD == SCREEN_METHOD_screen5 scor = actual_screen5(state, scn_fac); #elif SCREEN_METHOD == SCREEN_METHOD_chugunov2007 diff --git a/unit_test/_parameters b/unit_test/_parameters index 755f5354cb..2fd9224d2e 100644 --- a/unit_test/_parameters +++ b/unit_test/_parameters @@ -39,3 +39,6 @@ X32 real 0.0e0 X33 real 0.0e0 X34 real 0.0e0 X35 real 0.0e0 + +uniform_xn bool 0 + diff --git a/unit_test/burn_cell/_parameters b/unit_test/burn_cell/_parameters index 4ad72846d0..5a39427b8c 100644 --- a/unit_test/burn_cell/_parameters +++ b/unit_test/burn_cell/_parameters @@ -17,5 +17,3 @@ density real 1.e7 temperature real 3.e9 skip_initial_normalization bool 0 - -init_species_all_equal bool 0 diff --git a/unit_test/burn_cell/burn_cell.H b/unit_test/burn_cell/burn_cell.H index 43496f1a67..482fb580f9 100644 --- a/unit_test/burn_cell/burn_cell.H +++ b/unit_test/burn_cell/burn_cell.H @@ -21,14 +21,10 @@ void burn_cell_c() // Make sure user set all the mass fractions to values in the interval [0, 1] for (int n = 1; n <= NumSpec; ++n) { - if (unit_test_rp::init_species_all_equal) { - massfractions[n-1] = 1.0_rt / static_cast(NumSpec); - } else { - massfractions[n-1] = get_xn(n); + massfractions[n-1] = get_xn(n, unit_test_rp::uniform_xn); - if (massfractions[n-1] < 0 || massfractions[n-1] > 1) { - amrex::Error("mass fraction for " + short_spec_names_cxx[n-1] + " not initialized in the interval [0,1]!"); - } + if (massfractions[n-1] < 0 || massfractions[n-1] > 1) { + amrex::Error("mass fraction for " + short_spec_names_cxx[n-1] + " not initialized in the interval [0,1]!"); } } diff --git a/unit_test/burn_cell/ci-benchmarks/triple_alpha_plus_cago_FE_unit_test.out b/unit_test/burn_cell/ci-benchmarks/triple_alpha_plus_cago_FE_unit_test.out index d07b621b44..d985a36dc4 100644 --- a/unit_test/burn_cell/ci-benchmarks/triple_alpha_plus_cago_FE_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/triple_alpha_plus_cago_FE_unit_test.out @@ -33,7 +33,4 @@ omegadot(C12): 0.01497617246 omegadot(O16): 0.5494210118 omegadot(Fe56): 0 number of steps taken: 85669 -Unused ParmParse Variables: - [TOP]::unit_test.x4(nvals = 1) :: [0.0] - AMReX (23.07-7-g88f03408f18a) finalized diff --git a/unit_test/burn_cell/inputs_triple b/unit_test/burn_cell/inputs_triple index fef70305ca..7480f09a1c 100644 --- a/unit_test/burn_cell/inputs_triple +++ b/unit_test/burn_cell/inputs_triple @@ -23,4 +23,4 @@ unit_test.temperature = 3.e9 unit_test.X1 = 1.0 unit_test.X2 = 0.0 unit_test.X3 = 0.0 -unit_test.x4 = 0.0 +unit_test.X4 = 0.0 diff --git a/unit_test/burn_cell_sdc/burn_cell.H b/unit_test/burn_cell_sdc/burn_cell.H index 1583042eaa..f5e926e1c2 100644 --- a/unit_test/burn_cell_sdc/burn_cell.H +++ b/unit_test/burn_cell_sdc/burn_cell.H @@ -20,7 +20,7 @@ void burn_cell_c() // Make sure user set all the mass fractions to values in the interval [0, 1] for (int n = 1; n <= NumSpec; ++n) { - massfractions[n-1] = get_xn(n); + massfractions[n-1] = get_xn(n, unit_test_rp::uniform_xn); if (massfractions[n-1] < 0 || massfractions[n-1] > 1) { amrex::Error("mass fraction for " + short_spec_names_cxx[n-1] + " not initialized in the interval [0,1]!"); diff --git a/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out b/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out index bf346e3bde..bfc081b6e5 100644 --- a/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out +++ b/unit_test/eos_cell/ci-benchmarks/eos_gamma_law.out @@ -31,7 +31,4 @@ gam1 = 1.666666667 cs = 186127884.1 abar = 4 zbar = 2 -Unused ParmParse Variables: - [TOP]::unit_test.x4(nvals = 1) :: [0.0] - AMReX (23.06-31-g9c256b12b89f) finalized diff --git a/unit_test/eos_cell/ci-benchmarks/eos_helmholtz.out b/unit_test/eos_cell/ci-benchmarks/eos_helmholtz.out index 0277ec3084..69f2d64bb0 100644 --- a/unit_test/eos_cell/ci-benchmarks/eos_helmholtz.out +++ b/unit_test/eos_cell/ci-benchmarks/eos_helmholtz.out @@ -31,7 +31,4 @@ gam1 = 1.511401882 cs = 332889390.6 abar = 4 zbar = 2 -Unused ParmParse Variables: - [TOP]::unit_test.x4(nvals = 1) :: [0.0] - AMReX (22.10-75-g0a3ee1486cc3) finalized diff --git a/unit_test/eos_cell/eos_cell.H b/unit_test/eos_cell/eos_cell.H index e1d7a5c438..da033d7c20 100644 --- a/unit_test/eos_cell/eos_cell.H +++ b/unit_test/eos_cell/eos_cell.H @@ -24,7 +24,7 @@ void eos_cell_c() // Make sure user set all the mass fractions to values in the interval [0, 1] for (int n = 1; n <= NumSpec; ++n) { - massfractions[n-1] = get_xn(n); + massfractions[n-1] = get_xn(n, unit_test_rp::uniform_xn); if (massfractions[n-1] < 0 || massfractions[n-1] > 1) { amrex::Error("mass fraction for " + short_spec_names_cxx[n-1] + " not initialized in the interval [0,1]!"); diff --git a/unit_test/eos_cell/inputs_eos b/unit_test/eos_cell/inputs_eos index 13f9218081..118d6c74f8 100644 --- a/unit_test/eos_cell/inputs_eos +++ b/unit_test/eos_cell/inputs_eos @@ -7,6 +7,6 @@ unit_test.temperature = 1.e9 unit_test.X1 = 1.0 unit_test.X2 = 0.0 unit_test.X3 = 0.0 -unit_test.x4 = 0.0 +unit_test.X4 = 0.0 diff --git a/unit_test/jac_cell/README.md b/unit_test/jac_cell/README.md index 0d972f1f25..c09b05e00e 100644 --- a/unit_test/jac_cell/README.md +++ b/unit_test/jac_cell/README.md @@ -3,3 +3,9 @@ For a single thermodynamic condition, evaluate the Jacobian using both the analytic version from the network and the finite-difference numerical approximation. + +Note: most analytic Jacobians do not consider the derivative of the +screening factor with respect to composition. To get a more accurate +comparison of how the analytic and numerical terms compare, considering +just the reactive part, build with `SCREEN_METHOD=null`. + diff --git a/unit_test/jac_cell/_parameters b/unit_test/jac_cell/_parameters index 9047a84f8e..5894a90b97 100644 --- a/unit_test/jac_cell/_parameters +++ b/unit_test/jac_cell/_parameters @@ -1,19 +1,8 @@ @namespace: unit_test -run_prefix string "" - small_temp real 1.e5 small_dens real 1.e5 -# the final time to integrate to -tmax real 1.e-2 - -# first output time -- we will output in nsteps logarithmically spaced steps between [tfirst, tmax] -tfirst real 0.0 - -# number of steps (logarithmically spaced) -nsteps int 100 - density real 1.e7 temperature real 3.e9 diff --git a/unit_test/jac_cell/ci-benchmarks/jac_cell_aprox13.out b/unit_test/jac_cell/ci-benchmarks/jac_cell_aprox13.out index b85220eaee..8a3899192b 100644 --- a/unit_test/jac_cell/ci-benchmarks/jac_cell_aprox13.out +++ b/unit_test/jac_cell/ci-benchmarks/jac_cell_aprox13.out @@ -1,214 +1,215 @@ -Initializing AMReX (24.07-16-gdcb9cc0383dc)... -AMReX (24.07-16-gdcb9cc0383dc) initialized +Initializing AMReX (24.12-14-gb3f67385e62f)... +AMReX (24.12-14-gb3f67385e62f) initialized starting the single zone burn... -J( 1, 1) = -22.83252156 -22.76871089 -J( 1, 2) = -7074.184809 -7070.885902 -J( 1, 3) = -146332.0024 -146331.961 -J( 1, 4) = -2840218.246 -2840216.758 -J( 1, 5) = -1480087.416 -1480081.685 -J( 1, 6) = -150293.5967 -150288.8284 -J( 1, 7) = -155563.9558 -155558.7711 -J( 1, 8) = -90256.26203 -90254.79555 -J( 1, 9) = -4338.739586 -4338.971348 -J( 1, 10) = -1367.526773 -1364.641738 -J( 1, 11) = -754.4182788 -751.4739669 -J( 1, 12) = -164.8910475 -163.1314012 -J( 1, 13) = 0.177283612 8.045249939e-06 -J( 1, 14) = 2.036804679e-18 2.036804635e-18 + element numerical analytic +J( He4, He4) = -386901.0562 -389364.2481 +J( He4, C12) = 1538.83613 1527.484671 +J( He4, O16) = -11984.43367 -11604.67208 +J( He4, Ne20) = -226229.9236 -225570.3475 +J( He4, Mg24) = -118533.7928 -117650.0754 +J( He4, Si28) = -13102.92749 -12026.86162 +J( He4, S32) = -13739.13964 -12490.59603 +J( He4, Ar36) = -8735.824012 -7328.055338 +J( He4, Ca40) = -1917.505444 -359.8025782 +J( He4, Ti44) = -1801.580358 -100.6766494 +J( He4, Cr48) = -1902.430547 -63.38388515 +J( He4, Fe52) = -1987.160097 -13.86553734 +J( He4, Ni56) = -2104.486528 8.889625789e-06 +J( He4, e) = -3.302305762e-14 -3.302305672e-14 -J( 2, 1) = 22.83252156 22.76871089 -J( 2, 2) = -21223.06355 -21215.80687 -J( 2, 3) = -0.1431906097 0.01577663496 -J( 2, 4) = -0.1527366503 0 -J( 2, 5) = -0.1591006774 0 -J( 2, 6) = -0.1636464111 0 -J( 2, 7) = -0.1670557113 0 -J( 2, 8) = -0.1697073892 0 -J( 2, 9) = -0.1718287316 0 -J( 2, 10) = -0.1735643754 0 -J( 2, 11) = -0.1750107452 0 -J( 2, 12) = -0.1762345965 0 -J( 2, 13) = -0.177283612 0 -J( 2, 14) = -2.036804679e-18 -2.036804635e-18 +J( C12, He4) = -1544.888563 -1669.808542 +J( C12, C12) = -14163.60628 -14170.71908 +J( C12, O16) = -26.31037223 -15.90510669 +J( C12, Ne20) = -22.40516055 0 +J( C12, Mg24) = -31.64594981 0 +J( C12, Si28) = -39.31004354 0 +J( C12, S32) = -45.98896291 0 +J( C12, Ar36) = -52.01096713 0 +J( C12, Ca40) = -57.57313075 0 +J( C12, Ti44) = -62.80086484 0 +J( C12, Cr48) = -67.77757842 0 +J( C12, Fe52) = -72.56164943 0 +J( C12, Ni56) = -77.19392628 0 +J( C12, e) = -1.637120092e-15 -1.637119821e-15 -J( 3, 1) = -5.632962675e-25 -5.570411722e-25 -J( 3, 2) = 28286.69277 28286.69277 -J( 3, 3) = -585327.8861 -585327.8861 -J( 3, 4) = 6128.217777 6128.217777 -J( 3, 5) = -1.739488553e-26 0 -J( 3, 6) = -1.789188226e-26 0 -J( 3, 7) = -1.826462981e-26 0 -J( 3, 8) = -1.855454457e-26 0 -J( 3, 9) = -1.878647638e-26 0 -J( 3, 10) = -1.897623876e-26 0 -J( 3, 11) = -1.913437409e-26 0 -J( 3, 12) = -1.92681809e-26 0 -J( 3, 13) = -1.938287245e-26 0 -J( 3, 14) = -2.226890848e-43 -2.226890927e-43 +J( O16, He4) = -44251.69974 -44202.73796 +J( O16, C12) = 2167.279422 2205.366875 +J( O16, O16) = -46494.3376 -46450.60357 +J( O16, Ne20) = 6268.468387 6319.326423 +J( O16, Mg24) = -58.72063931 0 +J( O16, Si28) = -67.00603097 0 +J( O16, S32) = -75.55571467 0 +J( O16, Ar36) = -84.28093184 0 +J( O16, Ca40) = -93.12930016 0 +J( O16, Ti44) = -102.0668362 0 +J( O16, Cr48) = -111.0723391 0 +J( O16, Fe52) = -120.129117 0 +J( O16, Ni56) = -129.2264943 0 +J( O16, e) = 4.384415228e-16 4.384411885e-16 -J( 4, 1) = -1.366807088e-23 -1.347708423e-23 -J( 4, 2) = 4.753870727e-12 1.2621475e-25 -J( 4, 3) = 731659.8313 731659.8313 -J( 4, 4) = -14216404.34 -14216404.34 -J( 4, 5) = 0.003363940593 0.003363940593 -J( 4, 6) = -1.093051048e-24 0 -J( 4, 7) = -1.115822945e-24 0 -J( 4, 8) = -1.13353442e-24 0 -J( 4, 9) = -1.1477036e-24 0 -J( 4, 10) = -1.159296566e-24 0 -J( 4, 11) = -1.168957371e-24 0 -J( 4, 12) = -1.177131898e-24 0 -J( 4, 13) = -1.184138635e-24 0 -J( 4, 14) = -1.360452377e-41 -1.360452352e-41 +J(Ne20, He4) = -1071437.231 -1077714.229 +J(Ne20, C12) = 10295.4965 10403.37149 +J(Ne20, O16) = 56883.61948 58036.66646 +J(Ne20, Ne20) = -1145561.976 -1143650.054 +J(Ne20, Mg24) = -2527.616919 0.003495578555 +J(Ne20, Si28) = -3061.548719 0 +J(Ne20, S32) = -3544.317148 0 +J(Ne20, Ar36) = -3993.040141 0 +J(Ne20, Ca40) = -4417.866583 0 +J(Ne20, Ti44) = -4825.366803 0 +J(Ne20, Cr48) = -5219.853315 0 +J(Ne20, Fe52) = -5604.291196 0 +J(Ne20, Ni56) = -5980.914429 0 +J(Ne20, e) = -8.494068043e-14 -8.494067874e-14 -J( 5, 1) = 8.302657095e-24 8.170002764e-24 -J( 5, 2) = 3.153014201e-25 1.88760336e-28 -J( 5, 3) = 3.547140976e-25 1.88760336e-28 -J( 5, 4) = 17050492.88 17050492.88 -J( 5, 5) = -8880490.115 -8880490.115 -J( 5, 6) = 0.0001669112793 0.0001669112793 -J( 5, 7) = -5.810350496e-23 0 -J( 5, 8) = -5.809693618e-23 0 -J( 5, 9) = -5.809168116e-23 0 -J( 5, 10) = -5.808738159e-23 0 -J( 5, 11) = -5.808379862e-23 0 -J( 5, 12) = -5.808076688e-23 0 -J( 5, 13) = -5.807816824e-23 0 -J( 5, 14) = 5.045605543e-42 5.045605749e-42 +J(Mg24, He4) = 654844.2806 657000.6182 +J(Mg24, C12) = 435.4230384 15.9212509 +J(Mg24, O16) = 929.9962063 15.9212509 +J(Mg24, Ne20) = 1364203.922 1362901.075 +J(Mg24, Mg24) = -704261.6416 -705900.4607 +J(Mg24, Si28) = 1944.546334 0.0001724818945 +J(Mg24, S32) = 2231.306642 0 +J(Mg24, Ar36) = 2505.576486 0 +J(Mg24, Ca40) = 2770.981131 0 +J(Mg24, Ti44) = 3030.000681 0 +J(Mg24, Cr48) = 3284.193323 0 +J(Mg24, Fe52) = 3534.649522 0 +J(Mg24, Ni56) = 3782.237956 0 +J(Mg24, e) = 3.138100654e-14 3.138100793e-14 -J( 6, 1) = 9.417373366e-24 9.308549997e-24 -J( 6, 2) = 7.540626412e-25 2.20220392e-28 -J( 6, 3) = 6.930014139e-18 2.204043829e-28 -J( 6, 4) = 5.942221324e-23 0 -J( 6, 5) = 10360571.8 10360571.8 -J( 6, 6) = -1052021.799 -1052021.799 -J( 6, 7) = 2.618223472 2.618223472 -J( 6, 8) = 5.952275493e-23 0 -J( 6, 9) = 5.953532264e-23 0 -J( 6, 10) = 5.954560531e-23 0 -J( 6, 11) = 5.955417421e-23 0 -J( 6, 12) = 5.956142481e-23 0 -J( 6, 13) = 5.956763961e-23 0 -J( 6, 14) = 1.206687442e-41 1.206687398e-41 +J(Si28, He4) = 733740.6263 739362.5011 +J(Si28, C12) = -236.8292666 18.57479272 +J(Si28, O16) = 574.578026 18.59079981 +J(Si28, Ne20) = 1117.750234 0 +J(Si28, Mg24) = 825105.1884 823550.5326 +J(Si28, Si28) = -82267.7835 -84188.03177 +J(Si28, S32) = 2244.127686 2.723790242 +J(Si28, Ar36) = 2532.726899 0 +J(Si28, Ca40) = 2803.275313 0 +J(Si28, Ti44) = 3058.658852 0 +J(Si28, Cr48) = 3302.735385 0 +J(Si28, Fe52) = 3538.1484 0 +J(Si28, Ni56) = 3766.666407 0 +J(Si28, e) = 7.404959203e-14 7.404958884e-14 -J( 7, 1) = -4.433213398e-26 -4.216253314e-26 -J( 7, 2) = -1.713032535e-26 0 -J( 7, 3) = 9.755839735e-19 2.590167837e-32 -J( 7, 4) = -2.055639042e-26 0 -J( 7, 5) = -2.141290668e-26 0 -J( 7, 6) = 1202310.628 1202310.628 -J( 7, 7) = -1244476.153 -1244476.153 -J( 7, 8) = 10.56270301 10.56270301 -J( 7, 9) = -2.312593922e-26 0 -J( 7, 10) = -2.335953456e-26 0 -J( 7, 11) = -2.355419735e-26 0 -J( 7, 12) = -2.371891202e-26 0 -J( 7, 13) = -2.386009602e-26 0 -J( 7, 14) = -2.741277362e-43 -2.741278343e-43 +J( S32, He4) = -3580.328682 -3712.987952 +J( S32, C12) = 8.04608259 0 +J( S32, O16) = -10.40027404 0.002248083797 +J( S32, Ne20) = -22.96819788 0 +J( S32, Mg24) = -32.64411602 0 +J( S32, Si28) = 96174.2781 96214.89322 +J( S32, S32) = -99978.56511 -99930.99407 +J( S32, Ar36) = -42.78877369 11.02602735 +J( S32, Ca40) = -59.590923 0 +J( S32, Ti44) = -64.96182659 0 +J( S32, Cr48) = -70.09376854 0 +J( S32, Fe52) = -75.06398858 0 +J( S32, Ni56) = -79.83062948 0 +J( S32, e) = -1.735789154e-15 -1.735789442e-15 -J( 8, 1) = 5.935462112e-25 5.877272628e-25 -J( 8, 2) = 4.909905946e-26 0 -J( 8, 3) = 5.523644189e-26 0 -J( 8, 4) = 5.891887135e-26 0 -J( 8, 5) = 6.137382433e-26 0 -J( 8, 6) = 6.312736216e-26 0 -J( 8, 7) = 1400032.306 1400032.306 -J( 8, 8) = -812316.926 -812316.926 -J( 8, 9) = 1.524756764 1.524756764 -J( 8, 10) = 6.69532629e-26 0 -J( 8, 11) = 6.751120676e-26 0 -J( 8, 12) = 6.79833131e-26 0 -J( 8, 13) = 6.838797568e-26 0 -J( 8, 14) = 7.857068531e-43 7.857069044e-43 +J(Ar36, He4) = 46086.53162 46453.96399 +J(Ar36, C12) = -24.11889018 0 +J(Ar36, O16) = 26.26894988 0 +J(Ar36, Ne20) = 60.49242539 0 +J(Ar36, Mg24) = 86.63689247 0 +J(Ar36, Si28) = 108.1539016 0 +J(Ar36, S32) = 112545.6504 112418.8663 +J(Ar36, Ar36) = -65833.81509 -65977.3066 +J(Ar36, Ca40) = 160.4642177 1.609672543 +J(Ar36, Ti44) = 173.235213 0 +J(Ar36, Cr48) = 186.8760509 0 +J(Ar36, Fe52) = 199.9634644 0 +J(Ar36, Ni56) = 212.5928074 0 +J(Ar36, e) = 4.796861622e-15 4.796861932e-15 -J( 9, 1) = 8.75876269e-25 8.591697512e-25 -J( 9, 2) = 1.000399222e-25 0 -J( 9, 3) = 1.125449125e-25 0 -J( 9, 4) = 1.200479066e-25 0 -J( 9, 5) = 1.250499028e-25 0 -J( 9, 6) = 1.286227571e-25 0 -J( 9, 7) = 1.313023979e-25 0 -J( 9, 8) = 902561.1588 902561.1588 -J( 9, 9) = -43393.10183 -43393.10183 -J( 9, 10) = 136.1591093 136.1591093 -J( 9, 11) = 1.37554893e-25 0 -J( 9, 12) = 1.385168154e-25 0 -J( 9, 13) = 1.393413202e-25 0 -J( 9, 14) = 1.600887132e-42 1.600887039e-42 +J(Ca40, He4) = 68936.76349 69694.52161 +J(Ca40, C12) = -21.37672362 0 +J(Ca40, O16) = 92.60590241 0 +J(Ca40, Ne20) = 172.688813 0 +J(Ca40, Mg24) = 235.8250397 0 +J(Ca40, Si28) = 289.279374 0 +J(Ca40, S32) = 336.6818312 0 +J(Ca40, Ar36) = 73674.38243 73294.33591 +J(Ca40, Ca40) = -3181.015849 -3601.602832 +J(Ca40, Ti44) = 605.9215887 146.8474287 +J(Ca40, Cr48) = 496.0195798 0 +J(Ca40, Fe52) = 531.7806992 0 +J(Ca40, Ni56) = 566.6106385 0 +J(Ca40, e) = 1.005604098e-14 1.00560404e-14 -J( 10, 1) = 3.341465689e-26 3.256971428e-26 -J( 10, 2) = 3.653921958e-27 0 -J( 10, 3) = 4.110662202e-27 0 -J( 10, 4) = 4.384706349e-27 0 -J( 10, 5) = 4.567402447e-27 0 -J( 10, 6) = 4.69789966e-27 0 -J( 10, 7) = 4.79577257e-27 0 -J( 10, 8) = 4.871895944e-27 0 -J( 10, 9) = 47730.54842 47730.54842 -J( 10, 10) = -15310.60916 -15310.60916 -J( 10, 11) = 0.002383888334 0.002383888334 -J( 10, 12) = 5.059276557e-27 0 -J( 10, 13) = 5.089391298e-27 0 -J( 10, 14) = 5.847182318e-44 5.847181188e-44 +J(Ti44, He4) = 2668.743332 2690.820423 +J(Ti44, C12) = 1.238968178 0 +J(Ti44, O16) = 5.219659626 0 +J(Ti44, Ne20) = 8.176705476 0 +J(Ti44, Mg24) = 10.62185414 0 +J(Ti44, Si28) = 12.77460702 0 +J(Ti44, S32) = 14.74456604 0 +J(Ti44, Ar36) = 16.5926625 0 +J(Ti44, Ca40) = 3978.151195 3959.795738 +J(Ti44, Ti44) = -1410.451232 -1430.507487 +J(Ti44, Cr48) = 21.71303639 0.002593936486 +J(Ti44, Fe52) = 23.32886523 0 +J(Ti44, Ni56) = 24.91917548 0 +J(Ti44, e) = 3.036977677e-16 3.036977135e-16 -J( 11, 1) = 7.727968776e-27 7.52140159e-27 -J( 11, 2) = 1.166347238e-27 0 -J( 11, 3) = 1.312140643e-27 0 -J( 11, 4) = 1.399616686e-27 0 -J( 11, 5) = 1.457934048e-27 0 -J( 11, 6) = 1.499589306e-27 0 -J( 11, 7) = 1.53083075e-27 0 -J( 11, 8) = 1.555129651e-27 0 -J( 11, 9) = 1.574568772e-27 0 -J( 11, 10) = 16539.09179 16539.09179 -J( 11, 11) = -9017.692804 -9017.692804 -J( 11, 12) = 0.000548515294 0.000548515294 -J( 11, 13) = 1.624555082e-27 0 -J( 11, 14) = 1.866445159e-44 1.866444877e-44 +J(Cr48, He4) = 614.7614844 623.7272557 +J(Cr48, C12) = -0.1747033379 0 +J(Cr48, O16) = 1.201593305 0 +J(Cr48, Ne20) = 2.1753907 0 +J(Cr48, Mg24) = 2.947870284 0 +J(Cr48, Si28) = 3.605375763 0 +J(Cr48, S32) = 4.190958113 0 +J(Cr48, Ar36) = 4.728616507 0 +J(Cr48, Ca40) = 5.232734331 0 +J(Cr48, Ti44) = 1390.049151 1384.336707 +J(Cr48, Cr48) = -754.4384225 -760.6122813 +J(Cr48, Fe52) = 6.621805499 0.0006014485747 +J(Cr48, Ni56) = 7.057441608 0 +J(Cr48, e) = 1.19432887e-16 1.194328695e-16 -J( 12, 1) = 7.899733031e-27 7.648455578e-27 -J( 12, 2) = 1.384355059e-27 0 -J( 12, 3) = 1.557399442e-27 0 -J( 12, 4) = 1.661226071e-27 0 -J( 12, 5) = 1.730443824e-27 0 -J( 12, 6) = 1.779885076e-27 0 -J( 12, 7) = 1.816966015e-27 0 -J( 12, 8) = 1.845806746e-27 0 -J( 12, 9) = 1.86887933e-27 0 -J( 12, 10) = 1.887756899e-27 0 -J( 12, 11) = 9769.164387 9769.164387 -J( 12, 12) = -2120.709403 -2120.709403 -J( 12, 13) = 0.0001045882492 0.0001045882492 -J( 12, 14) = 2.215311799e-44 2.215311465e-44 +J(Fe52, He4) = 632.9421872 643.7409355 +J(Fe52, C12) = -0.1711193576 0 +J(Fe52, O16) = 1.500523358 0 +J(Fe52, Ne20) = 2.686529769 0 +J(Fe52, Mg24) = 3.629742825 0 +J(Fe52, Si28) = 4.434187596 0 +J(Fe52, S32) = 5.15193241 0 +J(Fe52, Ar36) = 5.811846257 0 +J(Fe52, Ca40) = 6.431285246 0 +J(Fe52, Ti44) = 7.02131108 0 +J(Fe52, Cr48) = 831.5828167 823.9935726 +J(Fe52, Fe52) = -172.1130683 -180.2532886 +J(Fe52, Ni56) = 8.67797769 0.0001155651353 +J(Fe52, e) = 1.4407613e-16 1.440761083e-16 -J( 13, 1) = 2.36704737e-27 2.283840256e-27 -J( 13, 2) = 4.518871185e-28 0 -J( 13, 3) = 5.083730083e-28 0 -J( 13, 4) = 5.422645422e-28 0 -J( 13, 5) = 5.648588981e-28 0 -J( 13, 6) = 5.809977238e-28 0 -J( 13, 7) = 5.93101843e-28 0 -J( 13, 8) = 6.02516158e-28 0 -J( 13, 9) = 6.1004761e-28 0 -J( 13, 10) = 6.16209707e-28 0 -J( 13, 11) = 6.213447879e-28 0 -J( 13, 12) = 2283.840256 2283.840256 -J( 13, 13) = -0.0001126334991 -0.0001126334991 -J( 13, 14) = 7.231315828e-45 7.231315425e-45 +J(Ni56, He4) = 190.5614891 194.1182245 +J(Ni56, C12) = -0.04871014481 0 +J(Ni56, O16) = 0.5045792589 0 +J(Ni56, Ne20) = 0.8977755497 0 +J(Ni56, Mg24) = 1.210925284 0 +J(Ni56, Si28) = 1.478328594 0 +J(Ni56, S32) = 1.717144845 0 +J(Ni56, Ar36) = 1.936906672 0 +J(Ni56, Ca40) = 2.143324048 0 +J(Ni56, Ti44) = 2.340035525 0 +J(Ni56, Cr48) = 2.52947455 0 +J(Ni56, Fe52) = 196.8315355 194.1182245 +J(Ni56, Ni56) = 2.892625537 -0.0001244547611 +J(Ni56, e) = 4.749774962e-17 4.749774706e-17 -J( 14, 1) = 1.335535433e+19 1.331810907e+19 -J( 14, 2) = 1.672678172e+22 1.221628179e+22 -J( 14, 3) = 1.764169333e+23 1.669883326e+23 -J( 14, 4) = 6.38519349e+24 6.382691446e+24 -J( 14, 5) = 3.575435081e+24 3.56443165e+24 -J( 14, 6) = 2.522960852e+23 2.518972785e+23 -J( 14, 7) = 2.567887496e+23 2.492143743e+23 -J( 14, 8) = 1.63734141e+23 1.532779936e+23 -J( 14, 9) = 1.46197563e+22 5.363178614e+21 -J( 14, 10) = 1.41303877e+22 2.541660459e+21 -J( 14, 11) = 7.680250114e+21 1.43978686e+21 -J( 14, 12) = 6.72190453e+21 3.147552442e+20 -J( 14, 13) = -1.03705583e+17 -1.56260658e+13 -J( 14, 14) = -1.191469501 -1.191474366 +J( e, He4) = 8.557355646e+23 8.612933787e+23 +J( e, C12) = 3.253492329e+21 3.293956754e+21 +J( e, O16) = 1.407824852e+22 1.326089094e+22 +J( e, Ne20) = 5.099953876e+23 5.085260898e+23 +J( e, Mg24) = 2.852905494e+23 2.833327761e+23 +J( e, Si28) = 2.258348807e+22 2.015807656e+22 +J( e, S32) = 2.278872757e+22 2.001064883e+22 +J( e, Ar36) = 1.55712523e+22 1.244521857e+22 +J( e, Ca40) = 3.918079028e+21 4.446572283e+20 +J( e, Ti44) = 3.973002198e+21 1.959966253e+20 +J( e, Cr48) = 4.23547637e+21 1.21440394e+20 +J( e, Fe52) = 4.413681758e+21 2.675297734e+19 +J( e, Ni56) = 4.682037199e+21 -1.725529485e+13 +J( e, e) = 74499.43446 74499.57032 -AMReX (24.07-16-gdcb9cc0383dc) finalized +AMReX (24.12-14-gb3f67385e62f) finalized diff --git a/unit_test/jac_cell/inputs b/unit_test/jac_cell/inputs new file mode 100644 index 0000000000..4da177146d --- /dev/null +++ b/unit_test/jac_cell/inputs @@ -0,0 +1,7 @@ +unit_test.small_temp = 1.e5 +unit_test.small_dens = 1.e5 + +unit_test.density = 1.e6 +unit_test.temperature = 3.e9 + +unit_test.uniform_xn = 1 diff --git a/unit_test/jac_cell/inputs_aprox13 b/unit_test/jac_cell/inputs_aprox13 deleted file mode 100644 index 5e138fc92d..0000000000 --- a/unit_test/jac_cell/inputs_aprox13 +++ /dev/null @@ -1,38 +0,0 @@ -unit_test.run_prefix = "react_aprox13_" - -unit_test.small_temp = 1.e5 -unit_test.small_dens = 1.e5 - -integrator.burner_verbose = 0 - -# Set which jacobian to use -# 1 = analytic jacobian -# 2 = numerical jacobian - -integrator.jacobian = 2 - -integrator.renormalize_abundances = 0 - -integrator.rtol_spec = 1.0e-6 -integrator.rtol_enuc = 1.0e-6 -integrator.atol_spec = 1.0e-6 -integrator.atol_enuc = 1.0e-6 - -unit_test.tmax = 1.e-2 - -unit_test.density = 1.e6 -unit_test.temperature = 3.e9 - -unit_test.X1 = 1.0 -unit_test.X2 = 0.0 -unit_test.X3 = 0.0 -unit_test.X4 = 0.0 -unit_test.X5 = 0.0 -unit_test.X6 = 0.0 -unit_test.X7 = 0.0 -unit_test.X8 = 0.0 -unit_test.X9 = 0.0 -unit_test.X10 = 0.0 -unit_test.X11 = 0.0 -unit_test.X12 = 0.0 -unit_test.X13 = 0.0 diff --git a/unit_test/jac_cell/jac_cell.H b/unit_test/jac_cell/jac_cell.H index fd49859883..54560a56dd 100644 --- a/unit_test/jac_cell/jac_cell.H +++ b/unit_test/jac_cell/jac_cell.H @@ -1,5 +1,5 @@ -#ifndef BURN_CELL_H -#define BURN_CELL_H +#ifndef JAC_CELL_H +#define JAC_CELL_H #include #include @@ -22,7 +22,7 @@ void jac_cell_c() // Make sure user set all the mass fractions to values in the interval [0, 1] for (int n = 1; n <= NumSpec; ++n) { - massfractions[n-1] = get_xn(n); + massfractions[n-1] = get_xn(n, unit_test_rp::uniform_xn); if (massfractions[n-1] < 0 || massfractions[n-1] > 1) { amrex::Error("mass fraction for " + short_spec_names_cxx[n-1] + " not initialized in the interval [0,1]!"); @@ -38,7 +38,9 @@ void jac_cell_c() burn_state.xn[n] = massfractions[n]; } - normalize_abundances_burn(burn_state); + if (! unit_test_rp::skip_initial_normalization) { + normalize_abundances_burn(burn_state); + } // get the energy @@ -65,9 +67,15 @@ void jac_cell_c() // output + // header + std::cout << std::setw(16) << "element" << std::setw(20) + << "numerical" << std::setw(20) << "analytic" << std::endl; + for (int ii = 1; ii <= neqs; ++ii) { + std::string ilabel = (ii < neqs) ? short_spec_names_cxx[ii-1] : "e"; for (int jj = 1; jj <= neqs; ++jj) { - std::cout << "J(" << std::setw(3) << ii << ", " << std::setw(3) << jj << ") = " + std::string jlabel = (jj < neqs) ? short_spec_names_cxx[jj-1] : "e"; + std::cout << "J(" << std::setw(4) << ilabel << ", " << std::setw(4) << jlabel << ") = " << std::setw(20) << jac_numerical(ii,jj) << " " << std::setw(20) << jac_analytic(ii,jj) << std::endl; } diff --git a/unit_test/test_nse_net/GNUmakefile b/unit_test/nse_net_cell/GNUmakefile similarity index 100% rename from unit_test/test_nse_net/GNUmakefile rename to unit_test/nse_net_cell/GNUmakefile diff --git a/unit_test/test_nse_net/Make.package b/unit_test/nse_net_cell/Make.package similarity index 100% rename from unit_test/test_nse_net/Make.package rename to unit_test/nse_net_cell/Make.package diff --git a/unit_test/test_nse_net/README.md b/unit_test/nse_net_cell/README.md similarity index 80% rename from unit_test/test_nse_net/README.md rename to unit_test/nse_net_cell/README.md index 098db51823..62f7308040 100644 --- a/unit_test/test_nse_net/README.md +++ b/unit_test/nse_net_cell/README.md @@ -1,6 +1,6 @@ -# test_nse_net +# `nse_net_cell` -`test_nse_net` finds the NSE mass fractions of a given state. +`nse_net_cell` finds the NSE mass fractions of a given state. The density, temperature, and composition are set in the inputs file. Then we update the NSE mass fraction to the current state to make sure we're in NSE. Then we @@ -9,4 +9,4 @@ we're in NSE or not. Upon completion, the new state is printed to the screen. And a statement is printed to see whether we're in NSE -or not. \ No newline at end of file +or not. diff --git a/unit_test/test_nse_net/_parameters b/unit_test/nse_net_cell/_parameters similarity index 100% rename from unit_test/test_nse_net/_parameters rename to unit_test/nse_net_cell/_parameters diff --git a/unit_test/test_nse_net/burn_cell.H b/unit_test/nse_net_cell/burn_cell.H similarity index 100% rename from unit_test/test_nse_net/burn_cell.H rename to unit_test/nse_net_cell/burn_cell.H diff --git a/unit_test/test_nse_net/ci-benchmarks/nse_net_unit_test.out b/unit_test/nse_net_cell/ci-benchmarks/nse_net_unit_test.out similarity index 100% rename from unit_test/test_nse_net/ci-benchmarks/nse_net_unit_test.out rename to unit_test/nse_net_cell/ci-benchmarks/nse_net_unit_test.out diff --git a/unit_test/test_nse_net/inputs_ase b/unit_test/nse_net_cell/inputs_ase similarity index 100% rename from unit_test/test_nse_net/inputs_ase rename to unit_test/nse_net_cell/inputs_ase diff --git a/unit_test/test_nse_net/main.cpp b/unit_test/nse_net_cell/main.cpp similarity index 100% rename from unit_test/test_nse_net/main.cpp rename to unit_test/nse_net_cell/main.cpp diff --git a/unit_test/test_nse_net/make_table/GNUmakefile b/unit_test/nse_net_cell/make_table/GNUmakefile similarity index 100% rename from unit_test/test_nse_net/make_table/GNUmakefile rename to unit_test/nse_net_cell/make_table/GNUmakefile diff --git a/unit_test/test_nse_net/make_table/Make.package b/unit_test/nse_net_cell/make_table/Make.package similarity index 100% rename from unit_test/test_nse_net/make_table/Make.package rename to unit_test/nse_net_cell/make_table/Make.package diff --git a/unit_test/test_nse_net/make_table/README.md b/unit_test/nse_net_cell/make_table/README.md similarity index 100% rename from unit_test/test_nse_net/make_table/README.md rename to unit_test/nse_net_cell/make_table/README.md diff --git a/unit_test/test_nse_net/make_table/_parameters b/unit_test/nse_net_cell/make_table/_parameters similarity index 100% rename from unit_test/test_nse_net/make_table/_parameters rename to unit_test/nse_net_cell/make_table/_parameters diff --git a/unit_test/test_nse_net/make_table/burn_cell.H b/unit_test/nse_net_cell/make_table/burn_cell.H similarity index 100% rename from unit_test/test_nse_net/make_table/burn_cell.H rename to unit_test/nse_net_cell/make_table/burn_cell.H diff --git a/unit_test/test_nse_net/make_table/ci-benchmarks/nse_net_make_table_unit_test.out b/unit_test/nse_net_cell/make_table/ci-benchmarks/nse_net_make_table_unit_test.out similarity index 100% rename from unit_test/test_nse_net/make_table/ci-benchmarks/nse_net_make_table_unit_test.out rename to unit_test/nse_net_cell/make_table/ci-benchmarks/nse_net_make_table_unit_test.out diff --git a/unit_test/test_nse_net/make_table/main.cpp b/unit_test/nse_net_cell/make_table/main.cpp similarity index 100% rename from unit_test/test_nse_net/make_table/main.cpp rename to unit_test/nse_net_cell/make_table/main.cpp diff --git a/unit_test/test_part_func/GNUmakefile b/unit_test/part_func_cell/GNUmakefile similarity index 100% rename from unit_test/test_part_func/GNUmakefile rename to unit_test/part_func_cell/GNUmakefile diff --git a/unit_test/test_part_func/Make.package b/unit_test/part_func_cell/Make.package similarity index 100% rename from unit_test/test_part_func/Make.package rename to unit_test/part_func_cell/Make.package diff --git a/unit_test/test_part_func/_parameters b/unit_test/part_func_cell/_parameters similarity index 100% rename from unit_test/test_part_func/_parameters rename to unit_test/part_func_cell/_parameters diff --git a/unit_test/test_part_func/ci-benchmarks/part_func.out b/unit_test/part_func_cell/ci-benchmarks/part_func.out similarity index 100% rename from unit_test/test_part_func/ci-benchmarks/part_func.out rename to unit_test/part_func_cell/ci-benchmarks/part_func.out diff --git a/unit_test/test_part_func/main.cpp b/unit_test/part_func_cell/main.cpp similarity index 100% rename from unit_test/test_part_func/main.cpp rename to unit_test/part_func_cell/main.cpp diff --git a/unit_test/test_part_func/pf_cell.H b/unit_test/part_func_cell/pf_cell.H similarity index 100% rename from unit_test/test_part_func/pf_cell.H rename to unit_test/part_func_cell/pf_cell.H diff --git a/unit_test/react_util.H b/unit_test/react_util.H index 16177f4a9c..7579b753e3 100644 --- a/unit_test/react_util.H +++ b/unit_test/react_util.H @@ -88,7 +88,7 @@ init_t setup_composition(const int nz) { AMREX_INLINE AMREX_GPU_HOST_DEVICE -void get_xn(const int k, const init_t cd, amrex::Real *xn_zone, int uniform_composition = 0) { +void get_xn(const int k, const init_t cd, amrex::Real *xn_zone, bool uniform_composition=false) { for (int n = 0; n < NumSpec; n++) { xn_zone[n] = 0.0_rt; @@ -167,10 +167,15 @@ void get_xn(const int k, const init_t cd, amrex::Real *xn_zone, int uniform_comp /// AMREX_INLINE AMREX_GPU_HOST_DEVICE -amrex::Real get_xn(const int index) { +amrex::Real get_xn(const int index, bool uniform_composition=false) { amrex::Real mass_fraction{}; + if (uniform_composition) { + mass_fraction = 1.0_rt / NumSpec; + return mass_fraction; + } + switch (index) { case 1: diff --git a/unit_test/test_react/README.md b/unit_test/test_react/README.md index cada1e5784..798b024def 100644 --- a/unit_test/test_react/README.md +++ b/unit_test/test_react/README.md @@ -3,40 +3,3 @@ This is a unit test that sets up a cube of data (rho, T, and X varying along dimensions) and calls the burner on it. You can specify the integrator via `INTEGRATOR_DIR` and the network via `NETWORK_DIR` in the `GNUmakefile` - -## CPU Status - -This table summarizes tests run with gfortran. - - -| network | VODE | VODE90 | BS | Rosenbrock | -|------------------------|-------------------------|-------------------------|-------------------------|---------------| -| aprox13 | works
(6:1517:2393) | works
(6:1517:23793) | works
(126:858:4678) | crashes | -| aprox19 | works
(27:152:8127) | works
(27:152:8127) | works
(72:297:2869) | runs too long | -| triple_alpha_plus_cago | works
(6:32:1950) | works
(6:32:1950) | works
(126:148:3044) | crashes | -| ignition_chamulak | works
(6:6:28) | works
(6:6:28) | works
(144:144:153) | works (252:252:252) | - - -## Running on GPUs with VODE - -To run a GPU test with the VODE integrator and aprox13, do: - -``` -make -j COMP=PGI USE_CUDA=TRUE AMREX_USE_CUDA=TRUE USE_GPU_PRAGMA=TRUE NETWORK_DIR=aprox13 EOS_DIR=helmholtz -``` - -Tested with PGI 18.10 and CUDA 9.2.148. - -To run in a different directory, the following files are necessary to -run in addition to the executable: - -- `helm_table.dat` (will be symlinked in at compile time) -- `inputs_aprox13` -- `probin.aprox13` -- `xin.aprox13` - -Then in an interactive session, (e.g. on Summit), do: - -``` -jsrun -n 1 -a 1 -g 1 ./[executable] inputs_aprox13 -``` diff --git a/unit_test/test_react/TODO b/unit_test/test_react/TODO deleted file mode 100644 index 02ceb37947..0000000000 --- a/unit_test/test_react/TODO +++ /dev/null @@ -1 +0,0 @@ -Note: this has no OpenMP or OpenACC directives yet \ No newline at end of file diff --git a/unit_test/test_react/_parameters b/unit_test/test_react/_parameters index 918582dd7b..6cfcdab221 100644 --- a/unit_test/test_react/_parameters +++ b/unit_test/test_react/_parameters @@ -5,8 +5,6 @@ dens_max real 1.e9 temp_min real 1.e6 temp_max real 1.e15 -uniform_xn int 0 - tmax real 0.1e0 small_temp real 1.e5 diff --git a/unit_test/test_react/test_network_combinations.py b/unit_test/test_react/test_network_combinations.py deleted file mode 100755 index eadd7fdc09..0000000000 --- a/unit_test/test_react/test_network_combinations.py +++ /dev/null @@ -1,198 +0,0 @@ -#!/usr/bin/env python3 - -from __future__ import print_function - -import itertools -import os -import subprocess -import sys - -executable = "main.Linux.gfortran.omp.exe" - -link_files = ["xin.rprox", - "gr0_3d.small", - "helm_table.dat"] - -link_files.append(executable) - -# this dictionary holds the parameters we want to set. For each key, -# we use a list to give all the possible values (even if there is only -# a single value). -aprox13_params = { - "dens_min": [1.e4, 1.e6], - "dens_max": [1.e8], - "temp_min": [5.e7], - "temp_max": [5.e9], - "tmax": [1.e-3], - "rtol_spec": [1.e-12, 1.e-8], - "atol_spec": [1.e-12, 1.e-8], - "jacobian": [1, 2] -} - -aprox19_params = { - "dens_min": [1.e5], - "dens_max": [5.e8], - "temp_min": [5.e7], - "temp_max": [5.e8], - "tmax": [1.e-9], - "rtol_spec": [1.e-8, 1.e-6], - "atol_spec": [1.e-8, 1.e-6], - "jacobian": [1, 2], - "call_eos_in_rhs": ["T", "F"] -} - -rprox_params = { - "dens_min": [1.e5], - "dens_max": [2.e6], - "temp_min": [1.e8], - "temp_max": [1.e9], - "tmax": [1.e-6], - "rtol_spec": [1.e-12, 1.e-8], - "atol_spec": [1.e-12, 1.e-8], - "jacobian": [1, 2], - "call_eos_in_rhs": ["T", "F"], - "renormalize_abundances": ["T", "F"] -} - - -params = rprox_params - -params_file = r""" -&PROBIN - - xin_file = "xin.rprox" - run_prefix = "react_rprox_" - - small_dens = 1.0 - - @@test-params@@ - -/ -""" - -# time out for a run, in seconds -TIMEOUT = 600 - -def run(command, stdin=False, outfile=None): - """ run a command in the unix shell """ - - sin = None - if stdin: sin = subprocess.PIPE - p0 = subprocess.Popen(command, stdin=sin, stdout=subprocess.PIPE, - stderr=subprocess.STDOUT, shell=True) - - stdout0 = p0.communicate(timeout=TIMEOUT) - if stdin: p0.stdin.close() - rc = p0.returncode - p0.stdout.close() - - if outfile is not None: - try: cf = open(outfile, "w") - except IOError: - sys.exit("ERROR: unable to open file for writing: {}".format(outfile)) - else: - for line in stdout0: - if line is not None: - cf.write(line.decode("ascii")) - cf.close() - - return stdout0[0], rc - - -def doit(): - - # itertools.product() will produce lists with every possible - # combination from the input lists. Here we use as the input lists - # the values from our dictionary. This magic comes from: - # http://stackoverflow.com/questions/3873654/combinations-from-dictionary-with-list-values-using-python - - combinations = [[{k: v} for (k, v) in zip(params.keys(), values)] - for values in itertools.product(*params.values())] - - - run_no = 0 - - top_dir = os.getcwd() - - outcomes = {} - - of = open("test-results.txt", "w") - - for c in combinations: - - cparams = {k: v for d in c for k, v in d.items()} - - # run this combination of test parameters - - # make the directory - run_no += 1 - odir = "{:02d}".format(run_no) - - print("running case {}...".format(run_no)) - - try: - os.mkdir(odir) - except: - sys.exit("unable to create directory") - - # copy the executable and support files - - for f in link_files: - try: - os.symlink(os.path.join(top_dir, f), - os.path.join(odir, os.path.basename(f))) - except: - sys.exit("unable to link file") - - # let's work in that directory - os.chdir(odir) - - # write the input file - infile = "inputs.{}".format(odir) - with open(infile, "w") as f: - for line in params_file.splitlines(): - if line.find("@@test-params@@") >= 0: - for k, v in cparams.items(): - f.write(" {} = {}\n".format(k, v)) - else: - f.write("{}\n".format(line)) - - # run - command = "./{} {}".format(executable, infile) - stdout, rc = run(command, outfile="run.out") - - os.chdir(top_dir) - - # log the result - if rc == 0: - result = "successful" - else: - result = "failed" - - outcomes[odir] = result - - # if we were successful, the last 4 lines are the summary of RHS evaluations - stats = "" - if rc == 0: - stats = stdout.splitlines()[-4:] - - # output summary - of.write("test {}: {}\n".format(odir, result)) - for k, v in sorted(cparams.items()): - of.write(" {} = {}\n".format(k, v)) - of.write("\n") - - for line in stats: - of.write("{}\n".format(line.decode("utf-8"))) - of.write("\n") - of.flush() - - of.close() - - for k, v in sorted(outcomes.items()): - print("{}: {}".format(k, v)) - - - -if __name__ == "__main__": - doit()