diff --git a/Exec/mhd_tests/LoopAdvection/problem_initialize.H b/Exec/mhd_tests/LoopAdvection/problem_initialize.H index 350fe3be27..fb229657cb 100644 --- a/Exec/mhd_tests/LoopAdvection/problem_initialize.H +++ b/Exec/mhd_tests/LoopAdvection/problem_initialize.H @@ -16,8 +16,8 @@ void problem_initialize () eos_state.rho = problem::rho_0; eos_state.p = problem::p_0; eos_state.T = 100000.0_rt; // initial guess - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = 0.0_rt; + for (auto & X : eos_state.xn) { + X = 0.0_rt; } eos_state.xn[0] = 1.0_rt; diff --git a/Exec/mhd_tests/LoopAdvection/problem_initialize_mhd_data.H b/Exec/mhd_tests/LoopAdvection/problem_initialize_mhd_data.H index 3241f0ab28..7c9716f44c 100644 --- a/Exec/mhd_tests/LoopAdvection/problem_initialize_mhd_data.H +++ b/Exec/mhd_tests/LoopAdvection/problem_initialize_mhd_data.H @@ -7,6 +7,8 @@ Real A_z (int i, int j, int k, const GeometryData& geomdata) { + amrex::ignore_unused(k); + // Compute A_z. This lives on edges, e.g., {A_z}_{i-1/2,j-1/2,k} // so it is centered only in the z-direction. diff --git a/Exec/mhd_tests/LoopAdvection/problem_initialize_state_data.H b/Exec/mhd_tests/LoopAdvection/problem_initialize_state_data.H index d3a6a25391..b7eab9ee6a 100644 --- a/Exec/mhd_tests/LoopAdvection/problem_initialize_state_data.H +++ b/Exec/mhd_tests/LoopAdvection/problem_initialize_state_data.H @@ -10,6 +10,8 @@ void problem_initialize_state_data (int i, int j, int k, const GeometryData& geomdata) { + amrex::ignore_unused(geomdata); + state(i,j,k,URHO) = problem::rho_0; state(i,j,k,UMX) = problem::rho_0 * problem::u_x; state(i,j,k,UMY) = problem::rho_0 * problem::u_y; diff --git a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc index 6316753695..9b49a701ed 100644 --- a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc +++ b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_simp_sdc @@ -1,5 +1,5 @@ ## Latest inputs file being used to reproduce initial conditions from Pakmor et al. 2022 -## with 12 km resolution using simplified sdc +## with 12.5 km resolution using simplified sdc ############################## CASTRO INPUTS ############################################### diff --git a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang index b56d9ba2a5..7977ab7787 100644 --- a/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang +++ b/Exec/science/wdmerger/tests/he_double_det/inputs_pakmor_strang @@ -1,5 +1,5 @@ ## Latest inputs file being used to reproduce initial conditions from Pakmor et al. 2022 -## with 12 km resolution using default strang-splitting +## with 12.5 km resolution using default strang-splitting ############################## CASTRO INPUTS ############################################### diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index 960d83acaf..c70fcd64eb 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -118,19 +118,6 @@ amrex::Array4 const& q_arr, amrex::Array4 const& qaux_arr); -/// -/// compute the flattening coefficient. This is 0 if we are in a shock and -/// 1 if we are not applying any flattening -/// -/// @param bx the box to operate over -/// @param q_arr the primitive variable state -/// @param flatn the flattening coefficient -/// @param pres_comp index into q_arr of the pressure component -/// - static void uflatten(const amrex::Box& bx, - amrex::Array4 const& q_arr, - amrex::Array4 const& flatn, - const int pres_comp); /// /// A multidimensional shock detection algorithm diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index fbcc8c6171..6eb352ac14 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -8,6 +8,7 @@ #include #include +#include using namespace amrex; @@ -117,19 +118,23 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) Array4 const flatn_arr = flatn.array(); if (first_order_hydro == 1) { - amrex::ParallelFor(obx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - flatn_arr(i,j,k) = 0.0; - }); + amrex::ParallelFor(obx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + flatn_arr(i,j,k) = 0.0; + }); } else if (use_flattening == 1) { - uflatten(obx, q_arr, flatn_arr, QPRES); + amrex::ParallelFor(obx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + flatn_arr(i,j,k) = hydro::flatten(i, j, k, q_arr, QPRES); + }); } else { - amrex::ParallelFor(obx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - flatn_arr(i,j,k) = 1.0; - }); + amrex::ParallelFor(obx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + flatn_arr(i,j,k) = 1.0; + }); } diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package index c5813e27c3..29c803c18a 100644 --- a/Source/hydro/Make.package +++ b/Source/hydro/Make.package @@ -10,7 +10,6 @@ CEXE_sources += Castro_mol.cpp CEXE_headers += advection_util.H CEXE_sources += advection_util.cpp CEXE_headers += flatten.H -CEXE_sources += flatten.cpp ifeq ($(USE_TRUE_SDC),TRUE) CEXE_sources += Castro_mol_hydro.cpp diff --git a/Source/hydro/flatten.cpp b/Source/hydro/flatten.cpp deleted file mode 100644 index 39318c9de0..0000000000 --- a/Source/hydro/flatten.cpp +++ /dev/null @@ -1,165 +0,0 @@ -#include - -#include - -#ifdef RADIATION -#include -#endif - -using namespace amrex; - -void -Castro::uflatten(const Box& bx, - Array4 const& q_arr, - Array4 const& flatn, const int pres_comp) { - - constexpr Real small_pres_flatn = 1.e-200_rt; - - // Knobs for detection of strong shock - constexpr Real shktst = 0.33_rt; - constexpr Real zcut1 = 0.75_rt; - constexpr Real zcut2 = 0.85_rt; - constexpr Real dzcut = 1.0_rt / (zcut2-zcut1); - - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - - // x-direction flattening coef - - Real dp = q_arr(i+1,j,k,pres_comp) - q_arr(i-1,j,k,pres_comp); - - int ishft = dp > 0.0_rt ? 1 : -1; - - Real denom = std::max(small_pres_flatn, std::abs(q_arr(i+2,j,k,pres_comp) - q_arr(i-2,j,k,pres_comp))); - Real zeta = std::abs(dp) / denom; - Real z = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt); - - Real tst = 0.0_rt; - if (q_arr(i-1,j,k,QU) - q_arr(i+1,j,k,QU) >= 0.0_rt) { - tst = 1.0_rt; - } - - Real tmp = std::min(q_arr(i+1,j,k,pres_comp), q_arr(i-1,j,k,pres_comp)); - - Real chi = 0.0_rt; - if (std::abs(dp) > shktst*tmp) { - chi = tst; - } - - - dp = q_arr(i+1-ishft,j,k,pres_comp) - q_arr(i-1-ishft,j,k,pres_comp); - - denom = std::max(small_pres_flatn, std::abs(q_arr(i+2-ishft,j,k,pres_comp)-q_arr(i-2-ishft,j,k,pres_comp))); - zeta = std::abs(dp) / denom; - Real z2 = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt); - - tst = 0.0_rt; - if (q_arr(i-1-ishft,j,k,QU) - q_arr(i+1-ishft,j,k,QU) >= 0.0_rt) { - tst = 1.0_rt; - } - - tmp = std::min(q_arr(i+1-ishft,j,k,pres_comp), q_arr(i-1-ishft,j,k,pres_comp)); - - Real chi2 = 0.0_rt; - if (std::abs(dp) > shktst*tmp) { - chi2 = tst; - } - - flatn(i,j,k) = 1.0_rt - std::max(chi2 * z2, chi * z); - - -#if AMREX_SPACEDIM >= 2 - // y-direction flattening coef - - dp = q_arr(i,j+1,k,pres_comp) - q_arr(i,j-1,k,pres_comp); - - ishft = dp > 0.0_rt ? 1 : -1; - - denom = std::max(small_pres_flatn, std::abs(q_arr(i,j+2,k,pres_comp) - q_arr(i,j-2,k,pres_comp))); - zeta = std::abs(dp) / denom; - z = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt); - - tst = 0.0_rt; - if (q_arr(i,j-1,k,QV) - q_arr(i,j+1,k,QV) >= 0.0_rt) { - tst = 1.0_rt; - } - - tmp = std::min(q_arr(i,j+1,k,pres_comp), q_arr(i,j-1,k,pres_comp)); - - chi = 0.0_rt; - if (std::abs(dp) > shktst*tmp) { - chi = tst; - } - - - dp = q_arr(i,j+1-ishft,k,pres_comp) - q_arr(i,j-1-ishft,k,pres_comp); - - denom = std::max(small_pres_flatn, std::abs(q_arr(i,j+2-ishft,k,pres_comp) - q_arr(i,j-2-ishft,k,pres_comp))); - zeta = std::abs(dp) / denom; - z2 = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt); - - tst = 0.0_rt; - if (q_arr(i,j-1-ishft,k,QV) - q_arr(i,j+1-ishft,k,QV) >= 0.0_rt) { - tst = 1.0_rt; - } - - tmp = std::min(q_arr(i,j+1-ishft,k,pres_comp), q_arr(i,j-1-ishft,k,pres_comp)); - - chi2 = 0.0_rt; - if (std::abs(dp) > shktst*tmp) { - chi2 = tst; - } - - flatn(i,j,k) = std::min(flatn(i,j,k), 1.0_rt - std::max(chi2 * z2, chi * z)); -#endif - - -#if AMREX_SPACEDIM == 3 - // z-direction flattening coef - - dp = q_arr(i,j,k+1,pres_comp) - q_arr(i,j,k-1,pres_comp); - - ishft = dp > 0.0_rt ? 1: -1; - - denom = std::max(small_pres_flatn, std::abs(q_arr(i,j,k+2,pres_comp) - q_arr(i,j,k-2,pres_comp))); - zeta = std::abs(dp) / denom; - z = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt); - - tst = 0.0_rt; - if (q_arr(i,j,k-1,QW) - q_arr(i,j,k+1,QW) >= 0.0_rt) { - tst = 1.0_rt; - } - - tmp = std::min(q_arr(i,j,k+1,pres_comp), q_arr(i,j,k-1,pres_comp)); - - chi = 0.0_rt; - if (std::abs(dp) > shktst*tmp) { - chi = tst; - } - - - dp = q_arr(i,j,k+1-ishft,pres_comp) - q_arr(i,j,k-1-ishft,pres_comp); - - denom = std::max(small_pres_flatn, std::abs(q_arr(i,j,k+2-ishft,pres_comp) - q_arr(i,j,k-2-ishft,pres_comp))); - zeta = std::abs(dp) / denom; - z2 = std::clamp(dzcut * (zeta - zcut1), 0.0_rt, 1.0_rt); - - tst = 0.0_rt; - if (q_arr(i,j,k-1-ishft,QW) - q_arr(i,j,k+1-ishft,QW) >= 0.0_rt) { - tst = 1.0_rt; - } - - tmp = std::min(q_arr(i,j,k+1-ishft,pres_comp), q_arr(i,j,k-1-ishft,pres_comp)); - - chi2 = 0.0_rt; - if (std::abs(dp) > shktst*tmp) { - chi2 = tst; - } - - flatn(i,j,k) = std::min(flatn(i,j,k), 1.0_rt - std::max(chi2 * z2, chi * z)); -#endif - - }); - -} diff --git a/Source/mhd/Castro_mhd.H b/Source/mhd/Castro_mhd.H index 12037acc75..3d72a01726 100644 --- a/Source/mhd/Castro_mhd.H +++ b/Source/mhd/Castro_mhd.H @@ -99,7 +99,7 @@ amrex::Array4 const& Ed2, const int d, const int d1, const int d2, const amrex::Real dt); - void + static void hlld(const amrex::Box& bx, amrex::Array4 const& qleft, amrex::Array4 const& qright, diff --git a/Source/mhd/Castro_mhd.cpp b/Source/mhd/Castro_mhd.cpp index 220a65ffb0..8d6af7594b 100644 --- a/Source/mhd/Castro_mhd.cpp +++ b/Source/mhd/Castro_mhd.cpp @@ -1,6 +1,7 @@ #include #include +#include using namespace amrex; @@ -196,10 +197,6 @@ Castro::construct_ctu_mhd_source(Real time, Real dt) auto flatn_arr = flatn.array(); auto elix_flatn = flatn.elixir(); - flatg.resize(bxi, 1); - auto flatg_arr = flatg.array(); - auto elix_flatg = flatg.elixir(); - if (use_flattening == 0) { amrex::ParallelFor(bxi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -209,13 +206,11 @@ Castro::construct_ctu_mhd_source(Real time, Real dt) } else { - uflatten(bxi, q_arr, flatn_arr, QPRES); - uflatten(bxi, q_arr, flatg_arr, QPTOT); - amrex::ParallelFor(bxi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - flatn_arr(i,j,k) = flatn_arr(i,j,k) * flatg_arr(i,j,k); + flatn_arr(i,j,k) = hydro::flatten(i, j, k, q_arr, QPRES) * + hydro::flatten(i, j, k, q_arr, QPTOT); }); } diff --git a/Source/mhd/mhd_eigen.H b/Source/mhd/mhd_eigen.H index 3fe6cfd919..7757a71e12 100644 --- a/Source/mhd/mhd_eigen.H +++ b/Source/mhd/mhd_eigen.H @@ -86,21 +86,24 @@ evecx(Array2D& leig, Real alf; Real als; - if (std::abs(cfx-csx) <= 1e-14){ + // if C_f = C_s, then we are degenerate, and set the alpha's accordingly. + // see Stone et a. 2008 comment just after Eql A17. + + if (std::abs(cfx - csx) <= 1e-14){ alf = 1.0_rt; als = 0.0_rt; - } - - if (as - csx < 0.0) { - alf = 0.0_rt; } else { - alf = std::sqrt((as - csx)/(cfx - csx)); - } + if (as - csx < 0.0) { + alf = 0.0_rt; + } else { + alf = std::sqrt((as - csx)/(cfx - csx)); + } - if (cfx - as < 0.0) { - als = 0.0_rt; - } else { - als = std::sqrt((cfx - as)/(cfx - csx)); + if (cfx - as < 0.0) { + als = 0.0_rt; + } else { + als = std::sqrt((cfx - as)/(cfx - csx)); + } } @@ -287,21 +290,24 @@ evecy(Array2D& leig, Real alf; Real als; - if (std::abs(cfy-csy) <= 1e-14){ + // if C_f = C_s, then we are degenerate, and set the alpha's accordingly. + // see Stone et a. 2008 comment just after Eql A17. + + if (std::abs(cfy - csy) <= 1e-14){ alf = 1.0_rt; als = 0.0_rt; - } - - if (as - csy < 0.0) { - alf = 0.0_rt; } else { - alf = std::sqrt((as - csy)/(cfy - csy)); - } + if (as - csy < 0.0) { + alf = 0.0_rt; + } else { + alf = std::sqrt((as - csy)/(cfy - csy)); + } - if (cfy - as < 0.0) { - als = 0.0_rt; - } else { - als = std::sqrt((cfy - as)/(cfy - csy)); + if (cfy - as < 0.0) { + als = 0.0_rt; + } else { + als = std::sqrt((cfy - as)/(cfy - csy)); + } } Real betx; @@ -485,21 +491,24 @@ evecz(Array2D& leig, Real alf; Real als; + // if C_f = C_s, then we are degenerate, and set the alpha's accordingly. + // see Stone et a. 2008 comment just after Eql A17. + if (std::abs(cfz-csz) <= 1e-14){ alf = 1.0_rt; als = 0.0_rt; - } - - if (as - csz < 0.0) { - alf = 0.0_rt; } else { - alf = std::sqrt((as - csz)/(cfz - csz)); - } + if (as - csz < 0.0) { + alf = 0.0_rt; + } else { + alf = std::sqrt((as - csz)/(cfz - csz)); + } - if (cfz - as < 0.0) { - als = 0.0_rt; - } else { - als = std::sqrt((cfz - as)/(cfz - csz)); + if (cfz - as < 0.0) { + als = 0.0_rt; + } else { + als = std::sqrt((cfz - as)/(cfz - csz)); + } } Real betx; diff --git a/Source/mhd/mhd_plm.cpp b/Source/mhd/mhd_plm.cpp index 2747812513..bf61c799f9 100644 --- a/Source/mhd/mhd_plm.cpp +++ b/Source/mhd/mhd_plm.cpp @@ -154,8 +154,8 @@ Castro::plm(const Box& bx, smhd[IEIGN_BTT] = q_zone(QW); // cross-talk of normal magnetic field direction - for (int n = 0; n < NEIGN; n++) { - smhd[n] = smhd[n] * (Bx(i+1,j,k) - Bx(i,j,k)) / dx[idir]; + for (auto & source : smhd) { + source *= (Bx(i+1,j,k) - Bx(i,j,k)) / dx[idir]; } } else if (idir == 1) { @@ -163,8 +163,8 @@ Castro::plm(const Box& bx, smhd[IEIGN_BTT] = q_zone(QW); // cross-talk of normal magnetic field direction - for (int n = 0; n < NEIGN; n++) { - smhd[n] = smhd[n] * (By(i,j+1,k) - By(i,j,k)) / dx[idir]; + for (auto & source : smhd) { + source *= (By(i,j+1,k) - By(i,j,k)) / dx[idir]; } } else { @@ -172,8 +172,8 @@ Castro::plm(const Box& bx, smhd[IEIGN_BTT] = q_zone(QV); // cross-talk of normal magnetic field direction - for (int n = 0; n < NEIGN; n++) { - smhd[n] = smhd[n] * (Bz(i,j,k+1) - Bz(i,j,k)) / dx[idir]; + for (auto & source : smhd) { + source *= (Bz(i,j,k+1) - Bz(i,j,k)) / dx[idir]; } } @@ -441,4 +441,3 @@ Castro::plm(const Box& bx, }); } - diff --git a/Source/mhd/mhd_ppm.cpp b/Source/mhd/mhd_ppm.cpp index f084a95a39..30cd9668d6 100644 --- a/Source/mhd/mhd_ppm.cpp +++ b/Source/mhd/mhd_ppm.cpp @@ -200,8 +200,8 @@ Castro::ppm_mhd(const Box& bx, smhd[IEIGN_BTT] = q_zone(QW); // cross-talk of normal magnetic field direction - for (int n = 0; n < NEIGN; n++) { - smhd[n] = smhd[n] * (Bx(i+1,j,k) - Bx(i,j,k)) / dx[idir]; + for (auto & source : smhd) { + source *= (Bx(i+1,j,k) - Bx(i,j,k)) / dx[idir]; } } else if (idir == 1) { @@ -209,8 +209,8 @@ Castro::ppm_mhd(const Box& bx, smhd[IEIGN_BTT] = q_zone(QW); // cross-talk of normal magnetic field direction - for (int n = 0; n < NEIGN; n++) { - smhd[n] = smhd[n] * (By(i,j+1,k) - By(i,j,k)) / dx[idir]; + for (auto & source : smhd) { + source *= (By(i,j+1,k) - By(i,j,k)) / dx[idir]; } } else { @@ -218,8 +218,8 @@ Castro::ppm_mhd(const Box& bx, smhd[IEIGN_BTT] = q_zone(QV); // cross-talk of normal magnetic field direction - for (int n = 0; n < NEIGN; n++) { - smhd[n] = smhd[n] * (Bz(i,j,k+1) - Bz(i,j,k)) / dx[idir]; + for (auto & source : smhd) { + source *= (Bz(i,j,k+1) - Bz(i,j,k)) / dx[idir]; } } @@ -505,4 +505,3 @@ Castro::ppm_mhd(const Box& bx, }); } -