Skip to content

Commit a59888b

Browse files
StaticObserverStaticObserver
committed
minor changes
1 parent 2b3a1f6 commit a59888b

1 file changed

Lines changed: 132 additions & 123 deletions

File tree

‎src/kernels/QED_process.hpp‎

Lines changed: 132 additions & 123 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<size_t*> offsets;
37+
const array_t<size_t*> offsets;
3838

3939
const real_t e_min;
4040
const real_t gamma_emit;
@@ -143,6 +143,7 @@ namespace kernel::QED{
143143
const real_t e_min;
144144
const real_t coeff;
145145
const real_t gamma_emit;
146+
const real_t gamma_min;
146147
const real_t rho;
147148

148149

@@ -155,6 +156,7 @@ namespace kernel::QED{
155156
const real_t e_min_,
156157
const real_t coeff_,
157158
const real_t gamma_emit_,
159+
const real_t gamma_min_,
158160
const real_t rho_,
159161
array_t<size_t*>& N_phs_)
160162
: ux1 { charge.ux1 }
@@ -164,6 +166,7 @@ namespace kernel::QED{
164166
, e_min { e_min_ }
165167
, coeff { coeff_ }
166168
, gamma_emit { gamma_emit_ }
169+
, gamma_min { gamma_min_ }
167170
, rho { rho_ }
168171
, N_phs { N_phs_ } {
169172
Kokkos::deep_copy(N_phs, 0);
@@ -181,6 +184,10 @@ namespace kernel::QED{
181184
return;
182185
}
183186
const real_t pp = math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p)));
187+
if (pp < gamma_min){
188+
N_phs(p) = 0;
189+
return;
190+
}
184191
const real_t zeta = e_min * rho * CUBE(gamma_emit / pp);
185192
auto N_ph = static_cast<size_t>(coeff * CDF(zeta) / SQR(pp));
186193

@@ -189,128 +196,7 @@ namespace kernel::QED{
189196
};
190197

191198

192-
template <Dimension D, Coord::type C>
193-
class AtomicCurvatureEmission_kernel{
194-
static_assert(D == Dim::_1D, "Curvature emission is only implemented in 1D");
195-
static_assert(C == Coord::Cart, "Curvature emission is only implemented in cartesian coordinates");
196-
197-
const array_t<real_t*> ux1, ux2, ux3;
198-
const array_t<real_t*> weight;
199-
const array_t<short*> tag;
200-
const array_t<int*> i1;
201-
const array_t<prtldx_t*> dx1;
202-
203-
array_t<real_t*> ux1_ph, ux2_ph, ux3_ph;
204-
array_t<real_t*> weight_ph;
205-
array_t<short*> tag_ph;
206-
array_t<real_t**> pld_ph;
207-
array_t<int*> i1_ph;
208-
array_t<prtldx_t*> dx1_ph;
209-
210-
const real_t e_min;
211-
const real_t gamma_emit;
212-
const real_t coeff;
213-
const real_t rho;
214-
const size_t N_max;
215-
const size_t npart_ph;
216-
random_number_pool_t random_pool;
217-
218-
array_t<size_t> n_inj { "n_inj" };
219-
220-
221-
public:
222-
AtomicCurvatureEmission_kernel(Particles<D, C>& charges,
223-
Particles<D, C>& photons,
224-
const real_t e_min_,
225-
const real_t gamma_emit_,
226-
const real_t coeff_,
227-
const real_t rho_,
228-
const size_t N_max,
229-
const size_t npart_ph_,
230-
random_number_pool_t& random_pool_)
231-
: ux1 { charges.ux1 }
232-
, ux2 { charges.ux2 }
233-
, ux3 { charges.ux3 }
234-
, weight { charges.weight }
235-
, tag { charges.tag }
236-
, i1 { charges.i1 }
237-
, dx1 { charges.dx1 }
238-
, ux1_ph { photons.ux1 }
239-
, ux2_ph { photons.ux2 }
240-
, ux3_ph { photons.ux3 }
241-
, weight_ph { photons.weight }
242-
, tag_ph { photons.tag }
243-
, pld_ph { photons.pld }
244-
, i1_ph { photons.i1 }
245-
, dx1_ph { photons.dx1 }
246-
, e_min { e_min_ }
247-
, gamma_emit { gamma_emit_ }
248-
, coeff { coeff_ }
249-
, rho { rho_ }
250-
, N_max { N_max }
251-
, npart_ph { npart_ph_ }
252-
, random_pool { random_pool_ } {
253-
Kokkos::deep_copy(n_inj, 0);
254-
}
255-
~AtomicCurvatureEmission_kernel() = default;
256-
257-
Inline auto CDF(real_t zeta_) const -> real_t{
258-
return math::exp(-zeta_);
259-
}
260-
261-
Inline auto inverseCDF(real_t u) const -> real_t{
262-
return -math::log(u);
263-
}
264-
265-
auto num_injected() const -> size_t{
266-
auto n_inj_h = Kokkos::create_mirror_view(n_inj);
267-
Kokkos::deep_copy(n_inj_h, n_inj);
268-
return n_inj_h();
269-
}
270-
271-
Inline void operator()(index_t p) const{
272-
if(tag(p) != ParticleTag::alive){
273-
if(tag(p) != ParticleTag::dead){
274-
raise::KernelError(HERE, "Invalid particle tag in pusher");
275-
}
276-
return;
277-
}
278-
const real_t pp = math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p)));
279-
const real_t zeta = e_min * rho * CUBE(gamma_emit / pp);
280-
auto N_ph = static_cast<size_t>(coeff * CDF(zeta) / SQR(pp));
281-
282-
if (N_ph < 1){
283-
return;
284-
}
285-
// weight correction
286-
auto w_crect = ONE;
287-
if (N_ph > N_max){
288-
w_crect *= N_ph / N_max;
289-
N_ph = N_max;
290-
}
291-
292-
auto offset = Kokkos::atomic_fetch_add(&n_inj(), N_ph);
293-
offset += npart_ph;
294-
if (offset + N_ph - 1 > ux1_ph.extent(0)){
295-
raise::KernelError(HERE, "Number of photons exceeds maxnpart.");
296-
}
297-
298-
for (short i = 0; i < N_ph; ++i){
299-
auto rand_gen = random_pool.get_state();
300-
301-
ux1_ph(offset + i) = SIGN(ux1(p)) * ONE;
302-
pld_ph(offset + i, 0) = inverseCDF(CDF(zeta) * Random<real_t>(rand_gen))
303-
* CUBE(pp / gamma_emit) / rho;
304-
pld_ph(offset + i, 1) = ZERO;
305-
i1_ph(offset + i) = i1(p);
306-
dx1_ph(offset + i) = dx1(p);
307-
weight_ph(offset + i) = w_crect * weight(p);
308-
tag_ph(offset + i) = ParticleTag::alive;
309199

310-
random_pool.free_state(rand_gen);
311-
}
312-
}
313-
};
314200

315201
template <Dimension D, Coord::type C>
316202
class PairCreation_kernel{
@@ -498,4 +384,127 @@ namespace kernel::QED{
498384
} // namespace QED
499385

500386

501-
#endif // QED_PROCESS_HPP
387+
#endif // QED_PROCESS_HPP
388+
389+
// template <Dimension D, Coord::type C>
390+
// class AtomicCurvatureEmission_kernel{
391+
// static_assert(D == Dim::_1D, "Curvature emission is only implemented in 1D");
392+
// static_assert(C == Coord::Cart, "Curvature emission is only implemented in cartesian coordinates");
393+
394+
// const array_t<real_t*> ux1, ux2, ux3;
395+
// const array_t<real_t*> weight;
396+
// const array_t<short*> tag;
397+
// const array_t<int*> i1;
398+
// const array_t<prtldx_t*> dx1;
399+
400+
// array_t<real_t*> ux1_ph, ux2_ph, ux3_ph;
401+
// array_t<real_t*> weight_ph;
402+
// array_t<short*> tag_ph;
403+
// array_t<real_t**> pld_ph;
404+
// array_t<int*> i1_ph;
405+
// array_t<prtldx_t*> dx1_ph;
406+
407+
// const real_t e_min;
408+
// const real_t gamma_emit;
409+
// const real_t coeff;
410+
// const real_t rho;
411+
// const size_t N_max;
412+
// const size_t npart_ph;
413+
// random_number_pool_t random_pool;
414+
415+
// array_t<size_t> n_inj { "n_inj" };
416+
417+
418+
// public:
419+
// AtomicCurvatureEmission_kernel(Particles<D, C>& charges,
420+
// Particles<D, C>& photons,
421+
// const real_t e_min_,
422+
// const real_t gamma_emit_,
423+
// const real_t coeff_,
424+
// const real_t rho_,
425+
// const size_t N_max,
426+
// const size_t npart_ph_,
427+
// random_number_pool_t& random_pool_)
428+
// : ux1 { charges.ux1 }
429+
// , ux2 { charges.ux2 }
430+
// , ux3 { charges.ux3 }
431+
// , weight { charges.weight }
432+
// , tag { charges.tag }
433+
// , i1 { charges.i1 }
434+
// , dx1 { charges.dx1 }
435+
// , ux1_ph { photons.ux1 }
436+
// , ux2_ph { photons.ux2 }
437+
// , ux3_ph { photons.ux3 }
438+
// , weight_ph { photons.weight }
439+
// , tag_ph { photons.tag }
440+
// , pld_ph { photons.pld }
441+
// , i1_ph { photons.i1 }
442+
// , dx1_ph { photons.dx1 }
443+
// , e_min { e_min_ }
444+
// , gamma_emit { gamma_emit_ }
445+
// , coeff { coeff_ }
446+
// , rho { rho_ }
447+
// , N_max { N_max }
448+
// , npart_ph { npart_ph_ }
449+
// , random_pool { random_pool_ } {
450+
// Kokkos::deep_copy(n_inj, 0);
451+
// }
452+
// ~AtomicCurvatureEmission_kernel() = default;
453+
454+
// Inline auto CDF(real_t zeta_) const -> real_t{
455+
// return math::exp(-zeta_);
456+
// }
457+
458+
// Inline auto inverseCDF(real_t u) const -> real_t{
459+
// return -math::log(u);
460+
// }
461+
462+
// auto num_injected() const -> size_t{
463+
// auto n_inj_h = Kokkos::create_mirror_view(n_inj);
464+
// Kokkos::deep_copy(n_inj_h, n_inj);
465+
// return n_inj_h();
466+
// }
467+
468+
// Inline void operator()(index_t p) const{
469+
// if(tag(p) != ParticleTag::alive){
470+
// if(tag(p) != ParticleTag::dead){
471+
// raise::KernelError(HERE, "Invalid particle tag in pusher");
472+
// }
473+
// return;
474+
// }
475+
// const real_t pp = math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p)));
476+
// const real_t zeta = e_min * rho * CUBE(gamma_emit / pp);
477+
// auto N_ph = static_cast<size_t>(coeff * CDF(zeta) / SQR(pp));
478+
479+
// if (N_ph < 1){
480+
// return;
481+
// }
482+
// // weight correction
483+
// auto w_crect = ONE;
484+
// if (N_ph > N_max){
485+
// w_crect *= N_ph / N_max;
486+
// N_ph = N_max;
487+
// }
488+
489+
// auto offset = Kokkos::atomic_fetch_add(&n_inj(), N_ph);
490+
// offset += npart_ph;
491+
// if (offset + N_ph - 1 > ux1_ph.extent(0)){
492+
// raise::KernelError(HERE, "Number of photons exceeds maxnpart.");
493+
// }
494+
495+
// for (short i = 0; i < N_ph; ++i){
496+
// auto rand_gen = random_pool.get_state();
497+
498+
// ux1_ph(offset + i) = SIGN(ux1(p)) * ONE;
499+
// pld_ph(offset + i, 0) = inverseCDF(CDF(zeta) * Random<real_t>(rand_gen))
500+
// * CUBE(pp / gamma_emit) / rho;
501+
// pld_ph(offset + i, 1) = ZERO;
502+
// i1_ph(offset + i) = i1(p);
503+
// dx1_ph(offset + i) = dx1(p);
504+
// weight_ph(offset + i) = w_crect * weight(p);
505+
// tag_ph(offset + i) = ParticleTag::alive;
506+
507+
// random_pool.free_state(rand_gen);
508+
// }
509+
// }
510+
// };

0 commit comments

Comments
 (0)