Merge pull request #263 from ECP-WarpX/convert_charge_deposition

Convert charge deposition

Author: RemiLehe

Source/Evolve/WarpXEvolveEM.cpp

@@ -483,11 +483,17 @@ WarpX::PushParticlesandDepose (int lev, Real cur_time)
                  Bfield_cax[lev][0].get(), Bfield_cax[lev][1].get(), Bfield_cax[lev][2].get(),
                  cur_time, dt[lev]);
 #ifdef WARPX_DIM_RZ
-    // This is called after all particles have deposited their current.
+    // This is called after all particles have deposited their current and charge.
     ApplyInverseVolumeScalingToCurrentDensity(current_fp[lev][0].get(), current_fp[lev][1].get(), current_fp[lev][2].get(), lev);
     if (current_buf[lev][0].get()) {
         ApplyInverseVolumeScalingToCurrentDensity(current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), lev-1);
     }
+    if (rho_fp[lev].get()) {
+        ApplyInverseVolumeScalingToChargeDensity(rho_fp[lev].get(), lev);
+        if (charge_buf[lev].get()) {
+            ApplyInverseVolumeScalingToChargeDensity(charge_buf[lev].get(), lev-1);
+        }
+    }
 #endif
 }
 

Source/FieldSolver/WarpXPushFieldsEM.cpp

@@ -671,4 +671,53 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
         });
     }
 }
+
+void
+WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev)
+{
+    const long ngRho = Rho->nGrow();
+    const std::array<Real,3>& dx = WarpX::CellSize(lev);
+    const Real dr = dx[0];
+
+    Box tilebox;
+
+    for ( MFIter mfi(*Rho, TilingIfNotGPU()); mfi.isValid(); ++mfi )
+    {
+
+        Array4<Real> const& Rho_arr = Rho->array(mfi);
+
+        tilebox = mfi.tilebox();
+        Box tb = convert(tilebox, IntVect::TheUnitVector());
+
+        // Lower corner of tile box physical domain
+        // Note that this is done before the tilebox.grow so that
+        // these do not include the guard cells.
+        const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, lev);
+        const Dim3 lo = lbound(tilebox);
+        const Real rmin = xyzmin[0];
+        const int irmin = lo.x;
+
+        // Rescale charge in r-z mode since the inverse volume factor was not
+        // included in the charge deposition.
+        amrex::ParallelFor(tb, Rho->nComp(),
+        [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp)
+        {
+            // Wrap the charge density deposited in the guard cells around
+            // to the cells above the axis.
+            // Rho is located on the boundary
+            if (rmin == 0. && 0 < i && i <= ngRho) {
+                Rho_arr(i,j,0,icomp) += Rho_arr(-i,j,0,icomp);
+            }
+
+            // Apply the inverse volume scaling
+            const amrex::Real r = std::abs(rmin + (i - irmin)*dr);
+            if (r == 0.) {
+                // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis
+                Rho_arr(i,j,0,icomp) /= (MathConst::pi*dr/3.);
+            } else {
+                Rho_arr(i,j,0,icomp) /= (2.*MathConst::pi*r);
+            }
+        });
+    }
+}
 #endif

Source/FortranInterface/WarpX_f.H

@@ -75,22 +75,6 @@ extern "C"
 {
 #endif
 
-	// Charge deposition
-	void warpx_charge_deposition(amrex::Real* rho,
-            const long* np,	const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp,	const amrex::Real* w,
-			const amrex::Real* q, const amrex::Real* xmin, const amrex::Real* ymin, const amrex::Real* zmin,
-			const amrex::Real* dx, const amrex::Real* dy, const amrex::Real* dz,
-			const long* nx, const long* ny, const long* nz,
-			const long* nxguard, const long* nyguard, const long* nzguard,
-			const long* nox, const long* noy,const long* noz,
-			const long* lvect, const long* charge_depo_algo);
-
-        // Charge deposition finalize for RZ
-        void warpx_charge_deposition_rz_volume_scaling(
-			amrex::Real* rho, const long* rho_ng, const int* rho_ntot,
-			const amrex::Real* rmin,
-			const amrex::Real* dr);
-
 	// Current deposition
 	void warpx_current_deposition(
 			amrex::Real* jx, const long* jx_ng, const int* jx_ntot,

Source/FortranInterface/WarpX_picsar.F90

@@ -7,7 +7,6 @@
 #ifdef WARPX_DIM_RZ
 
 #define WRPX_PXR_CURRENT_DEPOSITION       depose_jrjtjz_generic_rz
-#define WRPX_PXR_RZ_VOLUME_SCALING_RHO    apply_rz_volume_scaling_rho
 
 #else
 
@@ -49,149 +48,6 @@ module warpx_to_pxr_module
 
 contains
 
-! _________________________________________________________________
-!>
-!> @brief
-!> Main subroutine for the charge deposition
-!>
-!> @details
-!> This subroutines enable to controle the interpolation order
-!> via the parameters nox,noy,noz and the type of algorithm via
-!> the parameter charge_depo_algo
-!
-!> @param[inout] rho charge array
-!> @param[in] np number of particles
-!> @param[in] xp,yp,zp particle position arrays
-!> @param[in] w particle weight arrays
-!> @param[in] q particle species charge
-!> @param[in] xmin,ymin,zmin tile grid minimum position
-!> @param[in] dx,dy,dz space discretization steps
-!> @param[in] nx,ny,nz number of cells
-!> @param[in] nxguard,nyguard,nzguard number of guard cells
-!> @param[in] nox,noy,noz interpolation order
-!> @param[in] lvect vector length
-!> @param[in] charge_depo_algo algorithm choice for the charge deposition
-!>
-subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-   nxguard,nyguard,nzguard,nox,noy,noz,lvect,charge_depo_algo) &
-  bind(C, name="warpx_charge_deposition")
-
-  integer(c_long), intent(IN)                   :: np
-  integer(c_long), intent(IN)                   :: nx,ny,nz
-  integer(c_long), intent(IN)                   :: nxguard,nyguard,nzguard
-  integer(c_long), intent(IN)                   :: nox,noy,noz
-  real(amrex_real), intent(IN OUT)              :: rho(*)
-  real(amrex_real), intent(IN)                  :: q
-  real(amrex_real), intent(IN)                  :: dx,dy,dz
-  real(amrex_real), intent(IN)                  :: xmin,ymin,zmin
-  real(amrex_real), intent(IN),  dimension(np)  :: xp,yp,zp,w
-  integer(c_long), intent(IN)                   :: lvect
-  integer(c_long), intent(IN)                   :: charge_depo_algo
-
-
-  ! Dimension 3
-#if (AMREX_SPACEDIM==3)
-
-  SELECT CASE(charge_depo_algo)
-
-  ! Scalar classical charge deposition subroutines
-  CASE(1)
-    IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN
-
-      CALL depose_rho_scalar_1_1_1(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-    nxguard,nyguard,nzguard,lvect)
-
-    ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN
-
-      CALL depose_rho_scalar_2_2_2(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-    nxguard,nyguard,nzguard,lvect)
-
-    ELSE IF ((nox.eq.3).and.(noy.eq.3).and.(noz.eq.3)) THEN
-
-      CALL depose_rho_scalar_3_3_3(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-    nxguard,nyguard,nzguard,lvect)
-
-    ELSE
-      CALL pxr_depose_rho_n(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-                  nxguard,nyguard,nzguard,nox,noy,noz, &
-                  .TRUE._c_long,.FALSE._c_long)
-    ENDIF
-
-  ! Optimized subroutines
-  CASE DEFAULT
-
-    IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN
-      CALL depose_rho_vecHVv2_1_1_1(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-               nxguard,nyguard,nzguard,lvect)
-
-    ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN
-      CALL depose_rho_vecHVv2_2_2_2(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-                 nxguard,nyguard,nzguard,lvect)
-
-    ELSE
-      CALL pxr_depose_rho_n(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,&
-                  nxguard,nyguard,nzguard,nox,noy,noz, &
-                  .TRUE._c_long,.FALSE._c_long)
-    ENDIF
-  END SELECT
-
-  ! Dimension 2
-#elif (AMREX_SPACEDIM==2)
-
-#ifdef WARPX_DIM_RZ
-  logical(pxr_logical) :: l_2drz = .TRUE._c_long
-#else
-  logical(pxr_logical) :: l_2drz = .FALSE._c_long
-#endif
-
-  CALL pxr_depose_rho_n_2dxz(rho,np,xp,yp,zp,w,q,xmin,zmin,dx,dz,nx,nz,&
-       nxguard,nzguard,nox,noz, &
-       .TRUE._c_long, .FALSE._c_long, l_2drz, 0_c_long)
-
-#endif
-
- end subroutine warpx_charge_deposition
-
-  ! _________________________________________________________________
-  !>
-  !> @brief
-  !> Applies the inverse volume scaling for RZ charge deposition
-  !>
-  !> @details
-  !> The scaling is done for both single mode (FDTD) and
-  !> multi mode (spectral) (todo)
-  !
-  !> @param[inout] rho charge array
-  !> @param[in] rmin tile grid minimum radius
-  !> @param[in] dr radial space discretization steps
-  !> @param[in] nx,ny,nz number of cells
-  !> @param[in] nxguard,nyguard,nzguard number of guard cells
-  !>
-  subroutine warpx_charge_deposition_rz_volume_scaling(rho,rho_ng,rho_ntot,rmin,dr) &
-    bind(C, name="warpx_charge_deposition_rz_volume_scaling")
-
-    integer, intent(in) :: rho_ntot(AMREX_SPACEDIM)
-    integer(c_long), intent(in) :: rho_ng
-    real(amrex_real), intent(IN OUT):: rho(*)
-    real(amrex_real), intent(IN) :: rmin, dr
-
-#ifdef WARPX_DIM_RZ
-    integer(c_long) :: type_rz_depose = 1
-#endif
-
-    ! Compute the number of valid cells and guard cells
-    integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM)
-    rho_nvalid = rho_ntot - 2*rho_ng
-    rho_nguards = rho_ng
-
-#ifdef WARPX_DIM_RZ
-    CALL WRPX_PXR_RZ_VOLUME_SCALING_RHO(   &
-                 rho,rho_nguards,rho_nvalid, &
-                 rmin,dr,type_rz_depose)
-#endif
-
-  end subroutine warpx_charge_deposition_rz_volume_scaling
-
   ! _________________________________________________________________
   !>
   !> @brief

Source/Laser/LaserParticleContainer.cpp

@@ -453,7 +453,12 @@ LaserParticleContainer::Evolve (int lev,
             pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]);
             BL_PROFILE_VAR_STOP(blp_copy);
 
-            if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev);
+            if (rho) {
+                DepositCharge(pti, wp, rho, 0, 0, np_current, thread_num, lev, lev);
+                if (crho) {
+                    DepositCharge(pti, wp, crho, 0, np_current, np-np_current, thread_num, lev, lev-1);
+                }
+            }
 
             //
             // Particle Push
@@ -522,7 +527,12 @@ LaserParticleContainer::Evolve (int lev,
             pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]);
             BL_PROFILE_VAR_STOP(blp_copy);
 
-            if (rho) DepositCharge(pti, wp, rho, crho, 1, np_current, np, thread_num, lev);
+            if (rho) {
+                DepositCharge(pti, wp, rho, 1, 0, np_current, thread_num, lev, lev);
+                if (crho) {
+                    DepositCharge(pti, wp, crho, 1, np_current, np-np_current, thread_num, lev, lev-1);
+                }
+            }
 
             if (cost) {
                 const Box& tbx = pti.tilebox();

Source/Particles/Deposition/ChargeDeposition.H

@@ -0,0 +1,97 @@
+#ifndef CHARGEDEPOSITION_H_
+#define CHARGEDEPOSITION_H_
+
+#include "ShapeFactors.H"
+
+/* \brief Charge Deposition for thread thread_num
+ * /param xp, yp, zp   : Pointer to arrays of particle positions.
+ * \param wp           : Pointer to array of particle weights.
+ * \param rho_arr      : Array4 of charge density, either full array or tile.
+ * \param np_to_depose : 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 q            : species charge.
+ */
+template <int depos_order>
+void doChargeDepositionShapeN(const amrex::Real * const xp, 
+                              const amrex::Real * const yp, 
+                              const amrex::Real * const zp,
+                              const amrex::Real * const wp,
+                              const amrex::Array4<amrex::Real>& rho_arr,
+                              const long np_to_depose,
+                              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 dxi = 1.0/dx[0];
+    const amrex::Real dzi = 1.0/dx[2];
+#if (AMREX_SPACEDIM == 2)
+    const amrex::Real invvol = dxi*dzi;
+#elif (defined WARPX_DIM_3D)
+    const amrex::Real dyi = 1.0/dx[1];
+    const amrex::Real invvol = dxi*dyi*dzi;
+#endif
+
+    const amrex::Real xmin = xyzmin[0];
+    const amrex::Real ymin = xyzmin[1];
+    const amrex::Real zmin = xyzmin[2];
+
+    // Loop over particles and deposit into rho_arr
+    amrex::ParallelFor(
+        np_to_depose,
+        [=] AMREX_GPU_DEVICE (long ip) {
+            // --- Get particle quantities
+            const amrex::Real wq = q*wp[ip]*invvol;
+
+            // --- Compute shape factors
+            // x direction
+            // Get particle position in grid coordinates
+#if (defined WARPX_DIM_RZ)
+            const amrex::Real r = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]);
+            const amrex::Real x = (r - xmin)*dxi;
+#else
+            const amrex::Real x = (xp[ip] - xmin)*dxi;
+#endif
+            // Compute shape factors for node-centered quantities
+            amrex::Real AMREX_RESTRICT sx[depos_order + 1];
+            // i: leftmost grid point (node-centered) that the particle touches
+            const int i = compute_shape_factor<depos_order>(sx,  x);
+                     
+#if (defined WARPX_DIM_3D)
+            // y direction
+            const amrex::Real y = (yp[ip] - ymin)*dyi;
+            amrex::Real AMREX_RESTRICT sy[depos_order + 1];
+            const int j = compute_shape_factor<depos_order>(sy,  y);
+#endif
+            // z direction
+            const amrex::Real z = (zp[ip] - zmin)*dzi;
+            amrex::Real AMREX_RESTRICT sz[depos_order + 1];
+            const int k = compute_shape_factor<depos_order>(sz,  z);
+
+            // Deposit charge into rho_arr
+#if (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ)
+            for (int iz=0; iz<=depos_order; iz++){
+                for (int ix=0; ix<=depos_order; ix++){
+                    amrex::Gpu::Atomic::Add(
+                        &rho_arr(lo.x+i+ix, lo.y+k+iz, 0), 
+                        sx[ix]*sz[iz]*wq);
+                }
+            }
+#elif (defined WARPX_DIM_3D)
+            for (int iz=0; iz<=depos_order; iz++){
+                for (int iy=0; iy<=depos_order; iy++){
+                    for (int ix=0; ix<=depos_order; ix++){
+                        amrex::Gpu::Atomic::Add(
+                            &rho_arr(lo.x+i+ix, lo.y+j+iy, lo.z+k+iz),
+                            sx[ix]*sy[iy]*sz[iz]*wq);
+                    }
+                }
+            }
+#endif
+        }
+        );
+}
+
+#endif // CHARGEDEPOSITION_H_

Source/Particles/Deposition/Make.package

@@ -1,3 +1,4 @@
 CEXE_headers += CurrentDeposition.H
+CEXE_headers += ChargeDeposition.H
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Deposition
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/Deposition