From c52558d32206f25f5b6fa9d225fb35b12fc4cc92 Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Tue, 17 Feb 2026 16:34:17 +0100 Subject: [PATCH 1/8] 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 60a33e962c1304f222c647bada3aef251c5335bd Mon Sep 17 00:00:00 2001 From: Alexander Sinn Date: Fri, 20 Feb 2026 18:44:06 +0100 Subject: [PATCH 2/8] Bug search --- src/particles/collisions/Collision.H | 51 ++++- src/particles/collisions/Collision.cpp | 277 ++++++++++++++++++++++--- 2 files changed, 297 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(); } } From 74c41b99ef04c3c4d4b2cf60057f3f5df75fe879 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 23 Feb 2026 00:01:11 +0100 Subject: [PATCH 3/8] rebase lates updates --- src/Hipace.H | 3 ++ src/Hipace.cpp | 10 ++++ src/particles/collisions/Collision.H | 13 ++---- src/particles/collisions/Collision.cpp | 64 ++++++++++++++++++++------ 4 files changed, 68 insertions(+), 22 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index 4a38fb33ff..7c8d5c530f 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -23,6 +23,7 @@ #include "utils/MultiBuffer.H" #include "diagnostics/Diagnostic.H" #include "diagnostics/OpenPMDWriter.H" +#include "particles/collisions/Collision.H" #include #ifdef AMREX_USE_LINEAR_SOLVERS @@ -366,6 +367,8 @@ private: std::vector m_collision_names; /** 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, const amrex::Real time, const bool is_last_step); void FillBeamDiagnostics (const int step, const amrex::Real time, const bool is_last_step); diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 8e037b0ada..62126dfd13 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -269,6 +269,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); @@ -857,6 +863,10 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st // 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 52a5e6c709923e9a9a00ff322f1f350618717408 Mon Sep 17 00:00:00 2001 From: Emanuele Marchetti Date: Wed, 1 Apr 2026 14:40:44 +0200 Subject: [PATCH 4/8] WIP: impact ionization (electrons) --- src/particles/collisions/Collision.H | 8 +- src/particles/collisions/Collision.cpp | 366 +++++++++++++++--- src/particles/collisions/ImpactIonization.H | 361 +++++++++++++++++ .../collisions/ImpactIonizationSigma.H | 157 ++++++++ .../plasma/PlasmaParticleContainer.cpp | 7 +- src/utils/IonizationEnergiesTable.H | 2 +- 6 files changed, 840 insertions(+), 61 deletions(-) create mode 100644 src/particles/collisions/ImpactIonization.H create mode 100644 src/particles/collisions/ImpactIonizationSigma.H diff --git a/src/particles/collisions/Collision.H b/src/particles/collisions/Collision.H index 16cc4e249e..f3e5efa17e 100644 --- a/src/particles/collisions/Collision.H +++ b/src/particles/collisions/Collision.H @@ -8,8 +8,8 @@ #include /** - * \brief This class handles Coulomb collisions between 2 particle species - * (can be plasma-plasma or beam-plasma, the species can be the same) + * \brief This class handles inelastic Coulomb collisions between 2 plasma species, + * and eventually leads to an ionization event. */ class Collision { @@ -17,8 +17,10 @@ class Collision std::string m_inout_species2_name = ""; std::string m_out_species3_name = ""; std::string m_collision_type = ""; + std::string m_physical_element =""; bool m_has_collision_product = false; - amrex::Gpu::Buffer m_crosssection_data; + amrex::Gpu::Buffer m_binding_energies; + amrex::Gpu::Buffer m_ionization_energies; public: diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index 99d42ca7ab..206c418449 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -2,6 +2,51 @@ #include "Hipace.H" #include "particles/sorting/TileSort.H" #include "ShuffleFisherYates.H" +#include "utils/Constants.H" +#include "utils/IonizationEnergiesTable.H" +#include "utils/GPUUtil.H" +#include "ImpactIonization.H" +#include "ImpactIonizationSigma.H" + +using namespace amrex::literals; + +// Binding energies - reduced to H, N, He, Ar (eV) +// Innermost shell first, progressing outward +namespace { + constexpr amrex::Real bindingEnergies[28] = { + /* Z=1 (H) — 1 electron */ 13.5981, + /* Z=7 (N) — 7 electrons */ 403.78, 403.78, + 20.33, 20.33, + 14.534, 14.534, 14.534, + /* Z=2 (He) — 2 electrons */ 24.588, 24.588, + /* Z=18 (Ar) — 18 electrons */ 3206.2, 3206.2, + 324.2, 324.2, + 247.74, 247.74, 247.74, 247.74, 247.74, 247.74, + 29.24, 29.24, + 15.76, 15.76, 15.76, 15.76, 15.76, 15.76 + }; // Very structure specific to ImpactIonizationSigma.H + + constexpr amrex::Real ionizationEnergies[28] = { + /* Z=1 (H) */ 13.598434005136, + + /* Z=7 (N) : 1s, 1s, 2s, 2s, 2p, 2p, 2p */ + 667.04609, 552.06731, + 97.89013, 77.4735, + 47.4453, 29.60125, 14.53413, + + /* Z=2 (He) : 1s, 1s */ + 54.41776311, 24.587387936, + + /* Z=18 (Ar) : 1s,1s,2s,2s,2p x6,3s x2,3p x6 */ + 4426.2227, 4120.6655, + 918.374, 855.47, + 755.13, 685.47, 619.0, + 540.4, 479.76, 422.60, + 143.457, 124.41, + 91.290, 74.84, 59.58, + 40.735, 27.62967, 15.7596112 + }; +} void Collision::ReadParameters( @@ -13,12 +58,14 @@ Collision::ReadParameters( getWithParser(pp, "type", m_collision_type); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_collision_type == "a" || - m_collision_type == "b", + m_collision_type == "electron_impact" || + m_collision_type == "ion_impact", "Unknown collision type" ); - - m_has_collision_product = m_collision_type == "a"; + // The projectile can be either an electron or an ion + // Target particle is always a neutral/ion + // The user is asked to always input the target type as the second specie + m_has_collision_product = m_collision_type == "electron_impact" || m_collision_type == "ion_impact"; getWithParser(pp, "species1", m_inout_species1_name); getWithParser(pp, "species2", m_inout_species2_name); @@ -27,15 +74,23 @@ Collision::ReadParameters( getWithParser(pp, "species3", m_out_species3_name); } - - 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; + if (m_collision_type == "electron_impact") { + // Initialize binding energies for collision calculation + m_binding_energies.resize(28); + m_ionization_energies.resize(28); + for (int i=0; i<28 ; ++i) { + m_binding_energies[i] = bindingEnergies[i]; + m_ionization_energies[i] = ionizationEnergies[i]; } - m_crosssection_data.copyToDeviceAsync(); + m_binding_energies.copyToDeviceAsync(); + m_ionization_energies.copyToDeviceAsync(); + } + else { + } + // Plasma physical element name + amrex::ParmParse pp_s2(m_inout_species2_name); + getWithParser(pp_s2, "element", m_physical_element); } void @@ -43,46 +98,240 @@ Collision::doCollision ( int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma) { - if (m_collision_type == "a") { + if (m_collision_type == "electron_impact") { doCollisionA(lev, geom, multi_plasma); } else { doCollisionB(lev, geom, multi_plasma); } } - +// SI units are assumed (!!!) +// Different species colliding (electron impact) void Collision::doCollisionA ( int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma) { - auto p_crosssection_data = m_crosssection_data.data(); + constexpr amrex::Real c2 = PhysConstSI::c * PhysConstSI::c; + constexpr amrex::Real inv_c = 1.0_rt / PhysConstSI::c; + constexpr amrex::Real inv_c2 = 1.0_rt / (PhysConstSI::c * PhysConstSI::c); + + // auto p_crosssection_data = m_crosssection_data.data(); !!! + const amrex::Real* p_binding_energies = m_binding_energies.data(); + const amrex::Real* p_ionization_energies = m_ionization_energies.data(); + + auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name); + auto& species2 = multi_plasma.GetPlasma(m_inout_species2_name); + auto& species3 = multi_plasma.GetPlasma(m_out_species3_name); + + auto m1 = species1.GetMass(); + auto m2 = species2.GetMass(); + auto m3 = species3.GetMass(); + + const int ion_element_id = ion_map_ids[m_physical_element]; + const int ion_atomic_number = ion_atomic_numbers[ion_element_id]; + + const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2); + const amrex::Real dt = geom.CellSize(2) * inv_c; doCollisionImp(lev, geom, multi_plasma, [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, - amrex::Real cell_weight1, amrex::Real cell_weight2, + int N1, int N2, int icoll, 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]; + // Collision function + // clean from unused quantities + amrex::Real ux1 = ptd1.rdata(PlasmaIdx::ux_half_step)[i1]; + amrex::Real uy1 = ptd1.rdata(PlasmaIdx::uy_half_step)[i1]; + amrex::Real psi1 = ptd1.rdata(PlasmaIdx::psi_half_step)[i1]; + const amrex::Real w1 = ptd1.rdata(PlasmaIdx::w)[i1]; + + amrex::Real ux2 = ptd2.rdata(PlasmaIdx::ux_half_step)[i2]; + amrex::Real uy2 = ptd2.rdata(PlasmaIdx::uy_half_step)[i2]; + amrex::Real psi2 = ptd2.rdata(PlasmaIdx::psi_half_step)[i2]; + const amrex::Real w2 = ptd2.rdata(PlasmaIdx::w)[i2]; + int ion_lev2 = ptd2.idata(PlasmaIdx::ion_lev)[i2]; + + if (ion_lev2 > 0) { + // Already ionized, do not collide + return false; + } + + // particle's Lorentz factor + amrex::Real g1 = plasma_gamma(ux1, uy1, psi1, 1._rt / psi1, /* Assumes Aabssq == 0 */ 0._rt); + amrex::Real g2 = plasma_gamma(ux2, uy2, psi2, 1._rt / psi2, /* Assumes Aabssq == 0 */ 0._rt); + + // Convert from pseudo-potential to momentum + amrex::Real uz1 = plasma_uz(g1, psi1); + amrex::Real uz2 = plasma_uz(g2, psi2); - const int ion_lev = ptd1.idata(PlasmaIdx::ion_lev)[i1]; - const amrex::Real crosssection = p_crosssection_data[ion_lev]; + // 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 + const amrex::Real dt_fac = 0.5_rt * (g1/psi1 + g2/psi2); - return (crosssection + ux1) < ux2; + // Rescaling of the particle weights according to eq (22)-(23), + // Higginson et al., Journal of Computational Physics 413 (2020) + int N12 = amrex::max(N1,N2); + int D1; + int D2; + + if (N1>=N2) { + D1 = 1; + } else { + D1 = (N2/N1) + (icoll < (N2 % N1) ? 1 : 0); + } + + if (N2>=N1) { + D2 = 1; + } else { + D2 = (N1/N2) + (icoll < (N1 % N2) ? 1 : 0); + } + + amrex::Real w1r = w1 / amrex::max(D1,D2); + amrex::Real w2r = w2 / amrex::max(D1,D2); + + const amrex::Real diffx = amrex::Math::abs(ux1-ux2); + const amrex::Real diffy = amrex::Math::abs(uy1-uy2); + const amrex::Real diffz = amrex::Math::abs(uz1-uz2); + const amrex::Real diffm = std::sqrt(diffx*diffx+diffy*diffy+diffz*diffz); + const amrex::Real summ = std::sqrt(ux1*ux1+uy1*uy1+uz1*uz1) + std::sqrt(ux2*ux2+uy2*uy2+uz2*uz2); + // 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/summ < 1.0e-10 ) { + return false; + } + // Compute the particles relative velocity (lab frame) + const amrex::Real u1u2 = ux1*ux2+uy1*uy2+uz1*uz2; + const amrex::Real grel = g1*g2-u1u2; + const amrex::Real inv_grel2 = 1/(grel*grel); + const amrex::Real vrel = PhysConstSI::c * std::sqrt(1-inv_grel2); + + // Fetch the cross section + const amrex::Real Krel = (grel - 1.0) * m1 * c2 / PhysConstSI::q_e; // (eV) + auto cs_data = coll_ion::sigma_E(ion_atomic_number, ion_lev2, Krel, p_binding_energies, p_ionization_energies); + auto sigma = cs_data.sigma; + + // The ionization mechanism is based on the algorithm from Perez et al., Phys.Plasmas.19.083104 (2012) + // The weights are rescaled according to eq (22)-(23), + // Higginson et al., Journal of Computational Physics 413 (2020) + + // Ionization probability + const auto Pion = 1 - std::exp(-vrel * N12 * amrex::max(w1r,w2r) * inv_dV * sigma * dt*dt_fac); + // Get random numbers + auto r = amrex::Random(engine); + if (Pion > r) { + // Rejection method according to the adaptation of eq (14) in the ionization model, + // from Perez et al., Phys.Plasmas.19.083104 (2012) + r = amrex::Random(engine); + return ( w1r > r*amrex::max(w1r,w2r) ); + } else { + return false; + } }, [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, auto ptd3, int i3, + int N1, int N2, int icoll, 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; + // Ionization function + // Ionization occurs if the condition on Pion is satisfied + // A new electron is created following the rejection method + amrex::Real ux1 = ptd1.rdata(PlasmaIdx::ux_half_step)[i1]; + amrex::Real uy1 = ptd1.rdata(PlasmaIdx::uy_half_step)[i1]; + amrex::Real psi1 = ptd1.rdata(PlasmaIdx::psi_half_step)[i1]; + const amrex::Real w1 = ptd1.rdata(PlasmaIdx::w)[i1]; + const int ion_lev1 = ptd1.idata(PlasmaIdx::ion_lev)[i1]; + + amrex::Real ux2 = ptd2.rdata(PlasmaIdx::ux_half_step)[i2]; + amrex::Real uy2 = ptd2.rdata(PlasmaIdx::uy_half_step)[i2]; + amrex::Real psi2 = ptd2.rdata(PlasmaIdx::psi_half_step)[i2]; + const amrex::Real w2 = ptd2.rdata(PlasmaIdx::w)[i2]; + int ion_lev2 = ptd2.idata(PlasmaIdx::ion_lev)[i2]; + + // Initialize new electron properties + ptd3.id(i3) = 4; + ptd3.cpu(i3) = 0; + ptd3.rdata(PlasmaIdx::x)[i3] = ptd2.rdata(PlasmaIdx::x)[i2]; + ptd3.rdata(PlasmaIdx::y)[i3] = ptd2.rdata(PlasmaIdx::y)[i2]; + ptd3.rdata(PlasmaIdx::x_prev)[i3] = ptd2.rdata(PlasmaIdx::x_prev)[i2]; + ptd3.rdata(PlasmaIdx::y_prev)[i3] = ptd2.rdata(PlasmaIdx::y_prev)[i2]; + ptd3.rdata(PlasmaIdx::w)[i3] = ptd2.rdata(PlasmaIdx::w)[i2]; + ptd3.idata(PlasmaIdx::ion_lev)[i3] = 0; + // Quantities to be updated after collision + amrex::Real ux3 = 0._rt; + amrex::Real uy3 = 0._rt; + amrex::Real psi3 = 1._rt; + amrex::Real uz3 = 0._rt; + ptd3.rdata(PlasmaIdx::ux)[i3] = ux3; + ptd3.rdata(PlasmaIdx::uy)[i3] = uy3; + ptd3.rdata(PlasmaIdx::psi)[i3] = psi3; + ptd3.rdata(PlasmaIdx::ux_half_step)[i3] = ux3; + ptd3.rdata(PlasmaIdx::uy_half_step)[i3] = uy3; + ptd3.rdata(PlasmaIdx::psi_half_step)[i3] = psi3; + + // particle's Lorentz factor + amrex::Real g1 = plasma_gamma(ux1, uy1, psi1, 1._rt / psi1, /* Assumes Aabssq == 0 */ 0._rt); + amrex::Real g2 = plasma_gamma(ux2, uy2, psi2, 1._rt / psi2, /* Assumes Aabssq == 0 */ 0._rt); + + // Convert from pseudo-potential to momentum + amrex::Real uz1 = plasma_uz(g1, psi1); + amrex::Real uz2 = plasma_uz(g2, psi2); + + // Fetch the ionization energy + const amrex::Real u1u2 = ux1*ux2+uy1*uy2+uz1*uz2; + const amrex::Real grel = g1*g2-u1u2; + const amrex::Real Krel = (grel - 1.0) * m1 * c2 / PhysConstSI::q_e; // (eV) + auto cs_data = coll_ion::sigma_E(ion_atomic_number, ion_lev2, Krel, p_binding_energies, p_ionization_energies); + auto Eion_eV = cs_data.ionization_en; + + // Rescaling of the particle weights according to eq (22)-(23), + // Higginson et al., Journal of Computational Physics 413 (2020) + int N12 = amrex::max(N1,N2); + int D1; + int D2; + + if (N1>=N2) { + D1 = 1; + } else { + D1 = (N2/N1) + (icoll < (N2 % N1) ? 1 : 0); + } + + if (N2>=N1) { + D2 = 1; + } else { + D2 = (N1/N2) + (icoll < (N1 % N2) ? 1 : 0); + } + + amrex::Real w1r = w1 / amrex::max(D1,D2); + amrex::Real w2r = w2 / amrex::max(D1,D2); + + ImpactIonization( + ux1, uy1, uz1, g1, + ux2, uy2, uz2, g2, ion_lev2, + ux3, uy3, uz3, + m1, w1r, m2, w2r, m3, + Eion_eV, c2, inv_c, inv_c2, + engine + ); + + // Update particle properties + ptd1.rdata(PlasmaIdx::ux_half_step)[i1] = ux1; + ptd1.rdata(PlasmaIdx::uy_half_step)[i1] = uy1; + ptd1.rdata(PlasmaIdx::psi_half_step)[i1] = plasma_psi(ux1, uy1, uz1, /* Assumes Aabssq == 0 */ 0._rt); + + ptd2.rdata(PlasmaIdx::ux_half_step)[i2] = ux2; + ptd2.rdata(PlasmaIdx::uy_half_step)[i2] = uy2; + ptd2.rdata(PlasmaIdx::psi_half_step)[i2] = plasma_psi(ux2, uy2, uz2, /* Assumes Aabssq == 0 */ 0._rt); + ptd2.idata(PlasmaIdx::ion_lev)[i2] = ion_lev2; + + ptd3.rdata(PlasmaIdx::ux)[i3] = ux3; + ptd3.rdata(PlasmaIdx::uy)[i3] = uy3; + ptd3.rdata(PlasmaIdx::psi)[i3] = plasma_psi(ux3, uy3, uz3, /* Assumes Aabssq == 0 */ 0._rt); + ptd3.rdata(PlasmaIdx::ux_half_step)[i3] = ux3; + ptd3.rdata(PlasmaIdx::uy_half_step)[i3] = uy3; + ptd3.rdata(PlasmaIdx::psi_half_step)[i3] = plasma_psi(ux3, uy3, uz3, /* Assumes Aabssq == 0 */ 0._rt); }); } +// Same species colliding (Ion impact) void Collision::doCollisionB ( int lev, const amrex::Geometry& geom, @@ -90,12 +339,13 @@ Collision::doCollisionB ( { doCollisionImp(lev, geom, multi_plasma, [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, - amrex::Real cell_weight1, amrex::Real cell_weight2, + int N1, int N2, int icoll, amrex::RandomEngine const& engine) { return false; }, [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, auto ptd3, int i3, + int N1, int N2, int icoll, amrex::RandomEngine const& engine) { @@ -108,7 +358,7 @@ Collision::doCollisionImp ( int lev, const amrex::Geometry& geom, MultiPlasma& multi_plasma, F const& collision_function, - G const& ionizaiton_function) + G const& ionization_function) { // assume the two species are different for now auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name); @@ -122,6 +372,8 @@ Collision::doCollisionImp ( PlasmaBins bins1 = findParticlesInEachTile(geom.Domain(), 1, species1, geom); PlasmaBins bins2 = findParticlesInEachTile(geom.Domain(), 1, species2, geom); + // offset: start/end positions of particles for each cell + // perm: permutation array mapping cell entries to particle indices auto offset1 = bins1.offsetsPtr(); auto offset2 = bins2.offsetsPtr(); auto perm1 = bins1.permutationPtr(); @@ -134,9 +386,9 @@ Collision::doCollisionImp ( 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); + 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; @@ -151,6 +403,7 @@ Collision::doCollisionImp ( } ); + // loop over tiles for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { auto& ptile1 = species1.ParticlesAt(0, pti); @@ -201,6 +454,7 @@ Collision::doCollisionImp ( } ); + // loop over independent pairs 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); @@ -212,31 +466,33 @@ Collision::doCollisionImp ( const int icoll = ipair - p_num_ind_pairs[icell]; - const int n1 = offset1_stop - offset1_start; - const int n2 = offset2_stop - offset2_start; + 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]; + //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) { + 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, loc_cell_weight1, loc_cell_weight2, engine + ptd1, j1, ptd2, j2, + N1, N2, icoll, + engine ); - if (n2 < n1) { - idx1 += n2; + if (N2 < N1) { + idx1 += N2; if (make_new_particle) { - p_flag2[j2] = 1; + p_flag1[j1] = 1; } } else { - idx2 += n1; + idx2 += N1; if (make_new_particle) { - p_flag1[j1] = 1; + p_flag2[j2] = 1; } } } @@ -272,6 +528,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, @@ -285,29 +544,28 @@ Collision::doCollisionImp ( const int icoll = ipair - p_num_ind_pairs[icell]; - const int n1 = offset1_stop - offset1_start; - const int n2 = offset2_stop - offset2_start; + 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) { + 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[j1] : p_flag2[j2]; if (new_part_idx >= 0) { - ionizaiton_function( - ptd1, j1, - ptd2, j2, - ptd3, old_size3 + new_part_idx, + ionization_function( + ptd1, j1, ptd2, j2, ptd3, old_size3 + new_part_idx, + N1, N2, icoll, engine ); } - if (n2 < n1) { - idx1 += n2; + if (N2 < N1) { + idx1 += N2; } else { - idx2 += n1; + idx2 += N1; } } } diff --git a/src/particles/collisions/ImpactIonization.H b/src/particles/collisions/ImpactIonization.H new file mode 100644 index 0000000000..91c3d6c55a --- /dev/null +++ b/src/particles/collisions/ImpactIonization.H @@ -0,0 +1,361 @@ +/* This file has been derived from + * https://github.com/BLAST-WarpX/warpx/blob/development/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H + * + * Copyright 2023-2024 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Roelof Groenewald (TAE Technologies) + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef HIPACE_IMPACT_IONIZATION_H_ +#define HIPACE_IMPACT_IONIZATION_H_ + +#include "utils/Constants.H" + +#include + +#include // isnan() isinf() +#include // numeric_limits::min() + +/* \brief The implicit, non-relativistic assumption has not been + * changed + * Momentum and energy division after the ionization event + * is done as follows: + * Three numbers are generated that satisfy the triangle inequality. + * These numbers will be scaled to give the momentum + * for each particle (satisfying the triangle inequality + * ensures that the momentum vectors can be arranged + * such that the net linear momentum is zero). + * Species 2 is definitely a neutral/ion, so we first choose its + * proportion to be n_2 = min(E2 / Ecoll, 0.5 * R), where + * R is a random number between 0 and 1. + * The other two numbers (n_1 & n_3) are then found + * by choosing a random point on the ellipse characterized + * by the semi-major and semi-minor axes + * a = (1 - n_2) / 2.0 + * b = 0.5 * sqrt(1 - 2 * n_2) + * The numbers are found by randomly sampling an x value + * between -a and a, and finding the corresponding y value + * that falls on the ellipse: y^2 = b^2 - b^2/a^2 * x^2. + * Then n_1 = sqrt(y^2 + (x - n_2/2)^2) and + * n_3 = 1 - n_2 - n_1. + * Next, we need to find the number C, such that if + * p_i = C * n_i, we would have E_1 + E_2 + E_3 = E_out + * where E_i = p_i^2 / (2 m_i), i.e., + * C^2 * \sum_i [n_i^2 / (2 m_i)] = E_out + * After C is determined the momentum vectors are arranged + * such that the net linear momentum is zero. + * To see if there are nan or inf updated velocities, + * compile with USE_ASSERTION=TRUE. +*/ + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ImpactIonization ( + T_R& ux1, T_R& uy1, T_R& uz1, T_R const g1, + T_R& ux2, T_R& uy2, T_R& uz2, T_R const g2, int& ion_lev2, + T_R& ux3, T_R& uy3, T_R& uz3, + T_R m1, T_R const w1r, T_R m2, T_R const w2r, T_R const m3, + T_R const Eion_eV, T_R const c2, T_R const inv_c, T_R const inv_c2, + amrex::RandomEngine const& engine) +{ + using namespace amrex::literals; + + // Helper lambda to clamp a value between a lower and upper bound + auto clamp = [] AMREX_GPU_HOST_DEVICE (T_R x, T_R lo, T_R hi) -> T_R + { + return amrex::max(lo, amrex::min(hi, x)); + }; + // Tolerances for numerical checks + constexpr T_R rel_tol = T_R(1.0e-10); + constexpr T_R eps_n = T_R(1.0e-14); + constexpr T_R beta_tol = T_R(1.0e-12); + + T_R const Eion = Eion_eV * PhysConstSI::q_e; // [J] + + // ------------------------------------------------------------------ + // Early-out if relative motion is negligible + // ------------------------------------------------------------------ + T_R const diffx = ux1 - ux2; + T_R const diffy = uy1 - uy2; + T_R const diffz = uz1 - uz2; + T_R const diffm = std::sqrt(diffx*diffx + diffy*diffy + diffz*diffz); + T_R const summ = std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1) + std::sqrt(ux2*ux2 + uy2*uy2 + uz2*uz2); + // If g = u1 - u2 = 0, do not collide. + // Or if the relative difference is less than 1.0e-10. + if (summ <= T_R(0.0)) return; + if (diffm / summ < rel_tol) return; + + // ------------------------------------------------------------------ + // Initial lab-frame momenta + // ------------------------------------------------------------------ + T_R const px1 = ux1 * m1 * PhysConstSI::c; + T_R const py1 = uy1 * m1 * PhysConstSI::c; + T_R const pz1 = uz1 * m1 * PhysConstSI::c; + + T_R const px2 = ux2 * m2 * PhysConstSI::c; + T_R const py2 = uy2 * m2 * PhysConstSI::c; + T_R const pz2 = uz2 * m2 * PhysConstSI::c; + + // ------------------------------------------------------------------ + // COM velocity + // ------------------------------------------------------------------ + T_R const mass_g = m1 * g1 + m2 * g2; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mass_g > T_R(0.0), "Mass is zero in ImpactIonization"); + + T_R const vcx = (px1+px2) / mass_g; + T_R const vcy = (py1+py2) / mass_g; + T_R const vcz = (pz1+pz2) / mass_g; + T_R const vcms = vcx*vcx + vcy*vcy + vcz*vcz; + + T_R beta2 = vcms * inv_c2; + if (beta2 < T_R(0.0)) { + if (beta2 > -beta_tol) { + beta2 = T_R(0.0); + } else { + AMREX_DEVICE_PRINTF( + "ImpactIonization: beta < 0, beta2=%e vcms=%e\n", + beta2, vcms + ); + return; + } + } + + if (beta2 >= T_R(1.0) - beta_tol) { + AMREX_DEVICE_PRINTF( + "ImpactIonization: beta too close to 1, beta2=%e\n", + beta2 + ); + return; + } + + T_R const gc = T_R(1.0) / std::sqrt( T_R(1.0) - beta2 ); + + // Stable equivalent of (gc - 1)/vcms: + // (gamma - 1)/v^2 = gamma^2 / ((gamma + 1) c^2) + T_R const boost_v2 = gc*gc / ((gc + T_R(1.0)) * c2); + + // ------------------------------------------------------------------ + // COM-frame incoming particle 1 momentum + // ------------------------------------------------------------------ + T_R const vcDv1 = (vcx*ux1 + vcy*uy1 + vcz*uz1) * PhysConstSI::c / g1; + T_R const vcDv2 = (vcx*ux2 + vcy*uy2 + vcz*uz2) * PhysConstSI::c / g2; + + // Compute p1 star + T_R p1sx; + T_R p1sy; + T_R p1sz; + if ( vcms > std::numeric_limits::min() ) + { + T_R const lorentz_transform_factor = + ( boost_v2 * vcDv1 - gc ) * m1 * g1; + + p1sx = px1 + vcx*lorentz_transform_factor; + p1sy = py1 + vcy*lorentz_transform_factor; + p1sz = pz1 + vcz*lorentz_transform_factor; + } + else // If vcms = 0, don't do Lorentz-transform. + { + p1sx = px1; + p1sy = py1; + p1sz = pz1; + } + + // ------------------------------------------------------------------ + // COM-frame gammas + // ------------------------------------------------------------------ + T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2 )*gc*g1; + T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2 )*gc*g2; + // Redundant? + if (g1s < T_R(1.0) || g2s < T_R(1.0)) { + AMREX_DEVICE_PRINTF( + "ImpactIonization: invalid COM gammas, g1s=%e g2s=%e\n", + g1s, g2s + ); + return; + } + + // Ionization probability and final-state construction + // Available initial kinetic energy in the COM frame [J] + T_R const Ecm = + ( (g1s - T_R(1.0))*m1 + (g2s - T_R(1.0))*m2 ) * c2; + + // Kinetic energy available to the 3 outgoing particles in COM [J] + T_R const Eout = Ecm - Eion; + if (Eout <= T_R(0.0)) return; // Energy not sufficient for ionization + + // 3-body final state in the COM frame + // Particle 1 : scattered projectile (mass m1) + // Particle 2 : ion / neutral (mass m2) + // Particle 3 : secondary electron (mass m3) + // + // Momentum magnitudes satisfy triangle closure, + // scaled so that total kinetic energy = Eout. + + T_R const E_target = (g2s-T_R(1.0))*m2 * c2; + // Sample three normalized momentum magnitudes satisfying triangle inequality + T_R n_2 = amrex::min(E_target/Ecm, T_R(0.5)*amrex::Random(engine)); + n_2 = clamp(n_2, T_R(0.0), T_R(0.5) - eps_n); + + T_R const a = T_R(0.5) * (T_R(1.0) - n_2); + + T_R const b2_arg = T_R(1.0) - T_R(2.0)*n_2; + if (b2_arg < T_R(0.0)) { + AMREX_DEVICE_PRINTF( + "ImpactIonization: invalid ellipse, b2_arg=%e n2=%e E_target=%e Ecm=%e\n", + b2_arg, n_2, E_target, Ecm + ); + return; + } + + T_R const b = T_R(0.5) * std::sqrt(amrex::max(T_R(0.0), b2_arg)); + + T_R const x = (T_R(2.0)*amrex::Random(engine) - T_R(1.0)) * a; + T_R y2 = b*b - (b*b/(a*a))*x*x; + y2 = amrex::max(T_R(0.0), y2); + + T_R n1sq = y2 + x*x - x*n_2 + T_R(0.25)*n_2*n_2; + n1sq = amrex::max(T_R(0.0), n1sq); + + T_R const n_1 = std::sqrt(n1sq); + T_R const n_3 = T_R(1.0) - n_1 - n_2; + if (n_3 < T_R(0.0)) { + AMREX_DEVICE_PRINTF( + "ImpactIonization: invalid triangle, n1=%e n2=%e n3=%e\n", + n_1, n_2, n_3 + ); + return; + } + + // Common momentum scale from energy conservation (non-relativistic!) + T_R const C = std::sqrt(Eout / ( + n_1*n_1 / (T_R(2.0)*m1) + n_2*n_2 / (T_R(2.0)*m2) + n_3*n_3 / (T_R(2.0)*m3)) + ); + + // Momentum magnitudes in COM + T_R const p1f_mag = C * n_1; + T_R const p2f_mag = C * n_2; + T_R const p3f_mag = C * n_3; + + // Triangle angles from the law of cosines + T_R cos_alpha = (n_1*n_1 + n_3*n_3 - n_2*n_2) / (T_R(2.0)*n_1*n_3); + cos_alpha = clamp(cos_alpha, T_R(-1.0), T_R(1.0)); + + T_R cos_gamma = (n_1*n_1 + n_2*n_2 - n_3*n_3) / (T_R(2.0)*n_1*n_2); + cos_gamma = clamp(cos_gamma, T_R(-1.0), T_R(1.0)); + + T_R const sin_alpha = std::sqrt(amrex::max(T_R(0.0), T_R(1.0) - cos_alpha*cos_alpha)); + T_R const sin_gamma = std::sqrt(amrex::max(T_R(0.0), T_R(1.0) - cos_gamma*cos_gamma)); + + // Random orientation of the momentum triangle in 3D + T_R const Theta = amrex::Random(engine) * T_R(2.0) * MathConst::pi; + T_R const phi = amrex::Random(engine) * MathConst::pi; + + T_R const cos_theta = std::cos(Theta); + T_R const sin_theta = std::sin(Theta); + T_R const cos_phi = std::cos(phi); + T_R const sin_phi = std::sin(phi); + + // Outgoing momenta in the COM frame + T_R p1fsx = p1f_mag * cos_theta * cos_phi; + T_R p1fsy = p1f_mag * cos_theta * sin_phi; + T_R p1fsz = -p1f_mag * sin_theta; + + T_R p2fsx = p2f_mag * (-cos_alpha*cos_theta*cos_phi - sin_alpha*sin_phi); + T_R p2fsy = p2f_mag * (-cos_alpha*cos_theta*sin_phi + sin_alpha*cos_phi); + T_R p2fsz = p2f_mag * ( cos_alpha*sin_theta ); + + T_R p3fsx = p3f_mag * (-cos_gamma*cos_theta*cos_phi + sin_gamma*sin_phi); + T_R p3fsy = p3f_mag * (-cos_gamma*cos_theta*sin_phi - sin_gamma*cos_phi); + T_R p3fsz = p3f_mag * ( cos_gamma*sin_theta ); + + // Boost final-state momenta from COM back to lab + T_R const g1f_s = std::sqrt(T_R(1.0) + (p1fsx*p1fsx + p1fsy*p1fsy + p1fsz*p1fsz)/(m1*m1*c2)); + T_R const g2f_s = std::sqrt(T_R(1.0) + (p2fsx*p2fsx + p2fsy*p2fsy + p2fsz*p2fsz)/(m2*m2*c2)); + T_R const g3f_s = std::sqrt(T_R(1.0) + (p3fsx*p3fsx + p3fsy*p3fsy + p3fsz*p3fsz)/(m3*m3*c2)); + + T_R p1fx; T_R p1fy; T_R p1fz; + T_R p2fx; T_R p2fy; T_R p2fz; + T_R p3fx; T_R p3fy; T_R p3fz; + + if ( vcms > std::numeric_limits::min() ) + { + T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz; + T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz; + T_R const vcDp3fs = vcx*p3fsx + vcy*p3fsy + vcz*p3fsz; + + T_R const factor1 = boost_v2*vcDp1fs + m1*g1f_s*gc; + T_R const factor2 = boost_v2*vcDp2fs + m2*g2f_s*gc; + T_R const factor3 = boost_v2*vcDp3fs + m3*g3f_s*gc; + + p1fx = p1fsx + vcx*factor1; + p1fy = p1fsy + vcy*factor1; + p1fz = p1fsz + vcz*factor1; + + p2fx = p2fsx + vcx*factor2; + p2fy = p2fsy + vcy*factor2; + p2fz = p2fsz + vcz*factor2; + + p3fx = p3fsx + vcx*factor3; + p3fy = p3fsy + vcy*factor3; + p3fz = p3fsz + vcz*factor3; + } + else // If vcms = 0, don't do Lorentz-transform + { + p1fx = p1fsx; p1fy = p1fsy; p1fz = p1fsz; + p2fx = p2fsx; p2fy = p2fsy; p2fz = p2fsz; + p3fx = p3fsx; p3fy = p3fsy; p3fz = p3fsz; + } + + // Convert lab-frame momenta back to normalized momenta + // Rejection method according to the adaptation of eq (14) in the ionization model, + // from Perez et al., Phys.Plasmas.19.083104 (2012) + // The weights are rescaled according to eq (22)-(23), + // Higginson et al., Journal of Computational Physics 413 (2020) + auto r = amrex::Random(engine); + if ( w2r > r*amrex::max(w1r, w2r) ) + { + ux1 = inv_c * p1fx / m1; + uy1 = inv_c * p1fy / m1; + uz1 = inv_c * p1fz / m1; + + bool var_nan = + std::isnan(ux1) || std::isnan(uy1) || std::isnan(uz1); + + bool var_inf = + std::isinf(ux1) || std::isinf(uy1) || std::isinf(uz1); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!var_nan, "ImpactIonization: NaN detected in final-state velocities"); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!var_inf, "ImpactIonization: Inf detected in final-state velocities"); + } + + ion_lev2 += 1; + + ux2 = inv_c * p2fx / m2; + uy2 = inv_c * p2fy / m2; + uz2 = inv_c * p2fz / m2; + + ux3 = inv_c * p3fx / m3; + uy3 = inv_c * p3fy / m3; + uz3 = inv_c * p3fz / m3; + + bool var_nan = + std::isnan(ux1) || std::isnan(uy1) || std::isnan(uz1) || + std::isnan(ux2) || std::isnan(uy2) || std::isnan(uz2) || + std::isnan(ux3) || std::isnan(uy3) || std::isnan(uz3); + + bool var_inf = + std::isinf(ux1) || std::isinf(uy1) || std::isinf(uz1) || + std::isinf(ux2) || std::isinf(uy2) || std::isinf(uz2) || + std::isinf(ux3) || std::isinf(uy3) || std::isinf(uz3); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!var_nan, "ImpactIonization: NaN detected in final-state velocities"); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!var_inf, "ImpactIonization: Inf detected in final-state velocities"); +} + +#endif // HIPACE_IMPACT_IONIZATION_H_ \ No newline at end of file diff --git a/src/particles/collisions/ImpactIonizationSigma.H b/src/particles/collisions/ImpactIonizationSigma.H new file mode 100644 index 0000000000..c82543da40 --- /dev/null +++ b/src/particles/collisions/ImpactIonizationSigma.H @@ -0,0 +1,157 @@ +#ifndef HIPACE_IMPACT_IONIZATION_SIGMA_H_ +#define HIPACE_IMPACT_IONIZATION_SIGMA_H_ + +#include +#include +// #include +#include "utils/Constants.H" + +// Quantum degeneracy for each species (g = 2J + 1, ground-state term) +/* +AMREX_GPU_HOST_DEVICE inline amrex::Real g_deg[29] = { + // H (2 stages) + 2.0, 1.0, + // N (6 stages) + 4.0, 9.0, 6.0, 1.0, 2.0, 1.0, + // He (3 stages) + 1.0, 2.0, 1.0, + // Ar (18 stages) + 1.0, 4.0, 5.0, 4.0, 9.0, 6.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0, 1.0 +}; + +constexpr int num_atom_species = 4; + +constexpr amrex::Real X_min = 1e-20; // Threshold for considering a species to be present + +// Helper lambda to map species index to atomic species +AMREX_GPU_HOST_DEVICE +inline int get_atom_index(int mu) { + if (mu <= 1 && mu >= 0) return 0; // Hydrogen + if (mu <= 7 && mu > 1) return 1; // Nitrogen + if (mu <= 10 && mu > 7) return 2; // Helium + if (mu <= 14 && mu > 10) return 3; // Argon + return 0; //!!! +}; + +AMREX_GPU_HOST_DEVICE +const int zatom[num_atom_species] = {1, 7, 2, 18}; // Atomic numbers for H, N, He, Ar +*/ + +// Index mapping for reduced arrays +// Returns the starting index in the energy arrays for a given Z +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int get_energy_index(int Z) +{ + if (Z == 1) return 0; // H — 1 entry (indices 0) + if (Z == 7) return 1; // N — 7 entries (indices 1–7) + if (Z == 2) return 8; // He — 2 entries (indices 8–9) + if (Z == 18) return 10; // Ar — 18 entries (indices 10–27) + + return 0; //!!! +} + +namespace coll_ion +{ + using namespace amrex; + + inline constexpr Real alpha = 7.2973525643e-3; + inline constexpr Real a0 = 5.29177210544e-11; // m + inline constexpr Real normalization = 4.0 * MathConst::pi * a0*a0 * alpha*alpha*alpha*alpha; + + struct IonizationResult + { + Real sigma; + Real ionization_en; + }; + + // Cross section using RBEB model from Kim et al. (2000) + /** + * @brief Relativistic Binary-Encounter-Bethe (RBEB) electron-impact + * ionization cross section. + * + * Loops over all bound electrons of the ion (from innermost to the + * outermost occupied shell). Degenerate sub-shells sharing the same + * binding energy are grouped and counted together (factor `n`). + * + * @param atomic_number Atomic number Z of the target element. + * @param zstar Current ion charge (number of removed electrons). + * @param e Incident-electron kinetic energy [eV]. + * @return Ionization cross section [m²]. + */ + AMREX_GPU_HOST_DEVICE inline + IonizationResult sigma_E(int atomic_number, int zstar, Real e, const Real* p_binding_energies, const Real* p_ionization_energies) + { + Real sigma = 0.0; + Real weighted_en = 0.0; + + Real ionization_sum = 0.0; + int n = 0; // degeneracy counter + + int base_idx = get_energy_index(atomic_number); + const Real me_c2 = 510998.0; // electron rest energy [eV] + + // Convert incident energy to m_e c² units + const Real ep = e / me_c2; + + // Sum over all bound electrons (k = 0 … n_bound-1) + for (int k = 0; k < atomic_number - zstar; k++) { + + int idx_k = base_idx + k; + int idx_k_p = base_idx + (k + 1); + + // Add current electron to the running group + n++; + ionization_sum += p_ionization_energies[idx_k]; + + Real bp = p_binding_energies[idx_k] / me_c2; // binding energy in m_e c² units + + // Accumulate degenerate sub-shells (same binding energy) + if (k < atomic_number - zstar - 1 && + p_binding_energies[idx_k] == p_binding_energies[idx_k_p]) { + continue; + } + + Real e_red = ep / bp; // reduced energy T / B + + if (e_red > 1.0) { + // Relativistic β² factors + Real betae2 = 1.0 - 1.0 / ((1.0 + ep) * (1.0 + ep)); + Real betab2 = 1.0 - 1.0 / ((1.0 + bp) * (1.0 + bp)); + Real betau2 = betab2; // RBEB: U ≈ B + + // Pre-factor s0 = (r_e / 2c) * n / [B (βe² + βb² + βu²)] + Real s0 = normalization * n / + (2.0 * bp * (betae2 + betab2 + betau2)); + + Real ep2 = 1.0 / (1.0 + 0.5 * ep); + ep2 *= ep2; + + // RBEB kernel terms + Real a1c = (1.0 + 2.0 * ep) / (1.0 + e_red) * ep2; + Real a2c = (e_red - 1.0) * bp * bp * 0.5 * ep2; + Real a3 = std::log(betae2 / (1.0 - betae2)) + - betae2 + - std::log(2.0 * bp); + + Real sk = s0 * (0.5 * a3 * (1.0 - 1.0 / (e_red * e_red)) + + 1.0 - 1.0 / e_red + + a2c + - a1c * std::log(e_red)); + + sigma += sk; + + //is this correct? should degenerate electrons share one averaged ionization energy? + Real ionization_E = ionization_sum / static_cast(n); + weighted_en += sk * ionization_E; + } + + n = 0; // reset degeneracy for the next shell + ionization_sum = 0.0; + } + Real ionization_en = (sigma > 0.0) ? weighted_en / sigma : 0.0; + + return {sigma, ionization_en}; + } +} + +#endif // HIPACE_IMPACT_IONIZATION_SIGMA_H_ diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 23f085edb3..36a1354793 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -65,13 +65,14 @@ PlasmaParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_mass != 0, "The plasma particle mass must be specified"); bool ion_lev_specified = queryWithParser(pp, "initial_ion_level", m_init_ion_lev); - m_can_field_ionize = pp.contains("ionization_product"); + // m_can_field_ionize = pp.contains("ionization_product"); - queryWithParser(pp, "can_ionize", m_can_field_ionize); + // queryWithParser(pp, "can_ionize", m_can_field_ionize); + queryWithParser(pp, "can_ionize", m_can_ionize); m_can_laser_ionize = false; queryWithParser(pp, "can_laser_ionize", m_can_laser_ionize); - m_can_ionize = m_can_field_ionize || m_can_laser_ionize; + // m_can_ionize = m_can_field_ionize || m_can_laser_ionize; // !!! if(m_can_ionize) { m_neutralize_background = false; // change default diff --git a/src/utils/IonizationEnergiesTable.H b/src/utils/IonizationEnergiesTable.H index 71861a5c08..dfa45eca19 100644 --- a/src/utils/IonizationEnergiesTable.H +++ b/src/utils/IonizationEnergiesTable.H @@ -21,7 +21,7 @@ // // The Data written below is a reformatting of the data referenced form NIST. -std::map ion_map_ids = { +inline std::map ion_map_ids = { {"H", 0}, {"He", 1}, {"Li", 2}, From ac611f071836d43f96b388db9e4b5a6de04bf313 Mon Sep 17 00:00:00 2001 From: Emanuele Marchetti Date: Thu, 9 Apr 2026 14:31:51 +0200 Subject: [PATCH 5/8] WIP: bug search --- src/particles/collisions/Collision.cpp | 36 +++++++++++++++++++-- src/particles/collisions/ImpactIonization.H | 27 +++++++++------- 2 files changed, 49 insertions(+), 14 deletions(-) diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index 206c418449..cf1e5d0967 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -151,9 +151,15 @@ Collision::doCollisionA ( const amrex::Real w2 = ptd2.rdata(PlasmaIdx::w)[i2]; int ion_lev2 = ptd2.idata(PlasmaIdx::ion_lev)[i2]; - if (ion_lev2 > 0) { - // Already ionized, do not collide - return false; + // Already ionized, do not collide + if (ion_atomic_number == 1) { // H + if (ion_lev2 > 0) { + return false; + } + } else if (ion_atomic_number == 18) { // Ar + if (ion_lev2 > 3) { + return false; + } } // particle's Lorentz factor @@ -372,6 +378,9 @@ Collision::doCollisionImp ( PlasmaBins bins1 = findParticlesInEachTile(geom.Domain(), 1, species1, geom); PlasmaBins bins2 = findParticlesInEachTile(geom.Domain(), 1, species2, geom); + amrex::Print() << "Plasma bins 1 = " << bins1.numBins() << "\n"; + amrex::Print() << "Plasma bins 2 = " << bins2.numBins() << "\n"; + // offset: start/end positions of particles for each cell // perm: permutation array mapping cell entries to particle indices auto offset1 = bins1.offsetsPtr(); @@ -454,6 +463,7 @@ Collision::doCollisionImp ( } ); + amrex::Print() << "before first collision kernel\n" << std::flush; // DEBUG // loop over independent pairs amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ @@ -478,6 +488,19 @@ Collision::doCollisionImp ( const int j1 = perm1[offset1_start + idx1]; const int j2 = perm2[offset2_start + idx2]; + // DEBUG + if (j1 < 0 || j1 >= np1 || j2 < 0 || j2 >= np2) { + AMREX_DEVICE_PRINTF( + "OOB: tile np1=%d np2=%d j1=%d j2=%d icell=%d " + "offset1=[%d,%d) offset2=[%d,%d) idx1=%d idx2=%d\n", + np1, np2, j1, j2, icell, + offset1_start, offset1_stop, + offset2_start, offset2_stop, + idx1, idx2 + ); + return; + } + const bool make_new_particle = collision_function( ptd1, j1, ptd2, j2, N1, N2, icoll, @@ -499,6 +522,10 @@ Collision::doCollisionImp ( } ); + amrex::Print() << "before sync 1\n" << std::flush; // DEBUG + amrex::Gpu::streamSynchronize(); + amrex::Print() << "after sync 1\n" << std::flush; + if (!m_has_collision_product) { continue; } @@ -533,6 +560,7 @@ Collision::doCollisionImp ( ptd2 = ptile2.getParticleTileData(); auto ptd3 = ptile3.getParticleTileData(); + amrex::Print() << "before second collision kernel\n" << std::flush; // DEBUG 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); @@ -571,6 +599,8 @@ Collision::doCollisionImp ( } ); + amrex::Print() << "before sync 2\n" << std::flush; // DEBUG amrex::Gpu::streamSynchronize(); + amrex::Print() << "after sync 2\n" << std::flush; } } diff --git a/src/particles/collisions/ImpactIonization.H b/src/particles/collisions/ImpactIonization.H index 91c3d6c55a..484caffc3f 100644 --- a/src/particles/collisions/ImpactIonization.H +++ b/src/particles/collisions/ImpactIonization.H @@ -71,7 +71,7 @@ void ImpactIonization ( }; // Tolerances for numerical checks constexpr T_R rel_tol = T_R(1.0e-10); - constexpr T_R eps_n = T_R(1.0e-14); + constexpr T_R eps_n = T_R(1.0e-10); constexpr T_R beta_tol = T_R(1.0e-12); T_R const Eion = Eion_eV * PhysConstSI::q_e; // [J] @@ -222,18 +222,23 @@ void ImpactIonization ( T_R const n_1 = std::sqrt(n1sq); T_R const n_3 = T_R(1.0) - n_1 - n_2; - if (n_3 < T_R(0.0)) { - AMREX_DEVICE_PRINTF( - "ImpactIonization: invalid triangle, n1=%e n2=%e n3=%e\n", - n_1, n_2, n_3 - ); + + T_R const den_alpha = T_R(2.0)*n_1*n_3; + T_R const den_gamma = T_R(2.0)*n_1*n_2; + + if (n_1 <= eps_n || n_2 <= eps_n || n_3 <= eps_n || den_alpha <= eps_n || den_gamma <= eps_n) { + AMREX_DEVICE_PRINTF("ImpactIonization: degenerate triangle n1=%e n2=%e n3=%e\n", + n_1, n_2, n_3); return; } + T_R const C_den = n_1*n_1 / (T_R(2.0)*m1) + n_2*n_2 / (T_R(2.0)*m2) + n_3*n_3 / (T_R(2.0)*m3); + if (C_den <= eps_n) { + AMREX_DEVICE_PRINTF("ImpactIonization: bad energy_den=%e\n", C_den); + return; + } // Common momentum scale from energy conservation (non-relativistic!) - T_R const C = std::sqrt(Eout / ( - n_1*n_1 / (T_R(2.0)*m1) + n_2*n_2 / (T_R(2.0)*m2) + n_3*n_3 / (T_R(2.0)*m3)) - ); + T_R const C = std::sqrt(Eout / C_den); // Momentum magnitudes in COM T_R const p1f_mag = C * n_1; @@ -241,10 +246,10 @@ void ImpactIonization ( T_R const p3f_mag = C * n_3; // Triangle angles from the law of cosines - T_R cos_alpha = (n_1*n_1 + n_3*n_3 - n_2*n_2) / (T_R(2.0)*n_1*n_3); + T_R cos_alpha = (n_1*n_1 + n_3*n_3 - n_2*n_2) / den_alpha; cos_alpha = clamp(cos_alpha, T_R(-1.0), T_R(1.0)); - T_R cos_gamma = (n_1*n_1 + n_2*n_2 - n_3*n_3) / (T_R(2.0)*n_1*n_2); + T_R cos_gamma = (n_1*n_1 + n_2*n_2 - n_3*n_3) / den_gamma; cos_gamma = clamp(cos_gamma, T_R(-1.0), T_R(1.0)); T_R const sin_alpha = std::sqrt(amrex::max(T_R(0.0), T_R(1.0) - cos_alpha*cos_alpha)); From a022d5c7fe9ceeed4f10b1d92387c89b1bdceeb0 Mon Sep 17 00:00:00 2001 From: Emanuele Marchetti Date: Wed, 15 Apr 2026 15:13:30 +0200 Subject: [PATCH 6/8] Revised electron impact ionization: relativistic --- src/particles/collisions/Collision.cpp | 19 +-- src/particles/collisions/ImpactIonization.H | 121 ++++++++++++++------ 2 files changed, 97 insertions(+), 43 deletions(-) diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index cf1e5d0967..fe5e14bd7e 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -62,6 +62,7 @@ Collision::ReadParameters( m_collision_type == "ion_impact", "Unknown collision type" ); + // The projectile can be either an electron or an ion // Target particle is always a neutral/ion // The user is asked to always input the target type as the second specie @@ -86,7 +87,7 @@ Collision::ReadParameters( m_ionization_energies.copyToDeviceAsync(); } else { - + } // Plasma physical element name amrex::ParmParse pp_s2(m_inout_species2_name); @@ -130,6 +131,12 @@ Collision::doCollisionA ( const int ion_element_id = ion_map_ids[m_physical_element]; const int ion_atomic_number = ion_atomic_numbers[ion_element_id]; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ion_atomic_number == 1 || + ion_atomic_number == 18, + "The current implementation of electron-impact ionization only supports Hydrogen and Argon. Please check the input file and the physical element specified for species2." + ); + const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2); const amrex::Real dt = geom.CellSize(2) * inv_c; @@ -378,8 +385,8 @@ Collision::doCollisionImp ( PlasmaBins bins1 = findParticlesInEachTile(geom.Domain(), 1, species1, geom); PlasmaBins bins2 = findParticlesInEachTile(geom.Domain(), 1, species2, geom); - amrex::Print() << "Plasma bins 1 = " << bins1.numBins() << "\n"; - amrex::Print() << "Plasma bins 2 = " << bins2.numBins() << "\n"; + // amrex::Print() << "Plasma bins 1 = " << bins1.numBins() << "\n"; + // amrex::Print() << "Plasma bins 2 = " << bins2.numBins() << "\n"; // offset: start/end positions of particles for each cell // perm: permutation array mapping cell entries to particle indices @@ -463,7 +470,6 @@ Collision::doCollisionImp ( } ); - amrex::Print() << "before first collision kernel\n" << std::flush; // DEBUG // loop over independent pairs amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ @@ -522,9 +528,7 @@ Collision::doCollisionImp ( } ); - amrex::Print() << "before sync 1\n" << std::flush; // DEBUG amrex::Gpu::streamSynchronize(); - amrex::Print() << "after sync 1\n" << std::flush; if (!m_has_collision_product) { continue; @@ -560,7 +564,6 @@ Collision::doCollisionImp ( ptd2 = ptile2.getParticleTileData(); auto ptd3 = ptile3.getParticleTileData(); - amrex::Print() << "before second collision kernel\n" << std::flush; // DEBUG 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); @@ -599,8 +602,6 @@ Collision::doCollisionImp ( } ); - amrex::Print() << "before sync 2\n" << std::flush; // DEBUG amrex::Gpu::streamSynchronize(); - amrex::Print() << "after sync 2\n" << std::flush; } } diff --git a/src/particles/collisions/ImpactIonization.H b/src/particles/collisions/ImpactIonization.H index 484caffc3f..bc4e393c44 100644 --- a/src/particles/collisions/ImpactIonization.H +++ b/src/particles/collisions/ImpactIonization.H @@ -20,8 +20,8 @@ #include // isnan() isinf() #include // numeric_limits::min() -/* \brief The implicit, non-relativistic assumption has not been - * changed +/* \brief The implicit, non-relativistic assumption has been modified, + * and the ionization model is now relativistic. * Momentum and energy division after the ionization event * is done as follows: * Three numbers are generated that satisfy the triangle inequality. @@ -44,8 +44,8 @@ * n_3 = 1 - n_2 - n_1. * Next, we need to find the number C, such that if * p_i = C * n_i, we would have E_1 + E_2 + E_3 = E_out - * where E_i = p_i^2 / (2 m_i), i.e., - * C^2 * \sum_i [n_i^2 / (2 m_i)] = E_out + * where E_i = sqrt((p_i * c)^2 + (m_i * c^2)^2) - m_i * c^2. + * The root for C is found using the bisection method. * After C is determined the momentum vectors are arranged * such that the net linear momentum is zero. * To see if there are nan or inf updated velocities, @@ -139,23 +139,23 @@ void ImpactIonization ( T_R const boost_v2 = gc*gc / ((gc + T_R(1.0)) * c2); // ------------------------------------------------------------------ - // COM-frame incoming particle 1 momentum + // COM-frame boost of the initial momenta // ------------------------------------------------------------------ T_R const vcDv1 = (vcx*ux1 + vcy*uy1 + vcz*uz1) * PhysConstSI::c / g1; T_R const vcDv2 = (vcx*ux2 + vcy*uy2 + vcz*uz2) * PhysConstSI::c / g2; - // Compute p1 star + // Boost particle 1 to COM T_R p1sx; T_R p1sy; T_R p1sz; - if ( vcms > std::numeric_limits::min() ) + if ( vcms > T_R(0.0) ) { - T_R const lorentz_transform_factor = + T_R const lorentz_coeff_1 = ( boost_v2 * vcDv1 - gc ) * m1 * g1; - p1sx = px1 + vcx*lorentz_transform_factor; - p1sy = py1 + vcy*lorentz_transform_factor; - p1sz = pz1 + vcz*lorentz_transform_factor; + p1sx = px1 + vcx*lorentz_coeff_1; + p1sy = py1 + vcy*lorentz_coeff_1; + p1sz = pz1 + vcz*lorentz_coeff_1; } else // If vcms = 0, don't do Lorentz-transform. { @@ -164,24 +164,37 @@ void ImpactIonization ( p1sz = pz1; } - // ------------------------------------------------------------------ - // COM-frame gammas - // ------------------------------------------------------------------ - T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2 )*gc*g1; - T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2 )*gc*g2; - // Redundant? - if (g1s < T_R(1.0) || g2s < T_R(1.0)) { - AMREX_DEVICE_PRINTF( - "ImpactIonization: invalid COM gammas, g1s=%e g2s=%e\n", - g1s, g2s - ); - return; + // Boost particle 2 to COM + T_R p2sx; + T_R p2sy; + T_R p2sz; + if (vcms > T_R(0.0)) { + T_R const lorentz_coeff_2 = + ( boost_v2 * vcDv2 - gc ) * m2 * g2; + + p2sx = px2 + vcx*lorentz_coeff_2; + p2sy = py2 + vcy*lorentz_coeff_2; + p2sz = pz2 + vcz*lorentz_coeff_2; + } else { + p2sx = px2; + p2sy = py2; + p2sz = pz2; } + // ------------------------------------------------------------------ // Ionization probability and final-state construction + // ------------------------------------------------------------------ + T_R const p1s_sq = p1sx*p1sx + p1sy*p1sy + p1sz*p1sz; + T_R const p2s_sq = p2sx*p2sx + p2sy*p2sy + p2sz*p2sz; + + T_R const E1s = std::sqrt(p1s_sq*c2 + m1*m1*c2*c2); + T_R const E2s = std::sqrt(p2s_sq*c2 + m2*m2*c2*c2); + + T_R const K1s = E1s - m1*c2; + T_R const K2s = E2s - m2*c2; + // Available initial kinetic energy in the COM frame [J] - T_R const Ecm = - ( (g1s - T_R(1.0))*m1 + (g2s - T_R(1.0))*m2 ) * c2; + T_R const Ecm = K1s + K2s; // Kinetic energy available to the 3 outgoing particles in COM [J] T_R const Eout = Ecm - Eion; @@ -194,9 +207,8 @@ void ImpactIonization ( // // Momentum magnitudes satisfy triangle closure, // scaled so that total kinetic energy = Eout. - - T_R const E_target = (g2s-T_R(1.0))*m2 * c2; // Sample three normalized momentum magnitudes satisfying triangle inequality + T_R const E_target = K2s; T_R n_2 = amrex::min(E_target/Ecm, T_R(0.5)*amrex::Random(engine)); n_2 = clamp(n_2, T_R(0.0), T_R(0.5) - eps_n); @@ -232,13 +244,54 @@ void ImpactIonization ( return; } - T_R const C_den = n_1*n_1 / (T_R(2.0)*m1) + n_2*n_2 / (T_R(2.0)*m2) + n_3*n_3 / (T_R(2.0)*m3); - if (C_den <= eps_n) { - AMREX_DEVICE_PRINTF("ImpactIonization: bad energy_den=%e\n", C_den); - return; + // Conversion from non-relativistic to relativistic C + + // T_R const C_den = n_1*n_1 / (T_R(2.0)*m1) + n_2*n_2 / (T_R(2.0)*m2) + n_3*n_3 / (T_R(2.0)*m3); + // if (C_den <= eps_n) { + // AMREX_DEVICE_PRINTF("ImpactIonization: bad energy_den=%e\n", C_den); + // return; + // } + // T_R const C = std::sqrt(Eout / C_den); + + auto k_rel = [=] AMREX_GPU_HOST_DEVICE (T_R Crel) -> T_R + { + T_R const p1 = Crel * n_1; + T_R const p2 = Crel * n_2; + T_R const p3 = Crel * n_3; + + T_R const K1 = std::sqrt(p1*p1*c2 + m1*m1*c2*c2) - m1*c2; + T_R const K2 = std::sqrt(p2*p2*c2 + m2*m2*c2*c2) - m2*c2; + T_R const K3 = std::sqrt(p3*p3*c2 + m3*m3*c2*c2) - m3*c2; + + return K1 + K2 + K3; + }; + + // Root finding bisection method for C such that k_rel(C) = Eout + + T_R C_lo = T_R(0.0); + T_R C_hi = m1 * PhysConstSI::c; + // Grow upper bound until enough kinetic energy is reached + while (k_rel(C_hi) < Eout) { + C_hi *= T_R(2.0); + + // optional safety ceiling + if (C_hi > T_R(1.0e6) * m1 * PhysConstSI::c) { + AMREX_DEVICE_PRINTF("ImpactIonization: failed to bracket relativistic C\n"); + return; + } } - // Common momentum scale from energy conservation (non-relativistic!) - T_R const C = std::sqrt(Eout / C_den); + + // Is bisection fast enough? + for (int iter = 0; iter < 60; ++iter) { + T_R const C_mid = T_R(0.5) * (C_lo + C_hi); + if (k_rel(C_mid) < Eout) { + C_lo = C_mid; + } else { + C_hi = C_mid; + } + } + + T_R const C = T_R(0.5) * (C_lo + C_hi); // Momentum magnitudes in COM T_R const p1f_mag = C * n_1; @@ -286,7 +339,7 @@ void ImpactIonization ( T_R p2fx; T_R p2fy; T_R p2fz; T_R p3fx; T_R p3fy; T_R p3fz; - if ( vcms > std::numeric_limits::min() ) + if ( vcms > T_R(0.0) ) { T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz; T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz; From 93732499a9da9b4c98fa7c7ac9604ef74e6debfd Mon Sep 17 00:00:00 2001 From: Emanuele Marchetti Date: Sat, 25 Apr 2026 18:10:26 +0200 Subject: [PATCH 7/8] WIP: electron impact ionization --- docs/source/run/parameters.rst | 40 ++++++++++++-- src/Hipace.H | 5 +- src/Hipace.cpp | 12 ++-- src/particles/collisions/Collision.H | 6 +- src/particles/collisions/Collision.cpp | 55 ++++--------------- src/particles/plasma/MultiPlasma.cpp | 2 +- .../plasma/PlasmaParticleContainer.H | 1 + .../plasma/PlasmaParticleContainer.cpp | 13 +++-- 8 files changed, 69 insertions(+), 65 deletions(-) diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index d0c0b3abb0..c1ef383ea1 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -504,12 +504,15 @@ When both are specified, the per-species value is used. the specific ionization energy of each state. Options are: ``electron``, ``positron``, ``H``, ``D``, ``T``, ``He``, ``Li``, ``Be``, ``B``, …. -* ``.can_ionize`` (`bool`) optional (default `0`) - Whether this plasma can ionize. Can also be set to 1 by specifying ``.ionization_product``. +* ``.can_field_ionize`` (`bool`) optional (default `0`) + Whether this plasma can be ionized by fields, as for example from a driving particle beam. -* ``.can_laser_ionize`` (`bool`) optional (default `.can_ionize`) +* ``.can_laser_ionize`` (`bool`) optional (default `0`) Whether this plasma can be ionized by a laser. +* ``.can_impact_ionize`` (`bool`) optional (default `0`) + Whether this plasma can be ionized by impact ionization between plasma particles. + * ``.initial_ion_level`` (`int`) optional (default `-1`) The initial ionization state of the plasma. `0` for neutral gasses. If set, the plasma charge gets multiplied by this number. If the plasma species is not ionizable, @@ -517,7 +520,7 @@ When both are specified, the per-species value is used. * ``.ionization_product`` (`string`) optional (default "") Name of the plasma species that contains the new electrons that are produced - when this plasma gets ionized. Only needed if this plasma is ionizable. + when this plasma gets field or laser ionized. Only needed if this plasma is ionizable. * `` or plasmas.neutralize_background`` (`bool`) optional (default `1`) Whether to add a neutralizing background of immobile particles of opposite charge. @@ -1366,9 +1369,10 @@ WARNING: this module is in development. 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. +Please note that the impact ionization module is currently under development, and SI units are always required. * ``hipace.collisions`` (list of `strings`) optional - List of names of binary Coulomb collisions. + List of names of binary elastic Coulomb collisions. Each will represent collisions between 2 species. * ``.species`` (two `strings`) optional @@ -1380,6 +1384,32 @@ As collisions depend on the physical density, in normalized units `hipace.backgr Coulomb logarithm used for this collision. If not specified, the Coulomb logarithm is determined from the temperature in each cell. +* ``hipace.impact_ionization`` (list of `strings`) optional + List of names of binary inelastic Coulomb collisions. + Each will represent collisions between 2 plasma species. + Three species must be specified: the projectile, the target (ion or neutral atom), and the product (electron). + Currently only impact ionization between electrons and another species, which might be ions or neutral atoms, is implemented. + Please note that the impact ionization module is currently under development, and SI units are always required. + +* ``.type`` (string) optional + Type of impact ionization. Currently, only `electron_impact` is implemented. + +* ``.projectile`` (string) optional + The name of the projectile species for impact ionization. + The species must be defined according to `Plasma parameters`. + +* ``.target`` (string) optional + The name of the target species for impact ionization. + The species must be defined according to `Plasma parameters`. + +* ``.new_electron`` (string) optional + The name of the product species for impact ionization. + The species must be defined according to `Plasma parameters` and should be an electron species. + In electron impacts, initial and newly created electrons can be treated as on single species + simply by using the same name for the projectile and the new electrons. + However, to treat the new electrons as a separate species, a different name can be used. + In that case, the new electrons won't be contributing to further ionizations. + Radiation reaction ^^^^^^^^^^^^^^^^^^ diff --git a/src/Hipace.H b/src/Hipace.H index 7c8d5c530f..0722f18f75 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -367,8 +367,9 @@ private: std::vector m_collision_names; /** Vector of binary collisions */ amrex::Vector< CoulombCollision > m_all_collisions; - amrex::Vector m_collision_names2; - amrex::Vector m_all_collisions2; + /** Impact ionization */ + amrex::Vector m_impact_names; + amrex::Vector m_all_impacts; void InitDiagnostics (const int step, const amrex::Real time, const bool is_last_step); void FillBeamDiagnostics (const int step, const amrex::Real time, const bool is_last_step); diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 62126dfd13..6a6a8ddfb8 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -269,10 +269,10 @@ 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]); + queryWithParser(pph, "impact_ionization", m_impact_names); + for (int i = 0; i != m_impact_names.size(); ++i) { + m_all_impacts.emplace_back(Collision()); + m_all_impacts.back().ReadParameters(m_multi_plasma.m_names, m_impact_names[i]); } // external fields applied to the grid @@ -863,8 +863,8 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st // collisions for plasmas and beams doCoulombCollision(); - for (auto& collision : m_all_collisions2) { - collision.doCollision(0, m_slice_geom[0], m_multi_plasma); + for (auto& impact : m_all_impacts) { + impact.doCollision(0, m_slice_geom[0], m_multi_plasma); } // get minimum beam uz after push diff --git a/src/particles/collisions/Collision.H b/src/particles/collisions/Collision.H index f3e5efa17e..dc99067234 100644 --- a/src/particles/collisions/Collision.H +++ b/src/particles/collisions/Collision.H @@ -13,9 +13,9 @@ */ class Collision { - std::string m_inout_species1_name = ""; - std::string m_inout_species2_name = ""; - std::string m_out_species3_name = ""; + std::string m_projectile_name = ""; + std::string m_target_name = ""; + std::string m_new_electron_name = ""; std::string m_collision_type = ""; std::string m_physical_element =""; bool m_has_collision_product = false; diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index fe5e14bd7e..a98503b5c6 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -68,11 +68,11 @@ Collision::ReadParameters( // The user is asked to always input the target type as the second specie m_has_collision_product = m_collision_type == "electron_impact" || m_collision_type == "ion_impact"; - getWithParser(pp, "species1", m_inout_species1_name); - getWithParser(pp, "species2", m_inout_species2_name); + getWithParser(pp, "projectile", m_projectile_name); + getWithParser(pp, "target", m_target_name); if (m_has_collision_product) { - getWithParser(pp, "species3", m_out_species3_name); + getWithParser(pp, "new_electron", m_new_electron_name); } if (m_collision_type == "electron_impact") { @@ -90,7 +90,7 @@ Collision::ReadParameters( } // Plasma physical element name - amrex::ParmParse pp_s2(m_inout_species2_name); + amrex::ParmParse pp_s2(m_target_name); getWithParser(pp_s2, "element", m_physical_element); } @@ -120,9 +120,9 @@ Collision::doCollisionA ( const amrex::Real* p_binding_energies = m_binding_energies.data(); const amrex::Real* p_ionization_energies = m_ionization_energies.data(); - auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name); - auto& species2 = multi_plasma.GetPlasma(m_inout_species2_name); - auto& species3 = multi_plasma.GetPlasma(m_out_species3_name); + auto& species1 = multi_plasma.GetPlasma(m_projectile_name); + auto& species2 = multi_plasma.GetPlasma(m_target_name); + auto& species3 = multi_plasma.GetPlasma(m_new_electron_name); auto m1 = species1.GetMass(); auto m2 = species2.GetMass(); @@ -134,7 +134,7 @@ Collision::doCollisionA ( AMREX_ALWAYS_ASSERT_WITH_MESSAGE( ion_atomic_number == 1 || ion_atomic_number == 18, - "The current implementation of electron-impact ionization only supports Hydrogen and Argon. Please check the input file and the physical element specified for species2." + "The current implementation of electron-impact ionization only supports Hydrogen and Argon. Please check the input file and the physical element specified for target." ); const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2); @@ -374,8 +374,8 @@ Collision::doCollisionImp ( G const& ionization_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); + auto& species1 = multi_plasma.GetPlasma(m_projectile_name); + auto& species2 = multi_plasma.GetPlasma(m_target_name); for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) { amrex::removeInvalidParticles(species1.ParticlesAt(0, pti)); @@ -442,34 +442,6 @@ Collision::doCollisionImp ( 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; - } - } - ); - // loop over independent pairs amrex::ParallelForRNG(total_ind_pairs, [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ @@ -485,9 +457,6 @@ 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) { @@ -554,12 +523,12 @@ Collision::doCollisionImp ( amrex::Scan::Type::exclusive, amrex::Scan::retSum ); - auto& species3 = multi_plasma.GetPlasma(m_out_species3_name); + auto& species3 = multi_plasma.GetPlasma(m_new_electron_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 + // get new ptd after resize in case species3 is the same as species1 or species2 ptd1 = ptile1.getParticleTileData(); ptd2 = ptile2.getParticleTileData(); auto ptd3 = ptile3.getParticleTileData(); diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index f685f10ed1..3a1214fe36 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -54,7 +54,7 @@ void MultiPlasma::InitIonization (amrex::Vector gm) { for (auto& plasma : m_all_plasmas) { - if(plasma.m_can_ionize) { + if(plasma.m_can_field_ionize || plasma.m_can_laser_ionize) { for (int i=0; i Date: Fri, 22 May 2026 11:22:17 +0200 Subject: [PATCH 8/8] WIP: Add energy check --- src/particles/collisions/Collision.cpp | 65 ++++--- src/particles/collisions/ImpactIonization.H | 8 +- .../collisions/ImpactIonizationSigma.H | 158 +++++++++++++++++- .../plasma/PlasmaParticleContainer.cpp | 11 ++ 4 files changed, 211 insertions(+), 31 deletions(-) diff --git a/src/particles/collisions/Collision.cpp b/src/particles/collisions/Collision.cpp index a98503b5c6..4d89df75d5 100644 --- a/src/particles/collisions/Collision.cpp +++ b/src/particles/collisions/Collision.cpp @@ -26,18 +26,18 @@ namespace { 15.76, 15.76, 15.76, 15.76, 15.76, 15.76 }; // Very structure specific to ImpactIonizationSigma.H - constexpr amrex::Real ionizationEnergies[28] = { - /* Z=1 (H) */ 13.598434005136, +/* constexpr amrex::Real ionizationEnergies[28] = { + // Z=1 (H) 13.598434005136, - /* Z=7 (N) : 1s, 1s, 2s, 2s, 2p, 2p, 2p */ + // Z=7 (N) : 1s, 1s, 2s, 2s, 2p, 2p, 2p 667.04609, 552.06731, 97.89013, 77.4735, 47.4453, 29.60125, 14.53413, - /* Z=2 (He) : 1s, 1s */ + // Z=2 (He) : 1s, 1s 54.41776311, 24.587387936, - /* Z=18 (Ar) : 1s,1s,2s,2s,2p x6,3s x2,3p x6 */ + // Z=18 (Ar) : 1s,1s,2s,2s,2p x6,3s x2,3p x6 4426.2227, 4120.6655, 918.374, 855.47, 755.13, 685.47, 619.0, @@ -46,6 +46,22 @@ namespace { 91.290, 74.84, 59.58, 40.735, 27.62967, 15.7596112 }; +*/ + + constexpr amrex::Real ionizationEnergies[28] = { + /* H */ 13.598434005136, + + /* N: zstar = 0..6 */ + 14.53413, 29.60125, 47.4453, 77.4735, 97.89013, 552.06731, 667.04609, + + /* He: zstar = 0..1 */ + 24.587387936, 54.41776311, + + /* Ar: zstar = 0..17 */ + 15.7596112, 27.62967, 40.735, 59.58, 74.84, 91.290, + 124.41, 143.457, 422.60, 479.76, 540.4, 619.0, + 685.47, 755.13, 855.47, 918.374, 4120.6655, 4426.2227 + }; } void @@ -222,6 +238,15 @@ Collision::doCollisionA ( const amrex::Real Krel = (grel - 1.0) * m1 * c2 / PhysConstSI::q_e; // (eV) auto cs_data = coll_ion::sigma_E(ion_atomic_number, ion_lev2, Krel, p_binding_energies, p_ionization_energies); auto sigma = cs_data.sigma; + auto Eion_eV = cs_data.ionization_en; + + // Check if the collision energy is sufficient for ionization + bool impact = coll_ion::impact_energy ( + ux1, uy1, uz1, g1, + ux2, uy2, uz2, g2, + m1, m2, Eion_eV, c2, inv_c2 + ); + if (!impact) return false; // The ionization mechanism is based on the algorithm from Perez et al., Phys.Plasmas.19.083104 (2012) // The weights are rescaled according to eq (22)-(23), @@ -234,7 +259,8 @@ Collision::doCollisionA ( if (Pion > r) { // Rejection method according to the adaptation of eq (14) in the ionization model, // from Perez et al., Phys.Plasmas.19.083104 (2012) - r = amrex::Random(engine); + + r = amrex::Random(engine); return ( w1r > r*amrex::max(w1r,w2r) ); } else { return false; @@ -251,7 +277,6 @@ Collision::doCollisionA ( amrex::Real uy1 = ptd1.rdata(PlasmaIdx::uy_half_step)[i1]; amrex::Real psi1 = ptd1.rdata(PlasmaIdx::psi_half_step)[i1]; const amrex::Real w1 = ptd1.rdata(PlasmaIdx::w)[i1]; - const int ion_lev1 = ptd1.idata(PlasmaIdx::ion_lev)[i1]; amrex::Real ux2 = ptd2.rdata(PlasmaIdx::ux_half_step)[i2]; amrex::Real uy2 = ptd2.rdata(PlasmaIdx::uy_half_step)[i2]; @@ -259,26 +284,10 @@ Collision::doCollisionA ( const amrex::Real w2 = ptd2.rdata(PlasmaIdx::w)[i2]; int ion_lev2 = ptd2.idata(PlasmaIdx::ion_lev)[i2]; - // Initialize new electron properties - ptd3.id(i3) = 4; - ptd3.cpu(i3) = 0; - ptd3.rdata(PlasmaIdx::x)[i3] = ptd2.rdata(PlasmaIdx::x)[i2]; - ptd3.rdata(PlasmaIdx::y)[i3] = ptd2.rdata(PlasmaIdx::y)[i2]; - ptd3.rdata(PlasmaIdx::x_prev)[i3] = ptd2.rdata(PlasmaIdx::x_prev)[i2]; - ptd3.rdata(PlasmaIdx::y_prev)[i3] = ptd2.rdata(PlasmaIdx::y_prev)[i2]; - ptd3.rdata(PlasmaIdx::w)[i3] = ptd2.rdata(PlasmaIdx::w)[i2]; - ptd3.idata(PlasmaIdx::ion_lev)[i3] = 0; // Quantities to be updated after collision amrex::Real ux3 = 0._rt; amrex::Real uy3 = 0._rt; - amrex::Real psi3 = 1._rt; amrex::Real uz3 = 0._rt; - ptd3.rdata(PlasmaIdx::ux)[i3] = ux3; - ptd3.rdata(PlasmaIdx::uy)[i3] = uy3; - ptd3.rdata(PlasmaIdx::psi)[i3] = psi3; - ptd3.rdata(PlasmaIdx::ux_half_step)[i3] = ux3; - ptd3.rdata(PlasmaIdx::uy_half_step)[i3] = uy3; - ptd3.rdata(PlasmaIdx::psi_half_step)[i3] = psi3; // particle's Lorentz factor amrex::Real g1 = plasma_gamma(ux1, uy1, psi1, 1._rt / psi1, /* Assumes Aabssq == 0 */ 0._rt); @@ -335,6 +344,16 @@ Collision::doCollisionA ( ptd2.rdata(PlasmaIdx::psi_half_step)[i2] = plasma_psi(ux2, uy2, uz2, /* Assumes Aabssq == 0 */ 0._rt); ptd2.idata(PlasmaIdx::ion_lev)[i2] = ion_lev2; + // Initialize new electron properties + ptd3.id(i3) = 4; + ptd3.cpu(i3) = ptd2.cpu(i2); + ptd3.rdata(PlasmaIdx::x)[i3] = ptd2.rdata(PlasmaIdx::x)[i2]; + ptd3.rdata(PlasmaIdx::y)[i3] = ptd2.rdata(PlasmaIdx::y)[i2]; + ptd3.rdata(PlasmaIdx::x_prev)[i3] = ptd2.rdata(PlasmaIdx::x_prev)[i2]; + ptd3.rdata(PlasmaIdx::y_prev)[i3] = ptd2.rdata(PlasmaIdx::y_prev)[i2]; + ptd3.rdata(PlasmaIdx::w)[i3] = ptd2.rdata(PlasmaIdx::w)[i2]; + ptd3.idata(PlasmaIdx::ion_lev)[i3] = 0; + ptd3.rdata(PlasmaIdx::ux)[i3] = ux3; ptd3.rdata(PlasmaIdx::uy)[i3] = uy3; ptd3.rdata(PlasmaIdx::psi)[i3] = plasma_psi(ux3, uy3, uz3, /* Assumes Aabssq == 0 */ 0._rt); diff --git a/src/particles/collisions/ImpactIonization.H b/src/particles/collisions/ImpactIonization.H index bc4e393c44..332f2de690 100644 --- a/src/particles/collisions/ImpactIonization.H +++ b/src/particles/collisions/ImpactIonization.H @@ -198,7 +198,13 @@ void ImpactIonization ( // Kinetic energy available to the 3 outgoing particles in COM [J] T_R const Eout = Ecm - Eion; - if (Eout <= T_R(0.0)) return; // Energy not sufficient for ionization + if (Eout <= T_R(0.0)) { // Energy not sufficient for ionization + AMREX_DEVICE_PRINTF( + "ImpactIonization: energy not sufficient for ionization, E = %e = %e (Ecm) - %e (Eion)\n", + Eout, Ecm, Eion + ); + return; + } // Energy not sufficient for ionization // 3-body final state in the COM frame // Particle 1 : scattered projectile (mass m1) diff --git a/src/particles/collisions/ImpactIonizationSigma.H b/src/particles/collisions/ImpactIonizationSigma.H index c82543da40..5a08b5eba6 100644 --- a/src/particles/collisions/ImpactIonizationSigma.H +++ b/src/particles/collisions/ImpactIonizationSigma.H @@ -58,6 +58,149 @@ namespace coll_ion inline constexpr Real a0 = 5.29177210544e-11; // m inline constexpr Real normalization = 4.0 * MathConst::pi * a0*a0 * alpha*alpha*alpha*alpha; + // Check if the collision energy is sufficient for ionization + AMREX_GPU_HOST_DEVICE inline + bool impact_energy ( + Real const ux1, Real const uy1, Real const uz1, Real const g1, + Real const ux2, Real const uy2, Real const uz2, Real const g2, + Real const m1, Real const m2, Real const Eion_eV, Real const c2, Real const inv_c2) + { + using namespace amrex::literals; + + // Tolerances for numerical checks + constexpr Real rel_tol = 1e-10_rt; + constexpr Real beta_tol = 1e-12_rt; + + Real const Eion = Eion_eV * PhysConstSI::q_e; // [J] + + // ------------------------------------------------------------------ + // Early-out if relative motion is negligible + // ------------------------------------------------------------------ + Real const diffx = ux1 - ux2; + Real const diffy = uy1 - uy2; + Real const diffz = uz1 - uz2; + Real const diffm = std::sqrt(diffx*diffx + diffy*diffy + diffz*diffz); + Real const summ = std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1) + std::sqrt(ux2*ux2 + uy2*uy2 + uz2*uz2); + // If g = u1 - u2 = 0, do not collide. + // Or if the relative difference is less than 1.0e-10. + if (summ <= 0.0_rt) return false; + if (diffm / summ < rel_tol) return false; + + // ------------------------------------------------------------------ + // Initial lab-frame momenta + // ------------------------------------------------------------------ + Real const px1 = ux1 * m1 * PhysConstSI::c; + Real const py1 = uy1 * m1 * PhysConstSI::c; + Real const pz1 = uz1 * m1 * PhysConstSI::c; + + Real const px2 = ux2 * m2 * PhysConstSI::c; + Real const py2 = uy2 * m2 * PhysConstSI::c; + Real const pz2 = uz2 * m2 * PhysConstSI::c; + + // ------------------------------------------------------------------ + // COM velocity + // ------------------------------------------------------------------ + Real const mass_g = m1 * g1 + m2 * g2; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mass_g > 0.0_rt, "Mass is zero in ImpactIonizationSigma"); + + Real const vcx = (px1+px2) / mass_g; + Real const vcy = (py1+py2) / mass_g; + Real const vcz = (pz1+pz2) / mass_g; + Real const vcms = vcx*vcx + vcy*vcy + vcz*vcz; + + Real beta2 = vcms * inv_c2; + if (beta2 < 0.0_rt) { + if (beta2 > -beta_tol) { + beta2 = 0.0_rt; + } else { + AMREX_DEVICE_PRINTF( + "ImpactIonization: beta < 0, beta2=%e vcms=%e\n", + beta2, vcms + ); + return false; + } + } + + if (beta2 >= 1.0_rt - beta_tol) { + AMREX_DEVICE_PRINTF( + "ImpactIonization: beta too close to 1, beta2=%e\n", + beta2 + ); + return false; + } + + Real const gc = 1.0_rt / std::sqrt( 1.0_rt - beta2 ); + + // Stable equivalent of (gc - 1)/vcms: + // (gamma - 1)/v^2 = gamma^2 / ((gamma + 1) c^2) + Real const boost_v2 = gc*gc / ((gc + 1.0_rt) * c2); + + // ------------------------------------------------------------------ + // COM-frame boost of the initial momenta + // ------------------------------------------------------------------ + Real const vcDv1 = (vcx*ux1 + vcy*uy1 + vcz*uz1) * PhysConstSI::c / g1; + Real const vcDv2 = (vcx*ux2 + vcy*uy2 + vcz*uz2) * PhysConstSI::c / g2; + + // Boost particle 1 to COM + Real p1sx; + Real p1sy; + Real p1sz; + if ( vcms > 0.0_rt ) + { + Real const lorentz_coeff_1 = + ( boost_v2 * vcDv1 - gc ) * m1 * g1; + + p1sx = px1 + vcx*lorentz_coeff_1; + p1sy = py1 + vcy*lorentz_coeff_1; + p1sz = pz1 + vcz*lorentz_coeff_1; + } + else // If vcms = 0, don't do Lorentz-transform. + { + p1sx = px1; + p1sy = py1; + p1sz = pz1; + } + + // Boost particle 2 to COM + Real p2sx; + Real p2sy; + Real p2sz; + if (vcms > 0.0_rt) { + Real const lorentz_coeff_2 = + ( boost_v2 * vcDv2 - gc ) * m2 * g2; + + p2sx = px2 + vcx*lorentz_coeff_2; + p2sy = py2 + vcy*lorentz_coeff_2; + p2sz = pz2 + vcz*lorentz_coeff_2; + } else { + p2sx = px2; + p2sy = py2; + p2sz = pz2; + } + + // ------------------------------------------------------------------ + // Ionization probability and final-state construction + // ------------------------------------------------------------------ + Real const p1s_sq = p1sx*p1sx + p1sy*p1sy + p1sz*p1sz; + Real const p2s_sq = p2sx*p2sx + p2sy*p2sy + p2sz*p2sz; + + Real const E1s = std::sqrt(p1s_sq*c2 + m1*m1*c2*c2); + Real const E2s = std::sqrt(p2s_sq*c2 + m2*m2*c2*c2); + + Real const K1s = E1s - m1*c2; + Real const K2s = E2s - m2*c2; + + // Available initial kinetic energy in the COM frame [J] + Real const Ecm = K1s + K2s; + + // Kinetic energy available to the 3 outgoing particles in COM [J] + Real const Eout = Ecm - Eion; + if (Eout <= 0.0_rt) // Energy not sufficient for ionization + { return false; + } else { return true; + } + } + struct IonizationResult { Real sigma; @@ -82,9 +225,9 @@ namespace coll_ion IonizationResult sigma_E(int atomic_number, int zstar, Real e, const Real* p_binding_energies, const Real* p_ionization_energies) { Real sigma = 0.0; - Real weighted_en = 0.0; + Real ionization_en = 0.0; - Real ionization_sum = 0.0; + //Real ionization_sum = 0.0; int n = 0; // degeneracy counter int base_idx = get_energy_index(atomic_number); @@ -101,7 +244,7 @@ namespace coll_ion // Add current electron to the running group n++; - ionization_sum += p_ionization_energies[idx_k]; + //ionization_sum += p_ionization_energies[idx_k]; Real bp = p_binding_energies[idx_k] / me_c2; // binding energy in m_e c² units @@ -141,14 +284,15 @@ namespace coll_ion sigma += sk; //is this correct? should degenerate electrons share one averaged ionization energy? - Real ionization_E = ionization_sum / static_cast(n); - weighted_en += sk * ionization_E; + //Real ionization_E = ionization_sum / static_cast(n); + //weighted_en += sk * ionization_E; } n = 0; // reset degeneracy for the next shell - ionization_sum = 0.0; + //ionization_sum = 0.0; } - Real ionization_en = (sigma > 0.0) ? weighted_en / sigma : 0.0; + //Real ionization_en = (sigma > 0.0) ? weighted_en / sigma : 0.0; + ionization_en = p_ionization_energies[base_idx + zstar]; return {sigma, ionization_en}; } diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 5d857af30b..41528c02fe 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -13,6 +13,7 @@ #include "utils/DeprecatedInput.H" #include "utils/GPUUtil.H" #include "utils/InsituUtil.H" +#include "utils/IonizationEnergiesTable.H" #ifdef HIPACE_USE_OPENPMD # include #endif @@ -78,6 +79,16 @@ PlasmaParticleContainer::ReadParameters () "When specifying ionization_product, can_*_ionize must be set to 1 via field or laser"); if(m_can_ionize) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(element != "", + "For ionization, the element of the plasma must be specified. Please check the input file."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ion_map_ids.count(element) != 0, + "There are no ionization energies available for this element. " + "Please update src/utils/IonizationEnergiesTable.H using write_atomic_data_cpp.py"); + // Get atomic number and ionization energies from file + const int ion_element_id = ion_map_ids[element]; + const int ion_atomic_number = ion_atomic_numbers[ion_element_id]; + m_max_ion_lev = ion_atomic_number; + m_neutralize_background = false; // change default AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_init_ion_lev >= 0, "The initial ion level must be specified");