diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index b5cc679463..77b5c252f8 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -500,12 +500,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, @@ -513,7 +516,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. @@ -1307,9 +1310,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 @@ -1321,6 +1325,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 22b0aadac3..ce14bdc9d6 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 @@ -363,6 +364,9 @@ private: std::vector m_collision_names; /** Vector of binary collisions */ amrex::Vector< CoulombCollision > m_all_collisions; + /** 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 FillFieldDiagnostics (const int current_N_level, int islice); diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 717e53d4fa..782d2872d0 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -261,6 +261,12 @@ Hipace::ReadParameters () "be specified via 'hipace.background_density_SI'"); } + 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 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); @@ -844,6 +850,10 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st // collisions for plasmas and beams doCoulombCollision(); + for (auto& impact : m_all_impacts) { + impact.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/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..dc99067234 --- /dev/null +++ b/src/particles/collisions/Collision.H @@ -0,0 +1,74 @@ +#ifndef HIPACE_COLLISION_H_ +#define HIPACE_COLLISION_H_ + +#include "particles/plasma/MultiPlasma.H" + +#include + +#include + +/** + * \brief This class handles inelastic Coulomb collisions between 2 plasma species, + * and eventually leads to an ionization event. + */ +class Collision +{ + 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; + amrex::Gpu::Buffer m_binding_energies; + amrex::Gpu::Buffer m_ionization_energies; + +public: + + /** 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] geom corresponding geometry object + * \param[in,out] multi_plasma all plasmas + **/ + void doCollision ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma); + + /** + * \param[in] lev MR level + * \param[in] geom corresponding geometry object + * \param[in,out] multi_plasma all plasmas + * \param[in] collision_function return if a pair of particles collides + * \param[in] ionizaiton_function initialize product particle + **/ + template + void doCollisionImp ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionizaiton_function); + + /** + * \param[in] lev MR level + * \param[in] geom corresponding geometry object + * \param[in,out] multi_plasma all plasmas + **/ + void doCollisionA ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma); + + /** + * \param[in] lev MR level + * \param[in] geom corresponding geometry object + * \param[in,out] multi_plasma all plasmas + **/ + void doCollisionB ( + int lev, 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..a98503b5c6 --- /dev/null +++ b/src/particles/collisions/Collision.cpp @@ -0,0 +1,576 @@ +#include "Collision.H" +#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( + 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 == "electron_impact" || + 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 + m_has_collision_product = m_collision_type == "electron_impact" || m_collision_type == "ion_impact"; + + getWithParser(pp, "projectile", m_projectile_name); + getWithParser(pp, "target", m_target_name); + + if (m_has_collision_product) { + getWithParser(pp, "new_electron", m_new_electron_name); + } + + 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_binding_energies.copyToDeviceAsync(); + m_ionization_energies.copyToDeviceAsync(); + } + else { + + } + // Plasma physical element name + amrex::ParmParse pp_s2(m_target_name); + getWithParser(pp_s2, "element", m_physical_element); +} + +void +Collision::doCollision ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma) +{ + 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) +{ + 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_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(); + 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]; + + 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 target." + ); + + 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, + int N1, int N2, int icoll, + amrex::RandomEngine const& engine) + { + // 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]; + + // 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 + 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); + + // 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); + + // 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) + { + // 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, + MultiPlasma& multi_plasma) +{ + doCollisionImp(lev, geom, multi_plasma, + [=] AMREX_GPU_DEVICE (auto ptd1, int i1, auto ptd2, int i2, + 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) + { + + }); +} + +template +void +Collision::doCollisionImp ( + int lev, const amrex::Geometry& geom, + MultiPlasma& multi_plasma, + F const& collision_function, + G const& ionization_function) +{ + // assume the two species are different for now + 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)); + amrex::removeInvalidParticles(species2.ParticlesAt(0, pti)); + } + + 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(); + 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::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); + } + ); + + // loop over tiles + 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; + 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(); + + // 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); + + 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]; + + // 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, + engine + ); + + if (N2 < N1) { + idx1 += N2; + if (make_new_particle) { + p_flag1[j1] = 1; + } + } else { + idx2 += N1; + if (make_new_particle) { + p_flag2[j2] = 1; + } + } + } + } + ); + + amrex::Gpu::streamSynchronize(); + + 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_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 species2 + ptd1 = ptile1.getParticleTileData(); + ptd2 = ptile2.getParticleTileData(); + auto ptd3 = ptile3.getParticleTileData(); + + amrex::ParallelForRNG(total_ind_pairs, + [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){ + const int icell = amrex::bisect(p_num_ind_pairs, 0, num_cells, ipair); + + const int offset1_start = offset1[icell]; + const int offset1_stop = offset1[icell+1]; + const int 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_flag1[j1] : p_flag2[j2]; + + if (new_part_idx >= 0) { + ionization_function( + ptd1, j1, ptd2, j2, ptd3, old_size3 + new_part_idx, + N1, N2, icoll, + engine + ); + } + + if (N2 < N1) { + idx1 += N2; + } else { + idx2 += N1; + } + } + } + ); + + amrex::Gpu::streamSynchronize(); + } +} diff --git a/src/particles/collisions/ImpactIonization.H b/src/particles/collisions/ImpactIonization.H new file mode 100644 index 0000000000..bc4e393c44 --- /dev/null +++ b/src/particles/collisions/ImpactIonization.H @@ -0,0 +1,419 @@ +/* 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 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. + * 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 = 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, + * 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-10); + 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 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; + + // Boost particle 1 to COM + T_R p1sx; + T_R p1sy; + T_R p1sz; + if ( vcms > T_R(0.0) ) + { + T_R 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 + 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 = K1s + K2s; + + // 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. + // 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); + + 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; + + 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; + } + + // 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; + } + } + + // 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; + 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) / 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) / 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)); + 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 > T_R(0.0) ) + { + 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/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 ion_map_ids = { +inline std::map ion_map_ids = { {"H", 0}, {"He", 1}, {"Li", 2},