diff options
Diffstat (limited to 'Source/Particles')
-rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 2 | ||||
-rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 81 | ||||
-rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 1 | ||||
-rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 3 | ||||
-rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 30 |
5 files changed, 75 insertions, 42 deletions
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 9d39ec2f9..6475d1463 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -297,7 +297,7 @@ MultiParticleContainer::GetChargeDensity (int lev, bool local) std::unique_ptr<MultiFab> rho = allcontainers[0]->GetChargeDensity(lev, true); for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { std::unique_ptr<MultiFab> rhoi = allcontainers[i]->GetChargeDensity(lev, true); - MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow()); + MultiFab::Add(*rho, *rhoi, 0, 0, rho->nComp(), rho->nGrow()); } if (!local) { const Geometry& gm = allcontainers[0]->Geom(lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 43b46ec49..93a0ad7ea 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -41,16 +41,19 @@ NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac); for (int i_part=0; i_part<ref_num_ppc;i_part++) { - std::array<Real, 3> r; - plasma_injector->getPositionUnitBox(r, i_part, fac); + std::array<Real, 3> point; + plasma_injector->getPositionUnitBox(point, i_part, fac); + Real x = overlap_corner[0] + (iv[0] + point[0])*dx[0]; #if ( AMREX_SPACEDIM == 3 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; - Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1]; - Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2]; + Real y = overlap_corner[1] + (iv[1] + point[1])*dx[1]; + Real z = overlap_corner[2] + (iv[2] + point[2])*dx[2]; #elif ( AMREX_SPACEDIM == 2 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; Real y = 0; - Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1]; +#ifdef WARPX_RZ + Real z = overlap_corner[1] + (iv[1] + point[2])*dx[1]; +#else + Real z = overlap_corner[1] + (iv[1] + point[1])*dx[1]; +#endif #endif // If the new particle is not inside the tile box, // go to the next generated particle. @@ -448,16 +451,19 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac); for (int i_part=0; i_part<ref_num_ppc;i_part++) { - std::array<Real, 3> r; - plasma_injector->getPositionUnitBox(r, i_part, fac); + std::array<Real, 3> point; + plasma_injector->getPositionUnitBox(point, i_part, fac); + Real x = overlap_corner[0] + (iv[0] + point[0])*dx[0]; #if ( AMREX_SPACEDIM == 3 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; - Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1]; - Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2]; + Real y = overlap_corner[1] + (iv[1] + point[1])*dx[1]; + Real z = overlap_corner[2] + (iv[2] + point[2])*dx[2]; #elif ( AMREX_SPACEDIM == 2 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; Real y = 0; - Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1]; +#ifdef WARPX_RZ + Real z = overlap_corner[1] + (iv[1] + point[2])*dx[1]; +#else + Real z = overlap_corner[1] + (iv[1] + point[1])*dx[1]; +#endif #endif // If the new particle is not inside the tile box, // go to the next generated particle. @@ -473,11 +479,18 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) Real yb = y; #ifdef WARPX_RZ - // Replace the x and y, choosing the angle randomly. + // Replace the x and y, setting an angle theta. // These x and y are used to get the momentum and density - Real theta = 2.*MathConst::pi*amrex::Random(); - y = x*std::sin(theta); - x = x*std::cos(theta); + Real theta; + if (WarpX::nmodes == 1) { + // With only 1 mode, the angle doesn't matter so + // choose it randomly. + theta = 2.*MathConst::pi*amrex::Random(); + } else { + theta = 2.*MathConst::pi*point[1]; + } + x = xb*std::cos(theta); + y = xb*std::sin(theta); #endif Real dens; @@ -690,16 +703,19 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac); for (int i_part=0; i_part<ref_num_ppc;i_part++) { - std::array<Real, 3> r; - plasma_injector->getPositionUnitBox(r, i_part, fac); + std::array<Real, 3> point; + plasma_injector->getPositionUnitBox(point, i_part, fac); + Real x = overlap_corner[0] + (iv[0] + point[0])*dx[0]; #if ( AMREX_SPACEDIM == 3 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; - Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1]; - Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2]; + Real y = overlap_corner[1] + (iv[1] + point[1])*dx[1]; + Real z = overlap_corner[2] + (iv[2] + point[2])*dx[2]; #elif ( AMREX_SPACEDIM == 2 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; Real y = 0; - Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1]; +#ifdef WARPX_RZ + Real z = overlap_corner[1] + (iv[1] + point[2])*dx[1]; +#else + Real z = overlap_corner[1] + (iv[1] + point[1])*dx[1]; +#endif #endif // If the new particle is not inside the tile box, // go to the next generated particle. @@ -715,9 +731,16 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) Real yb = y; #ifdef WARPX_RZ - // Replace the x and y, choosing the angle randomly. + // Replace the x and y, setting an angle theta. // These x and y are used to get the momentum and density - Real theta = 2.*MathConst::pi*amrex::Random(); + Real theta; + if (WarpX::nmodes == 1) { + // With only 1 mode, the angle doesn't matter so + // choose it randomly. + theta = 2.*MathConst::pi*amrex::Random(); + } else { + theta = 2.*MathConst::pi*point[1]; + } x = xb*std::cos(theta); y = xb*std::sin(theta); #endif @@ -1133,6 +1156,7 @@ PhysicalParticleContainer::FieldGather (int lev, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), + &WarpX::nmodes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); @@ -1425,6 +1449,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*bxfab), BL_TO_FORTRAN_ANYD(*byfab), BL_TO_FORTRAN_ANYD(*bzfab), + &WarpX::nmodes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); @@ -1512,6 +1537,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*cbxfab), BL_TO_FORTRAN_ANYD(*cbyfab), BL_TO_FORTRAN_ANYD(*cbzfab), + &WarpX::nmodes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); } @@ -1853,6 +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, &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 2a3e8dd0d..d659b3854 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -426,6 +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, &ll4symtry, &l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 0e800bf1d..b2d86c587 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -42,6 +42,9 @@ namespace ParticleStringNames {"Bx", PIdx::Bx }, {"By", PIdx::By }, {"Bz", PIdx::Bz } +#ifdef WARPX_RZ + ,{"theta", PIdx::theta} +#endif }; } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 317f46fd4..a1517cb44 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -354,9 +354,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, tby.grow(ngJ); tbz.grow(ngJ); - local_jx[thread_num].resize(tbx); - local_jy[thread_num].resize(tby); - local_jz[thread_num].resize(tbz); + local_jx[thread_num].resize(tbx, jx->nComp()); + local_jy[thread_num].resize(tby, jy->nComp()); + local_jz[thread_num].resize(tbz, jz->nComp()); Real* jx_ptr = local_jx[thread_num].dataPtr(); Real* jy_ptr = local_jy[thread_num].dataPtr(); @@ -379,6 +379,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &np_to_depose, m_xp[thread_num].dataPtr() + offset, m_yp[thread_num].dataPtr() + offset, @@ -399,6 +400,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &xyzmin[0], &dx[0]); #endif BL_PROFILE_VAR_STOP(blp_pxr_cd); @@ -407,9 +409,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, BL_PROFILE_VAR_START(blp_accumulate); // CPU, tiling: atomicAdd local_jx into jx // (same for jx and jz) - (*jx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); - (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, 1); - (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, 1); + (*jx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, local_jx[thread_num].nComp()); + (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, local_jy[thread_num].nComp()); + (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, local_jz[thread_num].nComp()); BL_PROFILE_VAR_STOP(blp_accumulate); #endif } @@ -440,11 +442,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array<Real, 3>& xyzmin = xyzmin_tile; #ifdef AMREX_USE_GPU - data_ptr = (*rhomf)[pti].dataPtr(icomp); + data_ptr = (*rhomf)[pti].dataPtr(icomp*(rhomf->nComp()/2)); auto rholen = (*rhomf)[pti].length(); #else tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box); + local_rho[thread_num].resize(tile_box, rhomf->nComp()); data_ptr = local_rho[thread_num].dataPtr(); auto rholen = local_rho[thread_num].length(); @@ -483,7 +485,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); + (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp*(rhomf->nComp()/2), (rhomf->nComp()/2)); BL_PROFILE_VAR_STOP(blp_accumulate); #endif @@ -502,7 +504,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #else tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box); + local_rho[thread_num].resize(tile_box, crhomf->nComp()); data_ptr = local_rho[thread_num].dataPtr(); auto rholen = local_rho[thread_num].length(); @@ -543,7 +545,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); + (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp*(crhomf->nComp()/2), (crhomf->nComp()/2)); BL_PROFILE_VAR_STOP(blp_accumulate); #endif @@ -598,7 +600,7 @@ WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, BoxArray coarsened_fine_BA = fine_BA; coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, 1, 0); + MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); coarsened_fine_data.setVal(0.0); IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME @@ -629,7 +631,7 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) const int ng = WarpX::nox; - auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,1,ng)); + auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,WarpX::ncomps,ng)); rho->setVal(0.0); #ifdef _OPENMP @@ -657,7 +659,7 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); const std::array<Real, 3>& xyzmin = xyzmin_tile; tile_box.grow(ng); - rho_loc.resize(tile_box); + rho_loc.resize(tile_box, rho->nComp()); rho_loc = 0.0; data_ptr = rho_loc.dataPtr(); auto rholen = rho_loc.length(); |