From 979170a917a4ae5c921b6835e1f18e609b7b98a6 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 | 33 ++++++++++++++++++-- 3 files changed, 108 insertions(+), 3 deletions(-) diff --git a/Src/AmrCore/AMReX_Interp_2D_C.H b/Src/AmrCore/AMReX_Interp_2D_C.H index 40bc1b2ab0e..60ab3c4e06f 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.)/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.)/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..439ce86c219 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.)/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.)/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.)/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..7a78b3f305e 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -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