Skip to content

Commit 8ff0030

Browse files
haykhhaykh
authored
Merge pull request #136 from entity-toolkit/dev/uniform-maxw-injector
v1.2.3 Release Candidate
2 parents 63ce560 + 07fd04c commit 8ff0030

20 files changed

Lines changed: 757 additions & 116 deletions

File tree

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@ venv/
3939
*.png
4040
*.mov
4141
*.mp4
42+
!pgens/**/*.png
4243

4344
# Accidental files
4445
*.xc

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ set(PROJECT_NAME entity)
77

88
project(
99
${PROJECT_NAME}
10-
VERSION 1.2.1
10+
VERSION 1.2.3
1111
LANGUAGES CXX C)
1212
add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"")
1313
set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ")

‎input.example.toml‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -467,8 +467,9 @@
467467
# Field quantities to output
468468
# @type: array<string>
469469
# @default: ["B^2", "E^2", "ExB", "Rho", "T00"]
470-
# @enum: "B^2", "E^2", "ExB", "N", "Charge", "Rho", "T00", "T0i", "Tij"
471-
# @note: Same notation as for `output.fields.quantities`
470+
# @enum: "B^2", "E^2", "ExB", "N", "Npart", "Charge", "Rho", "T00", "T0i", "Tij"
471+
# @note: For particle moments, ...
472+
# @note: ... same notation is used as for `output.fields.quantities`
472473
quantities = ""
473474
# Custom (user-defined) stats
474475
# @type: array<string>

‎pgens/magnetosphere/magnetosphere.toml‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,9 @@
5858
pusher = "Boris,GCA"
5959

6060
[setup]
61-
Bsurf = 1.0
62-
period = 60.0
61+
Bsurf = 1.0
62+
field_geometry = "dipole" # can be "dipole" or "monopole" (default: dipole)
63+
period = 60.0
6364

6465
[output]
6566
format = "hdf5"

‎pgens/magnetosphere/pgen.hpp‎

Lines changed: 38 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -6,33 +6,58 @@
66

77
#include "arch/kokkos_aliases.h"
88
#include "arch/traits.h"
9+
#include "utils/numeric.h"
910

1011
#include "archetypes/problem_generator.h"
1112
#include "framework/domain/metadomain.h"
1213

14+
#include <string>
15+
1316
namespace user {
1417
using namespace ntt;
1518

19+
enum class FieldGeometry {
20+
dipole,
21+
monopole
22+
};
23+
1624
template <Dimension D>
1725
struct InitFields {
18-
InitFields(real_t bsurf, real_t rstar) : Bsurf { bsurf }, Rstar { rstar } {}
26+
InitFields(real_t bsurf, real_t rstar, const std::string& field_geometry)
27+
: Bsurf { bsurf }
28+
, Rstar { rstar }
29+
, field_geom { field_geometry == "monopole" ? FieldGeometry::monopole
30+
: FieldGeometry::dipole } {}
1931

2032
Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
21-
return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
33+
if (field_geom == FieldGeometry::monopole) {
34+
return Bsurf / SQR(x_Ph[0] / Rstar);
35+
} else {
36+
return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
37+
}
2238
}
2339

2440
Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t {
25-
return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
41+
if (field_geom == FieldGeometry::monopole) {
42+
return ZERO;
43+
} else {
44+
return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
45+
}
2646
}
2747

2848
private:
29-
const real_t Bsurf, Rstar;
49+
const real_t Bsurf, Rstar;
50+
const FieldGeometry field_geom;
3051
};
3152

3253
template <Dimension D>
3354
struct DriveFields : public InitFields<D> {
34-
DriveFields(real_t time, real_t bsurf, real_t rstar, real_t omega)
35-
: InitFields<D> { bsurf, rstar }
55+
DriveFields(real_t time,
56+
real_t bsurf,
57+
real_t rstar,
58+
real_t omega,
59+
const std::string& field_geometry)
60+
: InitFields<D> { bsurf, rstar, field_geometry }
3661
, time { time }
3762
, Omega { omega } {}
3863

@@ -73,25 +98,27 @@ namespace user {
7398
using arch::ProblemGenerator<S, M>::C;
7499
using arch::ProblemGenerator<S, M>::params;
75100

76-
const real_t Bsurf, Rstar, Omega;
77-
InitFields<D> init_flds;
101+
const real_t Bsurf, Rstar, Omega;
102+
const std::string field_geom;
103+
InitFields<D> init_flds;
78104

79105
inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
80106
: arch::ProblemGenerator<S, M>(p)
81107
, Bsurf { p.template get<real_t>("setup.Bsurf", ONE) }
82108
, Rstar { m.mesh().extent(in::x1).first }
83109
, Omega { static_cast<real_t>(constant::TWO_PI) /
84110
p.template get<real_t>("setup.period", ONE) }
85-
, init_flds { Bsurf, Rstar } {}
111+
, field_geom { p.template get<std::string>("setup.field_geometry", "dipole") }
112+
, init_flds { Bsurf, Rstar, field_geom } {}
86113

87114
inline PGen() {}
88115

89116
auto AtmFields(real_t time) const -> DriveFields<D> {
90-
return DriveFields<D> { time, Bsurf, Rstar, Omega };
117+
return DriveFields<D> { time, Bsurf, Rstar, Omega, field_geom };
91118
}
92119

93120
auto MatchFields(real_t) const -> InitFields<D> {
94-
return InitFields<D> { Bsurf, Rstar };
121+
return InitFields<D> { Bsurf, Rstar, field_geom };
95122
}
96123
};
97124

‎pgens/magnetosphere/sketch.png‎

190 KB
Loading

‎pgens/magnetosphere/sketch.py‎

Lines changed: 163 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,163 @@
1+
import matplotlib.pyplot as plt
2+
import matplotlib.patches as mpatches
3+
import myplotlib
4+
import numpy as np
5+
6+
plt.style.use("latex")
7+
8+
fig = plt.figure(dpi=600)
9+
ax = fig.add_subplot(111)
10+
11+
r0 = 0.15
12+
13+
ax.plot([0, 0], [r0, 1], c="k", lw=0.5)
14+
ax.plot([0, 0], [-r0, -1], c="k", lw=0.5)
15+
ax.add_patch(
16+
mpatches.Arc((0, 0), 2, 2, angle=0, theta1=-90, theta2=90, color="k", lw=0.5)
17+
)
18+
ax.add_patch(
19+
mpatches.Arc(
20+
(0, 0), 2 * r0, 2 * r0, angle=0, theta1=-90, theta2=90, color="k", lw=0.5
21+
)
22+
)
23+
ax.add_patch(
24+
mpatches.Arc(
25+
(0, 0),
26+
2 - 0.25,
27+
2 - 0.25,
28+
angle=0,
29+
theta1=-90,
30+
theta2=90,
31+
color="b",
32+
lw=0.5,
33+
ls="--",
34+
)
35+
)
36+
ax.add_patch(
37+
mpatches.Arc(
38+
(0, 0),
39+
2 - 0.1,
40+
2 - 0.1,
41+
angle=0,
42+
theta1=-90,
43+
theta2=90,
44+
color="r",
45+
lw=0.5,
46+
ls="--",
47+
)
48+
)
49+
ax.text(
50+
1.05 / np.sqrt(2), 1.05 / np.sqrt(2), "absorbing particle boundaries", c="r", size=5
51+
)
52+
ax.text(
53+
1.05 / np.sqrt(2),
54+
1.05 / np.sqrt(2) - 0.04,
55+
r"$\mathtt{grid.boundaries.absorb.ds}$",
56+
c="r",
57+
size=5,
58+
)
59+
ax.annotate(
60+
"",
61+
xy=(1.02 / np.sqrt(2), 1.02 / np.sqrt(2)),
62+
xytext=(0.93 / np.sqrt(2), 0.93 / np.sqrt(2)),
63+
arrowprops=dict(arrowstyle="<->", lw=0.5, color="r"),
64+
size=5,
65+
)
66+
ax.annotate(
67+
"",
68+
xy=(1.01 * np.cos(np.pi / 6), 1.01 * np.sin(np.pi / 6)),
69+
xytext=(0.87 * np.cos(np.pi / 6), 0.87 * np.sin(np.pi / 6)),
70+
arrowprops=dict(arrowstyle="<->", lw=0.5, color="b"),
71+
size=5,
72+
)
73+
ax.text(
74+
1.05 * np.cos(np.pi / 6),
75+
1.05 * np.sin(np.pi / 6),
76+
"matching field boundaries",
77+
c="b",
78+
size=5,
79+
)
80+
ax.text(
81+
1.05 * np.cos(np.pi / 6),
82+
1.05 * np.sin(np.pi / 6) - 0.04,
83+
r"$\mathtt{grid.boundaries.match.ds}$",
84+
c="b",
85+
size=5,
86+
)
87+
ax.add_patch(
88+
mpatches.Arc(
89+
(0, 0),
90+
0.4,
91+
0.4,
92+
angle=0,
93+
theta1=-90,
94+
theta2=90,
95+
color="g",
96+
lw=0.5,
97+
ls=":",
98+
)
99+
)
100+
ax.add_patch(
101+
mpatches.Arc(
102+
(0, 0),
103+
0.5,
104+
0.5,
105+
angle=0,
106+
theta1=-90,
107+
theta2=90,
108+
color="g",
109+
lw=0.5,
110+
ls=":",
111+
)
112+
)
113+
ax.annotate(
114+
"",
115+
xy=(0.27 / np.sqrt(2), 0.27 / np.sqrt(2)),
116+
xytext=(0.18 / np.sqrt(2), 0.18 / np.sqrt(2)),
117+
arrowprops=dict(arrowstyle="<->", lw=0.5, color="g"),
118+
size=5,
119+
)
120+
ax.text(
121+
0.27 / np.sqrt(2),
122+
0.27 / np.sqrt(2) + 0.04,
123+
"particle atmosphere injection",
124+
c="g",
125+
size=5,
126+
)
127+
ax.text(
128+
0.27 / np.sqrt(2),
129+
0.27 / np.sqrt(2),
130+
r"$\mathtt{grid.boundaries.atmosphere.height}$",
131+
c="g",
132+
size=5,
133+
)
134+
ax.annotate(
135+
"",
136+
xy=(0.13 / np.sqrt(2), -0.13 / np.sqrt(2)),
137+
xytext=(0.22 / np.sqrt(2), -0.22 / np.sqrt(2)),
138+
arrowprops=dict(arrowstyle="<->", lw=0.5, color="b"),
139+
size=5,
140+
)
141+
ax.text(
142+
0.23 / np.sqrt(2),
143+
-0.23 / np.sqrt(2),
144+
"buffer zone for resetting fields",
145+
c="b",
146+
size=5,
147+
)
148+
ax.text(
149+
0.23 / np.sqrt(2),
150+
-0.23 / np.sqrt(2)-0.04,
151+
r"size in cells = \# of filters",
152+
c="b",
153+
size=5,
154+
)
155+
ax.set(
156+
xlim=(-0.05, 1.05),
157+
ylim=(-1.05, 1.05),
158+
aspect=1,
159+
xticks=[],
160+
yticks=[],
161+
frame_on=False,
162+
)
163+
plt.savefig("sketch.png", bbox_inches="tight")

‎pgens/reconnection/pgen.hpp‎

Lines changed: 6 additions & 11 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 "archetypes/spatial_dist.h"
16+
#include "archetypes/utils.h"
1617
#include "framework/domain/metadomain.h"
1718

1819
#include "kernels/particle_moments.hpp"
@@ -205,17 +206,11 @@ namespace user {
205206

206207
inline void InitPrtls(Domain<S, M>& local_domain) {
207208
// background
208-
const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
209-
local_domain.random_pool,
210-
bg_temperature);
211-
const auto injector = arch::UniformInjector<S, M, arch::Maxwellian>(
212-
energy_dist,
213-
{ 1, 2 });
214-
arch::InjectUniform<S, M, arch::UniformInjector<S, M, arch::Maxwellian>>(
215-
params,
216-
local_domain,
217-
injector,
218-
ONE);
209+
arch::InjectUniformMaxwellian<S, M>(params,
210+
local_domain,
211+
ONE,
212+
bg_temperature,
213+
{ 1, 2 });
219214

220215
const auto sigma = params.template get<real_t>("scales.sigma0");
221216
const auto c_omp = params.template get<real_t>("scales.skindepth0");

‎pgens/reconnection/reconnection.toml‎

Lines changed: 11 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
[simulation]
22
name = "reconnection"
33
engine = "srpic"
4-
runtime = 10.0
4+
runtime = 10000.0
55

66
[simulation.domain]
77
decomposition = [-1, 2]
88

99
[grid]
10-
resolution = [512, 512]
11-
extent = [[-1.0, 1.0], [-1.0, 1.0]]
10+
resolution = [4096, 2048]
11+
extent = [[-500.0, 500.0], [-250.0, 250.0]]
1212

1313
[grid.metric]
1414
metric = "minkowski"
@@ -18,11 +18,11 @@
1818
particles = [["PERIODIC"], ["ABSORB", "ABSORB"]]
1919

2020
[grid.boundaries.match]
21-
ds = [[0.04], [0.1]]
21+
ds = [[10.0], [20.0]]
2222

2323
[scales]
24-
larmor0 = 2e-4
25-
skindepth0 = 2e-3
24+
larmor0 = 0.1
25+
skindepth0 = 1.0
2626

2727
[algorithms]
2828
current_filters = 8
@@ -49,23 +49,19 @@
4949
bg_B = 1.0
5050
bg_Bguide = 0.0
5151
bg_temperature = 1e-4
52-
inj_ypad = 0.25
53-
cs_width = 0.05
52+
inj_ypad = 50.0
53+
cs_width = 10.0
5454
cs_overdensity = 3.0
5555

5656
[output]
57-
format = "hdf5"
58-
interval_time = 0.1
57+
format = "bpfile"
58+
interval_time = 100.0
5959

6060
[output.fields]
6161
quantities = ["N_1", "N_2", "E", "B", "J"]
6262

6363
[output.particles]
64-
enable = false
64+
stride = 25
6565

6666
[output.spectra]
6767
enable = false
68-
69-
[diagnostics]
70-
colored_stdout = true
71-
interval = 10

‎pgens/reconnection/sketch.png‎

109 KB
Loading

0 commit comments

Comments
 (0)