diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 6762648aac..d3559c5bfd 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -279,6 +279,7 @@ public: bool m_can_field_ionize = false; /**< whether this plasma can ionize from the field */ bool m_can_laser_ionize = false; /**< whether this plasma can ionize from a laser */ int m_init_ion_lev = -1; /**< initial Ion level of each particle */ + int m_max_ion_lev = 0; /**< maximum ionization level */ std::string m_product_name = ""; /**< name of Ionization product plasma */ PlasmaParticleContainer* m_product_pc = nullptr; /**< Ionization product plasma */ /** to calculate Ionization probability with ADK formula */ diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index b7dbd3501b..ae1084b77e 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -418,6 +418,7 @@ IonizationModule (const int lev, amrex::Real* AMREX_RESTRICT adk_prefactor = m_adk_prefactor.data(); amrex::Real* AMREX_RESTRICT adk_exp_prefactor = m_adk_exp_prefactor.data(); amrex::Real* AMREX_RESTRICT adk_power = m_adk_power.data(); + const int max_ion_lev = m_max_ion_lev; long num_ions = ptile_ion.numParticles(); @@ -466,6 +467,9 @@ IonizationModule (const int lev, const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, /* Assumes Aabssq == 0 */ 0._rt); const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip]; + if (ion_lev_loc >= max_ion_lev) { + return; + } // gamma / (psi + 1) to complete dt for QSA amrex::Real w_dtau = gamma_psi * adk_prefactor[ion_lev_loc] * std::pow(Ep, adk_power[ion_lev_loc]) * @@ -611,6 +615,7 @@ LaserIonization (const int islice, amrex::Real* AMREX_RESTRICT laser_adk_prefactor = m_laser_adk_prefactor.data(); amrex::Real* AMREX_RESTRICT laser_dp_prefactor = m_laser_dp_prefactor.data(); amrex::Real* AMREX_RESTRICT laser_dp_second_prefactor = m_laser_dp_second_prefactor.data(); + const int max_ion_lev = m_max_ion_lev; long num_ions = ptile_ion.numParticles(); @@ -660,6 +665,9 @@ LaserIonization (const int islice, const amrex::Real gamma_psi = plasma_gamma_psi(ux, uy, 1._rt / psi, /* Assumes Aabssq == 0 */ 0._rt); const int ion_lev_loc = ptd_ion.idata(PlasmaIdx::ion_lev)[ip]; + if (ion_lev_loc >= max_ion_lev) { + return; + } // gamma / (psi + 1) to complete dt for QSA amrex::Real w_dtau_dc = gamma_psi * adk_prefactor[ion_lev_loc] * std::pow(Ep, adk_power[ion_lev_loc]) * diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp index 46e38a0fc4..4f6d955d0a 100644 --- a/src/particles/plasma/PlasmaParticleContainerInit.cpp +++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp @@ -426,6 +426,7 @@ InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_ // Get atomic number and ionization energies from file const int ion_element_id = ion_map_ids[physical_element]; const int ion_atomic_number = ion_atomic_numbers[ion_element_id]; + m_max_ion_lev = ion_atomic_number; amrex::Vector h_ionization_energies(ion_atomic_number); const int offset = ion_energy_offsets[ion_element_id]; for(int i=0; i