From 1c5f3116058c4ee84287f77342c3f3d27c79ca09 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Tue, 17 Feb 2026 16:34:17 +0100 Subject: [PATCH 1/6] Collision Skeleton --- src/particles/collisions/CMakeLists.txt | 1 + src/particles/collisions/Collision.H | 36 ++++++++++++++++ src/particles/collisions/Collision.cpp | 55 +++++++++++++++++++++++++ 3 files changed, 92 insertions(+) create mode 100644 src/particles/collisions/Collision.H create mode 100644 src/particles/collisions/Collision.cpp diff --git a/src/particles/collisions/CMakeLists.txt b/src/particles/collisions/CMakeLists.txt index 6a5cd3cf3c..5f37663fb8 100644 --- a/src/particles/collisions/CMakeLists.txt +++ b/src/particles/collisions/CMakeLists.txt @@ -8,4 +8,5 @@ target_sources(HiPACE PRIVATE CoulombCollision.cpp + Collision.cpp ) diff --git a/src/particles/collisions/Collision.H b/src/particles/collisions/Collision.H new file mode 100644 index 0000000000..3925439924 --- /dev/null +++ b/src/particles/collisions/Collision.H @@ -0,0 +1,36 @@ +#ifndef HIPACE_COLLISION_H_ +#define HIPACE_COLLISION_H_ + +#include "particles/plasma/MultiPlasma.H" + +#include + +/** + * \brief This class handles Coulomb collisions between 2 particle species + * (can be plasma-plasma or beam-plasma, the species can be the same) + */ +class Collision +{ +public: + int m_inout_species1_index = -1; + int m_inout_species2_index = -1; + int m_out_species3_index = -1; + std::string m_collision_type = ""; + + /** Read parameters from the input file */ + void ReadParameters ( + const std::vector& plasma_species_names, + std::string const collision_name); + + /** + * \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,out] multi_plasma all plasmas + **/ + static void doCollision ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma); +}; + +#endif // HIPACE_COLLISION_H_ diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp new file mode 100644 index 0000000000..bae1e2724d --- /dev/null +++ b/src/particles/collisions/Collision.cpp @@ -0,0 +1,55 @@ +#include "Collision.H" +#include "Hipace.H" + +void +Collision::ReadParameters( + const std::vector& plasma_species_names, + std::string const collision_name) +{ + amrex::ParmParse pp(collision_name); + + getWithParser(pp, "type", m_collision_type); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + m_collision_type == "a" || + m_collision_type == "b", + "Unknown collision type" + ); + + std::string name1 = ""; + getWithParser(pp, "species1", name1); + std::string name2 = ""; + getWithParser(pp, "species2", name1); + + std::string name3 = ""; + if (m_collision_type == "a") { + getWithParser(pp, "species3", name1); + } + + for (int i = 0; i < plasma_species_names.size(); ++i) { + if (plasma_species_names[i] == name1) { + m_inout_species1_index = i; + } + if (plasma_species_names[i] == name2) { + m_inout_species2_index = i; + } + if (plasma_species_names[i] == name3) { + m_out_species3_index = i; + } + } + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + m_inout_species1_index != -1, + "Collision: Unknown species1" + ); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + m_inout_species2_index != -1, + "Collision: Unknown species2" + ); + if (m_collision_type == "a") { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + m_out_species3_index != -1, + "Collision: Unknown species3" + ); + } +} From 44ba1b99bfe85d53e4de4bdef16b83560782826e Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Fri, 20 Feb 2026 18:44:06 +0100 Subject: [PATCH 2/6] Add doCollisionImp --- src/particles/collisions/Collision.H | 51 ++++- src/particles/collisions/Collision.cpp | 277 ++++++++++++++++++++++--- src/particles/plasma/MultiPlasma.H | 10 + 3 files changed, 307 insertions(+), 31 deletions(-) diff --git a/src/particles/collisions/Collision.H b/src/particles/collisions/Collision.H index 3925439924..a17b97f2e6 100644 --- a/src/particles/collisions/Collision.H +++ b/src/particles/collisions/Collision.H @@ -5,17 +5,22 @@ #include +#include + /** * \brief This class handles Coulomb collisions between 2 particle species * (can be plasma-plasma or beam-plasma, the species can be the same) */ class Collision { -public: - int m_inout_species1_index = -1; - int m_inout_species2_index = -1; - int m_out_species3_index = -1; + std::string m_inout_species1_name = ""; + std::string m_inout_species2_name = ""; + std::string m_out_species3_name = ""; std::string m_collision_type = ""; + bool m_has_collision_product = false; + amrex::Gpu::Buffer m_crosssection_data; + +public: /** Read parameters from the input file */ void ReadParameters ( @@ -28,7 +33,43 @@ public: * \param[in] geom corresponding geometry object * \param[in,out] multi_plasma all plasmas **/ - static void doCollision ( + void doCollision ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma); + + /** + * \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,out] multi_plasma all plasmas + * \param[in] collision_function return if a pair of particles collides + * \param[in] ionizaiton_function initialize product particle + **/ + template + void doCollisionImp ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionizaiton_function); + + + /** + * \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,out] multi_plasma all plasmas + **/ + void doCollisionA ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma); + + /** + * \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,out] multi_plasma all plasmas + **/ + void doCollisionB ( int lev, const amrex::Box& bx, const amrex::Geometry& geom, MultiPlasma& multi_plasma); }; diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index bae1e2724d..814e9f61ef 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -1,5 +1,7 @@ #include "Collision.H" #include "Hipace.H" +#include "particles/sorting/TileSort.H" +#include "ShuffleFisherYates.H" void Collision::ReadParameters( @@ -16,40 +18,263 @@ Collision::ReadParameters( "Unknown collision type" ); - std::string name1 = ""; - getWithParser(pp, "species1", name1); - std::string name2 = ""; - getWithParser(pp, "species2", name1); + m_has_collision_product = m_collision_type == "a"; - std::string name3 = ""; - if (m_collision_type == "a") { - getWithParser(pp, "species3", name1); + getWithParser(pp, "species1", m_inout_species1_name); + getWithParser(pp, "species2", m_inout_species2_name); + + if (m_has_collision_product) { + getWithParser(pp, "species3", m_out_species3_name); } - for (int i = 0; i < plasma_species_names.size(); ++i) { - if (plasma_species_names[i] == name1) { - m_inout_species1_index = i; - } - if (plasma_species_names[i] == name2) { - m_inout_species2_index = i; - } - if (plasma_species_names[i] == name3) { - m_out_species3_index = i; + + if (m_collision_type == "a") { + // Initialize crosssection data for collision calculation + m_crosssection_data.resize(10); + for (int i=0; i<10 ; ++i) { + m_crosssection_data[i] = i; } + m_crosssection_data.copyToDeviceAsync(); } +} - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_inout_species1_index != -1, - "Collision: Unknown species1" +void +Collision::doCollision ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma) +{ + if (m_collision_type == "a") { + doCollisionA(lev, bx, geom, multi_plasma); + } else { + doCollisionB(lev, bx, geom, multi_plasma); + } +} + +void +Collision::doCollisionA ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma) +{ + auto p_crosssection_data = m_crosssection_data.data(); + + doCollisionImp(lev, bx, geom, multi_plasma, + [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, + amrex::RandomEngine const& engine) + { + // do some calculation with the two particles to see if they should collide + // can modify the particles here + const amrex::Real ux1 = ptd1.rdata(PlasmaIdx::ux)[i1]; + const amrex::Real ux2 = ptd2.rdata(PlasmaIdx::ux)[i2]; + + const int ion_lev = ptd1.idata(PlasmaIdx::ion_lev)[i1]; + const amrex::Real crosssection = p_crosssection_data[ion_lev]; + + return (crosssection + ux1) < ux2; + }, + [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, auto ptd3, int i3, + amrex::RandomEngine const& engine) + { + // initialize the new particle + // can modify all three particles here + const amrex::Real ux1 = ptd1.rdata(PlasmaIdx::ux)[i1]; + const amrex::Real ux2 = ptd2.rdata(PlasmaIdx::ux)[i2]; + ptd3.rdata(PlasmaIdx::ux)[i3] = ux1 + ux2; + }); +} + +void +Collision::doCollisionB ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma) +{ + doCollisionImp(lev, bx, geom, multi_plasma, + [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, + amrex::RandomEngine const& engine) + { + return false; + }, + [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, auto ptd3, int i3, + amrex::RandomEngine const& engine) + { + + }); +} + +template +void +Collision::doCollisionImp ( + int lev, const amrex::Box& bx, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionizaiton_function) +{ + // assume the two species are different for now + auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name); + auto& species2 = multi_plasma.GetPlasma(m_inout_species2_name); + + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + amrex::removeInvalidParticles(species1.ParticlesAt(0, pti)); + amrex::removeInvalidParticles(species2.ParticlesAt(0, pti)); + } + + PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); + PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom); + + auto offset1 = bins1.offsetsPtr(); + auto offset2 = bins2.offsetsPtr(); + auto perm1 = bins1.permutationPtr(); + auto perm2 = bins2.permutationPtr(); + + const int num_cells = bins1.numBins(); + + amrex::Gpu::DeviceVector num_ind_pairs(num_cells, 0); + auto p_num_ind_pairs = num_ind_pairs.dataPtr(); + + const int total_ind_pairs = amrex::Scan::PrefixSum(num_cells, + [=] AMREX_GPU_DEVICE (int i) { + int n1 = offset1[i+1] - offset1[i]; + int n2 = offset2[i+1] - offset2[i]; + return std::min(n1, n2); + }, + [=] AMREX_GPU_DEVICE (int i, int s) { + p_num_ind_pairs[i] = s; + }, + amrex::Scan::Type::exclusive, amrex::Scan::retSum ); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_inout_species2_index != -1, - "Collision: Unknown species2" + + amrex::ParallelForRNG(num_cells, + [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) { + ShuffleFisherYates(perm1, offset1[i], offset1[i+1], engine); + ShuffleFisherYates(perm2, offset2[i], offset2[i+1], engine); + } ); - if (m_collision_type == "a") { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_out_species3_index != -1, - "Collision: Unknown species3" + + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + + auto& ptile1 = species1.ParticlesAt(0, pti); + auto& ptile2 = species2.ParticlesAt(0, pti); + + auto ptd1 = ptile1.getParticleTileData(); + auto ptd2 = ptile2.getParticleTileData(); + + const int np1 = ptile1.numParticles(); + const int np2 = ptile2.numParticles(); + + amrex::Gpu::DeviceVector flag1(np1, 0); + amrex::Gpu::DeviceVector flag2(np2, 0); + + auto p_flag1 = flag1.dataPtr(); + auto p_flag2 = flag2.dataPtr(); + + amrex::ParallelForRNG(total_ind_pairs, + [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + + const int offset1_start = offset1[icell]; + const int offset1_stop = offset1[icell+1]; + const int offset2_start = offset2[icell]; + const int offset2_stop = offset2[icell+1]; + + const int icoll = ipair - p_num_ind_pairs[icell]; + + const int n1 = offset1_stop - offset1_start; + const int n2 = offset2_stop - offset2_start; + + int idx1 = icoll; + int idx2 = icoll; + while (idx1 < n1 && idx2 < n2) { + const int j1 = perm1[offset1_start + idx1]; + const int j2 = perm2[offset2_start + idx2]; + + const bool make_new_particle = collision_function( + ptd1, j1, ptd2, j2, engine + ); + + if (n2 < n1) { + idx1 += n2; + if (make_new_particle) { + p_flag2[j2] = 1; + } + } else { + idx2 += n1; + if (make_new_particle) { + p_flag1[j1] = 1; + } + } + } + } + ); + + if (!m_has_collision_product) { + continue; + } + + const int num_new_particles = amrex::Scan::PrefixSum(np1+np2, + [=] AMREX_GPU_DEVICE (int i) { + if (i < np1) { + return p_flag1[i]; + } else { + i -= np1; + return p_flag2[i]; + } + }, + [=] AMREX_GPU_DEVICE (int i, int s) { + if (i < np1) { + p_flag1[i] = p_flag1[i] == 0 ? -1 : s; + } else { + i -= np1; + p_flag2[i] = p_flag2[i] == 0 ? -1 : s; + } + }, + amrex::Scan::Type::exclusive, amrex::Scan::retSum + ); + + auto& species3 = multi_plasma.GetPlasma(m_out_species3_name); + auto& ptile3 = species3.ParticlesAt(0, pti); + int old_size3 = ptile3.size(); + ptile3.resize(old_size3 + num_new_particles); + + auto ptd3 = ptile3.getParticleTileData(); + + amrex::ParallelForRNG(total_ind_pairs, + [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + + const int offset1_start = offset1[icell]; + const int offset1_stop = offset1[icell+1]; + const int offset2_start = offset2[icell]; + const int offset2_stop = offset2[icell+1]; + + const int icoll = ipair - p_num_ind_pairs[icell]; + + const int n1 = offset1_stop - offset1_start; + const int n2 = offset2_stop - offset2_start; + + int idx1 = icoll; + int idx2 = icoll; + while (idx1 < n1 && idx2 < n2) { + const int j1 = perm1[offset1_start + idx1]; + const int j2 = perm2[offset2_start + idx2]; + const int new_part_idx = n2 < n1 ? p_flag2[j2] : p_flag1[j1]; + + if (new_part_idx >= 0) { + ionizaiton_function( + ptd1, j1, + ptd2, j2, + ptd3, old_size3 + new_part_idx, + engine + ); + } + + if (n2 < n1) { + idx1 += n2; + } else { + idx2 += n1; + } + } + } ); + + amrex::Gpu::streamSynchronize(); } } diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index 12549d09d0..3212d35961 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -159,6 +159,16 @@ public: void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom, int max_step, amrex::Real max_time); + PlasmaParticleContainer& GetPlasma (const std::string& name) { + for (int i = 0; i < m_nplasmas; ++i) { + if (m_names[i] == name) { + return m_all_plasmas[i]; + } + } + amrex::Abort("Unknown plasma name " + name); + return m_all_plasmas[0]; + } + amrex::Vector m_all_plasmas; /**< contains all plasma containers */ int m_nplasmas = 0; /**< number of plasma containers */ amrex::Vector m_names {"no_plasma"}; /**< names of all plasma containers */ From 74605e9814b6e362cb0bfe7a35345545887ac388 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 23 Feb 2026 00:01:11 +0100 Subject: [PATCH 3/6] add cell weight --- src/Hipace.H | 4 ++ src/Hipace.cpp | 10 ++++ src/particles/collisions/Collision.H | 13 ++---- src/particles/collisions/Collision.cpp | 64 ++++++++++++++++++++------ 4 files changed, 69 insertions(+), 22 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index 240e357775..d94f696367 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -22,6 +22,7 @@ #include "utils/MultiBuffer.H" #include "diagnostics/Diagnostic.H" #include "diagnostics/OpenPMDWriter.H" +#include "particles/collisions/Collision.H" #include #ifdef AMREX_USE_LINEAR_SOLVERS @@ -351,6 +352,9 @@ private: /** Vector of binary collisions */ amrex::Vector< CoulombCollision > m_all_collisions; + amrex::Vector m_collision_names2; + amrex::Vector m_all_collisions2; + void InitDiagnostics (const int step); void FillFieldDiagnostics (const int current_N_level, int islice); void FillBeamDiagnostics (const int step); diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 5d551e7dab..dcf4b1c188 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -257,6 +257,12 @@ Hipace::ReadParameters () "be specified via 'hipace.background_density_SI'"); } + queryWithParser(pph, "collisions2", m_collision_names2); + for (int i = 0; i != m_collision_names2.size(); ++i) { + m_all_collisions2.emplace_back(Collision()); + m_all_collisions2.back().ReadParameters(m_multi_plasma.m_names, m_collision_names2[i]); + } + // external fields applied to the grid amrex::Array field_str = {"0", "0", "0", "0", "0"}; m_use_grid_external_fields = queryWithParser(pph, "grid_external_fields(x,y,z,t)", field_str); @@ -820,6 +826,10 @@ Hipace::SolveOneSlice (int islice, int step) // collisions for plasmas and beams doCoulombCollision(); + for (auto& collision : m_all_collisions2) { + collision.doCollision(0, m_slice_geom[0], m_multi_plasma); + } + // get minimum beam uz after push m_adaptive_time_step.GatherMinUzSlice(m_multi_beam, false); diff --git a/src/particles/collisions/Collision.H b/src/particles/collisions/Collision.H index a17b97f2e6..16cc4e249e 100644 --- a/src/particles/collisions/Collision.H +++ b/src/particles/collisions/Collision.H @@ -29,17 +29,15 @@ public: /** * \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,out] multi_plasma all plasmas **/ void doCollision ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma); /** * \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,out] multi_plasma all plasmas * \param[in] collision_function return if a pair of particles collides @@ -47,30 +45,27 @@ public: **/ template void doCollisionImp ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma, F const& collision_function, G const& ionizaiton_function); - /** * \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,out] multi_plasma all plasmas **/ void doCollisionA ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma); /** * \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,out] multi_plasma all plasmas **/ void doCollisionB ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma); }; diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index 814e9f61ef..99d42ca7ab 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -40,25 +40,26 @@ Collision::ReadParameters( void Collision::doCollision ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma) { if (m_collision_type == "a") { - doCollisionA(lev, bx, geom, multi_plasma); + doCollisionA(lev, geom, multi_plasma); } else { - doCollisionB(lev, bx, geom, multi_plasma); + doCollisionB(lev, geom, multi_plasma); } } void Collision::doCollisionA ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma) { auto p_crosssection_data = m_crosssection_data.data(); - doCollisionImp(lev, bx, geom, multi_plasma, + doCollisionImp(lev, geom, multi_plasma, [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, + amrex::Real cell_weight1, amrex::Real cell_weight2, amrex::RandomEngine const& engine) { // do some calculation with the two particles to see if they should collide @@ -84,11 +85,12 @@ Collision::doCollisionA ( void Collision::doCollisionB ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma) { - doCollisionImp(lev, bx, geom, multi_plasma, + doCollisionImp(lev, geom, multi_plasma, [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, + amrex::Real cell_weight1, amrex::Real cell_weight2, amrex::RandomEngine const& engine) { return false; @@ -103,7 +105,7 @@ Collision::doCollisionB ( template void Collision::doCollisionImp ( - int lev, const amrex::Box& bx, const amrex::Geometry& geom, + int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma, F const& collision_function, G const& ionizaiton_function) @@ -117,8 +119,8 @@ Collision::doCollisionImp ( amrex::removeInvalidParticles(species2.ParticlesAt(0, pti)); } - PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom); - PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom); + PlasmaBins bins1 = findParticlesInEachTile(geom.Domain(), 1, species1, geom); + PlasmaBins bins2 = findParticlesInEachTile(geom.Domain(), 1, species2, geom); auto offset1 = bins1.offsetsPtr(); auto offset2 = bins2.offsetsPtr(); @@ -160,12 +162,45 @@ Collision::doCollisionImp ( const int np1 = ptile1.numParticles(); const int np2 = ptile2.numParticles(); - amrex::Gpu::DeviceVector flag1(np1, 0); - amrex::Gpu::DeviceVector flag2(np2, 0); + amrex::Gpu::DeviceVector flag1; + amrex::Gpu::DeviceVector flag2; + + if (m_has_collision_product) { + flag1.resize(np1, 0); + flag2.resize(np2, 0); + } auto p_flag1 = flag1.dataPtr(); auto p_flag2 = flag2.dataPtr(); + amrex::Gpu::DeviceVector cell_weight1(num_cells); + amrex::Gpu::DeviceVector cell_weight2(num_cells); + auto p_cell_weight1 = cell_weight1.dataPtr(); + auto p_cell_weight2 = cell_weight2.dataPtr(); + + amrex::ParallelFor(2*num_cells, + [=] AMREX_GPU_DEVICE (int icell) { + if (icell < num_cells) { + auto start = offset1[icell]; + auto stop = offset1[icell+1]; + amrex::Real loc_cell_weight1 = 0; + for (int idx1 = start; idx1 < stop; ++idx1) { + loc_cell_weight1 += ptd1.rdata(PlasmaIdx::w)[perm1[idx1]]; + } + p_cell_weight1[icell] = loc_cell_weight1; + } else { + icell -= num_cells; + auto start = offset2[icell]; + auto stop = offset2[icell+1]; + amrex::Real loc_cell_weight2 = 0; + for (int idx2 = start; idx2 < stop; ++idx2) { + loc_cell_weight2 += ptd2.rdata(PlasmaIdx::w)[perm2[idx2]]; + } + p_cell_weight2[icell] = loc_cell_weight2; + } + } + ); + amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); @@ -180,6 +215,9 @@ Collision::doCollisionImp ( const int n1 = offset1_stop - offset1_start; const int n2 = offset2_stop - offset2_start; + const amrex::Real loc_cell_weight1 = p_cell_weight1[icell]; + const amrex::Real loc_cell_weight2 = p_cell_weight2[icell]; + int idx1 = icoll; int idx2 = icoll; while (idx1 < n1 && idx2 < n2) { @@ -187,7 +225,7 @@ Collision::doCollisionImp ( const int j2 = perm2[offset2_start + idx2]; const bool make_new_particle = collision_function( - ptd1, j1, ptd2, j2, engine + ptd1, j1, ptd2, j2, loc_cell_weight1, loc_cell_weight2, engine ); if (n2 < n1) { From bcc5e0c0ab9679c8da129902dc0a819c28959eb9 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Tue, 17 Mar 2026 09:56:49 +0100 Subject: [PATCH 4/6] add same species collisions --- src/particles/collisions/Collision.H | 28 ++++ src/particles/collisions/Collision.cpp | 212 ++++++++++++++++++++++++- 2 files changed, 236 insertions(+), 4 deletions(-) diff --git a/src/particles/collisions/Collision.H b/src/particles/collisions/Collision.H index 16cc4e249e..250e3a2bb7 100644 --- a/src/particles/collisions/Collision.H +++ b/src/particles/collisions/Collision.H @@ -50,6 +50,34 @@ public: F const& collision_function, G const& ionizaiton_function); + /** + * \param[in] lev MR level + * \param[in] geom corresponding geometry object + * \param[in,out] multi_plasma all plasmas + * \param[in] collision_function return if a pair of particles collides + * \param[in] ionizaiton_function initialize product particle + **/ + template + void doCollisionImpSameSpecies ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionizaiton_function); + + /** + * \param[in] lev MR level + * \param[in] geom corresponding geometry object + * \param[in,out] multi_plasma all plasmas + * \param[in] collision_function return if a pair of particles collides + * \param[in] ionizaiton_function initialize product particle + **/ + template + void doCollisionImpDifferentSpecies ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionizaiton_function); + /** * \param[in] lev MR level * \param[in] geom corresponding geometry object diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index 99d42ca7ab..ee913f09b9 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -110,7 +110,208 @@ Collision::doCollisionImp ( F const& collision_function, G const& ionizaiton_function) { - // assume the two species are different for now + if (m_inout_species1_name == m_inout_species2_name) { + doCollisionImpSameSpecies(lev, geom, multi_plasma, + collision_function, + ionizaiton_function); + } else { + doCollisionImpDifferentSpecies(lev, geom, multi_plasma, + collision_function, + ionizaiton_function); + } +} + +template +void +Collision::doCollisionImpSameSpecies ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionizaiton_function) +{ + auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name); + + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + amrex::removeInvalidParticles(species1.ParticlesAt(0, pti)); + } + + PlasmaBins bins1 = findParticlesInEachTile(geom.Domain(), 1, species1, geom); + + auto offset1 = bins1.offsetsPtr(); + auto perm1 = bins1.permutationPtr(); + + const int num_cells = bins1.numBins(); + + amrex::Gpu::DeviceVector num_ind_pairs(num_cells, 0); + auto p_num_ind_pairs = num_ind_pairs.dataPtr(); + + const int total_ind_pairs = amrex::Scan::PrefixSum(num_cells, + [=] AMREX_GPU_DEVICE (int i) { + int n = offset1[i+1] - offset1[i]; + return n / 2; + }, + [=] AMREX_GPU_DEVICE (int i, int s) { + p_num_ind_pairs[i] = s; + }, + amrex::Scan::Type::exclusive, amrex::Scan::retSum + ); + + amrex::ParallelForRNG(num_cells, + [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) { + ShuffleFisherYates(perm1, offset1[i], offset1[i+1], engine); + } + ); + + for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { + + auto& ptile1 = species1.ParticlesAt(0, pti); + + auto ptd1 = ptile1.getParticleTileData(); + + const int np1 = ptile1.numParticles(); + + amrex::Gpu::DeviceVector flag1; + + if (m_has_collision_product) { + flag1.resize(np1, 0); + } + + auto p_flag1 = flag1.dataPtr(); + + amrex::Gpu::DeviceVector cell_weight1(num_cells); + amrex::Gpu::DeviceVector cell_weight2(num_cells); + auto p_cell_weight1 = cell_weight1.dataPtr(); + auto p_cell_weight2 = cell_weight2.dataPtr(); + + amrex::ParallelFor(2*num_cells, + [=] AMREX_GPU_DEVICE (int icell) { + + if (icell < num_cells) { + auto start = offset1[icell]; + auto stop = offset1[icell+1]; + auto mid = start + (stop - start + 1) / 2; + amrex::Real loc_cell_weight1 = 0; + for (int idx1 = start; idx1 < mid; ++idx1) { + loc_cell_weight1 += ptd1.rdata(PlasmaIdx::w)[perm1[idx1]]; + } + p_cell_weight1[icell] = loc_cell_weight1; + } else { + icell -= num_cells; + auto start = offset1[icell]; + auto stop = offset1[icell+1]; + auto mid = start + (stop - start + 1) / 2; + amrex::Real loc_cell_weight2 = 0; + for (int idx2 = mid; idx2 < stop; ++idx2) { + loc_cell_weight2 += ptd1.rdata(PlasmaIdx::w)[perm1[idx2]]; + } + p_cell_weight2[icell] = loc_cell_weight2; + } + } + ); + + amrex::ParallelForRNG(total_ind_pairs, + [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + + const int offset1_start = offset1[icell]; + const int offset1_stop = offset1[icell+1]; + const int mid = offset1_start + (offset1_stop - offset1_start + 1) / 2; + + const int icoll = ipair - p_num_ind_pairs[icell]; + + const int n1 = mid - offset1_start; + const int n2 = offset1_stop - mid; + + const amrex::Real loc_cell_weight1 = p_cell_weight1[icell]; + const amrex::Real loc_cell_weight2 = p_cell_weight2[icell]; + + int idx1 = icoll; + int idx2 = icoll; + while (idx1 < n1) { + const int j1 = perm1[offset1_start + idx1]; + const int j2 = perm1[mid + idx2]; + + const bool make_new_particle = collision_function( + ptd1, j1, ptd1, j2, loc_cell_weight1, loc_cell_weight2, engine + ); + + idx1 += n2; + if (make_new_particle) { + p_flag1[j1] = 1; + } + } + } + ); + + if (!m_has_collision_product) { + continue; + } + + const int num_new_particles = amrex::Scan::PrefixSum(np1, + [=] AMREX_GPU_DEVICE (int i) { + return p_flag1[i]; + }, + [=] AMREX_GPU_DEVICE (int i, int s) { + p_flag1[i] = p_flag1[i] == 0 ? -1 : s; + }, + amrex::Scan::Type::exclusive, amrex::Scan::retSum + ); + + auto& species3 = multi_plasma.GetPlasma(m_out_species3_name); + auto& ptile3 = species3.ParticlesAt(0, pti); + int old_size3 = ptile3.size(); + ptile3.resize(old_size3 + num_new_particles); + + // get new ptd after resize in case species3 is the same as species1 or 2 + ptd1 = ptile1.getParticleTileData(); + auto ptd3 = ptile3.getParticleTileData(); + + amrex::ParallelForRNG(total_ind_pairs, + [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + + const int offset1_start = offset1[icell]; + const int offset1_stop = offset1[icell+1]; + const int mid = offset1_start + (offset1_stop - offset1_start + 1) / 2; + + const int icoll = ipair - p_num_ind_pairs[icell]; + + const int n1 = mid - offset1_start; + const int n2 = offset1_stop - mid; + + int idx1 = icoll; + int idx2 = icoll; + while (idx1 < n1) { + const int j1 = perm1[offset1_start + idx1]; + const int j2 = perm1[mid + idx2]; + const int new_part_idx = p_flag1[j1]; + + if (new_part_idx >= 0) { + ionizaiton_function( + ptd1, j1, + ptd1, j2, + ptd3, old_size3 + new_part_idx, + engine + ); + } + + idx1 += n2; + } + } + ); + + amrex::Gpu::streamSynchronize(); + } +} + +template +void +Collision::doCollisionImpDifferentSpecies ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionizaiton_function) +{ auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name); auto& species2 = multi_plasma.GetPlasma(m_inout_species2_name); @@ -231,12 +432,12 @@ Collision::doCollisionImp ( if (n2 < n1) { idx1 += n2; if (make_new_particle) { - p_flag2[j2] = 1; + p_flag1[j1] = 1; } } else { idx2 += n1; if (make_new_particle) { - p_flag1[j1] = 1; + p_flag2[j2] = 1; } } } @@ -272,6 +473,9 @@ Collision::doCollisionImp ( int old_size3 = ptile3.size(); ptile3.resize(old_size3 + num_new_particles); + // get new ptd after resize in case species3 is the same as species1 or 2 + ptd1 = ptile1.getParticleTileData(); + ptd2 = ptile2.getParticleTileData(); auto ptd3 = ptile3.getParticleTileData(); amrex::ParallelForRNG(total_ind_pairs, @@ -293,7 +497,7 @@ Collision::doCollisionImp ( while (idx1 < n1 && idx2 < n2) { const int j1 = perm1[offset1_start + idx1]; const int j2 = perm2[offset2_start + idx2]; - const int new_part_idx = n2 < n1 ? p_flag2[j2] : p_flag1[j1]; + const int new_part_idx = n2 < n1 ? p_flag1[j2] : p_flag2[j1]; if (new_part_idx >= 0) { ionizaiton_function( From 1c6a159c9c3f553ab36042e1b373abbeac7ffbd9 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Tue, 17 Mar 2026 12:17:30 +0100 Subject: [PATCH 5/6] clean up UpdateMomentumPerezElastic --- .../collisions/UpdateMomentumPerez2.H | 251 ++++++++++++++++++ 1 file changed, 251 insertions(+) create mode 100644 src/particles/collisions/UpdateMomentumPerez2.H diff --git a/src/particles/collisions/UpdateMomentumPerez2.H b/src/particles/collisions/UpdateMomentumPerez2.H new file mode 100644 index 0000000000..ea97e7b400 --- /dev/null +++ b/src/particles/collisions/UpdateMomentumPerez2.H @@ -0,0 +1,251 @@ +/* Copyright 2019 Yinjian Zhao + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC2_H_ +#define HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC2_H_ + +#include "utils/Constants.H" + +#include + +#include // isnan() isinf() +#include // numeric_limits::min() + +/* \brief Update particle velocities according to + * F. Perez et al., Phys.Plasmas.19.083104 (2012), + * which is based on Nanbu's method, PhysRevE.55.4642 (1997). + * @param[in] LmdD is max(Debye length, minimal interparticle distance). + * @param[in] L is the Coulomb log. A fixed L will be used if L > 0, + * otherwise L will be calculated based on the algorithm. + * To see if there are nan or inf updated velocities, + * compile with USE_ASSERTION=TRUE. +*/ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdateMomentumPerezElastic2 ( + amrex::Real& u1x, amrex::Real& u1y, amrex::Real& u1z, const amrex::Real g1, + amrex::Real& u2x, amrex::Real& u2y, amrex::Real& u2z, const amrex::Real g2, + const amrex::Real n1, const amrex::Real n2, const amrex::Real n12, + const amrex::Real q1, amrex::Real m1, const amrex::Real w1, + const amrex::Real q2, amrex::Real m2, const amrex::Real w2, + const amrex::Real dt, const amrex::Real L, const amrex::Real lmdD, + const amrex::Real inv_c_SI, const amrex::Real inv_c2_SI, const bool normalized_units, + amrex::RandomEngine const& engine) +{ + using namespace amrex::literals; + + const amrex::Real diffx = amrex::Math::abs(u1x - u2x); + const amrex::Real diffy = amrex::Math::abs(u1y - u2y); + const amrex::Real diffz = amrex::Math::abs(u1z - u2z); + const amrex::Real diffm = std::sqrt(diffx * diffx + diffy * diffy + diffz * diffz); + const amrex::Real summm = std::sqrt(u1x * u1x + u1y * u1y + u1z * u1z) + + std::sqrt(u2x * u2x + u2y * u2y + u2z * u2z); + // If g = u1 - u2 = 0, do not collide. + // Or if the relative difference is less than 1.0e-10. + if (diffm < std::numeric_limits::min() || diffm < 1.0e-10 * summm) { + return; + } + + // Transform to SI units (u is the normalized momentum here) + if (normalized_units) { + m1 *= PhysConstSI::m_e; + m2 *= PhysConstSI::m_e; + } + // Compute momenta + const amrex::Real p1x = u1x * m1 * PhysConstSI::c; + const amrex::Real p1y = u1y * m1 * PhysConstSI::c; + const amrex::Real p1z = u1z * m1 * PhysConstSI::c; + const amrex::Real p2x = u2x * m2 * PhysConstSI::c; + const amrex::Real p2y = u2y * m2 * PhysConstSI::c; + const amrex::Real p2z = u2z * m2 * PhysConstSI::c; + + // Compute center-of-mass (COM) velocity and gamma + const amrex::Real mass_g = m1 * g1 + m2 * g2; + const amrex::Real mass_g_inv = 1._rt / mass_g; + const amrex::Real vcx = (p1x + p2x) * mass_g_inv; + const amrex::Real vcy = (p1y + p2y) * mass_g_inv; + const amrex::Real vcz = (p1z + p2z) * mass_g_inv; + const amrex::Real vcms = vcx * vcx + vcy * vcy + vcz * vcz; + const amrex::Real gc_inv = std::sqrt(std::max(1._rt - vcms * inv_c2_SI, 0._rt)); + if (gc_inv == 0._rt) { + return; + } + const amrex::Real gc = 1._rt / gc_inv; + + // Compute vc dot v1 and v2 + const amrex::Real vcDv1 = (vcx * u1x + vcy * u1y + vcz * u1z) * PhysConstSI::c / g1; + const amrex::Real vcDv2 = (vcx * u2x + vcy * u2y + vcz * u2z) * PhysConstSI::c / g2; + + // Compute p1 star + amrex::Real p1sx = p1x; + amrex::Real p1sy = p1y; + amrex::Real p1sz = p1z; + // If vcms = 0, don't do Lorentz-transform. + if (vcms > std::numeric_limits::min()) { + const amrex::Real lorentz_transform_factor = ((gc - 1._rt) / vcms * vcDv1 - gc) * m1 * g1; + p1sx += vcx * lorentz_transform_factor; + p1sy += vcy * lorentz_transform_factor; + p1sz += vcz * lorentz_transform_factor; + } + const double p1sm = std::sqrt(static_cast(p1sx) * p1sx + p1sy * p1sy + p1sz * p1sz); + const double p1sm_inv_sq = 1. / (p1sm * p1sm * inv_c2_SI); + + // Compute gamma star + const amrex::Real g1s = (1._rt - vcDv1 * inv_c2_SI) * gc * g1; + const amrex::Real g2s = (1._rt - vcDv2 * inv_c2_SI) * gc * g2; + + // Compute the Coulomb log lnLmd + double lnLmd = L; + if (lnLmd <= 0._rt) { + // Compute b0 according to eq (22) from Perez et al., + // Phys.Plasmas.19.083104 (2012) Note: there is a typo in the equation, + // the last square is incorrect! See the SMILEI documentation: + // https://smileipic.github.io/Smilei/Understand/collisions.html + const double b0 = amrex::Math::abs(q1 * q2) * inv_c2_SI * + (1._rt / (4._rt * MathConst::pi * PhysConstSI::ep0)) + * gc * mass_g_inv * (m1 * g1s * m2 * g2s * p1sm_inv_sq + 1._rt); + + // Compute the minimal impact parameter + double bmin_inv = 1. / std::max(PhysConstSI::hbar * MathConst::pi * (1. / p1sm), b0); + + // Compute the Coulomb log lnLmd according to eq (23) from Perez et al., + // Phys.Plasmas.19.083104 (2012) + lnLmd = amrex::max(static_cast(2.0), + 0.5_rt * std::log(1._rt + lmdD * lmdD * bmin_inv * bmin_inv)); + } + + // Compute s according to eq (17) from Perez et al., Phys.Plasmas.19.083104 + // (2012) + const double tts = static_cast(m1) * g1s * m2 * g2s * p1sm_inv_sq + 1._rt; + const double tts2 = tts * tts; + const double charge_fac = normalized_units ? static_cast(PhysConstSI::q_e) + * PhysConstSI::q_e * PhysConstSI::q_e * PhysConstSI::q_e : 1.; + amrex::Real s = static_cast(n1) * n2 / n12 * dt * lnLmd * q1 * q1 * + q2 * q2 * charge_fac * inv_c2_SI * inv_c2_SI / (static_cast(4.0) * MathConst::pi * + PhysConstSI::ep0 * PhysConstSI::ep0 * m1 * g1 * m2 * g2) * gc * p1sm * mass_g_inv * tts2; + + const auto cbrt_n1 = std::cbrt(n1); + const auto cbrt_n2 = std::cbrt(n2); + const auto coeff = static_cast(std::pow(4.0 * MathConst::pi / 3.0, 1.0 / 3.0)); + const amrex::Real vrel = static_cast(mass_g) * p1sm / (m1 * g1s * m2 * g2s * gc); + const amrex::Real sp = static_cast(coeff) * n1 * n2 / n12 * dt * vrel * (m1 + m2) / + amrex::max(m1 * cbrt_n1 * cbrt_n1, m2 * cbrt_n2 * cbrt_n2); + + // Determine s + s = amrex::min(s, sp); + + // Get random numbers + amrex::Real r = amrex::Random(engine); + + // Compute scattering angle according to eq (10) and following from Perez et + // al., Phys.Plasmas.19.083104 (2012) + amrex::Real cosXs = 0; + if (s <= 0.1_rt) { + while (true) { + cosXs = 1._rt + s * std::log(r); + // Avoid the bug when r is too small such that cosXs < -1 + if (cosXs >= -1._rt) { + break; + } + r = amrex::Random(engine); + } + } else if (s > 0.1_rt && s <= 3._rt) { + const amrex::Real Ainv = static_cast( + 0.0056958_rt + + 0.9560202_rt * s + - 0.508139_rt * s * s + + 0.47913906_rt * s * s * s + - 0.12788975_rt * s * s * s * s + + 0.02389567_rt * s * s * s * s * s); + cosXs = Ainv * std::log(std::exp(-(1._rt / Ainv)) + 2._rt * r * std::sinh(1._rt / Ainv)); + } else if (s > 3._rt && s <= 6._rt) { + const amrex::Real A = 3._rt * std::exp(-s); + cosXs = 1._rt / A * std::log(std::exp(-A) + 2._rt * r * std::sinh(A)); + } else { + cosXs = 2._rt * r - 1._rt; + } + const amrex::Real sinXs = std::sqrt(1._rt - cosXs * cosXs); + + // Get random azimuthal angle + const amrex::Real phis = amrex::Random(engine) * 2._rt * MathConst::pi; + const amrex::Real cosphis = std::cos(phis); + const amrex::Real sinphis = std::sin(phis); + + // Compute post-collision momenta pfs in COM + amrex::Real p1fsx = 0; + amrex::Real p1fsy = 0; + amrex::Real p1fsz = 0; + // p1sp is the p1s perpendicular + double p1sp = std::sqrt(static_cast(p1sx) * p1sx + p1sy * p1sy); + // Make sure p1sp is not almost zero + if (p1sp > std::numeric_limits::min()) { + const double p1sp_inv = 1. / p1sp; + p1fsx = (p1sx * p1sz * p1sp_inv) * sinXs * cosphis + + (p1sy * p1sm * p1sp_inv) * sinXs * sinphis + (p1sx)*cosXs; + p1fsy = (p1sy * p1sz * p1sp_inv) * sinXs * cosphis + + (-p1sx * p1sm * p1sp_inv) * sinXs * sinphis + (p1sy)*cosXs; + p1fsz = (-p1sp) * sinXs * cosphis + + (0._rt) * sinXs * sinphis + (p1sz)*cosXs; + // Note a negative sign is different from + // Eq. (12) in Perez's paper, + // but they are the same due to the random nature of phis. + } else { + // If the previous p1sp is almost zero + // x->y y->z z->x + // This set is equivalent to the one in Nanbu's paper + p1sp = std::sqrt(static_cast(p1sy) * p1sy + p1sz * p1sz); + const double p1sp_inv = 1. / p1sp; + p1fsy = (p1sy * p1sx * p1sp_inv) * sinXs * cosphis + + (p1sz * p1sm * p1sp_inv) * sinXs * sinphis + (p1sy)*cosXs; + p1fsz = (p1sz * p1sx * p1sp_inv) * sinXs * cosphis + + (-p1sy * p1sm * p1sp_inv) * sinXs * sinphis + (p1sz)*cosXs; + p1fsx = (-p1sp) * sinXs * cosphis + + (0._rt) * sinXs * sinphis + (p1sx)*cosXs; + } + + const amrex::Real p2fsx = -p1fsx; + const amrex::Real p2fsy = -p1fsy; + const amrex::Real p2fsz = -p1fsz; + + // Transform from COM to lab frame + amrex::Real p1fx = p1fsx; + amrex::Real p2fx = p1fsy; + amrex::Real p1fy = p1fsz; + amrex::Real p2fy = p2fsx; + amrex::Real p1fz = p2fsy; + amrex::Real p2fz = p2fsz; + // If vcms = 0, don't do Lorentz-transform. + if (vcms > std::numeric_limits::min()) { + const amrex::Real vcDp1fs = vcx * p1fsx + vcy * p1fsy + vcz * p1fsz; + const amrex::Real vcDp2fs = vcx * p2fsx + vcy * p2fsy + vcz * p2fsz; + const amrex::Real factor = (gc - 1._rt) / vcms; + const amrex::Real factor1 = factor * vcDp1fs + m1 * g1s * gc; + const amrex::Real factor2 = factor * vcDp2fs + m2 * g2s * gc; + p1fx += vcx * factor1; + p1fy += vcy * factor1; + p1fz += vcz * factor1; + p2fx += vcx * factor2; + p2fy += vcy * factor2; + p2fz += vcz * factor2; + } + + // Rejection method according to eq (14) and following from Perez et al., + // Phys.Plasmas.19.083104 (2012) + const amrex::Real m1_inv = 1._rt / m1; + const amrex::Real m2_inv = 1._rt / m2; + if (w2 > amrex::Random(engine) * amrex::max(w1, w2)) { + u1x = inv_c_SI * p1fx * m1_inv; + u1y = inv_c_SI * p1fy * m1_inv; + u1z = inv_c_SI * p1fz * m1_inv; + } + if (w1 > amrex::Random(engine) * amrex::max(w1, w2)) { + u2x = inv_c_SI * p2fx * m2_inv; + u2y = inv_c_SI * p2fy * m2_inv; + u2z = inv_c_SI * p2fz * m2_inv; + } +} + +#endif // HIPACE_UPDATE_MOMENTUM_PEREZ_ELASTIC2_H_ From 482a0b8de95b82f1b5adb2d49f3d3ddb9cde8603 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Fri, 10 Apr 2026 13:09:39 +0200 Subject: [PATCH 6/6] fix bug --- src/particles/collisions/Collision.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index ee913f09b9..e7dc8c0502 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -211,7 +211,7 @@ Collision::doCollisionImpSameSpecies ( amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ - const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells-1, ipair); const int offset1_start = offset1[icell]; const int offset1_stop = offset1[icell+1]; @@ -268,7 +268,7 @@ Collision::doCollisionImpSameSpecies ( amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ - const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells-1, ipair); const int offset1_start = offset1[icell]; const int offset1_stop = offset1[icell+1]; @@ -404,7 +404,7 @@ Collision::doCollisionImpDifferentSpecies ( amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ - const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells-1, ipair); const int offset1_start = offset1[icell]; const int offset1_stop = offset1[icell+1]; @@ -480,7 +480,7 @@ Collision::doCollisionImpDifferentSpecies ( amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ - const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells-1, ipair); const int offset1_start = offset1[icell]; const int offset1_stop = offset1[icell+1]; @@ -497,7 +497,7 @@ Collision::doCollisionImpDifferentSpecies ( while (idx1 < n1 && idx2 < n2) { const int j1 = perm1[offset1_start + idx1]; const int j2 = perm2[offset2_start + idx2]; - const int new_part_idx = n2 < n1 ? p_flag1[j2] : p_flag2[j1]; + const int new_part_idx = n2 < n1 ? p_flag1[j1] : p_flag2[j2]; if (new_part_idx >= 0) { ionizaiton_function(