From 543feb3877fc08af6bed23534f37a74c198b67d8 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 10 Mar 2024 18:49:47 -0700 Subject: [PATCH 1/9] Add option to interpolate data on faces linearly in the tangential direction -- previously we only did piecewise constant within the face --- Src/AmrCore/AMReX_InterpFaceReg_2D_C.H | 6 +- Src/AmrCore/AMReX_InterpFaceReg_3D_C.H | 21 +- Src/AmrCore/AMReX_InterpFaceRegister.cpp | 12 +- Src/AmrCore/AMReX_Interp_C.H | 152 +++++++++- Src/AmrCore/AMReX_Interpolater.H | 126 ++++++++- Src/AmrCore/AMReX_Interpolater.cpp | 341 ++++++++++++++++++++++- 6 files changed, 639 insertions(+), 19 deletions(-) diff --git a/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H b/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H index 5b80958271d..9af617d4216 100644 --- a/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H +++ b/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H @@ -6,12 +6,12 @@ namespace amrex { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interp_face_reg (int i, int j, IntVect const& rr, Array4 const& fine, int scomp, Array4 const& crse, Array4 const& slope, int ncomp, - Box const& domface, int idim) + Box const& per_grown_domain, int idim) { int ic = amrex::coarsen(i,rr[0]); int jc = amrex::coarsen(j,rr[1]); if (idim == 0) { - if (jc == domface.smallEnd(1) || jc == domface.bigEnd(1)) { + if (jc == per_grown_domain.smallEnd(1) || jc == per_grown_domain.bigEnd(1)) { for (int n = 0; n < ncomp; ++n) { fine(i,j,0,n+scomp) = crse(ic,jc,0,n); } @@ -35,7 +35,7 @@ void interp_face_reg (int i, int j, IntVect const& rr, Array4 const& fine, } } } else { - if (ic == domface.smallEnd(0) || ic == domface.bigEnd(0)) { + if (ic == per_grown_domain.smallEnd(0) || ic == per_grown_domain.bigEnd(0)) { for (int n = 0; n < ncomp; ++n) { fine(i,j,0,n+scomp) = crse(ic,jc,0,n); } diff --git a/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H b/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H index 2df7fef055a..2d1728cb929 100644 --- a/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H +++ b/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H @@ -6,7 +6,7 @@ namespace amrex { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const& fine, int scomp, Array4 const& crse, Array4 const& slope, int ncomp, - Box const& domface, int idim) + Box const& per_grown_domain, int idim) { int ic = amrex::coarsen(i,rr[0]); int jc = amrex::coarsen(j,rr[1]); @@ -15,7 +15,8 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } - if (jc != domface.smallEnd(1) && jc != domface.bigEnd(1) && rr[1] > 1) { +#if 0 + if (jc != per_grown_domain.smallEnd(1) && jc != per_grown_domain.bigEnd(1) && rr[1] > 1) { Real sfy = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)); @@ -34,8 +35,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += yoff * slope(i,j,k,n) * sfy; } } - - if (kc != domface.smallEnd(2) && kc != domface.bigEnd(2) && rr[2] > 1) { + if (kc != per_grown_domain.smallEnd(2) && kc != per_grown_domain.bigEnd(2) && rr[2] > 1) { Real sfz = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)); @@ -54,11 +54,13 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += zoff * slope(i,j,k,n) * sfz; } } +#endif } else if (idim == 1) { for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } - if (ic != domface.smallEnd(0) && ic != domface.bigEnd(0) && rr[0] > 1) { +#if 0 + if (ic != per_grown_domain.smallEnd(0) && ic != per_grown_domain.bigEnd(0) && rr[0] > 1) { Real sfx = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)); @@ -78,7 +80,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const } } - if (kc != domface.smallEnd(2) && kc != domface.bigEnd(2) && rr[2] > 1) { + if (kc != per_grown_domain.smallEnd(2) && kc != per_grown_domain.bigEnd(2) && rr[2] > 1) { Real sfz = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)); @@ -97,11 +99,13 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += zoff * slope(i,j,k,n) * sfz; } } +#endif } else { for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } - if (ic != domface.smallEnd(0) && ic != domface.bigEnd(0) && rr[0] > 1) { +#if 0 + if (ic != per_grown_domain.smallEnd(0) && ic != per_grown_domain.bigEnd(0) && rr[0] > 1) { Real sfx = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)); @@ -121,7 +125,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const } } - if (jc != domface.smallEnd(1) && jc != domface.bigEnd(1) && rr[1] > 1) { + if (jc != per_grown_domain.smallEnd(1) && jc != per_grown_domain.bigEnd(1) && rr[1] > 1) { Real sfy = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)); @@ -140,6 +144,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += yoff * slope(i,j,k,n) * sfy; } } +#endif } } diff --git a/Src/AmrCore/AMReX_InterpFaceRegister.cpp b/Src/AmrCore/AMReX_InterpFaceRegister.cpp index ae5beb1985c..33de90c15b4 100644 --- a/Src/AmrCore/AMReX_InterpFaceRegister.cpp +++ b/Src/AmrCore/AMReX_InterpFaceRegister.cpp @@ -20,6 +20,14 @@ void InterpFaceRegister::define (BoxArray const& fba, DistributionMapping const& m_crse_geom = amrex::coarsen(m_fine_geom, m_ref_ratio); + // We don't need to worry about face-based domain because this is only used in the tangential interpolation + Box per_grown_domain = m_crse_geom.Domain(); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (m_crse_geom.isPeriodic(dim)) { + per_grown_domain.grow(dim,1); + } + } + constexpr int crse_fine_face = 1; constexpr int fine_fine_face = 0; // First, set the face mask to 1 (i.e., coarse/fine boundary), @@ -91,7 +99,7 @@ namespace { Array4 slope; Array4 crse; Array4 mask; - Box domface; + Box per_grown_domain; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box box() const noexcept { return Box(mask); } }; @@ -151,7 +159,7 @@ InterpFaceRegister::interp (Array const& fine, // NOL { if (tag.mask(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, tag.fine, scomp, tag.crse, - tag.slope, ncomp, tag.domface, idim); + tag.slope, ncomp, tag.per_grown_domain, idim); } }); } else diff --git a/Src/AmrCore/AMReX_Interp_C.H b/Src/AmrCore/AMReX_Interp_C.H index 86e7935a436..e209be2581d 100644 --- a/Src/AmrCore/AMReX_Interp_C.H +++ b/Src/AmrCore/AMReX_Interp_C.H @@ -15,8 +15,8 @@ namespace amrex { // // Fill fine values with piecewise-constant interpolation of coarse data. -// Operate only on faces that overlap--ie, only fill the fine faces that make up each -// coarse face, leave the in-between faces alone. +// Operate only on faces that overlap -- i.e., only fill the fine faces that +// make up each coarse face, leave the in-between faces alone. // template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void @@ -84,6 +84,154 @@ face_linear_face_interp_z (int fi, int fj, int fk, int n, Array4 const& fine, } } +// +// Fill fine values with tangential interpolation of coarse data. +// Operate only on faces that overlap -- i.e., only fill the fine faces that +// make up each coarse face, leave the in-between faces alone. +// +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, + Array4 const& crse, Array4 const& mask, + IntVect const& ratio, Box const& per_grown_domain, int dim) noexcept +{ + int ci = amrex::coarsen(i, ratio[0]); +#if (AMREX_SPACEDIM >= 2) + int cj = amrex::coarsen(j, ratio[1]); +#else + int cj = 0; +#endif +#if (AMREX_SPACEDIM == 3) + int ck = amrex::coarsen(k, ratio[2]); +#else + int ck = 0; +#endif + if (dim == 0 && ci*ratio[0] == i) { + // Check solve mask to ensure we don't overwrite valid fine data. + if (!mask || mask(ci, cj, ck, n)) { + fine(i, j, k, n) = crse(ci, cj, ck, n); +#if 0 +#if (AMREX_SPACEDIM >= 2) + if (cj != per_grown_domain.smallEnd(1) && cj != per_grown_domain.bigEnd(1) && ratio[1] > 1) { + Real sfy = Real(1.0); + Real dc = Real(0.5) * (crse(ci,cj+1,ck,n) - crse(ci,cj-1,ck,n)); + Real df = Real(2.0) * (crse(ci,cj+1,ck,n) - crse(ci,cj ,ck,n)); + Real db = Real(2.0) * (crse(ci,cj ,ck,n) - crse(ci,cj-1,ck,n)); + Real sy = (df*db >= Real(0.0)) ? + amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + if (dc != Real(0.0)) { + sfy = amrex::min(sfy, sy / dc); + } + Real slope = dc; + Real yoff = (static_cast(j - cj*ratio[1]) + Real(0.5)) / Real(ratio[1]) - Real(0.5); + fine(i,j,k,n) += yoff * slope * sfy; + } // jc +#endif +#if (AMREX_SPACEDIM == 3) + if (ck != per_grown_domain.smallEnd(2) && ck != per_grown_domain.bigEnd(2) && ratio[2] > 1) { + Real sfz = Real(1.0); + Real dc = Real(0.5) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck-1,n)); + Real df = Real(2.0) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck ,n)); + Real db = Real(2.0) * (crse(ci,cj,ck ,n) - crse(ci,cj,ck-1,n)); + Real sz = (df*db >= Real(0.0)) ? + amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc)); + if (dc != Real(0.0)) { + sfz = amrex::min(sfz, sz / dc); + } + Real slope = dc; + Real zoff = (static_cast(k - ck*ratio[2]) + Real(0.5)) / Real(ratio[2]) - Real(0.5); + fine(i,j,k,n) += zoff * slope * sfz; + } // ck +#endif +#endif + } // mask + } // dim + +#if (AMREX_SPACEDIM >= 2) + if (dim == 1 && cj*ratio[1] == j) { + // Check solve mask to ensure we don't overwrite valid fine data. + if (!mask || mask(ci, cj, ck, n)) { + fine(i, j, k, n) = crse(ci, cj, ck, n); + if (ci != per_grown_domain.smallEnd(0) && ci != per_grown_domain.bigEnd(0) && ratio[0] > 1) { + Real sfx = Real(1.0); + Real dc = Real(0.5) * (crse(ci+1,cj,ck,n) - crse(ci-1,cj,ck,n)); + Real df = Real(2.0) * (crse(ci+1,cj,ck,n) - crse(ci ,cj,ck,n)); + Real db = Real(2.0) * (crse(ci ,cj,ck,n) - crse(ci-1,cj,ck,n)); + Real sx = (df*db >= Real(0.0)) ? + amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + if (dc != Real(0.0)) { + sfx = amrex::min(sfx, sx / dc); + } + Real slope = dc; + Real xoff = (static_cast(i - ci*ratio[0]) + Real(0.5)) / Real(ratio[0]) - Real(0.5); + fine(i,j,k,n) += xoff * slope * sfx; + } // ci +#if (AMREX_SPACEDIM == 3) + if (ck != per_grown_domain.smallEnd(2) && ck != per_grown_domain.bigEnd(2) && ratio[2] > 1) { + Real sfz = Real(1.0); + Real dc = Real(0.5) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck-1,n)); + Real df = Real(2.0) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck ,n)); + Real db = Real(2.0) * (crse(ci,cj,ck ,n) - crse(ci,cj,ck-1,n)); + Real sz = (df*db >= Real(0.0)) ? + amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc)); + if (dc != Real(0.0)) { + sfz = amrex::min(sfz, sz / dc); + } + Real slope = dc; + Real zoff = (static_cast(k - ck*ratio[2]) + Real(0.5)) / Real(ratio[2]) - Real(0.5); + fine(i,j,k,n) += zoff * slope * sfz; + } // ck + } // mask +#endif // SPACEDIM >= 3 + } // dim == 1 +#endif // SPACEDIM >= 2 + +#if (AMREX_SPACEDIM == 3) + if (dim == 2 && ck*ratio[2] == k) { + // Check solve mask to ensure we don't overwrite valid fine data. + if (!mask || mask(ci, cj, ck, n)) { + fine(i, j, k, n) = crse(ci, cj, ck, n); +#if 0 + if (ci != per_grown_domain.smallEnd(0) && ci != per_grown_domain.bigEnd(0) && ratio[0] > 1) { + Real sfx = Real(1.0); + Real dc = Real(0.5) * (crse(ci+1,cj,ck,n) - crse(ci-1,cj,ck,n)); + Real df = Real(2.0) * (crse(ci+1,cj,ck,n) - crse(ci ,cj,ck,n)); + Real db = Real(2.0) * (crse(ci ,cj,ck,n) - crse(ci-1,cj,ck,n)); + Real sx = (df*db >= Real(0.0)) ? + amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + if (dc != Real(0.0)) { + sfx = amrex::min(sfx, sx / dc); + } + Real slope = dc; + Real xoff = (static_cast(i - ci*ratio[0]) + Real(0.5)) / Real(ratio[0]) - Real(0.5); + fine(i,j,k,n) += xoff * slope * sfx; + } // ci + if (cj != per_grown_domain.smallEnd(1) && cj != per_grown_domain.bigEnd(1) && ratio[1] > 1) { + Real sfy = Real(1.0); + Real dc = Real(0.5) * (crse(ci,cj+1,ck,n) - crse(ci,cj-1,ck,n)); + Real df = Real(2.0) * (crse(ci,cj+1,ck,n) - crse(ci,cj ,ck,n)); + Real db = Real(2.0) * (crse(ci,cj ,ck,n) - crse(ci,cj-1,ck,n)); + Real sy = (df*db >= Real(0.0)) ? + amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + if (dc != Real(0.0)) { + sfy = amrex::min(sfy, sy / dc); + } + Real slope = dc; + Real yoff = (static_cast(j - cj*ratio[1]) + Real(0.5)) / Real(ratio[1]) - Real(0.5); + fine(i,j,k,n) += yoff * slope * sfy; + } // cj +#endif + } // mask + } // dim == 2 +#endif +} + // // Do linear in dir, pc transverse to dir, leave alone the fine values // lining up with coarse edges--assume these have been set to hold the diff --git a/Src/AmrCore/AMReX_Interpolater.H b/Src/AmrCore/AMReX_Interpolater.H index e1210a83329..c13fb283a35 100644 --- a/Src/AmrCore/AMReX_Interpolater.H +++ b/Src/AmrCore/AMReX_Interpolater.H @@ -654,9 +654,9 @@ public: /** -* \brief Bilinear interpolation on face data. +* \brief Piecewise constant tangential interpolation / linear normal interpolation of face data. * -* Bilinear interpolation on data. +* Piecewise constant tangential interpolation / linear normal interpolation of face data. */ class FaceLinear : @@ -772,6 +772,127 @@ public: RunOn runon) override; +}; + +/** +* \brief Bilinear tangential interpolation / linear normal interpolation of face data. +* +* Bilinear tangential interpolation / linear normal interpolation of face data. +*/ +class FaceConservativeLinear + : + public Interpolater +{ +public: + /** + * \brief Returns coarsened box given fine box and refinement ratio. + * + * \param fine + * \param ratio + */ + Box CoarseBox (const Box& fine, int ratio) override; + + /** + * \brief Returns coarsened box given fine box and refinement ratio. + * + * \param fine + * \param ratio + */ + Box CoarseBox (const Box& fine, const IntVect& ratio) override; + + /** + * \brief Coarse to fine interpolation in space. + * + * \param crse + * \param crse_comp + * \param fine + * \param fine_comp + * \param ncomp + * \param fine_region + * \param ratio + * \param crse_geom + * \param fine_geom + * \param bcr + * \param actual_comp + * \param actual_state + */ + void interp (const FArrayBox& crse, + int crse_comp, + FArrayBox& fine, + int fine_comp, + int ncomp, + const Box& fine_region, + const IntVect& ratio, + const Geometry& crse_geom, + const Geometry& fine_geom, + Vector const& bcr, + int actual_comp, + int actual_state, + RunOn runon) override; + + /** + * \brief Coarse to fine interpolation in space for face-based data. + * + * \param crse + * \param crse_comp + * \param fine + * \param fine_comp + * \param ncomp + * \param fine_region + * \param ratio + * \param solve_mask + * \param crse_geom + * \param fine_geom + * \param bcr + * \param bccomp + * \param runon + */ + void interp_face (const FArrayBox& crse, + int crse_comp, + FArrayBox& fine, + int fine_comp, + int ncomp, + const Box& fine_region, + const IntVect& ratio, + const IArrayBox& solve_mask, + const Geometry& crse_geom, + const Geometry& fine_geom, + Vector const & bcr, + int bccomp, + RunOn runon) override; + + /** + * \brief Coarse to fine interpolation in space. + * + * \param crse + * \param crse_comp + * \param fine + * \param fine_comp + * \param ncomp + * \param fine_region + * \param ratio + * \param crse_geom + * \param fine_geom + * \param bcr + * \param actual_comp + * \param actual_state + */ + void interp_arr (Array const& crse, + int crse_comp, + Array const& fine, + int fine_comp, + int ncomp, + const Box& fine_region, + const IntVect& ratio, + Array const& solve_mask, + const Geometry& /*crse_geom*/, + const Geometry& /*fine_geom*/, + Vector > const& /*bcr*/, + int /*actual_comp*/, + int /*actual_state*/, + RunOn runon) override; + + }; /** @@ -836,6 +957,7 @@ extern AMREX_EXPORT PCInterp pc_interp; extern AMREX_EXPORT NodeBilinear node_bilinear_interp; extern AMREX_EXPORT FaceDivFree face_divfree_interp; extern AMREX_EXPORT FaceLinear face_linear_interp; +extern AMREX_EXPORT FaceConservativeLinear face_cons_linear_interp; extern AMREX_EXPORT CellConservativeLinear lincc_interp; extern AMREX_EXPORT CellConservativeLinear cell_cons_interp; extern AMREX_EXPORT CellBilinear cell_bilinear_interp; diff --git a/Src/AmrCore/AMReX_Interpolater.cpp b/Src/AmrCore/AMReX_Interpolater.cpp index 601b8b4b861..339a1f9dcd0 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -11,7 +11,7 @@ namespace amrex { /* - * PCInterp, NodeBilinear, FaceLinear, CellConservativeLinear, and + * PCInterp, NodeBilinear, FaceLinear, CellConservativeLinear and * CellBilinear are supported for all dimensions on cpu and gpu. * * CellConservativeProtected only works in 2D and 3D on cpu and gpu @@ -23,6 +23,8 @@ namespace amrex { * * CellConservativeQuartic only works with ref ratio of 2 on cpu and gpu. * + * FaceConservativeLinear works in 2D and 3D on cpu and gpu. + * * FaceDivFree works in 2D and 3D on cpu and gpu. * The algorithm is restricted to ref ratio of 2. */ @@ -33,6 +35,7 @@ namespace amrex { PCInterp pc_interp; NodeBilinear node_bilinear_interp; FaceLinear face_linear_interp; +FaceConservativeLinear face_cons_linear_interp; FaceDivFree face_divfree_interp; CellConservativeLinear lincc_interp; CellConservativeLinear cell_cons_interp(false); @@ -142,7 +145,14 @@ FaceLinear::interp (const FArrayBox& crse, RunOn runon) { // - // This version is called from InterpFromCoarseLevel + // This version is called from FillPatchInterp which is called by + // InterpFromCoarseLevel in AMReX_FillPatchUtil_I.H + // + // It assumes no existing fine values that need to be preserved (unlike interp_face below) + // + // Inside each call to face_linear_interp_* (in AMRex_Interp_*D_C.H), we do: + // * on fine faces which overlie crse faces, the fine value is set to the crse value (piecewise constant) + // * on fine faces which are between two crse faces, the fine value is set to the average of the crse values (linear) // BL_PROFILE("FaceLinear::interp()"); @@ -193,6 +203,18 @@ FaceLinear::interp_face (const FArrayBox& crse, const int /*bccomp*/, RunOn runon) { + // + // This version is called from InterpFace which is called from the version FillPatchTwoLevels_doit + // that takes a single MF (in AMReX_FillPatchUtil_I.H) + // + // It assumes there are existing fine values which we want to preserve (unlike interp above) + // + // We do the interpolation in two steps: + // 1) face_linear_face_interp_*: on fine faces which overlie crse faces, the fine value is set to the crse value (piecewise constant) ONLY IF + // there is not already fine data there + // 2) face_linear_interp_*: on fine faces which are between two crse faces, the fine value is set to the average of the values + // on the faces overlying -- this uses only the results of step 1, it does not take the crse values + // BL_PROFILE("FaceLinear::interp_face()"); AMREX_ASSERT(AMREX_D_TERM(fine_region.type(0),+fine_region.type(1),+fine_region.type(2)) == 1); @@ -283,6 +305,17 @@ void FaceLinear::interp_arr (Array const& crse, const int /*actual_state*/, const RunOn runon) { + // + // This version is called from FillPatchTwoLevels_doit (that takes an Array of MF*) in AMReX_FillPatchUtil_I.H + // + // It assumes there are existing fine values which we want to preserve (like face_interp, unlike interp above) + // + // We do the interpolation in two steps: + // 1) face_linear_face_interp_*: on fine faces which overlie crse faces, the fine value is set to the crse value (piecewise constant) ONLY IF + // there is not already fine data there + // 2) face_linear_interp_*: on fine faces which are between two crse faces, the fine value is set to the average of the values + // on the faces overlying -- this uses only the results of step 1, it does not take the crse values + // BL_PROFILE("FaceLinear::interp_arr()"); Array types; @@ -377,6 +410,310 @@ void FaceLinear::interp_arr (Array const& crse, }); } +#if (AMREX_SPACEDIM > 1) +Box +FaceConservativeLinear::CoarseBox (const Box& fine, int ratio) +{ + return CoarseBox(fine, IntVect(ratio)); +} + +Box +FaceConservativeLinear::CoarseBox (const Box& fine, const IntVect& ratio) +{ + IntVect ng(1); + for (int i = 0; i < AMREX_SPACEDIM; i++) { + if (fine.type(i) == IndexType::NODE or ratio[i] == 1) { + ng[i] = 0; + } + } + Box b = amrex::coarsen(fine,ratio); b.grow(ng); + + for (int i = 0; i < AMREX_SPACEDIM; i++) { + if (b.type(i) == IndexType::NODE) { + ng[i] = 0; + if (b.type(i) == IndexType::NODE && b.length(i) < 2) { + // Don't want degenerate boxes in nodal direction. + b.growHi(i,1); + } + } + } + return b; +} + +void +FaceConservativeLinear::interp (const FArrayBox& crse, + int crse_comp, + FArrayBox& fine, + int fine_comp, + int ncomp, + const Box& fine_region, + const IntVect& ratio, + const Geometry& crse_geom, + const Geometry& fine_geom, + Vector const& bcr, + int actual_comp, + int actual_state, + RunOn runon) +{ + // + // This version is called from FillPatchInterp which is called by + // InterpFromCoarseLevel in AMReX_FillPatchUtil_I.H + // + // It assumes no existing fine values that need to be preserved thus does not send a mask to interp_face + // + BL_PROFILE("FaceConservativeLinear::interp()"); + + AMREX_ASSERT(AMREX_D_TERM(fine_region.type(0),+fine_region.type(1),+fine_region.type(2)) == 1); + + // We intentionally do not allocate the mask so that all faces are filled from coarse values + IArrayBox dummy_mask; + int bccomp = 0; // This is also a dummy -- it's not used + interp_face(crse,crse_comp,fine,fine_comp,ncomp,fine_region,ratio,dummy_mask, + crse_geom,fine_geom,bcr,bccomp,runon); +} + +void +FaceConservativeLinear::interp_face (const FArrayBox& crse, + const int crse_comp, + FArrayBox& fine, + const int fine_comp, + const int ncomp, + const Box& fine_region, + const IntVect& ratio, + const IArrayBox& solve_mask, + const Geometry& crse_geom, + const Geometry& /*fine_geom */, + Vector const& /*bcr*/, + const int /*bccomp*/, + RunOn runon) +{ + // + // This version is called from InterpFace which is called from the version FillPatchTwoLevels_doit + // that takes a single MF (in AMReX_FillPatchUtil_I.H) + // + // It assumes there are existing fine values which we want to preserve (unlike interp above) + // + // We do the interpolation in two steps: + // 1) face_cons_linear_face_interp: on fine faces which overlie crse faces, slopes are computed (linear in 2d, bilinear in 3d) + // and the fine value is over-written ONLY IF there is not already fine data there (assuming the mask is used) + // 2) face_linear_interp_*: on fine faces which are between two crse faces, the fine value is set to the average of the values + // on the faces overlying -- this uses only the results of step 1 + // NOTE: we use the same routines as used by FaceLinear since this interpolation is only in the normal direction + // + BL_PROFILE("FaceConservativeLinear::interp_face()"); + + AMREX_ASSERT(AMREX_D_TERM(fine_region.type(0),+fine_region.type(1),+fine_region.type(2)) == 1); + Array4 const& fine_arr = fine.array(fine_comp); + Array4 const& crse_arr = crse.const_array(crse_comp); + Array4 mask_arr; + if (solve_mask.isAllocated()) { + mask_arr = solve_mask.const_array(); + } + + // We don't need to worry about face-based domain because this is only used in the tangential interpolation + Box per_grown_domain = crse_geom.Domain(); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (crse_geom.isPeriodic(dim)) { + per_grown_domain.grow(dim,1); + } + } + + // + // Fill fine ghost faces with interpolation of coarse data that is conservative linear + // in the tangential direction. + // Operate only on faces that overlap--ie, only fill the fine faces that make up each + // coarse face, leave the in-between faces alone. + // The mask ensures we do not overwrite valid fine cells. + // + if (fine_region.type(0) == IndexType::NODE) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_cons_linear_face_interp(i,j,k,n,fine_arr,crse_arr,mask_arr,ratio,per_grown_domain,0); + }); + } +#if (AMREX_SPACEDIM >= 2) + else if (fine_region.type(1) == IndexType::NODE) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_cons_linear_face_interp(i,j,k,n,fine_arr,crse_arr,mask_arr,ratio,per_grown_domain,1); + }); + } +#if (AMREX_SPACEDIM == 3) + else + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_cons_linear_face_interp(i,j,k,n,fine_arr,crse_arr,mask_arr,ratio,per_grown_domain,2); + }); + } +#endif +#endif + + // + // Interpolate unfilled grow cells using best data from + // surrounding faces of valid region, and pc-interpd data + // on fine faces overlaying coarse edges. + // + if (fine_region.type(0) == IndexType::NODE) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_linear_interp_x(i,j,k,n,fine_arr,ratio); + }); + } +#if (AMREX_SPACEDIM >= 2) + else if (fine_region.type(1) == IndexType::NODE) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_linear_interp_y(i,j,k,n,fine_arr,ratio); + }); + } +#if (AMREX_SPACEDIM == 3) + else + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_linear_interp_z(i,j,k,n,fine_arr,ratio); + }); + } +#endif +#endif +} + +void FaceConservativeLinear::interp_arr (Array const& crse, + const int crse_comp, + Array const& fine, + const int fine_comp, + const int ncomp, + const Box& fine_region, + const IntVect& ratio, + Array const& solve_mask, + const Geometry& crse_geom, + const Geometry& /*fine_geom*/, + Vector > const& /*bcr*/, + const int /*actual_comp*/, + const int /*actual_state*/, + const RunOn runon) +{ + // + // This version is called from FillPatchTwoLevels_doit (that takes an Array of MF*) in AMReX_FillPatchUtil_I.H + // + // It assumes there are existing fine values which we want to preserve (like face_interp, unlike interp above) + // + // We do the interpolation in two steps: + // 1) face_cons_linear_face_interp_*: on fine faces which overlie crse faces, we compute tangential slopes + // to compute the fine values (linear in 2d, bilinear in 3d) ONLY IF there is not already fine data there + // 2) face_cons_linear_interp_*: on fine faces which are between two crse faces, the fine value is set to the average of the values + // on the faces overlying -- this uses only the results of step 1, it does not take the crse values + // NOTE: here we use the same routines as used by FaceLinear since this interpolation is only in the normal direction + // + BL_PROFILE("FaceConservativeLinear::interp_arr()"); + + Array types; + for (int d=0; d, AMREX_SPACEDIM> crse_arr; + GpuArray, AMREX_SPACEDIM> fine_arr; + GpuArray, AMREX_SPACEDIM> mask_arr; + for (int d=0; dconst_array(crse_comp); + fine_arr[d] = fine[d]->array(fine_comp); + if (solve_mask[d] != nullptr) + { mask_arr[d] = solve_mask[d]->const_array(0); } + } + + // We don't need to worry about face-based domain because this is only used in the tangential interpolation + Box per_grown_domain = crse_geom.Domain(); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (crse_geom.isPeriodic(dim)) { + per_grown_domain.grow(dim,1); + } + } + + // + // Fill fine ghost faces with interpolation of coarse data that is conservative linear + // in the tangential direction. + // Operate only on faces that overlap--ie, only fill the fine faces that make up each + // coarse face, leave the in-between faces alone. + // The mask ensures we do not overwrite valid fine cells. + // + // Fuse the launches, 1 for each dimension, into a single launch. + AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM_FLAG(runon, + amrex::convert(fine_region,types[0]), bx0, + { + AMREX_LOOP_3D(bx0, i, j, k, + { + for (int n=0; n 1 + Box CellBilinear::CoarseBox (const Box& fine, int ratio) { From b1e5990cd33bc6553029ad599eab8161c2549bc5 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 11 Mar 2024 11:42:02 -0700 Subject: [PATCH 2/9] fix for compilation --- Src/AmrCore/AMReX_InterpFaceReg_3D_C.H | 6 ------ Src/AmrCore/AMReX_Interpolater.cpp | 6 +++--- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H b/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H index 2d1728cb929..921395be47f 100644 --- a/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H +++ b/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H @@ -15,7 +15,6 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } -#if 0 if (jc != per_grown_domain.smallEnd(1) && jc != per_grown_domain.bigEnd(1) && rr[1] > 1) { Real sfy = Real(1.0); for (int n = 0; n < ncomp; ++n) { @@ -54,12 +53,10 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += zoff * slope(i,j,k,n) * sfz; } } -#endif } else if (idim == 1) { for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } -#if 0 if (ic != per_grown_domain.smallEnd(0) && ic != per_grown_domain.bigEnd(0) && rr[0] > 1) { Real sfx = Real(1.0); for (int n = 0; n < ncomp; ++n) { @@ -99,12 +96,10 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += zoff * slope(i,j,k,n) * sfz; } } -#endif } else { for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } -#if 0 if (ic != per_grown_domain.smallEnd(0) && ic != per_grown_domain.bigEnd(0) && rr[0] > 1) { Real sfx = Real(1.0); for (int n = 0; n < ncomp; ++n) { @@ -144,7 +139,6 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += yoff * slope(i,j,k,n) * sfy; } } -#endif } } diff --git a/Src/AmrCore/AMReX_Interpolater.cpp b/Src/AmrCore/AMReX_Interpolater.cpp index 339a1f9dcd0..ec036b05399 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -422,7 +422,7 @@ FaceConservativeLinear::CoarseBox (const Box& fine, const IntVect& ratio) { IntVect ng(1); for (int i = 0; i < AMREX_SPACEDIM; i++) { - if (fine.type(i) == IndexType::NODE or ratio[i] == 1) { + if ( (fine.type(i) == IndexType::NODE) || (ratio[i] == 1) ) { ng[i] = 0; } } @@ -451,8 +451,8 @@ FaceConservativeLinear::interp (const FArrayBox& crse, const Geometry& crse_geom, const Geometry& fine_geom, Vector const& bcr, - int actual_comp, - int actual_state, + int /*actual_comp*/, + int /*actual_state*/, RunOn runon) { // From 8d447be956e3e43a579539a945a8b7a5d70efe0a Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 11 Mar 2024 11:57:56 -0700 Subject: [PATCH 3/9] remove if 0 accidentally left in --- Src/AmrCore/AMReX_Interp_C.H | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Src/AmrCore/AMReX_Interp_C.H b/Src/AmrCore/AMReX_Interp_C.H index e209be2581d..be03ce04c24 100644 --- a/Src/AmrCore/AMReX_Interp_C.H +++ b/Src/AmrCore/AMReX_Interp_C.H @@ -110,7 +110,6 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, // Check solve mask to ensure we don't overwrite valid fine data. if (!mask || mask(ci, cj, ck, n)) { fine(i, j, k, n) = crse(ci, cj, ck, n); -#if 0 #if (AMREX_SPACEDIM >= 2) if (cj != per_grown_domain.smallEnd(1) && cj != per_grown_domain.bigEnd(1) && ratio[1] > 1) { Real sfy = Real(1.0); @@ -127,7 +126,6 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, Real yoff = (static_cast(j - cj*ratio[1]) + Real(0.5)) / Real(ratio[1]) - Real(0.5); fine(i,j,k,n) += yoff * slope * sfy; } // jc -#endif #if (AMREX_SPACEDIM == 3) if (ck != per_grown_domain.smallEnd(2) && ck != per_grown_domain.bigEnd(2) && ratio[2] > 1) { Real sfz = Real(1.0); @@ -195,7 +193,6 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, // Check solve mask to ensure we don't overwrite valid fine data. if (!mask || mask(ci, cj, ck, n)) { fine(i, j, k, n) = crse(ci, cj, ck, n); -#if 0 if (ci != per_grown_domain.smallEnd(0) && ci != per_grown_domain.bigEnd(0) && ratio[0] > 1) { Real sfx = Real(1.0); Real dc = Real(0.5) * (crse(ci+1,cj,ck,n) - crse(ci-1,cj,ck,n)); @@ -226,7 +223,6 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, Real yoff = (static_cast(j - cj*ratio[1]) + Real(0.5)) / Real(ratio[1]) - Real(0.5); fine(i,j,k,n) += yoff * slope * sfy; } // cj -#endif } // mask } // dim == 2 #endif From a6466d287f65435532dcb5b8fedf92e129e46a24 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 11 Mar 2024 12:30:05 -0700 Subject: [PATCH 4/9] fix oops --- Src/AmrCore/AMReX_Interp_C.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Src/AmrCore/AMReX_Interp_C.H b/Src/AmrCore/AMReX_Interp_C.H index be03ce04c24..327f946219f 100644 --- a/Src/AmrCore/AMReX_Interp_C.H +++ b/Src/AmrCore/AMReX_Interp_C.H @@ -183,8 +183,8 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, Real zoff = (static_cast(k - ck*ratio[2]) + Real(0.5)) / Real(ratio[2]) - Real(0.5); fine(i,j,k,n) += zoff * slope * sfz; } // ck - } // mask #endif // SPACEDIM >= 3 + } // mask } // dim == 1 #endif // SPACEDIM >= 2 From 6d7b1b824072da275afcb1f1d4e7e25b418c1a70 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 11 Mar 2024 12:49:33 -0700 Subject: [PATCH 5/9] add ignore_unused --- Src/AmrCore/AMReX_Interp_C.H | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/Src/AmrCore/AMReX_Interp_C.H b/Src/AmrCore/AMReX_Interp_C.H index 327f946219f..0b2a2a74bb0 100644 --- a/Src/AmrCore/AMReX_Interp_C.H +++ b/Src/AmrCore/AMReX_Interp_C.H @@ -96,16 +96,20 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, IntVect const& ratio, Box const& per_grown_domain, int dim) noexcept { int ci = amrex::coarsen(i, ratio[0]); -#if (AMREX_SPACEDIM >= 2) - int cj = amrex::coarsen(j, ratio[1]); -#else + +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(per_grown_domain); int cj = 0; +#else + int cj = amrex::coarsen(j, ratio[1]); #endif + #if (AMREX_SPACEDIM == 3) int ck = amrex::coarsen(k, ratio[2]); #else int ck = 0; #endif + if (dim == 0 && ci*ratio[0] == i) { // Check solve mask to ensure we don't overwrite valid fine data. if (!mask || mask(ci, cj, ck, n)) { From 7f7aa67e248fd967ed2cdd5b85cb60c20f94d193 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 11 Mar 2024 15:25:14 -0700 Subject: [PATCH 6/9] allow FaceConservativeLinearInterpolater even if DIM == 1 --- Src/AmrCore/AMReX_Interpolater.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Src/AmrCore/AMReX_Interpolater.cpp b/Src/AmrCore/AMReX_Interpolater.cpp index ec036b05399..532730c06b1 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -410,7 +410,6 @@ void FaceLinear::interp_arr (Array const& crse, }); } -#if (AMREX_SPACEDIM > 1) Box FaceConservativeLinear::CoarseBox (const Box& fine, int ratio) { @@ -712,7 +711,6 @@ void FaceConservativeLinear::interp_arr (Array const }); }); } -#endif // AMREX_SPACEDIM > 1 Box CellBilinear::CoarseBox (const Box& fine, int ratio) From 61f0e38a990dec6449aad743b4603cd7962d935e Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 12 Mar 2024 09:31:55 -0700 Subject: [PATCH 7/9] fixes based on PR feedback --- Src/AmrCore/AMReX_InterpFaceRegister.H | 2 ++ Src/AmrCore/AMReX_InterpFaceRegister.cpp | 10 +++++----- Src/AmrCore/AMReX_Interp_C.H | 12 ++++++------ Src/AmrCore/AMReX_Interpolater.cpp | 1 - 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/Src/AmrCore/AMReX_InterpFaceRegister.H b/Src/AmrCore/AMReX_InterpFaceRegister.H index c54879bcaf6..42e38df29c0 100644 --- a/Src/AmrCore/AMReX_InterpFaceRegister.H +++ b/Src/AmrCore/AMReX_InterpFaceRegister.H @@ -74,6 +74,8 @@ private: Geometry m_crse_geom; + Box per_grown_domain; + Array m_fine_face_ba; Array m_crse_face_ba; Array m_face_mask; // crse/fine: 1, fine/fine & fine/physbc: 0 diff --git a/Src/AmrCore/AMReX_InterpFaceRegister.cpp b/Src/AmrCore/AMReX_InterpFaceRegister.cpp index 33de90c15b4..858b2c0292a 100644 --- a/Src/AmrCore/AMReX_InterpFaceRegister.cpp +++ b/Src/AmrCore/AMReX_InterpFaceRegister.cpp @@ -21,7 +21,7 @@ void InterpFaceRegister::define (BoxArray const& fba, DistributionMapping const& m_crse_geom = amrex::coarsen(m_fine_geom, m_ref_ratio); // We don't need to worry about face-based domain because this is only used in the tangential interpolation - Box per_grown_domain = m_crse_geom.Domain(); + per_grown_domain = m_crse_geom.Domain(); for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { if (m_crse_geom.isPeriodic(dim)) { per_grown_domain.grow(dim,1); @@ -47,6 +47,7 @@ void InterpFaceRegister::define (BoxArray const& fba, DistributionMapping const& #ifdef AMREX_USE_GPU Vector > tags; #endif + #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -99,7 +100,6 @@ namespace { Array4 slope; Array4 crse; Array4 mask; - Box per_grown_domain; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box box() const noexcept { return Box(mask); } }; @@ -159,7 +159,7 @@ InterpFaceRegister::interp (Array const& fine, // NOL { if (tag.mask(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, tag.fine, scomp, tag.crse, - tag.slope, ncomp, tag.per_grown_domain, idim); + tag.slope, ncomp, per_grown_domain, idim); } }); } else @@ -187,7 +187,7 @@ InterpFaceRegister::interp (Array const& fine, // NOL { if (mlo_arr(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, fine_arr, scomp, clo_arr, - slope_arr, ncomp, domlo, idim); + slope_arr, ncomp, per_grown_domain, idim); } }); slope.resize(hibx, ncomp); @@ -196,7 +196,7 @@ InterpFaceRegister::interp (Array const& fine, // NOL { if (mhi_arr(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, fine_arr, scomp, chi_arr, - slope_arr, ncomp, domhi, idim); + slope_arr, ncomp, per_grown_domain, idim); } }); } diff --git a/Src/AmrCore/AMReX_Interp_C.H b/Src/AmrCore/AMReX_Interp_C.H index 0b2a2a74bb0..01e920680c0 100644 --- a/Src/AmrCore/AMReX_Interp_C.H +++ b/Src/AmrCore/AMReX_Interp_C.H @@ -115,7 +115,7 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, if (!mask || mask(ci, cj, ck, n)) { fine(i, j, k, n) = crse(ci, cj, ck, n); #if (AMREX_SPACEDIM >= 2) - if (cj != per_grown_domain.smallEnd(1) && cj != per_grown_domain.bigEnd(1) && ratio[1] > 1) { + if (cj > per_grown_domain.smallEnd(1) && cj < per_grown_domain.bigEnd(1) && ratio[1] > 1) { Real sfy = Real(1.0); Real dc = Real(0.5) * (crse(ci,cj+1,ck,n) - crse(ci,cj-1,ck,n)); Real df = Real(2.0) * (crse(ci,cj+1,ck,n) - crse(ci,cj ,ck,n)); @@ -131,7 +131,7 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, fine(i,j,k,n) += yoff * slope * sfy; } // jc #if (AMREX_SPACEDIM == 3) - if (ck != per_grown_domain.smallEnd(2) && ck != per_grown_domain.bigEnd(2) && ratio[2] > 1) { + if (ck > per_grown_domain.smallEnd(2) && ck < per_grown_domain.bigEnd(2) && ratio[2] > 1) { Real sfz = Real(1.0); Real dc = Real(0.5) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck-1,n)); Real df = Real(2.0) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck ,n)); @@ -156,7 +156,7 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, // Check solve mask to ensure we don't overwrite valid fine data. if (!mask || mask(ci, cj, ck, n)) { fine(i, j, k, n) = crse(ci, cj, ck, n); - if (ci != per_grown_domain.smallEnd(0) && ci != per_grown_domain.bigEnd(0) && ratio[0] > 1) { + if (ci > per_grown_domain.smallEnd(0) && ci < per_grown_domain.bigEnd(0) && ratio[0] > 1) { Real sfx = Real(1.0); Real dc = Real(0.5) * (crse(ci+1,cj,ck,n) - crse(ci-1,cj,ck,n)); Real df = Real(2.0) * (crse(ci+1,cj,ck,n) - crse(ci ,cj,ck,n)); @@ -172,7 +172,7 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, fine(i,j,k,n) += xoff * slope * sfx; } // ci #if (AMREX_SPACEDIM == 3) - if (ck != per_grown_domain.smallEnd(2) && ck != per_grown_domain.bigEnd(2) && ratio[2] > 1) { + if (ck > per_grown_domain.smallEnd(2) && ck < per_grown_domain.bigEnd(2) && ratio[2] > 1) { Real sfz = Real(1.0); Real dc = Real(0.5) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck-1,n)); Real df = Real(2.0) * (crse(ci,cj,ck+1,n) - crse(ci,cj,ck ,n)); @@ -197,7 +197,7 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, // Check solve mask to ensure we don't overwrite valid fine data. if (!mask || mask(ci, cj, ck, n)) { fine(i, j, k, n) = crse(ci, cj, ck, n); - if (ci != per_grown_domain.smallEnd(0) && ci != per_grown_domain.bigEnd(0) && ratio[0] > 1) { + if (ci > per_grown_domain.smallEnd(0) && ci < per_grown_domain.bigEnd(0) && ratio[0] > 1) { Real sfx = Real(1.0); Real dc = Real(0.5) * (crse(ci+1,cj,ck,n) - crse(ci-1,cj,ck,n)); Real df = Real(2.0) * (crse(ci+1,cj,ck,n) - crse(ci ,cj,ck,n)); @@ -212,7 +212,7 @@ face_cons_linear_face_interp (int i, int j, int k, int n, Array4 const& fine, Real xoff = (static_cast(i - ci*ratio[0]) + Real(0.5)) / Real(ratio[0]) - Real(0.5); fine(i,j,k,n) += xoff * slope * sfx; } // ci - if (cj != per_grown_domain.smallEnd(1) && cj != per_grown_domain.bigEnd(1) && ratio[1] > 1) { + if (cj > per_grown_domain.smallEnd(1) && cj < per_grown_domain.bigEnd(1) && ratio[1] > 1) { Real sfy = Real(1.0); Real dc = Real(0.5) * (crse(ci,cj+1,ck,n) - crse(ci,cj-1,ck,n)); Real df = Real(2.0) * (crse(ci,cj+1,ck,n) - crse(ci,cj ,ck,n)); diff --git a/Src/AmrCore/AMReX_Interpolater.cpp b/Src/AmrCore/AMReX_Interpolater.cpp index 532730c06b1..b5e855feb13 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -429,7 +429,6 @@ FaceConservativeLinear::CoarseBox (const Box& fine, const IntVect& ratio) for (int i = 0; i < AMREX_SPACEDIM; i++) { if (b.type(i) == IndexType::NODE) { - ng[i] = 0; if (b.type(i) == IndexType::NODE && b.length(i) < 2) { // Don't want degenerate boxes in nodal direction. b.growHi(i,1); From ab41da34468f54df55cbf87f96c1c8012d7e2680 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 12 Mar 2024 09:43:47 -0700 Subject: [PATCH 8/9] fix --- Src/AmrCore/AMReX_InterpFaceRegister.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/Src/AmrCore/AMReX_InterpFaceRegister.cpp b/Src/AmrCore/AMReX_InterpFaceRegister.cpp index 858b2c0292a..13a5ca33675 100644 --- a/Src/AmrCore/AMReX_InterpFaceRegister.cpp +++ b/Src/AmrCore/AMReX_InterpFaceRegister.cpp @@ -100,6 +100,7 @@ namespace { Array4 slope; Array4 crse; Array4 mask; + Box per_grown_domain; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box box() const noexcept { return Box(mask); } }; @@ -131,10 +132,6 @@ InterpFaceRegister::interp (Array const& fine, // NOL IntVect rr = m_ref_ratio; - Box const& domain = m_crse_geom.growPeriodicDomain(1); - Box const& domlo = amrex::bdryLo(domain, idim); - Box const& domhi = amrex::bdryHi(domain, idim); - #ifdef AMREX_USE_GPU if (Gpu::inLaunchRegion()) { auto const& mlo_mf = m_face_mask[Orientation(idim,Orientation::low)]; @@ -151,15 +148,15 @@ InterpFaceRegister::interp (Array const& fine, // NOL auto const& chi_arr = chidata.const_array(mfi); auto const& mlo_arr = mlo_mf.const_array(mfi); auto const& mhi_arr = mhi_mf.const_array(mfi); - tags.push_back(IFRTag{fine_arr, slo_arr, clo_arr, mlo_arr, domlo}); - tags.push_back(IFRTag{fine_arr, shi_arr, chi_arr, mhi_arr, domhi}); + tags.push_back(IFRTag{fine_arr, slo_arr, clo_arr, mlo_arr, per_grown_domain}); + tags.push_back(IFRTag{fine_arr, shi_arr, chi_arr, mhi_arr, per_grown_domain}); } ParallelFor(tags, [=] AMREX_GPU_DEVICE (int i, int j, int k, IFRTag const& tag) noexcept { if (tag.mask(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, tag.fine, scomp, tag.crse, - tag.slope, ncomp, per_grown_domain, idim); + tag.slope, ncomp, tags.per_grown_domain, idim); } }); } else From c41f8c9775d400eb03a7dd05242aba63243f5c4b Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 12 Mar 2024 10:40:20 -0700 Subject: [PATCH 9/9] revert to original versions of AMReX_InterpFaceReg* --- Src/AmrCore/AMReX_InterpFaceReg_2D_C.H | 6 +++--- Src/AmrCore/AMReX_InterpFaceReg_3D_C.H | 15 +++++++------- Src/AmrCore/AMReX_InterpFaceRegister.H | 2 -- Src/AmrCore/AMReX_InterpFaceRegister.cpp | 25 ++++++++++-------------- 4 files changed, 21 insertions(+), 27 deletions(-) diff --git a/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H b/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H index 9af617d4216..5b80958271d 100644 --- a/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H +++ b/Src/AmrCore/AMReX_InterpFaceReg_2D_C.H @@ -6,12 +6,12 @@ namespace amrex { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interp_face_reg (int i, int j, IntVect const& rr, Array4 const& fine, int scomp, Array4 const& crse, Array4 const& slope, int ncomp, - Box const& per_grown_domain, int idim) + Box const& domface, int idim) { int ic = amrex::coarsen(i,rr[0]); int jc = amrex::coarsen(j,rr[1]); if (idim == 0) { - if (jc == per_grown_domain.smallEnd(1) || jc == per_grown_domain.bigEnd(1)) { + if (jc == domface.smallEnd(1) || jc == domface.bigEnd(1)) { for (int n = 0; n < ncomp; ++n) { fine(i,j,0,n+scomp) = crse(ic,jc,0,n); } @@ -35,7 +35,7 @@ void interp_face_reg (int i, int j, IntVect const& rr, Array4 const& fine, } } } else { - if (ic == per_grown_domain.smallEnd(0) || ic == per_grown_domain.bigEnd(0)) { + if (ic == domface.smallEnd(0) || ic == domface.bigEnd(0)) { for (int n = 0; n < ncomp; ++n) { fine(i,j,0,n+scomp) = crse(ic,jc,0,n); } diff --git a/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H b/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H index 921395be47f..2df7fef055a 100644 --- a/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H +++ b/Src/AmrCore/AMReX_InterpFaceReg_3D_C.H @@ -6,7 +6,7 @@ namespace amrex { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const& fine, int scomp, Array4 const& crse, Array4 const& slope, int ncomp, - Box const& per_grown_domain, int idim) + Box const& domface, int idim) { int ic = amrex::coarsen(i,rr[0]); int jc = amrex::coarsen(j,rr[1]); @@ -15,7 +15,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } - if (jc != per_grown_domain.smallEnd(1) && jc != per_grown_domain.bigEnd(1) && rr[1] > 1) { + if (jc != domface.smallEnd(1) && jc != domface.bigEnd(1) && rr[1] > 1) { Real sfy = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)); @@ -34,7 +34,8 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const fine(i,j,k,n+scomp) += yoff * slope(i,j,k,n) * sfy; } } - if (kc != per_grown_domain.smallEnd(2) && kc != per_grown_domain.bigEnd(2) && rr[2] > 1) { + + if (kc != domface.smallEnd(2) && kc != domface.bigEnd(2) && rr[2] > 1) { Real sfz = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)); @@ -57,7 +58,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } - if (ic != per_grown_domain.smallEnd(0) && ic != per_grown_domain.bigEnd(0) && rr[0] > 1) { + if (ic != domface.smallEnd(0) && ic != domface.bigEnd(0) && rr[0] > 1) { Real sfx = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)); @@ -77,7 +78,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const } } - if (kc != per_grown_domain.smallEnd(2) && kc != per_grown_domain.bigEnd(2) && rr[2] > 1) { + if (kc != domface.smallEnd(2) && kc != domface.bigEnd(2) && rr[2] > 1) { Real sfz = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)); @@ -100,7 +101,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const for (int n = 0; n < ncomp; ++n) { fine(i,j,k,n+scomp) = crse(ic,jc,kc,n); } - if (ic != per_grown_domain.smallEnd(0) && ic != per_grown_domain.bigEnd(0) && rr[0] > 1) { + if (ic != domface.smallEnd(0) && ic != domface.bigEnd(0) && rr[0] > 1) { Real sfx = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)); @@ -120,7 +121,7 @@ void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4 const } } - if (jc != per_grown_domain.smallEnd(1) && jc != per_grown_domain.bigEnd(1) && rr[1] > 1) { + if (jc != domface.smallEnd(1) && jc != domface.bigEnd(1) && rr[1] > 1) { Real sfy = Real(1.0); for (int n = 0; n < ncomp; ++n) { Real dc = Real(0.5) * (crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)); diff --git a/Src/AmrCore/AMReX_InterpFaceRegister.H b/Src/AmrCore/AMReX_InterpFaceRegister.H index 42e38df29c0..c54879bcaf6 100644 --- a/Src/AmrCore/AMReX_InterpFaceRegister.H +++ b/Src/AmrCore/AMReX_InterpFaceRegister.H @@ -74,8 +74,6 @@ private: Geometry m_crse_geom; - Box per_grown_domain; - Array m_fine_face_ba; Array m_crse_face_ba; Array m_face_mask; // crse/fine: 1, fine/fine & fine/physbc: 0 diff --git a/Src/AmrCore/AMReX_InterpFaceRegister.cpp b/Src/AmrCore/AMReX_InterpFaceRegister.cpp index 13a5ca33675..ae5beb1985c 100644 --- a/Src/AmrCore/AMReX_InterpFaceRegister.cpp +++ b/Src/AmrCore/AMReX_InterpFaceRegister.cpp @@ -20,14 +20,6 @@ void InterpFaceRegister::define (BoxArray const& fba, DistributionMapping const& m_crse_geom = amrex::coarsen(m_fine_geom, m_ref_ratio); - // We don't need to worry about face-based domain because this is only used in the tangential interpolation - per_grown_domain = m_crse_geom.Domain(); - for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { - if (m_crse_geom.isPeriodic(dim)) { - per_grown_domain.grow(dim,1); - } - } - constexpr int crse_fine_face = 1; constexpr int fine_fine_face = 0; // First, set the face mask to 1 (i.e., coarse/fine boundary), @@ -47,7 +39,6 @@ void InterpFaceRegister::define (BoxArray const& fba, DistributionMapping const& #ifdef AMREX_USE_GPU Vector > tags; #endif - #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -100,7 +91,7 @@ namespace { Array4 slope; Array4 crse; Array4 mask; - Box per_grown_domain; + Box domface; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Box box() const noexcept { return Box(mask); } }; @@ -132,6 +123,10 @@ InterpFaceRegister::interp (Array const& fine, // NOL IntVect rr = m_ref_ratio; + Box const& domain = m_crse_geom.growPeriodicDomain(1); + Box const& domlo = amrex::bdryLo(domain, idim); + Box const& domhi = amrex::bdryHi(domain, idim); + #ifdef AMREX_USE_GPU if (Gpu::inLaunchRegion()) { auto const& mlo_mf = m_face_mask[Orientation(idim,Orientation::low)]; @@ -148,15 +143,15 @@ InterpFaceRegister::interp (Array const& fine, // NOL auto const& chi_arr = chidata.const_array(mfi); auto const& mlo_arr = mlo_mf.const_array(mfi); auto const& mhi_arr = mhi_mf.const_array(mfi); - tags.push_back(IFRTag{fine_arr, slo_arr, clo_arr, mlo_arr, per_grown_domain}); - tags.push_back(IFRTag{fine_arr, shi_arr, chi_arr, mhi_arr, per_grown_domain}); + tags.push_back(IFRTag{fine_arr, slo_arr, clo_arr, mlo_arr, domlo}); + tags.push_back(IFRTag{fine_arr, shi_arr, chi_arr, mhi_arr, domhi}); } ParallelFor(tags, [=] AMREX_GPU_DEVICE (int i, int j, int k, IFRTag const& tag) noexcept { if (tag.mask(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, tag.fine, scomp, tag.crse, - tag.slope, ncomp, tags.per_grown_domain, idim); + tag.slope, ncomp, tag.domface, idim); } }); } else @@ -184,7 +179,7 @@ InterpFaceRegister::interp (Array const& fine, // NOL { if (mlo_arr(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, fine_arr, scomp, clo_arr, - slope_arr, ncomp, per_grown_domain, idim); + slope_arr, ncomp, domlo, idim); } }); slope.resize(hibx, ncomp); @@ -193,7 +188,7 @@ InterpFaceRegister::interp (Array const& fine, // NOL { if (mhi_arr(i,j,k)) { interp_face_reg(AMREX_D_DECL(i,j,k), rr, fine_arr, scomp, chi_arr, - slope_arr, ncomp, per_grown_domain, idim); + slope_arr, ncomp, domhi, idim); } }); }