Fix Bug in FaceConservativeLinear (#4630)

This fixes a bug when the fine box is not coarsenable in the face
direction.

Author: WeiqunZhang

Src/AmrCore/AMReX_Interp_C.H

@@ -283,6 +283,60 @@ face_linear_interp_z (int i, int j, int k, int n, amrex::Array4<amrex::Real> con
     }
 }
 
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+Real face_linear_interp_safe_x (int i, int j, int k, int n,
+                                amrex::Array4<amrex::Real> const& fine,
+                                IntVect const& ratio) noexcept
+{
+    const int ci = amrex::coarsen(i,ratio[0]);
+
+    if (i-ci*ratio[0] != 0)
+    {
+        Real const w = static_cast<Real>(i-ci*ratio[0]) * (Real(1.)/Real(ratio[0]));
+        int i1 = ci*ratio[0];
+        int i2 = (ci+1)*ratio[0];
+        return (Real(1.)-w) * fine(i1,j,k,n) + w * fine(i2,j,k,n);
+    } else {
+        return fine(i,j,k,n);
+    }
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+Real face_linear_interp_safe_y (int i, int j, int k, int n,
+                                amrex::Array4<amrex::Real> const& fine,
+                                IntVect const& ratio) noexcept
+{
+    const int cj = amrex::coarsen(j,ratio[1]);
+
+    if (j-cj*ratio[1] != 0)
+    {
+        Real const w = static_cast<Real>(j-cj*ratio[1]) * (Real(1.)/Real(ratio[1]));
+        int j1 = cj*ratio[1];
+        int j2 = (cj+1)*ratio[1];
+        return (Real(1.)-w) * fine(i,j1,k,n) + w * fine(i,j2,k,n);
+    } else {
+        return fine(i,j,k,n);
+    }
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+Real face_linear_interp_safe_z (int i, int j, int k, int n,
+                                amrex::Array4<amrex::Real> const& fine,
+                                IntVect const& ratio) noexcept
+{
+    const int ck = amrex::coarsen(k,ratio[2]);
+
+    if (k-ck*ratio[2] != 0)
+    {
+        Real const w = static_cast<Real>(k-ck*ratio[2]) * (Real(1.)/Real(ratio[2]));
+        int k1 = ck*ratio[2];
+        int k2 = (ck+1)*ratio[2];
+        return (Real(1.)-w) * fine(i,j,k1,n) + w * fine(i,j,k2,n);
+    } else {
+        return fine(i,j,k,n);
+    }
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void cell_quartic_interp_x (int i, int j, int k, int n, Array4<Real> const& fine,
                             Array4<Real const> const& crse) noexcept

Src/AmrCore/AMReX_Interpolater.cpp

@@ -516,6 +516,24 @@ FaceConservativeLinear::interp_face (const FArrayBox&       crse,
         }
     }
 
+    bool is_safe = true;
+    FArrayBox safe_fine;
+    Array4<Real> safe_fine_arr = fine_arr;
+    Box safe_fine_region = fine_region;
+    int facedir = 0;
+    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+        if (fine_region.type(idim) == IndexType::NODE) { facedir = idim; }
+    }
+    IntVect rrtmp(1);
+    rrtmp[facedir] = ratio[facedir];
+    if (! safe_fine_region.coarsenable(rrtmp)) {
+        is_safe = false;
+        safe_fine_region.coarsen(rrtmp);
+        safe_fine_region.refine(rrtmp);
+        safe_fine.resize(safe_fine_region, ncomp, The_Async_Arena());
+        safe_fine_arr = safe_fine.array();
+    }
+
     //
     // Fill fine ghost faces with interpolation of coarse data that is conservative linear
     //      in the tangential direction.
@@ -525,25 +543,25 @@ FaceConservativeLinear::interp_face (const FArrayBox&       crse,
     //
     if (fine_region.type(0) == IndexType::NODE)
     {
-        AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n,
+        AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,safe_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);
+            face_cons_linear_face_interp(i,j,k,n,safe_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,
+        AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,safe_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);
+            face_cons_linear_face_interp(i,j,k,n,safe_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,
+        AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,safe_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);
+            face_cons_linear_face_interp(i,j,k,n,safe_fine_arr,crse_arr,mask_arr,ratio,per_grown_domain,2);
         });
     }
 #endif
@@ -556,26 +574,47 @@ FaceConservativeLinear::interp_face (const FArrayBox&       crse,
     //
     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 (is_safe) {
+            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);
+            });
+        } else {
+            AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n,
+            {
+                fine_arr(i,j,k,n) = face_linear_interp_safe_x(i,j,k,n,safe_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 (is_safe) {
+            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);
+            });
+        } else {
+            AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n,
+            {
+                fine_arr(i,j,k,n) = face_linear_interp_safe_y(i,j,k,n,safe_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);
-        });
+        if (is_safe) {
+            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);
+            });
+        } else {
+            AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n,
+            {
+                fine_arr(i,j,k,n) = face_linear_interp_safe_z(i,j,k,n,safe_fine_arr,ratio);
+            });
+        }
     }
 #endif
 #endif