1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
/* Copyright 2019-2020 Andrew Myers, Axel Huebl,
* Maxence Thevenet
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef IONIZATION_H_
#define IONIZATION_H_
#include "Particles/Gather/FieldGather.H"
#include "Particles/Gather/GetExternalFields.H"
#include "Particles/Pusher/GetAndSetPosition.H"
#include "Particles/WarpXParticleContainer.H"
#include "Utils/WarpXConst.H"
#include <AMReX_Array.H>
#include <AMReX_Array4.H>
#include <AMReX_Dim3.H>
#include <AMReX_Extension.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_IndexType.H>
#include <AMReX_REAL.H>
#include <AMReX_Random.H>
#include <AMReX_BaseFwd.H>
#include <cmath>
struct IonizationFilterFunc
{
const amrex::Real* AMREX_RESTRICT m_ionization_energies;
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;
int comp;
int m_atomic_number;
GetParticlePosition m_get_position;
GetExternalEBField m_get_externalEB;
amrex::Array4<const amrex::Real> m_ex_arr;
amrex::Array4<const amrex::Real> m_ey_arr;
amrex::Array4<const amrex::Real> m_ez_arr;
amrex::Array4<const amrex::Real> m_bx_arr;
amrex::Array4<const amrex::Real> m_by_arr;
amrex::Array4<const amrex::Real> m_bz_arr;
amrex::IndexType m_ex_type;
amrex::IndexType m_ey_type;
amrex::IndexType m_ez_type;
amrex::IndexType m_bx_type;
amrex::IndexType m_by_type;
amrex::IndexType m_bz_type;
amrex::GpuArray<amrex::Real, 3> m_dx_arr;
amrex::GpuArray<amrex::Real, 3> m_xyzmin_arr;
bool m_galerkin_interpolation;
int m_nox;
int m_n_rz_azimuthal_modes;
amrex::Dim3 m_lo;
IonizationFilterFunc (const WarpXParIter& a_pti, int lev, amrex::IntVect ngEB,
amrex::FArrayBox const& exfab,
amrex::FArrayBox const& eyfab,
amrex::FArrayBox const& ezfab,
amrex::FArrayBox const& bxfab,
amrex::FArrayBox const& byfab,
amrex::FArrayBox const& bzfab,
const amrex::Real* AMREX_RESTRICT a_ionization_energies,
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,
int a_comp,
int a_atomic_number,
int a_offset = 0) noexcept;
template <typename PData>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
bool operator() (const PData& ptd, int i, amrex::RandomEngine const& engine) const noexcept
{
using namespace amrex::literals;
const int ion_lev = ptd.m_runtime_idata[comp][i];
if (ion_lev < m_atomic_number)
{
constexpr amrex::Real c = PhysConst::c;
constexpr amrex::Real c2_inv = amrex::Real(1.)/c/c;
// gather E and B
amrex::ParticleReal xp, yp, zp;
m_get_position(i, xp, yp, zp);
amrex::ParticleReal ex = 0._rt, ey = 0._rt, ez = 0._rt;
amrex::ParticleReal bx = 0._rt, by = 0._rt, bz = 0._rt;
m_get_externalEB(i, ex, ey, ez, bx, by, bz);
doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes,
m_nox, m_galerkin_interpolation);
// Compute electric field amplitude in the particle's frame of
// reference (particularly important when in boosted frame).
const amrex::ParticleReal ux = ptd.m_rdata[PIdx::ux][i];
const amrex::ParticleReal uy = ptd.m_rdata[PIdx::uy][i];
const amrex::ParticleReal uz = ptd.m_rdata[PIdx::uz][i];
const amrex::Real ga = static_cast<amrex::Real>(
std::sqrt(1. + (ux*ux + uy*uy + uz*uz) * c2_inv));
const amrex::Real E = std::sqrt(
- ( ux*ex + uy*ey + uz*ez ) * ( ux*ex + uy*ey + uz*ez ) * c2_inv
+ ( ga *ex + uy*bz - uz*by ) * ( ga *ex + uy*bz - uz*by )
+ ( ga *ey + uz*bx - ux*bz ) * ( ga *ey + uz*bx - ux*bz )
+ ( ga *ez + ux*by - uy*bx ) * ( ga *ez + ux*by - uy*bx )
);
// Compute probability of ionization p
const 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 );
const amrex::Real p = 1._rt - std::exp( - w_dtau );
const amrex::Real random_draw = amrex::Random(engine);
if (random_draw < p)
{
return true;
}
}
return false;
}
};
struct IonizationTransformFunc
{
template <typename DstData, typename SrcData>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void operator() (DstData& /*dst*/, SrcData& src,
int i_src, int /*i_dst*/,
amrex::RandomEngine const& /*engine*/) const noexcept
{
src.m_runtime_idata[0][i_src] += 1;
}
};
#endif
|
