Skip to content

Commit bda9732

Browse files
haykhhaykh
committed
compton emission tested
1 parent 8b0cc5c commit bda9732

7 files changed

Lines changed: 85 additions & 41 deletions

File tree

‎input.example.toml‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -248,6 +248,18 @@
248248
# @required
249249
photon_species = ""
250250

251+
# @inferred:
252+
# - nominal_probability
253+
# @brief: Nominal probability of the emission for a particle with `gamma * beta = 1`, charge-to-mass = `q0 / m0`
254+
# @type: float
255+
# @from: `.gamma_qed`, `.photon_weight`, `...drag.compton.gamma_rad`, `scales.omegaB0`, `algorithms.timestep.dt`
256+
# @value: `0.1 * omegaB0 * dt * (gamma_qed / gamma_rad)^2 / photon_weight`
257+
# - nominal_photon_energy
258+
# @brief: Nominal energy of the emitted photon for a particle with `gamma * beta = 1`, mass = `m0`
259+
# @type: float
260+
# @from: `.gamma_qed`
261+
# @value: `(1 / gamma_qed)^2`
262+
251263
[algorithms]
252264
# Number of current smoothing passes
253265
# @type: ushort [>= 0]

‎src/engines/reporter.cpp‎

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,7 @@ namespace ntt {
7979
fmt::formatVector(ndomains_per_dim).c_str(),
8080
ndomains);
8181
report += "\n";
82-
reporter::AddCategory(report, 4, "Fiducial parameters");
82+
reporter::AddCategory(report, 4, "Nominal parameters");
8383
reporter::AddParam(report,
8484
4,
8585
"Particles per cell [ppc0]",
@@ -116,6 +116,22 @@ namespace ntt {
116116
"Magnetization [sigma0]",
117117
"%.3e",
118118
params.template get<real_t>("scales.sigma0"));
119+
120+
reporter::AddCategory(report, 6, "Compton emission");
121+
if (params.contains("radiation.emission.compton.photon_species")) {
122+
reporter::AddParam(report,
123+
6,
124+
"Nominal probability",
125+
"%.3e",
126+
params.template get<real_t>(
127+
"radiation.emission.compton.nominal_probability"));
128+
reporter::AddParam(report,
129+
6,
130+
"Nominal photon energy",
131+
"%.3e",
132+
params.template get<real_t>(
133+
"radiation.emission.compton.nominal_photon_energy"));
134+
}
119135
return report;
120136
}
121137

‎src/framework/parameters/extra.cpp‎

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,8 @@ namespace ntt {
88
namespace params {
99

1010
void Extra::read(const std::map<std::string, bool>& extra,
11-
const toml::value& toml_data) {
11+
const toml::value& toml_data,
12+
const SimulationParams* const params) {
1213
if (extra.at("synchrotron_drag")) {
1314
synchrotron_gamma_rad = toml::find_or(toml_data,
1415
"radiation",
@@ -60,35 +61,43 @@ namespace ntt {
6061
}
6162

6263
if (extra.at("compton_emission")) {
63-
compton_gamma_rad = toml::find_or(toml_data,
64+
compton_gamma_rad = toml::find_or(toml_data,
6465
"radiation",
6566
"drag",
6667
"compton",
6768
"gamma_rad",
6869
defaults::compton::gamma_rad);
69-
compton_gamma_qed = toml::find_or(toml_data,
70+
compton_gamma_qed = toml::find_or(toml_data,
7071
"radiation",
7172
"emission",
7273
"compton",
7374
"gamma_qed",
7475
defaults::compton::gamma_qed);
75-
compton_energy_min = toml::find_or(toml_data,
76+
compton_energy_min = toml::find_or(toml_data,
7677
"radiation",
7778
"emission",
7879
"compton",
7980
"photon_energy_min",
8081
defaults::compton::energy_min);
81-
compton_photon_weight = toml::find_or(toml_data,
82+
compton_photon_weight = toml::find_or(toml_data,
8283
"radiation",
8384
"emission",
8485
"compton",
8586
"photon_weight",
8687
ONE);
87-
compton_photon_species = toml::find<spidx_t>(toml_data,
88+
compton_photon_species = toml::find<spidx_t>(toml_data,
8889
"radiation",
8990
"emission",
9091
"compton",
9192
"photon_species");
93+
compton_nominal_probability = params->template get<real_t>(
94+
"scales.omegaB0") *
95+
static_cast<real_t>(0.1) *
96+
params->template get<real_t>(
97+
"algorithms.timestep.dt") *
98+
SQR(compton_gamma_qed / compton_gamma_rad) /
99+
compton_photon_weight;
100+
compton_nominal_photon_energy = ONE / SQR(compton_gamma_qed);
92101
}
93102
}
94103

@@ -123,6 +132,10 @@ namespace ntt {
123132
compton_photon_weight);
124133
params->set("radiation.emission.compton.photon_species",
125134
compton_photon_species);
135+
params->set("radiation.emission.compton.nominal_probability",
136+
compton_nominal_probability);
137+
params->set("radiation.emission.compton.nominal_photon_energy",
138+
compton_nominal_photon_energy);
126139
}
127140
}
128141
} // namespace params

‎src/framework/parameters/extra.h‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,8 +39,12 @@ namespace ntt {
3939
real_t compton_gamma_qed;
4040
real_t compton_photon_weight;
4141
spidx_t compton_photon_species;
42+
real_t compton_nominal_probability;
43+
real_t compton_nominal_photon_energy;
4244

43-
void read(const std::map<std::string, bool>&, const toml::value&);
45+
void read(const std::map<std::string, bool>&,
46+
const toml::value&,
47+
const SimulationParams* const);
4448
void setParams(const std::map<std::string, bool>&, SimulationParams*) const;
4549
};
4650

‎src/framework/parameters/parameters.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -197,7 +197,7 @@ namespace ntt {
197197
boundaries_params.setParams(this);
198198

199199
/* [algorithms] --------------------------------------------------------- */
200-
ntt::params::Algorithms alg_params {};
200+
params::Algorithms alg_params {};
201201
std::map<std::string, bool> alg_extra_flags = {
202202
{ "gr", engine_enum == SimEngine::GRPIC },
203203
{ "gca", isPromised("algorithms.gca.e_ovr_b_max") },
@@ -206,15 +206,15 @@ namespace ntt {
206206
alg_params.setParams(alg_extra_flags, this);
207207

208208
/* extra physics ------------------------------------------------------ */
209-
ntt::params::Extra extra_params {};
209+
params::Extra extra_params {};
210210
std::map<std::string, bool> extra_extra_flags = {
211211
{ "synchrotron_drag",isPromised("radiation.drag.synchrotron.gamma_rad") },
212212
{ "compton_drag", isPromised("radiation.drag.compton.gamma_rad") },
213213
{ "synchrotron_emission",
214214
isPromised("radiation.emission.synchrotron.photon_species") },
215215
{ "compton_emission", isPromised("radiation.emission.compton.photon_species") }
216216
};
217-
extra_params.read(extra_extra_flags, toml_data);
217+
extra_params.read(extra_extra_flags, toml_data, this);
218218
extra_params.setParams(extra_extra_flags, this);
219219

220220
// @TODO: disabling stats for non-Cartesian

‎src/kernels/emission/compton.hpp‎

Lines changed: 27 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -54,21 +54,17 @@ namespace kernel {
5454
npart_t domain_idx,
5555
const SimulationParams& params,
5656
random_number_pool_t& random_pool)
57-
: emitted_photon_weight { params.get<real_t>(
57+
: emitted_photon_weight { params.template get<real_t>(
5858
"radiation.emission.compton.photon_weight") }
59-
, emitted_photon_min_energy { params.get<real_t>(
59+
, emitted_photon_min_energy { params.template get<real_t>(
6060
"radiation.emission.compton.photon_energy_min") }
6161
, nominal_probability { math::abs(species_charge / species_mass) *
62-
params.get<real_t>("scales.omegaB0") *
63-
static_cast<real_t>(0.1) * dt *
64-
SQR(params.get<real_t>(
65-
"radiation.emission.compton.gamma_qed") /
66-
params.get<real_t>(
67-
"radiation.drag.compton.gamma_rad")) /
68-
emitted_photon_weight }
69-
, nominal_photon_energy { species_mass /
70-
SQR(params.get<real_t>(
71-
"radiation.emission.compton.gamma_qed")) }
62+
params.template get<real_t>(
63+
"radiation.emission.compton.nominal_"
64+
"probability") }
65+
, nominal_photon_energy { species_mass * params.template get<real_t>(
66+
"radiation.emission.compton."
67+
"nominal_photon_energy") }
7268
, should_drag { static_cast<bool>(radiative_drag_flags &
7369
RadiativeDrag::COMPTON) }
7470
, photon_i1 { photon_species.i1 }
@@ -98,8 +94,15 @@ namespace kernel {
9894

9995
/**
10096
*
101-
* @brief Determine whether a photon is emitted and compute its energy and
102-
* the recoil on the emitting particle
97+
* @brief Determine whether a photon is emitted, the drag is applied,
98+
* and compute its energy and the recoil on the emitting particle
99+
*
100+
* @param x_Cd Position of the particle (code)
101+
* @param x_Ph Position of the particle (physical)
102+
* @param u_Ph Velocity of the particle (physical)
103+
* @param ep Interpolated electric field at the particle position (physical)
104+
* @param bp Interpolated magnetic field at the particle position (physical)
105+
* @param delta_u_Ph Output parameter for the recoil on the emitting particle (physical units)
103106
*
104107
* @note
105108
*
@@ -112,18 +115,18 @@ namespace kernel {
112115
* e_gamma = (gamma / gamma_QED)^2 * (m / m0)
113116
*
114117
* drag force [in units of m c]:
115-
* du / dt = - photon_weight * p_gamma * e_gamma * u_hat
118+
* du / dt = - photon_weight * p_gamma * e_gamma * u_hat / beta
116119
*
117120
* @returns Pair of booleans to indicate whether a particle should be emitted
118121
* and whether the emitting particle should experience a recoil (i.e. radiative drag)
119122
*
120123
*/
121-
Inline auto shouldEmit(const coord_t<M::PrtlDim>& xp_Cd,
122-
const coord_t<M::PrtlDim>& xp_Ph,
123-
const vec_t<Dim::_3D>& u_Ph,
124-
const vec_t<Dim::_3D>& ep,
125-
const vec_t<Dim::_3D>& bp,
126-
vec_t<Dim::_3D>& delta_u_Ph,
124+
Inline auto shouldEmit(const coord_t<M::PrtlDim>&,
125+
const coord_t<M::PrtlDim>&,
126+
const vec_t<Dim::_3D>& u_Ph,
127+
const vec_t<Dim::_3D>&,
128+
const vec_t<Dim::_3D>&,
129+
vec_t<Dim::_3D>& delta_u_Ph,
127130
Payload& payload) const -> Kokkos::pair<bool, bool> {
128131
const auto u_sqr = NORM_SQR(u_Ph[0], u_Ph[1], u_Ph[2]);
129132
const auto gamma_sqr = ONE + u_sqr;
@@ -132,8 +135,8 @@ namespace kernel {
132135

133136
payload.photon_energy = gamma_sqr * nominal_photon_energy;
134137

135-
const auto delta_u = -probability * payload.photon_energy /
136-
math::sqrt(u_sqr);
138+
const auto delta_u = -emitted_photon_weight * payload.photon_energy /
139+
(math::sqrt(u_sqr) * beta);
137140

138141
delta_u_Ph[0] = delta_u * u_Ph[0];
139142
delta_u_Ph[1] = delta_u * u_Ph[1];
@@ -142,11 +145,10 @@ namespace kernel {
142145
auto rand_gen = random_pool.get_state();
143146
const auto should_emit = Random<real_t>(rand_gen) < probability;
144147
random_pool.free_state(rand_gen);
145-
Kokkos::printf("probability: %e\n", probability);
146148

147149
return Kokkos::make_pair(
148150
should_emit and (payload.photon_energy >= emitted_photon_min_energy),
149-
should_drag);
151+
should_drag and should_emit);
150152
}
151153

152154
Inline void emit(const tuple_t<int, M::Dim>& xi_Cd,

‎src/kernels/particle_pusher_sr.hpp‎

Lines changed: 2 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -471,16 +471,13 @@ namespace kernel::sr {
471471
delta_u_Ph,
472472
payload);
473473

474-
const auto should_emit = emission_response.first;
475-
const auto should_drag = emission_response.second;
476-
477-
if (should_drag) {
474+
if (emission_response.second) {
478475
ux1(p) += delta_u_Ph[0];
479476
ux2(p) += delta_u_Ph[1];
480477
ux3(p) += delta_u_Ph[2];
481478
}
482479

483-
if (should_emit) {
480+
if (emission_response.first) {
484481
vec_t<Dim::_3D> direction { ZERO };
485482
const auto delta_u_Ph_mag = NORM(delta_u_Ph[0],
486483
delta_u_Ph[1],

0 commit comments

Comments
 (0)