From 10d2a964e5ce66453547c6c01621acc91e1a5951 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Sun, 6 Aug 2023 16:29:01 -0700 Subject: [PATCH] Fix Bug in FaceLinear Interpolater The bug was introduced in #2539. The issue was when FaceLinear was used by functions other than FillPatchTwoLevels, the assumption of the fine box was a refined version of the coarse box could be wrong. The bug did not affect FillPatchFromTwoLevels. --- Src/AmrCore/AMReX_Interp_2D_C.H | 30 +++++++++++++++++++ Src/AmrCore/AMReX_Interp_3D_C.H | 48 ++++++++++++++++++++++++++++++ Src/AmrCore/AMReX_Interpolater.cpp | 41 ++++++++++++++++++++----- 3 files changed, 112 insertions(+), 7 deletions(-) diff --git a/Src/AmrCore/AMReX_Interp_2D_C.H b/Src/AmrCore/AMReX_Interp_2D_C.H index 40bc1b2ab0e..20f8b1c3a0f 100644 --- a/Src/AmrCore/AMReX_Interp_2D_C.H +++ b/Src/AmrCore/AMReX_Interp_2D_C.H @@ -164,6 +164,36 @@ facediv_int (int ci, int cj, int /*ck*/, int nf, } +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_x (int i, int j, int /*k*/, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + const int ii = amrex::coarsen(i,ratio[0]); + const int jj = amrex::coarsen(j,ratio[1]); + if (i-ii*ratio[0] == 0) { + fine(i,j,0,n) = crse(ii,jj,0,n); + } else { + Real const w = static_cast(i-ii*ratio[0]) * (Real(1.)/Real(ratio[0])); + fine(i,j,0,n) = (Real(1.)-w) * crse(ii,jj,0,n) + w * crse(ii+1,jj,0,n); + } +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_y (int i, int j, int /*k*/, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + const int ii = amrex::coarsen(i,ratio[0]); + const int jj = amrex::coarsen(j,ratio[1]); + if (j-jj*ratio[1] == 0) { + fine(i,j,0,n) = crse(ii,jj,0,n); + } else { + Real const w = static_cast(j-jj*ratio[1]) * (Real(1.)/Real(ratio[1])); + fine(i,j,0,n) = (Real(1.)-w) * crse(ii,jj,0,n) + w * crse(ii,jj+1,0,n); + } +} + template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ccprotect_2d (int ic, int jc, int /*kc*/, int nvar, diff --git a/Src/AmrCore/AMReX_Interp_3D_C.H b/Src/AmrCore/AMReX_Interp_3D_C.H index 510e0dabdfe..daa19f6d5e3 100644 --- a/Src/AmrCore/AMReX_Interp_3D_C.H +++ b/Src/AmrCore/AMReX_Interp_3D_C.H @@ -332,6 +332,54 @@ facediv_int (int ci, int cj, int ck, int nf, + dz3/(8*dy*zspxs)*(v000+v021-v001-v020); } +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_x (int i, int j, int k, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + const int ii = amrex::coarsen(i,ratio[0]); + const int jj = amrex::coarsen(j,ratio[1]); + const int kk = amrex::coarsen(k,ratio[2]); + if (i-ii*ratio[0] == 0) { + fine(i,j,k,n) = crse(ii,jj,kk,n); + } else { + Real const w = static_cast(i-ii*ratio[0]) * (Real(1.)/Real(ratio[0])); + fine(i,j,k,n) = (Real(1.)-w) * crse(ii,jj,kk,n) + w * crse(ii+1,jj,kk,n); + } +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_y (int i, int j, int k, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + const int ii = amrex::coarsen(i,ratio[0]); + const int jj = amrex::coarsen(j,ratio[1]); + const int kk = amrex::coarsen(k,ratio[2]); + if (j-jj*ratio[1] == 0) { + fine(i,j,k,n) = crse(ii,jj,kk,n); + } else { + Real const w = static_cast(j-jj*ratio[1]) * (Real(1.)/Real(ratio[1])); + fine(i,j,k,n) = (Real(1.)-w) * crse(ii,jj,kk,n) + w * crse(ii,jj+1,kk,n); + } +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_z (int i, int j, int k, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + const int ii = amrex::coarsen(i,ratio[0]); + const int jj = amrex::coarsen(j,ratio[1]); + const int kk = amrex::coarsen(k,ratio[2]); + if (k-kk*ratio[2] == 0) { + fine(i,j,k,n) = crse(ii,jj,kk,n); + } else { + Real const w = static_cast(k-kk*ratio[2]) * (Real(1.)/Real(ratio[2])); + fine(i,j,k,n) = (Real(1.)-w) * crse(ii,jj,kk,n) + w * crse(ii,jj,kk+1,n); + } +} + template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ccprotect_3d (int ic, int jc, int kc, int nvar, diff --git a/Src/AmrCore/AMReX_Interpolater.cpp b/Src/AmrCore/AMReX_Interpolater.cpp index a5947f18aab..601b8b4b861 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -134,10 +134,10 @@ FaceLinear::interp (const FArrayBox& crse, int ncomp, const Box& fine_region, const IntVect& ratio, - const Geometry& crse_geom , - const Geometry& fine_geom , - Vector const& bcr, - int actual_comp, + const Geometry& /* crse_geom */, + const Geometry& /* fine_geom */, + Vector const& /*bcr*/, + int /* actual_comp*/, int /*actual_state*/, RunOn runon) { @@ -146,9 +146,36 @@ FaceLinear::interp (const FArrayBox& crse, // BL_PROFILE("FaceLinear::interp()"); - // pass unallocated IArrayBox for solve_mask, so all fine values get filled. - interp_face(crse, crse_comp, fine, fine_comp, ncomp, fine_region, - ratio, IArrayBox(), crse_geom, fine_geom, bcr, actual_comp, runon); + 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); + + 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,crse_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,crse_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,crse_arr,ratio); + }); + } +#endif +#endif } void