From 167f7e39aec230bf83bb38596087000d7622b8bd Mon Sep 17 00:00:00 2001 From: Ibrahim Najmudin Date: Wed, 26 Aug 2026 13:43:14 +0200 Subject: [PATCH 1/3] Added option collide_every, as an interval between applying the collision operator --- src/Hipace.H | 3 +- src/Hipace.cpp | 8 +- src/particles/collisions/CoulombCollision.H | 8 +- src/particles/collisions/CoulombCollision.cpp | 342 +++++++++--------- 4 files changed, 186 insertions(+), 175 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index 4a38fb33ff..12f7728761 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -147,8 +147,9 @@ public: /** * \brief does Coulomb collisions between plasmas and beams + * \param[in] islice slice number */ - void doCoulombCollision (); + void doCoulombCollision (int islice); /** * \brief add external fields to the field grid diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 4c81cdf2fc..6fe6eba8a5 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -858,7 +858,7 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st m_multi_beam.shiftSlippedParticles(islice, m_3D_geom[0]); // collisions for plasmas and beams - doCoulombCollision(); + doCoulombCollision(islice); // get minimum beam uz after push m_adaptive_time_step.GatherMinUzSlice(m_multi_beam, false); @@ -1286,7 +1286,7 @@ Hipace::AddGridExternalFields (const int lev, const int islice) } void -Hipace::doCoulombCollision () +Hipace::doCoulombCollision (int islice) { // collisions for all particles calculated on level 0 @@ -1301,7 +1301,7 @@ Hipace::doCoulombCollision () // TODO: enable tiling - CoulombCollision::doBeamPlasmaCoulombCollision( + CoulombCollision::doBeamPlasmaCoulombCollision( islice, lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2, m_all_collisions[i].m_CoulombLog, m_background_density_SI); } else { @@ -1311,7 +1311,7 @@ Hipace::doCoulombCollision () // TODO: enable tiling - CoulombCollision::doPlasmaPlasmaCoulombCollision( + CoulombCollision::doPlasmaPlasmaCoulombCollision( islice, m_all_collisions[i].m_collide_every, 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); } diff --git a/src/particles/collisions/CoulombCollision.H b/src/particles/collisions/CoulombCollision.H index 4d201e3a88..049b1a3f28 100644 --- a/src/particles/collisions/CoulombCollision.H +++ b/src/particles/collisions/CoulombCollision.H @@ -20,6 +20,7 @@ public: int m_nbeams {0}; bool m_isSameSpecies {false}; amrex::Real m_CoulombLog {-1.}; + int m_collide_every {1}; /** Read parameters from the input file */ void ReadParameters ( @@ -31,6 +32,8 @@ public: * \brief Perform Coulomb collisions of plasma species over longitudinal push by 1 cell. * Particles of both species are sorted per cell, paired, and collided pairwise. * + * \param[in] islice current slice number + * \param[in] collide_every interval between collision calculations * \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 @@ -41,7 +44,7 @@ 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 doPlasmaPlasmaCoulombCollision ( + static void doPlasmaPlasmaCoulombCollision ( int islice, int collide_every, 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); @@ -50,6 +53,7 @@ public: * \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] islice current slice number * \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 @@ -59,7 +63,7 @@ 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 doBeamPlasmaCoulombCollision ( + static void doBeamPlasmaCoulombCollision ( int islice, int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, amrex::Real background_density_SI); diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index 5d1466f73a..e7db80fb1a 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -25,6 +25,8 @@ CoulombCollision::ReadParameters( // default Coulomb log is -1, if < 0 (e.g. not specified), will be computed automatically pp.query("CoulombLog", m_CoulombLog); + // how often the collision operator should be applied - every m_collide_every'th slice + pp.query("collide_every", m_collide_every); for (int i=0; i<(int) beam_species_names.size(); i++) { if (beam_species_names[i] == collision_species[0]) m_nbeams += 1; @@ -58,183 +60,187 @@ CoulombCollision::ReadParameters( } void -CoulombCollision::doPlasmaPlasmaCoulombCollision ( +CoulombCollision::doPlasmaPlasmaCoulombCollision ( int islice, int collide_every, 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) { - HIPACE_PROFILE("CoulombCollision::doCoulombCollision()"); - AMREX_ALWAYS_ASSERT(lev == 0); - - if (species1.TotalNumberOfParticles(false, true) == 0 || - species2.TotalNumberOfParticles(false, true) == 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; - constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c; - constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c ); - - if ( is_same_species ) // species_1 == species_2 - { - // Logically particles per-cell, and return indices of particles in each cell - PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); - int const n_cells = bins1.numBins(); - - // Counter to check there is only 1 box - int count = 0; - for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { - - // Get particles SoA data - auto& ptile1 = pti.GetParticleTile(); - amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); - 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(); - PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); - PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); - amrex::Real q1 = species1.GetCharge(); - amrex::Real m1 = species1.GetMass(); - const bool can_ionize1 = species1.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 ? geom.CellSize(2)/wp - : geom.CellSize(2)/PhysConstSI::c; - - 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]` - PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; - PlasmaBins::index_type const cell_stop1 = offsets1[i_cell+1]; - PlasmaBins::index_type const cell_half1 = (cell_start1+cell_stop1)/2; - - if ( cell_stop1 - cell_start1 <= 1 ) return; - // Do not collide if there is only one particle in the cell - // shuffle - ShuffleFisherYates( - indices1, cell_start1, cell_half1, engine ); - - // TODO: FIX DT - // Call the function in order to perform collisions - ElasticCollisionPerez( - cell_start1, cell_half1, - cell_half1, cell_stop1, - indices1, indices1, - 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_SI, inv_c2_SI, - normalized_units, background_density_SI, is_same_species, false, engine ); - } - ); - count++; - } - AMREX_ALWAYS_ASSERT(count == 1); - - } else { - - // Logically particles per-cell, and return indices of particles in each cell - PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); - PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom); - - int const n_cells = bins1.numBins(); - - // Counter to check there is only 1 box - int count = 0; - for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { - - // Get particles SoA data for species 1 - auto& ptile1 = pti.GetParticleTile(); - amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); - 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(); - PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); - PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); - amrex::Real q1 = species1.GetCharge(); - amrex::Real m1 = species1.GetMass(); - const bool can_ionize1 = species1.m_can_ionize; - - // Get particles SoA data for species 2 - auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); - amrex::Real* const ux2 = ptile2.GetRealData(PlasmaIdx::ux_half_step).data(); - 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(); - 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 ? geom.CellSize(2)/wp - : geom.CellSize(2)/PhysConstSI::c; - // 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]` - PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; - PlasmaBins::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_lev1, ion_lev2, - 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, is_same_species, false, engine ); - } - ); - count++; + if (islice%collide_every == 0){ + HIPACE_PROFILE("CoulombCollision::doCoulombCollision()"); + AMREX_ALWAYS_ASSERT(lev == 0); + + if (species1.TotalNumberOfParticles(false, true) == 0 || + species2.TotalNumberOfParticles(false, true) == 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; + constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c; + constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c ); + + if ( is_same_species ) // species_1 == species_2 + { + // Logically particles per-cell, and return indices of particles in each cell + PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); + int const n_cells = bins1.numBins(); + + // Counter to check there is only 1 box + int count = 0; + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + + // Get particles SoA data + auto& ptile1 = pti.GetParticleTile(); + amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); + 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(); + PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); + PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); + amrex::Real q1 = species1.GetCharge(); + amrex::Real m1 = species1.GetMass(); + const bool can_ionize1 = species1.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)); + amrex::Real dt = normalized_units ? geom.CellSize(2)/wp + : geom.CellSize(2)/PhysConstSI::c; + dt = dt * collide_every; + + 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]` + PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; + PlasmaBins::index_type const cell_stop1 = offsets1[i_cell+1]; + PlasmaBins::index_type const cell_half1 = (cell_start1+cell_stop1)/2; + + if ( cell_stop1 - cell_start1 <= 1 ) return; + // Do not collide if there is only one particle in the cell + // shuffle + ShuffleFisherYates( + indices1, cell_start1, cell_half1, engine ); + + // TODO: FIX DT + // Call the function in order to perform collisions + ElasticCollisionPerez( + cell_start1, cell_half1, + cell_half1, cell_stop1, + indices1, indices1, + 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_SI, inv_c2_SI, + normalized_units, background_density_SI, is_same_species, false, engine ); + } + ); + count++; + } + AMREX_ALWAYS_ASSERT(count == 1); + + } else { + + // Logically particles per-cell, and return indices of particles in each cell + PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); + PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom); + + int const n_cells = bins1.numBins(); + + // Counter to check there is only 1 box + int count = 0; + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + + // Get particles SoA data for species 1 + auto& ptile1 = pti.GetParticleTile(); + amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); + 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(); + PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); + PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); + amrex::Real q1 = species1.GetCharge(); + amrex::Real m1 = species1.GetMass(); + const bool can_ionize1 = species1.m_can_ionize; + + // Get particles SoA data for species 2 + auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); + amrex::Real* const ux2 = ptile2.GetRealData(PlasmaIdx::ux_half_step).data(); + 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(); + 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)); + amrex::Real dt = normalized_units ? geom.CellSize(2)/wp + : geom.CellSize(2)/PhysConstSI::c; + dt = dt * collide_every; + // 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]` + PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; + PlasmaBins::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_lev1, ion_lev2, + 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, is_same_species, false, engine ); + } + ); + count++; + } + AMREX_ALWAYS_ASSERT(count == 1); } - AMREX_ALWAYS_ASSERT(count == 1); } } void -CoulombCollision::doBeamPlasmaCoulombCollision ( +CoulombCollision::doBeamPlasmaCoulombCollision ( int islice, int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, amrex::Real background_density_SI) From 2fba76701ed1b38711806e271c5a5452c75730f2 Mon Sep 17 00:00:00 2001 From: Ibrahim Najmudin Date: Thu, 10 Sep 2026 13:41:48 +0200 Subject: [PATCH 2/3] Finished addition of a collide_every option for coarsening collisions, to speed up simulation. --- src/Hipace.cpp | 2 +- src/particles/collisions/CoulombCollision.H | 3 +- src/particles/collisions/CoulombCollision.cpp | 344 +++++++++--------- 3 files changed, 176 insertions(+), 173 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 6fe6eba8a5..4608d9c24f 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -1301,7 +1301,7 @@ Hipace::doCoulombCollision (int islice) // TODO: enable tiling - CoulombCollision::doBeamPlasmaCoulombCollision( islice, + CoulombCollision::doBeamPlasmaCoulombCollision( islice, m_all_collisions[i].m_collide_every, lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2, m_all_collisions[i].m_CoulombLog, m_background_density_SI); } else { diff --git a/src/particles/collisions/CoulombCollision.H b/src/particles/collisions/CoulombCollision.H index 049b1a3f28..c2597c8a07 100644 --- a/src/particles/collisions/CoulombCollision.H +++ b/src/particles/collisions/CoulombCollision.H @@ -54,6 +54,7 @@ public: * Particles of both species are sorted per cell, paired, and collided pairwise. * * \param[in] islice current slice number + * \param[in] collide_every interval between collision operations * \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 @@ -63,7 +64,7 @@ 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 doBeamPlasmaCoulombCollision ( int islice, + static void doBeamPlasmaCoulombCollision ( int islice, int collide_every, int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, amrex::Real background_density_SI); diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index e7db80fb1a..be244be22c 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -65,186 +65,188 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision ( int islice, int collide_every PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog, amrex::Real background_density_SI) { - if (islice%collide_every == 0){ - HIPACE_PROFILE("CoulombCollision::doCoulombCollision()"); - AMREX_ALWAYS_ASSERT(lev == 0); - - if (species1.TotalNumberOfParticles(false, true) == 0 || - species2.TotalNumberOfParticles(false, true) == 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; - constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c; - constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c ); - - if ( is_same_species ) // species_1 == species_2 - { - // Logically particles per-cell, and return indices of particles in each cell - PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); - int const n_cells = bins1.numBins(); - - // Counter to check there is only 1 box - int count = 0; - for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { - - // Get particles SoA data - auto& ptile1 = pti.GetParticleTile(); - amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); - 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(); - PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); - PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); - amrex::Real q1 = species1.GetCharge(); - amrex::Real m1 = species1.GetMass(); - const bool can_ionize1 = species1.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)); - amrex::Real dt = normalized_units ? geom.CellSize(2)/wp - : geom.CellSize(2)/PhysConstSI::c; - dt = dt * collide_every; - - 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]` - PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; - PlasmaBins::index_type const cell_stop1 = offsets1[i_cell+1]; - PlasmaBins::index_type const cell_half1 = (cell_start1+cell_stop1)/2; - - if ( cell_stop1 - cell_start1 <= 1 ) return; - // Do not collide if there is only one particle in the cell - // shuffle - ShuffleFisherYates( - indices1, cell_start1, cell_half1, engine ); - - // TODO: FIX DT - // Call the function in order to perform collisions - ElasticCollisionPerez( - cell_start1, cell_half1, - cell_half1, cell_stop1, - indices1, indices1, - 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_SI, inv_c2_SI, - normalized_units, background_density_SI, is_same_species, false, engine ); - } - ); - count++; - } - AMREX_ALWAYS_ASSERT(count == 1); - - } else { - - // Logically particles per-cell, and return indices of particles in each cell - PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); - PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom); - - int const n_cells = bins1.numBins(); - - // Counter to check there is only 1 box - int count = 0; - for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { - - // Get particles SoA data for species 1 - auto& ptile1 = pti.GetParticleTile(); - amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); - 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(); - PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); - PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); - amrex::Real q1 = species1.GetCharge(); - amrex::Real m1 = species1.GetMass(); - const bool can_ionize1 = species1.m_can_ionize; - - // Get particles SoA data for species 2 - auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); - amrex::Real* const ux2 = ptile2.GetRealData(PlasmaIdx::ux_half_step).data(); - 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(); - 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)); - amrex::Real dt = normalized_units ? geom.CellSize(2)/wp - : geom.CellSize(2)/PhysConstSI::c; - dt = dt * collide_every; - // 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]` - PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; - PlasmaBins::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_lev1, ion_lev2, - 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, is_same_species, false, engine ); - } - ); - count++; - } - AMREX_ALWAYS_ASSERT(count == 1); + if (islice%collide_every != 0) {return;} + + HIPACE_PROFILE("CoulombCollision::doCoulombCollision()"); + AMREX_ALWAYS_ASSERT(lev == 0); + + if (species1.TotalNumberOfParticles(false, true) == 0 || + species2.TotalNumberOfParticles(false, true) == 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; + constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c; + constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c ); + + if ( is_same_species ) // species_1 == species_2 + { + // Logically particles per-cell, and return indices of particles in each cell + PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); + int const n_cells = bins1.numBins(); + + // Counter to check there is only 1 box + int count = 0; + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + + // Get particles SoA data + auto& ptile1 = pti.GetParticleTile(); + amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); + 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(); + PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); + PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); + amrex::Real q1 = species1.GetCharge(); + amrex::Real m1 = species1.GetMass(); + const bool can_ionize1 = species1.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)); + amrex::Real dt = normalized_units ? geom.CellSize(2)/wp + : geom.CellSize(2)/PhysConstSI::c; + dt = dt * collide_every; + + 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]` + PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; + PlasmaBins::index_type const cell_stop1 = offsets1[i_cell+1]; + PlasmaBins::index_type const cell_half1 = (cell_start1+cell_stop1)/2; + + if ( cell_stop1 - cell_start1 <= 1 ) return; + // Do not collide if there is only one particle in the cell + // shuffle + ShuffleFisherYates( + indices1, cell_start1, cell_half1, engine ); + + // TODO: FIX DT + // Call the function in order to perform collisions + ElasticCollisionPerez( + cell_start1, cell_half1, + cell_half1, cell_stop1, + indices1, indices1, + 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_SI, inv_c2_SI, + normalized_units, background_density_SI, is_same_species, false, engine ); + } + ); + count++; + } + AMREX_ALWAYS_ASSERT(count == 1); + + } else { + + // Logically particles per-cell, and return indices of particles in each cell + PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); + PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom); + + int const n_cells = bins1.numBins(); + + // Counter to check there is only 1 box + int count = 0; + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + + // Get particles SoA data for species 1 + auto& ptile1 = pti.GetParticleTile(); + amrex::Real* const ux1 = ptile1.GetRealData(PlasmaIdx::ux_half_step).data(); + 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(); + PlasmaBins::index_type * const indices1 = bins1.permutationPtr(); + PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr(); + amrex::Real q1 = species1.GetCharge(); + amrex::Real m1 = species1.GetMass(); + const bool can_ionize1 = species1.m_can_ionize; + + // Get particles SoA data for species 2 + auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex()); + amrex::Real* const ux2 = ptile2.GetRealData(PlasmaIdx::ux_half_step).data(); + 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(); + 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)); + amrex::Real dt = normalized_units ? geom.CellSize(2)/wp + : geom.CellSize(2)/PhysConstSI::c; + dt = dt * collide_every; + // 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]` + PlasmaBins::index_type const cell_start1 = offsets1[i_cell]; + PlasmaBins::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_lev1, ion_lev2, + 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, is_same_species, false, engine ); + } + ); + count++; } + AMREX_ALWAYS_ASSERT(count == 1); } } void -CoulombCollision::doBeamPlasmaCoulombCollision ( int islice, +CoulombCollision::doBeamPlasmaCoulombCollision ( int islice, int collide_every, int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog, amrex::Real background_density_SI) { + if (islice%collide_every != 0) {return;} + HIPACE_PROFILE("CoulombCollision::doBeamPlasmaCoulombCollision()"); AMREX_ALWAYS_ASSERT(lev == 0); From 0d8facbdee89d15e829739b59263d1f26aa3fcca Mon Sep 17 00:00:00 2001 From: Ibrahim Najmudin Date: Thu, 17 Sep 2026 14:43:52 +0200 Subject: [PATCH 3/3] Removed trailing whitespaces --- src/particles/collisions/CoulombCollision.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp index be244be22c..e29c7dd3b4 100644 --- a/src/particles/collisions/CoulombCollision.cpp +++ b/src/particles/collisions/CoulombCollision.cpp @@ -246,7 +246,7 @@ CoulombCollision::doBeamPlasmaCoulombCollision ( int islice, int collide_every, amrex::Real background_density_SI) { if (islice%collide_every != 0) {return;} - + HIPACE_PROFILE("CoulombCollision::doBeamPlasmaCoulombCollision()"); AMREX_ALWAYS_ASSERT(lev == 0);