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 aa51e8cd3c..0000000000 --- a/Source/hydro/flatten.cpp +++ /dev/null @@ -1,166 +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 = amrex::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 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1))); - - 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 = amrex::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 = amrex::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 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1))); - - 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 = amrex::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 - amrex::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 = amrex::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 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1))); - - 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 = amrex::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 = amrex::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 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1))); - - 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 = amrex::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) = amrex::min(flatn(i,j,k), 1.0_rt - amrex::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 = amrex::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 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1))); - - 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 = amrex::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 = amrex::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 = amrex::min(1.0_rt, amrex::max(0.0_rt, dzcut * (zeta - zcut1))); - - 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 = amrex::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) = amrex::min(flatn(i,j,k), 1.0_rt - amrex::max(chi2 * z2, chi * z)); -#endif - - }); - -} - 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); }); }