WarpX
Loading...
Searching...
No Matches
BremsstrahlungFunc.H
Go to the documentation of this file.
1/* Copyright 2024 David Grote
2 *
3 * This file is part of WarpX.
4 *
5 * License: BSD-3-Clause-LBNL
6 */
7
8#ifndef WARPX_BREMSSTRAHLUNG_FUNC_H_
9#define WARPX_BREMSSTRAHLUNG_FUNC_H_
10
16#include "Utils/ParticleUtils.H"
17#include "Utils/TextMsg.H"
18
19#include <AMReX_Algorithm.H>
20#include <AMReX_DenseBins.H>
21#include <AMReX_ParmParse.H>
22#include <AMReX_Random.H>
23#include <AMReX_REAL.H>
24#include <AMReX_Vector.H>
25
34{
35 // Define shortcuts for frequently-used type names
41 using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
42
43public:
44 /*
45 * \brief Default constructor of the BremsstrahlungFunc class.
46 */
47 BremsstrahlungFunc () = default;
48
56 BremsstrahlungFunc (std::string const& collision_name, MultiParticleContainer const * mypc,
57 bool /*isSameSpecies*/);
58
59 struct Executor {
83 index_type const I1s, index_type const I1e,
84 index_type const I2s, index_type const I2e,
85 index_type const* AMREX_RESTRICT I1,
86 index_type const* AMREX_RESTRICT I2,
87 SoaData_type const& soa_1, const SoaData_type& soa_2,
88 GetParticlePosition<PIdx> /*get_position_1*/, GetParticlePosition<PIdx> /*get_position_2*/,
89 amrex::ParticleReal const n1, amrex::ParticleReal const n2,
90 amrex::ParticleReal const /*T1*/, amrex::ParticleReal const /*T2*/,
91 amrex::Real const /*global_lamdb*/,
92 amrex::ParticleReal const /*q1*/, amrex::ParticleReal const /*q2*/,
93 amrex::ParticleReal const m1, amrex::ParticleReal const m2,
94 amrex::Real const dt, amrex::Real const /*dV*/, index_type coll_idx,
95 index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
96 index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
97 amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
98 amrex::ParticleReal* AMREX_RESTRICT p_product_data,
99 amrex::RandomEngine const& engine) const
100 {
101 amrex::ParticleReal * const AMREX_RESTRICT weight1 = soa_1.m_rdata[PIdx::w];
102 amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
103 amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
104 amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
105
106 amrex::ParticleReal * const AMREX_RESTRICT weight2 = soa_2.m_rdata[PIdx::w];
107 amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
108 amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
109 amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
110
111 // Number of macroparticles of each species
112 index_type const NI1 = I1e - I1s;
113 index_type const NI2 = I2e - I2s;
114
115 // Only loop over the number of particles in the first species (the electrons)
116 index_type const max_N = NI1;
117 index_type const min_N = NI2;
118
119 index_type pair_index = cell_start_pair + coll_idx;
120
121#if (defined WARPX_DIM_RZ)
122 amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
123 amrex::ParticleReal * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta];
124#endif
125
126 index_type i1 = I1s + coll_idx;
127 index_type const i2 = I2s + coll_idx;
128 // we will start from collision number = coll_idx and then add
129 // stride (smaller set size) until we do all collisions (larger set size)
130 // If there are more atoms than electrons, this loop executes once
131 // If there are more electrons than atoms, this loops over the electrons matched to the atom
132 for (index_type k = coll_idx; k < max_N; k += min_N)
133 {
134
135#if (defined WARPX_DIM_RZ)
136 /* In RZ geometry, macroparticles can collide with other macroparticles
137 * in the same *cylindrical* cell. For this reason, collisions between macroparticles
138 * are actually not local in space. In this case, the underlying assumption is that
139 * particles within the same cylindrical cell represent a cylindrically-symmetry
140 * momentum distribution function. Therefore, here, we temporarily rotate the
141 * momentum of one of the macroparticles in agreement with this cylindrical symmetry.
142 * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
143 * there is a corresponding assert statement at initialization.) */
144 amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
145 amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
146 u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
147 u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
148#endif
149
151 u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
152 u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
153 m1, m2, weight1[ I1[i1] ], weight2[ I2[i2] ],
154 n1, n2, dt, static_cast<int>(pair_index), p_mask,
155 p_pair_reaction_weight, p_product_data,
156 engine);
157
158#if (defined WARPX_DIM_RZ)
159 amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
160 u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
161 u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
162#endif
163
164 p_pair_indices_1[pair_index] = I1[i1];
165 p_pair_indices_2[pair_index] = I2[i2];
166 i1 += min_N;
167 pair_index += min_N;
168
169 }
170 }
171
172 /* \brief Calculate the cross section given the electron energy.
173 * The differential cross section is cut off below the plasma frequency since the
174 * low energy photons would all be immediately absorbed by the surrounding plasma.
175 *
176 * \param[in] E electron energy in Joules
177 * \param[in] n1 the electron number density
178 * \param[in] m1 the mass of the electrons
179 * \param[out] index the index for the electron energy grid
180 * \param[out] i0_cut the index below the cut off for the photon energy grid
181 * \param[out] koT1_cut the k of the photon energy cut off
182 * \param[out] kdsigdk_cut the ksigdk at the cut off
183 * \param[out] w0 the grid weighting at the cut off
184 * \result the cross section
185 */
187 amrex::ParticleReal
188 CalculateCrossSection (amrex::ParticleReal KErel_eV,
189 amrex::ParticleReal wpe,
190 int & index,
191 int & i0_cut,
192 amrex::ParticleReal & koT1_cut, amrex::ParticleReal & kdsigdk_cut,
193 amrex::ParticleReal & w0) const
194 {
195 using namespace amrex::literals;
196
197 amrex::ParticleReal const KEcut_eV = PhysConst::hbar*wpe/PhysConst::q_e;
198
199 if (KErel_eV < KEcut_eV) {
200 // Ignore low energy electrons which would only produce low energy photons
201 return 0._prt;
202 }
203
204 // find lo-index for koT1 cutoff (will typically be 1)
205 koT1_cut = std::max(KEcut_eV/KErel_eV, m_koT1_cut_default);
206 i0_cut = 0;
207 for (int i=1; i < nkoT1; i++) {
208 if (m_koT1_grid[i] > koT1_cut) {
209 break;
210 } else {
211 i0_cut = i;
212 }
213 }
214 if (i0_cut == nkoT1 - 1) {
215 return 0.0_prt;
216 }
217
218 // kdsigdk is linearly weighted in electron energy
219 if (KErel_eV >= m_KEgrid_eV[nKE-1]) {
220 index = nkoT1 - 2;
221 w0 = 0.0_prt;
222 } else {
223 // find index for electron energy interpolation
224 index = amrex::bisect(m_KEgrid_eV, 0, nKE, KErel_eV);
225
226 // compute interpolation weights for k*dsigma/dk
227 w0 = (m_KEgrid_eV[index+1] - KErel_eV)/(m_KEgrid_eV[index+1] - m_KEgrid_eV[index]);
228 }
229
230 // kdsigdk at the cutoff
231 amrex::ParticleReal const w00 = (m_koT1_grid[i0_cut+1] - koT1_cut)/(m_koT1_grid[i0_cut+1] - m_koT1_grid[i0_cut]);
232 amrex::ParticleReal const w01 = 1.0_prt - w00;
233 kdsigdk_cut = ( w00*(w0*m_kdsigdk[index*nkoT1 + i0_cut ] + (1.0_prt - w0)*m_kdsigdk[(index+1)*nkoT1 + i0_cut])
234 + w01*(w0*m_kdsigdk[index*nkoT1 + i0_cut+1] + (1.0_prt - w0)*m_kdsigdk[(index+1)*nkoT1 + i0_cut+1]));
235
236 amrex::ParticleReal sigma_total;
237 if (KEcut_eV/KErel_eV <= m_koT1_cut_default) {
238 // Use precalculated value
239 sigma_total = w0*m_sigma_total[index] + (1.0_prt - w0)*m_sigma_total[index+1];
240 } else {
241 // Integrate the distribution using trapezoidal rule
242 amrex::ParticleReal koT1_grid_im1 = koT1_cut;
243 amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut;
244 amrex::ParticleReal sigma = 0.0_prt;
245 for (int i=i0_cut+1; i < nkoT1; i++) {
246 amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index*nkoT1 + i] + (1.0_prt - w0)*m_kdsigdk[(index+1)*nkoT1 + i];
247 amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1);
248 amrex::ParticleReal const k_ave = (m_koT1_grid[i] + koT1_grid_im1)*0.5_prt;
249 amrex::ParticleReal const dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/k_ave;
250 sigma = sigma + dsigdk*dk;
251 koT1_grid_im1 = m_koT1_grid[i];
252 kdsigdk_im1 = kdsigdk_i;
253 }
254 sigma_total = sigma;
255 }
256
257 return sigma_total;
258 }
259
260 /* \brief Generate the photon energy from the differential cross section
261 *
262 * \param[in] index the index for the electron energy grid
263 * \param[in] i0_cut the index below the cut off for the photon energy grid
264 * \param[in] koT1_cut the k of the photon energy cut off
265 * \param[in] kdsigdk_cut the ksigdk at the cut off
266 * \param[in] w0 the grid weighting at the cut off
267 * \param[in] sigma_total the total cross section
268 * \param[in] random_number the uniformly spaced random number
269 * \result the photon energy
270 */
272 amrex::ParticleReal
273 Photon_energy(int index, int i0_cut, amrex::ParticleReal koT1_cut, amrex::ParticleReal kdsigdk_cut,
274 amrex::ParticleReal w0, amrex::ParticleReal sigma_total, amrex::ParticleReal random_number) const
275 {
276 using namespace amrex::literals;
277 amrex::ParticleReal koT1_grid_im1 = koT1_cut;
278 amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut;
279 amrex::ParticleReal sigma = 0._prt;
280 amrex::ParticleReal kdsigdk_i = kdsigdk_cut;
281 amrex::ParticleReal dk = m_koT1_grid[i0_cut+1] - koT1_grid_im1;
282 amrex::ParticleReal k_ave = koT1_cut;
283 for (int i=i0_cut+1; i < nkoT1; i++) {
284 dk = (m_koT1_grid[i] - koT1_grid_im1);
285 k_ave = (m_koT1_grid[i] + koT1_grid_im1)*0.5_prt;
286 kdsigdk_i = w0*m_kdsigdk[index*nkoT1 + i] + (1.0_prt - w0)*m_kdsigdk[(index+1)*nkoT1 + i];
287 amrex::ParticleReal const dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/k_ave;
288 amrex::ParticleReal const next_sigma = sigma + dsigdk*dk;
289 if (random_number*sigma_total > next_sigma) {
290 sigma = next_sigma;
291 koT1_grid_im1 = m_koT1_grid[i];
292 kdsigdk_im1 = kdsigdk_i;
293 } else {
294 break;
295 }
296 }
297
298 // k will be between k_im1 and k_i
299 amrex::ParticleReal const f_im1 = kdsigdk_im1/k_ave;
300 amrex::ParticleReal const f_i = kdsigdk_i/k_ave;
301 amrex::ParticleReal const x = (std::sqrt(f_im1*f_im1 + 2._prt*(f_i - f_im1)*(random_number*sigma_total - sigma)/dk)
302 - f_im1)/(f_i - f_im1);
303 amrex::ParticleReal const result = x*dk + koT1_grid_im1;
304
305 return result;
306 }
307
308 /* \brief Do the Bremsstrahlung calculation
309 *
310 * \param[in] u1x, u1y, u1z the gamma*velocity of the first species, the electrons
311 * \param[in] u2x, u2y, u2z the gamma*velocity of the second species, the atoms
312 * \param[in] m1, m2 the mass of the electrons and ions
313 * \param[in] weight1 the particle weight of the electrons
314 * \param[in] weight2 the particle weight of the target (unused)
315 * \param[in] n1 the number density of the electrons
316 * \param[in] n2 the number density of the target
317 * \param[in] dt the time step size
318 * \param[in] pair_index the index number of the pair colliding
319 * \param[out] p_mask flags whether the pair collided
320 * \param[out] p_pair_reaction_weight the particle weight of the generated photon
321 * \param[out] p_product_data the energy of the generated photon
322 * \param[in] engine the random number generating engine
323 */
325 void BremsstrahlungEvent (amrex::ParticleReal & u1x, amrex::ParticleReal & u1y,
326 amrex::ParticleReal & u1z, amrex::ParticleReal & u2x,
327 amrex::ParticleReal & u2y, amrex::ParticleReal & u2z,
328 amrex::ParticleReal m1, amrex::ParticleReal m2,
329 amrex::ParticleReal weight1, amrex::ParticleReal weight2,
330 amrex::ParticleReal n1, amrex::ParticleReal n2,
331 amrex::Real dt, int pair_index,
333 amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
334 amrex::ParticleReal* AMREX_RESTRICT p_product_data,
335 amrex::RandomEngine const& engine) const
336 {
337 using namespace amrex::literals;
338
339 constexpr auto c2 = PhysConst::c * PhysConst::c;
340
341 // compute gamma for particles 1 (electron) and 2 (ion)
342 amrex::ParticleReal const u1sq = (u1x*u1x + u1y*u1y + u1z*u1z);
343 amrex::ParticleReal const u2sq = (u2x*u2x + u2y*u2y + u2z*u2z);
344 amrex::ParticleReal const gamma1 = std::sqrt(1.0_prt + u1sq/c2);
345 amrex::ParticleReal const gamma2 = std::sqrt(1.0_prt + u2sq/c2);
346
347 // transform electron to frame where ion is at rest
348 amrex::ParticleReal u1x_rel = u1x;
349 amrex::ParticleReal u1y_rel = u1y;
350 amrex::ParticleReal u1z_rel = u1z;
351 ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, u2x, u2y, u2z);
352
353 // compute electron KE in frame where ion is at rest
354 amrex::ParticleReal const u1sq_rel = (u1x_rel*u1x_rel + u1y_rel*u1y_rel + u1z_rel*u1z_rel);
355 amrex::ParticleReal const gamma1_rel = std::sqrt(1.0_prt + u1sq_rel/c2);
356 amrex::ParticleReal const KE_eV = m1*u1sq_rel/(1.0_prt + gamma1_rel)/PhysConst::q_e;
357
358 amrex::ParticleReal const wpe = PhysConst::q_e*std::sqrt(n1/(m1*PhysConst::epsilon_0));
359
360 int i0_cut, index;
361 amrex::ParticleReal koT1_cut, kdsigdk_cut;
362 amrex::ParticleReal w0;
363 amrex::ParticleReal const sigma_total = CalculateCrossSection(KE_eV, wpe,
364 index, i0_cut, koT1_cut, kdsigdk_cut, w0);
365
366 if (sigma_total == 0._prt) { return; }
367
368 // determine if the pair collide and then do the collision
369 amrex::ParticleReal fmulti = m_multiplier; // >= 1.0
370 amrex::ParticleReal const gamma_factor = gamma1_rel/(gamma1*gamma2);
371 amrex::ParticleReal const v1 = std::sqrt(u1sq_rel)/gamma1_rel; // magnitude electron velocity
372 amrex::ParticleReal arg = fmulti*v1*sigma_total*n2*dt*gamma_factor;
373 if (arg > 1.0_prt) {
374#ifndef AMREX_USE_GPU
375 amrex::Print() << "Notice: arg = " << arg << " in interSpeciesBremsstrahlung" << "\n";
376 amrex::Print() << " v1 = " << v1 << "\n";
377 amrex::Print() << " n2 = " << n2 << "\n";
378 amrex::Print() << " sigma_total = " << sigma_total << "\n";
379#endif
380 arg = 1.0_prt;
381 fmulti = fmulti/arg;
382 if (fmulti < 1.0_prt) { fmulti = 1.0_prt; }
383 }
384 //q12 = 1.0 - exp(-arg)
385 amrex::ParticleReal const q12 = arg;
386
387 // Get a random number
388 amrex::ParticleReal const random_number = amrex::Random(engine);
389 if (random_number < q12) {
390
391 // compute weight for photon
392 amrex::ParticleReal const w_photon = weight1/fmulti;
393
394 // get energy of photon (hbar*omega)
395 amrex::ParticleReal const random_number2 = amrex::Random(engine);
396 amrex::ParticleReal const EphotonoT1 = Photon_energy(index, i0_cut, koT1_cut, kdsigdk_cut,
397 w0, sigma_total, random_number2);
398 amrex::ParticleReal Ephoton = EphotonoT1*KE_eV*PhysConst::q_e;
399
400 // limit Ephoton to KE in weighted center-of-momentum frame
401 // This is needed for physical solution that conserves both momentum and energy
402 amrex::ParticleReal const mime = (weight2*m2)/(weight1*m1);
403 amrex::ParticleReal const u1_rel = std::sqrt(u1sq_rel);
404 amrex::ParticleReal const gamma_m_1 = (u1sq_rel/c2)/(gamma1_rel + 1._prt); // gamma - 1
405 amrex::ParticleReal const gamma_m_u = 1._prt/(gamma1_rel + u1_rel/PhysConst::c); // gamma - u
406 amrex::ParticleReal const KE_cm = mime*gamma_m_1/(gamma_m_u + mime)*m1*c2;
407 Ephoton = std::min(Ephoton, KE_cm);
408
409 // The electron and ion momentum are updated in the same place where the photon is created,
410 // in PhotonCreationFunc.
411
412 p_product_data[pair_index] = Ephoton;
413 p_pair_reaction_weight[pair_index] = w_photon;
414 p_mask[pair_index] = 1;
415
416 }
417
418 }
419
423 amrex::ParticleReal m_multiplier;
424 amrex::ParticleReal m_koT1_cut_default;
425
426 int nkoT1;
427 int nKE;
428
429 amrex::ParticleReal * m_koT1_grid;
430 amrex::ParticleReal * m_KEgrid_eV;
431 amrex::ParticleReal * m_kdsigdk;
432 amrex::ParticleReal * m_sigma_total;
433
434 };
435
436 [[nodiscard]] Executor const& executor () const { return m_exe; }
437
439
440private:
441
443
445
446 // Cross section data
447
448 void UploadCrossSection (int Z);
449
450 int nkoT1;
451 int nKE;
452
453 // relative photon energy grid ( koT1 = Ephoton_eV/KErel_eV )
454 amrex::GpuArray<amrex::ParticleReal, 17> m_koT1_grid_h = {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};
455
456 // energy grid for electrons in eV
457 amrex::GpuArray<amrex::ParticleReal, 11> m_KEgrid_eV_h = {1.0e3, 2.0e3, 5.0e3, 1.0e4, 2.0e4, 5.0e4, 1.0e5, 2.0e5, 5.0e5, 1.0e6, 2.0e6};
458
459 // differential total cross section
462
467
468 // scaled energy-weighted differential total cross section for Bremsstrahlung (kdsigdk_n + Z*kdsigdk_e)
469 // for e + H from Table I of Seltzer and Berger, ATOMIC DATA AND NUCLEAR DATA TABLES 35,345-418 (1985).
470 // The values here are actualy beta1^2/Z^2*k*dsig/dk in units of mBarns = 1e-31 m^2 where Z is the
471 // atomic number and beta1 = sqrt(1.0 - 1/gamma1^2) with gamma1 = 1 + T1/(me*c^2). These vectors are
472 // converted to k*dsig/dk in units of m^2 on initialization.
473
474 std::map<int, std::vector<std::vector<amrex::ParticleReal>>> m_kdsigdk_map = {
475
476 // atomic number Z = 1
477 {1, {
478 {7.853, 7.849, 7.833, 7.800, 7.746, 7.446, 7.040, 6.586, 6.124, 5.664, 5.230, 4.841, 4.521, 4.400, 4.372, 4.362, 4.360}, // 1.0 keV
479 {8.805, 8.817, 8.801, 8.738, 8.638, 8.059, 7.377, 6.699, 6.052, 5.431, 4.839, 4.263, 3.682, 3.437, 3.374, 3.326, 3.320}, // 2.0 keV
480 {10.32, 10.25, 10.15, 9.976, 9.753, 8.703, 7.661, 6.728, 5.899, 5.148, 4.424, 3.697, 2.897, 2.446, 2.286, 2.161, 2.104}, // 5.0 keV
481 {11.55, 11.40, 11.19, 10.89, 10.55, 9.087, 7.795, 6.711, 5.776, 4.952, 4.172, 3.387, 2.509, 1.982, 1.762, 1.548, 1.439}, // 10 keV
482 {13.07, 12.61, 12.13, 11.62, 11.11, 9.370, 7.915, 6.674, 5.617, 4.734, 3.933, 3.141, 2.239, 1.684, 1.424, 1.129, 0.967}, // 20 keV
483 {14.90, 14.17, 13.42, 12.65, 11.92, 9.620, 7.858, 6.426, 5.295, 4.370, 3.575, 2.777, 1.940, 1.401, 1.112, 0.754, 0.551}, // 50 keV
484 {16.94, 15.72, 14.55, 13.47, 12.51, 9.686, 7.689, 6.107, 4.937, 3.951, 3.179, 2.438, 1.670, 1.180, 0.903, 0.548, 0.343}, // 100 keV
485 {20.31, 18.17, 16.25, 14.71, 13.46, 9.837, 7.509, 5.808, 4.595, 3.549, 2.694, 2.022, 1.362, 0.938, 0.697, 0.380, 0.197}, // 200 keV
486 {27.90, 23.39, 19.68, 17.28, 15.68, 10.84, 7.831, 5.945, 4.555, 3.400, 2.386, 1.634, 1.045, 0.706, 0.505, 0.244, 0.094}, // 500 keV
487 {35.57, 28.83, 23.45, 20.32, 18.41, 12.21, 8.960, 6.849, 5.241, 3.927, 2.763, 1.732, 1.038, 0.688, 0.483, 0.218, 0.065}, // 1.0 MeV
488 {43.42, 34.55, 27.56, 23.66, 21.45, 14.61, 11.03, 8.605, 6.744, 5.195, 3.809, 2.472, 1.306, 0.846, 0.587, 0.250, 0.057} // 2.0 MeV
489 }},
490
491 // atomic number Z = 2
492 {2, {
493 {7.167, 7.192, 7.206, 7.201, 7.181, 7.001, 6.726, 6.409, 6.077, 5.740, 5.422, 5.150, 4.939, 4.858, 4.837, 4.825, 4.821}, // 1.0 keV
494 {8.232, 8.239, 8.229, 8.187, 8.120, 7.713, 7.200, 6.666, 6.147, 5.646, 5.173, 4.729, 4.313, 4.145, 4.099, 4.064, 4.053}, // 2.0 keV
495 {9.678, 9.640, 9.570, 9.444, 9.276, 8.439, 7.568, 6.770, 6.055, 5.397, 4.770, 4.135, 3.510, 3.180, 3.070, 2.986, 2.951}, // 5.0 keV
496 {10.81, 10.71, 10.56, 10.33, 10.06, 8.834, 7.706, 6.747, 5.914, 5.166, 4.461, 3.762, 3.009, 2.590, 2.430, 2.283, 2.213}, // 10 keV
497 {12.18, 11.80, 11.40, 10.98, 10.55, 9.062, 7.784, 6.678, 5.721, 4.897, 4.148, 3.420, 2.605, 2.131, 1.927, 1.711, 1.600}, // 20 keV
498 {13.49, 12.92, 12.33, 11.71, 11.11, 9.139, 7.592, 6.336, 5.325, 4.473, 3.707, 2.950, 2.148, 1.657, 1.413, 1.130, 0.974}, // 50 keV
499 {14.54, 13.71, 12.88, 12.05, 11.28, 8.928, 7.241, 5.917, 4.892, 4.017, 3.265, 2.544, 1.797, 1.334, 1.093, 0.792, 0.624}, // 100 keV
500 {16.12, 14.79, 13.54, 12.44, 11.49, 8.747, 6.868, 5.470, 4.418, 3.513, 2.751, 2.092, 1.433, 1.023, 0.802, 0.525, 0.365}, // 200 keV
501 {19.94, 17.44, 16.27, 13.69, 12.50, 9.013, 6.783, 5.288, 4.140, 3.184, 2.345, 1.667, 1.080, 0.745, 0.556, 0.315, 0.178}, // 500 keV
502 {24.04, 20.39, 17.38, 15.42, 14.06, 9.865, 7.470, 5.820, 4.535, 3.479, 2.549, 1.714, 1.060, 0.713, 0.517, 0.266, 0.123}, // 1.0 MeV
503 {28.11, 23.46, 19.71, 17.43, 15.98, 11.46, 8.864, 7.025, 5.582, 4.375, 3.303, 2.268, 1.320, 0.866, 0.613, 0.292, 0.108} // 2.0 MeV
504 }},
505
506 // atomic number Z = 5
507 {5, {
508 {5.670, 5.717, 5.760, 5.793, 5.820, 5.871, 5.856, 5.797, 5.713, 5.610, 5.502, 5.411, 5.350, 5.323, 5.311, 5.299, 5.293}, // 1.0 keV
509 {6.814, 6.863, 6.903, 6.925, 6.932, 6.863, 6.698, 6.484, 6.250, 6.010, 5.789, 5.596, 5.438, 5.366, 5.340, 5.318, 5.306}, // 2.0 keV
510 {8.325, 8.350, 8.360, 8.324, 8.265, 7.870, 7.382, 6.895, 6.441, 6.020, 5.635, 5.285, 4.987, 4.880, 4.847, 4.814, 4.799}, // 5.0 keV
511 {9.457, 9.432, 9.377, 9.272, 9.128, 8.384, 7.614, 6.920, 6.310, 5.759, 5.259, 4.795, 4.372, 4.201, 4.150, 4.111, 4.097}, // 10 keV
512 {10.56, 10.44, 10.29, 10.07, 9.819, 8.684, 7.651, 6.784, 6.042, 5.381, 4.785, 4.222, 3.676, 3.414, 3.332, 3.275, 3.252}, // 20 keV
513 {12.00, 11.56, 11.11, 10.67, 10.23, 8.718, 7.455, 6.396, 5.515, 4.765, 4.097, 3.463, 2.785, 2.431, 2.302, 2.185, 2.136}, // 50 keV
514 {12.69, 12.07, 11.45, 10.86, 10.30, 8.452, 7.031, 5.897, 4.996, 4.223, 3.526, 2.867, 2.182, 1.807, 1.650, 1.490, 1.410}, // 100 keV
515 {13.49, 12.57, 11.70, 10.91, 10.22, 8.039, 6.503, 5.341, 4.421, 3.625, 2.914, 2.286, 1.651, 1.287, 1.121, 0.940, 0.844}, // 200 keV
516 {15.27, 13.74, 12.38, 11.33, 10.49, 7.918, 6.188, 4.945, 3.953, 3.121, 2.392, 1.765, 1.185, 0.865, 0.706, 0.521, 0.419}, // 500 keV
517 {17.23, 15.27, 13.58, 12.36, 11.41, 8.400, 6.550, 5.235, 4.166, 3.260, 2.464, 1.749, 1.131, 0.788, 0.617, 0.409, 0.293}, // 1.0 MeV
518 {19.34, 16.89, 14.86, 13.50, 12.53, 9.504, 7.545, 6.077, 4.899, 3.908, 3.027, 2.183, 1.368, 0.924, 0.693, 0.414, 0.256} // 2.0 MeV
519 }},
520
521 // atomic number Z = 6
522 {6, {
523 {5.336, 5.384, 5.427, 5.464, 5.494, 5.570, 5.585, 5.557, 5.501, 5.425, 5.337, 5.259, 5.206, 5.185, 5.175, 5.164, 5.158}, // 1.0 keV
524 {6.498, 6.553, 6.600, 6.630, 6.648, 6.628, 6.518, 6.359, 6.174, 5.976, 5.790, 5.619, 5.470, 5.400, 5.375, 5.355, 5.346}, // 2.0 keV
525 {8.028, 8.065, 8.084, 8.070, 8.031, 7.721, 7.315, 6.897, 6.502, 6.135, 5.801, 5.503, 5.251, 5.160, 5.129, 5.096, 5.080}, // 5.0 keV
526 {9.168, 9.159, 9.123, 9.041, 8.922, 8.281, 7.593, 6.964, 6.411, 5.913, 5.467, 5.064, 4.712, 4.577, 4.536, 4.500, 4.485}, // 10 keV
527 {10.27, 10.18, 10.05, 9.865, 9.638, 8.606, 7.644, 6.832, 6.140, 5.522, 4.973, 4.465, 3.994, 3.779, 3.715, 3.671, 3.654}, // 20 keV
528 {11.72, 11.31, 10.90, 10.49, 10.08, 8.651, 7.450, 6.436, 5.587, 4.862, 4.221, 3.622, 2.997, 2.684, 2.578, 2.489, 2.455}, // 50 keV
529 {12.38, 11.81, 11.24, 10.68, 10.16, 8.400, 7.026, 5.924, 5.045, 4.293, 3.611, 2.971, 2.314, 1.967, 1.830, 1.701, 1.638}, // 100 keV
530 {13.10, 12.26, 11.45, 10.72, 10.06, 7.971, 6.481, 5.348, 4.447, 3.668, 2.970, 2.351, 1.726, 1.378, 1.225, 1.068, 0.987}, // 200 keV
531 {14.65, 13.27, 12.04, 11.06, 10.26, 7.806, 6.132, 4.921, 3.953, 3.133, 2.417, 1.797, 1.222, 0.904, 0.756, 0.586, 0.494}, // 500 keV
532 {16.39, 14.64, 13.12, 11.99, 11.09, 8.243, 6.460, 5.179, 4.135, 3.248, 2.467, 1.769, 1.153, 0.814, 0.650, 0.453, 0.345}, // 1.0 MeV
533 {18.14, 16.07, 14.31, 13.09, 12.18, 9.254, 7.367, 5.974, 4.837, 3.869, 3.000, 2.184, 1.385, 0.942, 0.720, 0.453, 0.301} // 2.0 MeV
534 }}
535
536 };
537
538};
539
540#endif // WARPX_BREMSSTRAHLUNG_FUNC_H_
#define AMREX_RESTRICT
#define AMREX_INLINE
#define AMREX_GPU_DEVICE
amrex::Vector< amrex::ParticleReal > m_sigma_total_h
Definition BremsstrahlungFunc.H:461
int nKE
Definition BremsstrahlungFunc.H:451
int nkoT1
Definition BremsstrahlungFunc.H:450
Executor m_exe
Definition BremsstrahlungFunc.H:442
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_sigma_total_d
Definition BremsstrahlungFunc.H:466
void UploadCrossSection(int Z)
Definition BremsstrahlungFunc.cpp:50
WarpXParticleContainer::ParticleTileType ParticleTileType
Definition BremsstrahlungFunc.H:37
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_kdsigdk_d
Definition BremsstrahlungFunc.H:465
WarpXParticleContainer::ParticleType ParticleType
Definition BremsstrahlungFunc.H:36
std::map< int, std::vector< std::vector< amrex::ParticleReal > > > m_kdsigdk_map
Definition BremsstrahlungFunc.H:474
BremsstrahlungFunc()=default
amrex::GpuArray< amrex::ParticleReal, 17 > m_koT1_grid_h
Definition BremsstrahlungFunc.H:454
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_koT1_grid_d
Definition BremsstrahlungFunc.H:463
WarpXParticleContainer::ParticleTileType::ParticleTileDataType SoaData_type
Definition BremsstrahlungFunc.H:41
amrex::Gpu::DeviceVector< amrex::ParticleReal > m_KEgrid_eV_d
Definition BremsstrahlungFunc.H:464
bool use_global_debye_length()
Definition BremsstrahlungFunc.H:438
Executor const & executor() const
Definition BremsstrahlungFunc.H:436
amrex::DenseBins< ParticleTileDataType > ParticleBins
Definition BremsstrahlungFunc.H:39
bool m_use_global_debye_length
Definition BremsstrahlungFunc.H:444
amrex::Vector< amrex::ParticleReal > m_kdsigdk_h
Definition BremsstrahlungFunc.H:460
ParticleTileType::ParticleTileDataType ParticleTileDataType
Definition BremsstrahlungFunc.H:38
amrex::GpuArray< amrex::ParticleReal, 11 > m_KEgrid_eV_h
Definition BremsstrahlungFunc.H:457
ParticleBins::index_type index_type
Definition BremsstrahlungFunc.H:40
Definition MultiParticleContainer.H:68
AMREX_GPU_HOST_DEVICE AMREX_INLINE void doLorentzTransformWithU(amrex::ParticleReal &ux, amrex::ParticleReal &uy, amrex::ParticleReal &uz, amrex::ParticleReal const Ux, amrex::ParticleReal const Uy, amrex::ParticleReal const Uz)
Perform a Lorentz transformation of the given velocity to a frame moving with gamma*velocity (Ux,...
Definition ParticleUtils.H:119
constexpr auto c
vacuum speed of light [m/s]
Definition constant.H:149
constexpr auto epsilon_0
vacuum permittivity: dielectric permittivity of vacuum [F/m]
Definition constant.H:152
constexpr auto hbar
reduced Planck Constant = h / tau [J*s]
Definition constant.H:170
constexpr auto q_e
elementary charge [C]
Definition constant.H:158
PODVector< T, ArenaAllocator< T > > DeviceVector
__host__ __device__ T bisect(T lo, T hi, F f, T tol=1e-12, int max_iter=100)
Real Random()
Definition BremsstrahlungFunc.H:59
amrex::ParticleReal * m_kdsigdk
Definition BremsstrahlungFunc.H:431
int nkoT1
Definition BremsstrahlungFunc.H:426
AMREX_GPU_DEVICE AMREX_INLINE amrex::ParticleReal CalculateCrossSection(amrex::ParticleReal KErel_eV, amrex::ParticleReal wpe, int &index, int &i0_cut, amrex::ParticleReal &koT1_cut, amrex::ParticleReal &kdsigdk_cut, amrex::ParticleReal &w0) const
Definition BremsstrahlungFunc.H:188
amrex::ParticleReal m_multiplier
Definition BremsstrahlungFunc.H:423
int nKE
Definition BremsstrahlungFunc.H:427
amrex::ParticleReal * m_KEgrid_eV
Definition BremsstrahlungFunc.H:430
amrex::ParticleReal * m_sigma_total
Definition BremsstrahlungFunc.H:432
AMREX_GPU_DEVICE AMREX_INLINE void operator()(index_type const I1s, index_type const I1e, index_type const I2s, index_type const I2e, index_type const *AMREX_RESTRICT I1, index_type const *AMREX_RESTRICT I2, SoaData_type const &soa_1, const SoaData_type &soa_2, GetParticlePosition< PIdx >, GetParticlePosition< PIdx >, amrex::ParticleReal const n1, amrex::ParticleReal const n2, amrex::ParticleReal const, amrex::ParticleReal const, amrex::Real const, amrex::ParticleReal const, amrex::ParticleReal const, amrex::ParticleReal const m1, amrex::ParticleReal const m2, amrex::Real const dt, amrex::Real const, index_type coll_idx, index_type const cell_start_pair, index_type *AMREX_RESTRICT p_mask, index_type *AMREX_RESTRICT p_pair_indices_1, index_type *AMREX_RESTRICT p_pair_indices_2, amrex::ParticleReal *AMREX_RESTRICT p_pair_reaction_weight, amrex::ParticleReal *AMREX_RESTRICT p_product_data, amrex::RandomEngine const &engine) const
Executor of the BremsstrahlungFunc class. Produces Bremsstrahlung photons at the cell level....
Definition BremsstrahlungFunc.H:82
bool m_need_product_data
Definition BremsstrahlungFunc.H:422
bool m_computeSpeciesTemperatures
Definition BremsstrahlungFunc.H:421
AMREX_GPU_DEVICE AMREX_INLINE amrex::ParticleReal Photon_energy(int index, int i0_cut, amrex::ParticleReal koT1_cut, amrex::ParticleReal kdsigdk_cut, amrex::ParticleReal w0, amrex::ParticleReal sigma_total, amrex::ParticleReal random_number) const
Definition BremsstrahlungFunc.H:273
amrex::ParticleReal * m_koT1_grid
Definition BremsstrahlungFunc.H:429
AMREX_GPU_DEVICE AMREX_INLINE void BremsstrahlungEvent(amrex::ParticleReal &u1x, amrex::ParticleReal &u1y, amrex::ParticleReal &u1z, amrex::ParticleReal &u2x, amrex::ParticleReal &u2y, amrex::ParticleReal &u2z, amrex::ParticleReal m1, amrex::ParticleReal m2, amrex::ParticleReal weight1, amrex::ParticleReal weight2, amrex::ParticleReal n1, amrex::ParticleReal n2, amrex::Real dt, int pair_index, index_type *AMREX_RESTRICT p_mask, amrex::ParticleReal *AMREX_RESTRICT p_pair_reaction_weight, amrex::ParticleReal *AMREX_RESTRICT p_product_data, amrex::RandomEngine const &engine) const
Definition BremsstrahlungFunc.H:325
amrex::ParticleReal m_koT1_cut_default
Definition BremsstrahlungFunc.H:424
bool m_computeSpeciesDensities
Definition BremsstrahlungFunc.H:420
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel.
Definition GetAndSetPosition.H:75
@ uz
Definition WarpXParticleContainer.H:70
@ w
Definition WarpXParticleContainer.H:70
@ uy
Definition WarpXParticleContainer.H:70
@ ux
Definition WarpXParticleContainer.H:70
ParticleTileData< StorageParticleType, NArrayReal, NArrayInt > ParticleTileDataType