diff options
Diffstat (limited to 'Source/Parallelization')
| -rw-r--r-- | Source/Parallelization/WarpXComm.cpp | 118 | ||||
| -rw-r--r-- | Source/Parallelization/WarpXComm_K.H | 148 |
2 files changed, 184 insertions, 82 deletions
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index e5a912c8b..4131b7307 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -56,6 +56,9 @@ WarpX::UpdateAuxilaryDataStagToNodal () const amrex::IntVect& Ey_stag = Emf[0][1]->ixType().toIntVect(); const amrex::IntVect& Ez_stag = Emf[0][2]->ixType().toIntVect(); + // Destination MultiFab (aux) always has nodal index type when this function is called + const amrex::IntVect& dst_stag = amrex::IntVect::TheNodeVector(); + // For level 0, we only need to do the average. #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -82,44 +85,49 @@ WarpX::UpdateAuxilaryDataStagToNodal () Box bx = mfi.fabbox(); if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) { + #ifdef WARPX_USE_PSATD - const int fg_nox = WarpX::field_gathering_nox; - const int fg_noy = WarpX::field_gathering_noy; - const int fg_noz = WarpX::field_gathering_noz; - // Device vectors for Fornberg stencil coefficients used for finite-order centering - amrex::Real const * stencil_coeffs_x = WarpX::device_centering_stencil_coeffs_x.data(); - amrex::Real const * stencil_coeffs_y = WarpX::device_centering_stencil_coeffs_y.data(); - amrex::Real const * stencil_coeffs_z = WarpX::device_centering_stencil_coeffs_z.data(); + + // Order of finite-order centering of fields + const int fg_nox = WarpX::field_centering_nox; + const int fg_noy = WarpX::field_centering_noy; + const int fg_noz = WarpX::field_centering_noz; + + // Device vectors of stencil coefficients used for finite-order centering of fields + amrex::Real const * stencil_coeffs_x = WarpX::device_field_centering_stencil_coeffs_x.data(); + amrex::Real const * stencil_coeffs_y = WarpX::device_field_centering_stencil_coeffs_y.data(); + amrex::Real const * stencil_coeffs_z = WarpX::device_field_centering_stencil_coeffs_z.data(); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept { - warpx_interp<true>(j, k, l, bx_aux, bx_fp, Bx_stag, fg_nox, fg_noy, fg_noz, + warpx_interp<true>(j, k, l, bx_aux, bx_fp, dst_stag, Bx_stag, fg_nox, fg_noy, fg_noz, stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); - warpx_interp<true>(j, k, l, by_aux, by_fp, By_stag, fg_nox, fg_noy, fg_noz, + warpx_interp<true>(j, k, l, by_aux, by_fp, dst_stag, By_stag, fg_nox, fg_noy, fg_noz, stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); - warpx_interp<true>(j, k, l, bz_aux, bz_fp, Bz_stag, fg_nox, fg_noy, fg_noz, + warpx_interp<true>(j, k, l, bz_aux, bz_fp, dst_stag, Bz_stag, fg_nox, fg_noy, fg_noz, stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); - warpx_interp<true>(j, k, l, ex_aux, ex_fp, Ex_stag, fg_nox, fg_noy, fg_noz, + warpx_interp<true>(j, k, l, ex_aux, ex_fp, dst_stag, Ex_stag, fg_nox, fg_noy, fg_noz, stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); - warpx_interp<true>(j, k, l, ey_aux, ey_fp, Ey_stag, fg_nox, fg_noy, fg_noz, + warpx_interp<true>(j, k, l, ey_aux, ey_fp, dst_stag, Ey_stag, fg_nox, fg_noy, fg_noz, stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); - warpx_interp<true>(j, k, l, ez_aux, ez_fp, Ez_stag, fg_nox, fg_noy, fg_noz, + warpx_interp<true>(j, k, l, ez_aux, ez_fp, dst_stag, Ez_stag, fg_nox, fg_noy, fg_noz, stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); }); #endif } else { // FDTD amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept { - warpx_interp(j, k, l, bx_aux, bx_fp, Bx_stag); - warpx_interp(j, k, l, by_aux, by_fp, By_stag); - warpx_interp(j, k, l, bz_aux, bz_fp, Bz_stag); - warpx_interp(j, k, l, ex_aux, ex_fp, Ex_stag); - warpx_interp(j, k, l, ey_aux, ey_fp, Ey_stag); - warpx_interp(j, k, l, ez_aux, ez_fp, Ez_stag); + warpx_interp(j, k, l, bx_aux, bx_fp, dst_stag, Bx_stag); + warpx_interp(j, k, l, by_aux, by_fp, dst_stag, By_stag); + warpx_interp(j, k, l, bz_aux, bz_fp, dst_stag, Bz_stag); + warpx_interp(j, k, l, ex_aux, ex_fp, dst_stag, Ex_stag); + warpx_interp(j, k, l, ey_aux, ey_fp, dst_stag, Ey_stag); + warpx_interp(j, k, l, ez_aux, ez_fp, dst_stag, Ez_stag); }); } } @@ -360,6 +368,67 @@ WarpX::UpdateAuxilaryDataSameType () } } +void WarpX::UpdateCurrentNodalToStag (amrex::MultiFab& dst, amrex::MultiFab const& src) +{ + // If source and destination MultiFabs have the same index type, a simple copy is enough + // (for example, this happens with the current along y in 2D, which is always fully nodal) + if (dst.ixType() == src.ixType()) + { + amrex::MultiFab::Copy(dst, src, 0, 0, dst.nComp(), dst.nGrowVect()); + return; + } + + amrex::IntVect const& dst_stag = dst.ixType().toIntVect(); + + // Source MultiFab always has nodal index type when this function is called + amrex::IntVect const& src_stag = amrex::IntVect::TheNodeVector(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + + for (MFIter mfi(dst); mfi.isValid(); ++mfi) + { + // Loop over full box including ghost cells + // (input arrays will be padded with zeros beyond ghost cells + // for out-of-bound accesses due to large-stencil operations) + Box bx = mfi.fabbox(); + + amrex::Array4<amrex::Real const> const& src_arr = src.const_array(mfi); + amrex::Array4<amrex::Real> const& dst_arr = dst.array(mfi); + + if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) + { +#ifdef WARPX_USE_PSATD + + // Order of finite-order centering of currents + const int cc_nox = WarpX::current_centering_nox; + const int cc_noy = WarpX::current_centering_noy; + const int cc_noz = WarpX::current_centering_noz; + + // Device vectors of stencil coefficients used for finite-order centering of currents + amrex::Real const * stencil_coeffs_x = WarpX::device_current_centering_stencil_coeffs_x.data(); + amrex::Real const * stencil_coeffs_y = WarpX::device_current_centering_stencil_coeffs_y.data(); + amrex::Real const * stencil_coeffs_z = WarpX::device_current_centering_stencil_coeffs_z.data(); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp<true>(j, k, l, dst_arr, src_arr, dst_stag, src_stag, cc_nox, cc_noy, cc_noz, + stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); + }); +#endif + } + + else // FDTD + { + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp<false>(j, k, l, dst_arr, src_arr, dst_stag, src_stag); + }); + } + } +} + void WarpX::FillBoundaryB (IntVect ng) { @@ -710,6 +779,17 @@ WarpX::SyncCurrent () { WARPX_PROFILE("WarpX::SyncCurrent()"); + // If warpx.do_current_centering = 1, center currents from nodal grid to staggered grid + if (WarpX::do_current_centering) + { + for (int lev = 0; lev <= finest_level; lev++) + { + WarpX::UpdateCurrentNodalToStag(*current_fp[lev][0], *current_fp_nodal[lev][0]); + WarpX::UpdateCurrentNodalToStag(*current_fp[lev][1], *current_fp_nodal[lev][1]); + WarpX::UpdateCurrentNodalToStag(*current_fp[lev][2], *current_fp_nodal[lev][2]); + } + } + // Restrict fine patch current onto the coarse patch, before // summing the guard cells of the fine patch for (int lev = 1; lev <= finest_level; ++lev) diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index 64c5d0f7c..d1600ec7a 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -487,24 +487,24 @@ void warpx_interp_nd_efield_z (int j, int k, int l, } /** - * \brief Arbitrary-order interpolation function used to interpolate a given field from the Yee grid - * to the nodal grid, before gathering the field from the nodes to the particle positions - * (momentum-conserving field gathering). With the FDTD solver, this performs simple linear interpolation. + * \brief Arbitrary-order interpolation function used to center a given MultiFab between two grids + * with different staggerings. With the FDTD solver, this performs simple linear interpolation. * With the PSATD solver, this performs arbitrary-order interpolation based on the Fornberg coefficients. - * The result is stored in the output array \c arr_aux, allocated on the \c aux patch. + * The result is stored in the output array \c dst_arr. * * \param[in] j index along x of the output array * \param[in] k index along y (in 3D) or z (in 2D) of the output array * \param[in] l index along z (in 3D, \c l = 0 in 2D) of the output array - * \param[in,out] dst_arr output array allocated on the \c aux patch where interpolated values are stored - * \param[in] src_arr input array allocated on the fine patch - * \param[in] src_stag \c IndexType of the input array (Yee staggering) - * \param[in] nox order of finite-order interpolation along x - * \param[in] noy order of finite-order interpolation along y - * \param[in] noz order of finite-order interpolation along z - * \param[in] stencil_coeffs_x array of ordered Fornberg coefficients for finite-order interpolation stencil along x - * \param[in] stencil_coeffs_y array of ordered Fornberg coefficients for finite-order interpolation stencil along y - * \param[in] stencil_coeffs_z array of ordered Fornberg coefficients for finite-order interpolation stencil along z + * \param[in,out] dst_arr output array where interpolated values are stored + * \param[in] src_arr input array storing the values used for interpolation + * \param[in] dst_stag \c IndexType of the output array + * \param[in] src_stag \c IndexType of the input array + * \param[in] nox order of finite-order centering along x + * \param[in] noy order of finite-order centering along y + * \param[in] noz order of finite-order centering along z + * \param[in] stencil_coeffs_x array of ordered Fornberg coefficients for finite-order centering stencil along x + * \param[in] stencil_coeffs_y array of ordered Fornberg coefficients for finite-order centering stencil along y + * \param[in] stencil_coeffs_z array of ordered Fornberg coefficients for finite-order centering stencil along z */ template< bool IS_PSATD = false > AMREX_GPU_DEVICE AMREX_FORCE_INLINE @@ -513,6 +513,7 @@ void warpx_interp (const int j, const int l, amrex::Array4<amrex::Real > const& dst_arr, amrex::Array4<amrex::Real const> const& src_arr, + const amrex::IntVect& dst_stag, const amrex::IntVect& src_stag, const int nox = 2, const int noy = 2, @@ -540,11 +541,18 @@ void warpx_interp (const int j, amrex::ignore_unused(stencil_coeffs_x, stencil_coeffs_y, stencil_coeffs_z); #endif + // If dst_nodal = true , we are centering from a staggered grid to a nodal grid + // If dst_nodal = false, we are centering from a nodal grid to a staggered grid + const bool dst_nodal = (dst_stag == amrex::IntVect::TheNodeVector()); + + // See 1D examples below to understand the meaning of this integer shift + const int shift = (dst_nodal) ? 0 : 1; + // Staggering (s = 0 if cell-centered, s = 1 if nodal) - const int sj = src_stag[0]; - const int sk = src_stag[1]; + const int sj = (dst_nodal) ? src_stag[0] : dst_stag[0]; + const int sk = (dst_nodal) ? src_stag[1] : dst_stag[1]; #if (AMREX_SPACEDIM == 3) - const int sl = src_stag[2]; + const int sl = (dst_nodal) ? src_stag[2] : dst_stag[2]; #endif // Interpolate along j,k,l only if source MultiFab is staggered along j,k,l @@ -572,12 +580,12 @@ void warpx_interp (const int j, #endif // Min and max for interpolation loop along j - const int jmin = (interp_j) ? j - noj/2 : j; - const int jmax = (interp_j) ? j + noj/2 - 1 : j; + const int jmin = (interp_j) ? j - noj/2 + shift : j; + const int jmax = (interp_j) ? j + noj/2 + shift - 1 : j; // Min and max for interpolation loop along k - const int kmin = (interp_k) ? k - nok/2 : k; - const int kmax = (interp_k) ? k + nok/2 - 1 : k; + const int kmin = (interp_k) ? k - nok/2 + shift : k; + const int kmax = (interp_k) ? k + nok/2 + shift - 1 : k; // Min and max for interpolation loop along l #if (AMREX_SPACEDIM == 2) @@ -585,44 +593,73 @@ void warpx_interp (const int j, const int lmin = l; const int lmax = l; #elif (AMREX_SPACEDIM == 3) - const int lmin = (interp_l) ? l - nol/2 : l; - const int lmax = (interp_l) ? l + nol/2 - 1 : l; + const int lmin = (interp_l) ? l - nol/2 + shift : l; + const int lmax = (interp_l) ? l + nol/2 + shift - 1 : l; #endif + // Number of interpolation points + const int nj = jmax - jmin; + const int nk = kmax - kmin; + const int nl = lmax - lmin; + + // Example of 1D centering from nodal grid to nodal grid (simple copy): + // + // j + // --o-----o-----o-- result(j) = f(j) + // --o-----o-----o-- + // j-1 j j+1 + // + // Example of 1D linear centering from staggered grid to nodal grid: + // + // j + // --o-----o-----o-- result(j) = (f(j-1) + f(j)) / 2 + // -----x-----x----- + // j-1 j + // + // Example of 1D linear centering from nodal grid to staggered grid: + // (note the shift of +1 in the indices with respect to the case above, see variable "shift") + // + // j + // --x-----x-----x-- result(j) = (f(j) + f(j+1)) / 2 + // -----o-----o----- + // j j+1 + // + // Example of 1D finite-order centering from staggered grid to nodal grid: + // + // j + // --o-----o-----o-----o-----o-----o-----o-- result(j) = c_0 * (f(j-1) + f(j) ) / 2 + // -----x-----x-----x-----x-----x-----x----- + c_1 * (f(j-2) + f(j+1)) / 2 + // j-3 j-2 j-1 j j+1 j+2 + c_2 * (f(j-3) + f(j+2)) / 2 + // c_2 c_1 c_0 c_0 c_1 c_2 + ... + // + // Example of 1D finite-order centering from nodal grid to staggered grid: + // (note the shift of +1 in the indices with respect to the case above, see variable "shift") + // + // j + // --x-----x-----x-----x-----x-----x-----x-- result(j) = c_0 * (f(j) + f(j+1)) / 2 + // -----o-----o-----o-----o-----o-----o----- + c_1 * (f(j-1) + f(j+2)) / 2 + // j-2 j-1 j j+1 j+2 j+3 + c_2 * (f(j-2) + f(j+3)) / 2 + // c_2 c_1 c_0 c_0 c_1 c_2 + ... + amrex::Real res = 0.0_rt; - // FDTD (linear interpolation) - if( !IS_PSATD ) { - // Example of 1D linear interpolation from nodal grid to nodal grid: - // - // j - // --o-----o-----o-- result(j) = f(j) - // --o-----o-----o-- - // j-1 j j+1 - // - // Example of 1D linear interpolation from staggered grid to nodal grid: - // - // j - // --o-----o-----o-- result(j) = (f(j-1) + f(j)) / 2 - // -----x-----x----- - // j-1 j - - for (int ll = lmin; ll <= lmax; ll++) + if( !IS_PSATD ) // FDTD (linear interpolation) + { + for (int ll = 0; ll <= nl; ll++) { - for (int kk = kmin; kk <= kmax; kk++) + for (int kk = 0; kk <= nk; kk++) { - for (int jj = jmin; jj <= jmax; jj++) + for (int jj = 0; jj <= nj; jj++) { - res += src_arr_zeropad(jj,kk,ll); + res += src_arr_zeropad(jmin+jj,kmin+kk,lmin+ll); } } } - // PSATD (finite-order interpolation) - } else { - const int nj = jmax - jmin; - const int nk = kmax - kmin; - const int nl = lmax - lmin; + } + + else // PSATD (finite-order interpolation) + { amrex::Real cj = 1.0_rt; amrex::Real ck = 1.0_rt; amrex::Real cl = 1.0_rt; @@ -635,21 +672,6 @@ void warpx_interp (const int j, amrex::Real const* scl = stencil_coeffs_z; #endif - // Example of 1D finite-order interpolation from nodal grid to nodal grid: - // - // j - // --o-----o-----o-- result(j) = f(j) - // --o-----o-----o-- - // j-1 j j+1 - // - // Example of 1D finite-order interpolation from staggered grid to nodal grid: - // - // j - // --o-----o-----o-----o-----o-----o-----o-- result(j) = c_0 * (f(j-1) + f(j) ) / 2 - // -----x-----x-----x-----x-----x-----x----- + c_1 * (f(j-2) + f(j+1)) / 2 - // j-3 j-2 j-1 j j+1 j+2 + c_2 * (f(j-3) + f(j+2)) / 2 - // c_2 c_1 c_0 c_0 c_1 c_2 + ... - for (int ll = 0; ll <= nl; ll++) { #if (AMREX_SPACEDIM == 3) @@ -667,7 +689,7 @@ void warpx_interp (const int j, } } } - } // PSATD (finite-order interpolation) + } dst_arr(j,k,l) = wj * wk * wl * res; } |