Add precalculated cross section

Uses when the calculated cut off is below the default value.
ECP-WarpX · Jan 14, 2025 · 56d5042 · 56d5042
1 parent f2f34fc
commit 56d5042
Show file tree

Hide file tree

Showing 3 changed files with 52 additions and 15 deletions.
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
@@ -2207,6 +2207,10 @@ Details about the collision models can be found in the :ref:`theory section <mul
 * ``<collision_name>.create_photons`` (`integer`)
     Only for ``bremsstrahlung``. Whether photons will be created, defaults to 1 (true).
 
+* ``<collision_name>.koT1_cut`` (`float`)
+    Only for ``bremsstrahlung``. Minimum energy of the photons created.
+    This is relative to the electron energy, defaulting to 1.e-4.
+
 .. _running-cpp-parameters-numerics:
 
 Numerics and algorithms

diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H
@@ -207,7 +207,7 @@ public:
             }
 
             // find lo-index for koT1 cutoff (will typically be 1)
-            koT1_cut = std::max(KEcut_eV/KErel_eV, 1.0e-4_prt);
+            koT1_cut = std::max(KEcut_eV/KErel_eV, m_koT1_cut_default);
             i0_cut = 0;
             for (int i=1; i < nkoT1; i++) {
                 if (m_koT1_grid[i] > koT1_cut) {
@@ -242,19 +242,25 @@ public:
                 kdsigdk_cut *= koT1_cut/((m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt);
             }
 
-            // Integrate the distribution using trapezoidal rule
-            amrex::ParticleReal koT1_grid_im1 = koT1_cut;
-            amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut;
-            amrex::ParticleReal sigma = 0.0_prt;
-            for (int i=i0_cut+1; i < nkoT1; i++) {
-                amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1);
-                amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i];
-                amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt;
-                sigma = sigma + dsigdk*dk;
-                koT1_grid_im1 = m_koT1_grid[i];
-                kdsigdk_im1 = kdsigdk_i;
+            amrex::ParticleReal sigma_total;
+            if (KEcut_eV/KErel_eV <= m_koT1_cut_default) {
+                // Use precalculated value
+                sigma_total = w0*m_sigma_total[index] + (1.0_prt - w0)*m_sigma_total[index+1];
+            } else {
+                // Integrate the distribution using trapezoidal rule
+                amrex::ParticleReal koT1_grid_im1 = koT1_cut;
+                amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut;
+                amrex::ParticleReal sigma = 0.0_prt;
+                for (int i=i0_cut+1; i < nkoT1; i++) {
+                    amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1);
+                    amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i];
+                    amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt;
+                    sigma = sigma + dsigdk*dk;
+                    koT1_grid_im1 = m_koT1_grid[i];
+                    kdsigdk_im1 = kdsigdk_i;
+                }
+                sigma_total = sigma;
             }
-            amrex::ParticleReal const sigma_total = sigma; // total cross section [m^2]
 
             return sigma_total;
         }
@@ -443,6 +449,7 @@ public:
         bool m_create_photons = true;
         bool m_need_product_data = true;
         amrex::ParticleReal m_multiplier;
+        amrex::ParticleReal m_koT1_cut_default;
 
         // relative photon energy grid ( koT1 = Ephoton_eV/KErel_eV )
         amrex::GpuArray<amrex::ParticleReal, nkoT1> m_koT1_grid = {0., 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.99, 1.00};
@@ -452,6 +459,7 @@ public:
 
         // differential total cross section
         amrex::GpuArray<amrex::GpuArray<amrex::ParticleReal, nkoT1>, nKE> m_kdsigdk;
+        amrex::GpuArray<amrex::ParticleReal, nKE> m_sigma_total;
 
     };
 

diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp
@@ -41,6 +41,10 @@ BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, Multi
                                      "BremsstrahlungFunc: The multiplier must be greater than or equal to one");
     m_exe.m_multiplier = multiplier;
 
+    amrex::ParticleReal koT1_cut_default = 1.e-4;
+    pp_collision_name.query("koT1_cut", koT1_cut_default);
+    m_exe.m_koT1_cut_default = koT1_cut_default;
+
     // Fill in the m_kdsigdk array
     UploadCrossSection(Z);
 
@@ -57,6 +61,11 @@ BremsstrahlungFunc::UploadCrossSection (int Z)
 
     std::vector<std::vector<amrex::ParticleReal>> & kdsigdk = m_kdsigdk_map.at(Z);
 
+    // kdsigdk at the default cutoff
+    int i0_cut = 0;
+    amrex::ParticleReal const koT1_cut = m_exe.m_koT1_cut_default;
+    amrex::ParticleReal const w00 = (m_exe.m_koT1_grid[i0_cut+1] - koT1_cut)/(m_exe.m_koT1_grid[i0_cut+1] - m_exe.m_koT1_grid[i0_cut]);
+
     // Convert Seltzer and Berger energy-weighted differential cross section to units of [m^2]
     for (int iee=0; iee < Executor::nKE; iee++) {
         amrex::ParticleReal const E = m_exe.m_KEgrid_eV[iee]/m_e_eV;
@@ -65,9 +74,25 @@ BremsstrahlungFunc::UploadCrossSection (int Z)
         amrex::ParticleReal const betaSq = (E*E + 2._prt*E)/gamma/gamma;
         // 1.0e-31 converts mBarn to m**2
         amrex::ParticleReal const scale_factor = 1.0e-31_prt*Z*Z/betaSq;
-        for (int iep=0; iep < Executor::nkoT1; iep++) {
-            m_exe.m_kdsigdk[iee][iep] = kdsigdk[iee][iep]*scale_factor;
+        for (int i=0; i < Executor::nkoT1; i++) {
+            m_exe.m_kdsigdk[iee][i] = kdsigdk[iee][i]*scale_factor;
+        }
+
+        // Calculate the total cross section using the default k-cutoff
+        amrex::ParticleReal const kdsigdk_cut = w00*m_exe.m_kdsigdk[iee][i0_cut] + (1.0_prt - w00)*m_exe.m_kdsigdk[iee][i0_cut+1];
+        amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut*koT1_cut/((m_exe.m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt);
+        amrex::ParticleReal koT1_im1 = koT1_cut;
+        amrex::ParticleReal sigma = 0.0_prt;
+        for (int i=i0_cut+1; i < Executor::nkoT1; i++) {
+            amrex::ParticleReal const koT1_i = m_exe.m_koT1_grid[i];
+            amrex::ParticleReal const dk = (koT1_i - koT1_im1);
+            amrex::ParticleReal const kdsigdk_i = m_exe.m_kdsigdk[iee][i];
+            amrex::ParticleReal const dsigdk = (kdsigdk_i/koT1_i + kdsigdk_im1/koT1_im1)*0.5_prt;
+            sigma = sigma + dsigdk*dk;
+            koT1_im1 = koT1_i;
+            kdsigdk_im1 = kdsigdk_i;
         }
+        m_exe.m_sigma_total[iee] = sigma;
     }
 
 }