From 12ffd463b41c0819573ea2692d035914c47841ab Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Tue, 7 May 2019 15:19:09 -0700 Subject: Added code for RZ multimode for e and b advance, deposition, gather --- Source/Particles/RigidInjectedParticleContainer.cpp | 1 + 1 file changed, 1 insertion(+) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index a5acca281..77a45d146 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -427,6 +427,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), + &WarpX::nmodes, &ll4symtry, &l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); -- cgit v1.2.3 From 31edd3bad5d2747c7e045240b488d7fd101f1a30 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 12 Jul 2019 17:44:25 -0700 Subject: Changed nmodes to n_rz_azimuthal_modes --- Python/pywarpx/picmi.py | 2 +- Source/Diagnostics/FieldIO.cpp | 24 ++-- Source/Diagnostics/WarpXIO.cpp | 4 +- Source/Evolve/WarpXEvolveEM.cpp | 14 +- Source/FieldSolver/WarpXPushFieldsEM.cpp | 4 +- Source/FortranInterface/WarpX_f.H | 10 +- Source/FortranInterface/WarpX_picsar.F90 | 150 ++++++++++----------- Source/Particles/PhysicalParticleContainer.cpp | 12 +- .../Particles/RigidInjectedParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.cpp | 4 +- Source/WarpX.H | 2 +- Source/WarpX.cpp | 8 +- 12 files changed, 118 insertions(+), 118 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index b69a3ab21..67318c802 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -280,7 +280,7 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): pywarpx.geometry.is_periodic = '0 %d'%(self.bc_zmin=='periodic') # Is periodic? pywarpx.geometry.prob_lo = self.lower_bound # physical domain pywarpx.geometry.prob_hi = self.upper_bound - pywarpx.warpx.nmodes = self.n_azimuthal_modes + pywarpx.warpx.n_rz_azimuthal_modes = self.n_azimuthal_modes if self.moving_window_velocity is not None and np.any(np.not_equal(self.moving_window_velocity, 0.)): pywarpx.warpx.do_moving_window = 1 diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index c8eaa8ee0..135b7ff19 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -388,14 +388,14 @@ AverageAndPackScalarField( MultiFab& mf_avg, void AddToVarNames( Vector& varnames, std::string name, std::string suffix, - int nmodes ) { + int n_rz_azimuthal_modes ) { auto coords = {"x", "y", "z"}; auto coordsRZ = {"r", "theta", "z"}; for(auto coord:coords) varnames.push_back(name+coord+suffix); - if (nmodes > 1) { + if (n_rz_azimuthal_modes > 1) { // Note that the names are added in the same order as the fields // are packed in AverageAndPackVectorField - for (int i = 0 ; i < nmodes ; i++) { + for (int i = 0 ; i < n_rz_azimuthal_modes ; i++) { for(auto coord:coordsRZ) { varnames.push_back(name + coord + suffix + std::to_string(i) + "_real"); } @@ -415,12 +415,12 @@ WarpX::AverageAndPackFields ( Vector& varnames, amrex::Vector& mf_avg, const int ngrow) const { // Factor to account for quantities that have multiple components. - // If nmodes > 1, allow space for total field and the real and + // If n_rz_azimuthal_modes > 1, allow space for total field and the real and // imaginary part of each node. For now, also include the // imaginary part of mode 0 for code symmetry, even though // it is always zero. int modes_factor = 1; - if (nmodes > 1) modes_factor = 2*nmodes + 1; + if (n_rz_azimuthal_modes > 1) modes_factor = 2*n_rz_azimuthal_modes + 1; // Count how many different fields should be written (ncomp) const int ncomp = 0 @@ -450,17 +450,17 @@ WarpX::AverageAndPackFields ( Vector& varnames, int dcomp = 0; if (plot_J_field) { AverageAndPackVectorField(mf_avg[lev], current_fp[lev], dmap[lev], dcomp, ngrow); - if (lev == 0) AddToVarNames(varnames, "j", "", nmodes); + if (lev == 0) AddToVarNames(varnames, "j", "", n_rz_azimuthal_modes); dcomp += 3*modes_factor; } if (plot_E_field) { AverageAndPackVectorField(mf_avg[lev], Efield_aux[lev], dmap[lev], dcomp, ngrow); - if (lev == 0) AddToVarNames(varnames, "E", "", nmodes); + if (lev == 0) AddToVarNames(varnames, "E", "", n_rz_azimuthal_modes); dcomp += 3*modes_factor; } if (plot_B_field) { AverageAndPackVectorField(mf_avg[lev], Bfield_aux[lev], dmap[lev], dcomp, ngrow); - if (lev == 0) AddToVarNames(varnames, "B", "", nmodes); + if (lev == 0) AddToVarNames(varnames, "B", "", n_rz_azimuthal_modes); dcomp += 3*modes_factor; } @@ -569,10 +569,10 @@ WarpX::AverageAndPackFields ( Vector& varnames, if (plot_finepatch) { AverageAndPackVectorField( mf_avg[lev], Efield_fp[lev], dmap[lev], dcomp, ngrow ); - if (lev == 0) AddToVarNames(varnames, "E", "_fp", nmodes); + if (lev == 0) AddToVarNames(varnames, "E", "_fp", n_rz_azimuthal_modes); dcomp += 3*modes_factor; AverageAndPackVectorField( mf_avg[lev], Bfield_fp[lev], dmap[lev], dcomp, ngrow ); - if (lev == 0) AddToVarNames(varnames, "B", "_fp", nmodes); + if (lev == 0) AddToVarNames(varnames, "B", "_fp", n_rz_azimuthal_modes); dcomp += 3*modes_factor; } @@ -589,7 +589,7 @@ WarpX::AverageAndPackFields ( Vector& varnames, AverageAndPackVectorField( mf_avg[lev], E, dmap[lev], dcomp, ngrow ); } - if (lev == 0) AddToVarNames(varnames, "E", "_cp", nmodes); + if (lev == 0) AddToVarNames(varnames, "E", "_cp", n_rz_azimuthal_modes); dcomp += 3*modes_factor; // now the magnetic field @@ -603,7 +603,7 @@ WarpX::AverageAndPackFields ( Vector& varnames, std::array, 3> B = getInterpolatedB(lev); AverageAndPackVectorField( mf_avg[lev], B, dmap[lev], dcomp, ngrow ); } - if (lev == 0) AddToVarNames(varnames, "B", "_cp", nmodes); + if (lev == 0) AddToVarNames(varnames, "B", "_cp", n_rz_azimuthal_modes); dcomp += 3*modes_factor; } diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 214948b2b..bb6775680 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -416,12 +416,12 @@ WarpX::GetCellCenteredData() { Vector > cc(finest_level+1); // Factor to account for quantities that have multiple components. - // If nmodes > 1, allow space for total field and the real and + // If n_rz_azimuthal_modes > 1, allow space for total field and the real and // imaginary part of each node. For now, also include the // imaginary part of mode 0 for code symmetry, even though // it is always zero. int modes_factor = 1; - if (nmodes > 1) modes_factor = 2*nmodes + 1; + if (n_rz_azimuthal_modes > 1) modes_factor = 2*n_rz_azimuthal_modes + 1; for (int lev = 0; lev <= finest_level; ++lev) { diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index c5dbd73a7..110d07a6d 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -348,7 +348,7 @@ WarpX::OneStep_sub1 (Real curtime) RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); NodalSyncJ(fine_lev, PatchType::fine); - ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*nmodes); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*n_rz_azimuthal_modes); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -375,7 +375,7 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(coarse_lev, curtime); StoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, nmodes); + AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, n_rz_azimuthal_modes); EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); @@ -402,7 +402,7 @@ WarpX::OneStep_sub1 (Real curtime) RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); NodalSyncJ(fine_lev, PatchType::fine); - ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*nmodes); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*n_rz_azimuthal_modes); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -428,7 +428,7 @@ WarpX::OneStep_sub1 (Real curtime) // by only half a coarse step (second half) RestoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, nmodes, nmodes); + AddRhoFromFineLevelandSumBoundary(coarse_lev, n_rz_azimuthal_modes, n_rz_azimuthal_modes); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); @@ -501,12 +501,12 @@ WarpX::ComputeDt () // high enough (The simulation became unstable). Real multimode_coeffs[6] = { 0.2105, 1.0, 3.5234, 8.5104, 15.5059, 24.5037 }; Real multimode_alpha; - if (nmodes < 7) { + if (n_rz_azimuthal_modes < 7) { // Use the table of the coefficients - multimode_alpha = multimode_coeffs[nmodes-1]; + multimode_alpha = multimode_coeffs[n_rz_azimuthal_modes-1]; } else { // Use a realistic extrapolation - multimode_alpha = (nmodes - 1)*(nmodes - 1) - 0.4; + multimode_alpha = (n_rz_azimuthal_modes - 1)*(n_rz_azimuthal_modes - 1) - 0.4; } deltat = cfl * 1./( std::sqrt((1+multimode_alpha)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); #else diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index fc4fb902b..f08b44adb 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -109,7 +109,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) tbx.loVect(), tbx.hiVect(), tby.loVect(), tby.hiVect(), tbz.loVect(), tbz.hiVect(), - &nmodes, + &n_rz_azimuthal_modes, BL_TO_FORTRAN_3D((*Ex)[mfi]), BL_TO_FORTRAN_3D((*Ey)[mfi]), BL_TO_FORTRAN_3D((*Ez)[mfi]), @@ -272,7 +272,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.hiVect(), - &nmodes, + &n_rz_azimuthal_modes, BL_TO_FORTRAN_3D((*Ex)[mfi]), BL_TO_FORTRAN_3D((*Ey)[mfi]), BL_TO_FORTRAN_3D((*Ez)[mfi]), diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 8dcd4222b..f193f30c4 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -102,7 +102,7 @@ extern "C" amrex::Real* jx, const long* jx_ng, const int* jx_ntot, amrex::Real* jy, const long* jy_ng, const int* jy_ntot, amrex::Real* jz, const long* jz_ng, const int* jz_ntot, - const long* nmodes, + const long* n_rz_azimuthal_modes, const long* np, const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp, const amrex::Real* uxp, const amrex::Real* uyp,const amrex::Real* uzp, @@ -118,7 +118,7 @@ extern "C" amrex::Real* jx, const long* jx_ng, const int* jx_ntot, amrex::Real* jy, const long* jy_ng, const int* jy_ntot, amrex::Real* jz, const long* jz_ng, const int* jz_ntot, - const long* nmodes, + const long* n_rz_azimuthal_modes, const amrex::Real* rmin, const amrex::Real* dr); @@ -138,7 +138,7 @@ extern "C" const amrex::Real* bxg, const int* bxg_lo, const int* bxg_hi, const amrex::Real* byg, const int* byg_lo, const int* byg_hi, const amrex::Real* bzg, const int* bzg_lo, const int* bzg_hi, - const long* nmodes, + const long* n_rz_azimuthal_modes, const int* ll4symtry, const int* l_lower_order_in_v, const int* l_nodal, const long* lvect, const long* field_gathe_algo); @@ -197,7 +197,7 @@ extern "C" const int* xlo, const int* xhi, const int* ylo, const int* yhi, const int* zlo, const int* zhi, - const long* nmodes, + const long* n_rz_azimuthal_modes, BL_FORT_FAB_ARG_3D(ex), BL_FORT_FAB_ARG_3D(ey), BL_FORT_FAB_ARG_3D(ez), @@ -218,7 +218,7 @@ extern "C" const int* xlo, const int* xhi, const int* ylo, const int* yhi, const int* zlo, const int* zhi, - const long* nmodes, + const long* n_rz_azimuthal_modes, const BL_FORT_FAB_ARG_3D(ex), const BL_FORT_FAB_ARG_3D(ey), const BL_FORT_FAB_ARG_3D(ez), diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index e6e61b734..c94d8459f 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -87,7 +87,7 @@ contains ex,ey,ez,bx,by,bz,ixyzmin,xmin,ymin,zmin,dx,dy,dz,nox,noy,noz, & exg,exg_lo,exg_hi,eyg,eyg_lo,eyg_hi,ezg,ezg_lo,ezg_hi, & bxg,bxg_lo,bxg_hi,byg,byg_lo,byg_hi,bzg,bzg_lo,bzg_hi, & - nmodes, & + n_rz_azimuthal_modes, & ll4symtry,l_lower_order_in_v, l_nodal,& lvect,field_gathe_algo) & bind(C, name="warpx_geteb_energy_conserving") @@ -99,19 +99,19 @@ contains integer, intent(in) :: ixyzmin(AMREX_SPACEDIM) real(amrex_real), intent(in) :: xmin,ymin,zmin,dx,dy,dz integer(c_long), intent(in) :: field_gathe_algo - integer(c_long), intent(in) :: np,nmodes,nox,noy,noz + integer(c_long), intent(in) :: np,n_rz_azimuthal_modes,nox,noy,noz integer(c_int), intent(in) :: ll4symtry,l_lower_order_in_v, l_nodal integer(c_long),intent(in) :: lvect real(amrex_real), intent(in), dimension(np) :: xp,yp,zp real(amrex_real), intent(out), dimension(np) :: ex,ey,ez,bx,by,bz #ifdef WARPX_RZ ! The dimensions must be specified to allow the transpose - real(amrex_real),intent(in):: exg(exg_lo(1):exg_hi(1),exg_lo(2):exg_hi(2),2,nmodes) - real(amrex_real),intent(in):: eyg(eyg_lo(1):eyg_hi(1),eyg_lo(2):eyg_hi(2),2,nmodes) - real(amrex_real),intent(in):: ezg(ezg_lo(1):ezg_hi(1),ezg_lo(2):ezg_hi(2),2,nmodes) - real(amrex_real),intent(in):: bxg(bxg_lo(1):bxg_hi(1),bxg_lo(2):bxg_hi(2),2,nmodes) - real(amrex_real),intent(in):: byg(byg_lo(1):byg_hi(1),byg_lo(2):byg_hi(2),2,nmodes) - real(amrex_real),intent(in):: bzg(bzg_lo(1):bzg_hi(1),bzg_lo(2):bzg_hi(2),2,nmodes) + real(amrex_real),intent(in):: exg(exg_lo(1):exg_hi(1),exg_lo(2):exg_hi(2),2,n_rz_azimuthal_modes) + real(amrex_real),intent(in):: eyg(eyg_lo(1):eyg_hi(1),eyg_lo(2):eyg_hi(2),2,n_rz_azimuthal_modes) + real(amrex_real),intent(in):: ezg(ezg_lo(1):ezg_hi(1),ezg_lo(2):ezg_hi(2),2,n_rz_azimuthal_modes) + real(amrex_real),intent(in):: bxg(bxg_lo(1):bxg_hi(1),bxg_lo(2):bxg_hi(2),2,n_rz_azimuthal_modes) + real(amrex_real),intent(in):: byg(byg_lo(1):byg_hi(1),byg_lo(2):byg_hi(2),2,n_rz_azimuthal_modes) + real(amrex_real),intent(in):: bzg(bzg_lo(1):bzg_hi(1),bzg_lo(2):bzg_hi(2),2,n_rz_azimuthal_modes) #else ! These could be either 2d or 3d arrays real(amrex_real),intent(in):: exg(*), eyg(*), ezg(*), bxg(*), byg(*), bzg(*) @@ -148,14 +148,14 @@ contains bzg_nvalid = bzg_lo + bzg_hi - 2_c_long*ixyzmin + 1_c_long #ifdef WARPX_RZ - if (nmodes > 1) then - - allocate(erg_c(exg_lo(1):exg_hi(1),exg_lo(2):exg_hi(2),nmodes), & - etg_c(eyg_lo(1):eyg_hi(1),eyg_lo(2):eyg_hi(2),nmodes), & - ezg_c(ezg_lo(1):ezg_hi(1),ezg_lo(2):ezg_hi(2),nmodes), & - brg_c(bxg_lo(1):bxg_hi(1),bxg_lo(2):bxg_hi(2),nmodes), & - btg_c(byg_lo(1):byg_hi(1),byg_lo(2):byg_hi(2),nmodes), & - bzg_c(bzg_lo(1):bzg_hi(1),bzg_lo(2):bzg_hi(2),nmodes), stat=alloc_status) + if (n_rz_azimuthal_modes > 1) then + + allocate(erg_c(exg_lo(1):exg_hi(1),exg_lo(2):exg_hi(2),n_rz_azimuthal_modes), & + etg_c(eyg_lo(1):eyg_hi(1),eyg_lo(2):eyg_hi(2),n_rz_azimuthal_modes), & + ezg_c(ezg_lo(1):ezg_hi(1),ezg_lo(2):ezg_hi(2),n_rz_azimuthal_modes), & + brg_c(bxg_lo(1):bxg_hi(1),bxg_lo(2):bxg_hi(2),n_rz_azimuthal_modes), & + btg_c(byg_lo(1):byg_hi(1),byg_lo(2):byg_hi(2),n_rz_azimuthal_modes), & + bzg_c(bzg_lo(1):bzg_hi(1),bzg_lo(2):bzg_hi(2),n_rz_azimuthal_modes), stat=alloc_status) if (alloc_status /= 0) then print*,"Error: warpx_geteb_energy_conserving: complex arrays could not be allocated" stop @@ -173,7 +173,7 @@ contains bzg_c(:,:,:) = cmplx(bzg(:,:,1,:), bzg(:,:,2,:), amrex_real) call geteb2dcirc_energy_conserving_generic(np, xp, yp, zp, ex, ey, ez, bx, by, bz, & - xmin, zmin, dx, dz, nmodes, nox, noz, & + xmin, zmin, dx, dz, n_rz_azimuthal_modes, nox, noz, & pxr_l_lower_order_in_v, pxr_l_nodal, & erg_c, exg_nguards, exg_nvalid, & etg_c, eyg_nguards, eyg_nvalid, & @@ -372,22 +372,22 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> @param[in] charge_depo_algo algorithm choice for the charge deposition !> subroutine warpx_current_deposition( & - jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot,nmodes, & + jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot,n_rz_azimuthal_modes, & np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin,dt,dx,dy,dz,nox,noy,noz,& l_nodal,lvect,current_depo_algo) & bind(C, name="warpx_current_deposition") integer, intent(in) :: jx_ntot(AMREX_SPACEDIM), jy_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM) integer(c_long), intent(in) :: jx_ng, jy_ng, jz_ng - integer(c_long), intent(IN) :: nmodes + integer(c_long), intent(IN) :: n_rz_azimuthal_modes integer(c_long), intent(IN) :: np integer(c_long), intent(IN) :: nox,noy,noz integer(c_int), intent(in) :: l_nodal #ifdef WARPX_RZ - real(amrex_real), intent(IN OUT):: jx(jx_ntot(1),jx_ntot(2),2,nmodes) - real(amrex_real), intent(IN OUT):: jy(jy_ntot(1),jy_ntot(2),2,nmodes) - real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2,nmodes) + real(amrex_real), intent(IN OUT):: jx(jx_ntot(1),jx_ntot(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: jy(jy_ntot(1),jy_ntot(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2,n_rz_azimuthal_modes) #else real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) #endif @@ -432,11 +432,11 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! Dimension 2 #elif (AMREX_SPACEDIM==2) #ifdef WARPX_RZ - if (nmodes > 1) then + if (n_rz_azimuthal_modes > 1) then - allocate(jr_c(jx_ntot(1),jx_ntot(2),nmodes), & - jt_c(jy_ntot(1),jy_ntot(2),nmodes), & - jz_c(jz_ntot(1),jz_ntot(2),nmodes), stat=alloc_status) + allocate(jr_c(jx_ntot(1),jx_ntot(2),n_rz_azimuthal_modes), & + jt_c(jy_ntot(1),jy_ntot(2),n_rz_azimuthal_modes), & + jz_c(jz_ntot(1),jz_ntot(2),n_rz_azimuthal_modes), stat=alloc_status) if (alloc_status /= 0) then print*,"Error: warpx_current_deposition: complex arrays could not be allocated" stop @@ -450,7 +450,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jr_c,jx_nguards,jx_nvalid, & jt_c,jy_nguards,jy_nvalid, & jz_c,jz_nguards,jz_nvalid, & - nmodes, & + n_rz_azimuthal_modes, & np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, & xmin,zmin,dt,dx,dz, & nox,noz,l_particles_weight,type_rz_depose) @@ -500,15 +500,15 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> subroutine warpx_current_deposition_rz_volume_scaling( & jr,jr_ng,jr_ntot,jt,jt_ng,jt_ntot,jz,jz_ng,jz_ntot, & - nmodes,rmin,dr) & + n_rz_azimuthal_modes,rmin,dr) & bind(C, name="warpx_current_deposition_rz_volume_scaling") integer, intent(in) :: jr_ntot(AMREX_SPACEDIM), jt_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM) integer(c_long), intent(in) :: jr_ng, jt_ng, jz_ng - integer(c_long), intent(in) :: nmodes - real(amrex_real), intent(IN OUT):: jr(jr_ntot(1),jr_ntot(2),2,nmodes) - real(amrex_real), intent(IN OUT):: jt(jt_ntot(1),jt_ntot(2),2,nmodes) - real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2,nmodes) + integer(c_long), intent(in) :: n_rz_azimuthal_modes + real(amrex_real), intent(IN OUT):: jr(jr_ntot(1),jr_ntot(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: jt(jt_ntot(1),jt_ntot(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2,n_rz_azimuthal_modes) real(amrex_real), intent(IN) :: rmin, dr #ifdef WARPX_RZ @@ -528,11 +528,11 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jt_nguards = jt_ng jz_nguards = jz_ng - if (nmodes > 1) then + if (n_rz_azimuthal_modes > 1) then - allocate(jr_c(jr_ntot(1),jr_ntot(2),nmodes), & - jt_c(jt_ntot(1),jt_ntot(2),nmodes), & - jz_c(jz_ntot(1),jz_ntot(2),nmodes), stat=alloc_status) + allocate(jr_c(jr_ntot(1),jr_ntot(2),n_rz_azimuthal_modes), & + jt_c(jt_ntot(1),jt_ntot(2),n_rz_azimuthal_modes), & + jz_c(jz_ntot(1),jz_ntot(2),n_rz_azimuthal_modes), stat=alloc_status) if (alloc_status /= 0) then print*,"Error: warpx_current_deposition_rz_volume_scaling: complex arrays could not be allocated" stop @@ -550,7 +550,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jr_c, jr_nguards, jr_nvalid, & jt_c, jt_nguards, jt_nvalid, & jz_c, jz_nguards, jz_nvalid, & - nmodes, & + n_rz_azimuthal_modes, & rmin,dr, & type_rz_depose) @@ -710,7 +710,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> @param[in] dtsdx, dtsdy, dtsdz factors c**2 * dt/(dx, dy, dz) subroutine warpx_push_evec( & xlo, xhi, ylo, yhi, zlo, zhi, & - nmodes, & + n_rz_azimuthal_modes, & ex, exlo, exhi, & ey, eylo, eyhi, & ez, ezlo, ezhi, & @@ -730,19 +730,19 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jxlo(BL_SPACEDIM), jxhi(BL_SPACEDIM), jylo(BL_SPACEDIM), jyhi(BL_SPACEDIM), & jzlo(BL_SPACEDIM), jzhi(BL_SPACEDIM) - integer(c_long), intent(in) :: nmodes + integer(c_long), intent(in) :: n_rz_azimuthal_modes #ifdef WARPX_RZ ! The dimensions must be specified to allow the transpose - real(amrex_real), intent(IN OUT):: ex(exlo(1):exhi(1),exlo(2):exhi(2),2,nmodes) - real(amrex_real), intent(IN OUT):: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),2,nmodes) - real(amrex_real), intent(IN OUT):: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),2,nmodes) - real(amrex_real), intent(IN):: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),2,nmodes) - real(amrex_real), intent(IN):: by(bylo(1):byhi(1),bylo(2):byhi(2),2,nmodes) - real(amrex_real), intent(IN):: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2,nmodes) - real(amrex_real), intent(IN):: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2),2,nmodes) - real(amrex_real), intent(IN):: jy(jylo(1):jyhi(1),jylo(2):jyhi(2),2,nmodes) - real(amrex_real), intent(IN):: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2),2,nmodes) + real(amrex_real), intent(IN OUT):: ex(exlo(1):exhi(1),exlo(2):exhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: by(bylo(1):byhi(1),bylo(2):byhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: jy(jylo(1):jyhi(1),jylo(2):jyhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2),2,n_rz_azimuthal_modes) #else ! These can be either 2d or 3d real(amrex_real), intent(IN OUT):: ex(*) @@ -766,7 +766,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n complex(amrex_real), allocatable, dimension(:,:,:) :: jr_c, jt_c, jz_c integer :: alloc_status - if (nmodes == 1) then + if (n_rz_azimuthal_modes == 1) then #endif CALL WRPX_PXR_PUSH_EVEC(& xlo, xhi, ylo, yhi, zlo, zhi, & @@ -788,15 +788,15 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n #ifdef WARPX_RZ else - allocate(er_c(exlo(1):exhi(1),exlo(2):exhi(2),nmodes), & - et_c(eylo(1):eyhi(1),eylo(2):eyhi(2),nmodes), & - ez_c(ezlo(1):ezhi(1),ezlo(2):ezhi(2),nmodes), & - br_c(bxlo(1):bxhi(1),bxlo(2):bxhi(2),nmodes), & - bt_c(bylo(1):byhi(1),bylo(2):byhi(2),nmodes), & - bz_c(bzlo(1):bzhi(1),bzlo(2):bzhi(2),nmodes), & - jr_c(jxlo(1):jxhi(1),jxlo(2):jxhi(2),nmodes), & - jt_c(jylo(1):jyhi(1),jylo(2):jyhi(2),nmodes), & - jz_c(jzlo(1):jzhi(1),jzlo(2):jzhi(2),nmodes), stat=alloc_status) + allocate(er_c(exlo(1):exhi(1),exlo(2):exhi(2),n_rz_azimuthal_modes), & + et_c(eylo(1):eyhi(1),eylo(2):eyhi(2),n_rz_azimuthal_modes), & + ez_c(ezlo(1):ezhi(1),ezlo(2):ezhi(2),n_rz_azimuthal_modes), & + br_c(bxlo(1):bxhi(1),bxlo(2):bxhi(2),n_rz_azimuthal_modes), & + bt_c(bylo(1):byhi(1),bylo(2):byhi(2),n_rz_azimuthal_modes), & + bz_c(bzlo(1):bzhi(1),bzlo(2):bzhi(2),n_rz_azimuthal_modes), & + jr_c(jxlo(1):jxhi(1),jxlo(2):jxhi(2),n_rz_azimuthal_modes), & + jt_c(jylo(1):jyhi(1),jylo(2):jyhi(2),n_rz_azimuthal_modes), & + jz_c(jzlo(1):jzhi(1),jzlo(2):jzhi(2),n_rz_azimuthal_modes), stat=alloc_status) if (alloc_status /= 0) then print*,"Error: warpx_push_evec: complex arrays could not be allocated" stop @@ -818,7 +818,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n CALL pxrpush_emrz_evec_multimode(& xlo, xhi, ylo, yhi, zlo, zhi, & - nmodes, & + n_rz_azimuthal_modes, & er_c, exlo, exhi,& et_c, eylo, eyhi,& ez_c, ezlo, ezhi,& @@ -863,7 +863,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> @param[in] dtsdx, dtsdy, dtsdz factors 0.5 * dt/(dx, dy, dz) subroutine warpx_push_bvec( & xlo, xhi, ylo, yhi, zlo, zhi, & - nmodes, & + n_rz_azimuthal_modes, & ex, exlo, exhi, & ey, eylo, eyhi, & ez, ezlo, ezhi, & @@ -880,16 +880,16 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n bylo(BL_SPACEDIM), byhi(BL_SPACEDIM), bzlo(BL_SPACEDIM), bzhi(BL_SPACEDIM), & maxwell_fdtd_solver_id - integer(c_long), intent(in) :: nmodes + integer(c_long), intent(in) :: n_rz_azimuthal_modes #ifdef WARPX_RZ ! The dimensions must be specified to allow the transpose - real(amrex_real), intent(IN):: ex(exlo(1):exhi(1),exlo(2):exhi(2),2,nmodes) - real(amrex_real), intent(IN):: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),2,nmodes) - real(amrex_real), intent(IN):: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),2,nmodes) - real(amrex_real), intent(IN OUT):: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),2,nmodes) - real(amrex_real), intent(IN OUT):: by(bylo(1):byhi(1),bylo(2):byhi(2),2,nmodes) - real(amrex_real), intent(IN OUT):: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2,nmodes) + real(amrex_real), intent(IN):: ex(exlo(1):exhi(1),exlo(2):exhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN):: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: by(bylo(1):byhi(1),bylo(2):byhi(2),2,n_rz_azimuthal_modes) + real(amrex_real), intent(IN OUT):: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2,n_rz_azimuthal_modes) #else real(amrex_real), intent(IN):: ex(*) real(amrex_real), intent(IN):: ey(*) @@ -909,7 +909,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n integer :: alloc_status #endif - if (nmodes == 1) then + if (n_rz_azimuthal_modes == 1) then IF (maxwell_fdtd_solver_id .eq. 0) THEN ! Yee FDTD solver @@ -942,12 +942,12 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n #ifdef WARPX_RZ else - allocate(er_c(exlo(1):exhi(1),exlo(2):exhi(2),nmodes), & - et_c(eylo(1):eyhi(1),eylo(2):eyhi(2),nmodes), & - ez_c(ezlo(1):ezhi(1),ezlo(2):ezhi(2),nmodes), & - br_c(bxlo(1):bxhi(1),bxlo(2):bxhi(2),nmodes), & - bt_c(bylo(1):byhi(1),bylo(2):byhi(2),nmodes), & - bz_c(bzlo(1):bzhi(1),bzlo(2):bzhi(2),nmodes), stat=alloc_status) + allocate(er_c(exlo(1):exhi(1),exlo(2):exhi(2),n_rz_azimuthal_modes), & + et_c(eylo(1):eyhi(1),eylo(2):eyhi(2),n_rz_azimuthal_modes), & + ez_c(ezlo(1):ezhi(1),ezlo(2):ezhi(2),n_rz_azimuthal_modes), & + br_c(bxlo(1):bxhi(1),bxlo(2):bxhi(2),n_rz_azimuthal_modes), & + bt_c(bylo(1):byhi(1),bylo(2):byhi(2),n_rz_azimuthal_modes), & + bz_c(bzlo(1):bzhi(1),bzlo(2):bzhi(2),n_rz_azimuthal_modes), stat=alloc_status) if (alloc_status /= 0) then print*,"Error: warpx_push_bvec: complex arrays could not be allocated" stop @@ -966,7 +966,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n CALL pxrpush_emrz_bvec_multimode(& xlo, xhi, ylo, yhi, zlo, zhi, & - nmodes, & + n_rz_azimuthal_modes, & er_c, exlo, exhi,& et_c, eylo, eyhi,& ez_c, ezlo, ezhi,& diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 93a0ad7ea..aa9fbada7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -482,7 +482,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) // Replace the x and y, setting an angle theta. // These x and y are used to get the momentum and density Real theta; - if (WarpX::nmodes == 1) { + if (WarpX::n_rz_azimuthal_modes == 1) { // With only 1 mode, the angle doesn't matter so // choose it randomly. theta = 2.*MathConst::pi*amrex::Random(); @@ -734,7 +734,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) // Replace the x and y, setting an angle theta. // These x and y are used to get the momentum and density Real theta; - if (WarpX::nmodes == 1) { + if (WarpX::n_rz_azimuthal_modes == 1) { // With only 1 mode, the angle doesn't matter so // choose it randomly. theta = 2.*MathConst::pi*amrex::Random(); @@ -1156,7 +1156,7 @@ PhysicalParticleContainer::FieldGather (int lev, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), - &WarpX::nmodes, + &WarpX::n_rz_azimuthal_modes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); @@ -1449,7 +1449,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*bxfab), BL_TO_FORTRAN_ANYD(*byfab), BL_TO_FORTRAN_ANYD(*bzfab), - &WarpX::nmodes, + &WarpX::n_rz_azimuthal_modes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); @@ -1537,7 +1537,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*cbxfab), BL_TO_FORTRAN_ANYD(*cbyfab), BL_TO_FORTRAN_ANYD(*cbzfab), - &WarpX::nmodes, + &WarpX::n_rz_azimuthal_modes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); } @@ -1879,7 +1879,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), - &WarpX::nmodes, + &WarpX::n_rz_azimuthal_modes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index d659b3854..ff07527aa 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -426,7 +426,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), - &WarpX::nmodes, + &WarpX::n_rz_azimuthal_modes, &ll4symtry, &l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index a1517cb44..6ecf0fe98 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -379,7 +379,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), - &WarpX::nmodes, + &WarpX::n_rz_azimuthal_modes, &np_to_depose, m_xp[thread_num].dataPtr() + offset, m_yp[thread_num].dataPtr() + offset, @@ -400,7 +400,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), - &WarpX::nmodes, + &WarpX::n_rz_azimuthal_modes, &xyzmin[0], &dx[0]); #endif BL_PROFILE_VAR_STOP(blp_pxr_cd); diff --git a/Source/WarpX.H b/Source/WarpX.H index 5d3008972..581ca39ce 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -93,7 +93,7 @@ public: static long noz; // Number of modes for the RZ multimode version - static long nmodes; + static long n_rz_azimuthal_modes; static long ncomps; static bool use_fdtd_nci_corr; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 4ad2cf75d..a693952dd 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -43,7 +43,7 @@ long WarpX::field_gathering_algo; long WarpX::particle_pusher_algo; int WarpX::maxwell_fdtd_solver_id; -long WarpX::nmodes = 1; +long WarpX::n_rz_azimuthal_modes = 1; long WarpX::ncomps = 1; long WarpX::nox = 1; @@ -475,7 +475,7 @@ WarpX::ReadParameters () } // Only needs to be set with WARPX_RZ, otherwise defaults to 1. - pp.query("nmodes", nmodes); + pp.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes); } @@ -709,9 +709,9 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm { #if defined WARPX_RZ - if (nmodes > 1) { + if (n_rz_azimuthal_modes > 1) { // There is a real and imaginary component for each mode - ncomps = nmodes*2; + ncomps = n_rz_azimuthal_modes*2; } else { // With only mode 0, only reals are used ncomps = 1; -- cgit v1.2.3 From 69be3657b7410db305e102e89c366ebbc8be59d1 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 16 Aug 2019 16:37:23 -0700 Subject: Converted field gather for RZ multimode to C++ --- Source/Particles/Gather/FieldGather.H | 80 ++++++++++++++++++---- Source/Particles/PhysicalParticleContainer.H | 1 - Source/Particles/PhysicalParticleContainer.cpp | 20 +++--- .../Particles/RigidInjectedParticleContainer.cpp | 2 +- 4 files changed, 76 insertions(+), 27 deletions(-) (limited to 'Source/Particles/RigidInjectedParticleContainer.cpp') diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 8f5e8d4cf..219a9d055 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -2,9 +2,10 @@ #define FIELDGATHER_H_ #include "ShapeFactors.H" +#include /* \brief Field gather for particles handled by thread thread_num - * /param xp, yp, zp : Pointer to arrays of particle positions. + * \param xp, yp, zp : Pointer to arrays of particle positions. * \param Exp, Eyp, Ezp: Pointer to array of electric field on particles. * \param Bxp, Byp, Bzp: Pointer to array of magnetic field on particles. * \param ex_arr ey_arr: Array4 of current density, either full array or tile. @@ -15,6 +16,7 @@ * \param xyzmin : Physical lower bounds of domain. * \param lo : Index lower bounds of domain. * \param stagger_shift: 0 if nodal, 0.5 if staggered. + * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry */ template void doGatherShapeN(const amrex::Real * const xp, @@ -33,7 +35,8 @@ void doGatherShapeN(const amrex::Real * const xp, const std::array& dx, const std::array xyzmin, const amrex::Dim3 lo, - const amrex::Real stagger_shift) + const amrex::Real stagger_shift, + const long n_rz_azimuthal_modes) { const amrex::Real dxi = 1.0/dx[0]; const amrex::Real dzi = 1.0/dx[2]; @@ -56,8 +59,8 @@ void doGatherShapeN(const amrex::Real * const xp, // x direction // Get particle position #ifdef WARPX_DIM_RZ - const amrex::Real r = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); - const amrex::Real x = (r - xmin)*dxi; + const amrex::Real rp = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + const amrex::Real x = (rp - xmin)*dxi; #else const amrex::Real x = (xp[ip]-xmin)*dxi; #endif @@ -103,7 +106,7 @@ void doGatherShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ Eyp[ip] += sx[ix]*sz[iz]* - ey_arr(lo.x+j+ix, lo.y+l+iz, 0); + ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 0); } } // Gather field on particle Exp[i] from field on grid ex_arr @@ -111,9 +114,9 @@ void doGatherShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ Exp[ip] += sx0[ix]*sz[iz]* - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0); + ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 0); Bzp[ip] += sx0[ix]*sz[iz]* - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0); + bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 0); } } // Gather field on particle Ezp[i] from field on grid ez_arr @@ -121,30 +124,79 @@ void doGatherShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order; ix++){ Ezp[ip] += sx[ix]*sz0[iz]* - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0); + ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 0); Bxp[ip] += sx[ix]*sz0[iz]* - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0); + bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 0); } } // Gather field on particle Byp[i] from field on grid by_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ Byp[ip] += sx0[ix]*sz0[iz]* - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0); + by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 0); } } #ifdef WARPX_DIM_RZ - // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey + amrex::Real costheta; amrex::Real sintheta; - if (r > 0.) { - costheta = xp[ip]/r; - sintheta = yp[ip]/r; + if (rp > 0.) { + costheta = xp[ip]/rp; + sintheta = yp[ip]/rp; } else { costheta = 1.; sintheta = 0.; } + const Complex xy0 = Complex{costheta, -sintheta}; + Complex xy = xy0; + + for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { + + // Gather field on particle Eyp[i] from field on grid ey_arr + for (int iz=0; iz<=depos_order; iz++){ + for (int ix=0; ix<=depos_order; ix++){ + const amrex::Real dEy = (+ ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode)*xy.real() + - ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode+1)*xy.imag()); + Eyp[ip] += sx[ix]*sz[iz]*dEy; + } + } + // Gather field on particle Exp[i] from field on grid ex_arr + // Gather field on particle Bzp[i] from field on grid bz_arr + for (int iz=0; iz<=depos_order; iz++){ + for (int ix=0; ix<=depos_order-lower_in_v; ix++){ + const amrex::Real dEx = (+ ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.real() + - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode+1)*xy.imag()); + Exp[ip] += sx0[ix]*sz[iz]*dEx; + const amrex::Real dBz = (+ bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.real() + - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode+1)*xy.imag()); + Bzp[ip] += sx0[ix]*sz[iz]*dBz; + } + } + // Gather field on particle Ezp[i] from field on grid ez_arr + // Gather field on particle Bxp[i] from field on grid bx_arr + for (int iz=0; iz<=depos_order-lower_in_v; iz++){ + for (int ix=0; ix<=depos_order; ix++){ + const amrex::Real dEz = (+ ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.real() + - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode+1)*xy.imag()); + Ezp[ip] += sx[ix]*sz0[iz]*dEz; + const amrex::Real dBx = (+ bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.real() + - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode+1)*xy.imag()); + Bxp[ip] += sx[ix]*sz0[iz]*dBx; + } + } + // Gather field on particle Byp[i] from field on grid by_arr + for (int iz=0; iz<=depos_order-lower_in_v; iz++){ + for (int ix=0; ix<=depos_order-lower_in_v; ix++){ + const amrex::Real dBy = (+ by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode)*xy.real() + - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode+1)*xy.imag()); + Byp[ip] += sx0[ix]*sz0[iz]*dBy; + } + } + xy = xy*xy0; + } + + // Convert Exp and Eyp (which are actually Er and Etheta) to Ex and Ey const amrex::Real Exp_save = Exp[ip]; Exp[ip] = costheta*Exp[ip] - sintheta*Eyp[ip]; Eyp[ip] = costheta*Eyp[ip] + sintheta*Exp_save; diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ebd4eac2b..b80619733 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -53,7 +53,6 @@ public: amrex::FArrayBox const * byfab, amrex::FArrayBox const * bzfab, const int ngE, const int e_is_nodal, - const long n_rz_azimuthal_modes, const long offset, const long np_to_gather, int thread_num, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c44a9a8b7..ec360f2c4 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -899,7 +899,7 @@ PhysicalParticleContainer::FieldGather (int lev, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, - Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes, + Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); if (cost) { @@ -1176,7 +1176,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_START(blp_pxr_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, exfab, eyfab, ezfab, bxfab, byfab, bzfab, - Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes, + Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); if (np_gather < np) @@ -1250,7 +1250,6 @@ PhysicalParticleContainer::Evolve (int lev, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, cEx->nGrow(), e_is_nodal, - WarpX::n_rz_azimuthal_modes, nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1594,7 +1593,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, - Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes, + Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); // This wraps the momentum advance so that inheritors can modify the call. @@ -1835,7 +1834,6 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, - const long n_rz_azimuthal_modes, int thread_num, int lev, int gather_lev) @@ -1890,7 +1888,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doGatherShapeN<2,1>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1898,7 +1896,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doGatherShapeN<3,1>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1906,7 +1904,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } } else { if (WarpX::nox == 1){ @@ -1916,7 +1914,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doGatherShapeN<2,0>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1924,7 +1922,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doGatherShapeN<3,0>(xp, yp, zp, Exp.dataPtr() + offset, Eyp.dataPtr() + offset, @@ -1932,7 +1930,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, Byp.dataPtr() + offset, Bzp.dataPtr() + offset, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, np_to_gather, dx, - xyzmin, lo, stagger_shift); + xyzmin, lo, stagger_shift, WarpX::n_rz_azimuthal_modes); } } } diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 695955e15..a87168a75 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -426,7 +426,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, - Ex.nGrow(), e_is_nodal, WarpX::n_rz_azimuthal_modes, + Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); // Save the position and momenta, making copies -- cgit v1.2.3