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
|
/* Copyright 2021 Modern Electron
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef WARPX_PARTICLES_COLLISION_MCCPROCESS_H_
#define WARPX_PARTICLES_COLLISION_MCCPROCESS_H_
#include <AMReX_Math.H>
#include <AMReX_Vector.H>
#include <AMReX_RandomEngine.H>
#include <AMReX_GpuContainers.H>
enum class MCCProcessType {
INVALID,
ELASTIC,
BACK,
CHARGE_EXCHANGE,
EXCITATION,
IONIZATION,
};
class MCCProcess
{
public:
MCCProcess (
const std::string& scattering_process,
const std::string& cross_section_file,
amrex::ParticleReal energy
);
template <typename InputVector>
MCCProcess (
const std::string& scattering_process,
const InputVector&& energies,
const InputVector&& sigmas,
amrex::ParticleReal energy
);
MCCProcess (MCCProcess const&) = delete;
MCCProcess& operator= (MCCProcess const&) = delete;
MCCProcess (MCCProcess &&) = default;
MCCProcess& operator= (MCCProcess &&) = default;
/** Read the given cross-section data file to memory.
*
* @param cross_section_file the path to the file containing the cross-
* section data
* @param energies vector storing energy values in eV
* @param sigmas vector storing cross-section values
*
*/
static
void readCrossSectionFile (
std::string cross_section_file,
amrex::Vector<amrex::ParticleReal>& energies,
amrex::Gpu::HostVector<amrex::ParticleReal>& sigmas
);
static
void sanityCheckEnergyGrid (
const amrex::Vector<amrex::ParticleReal>& energies,
amrex::ParticleReal dE
);
struct Executor {
/** Get the collision cross-section using a simple linear interpolator. If
* the energy value is lower (higher) than the given energy range the
* first (last) cross-section value is used.
*
* @param E_coll collision energy in eV
*
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::ParticleReal getCrossSection (amrex::ParticleReal E_coll) const
{
if (E_coll < m_energy_lo) {
return m_sigma_lo;
} else if (E_coll > m_energy_hi) {
return m_sigma_hi;
} else {
using amrex::Math::floor;
using amrex::Math::ceil;
// calculate index of bounding energy pairs
amrex::ParticleReal temp = (E_coll - m_energy_lo) / m_dE;
const int idx_1 = static_cast<int>(floor(temp));
const int idx_2 = static_cast<int>(ceil(temp));
// linearly interpolate to the given energy value
temp -= idx_1;
return m_sigmas_data[idx_1] + (m_sigmas_data[idx_2] - m_sigmas_data[idx_1]) * temp;
}
}
amrex::ParticleReal* m_sigmas_data = nullptr;
amrex::ParticleReal m_energy_lo, m_energy_hi, m_sigma_lo, m_sigma_hi, m_dE;
amrex::ParticleReal m_energy_penalty;
MCCProcessType m_type;
};
Executor const& executor () const {
#ifdef AMREX_USE_GPU
return m_exe_d;
#else
return m_exe_h;
#endif
}
amrex::ParticleReal getCrossSection (amrex::ParticleReal E_coll) const
{
return m_exe_h.getCrossSection(E_coll);
}
amrex::ParticleReal getEnergyPenalty () const { return m_exe_h.m_energy_penalty; }
amrex::ParticleReal getMinEnergyInput () const { return m_exe_h.m_energy_lo; }
amrex::ParticleReal getMaxEnergyInput () const { return m_exe_h.m_energy_hi; }
amrex::ParticleReal getEnergyInputStep () const { return m_exe_h.m_dE; }
MCCProcessType type () const { return m_exe_h.m_type; }
private:
static
MCCProcessType parseProcessType(const std::string& process);
void init (const std::string& scattering_process, amrex::ParticleReal energy);
amrex::Vector<amrex::ParticleReal> m_energies;
#ifdef AMREX_USE_GPU
amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sigmas_d;
Executor m_exe_d;
#endif
amrex::Gpu::HostVector<amrex::ParticleReal> m_sigmas_h;
Executor m_exe_h;
int m_grid_size;
};
#endif // WARPX_PARTICLES_COLLISION_MCCPROCESS_H_
|
