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 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/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 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(); 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/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 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(); 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]); + + } }