From 15658267709e5a53a11ef6f945d1bfb283f994a6 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 24 Sep 2019 16:25:49 -0700 Subject: Current Synchronize: Port to GPU Port the current synchronize functions to GPU. --- Source/Parallelization/WarpXComm.cpp | 44 ++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 24 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 990d0f988..22fa4820e 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include @@ -348,36 +349,34 @@ WarpX::SyncCurrent () } } -/** \brief Fills the values of the current on the coarse patch by - * averaging the values of the current of the fine patch (on the same level). - */ void WarpX::SyncCurrent (const std::array& fine, - const std::array< amrex::MultiFab*,3>& crse, - int refinement_ratio) + const std::array< amrex::MultiFab*,3>& coarse, + int const refinement_ratio) { BL_ASSERT(refinement_ratio == 2); - const IntVect& ng = (fine[0]->nGrowVect() + 1) /refinement_ratio; + const IntVect& ng = (fine[0]->nGrowVect() + 1) / refinement_ratio; // add equivalent no. of guards to coarse patch #ifdef _OPEMP -#pragma omp parallel +#pragma omp parallel if (Gpu::notInLaunchRegion()) #endif { - FArrayBox ffab; - for (int idim = 0; idim < 3; ++idim) + FArrayBox ffab; // contiguous, temporary, copy of the tiled fine patch to read from + for (int idim = 0; idim < fine.size(); ++idim) // j-field components { - for (MFIter mfi(*crse[idim],true); mfi.isValid(); ++mfi) + for (MFIter mfi(*coarse[idim],true); mfi.isValid(); ++mfi) // OMP in-box decomposition of coarse into tilebox { - const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,refinement_ratio),1); - ffab.resize(fbx); - fbx &= (*fine[idim])[mfi].box(); - ffab.setVal(0.0); - ffab.copy((*fine[idim])[mfi], fbx, 0, fbx, 0, fine[idim]->nComp()); - WRPX_SYNC_CURRENT(bx.loVect(), bx.hiVect(), - BL_TO_FORTRAN_ANYD((*crse[idim])[mfi]), - BL_TO_FORTRAN_ANYD(ffab), - &idim); + const Box& bx = mfi.growntilebox(ng); // only grow to outer directions of tileboxes for filling guards + + auto const & arrFine = fine[idim]->const_array(mfi); + auto const & arrCoarse = coarse[idim]->array(mfi); + + if( idim == 0 ) + amrex::ParallelFor( bx, WarpxSyncCurrent<0>(arrFine, arrCoarse, refinement_ratio) ); + else if( idim == 1 ) + amrex::ParallelFor( bx, WarpxSyncCurrent<1>(arrFine, arrCoarse, refinement_ratio) ); + else if( idim == 2 ) + amrex::ParallelFor( bx, WarpxSyncCurrent<2>(arrFine, arrCoarse, refinement_ratio) ); } } } @@ -407,9 +406,6 @@ WarpX::SyncRho () } } -/** \brief Fills the values of the charge density on the coarse patch by - * averaging the values of the charge density of the fine patch (on the same level). - */ void WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int refinement_ratio) { @@ -418,7 +414,7 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int refinement_ratio) const int nc = fine.nComp(); #ifdef _OPEMP -#pragma omp parallel +#pragma omp parallel if (Gpu::notInLaunchRegion()) #endif { FArrayBox ffab; -- cgit v1.2.3 From d67030eb709d965672a9900132d0085922335cc1 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 25 Sep 2019 09:54:52 -0700 Subject: SyncCurrent: Disable Tiling for GPU Co-authored-by: Weiqun Zhang --- Source/Parallelization/WarpXComm.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 22fa4820e..40f8203a9 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -361,10 +361,10 @@ WarpX::SyncCurrent (const std::array& fine, #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif { - FArrayBox ffab; // contiguous, temporary, copy of the tiled fine patch to read from for (int idim = 0; idim < fine.size(); ++idim) // j-field components { - for (MFIter mfi(*coarse[idim],true); mfi.isValid(); ++mfi) // OMP in-box decomposition of coarse into tilebox + // OMP in-box decomposition of coarse into tilebox + for (MFIter mfi(*coarse[idim], TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); // only grow to outer directions of tileboxes for filling guards -- cgit v1.2.3 From c9577f8d200d99c40d15d5ff0d2fdacbddc2026f Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 27 Sep 2019 11:14:00 -0700 Subject: Rename, Profile & Remove TODO --- Source/Parallelization/CurrentSynchronize.H | 175 --------------------- .../InterpolateCurrentFineToCoarse.H | 175 +++++++++++++++++++++ Source/Parallelization/WarpXComm.cpp | 19 +-- Source/WarpX.H | 6 +- 4 files changed, 188 insertions(+), 187 deletions(-) delete mode 100644 Source/Parallelization/CurrentSynchronize.H create mode 100644 Source/Parallelization/InterpolateCurrentFineToCoarse.H (limited to 'Source/Parallelization/WarpXComm.cpp') diff --git a/Source/Parallelization/CurrentSynchronize.H b/Source/Parallelization/CurrentSynchronize.H deleted file mode 100644 index 5329ca242..000000000 --- a/Source/Parallelization/CurrentSynchronize.H +++ /dev/null @@ -1,175 +0,0 @@ -/* Copyright 2019 Axel Huebl, Weiqun Zhang - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ - -#ifndef WARPX_CURRENTSYNCHRONIZE_H -#define WARPX_CURRENTSYNCHRONIZE_H - -#include -#include -#include -#include -#include - -#include // std::move - - -/** Fill a current coarse patch with averaged values from a fine patch - * - * Fills the values of the current for a selected component on the coarse patch - * by averaging the values of the current of the fine patch. - * - * @tparam IDim j-field component on which the averaging is performed - */ -template< int IDim > -class WarpxSyncCurrent -{ -public: - /** Construct with fine and coarse patch and their refinement ratio - * - * @param[in] fine read-only fine patch - * @param[in,out] coarse overwritten coarse patch - * @param[in] refinement_ratio ratio between coarse and fine patch granularity - * (currently, only a value of is implemented) - */ - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - WarpxSyncCurrent( - amrex::Array4 const fine, - amrex::Array4 const coarse, - int const refinement_ratio - ) : m_fine(std::move(fine)), - m_coarse(std::move(coarse)), - m_refinement_ratio(std::move(refinement_ratio)) - { - //! @note constants and stencils in operator() implementation assume 2x refinement - BL_ASSERT(refinement_ratio == 2); - } - - AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void - operator()( - int const i, - int const j, - int const k - ) const noexcept // TODO rename to jkl - { - auto const & fine_unsafe = m_fine; // out-of-bounds access not secured with zero-values yet - auto const & coarse = m_coarse; // out-of-bounds access not secured but will also not occur - - // return zero for out-of-bounds accesses during interpolation - // this is efficiently used as a method to add neutral elements beyond guards in the average below - auto const fine = [fine_unsafe] AMREX_GPU_DEVICE (int const j, int const k, int const l) noexcept - { - return fine_unsafe.contains(j, k, l) ? fine_unsafe(j, k, l) : amrex::Real{0.}; - }; - - int const ii = i * m_refinement_ratio; - int const jj = j * m_refinement_ratio; - int const kk = k * m_refinement_ratio; -#if AMREX_SPACEDIM == 2 - if (IDim == 0) { - coarse(i, j, k) = 0.25 * ( - fine(ii, jj, kk) + fine(ii + 1, jj, kk) + - 0.5 * ( - fine(ii, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + - fine(ii, jj + 1, kk) + fine(ii + 1, jj + 1, kk) - ) - ); - } else if (IDim == 2) { - coarse(i, j, k) = 0.25 * ( - fine(ii, jj, kk) + fine(ii, jj + 1, kk) + - 0.5 * ( - fine(ii - 1, jj, kk) + fine(ii - 1, jj + 1, kk) + - fine(ii + 1, jj, kk) + fine(ii + 1, jj + 1, kk) - ) - ); - } else { - coarse(i, j, k) = 0.25 * ( - fine(ii, jj, kk) + - 0.5 * ( - fine(ii - 1, jj , kk) + fine(ii + 1, jj , kk) + - fine(ii , jj - 1, kk) + fine(ii , jj + 1, kk) - ) + - 0.25 * ( - fine(ii - 1, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + - fine(ii - 1, jj + 1, kk) + fine(ii + 1, jj + 1, kk) - ) - ); - } -#elif AMREX_SPACEDIM == 3 - if (IDim == 0) { - coarse(i,j,k) = 0.125 * ( - fine(ii , jj, kk) + - 0.5 * ( - fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) + - fine(ii , jj , kk-1) + fine(ii , jj , kk+1) - ) + - 0.25 * ( - fine(ii , jj-1, kk-1) + fine(ii , jj+1, kk-1) + - fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) - ) + - fine(ii+1, jj, kk) + - 0.5 * ( - fine(ii+1, jj-1, kk ) + fine(ii+1, jj+1, kk ) + - fine(ii+1, jj , kk-1) + fine(ii+1, jj , kk+1) - ) + - 0.25 * ( - fine(ii+1, jj-1, kk-1) + fine(ii+1, jj+1, kk-1) + - fine(ii+1, jj-1, kk+1) + fine(ii+1, jj+1, kk+1) - ) - ); - } else if (IDim == 1) { - coarse(i, j, k) = 0.125 * ( - fine(ii, jj , kk) + - 0.5 * ( - fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + - fine(ii , jj , kk-1) + fine(ii , jj , kk+1) - ) + - 0.25 * ( - fine(ii-1, jj , kk-1) + fine(ii+1, jj , kk-1) + - fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) - ) + - fine(ii, jj+1, kk) + - 0.5 * ( - fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) + - fine(ii , jj+1, kk-1) + fine(ii , jj+1, kk+1) - ) + - 0.25 * ( - fine(ii-1, jj+1, kk-1) + fine(ii+1, jj+1, kk-1) + - fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) - ) - ); - } else { - coarse(i, j, k) = 0.125 * ( - fine(ii, jj, kk ) + - 0.5 * ( - fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + - fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) - ) + - 0.25 * ( - fine(ii-1, jj-1, kk ) + fine(ii+1, jj-1, kk ) + - fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) - ) + - fine(ii, jj, kk+1) + - 0.5 * ( - fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) + - fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) - ) + - 0.25 * ( - fine(ii-1, jj-1, kk+1) + fine(ii+1, jj-1, kk+1) + - fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) - ) - ); - } -#endif - } -private: - amrex::Array4< amrex::Real const > m_fine; - amrex::Array4< amrex::Real > m_coarse; - int m_refinement_ratio; -}; - -#endif //WARPX_CURRENTSYNCHRONIZE_H diff --git a/Source/Parallelization/InterpolateCurrentFineToCoarse.H b/Source/Parallelization/InterpolateCurrentFineToCoarse.H new file mode 100644 index 000000000..148b725d0 --- /dev/null +++ b/Source/Parallelization/InterpolateCurrentFineToCoarse.H @@ -0,0 +1,175 @@ +/* Copyright 2019 Axel Huebl, Weiqun Zhang + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef INTERPOLATECURRENTFINETOCOARSE_H +#define INTERPOLATECURRENTFINETOCOARSE_H + +#include +#include +#include +#include +#include + +#include // std::move + + +/** Fill a current coarse patch with averaged values from a fine patch + * + * Fills the values of the current for a selected component on the coarse patch + * by averaging the values of the current of the fine patch. + * + * @tparam IDim j-field component on which the averaging is performed + */ +template< int IDim > +class InterpolateCurrentFineToCoarse +{ +public: + /** Construct with fine and coarse patch and their refinement ratio + * + * @param[in] fine read-only fine patch + * @param[in,out] coarse overwritten coarse patch + * @param[in] refinement_ratio ratio between coarse and fine patch granularity + * (currently, only a value of is implemented) + */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + InterpolateCurrentFineToCoarse( + amrex::Array4< amrex::Real const > const fine, + amrex::Array4< amrex::Real > const coarse, + int const refinement_ratio + ) : m_fine(std::move(fine)), + m_coarse(std::move(coarse)), + m_refinement_ratio(std::move(refinement_ratio)) + { + //! @note constants and stencils in operator() implementation assume 2x refinement + BL_ASSERT(refinement_ratio == 2); + } + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void + operator()( + int const i, + int const j, + int const k + ) const noexcept + { + auto const & fine_unsafe = m_fine; // out-of-bounds access not secured with zero-values yet + auto const & coarse = m_coarse; // out-of-bounds access not secured but will also not occur + + // return zero for out-of-bounds accesses during interpolation + // this is efficiently used as a method to add neutral elements beyond guards in the average below + auto const fine = [fine_unsafe] AMREX_GPU_DEVICE (int const j, int const k, int const l) noexcept + { + return fine_unsafe.contains(j, k, l) ? fine_unsafe(j, k, l) : amrex::Real{0.}; + }; + + int const ii = i * m_refinement_ratio; + int const jj = j * m_refinement_ratio; + int const kk = k * m_refinement_ratio; +#if AMREX_SPACEDIM == 2 + if (IDim == 0) { + coarse(i, j, k) = 0.25 * ( + fine(ii, jj, kk) + fine(ii + 1, jj, kk) + + 0.5 * ( + fine(ii, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + + fine(ii, jj + 1, kk) + fine(ii + 1, jj + 1, kk) + ) + ); + } else if (IDim == 2) { + coarse(i, j, k) = 0.25 * ( + fine(ii, jj, kk) + fine(ii, jj + 1, kk) + + 0.5 * ( + fine(ii - 1, jj, kk) + fine(ii - 1, jj + 1, kk) + + fine(ii + 1, jj, kk) + fine(ii + 1, jj + 1, kk) + ) + ); + } else { + coarse(i, j, k) = 0.25 * ( + fine(ii, jj, kk) + + 0.5 * ( + fine(ii - 1, jj , kk) + fine(ii + 1, jj , kk) + + fine(ii , jj - 1, kk) + fine(ii , jj + 1, kk) + ) + + 0.25 * ( + fine(ii - 1, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + + fine(ii - 1, jj + 1, kk) + fine(ii + 1, jj + 1, kk) + ) + ); + } +#elif AMREX_SPACEDIM == 3 + if (IDim == 0) { + coarse(i,j,k) = 0.125 * ( + fine(ii , jj, kk) + + 0.5 * ( + fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) + + fine(ii , jj , kk-1) + fine(ii , jj , kk+1) + ) + + 0.25 * ( + fine(ii , jj-1, kk-1) + fine(ii , jj+1, kk-1) + + fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) + ) + + fine(ii+1, jj, kk) + + 0.5 * ( + fine(ii+1, jj-1, kk ) + fine(ii+1, jj+1, kk ) + + fine(ii+1, jj , kk-1) + fine(ii+1, jj , kk+1) + ) + + 0.25 * ( + fine(ii+1, jj-1, kk-1) + fine(ii+1, jj+1, kk-1) + + fine(ii+1, jj-1, kk+1) + fine(ii+1, jj+1, kk+1) + ) + ); + } else if (IDim == 1) { + coarse(i, j, k) = 0.125 * ( + fine(ii, jj , kk) + + 0.5 * ( + fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + + fine(ii , jj , kk-1) + fine(ii , jj , kk+1) + ) + + 0.25 * ( + fine(ii-1, jj , kk-1) + fine(ii+1, jj , kk-1) + + fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) + ) + + fine(ii, jj+1, kk) + + 0.5 * ( + fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) + + fine(ii , jj+1, kk-1) + fine(ii , jj+1, kk+1) + ) + + 0.25 * ( + fine(ii-1, jj+1, kk-1) + fine(ii+1, jj+1, kk-1) + + fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) + ) + ); + } else { + coarse(i, j, k) = 0.125 * ( + fine(ii, jj, kk ) + + 0.5 * ( + fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + + fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) + ) + + 0.25 * ( + fine(ii-1, jj-1, kk ) + fine(ii+1, jj-1, kk ) + + fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) + ) + + fine(ii, jj, kk+1) + + 0.5 * ( + fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) + + fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) + ) + + 0.25 * ( + fine(ii-1, jj-1, kk+1) + fine(ii+1, jj-1, kk+1) + + fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) + ) + ); + } +#endif + } +private: + amrex::Array4< amrex::Real const > m_fine; + amrex::Array4< amrex::Real > m_coarse; + int m_refinement_ratio; +}; + +#endif // INTERPOLATECURRENTFINETOCOARSE_H diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 40f8203a9..4f870e79c 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -2,7 +2,7 @@ #include #include #include -#include +#include #include #include @@ -337,7 +337,7 @@ WarpX::SyncCurrent () std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, refinement_ratio[0]); + interpolateCurrentFineToCoarse(fine, crse, refinement_ratio[0]); } // For each level @@ -350,10 +350,11 @@ WarpX::SyncCurrent () } void -WarpX::SyncCurrent (const std::array& fine, - const std::array< amrex::MultiFab*,3>& coarse, - int const refinement_ratio) +WarpX::interpolateCurrentFineToCoarse ( std::array< amrex::MultiFab const *, 3 > const & fine, + std::array< amrex::MultiFab *, 3 > const & coarse, + int const refinement_ratio) { + BL_PROFILE("InterpolateCurrentFineToCoarse()"); BL_ASSERT(refinement_ratio == 2); const IntVect& ng = (fine[0]->nGrowVect() + 1) / refinement_ratio; // add equivalent no. of guards to coarse patch @@ -372,11 +373,11 @@ WarpX::SyncCurrent (const std::array& fine, auto const & arrCoarse = coarse[idim]->array(mfi); if( idim == 0 ) - amrex::ParallelFor( bx, WarpxSyncCurrent<0>(arrFine, arrCoarse, refinement_ratio) ); + amrex::ParallelFor( bx, InterpolateCurrentFineToCoarse<0>(arrFine, arrCoarse, refinement_ratio) ); else if( idim == 1 ) - amrex::ParallelFor( bx, WarpxSyncCurrent<1>(arrFine, arrCoarse, refinement_ratio) ); + amrex::ParallelFor( bx, InterpolateCurrentFineToCoarse<1>(arrFine, arrCoarse, refinement_ratio) ); else if( idim == 2 ) - amrex::ParallelFor( bx, WarpxSyncCurrent<2>(arrFine, arrCoarse, refinement_ratio) ); + amrex::ParallelFor( bx, InterpolateCurrentFineToCoarse<2>(arrFine, arrCoarse, refinement_ratio) ); } } } @@ -452,7 +453,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, refinement_ratio[0]); + interpolateCurrentFineToCoarse(fine, crse, refinement_ratio[0]); } void diff --git a/Source/WarpX.H b/Source/WarpX.H index 9b9fb045c..c59802427 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -428,9 +428,9 @@ private: * \param[out] coarse coarse patches to interpolate to * \param[in] refinement_ratio integer ratio between the two */ - void SyncCurrent (const std::array& fine, - const std::array< amrex::MultiFab*,3>& coarse, - int const refinement_ratio); + void interpolateCurrentFineToCoarse (std::array< amrex::MultiFab const *, 3 > const & fine, + std::array< amrex::MultiFab *, 3 > const & coarse, + int const refinement_ratio); /** \brief Fills the values of the charge density on the coarse patch by * averaging the values of the charge density of the fine patch (on the same level). -- cgit v1.2.3