Implement correction to ADK for high fields (ECP-WarpX#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 <[email protected]>
RevathiJambunathan · Jan 8, 2024 · 42e8c56 · 42e8c56
1 parent e944c7c
commit 42e8c56
Show file tree

Hide file tree

Showing 6 changed files with 51 additions and 2 deletions.
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
@@ -1167,6 +1167,11 @@ Particle initialization
 * ``<species>.do_field_ionization`` (`0` or `1`) optional (default `0`)
     Do field ionization for this species (using the ADK theory).
 
+* ``<species>.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) <https://doi.org/10.1103/PhysRevA.90.043410>`__.
+    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.
+
 * ``<species>.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
@@ -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<PIdx> 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 <typename PData>
@@ -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
@@ -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
@@ -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<Real> 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
@@ -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<amrex::Real> adk_power;
     amrex::Gpu::DeviceVector<amrex::Real> adk_prefactor;
     amrex::Gpu::DeviceVector<amrex::Real> adk_exp_prefactor;
+    /** for correction in Zhang et al., PRA 90, 043410 (2014). a1, a2, a3, Ecrit. */
+    amrex::Gpu::DeviceVector<amrex::Real> 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
@@ -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_