add option to not deposit Jx and Jy of the beam (#431)

doc/source/run/parameters.rst

@@ -47,6 +47,10 @@ General parameters
       requires
     | `hipace.beam_injection_cr = 8`.
 
+* ``hipace.do_beam_jx_jy_deposition`` (`bool`) optional (default `1`)
+    Using the default, the beam deposits all currents `Jx`, `Jy`, `Jz`. Using
+    `hipace.do_beam_jx_jy_deposition = 0` disables the transverse current deposition of the beams.
+
 Predictor-corrector loop parameters
 -----------------------------------
 

src/Hipace.H

@@ -209,6 +209,8 @@ public:
     /** Mixing factor between the transverse B field iterations in the predictor corrector loop
      */
     static amrex::Real m_predcorr_B_mixing_factor;
+    /** Whether the beams deposit Jx and Jy */
+    static bool m_do_beam_jx_jy_deposition;
     /** Whether to call amrex::Gpu::synchronize() around all profiler region */
     static bool m_do_device_synchronize;
     /** How much the box is coarsened for beam injection, to avoid exceeding max int in cell count.

src/Hipace.cpp

@@ -33,6 +33,7 @@ int Hipace::m_depos_order_z = 0;
 amrex::Real Hipace::m_predcorr_B_error_tolerance = 4e-2;
 int Hipace::m_predcorr_max_iterations = 30;
 amrex::Real Hipace::m_predcorr_B_mixing_factor = 0.05;
+bool Hipace::m_do_beam_jx_jy_deposition = true;
 bool Hipace::m_do_device_synchronize = false;
 int Hipace::m_beam_injection_cr = 1;
 amrex::Real Hipace::m_external_ExmBy_slope = 0.;
@@ -85,6 +86,7 @@ Hipace::Hipace () :
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_numprocs_x*m_numprocs_y*m_numprocs_z
                                      == amrex::ParallelDescriptor::NProcs(),
                                      "Check hipace.numprocs_x and hipace.numprocs_y");
+    pph.query("do_beam_jx_jy_deposition", m_do_beam_jx_jy_deposition);
     pph.query("do_device_synchronize", m_do_device_synchronize);
     pph.query("external_ExmBy_slope", m_external_ExmBy_slope);
     pph.query("external_Ez_slope", m_external_Ez_slope);
@@ -388,7 +390,8 @@ Hipace::SolveOneSlice (int islice, int lev, const int ibox,
     m_fields.SolvePoissonExmByAndEypBx(Geom(lev), m_comm_xy, lev);
 
     m_grid_current.DepositCurrentSlice(m_fields, geom[lev], lev, islice);
-    m_multi_beam.DepositCurrentSlice(m_fields, geom[lev], lev, islice, bx, bins, m_box_sorters, ibox);
+    m_multi_beam.DepositCurrentSlice(m_fields, geom[lev], lev, islice, bx, bins, m_box_sorters,
+                                     ibox, m_do_beam_jx_jy_deposition);
 
     j_slice.FillBoundary(Geom(lev).periodicity());
 

src/particles/MultiBeam.H

@@ -29,11 +29,13 @@ public:
      * \param[in] bins Vector (over species) of particles sorted by slices
      * \param[in] a_box_sorter_vec Vector (over species) of particles sorted by box
      * \param[in] ibox index of the current box
+     * \param[in] do_beam_jx_jy_deposition whether the beam deposits Jx and Jy
      */
     void DepositCurrentSlice (
         Fields& fields, const amrex::Geometry& geom, const int lev, int islice, const amrex::Box bx,
         amrex::Vector<amrex::DenseBins<BeamParticleContainer::ParticleType>> bins,
-        const amrex::Vector<BoxSorter>& a_box_sorter_vec, const int ibox);
+        const amrex::Vector<BoxSorter>& a_box_sorter_vec, const int ibox,
+        const bool do_beam_jx_jy_deposition);
 
     /** Loop over all beam species and build and return the indices of particles sorted per slice
      * \param[in] lev MR level

src/particles/MultiBeam.cpp

@@ -27,12 +27,14 @@ void
 MultiBeam::DepositCurrentSlice (
     Fields& fields, const amrex::Geometry& geom, const int lev, int islice, const amrex::Box bx,
     amrex::Vector<amrex::DenseBins<BeamParticleContainer::ParticleType>> bins,
-    const amrex::Vector<BoxSorter>& a_box_sorter_vec, const int ibox)
+    const amrex::Vector<BoxSorter>& a_box_sorter_vec, const int ibox,
+    const bool do_beam_jx_jy_deposition)
 
 {
     for (int i=0; i<m_nbeams; i++) {
         ::DepositCurrentSlice(m_all_beams[i], fields, geom, lev, islice, bx,
-                              a_box_sorter_vec[i].boxOffsetsPtr()[ibox], bins[i]);
+                              a_box_sorter_vec[i].boxOffsetsPtr()[ibox], bins[i],
+                              do_beam_jx_jy_deposition);
     }
 }
 

src/particles/deposition/BeamDepositCurrent.H

@@ -16,11 +16,13 @@
  * \param[in] bx current box to calculate in loop over longutidinal boxes
  * \param[in] offset offset to the particles on the current box
  * \param[in] bins beam particle container bins, to push only the beam particles on slice islice
+ * \param[in] do_beam_jx_jy_deposition whether the beams deposit Jx and Jy
  */
 void
 DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields, amrex::Geometry const& gm,
                      int const lev, const int islice, const amrex::Box bx, int const offset,
-                     amrex::DenseBins<BeamParticleContainer::ParticleType>& bins);
+                     amrex::DenseBins<BeamParticleContainer::ParticleType>& bins,
+                     const bool do_beam_jx_jy_deposition);
 
 
 #endif //  BEAMDEPOSITCURRENT_H_

src/particles/deposition/BeamDepositCurrent.cpp

@@ -11,7 +11,8 @@
 void
 DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields, amrex::Geometry const& gm,
                      int const lev ,const int islice, const amrex::Box bx, int const offset,
-                     amrex::DenseBins<BeamParticleContainer::ParticleType>& bins)
+                     amrex::DenseBins<BeamParticleContainer::ParticleType>& bins,
+                     const bool do_beam_jx_jy_deposition)
 {
     HIPACE_PROFILE("DepositCurrentSlice_BeamParticleContainer()");
     // Extract properties associated with physical size of the box
@@ -51,17 +52,17 @@ DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields, amrex::Geometr
 
     // Call deposition function in each box
     if        (Hipace::m_depos_order_xy == 0){
-        doDepositionShapeN<0, 0>( beam, jx_fab, jy_fab, jz_fab,
-                                  dx, xyzmin, lo, q, islice_local, bins, offset);
+        doDepositionShapeN<0, 0>( beam, jx_fab, jy_fab, jz_fab, dx, xyzmin, lo, q, islice_local,
+                                  bins, offset, do_beam_jx_jy_deposition);
     } else if (Hipace::m_depos_order_xy == 1){
-        doDepositionShapeN<1, 0>( beam, jx_fab, jy_fab, jz_fab,
-                                  dx, xyzmin, lo, q, islice_local, bins, offset);
+        doDepositionShapeN<1, 0>( beam, jx_fab, jy_fab, jz_fab, dx, xyzmin, lo, q, islice_local,
+                                  bins, offset, do_beam_jx_jy_deposition);
     } else if (Hipace::m_depos_order_xy == 2){
-        doDepositionShapeN<2, 0>( beam, jx_fab, jy_fab, jz_fab,
-                                  dx, xyzmin, lo, q, islice_local, bins, offset);
+        doDepositionShapeN<2, 0>( beam, jx_fab, jy_fab, jz_fab, dx, xyzmin, lo, q, islice_local,
+                                  bins, offset, do_beam_jx_jy_deposition);
     } else if (Hipace::m_depos_order_xy == 3){
-        doDepositionShapeN<3, 0>( beam, jx_fab, jy_fab, jz_fab,
-                                  dx, xyzmin, lo, q, islice_local, bins, offset);
+        doDepositionShapeN<3, 0>( beam, jx_fab, jy_fab, jz_fab, dx, xyzmin, lo, q, islice_local,
+                                  bins, offset, do_beam_jx_jy_deposition);
     } else {
         amrex::Abort("unknown deposition order");
     }

src/particles/deposition/BeamDepositCurrentInner.H

@@ -35,6 +35,7 @@
  * \param[in] islice Particles in slice islice will deposit
  * \param[in] bins Indices of particles arranged per slices.
  * \param[in] box_offset offset to particles on this box.
+ * \param[in] do_beam_jx_jy_deposition whether the beams deposit Jx and Jy
  */
 template <int depos_order_xy, int depos_order_z>
 void doDepositionShapeN (const BeamParticleContainer& ptile,
@@ -47,7 +48,8 @@ void doDepositionShapeN (const BeamParticleContainer& ptile,
                          amrex::Real const q,
                          int islice,
                          amrex::DenseBins<BeamParticleContainer::ParticleType>& bins,
-                         int box_offset)
+                         int box_offset,
+                         const bool do_beam_jx_jy_deposition)
 {
     using namespace amrex::literals;
 
@@ -153,12 +155,14 @@ void doDepositionShapeN (const BeamParticleContainer& ptile,
             for (int iz=0; iz<=depos_order_z; iz++){
                 for (int iy=0; iy<=depos_order_xy; iy++){
                     for (int ix=0; ix<=depos_order_xy; ix++){
-                        amrex::Gpu::Atomic::Add(
-                            &jx_arr(lo.x+j_cell+ix, lo.y+k_cell+iy, z_slice),
-                            sx_cell[ix]*sy_cell[iy]*sz_cell[iz]*wqx);
-                        amrex::Gpu::Atomic::Add(
-                            &jy_arr(lo.x+j_cell+ix, lo.y+k_cell+iy, z_slice),
-                            sx_cell[ix]*sy_cell[iy]*sz_cell[iz]*wqy);
+                        if (do_beam_jx_jy_deposition) {
+                            amrex::Gpu::Atomic::Add(
+                                &jx_arr(lo.x+j_cell+ix, lo.y+k_cell+iy, z_slice),
+                                sx_cell[ix]*sy_cell[iy]*sz_cell[iz]*wqx);
+                            amrex::Gpu::Atomic::Add(
+                                &jy_arr(lo.x+j_cell+ix, lo.y+k_cell+iy, z_slice),
+                                sx_cell[ix]*sy_cell[iy]*sz_cell[iz]*wqy);
+                        }
                         amrex::Gpu::Atomic::Add(
                             &jz_arr(lo.x+j_cell+ix, lo.y+k_cell+iy, z_slice),
                             sx_cell[ix]*sy_cell[iy]*sz_cell[iz]*wqz);