diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Evolve/WarpXEvolve.cpp | 15 | ||||
| -rw-r--r-- | Source/Initialization/WarpXInitData.cpp | 5 | ||||
| -rw-r--r-- | Source/Parallelization/WarpXComm.cpp | 118 | ||||
| -rw-r--r-- | Source/Parallelization/WarpXComm_K.H | 148 | ||||
| -rw-r--r-- | Source/WarpX.H | 79 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 200 |
6 files changed, 407 insertions, 158 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 39b51b261..e632409f6 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -21,7 +21,6 @@ #include <cmath> #include <limits> - using namespace amrex; void @@ -303,6 +302,7 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_particlescraper) warpx_py_particlescraper(); if (warpx_py_beforedeposition) warpx_py_beforedeposition(); PushParticlesandDepose(cur_time); + if (warpx_py_afterdeposition) warpx_py_afterdeposition(); // Synchronize J and rho @@ -317,7 +317,6 @@ WarpX::OneStep_nosub (Real cur_time) VayDeposition(); } - // At this point, J is up-to-date inside the domain, and E and B are // up-to-date including enough guard cells for first step of the field // solve. @@ -601,17 +600,27 @@ WarpX::PushParticlesandDepose (amrex::Real cur_time) void WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type) { + // If warpx.do_current_centering = 1, the current is deposited on the nodal MultiFab current_fp_nodal + // and then centered onto the staggered MultiFab current_fp + amrex::MultiFab* current_x = (WarpX::do_current_centering) ? current_fp_nodal[lev][0].get() + : current_fp[lev][0].get(); + amrex::MultiFab* current_y = (WarpX::do_current_centering) ? current_fp_nodal[lev][1].get() + : current_fp[lev][1].get(); + amrex::MultiFab* current_z = (WarpX::do_current_centering) ? current_fp_nodal[lev][2].get() + : current_fp[lev][2].get(); + mypc->Evolve(lev, *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2], *Efield_avg_aux[lev][0],*Efield_avg_aux[lev][1],*Efield_avg_aux[lev][2], *Bfield_avg_aux[lev][0],*Bfield_avg_aux[lev][1],*Bfield_avg_aux[lev][2], - *current_fp[lev][0],*current_fp[lev][1],*current_fp[lev][2], + *current_x, *current_y, *current_z, current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), rho_fp[lev].get(), charge_buf[lev].get(), Efield_cax[lev][0].get(), Efield_cax[lev][1].get(), Efield_cax[lev][2].get(), Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(), cur_time, dt[lev], a_dt_type); + #ifdef WARPX_DIM_RZ // This is called after all particles have deposited their current and charge. ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev); diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 1852d871f..f07995c84 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -316,6 +316,11 @@ WarpX::InitLevelData (int lev, Real /*time*/) } } + if (WarpX::do_current_centering) + { + current_fp_nodal[lev][i]->setVal(0.0); + } + if (B_ext_grid_s == "constant" || B_ext_grid_s == "default") { Bfield_fp[lev][i]->setVal(B_external_grid[i]); if (fft_do_time_averaging) { 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; } diff --git a/Source/WarpX.H b/Source/WarpX.H index 2c4b9c20e..4ba4c9d43 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -142,6 +142,9 @@ public: static amrex::Vector<int> particle_boundary_hi; + // If true, the current is deposited on a nodal grid and then centered onto a staggered grid + static bool do_current_centering; + // PSATD: If true (overwritten by the user in the input file), the current correction // defined in equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010 is applied bool current_correction = false; @@ -158,11 +161,15 @@ public: static long noy; static long noz; - // For momentum-conserving field gathering, order of interpolation from the - // staggered positions to the grid nodes - static int field_gathering_nox; - static int field_gathering_noy; - static int field_gathering_noz; + // Order of finite-order centering of fields (staggered to nodal) + static int field_centering_nox; + static int field_centering_noy; + static int field_centering_noz; + + // Order of finite-order centering of currents (nodal to staggered) + static int current_centering_nox; + static int current_centering_noy; + static int current_centering_noz; // Number of modes for the RZ multimode version static int n_rz_azimuthal_modes; @@ -441,6 +448,16 @@ public: void UpdateAuxilaryDataStagToNodal (); void UpdateAuxilaryDataSameType (); + /** + * \brief This function is called if \c warpx.do_current_centering = 1 and + * it centers the currents from a nodal grid to a staggered grid (Yee) using + * finite-order interpolation based on the Fornberg coefficients. + * + * \param[in,out] dst destination \c MultiFab where the results of the finite-order centering are stored + * \param[in] src source \c MultiFab that contains the values of the nodal current to be centered + */ + void UpdateCurrentNodalToStag (amrex::MultiFab& dst, amrex::MultiFab const& src); + // Fill boundary cells including coarse/fine boundaries void FillBoundaryB (amrex::IntVect ng); void FillBoundaryE (amrex::IntVect ng); @@ -605,13 +622,14 @@ public: void ApplyFilterandSumBoundaryRho (int lev, int glev, amrex::MultiFab& rho, int icomp, int ncomp); #ifdef WARPX_USE_PSATD - // Host and device vectors for ordered Fornberg stencil coefficients used for finite-order centering - amrex::Vector<amrex::Real> host_centering_stencil_coeffs_x; - amrex::Vector<amrex::Real> host_centering_stencil_coeffs_y; - amrex::Vector<amrex::Real> host_centering_stencil_coeffs_z; - amrex::Gpu::DeviceVector<amrex::Real> device_centering_stencil_coeffs_x; - amrex::Gpu::DeviceVector<amrex::Real> device_centering_stencil_coeffs_y; - amrex::Gpu::DeviceVector<amrex::Real> device_centering_stencil_coeffs_z; + // Device vectors of stencil coefficients used for finite-order centering of fields + amrex::Gpu::DeviceVector<amrex::Real> device_field_centering_stencil_coeffs_x; + amrex::Gpu::DeviceVector<amrex::Real> device_field_centering_stencil_coeffs_y; + amrex::Gpu::DeviceVector<amrex::Real> device_field_centering_stencil_coeffs_z; + // Device vectors of stencil coefficients used for finite-order centering of currents + amrex::Gpu::DeviceVector<amrex::Real> device_current_centering_stencil_coeffs_x; + amrex::Gpu::DeviceVector<amrex::Real> device_current_centering_stencil_coeffs_y; + amrex::Gpu::DeviceVector<amrex::Real> device_current_centering_stencil_coeffs_z; #endif protected: @@ -760,9 +778,37 @@ private: return gather_buffer_masks[lev].get(); } - void ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coeffs, - amrex::Vector<amrex::Real>& unordered_coeffs, - const int order); +#ifdef WARPX_USE_PSATD + /** + * \brief Re-orders the Fornberg coefficients so that they can be used more conveniently for + * finite-order centering operations. For example, for finite-order centering of order 6, + * the Fornberg coefficients \c (c_0,c_1,c_2) are re-ordered as \c (c_2,c_1,c_0,c_0,c_1,c_2). + * + * \param[in,out] ordered_coeffs host vector where the re-ordered Fornberg coefficients will be stored + * \param[in] unordered_coeffs host vector storing the original sequence of Fornberg coefficients + * \param[in] order order of the finite-order centering along a given direction + */ + void ReorderFornbergCoefficients (amrex::Vector<amrex::Real>& ordered_coeffs, + amrex::Vector<amrex::Real>& unordered_coeffs, + const int order); + /** + * \brief Allocates and initializes the stencil coefficients used for the finite-order centering + * of fields and currents, and stores them in the given device vectors. + * + * \param[in,out] device_centering_stencil_coeffs_x device vector where the stencil coefficients along x will be stored + * \param[in,out] device_centering_stencil_coeffs_y device vector where the stencil coefficients along y will be stored + * \param[in,out] device_centering_stencil_coeffs_z device vector where the stencil coefficients along z will be stored + * \param[in] centering_nox order of the finite-order centering along x + * \param[in] centering_noy order of the finite-order centering along y + * \param[in] centering_noz order of the finite-order centering along z + */ + void AllocateCenteringCoefficients (amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_x, + amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_y, + amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_z, + const int centering_nox, + const int centering_noy, + const int centering_noz); +#endif void AllocLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const amrex::IntVect& ngE, const amrex::IntVect& ngJ, @@ -805,6 +851,9 @@ private: // store fine patch amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_store; + // Nodal MultiFab for nodal current deposition if warpx.do_current_centering = 1 + amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>> current_fp_nodal; + // Coarse patch amrex::Vector< std::unique_ptr<amrex::MultiFab> > F_cp; amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_cp; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 92a2efb61..7417921ab 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -83,6 +83,8 @@ amrex::Vector<int> WarpX::field_boundary_hi(AMREX_SPACEDIM,0); amrex::Vector<int> WarpX::particle_boundary_lo(AMREX_SPACEDIM,0); amrex::Vector<int> WarpX::particle_boundary_hi(AMREX_SPACEDIM,0); +bool WarpX::do_current_centering = false; + int WarpX::n_rz_azimuthal_modes = 1; int WarpX::ncomps = 1; @@ -90,11 +92,15 @@ long WarpX::nox = 1; long WarpX::noy = 1; long WarpX::noz = 1; -// For momentum-conserving field gathering, order of interpolation from the -// staggered positions to the grid nodes -int WarpX::field_gathering_nox = 2; -int WarpX::field_gathering_noy = 2; -int WarpX::field_gathering_noz = 2; +// Order of finite-order centering of fields (staggered to nodal) +int WarpX::field_centering_nox = 2; +int WarpX::field_centering_noy = 2; +int WarpX::field_centering_noz = 2; + +// Order of finite-order centering of currents (nodal to staggered) +int WarpX::current_centering_nox = 2; +int WarpX::current_centering_noy = 2; +int WarpX::current_centering_noz = 2; bool WarpX::use_fdtd_nci_corr = false; bool WarpX::galerkin_interpolation = true; @@ -229,6 +235,11 @@ WarpX::WarpX () current_store.resize(nlevs_max); + if (do_current_centering) + { + current_fp_nodal.resize(nlevs_max); + } + F_cp.resize(nlevs_max); rho_cp.resize(nlevs_max); current_cp.resize(nlevs_max); @@ -630,6 +641,21 @@ WarpX::ReadParameters () // Only needs to be set with WARPX_DIM_RZ, otherwise defaults to 1 pp_warpx.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes); + + // If true, the current is deposited on a nodal grid and then interpolated onto a Yee grid + pp_warpx.query("do_current_centering", do_current_centering); + + // If do_nodal = 1, Maxwell's equations are solved on a nodal grid and + // the current should not be centered + if (do_nodal) + { + do_current_centering = false; + } + + if ((maxLevel() > 0) && do_current_centering) + { + amrex::Abort("\nFinite-order centering of currents is not implemented with mesh refinement"); + } } { @@ -695,64 +721,50 @@ WarpX::ReadParameters () pp_interpolation.query("noz", noz); #ifdef WARPX_USE_PSATD + if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { + // For momentum-conserving field gathering, read from input the order of // interpolation from the staggered positions to the grid nodes - if (field_gathering_algo == GatheringAlgo::MomentumConserving) { - pp_interpolation.query("field_gathering_nox", field_gathering_nox); - pp_interpolation.query("field_gathering_noy", field_gathering_noy); - pp_interpolation.query("field_gathering_noz", field_gathering_noz); + if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving) { + pp_interpolation.query("field_centering_nox", field_centering_nox); + pp_interpolation.query("field_centering_noy", field_centering_noy); + pp_interpolation.query("field_centering_noz", field_centering_noz); + } + + // Read order of finite-order centering of currents (nodal to staggered) + if (WarpX::do_current_centering) { + pp_interpolation.query("current_centering_nox", current_centering_nox); + pp_interpolation.query("current_centering_noy", current_centering_noy); + pp_interpolation.query("current_centering_noz", current_centering_noz); } - if (maxLevel() > 0) { + if (maxLevel() > 0) + { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - field_gathering_nox == 2 && field_gathering_noy == 2 && field_gathering_noz == 2, - "High-order interpolation (order > 2) is not implemented with mesh refinement"); + field_centering_nox == 2 && field_centering_noy == 2 && field_centering_noz == 2, + "High-order centering of fields (order > 2) is not implemented with mesh refinement"); } - // Host vectors for Fornberg stencil coefficients - amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_x; - amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_y; - amrex::Vector<amrex::Real> host_Fornberg_stencil_coeffs_z; - - host_Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(field_gathering_nox, false); - host_Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(field_gathering_noy, false); - host_Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(field_gathering_noz, false); - - // Host vectors for ordered Fornberg stencil coefficients used for finite-order centering - host_centering_stencil_coeffs_x.resize(field_gathering_nox); - host_centering_stencil_coeffs_y.resize(field_gathering_noy); - host_centering_stencil_coeffs_z.resize(field_gathering_noz); - - // Re-order Fornberg stencil coefficients - // example for order 6: (c_0,c_1,c_2) becomes (c_2,c_1,c_0,c_0,c_1,c_2) - ReorderFornbergCoefficients(host_centering_stencil_coeffs_x, - host_Fornberg_stencil_coeffs_x, field_gathering_nox); - ReorderFornbergCoefficients(host_centering_stencil_coeffs_y, - host_Fornberg_stencil_coeffs_y, field_gathering_noy); - ReorderFornbergCoefficients(host_centering_stencil_coeffs_z, - host_Fornberg_stencil_coeffs_z, field_gathering_noz); - - // Device vectors for ordered Fornberg stencil coefficients used for finite-order centering - device_centering_stencil_coeffs_x.resize(host_centering_stencil_coeffs_x.size()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, - host_centering_stencil_coeffs_x.begin(), - host_centering_stencil_coeffs_x.end(), - device_centering_stencil_coeffs_x.begin()); - - device_centering_stencil_coeffs_y.resize(host_centering_stencil_coeffs_y.size()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, - host_centering_stencil_coeffs_y.begin(), - host_centering_stencil_coeffs_y.end(), - device_centering_stencil_coeffs_y.begin()); - - device_centering_stencil_coeffs_z.resize(host_centering_stencil_coeffs_z.size()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, - host_centering_stencil_coeffs_z.begin(), - host_centering_stencil_coeffs_z.end(), - device_centering_stencil_coeffs_z.begin()); - - amrex::Gpu::synchronize(); + if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving) + { + AllocateCenteringCoefficients(device_field_centering_stencil_coeffs_x, + device_field_centering_stencil_coeffs_y, + device_field_centering_stencil_coeffs_z, + field_centering_nox, + field_centering_noy, + field_centering_noz); + } + + if (WarpX::do_current_centering) + { + AllocateCenteringCoefficients(device_current_centering_stencil_coeffs_x, + device_current_centering_stencil_coeffs_y, + device_current_centering_stencil_coeffs_z, + current_centering_nox, + current_centering_noy, + current_centering_noz); + } } #endif @@ -1031,6 +1043,11 @@ WarpX::ClearLevel (int lev) current_store[lev][i].reset(); + if (do_current_centering) + { + current_fp_nodal[lev][i].reset(); + } + current_cp[lev][i].reset(); Efield_cp [lev][i].reset(); Bfield_cp [lev][i].reset(); @@ -1209,6 +1226,14 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm current_fp[lev][1] = std::make_unique<MultiFab>(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ,tag("current_fp[y]")); current_fp[lev][2] = std::make_unique<MultiFab>(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ,tag("current_fp[z]")); + if (do_current_centering) + { + amrex::BoxArray const& nodal_ba = amrex::convert(ba, amrex::IntVect::TheNodeVector()); + current_fp_nodal[lev][0] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ); + current_fp_nodal[lev][1] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ); + current_fp_nodal[lev][2] = std::make_unique<MultiFab>(nodal_ba, dm, ncomps, ngJ); + } + Bfield_avg_fp[lev][0] = std::make_unique<MultiFab>(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[x]")); Bfield_avg_fp[lev][1] = std::make_unique<MultiFab>(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[y]")); Bfield_avg_fp[lev][2] = std::make_unique<MultiFab>(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE,tag("Bfield_avg_fp[z]")); @@ -1768,9 +1793,10 @@ WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_ma #endif } -void WarpX::ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coeffs, - amrex::Vector<amrex::Real>& unordered_coeffs, - const int order) +#ifdef WARPX_USE_PSATD +void WarpX::ReorderFornbergCoefficients (amrex::Vector<amrex::Real>& ordered_coeffs, + amrex::Vector<amrex::Real>& unordered_coeffs, + const int order) { const int n = order / 2; for (int i = 0; i < n; i++) { @@ -1781,6 +1807,64 @@ void WarpX::ReorderFornbergCoefficients(amrex::Vector<amrex::Real>& ordered_coef } } +void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_x, + amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_y, + amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_z, + const int centering_nox, + const int centering_noy, + const int centering_noz) +{ + // Vectors of Fornberg stencil coefficients + amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_x; + amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_y; + amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_z; + + // Host vectors of stencil coefficients used for finite-order centering + amrex::Vector<amrex::Real> host_centering_stencil_coeffs_x; + amrex::Vector<amrex::Real> host_centering_stencil_coeffs_y; + amrex::Vector<amrex::Real> host_centering_stencil_coeffs_z; + + Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(centering_nox, false); + Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(centering_noy, false); + Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(centering_noz, false); + + host_centering_stencil_coeffs_x.resize(centering_nox); + host_centering_stencil_coeffs_y.resize(centering_noy); + host_centering_stencil_coeffs_z.resize(centering_noz); + + // Re-order Fornberg stencil coefficients: + // example for order 6: (c_0,c_1,c_2) becomes (c_2,c_1,c_0,c_0,c_1,c_2) + ReorderFornbergCoefficients(host_centering_stencil_coeffs_x, + Fornberg_stencil_coeffs_x, centering_nox); + ReorderFornbergCoefficients(host_centering_stencil_coeffs_y, + Fornberg_stencil_coeffs_y, centering_noy); + ReorderFornbergCoefficients(host_centering_stencil_coeffs_z, + Fornberg_stencil_coeffs_z, centering_noz); + + // Device vectors of stencil coefficients used for finite-order centering + + device_centering_stencil_coeffs_x.resize(host_centering_stencil_coeffs_x.size()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_centering_stencil_coeffs_x.begin(), + host_centering_stencil_coeffs_x.end(), + device_centering_stencil_coeffs_x.begin()); + + device_centering_stencil_coeffs_y.resize(host_centering_stencil_coeffs_y.size()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_centering_stencil_coeffs_y.begin(), + host_centering_stencil_coeffs_y.end(), + device_centering_stencil_coeffs_y.begin()); + + device_centering_stencil_coeffs_z.resize(host_centering_stencil_coeffs_z.size()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_centering_stencil_coeffs_z.begin(), + host_centering_stencil_coeffs_z.end(), + device_centering_stencil_coeffs_z.begin()); + + amrex::Gpu::synchronize(); +} +#endif + const iMultiFab* WarpX::CurrentBufferMasks (int lev) { |