Skip to content

Commit 35bcce1

Browse files
JamesEdgeleyJamesEdgeley
committed
propagate normalisations to amjuel
1 parent ad8f549 commit 35bcce1

3 files changed

Lines changed: 73 additions & 52 deletions

File tree

‎src/ParticleSystems/AMJUEL.hpp‎

Lines changed: 18 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -87,24 +87,26 @@ class AMJUEL
8787
// Ionisation
8888
public:
8989
static inline const auto ionise_rate_data(
90+
double dens, double temp, double time,
9091
const std::string &filename = "data/H.4_2.1.5.csv")
9192
{
9293
const auto h4_2_1_5_coeffs =
9394
fetch_amjuel_coeffs<num_coeffs_n, num_coeffs_T>(filename);
9495

95-
return AMJUEL2DData<num_coeffs_T, num_coeffs_n>(
96-
1.0, norm::dens, norm::temp, norm::time, h4_2_1_5_coeffs);
96+
return AMJUEL2DData<num_coeffs_T, num_coeffs_n>(1.0, dens, temp, time,
97+
h4_2_1_5_coeffs);
9798
}
9899

99100
static inline const auto ionise_energy_data(
101+
double dens, double temp, double time, double vel,
100102
const std::string &filename = "data/H.10_2.1.8.csv")
101103
{
102104
const auto h10_2_1_5_coeffs =
103105
fetch_amjuel_coeffs<num_coeffs_np, num_coeffs_T>(filename);
104106

105107
return AMJUEL2DData<num_coeffs_T, num_coeffs_np>(
106-
norm::mass_amu * norm::vel * norm::vel, norm::dens, norm::temp,
107-
norm::time, h10_2_1_5_coeffs);
108+
norm::mass_amu * vel * vel, dens, temp, time,
109+
h10_2_1_5_coeffs);
108110
}
109111

110112
// Charge Exchange
@@ -132,45 +134,48 @@ class AMJUEL
132134

133135
public:
134136
static inline const auto cx_rate_data(
135-
double parent_mass, double child_mass,
137+
double parent_mass, double child_mass, double dens, double temp,
138+
double time, double vel,
136139
const std::string &filename = "data/H.3_3.1.8.csv")
137140
{
138141
const auto h3_3_1_8_coeffs =
139142
fetch_amjuel_coeffs<num_coeffs_E, num_coeffs_T>(filename);
140143
return AMJUEL2DDataH3<num_coeffs_T, num_coeffs_E, 2>(
141-
1.0, norm::dens, norm::temp / child_mass, norm::time, norm::vel,
142-
parent_mass, h3_3_1_8_coeffs);
144+
1.0, dens, temp / child_mass, time, vel, parent_mass,
145+
h3_3_1_8_coeffs);
143146
}
144147

145-
static inline const auto amjuel_fit_cross_section(double reduced_mass)
148+
static inline const auto amjuel_fit_cross_section(double reduced_mass,
149+
double vel)
146150
{
147151
static constexpr double amjuel_cs_norm = 1E-4;
148152

149153
return AMJUELFitCrossSection<h1_3_1_8_num_coeffs, h1_3_1_8_num_l_coeffs,
150154
h1_3_1_8_num_r_coeffs>(
151-
norm::vel, amjuel_cs_norm, reduced_mass, h1_3_1_8_coeffs,
155+
vel, amjuel_cs_norm, reduced_mass, h1_3_1_8_coeffs,
152156
h1_3_1_8_l_coeffs, h1_3_1_8_r_coeffs, h1_3_1_8_elabmin,
153157
h1_3_1_8_elabmax, h1_3_1_8_emax);
154158
}
155159

156160
// Recombination
157161
static inline const auto recomb_rate_data(
162+
double dens, double temp, double time,
158163
const std::string &filename = "data/H.4_2.1.8.csv")
159164
{
160165
const auto h4_2_1_8_coeffs =
161166
fetch_amjuel_coeffs<num_coeffs_T, num_coeffs_n>(filename);
162-
return AMJUEL2DData<num_coeffs_T, num_coeffs_n>(
163-
1.0, norm::dens, norm::temp, norm::time, h4_2_1_8_coeffs);
167+
return AMJUEL2DData<num_coeffs_T, num_coeffs_n>(1.0, dens, temp, time,
168+
h4_2_1_8_coeffs);
164169
}
165170

166171
static inline const auto recomb_energy_data(
172+
double dens, double temp, double time,double vel,
167173
const std::string &filename = "data/H.10_2.1.8.csv")
168174
{
169175
const auto h10_2_1_8_coeffs =
170176
fetch_amjuel_coeffs<num_coeffs_T, num_coeffs_np>(filename);
171177
return AMJUEL2DData<num_coeffs_T, num_coeffs_np>(
172-
norm::mass_amu * norm::vel * norm::vel, norm::dens, norm::temp,
173-
norm::time, h10_2_1_8_coeffs);
178+
norm::mass_amu * vel * vel, dens, temp, time, h10_2_1_8_coeffs);
174179
}
175180
};
176181
} // namespace PENKNIFE

‎src/ParticleSystems/ReactionSystem.cpp‎

Lines changed: 30 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -66,12 +66,16 @@ void ReactionSystem::set_up_reactions()
6666
if (this->vdim == 2)
6767
{
6868
reaction = ionise_reaction_amjuel<2>(
69-
this->sycl_target, target_species, electron_species);
69+
this->sycl_target, Nnorm, Tnorm, 1. / omega_c,
70+
mesh_length * omega_c, target_species,
71+
electron_species);
7072
}
7173
else if (this->vdim == 3)
7274
{
7375
reaction = ionise_reaction_amjuel<3>(
74-
this->sycl_target, target_species, electron_species);
76+
this->sycl_target, Nnorm, Tnorm, 1. / omega_c,
77+
mesh_length * omega_c, target_species,
78+
electron_species);
7579
}
7680
}
7781
}
@@ -96,28 +100,32 @@ void ReactionSystem::set_up_reactions()
96100
if (this->vdim == 2)
97101
{
98102
reaction = recombination_reaction_fixed<2>(
99-
this->sycl_target, rng_kernel, marker_species,
100-
electron_species, neutral_species, rate, energy_rate);
103+
this->sycl_target, rng_kernel, mesh_length * omega_c,
104+
marker_species, electron_species, neutral_species, rate,
105+
energy_rate);
101106
}
102107
else if (this->vdim == 3)
103108
{
104109
reaction = recombination_reaction_fixed<3>(
105-
this->sycl_target, rng_kernel, marker_species,
106-
electron_species, neutral_species, rate, energy_rate);
110+
this->sycl_target, rng_kernel, mesh_length * omega_c,
111+
marker_species, electron_species, neutral_species, rate,
112+
energy_rate);
107113
}
108114
}
109115
else if (std::get<2>(v).first == "AMJUEL")
110116
{
111117
if (this->vdim == 2)
112118
{
113119
reaction = recombination_reaction_amjuel<2>(
114-
this->sycl_target, rng_kernel, marker_species,
120+
this->sycl_target, rng_kernel, Nnorm, Tnorm,
121+
1. / omega_c, mesh_length * omega_c, marker_species,
115122
electron_species, neutral_species);
116123
}
117124
else if (this->vdim == 3)
118125
{
119126
reaction = recombination_reaction_amjuel<3>(
120-
this->sycl_target, rng_kernel, marker_species,
127+
this->sycl_target, rng_kernel, Nnorm, Tnorm,
128+
1. / omega_c, mesh_length * omega_c, marker_species,
121129
electron_species, neutral_species);
122130
}
123131
}
@@ -141,29 +149,33 @@ void ReactionSystem::set_up_reactions()
141149
if (this->vdim == 2)
142150
{
143151
reaction = cx_reaction_fixed<2>(
144-
this->sycl_target, rng_kernel, target_species,
145-
projectile_species, rate, cross_section);
152+
this->sycl_target, rng_kernel, mesh_length * omega_c,
153+
target_species, projectile_species, rate,
154+
cross_section);
146155
}
147156
else if (this->vdim == 3)
148157
{
149158
reaction = cx_reaction_fixed<3>(
150-
this->sycl_target, rng_kernel, target_species,
151-
projectile_species, rate, cross_section);
159+
this->sycl_target, rng_kernel, mesh_length * omega_c,
160+
target_species, projectile_species, rate,
161+
cross_section);
152162
}
153163
}
154164
else if (std::get<2>(v).first == "AMJUEL")
155165
{
156166
if (this->vdim == 2)
157167
{
158-
reaction = cx_reaction_amjuel<2>(this->sycl_target,
159-
rng_kernel, target_species,
160-
projectile_species);
168+
reaction = cx_reaction_amjuel<2>(
169+
this->sycl_target, rng_kernel, Nnorm, Tnorm,
170+
1. / omega_c, mesh_length * omega_c, target_species,
171+
projectile_species);
161172
}
162173
else if (this->vdim == 3)
163174
{
164-
reaction = cx_reaction_amjuel<3>(this->sycl_target,
165-
rng_kernel, target_species,
166-
projectile_species);
175+
reaction = cx_reaction_amjuel<3>(
176+
this->sycl_target, rng_kernel, Nnorm, Tnorm,
177+
1. / omega_c, mesh_length * omega_c, target_species,
178+
projectile_species);
167179
}
168180
}
169181
}

‎src/ParticleSystems/Reactions.hpp‎

Lines changed: 25 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -271,12 +271,13 @@ inline auto thermal_reflection(SYCLTargetSharedPtr sycl_target,
271271
}
272272

273273
template <size_t vdim>
274-
inline auto ionise_reaction_amjuel(SYCLTargetSharedPtr sycl_target,
274+
inline auto ionise_reaction_amjuel(SYCLTargetSharedPtr sycl_target, double dens,
275+
double temp, double time, double vel,
275276
const Species &target_species,
276277
const Species &electron_species)
277278
{
278-
auto ionise_rate_data = AMJUEL::ionise_rate_data();
279-
auto ionise_energy_data = AMJUEL::ionise_energy_data();
279+
auto ionise_rate_data = AMJUEL::ionise_rate_data(dens, temp, time);
280+
auto ionise_energy_data = AMJUEL::ionise_energy_data(dens, temp, time, vel);
280281

281282
auto ionise_reaction = std::make_shared<ElectronImpactIonisation<
282283
decltype(ionise_rate_data), decltype(ionise_energy_data), vdim>>(
@@ -305,19 +306,21 @@ inline auto ionise_reaction_fixed(SYCLTargetSharedPtr sycl_target,
305306
template <size_t vdim>
306307
inline auto cx_reaction_amjuel(
307308
SYCLTargetSharedPtr sycl_target,
308-
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel,
309-
const Species &parent_species, const Species &descendant_species)
309+
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel, double dens,
310+
double temp, double time, double vel, const Species &parent_species,
311+
const Species &descendant_species)
310312
{
311313
auto parent_mass = parent_species.get_mass();
312314
auto child_mass = descendant_species.get_mass();
313315
auto reduced_mass = (parent_mass * child_mass) / (parent_mass + child_mass);
314-
auto rate_data = AMJUEL::cx_rate_data(parent_mass, child_mass);
315-
auto cross_section = AMJUEL::amjuel_fit_cross_section(reduced_mass);
316+
auto rate_data =
317+
AMJUEL::cx_rate_data(parent_mass, child_mass, dens, temp, time, vel);
318+
auto cross_section = AMJUEL::amjuel_fit_cross_section(reduced_mass, vel);
316319

317320
auto data_calc_sampler =
318321
FilteredMaxwellianSampler<vdim, decltype(cross_section)>(
319322
(constants::temp_SI * constants::k_B) /
320-
(child_mass * norm::mass_amu_SI * norm::vel * norm::vel),
323+
(child_mass * norm::mass_amu_SI * vel * vel),
321324
cross_section, rng_kernel);
322325

323326
auto data_calculator =
@@ -346,7 +349,7 @@ inline auto cx_reaction_amjuel(
346349
template <size_t vdim>
347350
inline auto cx_reaction_fixed(
348351
SYCLTargetSharedPtr sycl_target,
349-
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel,
352+
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel, double vel,
350353
const Species &parent_species, const Species &descendant_species, REAL rate,
351354
REAL sigma)
352355
{
@@ -359,7 +362,7 @@ inline auto cx_reaction_fixed(
359362
auto data_calc_sampler =
360363
FilteredMaxwellianSampler<vdim, decltype(cross_section)>(
361364
(constants::temp_SI * constants::k_B) /
362-
(child_mass * norm::mass_amu_SI * norm::vel * norm::vel),
365+
(child_mass * norm::mass_amu_SI * vel * vel),
363366
cross_section, rng_kernel);
364367

365368
auto data_calculator =
@@ -388,27 +391,29 @@ inline auto cx_reaction_fixed(
388391
template <size_t vdim>
389392
inline auto recombination_reaction_amjuel(
390393
SYCLTargetSharedPtr sycl_target,
391-
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel,
392-
const Species &marker_species, const Species &electron_species,
393-
const Species &neutral_species)
394+
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel, double dens,
395+
double temp, double time, double vel, const Species &marker_species,
396+
const Species &electron_species, const Species &neutral_species)
394397
{
395-
auto recomb_data = AMJUEL::recomb_rate_data();
396-
auto recomb_energy_data = AMJUEL::recomb_energy_data();
398+
auto recomb_data = AMJUEL::recomb_rate_data(dens, temp, time);
399+
auto recomb_energy_data = AMJUEL::recomb_energy_data(dens, temp, time, vel);
397400

398401
auto constant_rate_cross_section = ConstantRateCrossSection(1.0);
399402
auto recomb_data_calc_sampler =
400403
FilteredMaxwellianSampler<vdim, decltype(constant_rate_cross_section)>(
401404
(constants::temp_SI * constants::k_B) /
402-
(marker_species.get_mass() * norm::mass_amu_SI * norm::vel *
403-
norm::vel),
405+
(marker_species.get_mass() * norm::mass_amu_SI * vel * vel),
404406
constant_rate_cross_section, rng_kernel);
405407
auto recomb_data_calc_obj =
406408
DataCalculator<decltype(recomb_energy_data),
407409
decltype(recomb_data_calc_sampler)>(
408410
recomb_energy_data, recomb_data_calc_sampler);
409411

412+
double potential_energy =
413+
13.6 * constants::e / (norm::mass_amu_SI * vel * vel);
414+
410415
auto recomb_reaction_kernel = RecombReactionKernels<vdim>(
411-
marker_species, electron_species, norm::potential_energy);
416+
marker_species, electron_species, potential_energy);
412417

413418
const int out_state = neutral_species.get_id();
414419
std::array<int, 1> recomb_out_states = {out_state};
@@ -425,7 +430,7 @@ inline auto recombination_reaction_amjuel(
425430
template <size_t vdim>
426431
inline auto recombination_reaction_fixed(
427432
SYCLTargetSharedPtr sycl_target,
428-
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel,
433+
std::shared_ptr<HostAtomicBlockKernelRNG<REAL>> rng_kernel, double vel,
429434
const Species &marker_species, const Species &electron_species,
430435
const Species &neutral_species, REAL rate, REAL energy_rate)
431436
{
@@ -436,8 +441,7 @@ inline auto recombination_reaction_fixed(
436441
auto recomb_data_calc_sampler =
437442
FilteredMaxwellianSampler<vdim, decltype(constant_rate_cross_section)>(
438443
(constants::temp_SI * constants::k_B) /
439-
(marker_species.get_mass() * norm::mass_amu_SI * norm::vel *
440-
norm::vel),
444+
(marker_species.get_mass() * norm::mass_amu_SI * vel * vel),
441445
constant_rate_cross_section, rng_kernel);
442446
auto recomb_data_calc_obj =
443447
DataCalculator<decltype(recomb_energy_data),

0 commit comments

Comments
 (0)