InterpFaceRegister

Src/AmrCore/AMReX_AmrCoreFwd.H

@@ -9,10 +9,13 @@ class AmrMesh;
 
 class FluxRegister;
 class Interpolater;
+class MFInterpolater;
 
 class TagBox;
 class TagBoxArray;
 
+class InterpFaceRegister;
+
 }
 
 #endif

Src/AmrCore/AMReX_InterpFaceReg_1D_C.H

@@ -0,0 +1,8 @@
+#ifndef AMREX_INTERP_FACE_REG_1D_C_H_
+#define AMREX_INTERP_FACE_REG_1D_C_H_
+
+namespace amrex {
+
+}
+
+#endif

Src/AmrCore/AMReX_InterpFaceReg_2D_C.H

@@ -0,0 +1,66 @@
+#ifndef AMREX_INTERP_FACE_REG_2D_C_H_
+#define AMREX_INTERP_FACE_REG_2D_C_H_
+
+namespace amrex {
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void interp_face_reg (int i, int j, IntVect const& rr, Array4<Real> const& fine, int scomp,
+                      Array4<Real const> const& crse, Array4<Real> const& slope, int ncomp,
+                      Box const& domface, int idim)
+{
+    int ic = amrex::coarsen(i,rr[0]);
+    int jc = amrex::coarsen(j,rr[1]);
+    if (idim == 0) {
+        if (jc == domface.smallEnd(1) || jc == domface.bigEnd(1)) {
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,0,n+scomp) = crse(ic,jc,0,n);
+            }
+        } else {
+            Real sfy = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic,jc+1,0,n) - crse(ic,jc-1,0,n));
+                Real df = Real(2.0) * (crse(ic,jc+1,0,n) - crse(ic,jc  ,0,n));
+                Real db = Real(2.0) * (crse(ic,jc  ,0,n) - crse(ic,jc-1,0,n));
+                Real sy = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sy = amrex::Math::copysign(Real(1.),dc)*amrex::min(sy,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfy = amrex::min(sfy, sy / dc);
+                }
+                slope(i,j,0,n) = dc;
+            }
+            Real yoff = (j - jc*rr[1] + Real(0.5)) / Real(rr[1]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,0,n+scomp) = crse(ic,jc,0,n) + yoff * slope(i,j,0,n) * sfy;
+            }
+        }
+    } else {
+        if (ic == domface.smallEnd(0) || ic == domface.bigEnd(0)) {
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,0,n+scomp) = crse(ic,jc,0,n);
+            }
+        } else {
+            Real sfx = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic+1,jc,0,n) - crse(ic-1,jc,0,n));
+                Real df = Real(2.0) * (crse(ic+1,jc,0,n) - crse(ic  ,jc,0,n));
+                Real db = Real(2.0) * (crse(ic  ,jc,0,n) - crse(ic-1,jc,0,n));
+                Real sx = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sx = amrex::Math::copysign(Real(1.),dc)*amrex::min(sx,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfx = amrex::min(sfx, sx / dc);
+                }
+                slope(i,j,0,n) = dc;
+            }
+            Real xoff = (i - ic*rr[0] + Real(0.5)) / Real(rr[0]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,0,n+scomp) = crse(ic,jc,0,n) + xoff * slope(i,j,0,n) * sfx;
+            }
+        }
+    }
+}
+
+}
+
+#endif

Src/AmrCore/AMReX_InterpFaceReg_3D_C.H

@@ -0,0 +1,151 @@
+#ifndef AMREX_INTERP_FACE_REG_3D_C_H_
+#define AMREX_INTERP_FACE_REG_3D_C_H_
+
+namespace amrex {
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void interp_face_reg (int i, int j, int k, IntVect const& rr, Array4<Real> const& fine, int scomp,
+                      Array4<Real const> const& crse, Array4<Real> const& slope, int ncomp,
+                      Box const& domface, int idim)
+{
+    int ic = amrex::coarsen(i,rr[0]);
+    int jc = amrex::coarsen(j,rr[1]);
+    int kc = amrex::coarsen(k,rr[2]);
+    if (idim == 0) {
+        if (jc == domface.smallEnd(1) || jc == domface.bigEnd(1)) {
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) = crse(ic,jc,kc,n);
+            }
+        } else {
+            Real sfy = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n));
+                Real df = Real(2.0) * (crse(ic,jc+1,kc,n) - crse(ic,jc  ,kc,n));
+                Real db = Real(2.0) * (crse(ic,jc  ,kc,n) - crse(ic,jc-1,kc,n));
+                Real sy = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sy = amrex::Math::copysign(Real(1.),dc)*amrex::min(sy,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfy = amrex::min(sfy, sy / dc);
+                }
+                slope(i,j,k,n) = dc;
+            }
+            Real yoff = (j - jc*rr[1] + Real(0.5)) / Real(rr[1]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) = crse(ic,jc,kc,n) + yoff * slope(i,j,k,n) * sfy;
+            }
+        }
+
+        if (kc != domface.smallEnd(2) && kc != domface.bigEnd(2)) {
+            Real sfz = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n));
+                Real df = Real(2.0) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc  ,n));
+                Real db = Real(2.0) * (crse(ic,jc,kc  ,n) - crse(ic,jc,kc-1,n));
+                Real sz = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sz = amrex::Math::copysign(Real(1.),dc)*amrex::min(sz,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfz = amrex::min(sfz, sz / dc);
+                }
+                slope(i,j,k,n) = dc;
+            }
+            Real zoff = (k - kc*rr[2] + Real(0.5)) / Real(rr[2]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) += zoff * slope(i,j,k,n) * sfz;
+            }
+        }
+    } else if (idim == 1) {
+        if (ic == domface.smallEnd(0) || ic == domface.bigEnd(0)) {
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) = crse(ic,jc,kc,n);
+            }
+        } else {
+            Real sfx = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n));
+                Real df = Real(2.0) * (crse(ic+1,jc,kc,n) - crse(ic  ,jc,kc,n));
+                Real db = Real(2.0) * (crse(ic  ,jc,kc,n) - crse(ic-1,jc,kc,n));
+                Real sx = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sx = amrex::Math::copysign(Real(1.),dc)*amrex::min(sx,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfx = amrex::min(sfx, sx / dc);
+                }
+                slope(i,j,k,n) = dc;
+            }
+            Real xoff = (i - ic*rr[0] + Real(0.5)) / Real(rr[0]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) = crse(ic,jc,kc,n) + xoff * slope(i,j,k,n) * sfx;
+            }
+        }
+
+        if (kc != domface.smallEnd(2) && kc != domface.bigEnd(2)) {
+            Real sfz = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n));
+                Real df = Real(2.0) * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc  ,n));
+                Real db = Real(2.0) * (crse(ic,jc,kc  ,n) - crse(ic,jc,kc-1,n));
+                Real sz = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sz = amrex::Math::copysign(Real(1.),dc)*amrex::min(sz,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfz = amrex::min(sfz, sz / dc);
+                }
+                slope(i,j,k,n) = dc;
+            }
+            Real zoff = (k - kc*rr[2] + Real(0.5)) / Real(rr[2]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) += zoff * slope(i,j,k,n) * sfz;
+            }
+        }
+    } else {
+        if (ic == domface.smallEnd(0) || ic == domface.bigEnd(0)) {
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) = crse(ic,jc,kc,n);
+            }
+        } else {
+            Real sfx = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n));
+                Real df = Real(2.0) * (crse(ic+1,jc,kc,n) - crse(ic  ,jc,kc,n));
+                Real db = Real(2.0) * (crse(ic  ,jc,kc,n) - crse(ic-1,jc,kc,n));
+                Real sx = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sx = amrex::Math::copysign(Real(1.),dc)*amrex::min(sx,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfx = amrex::min(sfx, sx / dc);
+                }
+                slope(i,j,k,n) = dc;
+            }
+            Real xoff = (i - ic*rr[0] + Real(0.5)) / Real(rr[0]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) = crse(ic,jc,kc,n) + xoff * slope(i,j,k,n) * sfx;
+            }
+        }
+
+        if (jc != domface.smallEnd(1) && jc != domface.bigEnd(1)) {
+            Real sfy = Real(1.0);
+            for (int n = 0; n < ncomp; ++n) {
+                Real dc = Real(0.5) * (crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n));
+                Real df = Real(2.0) * (crse(ic,jc+1,kc,n) - crse(ic,jc  ,kc,n));
+                Real db = Real(2.0) * (crse(ic,jc  ,kc,n) - crse(ic,jc-1,kc,n));
+                Real sy = (df*db >= Real(0.0)) ?
+                    amrex::min(amrex::Math::abs(df),amrex::Math::abs(db)) : Real(0.);
+                sy = amrex::Math::copysign(Real(1.),dc)*amrex::min(sy,amrex::Math::abs(dc));
+                if (dc != Real(0.0)) {
+                    sfy = amrex::min(sfy, sy / dc);
+                }
+                slope(i,j,k,n) = dc;
+            }
+            Real yoff = (j - jc*rr[1] + Real(0.5)) / Real(rr[1]) - Real(0.5);
+            for (int n = 0; n < ncomp; ++n) {
+                fine(i,j,k,n+scomp) += yoff * slope(i,j,k,n) * sfy;
+            }
+        }
+    }
+}
+
+}
+
+#endif

Src/AmrCore/AMReX_InterpFaceReg_C.H

@@ -0,0 +1,12 @@
+#ifndef AMREX_INTERP_FACE_REG_C_H_
+#define AMREX_INTERP_FACE_REG_C_H_
+
+#include <AMReX_Array4.H>
+
+#if (AMREX_SPACEDIM == 2)
+#include <AMReX_InterpFaceReg_2D_C.H>
+#elif (AMREX_SPACEDIM == 3)
+#include <AMReX_InterpFaceReg_3D_C.H>
+#endif
+
+#endif

Src/AmrCore/AMReX_InterpFaceRegister.H

@@ -0,0 +1,85 @@
+#ifndef AMREX_INTERP_FACE_REGISTER_H_
+#define AMREX_INTERP_FACE_REGISTER_H_
+#include <AMReX_Config.H>
+
+#include <AMReX_MultiFab.H>
+#include <AMReX_iMultiFab.H>
+#include <AMReX_Geometry.H>
+
+namespace amrex {
+
+/**
+ *  \brief InterpFaceRegister is a coarse/fine boundary register for
+ *  interpolation of face data at the coarse/fine boundadry.
+ */
+
+class InterpFaceRegister
+{
+public:
+
+    InterpFaceRegister () = default;
+
+    /**
+     * \brief InterpFaceRegister Constructor.
+     *
+     * \param fba  The fine level BoxArray
+     * \param fdm  The fine level DistributionMapping
+     * \param fgeom  The fine level Geometry
+     * \param ref_ratio  The refinement ratio
+     */
+    InterpFaceRegister (BoxArray const& fba, DistributionMapping const& fdm,
+                        Geometry const& fgeom, IntVect const& ref_ratio);
+
+    /**
+     * \brief Defines an InterpFaceRegister objecct.
+     *
+     * \param fba  The fine level BoxArray
+     * \param fdm  The fine level DistributionMapping
+     * \param fgeom  The fine level Geometry
+     * \param ref_ratio  The refinement ratio
+     */
+    void define (BoxArray const& fba, DistributionMapping const& fdm,
+                 Geometry const& fgeom, IntVect const& ref_ratio);
+
+    /**
+     * \brief Returns a coarse/fine boundary mask for a given face.
+     *
+     * This returns an iMultiFab at the coarse/fine boundary of a given
+     * face.  The data are only defined on one face of each box in the fine
+     * level BoxArray passed to the InterpFaceRegister constructor.  It has
+     * two possible values: 1 for coarse/fine boundary and 0 for fine/fine
+     * boundary including non-periodic physical domain face.
+     *
+     * \param face  The face
+     */
+    iMultiFab const& mask (Orientation face) const;
+
+    /**
+     * \brief Interpolates from coarse to fine data at coarse/fine
+     * boundaries.
+     *
+     * \param fine  Array of pointers to the fine data.
+     * \param crse  Array of const pointers to the coarse data.
+     * \param scomp  Starting component
+     * \param ncomp  Number of components
+     */
+    void interp (Array<MultiFab*, AMREX_SPACEDIM> const& fine,
+                 Array<MultiFab const*, AMREX_SPACEDIM> const& crse,
+                 int scomp, int ncomp);
+
+private:
+    BoxArray m_fine_ba;
+    DistributionMapping m_fine_dm;
+    Geometry m_fine_geom;
+    IntVect m_ref_ratio;
+
+    Geometry m_crse_geom;
+
+    Array<BoxArray, 2*AMREX_SPACEDIM> m_fine_face_ba;
+    Array<BoxArray, 2*AMREX_SPACEDIM> m_crse_face_ba;
+    Array<iMultiFab,2*AMREX_SPACEDIM> m_face_mask; // crse/fine: 1, fine/fine & fine/physbc: 0
+};
+
+}
+
+#endif