Skip to content

Commit 77a06e3

Browse files
StaticObserverStaticObserver
committed
fix errors
1 parent 6eb42be commit 77a06e3

4 files changed

Lines changed: 90 additions & 52 deletions

File tree

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

Lines changed: 39 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -15,67 +15,86 @@ namespace user {
1515

1616
template <Dimension D>
1717
struct InitFields {
18-
InitFields(real_t b0, real_t angle, real_t Omega, real_t d0, real_t rho0)
18+
InitFields(real_t b0, real_t angle, real_t Omega, real_t d0, real_t rho0, real_t r, real_t bp)
1919
: b0 { b0 }
2020
, angle { angle }
2121
, Omega { Omega }
2222
, skindepth0 { d0 }
23-
, larmor0 { rho0 } {}
23+
, larmor0 { rho0 }
24+
, R { r }
25+
, bp { bp } {}
2426

2527
Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
2628
return b0 * math::cos(angle);
2729
}
2830

2931
Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t {
32+
return bp;
33+
}
34+
35+
Inline auto bx3(const coord_t<D>& x_Ph) const -> real_t {
3036
return b0 * math::sin(angle);
3137
}
3238

3339
Inline auto ex1(const coord_t<D>& x_Ph) const -> real_t {
34-
return Omega * b0 * (math::sin(angle) * R + TWO * SQR(skindepth0) * x_Ph[0] / larmor0);
40+
return Omega * b0 * (-math::sin(angle) * R + TWO * SQR(skindepth0) * x_Ph[0] / larmor0);
41+
}
42+
43+
Inline auto ex3(const coord_t<D>& x_Ph) const -> real_t {
44+
return -Omega * bx1(x_Ph) * R;
3545
}
3646

3747
private:
38-
const real_t b0, angle, R, Omega, skindepth0, larmor0;
48+
const real_t b0, angle, R, Omega, skindepth0, larmor0, bp;
3949
};
4050

4151
template <Dimension D>
4252
struct MFields {
43-
MFields(real_t b0, real_t angle) : b0 { b0 }, angle { angle } {}
53+
MFields(real_t b0, real_t angle, real_t omega, real_t r, real_t bp)
54+
: b0 { b0 }
55+
, angle { angle }
56+
, Omega { omega }
57+
, R { r }
58+
, bp { bp } {}
4459

4560
Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
4661
return b0 * math::cos(angle);
4762
}
4863

4964
Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t {
65+
return bp;
66+
}
67+
68+
Inline auto bx3(const coord_t<D>& x_Ph) const -> real_t {
5069
return b0 * math::sin(angle);
5170
}
5271

72+
Inline auto ex3(const coord_t<D>& x_Ph) const -> real_t {
73+
return -Omega * bx1(x_Ph) * R;
74+
}
75+
5376

5477
private:
55-
const real_t b0, angle;
78+
const real_t b0, angle, Omega, R, bp;
5679
};
5780

5881
template <Dimension D>
5982
struct DriveFields : public InitFields<D> {
60-
DriveFields(real_t time, real_t b0, real_t angle, real_t omega, real_t r)
61-
: InitFields<D> { b0, angle}
83+
DriveFields(real_t time, real_t b0, real_t angle, real_t omega, real_t r, real_t d0, real_t rho0, real_t bp)
84+
: InitFields<D> { b0, angle, omega, d0, rho0, r, bp}
6285
, time { time }
6386
, Omega { omega }
6487
, R { r }{}
6588

6689
using InitFields<D>::bx1;
67-
using InitFields<D>::bx2;
68-
69-
Inline auto bx3(const coord_t<D>&) const -> real_t {
70-
return ZERO;
71-
}
90+
using InitFields<D>::bx3;
7291

7392
Inline auto ex1(const coord_t<D>& x_Ph) const -> real_t {
74-
return Omega * bx2(x_Ph) * R;
93+
return -Omega * bx3(x_Ph) * R;
7594
}
7695

7796
Inline auto ex2(const coord_t<D>& x_Ph) const -> real_t {
78-
return -Omega * bx1(x_Ph) * R;
97+
return Omega * bx1(x_Ph) * R;
7998
}
8099

81100
Inline auto ex3(const coord_t<D>&) const -> real_t {
@@ -100,7 +119,7 @@ namespace user {
100119
using arch::ProblemGenerator<S, M>::C;
101120
using arch::ProblemGenerator<S, M>::params;
102121

103-
const real_t B0, angle, R, Omega, skindepth0, larmor0;
122+
const real_t B0, angle, R, Omega, skindepth0, larmor0, bp;
104123
InitFields<D> init_flds;
105124

106125
inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
@@ -112,16 +131,17 @@ namespace user {
112131
, angle { p.template get<real_t>("setup.angle", ZERO) }
113132
, skindepth0 { p.template get<real_t>("scales.skindepth0") }
114133
, larmor0 { p.template get<real_t>("scales.larmor0") }
115-
, init_flds { B0, angle} {}
134+
, bp { p.template get<real_t>("setup.Bp") }
135+
, init_flds { B0, angle, Omega, skindepth0, larmor0, R, bp } {}
116136

117137
inline PGen() {}
118138

119139
auto AtmFields(real_t time) const -> DriveFields<D> {
120-
return DriveFields<D> { time, B0, angle, Omega, R, skindepth0, larmor0 };
140+
return DriveFields<D> { time, B0, angle, Omega, R, skindepth0, larmor0, bp };
121141
}
122142

123143
auto MatchFields(real_t) const -> MFields<D> {
124-
return MFields<D> { B0, angle };
144+
return MFields<D> { B0, angle, Omega, R, bp };
125145
}
126146
};
127147

‎src/engines/srpic.hpp‎

Lines changed: 47 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -89,20 +89,20 @@ namespace ntt {
8989
ParticleInjector(dom);
9090
}
9191

92-
if (fieldsolver_enabled) {
93-
timers.start("FieldSolver");
94-
Faraday(dom, HALF);
95-
timers.stop("FieldSolver");
92+
// if (fieldsolver_enabled) {
93+
// timers.start("FieldSolver");
94+
// Faraday(dom, HALF);
95+
// timers.stop("FieldSolver");
9696

97-
timers.start("Communications");
98-
m_metadomain.CommunicateFields(dom, Comm::B);
99-
timers.stop("Communications");
97+
// timers.start("Communications");
98+
// m_metadomain.CommunicateFields(dom, Comm::B);
99+
// timers.stop("Communications");
100100

101-
timers.start("FieldBoundaries");
102-
FieldBoundaries(dom, BC::B);
103-
timers.stop("FieldBoundaries");
104-
Kokkos::fence();
105-
}
101+
// timers.start("FieldBoundaries");
102+
// FieldBoundaries(dom, BC::B);
103+
// timers.stop("FieldBoundaries");
104+
// Kokkos::fence();
105+
// }
106106

107107
{
108108
timers.start("ParticlePusher");
@@ -131,28 +131,40 @@ namespace ntt {
131131
}
132132

133133
if (fieldsolver_enabled) {
134-
timers.start("FieldSolver");
135-
Faraday(dom, HALF);
136-
timers.stop("FieldSolver");
134+
// timers.start("FieldSolver");
135+
// Faraday(dom, HALF);
136+
// timers.stop("FieldSolver");
137137

138-
timers.start("Communications");
139-
m_metadomain.CommunicateFields(dom, Comm::B);
140-
timers.stop("Communications");
138+
// timers.start("Communications");
139+
// m_metadomain.CommunicateFields(dom, Comm::B);
140+
// timers.stop("Communications");
141141

142-
timers.start("FieldBoundaries");
143-
FieldBoundaries(dom, BC::B);
144-
timers.stop("FieldBoundaries");
142+
// timers.start("FieldBoundaries");
143+
// FieldBoundaries(dom, BC::B);
144+
// timers.stop("FieldBoundaries");
145145

146-
timers.start("FieldSolver");
147-
Ampere(dom, ONE);
148-
timers.stop("FieldSolver");
146+
// timers.start("FieldSolver");
147+
// Ampere(dom, ONE);
148+
// timers.stop("FieldSolver");
149149

150150
if (deposit_enabled) {
151151
timers.start("FieldSolver");
152152
CurrentsAmpere(dom);
153153
timers.stop("FieldSolver");
154154
}
155155

156+
auto& cur_view = dom.fields.cur; // Reference to the view
157+
real_t j_mean = 0.0;
158+
Kokkos::parallel_reduce(
159+
"MeanJ",
160+
dom.mesh.rangeActiveCells(),
161+
KOKKOS_LAMBDA(const std::size_t i, real_t& j_mean_local) {
162+
j_mean_local += cur_view(i, cur::jx1); // Use the captured view
163+
},
164+
j_mean);
165+
166+
std::cout << "Mean Jx1: " << j_mean / dom.mesh.n_active(in::x1) << std::endl;
167+
156168
timers.start("Communications");
157169
m_metadomain.CommunicateFields(dom, Comm::E | Comm::J);
158170
timers.stop("Communications");
@@ -525,20 +537,25 @@ namespace ntt {
525537
const auto q0 = m_params.template get<real_t>("scales.q0");
526538
const auto n0 = m_params.template get<real_t>("scales.n0");
527539
const auto B0 = m_params.template get<real_t>("scales.B0");
540+
const auto skin0 = m_params.template get<real_t>("scales.skindepth0");
541+
const auto larmor0 = m_params.template get<real_t>("scales.larmor0");
542+
const auto b0 = m_params.template get<real_t>("setup.B0");
543+
const auto omega = constant::TWO_PI / m_params.template get<real_t>("setup.period");
544+
const auto rho_ff = -TWO * b0 * omega * SQR(skin0) / larmor0;
528545
const auto coeff = -dt * q0 * n0 / B0;
529546
if constexpr (M::CoordType == Coord::Cart) {
530547
// minkowski case
531548
const auto V0 = m_params.template get<real_t>("scales.V0");
532-
const auto jff = m_params.template get<real_t>("setup.jff");
533-
534-
Kokkos::parallel_for(
549+
auto j0 = TWO * constant::TWO_PI * rho_ff * m_params.template get<real_t>("setup.j0");
550+
j0 /= domain.mesh.metric.dxMin();
551+
Kokkos::parallel_for(
535552
"Ampere",
536553
domain.mesh.rangeActiveCells(),
537554
kernel::mink::CurrentsAmpere_kernel<M::Dim>(domain.fields.em,
538555
domain.fields.cur,
539-
coeff / V0,
540-
ONE / n0,
541-
jff));
556+
coeff ,
557+
ONE / (n0 * V0),
558+
j0));
542559
} else {
543560
auto range = range_with_axis_BCs(domain);
544561
const auto ni2 = domain.mesh.n_active(in::x2);

‎src/kernels/ampere_mink.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -114,8 +114,9 @@ namespace kernel::mink {
114114
J(i1, cur::jx1) *= inv_n0;
115115
J(i1, cur::jx2) *= inv_n0;
116116
J(i1, cur::jx3) *= inv_n0;
117+
117118

118-
E(i1, em::ex1) += (J(i1, cur::jx1) - J0 * inv_n0) * coeff;
119+
E(i1, em::ex1) += (J(i1, cur::jx1) - J0) * coeff;
119120
E(i1, em::ex2) += J(i1, cur::jx2) * coeff;
120121
E(i1, em::ex3) += J(i1, cur::jx3) * coeff;
121122
} else {

‎src/kernels/particle_pusher_sr.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1051,7 +1051,7 @@ namespace kernel::sr {
10511051
Inline void boundaryConditions(index_t& p, coord_t<M::PrtlDim>& xp) const {
10521052
if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) {
10531053
auto invert_vel = false;
1054-
if (i1(p) < 0) {
1054+
if (i1(p) < -1) {
10551055
if (is_periodic_i1min) {
10561056
i1(p) += ni1;
10571057
i1_prev(p) += ni1;
@@ -1062,7 +1062,7 @@ namespace kernel::sr {
10621062
dx1(p) = ONE - dx1(p);
10631063
invert_vel = true;
10641064
}
1065-
} else if (i1(p) >= ni1) {
1065+
} else if (i1(p) >= ni1 + 1) {
10661066
if (is_periodic_i1max) {
10671067
i1(p) -= ni1;
10681068
i1_prev(p) -= ni1;

0 commit comments

Comments
 (0)