Fix Bug in FaceLinear Interpolater

The bug was introduced in AMReX-Codes#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

@@ -164,6 +164,36 @@ facediv_int (int ci, int cj, int /*ck*/, int nf,
 
 }
 
+template<typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
+face_linear_interp_x (int i, int j, int /*k*/, int n, Array4<T> const& fine,
+ Array4<T const> 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<Real>(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<typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
+face_linear_interp_y (int i, int j, int /*k*/, int n, Array4<T> const& fine,
+ Array4<T const> 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<Real>(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 <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void ccprotect_2d (int ic, int jc, int /*kc*/, int nvar,

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<typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
+face_linear_interp_x (int i, int j, int k, int n, Array4<T> const& fine,
+ Array4<T const> 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<Real>(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<typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
+face_linear_interp_y (int i, int j, int k, int n, Array4<T> const& fine,
+ Array4<T const> 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<Real>(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<typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
+face_linear_interp_z (int i, int j, int k, int n, Array4<T> const& fine,
+ Array4<T const> 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<Real>(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 <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void ccprotect_3d (int ic, int jc, int kc, int nvar,

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<Real> const& fine_arr = fine.array(fine_comp);
+ Array4<Real const> 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