Skip to content

Commit b864ded

Browse files
StaticObserverStaticObserver
committed
working version for photon decay
1 parent 9daf274 commit b864ded

2 files changed

Lines changed: 355 additions & 18 deletions

File tree

‎src/kernels/QED_process.hpp‎

Lines changed: 192 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ namespace kernel::QED{
3434
array_t<prtldx_t*> dx1_ph;
3535

3636
const array_t<size_t*> N_phs;
37-
const array_t<int*> offsets;
37+
const array_t<size_t*> offsets;
3838

3939
const real_t e_min;
4040
const real_t gamma_emit;
@@ -52,7 +52,7 @@ namespace kernel::QED{
5252
const real_t rho_,
5353
const size_t N_max,
5454
const array_t<size_t*>& N_phs_,
55-
const array_t<int*>& offsets_,
55+
const array_t<size_t*>& offsets_,
5656
random_number_pool_t& random_pool_)
5757
: ux1 { charges.ux1 }
5858
, ux2 { charges.ux2 }
@@ -119,6 +119,7 @@ namespace kernel::QED{
119119
ux1_ph(offsets(p) + i) = SIGN(ux1(p)) * ONE;
120120
pld_ph(offsets(p) + i, 0) = inverseCDF(CDF(zeta) * Random<real_t>(rand_gen))
121121
* CUBE(pp / gamma_emit) / rho;
122+
pld_ph(offsets(p) + i, 1) = ZERO;
122123
i1_ph(offsets(p) + i) = i1(p);
123124
dx1_ph(offsets(p) + i) = dx1(p);
124125
weight_ph(offsets(p) + i) = w_crect * weight(p);
@@ -298,8 +299,9 @@ namespace kernel::QED{
298299
auto rand_gen = random_pool.get_state();
299300

300301
ux1_ph(offset + i) = SIGN(ux1(p)) * ONE;
301-
pld_ph(offset + i, 0) = inverseCDF(CDF(zeta) * (Random<real_t>(rand_gen) + 1e-12))
302+
pld_ph(offset + i, 0) = inverseCDF(CDF(zeta) * Random<real_t>(rand_gen))
302303
* CUBE(pp / gamma_emit) / rho;
304+
pld_ph(offset + i, 1) = ZERO;
303305
i1_ph(offset + i) = i1(p);
304306
dx1_ph(offset + i) = dx1(p);
305307
weight_ph(offset + i) = w_crect * weight(p);
@@ -309,6 +311,193 @@ namespace kernel::QED{
309311
}
310312
}
311313
};
314+
315+
template <Dimension D, Coord::type C>
316+
class PairCreation_kernel{
317+
static_assert(D == Dim::_1D, "Pair creation is only implemented in 1D");
318+
static_assert(C == Coord::Cart, "Pair creation is only implemented in cartesian coordinates");
319+
320+
const array_t<real_t*> ux1_ph, ux2_ph, ux3_ph;
321+
const array_t<real_t**> pld_ph;
322+
const array_t<short*> tag_ph;
323+
const array_t<int*> i1_ph;
324+
const array_t<prtldx_t*> dx1_ph;
325+
326+
array_t<bool*> should_inj;
327+
328+
const real_t coeff1;
329+
const real_t coeff2;
330+
const real_t rho0;
331+
const real_t L;
332+
333+
const real_t dt;
334+
const size_t n_steps;
335+
336+
Inline auto Rho(const real_t x) const -> real_t{
337+
return rho0 * (ONE + 0.8 * x / L);
338+
}
339+
340+
341+
Inline auto integrate(const index_t p, const size_t n_steps) const -> real_t{
342+
const real_t dx { ux1_ph(p) * dt / n_steps };
343+
344+
const real_t rho { Rho(static_cast<real_t>(i1_ph(p)) + static_cast<real_t>(dx1_ph(p))) };
345+
346+
real_t sum { math::sin(pld_ph(p, 2)) * exp(-coeff2 / (math::sin(pld_ph(p, 2)) * pld_ph(p, 0)))
347+
+ math::sin(pld_ph(p, 2) + ux1_ph(p) * dt / rho) * exp(-coeff2 / (math::sin(pld_ph(p, 2) + ux1_ph(p) * dt / rho) * pld_ph(p, 0))) };
348+
349+
for (size_t i = 1; i < n_steps; i += 2){
350+
sum += FOUR * math::sin(pld_ph(p, 2) + i * dx / rho) * exp(-coeff2 / (math::sin(pld_ph(p, 2)+ i * dx / rho) * pld_ph(p, 0)));
351+
}
352+
353+
for (size_t i = 2; i < n_steps; i += 2){
354+
sum += TWO * math::sin(pld_ph(p, 2) + i * dx / rho) * exp(-coeff2 / (math::sin(pld_ph(p, 2) + i * dx / rho) * pld_ph(p, 0)));
355+
}
356+
357+
return dx / THREE * sum;
358+
359+
}
360+
361+
public:
362+
PairCreation_kernel(Particles<D, C>& photons,
363+
const real_t coeff1_,
364+
const real_t coeff2_,
365+
const real_t rho0_,
366+
const real_t L_,
367+
const real_t dt_,
368+
const size_t n_steps_,
369+
array_t<bool*>& should_inj_)
370+
: ux1_ph { photons.ux1 }
371+
, ux2_ph { photons.ux2 }
372+
, ux3_ph { photons.ux3 }
373+
, pld_ph { photons.pld }
374+
, tag_ph { photons.tag }
375+
, i1_ph { photons.i1 }
376+
, dx1_ph { photons.dx1 }
377+
, coeff1 { coeff1_ }
378+
, coeff2 { coeff2_ }
379+
, rho0 { rho0_ }
380+
, L { L_ }
381+
, dt { dt_ }
382+
, n_steps { n_steps_ }
383+
, should_inj { should_inj_ } {}
384+
~PairCreation_kernel() = default;
385+
386+
Inline void operator()(index_t p) const{
387+
if (tag_ph(p) != ParticleTag::alive){
388+
if (tag_ph(p) != ParticleTag::dead){
389+
raise::KernelError(HERE, "Invalid particle tag in pusher");
390+
}
391+
return;
392+
}
393+
pld_ph(p, 1) += coeff1 * integrate(p, n_steps);
394+
pld_ph(p, 2) += dt * ux1_ph(p) / Rho(static_cast<real_t>(i1_ph(p)) + static_cast<real_t>(dx1_ph(p)));
395+
if (pld_ph(p, 0) * math::sin(pld_ph(p, 1)) < TWO){
396+
tag_ph(p) = ParticleTag::dead;
397+
return;
398+
}
399+
if (pld_ph(p, 1) > ONE){
400+
should_inj(p) = true;
401+
}
402+
}
403+
404+
}; // class PairCreation_kernel
405+
406+
template <Dimension D, Coord::type C>
407+
class InjectPairs_kernel{
408+
static_assert(D == Dim::_1D, "Pair creation is only implemented in 1D");
409+
static_assert(C == Coord::Cart, "Pair creation is only implemented in cartesian coordinates");
410+
411+
const array_t<real_t*> ux1_ph, ux2_ph, ux3_ph;
412+
const array_t<real_t**> pld_ph;
413+
const array_t<int*> i1_ph;
414+
const array_t<prtldx_t*> dx1_ph;
415+
const array_t<short*> tag_ph;
416+
const array_t<real_t*> weight_ph;
417+
418+
array_t<real_t*> ux1_p, ux2_p, ux3_p;
419+
array_t<real_t*> weight_p;
420+
array_t<short*> tag_p;
421+
array_t<int*> i1_p;
422+
array_t<prtldx_t*> dx1_p;
423+
424+
array_t<real_t*> ux1_e, ux2_e, ux3_e;
425+
array_t<real_t*> weight_e;
426+
array_t<short*> tag_e;
427+
array_t<int*> i1_e;
428+
array_t<prtldx_t*> dx1_e;
429+
430+
array_t<bool*> should_inj;
431+
array_t<size_t*> offsets_e;
432+
array_t<size_t*> offsets_p;
433+
434+
public:
435+
436+
InjectPairs_kernel(Particles<D, C>& photon,
437+
Particles<D, C>& electron,
438+
Particles<D, C>& positron,
439+
array_t<bool*>& should_inj_,
440+
array_t<size_t*>& offsets_e_,
441+
array_t<size_t*>& offsets_p_)
442+
: ux1_ph { photon.ux1 }
443+
, ux2_ph { photon.ux2 }
444+
, ux3_ph { photon.ux3 }
445+
, pld_ph { photon.pld }
446+
, i1_ph { photon.i1 }
447+
, dx1_ph { photon.dx1 }
448+
, tag_ph { photon.tag }
449+
, weight_ph { photon.weight }
450+
, ux1_p { positron.ux1 }
451+
, ux2_p { positron.ux2 }
452+
, ux3_p { positron.ux3 }
453+
, i1_p { positron.i1 }
454+
, dx1_p { positron.dx1 }
455+
, tag_p { positron.tag }
456+
, weight_p { positron.weight }
457+
, ux1_e { electron.ux1 }
458+
, ux2_e { electron.ux2 }
459+
, ux3_e { electron.ux3 }
460+
, i1_e { electron.i1 }
461+
, dx1_e { electron.dx1 }
462+
, tag_e { electron.tag }
463+
, weight_e { electron.weight }
464+
, should_inj { should_inj_ }
465+
, offsets_e { offsets_e_ }
466+
, offsets_p { offsets_p_ } {}
467+
468+
~InjectPairs_kernel() = default;
469+
470+
Inline void operator()(index_t p) const{
471+
if (!should_inj(p)){
472+
return;
473+
}
474+
const real_t u = math::abs(math::cos(pld_ph(p, 2))) * math::cos(pld_ph(p, 2))
475+
* math::sqrt(SQR(pld_ph(p, 0)) - FOUR)
476+
/ math::sqrt(SQR(math::sin(pld_ph(p, 2) * pld_ph(p, 0)))
477+
+ FOUR * math::cos(pld_ph(p, 2)));
478+
479+
ux1_p(offsets_p(p)) = u;
480+
ux2_p(offsets_p(p)) = ZERO;
481+
ux3_p(offsets_p(p)) = ZERO;
482+
i1_p(offsets_p(p)) = i1_ph(p);
483+
dx1_p(offsets_p(p)) = dx1_ph(p);
484+
weight_p(offsets_p(p)) = weight_ph(p);
485+
tag_p(offsets_p(p)) = ParticleTag::alive;
486+
487+
ux1_e(offsets_e(p)) = u;
488+
ux2_e(offsets_e(p)) = ZERO;
489+
ux3_e(offsets_e(p)) = ZERO;
490+
i1_e(offsets_e(p)) = i1_ph(p);
491+
dx1_e(offsets_e(p)) = dx1_ph(p);
492+
weight_e(offsets_e(p)) = weight_ph(p);
493+
tag_e(offsets_e(p)) = ParticleTag::alive;
494+
495+
tag_ph(p) = ParticleTag::dead;
496+
}
497+
498+
}; // class InjectPairs_kernel
499+
500+
312501
} // namespace QED
313502

314503

0 commit comments

Comments
 (0)