Skip to content

Commit 5ba8da8

Browse files
StaticObserverStaticObserver
committed
working version
1 parent 5dbe653 commit 5ba8da8

3 files changed

Lines changed: 135 additions & 43 deletions

File tree

‎pgens/BZjet/accretion.toml‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
[simulation]
2+
name = "BZjet"
3+
engine = "grpic"
4+
runtime = 500.0
5+
6+
[grid]
7+
resolution = [256, 256]
8+
extent = [[0.98, 10.0]]
9+
10+
[grid.metric]
11+
metric = "qkerr_schild"
12+
qsph_r0 = 0.0
13+
qsph_h = 0.0
14+
ks_a = 0.99
15+
16+
[grid.boundaries]
17+
fields = [["MATCH"]]
18+
particles = [["ABSORB"]]
19+
20+
[grid.boundaries.absorb]
21+
ds = 1.0
22+
23+
[scales]
24+
larmor0 = 0.025
25+
skindepth0 = 0.5
26+
27+
[algorithms]
28+
current_filters = 4
29+
30+
[algorithms.gr]
31+
pusher_niter = 10
32+
pusher_eps = 1e-2
33+
34+
[algorithms.timestep]
35+
CFL = 0.5
36+
correction = 1.0
37+
38+
[algorithms.toggles]
39+
deposit = true
40+
fieldsolver = true
41+
42+
[particles]
43+
ppc0 = 4.0
44+
use_weights = true
45+
clear_interval = 100
46+
47+
[[particles.species]]
48+
label = "e-"
49+
mass = 1.0
50+
charge = -1.0
51+
maxnpart = 2e8
52+
pusher = "Boris"
53+
54+
[[particles.species]]
55+
label = "e+"
56+
mass = 1.0
57+
charge = 1.0
58+
maxnpart = 2e8
59+
pusher = "Boris"
60+
61+
[setup]
62+
sigma_max = 1000.0
63+
temperature = 0.01
64+
xi_min = [1.1, 0.0]
65+
xi_max = [6.0, 3.14159265]
66+
m_eps = 1.0
67+
inj_coeff = 0.5
68+
db_thr = 1e-2
69+
70+
[output]
71+
format = "hdf5"
72+
73+
[output.fields]
74+
interval_time = 1.0
75+
quantities = ["D", "B", "N_1", "N_2", "A"]
76+
77+
[output.particles]
78+
enable = false
79+
80+
[output.spectra]
81+
enable = false
82+
83+
[diagnostics]
84+
interval = 2
85+
colored_stdout = true
86+
blocking_timers = true

‎pgens/BZjet/pgen.hpp‎

Lines changed: 45 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ namespace user {
2424

2525
template <class M, Dimension D>
2626
struct InitFields {
27-
InitFields(M metric_) : metric { metric_ } {}
27+
InitFields(M metric_, real_t m_eps_) : metric { metric_ }, m_eps { m_eps_ } {}
2828

2929
Inline auto A_3(const coord_t<D>& x_Cd) const -> real_t {
3030
return HALF * (metric.template h_<3, 3>(x_Cd) +
@@ -51,33 +51,33 @@ namespace user {
5151
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
5252

5353
x0m[0] = xi[0];
54-
x0m[1] = xi[1] - HALF;
54+
x0m[1] = xi[1] - HALF * m_eps;
5555
x0p[0] = xi[0];
56-
x0p[1] = xi[1] + HALF;
56+
x0p[1] = xi[1] + HALF * m_eps;
5757

5858
real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
5959

6060
if (cmp::AlmostZero(x_Ph[1])) {
6161
return ONE;
6262
} else {
63-
return (A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP;
63+
return (A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP / m_eps;
6464
}
6565
}
6666

6767
Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t { // at ( i + HALF , j )
6868
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
6969
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
7070

71-
x0m[0] = xi[0] - HALF;
71+
x0m[0] = xi[0] - HALF * m_eps;
7272
x0m[1] = xi[1];
73-
x0p[0] = xi[0] + HALF;
73+
x0p[0] = xi[0] + HALF * m_eps;
7474
x0p[1] = xi[1];
7575

7676
real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
7777
if (cmp::AlmostZero(x_Ph[1])) {
7878
return ZERO;
7979
} else {
80-
return -(A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP;
80+
return -(A_3(x0p) - A_3(x0m)) * inv_sqrt_detH_ijP / m_eps;
8181
}
8282
}
8383

@@ -87,12 +87,12 @@ namespace user {
8787
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
8888

8989
x0m[0] = xi[0];
90-
x0m[1] = xi[1] - HALF;
90+
x0m[1] = xi[1] - HALF * m_eps;
9191
x0p[0] = xi[0];
92-
x0p[1] = xi[1] + HALF;
92+
x0p[1] = xi[1] + HALF * m_eps;
9393

9494
real_t inv_sqrt_detH_iPjP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
95-
return -(A_1(x0p) - A_1(x0m)) * inv_sqrt_detH_iPjP;
95+
return -(A_1(x0p) - A_1(x0m)) * inv_sqrt_detH_iPjP / m_eps;
9696
}
9797

9898
Inline auto dx1(const coord_t<D>& x_Ph) const -> real_t { // at ( i + HALF , j )
@@ -101,24 +101,24 @@ namespace user {
101101

102102
real_t alpha_iPj { metric.alpha({ xi[0], xi[1] }) };
103103
real_t inv_sqrt_detH_ij { ONE / metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
104-
real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
104+
real_t sqrt_detH_ij { metric.sqrt_det_h({ xi[0] - HALF, xi[1] }) };
105105
real_t beta_ij { metric.beta1({ xi[0] - HALF, xi[1] }) };
106106
real_t alpha_ij { metric.alpha({ xi[0] - HALF, xi[1] }) };
107107

108108
// D1 at ( i + HALF , j )
109-
x0m[0] = xi[0] - HALF;
109+
x0m[0] = xi[0] - HALF * m_eps;
110110
x0m[1] = xi[1];
111-
x0p[0] = xi[0] + HALF;
111+
x0p[0] = xi[0] + HALF * m_eps;
112112
x0p[1] = xi[1];
113-
real_t E1d { (A_0(x0p) - A_0(x0m)) };
113+
real_t E1d { (A_0(x0p) - A_0(x0m)) / m_eps };
114114
real_t D1d { E1d / alpha_iPj };
115115

116116
// D3 at ( i , j )
117-
x0m[0] = xi[0] - HALF - HALF;
117+
x0m[0] = xi[0] - HALF - HALF * m_eps;
118118
x0m[1] = xi[1];
119-
x0p[0] = xi[0] - HALF + HALF;
119+
x0p[0] = xi[0] - HALF + HALF * m_eps;
120120
x0p[1] = xi[1];
121-
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
121+
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij / m_eps};
122122

123123
real_t D1u { metric.template h<1, 1>({ xi[0], xi[1] }) * D1d +
124124
metric.template h<1, 3>({ xi[0], xi[1] }) * D3d };
@@ -130,16 +130,16 @@ namespace user {
130130
coord_t<D> xi { ZERO }, x0m { ZERO }, x0p { ZERO };
131131
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
132132
x0m[0] = xi[0];
133-
x0m[1] = xi[1] - HALF;
133+
x0m[1] = xi[1] - HALF * m_eps;
134134
x0p[0] = xi[0];
135-
x0p[1] = xi[1] + HALF;
135+
x0p[1] = xi[1] + HALF * m_eps;
136136
real_t inv_sqrt_detH_ijP { ONE / metric.sqrt_det_h({ xi[0], xi[1] }) };
137137
real_t sqrt_detH_ijP { metric.sqrt_det_h({ xi[0], xi[1] }) };
138138
real_t alpha_ijP { metric.alpha({ xi[0], xi[1] }) };
139139
real_t beta_ijP { metric.beta1({ xi[0], xi[1] }) };
140140

141-
real_t E2d { (A_0(x0p) - A_0(x0m)) };
142-
real_t D2d { E2d / alpha_ijP - (A_1(x0p) - A_1(x0m)) * beta_ijP / alpha_ijP };
141+
real_t E2d { (A_0(x0p) - A_0(x0m)) / m_eps };
142+
real_t D2d { E2d / alpha_ijP - (A_1(x0p) - A_1(x0m)) * beta_ijP / alpha_ijP / m_eps};
143143
real_t D2u { metric.template h<2, 2>({ xi[0], xi[1] }) * D2d };
144144

145145
return D2u;
@@ -155,18 +155,18 @@ namespace user {
155155
real_t alpha_iPj { metric.alpha({ xi[0] + HALF, xi[1] }) };
156156

157157
// D3 at ( i , j )
158-
x0m[0] = xi[0] - HALF;
158+
x0m[0] = xi[0] - HALF * m_eps;
159159
x0m[1] = xi[1];
160-
x0p[0] = xi[0] + HALF;
160+
x0p[0] = xi[0] + HALF * m_eps;
161161
x0p[1] = xi[1];
162-
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij };
162+
real_t D3d { (A_3(x0p) - A_3(x0m)) * beta_ij / alpha_ij / m_eps};
163163

164164
// D1 at ( i + HALF , j )
165-
x0m[0] = xi[0] + HALF - HALF;
165+
x0m[0] = xi[0] + HALF - HALF * m_eps;
166166
x0m[1] = xi[1];
167-
x0p[0] = xi[0] + HALF + HALF;
167+
x0p[0] = xi[0] + HALF + HALF * m_eps;
168168
x0p[1] = xi[1];
169-
real_t E1d { (A_0(x0p) - A_0(x0m)) };
169+
real_t E1d { (A_0(x0p) - A_0(x0m)) / m_eps };
170170
real_t D1d { E1d / alpha_iPj };
171171

172172
if (cmp::AlmostZero(x_Ph[1])) {
@@ -179,17 +179,16 @@ namespace user {
179179

180180
private:
181181
const M metric;
182+
const real_t m_eps;
182183
};
183184

184185
template <SimEngine::type S, class M>
185186
struct PointDistribution : public arch::SpatialDistribution<S, M> {
186187
PointDistribution(const std::vector<real_t>& xi_min,
187188
const std::vector<real_t>& xi_max,
188-
const real_t d0,
189-
const real_t rho0,
190189
const real_t sigma_thr,
190+
const real_t inj_coeff,
191191
const real_t db_thr,
192-
const real_t dens_thr,
193192
const SimulationParams& params,
194193
Domain<S, M>* domain_ptr)
195194
: arch::SpatialDistribution<S, M> { domain_ptr->mesh.metric }
@@ -201,7 +200,7 @@ namespace user {
201200
, d0 { params.template get<real_t>("scales.skindepth0") }
202201
, rho0 { params.template get<real_t>("scales.larmor0") }
203202
, inv_n0 { ONE / params.template get<real_t>("scales.n0") }
204-
, dens_thr { dens_thr } {
203+
, inj_coeff { inj_coeff } {
205204
std::copy(xi_min.begin(), xi_min.end(), x_min);
206205
std::copy(xi_max.begin(), xi_max.end(), x_max);
207206

@@ -258,16 +257,21 @@ namespace user {
258257
DOT(B_cntrv[0], B_cntrv[1], B_cntrv[2], B_cov[0], B_cov[1], B_cov[2]);
259258
const auto db = DOT(D_cntrv[0], D_cntrv[1], D_cntrv[2], B_cov[0], B_cov[1], B_cov[2]);
260259
const auto dens = density(i1, i2, 0);
261-
return (bsqr > sigma_thr * dens) && (db > db_thr * bsqr);
260+
if(cmp::AlmostZero(dens)) {
261+
return db * SIGN(db) > db_thr * bsqr;
262+
}else{
263+
return (bsqr > sigma_thr * dens) && (db * SIGN(db) > db_thr * bsqr);
264+
}
262265
}
263266
return false;
264267
}
265268

266-
Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t override {
267-
auto fill = false;
269+
Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t {
270+
auto fill = true;
268271
for (auto d = 0u; d < M::Dim; ++d) {
269272
fill &= x_Ph[d] > x_min[d] and x_Ph[d] < x_max[d] and sigma_crit(x_Ph);
270273
}
274+
coord_t<M::Dim> xi { ZERO };
271275
metric.template convert<Crd::Ph, Crd::Cd>(x_Ph, xi);
272276
const auto i1 = static_cast<int>(xi[0]) + static_cast<int>(N_GHOSTS);
273277
const auto i2 = static_cast<int>(xi[1]) + static_cast<int>(N_GHOSTS);
@@ -282,7 +286,7 @@ namespace user {
282286
const auto bsqr =
283287
DOT(B_cntrv[0], B_cntrv[1], B_cntrv[2], B_cov[0], B_cov[1], B_cov[2]);
284288
const auto db = DOT(D_cntrv[0], D_cntrv[1], D_cntrv[2], B_cov[0], B_cov[1], B_cov[2]);
285-
const real_t inj_n = dens_thr * db / SQRT(bsqr) * SQR(d0) / rho0;
289+
const real_t inj_n = inj_coeff * db * SIGN(db) / math::sqrt(bsqr) * SQR(d0) / rho0;
286290
return fill ? inj_n : ZERO;
287291
}
288292

@@ -291,7 +295,7 @@ namespace user {
291295
tuple_t<real_t, M::Dim> x_max;
292296
const real_t sigma_thr;
293297
const real_t db_thr;
294-
const real_t dens_thr;
298+
const real_t inj_coeff;
295299
const real_t inv_n0;
296300
const real_t d0;
297301
const real_t rho0;
@@ -318,7 +322,7 @@ namespace user {
318322

319323
const std::vector<real_t> xi_min;
320324
const std::vector<real_t> xi_max;
321-
const real_t sigma0, sigma_max, multiplicity, temperature, m_eps;
325+
const real_t sigma0, sigma_max, inj_coeff, db_thr, temperature, m_eps;
322326

323327
InitFields<M, D> init_flds;
324328
const Metadomain<S, M>* metadomain;
@@ -329,7 +333,8 @@ namespace user {
329333
, xi_max { p.template get<std::vector<real_t>>("setup.xi_max") }
330334
, sigma_max { p.template get<real_t>("setup.sigma_max") }
331335
, sigma0 { p.template get<real_t>("scales.sigma0") }
332-
, multiplicity { p.template get<real_t>("setup.multiplicity") }
336+
, inj_coeff { p.template get<real_t>("setup.inj_coeff") }
337+
, db_thr { p.template get<real_t>("setup.db_thr") }
333338
, temperature { p.template get<real_t>("setup.temperature") }
334339
, m_eps { p.template get<real_t>("setup.m_eps") }
335340
, init_flds { m.mesh().metric, m_eps }
@@ -342,7 +347,8 @@ namespace user {
342347
const auto spatial_dist = PointDistribution<S, M>(xi_min,
343348
xi_max,
344349
sigma_max / sigma0,
345-
multiplicity,
350+
inj_coeff,
351+
db_thr,
346352
params,
347353
&local_domain);
348354

‎pgens/accretion/accretion.toml‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111
metric = "qkerr_schild"
1212
qsph_r0 = 0.0
1313
qsph_h = 0.0
14-
ks_a = 0.95
14+
ks_a = 0.99
1515

1616
[grid.boundaries]
1717
fields = [["MATCH"]]
@@ -60,10 +60,10 @@
6060

6161
[setup]
6262
multiplicity = 1.0
63-
sigma_max = 1000.0
63+
sigma_max = 2000.0
6464
temperature = 0.01
65-
xi_min = [1.5, 0.0]
66-
xi_max = [4.0, 3.14159265]
65+
xi_min = [1.1, 0.0]
66+
xi_max = [6.0, 3.14159265]
6767
m_eps = 1.0
6868

6969
[output]

0 commit comments

Comments
 (0)