added general divergence preserving interpolation of face centered fields

Src/AmrCore/AMReX_Interp_1D_C.H

@@ -73,7 +73,7 @@ nodebilin_interp (Box const& bx, Array4<T> const& fine, const int fcomp, const i
 template<typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
 facediv_face_interp (int /*ci*/, int /*cj*/, int /*ck*/,
-                     int /*nc*/, int /*nf*/, int /*idir*/,
+                     int /*n*/, int /*idir*/,
                      Array4<T const> const& /*crse*/,
                      Array4<T> const& /*fine*/,
                      Array4<const int> const& /*mask*/,

Src/AmrCore/AMReX_Interp_2D_C.H

@@ -4,6 +4,7 @@
 
 #include <AMReX_Array.H>
 #include <AMReX_Geometry.H>
+#include <AMReX_LUSolver.H>
 
 namespace amrex {
 
@@ -97,14 +98,14 @@ nodebilin_interp (Box const& bx, Array4<T> const& fine, const int fcomp, const i
 template<typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
 facediv_face_interp (int ci, int cj, int /*ck*/,
-                     int nc, int nf, int idir,
+                     int n, int idir,
                      Array4<T const> const& crse,
                      Array4<T> const& fine,
                      Array4<int const> const& mask,
                      IntVect const& ratio) noexcept
 {
     if (mask) {
-        if (!mask(ci, cj, 0, nc))
+        if (!mask(ci, cj, 0, n))
             { return; }
     }
 
@@ -114,23 +115,25 @@ facediv_face_interp (int ci, int cj, int /*ck*/,
     switch (idir) {
         case 0:
         {
-            const Real neg = crse(ci, cj-1, 0, nc);
-            const Real cen = crse(ci, cj  , 0, nc);
-            const Real pos = crse(ci, cj+1, 0, nc);
+            const Real slpy = Real(0.5)*(crse(ci, cj+1, 0, n) - crse(ci, cj-1, 0, n));
+            const Real dys = Real(1.)/Real(ratio[1]);
 
-            fine(fi, fj  , 0, nf) = Real(0.125)*(8*cen + neg - pos);
-            fine(fi, fj+1, 0, nf) = Real(0.125)*(8*cen + pos - neg);
+            for (int j = 0; j < ratio[1]; ++j) {
+                Real ysloc = (Real(j)+Real(0.5))*dys - Real(0.5);
+                fine(fi, fj+j, 0, n) = crse(ci, cj, 0, n) + ysloc*slpy;
+            }
 
             break;
         }
         case 1:
         {
-            const Real neg = crse(ci-1, cj, 0, nc);
-            const Real cen = crse(ci  , cj, 0, nc);
-            const Real pos = crse(ci+1, cj, 0, nc);
+            const Real slpx = Real(0.5)*(crse(ci+1, cj, 0, n) - crse(ci-1, cj, 0, n));
+            const Real dxs = Real(1.)/Real(ratio[0]);
 
-            fine(fi  , fj, 0, nf) = Real(0.125)*(8*cen + neg - pos);
-            fine(fi+1, fj, 0, nf) = Real(0.125)*(8*cen + pos - neg);
+            for (int i = 0; i < ratio[0]; ++i) {
+                Real xsloc = (Real(i)+Real(0.5))*dxs - Real(0.5);
+                fine(fi+i, fj, 0, n) = crse(ci, cj, 0, n) + xsloc*slpx;
+            }
 
             break;
         }
@@ -140,7 +143,7 @@ facediv_face_interp (int ci, int cj, int /*ck*/,
 
 template<typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
-facediv_int (int ci, int cj, int /*ck*/, int nf,
+facediv_int (int ci, int cj, int /*ck*/, int n,
              GpuArray<Array4<T>, AMREX_SPACEDIM> const& fine,
              IntVect const& ratio,
              GpuArray<Real, AMREX_SPACEDIM> const& cellSize) noexcept
@@ -149,26 +152,167 @@ facediv_int (int ci, int cj, int /*ck*/, int nf,
     const int fj = cj*ratio[1];
 
     // References to fine exterior values.
-    const Real umm = fine[0](fi,   fj,   0, nf);
-    const Real ump = fine[0](fi,   fj+1, 0, nf);
-    const Real upm = fine[0](fi+2, fj,   0, nf);
-    const Real upp = fine[0](fi+2, fj+1, 0, nf);
+    const Real umm = fine[0](fi,   fj,   0, n);
+    const Real ump = fine[0](fi,   fj+1, 0, n);
+    const Real upm = fine[0](fi+2, fj,   0, n);
+    const Real upp = fine[0](fi+2, fj+1, 0, n);
 
-    const Real vmm = fine[1](fi,   fj,   0, nf);
-    const Real vmp = fine[1](fi+1, fj,   0, nf);
-    const Real vpm = fine[1](fi,   fj+2, 0, nf);
-    const Real vpp = fine[1](fi+1, fj+2, 0, nf);
+    const Real vmm = fine[1](fi,   fj,   0, n);
+    const Real vmp = fine[1](fi+1, fj,   0, n);
+    const Real vpm = fine[1](fi,   fj+2, 0, n);
+    const Real vpp = fine[1](fi+1, fj+2, 0, n);
 
     const Real dxdy = cellSize[0]/cellSize[1];
     const Real x_corr = Real(0.25)*dxdy * (vpp+vmm-vmp-vpm);
     const Real y_corr = Real(0.25)/dxdy * (upp+umm-ump-upm);
 
     // Calc fine faces on interior of coarse cells.
-    fine[0](fi+1,fj  ,0,nf) = Real(0.5)*(umm+upm) + x_corr;
-    fine[0](fi+1,fj+1,0,nf) = Real(0.5)*(ump+upp) + x_corr;
-    fine[1](fi,  fj+1,0,nf) = Real(0.5)*(vmm+vpm) + y_corr;
-    fine[1](fi+1,fj+1,0,nf) = Real(0.5)*(vmp+vpp) + y_corr;
+    fine[0](fi+1,fj  ,0,n) = Real(0.5)*(umm+upm) + x_corr;
+    fine[0](fi+1,fj+1,0,n) = Real(0.5)*(ump+upp) + x_corr;
+    fine[1](fi,  fj+1,0,n) = Real(0.5)*(vmm+vpm) + y_corr;
+    fine[1](fi+1,fj+1,0,n) = Real(0.5)*(vmp+vpp) + y_corr;
+
+}
+
+template <int i, std::enable_if_t<(i>=0 && i<4), int> = 0>
+constexpr IntVect facediv_get_refratio ()
+{
+    if constexpr (i == 0) {
+        return IntVect(2,2);
+    } else if constexpr (i == 1) {
+        return IntVect(4,2);
+    } else if constexpr (i == 2) {
+        return IntVect(2,4);
+    } else {
+        return IntVect(4,4);
+    }
+}
+
+template <int N>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+LUSolver<N,Real>
+facediv_init_lusolver (IntVect const& ratio,
+                       GpuArray<Real, AMREX_SPACEDIM> const& cellSize)
+{
+    Array2D<Real,0,N-1,0,N-1,Order::C> mat;
+
+    for (int irow = 0; irow < N; ++irow) {
+        for (int icol = 0; icol < N; ++icol) {
+            mat(irow,icol) = Real(0.0);
+        }
+    }
+
+    Real dxinv = Real(1.)/cellSize[0];
+    Real dyinv = Real(1.)/cellSize[1];
+    Real ddx2 = dxinv*dxinv;
+    Real ddy2 = dyinv*dyinv;
+
+    for (int j = 0; j < ratio[1]; ++j) {
+    for (int i = 0; i < ratio[0]; ++i) {
+        int irow = i + j*ratio[0];
+
+        if (i > 0) {
+            mat(irow,irow) += ddx2;
+            mat(irow,irow-1) += -ddx2;
+        }
+
+        if (i < ratio[0]-1) {
+            mat(irow,irow) += ddx2;
+            mat(irow,irow+1) += -ddx2;
+        }
+
+        if (j > 0) {
+            mat(irow,irow) += ddy2;
+            mat(irow,irow-ratio[0]) += -ddy2;
+        }
+
+        if (j < ratio[1]-1) {
+            mat(irow,irow) += ddy2;
+            mat(irow,irow+ratio[0]) += -ddy2;
+        }
+    }}
+
+    // Hack matrix to be nonsingular. This just picks a particular gauge
+    // with the last row being zero, doesn’t change the gradients we will
+    // apply.
+    mat(N-1,N-1) += ddx2;
+
+    return LUSolver<N,Real>{mat};
+}
+
+template<typename T, int N>
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
+facediv_int_gen (int ci, int cj, int /*ck*/, int nf,
+                 GpuArray<Array4<T>, AMREX_SPACEDIM> const& fine,
+                 IntVect const& ratio,
+                 GpuArray<Real, AMREX_SPACEDIM> const& cellSize,
+                 LUSolver<N,Real> const& lusolver) noexcept
+{
+    const int fi = ci*ratio[0];
+    const int fj = cj*ratio[1];
+
+    const Real dx = cellSize[0];
+    const Real dy = cellSize[1];
+    const Real dxinv = Real(1.)/dx;
+    const Real dyinv = Real(1.)/dy;
+
+    for (int j = 0; j < ratio[1]; ++j) {
+    for (int i = 1; i < ratio[0]; ++i) {
+        fine[0](fi+i, fj+j, 0, nf) =
+            ((ratio[0]-i)*fine[0](fi,           fj+j, 0, nf)
+           +          i *fine[0](fi+ratio[0], fj+j, 0, nf))/ratio[0];
+    }}
+
+    for (int j = 1; j < ratio[1]; ++j) {
+    for (int i = 0; i < ratio[0]; ++i) {
+        fine[1](fi+i, fj+j, 0, nf) =
+            ((ratio[1]-j)*fine[1](fi+i, fj,           0, nf)
+           +          j *fine[1](fi+i, fj+ratio[1], 0, nf))/ratio[1];
+    }}
+
+    Real uplus = Real(0.);
+    Real uminus = Real(0.);
+    for (int j = 0; j < ratio[1]; ++j) {
+        uplus  += fine[0](fi+ratio[0], fj+j, 0, nf);
+        uminus += fine[0](fi,          fj+j, 0, nf);
+    }
+
+    Real vplus = Real(0.);
+    Real vminus = Real(0.);
+    for (int i = 0; i < ratio[0]; ++i) {
+        vplus  += fine[1](fi+i, fj+ratio[1], 0, nf);
+        vminus += fine[1](fi+i, fj,          0, nf);
+    }
 
+    Real target_div = (uplus-uminus)/(N*dx) + (vplus-vminus)/(N*dy);
+
+    Array1D<Real,0,N-1> rhs;
+    Array1D<Real,0,N-1> phi;
+
+    for (int j = 0; j < ratio[1]; ++j) {
+    for (int i = 0; i < ratio[0]; ++i) {
+        int irow = i + j*ratio[0];
+        rhs(irow) =
+            (fine[0](fi+i+1, fj+j,   0, nf) - fine[0](fi+i,   fj+j,   0, nf))*dxinv
+          + (fine[1](fi+i,   fj+j+1, 0, nf) - fine[1](fi+i,   fj+j,   0, nf))*dyinv
+          - target_div;
+    }}
+
+    lusolver(phi.begin(), rhs.begin());
+
+    for (int j = 0; j < ratio[1]; ++j) {
+    for (int i = 1; i < ratio[0]; ++i) {
+        int plus  = i   + j*ratio[0];
+        int minus = i-1 + j*ratio[0];
+        fine[0](fi+i, fj+j, 0, nf) += (phi(plus) - phi(minus))*dxinv;
+    }}
+
+    for (int j = 1; j < ratio[1]; ++j) {
+    for (int i = 0; i < ratio[0]; ++i) {
+        int plus  = i + j   *ratio[0];
+        int minus = i + (j-1)*ratio[0];
+        fine[1](fi+i, fj+j, 0, nf) += (phi(plus) - phi(minus))*dyinv;
+    }}
 }
 
 template<typename T>