Skip to content

Commit bdb6354

Browse files
StaticObserverStaticObserver
committed
adding curvature drag; combine QED kernal into pgen
1 parent fefe8bf commit bdb6354

File tree

5 files changed

+194
-81
lines changed

5 files changed

+194
-81
lines changed

‎pgens/1d_polar_cap/pgen.hpp‎

Lines changed: 72 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
#include "archetypes/particle_injector.h"
1414
#include "archetypes/problem_generator.h"
1515
#include "framework/domain/metadomain.h"
16+
#include "kernels/QED_process.hpp"
1617

1718
namespace user {
1819
using namespace ntt;
@@ -141,7 +142,8 @@ namespace user {
141142
private:
142143
const real_t J0;
143144
};
144-
145+
146+
145147
template <SimEngine::type S, class M>
146148
struct PGen : public arch::ProblemGenerator<S, M> {
147149
// compatibility traits for the problem generator
@@ -156,12 +158,17 @@ namespace user {
156158
using arch::ProblemGenerator<S, M>::C;
157159
using arch::ProblemGenerator<S, M>::params;
158160

161+
159162
const real_t b0, skindepth0, larmor0;
160163
const real_t temp;
161164
const real_t j0;
162165
const real_t xsurf, ds;
163166
InitFields<D> init_flds;
164167
MagnetosphericCurrent<D> ext_current;
168+
169+
const kernel::QED::cdfTable cdf;
170+
random_number_pool_t random_pool(12345);
171+
const real_t e_min, gamma_emit, coeff, rho, N_max, dt;
165172

166173

167174
inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
@@ -172,6 +179,15 @@ namespace user {
172179
, temp { p.template get<real_t>("setup.temp") }
173180
, j0 { p.template get<real_t>("setup.j0") / p.template get<real_t>("scales.V0")}
174181
, ds { p.template get<real_t>("grid.boundaries.atmosphere.ds") }
182+
, dt { p.template get<real_t>("algorithms.timestep.dt") }
183+
, cdf { "cdf_table.txt", "inverse_cdf_table.txt" }
184+
, e_min { p.template get<real_t>("setup.e_min") }
185+
, gamma_emit { p.template get<real_t>("setup.gamma_emit") }
186+
, coeff { p.template get<real_t>("setup.coeff") }
187+
, rho { p.template get<real_t>("setup.rho") }
188+
, N_max { p.template get<real_t>("setup.N_max") }
189+
, coeff1 { coeff * 0.23 * constant::PI * b0 * constant::SQRT3 }
190+
, coeff2 { FOUR * TWO / THREE / b0 }
175191
, xsurf([&m, this]() {
176192
const auto min_buff = params.template get<unsigned short>("algorithms.current_filters") + 2;
177193
const auto buffer_ncells = min_buff > 5 ? min_buff : 5;
@@ -199,17 +215,6 @@ namespace user {
199215
const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
200216
local_domain.random_pool,
201217
temp);
202-
// const auto injector = arch::UniformInjector<S, M, arch::Maxwellian>(
203-
// energy_dist,
204-
// { 1, 2 }
205-
// );
206-
207-
// arch::InjectUniform<S, M, decltype(injector)>(
208-
// params,
209-
// local_domain,
210-
// injector,
211-
// ONE
212-
// );
213218
const auto spatial_dist = TargetDensityProfile<S, M>(
214219
local_domain.mesh.metric,
215220
params.template get<real_t>("grid.boundaries.atmosphere.density"),
@@ -243,6 +248,61 @@ namespace user {
243248
injector_extra_charge,
244249
TWO);
245250
}
251+
252+
void CustomPostStep(std::size_t, long double time, Domain<S, M>& local_domain) {
253+
kernel::QED::CurvatureEmission_kernel<D, C> curvature_emission1(local_domain.species[0],
254+
local_domain.species[2],
255+
e_min,
256+
gamma_emit,
257+
coeff * dt / skindepth0,
258+
rho,
259+
N_max,
260+
random_pool,
261+
cdf);
262+
Kokkos::parallel_for("CurvatureEmission",
263+
local_domain.species[0].rangeActiveParticles(),
264+
curvature_emission1);
265+
Kokkos::fence();
266+
auto n_injected = curvature_emission1.num_injected();
267+
local_domain.species[2].set_npart(local_domain.species[2].npart() + n_injected);
268+
kernel::QED::CurvatureEmission_kernel<D, C> curvature_emission2(local_domain.species[1],
269+
local_domain.species[2],
270+
e_min,
271+
gamma_emit,
272+
coeff * dt / skindepth0,
273+
rho,
274+
N_max,
275+
random_pool,
276+
cdf);
277+
Kokkos::parallel_for("CurvatureEmission",
278+
local_domain.species[1].rangeActiveParticles(),
279+
curvature_emission2);
280+
Kokkos::fence();
281+
n_injected = curvature_emission2.num_injected();
282+
local_domain.species[2].set_npart(local_domain.species[2].npart() + n_injected);
283+
284+
kernel::QED::PairCreation_kernel<D, C> pair_creation(local_domain.species[2],
285+
local_domain.species[1],
286+
local_domain.species[0]);
287+
288+
Kokkos::fence();
289+
290+
n_injected = pair_creation.num_injected();
291+
local_domain.species[1].set_npart(local_domain.species[1].npart() + n_injected);
292+
local_domain.species[0].set_npart(local_domain.species[0].npart() + n_injected);
293+
294+
Kokkos::parallel_for("PayloadUpdate",
295+
local_domain.species[2].rangeActiveParticles(),
296+
PayloadUpdate<D, C>(local_domain.species[2],
297+
coeff1 / skindepth0,
298+
coeff2,
299+
dt,
300+
rho,
301+
5));
302+
303+
}
304+
305+
246306
}; // PGen
247307

248308

‎src/engines/srpic.hpp‎

Lines changed: 12 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -328,6 +328,11 @@ namespace ntt {
328328
"scales.omegaB0") /
329329
(SQR(sync_grad) * species.mass())
330330
: ZERO;
331+
const auto has_curvature = (cooling == Cooling::CURVATURE);
332+
const auto cur_coeff = has_curvature
333+
? m_params.template get<real_t>(
334+
"algorithms.curvature.coeff")
335+
: ZERO;
331336

332337
// toggle to indicate whether pgen defines the external force
333338
bool has_extforce = false;
@@ -346,6 +351,9 @@ namespace ntt {
346351
if (cooling == Cooling::SYNCHROTRON) {
347352
cooling_tags = kernel::sr::Cooling::Synchrotron;
348353
}
354+
if (cooling == Cooling::CURVATURE) {
355+
cooling_tags = kernel::sr::Cooling::Curvature;
356+
}
349357
// clang-format off
350358
if (not has_atmosphere and not has_extforce) {
351359
Kokkos::parallel_for(
@@ -368,7 +376,7 @@ namespace ntt {
368376
domain.mesh.n_active(in::x2),
369377
domain.mesh.n_active(in::x3),
370378
domain.mesh.prtl_bc(),
371-
gca_larmor_max, gca_eovrb_max, sync_coeff
379+
gca_larmor_max, gca_eovrb_max, sync_coeff, cur_coeff
372380
));
373381
} else if (has_atmosphere and not has_extforce) {
374382
const auto force =
@@ -398,7 +406,7 @@ namespace ntt {
398406
domain.mesh.n_active(in::x2),
399407
domain.mesh.n_active(in::x3),
400408
domain.mesh.prtl_bc(),
401-
gca_larmor_max, gca_eovrb_max, sync_coeff
409+
gca_larmor_max, gca_eovrb_max, sync_coeff, cur_coeff
402410
));
403411
} else if (not has_atmosphere and has_extforce) {
404412
if constexpr (traits::has_member<traits::pgen::ext_force_t, pgen_t>::value) {
@@ -427,7 +435,7 @@ namespace ntt {
427435
domain.mesh.n_active(in::x2),
428436
domain.mesh.n_active(in::x3),
429437
domain.mesh.prtl_bc(),
430-
gca_larmor_max, gca_eovrb_max, sync_coeff
438+
gca_larmor_max, gca_eovrb_max, sync_coeff, cur_coeff
431439
));
432440
} else {
433441
raise::Error("External force not implemented", HERE);
@@ -459,7 +467,7 @@ namespace ntt {
459467
domain.mesh.n_active(in::x2),
460468
domain.mesh.n_active(in::x3),
461469
domain.mesh.prtl_bc(),
462-
gca_larmor_max, gca_eovrb_max, sync_coeff
470+
gca_larmor_max, gca_eovrb_max, sync_coeff, cur_coeff
463471
));
464472
} else {
465473
raise::Error("External force not implemented", HERE);

‎src/kernels/QED_process.hpp‎

Lines changed: 71 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -280,6 +280,68 @@ namespace kernel::QED{
280280
}
281281
}; // class CurvatureEmission_kernel
282282

283+
template <Dimension D, Coord::type C>
284+
struct PayloadUpdate {
285+
const real_t coeff1, coeff2, dt, rho;
286+
const array_t<real_t**> pld_ph;
287+
const array_t<real_t*> ux1_ph;
288+
const array_t<int*> i1_ph;
289+
const array_t<short*> tag_ph;
290+
array_t<prtldx_t*> dx1_ph;
291+
292+
const size_t n_steps;
293+
294+
295+
296+
Inline auto integrate(const index_t p, const size_t n_steps) const -> real_t{
297+
const real_t dx { ux1_ph(p) * dt / n_steps };
298+
299+
real_t sum { pld_ph(p, 2) * exp(-coeff2 / (pld_ph(p, 2) * pld_ph(p, 0)))
300+
+ (pld_ph(p, 2) + ux1_ph(p) * dt / rho) * exp(-coeff2 / ((pld_ph(p, 2) + ux1_ph(p) * dt / rho) * pld_ph(p, 0))) };
301+
302+
for (size_t i = 1; i < n_steps; i += 2){
303+
sum += FOUR * (pld_ph(p, 2) + i * dx / rho) * exp(-coeff2 / ((pld_ph(p, 2)+ i * dx / rho) * pld_ph(p, 0)));
304+
}
305+
306+
for (size_t i = 2; i < n_steps; i += 2){
307+
sum += TWO * (pld_ph(p, 2) + i * dx / rho) * exp(-coeff2 / ((pld_ph(p, 2) + i * dx / rho) * pld_ph(p, 0)));
308+
}
309+
310+
return dx / THREE * sum;
311+
312+
}
313+
314+
public:
315+
PayloadUpdate(Particles<D, C>& photon,
316+
const real_t coeff1_,
317+
const real_t coeff2_,
318+
const real_t dt_,
319+
const real_t rho_,
320+
const size_t n_steps_)
321+
: coeff1 { coeff1_ }
322+
, coeff2 { coeff2_ }
323+
, dt { dt_ }
324+
, rho { rho_ }
325+
, pld_ph { photon.pld }
326+
, ux1_ph { photon.ux1 }
327+
, i1_ph { photon.i1 }
328+
, tag_ph { photon.tag }
329+
, dx1_ph { photon.dx1 }
330+
, n_steps { n_steps_ } {}
331+
~PayloadUpdate() = default;
332+
333+
Inline void operator()(index_t p) const{
334+
if (tag_ph(p) != ParticleTag::alive){
335+
if (tag_ph(p) != ParticleTag::dead){
336+
raise::KernelError(HERE, "Invalid particle tag in pusher");
337+
}
338+
return;
339+
}
340+
pld_ph(p, 1) += coeff1 * integrate(p, n_steps);
341+
pld_ph(p, 2) += dt * ux1_ph(p) / rho;
342+
}
343+
};// class PayloadUpdate
344+
283345
template <Dimension D, Coord::type C>
284346
class PairCreation_kernel{
285347
static_assert(D == Dim::_1D, "Pair creation is only implemented in 1D");
@@ -295,7 +357,7 @@ namespace kernel::QED{
295357
array_t<real_t*> ux1_p, ux2_p, ux3_p;
296358
array_t<real_t*> weight_p;
297359
array_t<short*> tag_p;
298-
array_t<int*> i1_p;
360+
array_t<int*> i1_p;
299361
array_t<prtldx_t*> dx1_p;
300362
const size_t npart_p;
301363

@@ -306,58 +368,12 @@ namespace kernel::QED{
306368
array_t<prtldx_t*> dx1_e;
307369
const size_t npart_e;
308370

309-
const real_t coeff1;
310-
const real_t coeff2;
311-
const real_t rho0;
312-
const real_t L;
313-
314-
const real_t dt;
315-
const size_t n_steps;
316-
317371
array_t<size_t> n_inj { "n_inj" };
318372

319-
auto num_injected() const -> size_t{
320-
auto n_inj_h = Kokkos::create_mirror_view(n_inj);
321-
Kokkos::deep_copy(n_inj_h, n_inj);
322-
return n_inj_h();
323-
}
324-
325-
Inline auto Rho(const real_t x) const -> real_t{
326-
// return rho0 * (ONE + 0.8 * x / L);
327-
return rho0;
328-
}
329-
330-
331-
Inline auto integrate(const index_t p, const size_t n_steps) const -> real_t{
332-
const real_t dx { ux1_ph(p) * dt / n_steps };
333-
334-
const real_t rho { Rho(static_cast<real_t>(i1_ph(p)) + static_cast<real_t>(dx1_ph(p))) };
335-
336-
real_t sum { pld_ph(p, 2) * exp(-coeff2 / (pld_ph(p, 2) * pld_ph(p, 0)))
337-
+ (pld_ph(p, 2) + ux1_ph(p) * dt / rho) * exp(-coeff2 / ((pld_ph(p, 2) + ux1_ph(p) * dt / rho) * pld_ph(p, 0))) };
338-
339-
for (size_t i = 1; i < n_steps; i += 2){
340-
sum += FOUR * (pld_ph(p, 2) + i * dx / rho) * exp(-coeff2 / ((pld_ph(p, 2)+ i * dx / rho) * pld_ph(p, 0)));
341-
}
342-
343-
for (size_t i = 2; i < n_steps; i += 2){
344-
sum += TWO * (pld_ph(p, 2) + i * dx / rho) * exp(-coeff2 / ((pld_ph(p, 2) + i * dx / rho) * pld_ph(p, 0)));
345-
}
346-
347-
return dx / THREE * sum;
348-
349-
}
350-
351373
public:
352374
PairCreation_kernel(Particles<D, C>& photons,
353375
Particles<D, C>& positrons,
354-
Particles<D, C>& electrons,
355-
const real_t coeff1_,
356-
const real_t coeff2_,
357-
const real_t rho0_,
358-
const real_t L_,
359-
const real_t dt_,
360-
const size_t n_steps_)
376+
Particles<D, C>& electrons)
361377
: ux1_ph { photons.ux1 }
362378
, ux2_ph { photons.ux2 }
363379
, ux3_ph { photons.ux3 }
@@ -392,15 +408,20 @@ namespace kernel::QED{
392408
}
393409
~PairCreation_kernel() = default;
394410

411+
auto num_injected() const -> size_t{
412+
auto n_inj_h = Kokkos::create_mirror_view(n_inj);
413+
Kokkos::deep_copy(n_inj_h, n_inj);
414+
return n_inj_h();
415+
}
416+
417+
395418
Inline void operator()(index_t p) const{
396419
if (tag_ph(p) != ParticleTag::alive){
397420
if (tag_ph(p) != ParticleTag::dead){
398421
raise::KernelError(HERE, "Invalid particle tag in pusher");
399422
}
400423
return;
401424
}
402-
pld_ph(p, 1) += coeff1 * integrate(p, n_steps); // optical depth
403-
pld_ph(p, 2) += dt * ux1_ph(p) / Rho(static_cast<real_t>(i1_ph(p)) + static_cast<real_t>(dx1_ph(p))); // angle
404425
if (pld_ph(p, 0) * math::sin(pld_ph(p, 2)) < TWO){
405426
tag_ph(p) = ParticleTag::dead;
406427
return;

0 commit comments

Comments
 (0)