From 42e8c562345e0b67ee8026601b96d95879cf1aae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maxence=20Th=C3=A9venet?= Date: Mon, 8 Jan 2024 17:28:22 +0100 Subject: [PATCH] Implement correction to ADK for high fields (#4505) * Implement correction to ADK for high fields * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix Const Type and Typo * add documentation --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Axel Huebl --- Docs/source/usage/parameters.rst | 5 +++++ .../Particles/ElementaryProcess/Ionization.H | 13 +++++++++++- .../ElementaryProcess/Ionization.cpp | 4 ++++ .../Particles/PhysicalParticleContainer.cpp | 21 ++++++++++++++++++- Source/Particles/WarpXParticleContainer.H | 3 +++ .../Utils/Physics/IonizationEnergiesTable.H | 7 +++++++ 6 files changed, 51 insertions(+), 2 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index d3cf9b61aba..0597b21c393 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1167,6 +1167,11 @@ Particle initialization * ``.do_field_ionization`` (`0` or `1`) optional (default `0`) Do field ionization for this species (using the ADK theory). +* ``.do_adk_correction`` (`0` or `1`) optional (default `0`) + Whether to apply the correction to the ADK theory proposed by Zhang, Lan and Lu in `Q. Zhang et al. (Phys. Rev. A 90, 043410, 2014) `__. + If so, the probability of ionization is modified using an empirical model that should be more accurate in the regime of high electric fields. + Currently, this is only implemented for Hydrogen, although Argon is also available in the same reference. + * ``.physical_element`` (`string`) Only read if `do_field_ionization = 1`. Symbol of chemical element for this species. Example: for Helium, use ``physical_element = He``. diff --git a/Source/Particles/ElementaryProcess/Ionization.H b/Source/Particles/ElementaryProcess/Ionization.H index 573491645ea..61aeb19aea1 100644 --- a/Source/Particles/ElementaryProcess/Ionization.H +++ b/Source/Particles/ElementaryProcess/Ionization.H @@ -33,9 +33,11 @@ struct IonizationFilterFunc const amrex::Real* AMREX_RESTRICT m_adk_prefactor; const amrex::Real* AMREX_RESTRICT m_adk_exp_prefactor; const amrex::Real* AMREX_RESTRICT m_adk_power; + const amrex::Real* AMREX_RESTRICT m_adk_correction_factors; int comp; int m_atomic_number; + int m_do_adk_correction = 0; GetParticlePosition m_get_position; GetExternalEBField m_get_externalEB; @@ -82,8 +84,10 @@ struct IonizationFilterFunc const amrex::Real* AMREX_RESTRICT a_adk_prefactor, const amrex::Real* AMREX_RESTRICT a_adk_exp_prefactor, const amrex::Real* AMREX_RESTRICT a_adk_power, + const amrex::Real* AMREX_RESTRICT a_adk_correction_factors, int a_comp, int a_atomic_number, + int a_do_adk_correction, int a_offset = 0) noexcept; template @@ -133,9 +137,16 @@ struct IonizationFilterFunc ); // Compute probability of ionization p - const amrex::Real w_dtau = (E <= 0._rt) ? 0._rt : 1._rt/ ga * m_adk_prefactor[ion_lev] * + amrex::Real w_dtau = (E <= 0._rt) ? 0._rt : 1._rt/ ga * m_adk_prefactor[ion_lev] * std::pow(E, m_adk_power[ion_lev]) * std::exp( m_adk_exp_prefactor[ion_lev]/E ); + // if requested, do Zhang's correction of ADK + if (m_do_adk_correction) { + const amrex::Real r = E / m_adk_correction_factors[3]; + w_dtau *= std::exp(m_adk_correction_factors[0]*r*r+m_adk_correction_factors[1]*r+ + m_adk_correction_factors[2]); + } + const amrex::Real p = 1._rt - std::exp( - w_dtau ); const amrex::Real random_draw = amrex::Random(engine); diff --git a/Source/Particles/ElementaryProcess/Ionization.cpp b/Source/Particles/ElementaryProcess/Ionization.cpp index c3681a30cad..b7b91e4d4e3 100644 --- a/Source/Particles/ElementaryProcess/Ionization.cpp +++ b/Source/Particles/ElementaryProcess/Ionization.cpp @@ -30,15 +30,19 @@ IonizationFilterFunc::IonizationFilterFunc (const WarpXParIter& a_pti, int lev, const amrex::Real* const AMREX_RESTRICT a_adk_prefactor, const amrex::Real* const AMREX_RESTRICT a_adk_exp_prefactor, const amrex::Real* const AMREX_RESTRICT a_adk_power, + const amrex::Real* const AMREX_RESTRICT a_adk_correction_factors, int a_comp, int a_atomic_number, + int a_do_adk_correction, int a_offset) noexcept: m_ionization_energies{a_ionization_energies}, m_adk_prefactor{a_adk_prefactor}, m_adk_exp_prefactor{a_adk_exp_prefactor}, m_adk_power{a_adk_power}, + m_adk_correction_factors{a_adk_correction_factors}, comp{a_comp}, m_atomic_number{a_atomic_number}, + m_do_adk_correction{a_do_adk_correction}, m_Ex_external_particle{E_external_particle[0]}, m_Ey_external_particle{E_external_particle[1]}, m_Ez_external_particle{E_external_particle[2]}, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0fc718c7355..0d9180c090e 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -3174,10 +3174,15 @@ PhysicalParticleContainer::InitIonizationModule () "overriding user value and setting charge = q_e."); charge = PhysConst::q_e; } + utils::parser::queryWithParser(pp_species_name, "do_adk_correction", do_adk_correction); + utils::parser::queryWithParser( pp_species_name, "ionization_initial_level", ionization_initial_level); pp_species_name.get("ionization_product_species", ionization_product_name); pp_species_name.get("physical_element", physical_element); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + physical_element == "H" || !do_adk_correction, + "Correction to ADK by Zhang et al., PRA 90, 043410 (2014) only works with Hydrogen"); // Add runtime integer component for ionization level AddIntComp("ionizationLevel"); // Get atomic number and ionization energies from file @@ -3212,6 +3217,18 @@ PhysicalParticleContainer::InitIonizationModule () h_ionization_energies.begin(), h_ionization_energies.end(), ionization_energies.begin()); + adk_correction_factors.resize(4); + if (do_adk_correction) { + Vector h_correction_factors(4); + constexpr int offset_corr = 0; // hard-coded: only Hydrogen + for(int i=0; i<4; i++){ + h_correction_factors[i] = table_correction_factors[i+offset_corr]; + } + Gpu::copyAsync(Gpu::hostToDevice, + h_correction_factors.begin(), h_correction_factors.end(), + adk_correction_factors.begin()); + } + Real const* AMREX_RESTRICT p_ionization_energies = ionization_energies.data(); Real * AMREX_RESTRICT p_adk_power = adk_power.data(); Real * AMREX_RESTRICT p_adk_prefactor = adk_prefactor.data(); @@ -3249,8 +3266,10 @@ PhysicalParticleContainer::getIonizationFunc (const WarpXParIter& pti, adk_prefactor.dataPtr(), adk_exp_prefactor.dataPtr(), adk_power.dataPtr(), + adk_correction_factors.dataPtr(), particle_icomps["ionizationLevel"], - ion_atomic_number}; + ion_atomic_number, + do_adk_correction}; } PlasmaInjector* PhysicalParticleContainer::GetPlasmaInjector (int i) diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index e791ac25d22..33aa71d1c7d 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -445,6 +445,7 @@ protected: int do_continuous_injection = 0; int do_field_ionization = 0; + int do_adk_correction = 0; int ionization_product; std::string ionization_product_name; int ion_atomic_number; @@ -453,6 +454,8 @@ protected: amrex::Gpu::DeviceVector adk_power; amrex::Gpu::DeviceVector adk_prefactor; amrex::Gpu::DeviceVector adk_exp_prefactor; + /** for correction in Zhang et al., PRA 90, 043410 (2014). a1, a2, a3, Ecrit. */ + amrex::Gpu::DeviceVector adk_correction_factors; std::string physical_element; int do_resampling = 0; diff --git a/Source/Utils/Physics/IonizationEnergiesTable.H b/Source/Utils/Physics/IonizationEnergiesTable.H index f956e308dbe..b7c82e88247 100644 --- a/Source/Utils/Physics/IonizationEnergiesTable.H +++ b/Source/Utils/Physics/IonizationEnergiesTable.H @@ -1963,4 +1963,11 @@ namespace utils::physics } +constexpr int correction_factors_length = 4; +// a1, a2, a3, Ecrit from Zhang et al., PRA 90, 043410 (2014). +constexpr amrex::Real table_correction_factors[correction_factors_length]{ + // H + amrex::Real(0.11714), amrex::Real(-0.90933), amrex::Real(-0.06034), amrex::Real(32125000000.) +}; + #endif // #ifndef WARPX_UTILS_PHYSICS_IONIZATION_TABLE_H_