Skip to content

Commit

Permalink
Quartic interpolation for cell centered data (AMReX-Codes#2960)
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 23, 2022
1 parent c4b7982 commit 27ef106
Show file tree
Hide file tree
Showing 3 changed files with 210 additions and 0 deletions.
48 changes: 48 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,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
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
93 changes: 93 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,94 @@ 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());
amrex::ignore_unused(run_on_gpu);

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;
if (run_on_gpu) 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;
if (run_on_gpu) 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 27ef106

Please sign in to comment.