Quartic interpolation for cell centered data

New Interpolator for interpolation of cell centered data using a
fourth-degreee polynomial.  Note that the interpolation is not conservative
and does not do any slope limiting.

Author: WeiqunZhang

Src/AmrCore/AMReX_Interp_C.H

@@ -135,5 +135,53 @@ face_linear_interp_z (int i, int j, int k, int n, amrex::Array4<amrex::Real> con
     }
 }
 
+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
+{
+    constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
+                                      Real(0.92285156), Real(0.20507812),
+                                      Real(-0.02197266)};
+    int ii = amrex::coarsen(i,2);
+    int s = 2*(i-ii*2) - 1;  // if i == ii*2, s = -1; if i == ii*2+1, s = 1;
+    fine(i,j,k,n) = c[-2*s]*crse(ii-2,j,k,n)
+        +           c[  -s]*crse(ii-1,j,k,n)
+        +           c[   0]*crse(ii  ,j,k,n)
+        +           c[   s]*crse(ii+1,j,k,n)
+        +           c[ 2*s]*crse(ii+2,j,k,n);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void cell_quartic_interp_y (int i, int j, int k, int n, Array4<Real> const& fine,
+                            Array4<Real const> const& crse) noexcept
+{
+    constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
+                                      Real(0.92285156), Real(0.20507812),
+                                      Real(-0.02197266)};
+    int jj = amrex::coarsen(j,2);
+    int s = 2*(j-jj*2) - 1;  // if j == jj*2, s = -1; if j == jj*2+1, s = 1;
+    fine(i,j,k,n) = c[-2*s]*crse(i,jj-2,k,n)
+        +           c[  -s]*crse(i,jj-1,k,n)
+        +           c[   0]*crse(i,jj  ,k,n)
+        +           c[   s]*crse(i,jj+1,k,n)
+        +           c[ 2*s]*crse(i,jj+2,k,n);
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void cell_quartic_interp_z (int i, int j, int k, int n, Array4<Real> const& fine,
+                            Array4<Real const> const& crse) noexcept
+{
+    constexpr Array1D<Real,-2,2> c = {Real(0.01708984), Real(-0.12304688),
+                                      Real(0.92285156), Real(0.20507812),
+                                      Real(-0.02197266)};
+    int kk = amrex::coarsen(k,2);
+    int s = 2*(k-kk*2) - 1;  // if k == kk*2, s = -1; if k == kk*2+1, s = 1;
+    fine(i,j,k,n) = c[-2*s]*crse(i,j,kk-2,n)
+        +           c[  -s]*crse(i,j,kk-1,n)
+        +           c[   0]*crse(i,j,kk  ,n)
+        +           c[   s]*crse(i,j,kk+1,n)
+        +           c[ 2*s]*crse(i,j,kk+2,n);
+}
+
 }
 #endif

Src/AmrCore/AMReX_Interpolater.H

@@ -844,6 +844,74 @@ public:
 
 };
 
+/**
+* \brief Quartic interpolation on cell centered data.
+*
+* Quartic interpolation on cell centered data.
+*/
+
+class CellQuartic
+    :
+    public Interpolater
+{
+public:
+
+    /**
+    * \brief The constructor.
+    */
+    explicit CellQuartic ();
+
+    /**
+    * \brief The destructor.
+    */
+    virtual ~CellQuartic () override;
+
+    /**
+    * \brief Returns coarsened box given fine box and refinement ratio.
+    *
+    * \param fine
+    * \param ratio
+    */
+    virtual Box CoarseBox (const Box& fine, int ratio) override;
+
+    /**
+    * \brief Returns coarsened box given fine box and refinement ratio.
+    *
+    * \param fine
+    * \param ratio
+    */
+    virtual Box CoarseBox (const Box& fine, const IntVect& ratio) override;
+
+    /**
+    * \brief Coarse to fine interpolation in space.
+    *
+    * \param crse
+    * \param crse_comp
+    * \param fine
+    * \param fine_comp
+    * \param ncomp
+    * \param fine_region
+    * \param ratio
+    * \param crse_geom
+    * \param fine_geom
+    * \param bcr
+    * \param actual_comp
+    * \param actual_state
+    */
+    virtual void interp (const FArrayBox& crse,
+                         int              crse_comp,
+                         FArrayBox&       fine,
+                         int              fine_comp,
+                         int              ncomp,
+                         const Box&       fine_region,
+                         const IntVect&   ratio,
+                         const Geometry&  crse_geom,
+                         const Geometry&  fine_geom,
+                         Vector<BCRec> const& bcr,
+                         int              actual_comp,
+                         int              actual_state,
+                         RunOn            gpu_or_cpu) override;
+};
 
 //! CONSTRUCT A GLOBAL OBJECT OF EACH VERSION.
 extern AMREX_EXPORT PCInterp                  pc_interp;
@@ -856,6 +924,7 @@ extern AMREX_EXPORT CellBilinear              cell_bilinear_interp;
 extern AMREX_EXPORT CellConservativeProtected protected_interp;
 extern AMREX_EXPORT CellConservativeQuartic   quartic_interp;
 extern AMREX_EXPORT CellQuadratic             quadratic_interp;
+extern AMREX_EXPORT CellQuartic               cell_quartic_interp;
 
 }
 

Src/AmrCore/AMReX_Interpolater.cpp

@@ -18,6 +18,8 @@ namespace amrex {
  *
  * CellQuadratic only works in 2D and 3D on cpu and gpu.
  *
+ * CellQuartic works in 1D, 2D and 3D on cpu and gpu with ref ratio of 2
+ *
  * CellConservativeQuartic only works with ref ratio of 2 on cpu and gpu.
  *
  * FaceDivFree works in 2D and 3D on cpu and gpu.
@@ -37,6 +39,7 @@ CellConservativeProtected protected_interp;
 CellConservativeQuartic   quartic_interp;
 CellBilinear              cell_bilinear_interp;
 CellQuadratic             quadratic_interp;
+CellQuartic               cell_quartic_interp;
 
 NodeBilinear::~NodeBilinear () {}
 
@@ -988,4 +991,91 @@ FaceDivFree::interp_arr (Array<FArrayBox*, AMREX_SPACEDIM> const& crse,
     });
 }
 
+CellQuartic::CellQuartic () {}
+
+CellQuartic::~CellQuartic () {}
+
+Box
+CellQuartic::CoarseBox (const Box& fine, const IntVect& ratio)
+{
+    Box crse = amrex::coarsen(fine,ratio);
+    crse.grow(2);
+    return crse;
+}
+
+Box
+CellQuartic::CoarseBox (const Box& fine, int ratio)
+{
+    Box crse = amrex::coarsen(fine,ratio);
+    crse.grow(2);
+    return crse;
+}
+
+void
+CellQuartic::interp (const FArrayBox& crse,
+                     int              crse_comp,
+                     FArrayBox&       fine,
+                     int              fine_comp,
+                     int              ncomp,
+                     const Box&       fine_region,
+                     const IntVect&   ratio,
+                     const Geometry&  /*crse_geom*/,
+                     const Geometry&  /*fine_geom*/,
+                     Vector<BCRec> const&  /*bcr*/,
+                     int              /* actual_comp */,
+                     int              /* actual_state */,
+                     RunOn            runon)
+{
+    BL_PROFILE("CellQuartic::interp()");
+    amrex::ignore_unused(ratio);
+    AMREX_ASSERT(ratio == 2);
+
+    Box target_fine_region = fine_region & fine.box();
+
+    bool run_on_gpu = (runon == RunOn::Gpu && Gpu::inLaunchRegion());
+
+    Array4<Real const> const& crsearr = crse.const_array(crse_comp);
+    Array4<Real>       const& finearr = fine.array(fine_comp);
+
+#if (AMREX_SPACEDIM == 3)
+    Box bz = amrex::coarsen(target_fine_region, IntVect(2,2,1));
+    bz.grow(IntVect(2,2,0));
+    FArrayBox tmpz(bz, ncomp);
+    Elixir tmpz_eli = tmpz.elixir();
+    Array4<Real> const& tmpzarr = tmpz.array();
+    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon, bz, ncomp, i, j, k, n,
+    {
+        cell_quartic_interp_z(i,j,k,n,tmpzarr,crsearr);
+    });
+#endif
+
+#if (AMREX_SPACEDIM >= 2)
+    Box by = amrex::coarsen(target_fine_region, IntVect(AMREX_D_DECL(2,1,1)));
+    by.grow(IntVect(AMREX_D_DECL(2,0,0)));
+    FArrayBox tmpy(by, ncomp);
+    Elixir tmpy_eli = tmpy.elixir();
+    Array4<Real> const& tmpyarr = tmpy.array();
+#if (AMREX_SPACEDIM == 2)
+    Array4<Real const> srcarr = crsearr;
+#else
+    Array4<Real const> srcarr = tmpz.const_array();
+#endif
+    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon, by, ncomp, i, j, k, n,
+    {
+        cell_quartic_interp_y(i,j,k,n,tmpyarr,srcarr);
+    });
+#endif
+
+#if (AMREX_SPACEDIM == 1)
+    Array4<Real const> srcarr = crsearr;
+#else
+    srcarr = tmpy.const_array();
+#endif
+    AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon, target_fine_region, ncomp,
+                                           i, j, k, n,
+    {
+        cell_quartic_interp_x(i,j,k,n,finearr,srcarr);
+    });
+}
+
 }