aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Evolve/WarpXEvolveEM.cpp16
-rw-r--r--Source/Particles/Deposition/CurrentDeposition.H139
-rw-r--r--Source/Particles/Gather/FieldGather.H171
-rw-r--r--Source/Particles/Gather/Make.package3
-rw-r--r--Source/Particles/Make.package2
-rw-r--r--Source/Particles/MultiParticleContainer.H6
-rw-r--r--Source/Particles/MultiParticleContainer.cpp8
-rw-r--r--Source/Particles/PhysicalParticleContainer.H34
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp178
-rw-r--r--Source/Particles/ShapeFactors.H60
-rw-r--r--Source/Particles/WarpXParticleContainer.H6
-rw-r--r--Source/Particles/WarpXParticleContainer.cpp6
12 files changed, 493 insertions, 136 deletions
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp
index 32a4747db..f1f7c698c 100644
--- a/Source/Evolve/WarpXEvolveEM.cpp
+++ b/Source/Evolve/WarpXEvolveEM.cpp
@@ -189,9 +189,9 @@ WarpX::EvolveEM (int numsteps)
UpdateAuxilaryData();
for (int lev = 0; lev <= finest_level; ++lev) {
- mypc->FieldGather(lev,
- *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2],
- *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]);
+ mypc->FieldGatherFortran(lev,
+ *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2],
+ *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]);
}
last_plot_file_step = step+1;
@@ -241,11 +241,11 @@ WarpX::EvolveEM (int numsteps)
UpdateAuxilaryData();
for (int lev = 0; lev <= finest_level; ++lev) {
- mypc->FieldGather(lev,
- *Efield_aux[lev][0],*Efield_aux[lev][1],
- *Efield_aux[lev][2],
- *Bfield_aux[lev][0],*Bfield_aux[lev][1],
- *Bfield_aux[lev][2]);
+ mypc->FieldGatherFortran(lev,
+ *Efield_aux[lev][0],*Efield_aux[lev][1],
+ *Efield_aux[lev][2],
+ *Bfield_aux[lev][0],*Bfield_aux[lev][1],
+ *Bfield_aux[lev][2]);
}
if (write_plot_file)
diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H
index 1db477b1d..b4531d3b3 100644
--- a/Source/Particles/Deposition/CurrentDeposition.H
+++ b/Source/Particles/Deposition/CurrentDeposition.H
@@ -1,52 +1,7 @@
#ifndef CURRENTDEPOSITION_H_
#define CURRENTDEPOSITION_H_
-using namespace amrex;
-
-// Compute shape factor and return index of leftmost cell where
-// particle writes.
-// Specialized templates are defined below for orders 1, 2 and 3.
-template <int depos_order>
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-int compute_shape_factor(Real* const sx, Real xint);
-
-// Compute shape factor for order 1.
-template <>
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-int compute_shape_factor <1> (Real* const sx, Real xmid){
- int j = (int) xmid;
- Real xint = xmid-j;
- sx[0] = 1.0 - xint;
- sx[1] = xint;
- return j;
-}
-
-// Compute shape factor for order 2.
-template <>
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-int compute_shape_factor <2> (Real* const sx, Real xmid){
- int j = (int) (xmid+0.5);
- Real xint = xmid-j;
- sx[0] = 0.5*(0.5-xint)*(0.5-xint);
- sx[1] = 0.75-xint*xint;
- sx[2] = 0.5*(0.5+xint)*(0.5+xint);
- // index of the leftmost cell where particle deposits
- return j-1;
-}
-
-// Compute shape factor for order 3.
-template <>
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-int compute_shape_factor <3> (Real* const sx, Real xmid){
- int j = (int) xmid;
- Real xint = xmid-j;
- sx[0] = 1.0/6.0*(1.0-xint)*(1.0-xint)*(1.0-xint);
- sx[1] = 2.0/3.0-xint*xint*(1-xint/2.0);
- sx[2] = 2.0/3.0-(1-xint)*(1-xint)*(1.0-0.5*(1-xint));
- sx[3] = 1.0/6.0*xint*xint*xint;
- // index of the leftmost cell where particle deposits
- return j-1;
-}
+#include "ShapeFactors.H"
/* \brief Current Deposition for thread thread_num
* /param xp, yp, zp : Pointer to arrays of particle positions.
@@ -55,9 +10,7 @@ int compute_shape_factor <3> (Real* const sx, Real xmid){
* \param jx_arr : Array4 of current density, either full array or tile.
* \param jy_arr : Array4 of current density, either full array or tile.
* \param jz_arr : Array4 of current density, either full array or tile.
- * \param offset : Index of first particle for which current is deposited
* \param np_to_depose : Number of particles for which current is deposited.
- Particles [offset,offset+np_tp_depose] deposit current.
* \param dt : Time step for particle level
* \param dx : 3D cell size
* \param xyzmin : Physical lower bounds of domain.
@@ -66,79 +19,83 @@ int compute_shape_factor <3> (Real* const sx, Real xmid){
* /param q : species charge.
*/
template <int depos_order>
-void doDepositionShapeN(const Real * const xp, const Real * const yp, const Real * const zp,
- const Real * const wp, const Real * const uxp,
- const Real * const uyp, const Real * const uzp,
- const amrex::Array4<amrex::Real>& jx_arr,
- const amrex::Array4<amrex::Real>& jy_arr,
+void doDepositionShapeN(const amrex::Real * const xp,
+ const amrex::Real * const yp,
+ const amrex::Real * const zp,
+ const amrex::Real * const wp,
+ const amrex::Real * const uxp,
+ const amrex::Real * const uyp,
+ const amrex::Real * const uzp,
+ const amrex::Array4<amrex::Real>& jx_arr,
+ const amrex::Array4<amrex::Real>& jy_arr,
const amrex::Array4<amrex::Real>& jz_arr,
- const long offset, const long np_to_depose,
- const amrex::Real dt, const std::array<amrex::Real,3>& dx,
+ const long np_to_depose, const amrex::Real dt,
+ const std::array<amrex::Real,3>& dx,
const std::array<Real, 3> xyzmin,
- const Dim3 lo,
+ const amrex::Dim3 lo,
const amrex::Real stagger_shift,
const amrex::Real q)
{
- const Real dxi = 1.0/dx[0];
- const Real dzi = 1.0/dx[2];
- const Real dts2dx = 0.5*dt*dxi;
- const Real dts2dz = 0.5*dt*dzi;
+ const amrex::Real dxi = 1.0/dx[0];
+ const amrex::Real dzi = 1.0/dx[2];
+ const amrex::Real dts2dx = 0.5*dt*dxi;
+ const amrex::Real dts2dz = 0.5*dt*dzi;
#if (AMREX_SPACEDIM == 2)
- const Real invvol = dxi*dzi;
+ const amrex::Real invvol = dxi*dzi;
#else // (AMREX_SPACEDIM == 3)
- const Real dyi = 1.0/dx[1];
- const Real dts2dy = 0.5*dt*dyi;
- const Real invvol = dxi*dyi*dzi;
+ const amrex::Real dyi = 1.0/dx[1];
+ const amrex::Real dts2dy = 0.5*dt*dyi;
+ const amrex::Real invvol = dxi*dyi*dzi;
#endif
- const Real xmin = xyzmin[0];
- const Real ymin = xyzmin[1];
- const Real zmin = xyzmin[2];
- const Real clightsq = 1.0/PhysConst::c/PhysConst::c;
+ const amrex::Real xmin = xyzmin[0];
+ const amrex::Real ymin = xyzmin[1];
+ const amrex::Real zmin = xyzmin[2];
+ const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c;
// Loop over particles and deposit into jx_arr, jy_arr and jz_arr
ParallelFor( np_to_depose,
[=] AMREX_GPU_DEVICE (long ip) {
// --- Get particle quantities
- const Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
- + uyp[ip]*uyp[ip]*clightsq
- + uzp[ip]*uzp[ip]*clightsq);
- const Real wq = q*wp[ip];
- const Real vx = uxp[ip]*gaminv;
- const Real vy = uyp[ip]*gaminv;
- const Real vz = uzp[ip]*gaminv;
+ const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq
+ + uyp[ip]*uyp[ip]*clightsq
+ + uzp[ip]*uzp[ip]*clightsq);
+ const amrex::Real wq = q*wp[ip];
+ const amrex::Real vx = uxp[ip]*gaminv;
+ const amrex::Real vy = uyp[ip]*gaminv;
+ const amrex::Real vz = uzp[ip]*gaminv;
// wqx, wqy wqz are particle current in each direction
- const Real wqx = wq*invvol*vx;
- const Real wqy = wq*invvol*vy;
- const Real wqz = wq*invvol*vz;
+ const amrex::Real wqx = wq*invvol*vx;
+ const amrex::Real wqy = wq*invvol*vy;
+ const amrex::Real wqz = wq*invvol*vz;
// --- Compute shape factors
// x direction
// Get particle position after 1/2 push back in position
- const Real xmid = (xp[ip]-xmin)*dxi-dts2dx*vx;
+ const amrex::Real xmid = (xp[ip]-xmin)*dxi-dts2dx*vx;
// Compute shape factors for node-centered quantities
- Real AMREX_RESTRICT sx [depos_order + 1];
+ amrex::Real AMREX_RESTRICT sx [depos_order + 1];
// j: leftmost grid point (node-centered) that the particle touches
const int j = compute_shape_factor<depos_order>(sx, xmid);
// Compute shape factors for cell-centered quantities
- Real AMREX_RESTRICT sx0[depos_order + 1];
+ amrex::Real AMREX_RESTRICT sx0[depos_order + 1];
// j0: leftmost grid point (cell-centered) that the particle touches
const int j0 = compute_shape_factor<depos_order>(sx0, xmid-stagger_shift);
#if (AMREX_SPACEDIM == 3)
// y direction
- const Real ymid= (yp[ip]-ymin)*dyi-dts2dy*vy;
- Real AMREX_RESTRICT sy [depos_order + 1];
- const int k = compute_shape_factor<depos_order>(sy, ymid);
- Real AMREX_RESTRICT sy0[depos_order + 1];
- const int k0 = compute_shape_factor<depos_order>(sy0, ymid-stagger_shift);
+ const amrex::Real ymid= (yp[ip]-ymin)*dyi-dts2dy*vy;
+ amrex::Real AMREX_RESTRICT sy [depos_order + 1];
+ const int k = compute_shape_factor<depos_order>(sy, ymid);
+ amrex::Real AMREX_RESTRICT sy0[depos_order + 1];
+ const int k0 = compute_shape_factor<depos_order>(sy0, ymid-stagger_shift);
#endif
// z direction
- const Real zmid= (zp[ip]-zmin)*dzi-dts2dz*vz;
- Real AMREX_RESTRICT sz [depos_order + 1];
- const int l = compute_shape_factor<depos_order>(sz, zmid);
- Real AMREX_RESTRICT sz0[depos_order + 1];
- const int l0 = compute_shape_factor<depos_order>(sz0, zmid-stagger_shift);
+ const amrex::Real zmid= (zp[ip]-zmin)*dzi-dts2dz*vz;
+ amrex::Real AMREX_RESTRICT sz [depos_order + 1];
+ const int l = compute_shape_factor<depos_order>(sz, zmid);
+ amrex::Real AMREX_RESTRICT sz0[depos_order + 1];
+ const int l0 = compute_shape_factor<depos_order>(sz0, zmid-stagger_shift);
// Deposit current into jx_arr, jy_arr and jz_arr
#if (AMREX_SPACEDIM == 2)
diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H
new file mode 100644
index 000000000..789c7f775
--- /dev/null
+++ b/Source/Particles/Gather/FieldGather.H
@@ -0,0 +1,171 @@
+#ifndef FIELDGATHER_H_
+#define FIELDGATHER_H_
+
+#include "ShapeFactors.H"
+
+using namespace amrex;
+
+/* \brief Field gather for particles handled by thread thread_num
+ * /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.
+ * \param ez_arr bx_arr: Array4 of current density, either full array or tile.
+ * \param by_arr bz_arr: Array4 of current density, either full array or tile.
+ * \param np_to_gather : Number of particles for which current is deposited.
+ * \param dx : 3D cell size
+ * \param xyzmin : Physical lower bounds of domain.
+ * \param lo : Index lower bounds of domain.
+ * \param stagger_shift: 0 if nodal, 0.5 if staggered.
+ */
+template <int depos_order, int lower_in_v>
+void doGatherShapeN(const Real * const xp,
+ const Real * const yp,
+ const Real * const zp,
+ Real * const Exp, Real * const Eyp, Real * const Ezp,
+ Real * const Bxp, Real * const Byp, Real * const Bzp,
+ const amrex::Array4<const amrex::Real>& ex_arr,
+ const amrex::Array4<const amrex::Real>& ey_arr,
+ const amrex::Array4<const amrex::Real>& ez_arr,
+ const amrex::Array4<const amrex::Real>& bx_arr,
+ const amrex::Array4<const amrex::Real>& by_arr,
+ const amrex::Array4<const amrex::Real>& bz_arr,
+ const long np_to_gather,
+ const std::array<amrex::Real,3>& dx,
+ const std::array<Real, 3> xyzmin,
+ const Dim3 lo,
+ const amrex::Real stagger_shift)
+{
+ const Real dxi = 1.0/dx[0];
+ const Real dzi = 1.0/dx[2];
+#if (AMREX_SPACEDIM == 3)
+ const Real dyi = 1.0/dx[1];
+#endif
+
+ const Real xmin = xyzmin[0];
+ const Real ymin = xyzmin[1];
+ const Real zmin = xyzmin[2];
+
+ // Loop over particles and gather fields from {e,b}{x,y,z}_arr to {E,B}{xyz}p.
+ ParallelFor( np_to_gather,
+ [=] AMREX_GPU_DEVICE (long ip) {
+ // --- Compute shape factors
+ // x direction
+ // Get particle position after 1/2 push back in position
+ const Real xmid = (xp[ip]-xmin)*dxi;
+ // Compute shape factors for node-centered quantities
+ Real AMREX_RESTRICT sx [depos_order + 1];
+ // j: leftmost grid point (node-centered) that the particle touches
+ const int j = compute_shape_factor<depos_order>(sx, xmid);
+ // Compute shape factors for cell-centered quantities
+ Real AMREX_RESTRICT sx0[depos_order + 1 - lower_in_v];
+ // j0: leftmost grid point (cell-centered) that the particle touches
+ const int j0 = compute_shape_factor<depos_order - lower_in_v>(sx0, xmid-stagger_shift);
+
+#if (AMREX_SPACEDIM == 3)
+ // y direction
+ const Real ymid= (yp[ip]-ymin)*dyi;
+ Real AMREX_RESTRICT sy [depos_order + 1];
+ const int k = compute_shape_factor<depos_order>(sy, ymid);
+ Real AMREX_RESTRICT sy0[depos_order + 1 - lower_in_v];
+ const int k0 = compute_shape_factor<depos_order-lower_in_v>(sy0, ymid-stagger_shift);
+#endif
+ // z direction
+ const Real zmid= (zp[ip]-zmin)*dzi;
+ Real AMREX_RESTRICT sz [depos_order + 1];
+ const int l = compute_shape_factor<depos_order>(sz, zmid);
+ Real AMREX_RESTRICT sz0[depos_order + 1 - lower_in_v];
+ const int l0 = compute_shape_factor<depos_order - lower_in_v>(sz0, zmid-stagger_shift);
+
+ // Deposit current into jx_arr, jy_arr and jz_arr
+ Exp[ip] = 0;
+ Eyp[ip] = 0;
+ Ezp[ip] = 0;
+ Bxp[ip] = 0;
+ Byp[ip] = 0;
+ Bzp[ip] = 0;
+#if (AMREX_SPACEDIM == 2)
+ 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);
+ }
+ }
+ 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);
+ Bzp[ip] += sx0[ix]*sz[iz]*
+ bz_arr(lo.x+j0+ix, lo.y+l +iz, 0);
+ }
+ }
+ 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);
+ Bxp[ip] += sx[ix]*sz0[iz]*
+ bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0);
+ }
+ }
+ 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);
+ }
+ }
+#else // (AMREX_SPACEDIM == 3)
+ for (int iz=0; iz<=depos_order; iz++){
+ for (int iy=0; iy<=depos_order; iy++){
+ for (int ix=0; ix<=depos_order-lower_in_v; ix++){
+ Exp[ip] += sx0[ix]*sy[iy]*sz[iz]*
+ ex_arr(lo.x+j0+ix, lo.y+k+iy, lo.z+l+iz);
+ }
+ }
+ }
+ for (int iz=0; iz<=depos_order; iz++){
+ for (int iy=0; iy<=depos_order-lower_in_v; iy++){
+ for (int ix=0; ix<=depos_order; ix++){
+ Eyp[ip] += sx[ix]*sy0[iy]*sz[iz]*
+ ey_arr(lo.x+j+ix, lo.y+k0+iy, lo.z+l+iz);
+ }
+ }
+ }
+ for (int iz=0; iz<=depos_order-lower_in_v; iz++){
+ for (int iy=0; iy<=depos_order; iy++){
+ for (int ix=0; ix<=depos_order; ix++){
+ Ezp[ip] += sx[ix]*sy[iy]*sz0[iz]*
+ ez_arr(lo.x+j+ix, lo.y+k+iy, lo.z+l0+iz);
+ }
+ }
+ }
+
+ for (int iz=0; iz<=depos_order; iz++){
+ for (int iy=0; iy<=depos_order-lower_in_v; iy++){
+ for (int ix=0; ix<=depos_order-lower_in_v; ix++){
+ Bzp[ip] += sx0[ix]*sy0[iy]*sz[iz]*
+ bz_arr(lo.x+j0+ix, lo.y+k0+iy, lo.z+l+iz);
+ }
+ }
+ }
+ for (int iz=0; iz<=depos_order-lower_in_v; iz++){
+ for (int iy=0; iy<=depos_order; iy++){
+ for (int ix=0; ix<=depos_order-lower_in_v; ix++){
+ Byp[ip] += sx0[ix]*sy[iy]*sz0[iz]*
+ by_arr(lo.x+j0+ix, lo.y+k+iy, lo.z+l0+iz);
+ }
+ }
+ }
+ for (int iz=0; iz<=depos_order-lower_in_v; iz++){
+ for (int iy=0; iy<=depos_order-lower_in_v; iy++){
+ for (int ix=0; ix<=depos_order; ix++){
+ Bxp[ip] += sx[ix]*sy0[iy]*sz0[iz]*
+ bx_arr(lo.x+j+ix, lo.y+k0+iy, lo.z+l0+iz);
+ }
+ }
+ }
+#endif
+ }
+ );
+}
+
+#endif // FIELDGATHER_H_
diff --git a/Source/Particles/Gather/Make.package b/Source/Particles/Gather/Make.package
new file mode 100644
index 000000000..10abfcaaf
--- /dev/null
+++ b/Source/Particles/Gather/Make.package
@@ -0,0 +1,3 @@
+CEXE_headers += FieldGather.H
+INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Gather
+VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Gather
diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package
index 2038472a1..db90de1dc 100644
--- a/Source/Particles/Make.package
+++ b/Source/Particles/Make.package
@@ -9,9 +9,11 @@ CEXE_headers += MultiParticleContainer.H
CEXE_headers += WarpXParticleContainer.H
CEXE_headers += RigidInjectedParticleContainer.H
CEXE_headers += PhysicalParticleContainer.H
+CEXE_headers += ShapeFactors.H
include $(WARPX_HOME)/Source/Particles/Pusher/Make.package
include $(WARPX_HOME)/Source/Particles/Deposition/Make.package
+include $(WARPX_HOME)/Source/Particles/Gather/Make.package
INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles
diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H
index 869126fef..76e3c44bc 100644
--- a/Source/Particles/MultiParticleContainer.H
+++ b/Source/Particles/MultiParticleContainer.H
@@ -84,9 +84,9 @@ public:
/// Performs the field gather operation using the input fields E and B, for all the species
/// in the MultiParticleContainer. This is the electromagnetic version of the field gather.
///
- void FieldGather (int lev,
- const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez,
- const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz);
+ void FieldGatherFortran (int lev,
+ const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez,
+ const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz);
///
/// This evolves all the particles by one PIC time step, including current deposition, the
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index 9d39ec2f9..ab4400dad 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -239,12 +239,12 @@ MultiParticleContainer::sumParticleCharge (bool local)
#endif // WARPX_DO_ELECTROSTATIC
void
-MultiParticleContainer::FieldGather (int lev,
- const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
- const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
+MultiParticleContainer::FieldGatherFortran (int lev,
+ const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
+ const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
{
for (auto& pc : allcontainers) {
- pc->FieldGather(lev, Ex, Ey, Ez, Bx, By, Bz);
+ pc->FieldGatherFortran(lev, Ex, Ey, Ez, Bx, By, Bz);
}
}
diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H
index d55764682..ec9d7f98b 100644
--- a/Source/Particles/PhysicalParticleContainer.H
+++ b/Source/Particles/PhysicalParticleContainer.H
@@ -31,13 +31,33 @@ public:
amrex::Real t, amrex::Real dt) override;
#endif // WARPX_DO_ELECTROSTATIC
- virtual void FieldGather(int lev,
- const amrex::MultiFab& Ex,
- const amrex::MultiFab& Ey,
- const amrex::MultiFab& Ez,
- const amrex::MultiFab& Bx,
- const amrex::MultiFab& By,
- const amrex::MultiFab& Bz) final;
+ virtual void FieldGatherFortran(int lev,
+ const amrex::MultiFab& Ex,
+ const amrex::MultiFab& Ey,
+ const amrex::MultiFab& Ez,
+ const amrex::MultiFab& Bx,
+ const amrex::MultiFab& By,
+ const amrex::MultiFab& Bz) final;
+
+ void FieldGather(WarpXParIter& pti,
+ RealVector& Exp,
+ RealVector& Eyp,
+ RealVector& Ezp,
+ RealVector& Bxp,
+ RealVector& Byp,
+ RealVector& Bzp,
+ amrex::FArrayBox const * exfab,
+ amrex::FArrayBox const * eyfab,
+ amrex::FArrayBox const * ezfab,
+ amrex::FArrayBox const * bxfab,
+ amrex::FArrayBox const * byfab,
+ amrex::FArrayBox const * bzfab,
+ const int ngE, const int e_is_nodal,
+ const long offset,
+ const long np_to_gather,
+ int thread_num,
+ int lev,
+ int depos_lev);
virtual void Evolve (int lev,
const amrex::MultiFab& Ex,
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index d47a7b220..357a187fd 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -6,7 +6,7 @@
#include <WarpX.H>
#include <WarpXConst.H>
#include <WarpXWrappers.h>
-
+#include <FieldGather.H>
using namespace amrex;
@@ -1055,9 +1055,9 @@ PhysicalParticleContainer::EvolveES (const Vector<std::array<std::unique_ptr<Mul
#endif // WARPX_DO_ELECTROSTATIC
void
-PhysicalParticleContainer::FieldGather (int lev,
- const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
- const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
+PhysicalParticleContainer::FieldGatherFortran (int lev,
+ const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
+ const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
{
const std::array<Real,3>& dx = WarpX::CellSize(lev);
@@ -1395,19 +1395,19 @@ PhysicalParticleContainer::Evolve (int lev,
if (! do_not_push)
{
+ const long np_gather = (cEx) ? nfine_gather : np;
+
+ int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
+
//
// Field Gather of Aux Data (i.e., the full solution)
//
+ BL_PROFILE_VAR_START(blp_pxr_fg);
+#ifdef WARPX_RZ
const int ll4symtry = false;
long lvect_fieldgathe = 64;
-
const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev);
const int* ixyzmin_grid = box.loVect();
-
- const long np_gather = (cEx) ? nfine_gather : np;
-
- BL_PROFILE_VAR_START(blp_pxr_fg);
-
warpx_geteb_energy_conserving(
&np_gather,
m_xp[thread_num].dataPtr(),
@@ -1427,6 +1427,11 @@ PhysicalParticleContainer::Evolve (int lev,
BL_TO_FORTRAN_ANYD(*bzfab),
&ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal,
&lvect_fieldgathe, &WarpX::field_gathering_algo);
+#else
+ 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);
+#endif
if (np_gather < np)
{
@@ -1435,13 +1440,14 @@ PhysicalParticleContainer::Evolve (int lev,
const std::array<Real,3>& cxyzmin_grid = WarpX::LowerCorner(cbox, lev-1);
const int* cixyzmin_grid = cbox.loVect();
- const FArrayBox* cexfab = &(*cEx)[pti];
- const FArrayBox* ceyfab = &(*cEy)[pti];
- const FArrayBox* cezfab = &(*cEz)[pti];
- const FArrayBox* cbxfab = &(*cBx)[pti];
- const FArrayBox* cbyfab = &(*cBy)[pti];
- const FArrayBox* cbzfab = &(*cBz)[pti];
-
+ // Data on the grid
+ FArrayBox const* cexfab = &(*cEx)[pti];
+ FArrayBox const* ceyfab = &(*cEy)[pti];
+ FArrayBox const* cezfab = &(*cEz)[pti];
+ FArrayBox const* cbxfab = &(*cBx)[pti];
+ FArrayBox const* cbyfab = &(*cBy)[pti];
+ FArrayBox const* cbzfab = &(*cBz)[pti];
+
if (WarpX::use_fdtd_nci_corr)
{
#if (AMREX_SPACEDIM == 2)
@@ -1494,6 +1500,9 @@ PhysicalParticleContainer::Evolve (int lev,
#endif
}
+ // Field gather for particles in gather buffers
+#ifdef WARPX_RZ
+
long ncrse = np - nfine_gather;
warpx_geteb_energy_conserving(
&ncrse,
@@ -1514,6 +1523,15 @@ PhysicalParticleContainer::Evolve (int lev,
BL_TO_FORTRAN_ANYD(*cbzfab),
&ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal,
&lvect_fieldgathe, &WarpX::field_gathering_algo);
+#else
+ e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal();
+ FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+ cexfab, ceyfab, cezfab,
+ cbxfab, cbyfab, cbzfab,
+ cEx->nGrow(), e_is_nodal,
+ nfine_gather, np-nfine_gather,
+ thread_num, lev, lev-1);
+#endif
}
BL_PROFILE_VAR_STOP(blp_pxr_fg);
@@ -2112,3 +2130,129 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box)
const int lev=0;
AddPlasma(lev, injection_box);
}
+
+/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab,
+ * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp.
+ * \param Exp-Bzp: fields on particles.
+ * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti
+ * \param ngE: number of guard cells for E
+ * \param e_is_nodal: 0 if E is staggered, 1 if E is nodal
+ * \param offset: index of first particle for which fields are gathered
+ * \param np_to_gather: number of particles onto which fields are gathered
+ * \param thread_num: if using OpenMP, thread number
+ * \param lev: level on which particles are located
+ * \param gather_lev: level from which particles gather fields (lev-1) for
+ particles in buffers.
+ */
+void
+PhysicalParticleContainer::FieldGather(WarpXParIter& pti,
+ RealVector& Exp,
+ RealVector& Eyp,
+ RealVector& Ezp,
+ RealVector& Bxp,
+ RealVector& Byp,
+ RealVector& Bzp,
+ FArrayBox const * exfab,
+ FArrayBox const * eyfab,
+ FArrayBox const * ezfab,
+ FArrayBox const * bxfab,
+ FArrayBox const * byfab,
+ FArrayBox const * bzfab,
+ const int ngE, const int e_is_nodal,
+ const long offset,
+ const long np_to_gather,
+ int thread_num,
+ int lev,
+ int gather_lev)
+{
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) ||
+ (gather_lev==(lev )),
+ "Gather buffers only work for lev-1");
+
+ // If no particles, do not do anything
+ if (np_to_gather == 0) return;
+
+ const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0));
+ const Real stagger_shift = e_is_nodal ? 0.0 : 0.5;
+
+ // Get box =from which field is gathered.
+ Box box;
+ if (lev == gather_lev) {
+ box = pti.tilebox();
+ } else {
+ const IntVect& ref_ratio = WarpX::RefRatio(gather_lev);
+ box = amrex::coarsen(pti.tilebox(),ref_ratio);
+ }
+
+ box.grow(ngE);
+
+ const Array4<const Real>& ex_arr = exfab->array();
+ const Array4<const Real>& ey_arr = eyfab->array();
+ const Array4<const Real>& ez_arr = ezfab->array();
+ const Array4<const Real>& bx_arr = bxfab->array();
+ const Array4<const Real>& by_arr = byfab->array();
+ const Array4<const Real>& bz_arr = bzfab->array();
+
+ const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset;
+ const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset;
+ const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset;
+
+ // Lower corner of tile box physical domain
+ const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev);
+
+ const Dim3 lo = lbound(box);
+
+ if (WarpX::l_lower_order_in_v){
+ if (WarpX::nox == 1){
+ doGatherShapeN<1,1>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ 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);
+ } else if (WarpX::nox == 2){
+ doGatherShapeN<2,1>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ 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);
+ } else if (WarpX::nox == 3){
+ doGatherShapeN<3,1>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ 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);
+ }
+ } else {
+ if (WarpX::nox == 1){
+ doGatherShapeN<1,0>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ 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);
+ } else if (WarpX::nox == 2){
+ doGatherShapeN<2,0>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ 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);
+ } else if (WarpX::nox == 3){
+ doGatherShapeN<3,0>(xp, yp, zp,
+ Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
+ Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
+ 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);
+ }
+ }
+}
diff --git a/Source/Particles/ShapeFactors.H b/Source/Particles/ShapeFactors.H
new file mode 100644
index 000000000..6258f05af
--- /dev/null
+++ b/Source/Particles/ShapeFactors.H
@@ -0,0 +1,60 @@
+#ifndef SHAPEFACTORS_H_
+#define SHAPEFACTORS_H_
+
+using namespace amrex;
+
+// Compute shape factor and return index of leftmost cell where
+// particle writes.
+// Specialized templates are defined below for orders 1, 2 and 3.
+template <int depos_order>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int compute_shape_factor(Real* const sx, Real xint) {return 0;};
+
+// Compute shape factor for order 0.
+template <>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int compute_shape_factor <0> (Real* const sx, Real xmid){
+ int j = (int) (xmid+0.5);
+ sx[0] = 1.0;
+ return j;
+}
+
+// Compute shape factor for order 1.
+template <>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int compute_shape_factor <1> (Real* const sx, Real xmid){
+ int j = (int) xmid;
+ Real xint = xmid-j;
+ sx[0] = 1.0 - xint;
+ sx[1] = xint;
+ return j;
+}
+
+// Compute shape factor for order 2.
+template <>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int compute_shape_factor <2> (Real* const sx, Real xmid){
+ int j = (int) (xmid+0.5);
+ Real xint = xmid-j;
+ sx[0] = 0.5*(0.5-xint)*(0.5-xint);
+ sx[1] = 0.75-xint*xint;
+ sx[2] = 0.5*(0.5+xint)*(0.5+xint);
+ // index of the leftmost cell where particle deposits
+ return j-1;
+}
+
+// Compute shape factor for order 3.
+template <>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int compute_shape_factor <3> (Real* const sx, Real xmid){
+ int j = (int) xmid;
+ Real xint = xmid-j;
+ sx[0] = 1.0/6.0*(1.0-xint)*(1.0-xint)*(1.0-xint);
+ sx[1] = 2.0/3.0-xint*xint*(1-xint/2.0);
+ sx[2] = 2.0/3.0-(1-xint)*(1-xint)*(1.0-0.5*(1-xint));
+ sx[3] = 1.0/6.0*xint*xint*xint;
+ // index of the leftmost cell where particle deposits
+ return j-1;
+}
+
+#endif // SHAPEFACTORS_H_
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index 662b2e1b8..1b2a06171 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -103,9 +103,9 @@ public:
virtual void FieldGatherES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E,
const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks) {}
- virtual void FieldGather (int lev,
- const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez,
- const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) {}
+ virtual void FieldGatherFortran (int lev,
+ const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez,
+ const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) {}
#ifdef WARPX_DO_ELECTROSTATIC
virtual void EvolveES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E,
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index a20f0035e..89f233b2c 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -532,17 +532,17 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti,
if (WarpX::nox == 1){
doDepositionShapeN<1>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(),
uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr,
- jz_arr, offset, np_to_depose, dt, dx,
+ jz_arr, np_to_depose, dt, dx,
xyzmin, lo, stagger_shift, q);
} else if (WarpX::nox == 2){
doDepositionShapeN<2>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(),
uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr,
- jz_arr, offset, np_to_depose, dt, dx,
+ jz_arr, np_to_depose, dt, dx,
xyzmin, lo, stagger_shift, q);
} else if (WarpX::nox == 3){
doDepositionShapeN<3>(xp, yp, zp, wp.dataPtr(), uxp.dataPtr(),
uyp.dataPtr(), uzp.dataPtr(), jx_arr, jy_arr,
- jz_arr, offset, np_to_depose, dt, dx,
+ jz_arr, np_to_depose, dt, dx,
xyzmin, lo, stagger_shift, q);
}
}
generated by cgit v1.2.3 (git 2.46.0) at 2025-09-01 18:05:18 +0000