diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index 5ce876eb83..f3ce00abea 100644 --- a/.github/workflows/linux.yml +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 6d71602622..10974bc452 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/cmake/HiPACEFunctions.cmake b/cmake/HiPACEFunctions.cmake index e4f07dc998..98af4a93d1 100644 --- a/cmake/HiPACEFunctions.cmake +++ b/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() diff --git a/docs/source/building/building.rst b/docs/source/building/building.rst index 2b6856829a..00e1ddb1eb 100644 --- a/docs/source/building/building.rst +++ b/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 `__. diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 8e037b0ada..3bb6f3247f 100644 --- a/src/Hipace.cpp +++ b/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 10. && m_predcorr_B_error_tolerance > 0.) { amrex::Print() << "WARNING: Predictor corrector loop may have diverged!\n" diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index 5d1466f73a..2e1c65cdfc 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/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 ); diff --git a/src/particles/deposition/PlasmaDepositCurrent.cpp b/src/particles/deposition/PlasmaDepositCurrent.cpp index 418bad9e8f..e9037d1e17 100644 --- a/src/particles/deposition/PlasmaDepositCurrent.cpp +++ b/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 gpu_n_qsa_violation{}; int* AMREX_RESTRICT p_n_qsa_violation = nullptr; @@ -97,6 +94,9 @@ 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."); @@ -104,13 +104,9 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, amrex::AnyCTO( // use compile-time options amrex::TypeList< - amrex::CompileTimeOptions<0, 1, 2, 3>, // depos_order - amrex::CompileTimeOptions, // can_ionize - amrex::CompileTimeOptions // 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 @@ -118,28 +114,17 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, [&](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( - int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{aabs}, - amrex::GpuArray{jx, jy, jz, rho, chi, rhomjz, n}); - } else { - SharedMemoryDeposition( - int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{}, - amrex::GpuArray{jx, jy, jz, rho, chi, rhomjz, n}); - } + SharedMemoryDeposition( + int(pti.numParticles()), is_valid, get_cell, deposit, isl_fab.array(), + isl_fab.box(), pti.GetParticleTile().getParticleTileData(), + amrex::GpuArray{}, + amrex::GpuArray{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(xmid, 0); - auto [shape_y, j] = shape_factor(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 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(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 diff --git a/src/particles/deposition/TemperatureDeposition.cpp b/src/particles/deposition/TemperatureDeposition.cpp index 603c21f82e..b1e1e8415d 100644 --- a/src/particles/deposition/TemperatureDeposition.cpp +++ b/src/particles/deposition/TemperatureDeposition.cpp @@ -51,23 +51,15 @@ DepositTemperature (PlasmaParticleContainer& plasma, // Extract box properties const amrex::Real dx_inv = gm[lev].InvCellSize(0); const amrex::Real dy_inv = gm[lev].InvCellSize(1); - // extract laser properties and boolean for the presence of the laser and for ionization - const PhysConst pc = get_phys_const(); - const int aabs = Hipace::m_use_laser ? Comps[WhichSlice::This]["aabs"] : -1; - const amrex::Real laser_norm = (plasma.m_charge/pc.q_e) * (pc.m_e/plasma.m_mass) - * (plasma.m_charge/pc.q_e) * (pc.m_e/plasma.m_mass); + const bool use_laser = Hipace::m_use_laser; // Loop over particles amrex::AnyCTO( // use compile-time options amrex::TypeList< - amrex::CompileTimeOptions<0, 1, 2, 3>, // depos_order - amrex::CompileTimeOptions, // can_ionize - amrex::CompileTimeOptions // use_laser + amrex::CompileTimeOptions<0, 1, 2, 3> // depos_order >{}, { - Hipace::m_temperature_depos_order, - plasma.m_can_ionize, - Hipace::m_use_laser + Hipace::m_temperature_depos_order }, // call deposition function // The three functions passed as arguments to this lambda @@ -75,29 +67,18 @@ DepositTemperature (PlasmaParticleContainer& plasma, [&](auto is_valid, auto get_start_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( - int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{aabs}, - amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); - } else { - SharedMemoryDeposition( - int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), - isl_fab.box(), pti.GetParticleTile().getParticleTileData(), - amrex::GpuArray{}, - amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); - } + SharedMemoryDeposition( + int(pti.numParticles()), is_valid, get_start_cell, deposit, isl_fab.array(), + isl_fab.box(), pti.GetParticleTile().getParticleTileData(), + amrex::GpuArray{}, + amrex::GpuArray{w, ux, uy, uz, uxsq, uysq, uzsq}); }, // 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 on or below their according MR level return ptd.id(ip).is_valid() && (lev == 0 || ptd.cpu(ip) >= lev); @@ -105,9 +86,7 @@ DepositTemperature (PlasmaParticleContainer& plasma, // get_start_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); @@ -118,30 +97,21 @@ DepositTemperature (PlasmaParticleContainer& plasma, auto [shape_x, i] = shape_factor(xmid, 0); auto [shape_y, j] = shape_factor(ymid, 0); - return {i-1, j-1}; + return {i, j}; }, // do_deposit // deposit of weight, momentum (ux, uy, uz) and their squares (uxsq, uysq, uzsq) [=] AMREX_GPU_DEVICE (int ip, auto ptd, Array3 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 xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); amrex::Real Aabssqp = 0._rt; if (use_laser) { - amrex::Real laser_norm_ion = laser_norm; - if (can_ionize) { - laser_norm_ion *= - ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; - } - doLaserGatherShapeN<2>(xp, yp, Aabssqp, arr, cache_idx[0], - dx_inv, dy_inv, x_pos_offset, y_pos_offset); - Aabssqp *= laser_norm_ion; + Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; } const amrex::Real uxp = ptd.rdata(PlasmaIdx::ux)[ip]; diff --git a/src/particles/particles_utils/FieldGather.H b/src/particles/particles_utils/FieldGather.H index 40d9c39a9a..462aa96558 100644 --- a/src/particles/particles_utils/FieldGather.H +++ b/src/particles/particles_utils/FieldGather.H @@ -251,8 +251,7 @@ void doLaserGatherShapeN (const amrex::Real xp, auto [shape_y, j] = shape_factor(y, iy); auto [shape_x, i] = shape_factor(x, ix); - const amrex::Real x00y00 = slice_arr(i, j, aabs_comp); - Aabssqp += shape_x * shape_y * x00y00; + Aabssqp += shape_x * shape_y * slice_arr(i, j, aabs_comp); } } } diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index 74ebb3ffa6..897c1aa2bc 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -83,7 +83,7 @@ public: /** \brief Gather field values and push particles * - * \param[in,out] fields the general field class, modified by this function + * \param[in] fields the general field class * \param[in] gm Geometry of the simulation, to get the cell size etc. * \param[in] temp_slice if true, the temporary data (x_temp, ...) will be used * \param[in] lev MR level @@ -93,6 +93,17 @@ public: const Fields & fields, amrex::Vector const& gm, bool temp_slice, int lev, int const current_N_level); + /** \brief Reset positions after push to temp slice */ + void ResetPositions (); + + /** \brief Gather Laser field to the particles + * + * \param[in] lev MR level + * \param[in] gm Geometry of the simulation, to get the cell size etc. + * \param[in] fields the general field class + */ + void GatherLaser (int lev, amrex::Geometry const& gm, const Fields & fields); + /** \brief Loop over plasma species and deposit their neutralizing background, if needed * * \param[in,out] fields the general field class, modified by this function diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index f685f10ed1..f6816a21cd 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -133,6 +133,22 @@ MultiPlasma::AdvanceParticles ( } } +void +MultiPlasma::ResetPositions () +{ + for (int i=0; i @@ -184,6 +203,19 @@ public: const MultiLaser& laser, const amrex::Real background_density_SI); + /** Reset positions after push to temp slice */ + void ResetPositions (); + + /** Gather Laser field to the particles + * + * \param[in] lev MR level + * \param[in] geom Geometry of the simulation, to get the cell size + * \param[in] fields the general field class + */ + void GatherLaser (const int lev, + const amrex::Geometry& geom, + const Fields& fields); + /** Reorder particles to speed-up current deposition * \param[in] islice zeta slice index */ @@ -272,6 +304,12 @@ public: int m_n_subcycles = 1; /**< number of subcycles in the plasma particle push */ bool m_do_push = true; /**< whether the plasma pusher is enabled */ int m_ab5_permutation = 0; /**< index offset of the first ab5 force term */ + /** whether to use the AB5 or leapfrog pusher */ + bool m_use_ab5_push = false; + PlasmaComps m_comps; + bool m_is_on_temp_slice = false; + /** If the runtime attributes of the AMReX particle container were already set up */ + bool m_components_allocated = false; // ionization: @@ -322,8 +360,6 @@ private: amrex::Vector m_insitu_sum_idata; /** Prefix/path for the output files */ std::string m_insitu_file_prefix = ""; - /** If the runtime attributes of the AMReX particle container were already set up */ - bool m_components_allocated = false; }; /** \brief Iterator over boxes in a particle container */ diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index c7e2ab6b52..b81a2b7c04 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -12,6 +12,7 @@ #include "utils/AtomicWeightTable.H" #include "utils/DeprecatedInput.H" #include "utils/GPUUtil.H" +#include "utils/OMPUtil.H" #include "utils/InsituUtil.H" #ifdef HIPACE_USE_OPENPMD # include @@ -49,14 +50,13 @@ PlasmaParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mass_Da != 0, "Unknown Element"); } + queryWithParserAlt(pp, "use_ab5_push", m_use_ab5_push, pp_alt); queryWithParserAlt(pp, "n_subcycles", m_n_subcycles, pp_alt); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_n_subcycles >= 1, - "n_subcycles must be larger or equal to 1 sub-cycle (default is 1)"); -#ifdef HIPACE_USE_AB5_PUSH - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_n_subcycles == 1, - "Plasma subcycling only implemeted for leapfrog pusher!" - "Please set plasmas.n_subcycles = 1"); -#endif + "n_subcycles must be larger or equal to 1 sub-cycle (default is 1)"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_use_ab5_push || m_n_subcycles == 1, + "Plasma subcycling only implemeted for leapfrog pusher!" + "Please set plasmas.n_subcycles = 1"); queryWithParser(pp, "mass_Da", mass_Da); if(mass_Da != 0) { m_mass = phys_const.m_p * mass_Da / 1.007276466621; @@ -187,12 +187,46 @@ PlasmaParticleContainer::InitData (const amrex::Vector& geom3d) SetArena(amrex::The_Arena()); + m_comps.use_laser = Hipace::m_use_laser; + m_comps.use_temp_slice = Hipace::GetInstance().m_multi_beam.AnySpeciesSalame() + || !Hipace::m_explicit; + m_comps.use_ab5_push = m_use_ab5_push; + m_comps.use_ion_level = m_can_ionize; + + int num_real_comps = 0; + for (int j = 0; j < PlasmaIdx::real_nattribs; ++j) { AddRealComp(); + ++num_real_comps; } - for (int j = 0; j < PlasmaIdx::int_nattribs; ++j) { - AddIntComp(); + if (m_comps.use_laser) { + for (int j = 0; j < PlasmaIdx::real_nattribs_laser; ++j) { + AddRealComp(); + ++num_real_comps; + } + } + + if (m_comps.use_temp_slice) { + m_comps.offset_temp = num_real_comps; + for (int j = 0; j < PlasmaIdx::real_nattribs_temp_slice; ++j) { + AddRealComp(); + ++num_real_comps; + } + } + + if (m_comps.use_ab5_push) { + m_comps.offset_ab5 = num_real_comps; + for (int j = 0; j < PlasmaIdx::real_nattribs_ab5_push; ++j) { + AddRealComp(); + ++num_real_comps; + } + } + + if (m_comps.use_ion_level) { + for (int j = 0; j < PlasmaIdx::int_nattribs_ion_level; ++j) { + AddIntComp(); + } } reserveData(); @@ -345,9 +379,11 @@ PlasmaParticleContainer::TagByLevel (const int current_N_level, { auto& ptile = pti.GetParticleTile(); const amrex::Real * const AMREX_RESTRICT pos_x = to_prev ? - ptile.GetRealData(PlasmaIdx::x_prev).data() : ptile.GetRealData(PlasmaIdx::x).data(); + ptile.GetRealData(PlasmaIdx::x_prev + m_comps.offset_temp).data() + : ptile.GetRealData(PlasmaIdx::x).data(); const amrex::Real * const AMREX_RESTRICT pos_y = to_prev ? - ptile.GetRealData(PlasmaIdx::y_prev).data() : ptile.GetRealData(PlasmaIdx::y).data(); + ptile.GetRealData(PlasmaIdx::y_prev + m_comps.offset_temp).data() + : ptile.GetRealData(PlasmaIdx::y).data(); auto * AMREX_RESTRICT idcpup = ptile.GetIdCPUData().data(); const int lev1_idx = std::min(1, current_N_level-1); @@ -438,7 +474,7 @@ IonizationModule (const int lev, const int max_ion_lev = m_max_ion_lev; long num_ions = ptile_ion.numParticles(); - + const bool use_laser = Hipace::m_use_laser; // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time // and calculates ionization probability. If ionization occurs, it increments @@ -459,9 +495,8 @@ IonizationModule (const int lev, if (!ptd_ion.id(ip).is_valid() || ptd_ion.cpu(ip) != lev) return; - // Avoid temp slice - const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y)[ip]; // Define field at particle position reals amrex::ParticleReal ExmByp = 0., EypBxp = 0., Ezp = 0.; @@ -481,8 +516,11 @@ IonizationModule (const int lev, const amrex::Real psi = ptd_ion.rdata(PlasmaIdx::psi_half_step)[ip]; // Compute probability of ionization p - const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, - /* Assumes Aabssq == 0 */ 0._rt); + amrex::Real Aabssqp = 0._rt; + if (use_laser) { + Aabssqp = ptd_ion.rdata(PlasmaIdx::aabssq)[ip]; + } + const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, Aabssqp); const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip]; if (ion_lev_loc >= max_ion_lev) { return; @@ -519,43 +557,53 @@ IonizationModule (const int lev, // Load electron after resize auto ptd_elec = ptile_elec.getParticleTileData(); - const int init_ion_lev = m_product_pc->m_init_ion_lev; - amrex::Gpu::DeviceScalar ip_elec(0); uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); + auto comps = m_product_pc->m_comps; + AMREX_ALWAYS_ASSERT(comps.use_temp_slice == m_comps.use_temp_slice); + AMREX_ALWAYS_ASSERT(comps.offset_temp == m_comps.offset_temp); + AMREX_ALWAYS_ASSERT(comps.use_ion_level == false); + // This kernel adds the new ionized electrons to the Plasma Particle Container amrex::ParallelFor(num_ions, [=] AMREX_GPU_DEVICE (long ip) { - if(p_ion_mask[ip] != 0) { + if (p_ion_mask[ip] != 0) { const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` const long pidx = pid + old_size; // Copy ion data to new electron // Set the ionized electron ID to 2 (valid/invalid) for the ionized electrons + // Later we could consider adding a finite temperature to the ionized electrons ptd_elec.id(pidx) = 2; ptd_elec.cpu(pidx) = lev; // current level - ptd_elec.rdata(PlasmaIdx::x )[pidx] = ptd_ion.rdata(PlasmaIdx::x)[ip]; - ptd_elec.rdata(PlasmaIdx::y )[pidx] = ptd_ion.rdata(PlasmaIdx::y)[ip]; - - ptd_elec.rdata(PlasmaIdx::w )[pidx] = ptd_ion.rdata(PlasmaIdx::w)[ip]; - ptd_elec.rdata(PlasmaIdx::ux )[pidx] = 0._rt; - ptd_elec.rdata(PlasmaIdx::uy )[pidx] = 0._rt; - // Later we could consider adding a finite temperature to the ionized electrons - ptd_elec.rdata(PlasmaIdx::psi )[pidx] = 1._rt; // Assumes Aabssq == 0 - ptd_elec.rdata(PlasmaIdx::x_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - ptd_elec.rdata(PlasmaIdx::y_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + ptd_elec.rdata(PlasmaIdx::x)[pidx] = ptd_ion.rdata(PlasmaIdx::x)[ip]; + ptd_elec.rdata(PlasmaIdx::y)[pidx] = ptd_ion.rdata(PlasmaIdx::y)[ip]; + ptd_elec.rdata(PlasmaIdx::w)[pidx] = ptd_ion.rdata(PlasmaIdx::w)[ip]; + ptd_elec.rdata(PlasmaIdx::ux)[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::uy)[pidx] = 0._rt; + ptd_elec.rdata(PlasmaIdx::psi)[pidx] = 1._rt; // Assumes Aabssq == 0 ptd_elec.rdata(PlasmaIdx::ux_half_step )[pidx] = 0._rt; ptd_elec.rdata(PlasmaIdx::uy_half_step )[pidx] = 0._rt; ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = 1._rt; -#ifdef HIPACE_USE_AB5_PUSH - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd_elec.rdata(iforce)[pidx] = 0._rt; + + if (comps.use_laser) { + ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; + } + + if (comps.use_temp_slice) { + ptd_elec.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + ptd_elec.rdata(PlasmaIdx::y_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd_elec.rdata(iforce + comps.offset_ab5)[pidx] = 0._rt; + } } -#endif - ptd_elec.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; } }); @@ -650,9 +698,8 @@ LaserIonization (const int islice, [=] AMREX_GPU_DEVICE (long ip, const amrex::RandomEngine& engine, auto depos_order_xy) { - // Avoid temp slice - const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y)[ip]; if (!ptd_ion.id(ip).is_valid() || !laser_bounds.contains(xp, yp)) return; @@ -675,9 +722,9 @@ LaserIonization (const int islice, const amrex::Real uy = ptd_ion.rdata(PlasmaIdx::uy_half_step)[ip]; const amrex::Real psi = ptd_ion.rdata(PlasmaIdx::psi_half_step)[ip]; + // const amrex::Real Aabssqp = ptd_ion.rdata(PlasmaIdx::aabssq)[ip]; // Compute probability of ionization p - const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, - /* Assumes Aabssq == 0 */ 0._rt); + const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, 0._rt); // TODO Add Aabssqp const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip]; if (ion_lev_loc >= max_ion_lev) { return; @@ -718,11 +765,14 @@ LaserIonization (const int islice, // Load electron after resize auto ptd_elec = ptile_elec.getParticleTileData(); - const int init_ion_lev = m_product_pc->m_init_ion_lev; - amrex::Gpu::DeviceScalar ip_elec(0); uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr(); + auto comps = m_product_pc->m_comps; + AMREX_ALWAYS_ASSERT(comps.use_temp_slice == m_comps.use_temp_slice); + AMREX_ALWAYS_ASSERT(comps.offset_temp == m_comps.offset_temp); + AMREX_ALWAYS_ASSERT(comps.use_ion_level == false); + // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time. // It calculates the momentum of ionized electrons based on equations (B8) and (B9) // from the F. Massimo (2020) article and equation (12) from the P. Tomassini (2021) article. @@ -740,11 +790,10 @@ LaserIonization (const int islice, [=] AMREX_GPU_DEVICE (long ip, const amrex::RandomEngine& engine, auto depos_order_xy) { - if(p_ion_mask[ip] != 0) { + if (p_ion_mask[ip] != 0) { - // Avoid temp slice - const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd_ion.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd_ion.rdata(PlasmaIdx::y)[ip]; if (!ptd_ion.id(ip).is_valid() || !laser_bounds.contains(xp, yp)) return; @@ -792,29 +841,38 @@ LaserIonization (const int islice, const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec` const long pidx = pid + old_size; + // No normalization required since these are always electrons const amrex::Real psi = plasma_psi(ux, uy, uz, amrex::abs(A*A)); // Copy ion data to new electron // Set the ionized electron ID to 2 (valid/invalid) for the ionized electrons ptd_elec.id(pidx) = 2; ptd_elec.cpu(pidx) = ptd_ion.cpu(ip); // current level - ptd_elec.rdata(PlasmaIdx::x )[pidx] = ptd_ion.rdata(PlasmaIdx::x)[ip]; - ptd_elec.rdata(PlasmaIdx::y )[pidx] = ptd_ion.rdata(PlasmaIdx::y)[ip]; - ptd_elec.rdata(PlasmaIdx::w )[pidx] = ptd_ion.rdata(PlasmaIdx::w)[ip]; - ptd_elec.rdata(PlasmaIdx::ux )[pidx] = ux; - ptd_elec.rdata(PlasmaIdx::uy )[pidx] = uy; - ptd_elec.rdata(PlasmaIdx::psi )[pidx] = psi; - ptd_elec.rdata(PlasmaIdx::x_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::x_prev)[ip]; - ptd_elec.rdata(PlasmaIdx::y_prev )[pidx] = ptd_ion.rdata(PlasmaIdx::y_prev)[ip]; + ptd_elec.rdata(PlasmaIdx::x)[pidx] = ptd_ion.rdata(PlasmaIdx::x)[ip]; + ptd_elec.rdata(PlasmaIdx::y)[pidx] = ptd_ion.rdata(PlasmaIdx::y)[ip]; + ptd_elec.rdata(PlasmaIdx::w)[pidx] = ptd_ion.rdata(PlasmaIdx::w)[ip]; + ptd_elec.rdata(PlasmaIdx::ux)[pidx] = ux; + ptd_elec.rdata(PlasmaIdx::uy)[pidx] = uy; + ptd_elec.rdata(PlasmaIdx::psi)[pidx] = psi; ptd_elec.rdata(PlasmaIdx::ux_half_step )[pidx] = ux; ptd_elec.rdata(PlasmaIdx::uy_half_step )[pidx] = uy; ptd_elec.rdata(PlasmaIdx::psi_half_step)[pidx] = psi; -#ifdef HIPACE_USE_AB5_PUSH - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd_elec.rdata(iforce)[pidx] = 0._rt; + + if (comps.use_laser) { + ptd_elec.rdata(PlasmaIdx::aabssq)[pidx] = amrex::abs(A*A); + } + + if (comps.use_temp_slice) { + ptd_elec.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + ptd_elec.rdata(PlasmaIdx::y_prev + comps.offset_temp)[pidx] = + ptd_ion.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd_elec.rdata(iforce + comps.offset_ab5)[pidx] = 0._rt; + } } -#endif - ptd_elec.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; } }); @@ -823,6 +881,103 @@ LaserIonization (const int islice, } } +void +PlasmaParticleContainer:: +ResetPositions () +{ + HIPACE_PROFILE("PlasmaParticleContainer::ResetPositions()"); + + AMREX_ALWAYS_ASSERT(m_is_on_temp_slice); + m_is_on_temp_slice = false; + + // Loop over particle boxes + for (PlasmaParticleIterator pti(*this); pti.isValid(); ++pti) + { + // loading the data + const auto ptd = pti.GetParticleTile().getParticleTileData(); + const auto comps = m_comps; + AMREX_ALWAYS_ASSERT(m_comps.use_temp_slice); + + // Use OMP ParallelFor to use multiple threads when running on CPU + omp::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (int ip) { + if (!ptd.id(ip).is_valid()) return; + + ptd.rdata(PlasmaIdx::x)[ip] = ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + ptd.rdata(PlasmaIdx::y)[ip] = ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; + }); + } +} + +void +PlasmaParticleContainer:: +GatherLaser (const int lev, + const amrex::Geometry& geom, + const Fields& fields) +{ + if (!Hipace::m_use_laser) return; + HIPACE_PROFILE("PlasmaParticleContainer::GatherLaser()"); + + using namespace amrex::literals; + + const PhysConst phys_const = get_phys_const(); + + // Loop over particle boxes + for (PlasmaParticleIterator pti(*this); pti.isValid(); ++pti) + { + // Extract field array from FabArray + const amrex::FArrayBox& slice_fab = fields.getSlices(lev)[pti]; + Array3 const slice_arr = slice_fab.const_array(); + const int aabs_comp = Comps[WhichSlice::This]["aabs"]; + + // Extract properties associated with physical size of the box + const amrex::Real dx_inv = geom.InvCellSize(0); + const amrex::Real dy_inv = geom.InvCellSize(1); + + // Offset for converting positions to indexes + amrex::Real const x_pos_offset = GetPosOffset(0, geom, slice_fab.box()); + const amrex::Real y_pos_offset = GetPosOffset(1, geom, slice_fab.box()); + + // loading the data + const auto ptd = pti.GetParticleTile().getParticleTileData(); + const bool can_ionize = m_can_ionize; + + const amrex::Real laser_norm_qm = (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass) + * (m_charge/phys_const.q_e) * (phys_const.m_e/m_mass); + + // Use OMP ParallelFor to use multiple threads when running on CPU + omp::ParallelFor( + amrex::TypeList< + amrex::CompileTimeOptions<0, 1, 2, 3> + >{}, { + Hipace::m_depos_order_xy + }, + pti.numParticles(), + [=] AMREX_GPU_DEVICE (int ip, auto depos_order) { + // only push plasma particles on their according MR level + if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; + + amrex::Real laser_norm_qm_ion = laser_norm_qm; + if (can_ionize) { + const amrex::Real p_ion_lev = amrex::Real(ptd.idata(PlasmaIdx::ion_lev)[ip]); + laser_norm_qm_ion *= p_ion_lev * p_ion_lev; + } + + const amrex::Real xp = ptd.rdata(PlasmaIdx::x)[ip]; + const amrex::Real yp = ptd.rdata(PlasmaIdx::y)[ip]; + + amrex::Real Aabssqp = 0._rt; + doLaserGatherShapeN(xp, yp, + Aabssqp, slice_arr, aabs_comp, + dx_inv, dy_inv, x_pos_offset, y_pos_offset); + Aabssqp *= laser_norm_qm_ion; + + ptd.rdata(PlasmaIdx::aabssq)[ip] = Aabssqp; + }); + } +} + void PlasmaParticleContainer::InSituComputeDiags (int islice) { @@ -840,22 +995,8 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) { // Loading the data const auto ptd = pti.GetParticleTile().getParticleTileData(); - amrex::Long const num_particles = pti.numParticles(); - - const PhysConst pc = get_phys_const(); const bool use_laser = Hipace::m_use_laser; - const amrex::Geometry& gm = Hipace::GetInstance().m_3D_geom[0]; - const int aabs_comp = Hipace::m_use_laser ? Comps[WhichSlice::This]["aabs"] : -1; - amrex::FArrayBox& isl_fab = Hipace::GetInstance().m_fields.getSlices(0)[pti]; - Array3 arr = isl_fab.array(); - const amrex::Real x_pos_offset = GetPosOffset(0, gm, isl_fab.box()); - const amrex::Real y_pos_offset = GetPosOffset(1, gm, isl_fab.box()); - const amrex::Real dx_inv = gm.InvCellSize(0); - const amrex::Real dy_inv = gm.InvCellSize(1); - const bool can_ionize = m_can_ionize; - const amrex::Real laser_norm = (m_charge/pc.q_e) * (pc.m_e/m_mass) - * (m_charge/pc.q_e) * (pc.m_e/m_mass); amrex::TypeMultiplier reduce_op; amrex::TypeMultiplier reduce_data(reduce_op); @@ -877,14 +1018,7 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) amrex::Real Aabssqp = 0._rt; if (use_laser) { - amrex::Real laser_norm_ion = laser_norm; - if (can_ionize) { - laser_norm_ion *= - ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; - } - doLaserGatherShapeN<2>(x, y, Aabssqp, arr, aabs_comp, - dx_inv, dy_inv, x_pos_offset, y_pos_offset); - Aabssqp *= laser_norm_ion; + Aabssqp = ptd.rdata(PlasmaIdx::aabssq)[ip]; } // Particle's Lorentz factor diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index 02b6e2ca1d..9d171287e3 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -298,6 +298,8 @@ InitParticles (const amrex::RealVect& a_u_std, AMREX_ALWAYS_ASSERT(total_non_mirrored_particles == current_size); + auto comps = m_comps; + amrex::ParallelForRNG(current_size, [=] AMREX_GPU_DEVICE (unsigned int pidx, const amrex::RandomEngine& engine) { const amrex::Real x = ptd.rdata(PlasmaIdx::x)[pidx]; @@ -308,20 +310,30 @@ InitParticles (const amrex::RealVect& a_u_std, ptd.rdata(PlasmaIdx::ux)[pidx] = u[0]; ptd.rdata(PlasmaIdx::uy)[pidx] = u[1]; - ptd.rdata(PlasmaIdx::psi)[pidx] = plasma_psi(u[0], u[1], u[2], - /* Assumes Aabssq == 0 */ 0._rt); - ptd.rdata(PlasmaIdx::x_prev)[pidx] = x; - ptd.rdata(PlasmaIdx::y_prev)[pidx] = y; + ptd.rdata(PlasmaIdx::psi)[pidx] = + plasma_psi(u[0], u[1], u[2], /* Assumes Aabssq == 0 */ 0._rt); ptd.rdata(PlasmaIdx::ux_half_step)[pidx] = u[0]; ptd.rdata(PlasmaIdx::uy_half_step)[pidx] = u[1]; ptd.rdata(PlasmaIdx::psi_half_step)[pidx] = ptd.rdata(PlasmaIdx::psi)[pidx]; -#ifdef HIPACE_USE_AB5_PUSH - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd.rdata(iforce)[pidx] = 0._rt; + + if (comps.use_laser) { + ptd.rdata(PlasmaIdx::aabssq)[pidx] = 0._rt; + } + + if (comps.use_temp_slice) { + ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[pidx] = x; + ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[pidx] = y; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd.rdata(iforce + comps.offset_ab5)[pidx] = 0._rt; + } + } + + if (comps.use_ion_level) { + ptd.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; } -#endif - ptd.idata(PlasmaIdx::ion_lev)[pidx] = init_ion_lev; }); if (m_do_symmetrize) { @@ -349,27 +361,38 @@ InitParticles (const amrex::RealVect& a_u_std, ptd.cpu(midx) = 0; // level 0 ptd.rdata(PlasmaIdx::x)[midx] = x_arr[imirror]; ptd.rdata(PlasmaIdx::y)[midx] = y_arr[imirror]; - ptd.rdata(PlasmaIdx::w)[midx] = ptd.rdata(PlasmaIdx::w)[pidx]; - ptd.rdata(PlasmaIdx::ux)[midx] = ptd.rdata(PlasmaIdx::ux)[pidx] * ux_arr[imirror]; - ptd.rdata(PlasmaIdx::uy)[midx] = ptd.rdata(PlasmaIdx::uy)[pidx] * uy_arr[imirror]; - ptd.rdata(PlasmaIdx::psi)[midx] = ptd.rdata(PlasmaIdx::psi)[pidx]; - ptd.rdata(PlasmaIdx::x_prev)[midx] = x_arr[imirror]; - ptd.rdata(PlasmaIdx::y_prev)[midx] = y_arr[imirror]; + ptd.rdata(PlasmaIdx::ux)[midx] = + ptd.rdata(PlasmaIdx::ux)[pidx] * ux_arr[imirror]; + ptd.rdata(PlasmaIdx::uy)[midx] = + ptd.rdata(PlasmaIdx::uy)[pidx] * uy_arr[imirror]; + ptd.rdata(PlasmaIdx::psi)[midx] = + ptd.rdata(PlasmaIdx::psi)[pidx]; ptd.rdata(PlasmaIdx::ux_half_step)[midx] = ptd.rdata(PlasmaIdx::ux_half_step)[pidx] * ux_arr[imirror]; ptd.rdata(PlasmaIdx::uy_half_step)[midx] = ptd.rdata(PlasmaIdx::uy_half_step)[pidx] * uy_arr[imirror]; ptd.rdata(PlasmaIdx::psi_half_step)[midx] = ptd.rdata(PlasmaIdx::psi_half_step)[pidx]; -#ifdef HIPACE_USE_AB5_PUSH - HIPACE_LOOP_UNROLL - for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { - ptd.rdata(iforce)[midx] = 0._rt; + + if (comps.use_laser) { + ptd.rdata(PlasmaIdx::aabssq)[midx] = 0._rt; } -#endif - ptd.idata(PlasmaIdx::ion_lev)[midx] = ptd.idata(PlasmaIdx::ion_lev)[pidx]; + if (comps.use_temp_slice) { + ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[midx] = x_arr[imirror]; + ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[midx] = y_arr[imirror]; + } + + if (comps.use_ab5_push) { + for (int iforce = PlasmaIdx::Fx1; iforce <= PlasmaIdx::Fpsi5; ++iforce) { + ptd.rdata(iforce + comps.offset_ab5)[midx] = 0._rt; + } + } + + if (comps.use_ion_level) { + ptd.idata(PlasmaIdx::ion_lev)[midx] = init_ion_lev; + } } }); } diff --git a/src/particles/pusher/PlasmaParticleAdvance.cpp b/src/particles/pusher/PlasmaParticleAdvance.cpp index 3c5256c2a5..cf246dc366 100644 --- a/src/particles/pusher/PlasmaParticleAdvance.cpp +++ b/src/particles/pusher/PlasmaParticleAdvance.cpp @@ -79,17 +79,26 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, const amrex::Real clight_inv = 1._rt/phys_const.c; const amrex::Real charge_mass_clight_ratio = plasma.m_charge/(plasma.m_mass * phys_const.c); + const auto comps = plasma.m_comps; + const bool read_from_prev = plasma.m_is_on_temp_slice; + AMREX_ALWAYS_ASSERT(!(read_from_prev || temp_slice) || comps.use_temp_slice); + // Use OMP ParallelFor to use multiple threads when running on CPU omp::ParallelFor( amrex::TypeList< amrex::CompileTimeOptions<0, 1, 2, 3>, + amrex::CompileTimeOptions, + amrex::CompileTimeOptions, amrex::CompileTimeOptions >{}, { Hipace::m_depos_order_xy, - Hipace::m_use_laser + Hipace::m_use_laser, + plasma.m_use_ab5_push, + n_subcycles > 1 }, int(pti.numParticles()), // int ParallelFor is 3-5% faster than amrex::Long version - [=] AMREX_GPU_DEVICE (int ip, auto depos_order, auto use_laser) { + [=] AMREX_GPU_DEVICE (int ip, auto depos_order, auto use_laser, + auto use_ab5_push, auto use_subcycling) { // only push plasma particles on their according MR level if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; @@ -101,98 +110,131 @@ AdvancePlasmaParticles (PlasmaParticleContainer& plasma, const Fields & fields, amrex::Real q_mass_clight_ratio = charge_mass_clight_ratio; amrex::Real laser_norm_ion = laser_norm; if (can_ionize) { - q_mass_clight_ratio *= ptd.idata(PlasmaIdx::ion_lev)[ip]; - laser_norm_ion *= - ptd.idata(PlasmaIdx::ion_lev)[ip] * ptd.idata(PlasmaIdx::ion_lev)[ip]; + const amrex::Real p_ion_lev = amrex::Real(ptd.idata(PlasmaIdx::ion_lev)[ip]); + q_mass_clight_ratio *= p_ion_lev; + laser_norm_ion *= p_ion_lev * p_ion_lev; + } + + amrex::Real xp = 0._rt; + amrex::Real yp = 0._rt; + + if (!read_from_prev) { + xp = ptd.rdata(PlasmaIdx::x)[ip]; + yp = ptd.rdata(PlasmaIdx::y)[ip]; + if (temp_slice) { + // first temp slice + ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip] = xp; + ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip] = yp; + } + } else { + xp = ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + yp = ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; } - for (int i = 0; i < n_subcycles; i++) { + amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; + amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; + amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; + + amrex::Real ux_half = ux; + amrex::Real uy_half = uy; + amrex::Real psi_half = psi; - amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev)[ip]; - amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev)[ip]; + for (int i = 0; i < (use_subcycling ? n_subcycles : 1); i++) { - if (lev == 0 || lev_bounds.contains(xp, yp)) { + if (i == 0 || lev == 0 || lev_bounds.contains(xp, yp)) { ExmByp = 0._rt, EypBxp = 0._rt, Ezp = 0._rt; Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt; - Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt; doGatherShapeN(xp, yp, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp, slice_arr, psi_comp, ez_comp, bx_comp, by_comp, bz_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); + ExmByp *= clight_inv; + EypBxp *= clight_inv; + Ezp *= clight_inv; + if (use_laser.value) { + Aabssqp = 0._rt, AabssqDxp = 0._rt, AabssqDyp = 0._rt; + doLaserGatherShapeN(xp, yp, Aabssqp, AabssqDxp, AabssqDyp, slice_arr, aabs_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); - } - ExmByp *= clight_inv; - EypBxp *= clight_inv; - Ezp *= clight_inv; - Aabssqp *= 0.5_rt * laser_norm_ion; - AabssqDxp *= 0.25_rt * laser_norm_ion; - AabssqDyp *= 0.25_rt * laser_norm_ion; + Aabssqp *= 0.5_rt * laser_norm_ion; + AabssqDxp *= 0.25_rt * laser_norm_ion; + AabssqDyp *= 0.25_rt * laser_norm_ion; + } } -#ifndef HIPACE_USE_AB5_PUSH + if (!use_ab5_push.value) { - constexpr int nsub = 4; - const amrex::Real sdz = dz/nsub; + constexpr int nsub = 4; + const amrex::Real sdz = dz/nsub; - amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; - amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; - amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; + // full push in momentum + // from t-1/2 to t+1/2 + // using the fields at t + for (int isub=0; isub= 5) { + p -= 5; + } + p += comps.offset_ab5; + xp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fx1 + p)[ip]; + yp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fy1 + p)[ip]; + ux += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fux1 + p)[ip]; + uy += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fuy1 + p)[ip]; + psi += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fpsi1 + p)[ip]; + } - ux += sdz*dz_ux + 0.5_rt*sdz*sdz*dz_ux_dual.epsilon; - uy += sdz*dz_uy + 0.5_rt*sdz*sdz*dz_uy_dual.epsilon; - psi += sdz*dz_psi + 0.5_rt*sdz*sdz*dz_psi_dual.epsilon; + if (enforceBC(ptd, ip, xp, yp, ux, uy, PlasmaIdx::w)) return; + ux_half = ux; + uy_half = uy; + psi_half = psi; } - ptd.rdata(PlasmaIdx::ux)[ip] = ux; - ptd.rdata(PlasmaIdx::uy)[ip] = uy; - ptd.rdata(PlasmaIdx::psi)[ip] = psi; -#else - amrex::Real ux = ptd.rdata(PlasmaIdx::ux_half_step)[ip]; - amrex::Real uy = ptd.rdata(PlasmaIdx::uy_half_step)[ip]; - amrex::Real psi = ptd.rdata(PlasmaIdx::psi_half_step)[ip]; - const amrex::Real psi_inv = 1._rt/psi; - - auto [dz_ux, dz_uy, dz_psi] = PlasmaMomentumPush( - ux, uy, psi_inv, ExmByp, EypBxp, Ezp, Bxp, Byp, Bzp, - Aabssqp, AabssqDxp, AabssqDyp, q_mass_clight_ratio); - - ptd.rdata(PlasmaIdx::Fx1 + ab5_permutation)[ip] = ux * psi_inv; - ptd.rdata(PlasmaIdx::Fy1 + ab5_permutation)[ip] = uy * psi_inv; - ptd.rdata(PlasmaIdx::Fux1 + ab5_permutation)[ip] = dz_ux; - ptd.rdata(PlasmaIdx::Fuy1 + ab5_permutation)[ip] = dz_uy; - ptd.rdata(PlasmaIdx::Fpsi1 + ab5_permutation)[ip] = dz_psi; - - const amrex::Real ab5_coeffs[5] = { - ( 1901._rt / 720._rt ) * dz, // a1 times dz - ( -1387._rt / 360._rt ) * dz, // a2 times dz - ( 109._rt / 30._rt ) * dz, // a3 times dz - ( -637._rt / 360._rt ) * dz, // a4 times dz - ( 251._rt / 720._rt ) * dz // a5 times dz - }; - - HIPACE_LOOP_UNROLL - for (int iab=0; iab<5; ++iab) { - int p = ab5_permutation + iab; - if (p >= 5) { - p -= 5; - } - xp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fx1 + p)[ip]; - yp += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fy1 + p)[ip]; - ux += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fux1 + p)[ip]; - uy += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fuy1 + p)[ip]; - psi += ab5_coeffs[iab] * ptd.rdata(PlasmaIdx::Fpsi1 + p)[ip]; - } + } // loop over subcycles - if (enforceBC(ptd, ip, xp, yp, ux, uy, PlasmaIdx::w)) return; - ptd.pos(0, ip) = xp; - ptd.pos(1, ip) = yp; - - if (!temp_slice) { - // update values of the last non temp slice - // the next push always starts from these - ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux; - ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy; - ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi; - ptd.rdata(PlasmaIdx::x_prev)[ip] = xp; - ptd.rdata(PlasmaIdx::y_prev)[ip] = yp; - } + ptd.pos(0, ip) = xp; + ptd.pos(1, ip) = yp; - ptd.rdata(PlasmaIdx::ux)[ip] = ux; - ptd.rdata(PlasmaIdx::uy)[ip] = uy; - ptd.rdata(PlasmaIdx::psi)[ip] = psi; -#endif - } // loop over subcycles + ptd.rdata(PlasmaIdx::ux)[ip] = ux; + ptd.rdata(PlasmaIdx::uy)[ip] = uy; + ptd.rdata(PlasmaIdx::psi)[ip] = psi; + + if (!temp_slice) { + // update values of the last non temp slice + // the next push always starts from these + ptd.rdata(PlasmaIdx::ux_half_step)[ip] = ux_half; + ptd.rdata(PlasmaIdx::uy_half_step)[ip] = uy_half; + ptd.rdata(PlasmaIdx::psi_half_step)[ip] = psi_half; + } }); } -#ifdef HIPACE_USE_AB5_PUSH - if (!temp_slice && lev == current_N_level - 1) { + if (lev == current_N_level - 1) { + plasma.m_is_on_temp_slice = temp_slice; + } + + if (plasma.m_use_ab5_push && !temp_slice && lev == current_N_level - 1) { plasma.m_ab5_permutation = (plasma.m_ab5_permutation + 4) % 5; } -#endif } diff --git a/src/salame/Salame.cpp b/src/salame/Salame.cpp index abaefe8b75..ad9c8b54ce 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -182,12 +182,14 @@ SalameModule (Hipace* hipace, const int n_iter, const bool do_advance, int& last hipace->ExplicitMGSolveBxBy(lev, WhichSlice::This); } - } - if (hipace->m_N_level > 1) { - // tag to prev slice for push - hipace->m_multi_plasma.TagByLevel(current_N_level, hipace->m_3D_geom, true); + if (hipace->m_N_level > 1) { + // tag to prev slice for next push + hipace->m_multi_plasma.TagByLevel(current_N_level, hipace->m_3D_geom, true); + } } + + hipace->m_multi_plasma.ResetPositions(); } @@ -234,11 +236,20 @@ SalameGetJxJyFromBxBy (Hipace* hipace, const int lev) amrex::MultiFab& slicemf = hipace->m_fields.getSlices(lev); -#ifdef HIPACE_USE_AB5_PUSH - const amrex::Real dz = ( 1901._rt / 720._rt ) * hipace->m_3D_geom[lev].CellSize(Direction::z); -#else - const amrex::Real dz = 1.5_rt * hipace->m_3D_geom[lev].CellSize(Direction::z); -#endif + bool use_ab5_pusher = false; + for (int i=0; im_multi_plasma.GetNPlasmas(); ++i) { + if (i == 0) { + use_ab5_pusher = hipace->m_multi_plasma.m_all_plasmas[i].m_use_ab5_push; + } else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + hipace->m_multi_plasma.m_all_plasmas[i].m_use_ab5_push == use_ab5_pusher, + "All plasmas must use the same pusher when using SALAME" + ); + } + } + + const amrex::Real dz = (use_ab5_pusher ? 1901._rt / 720._rt : 1.5_rt) + * hipace->m_3D_geom[lev].CellSize(Direction::z); for ( amrex::MFIter mfi(slicemf, DfltMfiTlng); mfi.isValid(); ++mfi ){ @@ -280,7 +291,8 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) const amrex::Real dx_inv = gm.InvCellSize(0); const amrex::Real dy_inv = gm.InvCellSize(1); - const amrex::Real dz = gm.CellSize(2); + const amrex::Real dz = (plasma.m_use_ab5_push ? 1901._rt / 720._rt : 1.5_rt) + * gm.CellSize(2); // Offset for converting positions to indexes amrex::Real const x_pos_offset = GetPosOffset(0, gm, slice_fab.box()); @@ -290,7 +302,8 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) const amrex::Real charge_mass_c_ratio = plasma.m_charge / (plasma.m_mass * get_phys_const().c); - const bool can_ionize = plasma.m_can_ionize; + auto comps = plasma.m_comps; + AMREX_ALWAYS_ASSERT(comps.use_temp_slice); omp::ParallelFor( amrex::TypeList>{}, @@ -300,8 +313,8 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) // only push plasma particles on their according MR level if (!ptd.id(ip).is_valid() || ptd.cpu(ip) != lev) return; - const amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev)[ip]; - const amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev)[ip]; + const amrex::Real xp = ptd.rdata(PlasmaIdx::x_prev + comps.offset_temp)[ip]; + const amrex::Real yp = ptd.rdata(PlasmaIdx::y_prev + comps.offset_temp)[ip]; amrex::Real Bxp = 0._rt; amrex::Real Byp = 0._rt; @@ -310,17 +323,12 @@ SalameOnlyAdvancePlasma (Hipace* hipace, const int lev) doBxByGatherShapeN(xp, yp, Bxp, Byp, slice_arr, bx_comp, by_comp, dx_inv, dy_inv, x_pos_offset, y_pos_offset); - const amrex::Real q_m_c_ratio = can_ionize ? + const amrex::Real q_m_c_ratio = comps.use_ion_level ? ptd.idata(PlasmaIdx::ion_lev)[ip] * charge_mass_c_ratio : charge_mass_c_ratio; -#ifdef HIPACE_USE_AB5_PUSH - ptd.rdata(PlasmaIdx::ux)[ip] = ( 1901._rt / 720._rt )*dz * q_m_c_ratio * Byp; - ptd.rdata(PlasmaIdx::uy)[ip] = -( 1901._rt / 720._rt )*dz * q_m_c_ratio * Bxp; -#else - ptd.rdata(PlasmaIdx::ux)[ip] = 1.5_rt*dz * q_m_c_ratio * Byp; - ptd.rdata(PlasmaIdx::uy)[ip] = -1.5_rt*dz * q_m_c_ratio * Bxp; -#endif + ptd.rdata(PlasmaIdx::ux)[ip] = dz * q_m_c_ratio * Byp; + ptd.rdata(PlasmaIdx::uy)[ip] = -dz * q_m_c_ratio * Bxp; }); } diff --git a/tests/beam_in_vacuum.SI.Serial.sh b/tests/beam_in_vacuum.SI.Serial.sh index d667aa783a..c9bd631385 100755 --- a/tests/beam_in_vacuum.SI.Serial.sh +++ b/tests/beam_in_vacuum.SI.Serial.sh @@ -30,6 +30,7 @@ $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \ hipace.tile_size = 8 \ hipace.depos_order_xy=0 \ diagnostic.field_data = all rho \ + plasmas.use_ab5_push = 1 \ hipace.file_prefix=$TEST_NAME # Compare the result with theory diff --git a/tests/beam_in_vacuum.normalized.Serial.sh b/tests/beam_in_vacuum.normalized.Serial.sh index bc339714a5..8ae53d8154 100755 --- a/tests/beam_in_vacuum.normalized.Serial.sh +++ b/tests/beam_in_vacuum.normalized.Serial.sh @@ -30,6 +30,7 @@ $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \ hipace.tile_size = 8 \ hipace.depos_order_xy=0 \ diagnostic.field_data = all rho \ + plasmas.use_ab5_push = 1 \ hipace.file_prefix=$TEST_NAME # Compare the result with theory diff --git a/tests/blowout_wake.Serial.sh b/tests/blowout_wake.Serial.sh index 0c9d46b408..94964bfe65 100755 --- a/tests/blowout_wake.Serial.sh +++ b/tests/blowout_wake.Serial.sh @@ -28,6 +28,7 @@ HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests # Run the simulation $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \ hipace.tile_size = 8 \ + plasmas.use_ab5_push = 1 \ hipace.file_prefix=$TEST_NAME # Compare the results with checksum benchmark diff --git a/tests/checksum/benchmarks_json/laser_ionization.1Rank.json b/tests/checksum/benchmarks_json/laser_ionization.1Rank.json index e2b263b9d6..ce0d5e1cb6 100644 --- a/tests/checksum/benchmarks_json/laser_ionization.1Rank.json +++ b/tests/checksum/benchmarks_json/laser_ionization.1Rank.json @@ -31,11 +31,11 @@ "uy^2_ion": 0.0, "uy_elec": 0.0, "uy_ion": 0.0, - "uz^2_elec": 1.6568753834647e-06, - "uz^2_ion": 1.5025085926483e-20, - "uz_elec": 0.15044101093153, - "uz_ion": 4.8989912126984e-09, - "w_elec": 483046484939.32, + "uz^2_elec": 1.6613402997284e-06, + "uz^2_ion": 1.5330025597749e-20, + "uz_elec": 0.15061530702838, + "uz_ion": 4.9506114763176e-09, + "w_elec": 483046486891.44, "w_ion": 1400000000000.1, "|a^2|": 0.46877495720552 }