diff options
Diffstat (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp')
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 360 |
1 files changed, 297 insertions, 63 deletions
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index f08b44adb..21000424d 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -6,6 +6,7 @@ #include <WarpXConst.H> #include <WarpX_f.H> #include <WarpX_K.H> +#include <WarpX_FDTD.H> #ifdef WARPX_USE_PY #include <WarpX_py.H> #endif @@ -16,6 +17,75 @@ using namespace amrex; +#ifdef WARPX_USE_PSATD +namespace { + void + PushPSATDSinglePatch ( + SpectralSolver& solver, + std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, + std::array<std::unique_ptr<amrex::MultiFab>,3>& Bfield, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + std::unique_ptr<amrex::MultiFab>& rho ) { + + using Idx = SpectralFieldIndex; + + // Perform forward Fourier transform + solver.ForwardTransform(*Efield[0], Idx::Ex); + solver.ForwardTransform(*Efield[1], Idx::Ey); + solver.ForwardTransform(*Efield[2], Idx::Ez); + solver.ForwardTransform(*Bfield[0], Idx::Bx); + solver.ForwardTransform(*Bfield[1], Idx::By); + solver.ForwardTransform(*Bfield[2], Idx::Bz); + solver.ForwardTransform(*current[0], Idx::Jx); + solver.ForwardTransform(*current[1], Idx::Jy); + solver.ForwardTransform(*current[2], Idx::Jz); + solver.ForwardTransform(*rho, Idx::rho_old, 0); + solver.ForwardTransform(*rho, Idx::rho_new, 1); + // Advance fields in spectral space + solver.pushSpectralFields(); + // Perform backward Fourier Transform + solver.BackwardTransform(*Efield[0], Idx::Ex); + solver.BackwardTransform(*Efield[1], Idx::Ey); + solver.BackwardTransform(*Efield[2], Idx::Ez); + solver.BackwardTransform(*Bfield[0], Idx::Bx); + solver.BackwardTransform(*Bfield[1], Idx::By); + solver.BackwardTransform(*Bfield[2], Idx::Bz); + } +} + +void +WarpX::PushPSATD (amrex::Real a_dt) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent"); + if (fft_hybrid_mpi_decomposition){ +#ifndef AMREX_USE_CUDA // Only available on CPU + PushPSATD_hybridFFT(lev, a_dt); +#endif + } else { + PushPSATD_localFFT(lev, a_dt); + } + + // Evolve the fields in the PML boxes + if (do_pml && pml[lev]->ok()) { + pml[lev]->PushPSATD(); + } + } +} + +void +WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) +{ + // Update the fields on the fine and coarse patch + PushPSATDSinglePatch( *spectral_solver_fp[lev], + Efield_fp[lev], Bfield_fp[lev], current_fp[lev], rho_fp[lev] ); + if (spectral_solver_cp[lev]) { + PushPSATDSinglePatch( *spectral_solver_cp[lev], + Efield_cp[lev], Bfield_cp[lev], current_cp[lev], rho_cp[lev] ); + } +} +#endif + void WarpX::EvolveB (Real a_dt) { @@ -40,7 +110,8 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; + const Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; + const Real dxinv = 1./dx[0]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; if (patch_type == PatchType::fine) @@ -81,13 +152,13 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) const Box& tby = mfi.tilebox(By_nodal_flag); const Box& tbz = mfi.tilebox(Bz_nodal_flag); + auto const& Bxfab = Bx->array(mfi); + auto const& Byfab = By->array(mfi); + auto const& Bzfab = Bz->array(mfi); + auto const& Exfab = Ex->array(mfi); + auto const& Eyfab = Ey->array(mfi); + auto const& Ezfab = Ez->array(mfi); if (do_nodal) { - auto const& Bxfab = Bx->array(mfi); - auto const& Byfab = By->array(mfi); - auto const& Bzfab = Bz->array(mfi); - auto const& Exfab = Ex->array(mfi); - auto const& Eyfab = Ey->array(mfi); - auto const& Ezfab = Ez->array(mfi); amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int j, int k, int l) { @@ -103,22 +174,54 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); }); - } else { - // Call picsar routine for each tile - warpx_push_bvec( - tbx.loVect(), tbx.hiVect(), - tby.loVect(), tby.hiVect(), - tbz.loVect(), tbz.hiVect(), - &n_rz_azimuthal_modes, - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - &dtsdx, &dtsdy, &dtsdz, - &xmin, &dx[0], - &WarpX::maxwell_fdtd_solver_id); + } else if (WarpX::maxwell_fdtd_solver_id == 0) { + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); + }); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); + }); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin); + }); + } else if (WarpX::maxwell_fdtd_solver_id == 1) { + Real betaxy, betaxz, betayx, betayz, betazx, betazy; + Real gammax, gammay, gammaz; + Real alphax, alphay, alphaz; + warpx_calculate_ckc_coefficients(dtsdx, dtsdy, dtsdz, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + amrex::ParallelFor(tbx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bx_ckc(j,k,l,Bxfab,Eyfab,Ezfab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tby, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_by_ckc(j,k,l,Byfab,Exfab,Ezfab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_bz_ckc(j,k,l,Bzfab,Exfab,Eyfab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); } if (cost) { @@ -190,9 +293,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt; const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt; - int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2]; + const Real dtsdx_c2 = c2dt/dx[0], dtsdy_c2 = c2dt/dx[1], dtsdz_c2 = c2dt/dx[2]; + const Real dxinv = 1./dx[0]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F; if (patch_type == PatchType::fine) @@ -241,16 +345,17 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const Box& tey = mfi.tilebox(Ey_nodal_flag); const Box& tez = mfi.tilebox(Ez_nodal_flag); + auto const& Exfab = Ex->array(mfi); + auto const& Eyfab = Ey->array(mfi); + auto const& Ezfab = Ez->array(mfi); + auto const& Bxfab = Bx->array(mfi); + auto const& Byfab = By->array(mfi); + auto const& Bzfab = Bz->array(mfi); + auto const& jxfab = jx->array(mfi); + auto const& jyfab = jy->array(mfi); + auto const& jzfab = jz->array(mfi); + if (do_nodal) { - auto const& Exfab = Ex->array(mfi); - auto const& Eyfab = Ey->array(mfi); - auto const& Ezfab = Ez->array(mfi); - auto const& Bxfab = Bx->array(mfi); - auto const& Byfab = By->array(mfi); - auto const& Bzfab = Bz->array(mfi); - auto const& jxfab = jx->array(mfi); - auto const& jyfab = jy->array(mfi); - auto const& jzfab = jz->array(mfi); amrex::ParallelFor(tex, [=] AMREX_GPU_DEVICE (int j, int k, int l) { @@ -267,38 +372,76 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2); }); } else { - // Call picsar routine for each tile - warpx_push_evec( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - &n_rz_azimuthal_modes, - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*Bx)[mfi]), - BL_TO_FORTRAN_3D((*By)[mfi]), - BL_TO_FORTRAN_3D((*Bz)[mfi]), - BL_TO_FORTRAN_3D((*jx)[mfi]), - BL_TO_FORTRAN_3D((*jy)[mfi]), - BL_TO_FORTRAN_3D((*jz)[mfi]), - &mu_c2_dt, - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, - &xmin, &dx[0]); + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); + }); + amrex::ParallelFor(tey, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin); + }); + amrex::ParallelFor(tez, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin); + }); } if (F) { - warpx_push_evec_f( - tex.loVect(), tex.hiVect(), - tey.loVect(), tey.hiVect(), - tez.loVect(), tez.hiVect(), - BL_TO_FORTRAN_3D((*Ex)[mfi]), - BL_TO_FORTRAN_3D((*Ey)[mfi]), - BL_TO_FORTRAN_3D((*Ez)[mfi]), - BL_TO_FORTRAN_3D((*F)[mfi]), - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, - &WarpX::maxwell_fdtd_solver_id); + auto const& Ffab = F->array(mfi); + if (WarpX::maxwell_fdtd_solver_id == 0) { + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_f_yee(j,k,l,Exfab,Ffab,dtsdx_c2); + }); + amrex::ParallelFor(tey, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ey_f_yee(j,k,l,Eyfab,Ffab,dtsdy_c2); + }); + amrex::ParallelFor(tez, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ez_f_yee(j,k,l,Ezfab,Ffab,dtsdz_c2); + }); + } + else if (WarpX::maxwell_fdtd_solver_id == 1) { + Real betaxy, betaxz, betayx, betayz, betazx, betazy; + Real gammax, gammay, gammaz; + Real alphax, alphay, alphaz; + warpx_calculate_ckc_coefficients(dtsdx_c2, dtsdy_c2, dtsdz_c2, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_f_ckc(j,k,l,Exfab,Ffab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tey, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ey_f_ckc(j,k,l,Eyfab,Ffab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + amrex::ParallelFor(tez, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ez_f_ckc(j,k,l,Ezfab,Ffab, + betaxy, betaxz, betayx, betayz, betazx, betazy, + gammax, gammay, gammaz, + alphax, alphay, alphaz); + }); + } } if (cost) { @@ -388,7 +531,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c; - int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; + const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const auto& dx = WarpX::CellSize(patch_level); const std::array<Real,3> dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]}; @@ -438,3 +581,94 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) } } +#ifdef WARPX_DIM_RZ +// This scales the current by the inverse volume and wraps around the depostion at negative radius. +// It is faster to apply this on the grid than to do it particle by particle. +// It is put here since there isn't another nice place for it. +void +WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, MultiFab* Jz, int lev) +{ + const long ngJ = Jx->nGrow(); + const std::array<Real,3>& dx = WarpX::CellSize(lev); + const Real dr = dx[0]; + + Box tilebox; + + for ( MFIter mfi(*Jx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + + Array4<Real> const& Jr_arr = Jx->array(mfi); + Array4<Real> const& Jt_arr = Jy->array(mfi); + Array4<Real> const& Jz_arr = Jz->array(mfi); + + tilebox = mfi.tilebox(); + Box tbr = convert(tilebox, WarpX::jx_nodal_flag); + Box tbt = convert(tilebox, WarpX::jy_nodal_flag); + Box tbz = convert(tilebox, WarpX::jz_nodal_flag); + + // Lower corner of tile box physical domain + // Note that this is done before the tilebox.grow so that + // these do not include the guard cells. + const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, lev); + const Dim3 lo = lbound(tilebox); + const Real rmin = xyzmin[0]; + const int irmin = lo.x; + + // Rescale current in r-z mode since the inverse volume factor was not + // included in the current deposition. + amrex::ParallelFor(tbr, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Note that Jr(i==0) is at 1/2 dr. + if (rmin == 0. && 0 <= i && i < ngJ) { + Jr_arr(i,j,0) -= Jr_arr(-1-i,j,0); + } + // Apply the inverse volume scaling + // Since Jr is not cell centered in r, no need for distinction + // between on axis and off-axis factors + const amrex::Real r = std::abs(rmin + (i - irmin + 0.5)*dr); + Jr_arr(i,j,0) /= (2.*MathConst::pi*r); + }); + amrex::ParallelFor(tbt, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Jt is located on the boundary + if (rmin == 0. && 0 < i && i <= ngJ) { + Jt_arr(i,j,0) += Jt_arr(-i,j,0); + } + + // Apply the inverse volume scaling + // Jt is forced to zero on axis. + const amrex::Real r = std::abs(rmin + (i - irmin)*dr); + if (r == 0.) { + Jt_arr(i,j,0) = 0.; + } else { + Jt_arr(i,j,0) /= (2.*MathConst::pi*r); + } + }); + amrex::ParallelFor(tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Wrap the current density deposited in the guard cells around + // to the cells above the axis. + // Jz is located on the boundary + if (rmin == 0. && 0 < i && i <= ngJ) { + Jz_arr(i,j,0) += Jz_arr(-i,j,0); + } + + // Apply the inverse volume scaling + const amrex::Real r = std::abs(rmin + (i - irmin)*dr); + if (r == 0.) { + // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis + Jz_arr(i,j,0) /= (MathConst::pi*dr/3.); + } else { + Jz_arr(i,j,0) /= (2.*MathConst::pi*r); + } + }); + } +} +#endif |