From e58453c06f8058318eff110389f36078e8a0dd55 Mon Sep 17 00:00:00 2001 From: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com> Date: Thu, 6 Jul 2023 15:30:33 +0200 Subject: [PATCH] Add beam-plasma collisions (#983) * prepare functions * add basic framework * add transverse beam particle binning per slice * add beam collisions * fix doc * fix CI * fix automatic temperature calculation for beams * add benchmark CI test for collisions with beam * reset benchmark for new CI test * cleaning and doc * fix doc * Apply suggestions from code review * Update docs/source/run/parameters.rst * Update docs/source/run/parameters.rst --------- Co-authored-by: Tools --- CMakeLists.txt | 6 + docs/source/run/parameters.rst | 13 +- src/Hipace.H | 13 ++ src/Hipace.cpp | 49 +++++- src/particles/beam/BeamParticleContainer.H | 4 + src/particles/beam/MultiBeam.H | 4 +- src/particles/collisions/ComputeTemperature.H | 10 +- src/particles/collisions/CoulombCollision.H | 29 +++- src/particles/collisions/CoulombCollision.cpp | 150 ++++++++++++++++-- .../collisions/ElasticCollisionPerez.H | 28 ++-- src/particles/plasma/MultiPlasma.H | 19 +-- src/particles/plasma/MultiPlasma.cpp | 31 +--- src/particles/sorting/TileSort.H | 16 ++ src/particles/sorting/TileSort.cpp | 52 ++++++ .../collisions_beam.SI.1Rank.json | 32 ++++ tests/checksum/reset_all_benchmarks.sh | 12 ++ tests/collisions.SI.1Rank.sh | 2 +- tests/collisions_beam.SI.1Rank.sh | 40 +++++ 18 files changed, 423 insertions(+), 87 deletions(-) create mode 100644 tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json create mode 100755 tests/collisions_beam.SI.1Rank.sh diff --git a/CMakeLists.txt b/CMakeLists.txt index b7091e3622..a67e2db460 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -345,6 +345,12 @@ if(BUILD_TESTING) WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} ) + add_test(NAME collisions_beam.SI.1Rank + COMMAND bash ${HiPACE_SOURCE_DIR}/tests/collisions_beam.SI.1Rank.sh + $ ${HiPACE_SOURCE_DIR} + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} + ) + add_test(NAME ion_motion.SI.1Rank COMMAND bash ${HiPACE_SOURCE_DIR}/tests/ion_motion.SI.1Rank.sh $ ${HiPACE_SOURCE_DIR} diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index 66842a089f..dbb863e586 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -810,15 +810,18 @@ Binary collisions WARNING: this module is in development. -HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX, between plasma species. +HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX, +for collisions between plasma-plasma and beam-plasma. As collisions depend on the physical density, in normalized units `hipace.background_density_SI` must be specified. -* ``plasmas.collisions`` (list of `strings`) optional - List of names of types binary Coulomb collisions. - Each will represent collisions between 2 plasma species (potentially the same). +* ``hipace.collisions`` (list of `strings`) optional + List of names of binary Coulomb collisions. + Each will represent collisions between 2 species. * ``.species`` (two `strings`) optional - The name of the two plasma species for which collisions should be included. + The name of the two species for which collisions should be included. + This can either be plasma-plasma or beam-plasma collisions. For plasma-plasma collisions, the species can be the same to model collisions within a species. + The names must be in `plasmas.names` or `beams.names` (for beam-plasma collisions). * ``.CoulombLog`` (`float`) optional (default `-1.`) Coulomb logarithm used for this collision. diff --git a/src/Hipace.H b/src/Hipace.H index 29cfca7d45..367f3d1bb0 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -162,6 +162,12 @@ public: /** \brief reset all laser slices to 0 */ void ResetLaser (); + /** + * \brief does Coulomb collisions between plasmas and beams + * \param[in] islice_local slice on which collisions are calculated, needed for beam binning + */ + void doCoulombCollision (const int islice_local); + /** \brief Returns true on the head rank, otherwise false. */ @@ -310,6 +316,8 @@ public: inline static amrex::Real m_background_density_SI = 0.; /** Time step for the beam evolution */ amrex::Real m_dt = 0.0; + /** Number of binary collisions */ + inline static int m_ncollisions = 0; /** Adaptive time step instance */ AdaptiveTimeStep m_adaptive_time_step; /** Laser instance (soon to be multi laser container) */ @@ -368,6 +376,11 @@ private: /** Diagnostics */ Diagnostic m_diags; + /** User-input names of the binary collisions to be used */ + std::vector m_collision_names; + /** Vector of binary collisions */ + amrex::Vector< CoulombCollision > m_all_collisions; + /** \brief resizes the diagnostic fab to the correct box in a loop over boxes * * \param[in] it index of box to be resized to diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 79e06e3649..552288211e 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -170,6 +170,18 @@ Hipace::Hipace () : MakeGeometry(); m_use_laser = m_multi_laser.m_use_laser; + + queryWithParser(pph, "collisions", m_collision_names); + /** Initialize the collision objects */ + m_ncollisions = m_collision_names.size(); + for (int i = 0; i < m_ncollisions; ++i) { + m_all_collisions.emplace_back(CoulombCollision(m_multi_plasma.m_names, m_multi_beam.m_names, m_collision_names[i])); + } + if (m_normalized_units && m_ncollisions > 0) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_background_density_SI!=0, + "For collisions with normalized units, a background plasma density must " + "be specified via 'hipace.background_density_SI'"); + } } Hipace::~Hipace () @@ -618,8 +630,8 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step) m_external_E_uniform, m_external_B_uniform, m_external_E_slope, m_external_B_slope); - // collisions for all particles calculated on level 0 - m_multi_plasma.doCoulombCollision(0, m_slice_geom[0].Domain(), m_slice_geom[0]); + // collisions for plasmas and beams + doCoulombCollision(islice_local); // Advance laser slice by 1 step and store result to 3D array // no MR for laser @@ -1350,6 +1362,39 @@ Hipace::NotifyFinish (const int it, bool only_ghost, bool only_time) #endif } +void +Hipace::doCoulombCollision (const int islice_local) +{ + + // collisions for all particles calculated on level 0 + const int lev = 0; + + for (int i = 0; i < m_ncollisions; ++i) + { + if (m_all_collisions[i].m_nbeams == 1) { + // do beam-plasma collisions + auto& species1 = m_multi_beam.m_all_beams[ m_all_collisions[i].m_species1_index ]; + auto& species2 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species2_index ]; + + // TODO: enable tiling + + CoulombCollision::doBeamPlasmaCoulombCollision( + lev, m_slice_geom[0].Domain(), m_slice_geom[0], islice_local, species1, species2, + m_all_collisions[i].m_CoulombLog, m_background_density_SI); + } else { + // do plasma-plasma collisions + auto& species1 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species1_index ]; + auto& species2 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species2_index ]; + + // TODO: enable tiling + + CoulombCollision::doPlasmaPlasmaCoulombCollision( + lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2, m_all_collisions[i].m_isSameSpecies, + m_all_collisions[i].m_CoulombLog, m_background_density_SI); + } + } +} + void Hipace::ResizeFDiagFAB (const int it, const int step) { diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H index 0e87f192ed..82888fb7c7 100644 --- a/src/particles/beam/BeamParticleContainer.H +++ b/src/particles/beam/BeamParticleContainer.H @@ -157,6 +157,10 @@ public: std::string get_name () const {return m_name;} amrex::Real m_charge; /**< charge of each particle of this species */ amrex::Real m_mass; /**< mass of each particle of this species */ + /** Returns elementary charge q_e (or -q_e for electrons). */ + amrex::Real GetCharge () const {return m_charge;} + /** Returns mass of physical species */ + amrex::Real GetMass () const {return m_mass;} bool m_do_z_push {true}; /**< Pushing beam particles in z direction */ int m_n_subcycles {10}; /**< Number of sub-cycles in the beam pusher */ bool m_do_radiation_reaction {false}; /**< whether to calculate radiation losses */ diff --git a/src/particles/beam/MultiBeam.H b/src/particles/beam/MultiBeam.H index ebd1b5a4f9..13d0d2b2cf 100644 --- a/src/particles/beam/MultiBeam.H +++ b/src/particles/beam/MultiBeam.H @@ -193,11 +193,11 @@ public: /** \brief returns if any beam uses the SALAME algorithm */ bool AnySpeciesSalame (); + amrex::Vector m_names {"no_beam"}; /**< names of all beam containers */ + amrex::Vector m_all_beams; /**< contains all beam containers */ private: - amrex::Vector m_all_beams; /**< contains all beam containers */ - amrex::Vector m_names {"no_beam"}; /**< names of all beam containers */ int m_nbeams {0}; /**< number of beam containers */ /** number of real particles per beam, as opposed to ghost particles */ amrex::Vector m_n_real_particles; diff --git a/src/particles/collisions/ComputeTemperature.H b/src/particles/collisions/ComputeTemperature.H index 107ff2e027..413a82d475 100644 --- a/src/particles/collisions/ComputeTemperature.H +++ b/src/particles/collisions/ComputeTemperature.H @@ -13,7 +13,7 @@ AMREX_GPU_HOST_DEVICE T_R ComputeTemperature ( T_index const Is, T_index const Ie, T_index const * AMREX_RESTRICT I, T_R const * AMREX_RESTRICT ux, T_R const * AMREX_RESTRICT uy, T_R const * AMREX_RESTRICT psi, - T_R const m, T_R const clight, T_R const inv_c2 ) + T_R const m, T_R const clight, T_R const inv_c2, bool is_beam_coll ) { using namespace amrex::literals; @@ -27,9 +27,11 @@ T_R ComputeTemperature ( for (int i = Is; i < (int) Ie; ++i) { // particle's Lorentz factor - gm = (1.0_rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]])*inv_c2 - + psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] ); - const amrex::Real uz = clight * (gm - psi[I[i]]); + gm = is_beam_coll ? std::sqrt( 1._rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]] + + psi[I[i]]*psi[I[i]])*inv_c2 ) + : (1.0_rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]])*inv_c2 + + psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] ); + const amrex::Real uz = is_beam_coll ? psi[I[i]] : clight * (gm - psi[I[i]]); us = ( ux[ I[i] ] * ux[ I[i] ] + uy[ I[i] ] * uy[ I[i] ] + uz * uz); diff --git a/src/particles/collisions/CoulombCollision.H b/src/particles/collisions/CoulombCollision.H index 675bbca18b..4c38d85415 100644 --- a/src/particles/collisions/CoulombCollision.H +++ b/src/particles/collisions/CoulombCollision.H @@ -2,25 +2,29 @@ #define HIPACE_COULOMB_COLLISION_H_ #include "particles/plasma/PlasmaParticleContainer.H" +#include "particles/beam/BeamParticleContainer.H" #include #include #include /** - * \brief This class handles Coulomb collisions between 2 plasma species (can be the same). + * \brief This class handles Coulomb collisions between 2 particle species + * (can be plasma-plasma or beam-plasma, the species can be the same) */ class CoulombCollision { public: int m_species1_index {-1}; int m_species2_index {-1}; + int m_nbeams {0}; bool m_isSameSpecies {false}; amrex::Real m_CoulombLog {-1.}; /** Constructor */ CoulombCollision( - const std::vector& species_names, + const std::vector& plasma_species_names, + const std::vector& beam_species_names, std::string const collision_name); /** @@ -37,11 +41,30 @@ public: * Coulomb logarithm is deduced from the plasma temperature, measured in each cell. * \param[in] background_density_SI background plasma density (only needed for normalized units) **/ - static void doCoulombCollision ( + static void doPlasmaPlasmaCoulombCollision ( int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1, PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog, amrex::Real background_density_SI); + /** + * \brief Perform Coulomb collisions of a beam with a plasma species over a push by one beam time step + * Particles of both species are sorted per cell, paired, and collided pairwise. + * + * \param[in] lev MR level + * \param[in] bx transverse box (plasma particles will be sorted per-cell on this box) + * \param[in] geom corresponding geometry object + * \param[in] islice_local local slice on which collisions are calculated + * \param[in,out] species1 beam species + * \param[in,out] species2 plasma species + * \param[in] CoulombLog Value of the Coulomb logarithm used for the collisions. If <0, the + * Coulomb logarithm is deduced from the plasma temperature, measured in each cell. + * \param[in] background_density_SI background plasma density (only needed for normalized units) + **/ + static void doBeamPlasmaCoulombCollision ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, int islice_local, + BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, + amrex::Real background_density_SI); + }; #endif // HIPACE_COULOMB_COLLISION_H_ diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index dd905672a1..b71a4a52eb 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -6,7 +6,8 @@ #include "utils/HipaceProfilerWrapper.H" CoulombCollision::CoulombCollision( - const std::vector& species_names, + const std::vector& plasma_species_names, + const std::vector& beam_species_names, std::string const collision_name) { using namespace amrex::literals; @@ -24,20 +25,39 @@ CoulombCollision::CoulombCollision( // default Coulomb log is -1, if < 0 (e.g. not specified), will be computed automatically pp.query("CoulombLog", m_CoulombLog); - for (int i=0; i<(int) species_names.size(); i++) - { - if (species_names[i] == collision_species[0]) m_species1_index = i; - if (species_names[i] == collision_species[1]) m_species2_index = i; + for (int i=0; i<(int) beam_species_names.size(); i++) { + if (beam_species_names[i] == collision_species[0]) m_nbeams += 1; + if (beam_species_names[i] == collision_species[1]) m_nbeams += 1; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + m_nbeams <= 1, + ".species must contain at maximum one beam species name!"); + } + if (m_nbeams == 1) { + for (int i=0; i<(int) beam_species_names.size(); i++) { + if (beam_species_names[i] == collision_species[0] || + beam_species_names[i] == collision_species[1]) m_species1_index = i; + } + for (int i=0; i<(int) plasma_species_names.size(); i++) { + if (plasma_species_names[i] == collision_species[0] || + plasma_species_names[i] == collision_species[1]) m_species2_index = i; + } + } else { + for (int i=0; i<(int) plasma_species_names.size(); i++) { + if (plasma_species_names[i] == collision_species[0]) m_species1_index = i; + if (plasma_species_names[i] == collision_species[1]) m_species2_index = i; + } } + m_isSameSpecies = collision_species[0] == collision_species[1] ? true : false; AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_species1_index >= 0 && m_species2_index >= 0, - ".species must contain exactly the name of 2 species in plasma.names" + ".species must contain exactly the name of 2 species in plasmas.names " + "or one in plasmas.names and one in beams.names" ); } void -CoulombCollision::doCoulombCollision ( +CoulombCollision::doPlasmaPlasmaCoulombCollision ( int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1, PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog, amrex::Real background_density_SI) @@ -114,7 +134,7 @@ CoulombCollision::doCoulombCollision ( ux1, uy1, psi1, ux1, uy1, psi1, w1, w1, ion_lev1, ion_lev1, q1, q1, m1, m1, -1.0_rt, -1.0_rt, can_ionize1, can_ionize1, dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI, - normalized_units, background_density_SI, is_same_species, engine ); + normalized_units, background_density_SI, is_same_species, false, engine ); } ); count++; @@ -205,7 +225,7 @@ CoulombCollision::doCoulombCollision ( ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev1, ion_lev2, q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2, dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI, - normalized_units, background_density_SI, is_same_species, engine ); + normalized_units, background_density_SI, is_same_species, false, engine ); } ); count++; @@ -213,3 +233,115 @@ CoulombCollision::doCoulombCollision ( AMREX_ALWAYS_ASSERT(count == 1); } } + +void +CoulombCollision::doBeamPlasmaCoulombCollision ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, int islice_local, + BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, + amrex::Real background_density_SI) +{ + HIPACE_PROFILE("CoulombCollision::doBeamPlasmaCoulombCollision()"); + AMREX_ALWAYS_ASSERT(lev == 0); + + if (species1.TotalNumberOfParticles() == 0 || species2.TotalNumberOfParticles() == 0) return; + + using namespace amrex::literals; + const PhysConst cst = get_phys_const(); + bool normalized_units = Hipace::m_normalized_units; + + const amrex::Real clight = cst.c; + const amrex::Real inv_c = 1.0_rt / cst.c; + const amrex::Real inv_c2 = 1.0_rt / ( cst.c * cst.c ); + constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c; + constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c ); + + // Logically particles per-cell, and return indices of particles in each cell + BeamBins bins1 = findBeamParticlesInEachTile(bx, 1, islice_local, species1, geom); + PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom); + + int const n_cells = bins2.numBins(); + + // Counter to check there is only 1 box + int count = 0; + for (PlasmaParticleIterator pti(species2); pti.isValid(); ++pti) { + + // // Get particles SoA data for species 1 + const int box_offset = species1.m_box_sorter.boxOffsetsPtr()[species1.m_ibox]; + auto& soa1 = species1.GetStructOfArrays(); + amrex::Real* const ux1 = soa1.GetRealData(BeamIdx::ux).data() + box_offset; + amrex::Real* const uy1 = soa1.GetRealData(BeamIdx::uy).data() + box_offset; + amrex::Real* const psi1 = soa1.GetRealData(BeamIdx::uz).data() + box_offset; + const amrex::Real* const w1 = soa1.GetRealData(BeamIdx::w).data() + box_offset; + BeamBins::index_type * const indices1 = bins1.permutationPtr(); + BeamBins::index_type const * const offsets1 = bins1.offsetsPtr(); + amrex::Real q1 = species1.GetCharge(); + amrex::Real m1 = species1.GetMass(); + constexpr bool can_ionize1 = false; + + // Get particles SoA data for species 2 + //auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); + auto& soa2 = pti.GetStructOfArrays(); + amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux_half_step).data(); + amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy_half_step).data(); + amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi_half_step).data(); + const amrex::Real* const w2 = soa2.GetRealData(PlasmaIdx::w).data(); + const int* const ion_lev2 = soa2.GetIntData(PlasmaIdx::ion_lev).data(); + PlasmaBins::index_type * const indices2 = bins2.permutationPtr(); + PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr(); + amrex::Real q2 = species2.GetCharge(); + amrex::Real m2 = species2.GetMass(); + const bool can_ionize2 = species2.m_can_ionize; + + // volume is used to calculate density, but weights already represent density in normalized units + const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2); + // static_cast to avoid precision problems in FP32 + const amrex::Real wp = std::sqrt(static_cast(background_density_SI) * + PhysConstSI::q_e*PhysConstSI::q_e / + (PhysConstSI::ep0*PhysConstSI::m_e)); + const amrex::Real dt = normalized_units ? Hipace::GetInstance().m_dt/wp + : Hipace::GetInstance().m_dt; + + // Extract particles in the tile that `mfi` points to + // ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi); + // ParticleTileType& ptile_2 = species_2->ParticlesAt(lev, mfi); + // Loop over cells, and collide the particles in each cell + + // Loop over cells + amrex::ParallelForRNG( + n_cells, + [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept + { + // The particles from species1 that are in the cell `i_cell` are + // given by the `indices_1[cell_start_1:cell_stop_1]` + BeamBins::index_type const cell_start1 = offsets1[i_cell]; + BeamBins::index_type const cell_stop1 = offsets1[i_cell+1]; + // Same for species 2 + PlasmaBins::index_type const cell_start2 = offsets2[i_cell]; + PlasmaBins::index_type const cell_stop2 = offsets2[i_cell+1]; + + // ux from species1 can be accessed like this: + // ux_1[ indices_1[i] ], where i is between + // cell_start_1 (inclusive) and cell_start_2 (exclusive) + + // Do not collide if one species is missing in the cell + if ( cell_stop1 - cell_start1 < 1 || + cell_stop2 - cell_start2 < 1 ) return; + // shuffle + ShuffleFisherYates(indices1, cell_start1, cell_stop1, engine); + ShuffleFisherYates(indices2, cell_start2, cell_stop2, engine); + + // TODO: FIX DT. + // Call the function in order to perform collisions + 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 + q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2, + dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI, + normalized_units, background_density_SI, false, true, engine ); + } + ); + count++; + } + AMREX_ALWAYS_ASSERT(count == 1); +} diff --git a/src/particles/collisions/ElasticCollisionPerez.H b/src/particles/collisions/ElasticCollisionPerez.H index 0dfa68c581..cb33f0df54 100644 --- a/src/particles/collisions/ElasticCollisionPerez.H +++ b/src/particles/collisions/ElasticCollisionPerez.H @@ -50,6 +50,7 @@ * @param[in] normalized_units whether normalized units are used * @param[in] background_density_SI background plasma density (only needed for normalized units) * @param[in] is_same_species whether the collisions happen within the same species + * @param[in] is_beam_coll whether species1 is a beam * @param[in] engine AMReX engine for the random number generator. */ @@ -69,7 +70,7 @@ void ElasticCollisionPerez ( const bool can_ionize1, const bool can_ionize2, T_R const dt, T_R const L, T_R const inv_dV, T_R const clight, T_R const inv_c, T_R const inv_c_SI, T_R const inv_c2, T_R const inv_c2_SI, const bool normalized_units, - const amrex::Real background_density_SI, const bool is_same_species, + const amrex::Real background_density_SI, const bool is_same_species, bool is_beam_coll, amrex::RandomEngine const& engine) { using namespace amrex::literals; @@ -81,12 +82,13 @@ void ElasticCollisionPerez ( T_R T1t; T_R T2t; if ( T1 <= T_R(0.0) && L <= T_R(0.0) ) { - T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,psi1,m1,clight,inv_c2); + T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,psi1,m1,clight,inv_c2,is_beam_coll); } else { T1t = T1; } if ( T2 <= T_R(0.0) && L <= T_R(0.0) ) { - T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,psi2,m2,clight,inv_c2); + // second species is always plasma thus is_beam_coll is always false + T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,psi2,m2,clight,inv_c2,false); } else { T2t = T2; } @@ -151,24 +153,22 @@ void ElasticCollisionPerez ( if (can_ionize1) q1 *= ion_lev1[I1[i1]]; if (can_ionize2) q2 *= ion_lev2[I2[i2]]; // particle's Lorentz factor - amrex::Real g1 = ( - 1.0_rt + - u1x[I1[i1]]*u1x[I1[i1]]*inv_c2 + u1y[I1[i1]]*u1y[I1[i1]]*inv_c2 + - psi1[I1[i1]]*psi1[I1[i1]]) / (2.0_rt * psi1[I1[i1]] ); + amrex::Real g1 = is_beam_coll ? std::sqrt( 1._rt + + (u1x[I1[i1]]*u1x[I1[i1]] + u1y[I1[i1]]*u1y[I1[i1]] + psi1[I1[i1]]*psi1[I1[i1]])*inv_c2 ) + : (1.0_rt + u1x[I1[i1]]*u1x[I1[i1]]*inv_c2 + u1y[I1[i1]]*u1y[I1[i1]]*inv_c2 + + psi1[I1[i1]]*psi1[I1[i1]]) / (2.0_rt * psi1[I1[i1]] ); // particle's Lorentz factor - amrex::Real g2 = ( - 1.0_rt + - u2x[I2[i2]]*u2x[I2[i2]]*inv_c2 + u2y[I2[i2]]*u2y[I2[i2]]*inv_c2 + - psi2[I2[i2]]*psi2[I2[i2]]) / (2.0_rt * psi2[I2[i2]] ); + amrex::Real g2 = (1.0_rt + u2x[I2[i2]]*u2x[I2[i2]]*inv_c2 + u2y[I2[i2]]*u2y[I2[i2]]*inv_c2 + + psi2[I2[i2]]*psi2[I2[i2]]) / (2.0_rt * psi2[I2[i2]] ); // Convert from pseudo-potential to momentum - amrex::Real u1z = clight * (g1 - psi1[I1[i1]]); + amrex::Real u1z = is_beam_coll ? psi1[I1[i1]] : clight * (g1 - psi1[I1[i1]]); amrex::Real u2z = clight * (g2 - psi2[I2[i2]]); // In the longitudinal push of plasma particles, the dt is different for each particle. // The dt applied for collision probability is the average (in the lab frame) of these // dts. This is NOT clean. TODO FIXME. - const amrex::Real dt_fac = 0.5_rt * (g1/psi1[I1[i1]] + g2/psi2[I2[i2]]); + const amrex::Real dt_fac = is_beam_coll ? 1.0_rt : 0.5_rt * (g1/psi1[I1[i1]] + g2/psi2[I2[i2]]); UpdateMomentumPerezElastic( u1x[ I1[i1] ], u1y[ I1[i1] ], u1z, g1, u2x[ I2[i2] ], u2y[ I2[i2] ], u2z, g2, @@ -177,7 +177,7 @@ void ElasticCollisionPerez ( g1 = std::sqrt( T_R(1.0) + (u1x[I1[i1]]*u1x[I1[i1]]+u1y[I1[i1]]*u1y[I1[i1]]+u1z*u1z)*inv_c2 ); - psi1[I1[i1]] = g1 - u1z*inv_c; + psi1[I1[i1]] = is_beam_coll ? u1z : g1 - u1z*inv_c; g2 = std::sqrt( T_R(1.0) + (u2x[I2[i2]]*u2x[I2[i2]]+u2y[I2[i2]]*u2y[I2[i2]]+u2z*u2z)*inv_c2 ); psi2[I2[i2]] = g2 - u2z*inv_c; diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index 2ac504ea80..1616f8c2f9 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -122,16 +122,6 @@ public: /** returns number of plasma species */ int GetNPlasmas() const {return m_nplasmas;} - /** Perform binary elastic Coulomb collision, inter- and/or intra-species. - * - * The algorithm implemented is that of [Perez et al., Phys. Plasmas 19, 083104 (2012)] - * - * \param[in] lev MR level - * \param[in] bx box on which plasma particles are sorted per-cell and collide together - * \param[in] geom Corresponding gemetry - */ - void doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom); - /** Reorder particles to speed-up current deposition * \param[in] islice zeta slice index */ @@ -166,22 +156,15 @@ public: int m_sort_bin_size {32}; /**< Tile size to sort plasma particles */ - /** Number of binary collisions */ - int m_ncollisions = 0; - amrex::Vector m_all_plasmas; /**< contains all plasma containers */ int m_nplasmas = 0; /**< number of plasma containers */ amrex::Vector m_all_bins; /**< Logical tile bins for all plasma containers */ -private: amrex::Vector m_names {"no_plasma"}; /**< names of all plasma containers */ +private: /** Background (hypothetical) density, used to compute the adaptive time step */ amrex::Real m_adaptive_density = 0.; - /** User-input names of the binary collisions to be used */ - std::vector m_collision_names; - /** Vector of binary collisions */ - amrex::Vector< CoulombCollision > m_all_collisions; }; #endif // MULTIPLASMA_H_ diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index 46233015a2..a6954be746 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -21,8 +21,9 @@ MultiPlasma::MultiPlasma () queryWithParser(pp, "names", m_names); queryWithParser(pp, "adaptive_density", m_adaptive_density); queryWithParser(pp, "sort_bin_size", m_sort_bin_size); - queryWithParser(pp, "collisions", m_collision_names); + DeprecatedInput ("plasmas", "collisions", + "hipace.collisions", "", true); DeprecatedInput ("plasmas", "background_density_SI", "hipace.background_density_SI", "", true); @@ -33,16 +34,6 @@ MultiPlasma::MultiPlasma () m_all_plasmas.emplace_back(PlasmaParticleContainer(m_names[i])); } - /** Initialize the collision objects */ - m_ncollisions = m_collision_names.size(); - for (int i = 0; i < m_ncollisions; ++i) { - m_all_collisions.emplace_back(CoulombCollision(m_names, m_collision_names[i])); - } - if (Hipace::m_normalized_units && m_ncollisions > 0) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Hipace::m_background_density_SI!=0, - "For collisions with normalized units, a background plasma density must " - "be specified via 'hipace.background_density_SI'"); - } } void @@ -168,24 +159,6 @@ MultiPlasma::TileSort (amrex::Box bx, amrex::Geometry geom) } } -void -MultiPlasma::doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom) -{ - HIPACE_PROFILE("MultiPlasma::doCoulombCollision()"); - for (int i = 0; i < m_ncollisions; ++i) - { - AMREX_ALWAYS_ASSERT(lev == 0); - auto& species1 = m_all_plasmas[ m_all_collisions[i].m_species1_index ]; - auto& species2 = m_all_plasmas[ m_all_collisions[i].m_species2_index ]; - - // TODO: enable tiling - - CoulombCollision::doCoulombCollision( - lev, bx, geom, species1, species2, m_all_collisions[i].m_isSameSpecies, - m_all_collisions[i].m_CoulombLog, Hipace::m_background_density_SI); - } -} - void MultiPlasma::ReorderParticles (const int islice) { diff --git a/src/particles/sorting/TileSort.H b/src/particles/sorting/TileSort.H index 9137806881..e7254d2467 100644 --- a/src/particles/sorting/TileSort.H +++ b/src/particles/sorting/TileSort.H @@ -9,6 +9,7 @@ #define HIPACE_TILESORT_H_ #include "particles/plasma/PlasmaParticleContainer.H" +#include "particles/beam/BeamParticleContainer.H" #include @@ -28,4 +29,19 @@ findParticlesInEachTile ( amrex::Box bx, int bin_size, PlasmaParticleContainer& plasma, const amrex::Geometry& geom); +/** \brief Find beam particles in each bin, and return collections of indices per bin (tile). + * + * Note that this does *not* rearrange particle arrays + * + * \param[in] bx 3d box in which particles are sorted per slice + * \param[in] bin_size number of cells per tile (square) + * \param[in] islice_local slice of which the beam particles should be sorted + * \param[in] beam beam particle container + * \param[in] geom Geometry + */ +BeamBins +findBeamParticlesInEachTile ( + amrex::Box bx, int bin_size, int islice_local, + BeamParticleContainer& beam, const amrex::Geometry& geom); + #endif // HIPACE_TILESORT_H_ diff --git a/src/particles/sorting/TileSort.cpp b/src/particles/sorting/TileSort.cpp index f1abb2bf29..85d16343f4 100644 --- a/src/particles/sorting/TileSort.cpp +++ b/src/particles/sorting/TileSort.cpp @@ -53,3 +53,55 @@ findParticlesInEachTile ( AMREX_ALWAYS_ASSERT(count <= 1); return bins; } + +BeamBins +findBeamParticlesInEachTile ( + amrex::Box bx, int bin_size, int islice_local, + BeamParticleContainer& beam, const amrex::Geometry& geom) +{ + HIPACE_PROFILE("findBeamParticlesInEachTile()"); + + // Tile box: only 1 cell longitudinally, same as bx transversally. + const amrex::Box tcbx = bx.coarsen(bin_size); + const amrex::Box cbx = {{tcbx.smallEnd(0),tcbx.smallEnd(1),0}, {tcbx.bigEnd(0),tcbx.bigEnd(1),0}}; + + BeamBins bins; + + const int box_offset = beam.m_box_sorter.boxOffsetsPtr()[beam.m_ibox]; + + BeamBins::index_type const * const indices = beam.m_slice_bins.permutationPtr(); + BeamBins::index_type const * const slice_offsets = beam.m_slice_bins.offsetsPtrCpu(); + BeamBins::index_type const + cell_start = slice_offsets[islice_local], cell_stop = slice_offsets[islice_local+1]; + // The particles that are in slice islice_local are + // given by the indices[cell_start:cell_stop] + + int const num_particles = cell_stop-cell_start; + + // Extract box properties + const auto lo = lbound(cbx); + const auto dxi = amrex::GpuArray({ + geom.InvCellSizeArray()[0]/bin_size, + geom.InvCellSizeArray()[1]/bin_size, + 1.}); + const auto plo = geom.ProbLoArray(); + + // Find the particles that are in each slice and return collections of indices per slice. + bins.build( + num_particles, beam.getParticleTileData(), cbx, + // Pass lambda function that returns the slice index + [=] AMREX_GPU_DEVICE (const BeamParticleContainer::ParticleTileDataType& pdt, + unsigned int idx) + noexcept -> amrex::IntVect + { + auto p = pdt[indices[cell_start+idx] + box_offset]; + + return amrex::IntVect( + AMREX_D_DECL( + static_cast((p.pos(0)-plo[0])*dxi[0]-lo.x), + static_cast((p.pos(1)-plo[1])*dxi[1]-lo.y), + 0)); + }); + + return bins; +} diff --git a/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json b/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json new file mode 100644 index 0000000000..95d033bfb2 --- /dev/null +++ b/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json @@ -0,0 +1,32 @@ +{ + "lev=0": { + "Bx": 185897.4521776, + "By": 185897.45217553, + "Bz": 5395.5503546258, + "ExmBy": 84161261851900.0, + "EypBx": 84161261852289.0, + "Ez": 71705050369317.0, + "Psi": 1089585554.3813, + "Sx": 5219225582596500.0, + "Sy": 5219225582743600.0, + "chi": 4056826388137200.0, + "jx": 1.4609497758209e+16, + "jx_beam": 2.1124728260007e-06, + "jy": 1.4609497758179e+16, + "jy_beam": 4.2513859542859e-06, + "jz_beam": 1.0840992316526e+16, + "rhomjz": 178176970.1344 + }, + "beam": { + "charge": 1.1933011570032e-15, + "id": 27740076, + "mass": 6.7846689808772e-27, + "x": 0.03871, + "y": 0.03871, + "z": 0.2189712, + "ux": 0.0, + "uy": 0.0, + "uz": 14896000.0, + "w": 1692775082.5317 + } +} diff --git a/tests/checksum/reset_all_benchmarks.sh b/tests/checksum/reset_all_benchmarks.sh index ccf8647931..9a4ca5812e 100755 --- a/tests/checksum/reset_all_benchmarks.sh +++ b/tests/checksum/reset_all_benchmarks.sh @@ -367,6 +367,18 @@ then --test-name collisions.SI.1Rank fi +# collisions_beam.SI.1Rank +if [[ $all_tests = true ]] || [[ $one_test_name = "collisions_beam.SI.1Rank" ]] +then + cd $build_dir + ctest --output-on-failure -R collisions_beam.SI.1Rank \ + || echo "ctest command failed, maybe just because checksums are different. Keep going" + cd $checksum_dir + ./checksumAPI.py --reset-benchmark \ + --file_name ${build_dir}/bin/collisions_beam.SI.1Rank \ + --test-name collisions_beam.SI.1Rank +fi + # beam_in_vacuum_open_boundary.normalized.1Rank if [[ $all_tests = true ]] || [[ $one_test_name = "beam_in_vacuum_open_boundary.normalized.1Rank" ]] then diff --git a/tests/collisions.SI.1Rank.sh b/tests/collisions.SI.1Rank.sh index 6aaa80dedf..6ed72e0111 100755 --- a/tests/collisions.SI.1Rank.sh +++ b/tests/collisions.SI.1Rank.sh @@ -29,7 +29,7 @@ TEST_NAME="${FILE_NAME%.*}" OMP_NUM_THREADS=1 mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \ hipace.do_tiling = 0 \ hipace.file_prefix=$TEST_NAME \ - plasmas.collisions = collision1 \ + hipace.collisions = collision1 \ collision1.species = plasma plasma # Compare the results with checksum benchmark diff --git a/tests/collisions_beam.SI.1Rank.sh b/tests/collisions_beam.SI.1Rank.sh new file mode 100755 index 0000000000..39c43d9cf6 --- /dev/null +++ b/tests/collisions_beam.SI.1Rank.sh @@ -0,0 +1,40 @@ +#! /usr/bin/env bash + +# Copyright 2022 +# +# This file is part of HiPACE++. +# +# Authors: MaxThevenet +# License: BSD-3-Clause-LBNL + + +# This file is part of the HiPACE++ test suite. +# It runs a Hipace simulation with plasma collisions and compare with benchmark +# in normalized units + +# abort on first encounted error +set -eu -o pipefail + +# Read input parameters +HIPACE_EXECUTABLE=$1 +HIPACE_SOURCE_DIR=$2 + +HIPACE_EXAMPLE_DIR=${HIPACE_SOURCE_DIR}/examples/blowout_wake +HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests + +FILE_NAME=`basename "$0"` +TEST_NAME="${FILE_NAME%.*}" + +# Run the simulation +OMP_NUM_THREADS=1 mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \ + hipace.do_tiling = 0 \ + hipace.file_prefix=$TEST_NAME \ + hipace.collisions = collision1 \ + collision1.species = beam plasma + +# Compare the results with checksum benchmark +$HIPACE_TEST_DIR/checksum/checksumAPI.py \ + --evaluate \ + --file_name $TEST_NAME \ + --test-name $TEST_NAME \ + --skip "{'beam': 'id'}"