Skip to content

Commit 470b202

Browse files
StaticObserverStaticObserver
committed
bug fix
1 parent bdb6354 commit 470b202

7 files changed

Lines changed: 112 additions & 113 deletions

File tree

‎pgens/1d_polar_cap/pgen.hpp‎

Lines changed: 68 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -31,16 +31,12 @@ namespace user {
3131
}
3232

3333
Inline auto ex1(const coord_t<D>& x_Ph) const -> real_t {
34-
if (x_Ph[0] < xsurf + ds ){
34+
if (x_Ph[0] < xsurf ){
3535
return ZERO;
36-
}else if (x_Ph[0] > xsurf + 1.33 * ds){
37-
return -coeff * (x_Ph[0] - xsurf - ds * (ONE + 0.03 * math::log(1.01) - 0.33 * math::log(0.01 + math::exp(-11.0 / ds))));
38-
3936
}else{
40-
return -coeff * (x_Ph[0] - xsurf - ds
37+
return -coeff * (x_Ph[0] - xsurf
4138
+ 0.03 * ds * math::log(0.01 + math::exp((xsurf + 1.0 * ds - x_Ph[0]) / 0.03 / ds))
42-
- 0.03 * ds * math::log(0.01 + ONE));
43-
//return -coeff * (x_Ph[0] - xsurf - ds);
39+
- 0.03 * ds * math::log(0.01 + math::exp(1.0/0.03)));
4440
}
4541
}
4642

@@ -109,14 +105,10 @@ namespace user {
109105
, ds { ds_ } {}
110106

111107
Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t {
112-
if ( x_Ph[0] < xsurf or x_Ph[0] > 1.33 * ds + xsurf){
108+
if ( x_Ph[0] < xsurf){
113109
return ZERO;
114110
}else{
115-
if (x_Ph[0] < ds + xsurf){
116-
return ONE;
117-
}else{
118-
return ONE - 0.01 / (0.01 + math::exp(-(x_Ph[0] - xsurf - 1.0 * ds) / 0.03 / ds));
119-
}
111+
return ONE - 0.01 / (0.01 + math::exp(-(x_Ph[0] - xsurf - 1.0 * ds) / 0.03 / ds));
120112
}
121113
}
122114
}; // ExtraCharge
@@ -167,8 +159,11 @@ namespace user {
167159
MagnetosphericCurrent<D> ext_current;
168160

169161
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;
162+
random_number_pool_t random_pool;
163+
const real_t e_min, gamma_emit, coeff, rho, dt, coeff1, coeff2;
164+
const size_t N_max;
165+
166+
const bool QED_on;
172167

173168

174169
inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
@@ -185,18 +180,18 @@ namespace user {
185180
, gamma_emit { p.template get<real_t>("setup.gamma_emit") }
186181
, coeff { p.template get<real_t>("setup.coeff") }
187182
, rho { p.template get<real_t>("setup.rho") }
188-
, N_max { p.template get<real_t>("setup.N_max") }
183+
, N_max { static_cast<size_t>(p.template get<int>("setup.N_max")) }
189184
, coeff1 { coeff * 0.23 * constant::PI * b0 * constant::SQRT3 }
190185
, coeff2 { FOUR * TWO / THREE / b0 }
186+
, random_pool { 12345 }
191187
, xsurf([&m, this]() {
192188
const auto min_buff = params.template get<unsigned short>("algorithms.current_filters") + 2;
193189
const auto buffer_ncells = min_buff > 5 ? min_buff : 5;
194190
return m.mesh().metric.template convert<1, Crd::Cd, Crd::Ph>(static_cast<real_t>(buffer_ncells));
195191
}())
196-
//, l_atm { ZERO }
197-
// , init_flds(b0, TWO * FOUR * constant::PI * b0 * Omega * SQR(skindepth0) / larmor0, l_atm, ds)
198192
, init_flds(b0, larmor0 / SQR(skindepth0), xsurf, ds)
199193
, ext_current(j0)
194+
, QED_on { p.template get<bool>("setup.QED") }
200195
{}
201196

202197
inline PGen() {}
@@ -230,7 +225,8 @@ namespace user {
230225
params,
231226
local_domain,
232227
injector,
233-
ONE);
228+
ONE,
229+
true);
234230

235231
const auto extra_charge = ExtraCharge<S, M>(
236232
local_domain.mesh.metric,
@@ -246,61 +242,64 @@ namespace user {
246242
params,
247243
local_domain,
248244
injector_extra_charge,
249-
TWO);
245+
TWO,
246+
true);
250247
}
251248

252249
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));
250+
if (QED_on) {
251+
kernel::QED::CurvatureEmission_kernel<D, M::CoordType> curvature_emission1(local_domain.species[0],
252+
local_domain.species[2],
253+
e_min,
254+
gamma_emit,
255+
coeff * dt / skindepth0,
256+
rho,
257+
N_max,
258+
random_pool,
259+
cdf);
260+
Kokkos::parallel_for("CurvatureEmission",
261+
local_domain.species[0].rangeActiveParticles(),
262+
curvature_emission1);
263+
Kokkos::fence();
264+
auto n_injected = curvature_emission1.num_injected();
265+
local_domain.species[2].set_npart(local_domain.species[2].npart() + n_injected);
266+
kernel::QED::CurvatureEmission_kernel<D, M::CoordType> curvature_emission2(local_domain.species[1],
267+
local_domain.species[2],
268+
e_min,
269+
gamma_emit,
270+
coeff * dt / skindepth0,
271+
rho,
272+
N_max,
273+
random_pool,
274+
cdf);
275+
Kokkos::parallel_for("CurvatureEmission",
276+
local_domain.species[1].rangeActiveParticles(),
277+
curvature_emission2);
278+
Kokkos::fence();
279+
n_injected = curvature_emission2.num_injected();
280+
local_domain.species[2].set_npart(local_domain.species[2].npart() + n_injected);
281+
282+
kernel::QED::PairCreation_kernel<D, M::CoordType> pair_creation(local_domain.species[2],
283+
local_domain.species[1],
284+
local_domain.species[0]);
285+
286+
Kokkos::fence();
287+
288+
n_injected = pair_creation.num_injected();
289+
local_domain.species[1].set_npart(local_domain.species[1].npart() + n_injected);
290+
local_domain.species[0].set_npart(local_domain.species[0].npart() + n_injected);
291+
292+
Kokkos::parallel_for("PayloadUpdate",
293+
local_domain.species[2].rangeActiveParticles(),
294+
kernel::QED::PayloadUpdate<D, M::CoordType>(local_domain.species[2],
295+
coeff1 / skindepth0,
296+
coeff2,
297+
dt,
298+
rho,
299+
5));
302300

303301
}
302+
}
304303

305304

306305
}; // PGen

‎src/archetypes/particle_injector.h‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -447,9 +447,9 @@ namespace arch {
447447
raise::ErrorIf((M::CoordType != Coord::Cart) && (not use_weights),
448448
"Weights must be used for non-Cartesian coordinates",
449449
HERE);
450-
raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights,
451-
"Weights should not be used for Cartesian coordinates",
452-
HERE);
450+
// raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights,
451+
// "Weights should not be used for Cartesian coordinates",
452+
// HERE);
453453
raise::ErrorIf(params.template get<bool>("particles.use_weights") != use_weights,
454454
"Weights must be enabled from the input file to use them in "
455455
"the injector",
@@ -558,9 +558,9 @@ namespace arch {
558558
raise::ErrorIf((M::CoordType != Coord::Cart) && (not use_weights),
559559
"Weights must be used for non-Cartesian coordinates",
560560
HERE);
561-
raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights,
562-
"Weights should not be used for Cartesian coordinates",
563-
HERE);
561+
// raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights,
562+
// "Weights should not be used for Cartesian coordinates",
563+
// HERE);
564564
raise::ErrorIf(
565565
params.template get<bool>("particles.use_weights") != use_weights,
566566
"Weights must be enabled from the input file to use them in "
@@ -671,9 +671,9 @@ namespace arch {
671671
raise::ErrorIf((M::CoordType != Coord::Cart) && (not use_weights),
672672
"Weights must be used for non-Cartesian coordinates",
673673
HERE);
674-
raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights,
675-
"Weights should not be used for Cartesian coordinates",
676-
HERE);
674+
// raise::ErrorIf((M::CoordType == Coord::Cart) && use_weights,
675+
// "Weights should not be used for Cartesian coordinates",
676+
// HERE);
677677
raise::ErrorIf(
678678
params.template get<bool>("particles.use_weights") and not use_weights,
679679
"Weights are enabled in the input but not enabled in the injector",

‎src/engines/srpic.hpp‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -330,8 +330,9 @@ namespace ntt {
330330
: ZERO;
331331
const auto has_curvature = (cooling == Cooling::CURVATURE);
332332
const auto cur_coeff = has_curvature
333-
? m_params.template get<real_t>(
334-
"algorithms.curvature.coeff")
333+
? -m_params.template get<real_t>(
334+
"scales.q0") * dt / real_t(6.0) / constant::PI
335+
/ SQR(m_params.template get<real_t>("setup.rho0"))
335336
: ZERO;
336337

337338
// toggle to indicate whether pgen defines the external force
@@ -1196,7 +1197,7 @@ namespace ntt {
11961197
Kokkos::deep_copy(domain.fields.bckp, ZERO);
11971198
auto scatter_bckp = Kokkos::Experimental::create_scatter_view(
11981199
domain.fields.bckp);
1199-
const auto use_weights = M::CoordType != Coord::Cart;
1200+
const auto use_weights = true;
12001201
const auto ni2 = domain.mesh.n_active(in::x2);
12011202
const auto inv_n0 = ONE / m_params.template get<real_t>("scales.n0");
12021203

‎src/global/enums.h‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -265,14 +265,15 @@ namespace ntt {
265265
enum type : uint8_t {
266266
INVALID = 0,
267267
SYNCHROTRON = 1,
268-
NONE = 2,
268+
CURVATURE = 2,
269+
NONE = 3,
269270
};
270271

271272
constexpr Cooling(uint8_t c) : enums_hidden::BaseEnum<Cooling> { c } {}
272273

273-
static constexpr type variants[] = { SYNCHROTRON, NONE };
274-
static constexpr const char* lookup[] = { "synchrotron", "none" };
275-
static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]);
274+
static constexpr type variants[] = { SYNCHROTRON, CURVATURE, NONE };
275+
static constexpr const char* lookup[] = { "synchrotron", "curvature", "none" };
276+
static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]);
276277
};
277278

278279
struct FldsID : public enums_hidden::BaseEnum<FldsID> {

‎src/kernels/QED_process.hpp‎

Lines changed: 6 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -219,7 +219,7 @@ namespace kernel::QED{
219219
, dx1_ph { photons.dx1 }
220220
, e_min { e_min_ }
221221
, gamma_emit { gamma_emit_ }
222-
, coeff { coeff_ }
222+
, coeff { coeff_ * 5.236}
223223
, rho { rho_ }
224224
, N_max { N_max }
225225
, npart_ph { photons.npart() }
@@ -245,7 +245,7 @@ namespace kernel::QED{
245245
}
246246
const real_t pp = math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p)));
247247
const real_t zeta = e_min * rho * CUBE(gamma_emit / pp);
248-
auto N_ph = static_cast<size_t>(coeff * cdf.CDF(zeta) / SQR(pp));
248+
auto N_ph = static_cast<size_t>(coeff * cdf.CDF(zeta) * pp / CUBE(gamma_emit) / rho);
249249

250250
if (N_ph < 1){
251251
return;
@@ -397,13 +397,7 @@ namespace kernel::QED{
397397
, dx1_e { electrons.dx1 }
398398
, tag_e { electrons.tag }
399399
, weight_e { electrons.weight }
400-
, npart_e { electrons.npart() }
401-
, coeff1 { coeff1_ }
402-
, coeff2 { coeff2_ }
403-
, rho0 { rho0_ }
404-
, L { L_ }
405-
, dt { dt_ }
406-
, n_steps { n_steps_ }{
400+
, npart_e { electrons.npart() }{
407401
Kokkos::deep_copy(n_inj, 0);
408402
}
409403
~PairCreation_kernel() = default;
@@ -426,11 +420,10 @@ namespace kernel::QED{
426420
tag_ph(p) = ParticleTag::dead;
427421
return;
428422
}
429-
auto offset_p = Kokkos::atomic_fetch_add(&n_inj(), 1);
430-
auto offset_e = Kokkos::atomic_fetch_add(&n_inj(), 1);
423+
auto offset = Kokkos::atomic_fetch_add(&n_inj(), 1);
431424

432-
offset_p += npart_p;
433-
offset_e += npart_e;
425+
auto offset_p = offset + npart_p;
426+
auto offset_e = offset + npart_e;
434427

435428
if (pld_ph(p, 1) > ONE){
436429
const real_t u = math::sqrt((SQR(pld_ph(p, 0)) - FOUR) / (SQR(pld_ph(p, 0) * pld_ph(p, 2)) + FOUR));

‎src/kernels/particle_pusher_sr.hpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -461,12 +461,12 @@ namespace kernel::sr {
461461

462462
Inline void curvatureDrag(index_t& p,
463463
vec_t<Dim::_3D>& u_prime) const {
464-
real_t gamma_prime_sqr = ONE / (ONE + NORM_SQR(u_prime[0],
465-
u_prime[1],
466-
u_prime[2]));
467-
ux1(p) += coeff_cur * SQR(gamma_prime_sqr) * u_prime[0];
468-
ux2(p) += coeff_cur * SQR(gamma_prime_sqr) * u_prime[1];
469-
ux3(p) += coeff_cur * SQR(gamma_prime_sqr) * u_prime[2];
464+
real_t gamma_prime_sqr = math::sqrt(ONE + NORM_SQR(u_prime[0],
465+
u_prime[1],
466+
u_prime[2]));
467+
ux1(p) += coeff_cur * CUBE(gamma_prime_sqr) * u_prime[0];
468+
ux2(p) += coeff_cur * CUBE(gamma_prime_sqr) * u_prime[1];
469+
ux3(p) += coeff_cur * CUBE(gamma_prime_sqr) * u_prime[2];
470470
}
471471

472472
Inline void operator()(index_t p) const {

0 commit comments

Comments
 (0)