diff --git a/.clang-tidy b/.clang-tidy index 391149e1da..df06bf21df 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -7,12 +7,12 @@ Checks: > clang-diagnostic-*, cppcoreguidelines-*, -cppcoreguidelines-avoid-c-arrays, + -cppcoreguidelines-avoid-do-while, -cppcoreguidelines-avoid-magic-numbers, -cppcoreguidelines-avoid-non-const-global-variables, -cppcoreguidelines-init-variables, -cppcoreguidelines-interfaces-global-init, -cppcoreguidelines-macro-usage, - -cppcoreguidelines-no-malloc, -cppcoreguidelines-non-private-member-variables-in-classes, -cppcoreguidelines-owning-memory, -cppcoreguidelines-pro-*, @@ -20,10 +20,10 @@ Checks: > -misc-const-correctness, -misc-include-cleaner, -misc-non-private-member-variables-in-classes, + -misc-no-recursion, modernize-*, -modernize-avoid-c-arrays, -modernize-use-trailing-return-type, - -modernize-use-using, performance-*, -performance-avoid-endl, portability-*, diff --git a/.github/workflows/c-linter.yml b/.github/workflows/c-linter.yml index d575a2f236..a779df810a 100644 --- a/.github/workflows/c-linter.yml +++ b/.github/workflows/c-linter.yml @@ -2,6 +2,10 @@ name: cpp-linter on: [pull_request] +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }} + cancel-in-progress: true + jobs: cpp-linter: runs-on: ubuntu-latest diff --git a/.github/workflows/check_pr_branch.yml b/.github/workflows/check_pr_branch.yml new file mode 100644 index 0000000000..3543088221 --- /dev/null +++ b/.github/workflows/check_pr_branch.yml @@ -0,0 +1,20 @@ +name: check PR branch + +on: + pull_request: + types: + - opened + - synchronize + - reopened + - edited + +jobs: + check-PR-branch: + runs-on: ubuntu-latest + steps: + - name: PRs should not target main + run: | + if [[ "${{ github.base_ref }}" == "main" ]]; then + echo 'Pull requests must not be made against main. Please target development instead.' + exit 1 + fi diff --git a/.github/workflows/clang-tidy.yml b/.github/workflows/clang-tidy.yml new file mode 100644 index 0000000000..f712591a4b --- /dev/null +++ b/.github/workflows/clang-tidy.yml @@ -0,0 +1,36 @@ +name: "clang-tidy" + +on: [pull_request] + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }} + cancel-in-progress: true + +jobs: + clang_tidy: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Get submodules + run: | + git submodule update --init + cd external/Microphysics + git fetch; git checkout development + echo "MICROPHYSICS_HOME=$(pwd)" >> $GITHUB_ENV + cd ../amrex + git fetch; git checkout development + echo "AMREX_HOME=$(pwd)" >> $GITHUB_ENV + cd ../.. + + - name: Install dependencies + run: | + .github/workflows/dependencies_clang-tidy-apt-llvm.sh 17 + + - name: Compile flame_wave + run: | + echo $AMREX_HOME + echo $MICROPHYSICS_HOME + cd Exec/science/subchandra + make USE_MPI=FALSE USE_CLANG_TIDY=TRUE CLANG_TIDY=clang-tidy-17 CLANG_TIDY_WARN_ERROR=TRUE -j 4 + diff --git a/.github/workflows/compiler-warnings.yml b/.github/workflows/compiler-warnings.yml index e8092be3d3..bdc0c627ac 100644 --- a/.github/workflows/compiler-warnings.yml +++ b/.github/workflows/compiler-warnings.yml @@ -1,6 +1,11 @@ name: "compiler warnings" on: [pull_request] + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }} + cancel-in-progress: true + jobs: compiler_warnings: runs-on: ubuntu-latest diff --git a/.github/workflows/dependencies_clang-tidy-apt-llvm.sh b/.github/workflows/dependencies_clang-tidy-apt-llvm.sh new file mode 100755 index 0000000000..b64cd619eb --- /dev/null +++ b/.github/workflows/dependencies_clang-tidy-apt-llvm.sh @@ -0,0 +1,22 @@ +#!/usr/bin/env bash + +set -eu -o pipefail + +# `man apt.conf`: +# Number of retries to perform. If this is non-zero APT will retry +# failed files the given number of times. +echo 'Acquire::Retries "3";' | sudo tee /etc/apt/apt.conf.d/80-retries + +if [[ ! -f /etc/apt/trusted.gpg.d/apt.llvm.org.asc ]]; then + wget -qO- https://apt.llvm.org/llvm-snapshot.gpg.key | sudo tee /etc/apt/trusted.gpg.d/apt.llvm.org.asc +fi + +source /etc/os-release # set UBUNTU_CODENAME + +sudo add-apt-repository "deb http://apt.llvm.org/${UBUNTU_CODENAME}/ llvm-toolchain-${UBUNTU_CODENAME} main" +sudo add-apt-repository "deb http://apt.llvm.org/${UBUNTU_CODENAME}/ llvm-toolchain-${UBUNTU_CODENAME}-$1 main" + +sudo apt-get update + +sudo apt-get install -y --no-install-recommends \ + clang-tidy-$1 libomp-$1-dev diff --git a/.github/workflows/detonation-sdc-compare.yml b/.github/workflows/detonation-sdc-compare.yml index b8fe7e8568..46038aa534 100644 --- a/.github/workflows/detonation-sdc-compare.yml +++ b/.github/workflows/detonation-sdc-compare.yml @@ -1,6 +1,11 @@ name: detonation-sdc on: [pull_request] + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }} + cancel-in-progress: true + jobs: detonation-sdc: runs-on: ubuntu-latest diff --git a/.github/workflows/flame_wave-compare.yml b/.github/workflows/flame_wave-compare.yml index 46ad5b4285..11d21675a7 100644 --- a/.github/workflows/flame_wave-compare.yml +++ b/.github/workflows/flame_wave-compare.yml @@ -1,6 +1,11 @@ name: flame_wave on: [pull_request] + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }} + cancel-in-progress: true + jobs: flame_wave-2d: runs-on: ubuntu-latest @@ -31,7 +36,7 @@ jobs: - name: Run flame_wave run: | cd Exec/science/flame_wave - ./Castro2d.gnu.DEBUG.ex inputs_2d.testsuite max_step=2 castro.sum_interval=1 amr.plot_files_output=0 amr.checkpoint_files_output=0 + ./Castro2d.gnu.DEBUG.ex inputs_2d.testsuite max_step=2 castro.sum_interval=1 amr.plot_files_output=1 amr.checkpoint_files_output=0 - name: Check grid_diag.out run: | @@ -42,3 +47,10 @@ jobs: run: | cd Exec/science/flame_wave diff species_diag.out ci-benchmarks/species_diag.out + + - name: Check job_info parameters + run: | + cd Exec/science/flame_wave + line_no=`grep -n "Inputs File Parameters" plt00000/job_info | cut -d ":" -f 1` + tail -n +$(($line_no -1)) plt00000/job_info > tmp.txt + diff tmp.txt ci-benchmarks/job_info_params.txt diff --git a/.github/workflows/fsanitizer.yml b/.github/workflows/fsanitizer.yml new file mode 100644 index 0000000000..7267a62488 --- /dev/null +++ b/.github/workflows/fsanitizer.yml @@ -0,0 +1,37 @@ +name: fsanitizer + +on: [push, pull_request] + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-fsanitizer + cancel-in-progress: true + +jobs: + fsanitizer: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Get submodules + run: | + git submodule update --init + cd external/Microphysics + git fetch; git checkout development + cd ../amrex + git fetch; git checkout development + cd ../.. + + - name: Install dependencies + run: | + sudo apt-get update -y -qq + sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0 + + - name: Compile subchandra + run: | + cd Exec/science/subchandra + make FSANITIZER=TRUE USE_MPI=FALSE -j 4 + + - name: Run subchandra + run: | + cd Exec/science/subchandra + ./Castro2d.gnu.SMPLSDC.ex inputs_2d.N14.coarse amr.max_level=1 max_step=1 amr.n_cell=320 640 diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml index d5b3995fa3..5ea459fdc7 100644 --- a/.github/workflows/gh-pages.yml +++ b/.github/workflows/gh-pages.yml @@ -44,7 +44,7 @@ jobs: run: ./deploy_docs_action.sh - name: Deploy - uses: peaceiris/actions-gh-pages@v3 + uses: peaceiris/actions-gh-pages@v4 with: github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: ./out diff --git a/.github/workflows/gpu_action.yml b/.github/workflows/gpu_action.yml index 95d65496b0..3c6af431e1 100644 --- a/.github/workflows/gpu_action.yml +++ b/.github/workflows/gpu_action.yml @@ -1,6 +1,11 @@ name: GPU compilation on: [pull_request] + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }} + cancel-in-progress: true + jobs: gpu-compilation: runs-on: ubuntu-latest diff --git a/.github/workflows/style/check_trailing_whitespaces.sh b/.github/workflows/style/check_trailing_whitespaces.sh index b33b62cf5c..c7f921a8df 100755 --- a/.github/workflows/style/check_trailing_whitespaces.sh +++ b/.github/workflows/style/check_trailing_whitespaces.sh @@ -6,6 +6,7 @@ find . -type d \( -name .git \ -o -path ./paper \ -o -name build -o -name install \ -o -name tmp_build_dir -o -name tmp_install_dir \ + -o -name ci-benchmarks \ -o -path ./util/gcem \ \) -prune -o \ -type f \( \( -name "*.H" -o -name "*.h" -o -name "*.hh" -o -name "*.hpp" \ diff --git a/.github/workflows/wdmerger_collision-compare.yml b/.github/workflows/wdmerger_collision-compare.yml index 99551f09dd..f8707b0ad6 100644 --- a/.github/workflows/wdmerger_collision-compare.yml +++ b/.github/workflows/wdmerger_collision-compare.yml @@ -1,6 +1,11 @@ name: wdmerger_collision on: [pull_request] + +concurrency: + group: ${{ github.ref }}-${{ github.head_ref }}-${{ github.workflow }} + cancel-in-progress: true + jobs: wdmerger_collision-2d: runs-on: ubuntu-latest @@ -43,5 +48,5 @@ jobs: - name: Check the extrema run: | cd Exec/science/wdmerger - ../../../external/amrex/Tools/Plotfile/fextrema.gnu.ex plt00095 > fextrema.out + ../../../external/amrex/Tools/Plotfile/fextrema.gnu.ex plt00086 > fextrema.out diff fextrema.out ci-benchmarks/wdmerger_collision_2D.out diff --git a/CHANGES.md b/CHANGES.md index 7aeb11c743..b267ba4db1 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,42 @@ +# 24.05 + + * Changed how the shock flag is computed. It is now computed once, + at the start of a timestep. It also takes into account sources + such that the pressure jump is only the pressure that is available + to drive a shock. (#2726, #2728) + + * Fixed a boundary condition check in gravity (#2833) + + * Some coverity and clang-tidy fixes (#2810, #2811, #2835, #2831, + #2832) + + * Fix the wdmerger initialization when the two stars are overlapping + (#2829) + + * Fixed the Castro retry mechanism when amr.subcycling_mode = None + (#2821, #2824, #2826) + + * Added more amr parameters to the job_info file (#2828) + + * Added OpenMP to the SDC burn loop (#2770) + +# 24.04 + + * Some clang-tidy fixes (#2779, #2780, #2781, #2784, #2786, #2787, #2790, + #2791, #2792, #2793, #2797, #2798, #2799, #2800, #2801, #2804) + + * Fix species initialization in the Detonation problem (#2806) + and let it work with NSE (#2765) + + * Sync up with AMReX changes (#2794) + + * wdmerger now reports composition of initial stars (#2767) + + * flame_wave now checks if the `atm_delta` is too small (#2782) + and stores X(ash) (#2773) + + * a bounds issue in the true SDC integration was fixed (#2775) + # 24.03 * Documentation updates (#2742, #2752, #2753) diff --git a/Docs/source/debugging.rst b/Docs/source/debugging.rst index 4baf52c717..3341984d51 100644 --- a/Docs/source/debugging.rst +++ b/Docs/source/debugging.rst @@ -85,6 +85,12 @@ You can also use ask it to fix errors automatically via: and you can treat warnings as errors by adding ``CLANG_TIDY_WARN_ERROR=TRUE``. +.. note:: + + Building a Castro problem with ``clang-tidy`` will suppress the + checks in AMReX and Microphysics sources and headers. This is set by + the parameter ``CLANG_TIDY_IGNORE_SOURCES`` in ``Make.Castro``, and + the ``HeaderFilterRegex`` whitelist in ``.clang-tidy``. Thread sanitizer @@ -121,5 +127,3 @@ to operate over. Physics issues ============== - - diff --git a/Docs/source/io.rst b/Docs/source/io.rst index e9dc10993c..1e6191f6f4 100644 --- a/Docs/source/io.rst +++ b/Docs/source/io.rst @@ -304,6 +304,8 @@ Derived variables +-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+ | variable name | description | derive routine | units | +===================================+===================================================+=============================+=========================================+ +| ``abar`` | Mean atomic mass | ``derabar`` | :math:`\amu` | ++-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+ | ``angular_momentum_x``, | Angular momentum / volume in the x, y, or z dir | ``derangmomx``, | :math:`{\rm g~cm^{-1}~s^{-1}}` | | ``angular_momentum_y``, | computed as :math:`[(\rho \ub) \times {\bf r}]_n` | ``derangmomy``, | | | ``angular_momentum_z`` | where :math:`{\bf r}` is the distance from | ``derangmomz`` | | @@ -328,6 +330,8 @@ Derived variables | | :math:`s = s(\rho, e, X_k)`, where `e` is | | | | | computed from :math:`(\rho e)` | | | +-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+ +| ``enuc`` | Nuclear energy generation rate / gram | ``derenuc`` | :math:`{\rm erg~g^{-1}~s^{-1}}` | ++-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+ | ``Ertot`` | Total radiation energy density | ``derertot`` | | | | (for multigroup radiation problems) | | | +-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+ @@ -391,8 +395,6 @@ Derived variables | ``y_velocity``, | :math:`\ub = (\rho \ub)/\rho` | | | | ``z_velocity`` | | | | +-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+ -| ``enuc`` | Nuclear energy generation rate / gram | ``derenuc`` | :math:`{\rm erg~g^{-1}~s^{-1}}` | -+-----------------------------------+---------------------------------------------------+-----------------------------+-----------------------------------------+ problem-specific plotfile variables ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/Docs/source/reactions.rst b/Docs/source/reactions.rst index ed23b9d746..25df19d9d9 100644 --- a/Docs/source/reactions.rst +++ b/Docs/source/reactions.rst @@ -11,10 +11,17 @@ and it evolves those components through a nuclear burning step. All of the reaction networks that Castro uses are provided by the `Microphysics repository `_. +.. index:: USE_REACT + .. note:: - An arbitrary reaction network can be created for Castro via the - `pynucastro library `_. + To enable reactions in a simulation, you must compile with + + :: + + USE_REACT=TRUE + + This can be set in your ``GNUmakefile``. Microphysics comes with a ``general_null`` network. This is a bare interface for a @@ -24,6 +31,16 @@ are accepted. It contains several sets of isotopes; for example, isotopes needed to represent the triple-\ :math:`\alpha` reaction converting helium into carbon, as well as oxygen and iron. +Other reaction networks can be found in ``Microphysics/networks``. To select +a network, you set the make variable ``NETWORK_DIR`` to the name of the network +directory. + +.. note:: + + An arbitrary reaction network can be created for Castro via the + `pynucastro library `_. + + The main interface for burning is in ``Microphysics/interfaces/burner.H``: .. code:: c++ @@ -46,7 +63,7 @@ details. Controlling burning =================== -.. index:: castro.react_T_min, castro.react_T_max, castro.react_rho_min, castro.react_rho_max +.. index:: castro.react_T_min, castro.react_T_max, castro.react_rho_min, castro.react_rho_max, castro.do_react There are a number of reactions-related parameters that can be set at runtime in the inputs file. Reactions are enabled by setting:: @@ -63,7 +80,11 @@ reactions to occur in a zone using the parameters: * ``castro.react_rho_min`` and ``castro.react_rho_max`` for density -.. index:: castro.disable_shock_burning, USE_SHOCK_VAR + +Burning in Shocks +----------------- + +.. index:: USE_SHOCK_VAR, castro.diable_shock_burning, castro.shock_detection_threshold, castro.shock_detection_include_sources Burning can also be disabled inside shocks. This requires that the code be compiled with:: @@ -72,12 +93,27 @@ compiled with:: in the ``GNUmakefile``. This will allocate storage for a shock flag in the conserved state array. This flag is computed via a multidimensional shock detection algorithm -that looks for compression (:math:`\nabla \cdot \ub < 0`) along with a pressure jump -in the direction of compression. The runtime parameter:: +described in :cite:`doubledet2024`. A zone is tagged as a shock if the following +conditions are true: + +.. math:: + + \begin{align*} + \nabla \cdot \ub &< 0 \\ + \frac{|(\nabla p - \rho {\bf g}) \cdot \ub|}{p |\ub_\mathrm{cell}|} > &f_\mathrm{shock} + \end{align*} + +This requires that there is compression and that the pressure jump (excluding +the part of the pressure that balances gravity) is large. The runtime parameter + +:: castro.disable_shock_burning = 1 -will skip reactions in a zone where we've detected a shock. +will skip reactions in a zone where we've detected a shock. The runtime parameters +``castro.shock_detection_threshold`` and ``castro.shock_detection_include_sources`` +will set the value of $f_\mathrm{shock}$ and whether to subtract $\rho {\bf g}$ +from the pressure gradient. .. note:: diff --git a/Docs/source/refs.bib b/Docs/source/refs.bib index 612f332637..7683e16623 100644 --- a/Docs/source/refs.bib +++ b/Docs/source/refs.bib @@ -1119,4 +1119,21 @@ @article{castro_simple_sdc title = {An Improved Method for Coupling Hydrodynamics with Astrophysical Reaction Networks}, journal = {The Astrophysical Journal}, abstract = {Reacting astrophysical flows can be challenging to model, because of the difficulty in accurately coupling hydrodynamics and reactions. This can be particularly acute during explosive burning or at high temperatures where nuclear statistical equilibrium is established. We develop a new approach, based on the ideas of spectral deferred corrections (SDC) coupling of explicit hydrodynamics and stiff reaction sources as an alternative to operator splitting, that is simpler than the more comprehensive SDC approach we demonstrated previously. We apply the new method to a double-detonation problem with a moderately sized astrophysical nuclear reaction network and explore the time step size and reaction network tolerances, to show that the simplified-SDC approach provides improved coupling with decreased computational expense compared to traditional Strang operator splitting. This is all done in the framework of the Castro hydrodynamics code, and all algorithm implementations are freely available.} -} \ No newline at end of file +} + +@ARTICLE{doubledet2024, + author = {{Zingale}, Michael and {Chen}, Zhi and {Rasmussen}, Melissa and {Polin}, Abigail and {Katz}, Max and {Smith Clark}, Alexander and {Johnson}, Eric T.}, + title = "{Sensitivity of Simulations of Double Detonation Type Ia Supernova to Integration Methodology}", + journal = {arXiv e-prints}, + keywords = {Astrophysics - High Energy Astrophysical Phenomena, Astrophysics - Solar and Stellar Astrophysics}, + year = 2023, + month = sep, + eid = {arXiv:2309.01802}, + pages = {arXiv:2309.01802}, + doi = {10.48550/arXiv.2309.01802}, +archivePrefix = {arXiv}, + eprint = {2309.01802}, + primaryClass = {astro-ph.HE}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2023arXiv230901802Z}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} diff --git a/Exec/Make.Castro b/Exec/Make.Castro index 3a8ffa20c8..342f5fe97d 100644 --- a/Exec/Make.Castro +++ b/Exec/Make.Castro @@ -357,7 +357,7 @@ include $(AMREX_HOME)/Tools/GNUMake/Make.rules clean:: $(SILENT) $(RM) extern.F90 prob_params_auto.F90 extern_parameters.H extern_parameters_F.H extern_parameters.cpp $(SILENT) $(RM) AMReX_buildInfo.cpp - $(SILENT) $(RM) $(CASTRO_AUTO_SOURCE_DIR)/*.H $(CASTRO_AUTO_SOURCE_DIR)/*.[fF]90 + $(SILENT) $(RM) $(CASTRO_AUTO_SOURCE_DIR)/*.H $(CASTRO_AUTO_SOURCE_DIR)/*.[fF]90 $(CASTRO_AUTO_SOURCE_DIR)/*.cpp # these files are now created directly in the CASTRO_AUTO_SOURCE_DIR so eventually we # can get rid of explicitly removing them (this is for backwards compatibility) diff --git a/Exec/Make.auto_source b/Exec/Make.auto_source index 178e0f43e2..596542809d 100644 --- a/Exec/Make.auto_source +++ b/Exec/Make.auto_source @@ -67,6 +67,8 @@ AUTO_BUILD_SOURCES += $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H CPP_PARAMETERS := $(TOP)/Source/driver/_cpp_parameters +$(CASTRO_AUTO_SOURCE_DIR)/runtime_params.cpp: $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H + $(CASTRO_AUTO_SOURCE_DIR)/castro_params.H: $(CPP_PARAMETERS) @if [ ! -d $(CASTRO_AUTO_SOURCE_DIR) ]; then mkdir -p $(CASTRO_AUTO_SOURCE_DIR); fi PYTHONPATH=$(MICROPHYSICS_HOME)/util/build_scripts $(TOP)/Source/driver/parse_castro_params.py -o $(CASTRO_AUTO_SOURCE_DIR) $(CPP_PARAMETERS) diff --git a/Exec/gravity_tests/hse_convergence/_prob_params b/Exec/gravity_tests/hse_convergence/_prob_params index 07cd760aaf..f2a36e386b 100644 --- a/Exec/gravity_tests/hse_convergence/_prob_params +++ b/Exec/gravity_tests/hse_convergence/_prob_params @@ -5,4 +5,4 @@ temp_base real 1.0_rt y pert_width real 1.0_rt y -do_pert logical .true. y +do_pert bool true y diff --git a/Exec/hydro_tests/Sedov/ci-benchmarks/grid_diag.out b/Exec/hydro_tests/Sedov/ci-benchmarks/grid_diag.out index d27e6e4892..61946ff5b3 100644 --- a/Exec/hydro_tests/Sedov/ci-benchmarks/grid_diag.out +++ b/Exec/hydro_tests/Sedov/ci-benchmarks/grid_diag.out @@ -1,6 +1,6 @@ # COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 # TIMESTEP TIME MASS XMOM YMOM ZMOM ANG. MOM. X ANG. MOM. Y ANG. MOM. Z KIN. ENERGY INT. ENERGY GAS ENERGY CENTER OF MASS X-LOC CENTER OF MASS Y-LOC CENTER OF MASS Z-LOC CENTER OF MASS X-VEL CENTER OF MASS Y-VEL CENTER OF MASS Z-VEL MAXIMUM TEMPERATURE MAXIMUM DENSITY - 0 0.0000000000000000 1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.9540123498075317e-01 9.9540123498075317e-01 5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.6844909478137644e-03 1.0000000000000000e+00 + 0 0.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 9.9540123498075317e-01 9.9540123498075317e-01 5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.6844909478137644e-03 1.0000000000000000e+00 1 3.2340315269509180e-07 1.0000000000000000e+00 -1.2207196712653843e-22 -2.4091738146955404e-21 1.5915318696162594e-21 -6.8657515920158759e-26 -6.1490168351429670e-25 3.2333634833518791e-26 8.9703654836411012e-05 9.9531153132591343e-01 9.9540123498075073e-01 5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01 -1.2207196712653843e-22 -2.4091738146955404e-21 1.5915318696162594e-21 1.6858365601187205e-03 1.0027049946567383e+00 2 6.7914662065969288e-07 1.0000000000000000e+00 -4.8211172704689422e-22 -1.0918949986942682e-21 7.0508219992330420e-22 -1.8798232765775811e-24 8.0820625471129870e-25 -1.7707803610263417e-24 3.8423927161269466e-04 9.9501699570913438e-01 9.9540123498074762e-01 5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01 -4.8211172704689422e-22 -1.0918949986942682e-21 7.0508219992330420e-22 1.6864495062059968e-03 1.0056668439067880e+00 3 1.0704644354207541e-06 1.0000000000000000e+00 -2.9813388912628283e-22 -1.7291596040867426e-21 2.2238460006813763e-21 4.4807723826432593e-24 3.0419916181484195e-26 -1.1783993498209245e-25 9.2349575408761178e-04 9.9447773922665816e-01 9.9540123498074684e-01 5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01 -2.9813388912628283e-22 -1.7291596040867426e-21 2.2238460006813763e-21 1.6861149214870039e-03 1.0089335120072900e+00 diff --git a/Exec/hydro_tests/Sedov/problem_initialize.H b/Exec/hydro_tests/Sedov/problem_initialize.H index b056e41349..fb0eb76814 100644 --- a/Exec/hydro_tests/Sedov/problem_initialize.H +++ b/Exec/hydro_tests/Sedov/problem_initialize.H @@ -22,10 +22,7 @@ void problem_initialize () } } - Real xn_zone[NumSpec]; - for (int n = 0; n < NumSpec; ++n) { - xn_zone[n] = 0.0_rt; - } + Real xn_zone[NumSpec] = {0.0}; xn_zone[0] = 1.0_rt; eos_t eos_state; @@ -96,17 +93,20 @@ void problem_initialize () vctr = (4.0_rt / 3.0_rt) * M_PI * problem::r_init * problem::r_init * problem::r_init; #else - amrex::Abort("Sedov problem unsupported in 3D axisymmetric geometry."); - #endif } else if (coord_type == 2) { + +#if AMREX_SPACEDIM == 1 // Must have AMREX_SPACEDIM == 1 for this coord_type. vctr = (4.0_rt / 3.0_rt) * M_PI * problem::r_init * problem::r_init * problem::r_init; +#else + amrex::Abort("Sedov problem unsupported in 2-D or 3-D spherical geometry."); +#endif } diff --git a/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out b/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out index b23a288448..5d48f26c1b 100644 --- a/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out +++ b/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out @@ -1,79 +1,79 @@ plotfile = det_x_plt00040 - time = 5.1558159133685604e-06 + time = 5.1558159140336702e-06 variables minimum value maximum value - density 185259613.21 216639295.98 - xmom -76364812859 2.9487648003e+16 + density 185259659.3 216625503.34 + xmom -13760057465 2.9527146771e+16 ymom 0 0 zmom 0 0 - rho_E 1.3062473824e+26 2.7890630574e+26 - rho_e 1.3062473824e+26 2.7748033958e+26 - Temp 50000004.458 7845854782 - rho_H1 2.1207015241e-22 0.020000144873 - rho_He3 0.0017224909604 0.021023115362 - rho_He4 94365290.111 200001891.25 - rho_C12 0.020000000212 21243051.967 - rho_N14 1.9999977239e-22 0.020000190213 - rho_O16 0.019999999996 19209.485639 - rho_Ne20 0.019999999996 9734.6358293 - rho_Mg24 0.019999999996 23306.635461 - rho_Si28 0.019999999996 2013766.0467 - rho_S32 0.019999999996 1654854.387 - rho_Ar36 0.019999999996 821077.40202 - rho_Ca40 0.019999999996 724428.55057 - rho_Ti44 0.019999999996 34137.672926 - rho_Cr48 0.019999999996 78094.207516 - rho_Fe52 0.019999999996 277975.05929 - rho_Fe54 0.019999999996 94999519.274 - rho_Ni56 0.019999999996 2236723.6985 - rho_n 2.1207015241e-22 234580.251 - rho_p 0.019999995434 3609846.1745 - rho_enuc -4.6916352344e+29 3.5793978339e+32 - pressure 5.5236728673e+25 1.1610546577e+26 - kineng 0 2.0430835425e+24 - soundspeed 612864631.35 895226082.4 - Gamma_1 1.359975607 1.3820247131 - MachNumber 0 0.16100811998 - uplusc 612864631.35 999536286.21 - uminusc -895226453.03 -612859846.92 - entropy 98214771.47 336273518.24 + rho_E 1.3062473821e+26 2.7889438781e+26 + rho_e 1.3062473821e+26 2.7747059271e+26 + Temp 50000000.026 7845854559 + rho_H1 2.1203665799e-22 0.020000114043 + rho_He3 0.0017224849005 0.021021778892 + rho_He4 94384549.451 200001582.65 + rho_C12 0.020000000216 21181006.118 + rho_N14 1.9999997268e-22 0.020000159166 + rho_O16 0.02 19209.472033 + rho_Ne20 0.02 3403.302453 + rho_Mg24 0.02 23330.215568 + rho_Si28 0.02 2016074.1508 + rho_S32 0.02 1656865.76 + rho_Ar36 0.02 822126.88981 + rho_Ca40 0.02 725125.60248 + rho_Ti44 0.02 34184.806642 + rho_Cr48 0.02 78126.082381 + rho_Fe52 0.02 278040.09021 + rho_Fe54 0.02 94986843.638 + rho_Ni56 0.02 2239055.987 + rho_n 0.020000004562 234580.08802 + rho_p 0.019999995438 3609282.9073 + rho_enuc -4.6902991452e+29 3.5778382148e+32 + pressure 5.5236728651e+25 1.1610543806e+26 + kineng 0 2.047804488e+24 + soundspeed 612864631.21 895226062.8 + Gamma_1 1.3599756296 1.3819790602 + MachNumber 0 0.16114515958 + uplusc 612864631.21 999634014.46 + uminusc -895226128.14 -612860627.32 + entropy 98214767.758 336273492.62 magvort 0 0 - divu -97902.982158 33744.686211 - eint_E 6.5312369121e+17 1.3804344246e+18 - eint_e 6.5312369121e+17 1.3804344246e+18 - logden 8.267780753 8.3357372357 - StateErr_0 185259613.21 216639295.98 - StateErr_1 50000004.458 7845854782 - StateErr_2 1e-30 9.9999779304e-11 - X(H1) 1e-30 9.9999779304e-11 - X(He3) 8.9859771063e-12 9.9999601234e-11 - X(He4) 0.48406305162 0.9999999982 - X(C12) 1.0000000106e-10 0.10016992833 - X(N14) 1e-30 1.0000000525e-10 - X(O16) 9.999999998e-11 9.6047401986e-05 - X(Ne20) 9.999999998e-11 4.5902903915e-05 - X(Mg24) 9.999999998e-11 0.00011109993006 - X(Si28) 9.999999998e-11 0.010051430746 - X(S32) 9.999999998e-11 0.0083423100958 - X(Ar36) 9.999999998e-11 0.004162785426 - X(Ca40) 9.999999998e-11 0.0037083550771 - X(Ti44) 9.999999998e-11 0.00017290662768 - X(Cr48) 9.999999998e-11 0.00040059772358 - X(Fe52) 9.999999998e-11 0.001337522714 - X(Fe54) 9.999999998e-11 0.46488979858 - X(Ni56) 9.999999998e-11 0.010324644421 - X(n) 1e-30 0.0011729007238 - X(p) 9.9999977087e-11 0.017592929732 - abar 4.000000001 6.7312458561 - Ye 0.49998677577 0.50001556472 - x_velocity -381.82392939 138572160.27 + divu -97879.954761 33763.419428 + eint_E 6.5312369103e+17 1.3804343698e+18 + eint_e 6.5312369103e+17 1.3804343698e+18 + logden 8.267780861 8.3357095848 + StateErr_0 185259659.3 216625503.34 + StateErr_1 50000000.026 7845854559 + StateErr_2 1e-30 9.9999779324e-11 + X(H1) 1e-30 9.9999779324e-11 + X(He3) 8.9862145938e-12 9.9999601254e-11 + X(He4) 0.48407561714 0.9999999982 + X(C12) 1.0000000108e-10 0.099893133188 + X(N14) 1e-30 1.0000000432e-10 + X(O16) 1e-10 9.6047335509e-05 + X(Ne20) 1e-10 1.6050538078e-05 + X(Mg24) 1e-10 0.00011120690546 + X(Si28) 1e-10 0.010059247336 + X(S32) 1e-10 0.0083490446432 + X(Ar36) 1e-10 0.0041662673673 + X(Ca40) 1e-10 0.0037091533301 + X(Ti44) 1e-10 0.00017295101573 + X(Cr48) 1e-10 0.00040068985616 + X(Fe52) 1e-10 0.0013377863567 + X(Fe54) 1e-10 0.46486808681 + X(Ni56) 1e-10 0.010336068249 + X(n) 1.0000002281e-10 0.001172900139 + X(p) 9.9999977027e-11 0.017594951746 + abar 4.000000001 6.7310197958 + Ye 0.49998669648 0.50001557929 + x_velocity -68.800275809 138706560.7 y_velocity 0 0 z_velocity 0 0 - t_sound_t_enuc 3.4412406672e-13 0.97497430527 - enuc -2.4944881551e+21 1.6522384906e+24 - magvel 0 138572160.27 - radvel -381.82392939 138572160.27 - circvel 0 2 - magmom 0 2.9487648003e+16 + t_sound_t_enuc 3.4412406685e-13 0.97470686938 + enuc -2.4937801624e+21 1.6516237283e+24 + magvel 0 138706560.7 + radvel -68.800275809 138706560.7 + circvel 0 2.8284271247 + magmom 0 2.9527146771e+16 angular_momentum_x 0 0 angular_momentum_y 0 0 angular_momentum_z 0 0 diff --git a/Exec/science/Detonation/problem_initialize.H b/Exec/science/Detonation/problem_initialize.H index 58205bbf24..c625a16440 100644 --- a/Exec/science/Detonation/problem_initialize.H +++ b/Exec/science/Detonation/problem_initialize.H @@ -55,10 +55,15 @@ void problem_initialize () if (problem::in14 >= 0) { problem::xn[problem::in14] = amrex::max(problem::nfrac, problem::smallx); - problem::xn[problem::ihe4] = 1.0_rt - problem::cfrac - problem::nfrac - problem::ofrac - (NumSpec - 3) * problem::smallx; + problem::xn[problem::ihe4] = 1.0_rt - problem::xn[problem::ic12] + - problem::xn[problem::in14] + - problem::xn[problem::io16] + - (NumSpec - 4) * problem::smallx; } else { - problem::xn[problem::ihe4] = 1.0_rt - problem::cfrac - problem::ofrac - (NumSpec - 2) * problem::smallx; + problem::xn[problem::ihe4] = 1.0_rt - problem::xn[problem::ic12] + - problem::xn[problem::io16] + - (NumSpec - 3) * problem::smallx; } // Set the ambient material diff --git a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out index 2f054f110b..973d3f2aba 100644 --- a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out @@ -1,5 +1,5 @@ # COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 # TIMESTEP TIME MASS XMOM YMOM ZMOM ANG. MOM. X ANG. MOM. Y ANG. MOM. Z KIN. ENERGY INT. ENERGY GAS ENERGY GRAV. ENERGY TOTAL ENERGY CENTER OF MASS X-LOC CENTER OF MASS Y-LOC CENTER OF MASS Z-LOC CENTER OF MASS X-VEL CENTER OF MASS Y-VEL CENTER OF MASS Z-VEL MAXIMUM TEMPERATURE MAXIMUM DENSITY MAXIMUM T_S / T_E - 0 0.0000000000000000 1.4902024346464787e+20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.9583337617665660e+37 2.9583337617665660e+37 0.0000000000000000e+00 2.9583337617665660e+37 2.7306640614148750e+04 6.8067773166795473e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000156899e+08 3.1017131702472486e+07 0.0000000000000000e+00 - 1 5.9694928185132600e-08 1.4902038799260451e+20 -9.8852814000373391e+18 7.0758681407801342e+23 -3.7077151119947915e+15 -7.8559831403383900e+18 7.9789428550115738e+19 1.9321854546506743e+28 2.7622571688115688e+29 2.9583392229687520e+37 2.9583392505913144e+37 0.0000000000000000e+00 2.9583392505913144e+37 2.7306640614147251e+04 6.8067725884159722e+02 0.0000000000000000e+00 -6.6335093695554737e-02 4.7482550784469113e+03 -2.4880589575292177e-05 1.0000972452758998e+08 3.1017431302318867e+07 2.1834522208251888e-12 - 2 1.2535934918877847e-07 1.4902056479398969e+20 -2.0759055857352421e+19 1.4862415462814437e+24 -1.6351019861504854e+16 -3.4644556500407476e+19 3.5187152219157529e+20 4.0584338401784253e+28 5.5823502529036681e+29 2.9583452181108315e+37 2.9583452739343125e+37 0.0000000000000000e+00 2.9583452739343125e+37 2.7306640614140484e+04 6.8067697458182033e+02 0.0000000000000000e+00 -1.3930329606554864e-01 9.9733989623248763e+03 -1.0972324446702355e-04 1.0003863648118678e+08 3.1017667549506564e+07 1.8217100523183610e-12 + 0 0.0000000000000000e+00 1.4700685690736437e+20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.9163021935541997e+37 2.9163021935541997e+37 0.0000000000000000e+00 2.9163021935541997e+37 2.7306641645125415e+04 6.8087525965685256e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.3928678114307070e+09 3.0715027784567930e+07 0.0000000000000000e+00 + 1 5.9806047036494849e-08 1.4700700110495903e+20 3.3564221294760178e+23 7.2876252590129013e+23 1.2612378576902532e+20 3.1318392026636616e+23 -2.7141686393328650e+24 1.9066097490832445e+28 7.8605333893026802e+29 2.9163081863790971e+37 2.9163082649844082e+37 0.0000000000000000e+00 2.9163082649844082e+37 2.7306641717925890e+04 6.8087482283157362e+02 0.0000000000000000e+00 2.2831716205676648e+03 4.9573321027137572e+03 8.5794407627549896e-01 1.3929561430929687e+09 3.0715326643214259e+07 3.3873763068343888e-04 + 2 1.2559269877663918e-07 1.4700717750115456e+20 7.0484895627215180e+23 1.5304335700608418e+24 5.5619062084401319e+20 1.3810901029080806e+24 -1.1969287510263666e+25 4.0039912973182075e+28 3.0853296206485770e+30 2.9163145904134766e+37 2.9163148989464290e+37 0.0000000000000000e+00 2.9163148989464290e+37 2.7306641954244220e+04 6.8087461718998884e+02 0.0000000000000000e+00 4.7946567525018709e+03 1.0410604407725754e+04 3.7834249340624542e+00 1.3933944931271181e+09 3.0715562157661017e+07 3.3857586853705914e-04 diff --git a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt new file mode 100644 index 0000000000..bd1e776e0c --- /dev/null +++ b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt @@ -0,0 +1,258 @@ +============================================================================== + Inputs File Parameters +============================================================================== +[*] castro.diffuse_temp = 1 +[*] castro.diffuse_cutoff_density = 20000 +[*] castro.diffuse_cutoff_density_hi = 50000 +[*] castro.diffuse_cond_scale_fac = 10 + castro.use_point_mass = 0 + castro.point_mass = 0 + castro.point_mass_fix_solution = 0 + castro.gw_dist = 0 + castro.point_mass_offset_is_true = 0 + castro.point_mass_location_offset = 0 +[*] castro.rotational_period = 0.001 +[*] castro.rotation_include_centrifugal = 0 + castro.rotation_include_coriolis = 1 + castro.rot_source_type = 4 + castro.implicit_rotation_update = 1 +[*] castro.rot_axis = 2 + castro.sponge_lower_radius = -1 + castro.sponge_upper_radius = -1 +[*] castro.sponge_lower_density = 1 +[*] castro.sponge_upper_density = 100 + castro.sponge_lower_pressure = -1 + castro.sponge_upper_pressure = -1 + castro.sponge_lower_factor = 0 + castro.sponge_upper_factor = 1 + castro.sponge_target_x_velocity = 0 + castro.sponge_target_y_velocity = 0 + castro.sponge_target_z_velocity = 0 +[*] castro.sponge_timescale = 1e-07 + castro.state_interp_order = 1 + castro.lin_limit_state_interp = 0 + castro.do_reflux = 1 + castro.update_sources_after_reflux = 1 + castro.allow_non_unit_aspect_zones = 0 + castro.difmag = 0.1 +[*] castro.small_dens = 1e-05 +[*] castro.small_temp = 1e+06 +[*] castro.small_pres = 2.98096e+09 +[*] castro.small_ener = 8.25345e+14 +[*] castro.do_hydro = 1 + castro.time_integration_method = 0 + castro.limit_fourth_order = 1 + castro.initialization_is_cell_average = 0 + castro.use_reconstructed_gamma1 = 0 + castro.add_ext_src = 0 + castro.hybrid_hydro = 0 + castro.ppm_type = 1 + castro.ppm_do_limiting = 1 + castro.mhd_limit_characteristic = 1 + castro.ppm_temp_fix = 0 + castro.plm_iorder = 2 + castro.plm_limiter = 2 + castro.hybrid_riemann = 0 + castro.riemann_solver = 0 + castro.cg_maxiter = 12 + castro.cg_tol = 1e-05 + castro.cg_blend = 2 + castro.use_flattening = 1 + castro.transverse_use_eos = 0 + castro.transverse_reset_density = 1 + castro.transverse_reset_rhoe = 0 + castro.dual_energy_eta1 = 1 + castro.dual_energy_eta2 = 0.0001 + castro.use_pslope = 0 + castro.pslope_cutoff_density = -1e+20 + castro.limit_fluxes_on_small_dens = 0 + castro.speed_limit = 0 +[*] castro.do_sponge = 1 + castro.sponge_implicit = 1 + castro.ext_src_implicit = 0 + castro.source_term_predictor = 0 + castro.first_order_hydro = 0 + castro.xl_ext_bc_type = -1 + castro.xr_ext_bc_type = -1 +[*] castro.yl_ext_bc_type = 1 + castro.yr_ext_bc_type = -1 + castro.zl_ext_bc_type = -1 + castro.zr_ext_bc_type = -1 + castro.hse_zero_vels = 0 +[*] castro.hse_interp_temp = 1 + castro.hse_fixed_temp = -1e+200 +[*] castro.hse_reflect_vels = 1 +[*] castro.fill_ambient_bc = 1 +[*] castro.ambient_fill_dir = 1 +[*] castro.ambient_outflow_vel = 1 + castro.clamp_ambient_temp = 0 + castro.ambient_safety_factor = 1.1 + castro.ambient_density = -1e+200 + castro.ambient_temp = -1e+200 + castro.ambient_energy = -1e+200 + castro.sdc_order = 2 + castro.sdc_quadrature = 0 + castro.sdc_extra = 0 + castro.sdc_solver = 1 + castro.use_axisymmetric_geom_source = 1 + castro.add_sdc_react_source_to_advection = 1 + castro.hydro_memory_footprint_ratio = -1 + castro.fixed_dt = -1 + castro.initial_dt = -1 + castro.dt_cutoff = 1e-12 + castro.max_dt = 1e+200 + castro.cfl = 0.8 +[*] castro.init_shrink = 0.1 + castro.change_max = 1.1 + castro.check_dt_before_advance = 1 + castro.check_dt_after_advance = 1 + castro.plot_per_is_exact = 0 + castro.small_plot_per_is_exact = 0 + castro.use_retry = 1 + castro.retry_subcycle_factor = 0.5 + castro.retry_small_density_cutoff = -1e+200 + castro.abundance_failure_tolerance = 0.01 + castro.abundance_failure_rho_cutoff = -1e+200 + castro.use_post_step_regrid = 0 + castro.max_subcycles = 10 + castro.sdc_iters = 2 + castro.stopping_criterion_field = + castro.stopping_criterion_value = 1e+200 + castro.dtnuc_e = 1e+200 + castro.dtnuc_X = 1e+200 + castro.dtnuc_X_threshold = 0.001 +[*] castro.do_react = 1 +[*] castro.react_T_min = 6e+07 + castro.react_T_max = 1e+200 +[*] castro.react_rho_min = 100 +[*] castro.react_rho_max = 5e+07 + castro.disable_shock_burning = 0 + castro.shock_detection_threshold = 0.666667 + castro.shock_detection_include_sources = 1 + castro.T_guess = 1e+08 + castro.drive_initial_convection = 0 + castro.drive_initial_convection_tmax = 1e+200 + castro.drive_initial_convection_reinit_period = 1e+200 +[*] castro.do_grav = 1 + castro.moving_center = 0 +[*] castro.grav_source_type = 2 +[*] castro.do_rotation = 1 + castro.bndry_func_thread_safe = 1 + castro.grown_factor = 1 + castro.star_at_center = -1 + castro.do_scf_initial_model = 0 + castro.scf_maximum_density = -1e+06 + castro.scf_equatorial_radius = -1e+09 + castro.scf_polar_radius = -1e+09 + castro.scf_relax_tol = 0.001 + castro.scf_max_iterations = 30 + castro.do_special_tagging = 0 + castro.max_tagging_radius = 10 +[*] castro.verbose = 1 + castro.dump_old = 0 + castro.domain_is_plane_parallel = 0 +[*] castro.print_update_diagnostics = 1 +[*] castro.sum_interval = 1 + castro.sum_per = -1 + castro.job_name = Castro + castro.output_at_completion = 1 + castro.reset_checkpoint_time = -1e+200 + castro.reset_checkpoint_step = -1 + castro.store_omegadot = 0 + castro.store_burn_weights = 0 + castro.abort_on_invalid_params = 0 + castro.do_radiation = -1 +[*] gravity.gravity_type = ConstantGrav +[*] gravity.const_grav = -1.5e+14 + gravity.direct_sum_bcs = 0 + gravity.drdxfac = 1 + gravity.lnum = 0 + gravity.verbose = 0 + gravity.no_sync = 0 + gravity.do_composite_phi_correction = 1 + gravity.max_solve_level = 14 + gravity.get_g_from_phi = 0 + gravity.mlmg_max_fmg_iter = 0 + gravity.mlmg_agglomeration = 1 + gravity.mlmg_consolidation = 1 + gravity.mlmg_nsolve = 0 + diffusion.verbose = 0 + diffusion.mlmg_maxorder = 4 +[*] problem.dtemp = 1.2e+09 +[*] problem.x_half_max = 20480 +[*] problem.x_half_width = 2048 +[*] problem.X_min = 0.01 + problem.tag_by_density = 1 +[*] problem.cutoff_density = 25000 + problem.refine_height = 3600 +[*] problem.T_hi = 2e+08 + problem.T_star = 1e+08 +[*] problem.T_lo = 8e+06 +[*] problem.dens_base = 3.43e+06 +[*] problem.H_star = 2000 +[*] problem.atm_delta = 100 + problem.fuel1_name = helium-4 + problem.fuel2_name = + problem.fuel3_name = + problem.fuel4_name = +[*] problem.ash1_name = nickel-56 + problem.ash2_name = + problem.ash3_name = + problem.fuel1_frac = 1 + problem.fuel2_frac = 0 + problem.fuel3_frac = 0 + problem.fuel4_frac = 0 + problem.ash1_frac = 1 + problem.ash2_frac = 0 + problem.ash3_frac = 0 + problem.low_density_cutoff = 0.0001 + problem.smallx = 1e-10 +[*] problem.x_refine_distance = 61440 +[*] problem.max_hse_tagging_level = 3 +[*] problem.max_base_tagging_level = 2 + integrator.X_reject_buffer = 1 + integrator.call_eos_in_rhs = 1 + integrator.integrate_energy = 1 + integrator.jacobian = 1 + integrator.burner_verbose = 0 +[*] integrator.rtol_spec = 1e-06 + integrator.rtol_enuc = 1e-06 +[*] integrator.atol_spec = 1e-06 + integrator.atol_enuc = 1e-06 + integrator.renormalize_abundances = 0 + integrator.SMALL_X_SAFE = 1e-30 + integrator.MAX_TEMP = 1e+11 +[*] integrator.react_boost = 10 + integrator.ode_max_steps = 150000 + integrator.ode_max_dt = 1e+30 + integrator.use_jacobian_caching = 1 + integrator.nonaka_i = 0 + integrator.nonaka_j = 0 + integrator.nonaka_k = 0 + integrator.nonaka_level = 0 + integrator.nonaka_file = nonaka_plot.dat + integrator.use_burn_retry = 0 + integrator.retry_swap_jacobian = 1 + integrator.retry_rtol_spec = 1e-12 + integrator.retry_rtol_enuc = 1e-06 + integrator.retry_atol_spec = 1e-08 + integrator.retry_atol_enuc = 1e-06 + integrator.do_species_clip = 1 + integrator.use_number_densities = 0 + integrator.sdc_burn_tol_factor = 1 + integrator.scale_system = 0 + integrator.nse_iters = 3 + integrator.nse_deriv_dt_factor = 0.05 + integrator.nse_include_enu_weak = 1 + integrator.linalg_do_pivoting = 1 + screening.enable_chabrier1998_quantum_corr = 0 + integrator.subtract_internal_energy = 0 + eos.use_eos_coulomb = 1 + eos.eos_input_is_constant = 1 + eos.eos_ttol = 1e-08 + eos.eos_dtol = 1e-08 + eos.prad_limiter_rho_c = -1 + eos.prad_limiter_delta_rho = -1 + network.small_x = 1e-30 +[*] network.use_tables = 1 + network.use_c12ag_deboer17 = 0 diff --git a/Exec/science/flame_wave/ci-benchmarks/species_diag.out b/Exec/science/flame_wave/ci-benchmarks/species_diag.out index b56e8bb348..6baaafc600 100644 --- a/Exec/science/flame_wave/ci-benchmarks/species_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/species_diag.out @@ -1,5 +1,5 @@ # COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 # TIMESTEP TIME Mass He4 Mass C12 Mass O16 Mass Ne20 Mass Mg24 Mass Si28 Mass S32 Mass Ar36 Mass Ca40 Mass Ti44 Mass Cr48 Mass Fe52 Mass Ni56 - 0 0.0000000000000000 3.5899765861250859e-15 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.1354824912929246e-14 - 1 5.9694928185132600e-08 3.5899765861251656e-15 7.4944874211898590e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.1354897598484802e-14 - 2 1.2535934918877847e-07 3.5899765861253060e-15 7.4944963166333664e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.1354986514892597e-14 + 0 0.0000000000000000e+00 2.8142378112400895e-15 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.1117997526547418e-14 + 1 5.9806047036494849e-08 2.8142323856670296e-15 5.4329394513687282e-21 7.4194476167127106e-24 7.3934515289195991e-24 7.3935541950036108e-24 7.3933935990475538e-24 7.3932321468746301e-24 7.3932312019180565e-24 7.3932308328287436e-24 7.3932307869393029e-24 7.3932307850741563e-24 7.3932307849933102e-24 7.1118070045957091e-14 + 2 1.2559269877663918e-07 2.8142264165380905e-15 1.1401993461614438e-20 7.4929577845086585e-24 7.3944696393124482e-24 7.3939305113509452e-24 7.3935838899709814e-24 7.3932425308155010e-24 7.3932405326302037e-24 7.3932397568503172e-24 7.3932396603620439e-24 7.3932396564405536e-24 7.3932396562705448e-24 7.1118158758588142e-14 diff --git a/Exec/science/flame_wave/initial_model.H b/Exec/science/flame_wave/initial_model.H index 7deae43d94..27d7f2b5ca 100644 --- a/Exec/science/flame_wave/initial_model.H +++ b/Exec/science/flame_wave/initial_model.H @@ -236,6 +236,9 @@ generate_initial_model(const int npts_model, const Real xmin, const Real xmax, for (int i = index_base+1; i < npts_model; i++) { if ((model::profile(model_num).r(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) && !flipped) { + if (i == index_base + 1) { + amrex::Error("atm_delta is too small for the grid spacing (at least one transition zone is required)"); + } isentropic = true; flipped = true; diff --git a/Exec/science/flame_wave/inputs_2d.testsuite b/Exec/science/flame_wave/inputs_2d.testsuite index 347d66e14c..f4f04d5aed 100644 --- a/Exec/science/flame_wave/inputs_2d.testsuite +++ b/Exec/science/flame_wave/inputs_2d.testsuite @@ -102,7 +102,7 @@ problem.T_hi = 2.e8 problem.T_lo = 8.e6 problem.H_star = 2000.e0 -problem.atm_delta = 50.0 +problem.atm_delta = 100.0 problem.fuel1_name = "helium-4" problem.fuel1_frac = 1.0e0 diff --git a/Exec/science/flame_wave/inputs_3d.testsuite.gpu b/Exec/science/flame_wave/inputs_3d.testsuite.gpu index 030bef8aba..2024847e4e 100644 --- a/Exec/science/flame_wave/inputs_3d.testsuite.gpu +++ b/Exec/science/flame_wave/inputs_3d.testsuite.gpu @@ -100,7 +100,7 @@ problem.T_hi = 2.e8 problem.T_lo = 8.e6 problem.H_star = 2000.e0 -problem.atm_delta = 50.0 +problem.atm_delta = 100.0 problem.fuel1_name = "helium-4" problem.fuel1_frac = 1.0e0 diff --git a/Exec/science/massive_star/analysis/massive_star_slice_3d.py b/Exec/science/massive_star/analysis/massive_star_slice_3d.py new file mode 100755 index 0000000000..4de170ec4a --- /dev/null +++ b/Exec/science/massive_star/analysis/massive_star_slice_3d.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python3 + +import os +import sys + +import matplotlib +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import ImageGrid + +import yt +from yt.frontends.boxlib.api import CastroDataset + +matplotlib.use('agg') + + +plotfile = sys.argv[1] +ds = CastroDataset(plotfile) + +t_drive = 0.0 +if "[*] castro.drive_initial_convection_tmax" in ds.parameters: + t_drive = ds.parameters["[*] castro.drive_initial_convection_tmax"] +elif "castro.drive_initial_convection_tmax" in ds.parameters: + t_drive = ds.parameters["castro.drive_initial_convection_tmax"] +print(t_drive) + +domain_frac = 0.2 + +xmin = ds.domain_left_edge[0] +xmax = ds.domain_right_edge[0] +xctr = 0.5 * (xmin + xmax) +L_x = xmax - xmin +xmin = xctr - 0.5 * domain_frac * L_x +xmax = xctr + 0.5 * domain_frac * L_x +L_x = xmax - xmin + +ymin = ds.domain_left_edge[1] +ymax = ds.domain_right_edge[1] +yctr = 0.5 * (ymin + ymax) +L_y = ymax - ymin +ymin = yctr - 0.5 * domain_frac * L_y +ymax = yctr + 0.5 * domain_frac * L_y +L_y = ymax - ymin + +zmin = ds.domain_left_edge[2] +zmax = ds.domain_right_edge[2] +zctr = 0.5 * (zmin + zmax) +L_z = zmax - zmin +zmin = zctr - 0.5 * domain_frac * L_z +zmax = zctr + 0.5 * domain_frac * L_z +L_z = zmax - zmin + +fig = plt.figure() + + +fields = ["MachNumber", "magvort", "abar", "enuc"] + +grid = ImageGrid(fig, 111, nrows_ncols=(2, (len(fields)+1)//2), + axes_pad=0.75, cbar_pad=0.05, label_mode="L", cbar_mode="each") + + +for i, f in enumerate(fields): + + sp = yt.SlicePlot(ds, "y", f, + center=[xctr, yctr, zctr], width=[L_x, L_z], + fontsize="12") + + sp.set_buff_size((2400,2400)) + sp.swap_axes() + + if f == "Ye": + sp.set_zlim(f, 0.46, 0.5) + sp.set_log(f, False) + sp.set_cmap(f, "magma_r") + elif f == "abar": + sp.set_log(f, False) + sp.set_cmap(f, "viridis") + elif f == "enuc": + sp.set_log(f, True, linthresh=1.e12) + sp.set_zlim(f, -1.e20, 1.e20) + sp.set_cmap(f, "bwr") + elif f == "MachNumber": + sp.set_zlim(f, 1.e-4, 0.3) + sp.set_cmap(f, "plasma") + elif f == "magvel": + sp.set_zlim(f, 100.0, 2.e7) + sp.set_cmap(f, "viridis") + elif f == "magvort": + sp.set_cmap(f, "magma") + sp.set_zlim(f, 1.e-2, 5) + + sp.set_axes_unit("km") + + plot = sp.plots[f] + plot.figure = fig + plot.axes = grid[i].axes + plot.cax = grid.cbar_axes[i] + if i < len(fields)-1: + grid[i].axes.xaxis.offsetText.set_visible(False) + + sp._setup_plots() + +fig.text(0.02, 0.02, f"$t - \\tau_\\mathrm{{drive}}$ = {float(ds.current_time) - t_drive:6.1f} s", + transform=fig.transFigure) + +fig.set_size_inches(11, 10) +plt.tight_layout() +plt.savefig(f"{os.path.basename(plotfile)}_slice.png") diff --git a/Exec/science/massive_star/inputs_3d.nse b/Exec/science/massive_star/inputs_3d.nse index d7e2cdd48e..08f416b470 100644 --- a/Exec/science/massive_star/inputs_3d.nse +++ b/Exec/science/massive_star/inputs_3d.nse @@ -3,10 +3,10 @@ amr.plot_files_output = 1 amr.checkpoint_files_output = 1 -max_step = 50 +max_step = 5000000 stop_time = 360000 -geometry.is_periodic = 0 0 +geometry.is_periodic = 0 0 0 geometry.coord_sys = 0 # r-z coordinates geometry.prob_lo = 0. 0. 0. @@ -14,7 +14,7 @@ geometry.prob_hi = 1.6384e10 1.6384e10 1.6384e10 amr.n_cell = 512 512 512 -amr.max_level = 3 # maximum level number allowed +amr.max_level = 2 # maximum level number allowed castro.lo_bc = 2 2 2 castro.hi_bc = 2 2 2 @@ -50,21 +50,21 @@ castro.sponge_upper_density = 1.e3 castro.sponge_lower_density = 1.e2 castro.sponge_timescale = 1.e-3 -castro.cfl = 0.4 # cfl number for hyperbolic system -castro.init_shrink = 0.01 # scale back initial timestep by this factor +castro.cfl = 0.5 # cfl number for hyperbolic system +castro.init_shrink = 0.1 # scale back initial timestep by this factor castro.change_max = 1.2 # factor by which dt is allowed to change each timestep castro.sum_interval = 10 # timesteps between computing and printing volume averages #castro.dtnuc_e = 0.25 #castro.dtnuc_X = 0.25 -amr.ref_ratio = 4 2 2 2 # refinement ratio +amr.ref_ratio = 4 4 2 2 2 # refinement ratio amr.regrid_int = 10000 # how often to regrid -amr.n_error_buf = 4 2 2 2 # number of buffer cells in error est +amr.n_error_buf = 4 4 2 2 2 # number of buffer cells in error est amr.grid_eff = 0.7 # what constitutes an efficient grid amr.check_file = massive_star_chk # root name of checkpoint file -amr.check_int = 200 # number of timesteps between checkpoints +amr.check_int = 10 # number of timesteps between checkpoints amr.plot_file = massive_star_plt # root name of plot file amr.plot_per = 5.0 @@ -76,10 +76,6 @@ amr.small_plot_per = 0.5 amr.small_plot_vars = density Temp in_nse amr.derive_small_plot_vars = abar Ye enuc MachNumber magvel magvort -#amr.checkpoint_files_output = 0 -#amr.plot_files_output = 0 -#amr.smallplot_files_output = 0 - fab.format = NATIVE_32 castro.plot_per_is_exact = 0 @@ -88,6 +84,8 @@ castro.plot_per_is_exact = 0 amr.max_grid_size = 64 # maximum grid size allowed -- used to control parallelism amr.blocking_factor = 32 # block factor in grid generation +amr.subcycling_mode = None + amr.v = 1 # control verbosity in Amr.cpp castro.v = 1 # control verbosity in Castro.cpp @@ -116,17 +114,13 @@ castro.drive_initial_convection_tmax = 50 # refinement -amr.refinement_indicators = denerr denerr2 denerr3 +amr.refinement_indicators = denerr denerr3 amr.refine.denerr.max_level = 1 amr.refine.denerr.value_greater = 2.e3 amr.refine.denerr.field_name = density -amr.refine.denerr2.max_level = 2 -amr.refine.denerr2.value_greater = 3.e4 -amr.refine.denerr2.field_name = density - -amr.refine.denerr3.max_level = 3 +amr.refine.denerr3.max_level = 2 amr.refine.denerr3.value_greater = 3.e5 amr.refine.denerr3.field_name = density @@ -146,7 +140,15 @@ network.Si_nse = 0.02 network.C_nse = 1.0 network.O_nse = 1.0 -integrator.ode_max_steps = 500000 +integrator.ode_max_steps = 5000 + +integrator.use_burn_retry = 1 +integrator.retry_swap_jacobian = 1 + +integrator.retry_rtol_spec = 1.e-5 +integrator.retry_atol_spec = 1.e-5 +integrator.retry_rtol_enuc = 1.e-5 +integrator.retry_atol_enuc = 1.e-5 network.small_x = 1.e-10 diff --git a/Exec/science/wdmerger/_prob_params b/Exec/science/wdmerger/_prob_params index f4c0317d85..f1b0c49239 100644 --- a/Exec/science/wdmerger/_prob_params +++ b/Exec/science/wdmerger/_prob_params @@ -98,7 +98,7 @@ onemg_wd_mg_frac real 0.05e0_rt y # Parameters for interpolation from 1D model to 3D model # Default to interpolation that preserves temperature; otherwise, use pressure -interp_temp logical .true. y +interp_temp bool true y # Number of sub-grid-scale zones to use nsub integer 1 y diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 6bdddfb49a..4f6e0bf29a 100644 --- a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out +++ b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out @@ -1,29 +1,29 @@ - plotfile = plt00095 + plotfile = plt00086 time = 1.25 variables minimum value maximum value - density 8.7136170328e-05 13347626.831 - xmom -4.4623942798e+14 1.4982725823e+15 - ymom -1.8879302767e+15 1.980304858e+15 + density 8.693611703e-05 19565534.299 + xmom -5.4964100651e+14 1.3559128302e+14 + ymom -2.5530096328e+15 2.5530122744e+15 zmom 0 0 - rho_E 7.7948360809e+11 5.64675518e+24 - rho_e 7.4322366401e+11 5.6244082614e+24 - Temp 100000 3969215375.1 - rho_He4 8.7136170328e-17 1473.6068478 - rho_C12 3.4854468131e-05 5028784.2519 - rho_O16 5.2281702196e-05 7774667.7686 - rho_Ne20 8.7136170328e-17 448739.14288 - rho_Mg24 8.7136170328e-17 1131012.4412 - rho_Si28 8.7136170328e-17 4244213.8084 - rho_S32 8.7136170328e-17 2178795.6667 - rho_Ar36 8.7136170328e-17 496791.83986 - rho_Ca40 8.7136170328e-17 381536.89513 - rho_Ti44 8.7136170328e-17 1576.1652647 - rho_Cr48 8.7136170328e-17 1467.4931808 - rho_Fe52 8.7136170328e-17 14841.830382 - rho_Ni56 8.7136170328e-17 182905.91385 - phiGrav -4.6146743474e+17 -2.2055817929e+16 - grav_x -461199998.92 -48603.52543 - grav_y -444710966.19 392312565.37 + rho_E 7.4982062146e+11 5.0669247218e+24 + rho_e 7.1077581849e+11 5.0640768325e+24 + Temp 242288.68588 1409652233.6 + rho_He4 8.693611703e-17 3.5999033031 + rho_C12 3.4774446812e-05 7825956.6934 + rho_O16 5.2161670217e-05 11739149.75 + rho_Ne20 8.693611703e-17 181951.05725 + rho_Mg24 8.693611703e-17 1192.7969771 + rho_Si28 8.693611703e-17 6.6913703161 + rho_S32 8.693611703e-17 0.00019493291757 + rho_Ar36 8.693611703e-17 1.9565534609e-05 + rho_Ca40 8.693611703e-17 1.9565534331e-05 + rho_Ti44 8.693611703e-17 1.9565534308e-05 + rho_Cr48 8.693611703e-17 1.9565534308e-05 + rho_Fe52 8.693611703e-17 1.9565534308e-05 + rho_Ni56 8.693611703e-17 1.9565534308e-05 + phiGrav -5.8708033902e+17 -2.3375498341e+16 + grav_x -684991644 -51428.243166 + grav_y -739606241.84 739606820.44 grav_z 0 0 - rho_enuc -7.4954583352e+21 6.7150736332e+26 + rho_enuc -8.1740856019e+12 7.6429034917e+23 diff --git a/Exec/science/wdmerger/problem_initialize_state_data.H b/Exec/science/wdmerger/problem_initialize_state_data.H index 1b37fde218..3b19019039 100644 --- a/Exec/science/wdmerger/problem_initialize_state_data.H +++ b/Exec/science/wdmerger/problem_initialize_state_data.H @@ -39,7 +39,10 @@ void problem_initialize_state_data (int i, int j, int k, // things right on the coarse levels. So we can still use the // interpolation scheme, because it handles this special case // for us by simply using the central zone of the model; we - // just need to make sure we catch it. + // just need to make sure we catch it. Finally, sometimes + // the stars are so close that the interpolation will overlap. + // in this case, look at which model has the highest density + // at that location and use that model. const Real* dx = geomdata.CellSize(); @@ -53,43 +56,59 @@ void problem_initialize_state_data (int i, int j, int k, eos_t zone_state; - if (problem::mass_P > 0.0_rt && (dist_P < problem::radius_P || - (problem::radius_P <= max_dx && dist_P < max_dx))) { - Real pos[3] = {loc[0] - problem::center_P_initial[0], - loc[1] - problem::center_P_initial[1], - loc[2] - problem::center_P_initial[2]}; + bool P_star_test = problem::mass_P > 0.0_rt && + (dist_P < problem::radius_P || + (problem::radius_P <= max_dx && dist_P < max_dx)); - zone_state.rho = interpolate_3d(pos, dx, model::idens, problem::nsub, 0); - zone_state.T = interpolate_3d(pos, dx, model::itemp, problem::nsub, 0); - for (int n = 0; n < NumSpec; ++n) { - zone_state.xn[n] = interpolate_3d(pos, dx, model::ispec + n, problem::nsub, 0); + bool S_star_test = problem::mass_S > 0.0_rt && + (dist_S < problem::radius_S || + (problem::radius_S <= max_dx && dist_S < max_dx)); + + double rho_P{0.0}; + double rho_S{0.0}; + + if (P_star_test || S_star_test) { + + Real pos_P[3] = {loc[0] - problem::center_P_initial[0], + loc[1] - problem::center_P_initial[1], + loc[2] - problem::center_P_initial[2]}; + + if (problem::mass_P > 0.0_rt) { + rho_P = interpolate_3d(pos_P, dx, model::idens, problem::nsub, 0); } -#ifdef AUX_THERMO - set_aux_comp_from_X(zone_state); -#endif - eos(eos_input_rt, zone_state); + Real pos_S[3] = {loc[0] - problem::center_S_initial[0], + loc[1] - problem::center_S_initial[1], + loc[2] - problem::center_S_initial[2]}; - } - else if (problem::mass_S > 0.0_rt && (dist_S < problem::radius_S || - (problem::radius_S <= max_dx && dist_S < max_dx))) { - Real pos[3] = {loc[0] - problem::center_S_initial[0], - loc[1] - problem::center_S_initial[1], - loc[2] - problem::center_S_initial[2]}; - - zone_state.rho = interpolate_3d(pos, dx, model::idens, problem::nsub, 1); - zone_state.T = interpolate_3d(pos, dx, model::itemp, problem::nsub, 1); - for (int n = 0; n < NumSpec; ++n) { - zone_state.xn[n] = interpolate_3d(pos, dx, model::ispec + n, problem::nsub, 1); + if (problem::mass_S > 0.0_rt) { + rho_S = interpolate_3d(pos_S, dx, model::idens, problem::nsub, 1); } + + if (rho_P > rho_S) { + // use the primary star initialization + zone_state.rho = rho_P; + zone_state.T = interpolate_3d(pos_P, dx, model::itemp, problem::nsub, 0); + for (int n = 0; n < NumSpec; ++n) { + zone_state.xn[n] = interpolate_3d(pos_P, dx, model::ispec + n, problem::nsub, 0); + } + + } else { + // use the secondary star initialization + zone_state.rho = rho_S; + zone_state.T = interpolate_3d(pos_S, dx, model::itemp, problem::nsub, 1); + for (int n = 0; n < NumSpec; ++n) { + zone_state.xn[n] = interpolate_3d(pos_S, dx, model::ispec + n, problem::nsub, 1); + } + } + #ifdef AUX_THERMO set_aux_comp_from_X(zone_state); #endif eos(eos_input_rt, zone_state); - } - else { + } else { zone_state.rho = ambient::ambient_state[URHO]; zone_state.T = ambient::ambient_state[UTEMP]; @@ -131,17 +150,17 @@ void problem_initialize_state_data (int i, int j, int k, if (problem::problem != 1) { - if (dist_P < problem::radius_P) { - state(i,j,k,UMX) += problem::vel_P[0] * state(i,j,k,URHO); - state(i,j,k,UMY) += problem::vel_P[1] * state(i,j,k,URHO); - state(i,j,k,UMZ) += problem::vel_P[2] * state(i,j,k,URHO); - } - else if (dist_S < problem::radius_S) { - state(i,j,k,UMX) += problem::vel_S[0] * state(i,j,k,URHO); - state(i,j,k,UMY) += problem::vel_S[1] * state(i,j,k,URHO); - state(i,j,k,UMZ) += problem::vel_S[2] * state(i,j,k,URHO); + if (P_star_test || S_star_test) { + if (rho_P > rho_S && problem::mass_P > 0.0_rt && dist_P < problem::radius_P) { + state(i,j,k,UMX) += problem::vel_P[0] * state(i,j,k,URHO); + state(i,j,k,UMY) += problem::vel_P[1] * state(i,j,k,URHO); + state(i,j,k,UMZ) += problem::vel_P[2] * state(i,j,k,URHO); + } else if (problem::mass_S > 0.0_rt && dist_S < problem::radius_S) { + state(i,j,k,UMX) += problem::vel_S[0] * state(i,j,k,URHO); + state(i,j,k,UMY) += problem::vel_S[1] * state(i,j,k,URHO); + state(i,j,k,UMZ) += problem::vel_S[2] * state(i,j,k,URHO); + } } - } // If we're in the inertial reference frame, use rigid body rotation with velocity omega x r. diff --git a/Exec/science/wdmerger/tests/inputs_scaling b/Exec/science/wdmerger/tests/inputs_scaling new file mode 100644 index 0000000000..cd48f143b3 --- /dev/null +++ b/Exec/science/wdmerger/tests/inputs_scaling @@ -0,0 +1,459 @@ +## Latest inputs file being used to reproduce initial conditions from Pakmor et al. 2022 +## with 100 km resolution using for scaling + +############################## CASTRO INPUTS ############################################### + +############################################################################################ +# Geometry +############################################################################################ + +# Non-periodic boundary conditions +geometry.is_periodic = 0 0 0 + +#if AMREX_SPACEDIM == 3 + + # Cartesian coordinate system + geometry.coord_sys = 0 + + # Lower boundary limits in physical space + geometry.prob_lo = -5.12e9 -5.12e9 -5.12e9 + +#elif AMREX_SPACEDIM == 2 + + # Cartesian coordinate system + geometry.coord_sys = 1 + + # Lower boundary limits in physical space + geometry.prob_lo = 0.0e0 -5.12e9 + +#endif + +# Upper boundary limits in physical space +geometry.prob_hi = 5.12e9 5.12e9 5.12e9 + +############################################################################################ +# Boundary conditions +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +############################################################################################ + +#if AMREX_SPACEDIM == 3 + + # Boundary conditions on lo x, y, and z edges + castro.lo_bc = 2 2 2 + +#elif AMREX_SPACEDIM == 2 + + # Boundary conditions on lo x, y, and z edges + castro.lo_bc = 3 2 + +#endif + +# Boundary conditions on hi x, y, and z edges +castro.hi_bc = 2 2 2 + +############################################################################################ +# Timestepping +############################################################################################ + +# Maximum number of level 0 steps +max_step = 10000000 + +# Simulation end time +stop_time = 0.1 + +# CFL number for hyperbolic system +castro.cfl = 0.5 + +# Fixed level 0 timestep; unused if < 0 +castro.fixed_dt = -1.0 + +# Scale back initial timestep by this factor +castro.init_shrink = 0.01 + +# Factor by which dt is allowed to change each timestep +castro.change_max = 1.3 + +# If we regrid on Level 0, compute a new timestep afterward +amr.compute_new_dt_on_regrid = 1 + +# Use a retry if an advance violated our stability criteria +castro.use_retry = 1 + +# Skip retries for small density if the starting density was less than this threshold +castro.retry_small_density_cutoff = 1.0e0 + +# Don't abort for invalid X if the zone density is less than this threshold +castro.abundance_failure_rho_cutoff = 1.0e0 + +# Maximum number of subcycles +# Default is 10, 16 is recommended value +castro.max_subcycles = 16 + + +############################################################################################ +# Resolution, gridding and AMR +############################################################################################ + +#if AMREX_SPACEDIM == 3 + + # Number of cells on the coarse grid + amr.n_cell = 256 256 256 + +#elif AMREX_SPACEDIM == 2 + + # Number of cells on the coarse grid + amr.n_cell = 128 256 + +#endif + +# Maximum level number allowed +amr.max_level = 1 + +# Refinement ratio +amr.ref_ratio = 4 + +# How many coarse timesteps between regridding +amr.regrid_int = 2 + +# Allow special regrids based on stability criteria +castro.use_post_step_regrid = 0 + +# Number of buffer cells in error estimation +amr.n_error_buf = 2 2 2 2 2 2 2 2 2 2 + +# Maximum grid size at each level +amr.max_grid_size = 64 64 64 64 64 64 64 64 + +# Grid sizes must be a multiple of blocking factor +amr.blocking_factor = 32 + +# What constitutes an efficient grid +amr.grid_eff = 0.9 + +# Order of reconstruction for interpolation +castro.state_interp_order = 0 + +# Limiting on state data interpolation (preserve linear combinations) +castro.lin_limit_state_interp = 1 + +# Add refinement indicators +amr.refinement_indicators = density density2 temperature + +# Density refinement criterion +amr.refine.density.value_greater = 1.0e0 +amr.refine.density.field_name = density +amr.refine.density.max_level = 1 + +# Density2 refinement criterion +amr.refine.density2.value_greater = 1.0e5 +amr.refine.density2.field_name = density +amr.refine.density2.max_level = 20 + +# Temperature refinement criterion +amr.refine.temperature.value_greater = 5.0e8 +amr.refine.temperature.field_name = Temp +amr.refine.temperature.max_level = 3 + +# Avoid tagging near the domain boundary +castro.max_tagging_radius = 0.75e0 + +# Whether or not to use AMR subcycling +amr.subcycling_mode = None + +############################################################################################ +# Physics to include +############################################################################################ + +# Whether or not to do hydrodynamics +castro.do_hydro = 1 + +# Whether or not to do gravity +castro.do_grav = 1 + +# Whether or not to do reactions +castro.do_react = 1 + +# Whether or not to apply the sponge +castro.do_sponge = 1 + +# Whether or not to apply external source terms +castro.add_ext_src = 1 +castro.ext_src_implicit = 1 + +# Whether or not to include the rotation source term +castro.do_rotation = 1 + +############################################################################################ +# PPM/Hydro options +############################################################################################ + +# Piecewise parabolic with the original limiters (0 is piecewise linear; 2 is new limiters) +castro.ppm_type = 1 + +# Use the EOS in calculation of the edge states going into the Riemann solver +castro.ppm_temp_fix = 0 + +# Which Riemann solver to use. +# 0 = Colella, Glaz, and Ferguson (cheaper, less accurate) +# 1 = Colella and Glaz 1985 (more expensive, more accurate) +# 2 = HLL +castro.riemann_solver = 0 + +# For the CG Riemann solver, we need to tell the solver not to quit when +# the iterations don't converge, but instead to do additional bisection iteration. +castro.cg_blend = 2 + +# Use a lagged predictor estimate of the source terms in the hydro +castro.source_term_predictor = 1 + +# Whether to use the hybrid advection technique that conserves angular momentum +castro.hybrid_hydro = 0 + +# Reset (rho*e) if it goes negative in the transverse terms +castro.transverse_reset_rhoe = 1 + +# Reset rho if it goes negative in the transverse terms +castro.transverse_reset_density = 1 + +# Explicitly limit fluxes to avoid hitting a negative density +castro.limit_fluxes_on_small_dens = 1 + +# Set global simulation speed limit +castro.speed_limit = 1.498962290e9 # 0.05c + +############################################################################################ +# Thermodynamics +############################################################################################ + +# Minimum allowable temperature (K) +castro.small_temp = 1.e5 + +# Minimum allowable density (g / cm**3) +castro.small_dens = 1.e-5 + +# Threshold for when to use the internal energy in calculating pressure +castro.dual_energy_eta1 = 1.0e-3 + +# Threshold for when to use (E - K) in updating internal energy +castro.dual_energy_eta2 = 1.0e-4 + +# Use Coulomb corrections in Helmholtz EOS +eos.use_eos_coulomb = 1 + +# Keep EOS inputs constant after EOS evaluation +eos.eos_input_is_constant = 1 + +# Ambient temperature (K) +castro.ambient_temp = 1.0e7 + +# Ambient density (g / cm**3) +castro.ambient_density = 1.0e-4 + +# Clamp temperature in ambient zones to its initial value +castro.clamp_ambient_temp = 1 + +############################################################################################ +# Reactions/Network +############################################################################################ + +# Limit timestep based on nuclear burning considerations (changes in internal energy) +castro.dtnuc_e = 1.e200 + +#Limit timestep based on nuclear burning considerations (changes in species) +castro.dtnuc_X = 1.e200 + +# Minimum temperature for allowing nuclear burning +castro.react_T_min = 1.0e8 + +# Maximum temperature for allowing nuclear burning +castro.react_T_max = 1.0e12 + +# Minimum density for allowing nuclear burning +castro.react_rho_min = 1.0e6 + +# Maximum density for allowing nuclear burning +castro.react_rho_max = 1.0e12 + +# Smallest allowable mass fraction +network.small_x = 1.0e-12 + +# Evaluate the RHS during the burn +integrator.call_eos_in_rhs = 1 + +# Integration tolerances +integrator.rtol_spec = 1.0e-6 +integrator.atol_spec = 1.0e-6 + +integrator.rtol_enuc = 1.0e-6 +integrator.atol_enuc = 1.0e-6 + +# Do not abort or retry on a failed burn (Castro will handle this) +integrator.abort_on_failure = 0 + +# Renormalize abundances during the burn +integrator.renormalize_abundances = 1 + +# Maximum temperature allowed in the burn +integrator.MAX_TEMP = 1.0e10 + +# Use tabular rate evaluation when available +network.use_tables = 1 + +############################################################################################ +# Gravity +############################################################################################ + +# Full self-gravity with the Poisson equation +gravity.gravity_type = PoissonGrav + +# Multipole expansion includes terms up to r**(-max_multipole_order) +gravity.max_multipole_order = 6 + +# Tolerance for multigrid solver for phi solves +gravity.abs_tol = 1.e-10 + +# Use sync solve for gravity after refluxing +gravity.no_sync = 0 + +# Disable the use of the lagged composite correction for the potential +gravity.do_composite_phi_correction = 0 + +############################################################################################ +# Rotation +############################################################################################ + +# Rotational period of the rotating reference frame +castro.rotational_period = 100.0 + +############################################################################################ +# Sponge +############################################################################################ + +castro.sponge_lower_density = 1.0e0 +castro.sponge_upper_density = 1.0e0 +castro.sponge_timescale = 0.01e0 + +############################################################################################ +# Load balancing +############################################################################################ + +# Choice of load balancing strategy to use +DistributionMapping.strategy = KNAPSACK + +# Efficiency demanded from the knapsack algorithm +DistributionMapping.efficiency = 0.9 + +############################################################################################ +# Diagnostics and I/O +############################################################################################ + +# Timesteps between computing and printing volume averaged diagnostic quantities +castro.sum_interval = 0 + +# Simulation time between computing and printing volume averaged diagnostic quantities +castro.sum_per = -1.0 + +# Gravitational wave strain observation distance +castro.gw_dist = 10.0 + +# Name the job +castro.job_name = wdmerger + +# Whether or not to output plotfiles +amr.plot_files_output = 0 + +# Whether or not to output checkpoints +amr.checkpoint_files_output = 0 + +# Root name of checkpoint files +amr.check_file = chk + +# We want to store the 'old' state data in checkpoints +castro.dump_old = 1 + +# Simulation time between checkpoints +amr.check_per = 0.1 + +# Number of timesteps between checkpoints +amr.check_int = 20 + +# Root name of plot files +amr.plot_file = plt + +# Simulation time between plotfiles +amr.plot_per = -1 + +# Number of timesteps between plotfiles +#amr.plot_int = -1 + +# Root name of small plot files +amr.small_plot_file = smallplt + +# Simulation time between small plotfiles +amr.small_plot_per = -1 + +# Number of timesteps between small plotfiles +amr.small_plot_int = -1 + +# Do not write plotfiles when we dump checkpoints +amr.write_plotfile_with_checkpoint = 0 + +# Write final checkpoint/plotfile +castro.output_at_completion = 1 + +# Do not write a plotfile or checkpoint on restart +amr.plotfile_on_restart = 1 +amr.checkpoint_on_restart = 1 + +# Control verbosity in Amr.cpp +amr.v = 1 + +# Control verbosity in Castro.cpp +castro.v = 1 + +# Control verbosity in Gravity.cpp +gravity.v = 1 + +# State variables to add to plot files +amr.plot_vars = NONE + +# Derived variables to add to plot files +amr.derive_plot_vars = NONE + +# State variables to add to small plot files +amr.small_plot_vars = NONE + +# Derived variables to add to small plot files +amr.derive_small_plot_vars = NONE + +# Name of the diagnostic sum output files +#amr.data_log = NONE + +############################################################################################ +# Problem parameters +############################################################################################ + +problem.mass_P = 1.05 +problem.mass_S = 0.70 + +problem.co_wd_c_frac = 0.5 +problem.co_wd_o_frac = 0.5 +problem.co_wd_he_shell_mass = 0.03 + +problem.max_co_wd_mass = 1.50 + +problem.nsub = 4 + +problem.problem = 1 + +problem.roche_radius_factor = 1.0e0 + +problem.interp_temp = 1 + +problem.relaxation_damping_factor = -1.0e-1 +problem.relaxation_density_cutoff = 1.0e3 +problem.relaxation_cutoff_time = -1.0e0 + +problem.stellar_temp = 5.0e5 diff --git a/Exec/science/wdmerger/wdmerger_util.H b/Exec/science/wdmerger/wdmerger_util.H index 1b61f934d3..d2c28b3bfa 100644 --- a/Exec/science/wdmerger/wdmerger_util.H +++ b/Exec/science/wdmerger/wdmerger_util.H @@ -16,7 +16,7 @@ void kepler_third_law (Real radius_1, Real mass_1, Real radius_2, Real mass_2, Real& period, Real eccentricity, Real phi, Real& a, Real& r_1, Real& r_2, Real& v_1r, Real& v_2r, Real& v_1p, Real& v_2p); -void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec], Real envelope_comp[NumSpec]); +void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec], Real envelope_comp[NumSpec], const std::string& star_type); void ensure_primary_mass_larger (); diff --git a/Exec/science/wdmerger/wdmerger_util.cpp b/Exec/science/wdmerger/wdmerger_util.cpp index d2a5b46193..65b09f89e5 100644 --- a/Exec/science/wdmerger/wdmerger_util.cpp +++ b/Exec/science/wdmerger/wdmerger_util.cpp @@ -73,7 +73,7 @@ void kepler_third_law (Real radius_1, Real mass_1, Real radius_2, Real mass_2, // Given a WD mass, set its core and envelope composition. -void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec], Real envelope_comp[NumSpec]) +void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec], Real envelope_comp[NumSpec], const std::string& star_type) { int iHe4 = network_spec_index("helium-4"); int iC12 = network_spec_index("carbon-12"); @@ -98,6 +98,8 @@ void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec] core_comp[iHe4] = 1.0_rt; + amrex::Print() << "Created a pure He " << star_type << "." << std::endl; + for (int n = 0; n < NumSpec; ++n) { envelope_comp[n] = core_comp[n]; } @@ -117,6 +119,10 @@ void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec] envelope_mass = problem::hybrid_wd_he_shell_mass; + amrex::Print()<< "Creating " << star_type << "with CO core with mass fractions C = "<< problem::hybrid_wd_c_frac << "and O = " + << problem::hybrid_wd_o_frac <<" and a He shell of solar mass =" << problem::hybrid_wd_he_shell_mass << "." << std::endl; + + if (envelope_mass > 0.0_rt) { if (iHe4 < 0) { amrex::Error("Must have He4 in the nuclear network."); @@ -128,7 +134,6 @@ void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec] envelope_comp[n] = core_comp[n]; } } - } else if (mass >= problem::max_hybrid_wd_mass && mass < problem::max_co_wd_mass) { @@ -144,6 +149,9 @@ void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec] envelope_mass = problem::co_wd_he_shell_mass; + amrex::Print()<<"Creating " << star_type << " with CO core with mass fractions C = "< 0.0_rt) { if (iHe4 < 0) { amrex::Error("Must have He4 in the nuclear network."); @@ -173,6 +181,8 @@ void set_wd_composition (Real mass, Real& envelope_mass, Real core_comp[NumSpec] core_comp[iNe20] = problem::onemg_wd_ne_frac; core_comp[iMg24] = problem::onemg_wd_mg_frac; + amrex::Print()<<"Creating an ONeMg " << star_type << "." < #include diff --git a/Source/diffusion/Make.package b/Source/diffusion/Make.package index d39ddec718..be1fe61fbe 100644 --- a/Source/diffusion/Make.package +++ b/Source/diffusion/Make.package @@ -8,5 +8,4 @@ CEXE_sources += diffusion_util.cpp CEXE_sources += Castro_diffusion.cpp CEXE_sources += Diffusion.cpp -CEXE_sources += diffusion_params.cpp CEXE_headers += Diffusion.H diff --git a/Source/diffusion/diffusion_params.cpp b/Source/diffusion/diffusion_params.cpp deleted file mode 100644 index 2518c14740..0000000000 --- a/Source/diffusion/diffusion_params.cpp +++ /dev/null @@ -1,7 +0,0 @@ -#include - -#include - -#include - -#include diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index cbd002f889..574456bc53 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -112,10 +112,10 @@ enum int_method { CornerTransportUpwind = 0, // why an advance failed. struct advance_status { - bool success; - Real suggested_dt; - std::string reason; - advance_status() : success(true), suggested_dt(0.0_rt), reason("") {} + bool success{true}; + Real suggested_dt{}; + std::string reason{}; + advance_status() = default; }; /// diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 702385e75c..4ceac43446 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -3551,12 +3551,14 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time int boundary_buf = n_error_buf[dim] + blocking_factor[dim] / ref_ratio[dim]; - if ((physbc_lo[dim] != amrex::PhysBCType::symmetry && physbc_lo[dim] != amrex::PhysBCType::interior) && + if ((physbc_lo[dim] != amrex::PhysBCType::symmetry && + physbc_lo[dim] != amrex::PhysBCType::interior) && (idx[dim] <= domlo[dim] + boundary_buf)) { outer_boundary_test[dim] = true; } - if ((physbc_hi[dim] != amrex::PhysBCType::symmetry && physbc_lo[dim] != amrex::PhysBCType::interior) && + if ((physbc_hi[dim] != amrex::PhysBCType::symmetry && + physbc_hi[dim] != amrex::PhysBCType::interior) && (idx[dim] >= domhi[dim] - boundary_buf)) { outer_boundary_test[dim] = true; } diff --git a/Source/driver/Castro_advance_ctu.cpp b/Source/driver/Castro_advance_ctu.cpp index 911a2780e9..af387a92d2 100644 --- a/Source/driver/Castro_advance_ctu.cpp +++ b/Source/driver/Castro_advance_ctu.cpp @@ -12,13 +12,16 @@ using namespace amrex; advance_status -Castro::do_advance_ctu (Real time, Real dt) +Castro::do_advance_ctu (Real time, Real dt) // NOLINT(readability-convert-member-functions-to-static) { // this routine will advance the old state data (called Sborder here) // to the new time, for a single level. The new data is called // S_new here. The update includes reactions (if we are not doing // SDC), hydro, and the source terms. + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + BL_PROFILE("Castro::do_advance_ctu()"); advance_status status {}; @@ -154,11 +157,18 @@ Castro::retry_advance_ctu(Real dt, const advance_status& status) if (do_retry) { - if (status.suggested_dt > 0.0_rt && status.suggested_dt < dt) { - dt_subcycle = status.suggested_dt; + int max_level_to_advance = level; + + if (parent->subcyclingMode() == "None" && level == 0) { + max_level_to_advance = parent->finestLevel(); } - else { - dt_subcycle = std::min(dt, dt_subcycle) * retry_subcycle_factor; + + for (int lev = level; lev <= max_level_to_advance; ++lev) { + if (status.suggested_dt > 0.0_rt && status.suggested_dt < dt) { + getLevel(lev).dt_subcycle = status.suggested_dt; + } else { + getLevel(lev).dt_subcycle = std::min(dt, getLevel(lev).dt_subcycle) * retry_subcycle_factor; + } } if (verbose && ParallelDescriptor::IOProcessor()) { @@ -175,51 +185,56 @@ Castro::retry_advance_ctu(Real dt, const advance_status& status) // be useful to us at the end of the timestep when we need // to restore the original old data. - save_data_for_retry(); + for (int lev = level; lev <= max_level_to_advance; ++lev) { + getLevel(lev).save_data_for_retry(); - // Clear the contribution to the fluxes from this step. + // Clear the contribution to the fluxes from this step. - for (int dir = 0; dir < 3; ++dir) { - fluxes[dir]->setVal(0.0); - } + for (int dir = 0; dir < 3; ++dir) { + getLevel(lev).fluxes[dir]->setVal(0.0); + } - for (int dir = 0; dir < 3; ++dir) { - mass_fluxes[dir]->setVal(0.0); - } + for (int dir = 0; dir < 3; ++dir) { + getLevel(lev).mass_fluxes[dir]->setVal(0.0); + } #if (AMREX_SPACEDIM <= 2) - if (!Geom().IsCartesian()) { - P_radial.setVal(0.0); - } + if (!Geom().IsCartesian()) { + getLevel(lev).P_radial.setVal(0.0); + } #endif #ifdef RADIATION - if (Radiation::rad_hydro_combined) { - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - rad_fluxes[dir]->setVal(0.0); - } - } + if (Radiation::rad_hydro_combined) { + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + getLevel(lev).rad_fluxes[dir]->setVal(0.0); + } + } #endif #ifdef REACTIONS - burn_weights.setVal(0.0); + if (castro::store_burn_weights) { + getLevel(lev).burn_weights.setVal(0.0); + } #endif - // For simplified SDC, we'll have garbage data if we - // attempt to use the lagged source terms (both reacting - // and non-reacting) from the last timestep, since that - // advance failed and we don't know if we can trust it. - // So we zero out both source term correctors. + // For simplified SDC, we'll have garbage data if we + // attempt to use the lagged source terms (both reacting + // and non-reacting) from the last timestep, since that + // advance failed and we don't know if we can trust it. + // So we zero out both source term correctors. - if (time_integration_method == SimplifiedSpectralDeferredCorrections) { - source_corrector.setVal(0.0, source_corrector.nGrow()); + if (time_integration_method == SimplifiedSpectralDeferredCorrections) { + getLevel(lev).source_corrector.setVal(0.0, getLevel(lev).source_corrector.nGrow()); #ifdef SIMPLIFIED_SDC #ifdef REACTIONS - MultiFab& SDC_react_new = get_new_data(Simplified_SDC_React_Type); - SDC_react_new.setVal(0.0, SDC_react_new.nGrow()); + MultiFab& SDC_react_new = getLevel(lev).get_new_data(Simplified_SDC_React_Type); + SDC_react_new.setVal(0.0, SDC_react_new.nGrow()); #endif #endif + } + } } @@ -306,7 +321,9 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, if (num_subcycles_remaining > max_subcycles) { amrex::Print() << std::endl - << " The subcycle mechanism requested " << num_subcycles_remaining << " subcycled timesteps, which is larger than the maximum of " << max_subcycles << "." << std::endl + << " The subcycle mechanism requested " << num_subcycles_remaining + << " subcycled timesteps, which is larger than the maximum of " + << max_subcycles << "." << std::endl << " If you would like to override this, increase the parameter castro.max_subcycles." << std::endl; amrex::Abort("Error: too many subcycles."); } @@ -315,9 +332,12 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, if (verbose && ParallelDescriptor::IOProcessor()) { std::cout << std::endl; - std::cout << Font::Bold << FGColor::Green << " Beginning subcycle " << sub_iteration + 1 << " starting at time " << subcycle_time + std::cout << Font::Bold << FGColor::Green + << " Beginning subcycle " << sub_iteration + 1 + << " starting at time " << subcycle_time << " with dt = " << dt_subcycle << ResetDisplay << std::endl; - std::cout << " Estimated number of subcycles remaining: " << num_subcycles_remaining << std::endl << std::endl; + std::cout << " Estimated number of subcycles remaining: " + << num_subcycles_remaining << std::endl << std::endl; } // Swap the time levels. Only do this after the first iteration; @@ -326,25 +346,26 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, if (do_swap) { - swap_state_time_levels(0.0); + for (int lev = level; lev <= max_level_to_advance; ++lev) { + getLevel(lev).swap_state_time_levels(0.0); #ifdef GRAVITY - if (do_grav) { - gravity->swapTimeLevels(level); - } + if (do_grav) { + getLevel(lev).gravity->swapTimeLevels(lev); + } #endif + } - } - else { - + } else { do_swap = true; - } // Set the relevant time levels. - for (int k = 0; k < num_state_type; k++) { - state[k].setTimeLevel(subcycle_time + dt_subcycle, dt_subcycle, 0.0); + for (int lev = level; lev <= max_level_to_advance; ++lev) { + for (int k = 0; k < num_state_type; k++) { + getLevel(lev).state[k].setTimeLevel(subcycle_time + dt_subcycle, dt_subcycle, 0.0); + } } // Do the advance and construct the relevant source terms. For CTU this @@ -455,7 +476,9 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, // Record the number of subcycles we took for diagnostic purposes. - num_subcycles_taken = sub_iteration; + for (int lev = level; lev <= max_level_to_advance; ++lev) { + getLevel(lev).num_subcycles_taken = sub_iteration; + } if (sub_iteration > 1) { @@ -465,14 +488,14 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, // we still have the last iteration's old data if we need // it later. - for (int k = 0; k < num_state_type; k++) { - - if (prev_state[k]->hasOldData()) { - state[k].replaceOldData(*prev_state[k]); + for (int lev = level; lev <= max_level_to_advance; ++lev) { + for (int k = 0; k < num_state_type; k++) { + if (getLevel(lev).prev_state[k]->hasOldData()) { + getLevel(lev).state[k].replaceOldData(*getLevel(lev).prev_state[k]); + } + getLevel(lev).state[k].setTimeLevel(time + dt, dt, 0.0); + getLevel(lev).prev_state[k]->setTimeLevel(time + dt, dt_subcycle, 0.0); } - state[k].setTimeLevel(time + dt, dt, 0.0); - prev_state[k]->setTimeLevel(time + dt, dt_subcycle, 0.0); - } // If we took more than one step and are going to do a reflux, @@ -484,8 +507,10 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, // reflux immediately following this, skip this if we're on the // finest level and this is not the last iteration. - if (!(amr_iteration < amr_ncycle && level == parent->finestLevel())) { - keep_prev_state = true; + for (int lev = level; lev <= max_level_to_advance; ++lev) { + if (!(amr_iteration < amr_ncycle && lev == parent->finestLevel())) { + getLevel(lev).keep_prev_state = true; + } } } diff --git a/Source/driver/Castro_advance_sdc.cpp b/Source/driver/Castro_advance_sdc.cpp index 83b226c7bc..01c5f653be 100644 --- a/Source/driver/Castro_advance_sdc.cpp +++ b/Source/driver/Castro_advance_sdc.cpp @@ -30,6 +30,9 @@ Castro::do_advance_sdc (Real time, // this is the new "formal" SDC integration routine. + amrex::ignore_unused(amr_iteration); + amrex::ignore_unused(amr_ncycle); + // this does the entire update in time for 1 SDC iteration. BL_PROFILE("Castro::do_advance_sdc()"); diff --git a/Source/driver/Castro_io.H b/Source/driver/Castro_io.H index d630db0ade..d47c74019f 100644 --- a/Source/driver/Castro_io.H +++ b/Source/driver/Castro_io.H @@ -1 +1,6 @@ +#ifndef CASTRO_IO_H +#define CASTRO_IO_H + extern std::string inputs_name; + +#endif diff --git a/Source/driver/Castro_io.cpp b/Source/driver/Castro_io.cpp index b002bf210e..75b5a45d54 100644 --- a/Source/driver/Castro_io.cpp +++ b/Source/driver/Castro_io.cpp @@ -574,16 +574,11 @@ Castro::writeJobInfo (const std::string& dir, const Real io_time) jobInfoFile << " Plotfile Information\n"; jobInfoFile << PrettyLine; - time_t now = time(nullptr); + const std::time_t now = time(nullptr); + jobInfoFile << "output date / time: " + << std::put_time(std::localtime(&now), "%c\n") << "\n"; - // Convert now to tm struct for local timezone - tm* localtm = localtime(&now); - jobInfoFile << "output date / time: " << asctime(localtm); - - char currentDir[FILENAME_MAX]; - if (getcwd(currentDir, FILENAME_MAX) != nullptr) { - jobInfoFile << "output dir: " << currentDir << "\n"; - } + jobInfoFile << "output dir: " << amrex::FileSystem::CurrentPath() << "\n"; jobInfoFile << "I/O time (s): " << io_time << "\n"; @@ -758,8 +753,39 @@ Castro::writeJobInfo (const std::string& dir, const Real io_time) } jobInfoFile << "\n"; + jobInfoFile << " amr.n_error_buf: "; + for (int lev = 1; lev <= max_level; lev++) { + int errbuf = parent->nErrorBuf(lev-1); + jobInfoFile << errbuf << " "; + } + jobInfoFile << "\n"; + + jobInfoFile << " amr.regrid_int: "; + for (int lev = 1; lev <= max_level; lev++) { + int regridint = parent->regridInt(lev-1); + jobInfoFile << regridint << " "; + } + jobInfoFile << "\n"; + + jobInfoFile << " amr.blocking_factor: "; + for (int lev = 1; lev <= max_level; lev++) { + IntVect bf = parent->blockingFactor(lev-1); + jobInfoFile << bf[0] << " "; + } + jobInfoFile << "\n"; + + jobInfoFile << " amr.max_grid_size: "; + for (int lev = 1; lev <= max_level; lev++) { + IntVect mgs = parent->maxGridSize(lev-1); + jobInfoFile << mgs[0] << " "; + } jobInfoFile << "\n\n"; + jobInfoFile << " amr.subcycling_mode: " << parent->subcyclingMode(); + + jobInfoFile << "\n\n"; + + // species info int mlen = 20; diff --git a/Source/driver/Castro_setup.cpp b/Source/driver/Castro_setup.cpp index c4a5327e64..0d3c8190d3 100644 --- a/Source/driver/Castro_setup.cpp +++ b/Source/driver/Castro_setup.cpp @@ -20,137 +20,153 @@ using std::string; using namespace amrex; -static Box the_same_box (const Box& b) { return b; } -static Box grow_box_by_one (const Box& b) { return amrex::grow(b,1); } +using BndryFunc = StateDescriptor::BndryFunc; -typedef StateDescriptor::BndryFunc BndryFunc; +namespace { -// -// Components are: -// Interior, Inflow, Outflow, Symmetry, SlipWall, NoSlipWall -// -static int scalar_bc[] = - { - amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_even, amrex::BCType::reflect_even, amrex::BCType::reflect_even - }; - -static int norm_vel_bc[] = - { - amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_odd, amrex::BCType::reflect_odd, amrex::BCType::reflect_odd - }; + Box the_same_box (const Box& b) { return b; } + Box grow_box_by_one (const Box& b) { return amrex::grow(b,1); } -static int tang_vel_bc[] = - { - amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_even, amrex::BCType::reflect_even, amrex::BCType::reflect_even - }; + // + // Components are: + // Interior, Inflow, Outflow, Symmetry, SlipWall, NoSlipWall + // + int scalar_bc[] = + { + amrex::BCType::int_dir, + amrex::BCType::ext_dir, + amrex::BCType::foextrap, + amrex::BCType::reflect_even, + amrex::BCType::reflect_even, + amrex::BCType::reflect_even + }; + + int norm_vel_bc[] = + { + amrex::BCType::int_dir, + amrex::BCType::ext_dir, + amrex::BCType::foextrap, + amrex::BCType::reflect_odd, + amrex::BCType::reflect_odd, + amrex::BCType::reflect_odd + }; + + int tang_vel_bc[] = + { + amrex::BCType::int_dir, + amrex::BCType::ext_dir, + amrex::BCType::foextrap, + amrex::BCType::reflect_even, + amrex::BCType::reflect_even, + amrex::BCType::reflect_even + }; #ifdef MHD -static int mag_field_bc[] = -{ - amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_even, amrex::BCType::foextrap, amrex::BCType::hoextrap -}; + int mag_field_bc[] = + { + amrex::BCType::int_dir, + amrex::BCType::ext_dir, + amrex::BCType::foextrap, + amrex::BCType::reflect_even, + amrex::BCType::foextrap, + amrex::BCType::hoextrap + }; #endif -static -void -set_scalar_bc (BCRec& bc, const BCRec& phys_bc) -{ - const int* lo_bc = phys_bc.lo(); - const int* hi_bc = phys_bc.hi(); - for (int i = 0; i < AMREX_SPACEDIM; i++) + void + set_scalar_bc (BCRec& bc, const BCRec& phys_bc) { - bc.setLo(i,scalar_bc[lo_bc[i]]); - bc.setHi(i,scalar_bc[hi_bc[i]]); + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + for (int i = 0; i < AMREX_SPACEDIM; i++) + { + bc.setLo(i,scalar_bc[lo_bc[i]]); + bc.setHi(i,scalar_bc[hi_bc[i]]); + } } -} -static -void -set_x_vel_bc(BCRec& bc, const BCRec& phys_bc) -{ - const int* lo_bc = phys_bc.lo(); - const int* hi_bc = phys_bc.hi(); - bc.setLo(0,norm_vel_bc[lo_bc[0]]); - bc.setHi(0,norm_vel_bc[hi_bc[0]]); + void + set_x_vel_bc(BCRec& bc, const BCRec& phys_bc) + { + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + bc.setLo(0,norm_vel_bc[lo_bc[0]]); + bc.setHi(0,norm_vel_bc[hi_bc[0]]); #if (AMREX_SPACEDIM >= 2) - bc.setLo(1,tang_vel_bc[lo_bc[1]]); - bc.setHi(1,tang_vel_bc[hi_bc[1]]); + bc.setLo(1,tang_vel_bc[lo_bc[1]]); + bc.setHi(1,tang_vel_bc[hi_bc[1]]); #endif #if (AMREX_SPACEDIM == 3) - bc.setLo(2,tang_vel_bc[lo_bc[2]]); - bc.setHi(2,tang_vel_bc[hi_bc[2]]); + bc.setLo(2,tang_vel_bc[lo_bc[2]]); + bc.setHi(2,tang_vel_bc[hi_bc[2]]); #endif -} + } -static -void -set_y_vel_bc(BCRec& bc, const BCRec& phys_bc) -{ - const int* lo_bc = phys_bc.lo(); - const int* hi_bc = phys_bc.hi(); - bc.setLo(0,tang_vel_bc[lo_bc[0]]); - bc.setHi(0,tang_vel_bc[hi_bc[0]]); + void + set_y_vel_bc(BCRec& bc, const BCRec& phys_bc) + { + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + bc.setLo(0,tang_vel_bc[lo_bc[0]]); + bc.setHi(0,tang_vel_bc[hi_bc[0]]); #if (AMREX_SPACEDIM >= 2) - bc.setLo(1,norm_vel_bc[lo_bc[1]]); - bc.setHi(1,norm_vel_bc[hi_bc[1]]); + bc.setLo(1,norm_vel_bc[lo_bc[1]]); + bc.setHi(1,norm_vel_bc[hi_bc[1]]); #endif #if (AMREX_SPACEDIM == 3) - bc.setLo(2,tang_vel_bc[lo_bc[2]]); - bc.setHi(2,tang_vel_bc[hi_bc[2]]); + bc.setLo(2,tang_vel_bc[lo_bc[2]]); + bc.setHi(2,tang_vel_bc[hi_bc[2]]); #endif -} + } -static -void -set_z_vel_bc(BCRec& bc, const BCRec& phys_bc) -{ - const int* lo_bc = phys_bc.lo(); - const int* hi_bc = phys_bc.hi(); - bc.setLo(0,tang_vel_bc[lo_bc[0]]); - bc.setHi(0,tang_vel_bc[hi_bc[0]]); + void + set_z_vel_bc(BCRec& bc, const BCRec& phys_bc) + { + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + bc.setLo(0,tang_vel_bc[lo_bc[0]]); + bc.setHi(0,tang_vel_bc[hi_bc[0]]); #if (AMREX_SPACEDIM >= 2) - bc.setLo(1,tang_vel_bc[lo_bc[1]]); - bc.setHi(1,tang_vel_bc[hi_bc[1]]); + bc.setLo(1,tang_vel_bc[lo_bc[1]]); + bc.setHi(1,tang_vel_bc[hi_bc[1]]); #endif #if (AMREX_SPACEDIM == 3) - bc.setLo(2,norm_vel_bc[lo_bc[2]]); - bc.setHi(2,norm_vel_bc[hi_bc[2]]); + bc.setLo(2,norm_vel_bc[lo_bc[2]]); + bc.setHi(2,norm_vel_bc[hi_bc[2]]); #endif -} - + } #ifdef MHD -static -void -set_mag_field_bc(BCRec& bc, const BCRec& phys_bc) -{ - const int* lo_bc = phys_bc.lo(); - const int* hi_bc = phys_bc.hi(); - for (int i = 0; i < AMREX_SPACEDIM; i++) + void + set_mag_field_bc(BCRec& bc, const BCRec& phys_bc) { - bc.setLo(i, mag_field_bc[lo_bc[i]]); - bc.setHi(i, mag_field_bc[hi_bc[i]]); + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + for (int i = 0; i < AMREX_SPACEDIM; i++) + { + bc.setLo(i, mag_field_bc[lo_bc[i]]); + bc.setHi(i, mag_field_bc[hi_bc[i]]); + } } -} #endif -// In some cases we want to replace inflow boundaries with -// first-order extrapolation boundaries. This is intended to -// be used for state data that the user is not going to -// provide inflow boundary conditions for, like gravity -// and reactions, and it works in conjunction with the -// generic_fill boundary routine. + // In some cases we want to replace inflow boundaries with + // first-order extrapolation boundaries. This is intended to + // be used for state data that the user is not going to + // provide inflow boundary conditions for, like gravity + // and reactions, and it works in conjunction with the + // generic_fill boundary routine. -static -void -replace_inflow_bc (BCRec& bc) -{ - for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - if (bc.lo(dir) == amrex::BCType::ext_dir) { - bc.setLo(dir, amrex::BCType::foextrap); - } - if (bc.hi(dir) == amrex::BCType::ext_dir) { - bc.setHi(dir, amrex::BCType::foextrap); + void + replace_inflow_bc (BCRec& bc) + { + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + if (bc.lo(dir) == amrex::BCType::ext_dir) { + bc.setLo(dir, amrex::BCType::foextrap); + } + if (bc.hi(dir) == amrex::BCType::ext_dir) { + bc.setHi(dir, amrex::BCType::foextrap); + } } } } diff --git a/Source/driver/Make.package b/Source/driver/Make.package index c4ca1e4313..85febd5698 100644 --- a/Source/driver/Make.package +++ b/Source/driver/Make.package @@ -1,7 +1,7 @@ # these are the files that should be needed for any Castro build CEXE_sources += Castro.cpp -CEXE_sources += runparams_defaults.cpp +CEXE_sources += runtime_params.cpp CEXE_sources += Castro_advance.cpp CEXE_sources += Castro_advance_ctu.cpp ifeq ($(USE_TRUE_SDC), TRUE) diff --git a/Source/driver/castro_limits.H b/Source/driver/castro_limits.H index 0111b43432..31ddeb93a9 100644 --- a/Source/driver/castro_limits.H +++ b/Source/driver/castro_limits.H @@ -1,6 +1,6 @@ #ifndef CASTRO_LIMITS_H #define CASTRO_LIMITS_H -#define MAX_LEV 15 +constexpr int MAX_LEV{15}; #endif diff --git a/Source/driver/parse_castro_params.py b/Source/driver/parse_castro_params.py index aedeea6624..89f91a8ea7 100755 --- a/Source/driver/parse_castro_params.py +++ b/Source/driver/parse_castro_params.py @@ -1,9 +1,9 @@ #!/usr/bin/env python3 -""" -This script parses the list of C++ runtime parameters and writes the -necessary header files and Fortran routines to make them available -in Castro's C++ routines. +"""This script parses the list of C++ runtime parameters and writes +the necessary header and source files to make them available in +Castro's C++ routines. They are available in 2 ways: as global +parameter and in the form of a single struct. parameters have the format: @@ -41,9 +41,6 @@ -- name_params.H (for castro, included in Castro.H): sets up the namespace and extern parameters - -- name_declares.H (for castro, included in Castro.cpp): - declares the runtime parameters - -- name_queries.H (for castro, included in Castro.cpp): does the parmparse query to override the default in C++ @@ -51,6 +48,9 @@ this tests the current value against the default and outputs into a file + -- runtime_params.cpp + has the actual definition of the variables (without extern) + """ import argparse @@ -129,12 +129,12 @@ def read_param_file(infile): return params -def write_headers(params, out_directory, struct_name): +def write_headers_and_source(params, out_directory, struct_name): # output # find all the namespaces - namespaces = {q.namespace for q in params} + namespaces = sorted({q.namespace for q in params}) for nm in namespaces: @@ -142,32 +142,6 @@ def write_headers(params, out_directory, struct_name): # sort by repr since None may be present ifdefs = sorted({q.ifdef for q in params_nm}, key=repr) - # write name_declares.H - try: - cd = open(f"{out_directory}/{nm}_declares.H", "w", encoding="UTF-8") - except OSError: - sys.exit(f"unable to open {nm}_declares.H for writing") - - cd.write(CWARNING) - cd.write(f"#ifndef {nm.upper()}_DECLARES_H\n") - cd.write(f"#define {nm.upper()}_DECLARES_H\n") - - cd.write("\n") - cd.write(f"namespace {nm} {{\n") - - for ifdef in ifdefs: - if ifdef is None: - for p in [q for q in params_nm if q.ifdef is None]: - cd.write(p.get_declare_string()) - else: - cd.write(f"#ifdef {ifdef}\n") - for p in [q for q in params_nm if q.ifdef == ifdef]: - cd.write(p.get_declare_string()) - cd.write("#endif\n") - cd.write("}\n\n") - cd.write("#endif\n") - cd.close() - # write name_params.H try: cp = open(f"{out_directory}/{nm}_params.H", "w", encoding="UTF-8") @@ -238,7 +212,44 @@ def write_headers(params, out_directory, struct_name): jo.close() + # write a single C++ source file that actually defines the parameters + # (one file for all namespaces) + try: + pf = open(f"{out_directory}/runtime_params.cpp", "w", encoding="UTF-8") + except OSError: + sys.exit(f"unable to open runtime_params.cpp") + + pf.write("#include \n") + pf.write("#include \n") + pf.write("#include \n\n") + + for nm in namespaces: + pf.write(f"#include <{nm}_params.H>\n") + pf.write("\n") + + for nm in namespaces: + params_nm = [q for q in params if q.namespace == nm] + # sort by repr since None may be present + ifdefs = sorted({q.ifdef for q in params_nm}, key=repr) + + pf.write(f"namespace {nm} {{\n") + + for ifdef in ifdefs: + if ifdef is None: + for p in [q for q in params_nm if q.ifdef is None]: + pf.write(p.get_declare_string()) + else: + pf.write(f"#ifdef {ifdef}\n") + for p in [q for q in params_nm if q.ifdef == ifdef]: + pf.write(p.get_declare_string()) + pf.write("#endif\n") + pf.write("}\n\n") + + pf.close() + # now write a single file that contains all of the parameter structs + # to minimize padding, we want to sort on type + try: sf = open(f"{out_directory}/{struct_name}_type.H", "w", encoding="UTF-8") except OSError: @@ -255,20 +266,25 @@ def write_headers(params, out_directory, struct_name): params_nm = [q for q in params if q.namespace == nm] # sort by repr since None may be present ifdefs = sorted({q.ifdef for q in params_nm}, key=repr) - sf.write(f"struct {nm}_t {{\n") print("namespace = ", nm) for ifdef in ifdefs: + params_if = [q for q in params_nm if q.ifdef == ifdef] + types = sorted(set(q.dtype for q in params_if)) + if ifdef is None: - for p in [q for q in params_nm if q.ifdef is None]: - sf.write(p.get_struct_entry()) + for tt in types: + params_type = [q for q in params_if if q.dtype == tt] + for p in params_type: + sf.write(p.get_struct_entry()) else: sf.write(f"#ifdef {ifdef}\n") - for p in [q for q in params_nm if q.ifdef == ifdef]: - sf.write(p.get_struct_entry()) + for tt in types: + params_type = [q for q in params_if if q.dtype == tt] + for p in params_type: + sf.write(p.get_struct_entry()) sf.write("#endif\n") - sf.write("};\n\n") # now the parent struct @@ -295,7 +311,7 @@ def main(): args = parser.parse_args() p = read_param_file(args.input_file[0]) - write_headers(p, args.o, args.s) + write_headers_and_source(p, args.o, args.s) if __name__ == "__main__": main() diff --git a/Source/driver/runparams_defaults.cpp b/Source/driver/runparams_defaults.cpp deleted file mode 100644 index c9f5ed1de6..0000000000 --- a/Source/driver/runparams_defaults.cpp +++ /dev/null @@ -1,6 +0,0 @@ -#include - -#include - -#include - diff --git a/Source/driver/sum_integrated_quantities.cpp b/Source/driver/sum_integrated_quantities.cpp index 5a41f9242c..a92c65d885 100644 --- a/Source/driver/sum_integrated_quantities.cpp +++ b/Source/driver/sum_integrated_quantities.cpp @@ -375,13 +375,9 @@ Castro::sum_integrated_quantities () data_log1 << std::setw(intwidth) << timestep; - if (time == 0.0_rt) { - data_log1 << std::fixed; - } - else if (time < 1.e-4_rt || time > 1.e4_rt) { + if (time < 1.e-4_rt || time > 1.e4_rt) { data_log1 << std::scientific; - } - else { + } else { data_log1 << std::fixed; } @@ -515,10 +511,7 @@ Castro::sum_integrated_quantities () log << std::setw(intwidth) << timestep; - if (time == 0.0_rt) { - log << std::fixed; - } - else if (time < 1.e-4_rt || time > 1.e4_rt) { + if (time < 1.e-4_rt || time > 1.e4_rt) { log << std::scientific; } else { @@ -616,13 +609,9 @@ Castro::sum_integrated_quantities () log << std::setw(intwidth) << timestep; - if (time == 0.0_rt) { - log << std::fixed; - } - else if (time < 1.e-4_rt || time > 1.e4_rt) { + if (time < 1.e-4_rt || time > 1.e4_rt) { log << std::scientific; - } - else { + } else { log << std::fixed; } @@ -717,13 +706,9 @@ Castro::sum_integrated_quantities () log << std::setw(intwidth) << timestep; - if (time == 0.0_rt) { - log << std::fixed; - } - else if (time < 1.e-4_rt || time > 1.e4_rt) { + if (time < 1.e-4_rt || time > 1.e4_rt) { log << std::scientific; - } - else { + } else { log << std::fixed; } diff --git a/Source/gravity/Castro_gravity.H b/Source/gravity/Castro_gravity.H index b9caa2a1bd..3943d55f90 100644 --- a/Source/gravity/Castro_gravity.H +++ b/Source/gravity/Castro_gravity.H @@ -1,3 +1,6 @@ +#ifndef CASTRO_GRAVITY_H +#define CASTRO_GRAVITY_H + /// /// Construct gravitational field at old timestep /// @@ -35,3 +38,4 @@ /// void construct_new_gravity_source(amrex::MultiFab& source, amrex::MultiFab& state_old, amrex::MultiFab& state_new, amrex::Real time, amrex::Real dt); +#endif diff --git a/Source/gravity/Castro_gravity.cpp b/Source/gravity/Castro_gravity.cpp index ed8a2406e0..bf9dc4f6b7 100644 --- a/Source/gravity/Castro_gravity.cpp +++ b/Source/gravity/Castro_gravity.cpp @@ -388,7 +388,7 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in, Real SrE{}; - if (castro::grav_source_type == 1 || castro::grav_source_type == 2) { + if (castro::grav_source_type == 1 || castro::grav_source_type == 2) { // NOLINT(bugprone-branch-clone) // Src = rho u dot g, evaluated with all quantities at t^n diff --git a/Source/gravity/Gravity.H b/Source/gravity/Gravity.H index af7ff177fa..0f6aba88ac 100644 --- a/Source/gravity/Gravity.H +++ b/Source/gravity/Gravity.H @@ -127,7 +127,7 @@ public: /// @param time Current time /// @param multi_level boolean, do we iterate over all levels or mask off fine grids? /// - void set_mass_offset(amrex::Real time, bool multi_level=true); + void set_mass_offset(amrex::Real time, bool multi_level=true) const; /// /// Return ``grad_phi_prev`` at given level @@ -242,17 +242,17 @@ public: /// /// Compute the difference between level and composite solves /// -/// @param level level index -/// @param comp_phi MultiFab containing computed phi -/// @param comp_gphi Vector of MultiFabs containing computed grad phi -/// @param cml_phi MultiFab, computed minus level phi -/// @param cml_gphi Vector of MultiFabs, computed minus level grad phi +/// @param level level index +/// @param comp_phi MultiFab containing computed phi +/// @param comp_gphi Vector of MultiFabs containing computed grad phi +/// @param comp_minus_level_phi MultiFab, computed minus level phi +/// @param comp_minus_level_grad_phi Vector of MultiFabs, computed minus level grad phi /// void create_comp_minus_level_grad_phi(int level, amrex::MultiFab& comp_phi, const amrex::Vector& comp_gphi, - amrex::MultiFab& cml_phi, - amrex::Vector >& cml_gphi); + amrex::MultiFab& comp_minus_level_phi, + amrex::Vector >& comp_minus_level_grad_phi); /// @@ -348,7 +348,7 @@ public: /// @param radial_grav Radial gravity /// @param grav_vector Gravity vector /// - void interpolate_monopole_grav(int level, RealVector& radial_grav, amrex::MultiFab& grav_vector); + void interpolate_monopole_grav(int level, RealVector& radial_grav, amrex::MultiFab& grav_vector) const; /// /// Integrate radially outward to find radial mass distribution @@ -368,7 +368,7 @@ public: #ifdef GR_GRAV RealVector& radial_pres, #endif - int n1d, int level); + int n1d, int level) const; /// /// Implement multipole boundary conditions @@ -383,7 +383,7 @@ public: /// /// Initialize multipole gravity /// - void init_multipole_grav(); + void init_multipole_grav() const; #if (AMREX_SPACEDIM == 3) @@ -480,7 +480,7 @@ public: /// @param Rhs MultiFab /// @param coeffs Vector of MultiFabs /// - void applyMetricTerms(int level,amrex::MultiFab& Rhs, const amrex::Vector& coeffs); + void applyMetricTerms(int level,amrex::MultiFab& Rhs, const amrex::Vector& coeffs) const; /// /// @@ -488,7 +488,7 @@ public: /// @param level Index of level /// @param cc Cell-centered data /// - void unweight_cc(int level,amrex::MultiFab& cc); + void unweight_cc(int level,amrex::MultiFab& cc) const; /// /// @@ -496,7 +496,7 @@ public: /// @param level Index of level /// @param edges Edge-based data /// - void unweight_edges(int level, const amrex::Vector& edges); + void unweight_edges(int level, const amrex::Vector& edges) const; #endif /// @@ -504,7 +504,7 @@ public: /// @param phi Gravitational potential /// @param grav_vector Gravity vector /// - void add_pointmass_to_gravity (int level, amrex::MultiFab& phi, amrex::MultiFab& grav_vector); + void add_pointmass_to_gravity (int level, amrex::MultiFab& phi, amrex::MultiFab& grav_vector) const; /// /// Get the rhs @@ -546,7 +546,7 @@ public: const amrex::Vector >& grad_phi, const amrex::Vector& res, const amrex::MultiFab* const crse_bcdata, - amrex::Real rel_eps, amrex::Real abs_eps); + amrex::Real rel_eps, amrex::Real abs_eps) const; /// diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index b2f0f90d4e..72a2565935 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -39,7 +39,7 @@ Real Gravity::mass_offset = 0.0; // ************************************************************************************** // -static Real Ggravity = 0.; +const Real Ggravity = 4.0 * M_PI * C::Gconst; /// /// Multipole gravity data @@ -264,13 +264,6 @@ Gravity::read_params () } - Ggravity = 4.0 * M_PI * C::Gconst; - if (gravity::verbose > 1 && ParallelDescriptor::IOProcessor()) - { - std::cout << "Getting Gconst from constants: " << C::Gconst << std::endl; - std::cout << "Using " << Ggravity << " for 4 pi G in Gravity.cpp " << std::endl; - } - done = true; } } @@ -1325,7 +1318,7 @@ Gravity::test_composite_phi (int crse_level) } void -Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab& grav_vector) +Gravity::interpolate_monopole_grav(int level, RealVector& radial_grav, MultiFab& grav_vector) const { BL_PROFILE("Gravity::interpolate_monopole_grav()"); @@ -1441,7 +1434,7 @@ Gravity::compute_radial_mass(const Box& bx, #ifdef GR_GRAV RealVector& radial_pres_local, #endif - int n1d, int level) + int n1d, int level) const { const Geometry& geom = parent->Geom(level); @@ -1606,7 +1599,7 @@ Gravity::compute_radial_mass(const Box& bx, } void -Gravity::init_multipole_grav() +Gravity::init_multipole_grav() const { if (gravity::lnum < 0) { amrex::Abort("lnum negative"); @@ -1698,7 +1691,7 @@ Gravity::init_multipole_grav() multipole::parity_q0(l) = 1.0_rt; if (l % 2 != 0) { - if (AMREX_SPACEDIM == 3 && (multipole::doReflectionLo(2) || multipole::doReflectionHi(2))) { + if (AMREX_SPACEDIM == 3 && (multipole::doReflectionLo(2) || multipole::doReflectionHi(2))) { // NOLINT(bugprone-branch-clone) multipole::parity_q0(l) = 0.0_rt; } else if (AMREX_SPACEDIM == 2 && parent->Geom(0).Coord() == 1) { @@ -1732,7 +1725,7 @@ Gravity::init_multipole_grav() if (AMREX_SPACEDIM == 3) { multipole::parity_qC_qS(l,m) = 1.0_rt; } - else if (AMREX_SPACEDIM == 2 && parent->Geom(0).Coord() == 1) { + else if (AMREX_SPACEDIM == 2 && parent->Geom(0).Coord() == 1) { // NOLINT(bugprone-branch-clone) multipole::parity_qC_qS(l,m) = 0.0_rt; } else if (AMREX_SPACEDIM == 1 && parent->Geom(0).Coord() == 2) { @@ -2854,7 +2847,7 @@ Gravity::fill_direct_sum_BCs(int crse_level, int fine_level, const Vector& coeffs) +Gravity::applyMetricTerms(int level, MultiFab& Rhs, const Vector& coeffs) const { BL_PROFILE("Gravity::applyMetricTerms()"); @@ -2884,7 +2877,7 @@ Gravity::applyMetricTerms(int level, MultiFab& Rhs, const Vector& coe } void -Gravity::unweight_cc(int level, MultiFab& cc) +Gravity::unweight_cc(int level, MultiFab& cc) const { BL_PROFILE("Gravity::unweight_cc()"); @@ -2903,7 +2896,7 @@ Gravity::unweight_cc(int level, MultiFab& cc) } void -Gravity::unweight_edges(int level, const Vector& edges) +Gravity::unweight_edges(int level, const Vector& edges) const { BL_PROFILE("Gravity::unweight_edges()"); @@ -2956,7 +2949,7 @@ Gravity::make_mg_bc () } void -Gravity::set_mass_offset (Real time, bool multi_level) +Gravity::set_mass_offset (Real time, bool multi_level) const { BL_PROFILE("Gravity::set_mass_offset()"); @@ -2973,9 +2966,9 @@ Gravity::set_mass_offset (Real time, bool multi_level) { for (int lev = 0; lev <= parent->finestLevel(); lev++) { auto* cs = dynamic_cast(&parent->getLevel(lev)); - if (cs != nullptr) { - mass_offset += cs->volWgtSum("density", time); - } else { + if (cs != nullptr) { + mass_offset += cs->volWgtSum("density", time); + } else { amrex::Abort("unable to access volWgtSum"); } } @@ -2983,7 +2976,11 @@ Gravity::set_mass_offset (Real time, bool multi_level) else { auto* cs = dynamic_cast(&parent->getLevel(0)); - mass_offset = cs->volWgtSum("density", time, false, false); // do not mask off fine grids + if (cs != nullptr) { + mass_offset = cs->volWgtSum("density", time, false, false); // do not mask off fine grids + } else { + amrex::Abort("unable to access volWgtSum"); + } } mass_offset = mass_offset / geom.ProbSize(); @@ -3006,7 +3003,7 @@ Gravity::set_mass_offset (Real time, bool multi_level) } void -Gravity::add_pointmass_to_gravity (int level, MultiFab& phi, MultiFab& grav_vector) +Gravity::add_pointmass_to_gravity (int level, MultiFab& phi, MultiFab& grav_vector) const { BL_PROFILE("Gravity::add_pointmass_to_gravity()"); @@ -3101,8 +3098,7 @@ Gravity::make_radial_gravity(int level, Real time, RealVector& radial_grav) // Create MultiFab with NUM_STATE components and no ghost cells MultiFab S(grids[lev],dmap[lev],NUM_STATE,0); - if ( eps == 0.0 ) - { + if ( eps == 0.0 ) { // NOLINT(bugprone-branch-clone,-warnings-as-errors) // Old and new time are identical; this should only happen if // dt is smaller than roundoff compared to the current time, // in which case we're probably in trouble anyway, @@ -3720,7 +3716,7 @@ Gravity::actual_solve_with_mlmg (int crse_level, int fine_level, const amrex::Vector >& grad_phi, const amrex::Vector& res, const amrex::MultiFab* const crse_bcdata, - amrex::Real rel_eps, amrex::Real abs_eps) + amrex::Real rel_eps, amrex::Real abs_eps) const { BL_PROFILE("Gravity::actual_solve_with_mlmg()"); diff --git a/Source/gravity/Make.package b/Source/gravity/Make.package index f94894307f..6708aab63d 100644 --- a/Source/gravity/Make.package +++ b/Source/gravity/Make.package @@ -2,7 +2,6 @@ # this is included if USE_GRAV = TRUE CEXE_sources += Gravity.cpp -CEXE_sources += gravity_params.cpp CEXE_headers += Gravity.H CEXE_headers += Gravity_util.H CEXE_headers += Castro_gravity.H diff --git a/Source/gravity/gravity_params.cpp b/Source/gravity/gravity_params.cpp deleted file mode 100644 index c80baeadc7..0000000000 --- a/Source/gravity/gravity_params.cpp +++ /dev/null @@ -1,9 +0,0 @@ -#include - -#include - -#include - -#include - - diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index b98293b7c0..38e18b0cac 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -9,9 +9,6 @@ using namespace amrex; void Castro::consup_hydro(const Box& bx, -#ifdef SHOCK_VAR - Array4 const& shk, -#endif Array4 const& U_new, Array4 const& flux0, Array4 const& qx, @@ -68,11 +65,6 @@ Castro::consup_hydro(const Box& bx, U_new(i,j,k,n) = U_new(i,j,k,n) - dt * pdu; -#ifdef SHOCK_VAR - } else if (n == USHK) { - U_new(i,j,k,USHK) = shk(i,j,k); -#endif - } else if (n == UMX) { // Add gradp term to momentum equation -- only for axisymmetric // coords (and only for the radial flux). diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index c1a29ae0c2..b7db069e41 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -14,8 +14,11 @@ using namespace amrex; advance_status -Castro::construct_ctu_hydro_source(Real time, Real dt) +Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-convert-member-functions-to-static) { + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + advance_status status {}; if (!do_hydro) { @@ -1190,9 +1193,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #endif consup_hydro(bx, -#ifdef SHOCK_VAR - shk_arr, -#endif update_arr, flx_arr, qx_arr, #if AMREX_SPACEDIM >= 2 diff --git a/Source/hydro/Castro_hybrid.cpp b/Source/hydro/Castro_hybrid.cpp index afff6949b3..b1bdabc4f9 100644 --- a/Source/hydro/Castro_hybrid.cpp +++ b/Source/hydro/Castro_hybrid.cpp @@ -84,7 +84,7 @@ Castro::construct_new_hybrid_source(MultiFab& source, MultiFab& state_old, Multi void -Castro::fill_hybrid_hydro_source(MultiFab& sources, MultiFab& state_in, Real mult_factor) +Castro::fill_hybrid_hydro_source(MultiFab& sources, const MultiFab& state_in, Real mult_factor) { BL_PROFILE("Castro::fill_hybrid_hydro_source()"); @@ -97,7 +97,7 @@ Castro::fill_hybrid_hydro_source(MultiFab& sources, MultiFab& state_in, Real mul { const Box& bx = mfi.tilebox(); - auto u = state_in.array(mfi); + const auto u = state_in.array(mfi); auto src = sources.array(mfi); amrex::ParallelFor(bx, diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index 94093c6392..776aee48cc 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -135,13 +135,14 @@ /// /// A multidimensional shock detection algorithm /// -/// @param bx the box to operate over -/// @param q_arr the primitive variable state -/// @param shk the shock flag (1 = shock, 0 = no shock) +/// @param bx the box to operate over +/// @param q_arr the primitive variable state +/// @param U_scr_arr the conservative state sources +/// @param shk the shock flag (1 = shock, 0 = no shock) /// void shock(const amrex::Box& bx, amrex::Array4 const& q_arr, - amrex::Array4 const& q_src_arr, + amrex::Array4 const& U_src_arr, amrex::Array4 const& shk); /// @@ -710,9 +711,6 @@ /// @param dt timestep /// void consup_hydro(const amrex::Box& bx, -#ifdef SHOCK_VAR - amrex::Array4 const& shk, -#endif amrex::Array4 const& U_new, amrex::Array4 const& flux0, amrex::Array4 const& qx, @@ -856,7 +854,7 @@ /// @param source MultiFab to save source to /// @param mult_factor /// - void fill_hybrid_hydro_source(amrex::MultiFab& state, amrex::MultiFab& source, const amrex::Real mult_factor); + void fill_hybrid_hydro_source(amrex::MultiFab& source, const amrex::MultiFab& state, const amrex::Real mult_factor); /// diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 23a004c719..eb6314769c 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -199,7 +199,7 @@ Castro::mol_ppm_reconstruct(const Box& bx, void -Castro::mol_consup(const Box& bx, +Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-functions-to-static) #ifdef SHOCK_VAR Array4 const& shk, #endif diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 4d7e0f32db..fbcc8c6171 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -21,6 +21,8 @@ void Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) { + amrex::ignore_unused(time); + #ifdef RADIATION amrex::Abort("Error: radiation not supported for the MOL hydro source term"); #else @@ -648,6 +650,8 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) } Array4 pradial_fab = pradial.array(); +#endif +#if AMREX_SPACEDIM == 1 Array4 const qex_arr = qe[0].array(); #endif diff --git a/Source/hydro/advection_util.H b/Source/hydro/advection_util.H index 1669869dc4..cde9eaaba6 100644 --- a/Source/hydro/advection_util.H +++ b/Source/hydro/advection_util.H @@ -23,6 +23,9 @@ void src_to_prim(int i, int j, int k, const Real dt, #endif Array4 const& srcQ) { + + amrex::ignore_unused(dt); + for (int n = 0; n < NQSRC; ++n) { srcQ(i,j,k,n) = 0.0_rt; } diff --git a/Source/hydro/slope.H b/Source/hydro/slope.H index 212f00c123..9a8110f157 100644 --- a/Source/hydro/slope.H +++ b/Source/hydro/slope.H @@ -1,3 +1,6 @@ +#ifndef SLOPE_H +#define SLOPE_H + #include #ifdef RADIATION @@ -239,3 +242,4 @@ pslope(const Real* rho, const Real* p, const Real* src, const Real flatn, } } +#endif diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index 9c85de9df1..acbc4947df 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -64,7 +64,7 @@ Castro::trans_single(const Box& bx, void -Castro::actual_trans_single(const Box& bx, +Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member-functions-to-static) int idir_t, int idir_n, int d, Array4 const& q_arr, Array4 const& qo_arr, @@ -521,7 +521,7 @@ Castro::trans_final(const Box& bx, void -Castro::actual_trans_final(const Box& bx, +Castro::actual_trans_final(const Box& bx, // NOLINT(readability-convert-member-functions-to-static) int idir_n, int idir_t1, int idir_t2, int d, Array4 const& q_arr, Array4 const& qo_arr, diff --git a/Source/mhd/Castro_mhd.H b/Source/mhd/Castro_mhd.H index 3451fefe37..99821b31c7 100644 --- a/Source/mhd/Castro_mhd.H +++ b/Source/mhd/Castro_mhd.H @@ -1,3 +1,5 @@ +#ifndef CASTRO_MHD_H +#define CASTRO_MHD_H advance_status construct_ctu_mhd_source(amrex::Real time, amrex::Real dt); @@ -31,7 +33,7 @@ void consup_mhd(const amrex::Box& bx, const Real dt, - amrex::Array4 const& update, + amrex::Array4 const& U_new, amrex::Array4 const& flux0, amrex::Array4 const& flux1, amrex::Array4 const& flux2); @@ -105,3 +107,4 @@ const int dir); +#endif diff --git a/Source/mhd/mhd_util.H b/Source/mhd/mhd_util.H index 7b7cce3eba..7f9097e4f5 100644 --- a/Source/mhd/mhd_util.H +++ b/Source/mhd/mhd_util.H @@ -23,7 +23,7 @@ eos_soundspeed_mhd(Real& c, Real as, Real ca, Real bd) { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -qflux(Real* qflx, Real* flx, Real* q_zone) { +qflux(Real* qflx, const Real* flx, const Real* q_zone) { // Calculate the C to P Jacobian applied to the fluxes @@ -78,7 +78,7 @@ qflux(Real* qflx, Real* flx, Real* q_zone) { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -electric(Real* q_zone, Real& E_zone, const int comp) { +electric(const Real* q_zone, Real& E_zone, const int comp) { // this takes the cell-center primitive state, q_zone, and computes the cell-center // electric field, E_zone, using Faraday's law: @@ -99,7 +99,7 @@ electric(Real* q_zone, Real& E_zone, const int comp) { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -ConsToPrim(Real* q_zone, Real* U_zone) { +ConsToPrim(Real* q_zone, const Real* U_zone) { // calculate the primitive variables from the conserved // U has NUM_STATE+3 components diff --git a/Source/problems/Problem.H b/Source/problems/Problem.H index a32930b169..0a385762c9 100644 --- a/Source/problems/Problem.H +++ b/Source/problems/Problem.H @@ -1,3 +1,6 @@ +#ifndef PROBLEM_H +#define PROBLEM_H + /* problem-specific Castro:: declarations go here */ #ifdef DO_PROBLEM_POST_TIMESTEP @@ -11,3 +14,5 @@ void problem_post_restart(); #ifdef DO_PROBLEM_POST_INIT void problem_post_init(); #endif + +#endif diff --git a/Source/problems/Problem_Derive.H b/Source/problems/Problem_Derive.H index 0db4669c53..0cfcfe6a31 100644 --- a/Source/problems/Problem_Derive.H +++ b/Source/problems/Problem_Derive.H @@ -1 +1,6 @@ +#ifndef PROBLEM_DERIVE_H +#define PROBLEM_DERIVE_H + /* problem-specific stuff goes here */ + +#endif diff --git a/Source/problems/Problem_Derives.H b/Source/problems/Problem_Derives.H index ba416c112e..e09391c664 100644 --- a/Source/problems/Problem_Derives.H +++ b/Source/problems/Problem_Derives.H @@ -1 +1,6 @@ +#ifndef PROBLEM_DERIVES_H +#define PROBLEM_DERIVES_H + // problem-specific derived variables would be put here + +#endif diff --git a/Source/problems/ambient_fill.cpp b/Source/problems/ambient_fill.cpp index c5cb7cc701..fe16818df4 100644 --- a/Source/problems/ambient_fill.cpp +++ b/Source/problems/ambient_fill.cpp @@ -6,9 +6,10 @@ void ambient_denfill(const Box& bx, Array4 const& state, Geometry const& geom, const Vector& bcr) { + // Make a copy of the BC that can be passed by value to the ParallelFor. - BCRec bc = bcr[0]; + BCRec bc = bcr[URHO]; const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); @@ -62,12 +63,11 @@ void ambient_fill(const Box& bx, Array4 const& state, Geometry const& geom, const Vector& bcr) { - // Copy BCs to an Array1D so they can be passed by value to the ParallelFor. - Array1D bcs; - for (int n = 0; n < NUM_STATE; ++n) { - bcs(n) = bcr[n]; - } + // Make a copy of the BC that can be passed by value to the ParallelFor. + // even though this fills all components, we only check the density BCs + + BCRec bc = bcr[URHO]; const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); @@ -76,22 +76,22 @@ ambient_fill(const Box& bx, Array4 const& state, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { bool ambient_x_lo = (castro::ambient_fill_dir == 0 || castro::ambient_fill_dir == -1) && - (bcs(URHO).lo(0) == amrex::BCType::foextrap || bcs(URHO).lo(0) == amrex::BCType::hoextrap); + (bc.lo(0) == amrex::BCType::foextrap || bc.lo(0) == amrex::BCType::hoextrap); bool ambient_x_hi = (castro::ambient_fill_dir == 0 || castro::ambient_fill_dir == -1) && - (bcs(URHO).hi(0) == amrex::BCType::foextrap || bcs(URHO).hi(0) == amrex::BCType::hoextrap); + (bc.hi(0) == amrex::BCType::foextrap || bc.hi(0) == amrex::BCType::hoextrap); #if AMREX_SPACEDIM >= 2 bool ambient_y_lo = (castro::ambient_fill_dir == 1 || castro::ambient_fill_dir == -1) && - (bcs(URHO).lo(1) == amrex::BCType::foextrap || bcs(URHO).lo(1) == amrex::BCType::hoextrap); + (bc.lo(1) == amrex::BCType::foextrap || bc.lo(1) == amrex::BCType::hoextrap); bool ambient_y_hi = (castro::ambient_fill_dir == 1 || castro::ambient_fill_dir == -1) && - (bcs(URHO).hi(1) == amrex::BCType::foextrap || bcs(URHO).hi(1) == amrex::BCType::hoextrap); + (bc.hi(1) == amrex::BCType::foextrap || bc.hi(1) == amrex::BCType::hoextrap); #endif #if AMREX_SPACEDIM == 3 bool ambient_z_lo = (castro::ambient_fill_dir == 2 || castro::ambient_fill_dir == -1) && - (bcs(URHO).lo(2) == amrex::BCType::foextrap || bcs(URHO).lo(2) == amrex::BCType::hoextrap); + (bc.lo(2) == amrex::BCType::foextrap || bc.lo(2) == amrex::BCType::hoextrap); bool ambient_z_hi = (castro::ambient_fill_dir == 2 || castro::ambient_fill_dir == -1) && - (bcs(URHO).hi(2) == amrex::BCType::foextrap || bcs(URHO).hi(2) == amrex::BCType::hoextrap); + (bc.hi(2) == amrex::BCType::foextrap || bc.hi(2) == amrex::BCType::hoextrap); #endif if (castro::fill_ambient_bc == 1) { diff --git a/Source/radiation/Make.package b/Source/radiation/Make.package index 5a0023f4da..dfaa525b1b 100644 --- a/Source/radiation/Make.package +++ b/Source/radiation/Make.package @@ -5,7 +5,6 @@ CEXE_sources += HypreExtMultiABec.cpp CEXE_sources += HypreMultiABec.cpp CEXE_sources += HypreABec.cpp CEXE_sources += Radiation.cpp -CEXE_sources += radiation_params.cpp CEXE_sources += RadSolve.cpp CEXE_sources += RadBndry.cpp CEXE_sources += RadMultiGroup.cpp diff --git a/Source/radiation/Radiation.cpp b/Source/radiation/Radiation.cpp index 81a69d64bd..861c14a6d0 100644 --- a/Source/radiation/Radiation.cpp +++ b/Source/radiation/Radiation.cpp @@ -5,8 +5,6 @@ #include #include -#include - #include #include diff --git a/Source/radiation/radiation_params.cpp b/Source/radiation/radiation_params.cpp deleted file mode 100644 index 5c10888946..0000000000 --- a/Source/radiation/radiation_params.cpp +++ /dev/null @@ -1,9 +0,0 @@ -#include - -#include - -#include - -#include - -#include diff --git a/Source/reactions/Castro_react.H b/Source/reactions/Castro_react.H index 815e0a2774..76aa956761 100644 --- a/Source/reactions/Castro_react.H +++ b/Source/reactions/Castro_react.H @@ -1,3 +1,6 @@ +#ifndef CASTRO_REACT_H +#define CASTRO_REACT_H + /// /// Perform any old-time reactions. /// @@ -36,3 +39,5 @@ /// @param State State MultiFab /// static bool valid_zones_to_burn(amrex::MultiFab& State); + +#endif diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index 73c5a3c873..f226cd62e5 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -12,8 +12,8 @@ using namespace amrex; #ifndef TRUE_SDC advance_status -Castro::do_old_reactions (Real time, Real dt) -{ +Castro::do_old_reactions (Real time, Real dt) { // NOLINT(readability-convert-member-functions-to-static) + amrex::ignore_unused(time); amrex::ignore_unused(dt); @@ -203,7 +203,8 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra auto U = s.array(mfi); auto reactions = r.array(mfi); auto weights = store_burn_weights ? burn_weights.array(mfi) : Array4{}; - const auto mask = mask_covered_zones ? mask_mf.array(mfi) : Array4{}; + Array4 empty_arr{}; + const auto& mask = mask_covered_zones ? mask_mf.array(mfi) : empty_arr; const auto dx = geom.CellSizeArray(); #ifdef MODEL_PARSER @@ -213,7 +214,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) #else - LoopOnCpu(bx, [&] (int i, int j, int k) mutable + LoopOnCpu(bx, [&] (int i, int j, int k) #endif { @@ -533,7 +534,9 @@ Castro::react_state(Real time, Real dt) #endif int num_failed = 0; - // why no omp here? +#ifdef _OPENMP +#pragma omp parallel reduction(+:num_failed) +#endif for (MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); @@ -548,7 +551,8 @@ Castro::react_state(Real time, Real dt) auto I = SDC_react.array(mfi); auto react_src = reactions.array(mfi); auto weights = store_burn_weights ? burn_weights.array(mfi) : Array4{}; - const auto mask = mask_covered_zones ? mask_mf.array(mfi) : Array4{}; + Array4 empty_arr{}; + const auto& mask = mask_covered_zones ? mask_mf.array(mfi) : empty_arr; int lsdc_iteration = sdc_iteration; @@ -558,7 +562,7 @@ Castro::react_state(Real time, Real dt) #if defined(AMREX_USE_GPU) ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) #else - LoopOnCpu(bx, [&] (int i, int j, int k) mutable + LoopOnCpu(bx, [&] (int i, int j, int k) #endif { burn_t burn_state; diff --git a/Source/rotation/Castro_rotation.H b/Source/rotation/Castro_rotation.H index 0fec911450..af42388cc3 100644 --- a/Source/rotation/Castro_rotation.H +++ b/Source/rotation/Castro_rotation.H @@ -1,3 +1,5 @@ +#ifndef CASTRO_ROTATION_H +#define CASTRO_ROTATION_H /// /// Construct rotation source terms at old time @@ -80,3 +82,4 @@ fill_rotational_psi(const Box& bx, Array4 const& psi, const Real time); +#endif diff --git a/Source/sdc/Castro_sdc.cpp b/Source/sdc/Castro_sdc.cpp index 3172ef7a6b..00ba06c778 100644 --- a/Source/sdc/Castro_sdc.cpp +++ b/Source/sdc/Castro_sdc.cpp @@ -1,5 +1,4 @@ #include -// #include #include using namespace amrex; diff --git a/Source/sources/Castro_sources.H b/Source/sources/Castro_sources.H index 8dc34db4b3..2381d4e82d 100644 --- a/Source/sources/Castro_sources.H +++ b/Source/sources/Castro_sources.H @@ -1,3 +1,6 @@ +#ifndef CASTRO_SOURCES_H +#define CASTRO_SOURCES_H + /// /// Returns true if flag corresponding to source type ``src`` is set. /// @@ -54,7 +57,7 @@ /// @param dt the timestep to advance (e.g., go from time to /// time + dt) /// - advance_status do_old_sources (amrex::Real time, amrex::Real dt); + advance_status do_old_sources (amrex::Real time, amrex::Real dt, bool apply_to_state = true); /// @@ -146,7 +149,7 @@ /// /// @param update Real Vector of changes /// - static void print_source_change(amrex::Vector update); + static void print_source_change(const amrex::Vector& update); /// @@ -378,3 +381,5 @@ /// advance_status post_advance_operators (amrex::Real time, amrex::Real dt); + +#endif diff --git a/Source/sources/Castro_sources.cpp b/Source/sources/Castro_sources.cpp index f99aac52dc..9a0be7b695 100644 --- a/Source/sources/Castro_sources.cpp +++ b/Source/sources/Castro_sources.cpp @@ -151,7 +151,7 @@ Castro::do_old_sources( // Optionally print out diagnostic information about how much // these source terms changed the state. - if (print_update_diagnostics) { + if (apply_to_state && print_update_diagnostics) { bool is_new = false; print_all_source_changes(dt, is_new); } @@ -174,7 +174,7 @@ Castro::do_old_sources( } advance_status -Castro::do_old_sources (Real time, Real dt) +Castro::do_old_sources (Real time, Real dt, bool apply_to_state) { advance_status status {}; @@ -192,7 +192,8 @@ Castro::do_old_sources (Real time, Real dt) #ifdef MHD Bx_old, By_old, Bz_old, #endif - old_source, Sborder, S_new, time, dt); + old_source, Sborder, S_new, time, dt, + apply_to_state); return status; } @@ -459,7 +460,7 @@ Castro::evaluate_source_change(const MultiFab& source, Real dt, bool local) // interested in printing changes to energy, mass, etc. void -Castro::print_source_change(Vector update) +Castro::print_source_change(const Vector& update) { if (ParallelDescriptor::IOProcessor()) { @@ -522,13 +523,87 @@ Castro::print_all_source_changes(Real dt, bool is_new) // and the hydro advance. advance_status -Castro::pre_advance_operators (Real time, Real dt) +Castro::pre_advance_operators (Real time, Real dt) // NOLINT(readability-convert-member-functions-to-static) { amrex::ignore_unused(time); amrex::ignore_unused(dt); advance_status status {}; + // If we are using gravity, solve for the potential and + // gravitational field. note: since reactions don't change + // density, we can do this before or after the burn. + +#ifdef GRAVITY + construct_old_gravity(time); +#endif + +#ifdef SHOCK_VAR + // we want to compute the shock flag that will be used + // (optionally) in disabling reactions in shocks. We compute this + // only once per timestep using the time-level n data. + + // We need to compute the old sources -- note that we will + // recompute the old sources after the burn, so this is done here + // only for evaluating the shock flag. + + const bool apply_to_state{false}; + do_old_sources(time, dt, apply_to_state); + + MultiFab& S_old = get_old_data(State_Type); + MultiFab& old_source = get_old_data(Source_Type); + + FArrayBox shk(The_Async_Arena()); + FArrayBox q(The_Async_Arena()), qaux(The_Async_Arena()); + + for (MFIter mfi(Sborder, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + const Box& obx = mfi.growntilebox(1); + + shk.resize(bx, 1); +#ifdef RADIATION + q.resize(obx, NQ); +#else + q.resize(obx, NQTHERM); +#endif + qaux.resize(obx, NQAUX); + + Array4 const shk_arr = shk.array(); + Array4 const q_arr = q.array(); + Array4 const qaux_arr = qaux.array(); + + Array4 const Sborder_old_arr = Sborder.array(mfi); + Array4 const S_old_arr = S_old.array(mfi); + Array4 const old_src_arr = old_source.array(mfi); + + ctoprim(obx, time, Sborder_old_arr, q_arr, qaux_arr); + + shock(bx, q_arr, old_src_arr, shk_arr); + + // now store it in S_old -- we'll fillpatch into Sborder in a bit + + // Note: we still compute the shock flag in the hydro for the + // hybrid-Riemann (for now) since that version will have seen the + // effect of the burning in the first dt), but that version is + // never stored in State_Type + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + S_old_arr(i,j,k,USHK) = shk_arr(i,j,k); + }); + + } + + // we only computed the shock flag on the interior, but the first + // burn needs ghost cells, so FillPatch just the shock flag + + if (Sborder.nGrow() > 0) { + AmrLevel::FillPatch(*this, Sborder, Sborder.nGrow(), time, State_Type, USHK, 1, USHK); + } + +#endif + // If we are Strang splitting the reactions, do the old-time contribution now. #ifndef TRUE_SDC @@ -541,11 +616,6 @@ Castro::pre_advance_operators (Real time, Real dt) #endif #endif - // If we are using gravity, solve for the potential and gravitational field. - -#ifdef GRAVITY - construct_old_gravity(time); -#endif // Initialize the new-time data. This copy needs to come after all Strang-split operators. @@ -560,7 +630,7 @@ Castro::pre_advance_operators (Real time, Real dt) // but before the hydro advance. advance_status -Castro::pre_hydro_operators (Real time, Real dt) +Castro::pre_hydro_operators (Real time, Real dt) // NOLINT(readability-convert-member-functions-to-static) { amrex::ignore_unused(time); amrex::ignore_unused(dt); @@ -585,7 +655,7 @@ Castro::pre_hydro_operators (Real time, Real dt) // but before the corrector sources. advance_status -Castro::post_hydro_operators (Real time, Real dt) +Castro::post_hydro_operators (Real time, Real dt) // NOLINT(readability-convert-member-functions-to-static) { amrex::ignore_unused(time); amrex::ignore_unused(dt); @@ -602,8 +672,11 @@ Castro::post_hydro_operators (Real time, Real dt) // Perform all operations that occur after the corrector sources. advance_status -Castro::post_advance_operators (Real time, Real dt) +Castro::post_advance_operators (Real time, Real dt) // NOLINT(readability-convert-member-functions-to-static) { + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + advance_status status {}; #ifndef TRUE_SDC diff --git a/Util/scripts/write_probdata.py b/Util/scripts/write_probdata.py index 68cb716aac..2cd7cdd398 100755 --- a/Util/scripts/write_probdata.py +++ b/Util/scripts/write_probdata.py @@ -237,8 +237,8 @@ def write_probin(prob_param_files, cxx_prefix): size = p.size if (size == "nspec"): size = "NumSpec" - fout.write(f" for (int n = 0; n < {size}; n++) {{\n") - fout.write(f" problem::{p.name}[n] = {p.default_format(lang='C++')};\n") + fout.write(f" for (auto & e : problem::{p.name}) {{\n") + fout.write(f" e = {p.default_format(lang='C++')};\n") fout.write(f" }}\n") else: fout.write(f" {p.get_default_string()}") diff --git a/external/Microphysics b/external/Microphysics index 02d46a67cf..c89e33409e 160000 --- a/external/Microphysics +++ b/external/Microphysics @@ -1 +1 @@ -Subproject commit 02d46a67cf5852509cc35db7850e5f87da495679 +Subproject commit c89e33409e8c6bd35764dbb0f457a10e72a5af7b diff --git a/external/amrex b/external/amrex index d737886d57..76d09f554c 160000 --- a/external/amrex +++ b/external/amrex @@ -1 +1 @@ -Subproject commit d737886d574d4f1c0cf76323337b666ecd5bb4e0 +Subproject commit 76d09f554c3e5b8b4ecd4b8bf0c39d0fefddffa4