From cc8bbef91623e2c0837076148bf4e6651c2c912c Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 20 Aug 2024 08:24:41 -0700 Subject: [PATCH] Move more functionality into AMReX-Hydro (#1202) --------- Co-authored-by: Marc Henry de Frahan --- amr-wind/convection/CMakeLists.txt | 1 - amr-wind/convection/Godunov.H | 22 - amr-wind/convection/incflo_convection_K.H | 537 ---------------- .../convection/incflo_godunov_advection.cpp | 548 ----------------- amr-wind/convection/incflo_godunov_minmod.H | 207 ------- amr-wind/convection/incflo_godunov_ppm.H | 575 ------------------ amr-wind/convection/incflo_godunov_upwind.H | 30 - amr-wind/convection/incflo_mol_fluxes.cpp | 47 +- .../equation_systems/icns/icns_advection.H | 14 +- .../equation_systems/vof/vof_momentum_flux.H | 87 ++- unit_tests/multiphase/test_mflux_schemes.cpp | 103 ++-- 11 files changed, 162 insertions(+), 2009 deletions(-) delete mode 100644 amr-wind/convection/incflo_convection_K.H delete mode 100644 amr-wind/convection/incflo_godunov_advection.cpp delete mode 100644 amr-wind/convection/incflo_godunov_minmod.H delete mode 100644 amr-wind/convection/incflo_godunov_ppm.H delete mode 100644 amr-wind/convection/incflo_godunov_upwind.H diff --git a/amr-wind/convection/CMakeLists.txt b/amr-wind/convection/CMakeLists.txt index 82c968f7d5..5e2d510cde 100644 --- a/amr-wind/convection/CMakeLists.txt +++ b/amr-wind/convection/CMakeLists.txt @@ -1,6 +1,5 @@ target_sources(${amr_wind_lib_name} PRIVATE #C++ - incflo_godunov_advection.cpp incflo_mol_fluxes.cpp ) diff --git a/amr-wind/convection/Godunov.H b/amr-wind/convection/Godunov.H index e5e98af000..dcbf32bdbb 100644 --- a/amr-wind/convection/Godunov.H +++ b/amr-wind/convection/Godunov.H @@ -6,32 +6,10 @@ #ifndef Godunov_H #define Godunov_H -#include "amr-wind/core/FieldRepo.H" - namespace godunov { enum class scheme { PLM, PPM, PPM_NOLIM, BDS, WENO_JS, WENOZ, MINMOD, UPWIND }; -void compute_fluxes( - int lev, - amrex::Box const& bx, - int ncomp, - amrex::Array4 const& fx, - amrex::Array4 const& fy, - amrex::Array4 const& fz, - amrex::Array4 const& q, - amrex::Array4 const& umac, - amrex::Array4 const& vmac, - amrex::Array4 const& wmac, - amrex::Array4 const& fq, - amrex::BCRec const* pbc, - int const* iconserv, - amrex::Real* p, - amrex::Vector geom, - amrex::Real dt, - godunov::scheme godunov_scheme, - bool godunov_use_forces_in_trans); - } // namespace godunov #endif /* Godunov_H */ diff --git a/amr-wind/convection/incflo_convection_K.H b/amr-wind/convection/incflo_convection_K.H deleted file mode 100644 index 30a87d752f..0000000000 --- a/amr-wind/convection/incflo_convection_K.H +++ /dev/null @@ -1,537 +0,0 @@ -#ifndef INCFLO_CONVECTION_K_H_ -#define INCFLO_CONVECTION_K_H_ - -#include - -namespace { - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_xslope( - int i, - int j, - int k, - int n, - amrex::Array4 const& vcc) noexcept -{ - amrex::Real dl = 2.0 * (vcc(i, j, k, n) - vcc(i - 1, j, k, n)); - amrex::Real dr = 2.0 * (vcc(i + 1, j, k, n) - vcc(i, j, k, n)); - amrex::Real dc = 0.5 * (vcc(i + 1, j, k, n) - vcc(i - 1, j, k, n)); - amrex::Real slope = amrex::min(std::abs(dl), std::abs(dc), std::abs(dr)); - slope = (dr * dl > 0.0) ? slope : 0.0; - return (dc > 0.0) ? slope : -slope; -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_ho_xslope( - int i, - int j, - int k, - int n, - amrex::Array4 const& q) noexcept -{ - amrex::Real dlft, drgt, dcen, dfm, dfp, dlim, dsgn, dtemp; - amrex::Real qm, qp, qi; - qi = q(i, j, k, n); - qm = q(i - 1, j, k, n); - qp = q(i + 1, j, k, n); - - dlft = qm - q(i - 2, j, k, n); - drgt = qi - qm; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfm = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qp - qi; - drgt = q(i + 2, j, k, n) - qp; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfp = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qi - qm; - drgt = qp - qi; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - - return dsgn * amrex::min(dlim, std::abs(dtemp)); -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_xslope_extdir( - int i, - int j, - int k, - int n, - amrex::Array4 const& vcc, - bool edlo, - bool edhi, - int domlo, - int domhi) noexcept -{ - amrex::Real dl = 2.0 * (vcc(i, j, k, n) - vcc(i - 1, j, k, n)); - amrex::Real dr = 2.0 * (vcc(i + 1, j, k, n) - vcc(i, j, k, n)); - amrex::Real dc = 0.5 * (vcc(i + 1, j, k, n) - vcc(i - 1, j, k, n)); - if (edlo and i == domlo) { - dc = (vcc(i + 1, j, k, n) + 3.0 * vcc(i, j, k, n) - - 4.0 * vcc(i - 1, j, k, n)) / - 3.0; - } else if (edhi and i == domhi) { - dc = (4.0 * vcc(i + 1, j, k, n) - 3.0 * vcc(i, j, k, n) - - vcc(i - 1, j, k, n)) / - 3.0; - } - amrex::Real slope = amrex::min(std::abs(dl), std::abs(dc), std::abs(dr)); - slope = (dr * dl > 0.0) ? slope : 0.0; - return (dc > 0.0) ? slope : -slope; -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_ho_xslope_extdir( - int i, - int j, - int k, - int n, - amrex::Array4 const& q, - bool edlo, - bool edhi, - int domlo, - int domhi) noexcept -{ - amrex::Real dlft, drgt, dcen, dfm, dfp, dlim, dsgn, dtemp, dlimsh, dsgnsh; - amrex::Real qm, qp, qi; - qi = q(i, j, k, n); - qm = q(i - 1, j, k, n); - qp = q(i + 1, j, k, n); - - dlft = qm - q(i - 2, j, k, n); - drgt = qi - qm; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfm = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qp - qi; - drgt = q(i + 2, j, k, n) - qp; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfp = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qi - qm; - drgt = qp - qi; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - - if (edlo and i == domlo) { - dtemp = -16. / 15. * q(i - 1, j, k, n) + .5 * q(i, j, k, n) + - 2. / 3. * q(i + 1, j, k, n) - 0.1 * q(i + 2, j, k, n); - dlft = 2. * (q(i, j, k, n) - q(i - 1, j, k, n)); - drgt = 2. * (q(i + 1, j, k, n) - q(i, j, k, n)); - dlim = (dlft * drgt >= 0.0) ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgn = std::copysign(1.e0, dtemp); - } else if (edlo and i == domlo + 1) { - dfm = -16. / 15. * q(domlo - 1, j, k, n) + .5 * q(domlo, j, k, n) + - 2. / 3. * q(domlo + 1, j, k, n) - 0.1 * q(domlo + 2, j, k, n); - dlft = 2. * (q(domlo, j, k, n) - q(domlo - 1, j, k, n)); - drgt = 2. * (q(domlo + 1, j, k, n) - q(domlo, j, k, n)); - dlimsh = (dlft * drgt >= 0.0) - ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgnsh = std::copysign(1.e0, dfm); - dfm = dsgnsh * amrex::min(dlimsh, std::abs(dfm)); - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - } - - if (edhi and i == domhi) { - dtemp = 16. / 15. * q(i + 1, j, k, n) - .5 * q(i, j, k, n) - - 2. / 3. * q(i - 1, j, k, n) + 0.1 * q(i - 2, j, k, n); - dlft = 2. * (q(i, j, k, n) - q(i - 1, j, k, n)); - drgt = 2. * (q(i + 1, j, k, n) - q(i, j, k, n)); - dlim = (dlft * drgt >= 0.0) ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgn = std::copysign(1.e0, dtemp); - } else if (edhi and i == domhi - 1) { - dfp = 16. / 15. * q(domhi + 1, j, k, n) - .5 * q(domhi, j, k, n) - - 2. / 3. * q(domhi - 1, j, k, n) + 0.1 * q(domhi - 2, j, k, n); - dlft = 2. * (q(domhi, j, k, n) - q(domhi - 1, j, k, n)); - drgt = 2. * (q(domhi + 1, j, k, n) - q(domhi, j, k, n)); - dlimsh = (dlft * drgt >= 0.0) - ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgnsh = std::copysign(1.e0, dfp); - dfp = dsgnsh * amrex::min(dlimsh, std::abs(dfp)); - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - } - - return dsgn * amrex::min(dlim, std::abs(dtemp)); -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_yslope( - int i, - int j, - int k, - int n, - amrex::Array4 const& vcc) noexcept -{ - amrex::Real dl = 2.0 * (vcc(i, j, k, n) - vcc(i, j - 1, k, n)); - amrex::Real dr = 2.0 * (vcc(i, j + 1, k, n) - vcc(i, j, k, n)); - amrex::Real dc = 0.5 * (vcc(i, j + 1, k, n) - vcc(i, j - 1, k, n)); - amrex::Real slope = amrex::min(std::abs(dl), std::abs(dc), std::abs(dr)); - slope = (dr * dl > 0.0) ? slope : 0.0; - return (dc > 0.0) ? slope : -slope; -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_ho_yslope( - int i, - int j, - int k, - int n, - amrex::Array4 const& q) noexcept -{ - amrex::Real dlft, drgt, dcen, dfm, dfp, dlim, dsgn, dtemp; - amrex::Real qm, qp, qj; - qj = q(i, j, k, n); - qm = q(i, j - 1, k, n); - qp = q(i, j + 1, k, n); - - dlft = qm - q(i, j - 2, k, n); - drgt = qj - qm; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfm = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qp - qj; - drgt = q(i, j + 2, k, n) - qp; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfp = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qj - qm; - drgt = qp - qj; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - return dsgn * amrex::min(dlim, std::abs(dtemp)); -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_yslope_extdir( - int i, - int j, - int k, - int n, - amrex::Array4 const& vcc, - bool edlo, - bool edhi, - int domlo, - int domhi) noexcept -{ - amrex::Real dl = 2.0 * (vcc(i, j, k, n) - vcc(i, j - 1, k, n)); - amrex::Real dr = 2.0 * (vcc(i, j + 1, k, n) - vcc(i, j, k, n)); - amrex::Real dc = 0.5 * (vcc(i, j + 1, k, n) - vcc(i, j - 1, k, n)); - if (edlo and j == domlo) { - dc = (vcc(i, j + 1, k, n) + 3.0 * vcc(i, j, k, n) - - 4.0 * vcc(i, j - 1, k, n)) / - 3.0; - } else if (edhi and j == domhi) { - dc = (4.0 * vcc(i, j + 1, k, n) - 3.0 * vcc(i, j, k, n) - - vcc(i, j - 1, k, n)) / - 3.0; - } - amrex::Real slope = amrex::min(std::abs(dl), std::abs(dc), std::abs(dr)); - slope = (dr * dl > 0.0) ? slope : 0.0; - return (dc > 0.0) ? slope : -slope; -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_ho_yslope_extdir( - int i, - int j, - int k, - int n, - amrex::Array4 const& q, - bool edlo, - bool edhi, - int domlo, - int domhi) noexcept -{ - amrex::Real dlft, drgt, dcen, dfm, dfp, dlim, dsgn, dtemp, dlimsh, dsgnsh; - amrex::Real qm, qp, qj; - qj = q(i, j, k, n); - qm = q(i, j - 1, k, n); - qp = q(i, j + 1, k, n); - - dlft = qm - q(i, j - 2, k, n); - drgt = qj - qm; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfm = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qp - qj; - drgt = q(i, j + 2, k, n) - qp; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfp = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qj - qm; - drgt = qp - qj; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - - if (edlo and j == domlo) { - dtemp = -16. / 15. * q(i, j - 1, k, n) + .5 * q(i, j, k, n) + - 2. / 3. * q(i, j + 1, k, n) - 0.1 * q(i, j + 2, k, n); - dlft = 2. * (q(i, j, k, n) - q(i, j - 1, k, n)); - drgt = 2. * (q(i, j + 1, k, n) - q(i, j, k, n)); - dlim = (dlft * drgt >= 0.0) ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgn = std::copysign(1.e0, dtemp); - } else if (edlo and j == domlo + 1) { - dfm = -16. / 15. * q(i, domlo - 1, k, n) + .5 * q(i, domlo, k, n) + - 2. / 3. * q(i, domlo + 1, k, n) - 0.1 * q(i, domlo + 2, k, n); - dlft = 2. * (q(i, domlo, k, n) - q(i, domlo - 1, k, n)); - drgt = 2. * (q(i, domlo + 1, k, n) - q(i, domlo, k, n)); - dlimsh = (dlft * drgt >= 0.0) - ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgnsh = std::copysign(1.e0, dfm); - dfm = dsgnsh * amrex::min(dlimsh, std::abs(dfm)); - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - } - - if (edhi and j == domhi) { - dtemp = 16. / 15. * q(i, j + 1, k, n) - .5 * q(i, j, k, n) - - 2. / 3. * q(i, j - 1, k, n) + 0.1 * q(i, j - 2, k, n); - dlft = 2. * (q(i, j, k, n) - q(i, j - 1, k, n)); - drgt = 2. * (q(i, j + 1, k, n) - q(i, j, k, n)); - dlim = (dlft * drgt >= 0.0) ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgn = std::copysign(1.e0, dtemp); - } else if (edhi and j == domhi - 1) { - dfp = 16. / 15. * q(i, domhi + 1, k, n) - .5 * q(i, domhi, k, n) - - 2. / 3. * q(i, domhi - 1, k, n) + 0.1 * q(i, domhi - 2, k, n); - dlft = 2. * (q(i, domhi, k, n) - q(i, domhi - 1, k, n)); - drgt = 2. * (q(i, domhi + 1, k, n) - q(i, domhi, k, n)); - dlimsh = (dlft * drgt >= 0.0) - ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgnsh = std::copysign(1.e0, dfp); - dfp = dsgnsh * amrex::min(dlimsh, std::abs(dfp)); - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - } - - return dsgn * amrex::min(dlim, std::abs(dtemp)); -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_zslope( - int i, - int j, - int k, - int n, - amrex::Array4 const& vcc) noexcept -{ - amrex::Real dl = 2.0 * (vcc(i, j, k, n) - vcc(i, j, k - 1, n)); - amrex::Real dr = 2.0 * (vcc(i, j, k + 1, n) - vcc(i, j, k, n)); - amrex::Real dc = 0.5 * (vcc(i, j, k + 1, n) - vcc(i, j, k - 1, n)); - amrex::Real slope = amrex::min(std::abs(dl), std::abs(dc), std::abs(dr)); - slope = (dr * dl > 0.0) ? slope : 0.0; - return (dc > 0.0) ? slope : -slope; -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_ho_zslope( - int i, - int j, - int k, - int n, - amrex::Array4 const& q) noexcept -{ - amrex::Real dlft, drgt, dcen, dfm, dfp, dlim, dsgn, dtemp; - amrex::Real qm, qp, qk; - qk = q(i, j, k, n); - qm = q(i, j, k - 1, n); - qp = q(i, j, k + 1, n); - - dlft = qm - q(i, j, k - 2, n); - drgt = qk - qm; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfm = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qp - qk; - drgt = q(i, j, k + 2, n) - qp; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfp = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qk - qm; - drgt = qp - qk; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - return dsgn * amrex::min(dlim, std::abs(dtemp)); -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_zslope_extdir( - int i, - int j, - int k, - int n, - amrex::Array4 const& vcc, - bool edlo, - bool edhi, - int domlo, - int domhi) noexcept -{ - amrex::Real dl = 2.0 * (vcc(i, j, k, n) - vcc(i, j, k - 1, n)); - amrex::Real dr = 2.0 * (vcc(i, j, k + 1, n) - vcc(i, j, k, n)); - amrex::Real dc = 0.5 * (vcc(i, j, k + 1, n) - vcc(i, j, k - 1, n)); - if (edlo and k == domlo) { - dc = (vcc(i, j, k + 1, n) + 3.0 * vcc(i, j, k, n) - - 4.0 * vcc(i, j, k - 1, n)) / - 3.0; - } else if (edhi and k == domhi) { - dc = (4.0 * vcc(i, j, k + 1, n) - 3.0 * vcc(i, j, k, n) - - vcc(i, j, k - 1, n)) / - 3.0; - } - amrex::Real slope = amrex::min(std::abs(dl), std::abs(dc), std::abs(dr)); - slope = (dr * dl > 0.0) ? slope : 0.0; - return (dc > 0.0) ? slope : -slope; -} - -AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real incflo_ho_zslope_extdir( - int i, - int j, - int k, - int n, - amrex::Array4 const& q, - bool edlo, - bool edhi, - int domlo, - int domhi) noexcept -{ - amrex::Real dlft, drgt, dcen, dfm, dfp, dlim, dsgn, dtemp, dlimsh, dsgnsh; - amrex::Real qm, qp, qk; - qk = q(i, j, k, n); - qm = q(i, j, k - 1, n); - qp = q(i, j, k + 1, n); - - dlft = qm - q(i, j, k - 2, n); - drgt = qk - qm; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfm = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qp - qk; - drgt = q(i, j, k + 2, n) - qp; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dfp = dsgn * amrex::min(dlim, std::abs(dcen)); - - dlft = qk - qm; - drgt = qp - qk; - dcen = 0.5 * (dlft + drgt); - dsgn = std::copysign(1.e0, dcen); - dlim = (dlft * drgt >= 0.0) - ? 2.0 * amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - - if (edlo and k == domlo) { - dtemp = -16. / 15. * q(i, j, k - 1, n) + .5 * q(i, j, k, n) + - 2. / 3. * q(i, j, k + 1, n) - 0.1 * q(i, j, k + 2, n); - dlft = 2. * (q(i, j, k, n) - q(i, j, k - 1, n)); - drgt = 2. * (q(i, j, k + 1, n) - q(i, j, k, n)); - dlim = (dlft * drgt >= 0.0) ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgn = std::copysign(1.e0, dtemp); - } else if (edlo and k == domlo + 1) { - dfm = -16. / 15. * q(i, j, domlo - 1, n) + .5 * q(i, j, domlo, n) + - 2. / 3. * q(i, j, domlo + 1, n) - 0.1 * q(i, j, domlo + 2, n); - dlft = 2. * (q(i, j, domlo, n) - q(i, j, domlo - 1, n)); - drgt = 2. * (q(i, j, domlo + 1, n) - q(i, j, domlo, n)); - dlimsh = (dlft * drgt >= 0.0) - ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgnsh = std::copysign(1.e0, dfm); - dfm = dsgnsh * amrex::min(dlimsh, std::abs(dfm)); - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - } - - if (edhi and k == domhi) { - dtemp = 16. / 15. * q(i, j, k + 1, n) - .5 * q(i, j, k, n) - - 2. / 3. * q(i, j, k - 1, n) + 0.1 * q(i, j, k - 2, n); - dlft = 2. * (q(i, j, k, n) - q(i, j, k - 1, n)); - drgt = 2. * (q(i, j, k + 1, n) - q(i, j, k, n)); - dlim = (dlft * drgt >= 0.0) ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgn = std::copysign(1.e0, dtemp); - } else if (edhi and k == domhi - 1) { - dfp = 16. / 15. * q(i, j, domhi + 1, n) - .5 * q(i, j, domhi, n) - - 2. / 3. * q(i, j, domhi - 1, n) + 0.1 * q(i, j, domhi - 2, n); - dlft = 2. * (q(i, j, domhi, n) - q(i, j, domhi - 1, n)); - drgt = 2. * (q(i, j, domhi + 1, n) - q(i, j, domhi, n)); - dlimsh = (dlft * drgt >= 0.0) - ? amrex::min(std::abs(dlft), std::abs(drgt)) - : 0.0; - dsgnsh = std::copysign(1.e0, dfp); - dfp = dsgnsh * amrex::min(dlimsh, std::abs(dfp)); - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); - } - return dsgn * amrex::min(dlim, std::abs(dtemp)); -} - -} // namespace -#endif diff --git a/amr-wind/convection/incflo_godunov_advection.cpp b/amr-wind/convection/incflo_godunov_advection.cpp deleted file mode 100644 index 44dc773919..0000000000 --- a/amr-wind/convection/incflo_godunov_advection.cpp +++ /dev/null @@ -1,548 +0,0 @@ -#include "amr-wind/convection/incflo_godunov_ppm.H" -#include "amr-wind/convection/incflo_godunov_minmod.H" -#include "amr-wind/convection/incflo_godunov_upwind.H" -#include "amr-wind/convection/Godunov.H" -#include - -using namespace amrex; - -void godunov::compute_fluxes( - int lev, - Box const& bx, - int ncomp, - Array4 const& fx, - Array4 const& fy, - Array4 const& fz, - Array4 const& q, - Array4 const& umac, - Array4 const& vmac, - Array4 const& wmac, - Array4 const& fq, - BCRec const* pbc, - int const* iconserv, - Real* p, - Vector geom, - Real dt, - godunov::scheme godunov_scheme, - bool use_forces_in_trans) -{ - - /* iconserv functionality: (would be better served by two flags) - * ------------------------------------------------------------------------ - * == 0 : non-conservative formulation of interpolation, fluxes not - * multiplied by MAC velocity - * == 1 : conservative formulation of interpolation, fluxes include factor - * of MAC velocity - * (!= 1 && != 0) : conservative formulation of interpolation, fluxes not - * multiplied by MAC velocity */ - - BL_PROFILE("amr-wind::godunov::compute_fluxes"); - Box const& xbx = amrex::surroundingNodes(bx, 0); - Box const& ybx = amrex::surroundingNodes(bx, 1); - Box const& zbx = amrex::surroundingNodes(bx, 2); - Box const& bxg1 = amrex::grow(bx, 1); - Box xebox = Box(xbx).grow(1, 1).grow(2, 1); - Box yebox = Box(ybx).grow(0, 1).grow(2, 1); - Box zebox = Box(zbx).grow(0, 1).grow(1, 1); - - const Real dx = geom[lev].CellSize(0); - const Real dy = geom[lev].CellSize(1); - const Real dz = geom[lev].CellSize(2); - Real l_dt = dt; - Real dtdx = l_dt / dx; - Real dtdy = l_dt / dy; - Real dtdz = l_dt / dz; - - Box const& domain = geom[lev].Domain(); - const auto dlo = amrex::lbound(domain); - const auto dhi = amrex::ubound(domain); - - Array4 Imx = makeArray4(p, bxg1, ncomp); - p += Imx.size(); - Array4 Ipx = makeArray4(p, bxg1, ncomp); - p += Ipx.size(); - Array4 Imy = makeArray4(p, bxg1, ncomp); - p += Imy.size(); - Array4 Ipy = makeArray4(p, bxg1, ncomp); - p += Ipy.size(); - Array4 Imz = makeArray4(p, bxg1, ncomp); - p += Imz.size(); - Array4 Ipz = makeArray4(p, bxg1, ncomp); - p += Ipz.size(); - Array4 xlo = makeArray4(p, xebox, ncomp); - p += xlo.size(); - Array4 xhi = makeArray4(p, xebox, ncomp); - p += xhi.size(); - Array4 ylo = makeArray4(p, yebox, ncomp); - p += ylo.size(); - Array4 yhi = makeArray4(p, yebox, ncomp); - p += yhi.size(); - Array4 zlo = makeArray4(p, zebox, ncomp); - p += zlo.size(); - Array4 zhi = makeArray4(p, zebox, ncomp); - p += zhi.size(); - Array4 xyzlo = makeArray4(p, bxg1, ncomp); - p += xyzlo.size(); - Array4 xyzhi = makeArray4(p, bxg1, ncomp); - - // Use PPM to generate Im and Ip */ - switch (godunov_scheme) { - case godunov::scheme::MINMOD: { - amrex::ParallelFor( - bxg1, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Godunov_minmod_fpu_x( - i, j, k, n, l_dt, dx, Imx(i, j, k, n), Ipx(i, j, k, n), q, - umac, pbc[n], dlo.x, dhi.x); - Godunov_minmod_fpu_y( - i, j, k, n, l_dt, dy, Imy(i, j, k, n), Ipy(i, j, k, n), q, - vmac, pbc[n], dlo.y, dhi.y); - Godunov_minmod_fpu_z( - i, j, k, n, l_dt, dz, Imz(i, j, k, n), Ipz(i, j, k, n), q, - wmac, pbc[n], dlo.z, dhi.z); - }); - break; - } - case godunov::scheme::UPWIND: { - amrex::ParallelFor( - bxg1, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Godunov_upwind_fpu( - i, j, k, n, Imx(i, j, k, n), Ipx(i, j, k, n), q); - Godunov_upwind_fpu( - i, j, k, n, Imy(i, j, k, n), Ipy(i, j, k, n), q); - Godunov_upwind_fpu( - i, j, k, n, Imz(i, j, k, n), Ipz(i, j, k, n), q); - }); - break; - } - default: { - amrex::Abort( - "This code path only used in multiphase simulations with MINMOD " - "and UPWIND"); - } - } - - amrex::ParallelFor( - xebox, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - constexpr Real small_vel = 1.e-8; - - Real uad = umac(i, j, k); - Real fux = (std::abs(uad) < small_vel) ? 0. : 1.; - bool uval = uad >= 0.; - // divu = 0 - // Real cons1 = (iconserv[n]) ? -0.5*l_dt*q(i-1,j,k,n)*divu(i-1,j,k) - // : 0.; Real cons2 = (iconserv[n]) ? -0.5*l_dt*q(i ,j,k,n)*divu(i - // ,j,k) : 0.; - Real lo = Ipx(i - 1, j, k, n); // + cons1; - Real hi = Imx(i, j, k, n); // + cons2; - if (use_forces_in_trans && fq) { - lo += 0.5 * l_dt * fq(i - 1, j, k, n); - hi += 0.5 * l_dt * fq(i, j, k, n); - } - - auto bc = pbc[n]; - - Godunov_trans_xbc( - i, j, k, n, q, lo, hi, uad, bc.lo(0), bc.hi(0), dlo.x, dhi.x); - xlo(i, j, k, n) = lo; - xhi(i, j, k, n) = hi; - Real st = (uval) ? lo : hi; - Imx(i, j, k, n) = fux * st + (1. - fux) * 0.5 * (hi + lo); - }, - yebox, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - constexpr Real small_vel = 1.e-8; - - Real vad = vmac(i, j, k); - Real fuy = (std::abs(vad) < small_vel) ? 0. : 1.; - bool vval = vad >= 0.; - // divu = 0 - // Real cons1 = (iconserv[n]) ? -0.5*l_dt*q(i,j-1,k,n)*divu(i,j-1,k) - // : 0.; Real cons2 = (iconserv[n]) ? -0.5*l_dt*q(i,j ,k,n)*divu(i,j - // ,k) : 0.; - Real lo = Ipy(i, j - 1, k, n); // + cons1; - Real hi = Imy(i, j, k, n); // + cons2; - if (use_forces_in_trans && fq) { - lo += 0.5 * l_dt * fq(i, j - 1, k, n); - hi += 0.5 * l_dt * fq(i, j, k, n); - } - - auto bc = pbc[n]; - - Godunov_trans_ybc( - i, j, k, n, q, lo, hi, vad, bc.lo(1), bc.hi(1), dlo.y, dhi.y); - - ylo(i, j, k, n) = lo; - yhi(i, j, k, n) = hi; - Real st = (vval) ? lo : hi; - Imy(i, j, k, n) = fuy * st + (1. - fuy) * 0.5 * (hi + lo); - }, - zebox, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - constexpr Real small_vel = 1.e-8; - - Real wad = wmac(i, j, k); - Real fuz = (std::abs(wad) < small_vel) ? 0. : 1.; - bool wval = wad >= 0.; - auto bc = pbc[n]; - // divu = 0 - // Real cons1 = (iconserv[n]) ? -0.5*l_dt*q(i,j,k-1,n)*divu(i,j,k-1) - // : 0.; Real cons2 = (iconserv[n]) ? -0.5*l_dt*q(i,j,k - // ,n)*divu(i,j,k ) : 0.; - Real lo = Ipz(i, j, k - 1, n); // + cons1; - Real hi = Imz(i, j, k, n); // + cons2; - if (use_forces_in_trans && fq) { - lo += 0.5 * l_dt * fq(i, j, k - 1, n); - hi += 0.5 * l_dt * fq(i, j, k, n); - } - - Godunov_trans_zbc( - i, j, k, n, q, lo, hi, wad, bc.lo(2), bc.hi(2), dlo.z, dhi.z); - - zlo(i, j, k, n) = lo; - zhi(i, j, k, n) = hi; - Real st = (wval) ? lo : hi; - Imz(i, j, k, n) = fuz * st + (1. - fuz) * 0.5 * (hi + lo); - }); - - Array4 xedge = Imx; - Array4 yedge = Imy; - Array4 zedge = Imz; - - // We can reuse the space in Ipx, Ipy and Ipz. - - // - // x-direction - // - Box const& xbxtmp = amrex::grow(bx, 0, 1); - Array4 yzlo = - makeArray4(xyzlo.dataPtr(), amrex::surroundingNodes(xbxtmp, 1), ncomp); - Array4 zylo = - makeArray4(xyzhi.dataPtr(), amrex::surroundingNodes(xbxtmp, 2), ncomp); - amrex::ParallelFor( - Box(zylo), ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; - Real l_zylo, l_zyhi; - Godunov_corner_couple_zy( - l_zylo, l_zyhi, i, j, k, n, l_dt, dy, iconserv[n] != 0, - zlo(i, j, k, n), zhi(i, j, k, n), q, vmac, yedge); - - Real wad = wmac(i, j, k); - Godunov_trans_zbc( - i, j, k, n, q, l_zylo, l_zyhi, wad, bc.lo(2), bc.hi(2), dlo.z, - dhi.z); - - constexpr Real small_vel = 1.e-8; - - Real st = (wad >= 0.) ? l_zylo : l_zyhi; - Real fu = (std::abs(wad) < small_vel) ? 0.0 : 1.0; - zylo(i, j, k, n) = fu * st + (1.0 - fu) * 0.5 * (l_zyhi + l_zylo); - }, - Box(yzlo), ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; - Real l_yzlo, l_yzhi; - Godunov_corner_couple_yz( - l_yzlo, l_yzhi, i, j, k, n, l_dt, dz, iconserv[n] != 0, - ylo(i, j, k, n), yhi(i, j, k, n), q, wmac, zedge); - - Real vad = vmac(i, j, k); - Godunov_trans_ybc( - i, j, k, n, q, l_yzlo, l_yzhi, vad, bc.lo(1), bc.hi(1), dlo.y, - dhi.y); - - constexpr Real small_vel = 1.e-8; - - Real st = (vad >= 0.) ? l_yzlo : l_yzhi; - Real fu = (std::abs(vad) < small_vel) ? 0.0 : 1.0; - yzlo(i, j, k, n) = fu * st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); - }); - // - - amrex::ParallelFor( - xbx, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Real stl, sth; - constexpr Real small_vel = 1.e-8; - if (iconserv[n] != 0) { - stl = xlo(i, j, k, n) - - (0.5 * dtdy) * - (yzlo(i - 1, j + 1, k, n) * vmac(i - 1, j + 1, k) - - yzlo(i - 1, j, k, n) * vmac(i - 1, j, k)) - - (0.5 * dtdz) * - (zylo(i - 1, j, k + 1, n) * wmac(i - 1, j, k + 1) - - zylo(i - 1, j, k, n) * wmac(i - 1, j, k)) + - (0.5 * dtdy) * q(i - 1, j, k, n) * - (vmac(i - 1, j + 1, k) - vmac(i - 1, j, k)) + - (0.5 * dtdz) * q(i - 1, j, k, n) * - (wmac(i - 1, j, k + 1) - wmac(i - 1, j, k)); - - sth = xhi(i, j, k, n) - - (0.5 * dtdy) * (yzlo(i, j + 1, k, n) * vmac(i, j + 1, k) - - yzlo(i, j, k, n) * vmac(i, j, k)) - - (0.5 * dtdz) * (zylo(i, j, k + 1, n) * wmac(i, j, k + 1) - - zylo(i, j, k, n) * wmac(i, j, k)) + - (0.5 * dtdy) * q(i, j, k, n) * - (vmac(i, j + 1, k) - vmac(i, j, k)) + - (0.5 * dtdz) * q(i, j, k, n) * - (wmac(i, j, k + 1) - wmac(i, j, k)); - - } else { - stl = xlo(i, j, k, n) - - (0.25 * dtdy) * - (vmac(i - 1, j + 1, k) + vmac(i - 1, j, k)) * - (yzlo(i - 1, j + 1, k, n) - yzlo(i - 1, j, k, n)) - - (0.25 * dtdz) * - (wmac(i - 1, j, k + 1) + wmac(i - 1, j, k)) * - (zylo(i - 1, j, k + 1, n) - zylo(i - 1, j, k, n)); - - sth = xhi(i, j, k, n) - - (0.25 * dtdy) * (vmac(i, j + 1, k) + vmac(i, j, k)) * - (yzlo(i, j + 1, k, n) - yzlo(i, j, k, n)) - - (0.25 * dtdz) * (wmac(i, j, k + 1) + wmac(i, j, k)) * - (zylo(i, j, k + 1, n) - zylo(i, j, k, n)); - } - - stl += (!use_forces_in_trans && fq) - ? 0.5 * l_dt * fq(i - 1, j, k, n) - : 0.; - sth += - (!use_forces_in_trans && fq) ? 0.5 * l_dt * fq(i, j, k, n) : 0.; - - auto bc = pbc[n]; - Godunov_cc_xbc_lo(i, j, k, n, q, stl, sth, umac, bc.lo(0), dlo.x); - Godunov_cc_xbc_hi(i, j, k, n, q, stl, sth, umac, bc.hi(0), dhi.x); - - Real qx = (umac(i, j, k) >= 0.) ? stl : sth; - qx = (std::abs(umac(i, j, k)) < small_vel) ? 0.5 * (stl + sth) : qx; - - if (iconserv[n] == 1) { - fx(i, j, k, n) = umac(i, j, k) * qx; - } else { - fx(i, j, k, n) = qx; - } - }); - - // - // y-direction - // - Box const& ybxtmp = amrex::grow(bx, 1, 1); - Array4 xzlo = - makeArray4(xyzlo.dataPtr(), amrex::surroundingNodes(ybxtmp, 0), ncomp); - Array4 zxlo = - makeArray4(xyzhi.dataPtr(), amrex::surroundingNodes(ybxtmp, 2), ncomp); - amrex::ParallelFor( - Box(xzlo), ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; - Real l_xzlo, l_xzhi; - Godunov_corner_couple_xz( - l_xzlo, l_xzhi, i, j, k, n, l_dt, dz, iconserv[n] != 0, - xlo(i, j, k, n), xhi(i, j, k, n), q, wmac, zedge); - - Real uad = umac(i, j, k); - Godunov_trans_xbc( - i, j, k, n, q, l_xzlo, l_xzhi, uad, bc.lo(0), bc.hi(0), dlo.x, - dhi.x); - - constexpr Real small_vel = 1.e-8; - - Real st = (uad >= 0.) ? l_xzlo : l_xzhi; - Real fu = (std::abs(uad) < small_vel) ? 0.0 : 1.0; - xzlo(i, j, k, n) = fu * st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); - }, - Box(zxlo), ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; - Real l_zxlo, l_zxhi; - Godunov_corner_couple_zx( - l_zxlo, l_zxhi, i, j, k, n, l_dt, dx, iconserv[n] != 0, - zlo(i, j, k, n), zhi(i, j, k, n), q, umac, xedge); - - Real wad = wmac(i, j, k); - Godunov_trans_zbc( - i, j, k, n, q, l_zxlo, l_zxhi, wad, bc.lo(2), bc.hi(2), dlo.z, - dhi.z); - - constexpr Real small_vel = 1.e-8; - - Real st = (wad >= 0.) ? l_zxlo : l_zxhi; - Real fu = (std::abs(wad) < small_vel) ? 0.0 : 1.0; - zxlo(i, j, k, n) = fu * st + (1.0 - fu) * 0.5 * (l_zxhi + l_zxlo); - }); - - amrex::ParallelFor( - ybx, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Real stl, sth; - constexpr Real small_vel = 1.e-8; - - if (iconserv[n] != 0) { - stl = ylo(i, j, k, n) - - (0.5 * dtdx) * - (xzlo(i + 1, j - 1, k, n) * umac(i + 1, j - 1, k) - - xzlo(i, j - 1, k, n) * umac(i, j - 1, k)) - - (0.5 * dtdz) * - (zxlo(i, j - 1, k + 1, n) * wmac(i, j - 1, k + 1) - - zxlo(i, j - 1, k, n) * wmac(i, j - 1, k)) + - (0.5 * dtdx) * q(i, j - 1, k, n) * - (umac(i + 1, j - 1, k) - umac(i, j - 1, k)) + - (0.5 * dtdz) * q(i, j - 1, k, n) * - (wmac(i, j - 1, k + 1) - wmac(i, j - 1, k)); - - sth = yhi(i, j, k, n) - - (0.5 * dtdx) * (xzlo(i + 1, j, k, n) * umac(i + 1, j, k) - - xzlo(i, j, k, n) * umac(i, j, k)) - - (0.5 * dtdz) * (zxlo(i, j, k + 1, n) * wmac(i, j, k + 1) - - zxlo(i, j, k, n) * wmac(i, j, k)) + - (0.5 * dtdx) * q(i, j, k, n) * - (umac(i + 1, j, k) - umac(i, j, k)) + - (0.5 * dtdz) * q(i, j, k, n) * - (wmac(i, j, k + 1) - wmac(i, j, k)); - } else { - stl = ylo(i, j, k, n) - - (0.25 * dtdx) * - (umac(i + 1, j - 1, k) + umac(i, j - 1, k)) * - (xzlo(i + 1, j - 1, k, n) - xzlo(i, j - 1, k, n)) - - (0.25 * dtdz) * - (wmac(i, j - 1, k + 1) + wmac(i, j - 1, k)) * - (zxlo(i, j - 1, k + 1, n) - zxlo(i, j - 1, k, n)); - - sth = yhi(i, j, k, n) - - (0.25 * dtdx) * (umac(i + 1, j, k) + umac(i, j, k)) * - (xzlo(i + 1, j, k, n) - xzlo(i, j, k, n)) - - (0.25 * dtdz) * (wmac(i, j, k + 1) + wmac(i, j, k)) * - (zxlo(i, j, k + 1, n) - zxlo(i, j, k, n)); - } - - stl += (!use_forces_in_trans && fq) - ? 0.5 * l_dt * fq(i, j - 1, k, n) - : 0.; - sth += - (!use_forces_in_trans && fq) ? 0.5 * l_dt * fq(i, j, k, n) : 0.; - - auto bc = pbc[n]; - Godunov_cc_ybc_lo(i, j, k, n, q, stl, sth, vmac, bc.lo(1), dlo.y); - Godunov_cc_ybc_hi(i, j, k, n, q, stl, sth, vmac, bc.hi(1), dhi.y); - - Real qy = (vmac(i, j, k) >= 0.) ? stl : sth; - qy = (std::abs(vmac(i, j, k)) < small_vel) ? 0.5 * (stl + sth) : qy; - - if (iconserv[n] == 1) { - fy(i, j, k, n) = vmac(i, j, k) * qy; - } else { - fy(i, j, k, n) = qy; - } - }); - - // - // z-direction - // - Box const& zbxtmp = amrex::grow(bx, 2, 1); - Array4 xylo = - makeArray4(xyzlo.dataPtr(), amrex::surroundingNodes(zbxtmp, 0), ncomp); - Array4 yxlo = - makeArray4(xyzhi.dataPtr(), amrex::surroundingNodes(zbxtmp, 1), ncomp); - amrex::ParallelFor( - Box(xylo), ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; - Real l_xylo, l_xyhi; - Godunov_corner_couple_xy( - l_xylo, l_xyhi, i, j, k, n, l_dt, dy, iconserv[n] != 0, - xlo(i, j, k, n), xhi(i, j, k, n), q, vmac, yedge); - - Real uad = umac(i, j, k); - Godunov_trans_xbc( - i, j, k, n, q, l_xylo, l_xyhi, uad, bc.lo(0), bc.hi(0), dlo.x, - dhi.x); - - constexpr Real small_vel = 1.e-8; - - Real st = (uad >= 0.) ? l_xylo : l_xyhi; - Real fu = (std::abs(uad) < small_vel) ? 0.0 : 1.0; - xylo(i, j, k, n) = fu * st + (1.0 - fu) * 0.5 * (l_xyhi + l_xylo); - }, - Box(yxlo), ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; - Real l_yxlo, l_yxhi; - Godunov_corner_couple_yx( - l_yxlo, l_yxhi, i, j, k, n, l_dt, dx, iconserv[n] != 0, - ylo(i, j, k, n), yhi(i, j, k, n), q, umac, xedge); - - Real vad = vmac(i, j, k); - Godunov_trans_ybc( - i, j, k, n, q, l_yxlo, l_yxhi, vad, bc.lo(1), bc.hi(1), dlo.y, - dhi.y); - - constexpr Real small_vel = 1.e-8; - - Real st = (vad >= 0.) ? l_yxlo : l_yxhi; - Real fu = (std::abs(vad) < small_vel) ? 0.0 : 1.0; - yxlo(i, j, k, n) = fu * st + (1.0 - fu) * 0.5 * (l_yxhi + l_yxlo); - }); - - amrex::ParallelFor( - zbx, ncomp, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Real stl, sth; - constexpr Real small_vel = 1.e-8; - if (iconserv[n] != 0) { - stl = zlo(i, j, k, n) - - (0.5 * dtdx) * - (xylo(i + 1, j, k - 1, n) * umac(i + 1, j, k - 1) - - xylo(i, j, k - 1, n) * umac(i, j, k - 1)) - - (0.5 * dtdy) * - (yxlo(i, j + 1, k - 1, n) * vmac(i, j + 1, k - 1) - - yxlo(i, j, k - 1, n) * vmac(i, j, k - 1)) + - (0.5 * dtdx) * q(i, j, k - 1, n) * - (umac(i + 1, j, k - 1) - umac(i, j, k - 1)) + - (0.5 * dtdy) * q(i, j, k - 1, n) * - (vmac(i, j + 1, k - 1) - vmac(i, j, k - 1)); - - sth = zhi(i, j, k, n) - - (0.5 * dtdx) * (xylo(i + 1, j, k, n) * umac(i + 1, j, k) - - xylo(i, j, k, n) * umac(i, j, k)) - - (0.5 * dtdy) * (yxlo(i, j + 1, k, n) * vmac(i, j + 1, k) - - yxlo(i, j, k, n) * vmac(i, j, k)) + - (0.5 * dtdx) * q(i, j, k, n) * - (umac(i + 1, j, k) - umac(i, j, k)) + - (0.5 * dtdy) * q(i, j, k, n) * - (vmac(i, j + 1, k) - vmac(i, j, k)); - } else { - stl = zlo(i, j, k, n) - - (0.25 * dtdx) * - (umac(i + 1, j, k - 1) + umac(i, j, k - 1)) * - (xylo(i + 1, j, k - 1, n) - xylo(i, j, k - 1, n)) - - (0.25 * dtdy) * - (vmac(i, j + 1, k - 1) + vmac(i, j, k - 1)) * - (yxlo(i, j + 1, k - 1, n) - yxlo(i, j, k - 1, n)); - - sth = zhi(i, j, k, n) - - (0.25 * dtdx) * (umac(i + 1, j, k) + umac(i, j, k)) * - (xylo(i + 1, j, k, n) - xylo(i, j, k, n)) - - (0.25 * dtdy) * (vmac(i, j + 1, k) + vmac(i, j, k)) * - (yxlo(i, j + 1, k, n) - yxlo(i, j, k, n)); - } - - stl += (!use_forces_in_trans && fq) - ? 0.5 * l_dt * fq(i, j, k - 1, n) - : 0.; - sth += - (!use_forces_in_trans && fq) ? 0.5 * l_dt * fq(i, j, k, n) : 0.; - - auto bc = pbc[n]; - Godunov_cc_zbc_lo(i, j, k, n, q, stl, sth, wmac, bc.lo(2), dlo.z); - Godunov_cc_zbc_hi(i, j, k, n, q, stl, sth, wmac, bc.hi(2), dhi.z); - - Real qz = (wmac(i, j, k) >= 0.) ? stl : sth; - qz = (std::abs(wmac(i, j, k)) < small_vel) ? 0.5 * (stl + sth) : qz; - - if (iconserv[n] == 1) { - fz(i, j, k, n) = wmac(i, j, k) * qz; - } else { - fz(i, j, k, n) = qz; - } - }); -} diff --git a/amr-wind/convection/incflo_godunov_minmod.H b/amr-wind/convection/incflo_godunov_minmod.H deleted file mode 100644 index 499d5ebf03..0000000000 --- a/amr-wind/convection/incflo_godunov_minmod.H +++ /dev/null @@ -1,207 +0,0 @@ -#ifndef GODUNOV_MINMOD_H -#define GODUNOV_MINMOD_H - -#include -#include - -/* This header file contains the inlined __host__ __device__ functions required - for the scalar advection routines for 3D Godunov. It also contains function - declarations for controlling host functions. */ - -namespace { - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void minmod( - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, - amrex::Real& dsm, - amrex::Real& dsp) -{ - // Calculate gradients on both sides - dsp = sp1 - s0; - dsm = s0 - sm1; - - if (!(dsp * dsm < 0.0)) { - // Select the smaller slope if same sign - if (std::abs(dsp) < std::abs(dsm)) { - dsm = dsp; - } else { - dsp = dsm; - } - } else { - // Set to zero if opposite sign - dsp = 0.0; - dsm = 0.0; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_minmod_bc( - const int n, - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, - amrex::Real& dsm, - amrex::Real& dsp, - const int bclo, - const int bchi, - const int domlo, - const int domhi) -{ - using namespace amrex; - - if (bclo == BCType::ext_dir || bclo == BCType::hoextrap) { - if (n == domlo) { - // Ensure that left-side slope is used unchanged - dsm = s0 - sm1; - } - } - - if (bchi == BCType::ext_dir || bchi == BCType::hoextrap) { - if (n == domhi) { - // Ensure that the right-side slope is used unchanged - dsp = sp1 - s0; - } - } -} - -// This version is called after the MAC projection, when we use the -// MAC-projected velocity -// for upwinding -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_minmod_fpu_x( - const int i, - const int j, - const int k, - const int n, - const amrex::Real dt, - const amrex::Real dx, - amrex::Real& Im, - amrex::Real& Ip, - const amrex::Array4& S, - const amrex::Array4& vel_edge, - const amrex::BCRec& bc, - const int domlo, - const int domhi) -{ - - using namespace amrex; - - constexpr amrex::Real small_vel = 1e-8; - - amrex::Real sm1 = S(i - 1, j, k, n); - amrex::Real s0 = S(i, j, k, n); - amrex::Real sp1 = S(i + 1, j, k, n); - amrex::Real dsp = 0.0; - amrex::Real dsm = 0.0; - minmod(sm1, s0, sp1, dsm, dsp); - Godunov_minmod_bc( - i, sm1, s0, sp1, dsm, dsp, bc.lo(0), bc.hi(0), domlo, domhi); - - amrex::Real sigmap = std::abs(vel_edge(i + 1, j, k)) * dt / dx; - amrex::Real sigmam = std::abs(vel_edge(i, j, k)) * dt / dx; - - if (vel_edge(i + 1, j, k) > small_vel) { - Ip = s0 + 0.5 * (1.0 - sigmap) * dsp; - } else { - Ip = S(i, j, k, n); - } - - if (vel_edge(i, j, k) < -small_vel) { - Im = s0 - 0.5 * (1.0 - sigmam) * dsm; - } else { - Im = S(i, j, k, n); - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_minmod_fpu_y( - const int i, - const int j, - const int k, - const int n, - const amrex::Real dt, - const amrex::Real dx, - amrex::Real& Im, - amrex::Real& Ip, - const amrex::Array4& S, - const amrex::Array4& vel_edge, - const amrex::BCRec& bc, - const int domlo, - const int domhi) -{ - - using namespace amrex; - - constexpr amrex::Real small_vel = 1e-8; - - amrex::Real sm1 = S(i, j - 1, k, n); - amrex::Real s0 = S(i, j, k, n); - amrex::Real sp1 = S(i, j + 1, k, n); - amrex::Real dsp = 0.0; - amrex::Real dsm = 0.0; - minmod(sm1, s0, sp1, dsm, dsp); - Godunov_minmod_bc( - j, sm1, s0, sp1, dsm, dsp, bc.lo(1), bc.hi(1), domlo, domhi); - - amrex::Real sigmap = std::abs(vel_edge(i, j + 1, k)) * dt / dx; - amrex::Real sigmam = std::abs(vel_edge(i, j, k)) * dt / dx; - - if (vel_edge(i, j + 1, k) > small_vel) { - Ip = s0 + 0.5 * (1.0 - sigmap) * dsp; - } else { - Ip = S(i, j, k, n); - } - - if (vel_edge(i, j, k) < -small_vel) { - Im = s0 - 0.5 * (1.0 - sigmam) * dsm; - } else { - Im = S(i, j, k, n); - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_minmod_fpu_z( - const int i, - const int j, - const int k, - const int n, - const amrex::Real dt, - const amrex::Real dx, - amrex::Real& Im, - amrex::Real& Ip, - const amrex::Array4& S, - const amrex::Array4& vel_edge, - const amrex::BCRec& bc, - const int domlo, - const int domhi) -{ - - using namespace amrex; - - constexpr amrex::Real small_vel = 1e-8; - - amrex::Real sm1 = S(i, j, k - 1, n); - amrex::Real s0 = S(i, j, k, n); - amrex::Real sp1 = S(i, j, k + 1, n); - amrex::Real dsp = 0.0; - amrex::Real dsm = 0.0; - minmod(sm1, s0, sp1, dsm, dsp); - Godunov_minmod_bc( - k, sm1, s0, sp1, dsm, dsp, bc.lo(2), bc.hi(2), domlo, domhi); - - amrex::Real sigmap = std::abs(vel_edge(i, j, k + 1)) * dt / dx; - amrex::Real sigmam = std::abs(vel_edge(i, j, k)) * dt / dx; - - if (vel_edge(i, j, k + 1) > small_vel) { - Ip = s0 + 0.5 * (1.0 - sigmap) * dsp; - } else { - Ip = S(i, j, k, n); - } - - if (vel_edge(i, j, k) < -small_vel) { - Im = s0 - 0.5 * (1.0 - sigmam) * dsm; - } else { - Im = S(i, j, k, n); - } -} - -} // namespace - -#endif diff --git a/amr-wind/convection/incflo_godunov_ppm.H b/amr-wind/convection/incflo_godunov_ppm.H deleted file mode 100644 index 87f0c9e743..0000000000 --- a/amr-wind/convection/incflo_godunov_ppm.H +++ /dev/null @@ -1,575 +0,0 @@ -#ifndef GODUNOV_PPM_H -#define GODUNOV_PPM_H - -#include -#include - -/* This header file contains the inlined __host__ __device__ functions required - for the scalar advection routines for 3D Godunov. It also contains function - declarations for controlling host functions. */ - -namespace { - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_trans_xbc( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - amrex::Real& /* uedge */, - const int bclo, - const int bchi, - const int domlo, - const int domhi) -{ - using namespace amrex; - - // Low X - if (i <= domlo) { - if (bclo == BCType::ext_dir) { - // IAMR does this but it breaks lo/hi symmetry - // Real st = (uedge <= small_vel) ? hi : s(domlo-1,j,k,n); - // So here we do something simpler... - Real st = s(domlo - 1, j, k, n); - lo = st; - hi = st; - } - - else if ( - bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even) { - lo = hi; - - } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = 0.; - } - } - - // High X - if (i > domhi) { - if (bchi == BCType::ext_dir) { - // IAMR does this but it breaks lo/hi symmetry - // Real st = (uedge >= -small_vel)? lo : s(domhi+1,j,k,n); - // So here we do something simpler... - Real st = s(domhi + 1, j, k, n); - lo = st; - hi = st; - } - - else if ( - bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even) { - hi = lo; - - } else if (bchi == BCType::reflect_odd) { - hi = 0.; - lo = 0.; - } - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_trans_ybc( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - amrex::Real /* vedge */, - const int bclo, - const int bchi, - const int domlo, - const int domhi) -{ - using namespace amrex; - - // Low Y - if (j <= domlo) { - if (bclo == BCType::ext_dir) { - // IAMR does this but it breaks lo/hi symmetry - // Real st = (vedge <= small_vel) ? hi : s(i,domlo-1,k,n); - // So here we do something simpler... - Real st = s(i, domlo - 1, k, n); - lo = st; - hi = st; - } - - else if ( - bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even) { - lo = hi; - - } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = 0.; - } - } - - // High Y - if (j > domhi) { - if (bchi == BCType::ext_dir) { - // IAMR does this but it breaks lo/hi symmetry - // Real st = (vedge >= -small_vel)? lo : s(i,domhi+1,k,n); - // So here we do something simpler... - Real st = s(i, domhi + 1, k, n); - lo = st; - hi = st; - } - - else if ( - bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even) { - hi = lo; - - } else if (bchi == BCType::reflect_odd) { - hi = 0.; - lo = 0.; - } - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_trans_zbc( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - amrex::Real /* wedge */, - const int bclo, - const int bchi, - const int domlo, - const int domhi) -{ - using namespace amrex; - - // Low Z - if (k <= domlo) { - if (bclo == BCType::ext_dir) { - // IAMR does this but it breaks lo/hi symmetry - // Real st = (wedge <= small_vel) ? hi : s(i,j,domlo-1,n); - // So here we do something simpler... - Real st = s(i, j, domlo - 1, n); - lo = st; - hi = st; - } - - else if ( - bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even) { - lo = hi; - - } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = 0.; - } - } - - // High Z - if (k > domhi) { - if (bchi == BCType::ext_dir) { - // IAMR does this but it breaks lo/hi symmetry - // Real st = (wedge >= -small_vel)? lo : s(i,j,domhi+1,n); - // So here we do something simpler... - Real st = s(i, j, domhi + 1, n); - lo = st; - hi = st; - } else if ( - bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even) { - hi = lo; - - } else if (bchi == BCType::reflect_odd) { - hi = 0.; - lo = 0.; - } - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_cc_xbc_lo( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - const amrex::Array4& uedge, - const int bclo, - const int domlo) -{ - using namespace amrex; - - if (i == domlo) { - if (bclo == BCType::ext_dir && uedge(i, j, k) >= 0.) { - hi = s(domlo - 1, j, k, n); - lo = hi; - } else if ( - bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even || - (bclo == BCType::ext_dir && uedge(i, j, k) < 0)) { - lo = hi; - - } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = hi; - } - } else { - return; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_cc_xbc_hi( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - const amrex::Array4& uedge, - const int bchi, - const int domhi) -{ - using namespace amrex; - - if (i == domhi + 1) { - if (bchi == BCType::ext_dir && uedge(i, j, k) <= 0.) { - lo = s(domhi + 1, j, k, n); - hi = lo; - } else if ( - bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even || - (bchi == BCType::ext_dir && uedge(i, j, k) > 0)) { - hi = lo; - - } else if (bchi == BCType::reflect_odd) { - lo = 0.; - hi = lo; - } - } else { - return; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_cc_ybc_lo( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - const amrex::Array4& vedge, - const int bclo, - const int domlo) -{ - using namespace amrex; - - if (j == domlo) { - if (bclo == BCType::ext_dir && vedge(i, j, k) >= 0.) { - hi = s(i, domlo - 1, k, n); - lo = hi; - } else if ( - bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even || - (bclo == BCType::ext_dir && vedge(i, j, k) < 0)) { - lo = hi; - - } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = hi; - } - } else { - return; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_cc_ybc_hi( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - const amrex::Array4& vedge, - const int bchi, - const int domhi) -{ - using namespace amrex; - - if (j == domhi + 1) { - if (bchi == BCType::ext_dir && vedge(i, j, k) <= 0.) { - lo = s(i, domhi + 1, k, n); - hi = lo; - } else if ( - bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even || - (bchi == BCType::ext_dir && vedge(i, j, k) > 0)) { - hi = lo; - - } else if (bchi == BCType::reflect_odd) { - lo = 0.; - hi = lo; - } - } else { - return; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_cc_zbc_lo( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - const amrex::Array4& wedge, - const int bclo, - const int domlo) -{ - using namespace amrex; - - if (k == domlo) { - if (bclo == BCType::ext_dir && wedge(i, j, k) >= 0.) { - hi = s(i, j, domlo - 1, n); - lo = hi; - } - - else if ( - bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even || - (bclo == BCType::ext_dir && wedge(i, j, k) < 0)) { - lo = hi; - - } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = hi; - } - } else { - return; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_cc_zbc_hi( - const int i, - const int j, - const int k, - const int n, - const amrex::Array4& s, - amrex::Real& lo, - amrex::Real& hi, - const amrex::Array4& wedge, - const int bchi, - const int domhi) -{ - using namespace amrex; - - if (k == domhi + 1) { - if (bchi == BCType::ext_dir && wedge(i, j, k) <= 0.) { - lo = s(i, j, domhi + 1, n); - hi = lo; - } else if ( - bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even || - (bchi == BCType::ext_dir && wedge(i, j, k) > 0)) { - hi = lo; - - } else if (bchi == BCType::reflect_odd) { - lo = 0.; - hi = lo; - } - } else { - return; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_corner_couple_yx( - amrex::Real& lo1, - amrex::Real& hi1, - int i, - int j, - int k, - int n, - amrex::Real dt, - amrex::Real dx, - bool /*iconserv*/, - amrex::Real lo, - amrex::Real hi, - amrex::Array4 const& s, - amrex::Array4 const& mac, - amrex::Array4 const& state) -{ - lo1 = lo - - dt / (amrex::Real(3.0) * dx) * - (state(i + 1, j - 1, k, n) * mac(i + 1, j - 1, k) - - state(i, j - 1, k, n) * mac(i, j - 1, k)) + - dt / (amrex::Real(3.0) * dx) * s(i, j - 1, k, n) * - (mac(i + 1, j - 1, k) - mac(i, j - 1, k)); - hi1 = hi - - dt / (amrex::Real(3.0) * dx) * - (state(i + 1, j, k, n) * mac(i + 1, j, k) - - state(i, j, k, n) * mac(i, j, k)) + - dt / (amrex::Real(3.0) * dx) * s(i, j, k, n) * - (mac(i + 1, j, k) - mac(i, j, k)); -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_corner_couple_zx( - amrex::Real& lo1, - amrex::Real& hi1, - int i, - int j, - int k, - int n, - amrex::Real dt, - amrex::Real dx, - bool /*iconserv*/, - amrex::Real lo, - amrex::Real hi, - amrex::Array4 const& s, - amrex::Array4 const& mac, - amrex::Array4 const& state) -{ - lo1 = lo - - dt / (amrex::Real(3.0) * dx) * - (state(i + 1, j, k - 1, n) * mac(i + 1, j, k - 1) - - state(i, j, k - 1, n) * mac(i, j, k - 1)) + - dt / (amrex::Real(3.0) * dx) * s(i, j, k - 1, n) * - (mac(i + 1, j, k - 1) - mac(i, j, k - 1)); - hi1 = hi - - dt / (amrex::Real(3.0) * dx) * - (state(i + 1, j, k, n) * mac(i + 1, j, k) - - state(i, j, k, n) * mac(i, j, k)) + - dt / (amrex::Real(3.0) * dx) * s(i, j, k, n) * - (mac(i + 1, j, k) - mac(i, j, k)); -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_corner_couple_xy( - amrex::Real& lo1, - amrex::Real& hi1, - int i, - int j, - int k, - int n, - amrex::Real dt, - amrex::Real dy, - bool /*iconserv*/, - amrex::Real lo, - amrex::Real hi, - amrex::Array4 const& s, - amrex::Array4 const& mac, - amrex::Array4 const& state) -{ - lo1 = lo - - dt / (amrex::Real(3.0) * dy) * - (state(i - 1, j + 1, k, n) * mac(i - 1, j + 1, k) - - state(i - 1, j, k, n) * mac(i - 1, j, k)) + - dt / (amrex::Real(3.0) * dy) * s(i - 1, j, k, n) * - (mac(i - 1, j + 1, k) - mac(i - 1, j, k)); - hi1 = hi - - dt / (amrex::Real(3.0) * dy) * - (state(i, j + 1, k, n) * mac(i, j + 1, k) - - state(i, j, k, n) * mac(i, j, k)) + - dt / (amrex::Real(3.0) * dy) * s(i, j, k, n) * - (mac(i, j + 1, k) - mac(i, j, k)); -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_corner_couple_zy( - amrex::Real& lo1, - amrex::Real& hi1, - int i, - int j, - int k, - int n, - amrex::Real dt, - amrex::Real dy, - bool /*iconserv*/, - amrex::Real lo, - amrex::Real hi, - amrex::Array4 const& s, - amrex::Array4 const& mac, - amrex::Array4 const& state) -{ - lo1 = lo - - dt / (amrex::Real(3.0) * dy) * - (state(i, j + 1, k - 1, n) * mac(i, j + 1, k - 1) - - state(i, j, k - 1, n) * mac(i, j, k - 1)) + - dt / (amrex::Real(3.0) * dy) * s(i, j, k - 1, n) * - (mac(i, j + 1, k - 1) - mac(i, j, k - 1)); - hi1 = hi - - dt / (amrex::Real(3.0) * dy) * - (state(i, j + 1, k, n) * mac(i, j + 1, k) - - state(i, j, k, n) * mac(i, j, k)) + - dt / (amrex::Real(3.0) * dy) * s(i, j, k, n) * - (mac(i, j + 1, k) - mac(i, j, k)); -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_corner_couple_xz( - amrex::Real& lo1, - amrex::Real& hi1, - int i, - int j, - int k, - int n, - amrex::Real dt, - amrex::Real dz, - bool /*iconserv*/, - amrex::Real lo, - amrex::Real hi, - amrex::Array4 const& s, - amrex::Array4 const& mac, - amrex::Array4 const& state) -{ - lo1 = lo - - dt / (amrex::Real(3.0) * dz) * - (state(i - 1, j, k + 1, n) * mac(i - 1, j, k + 1) - - state(i - 1, j, k, n) * mac(i - 1, j, k)) + - dt / (amrex::Real(3.0) * dz) * s(i - 1, j, k, n) * - (mac(i - 1, j, k + 1) - mac(i - 1, j, k)); - hi1 = hi - - dt / (amrex::Real(3.0) * dz) * - (state(i, j, k + 1, n) * mac(i, j, k + 1) - - state(i, j, k, n) * mac(i, j, k)) + - dt / (amrex::Real(3.0) * dz) * s(i, j, k, n) * - (mac(i, j, k + 1) - mac(i, j, k)); -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_corner_couple_yz( - amrex::Real& lo1, - amrex::Real& hi1, - int i, - int j, - int k, - int n, - amrex::Real dt, - amrex::Real dz, - bool /*iconserv*/, - amrex::Real lo, - amrex::Real hi, - amrex::Array4 const& s, - amrex::Array4 const& mac, - amrex::Array4 const& state) -{ - lo1 = lo - - dt / (amrex::Real(3.0) * dz) * - (state(i, j - 1, k + 1, n) * mac(i, j - 1, k + 1) - - state(i, j - 1, k, n) * mac(i, j - 1, k)) + - dt / (amrex::Real(3.0) * dz) * s(i, j - 1, k, n) * - (mac(i, j - 1, k + 1) - mac(i, j - 1, k)); - hi1 = hi - - dt / (amrex::Real(3.0) * dz) * - (state(i, j, k + 1, n) * mac(i, j, k + 1) - - state(i, j, k, n) * mac(i, j, k)) + - dt / (amrex::Real(3.0) * dz) * s(i, j, k, n) * - (mac(i, j, k + 1) - mac(i, j, k)); -} - -} // namespace -#endif diff --git a/amr-wind/convection/incflo_godunov_upwind.H b/amr-wind/convection/incflo_godunov_upwind.H deleted file mode 100644 index d9a768593d..0000000000 --- a/amr-wind/convection/incflo_godunov_upwind.H +++ /dev/null @@ -1,30 +0,0 @@ -#ifndef GODUNOV_UPWIND_H -#define GODUNOV_UPWIND_H - -#include -#include - -/* This header file contains the inlined __host__ __device__ functions required - for the scalar advection routines for 3D Godunov. It also contains function - declarations for controlling host functions. */ - -namespace { - -// This version is called after the MAC projection, when we use the -// MAC-projected velocity for upwinding -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_upwind_fpu( - const int i, - const int j, - const int k, - const int n, - amrex::Real& Im, - amrex::Real& Ip, - const amrex::Array4& S) -{ - Ip = S(i, j, k, n); - Im = S(i, j, k, n); -} - -} // namespace - -#endif diff --git a/amr-wind/convection/incflo_mol_fluxes.cpp b/amr-wind/convection/incflo_mol_fluxes.cpp index 787c090254..0ee44417b8 100644 --- a/amr-wind/convection/incflo_mol_fluxes.cpp +++ b/amr-wind/convection/incflo_mol_fluxes.cpp @@ -1,5 +1,5 @@ #include -#include "amr-wind/convection/incflo_convection_K.H" +#include #include "amr-wind/convection/MOL.H" #include "amr-wind/utilities/bc_ops.H" @@ -81,13 +81,13 @@ void mol::compute_convective_fluxes( } else { Real qpls = q(i, j, k, n) - - 0.5 * incflo_xslope_extdir( - i, j, k, n, q, extdir_or_ho_ilo, + 0.5 * amrex_calc_xslope_extdir( + i, j, k, n, 2, q, extdir_or_ho_ilo, extdir_or_ho_ihi, domain_ilo, domain_ihi); Real qmns = q(i - 1, j, k, n) + - 0.5 * incflo_xslope_extdir( - i - 1, j, k, n, q, extdir_or_ho_ilo, + 0.5 * amrex_calc_xslope_extdir( + i - 1, j, k, n, 2, q, extdir_or_ho_ilo, extdir_or_ho_ihi, domain_ilo, domain_ihi); if (umac(i, j, k) > small_vel) { qs = qmns; @@ -104,9 +104,10 @@ void mol::compute_convective_fluxes( xbx, ncomp, [q, umac, fx] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Real qpls = q(i, j, k, n) - 0.5 * incflo_xslope(i, j, k, n, q); - Real qmns = - q(i - 1, j, k, n) + 0.5 * incflo_xslope(i - 1, j, k, n, q); + Real qpls = + q(i, j, k, n) - 0.5 * amrex_calc_xslope(i, j, k, n, 2, q); + Real qmns = q(i - 1, j, k, n) + + 0.5 * amrex_calc_xslope(i - 1, j, k, n, 2, q); Real qs; if (umac(i, j, k) > small_vel) { qs = qmns; @@ -143,13 +144,13 @@ void mol::compute_convective_fluxes( } else { Real qpls = q(i, j, k, n) - - 0.5 * incflo_yslope_extdir( - i, j, k, n, q, extdir_or_ho_jlo, + 0.5 * amrex_calc_yslope_extdir( + i, j, k, n, 2, q, extdir_or_ho_jlo, extdir_or_ho_jhi, domain_jlo, domain_jhi); Real qmns = q(i, j - 1, k, n) + - 0.5 * incflo_yslope_extdir( - i, j - 1, k, n, q, extdir_or_ho_jlo, + 0.5 * amrex_calc_yslope_extdir( + i, j - 1, k, n, 2, q, extdir_or_ho_jlo, extdir_or_ho_jhi, domain_jlo, domain_jhi); if (vmac(i, j, k) > small_vel) { qs = qmns; @@ -166,9 +167,10 @@ void mol::compute_convective_fluxes( ybx, ncomp, [q, vmac, fy] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Real qpls = q(i, j, k, n) - 0.5 * incflo_yslope(i, j, k, n, q); - Real qmns = - q(i, j - 1, k, n) + 0.5 * incflo_yslope(i, j - 1, k, n, q); + Real qpls = + q(i, j, k, n) - 0.5 * amrex_calc_yslope(i, j, k, n, 2, q); + Real qmns = q(i, j - 1, k, n) + + 0.5 * amrex_calc_yslope(i, j - 1, k, n, 2, q); Real qs; if (vmac(i, j, k) > small_vel) { qs = qmns; @@ -205,13 +207,13 @@ void mol::compute_convective_fluxes( } else { Real qpls = q(i, j, k, n) - - 0.5 * incflo_zslope_extdir( - i, j, k, n, q, extdir_or_ho_klo, + 0.5 * amrex_calc_zslope_extdir( + i, j, k, n, 2, q, extdir_or_ho_klo, extdir_or_ho_khi, domain_klo, domain_khi); Real qmns = q(i, j, k - 1, n) + - 0.5 * incflo_zslope_extdir( - i, j, k - 1, n, q, extdir_or_ho_klo, + 0.5 * amrex_calc_zslope_extdir( + i, j, k - 1, n, 2, q, extdir_or_ho_klo, extdir_or_ho_khi, domain_klo, domain_khi); if (wmac(i, j, k) > small_vel) { qs = qmns; @@ -228,9 +230,10 @@ void mol::compute_convective_fluxes( zbx, ncomp, [q, wmac, fz] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - Real qpls = q(i, j, k, n) - 0.5 * incflo_zslope(i, j, k, n, q); - Real qmns = - q(i, j, k - 1, n) + 0.5 * incflo_zslope(i, j, k - 1, n, q); + Real qpls = + q(i, j, k, n) - 0.5 * amrex_calc_zslope(i, j, k, n, 2, q); + Real qmns = q(i, j, k - 1, n) + + 0.5 * amrex_calc_zslope(i, j, k - 1, n, 2, q); Real qs; if (wmac(i, j, k) > small_vel) { qs = qmns; diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index c1a1f4e21f..a72d6d24d0 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -187,11 +187,10 @@ struct AdvectionOp for (int lev = 0; lev < repo.num_active_levels(); ++lev) { HydroUtils::ExtrapVelToFaces( dof_field(lev), src_term(lev), u_mac(lev), v_mac(lev), - w_mac(lev), dof_field.bcrec(), - dof_field.bcrec_device().data(), repo.mesh().Geom(lev), dt, - godunov_use_ppm, godunov_use_forces_in_trans, - premac_advection_type, limiter_type, - m_allow_inflow_on_outflow); + w_mac(lev), dof_field.bcrec(), bcrec_device.data(), + repo.mesh().Geom(lev), dt, godunov_use_ppm, + godunov_use_forces_in_trans, premac_advection_type, + limiter_type, m_allow_inflow_on_outflow); } } else { amrex::Abort("Invalid godunov scheme"); @@ -360,8 +359,9 @@ struct AdvectionOp multiphase::hybrid_fluxes( repo, ICNS::ndim, iconserv, (*flux_x), (*flux_y), (*flux_z), dof_field, src_term, rho_o, u_mac, v_mac, w_mac, - dof_field.bcrec_device().data(), rho_o.bcrec_device().data(), - dt, mflux_scheme, godunov_use_forces_in_trans); + dof_field.bcrec(), dof_field.bcrec_device().data(), + rho_o.bcrec(), rho_o.bcrec_device().data(), dt, mflux_scheme, + godunov_use_forces_in_trans); } amrex::Vector> fluxes( diff --git a/amr-wind/equation_systems/vof/vof_momentum_flux.H b/amr-wind/equation_systems/vof/vof_momentum_flux.H index c66131eac2..3c66aef57a 100644 --- a/amr-wind/equation_systems/vof/vof_momentum_flux.H +++ b/amr-wind/equation_systems/vof/vof_momentum_flux.H @@ -18,8 +18,10 @@ static void hybrid_fluxes( const Field& u_mac, const Field& v_mac, const Field& w_mac, - amrex::BCRec const* velbc, - amrex::BCRec const* rhobc, + amrex::Vector const& velbc, + amrex::BCRec const* velbc_d, + amrex::Vector const& rhobc, + amrex::BCRec const* rhobc_d, const amrex::Real dt, godunov::scheme mflux_scheme, bool use_forces_in_trans) @@ -42,6 +44,13 @@ static void hybrid_fluxes( auto ftmp_z = repo.create_scratch_field(ncomp, 0, amr_wind::FieldLoc::ZFACE); + auto face_x = + repo.create_scratch_field(ncomp, 0, amr_wind::FieldLoc::XFACE); + auto face_y = + repo.create_scratch_field(ncomp, 0, amr_wind::FieldLoc::YFACE); + auto face_z = + repo.create_scratch_field(ncomp, 0, amr_wind::FieldLoc::ZFACE); + // Create iconserv = 0 array for density, avoid multiplication by face vel amrex::Gpu::DeviceVector idnsty; idnsty.resize(1, -1); @@ -59,9 +68,6 @@ static void hybrid_fluxes( fvm::Godunov::nghost_src); amrex::MultiFab::Copy( fq, src_term(lev), 0, 0, ncomp, fvm::Godunov::nghost_src); - amrex::MultiFab frho( - src_term(lev).boxArray(), src_term(lev).DistributionMap(), 1, - fvm::Godunov::nghost_src); amrex::MultiFab::Multiply( q, rho_o(lev), 0, 0, 1, fvm::Godunov::nghost_state); @@ -75,8 +81,26 @@ static void hybrid_fluxes( fq, rho_o(lev), 0, 1, 1, fvm::Godunov::nghost_src); amrex::MultiFab::Multiply( fq, rho_o(lev), 0, 2, 1, fvm::Godunov::nghost_src); + + amrex::MultiFab frho( + src_term(lev).boxArray(), src_term(lev).DistributionMap(), 1, + fvm::Godunov::nghost_src); frho.setVal(0.0); + std::string advection_type = "Godunov"; + int limiter_type = PPM::UPWIND; // default + if (mflux_scheme == godunov::scheme::MINMOD) { + limiter_type = PPM::MINMOD; + } else if (mflux_scheme != godunov::scheme::UPWIND) { + amrex::Abort("Unknown limiter type in vof_momentum_flux.H"); + } + + const bool known_edge_state = false; + const bool godunov_use_ppm = true; + const bool is_velocity = false; + const bool allow_inflow_on_outflow = true; + const bool fluxes_are_area_weighted = false; + amrex::MFItInfo mfi_info; if (amrex::Gpu::notInLaunchRegion()) { mfi_info.EnableTiling(amrex::IntVect(1024, 1024, 1024)) @@ -88,12 +112,7 @@ static void hybrid_fluxes( for (amrex::MFIter mfi(dof_field(lev), mfi_info); mfi.isValid(); ++mfi) { const auto& bx = mfi.tilebox(); - const auto& bxg1 = amrex::grow(bx, 1); - const auto& xbx = amrex::surroundingNodes(bx, 0); - const auto& ybx = amrex::surroundingNodes(bx, 1); - const auto& zbx = amrex::surroundingNodes(bx, 2); - // Memory for working arrays in godunov::compute_fluxes - amrex::FArrayBox tmpfab(amrex::grow(bx, 1), ncomp * 14); + // Arrays for easy reference (w for working array) auto F_x = flux_x(lev).array(mfi); auto F_y = flux_y(lev).array(mfi); @@ -101,17 +120,27 @@ static void hybrid_fluxes( auto Fw_x = (*ftmp_x)(lev).array(mfi); auto Fw_y = (*ftmp_y)(lev).array(mfi); auto Fw_z = (*ftmp_z)(lev).array(mfi); - // At this point in the solve, these are advected - auto ar_x = advrho_x(lev).const_array(mfi); - auto ar_y = advrho_y(lev).const_array(mfi); - auto ar_z = advrho_z(lev).const_array(mfi); + + // Temporary divu + amrex::FArrayBox divufab(amrex::grow(bx, 1), 1); + divufab.setVal(0.0); + const auto& divu = divufab.array(); + // Compute momentum flux quantities (multiplied by face velocity) - godunov::compute_fluxes( - lev, bx, ncomp, Fw_x, Fw_y, Fw_z, q.const_array(mfi), + HydroUtils::ComputeFluxesOnBoxFromState( + bx, ncomp, mfi, q.const_array(mfi), Fw_x, Fw_y, Fw_z, + (*face_x)(lev).array(mfi), (*face_y)(lev).array(mfi), + (*face_z)(lev).array(mfi), known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), - w_mac(lev).const_array(mfi), fq.const_array(mfi), velbc, - iconserv.data(), tmpfab.dataPtr(), geom, dt, mflux_scheme, - use_forces_in_trans); + w_mac(lev).const_array(mfi), divu, fq.const_array(mfi), + geom[lev], dt, velbc, velbc_d, iconserv.data(), godunov_use_ppm, + use_forces_in_trans, is_velocity, fluxes_are_area_weighted, + advection_type, limiter_type, allow_inflow_on_outflow); + + const auto& bxg1 = amrex::grow(bx, 1); + const auto& xbx = amrex::surroundingNodes(bx, 0); + const auto& ybx = amrex::surroundingNodes(bx, 1); + const auto& zbx = amrex::surroundingNodes(bx, 2); // Where interface is present, replace current momentum flux // quantities with those from mflux_scheme @@ -150,12 +179,20 @@ static void hybrid_fluxes( }); // Compute density flux quantities (not multiplied by face velocity) - godunov::compute_fluxes( - lev, bx, 1, Fw_x, Fw_y, Fw_z, rho_o(lev).const_array(mfi), + HydroUtils::ComputeFluxesOnBoxFromState( + bx, 1, mfi, rho_o(lev).const_array(mfi), Fw_x, Fw_y, Fw_z, + (*face_x)(lev).array(mfi), (*face_y)(lev).array(mfi), + (*face_z)(lev).array(mfi), known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), - w_mac(lev).const_array(mfi), frho.const_array(mfi), rhobc, - idnsty.data(), tmpfab.dataPtr(), geom, dt, mflux_scheme, - use_forces_in_trans); + w_mac(lev).const_array(mfi), divu, frho.const_array(mfi), + geom[lev], dt, rhobc, rhobc_d, idnsty.data(), godunov_use_ppm, + use_forces_in_trans, is_velocity, fluxes_are_area_weighted, + advection_type, limiter_type, allow_inflow_on_outflow); + + // At this point in the solve, these are advected + auto ar_x = advrho_x(lev).const_array(mfi); + auto ar_y = advrho_y(lev).const_array(mfi); + auto ar_z = advrho_z(lev).const_array(mfi); // When interface is present, divide by interpolated density to get // Favre-averaged flux value. Multiply all fluxes with advected diff --git a/unit_tests/multiphase/test_mflux_schemes.cpp b/unit_tests/multiphase/test_mflux_schemes.cpp index fde1e46a5c..fa2977d7b5 100644 --- a/unit_tests/multiphase/test_mflux_schemes.cpp +++ b/unit_tests/multiphase/test_mflux_schemes.cpp @@ -1,8 +1,6 @@ #include "aw_test_utils/MeshTest.H" - -#include "amr-wind/convection/incflo_godunov_minmod.H" -#include "amr-wind/convection/incflo_godunov_upwind.H" +#include "hydro_godunov_ppm.H" namespace amr_wind_tests { @@ -62,24 +60,37 @@ void init_scalar_uniform(amr_wind::Field& fld, amrex::Real cst) void get_output_upwind( amr_wind::Field& fld, + amr_wind::Field& mac_fld, + amrex::Real dt, int ii, int jj, int kk, amrex::Real& Im, amrex::Real& Ip) { + const int nlevels = fld.repo().num_active_levels(); amrex::Gpu::DeviceVector dout(2, 0.0); auto* dout_ptr = dout.data(); - const int nlevels = fld.repo().num_active_levels(); + const auto* pbc = fld.bcrec_device().data(); + for (int lev = 0; lev < nlevels; ++lev) { + const auto& dx = fld.repo().mesh().Geom(lev).CellSizeArray(); + amrex::Box const& domain = fld.repo().mesh().Geom(lev).Domain(); + const auto dlo = amrex::lbound(domain); + const auto dhi = amrex::ubound(domain); + + auto limiter = PPM::upwind(); for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) { auto bx = mfi.validbox(); const auto& farr = fld(lev).const_array(mfi); + const auto& vel_mac = mac_fld(lev).const_array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { amrex::Real im_tmp, ip_tmp; - Godunov_upwind_fpu(i, j, k, 0, im_tmp, ip_tmp, farr); + PPM::PredictStateOnXFace( + i, j, k, 0, dt, dx[0], im_tmp, ip_tmp, farr, vel_mac, + pbc[0], dlo.x, dhi.x, limiter, PPM::UPWIND); if (i == ii && j == jj && k == kk) { dout_ptr[0] = im_tmp; dout_ptr[1] = ip_tmp; @@ -107,6 +118,9 @@ void get_output_minmod( amrex::Gpu::DeviceVector dout(2, 0.0); auto* dout_ptr = dout.data(); const auto* pbc = fld.bcrec_device().data(); + + auto limiter = PPM::minmod(); + for (int lev = 0; lev < nlevels; ++lev) { const auto& dx = fld.repo().mesh().Geom(lev).CellSizeArray(); amrex::Box const& domain = fld.repo().mesh().Geom(lev).Domain(); @@ -122,17 +136,17 @@ void get_output_minmod( amrex::Real im_tmp = 0.0; amrex::Real ip_tmp = 0.0; if (dir == 0) { - Godunov_minmod_fpu_x( + PPM::PredictStateOnXFace( i, j, k, 0, dt, dx[0], im_tmp, ip_tmp, farr, vel_mac, - pbc[0], dlo.x, dhi.x); + pbc[0], dlo.x, dhi.x, limiter, PPM::MINMOD); } else if (dir == 1) { - Godunov_minmod_fpu_y( + PPM::PredictStateOnYFace( i, j, k, 0, dt, dx[1], im_tmp, ip_tmp, farr, vel_mac, - pbc[0], dlo.y, dhi.y); + pbc[0], dlo.y, dhi.y, limiter, PPM::MINMOD); } else if (dir == 2) { - Godunov_minmod_fpu_z( + PPM::PredictStateOnZFace( i, j, k, 0, dt, dx[2], im_tmp, ip_tmp, farr, vel_mac, - pbc[0], dlo.z, dhi.z); + pbc[0], dlo.z, dhi.z, limiter, PPM::MINMOD); } if (i == ii && j == jj && k == kk) { dout_ptr[0] = im_tmp; @@ -185,7 +199,20 @@ TEST_F(MFluxSchemeTest, upwind) amrex::Real Ip = 0.0; // Initialize field variable - auto& sc = repo.declare_field("scalar", 1, 1); + auto& sc = repo.declare_field("scalar", 1, 2); + + // Parameters to be referenced during test + amrex::Real adv_vel = 1.2; + amrex::Real dt = 0.25; + + // Initialize mac velocity + repo.declare_face_normal_field({"umac", "vmac", "wmac"}, 1, 1, 1); + auto& umac = repo.get_field("umac"); + auto& vmac = repo.get_field("vmac"); + auto& wmac = repo.get_field("wmac"); + umac.setVal(adv_vel); + vmac.setVal(adv_vel); + wmac.setVal(-adv_vel); /* -- Constant field portion of test -- */ // Set up field @@ -195,32 +222,38 @@ TEST_F(MFluxSchemeTest, upwind) int i = m_nx / 2; int j = m_nx / 2; int k = m_nx / 2; - get_output_upwind(sc, i, j, k, Im, Ip); + get_output_upwind(sc, umac, dt, i, j, k, Im, Ip); // Check values EXPECT_NEAR(sc_cst, Im, m_tol); EXPECT_NEAR(sc_cst, Ip, m_tol); - /* -- Increasing in slope: x, y, z -- */ - for (int n = 0; n < 3; ++n) { - // Set up field - init_scalar_increasing(sc, n); - // Compute interpolated quantities at each face - get_output_upwind(sc, i, j, k, Im, Ip); - // Check values - EXPECT_NEAR(i * i, Im, m_tol); - EXPECT_NEAR(i * i, Ip, m_tol); - } - /* -- Change in sign of slope: x, y, z -- */ - for (int n = 0; n < 3; ++n) { - // Set up field - init_scalar_slopechange(sc, n, i); - // Compute interpolated quantities at each face - get_output_upwind(sc, i, j, k, Im, Ip); - // Check values - EXPECT_NEAR(0, Im, m_tol); - EXPECT_NEAR(0, Ip, m_tol); - } + // Set up field + int n = 0; + init_scalar_slopechange(sc, n, i); + // Compute interpolated quantities at each face + get_output_upwind(sc, umac, dt, i, j, k, Im, Ip); + // Check values + EXPECT_NEAR(0, Im, m_tol); + EXPECT_NEAR(0, Ip, m_tol); + + // Set up field + n = 1; + init_scalar_slopechange(sc, n, i); + // Compute interpolated quantities at each face + get_output_upwind(sc, vmac, dt, i, j, k, Im, Ip); + // Check values + EXPECT_NEAR(0, Im, m_tol); + EXPECT_NEAR(0, Ip, m_tol); + + // Set up field + n = 2; + init_scalar_slopechange(sc, n, i); + // Compute interpolated quantities at each face + get_output_upwind(sc, wmac, dt, i, j, k, Im, Ip); + // Check values + EXPECT_NEAR(0, Im, m_tol); + EXPECT_NEAR(0, Ip, m_tol); } TEST_F(MFluxSchemeTest, minmod) @@ -237,7 +270,7 @@ TEST_F(MFluxSchemeTest, minmod) amrex::Real Ip = 0.0; // Initialize field variable - auto& sc = repo.declare_field("scalar", 1, 1); + auto& sc = repo.declare_field("scalar", 1, 2); // Initialize mac velocity repo.declare_face_normal_field({"umac", "vmac", "wmac"}, 1, 1, 1); @@ -311,7 +344,7 @@ TEST_F(MFluxSchemeTest, minmodbdy) amrex::Real Ip = 0.0; // Initialize field variable - auto& sc = repo.declare_field("scalar", 1, 1); + auto& sc = repo.declare_field("scalar", 1, 2); // Initialize mac velocity repo.declare_face_normal_field({"umac", "vmac", "wmac"}, 1, 1, 1);