diff options
Diffstat (limited to 'Source/Utils')
| -rw-r--r-- | Source/Utils/Average.H | 176 | ||||
| -rw-r--r-- | Source/Utils/Average.cpp | 109 |
2 files changed, 202 insertions, 83 deletions
diff --git a/Source/Utils/Average.H b/Source/Utils/Average.H index bc7d81e66..fb1460413 100644 --- a/Source/Utils/Average.H +++ b/Source/Utils/Average.H @@ -8,87 +8,121 @@ namespace Average{ using namespace amrex; /** - * \brief Returns the cell-centered average of the floating point data contained - * in the input - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * \c mf_in_arr. + * \brief Interpolates the floating point data contained in the source Array4 + * \c mf_fp_arr, extracted from a fine MultiFab, by averaging over either + * 1 point or 2 equally distant points. * - * \param[in] mf_in_arr floating point data to be averaged - * \param[in] stag staggering (index type) of the data - * \param[in] i index along x to access the - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * \c mf_in_arr containing the data to be averaged - * \param[in] j index along y to access the - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * \c mf_in_arr containing the data to be averaged - * \param[in] k index along z to access the - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * \c mf_in_arr containing the data to be averaged - * \param[in] comp index along the fourth component of the - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * \c mf_in_arr containing the data to be averaged - * \return averaged field at cell (i,j,k) + * \param[in] arr_src floating point data to be interpolated + * \param[in] sf staggering of the source fine MultiFab + * \param[in] sc staggering of the destination coarsened MultiFab + * \param[in] cr coarsening ratio along each spatial direction + * \param[in] np number of points to loop over for interpolation + * \param[in] i index along x of the coarsened Array4 to be filled + * \param[in] j index along y of the coarsened Array4 to be filled + * \param[in] k index along z of the coarsened Array4 to be filled + * \param[in] comp index along the fourth component of the Array4 \c arr_src + * containing the data to be interpolated + * + * \return interpolated field at cell (i,j,k) of a coarsened Array4 */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - Real ToCellCenter ( Array4<Real const> const& mf_in_arr, - const IntVect stag, - const int i, - const int j, - const int k, - const int comp=0 ) + Real CoarsenAndInterpolateKernel ( Array4<Real const> const& arr_src, + int const* const AMREX_RESTRICT sf, + int const* const AMREX_RESTRICT sc, + int const* const AMREX_RESTRICT cr, + int const* const AMREX_RESTRICT np, + const int i, + const int j, + const int k, + const int comp ) { - const int sx = stag[0]; - const int sy = stag[1]; -#if (AMREX_SPACEDIM == 2) - constexpr int sz = 0; -#elif (AMREX_SPACEDIM == 3) - const int sz = stag[2]; -#endif - return 0.125_rt * ( mf_in_arr(i ,j ,k ,comp) - + mf_in_arr(i+sx,j ,k ,comp) - + mf_in_arr(i ,j+sy,k ,comp) - + mf_in_arr(i ,j ,k+sz,comp) - + mf_in_arr(i+sx,j+sy,k ,comp) - + mf_in_arr(i ,j+sy,k+sz,comp) - + mf_in_arr(i+sx,j ,k+sz,comp) - + mf_in_arr(i+sx,j+sy,k+sz,comp) ); + // Indices of source coarse array + const int ic[3] = { i, j, k }; + + // Index of first point of source Array4 from which interpolation is done + int idx_min[3]; + + // Compute number of points to loop over (either 1 or 2) + // and starting indices of fine array in each direction + for ( int l = 0; l < 3; ++l ) { + if ( cr[l] == 1 ) idx_min[l] = ic[l]-sc[l]*(1-sf[l]); // no coarsening + else idx_min[l] = ic[l]*cr[l]+static_cast<int>(cr[l]/2)*(1-sc[l])-(1-sf[l]); + } + + // Auxiliary integer variables + const int numx = np[0]; + const int numy = np[1]; + const int numz = np[2]; + const int imin = idx_min[0]; + const int jmin = idx_min[1]; + const int kmin = idx_min[2]; + + // Interpolate over points computed above + Real c = 0.0_rt; + for (int kref = 0; kref < numz; ++kref) { + for (int jref = 0; jref < numy; ++jref) { + for (int iref = 0; iref < numx; ++iref) { + c += arr_src(imin+iref,jmin+jref,kmin+kref,comp); + } + } + } + return c / (numx*numy*numz); }; /** - * \brief Stores the cell-centered average of the floating point data contained - * in the input - * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a> - * \c mf_in into the output - * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a> - * \c mf_out. + * \brief Loops over the boxes of the coarsened MultiFab \c mf_dst and fills + * them by interpolating the data contained in the fine MultiFab \c mf_src. + * + * \param[in,out] mf_dst coarsened MultiFab containing the floating point data + * to be filled by interpolating the source fine MultiFab + * \param[in] mf_src fine MultiFab containing the floating point data to be interpolated + * \param[in] dcomp offset for the fourth component of the coarsened Array4 + * object, extracted from its MultiFab, where the interpolated + * values will be stored + * \param[in] scomp offset for the fourth component of the fine Array4 + * object, extracted from its MultiFab, containing the + * data to be interpolated + * \param[in] ncomp number of components to loop over for the coarsened + * Array4 extracted from the coarsened MultiFab \c mf_dst + * \param[in] ngrow number of guard cells to fill + * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src + * and the coarsened MultiFab \c mf_dst + */ + void CoarsenAndInterpolateLoop ( MultiFab& mf_dst, + const MultiFab& mf_src, + const int dcomp, + const int scomp, + const int ncomp, + const int ngrow, + const IntVect crse_ratio=IntVect(1) ); + + /** + * \brief Stores in the coarsened MultiFab \c mf_dst the values obtained by + * interpolating the data contained in the fine MultiFab \c mf_src. * - * \param[in,out] mf_out <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a> - * containing the floating point data to be filled with cell-centered averages - * \param[in] mf_in <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a> - * containing the floating point data to be averaged - * \param[in] dcomp offset for the fourth component of the - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * object, extracted from its - * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a> - * , where the cell-centered averages will be stored - * \param[in] scomp optional offset for the fourth component of the - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * object, extracted from its - * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a> - * , containing the data to be averaged - * \param[in] ncomp optional number of components to loop over for both the - * input and the output - * <a href="https://amrex-codes.github.io/amrex/doxygen/structamrex_1_1Array4.html">Array4</a> - * objects, extracted from their respective - * <a href="https://amrex-codes.github.io/amrex/doxygen/classamrex_1_1MultiFab.html">MultiFab</a> + * \param[in,out] mf_dst coarsened MultiFab containing the floating point data + * to be filled by interpolating the fine MultiFab \c mf_src + * \param[in] mf_src fine MultiFab containing the floating point data to be interpolated + * \param[in] dcomp offset for the fourth component of the coarsened Array4 + * object, extracted from its MultiFab, where the interpolated + * values will be stored + * \param[in] scomp offset for the fourth component of the fine Array4 + * object, extracted from its MultiFab, containing the + * data to be interpolated + * \param[in] ncomp number of components to loop over for the coarsened + * Array4 extracted from the coarsened MultiFab \c mf_dst + * \param[in] ngrow number of guard cells to fill + * \param[in] crse_ratio coarsening ratio between the fine MultiFab \c mf_src + * and the coarsened MultiFab \c mf_dst */ - void ToCellCenter ( MultiFab& mf_out, - const MultiFab& mf_in, - const int dcomp, - const int ngrow, - const int scomp=0, - const int ncomp=1 ); + void CoarsenAndInterpolate ( MultiFab& mf_dst, + const MultiFab& mf_src, + const int dcomp, + const int scomp, + const int ncomp, + const int ngrow, + const IntVect crse_ratio=IntVect(1) ); } #endif // WARPX_AVERAGE_H_ diff --git a/Source/Utils/Average.cpp b/Source/Utils/Average.cpp index b98f1e948..ba03bf5b8 100644 --- a/Source/Utils/Average.cpp +++ b/Source/Utils/Average.cpp @@ -3,27 +3,112 @@ using namespace amrex; void -Average::ToCellCenter ( MultiFab& mf_out, - const MultiFab& mf_in, - const int dcomp, - const int ngrow, - const int scomp, - const int ncomp ) +Average::CoarsenAndInterpolateLoop ( MultiFab& mf_dst, + const MultiFab& mf_src, + const int dcomp, + const int scomp, + const int ncomp, + const int ngrow, + const IntVect crse_ratio ) { + // Staggering of source fine MultiFab and destination coarse MultiFab + const IntVect stag_src = mf_src.boxArray().ixType().toIntVect(); + const IntVect stag_dst = mf_dst.boxArray().ixType().toIntVect(); + + if ( crse_ratio > IntVect(1) ) AMREX_ALWAYS_ASSERT_WITH_MESSAGE( ngrow == 0, + "option of filling guard cells of destination MultiFab with coarsening not supported for this interpolation" ); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( mf_src.nGrowVect() >= stag_dst-stag_src+IntVect(ngrow), + "source fine MultiFab does not have enough guard cells for this interpolation" ); + + // Auxiliary integer arrays (always 3D) + Gpu::ManagedVector<int> sf_gpuarr, sc_gpuarr, cr_gpuarr, np_gpuarr; + sf_gpuarr.resize( 3 ); // staggering of source fine MultiFab + sc_gpuarr.resize( 3 ); // staggering of destination coarse MultiFab + cr_gpuarr.resize( 3 ); // coarsening ratio + np_gpuarr.resize( 3 ); // number of points to loop over for interpolation + + sf_gpuarr[0] = stag_src[0]; + sf_gpuarr[1] = stag_src[1]; +#if (AMREX_SPACEDIM == 2) + sf_gpuarr[2] = 0; +#elif (AMREX_SPACEDIM == 3) + sf_gpuarr[2] = stag_src[2]; +#endif + + sc_gpuarr[0] = stag_dst[0]; + sc_gpuarr[1] = stag_dst[1]; +#if (AMREX_SPACEDIM == 2) + sc_gpuarr[2] = 0; +#elif (AMREX_SPACEDIM == 3) + sc_gpuarr[2] = stag_dst[2]; +#endif + + cr_gpuarr[0] = crse_ratio[0]; + cr_gpuarr[1] = crse_ratio[1]; +#if (AMREX_SPACEDIM == 2) + cr_gpuarr[2] = 1; +#elif (AMREX_SPACEDIM == 3) + cr_gpuarr[2] = crse_ratio[2]; +#endif + + int const* const AMREX_RESTRICT sf = sf_gpuarr.data(); + int const* const AMREX_RESTRICT sc = sc_gpuarr.data(); + int const* const AMREX_RESTRICT cr = cr_gpuarr.data(); + int * const AMREX_RESTRICT np = np_gpuarr.data(); + + // Compute number of points to loop over (either 1 or 2) in each direction + for ( int l = 0; l < 3; ++l ) { + if ( cr[l] == 1 ) np[l] = 1+abs(sf[l]-sc[l]); // no coarsening + else np[l] = 2-sf[l]; + } + #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif // Loop over boxes (or tiles if not on GPU) - for (MFIter mfi( mf_out, TilingIfNotGPU() ); mfi.isValid(); ++mfi) + for (MFIter mfi( mf_dst, TilingIfNotGPU() ); mfi.isValid(); ++mfi) { - const Box bx = mfi.growntilebox( ngrow ); - Array4<Real> const& mf_out_arr = mf_out.array( mfi ); - Array4<Real const> const& mf_in_arr = mf_in.const_array( mfi ); - const IntVect stag = mf_in.boxArray().ixType().toIntVect(); + // Tiles defined at the coarse level + const Box& bx = mfi.growntilebox( ngrow ); + Array4<Real> const& arr_dst = mf_dst.array( mfi ); + Array4<Real const> const& arr_src = mf_src.const_array( mfi ); ParallelFor( bx, ncomp, [=] AMREX_GPU_DEVICE( int i, int j, int k, int n ) { - mf_out_arr(i,j,k,n+dcomp) = Average::ToCellCenter( mf_in_arr, stag, i, j, k, n+scomp ); + arr_dst(i,j,k,n+dcomp) = Average::CoarsenAndInterpolateKernel( + arr_src, sf, sc, cr, np, i, j, k, n+scomp ); } ); } } + +void +Average::CoarsenAndInterpolate ( MultiFab& mf_dst, + const MultiFab& mf_src, + const int dcomp, + const int scomp, + const int ncomp, + const int ngrow, + const IntVect crse_ratio ) +{ + BL_PROFILE( "Average::CoarsenAndInterpolate" ); + + // Convert BoxArray of source MultiFab to staggering of destination MultiFab and coarsen it (if possible) + BoxArray ba_tmp = amrex::convert( mf_src.boxArray(), mf_dst.ixType().toIntVect() ); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( ba_tmp.coarsenable( crse_ratio ), + "source MultiFab converted to staggering of destination MultiFab is not coarsenable" ); + ba_tmp.coarsen( crse_ratio ); + + if ( ba_tmp == mf_dst.boxArray() and mf_src.DistributionMap() == mf_dst.DistributionMap() ) + Average::CoarsenAndInterpolateLoop( mf_dst, mf_src, dcomp, scomp, ncomp, ngrow, crse_ratio ); + else + { + // Cannot coarsen into MultiFab with different BoxArray or DistributionMapping: + // 1) create temporary MultiFab on coarsened version of source BoxArray with same DistributionMapping + MultiFab mf_tmp( ba_tmp, mf_src.DistributionMap(), ncomp, 0, MFInfo(), FArrayBoxFactory() ); + // 2) interpolate from mf_src to mf_tmp (start writing into component 0) + Average::CoarsenAndInterpolateLoop( mf_tmp, mf_src, 0, scomp, ncomp, ngrow, crse_ratio ); + // 3) copy from mf_tmp to mf_dst (with different BoxArray or DistributionMapping) + mf_dst.copy( mf_tmp, 0, dcomp, ncomp ); + } +} |