Skip to content

Commit d86bf7d

Browse files
StaticObserverStaticObserver
committed
pgen for polar cap 1d
1 parent 65e6d77 commit d86bf7d

3 files changed

Lines changed: 222 additions & 0 deletions

File tree

‎setups/srpic/1d_polar_cap/1d_polar_cap.toml‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
[simulation]
2+
name = "1d_polar_cap"
3+
engine = "srpic"
4+
runtime = 60.0
5+
6+
[grid]
7+
resolution = 2048
8+
extent = [[0.0, 100.0]]
9+
10+
[grid.metric]
11+
metric = "minkowski"
12+
13+
[grid.boundaries]
14+
fields = [["ATMOSPHERE", "ABSORB"]]
15+
particles = [["ATMOSPHERE", "ABSORB"]]
16+
17+
[grid.boundaries.absorb]
18+
ds = 1.0
19+
20+
[grid.boundaries.atmosphere]
21+
temperature = 0.1
22+
density = 10.0
23+
height = 0.02
24+
species = [1, 2]
25+
ds = 2.0
26+
27+
[scales]
28+
larmor0 = 2e-5
29+
skindepth0 = 0.01
30+
31+
[algorithms]
32+
current_filters = 4
33+
34+
[algorithms.timestep]
35+
CFL = 0.5
36+
37+
[algorithms.gca]
38+
e_ovr_b_max = 0.9
39+
larmor_max = 1.0
40+
41+
[particles]
42+
ppc0 = 5.0
43+
use_weights = true
44+
clear_interval = 100
45+
46+
[[particles.species]]
47+
label = "e-"
48+
mass = 1.0
49+
charge = -1.0
50+
maxnpart = 1e7
51+
pusher = "Boris,GCA"
52+
53+
[[particles.species]]
54+
label = "e+"
55+
mass = 1.0
56+
charge = 1.0
57+
maxnpart = 1e7
58+
pusher = "Boris,GCA"
59+
60+
[setup]
61+
Bsurf = 1.0
62+
period = 1.0
63+
64+
[output]
65+
format = "hdf5"
66+
67+
[output.fields]
68+
interval_time = 1.0
69+
quantities = ["N_1", "N_2", "E", "B", "T00"]
70+
71+
[output.particles]
72+
enable = false
73+
74+
[output.spectra]
75+
enable = true
76+
77+
[diagnostics]
78+
interval = 50
79+
colored_stdout = true

‎setups/srpic/1d_polar_cap/pgen.hpp‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
#ifndef PROBLEM_GENERATOR_H
2+
#define PROBLEM_GENERATOR_H
3+
4+
#include "enums.h"
5+
#include "global.h"
6+
7+
#include "arch/kokkos_aliases.h"
8+
#include "arch/traits.h"
9+
10+
#include "archetypes/problem_generator.h"
11+
#include "framework/domain/metadomain.h"
12+
13+
namespace user {
14+
using namespace ntt;
15+
16+
template <Dimension D>
17+
struct InitFields {
18+
InitFields(real_t b0, real_t angle, real_t rho) : b0 { b0 }, angle { angle }, rho { rho } {}
19+
20+
Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
21+
return b0 * math::cos(angle);
22+
}
23+
24+
Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t {
25+
return ZERO;
26+
}
27+
28+
Inline auto ex1(const coord_t<D>& x_Ph) const -> real_t {
29+
return rho * x_Ph[0];
30+
}
31+
32+
private:
33+
const real_t b0, angle, rho;
34+
};
35+
36+
template <Dimension D>
37+
struct MFields {
38+
MFields(real_t b0, real_t angle) : b0 { b0 }, angle { angle } {}
39+
40+
Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
41+
return b0 * math::cos(angle);
42+
}
43+
44+
Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t {
45+
return b0 * math::sin(angle);
46+
}
47+
48+
private:
49+
const real_t b0, angle;
50+
};
51+
52+
template <Dimension D>
53+
struct DriveFields : public InitFields<D> {
54+
DriveFields(real_t time, real_t b0, real_t angle, real_t omega, real_t r, real_t theta)
55+
: InitFields<D> { b0, angle }
56+
, time { time }
57+
, Omega { omega }
58+
, R { r }
59+
, theta { theta } {}
60+
61+
using InitFields<D>::bx1;
62+
using InitFields<D>::bx2;
63+
64+
Inline auto bx3(const coord_t<D>&) const -> real_t {
65+
return ZERO;
66+
}
67+
68+
Inline auto ex1(const coord_t<D>& x_Ph) const -> real_t {
69+
return Omega * bx2(x_Ph) * R * math::sin(theta);
70+
}
71+
72+
Inline auto ex2(const coord_t<D>& x_Ph) const -> real_t {
73+
return -Omega * bx1(x_Ph) * R * math::sin(theta);
74+
}
75+
76+
Inline auto ex3(const coord_t<D>&) const -> real_t {
77+
return ZERO;
78+
}
79+
80+
private:
81+
const real_t time, Omega, R, theta;
82+
};
83+
84+
template <SimEngine::type S, class M>
85+
struct PGen : public arch::ProblemGenerator<S, M> {
86+
// compatibility traits for the problem generator
87+
static constexpr auto engines { traits::compatible_with<SimEngine::SRPIC>::value };
88+
static constexpr auto metrics {
89+
traits::compatible_with<Metric::Minkowski>::value
90+
};
91+
static constexpr auto dimensions { traits::compatible_with<Dim::_1D>::value };
92+
93+
// for easy access to variables in the child class
94+
using arch::ProblemGenerator<S, M>::D;
95+
using arch::ProblemGenerator<S, M>::C;
96+
using arch::ProblemGenerator<S, M>::params;
97+
98+
const real_t B0, angle, R, Omega, rho, theta, skin0, larmor0;
99+
InitFields<D> init_flds;
100+
101+
inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
102+
: arch::ProblemGenerator<S, M>(p)
103+
, B0 { p.template get<real_t>("setup.B0", ONE) }
104+
, R { p.template get<real_t>("setup.R") }
105+
, Omega { static_cast<real_t>(constant::TWO_PI) /
106+
p.template get<real_t>("setup.period", ONE) }
107+
, angle { p.template get<real_t>("setup.angle", ZERO) }
108+
, theta { p.template get<real_t>("setup.theta", ZERO) }
109+
, skin0 { p.template get<real_t>("scales.skindepth0") }
110+
, larmor0 { p.template get<real_t>("scales.larmor0") }
111+
, rho { TWO * B0 * Omega * SQR(skin0) / larmor0}
112+
, init_flds { B0, angle, rho} {}
113+
114+
inline PGen() {}
115+
116+
auto AtmFields(real_t time) const -> DriveFields<D> {
117+
return DriveFields<D> { time, B0, angle, Omega, R, theta};
118+
}
119+
120+
auto MatchFields(real_t) const -> MFields<D> {
121+
return MFields<D> { B0, angle };
122+
}
123+
};
124+
125+
} // namespace user
126+
127+
#endif

‎src/engines/srpic.hpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,8 @@ namespace ntt {
8181
"algorithms.toggles.deposit");
8282
const auto clear_interval = m_params.template get<std::size_t>(
8383
"particles.clear_interval");
84+
const auto force_free = m_params.template get<bool>(
85+
"setup.force_free");
8486

8587
if (step == 0) {
8688
// communicate fields and apply BCs on the first timestep
@@ -113,6 +115,20 @@ namespace ntt {
113115
timers.start("CurrentDeposit");
114116
Kokkos::deep_copy(dom.fields.cur, ZERO);
115117
CurrentsDeposit(dom);
118+
if constexpr (force_free) {
119+
auto jff = m_params.template get<real_t>("setup.jff");
120+
auto skin0 = m_params.template get<real_t>("scales.skin0");
121+
auto larmor0 = m_params.template get<real_t>("scales.larmor0");
122+
auto B0 = m_params.template get<real_t>("setup.B0");
123+
auto Omega = static_cast<real_t>(constant::TWO_PI) / m_params.template get<real_t>("setup.period", ONE)
124+
jff *= TWO * B0 * Omega * SQR(skin0) / larmor0;
125+
Kokkos::parallel_for("FFCurrentCorrection",
126+
dom.mesh.ExtentToRange(dom.mesh.extent(in::x1), {false, true}),
127+
[](const size_t i){
128+
dom.fields.cur.access(i, cur::jx1) -= jff;
129+
}
130+
)
131+
}
116132
timers.stop("CurrentDeposit");
117133

118134
timers.start("Communications");

0 commit comments

Comments
 (0)