From 56d504264852d3c402edc7aa9951f9c4ab48601c Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 15:31:39 -0800 Subject: [PATCH] Add precalculated cross section Uses when the calculated cut off is below the default value. --- Docs/source/usage/parameters.rst | 4 +++ .../Bremsstrahlung/BremsstrahlungFunc.H | 34 ++++++++++++------- .../Bremsstrahlung/BremsstrahlungFunc.cpp | 29 ++++++++++++++-- 3 files changed, 52 insertions(+), 15 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index a91cb03c815..a5193d7f163 100644 --- 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 .create_photons`` (`integer`) Only for ``bremsstrahlung``. Whether photons will be created, defaults to 1 (true). +* ``.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 index d36eaceaebd..bd9ac7fd38b 100644 --- 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 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, nKE> m_kdsigdk; + amrex::GpuArray m_sigma_total; }; diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 6cb8ce3b93e..4eff0bcf84f 100644 --- 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> & 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; } }