Skip to content

Commit

Permalink
Quartic interpolation for cell centered data
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Sep 22, 2022
1 parent 3e5cc77 commit 43d415f
Show file tree
Hide file tree
Showing 3 changed files with 210 additions and 0 deletions.
51 changes: 51 additions & 0 deletions Src/AmrCore/AMReX_Interp_C.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,5 +135,56 @@ 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 = (ii * 2 == i) ? -1 : 1;
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 = (jj * 2 == j) ? -1 : 1;
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 = (kk * 2 == k) ? -1 : 1;
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
69 changes: 69 additions & 0 deletions Src/AmrCore/AMReX_Interpolater.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;

}

Expand Down
90 changes: 90 additions & 0 deletions Src/AmrCore/AMReX_Interpolater.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -37,6 +39,7 @@ CellConservativeProtected protected_interp;
CellConservativeQuartic quartic_interp;
CellBilinear cell_bilinear_interp;
CellQuadratic quadratic_interp;
CellQuartic cell_quartic_interp;

NodeBilinear::~NodeBilinear () {}

Expand Down Expand Up @@ -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);
});
}

}

0 comments on commit 43d415f

Please sign in to comment.