aboutsummaryrefslogtreecommitdiff
path: root/Source/Particles
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Particles')
-rw-r--r--Source/Particles/Deposition/CurrentDeposition.H109
-rw-r--r--Source/Particles/Gather/FieldGather.H80
-rw-r--r--Source/Particles/MultiParticleContainer.cpp2
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp42
-rw-r--r--Source/Particles/RigidInjectedParticleContainer.cpp3
-rw-r--r--Source/Particles/WarpXParticleContainer.H3
-rw-r--r--Source/Particles/WarpXParticleContainer.cpp46
7 files changed, 214 insertions, 71 deletions
diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H
index e8e7fab0a..608c0109a 100644
--- a/Source/Particles/Deposition/CurrentDeposition.H
+++ b/Source/Particles/Deposition/CurrentDeposition.H
@@ -2,6 +2,7 @@
#define CURRENTDEPOSITION_H_
#include "ShapeFactors.H"
+#include <WarpX_Complex.H>
/* \brief Current Deposition for thread thread_num
* /param xp, yp, zp : Pointer to arrays of particle positions.
@@ -169,7 +170,7 @@ void doDepositionShapeN(const amrex::Real * const xp,
}
/* \brief Esirkepov Current Deposition for thread thread_num
- * /param xp, yp, zp : Pointer to arrays of particle positions.
+ * \param xp, yp, zp : Pointer to arrays of particle positions.
* \param wp : Pointer to array of particle weights.
* \param uxp uyp uzp : Pointer to arrays of particle momentum.
* \param ion_lev : Pointer to array of particle ionization level. This is
@@ -184,7 +185,8 @@ void doDepositionShapeN(const amrex::Real * const xp,
* \param dx : 3D cell size
* \param xyzmin : Physical lower bounds of domain.
* \param lo : Index lower bounds of domain.
- * /param q : species charge.
+ * \param q : species charge.
+ * \param n_rz_azimuthal_modes: Number of azimuthal modes when using RZ geometry
*/
template <int depos_order>
void doEsirkepovDepositionShapeN (const amrex::Real * const xp,
@@ -203,7 +205,8 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp,
const std::array<amrex::Real,3>& dx,
const std::array<amrex::Real, 3> xyzmin,
const amrex::Dim3 lo,
- const amrex::Real q)
+ const amrex::Real q,
+ const long n_rz_azimuthal_modes)
{
// Whether ion_lev is a null pointer (do_ionization=0) or a real pointer
// (do_ionization=1)
@@ -230,6 +233,10 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp,
const amrex::Real invvol = 1.0/(dx[0]*dx[2]);
#endif
+#if (defined WARPX_DIM_RZ)
+ const Complex I = Complex{0., 1.};
+#endif
+
const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c;
// Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr
@@ -255,11 +262,42 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp,
// computes current and old position in grid units
#if (defined WARPX_DIM_RZ)
- const amrex::Real r_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]);
- const amrex::Real r_old = std::sqrt((xp[ip] - dt*uxp[ip]*gaminv)*(xp[ip] - dt*uxp[ip]*gaminv) +
- (yp[ip] - dt*uyp[ip]*gaminv)*(yp[ip] - dt*uyp[ip]*gaminv));
- const amrex::Real x_new = (r_new - xmin)*dxi;
- const amrex::Real x_old = (r_old - xmin)*dxi;
+ const amrex::Real xp_mid = xp[ip] - 0.5*dt*uxp[ip]*gaminv;
+ const amrex::Real yp_mid = yp[ip] - 0.5*dt*uyp[ip]*gaminv;
+ const amrex::Real xp_old = xp[ip] - dt*uxp[ip]*gaminv;
+ const amrex::Real yp_old = yp[ip] - dt*uyp[ip]*gaminv;
+ const amrex::Real rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]);
+ const amrex::Real rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid);
+ const amrex::Real rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
+ amrex::Real costheta_new, sintheta_new;
+ if (rp_new > 0.) {
+ costheta_new = xp[ip]/rp_new;
+ sintheta_new = yp[ip]/rp_new;
+ } else {
+ costheta_new = 1.;
+ sintheta_new = 0.;
+ }
+ amrex::Real costheta_mid, sintheta_mid;
+ if (rp_mid > 0.) {
+ costheta_mid = xp_mid/rp_mid;
+ sintheta_mid = yp_mid/rp_mid;
+ } else {
+ costheta_mid = 1.;
+ sintheta_mid = 0.;
+ }
+ amrex::Real costheta_old, sintheta_old;
+ if (rp_old > 0.) {
+ costheta_old = xp_old/rp_old;
+ sintheta_old = yp_old/rp_old;
+ } else {
+ costheta_old = 1.;
+ sintheta_old = 0.;
+ }
+ const Complex xy_new0 = Complex{costheta_new, sintheta_new};
+ const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
+ const Complex xy_old0 = Complex{costheta_old, sintheta_old};
+ const amrex::Real x_new = (rp_new - xmin)*dxi;
+ const amrex::Real x_old = (rp_old - xmin)*dxi;
#else
const amrex::Real x_new = (xp[ip] - xmin)*dxi;
const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv;
@@ -272,16 +310,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp,
const amrex::Real z_old = z_new - dtsdz0*uzp[ip]*gaminv;
#if (defined WARPX_DIM_RZ)
- amrex::Real costheta;
- amrex::Real sintheta;
- if (r_new > 0.) {
- costheta = xp[ip]/r_new;
- sintheta = yp[ip]/r_new;
- } else {
- costheta = 1.;
- sintheta = 0.;
- }
- const amrex::Real vy = (-uxp[ip]*sintheta + uyp[ip]*costheta)*gaminv;
+ const amrex::Real vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv;
#elif (defined WARPX_DIM_2D)
const amrex::Real vy = uyp[ip]*gaminv;
#endif
@@ -363,25 +392,61 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp,
amrex::Real sdxi = 0.;
for (int i=dil; i<=depos_order+1-diu; i++) {
sdxi += wqx*(sx_old[i] - sx_new[i])*(sz_new[k] + 0.5*(sz_old[k] - sz_new[k]));
- amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0), sdxi);
+ amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi);
+#if (defined WARPX_DIM_RZ)
+ Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
+ for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
+ // The factor 2 comes from the normalization of the modes
+ const Complex djr_cmplx = 2.*sdxi*xy_mid;
+ amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real());
+ amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag());
+ xy_mid = xy_mid*xy_mid0;
+ }
+#endif
}
}
for (int k=dkl; k<=depos_order+2-dku; k++) {
for (int i=dil; i<=depos_order+2-diu; i++) {
const amrex::Real sdyj = wq*vy*invvol*((sz_new[k] + 0.5*(sz_old[k] - sz_new[k]))*sx_new[i] +
(0.5*sz_new[k] + 1./3.*(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i]));
- amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0), sdyj);
+ amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj);
+#if (defined WARPX_DIM_RZ)
+ Complex xy_new = xy_new0;
+ Complex xy_mid = xy_mid0;
+ Complex xy_old = xy_old0;
+ // Throughout the following loop, xy_ takes the value e^{i m theta_}
+ for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
+ // The factor 2 comes from the normalization of the modes
+ // The minus sign comes from the different convention with respect to Davidson et al.
+ const Complex djt_cmplx = -2.*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(double)imode*
+ (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old));
+ amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real());
+ amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag());
+ xy_new = xy_new*xy_new0;
+ xy_mid = xy_mid*xy_mid0;
+ xy_old = xy_old*xy_old0;
+ }
+#endif
}
}
for (int i=dil; i<=depos_order+2-diu; i++) {
amrex::Real sdzk = 0.;
for (int k=dkl; k<=depos_order+1-dku; k++) {
sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5*(sx_old[i] - sx_new[i]));
- amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0), sdzk);
+ amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk);
+#if (defined WARPX_DIM_RZ)
+ Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta}
+ for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) {
+ // The factor 2 comes from the normalization of the modes
+ const Complex djz_cmplx = 2.*sdzk*xy_mid;
+ amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real());
+ amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag());
+ xy_mid = xy_mid*xy_mid0;
+ }
+#endif
}
}
-
#endif
}
);
diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H
index 1978e8141..d8d1d78ef 100644
--- a/Source/Particles/Gather/FieldGather.H
+++ b/Source/Particles/Gather/FieldGather.H
@@ -2,9 +2,10 @@
#define FIELDGATHER_H_
#include "ShapeFactors.H"
+#include <WarpX_Complex.H>
/* \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 <int depos_order, int lower_in_v>
void doGatherShapeN(const amrex::Real * const xp,
@@ -33,7 +35,8 @@ void doGatherShapeN(const amrex::Real * const xp,
const std::array<amrex::Real, 3>& dx,
const std::array<amrex::Real, 3> 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-1)*xy.real()
+ - ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode)*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-1)*xy.real()
+ - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*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-1)*xy.real()
+ - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*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-1)*xy.real()
+ - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*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-1)*xy.real()
+ - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*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-1)*xy.real()
+ - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode)*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/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index 143f4b7f9..7803bdae1 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -257,7 +257,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 73acd60fa..1513297a1 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -160,12 +160,12 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
npart /= 4;
}
for (long i = 0; i < npart; ++i) {
-#if ( AMREX_SPACEDIM == 3 | WARPX_DIM_RZ)
+#if (defined WARPX_DIM_3D) || (WARPX_DIM_RZ)
Real weight = q_tot/npart/charge;
Real x = distx(mt);
Real y = disty(mt);
Real z = distz(mt);
-#elif ( AMREX_SPACEDIM == 2 )
+#elif (defined WARPX_DIM_2D)
Real weight = q_tot/npart/charge/y_rms;
Real x = distx(mt);
Real y = 0.;
@@ -332,6 +332,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
Real density_max = plasma_injector->density_max;
#ifdef WARPX_DIM_RZ
+ const long nmodes = WarpX::n_rz_azimuthal_modes;
bool radially_weighted = plasma_injector->radially_weighted;
#endif
@@ -489,7 +490,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
#else
Real x = overlap_corner[0] + (iv[0]+r.x)*dx[0];
Real y = 0.0;
+#ifdef WARPX_DIM_2D
Real z = overlap_corner[1] + (iv[1]+r.y)*dx[1];
+#elif defined WARPX_DIM_RZ
+ // Note that for RZ, r.y will be theta
+ Real z = overlap_corner[1] + (iv[1]+r.z)*dx[1];
+#endif
#endif
#if (AMREX_SPACEDIM == 3)
@@ -510,9 +516,16 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
Real yb = y;
#ifdef WARPX_DIM_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 (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*r.y;
+ }
x = xb*std::cos(theta);
y = xb*std::sin(theta);
#endif
@@ -903,7 +916,8 @@ 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, 0, np, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal,
+ 0, np, thread_num, lev, lev);
if (cost) {
const Box& tbx = pti.tilebox();
@@ -1201,7 +1215,8 @@ 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, 0, np_gather, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal,
+ 0, np_gather, thread_num, lev, lev);
if (np_gather < np)
{
@@ -1646,7 +1661,8 @@ 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, 0, np, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal,
+ 0, np, thread_num, lev, lev);
// This wraps the momentum advance so that inheritors can modify the call.
// Extract pointers to the different particle quantities
@@ -1939,7 +1955,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,
@@ -1947,7 +1963,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,
@@ -1955,7 +1971,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){
@@ -1965,7 +1981,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,
@@ -1973,7 +1989,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,
@@ -1981,7 +1997,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 0c94d1e33..4893b3294 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -409,7 +409,8 @@ 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, 0, np, thread_num, lev, lev);
+ Ex.nGrow(), e_is_nodal,
+ 0, np, thread_num, lev, lev);
// Save the position and momenta, making copies
auto uxp_save = uxp;
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index c4bb5e410..bc46a116b 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -51,6 +51,9 @@ namespace ParticleStringNames
{"Bx", PIdx::Bx },
{"By", PIdx::By },
{"Bz", PIdx::Bz }
+#ifdef WARPX_DIM_RZ
+ ,{"theta", PIdx::theta}
+#endif
};
}
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index b34e71178..8d2499a35 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -337,9 +337,9 @@ WarpXParticleContainer::DepositCurrentFortran(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();
@@ -362,6 +362,7 @@ WarpXParticleContainer::DepositCurrentFortran(WarpXParIter& pti,
jx_ptr, &ngJ, jxntot.getVect(),
jy_ptr, &ngJ, jyntot.getVect(),
jz_ptr, &ngJ, jzntot.getVect(),
+ &WarpX::n_rz_azimuthal_modes,
&np_to_depose,
m_xp[thread_num].dataPtr() + offset,
m_yp[thread_num].dataPtr() + offset,
@@ -382,9 +383,9 @@ WarpXParticleContainer::DepositCurrentFortran(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
}
@@ -460,9 +461,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());
// local_jx[thread_num] is set to zero
local_jx[thread_num].setVal(0.0);
@@ -496,17 +497,20 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti,
doEsirkepovDepositionShapeN<1>(
xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset,
uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev,
- jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q);
+ jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q,
+ WarpX::n_rz_azimuthal_modes);
} else if (WarpX::nox == 2){
doEsirkepovDepositionShapeN<2>(
xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset,
uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev,
- jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q);
+ jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q,
+ WarpX::n_rz_azimuthal_modes);
} else if (WarpX::nox == 3){
doEsirkepovDepositionShapeN<3>(
xp, yp, zp, wp.dataPtr() + offset, uxp.dataPtr() + offset,
uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev,
- jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q);
+ jx_arr, jy_arr, jz_arr, np_to_depose, dt, dx, xyzmin, lo, q,
+ WarpX::n_rz_azimuthal_modes);
}
} else {
if (WarpX::nox == 1){
@@ -534,9 +538,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, jx->nComp());
+ (*jy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, jy->nComp());
+ (*jz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, jz->nComp());
BL_PROFILE_VAR_STOP(blp_accumulate);
#endif
}
@@ -593,15 +597,17 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp,
tilebox.grow(ngRho);
+ const int nc = (rho->nComp() == 1 ? 1 : rho->nComp()/2);
+
#ifdef AMREX_USE_GPU
// No tiling on GPU: rho_arr points to the full rho array.
- MultiFab rhoi(*rho, amrex::make_alias, icomp, 1);
+ MultiFab rhoi(*rho, amrex::make_alias, icomp*nc, nc);
Array4<Real> const& rho_arr = rhoi.array(pti);
#else
// Tiling is on: rho_arr points to local_rho[thread_num]
const Box tb = amrex::convert(tilebox, IntVect::TheUnitVector());
- local_rho[thread_num].resize(tb);
+ local_rho[thread_num].resize(tb, nc);
// local_rho[thread_num] is set to zero
local_rho[thread_num].setVal(0.0);
@@ -637,7 +643,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp,
#ifndef AMREX_USE_GPU
BL_PROFILE_VAR_START(blp_accumulate);
- (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp, 1);
+ (*rho)[pti].atomicAdd(local_rho[thread_num], tb, tb, 0, icomp*nc, nc);
BL_PROFILE_VAR_STOP(blp_accumulate);
#endif
@@ -691,7 +697,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
@@ -720,7 +726,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
generated by cgit v1.2.3 (git 2.46.0) at 2025-10-02 10:23:53 +0000