From 0d45d8c407236d5e542b6a57e348def385f488a8 Mon Sep 17 00:00:00 2001 From: WeiqunZhang Date: Thu, 9 Jul 2020 09:04:56 -0700 Subject: [PATCH] FaceLinear Interpolater for face data (#1118) It's linear in the nodal direction and constant in others. This preserves the divergence, but it's only first-order. --- Src/AmrCore/AMReX_Interp_1D_C.H | 10 ++++ Src/AmrCore/AMReX_Interp_2D_C.H | 22 +++++++++ Src/AmrCore/AMReX_Interp_3D_C.H | 36 +++++++++++++++ Src/AmrCore/AMReX_Interpolater.H | 66 ++++++++++++++++++++++++++ Src/AmrCore/AMReX_Interpolater.cpp | 74 +++++++++++++++++++++++++++++- 5 files changed, 207 insertions(+), 1 deletion(-) diff --git a/Src/AmrCore/AMReX_Interp_1D_C.H b/Src/AmrCore/AMReX_Interp_1D_C.H index 21420e79715..146f833a0d3 100644 --- a/Src/AmrCore/AMReX_Interp_1D_C.H +++ b/Src/AmrCore/AMReX_Interp_1D_C.H @@ -296,6 +296,16 @@ nodebilin_interp (Box const& bx, Array4 const& fine, const int fcomp, const i } } +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_x (int i, int /*j*/, int /*k*/, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + int ii = amrex::coarsen(i,ratio[0]); + Real const w = static_cast(i-ii*ratio[0]) * (1.0_rt/ratio[0]); + fine(i,0,0,n) = (1.0_rt-w) * crse(ii,0,0,n) + w * crse(ii+1,0,0,n); +} + } #endif diff --git a/Src/AmrCore/AMReX_Interp_2D_C.H b/Src/AmrCore/AMReX_Interp_2D_C.H index 7d8396aa01e..16a0a1f6319 100644 --- a/Src/AmrCore/AMReX_Interp_2D_C.H +++ b/Src/AmrCore/AMReX_Interp_2D_C.H @@ -428,6 +428,28 @@ nodebilin_interp (Box const& bx, Array4 const& fine, const int fcomp, const i } } +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_x (int i, int j, int /*k*/, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + int ii = amrex::coarsen(i,ratio[0]); + int jj = amrex::coarsen(j,ratio[1]); + Real const w = static_cast(i-ii*ratio[0]) * (1.0_rt/ratio[0]); + fine(i,j,0,n) = (1.0_rt-w) * crse(ii,jj,0,n) + w * crse(ii+1,jj,0,n); +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_y (int i, int j, int /*k*/, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + int ii = amrex::coarsen(i,ratio[0]); + int jj = amrex::coarsen(j,ratio[1]); + Real const w = static_cast(j-jj*ratio[1]) * (1.0_rt/ratio[1]); + fine(i,j,0,n) = (1.0_rt-w) * crse(ii,jj,0,n) + w * crse(ii,jj+1,0,n); +} + } #endif diff --git a/Src/AmrCore/AMReX_Interp_3D_C.H b/Src/AmrCore/AMReX_Interp_3D_C.H index c3450f2bd0d..8badd0d1101 100644 --- a/Src/AmrCore/AMReX_Interp_3D_C.H +++ b/Src/AmrCore/AMReX_Interp_3D_C.H @@ -566,6 +566,42 @@ nodebilin_interp (Box const& bx, Array4 const& fine, const int fcomp, const i } } +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_x (int i, int j, int k, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + int ii = amrex::coarsen(i,ratio[0]); + int jj = amrex::coarsen(j,ratio[1]); + int kk = amrex::coarsen(k,ratio[2]); + Real const w = static_cast(i-ii*ratio[0]) * (1.0_rt/ratio[0]); + fine(i,j,k,n) = (1.0_rt-w) * crse(ii,jj,kk,n) + w * crse(ii+1,jj,kk,n); +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_y (int i, int j, int k, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + int ii = amrex::coarsen(i,ratio[0]); + int jj = amrex::coarsen(j,ratio[1]); + int kk = amrex::coarsen(k,ratio[2]); + Real const w = static_cast(j-jj*ratio[1]) * (1.0_rt/ratio[1]); + fine(i,j,k,n) = (1.0_rt-w) * crse(ii,jj,kk,n) + w * crse(ii,jj+1,kk,n); +} + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +face_linear_interp_z (int i, int j, int k, int n, Array4 const& fine, + Array4 const& crse, IntVect const& ratio) noexcept +{ + int ii = amrex::coarsen(i,ratio[0]); + int jj = amrex::coarsen(j,ratio[1]); + int kk = amrex::coarsen(k,ratio[2]); + Real const w = static_cast(k-kk*ratio[2]) * (1.0_rt/ratio[2]); + fine(i,j,k,n) = (1.0_rt-w) * crse(ii,jj,kk,n) + w * crse(ii,jj,kk+1,n); +} + } #endif diff --git a/Src/AmrCore/AMReX_Interpolater.H b/Src/AmrCore/AMReX_Interpolater.H index bd2c7f58f17..0421669a5f5 100644 --- a/Src/AmrCore/AMReX_Interpolater.H +++ b/Src/AmrCore/AMReX_Interpolater.H @@ -646,10 +646,76 @@ public: }; #endif +/** +* \brief Bilinear interpolation on face data. +* +* Bilinear interpolation on data. +*/ + +class FaceLinear + : + public Interpolater +{ +public: + + /** + * \brief The destructor. + */ + virtual ~FaceLinear () 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 const& bcr, + int actual_comp, + int actual_state, + RunOn gpu_or_cpu) override; +}; //! CONSTRUCT A GLOBAL OBJECT OF EACH VERSION. extern PCInterp pc_interp; extern NodeBilinear node_bilinear_interp; +extern FaceLinear face_linear_interp; extern CellConservativeLinear lincc_interp; extern CellConservativeLinear cell_cons_interp; diff --git a/Src/AmrCore/AMReX_Interpolater.cpp b/Src/AmrCore/AMReX_Interpolater.cpp index 4da717bdc0a..6b57b4b1683 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -13,7 +13,8 @@ namespace amrex { // -// PCInterp, NodeBilinear, and CellConservativeLinear are supported for all dimensions on cpu and gpu. +// PCInterp, NodeBilinear, FaceLinear, and CellConservativeLinear are supported for all dimensions +// on cpu and gpu. // // CellConsertiveProtected only works in 2D and 3D on cpu. // @@ -29,6 +30,7 @@ namespace amrex { // PCInterp pc_interp; NodeBilinear node_bilinear_interp; +FaceLinear face_linear_interp; CellConservativeLinear lincc_interp; CellConservativeLinear cell_cons_interp(0); @@ -141,6 +143,76 @@ NodeBilinear::interp (const FArrayBox& crse, }); } +Box +FaceLinear::CoarseBox (const Box& fine, int ratio) +{ + return CoarseBox(fine, IntVect(ratio)); +} + +Box +FaceLinear::CoarseBox (const Box& fine, const IntVect& ratio) +{ + Box b = amrex::coarsen(fine,ratio); + for (int i = 0; i < AMREX_SPACEDIM; i++) { + if (b.type(i) == IndexType::NODE && b.length(i) < 2) { + // Don't want degenerate boxes in nodal direction. + b.growHi(i,1); + } + } + return b; +} + +void +FaceLinear::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 const& /*bcr*/, + int /*actual_comp*/, + int /*actual_state*/, + RunOn runon) +{ + BL_PROFILE("FaceLinear::interp()"); + + AMREX_ASSERT(AMREX_D_TERM(fine_region.type(0),+fine_region.type(1),+fine_region.type(2)) == 1); + + Array4 const& fine_arr = fine.array(fine_comp); + Array4 const& crse_arr = crse.const_array(crse_comp); + + if (fine_region.type(0) == IndexType::NODE) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_linear_interp_x(i,j,k,n,fine_arr,crse_arr,ratio); + }); + } +#if (AMREX_SPACEDIM >= 2) + else if (fine_region.type(1) == IndexType::NODE) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_linear_interp_y(i,j,k,n,fine_arr,crse_arr,ratio); + }); + } +#if (AMREX_SPACEDIM == 3) + else + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon,fine_region,ncomp,i,j,k,n, + { + face_linear_interp_z(i,j,k,n,fine_arr,crse_arr,ratio); + }); + } +#endif +#endif +} + +FaceLinear::~FaceLinear () {} + #ifndef BL_NO_FORT CellBilinear::~CellBilinear () {}