From 7a3cec3076d19b8a4dab0aa6c7d80078ef45d7ec Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 25 Jul 2019 11:39:13 -0700 Subject: Correct minor PML bugs --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 3 +++ 1 file changed, 3 insertions(+) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 4fce4717b..d756fca24 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -31,6 +31,9 @@ WarpX::PushPSATD (amrex::Real a_dt) } else { PushPSATD_localFFT(lev, a_dt); } + + // Evolve the fields in the PML boxes + if (do_pml) pml[lev]->PushPSATD(); } } -- cgit v1.2.3 From c40c2e6b85955bdb0ee1bca6401139ce9e46e364 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 26 Jul 2019 16:49:01 -0700 Subject: Only push the PML fields if they have been allocated --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index d756fca24..4cc8e675a 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -33,7 +33,9 @@ WarpX::PushPSATD (amrex::Real a_dt) } // Evolve the fields in the PML boxes - if (do_pml) pml[lev]->PushPSATD(); + if (do_pml && pml[lev]->ok()) { + pml[lev]->PushPSATD(); + } } } -- cgit v1.2.3 From f4fe70dd5e7bea68f9b9fbb9ff4dbbf5afdff17c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 30 Jul 2019 16:25:48 -0700 Subject: Apply spectral solver to the fine and coarse patch --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 60 ++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 23 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 4fce4717b..95e3e78b3 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -34,35 +34,50 @@ WarpX::PushPSATD (amrex::Real a_dt) } } -void WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) -{ - auto& solver = *spectral_solver_fp[lev]; +void +PushPSATDSinglePatch ( + SpectralSolver& solver, + std::array,3>& Efield, + std::array,3>& Bfield, + std::array,3>& current, + std::unique_ptr& rho ) { - // Perform forward Fourier transform - solver.ForwardTransform(*Efield_fp[lev][0], SpectralFieldIndex::Ex); - solver.ForwardTransform(*Efield_fp[lev][1], SpectralFieldIndex::Ey); - solver.ForwardTransform(*Efield_fp[lev][2], SpectralFieldIndex::Ez); - solver.ForwardTransform(*Bfield_fp[lev][0], SpectralFieldIndex::Bx); - solver.ForwardTransform(*Bfield_fp[lev][1], SpectralFieldIndex::By); - solver.ForwardTransform(*Bfield_fp[lev][2], SpectralFieldIndex::Bz); - solver.ForwardTransform(*current_fp[lev][0], SpectralFieldIndex::Jx); - solver.ForwardTransform(*current_fp[lev][1], SpectralFieldIndex::Jy); - solver.ForwardTransform(*current_fp[lev][2], SpectralFieldIndex::Jz); - solver.ForwardTransform(*rho_fp[lev], SpectralFieldIndex::rho_old, 0); - solver.ForwardTransform(*rho_fp[lev], SpectralFieldIndex::rho_new, 1); + 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_fp[lev][0], SpectralFieldIndex::Ex); - solver.BackwardTransform(*Efield_fp[lev][1], SpectralFieldIndex::Ey); - solver.BackwardTransform(*Efield_fp[lev][2], SpectralFieldIndex::Ez); - solver.BackwardTransform(*Bfield_fp[lev][0], SpectralFieldIndex::Bx); - solver.BackwardTransform(*Bfield_fp[lev][1], SpectralFieldIndex::By); - solver.BackwardTransform(*Bfield_fp[lev][2], SpectralFieldIndex::Bz); + 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_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 @@ -559,4 +574,3 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) } } } - -- cgit v1.2.3 From c5b24255c3ad5cad34823e9234f71d4a2744c471 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 31 Jul 2019 10:16:32 -0700 Subject: Add anonymous namespace --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 67 ++++++++++++++++---------------- 1 file changed, 34 insertions(+), 33 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 95e3e78b3..475f7a1dd 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -18,6 +18,40 @@ using namespace amrex; #ifdef WARPX_USE_PSATD +namespace { + void + PushPSATDSinglePatch ( + SpectralSolver& solver, + std::array,3>& Efield, + std::array,3>& Bfield, + std::array,3>& current, + std::unique_ptr& 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) @@ -34,39 +68,6 @@ WarpX::PushPSATD (amrex::Real a_dt) } } -void -PushPSATDSinglePatch ( - SpectralSolver& solver, - std::array,3>& Efield, - std::array,3>& Bfield, - std::array,3>& current, - std::unique_ptr& 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_localFFT (int lev, amrex::Real /* dt */) { -- cgit v1.2.3 From f5207b1bb5bc9dd9b0f2b7d4170bcf67fc14c8a2 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 1 Aug 2019 17:12:29 -0700 Subject: Moved ApplyInverseVolumeScalingToCurrentDensity to WarpX and is called after all particles are deposited. --- Source/Evolve/WarpXEvolveEM.cpp | 7 +++ Source/FieldSolver/WarpXPushFieldsEM.cpp | 96 +++++++++++++++++++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 8 +-- Source/WarpX.H | 7 +++ 4 files changed, 111 insertions(+), 7 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 32a4747db..57a0c44c0 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -481,6 +481,13 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time) 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]); +#ifdef WARPX_DIM_RZ + // This is called after all particles have deposited their current. + ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev); + if (current_buf[lev][0].get()) { + ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1); + } +#endif } void diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 4fce4717b..892c08aa2 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -560,3 +560,99 @@ 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& dx = WarpX::CellSize(lev); + const Real dr = dx[0]; + + Box tilebox; + + for ( MFIter mfi(*Jx, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + + Array4 const& Jr_arr = Jx->array(mfi); + Array4 const& Jt_arr = Jy->array(mfi); + Array4 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& xyzmin = WarpX::LowerCorner(tilebox, lev); + const Dim3 lo = lbound(tilebox); + const Real rmin = xyzmin[0]; + const int irmin = lo.x; + + tilebox.grow(ngJ); + + // 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) + { + /* if (j == 5) std::cout << "Jz " << i << " " << Jz_arr(i,j,0) << "\n"; */ + // 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.); + // Alternative (when order in r is limited one, which is not implemented) + // Jz_arr(i,j,0) /= (MathConst::pi*dr/4.); + } else { + Jz_arr(i,j,0) /= (2.*MathConst::pi*r); + } + }); + } +} +#endif diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 787c0e397..0e76353a4 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -504,7 +504,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow - const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev);; + const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); // xyzmin is built on pti.tilebox(), so it does // not include staggering, so the stagger_shift has to be done by hand. // Alternatively, we could define xyzminx from tbx (and the same for 3 @@ -548,12 +548,6 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, } } -#ifdef WARPX_DIM_RZ - ApplyInverseVolumeScalingToCurrentDensity(jx_arr, jy_arr, jz_arr, - tbx, tby, tbz, - xyzmin[0], lo.x, dx[0], ngJ); -#endif - #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); // CPU, tiling: atomicAdd local_jx into jx diff --git a/Source/WarpX.H b/Source/WarpX.H index a25eef9e4..dde2278d7 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -178,6 +178,13 @@ public: void EvolveE (int lev, PatchType patch_type, amrex::Real dt); void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type); +#ifdef WARPX_DIM_RZ + void ApplyInverseVolumeScalingToCurrentDensity(amrex::MultiFab* Jx, + amrex::MultiFab* Jy, + amrex::MultiFab* Jz, + int lev); +#endif + void DampPML (); void DampPML (int lev); void DampPML (int lev, PatchType patch_type); -- cgit v1.2.3 From d92404807b4e916fe49021664c425dec0a72b9ea Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 5 Aug 2019 14:59:21 -0700 Subject: Minor fixes for RZ conversion --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 5 ----- 1 file changed, 5 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 787e61f11..f8e8a639c 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -609,8 +609,6 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu const Real rmin = xyzmin[0]; const int irmin = lo.x; - tilebox.grow(ngJ); - // Rescale current in r-z mode since the inverse volume factor was not // included in the current deposition. amrex::ParallelFor(tbr, @@ -650,7 +648,6 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - /* if (j == 5) std::cout << "Jz " << i << " " << Jz_arr(i,j,0) << "\n"; */ // Wrap the current density deposited in the guard cells around // to the cells above the axis. // Jz is located on the boundary @@ -663,8 +660,6 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu if (r == 0.) { // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis Jz_arr(i,j,0) /= (MathConst::pi*dr/3.); - // Alternative (when order in r is limited one, which is not implemented) - // Jz_arr(i,j,0) /= (MathConst::pi*dr/4.); } else { Jz_arr(i,j,0) /= (2.*MathConst::pi*r); } -- cgit v1.2.3 From 9fa3a16002c585720e1741b5fdbe5e6d5b974af9 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 8 Aug 2019 16:27:54 -0700 Subject: Implemented ApplyInverseVolumeScalingToChargeDensity --- Source/Evolve/WarpXEvolveEM.cpp | 10 +++++- Source/FieldSolver/WarpXPushFieldsEM.cpp | 49 +++++++++++++++++++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 5 --- Source/WarpX.H | 4 +++ 4 files changed, 62 insertions(+), 6 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index a097d8814..e05fde03e 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -483,11 +483,19 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time) Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(), cur_time, dt[lev]); #ifdef WARPX_DIM_RZ - // This is called after all particles have deposited their current. + // 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); if (current_buf[lev][0].get()) { ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1); } + if (rho_fp[lev].get()) { + for (int icomp = 0 ; icomp < rho_fp[lev].get()->nComp() ; icomp++) { + ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), icomp, lev); + if (charge_buf[lev].get()) { + ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), icomp, lev-1); + } + } + } #endif } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 21000424d..f199e0660 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -671,4 +671,53 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu }); } } + +void +WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int icomp, int lev) +{ + const long ngRho = Rho->nGrow(); + const std::array& dx = WarpX::CellSize(lev); + const Real dr = dx[0]; + + Box tilebox; + + for ( MFIter mfi(*Rho, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + + Array4 const& Rho_arr = Rho->array(mfi); + + tilebox = mfi.tilebox(); + Box tb = convert(tilebox, IntVect::TheUnitVector()); + + // 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& xyzmin = WarpX::LowerCorner(tilebox, lev); + const Dim3 lo = lbound(tilebox); + const Real rmin = xyzmin[0]; + const int irmin = lo.x; + + // Rescale charge in r-z mode since the inverse volume factor was not + // included in the charge deposition. + amrex::ParallelFor(tb, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Wrap the charge density deposited in the guard cells around + // to the cells above the axis. + // Rho is located on the boundary + if (rmin == 0. && 0 < i && i <= ngRho) { + Rho_arr(i,j,0,icomp) += Rho_arr(-i,j,0,icomp); + } + + // 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 + Rho_arr(i,j,0,icomp) /= (MathConst::pi*dr/3.); + } else { + Rho_arr(i,j,0,icomp) /= (2.*MathConst::pi*r); + } + }); + } +} #endif diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index ec0f78564..32abae3ad 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -623,11 +623,6 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, doChargeDepositionShapeN<3>(xp, yp, zp, wp.dataPtr()+offset, rho_arr, np_to_depose, dx, xyzmin, lo, q); } -#ifdef WARPX_DIM_RZ - warpx_charge_deposition_rz_volume_scaling( - data_ptr, &ngRho, rholen.getVect(), - &xyzmin[0], &dx[0]); -#endif BL_PROFILE_VAR_STOP(blp_ppc_chd); #ifndef AMREX_USE_GPU diff --git a/Source/WarpX.H b/Source/WarpX.H index d8b1b9ce0..67d2366d4 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -183,6 +183,10 @@ public: amrex::MultiFab* Jy, amrex::MultiFab* Jz, int lev); + + void ApplyInverseVolumeScalingToChargeDensity(amrex::MultiFab* Rho, + int icomp, + int lev); #endif void DampPML (); -- cgit v1.2.3 From 006c27ae9eb5ed31fdb4824b1c7e602712d5132f Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 8 Aug 2019 16:56:27 -0700 Subject: Moved loop over ncomps into ApplyInverseVolumeScalingToChargeDensity --- Source/Evolve/WarpXEvolveEM.cpp | 8 +++----- Source/FieldSolver/WarpXPushFieldsEM.cpp | 6 +++--- Source/WarpX.H | 1 - 3 files changed, 6 insertions(+), 9 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index e05fde03e..170a8f604 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -489,11 +489,9 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time) ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1); } if (rho_fp[lev].get()) { - for (int icomp = 0 ; icomp < rho_fp[lev].get()->nComp() ; icomp++) { - ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), icomp, lev); - if (charge_buf[lev].get()) { - ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), icomp, lev-1); - } + ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), lev); + if (charge_buf[lev].get()) { + ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), lev-1); } } #endif diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index f199e0660..1df05bc0f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -673,7 +673,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu } void -WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int icomp, int lev) +WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev) { const long ngRho = Rho->nGrow(); const std::array& dx = WarpX::CellSize(lev); @@ -699,8 +699,8 @@ WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int icomp, int l // Rescale charge in r-z mode since the inverse volume factor was not // included in the charge deposition. - amrex::ParallelFor(tb, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + amrex::ParallelFor(tb, Rho->nComp(), + [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) { // Wrap the charge density deposited in the guard cells around // to the cells above the axis. diff --git a/Source/WarpX.H b/Source/WarpX.H index 67d2366d4..8d910dd4a 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -185,7 +185,6 @@ public: int lev); void ApplyInverseVolumeScalingToChargeDensity(amrex::MultiFab* Rho, - int icomp, int lev); #endif -- cgit v1.2.3