Allocate Plasma Components Dynamically

.github/workflows/linux.yml

@@ -73,8 +73,7 @@ jobs:
         pip install -U -e ./tools
         cmake -S . -B build             \
             -DHiPACE_MPI=OFF            \
-            -DCMAKE_VERBOSE_MAKEFILE=ON \
-            -DHiPACE_PUSHER=AB5
+            -DCMAKE_VERBOSE_MAKEFILE=ON
         cmake --build build -j 2
     - name: Run Tests
       run: ctest --test-dir build --output-on-failure

CMakeLists.txt

@@ -77,13 +77,6 @@ option(HiPACE_MPI            "Multi-node support (message-passing)"    ON)
 option(HiPACE_OPENPMD        "openPMD I/O (HDF5, ADIOS)"               ON)
 option(HiPACE_openpmd_mpi    "parallel version of openPMD I/O"         ${HiPACE_MPI})
 
-set(HiPACE_PUSHER_VALUES LEAPFROG AB5)
-set(HiPACE_PUSHER LEAPFROG CACHE STRING "Plasma pusher (LEAPFROG/AB5)")
-set_property(CACHE HiPACE_PUSHER PROPERTY STRINGS ${HiPACE_PUSHER_VALUES})
-if(NOT HiPACE_PUSHER IN_LIST HiPACE_PUSHER_VALUES)
-    message(FATAL_ERROR "HiPACE_PUSHER (${HiPACE_PUSHER}) must be one of ${HiPACE_PUSHER_VALUES}")
-endif()
-
 set(HiPACE_PRECISION_VALUES SINGLE DOUBLE)
 set(HiPACE_PRECISION DOUBLE CACHE STRING "Floating point precision (SINGLE/DOUBLE)")
 set_property(CACHE HiPACE_PRECISION PROPERTY STRINGS ${HiPACE_PRECISION_VALUES})
@@ -170,10 +163,6 @@ if(HiPACE_OPENPMD)
     target_link_libraries(HiPACE PUBLIC openPMD::openPMD)
 endif()
 
-if(HiPACE_PUSHER STREQUAL "AB5")
-    target_compile_definitions(HiPACE PUBLIC HIPACE_USE_AB5_PUSH)
-endif()
-
 if(AMReX_LINEAR_SOLVERS)
     target_compile_definitions(HiPACE PUBLIC AMREX_USE_LINEAR_SOLVERS)
 endif()

cmake/HiPACEFunctions.cmake

@@ -194,12 +194,6 @@ function(set_hipace_binary_name)
         set_property(TARGET HiPACE APPEND_STRING PROPERTY OUTPUT_NAME ".DEBUG")
     endif()
 
-    if(HiPACE_PUSHER STREQUAL "AB5")
-        set_property(TARGET HiPACE APPEND_STRING PROPERTY OUTPUT_NAME ".AB5")
-    else()
-        set_property(TARGET HiPACE APPEND_STRING PROPERTY OUTPUT_NAME ".LF")
-    endif()
-
     # alias to the latest build, because using the full name is often confusing
     add_custom_command(TARGET HiPACE POST_BUILD
         COMMAND ${CMAKE_COMMAND} -E create_symlink
@@ -274,6 +268,5 @@ function(hipace_print_summary)
     message("    MPI: ${HiPACE_MPI}")
     message("    OPENPMD: ${HiPACE_OPENPMD}")
     message("    PRECISION: ${HiPACE_PRECISION}")
-    message("    PUSHER: ${HiPACE_PUSHER}")
     message("")
 endfunction()

docs/source/building/building.rst

@@ -190,7 +190,6 @@ or by providing arguments to the CMake call
  ``HiPACE_MPI``                **ON**/OFF                                Multi-node support (message-passing)
  ``HiPACE_PRECISION``          SINGLE/**DOUBLE**                         Floating point precision (single/double)
  ``HiPACE_OPENPMD``            **ON**/OFF                                openPMD I/O (HDF5, ADIOS2)
- ``HiPACE_PUSHER``             **LEAPFROG**/AB5                          Use leapfrog or fifth-order Adams-Bashforth plasma pusher
 =============================  ========================================  =========================================================
 
 HiPACE++ can be configured in further detail with options from AMReX, which are documented in the `AMReX manual <https://amrex-codes.github.io/amrex/docs_html/BuildingAMReX.html#customization-options>`__.

src/Hipace.cpp

@@ -293,11 +293,6 @@ Hipace::InitData ()
     amrex::Print() << "using CUDA version " << __CUDACC_VER_MAJOR__ << "." << __CUDACC_VER_MINOR__
                    << "." << __CUDACC_VER_BUILD__ << "\n";
 #endif
-#ifdef HIPACE_USE_AB5_PUSH
-    amrex::Print() << "using the Adams-Bashforth plasma particle pusher\n";
-#else
-    amrex::Print() << "using the leapfrog plasma particle pusher\n";
-#endif
 
     m_multi_laser.InitData();
 
@@ -713,6 +708,11 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st
     // write laser aabs into fields MultiFab
     m_multi_laser.UpdateLaserAabs(islice, current_N_level, m_fields, m_3D_geom);
 
+    // interpolate laser aabs to plasma particles
+    for (int lev=0; lev<current_N_level; ++lev) {
+        m_multi_plasma.GatherLaser(lev, m_3D_geom[lev], m_fields);
+    }
+
     // has to be after aabs writing
     m_multi_plasma.InSituComputeDiags(step, islice, m_physical_time, is_last_step);
 
@@ -1189,6 +1189,8 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice, const int current_N
         relative_Bfield_error_prev_iter = relative_Bfield_error;
     } // end of predictor corrector loop
 
+    m_multi_plasma.ResetPositions();
+
     if (relative_Bfield_error > 10. && m_predcorr_B_error_tolerance > 0.)
     {
         amrex::Print() << "WARNING: Predictor corrector loop may have diverged!\n"

src/particles/collisions/CoulombCollision.cpp

@@ -93,7 +93,8 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
             amrex::Real* const uy1 = ptile1.GetRealData(PlasmaIdx::uy_half_step).data();
             amrex::Real* const psi1 = ptile1.GetRealData(PlasmaIdx::psi_half_step).data();
             const amrex::Real* const w1 = ptile1.GetRealData(PlasmaIdx::w).data();
-            const int* const ion_lev1 = ptile1.GetIntData(PlasmaIdx::ion_lev).data();
+            const int* const ion_lev1 = species1.m_can_ionize ?
+                ptile1.GetIntData(PlasmaIdx::ion_lev).data() : nullptr;
             PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
             PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
             amrex::Real q1 = species1.GetCharge();
@@ -159,7 +160,8 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
             amrex::Real* const uy1 = ptile1.GetRealData(PlasmaIdx::uy_half_step).data();
             amrex::Real* const psi1 = ptile1.GetRealData(PlasmaIdx::psi_half_step).data();
             const amrex::Real* const w1 = ptile1.GetRealData(PlasmaIdx::w).data();
-            const int* const ion_lev1 = ptile1.GetIntData(PlasmaIdx::ion_lev).data();
+            const int* const ion_lev1 = species1.m_can_ionize ?
+                ptile1.GetIntData(PlasmaIdx::ion_lev).data() : nullptr;
             PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
             PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
             amrex::Real q1 = species1.GetCharge();
@@ -172,7 +174,8 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
             amrex::Real* const uy2 = ptile2.GetRealData(PlasmaIdx::uy_half_step).data();
             amrex::Real* const psi2= ptile2.GetRealData(PlasmaIdx::psi_half_step).data();
             const amrex::Real* const w2 = ptile2.GetRealData(PlasmaIdx::w).data();
-            const int* const ion_lev2 = ptile2.GetIntData(PlasmaIdx::ion_lev).data();
+            const int* const ion_lev2 = species2.m_can_ionize ?
+                ptile2.GetIntData(PlasmaIdx::ion_lev).data() : nullptr;
             PlasmaBins::index_type * const indices2 = bins2.permutationPtr();
             PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr();
             amrex::Real q2 = species2.GetCharge();
@@ -282,7 +285,8 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
         amrex::Real* const uy2 = ptile2.GetRealData(PlasmaIdx::uy_half_step).data();
         amrex::Real* const psi2= ptile2.GetRealData(PlasmaIdx::psi_half_step).data();
         const amrex::Real* const w2 = ptile2.GetRealData(PlasmaIdx::w).data();
-        const int* const ion_lev2 = ptile2.GetIntData(PlasmaIdx::ion_lev).data();
+        const int* const ion_lev2 = species2.m_can_ionize ?
+            ptile2.GetIntData(PlasmaIdx::ion_lev).data() : nullptr;
         PlasmaBins::index_type * const indices2 = bins2.permutationPtr();
         PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr();
         amrex::Real q2 = species2.GetCharge();
@@ -332,7 +336,7 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
                 ElasticCollisionPerez(
                     cell_start1, cell_stop1, cell_start2, cell_stop2,
                     indices1, indices2,
-                    ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev2, ion_lev2, // passing ion_lev2 for beam particles, will never be used
+                    ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, nullptr, ion_lev2, // passing nullptr for beam particles, will never be used
                     q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2,
                     dt, CoulombLog, inv_dV, clight, inv_c_SI, inv_c2_SI,
                     normalized_units, background_density_SI, false, true, engine );

src/particles/deposition/PlasmaDepositCurrent.cpp

@@ -64,7 +64,6 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
         const int    chi = deposit_chi    ? Comps[which_slice]["chi"]    : -1;
         const int rhomjz = deposit_rhomjz ? Comps[which_slice]["rhomjz"] : -1;
         const int      n = deposit_n      ? Comps[which_slice][n_str]    : -1;
-        const int   aabs = Hipace::m_use_laser ? Comps[WhichSlice::This]["aabs"] : -1;
 
         // Offset for converting positions to indexes
         const amrex::Real x_pos_offset = GetPosOffset(0, gm[lev], isl_fab.box());
@@ -84,8 +83,6 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
         const amrex::Real clight = pc.c;
         const amrex::Real charge_invvol = charge * invvol;
         const amrex::Real charge_mu0_mass_ratio = charge * pc.mu0 / mass;
-        const amrex::Real laser_norm = (charge/pc.q_e) * (pc.m_e/mass)
-                                     * (charge/pc.q_e) * (pc.m_e/mass);
 
         amrex::Gpu::DeviceScalar<int> gpu_n_qsa_violation{};
         int* AMREX_RESTRICT p_n_qsa_violation = nullptr;
@@ -97,49 +94,37 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
                 &n_qsa_violation, &n_qsa_violation + 1, p_n_qsa_violation);
         }
 
+        const bool use_laser = Hipace::m_use_laser;
+        const bool can_ionize = plasma.m_can_ionize;
+
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(isl_fab.box().ixType().cellCentered(),
             "jx, jy, jz, and rho must be cell centered in all directions.");
 
         // Loop over particles and deposit into jx_fab, jy_fab, jz_fab, and rho_fab
         amrex::AnyCTO(
             // use compile-time options
             amrex::TypeList<
-                amrex::CompileTimeOptions<0, 1, 2, 3>,  // depos_order
-                amrex::CompileTimeOptions<false, true>, // can_ionize
-                amrex::CompileTimeOptions<false, true>  // use_laser
+                amrex::CompileTimeOptions<0, 1, 2, 3>   // depos_order
             >{}, {
-                Hipace::m_depos_order_xy,
-                plasma.m_can_ionize,
-                Hipace::m_use_laser
+                Hipace::m_depos_order_xy
             },
             // call deposition function
             // The three functions passed as arguments to this lambda
             // are defined below as the next arguments.
             [&](auto is_valid, auto get_cell, auto deposit){
                 constexpr auto ctos = deposit.GetOptions();
                 constexpr int depos_order = ctos[0];
-                constexpr int use_laser = ctos[2];
                 constexpr int stencil_size = depos_order + 1;
-                if constexpr (use_laser) {
-                    SharedMemoryDeposition<stencil_size, stencil_size, true>(
-                        int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(),
-                        isl_fab.box(), pti.GetParticleTile().getParticleTileData(),
-                        amrex::GpuArray<int, 1>{aabs},
-                        amrex::GpuArray<int, 7>{jx, jy, jz, rho, chi, rhomjz, n});
-                } else {
-                    SharedMemoryDeposition<stencil_size, stencil_size, true>(
-                        int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(),
-                        isl_fab.box(), pti.GetParticleTile().getParticleTileData(),
-                        amrex::GpuArray<int, 0>{},
-                        amrex::GpuArray<int, 7>{jx, jy, jz, rho, chi, rhomjz, n});
-                }
+                SharedMemoryDeposition<stencil_size, stencil_size, true>(
+                    int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(),
+                    isl_fab.box(), pti.GetParticleTile().getParticleTileData(),
+                    amrex::GpuArray<int, 0>{},
+                    amrex::GpuArray<int, 7>{jx, jy, jz, rho, chi, rhomjz, n});
             },
             // is_valid
             // return whether the particle is valid and should deposit
             [=] AMREX_GPU_DEVICE (int ip, auto ptd,
-                                  auto /*depos_order*/,
-                                  auto can_ionize,
-                                  auto /*use_laser*/)
+                                  auto /*depos_order*/)
             {
                 // only deposit plasma currents on or below their according MR level
                 return ptd.id(ip).is_valid() &&
@@ -150,9 +135,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
             // get_cell
             // return the lowest cell index that the particle deposits into
             [=] AMREX_GPU_DEVICE (int ip, auto ptd,
-                                  auto depos_order,
-                                  auto /*can_ionize*/,
-                                  auto /*use_laser*/) -> amrex::IntVectND<2>
+                                  auto depos_order) -> amrex::IntVectND<2>
             {
                 const amrex::Real xp = ptd.pos(0, ip);
                 const amrex::Real yp = ptd.pos(1, ip);
@@ -161,7 +144,6 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
                 const amrex::Real ymid = (yp - y_pos_offset) * dy_inv;
 
                 auto [shape_x, i] = shape_factor<depos_order>(xmid, 0);
-
                 auto [shape_y, j] = shape_factor<depos_order>(ymid, 0);
 
                 return {i, j};
@@ -170,10 +152,8 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
             // deposit the charge / current of one particle
             [=] AMREX_GPU_DEVICE (int ip, auto ptd,
                                   Array3<amrex::Real> arr,
-                                  auto cache_idx, auto depos_idx,
-                                  auto depos_order,
-                                  auto can_ionize,
-                                  auto use_laser) noexcept
+                                  auto /*cache_idx*/, auto depos_idx,
+                                  auto depos_order) noexcept
             {
                 const amrex::Real psi_inv = 1._rt/ptd.rdata(PlasmaIdx::psi)[ip];
                 const amrex::Real xp = ptd.pos(0, ip);
@@ -184,22 +164,18 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
                 // calculate charge of the plasma particles
                 amrex::Real q_invvol = charge_invvol * w;
                 amrex::Real q_mu0_mass_ratio = charge_mu0_mass_ratio;
-                [[maybe_unused]] amrex::Real laser_norm_ion = laser_norm;
-                if constexpr (can_ionize) {
+                if (can_ionize) {
                     const amrex::Real p_ion_lev = amrex::Real(ptd.idata(PlasmaIdx::ion_lev)[ip]);
                     q_invvol *= p_ion_lev;
                     q_mu0_mass_ratio *= p_ion_lev;
-                    laser_norm_ion *= p_ion_lev * p_ion_lev;
                 }
 
                 const amrex::Real xmid = (xp - x_pos_offset) * dx_inv;
                 const amrex::Real ymid = (yp - y_pos_offset) * dy_inv;
 
                 amrex::Real Aabssqp = 0._rt;
-                if constexpr (use_laser) {
-                    doLaserGatherShapeN<depos_order>(xp, yp, Aabssqp, arr, cache_idx[0],
-                                                    dx_inv, dy_inv, x_pos_offset, y_pos_offset);
-                    Aabssqp *= laser_norm_ion;
+                if (use_laser) {
+                    Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip];
                 }
 
                 // calculate gamma/psi for plasma particles