From 73245d6c1a60509b6e7ce0ec6b7cb8c83a39ad70 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 8 Dec 2023 12:41:53 -0500 Subject: [PATCH 1/9] fix spelling (#2672) --- Source/hydro/trace_ppm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index c86a5a38b3..dcddd09dee 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -288,7 +288,7 @@ Castro::trace_ppm(const Box& bx, ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho); } - // normal velcoity + // normal velocity Real Ip_src_un_0 = 0.0_rt; Real Im_src_un_0 = 0.0_rt; From 058d0f0cca05ac212b91a037aabba06ddac49bf2 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 8 Dec 2023 12:46:33 -0500 Subject: [PATCH 2/9] add a bit more flexibility to the massive star plotting script (#2671) This can now do a vertical orientation and also pass in the width for zooming via the commandline --- .../analysis/massive_star_multi.py | 60 +++++++++++++------ 1 file changed, 41 insertions(+), 19 deletions(-) diff --git a/Exec/science/massive_star/analysis/massive_star_multi.py b/Exec/science/massive_star/analysis/massive_star_multi.py index 08f4e9a25f..29f6bf3829 100755 --- a/Exec/science/massive_star/analysis/massive_star_multi.py +++ b/Exec/science/massive_star/analysis/massive_star_multi.py @@ -1,22 +1,24 @@ #!/usr/bin/env python3 +import argparse +import os + import matplotlib -matplotlib.use('agg') +import matplotlib.pyplot as plt -import os -import sys import yt -import matplotlib.pyplot as plt -import numpy as np -from functools import reduce +from yt.frontends.boxlib.api import CastroDataset +# assume that our data is in CGS +from yt.units import cm + +matplotlib.use('agg') -from mpl_toolkits.axes_grid1 import ImageGrid -# assume that our data is in CGS -from yt.units import cm, amu -from yt.frontends.boxlib.api import CastroDataset -def make_plot(plotfile, prefix="plot", size=(19.2, 10.8), cbar_location="right"): + +def make_plot(plotfile, prefix="plot", width_frac=0.1, + layout=(1, 4), + size=(19.2, 10.8), cbar_location="right"): ds = CastroDataset(plotfile) @@ -36,8 +38,6 @@ def make_plot(plotfile, prefix="plot", size=(19.2, 10.8), cbar_location="right") fig = plt.figure() fig.set_size_inches(size) - width_frac = 0.1 - fields = ["MachNumber", "magvort", "abar", "enuc"] sp = yt.SlicePlot(ds, "theta", fields, @@ -75,17 +75,39 @@ def make_plot(plotfile, prefix="plot", size=(19.2, 10.8), cbar_location="right") sp.set_axes_unit("cm") - fig = sp.export_to_mpl_figure((1, len(fields)), cbar_location=cbar_location, cbar_pad="5%") + fig = sp.export_to_mpl_figure((layout[0], layout[1]), axes_pad=(1.0, 0.4), + cbar_location=cbar_location, cbar_pad="2%") fig.subplots_adjust(left=0.05, right=0.95, bottom=0.025, top=0.975) - fig.text(0.02, 0.02, "time = {:8.5f} s".format(float(ds.current_time)), transform=fig.transFigure) + fig.text(0.02, 0.02, f"time = {float(ds.current_time):8.2f} s", + transform=fig.transFigure) fig.set_size_inches(size) + fig.tight_layout() + extra = "" + if layout[0] >= layout[1]: + extra = "_vertical" - fig.savefig(f"{prefix}_{os.path.basename(plotfile)}_slice.png", pad_inches=0) + fig.savefig(f"{prefix}_{os.path.basename(plotfile)}_w{width_frac:04.2f}{extra}.pdf", pad_inches=0) if __name__ == "__main__": - plotfile = sys.argv[1] - - make_plot(plotfile, "all") + p = argparse.ArgumentParser() + p.add_argument("--vertical", action="store_true", + help="plot 2x2 or 1x4") + p.add_argument("--width_fraction", type=float, default=0.1, + help="fraction of domain to show") + p.add_argument("plotfile", type=str, nargs=1, + help="plotfile to plot") + + args = p.parse_args() + plotfile = args.plotfile[0] + + if args.vertical: + size = (7.5, 11.0) + layout = (2, 2) + else: + size = (19.2, 8.5) + layout = (1, 4) + + make_plot(plotfile, "all", layout=layout, width_frac=args.width_fraction, size=size) From 9cf0cb6b3e31a382e73dac7d94875c2e099b1e0b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 8 Dec 2023 13:16:51 -0500 Subject: [PATCH 3/9] sync to the new tabular NSE interfacet (#2673) this syncs up with the Microphysics changes in AMReX-Astro/Microphysics#1404 --- .../nse_test/problem_initialize_state_data.H | 26 ++++++++++--------- .../problem_initialize_state_data.H | 22 +++++++--------- 2 files changed, 23 insertions(+), 25 deletions(-) diff --git a/Exec/reacting_tests/nse_test/problem_initialize_state_data.H b/Exec/reacting_tests/nse_test/problem_initialize_state_data.H index 43b6666355..884a606411 100644 --- a/Exec/reacting_tests/nse_test/problem_initialize_state_data.H +++ b/Exec/reacting_tests/nse_test/problem_initialize_state_data.H @@ -5,6 +5,7 @@ #include #include #if defined(NSE_TABLE) +#include #include #elif defined(NSE_NET) #include @@ -64,14 +65,15 @@ void problem_initialize_state_data (int i, int j, int k, burn_state.y_e = ye; #if defined(NSE_TABLE) - Real abar; - Real dq; - Real dyedt; - Real dabardt; - Real dbeadt; - Real e_nu; - - nse_interp(T, problem::rho0, ye, abar, dq, dyedt, dabardt, dbeadt, e_nu, burn_state.xn); + nse_table_t nse_state; + nse_state.T = T; + nse_state.rho = problem::rho0; + nse_state.Ye = ye; + nse_interp(nse_state); + + for (int n = 0; n < NumSpec; ++n) { + burn_state.xn[n] = nse_state.X[n]; + } #elif defined(NSE_NET) Real eps = 1.e-10_rt; bool input_ye_is_valid = true; @@ -105,16 +107,16 @@ void problem_initialize_state_data (int i, int j, int k, for (int n = 0; n < NumSpec; n++) { burn_state.xn[n] /= sumX; } - + #ifdef NSE_NET state(i,j,k,UMUP) = burn_state.mu_p; state(i,j,k,UMUN) = burn_state.mu_n; #endif - + #ifdef AUX_THERMO burn_state.aux[AuxZero::iye] = ye; - burn_state.aux[AuxZero::iabar] = abar; - burn_state.aux[AuxZero::ibea] = dq; + burn_state.aux[AuxZero::iabar] = nse_state.abar; + burn_state.aux[AuxZero::ibea] = nse_state.bea; #endif eos(eos_input_rt, burn_state); diff --git a/Exec/science/massive_star/problem_initialize_state_data.H b/Exec/science/massive_star/problem_initialize_state_data.H index 5850ba23b4..c34101fd76 100644 --- a/Exec/science/massive_star/problem_initialize_state_data.H +++ b/Exec/science/massive_star/problem_initialize_state_data.H @@ -5,6 +5,7 @@ #include #include #ifdef NSE_TABLE +#include #include #endif @@ -80,21 +81,16 @@ void problem_initialize_state_data (int i, int j, int k, if (nse_check) { - Real abar; - Real dq; - Real dyedt; - Real dabardt; - Real dbeadt; - Real e_nu; - Real xn[NumSpec]; + nse_table_t nse_state; + nse_state.T = state(i,j,k,UTEMP); + nse_state.rho = state(i,j,k,URHO); + nse_state.Ye = state(i,j,k,UFX+AuxZero::iye); + nse_interp(nse_state); - nse_interp(state(i,j,k,UTEMP), state(i,j,k,URHO), - state(i,j,k,UFX+AuxZero::iye), abar, dq, dyedt, dabardt, dbeadt, e_nu, xn); - - state(i,j,k,UFX+AuxZero::iabar) = abar; - state(i,j,k,UFX+AuxZero::ibea) = dq; + state(i,j,k,UFX+AuxZero::iabar) = nse_state.abar; + state(i,j,k,UFX+AuxZero::ibea) = nse_state.bea; for (int n = 0; n < NumSpec; n++) { - state(i,j,k,UFS+n) = amrex::max(xn[n], small_x); + state(i,j,k,UFS+n) = amrex::max(nse_state.X[n], small_x); } // we just got new X's, so we need to re-renormalize them From 72789c07ab2949088ee327e73a2cc4c944c649ef Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 13:59:38 -0500 Subject: [PATCH 4/9] sync up with amrex hip script for ROCm 6.0 (#2680) --- .github/workflows/dependencies_hip.sh | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/.github/workflows/dependencies_hip.sh b/.github/workflows/dependencies_hip.sh index e5fb425b67..36df2f384b 100755 --- a/.github/workflows/dependencies_hip.sh +++ b/.github/workflows/dependencies_hip.sh @@ -1,17 +1,26 @@ #!/usr/bin/env bash # -# Copyright 2020-2022 The AMReX Community +# Copyright 2020 The AMReX Community # # License: BSD-3-Clause-LBNL # Authors: Axel Huebl +# search recursive inside a folder if a file contains tabs +# +# @result 0 if no files are found, else 1 +# + 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 # Ref.: https://rocmdocs.amd.com/en/latest/Installation_Guide/Installation-Guide.html#ubuntu curl -O https://repo.radeon.com/rocm/rocm.gpg.key sudo apt-key add rocm.gpg.key -echo 'deb [arch=amd64] https://repo.radeon.com/rocm/apt/debian/ ubuntu main' \ +echo "deb [arch=amd64] https://repo.radeon.com/rocm/apt/${1-debian}/ ubuntu main" \ | sudo tee /etc/apt/sources.list.d/rocm.list echo 'export PATH=/opt/rocm/llvm/bin:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin:$PATH' \ | sudo tee -a /etc/profile.d/rocm.sh @@ -34,12 +43,14 @@ sudo apt-get install -y --no-install-recommends \ roctracer-dev \ rocprofiler-dev \ rocrand-dev \ - rocprim-dev + rocprim-dev \ + hiprand-dev # activate # source /etc/profile.d/rocm.sh hipcc --version +hipconfig --full which clang which clang++ which flang From 8ee041fdca954c27fdfe5a793ccc18c8dc95c148 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 14:19:23 -0500 Subject: [PATCH 5/9] allow for optional passive sources to the conserved state (#2678) Space will be allocated for the passive scalars in the conserved state if we set USE_SPECIES_SOURCES=TRUE in the makefile. This is off by default, as most applications don't need it. --- Exec/Make.Castro | 4 ++++ Source/driver/_variables | 6 +++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/Exec/Make.Castro b/Exec/Make.Castro index cfd7e187a9..79e747b5c8 100644 --- a/Exec/Make.Castro +++ b/Exec/Make.Castro @@ -270,6 +270,10 @@ ifeq ($(USE_ROTATION), TRUE) DEFINES += -DROTATION endif +ifeq ($(USE_SPECIES_SOURCES), TRUE) + DEFINES += -DCONS_SPECIES_HAVE_SOURCES +endif + ifeq ($(USE_PARTICLES), TRUE) Bdirs += Source/particles endif diff --git a/Source/driver/_variables b/Source/driver/_variables index 36d05dfc06..b6408c82a8 100644 --- a/Source/driver/_variables +++ b/Source/driver/_variables @@ -16,9 +16,9 @@ energy-density UEDEN NSRC 1 None internal-energy UEINT NSRC 1 None temperature UTEMP NSRC 1 None - advected UFA None NumAdv None - species UFS None NumSpec None - auxiliary UFX None NumAux None + advected UFA [(NSRC, CONS_SPECIES_HAVE_SOURCES)] NumAdv None + species UFS [(NSRC, CONS_SPECIES_HAVE_SOURCES)] NumSpec None + auxiliary UFX [(NSRC, CONS_SPECIES_HAVE_SOURCES)] NumAux None shock USHK None 1 SHOCK_VAR mu_p UMUP None 1 NSE_NET mu_n UMUN None 1 NSE_NET From 5731a8fdd2380eaa02a6ed781851cd85c08236d9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 14:19:42 -0500 Subject: [PATCH 6/9] remove elixirs in Source/driver/ (#2677) --- Source/driver/Castro_advance_sdc.cpp | 9 +++------ Source/driver/Derive.cpp | 3 +-- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/Source/driver/Castro_advance_sdc.cpp b/Source/driver/Castro_advance_sdc.cpp index e6bec0804c..2d0463aa1f 100644 --- a/Source/driver/Castro_advance_sdc.cpp +++ b/Source/driver/Castro_advance_sdc.cpp @@ -122,8 +122,7 @@ Castro::do_advance_sdc (Real time, for (MFIter mfi(S_new); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - tmp.resize(bx, 1); - Elixir elix_tmp = tmp.elixir(); + tmp.resize(bx, 1, The_Async_Arena()); auto tmp_arr = tmp.array(); make_fourth_in_place(bx, old_source.array(mfi), tmp_arr, domain_lo, domain_hi); @@ -305,8 +304,7 @@ Castro::do_advance_sdc (Real time, // pass in the reaction source at centers (Sburn_arr), including // one ghost cell and derive everything that is needed including // 1 ghost cell - R_center.resize(obx, R_new.nComp()); - Elixir elix_r_center = R_center.elixir(); + R_center.resize(obx, R_new.nComp(), The_Async_Arena()); auto const R_center_arr = R_center.array(); Array4 const Sburn_arr = Sburn.array(mfi); @@ -315,8 +313,7 @@ Castro::do_advance_sdc (Real time, ca_store_reaction_state(obx, Sburn_arr, R_center_arr); // convert R_new from centers to averages in place - tmp.resize(bx, 1); - Elixir elix_tmp = tmp.elixir(); + tmp.resize(bx, 1, The_Async_Arena()); auto const tmp_arr = tmp.array(); make_fourth_in_place(bx, R_center_arr, tmp_arr, domain_lo, domain_hi); diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index 6a6f6207ad..3064fc3dd2 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -352,8 +352,7 @@ extern "C" const Box& obx = amrex::grow(bx, 1); FArrayBox coeff_cc; - coeff_cc.resize(obx, 1); - Elixir elix_coeff_cc = coeff_cc.elixir(); + coeff_cc.resize(obx, 1, The_Async_Arena()); Array4 const coeff_arr = coeff_cc.array(); auto const dat = datfab.array(); From 27d602e534f67a4dcfca6837a01537da843245fd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 14:20:05 -0500 Subject: [PATCH 7/9] remove elixirs in Castro_ctu_hydro (#2675) --- Source/hydro/Castro_ctu_hydro.cpp | 103 ++++++++++-------------------- 1 file changed, 35 insertions(+), 68 deletions(-) diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index a3e5b28a9c..985527f437 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -135,45 +135,54 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) int priv_nstep_fsp = -1; #endif - // Declare local storage now. This should be done outside the MFIter loop, - // and then we will resize the Fabs in each MFIter loop iteration. Then, - // we apply an Elixir to ensure that their memory is saved until it is no - // longer needed (only relevant for the asynchronous case, usually on GPUs). - - FArrayBox shk; - FArrayBox q, qaux; - FArrayBox rho_inv; - FArrayBox src_q; - FArrayBox qxm, qxp; + // Declare local storage now. This should be done outside the + // MFIter loop, and then we will resize the Fabs in each MFIter + // loop iteration. We use the async arenato ensure that their + // memory is saved until it is no longer needed (only relevant for + // the asynchronous case, usually on GPUs). + + FArrayBox shk(The_Async_Arena()); + FArrayBox q(The_Async_Arena()), qaux(The_Async_Arena()); + FArrayBox rho_inv(The_Async_Arena()); + FArrayBox src_q(The_Async_Arena()); + FArrayBox qxm(The_Async_Arena()), qxp(The_Async_Arena()); #if AMREX_SPACEDIM >= 2 - FArrayBox qym, qyp; + FArrayBox qym(The_Async_Arena()), qyp(The_Async_Arena()); #endif #if AMREX_SPACEDIM == 3 - FArrayBox qzm, qzp; + FArrayBox qzm(The_Async_Arena()), qzp(The_Async_Arena()); #endif - FArrayBox div; + FArrayBox div(The_Async_Arena()); #if AMREX_SPACEDIM >= 2 - FArrayBox ftmp1, ftmp2; + FArrayBox ftmp1(The_Async_Arena()), ftmp2(The_Async_Arena()); #ifdef RADIATION - FArrayBox rftmp1, rftmp2; + FArrayBox rftmp1(The_Async_Arena()), rftmp2(The_Async_Arena()); #endif - FArrayBox qgdnvtmp1, qgdnvtmp2; - FArrayBox ql, qr; + FArrayBox qgdnvtmp1(The_Async_Arena()), qgdnvtmp2(The_Async_Arena()); + FArrayBox ql(The_Async_Arena()), qr(The_Async_Arena()); #endif - FArrayBox flux[AMREX_SPACEDIM], qe[AMREX_SPACEDIM]; + Vector flux, qe; + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + flux.push_back(FArrayBox(The_Async_Arena())); + qe.push_back(FArrayBox(The_Async_Arena())); + } + #ifdef RADIATION - FArrayBox rad_flux[AMREX_SPACEDIM]; + Vector rad_flux; + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + rad_flux.push_back(FArrayBox(The_Async_Arena())); + } #endif #if AMREX_SPACEDIM <= 2 - FArrayBox pradial; + FArrayBox pradial(The_Async_Arena()); #endif #if AMREX_SPACEDIM == 3 - FArrayBox qmyx, qpyx; - FArrayBox qmzx, qpzx; - FArrayBox qmxy, qpxy; - FArrayBox qmzy, qpzy; - FArrayBox qmxz, qpxz; - FArrayBox qmyz, qpyz; + FArrayBox qmyx(The_Async_Arena()), qpyx(The_Async_Arena()); + FArrayBox qmzx(The_Async_Arena()), qpzx(The_Async_Arena()); + FArrayBox qmxy(The_Async_Arena()), qpxy(The_Async_Arena()); + FArrayBox qmzy(The_Async_Arena()), qpzy(The_Async_Arena()); + FArrayBox qmxz(The_Async_Arena()), qpxz(The_Async_Arena()); + FArrayBox qmyz(The_Async_Arena()), qpyz(The_Async_Arena()); #endif MultiFab& old_source = get_old_data(Source_Type); @@ -198,19 +207,16 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // primitive versions on demand as needed q.resize(qbx, NQTHERM); #endif - Elixir elix_q = q.elixir(); fab_size += q.nBytes(); Array4 const q_arr = q.array(); qaux.resize(qbx, NQAUX); - Elixir elix_qaux = qaux.elixir(); fab_size += qaux.nBytes(); Array4 const qaux_arr = qaux.array(); Array4 const U_old_arr = Sborder.array(mfi); rho_inv.resize(qbx3, 1); - Elixir elix_rho_inv = rho_inv.elixir(); fab_size += rho_inv.nBytes(); Array4 const rho_inv_arr = rho_inv.array(); @@ -248,7 +254,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #endif shk.resize(obx, 1); - Elixir elix_shk = shk.elixir(); fab_size += shk.nBytes(); Array4 const shk_arr = shk.array(); @@ -276,7 +281,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // get the primitive variable hydro sources src_q.resize(qbx3, NQSRC); - Elixir elix_src_q = src_q.elixir(); fab_size += src_q.nBytes(); Array4 const src_q_arr = src_q.array(); @@ -292,11 +296,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // work on the interface states qxm.resize(obx, NQ); - Elixir elix_qxm = qxm.elixir(); fab_size += qxm.nBytes(); qxp.resize(obx, NQ); - Elixir elix_qxp = qxp.elixir(); fab_size += qxp.nBytes(); Array4 const qxm_arr = qxm.array(); @@ -304,11 +306,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #if AMREX_SPACEDIM >= 2 qym.resize(obx, NQ); - Elixir elix_qym = qym.elixir(); fab_size += qym.nBytes(); qyp.resize(obx, NQ); - Elixir elix_qyp = qyp.elixir(); fab_size += qyp.nBytes(); Array4 const qym_arr = qym.array(); @@ -318,11 +318,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #if AMREX_SPACEDIM == 3 qzm.resize(obx, NQ); - Elixir elix_qzm = qzm.elixir(); fab_size += qzm.nBytes(); qzp.resize(obx, NQ); - Elixir elix_qzp = qzp.elixir(); fab_size += qzp.nBytes(); Array4 const qzm_arr = qzm.array(); @@ -387,7 +385,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) } div.resize(obx, 1); - Elixir elix_div = div.elixir(); fab_size += div.nBytes(); auto div_arr = div.array(); @@ -395,36 +392,30 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) divu(obx, q_arr, div_arr); flux[0].resize(gxbx, NUM_STATE); - Elixir elix_flux_x = flux[0].elixir(); fab_size += flux[0].nBytes(); Array4 const flux0_arr = (flux[0]).array(); qe[0].resize(gxbx, NGDNV); - Elixir elix_qe_x = qe[0].elixir(); auto qex_arr = qe[0].array(); fab_size += qe[0].nBytes(); #ifdef RADIATION rad_flux[0].resize(gxbx, Radiation::nGroups); - Elixir elix_rad_flux_x = rad_flux[0].elixir(); fab_size += rad_flux[0].nBytes(); auto rad_flux0_arr = (rad_flux[0]).array(); #endif #if AMREX_SPACEDIM >= 2 flux[1].resize(gybx, NUM_STATE); - Elixir elix_flux_y = flux[1].elixir(); fab_size += flux[1].nBytes(); Array4 const flux1_arr = (flux[1]).array(); qe[1].resize(gybx, NGDNV); - Elixir elix_qe_y = qe[1].elixir(); auto qey_arr = qe[1].array(); fab_size += qe[1].nBytes(); #ifdef RADIATION rad_flux[1].resize(gybx, Radiation::nGroups); - Elixir elix_rad_flux_y = rad_flux[1].elixir(); fab_size += rad_flux[1].nBytes(); auto const rad_flux1_arr = (rad_flux[1]).array(); #endif @@ -432,18 +423,15 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #if AMREX_SPACEDIM == 3 flux[2].resize(gzbx, NUM_STATE); - Elixir elix_flux_z = flux[2].elixir(); fab_size += flux[2].nBytes(); Array4 const flux2_arr = (flux[2]).array(); qe[2].resize(gzbx, NGDNV); - Elixir elix_qe_z = qe[2].elixir(); auto qez_arr = qe[2].array(); fab_size += qe[2].nBytes(); #ifdef RADIATION rad_flux[2].resize(gzbx, Radiation::nGroups); - Elixir elix_rad_flux_z = rad_flux[2].elixir(); fab_size += rad_flux[2].nBytes(); auto const rad_flux2_arr = (rad_flux[2]).array(); #endif @@ -453,7 +441,6 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) if (!Geom().IsCartesian()) { pradial.resize(xbx, 1); } - Elixir elix_pradial = pradial.elixir(); fab_size += pradial.nBytes(); #endif @@ -491,46 +478,38 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #if AMREX_SPACEDIM >= 2 ftmp1.resize(obx, NUM_STATE); - Elixir elix_ftmp1 = ftmp1.elixir(); auto ftmp1_arr = ftmp1.array(); fab_size += ftmp1.nBytes(); ftmp2.resize(obx, NUM_STATE); - Elixir elix_ftmp2 = ftmp2.elixir(); auto ftmp2_arr = ftmp2.array(); fab_size += ftmp2.nBytes(); #ifdef RADIATION rftmp1.resize(obx, Radiation::nGroups); - Elixir elix_rftmp1 = rftmp1.elixir(); auto rftmp1_arr = rftmp1.array(); fab_size += rftmp1.nBytes(); rftmp2.resize(obx, Radiation::nGroups); - Elixir elix_rftmp2 = rftmp2.elixir(); auto rftmp2_arr = rftmp2.array(); fab_size += rftmp2.nBytes(); #endif qgdnvtmp1.resize(obx, NGDNV); - Elixir elix_qgdnvtmp1 = qgdnvtmp1.elixir(); auto qgdnvtmp1_arr = qgdnvtmp1.array(); fab_size += qgdnvtmp1.nBytes(); #if AMREX_SPACEDIM == 3 qgdnvtmp2.resize(obx, NGDNV); - Elixir elix_qgdnvtmp2 = qgdnvtmp2.elixir(); auto qgdnvtmp2_arr = qgdnvtmp2.array(); fab_size += qgdnvtmp2.nBytes(); #endif ql.resize(obx, NQ); - Elixir elix_ql = ql.elixir(); auto ql_arr = ql.array(); fab_size += ql.nBytes(); qr.resize(obx, NQ); - Elixir elix_qr = qr.elixir(); auto qr_arr = qr.array(); fab_size += qr.nBytes(); #endif @@ -699,12 +678,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) const Box& tyxbx = amrex::grow(ybx, IntVect(AMREX_D_DECL(0,0,1))); qmyx.resize(tyxbx, NQ); - Elixir elix_qmyx = qmyx.elixir(); auto qmyx_arr = qmyx.array(); fab_size += qmyx.nBytes(); qpyx.resize(tyxbx, NQ); - Elixir elix_qpyx = qpyx.elixir(); auto qpyx_arr = qpyx.array(); fab_size += qpyx.nBytes(); @@ -730,12 +707,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) const Box& tzxbx = amrex::grow(zbx, IntVect(AMREX_D_DECL(0,1,0))); qmzx.resize(tzxbx, NQ); - Elixir elix_qmzx = qmzx.elixir(); auto qmzx_arr = qmzx.array(); fab_size += qmzx.nBytes(); qpzx.resize(tzxbx, NQ); - Elixir elix_qpzx = qpzx.elixir(); auto qpzx_arr = qpzx.array(); fab_size += qpzx.nBytes(); @@ -775,12 +750,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) const Box& txybx = amrex::grow(xbx, IntVect(AMREX_D_DECL(0,0,1))); qmxy.resize(txybx, NQ); - Elixir elix_qmxy = qmxy.elixir(); auto qmxy_arr = qmxy.array(); fab_size += qmxy.nBytes(); qpxy.resize(txybx, NQ); - Elixir elix_qpxy = qpxy.elixir(); auto qpxy_arr = qpxy.array(); fab_size += qpxy.nBytes(); @@ -806,12 +779,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) const Box& tzybx = amrex::grow(zbx, IntVect(AMREX_D_DECL(1,0,0))); qmzy.resize(tzybx, NQ); - Elixir elix_qmzy = qmzy.elixir(); auto qmzy_arr = qmzy.array(); fab_size += qmzy.nBytes(); qpzy.resize(tzybx, NQ); - Elixir elix_qpzy = qpzy.elixir(); auto qpzy_arr = qpzy.array(); fab_size += qpzy.nBytes(); @@ -854,12 +825,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) const Box& txzbx = amrex::grow(xbx, IntVect(AMREX_D_DECL(0,1,0))); qmxz.resize(txzbx, NQ); - Elixir elix_qmxz = qmxz.elixir(); auto qmxz_arr = qmxz.array(); fab_size += qmxz.nBytes(); qpxz.resize(txzbx, NQ); - Elixir elix_qpxz = qpxz.elixir(); auto qpxz_arr = qpxz.array(); fab_size += qpxz.nBytes(); @@ -885,12 +854,10 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) const Box& tyzbx = amrex::grow(ybx, IntVect(AMREX_D_DECL(1,0,0))); qmyz.resize(tyzbx, NQ); - Elixir elix_qmyz = qmyz.elixir(); auto qmyz_arr = qmyz.array(); fab_size += qmyz.nBytes(); qpyz.resize(tyzbx, NQ); - Elixir elix_qpyz = qpyz.elixir(); auto qpyz_arr = qpyz.array(); fab_size += qpyz.nBytes(); From 1613d3dbefa0b1dd861aa421bb958418dd3dda08 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 14:20:24 -0500 Subject: [PATCH 8/9] remove Elixir in Castro.cpp (#2674) we now use The_Async_Arena() --- Source/driver/Castro.cpp | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 7721244b93..8c6724d621 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -1260,8 +1260,7 @@ Castro::initData () { const Box& box = mfi.validbox(); - tmp.resize(box, 1); - Elixir elix_tmp = tmp.elixir(); + tmp.resize(box, 1, The_Async_Arena()); auto tmp_arr = tmp.array(); make_fourth_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi); @@ -1287,8 +1286,7 @@ Castro::initData () { const Box& box = mfi.growntilebox(2); - tmp.resize(box, 1); - Elixir elix_tmp = tmp.elixir(); + tmp.resize(box, 1, The_Async_Arena()); auto tmp_arr = tmp.array(); make_cell_center_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi); @@ -1336,8 +1334,7 @@ Castro::initData () { const Box& box = mfi.validbox(); - tmp.resize(box, 1); - Elixir elix_tmp = tmp.elixir(); + tmp.resize(box, 1, The_Async_Arena()); auto tmp_arr = tmp.array(); make_fourth_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi); @@ -3950,8 +3947,7 @@ Castro::computeTemp( compute_lap_term(bx0, Stemp.array(mfi), Eint_lap.array(mfi), UEINT, domain_lo, domain_hi); - tmp.resize(bx, 1); - Elixir elix_tmp = tmp.elixir(); + tmp.resize(bx, 1, The_Async_Arena()); auto tmp_arr = tmp.array(); make_cell_center_in_place(bx, Stemp.array(mfi), tmp_arr, domain_lo, domain_hi); @@ -4067,8 +4063,7 @@ Castro::computeTemp( const Box& bx = mfi.tilebox(); - tmp.resize(bx, 1); - Elixir elix_tmp = tmp.elixir(); + tmp.resize(bx, 1, The_Async_Arena()); auto tmp_arr = tmp.array(); // only temperature From e65b3b9c2b50c38d0c2ed381516a24e6867df9cd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 14:21:24 -0500 Subject: [PATCH 9/9] add a no-limiting option to PPM (#2670) This is useful for convergence testing. --- Source/driver/_cpp_parameters | 4 +- Source/hydro/ppm.H | 117 +++++++++++++++++++--------------- 2 files changed, 67 insertions(+), 54 deletions(-) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 774d6916ac..4d1dde2581 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -84,9 +84,11 @@ hybrid_hydro int 0 # reconstruction type: # 0: piecewise linear; # 1: classic Colella \& Woodward ppm; -# 2: extrema-preserving ppm ppm_type int 1 +# do we limit the ppm parabola? +ppm_do_limiting int 1 + # For MHD + PLM, do we limit on characteristic or primitive variables mhd_limit_characteristic int 1 diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 9070a48bc1..1fe14199fb 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -57,84 +57,95 @@ ppm_reconstruct(const Real* s, const Real flatn, Real& sm, Real& sp) { - // Compute van Leer slopes + if (ppm_do_limiting) { - Real dsl = 2.0_rt * (s[im1] - s[im2]); - Real dsr = 2.0_rt * (s[i0] - s[im1]); + // Compute van Leer slopes - Real dsvl_l = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[i0] - s[im2]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); - } + Real dsl = 2.0_rt * (s[im1] - s[im2]); + Real dsr = 2.0_rt * (s[i0] - s[im1]); - dsl = 2.0_rt * (s[i0] - s[im1]); - dsr = 2.0_rt * (s[ip1] - s[i0]); + Real dsvl_l = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[i0] - s[im2]); + dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); + } - Real dsvl_r = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); - } + dsl = 2.0_rt * (s[i0] - s[im1]); + dsr = 2.0_rt * (s[ip1] - s[i0]); - // Interpolate s to edges + Real dsvl_r = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[ip1] - s[im1]); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); + } - sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + // Interpolate s to edges - // Make sure sedge lies in between adjacent cell-centered values + sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); - sm = amrex::max(sm, amrex::min(s[i0], s[im1])); - sm = amrex::min(sm, amrex::max(s[i0], s[im1])); + // Make sure sedge lies in between adjacent cell-centered values + sm = amrex::max(sm, amrex::min(s[i0], s[im1])); + sm = amrex::min(sm, amrex::max(s[i0], s[im1])); - // Compute van Leer slopes - dsl = 2.0_rt * (s[i0] - s[im1]); - dsr = 2.0_rt * (s[ip1] - s[i0]); + // Compute van Leer slopes - dsvl_l = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr))); - } + dsl = 2.0_rt * (s[i0] - s[im1]); + dsr = 2.0_rt * (s[ip1] - s[i0]); - dsl = 2.0_rt * (s[ip1] - s[i0]); - dsr = 2.0_rt * (s[ip2] - s[ip1]); + dsvl_l = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[ip1] - s[im1]); + dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr))); + } - dsvl_r = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip2] - s[i0]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr))); - } + dsl = 2.0_rt * (s[ip1] - s[i0]); + dsr = 2.0_rt * (s[ip2] - s[ip1]); - // Interpolate s to edges + dsvl_r = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[ip2] - s[i0]); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), amrex::min(std::abs(dsl),std::abs(dsr))); + } - sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + // Interpolate s to edges - // Make sure sedge lies in between adjacent cell-centered values + sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); - sp = amrex::max(sp, amrex::min(s[ip1], s[i0])); - sp = amrex::min(sp, amrex::max(s[ip1], s[i0])); + // Make sure sedge lies in between adjacent cell-centered values + sp = amrex::max(sp, amrex::min(s[ip1], s[i0])); + sp = amrex::min(sp, amrex::max(s[ip1], s[i0])); - // Flatten the parabola - sm = flatn * sm + (1.0_rt - flatn) * s[i0]; - sp = flatn * sp + (1.0_rt - flatn) * s[i0]; + // Flatten the parabola - // Modify using quadratic limiters -- note this version of the limiting comes - // from Colella and Sekora (2008), not the original PPM paper. + sm = flatn * sm + (1.0_rt - flatn) * s[i0]; + sp = flatn * sp + (1.0_rt - flatn) * s[i0]; - if ((sp - s[i0]) * (s[i0] - sm) <= 0.0_rt) { - sp = s[i0]; - sm = s[i0]; + // Modify using quadratic limiters -- note this version of the limiting comes + // from Colella and Sekora (2008), not the original PPM paper. - } else if (std::abs(sp - s[i0]) >= 2.0_rt * std::abs(sm - s[i0])) { - sp = 3.0_rt * s[i0] - 2.0_rt * sm; + if ((sp - s[i0]) * (s[i0] - sm) <= 0.0_rt) { + sp = s[i0]; + sm = s[i0]; - } else if (std::abs(sm - s[i0]) >= 2.0_rt * std::abs(sp - s[i0])) { - sm = 3.0_rt * s[i0] - 2.0_rt * sp; - } + } else if (std::abs(sp - s[i0]) >= 2.0_rt * std::abs(sm - s[i0])) { + sp = 3.0_rt * s[i0] - 2.0_rt * sm; + + } else if (std::abs(sm - s[i0]) >= 2.0_rt * std::abs(sp - s[i0])) { + sm = 3.0_rt * s[i0] - 2.0_rt * sp; + } + + } else { + + // unlimited PPM reconstruction (Eq. 1.9 in the PPM paper) + + sm = (7.0_rt/12.0_rt) * (s[i0] + s[im1]) - (1.0_rt/12.0_rt) * (s[im2] + s[ip1]); + sp = (7.0_rt/12.0_rt) * (s[ip1] + s[i0]) - (1.0_rt/12.0_rt) * (s[im1] + s[ip2]); + + } }