Skip to content

Commit 068f524

Browse files
haykhhaykh
authored
Merge pull request #70 from entity-toolkit/dev/bulk
Add bulk velocity as moments to output
2 parents 330b173 + a22cd96 commit 068f524

File tree

9 files changed

+261
-39
lines changed

9 files changed

+261
-39
lines changed

‎README.md‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,10 +24,14 @@ Our [detailed documentation](https://entity-toolkit.github.io/) includes everyth
2424

2525
## Contributors (alphabetical)
2626

27+
🎸 __Ludwig Böss__ {[@LudwigBoess](https://github.com/LudwigBoess): PIC, framework}
28+
2729
👀 __Yangyang Cai__ {[@StaticObserver](https://github.com/StaticObserver): GRPIC}
2830

2931
🍵 __Benjamin Crinquand__ {[@bcrinquand](https://github.com/bcrinquand): GRPIC, cubed-sphere}
3032

33+
🚂 __Evgeny Gorbunov__ {[@Alcauchy](https://github.com/Alcauchy): PIC, framework}
34+
3135
:radio: __Siddhant Solanki__ {[@sidruns30](https://github.com/sidruns30): framework}
3236

3337
🤷 __Arno Vanthieghem__ {[@vanthieg](https://github.com/vanthieg): framework, PIC}

‎TASKLIST.md‎

Lines changed: 0 additions & 9 deletions
This file was deleted.

‎input.example.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -340,7 +340,7 @@
340340
# Field quantities to output:
341341
# @type: array of strings
342342
# @valid: fields: "E", "B", "J", "divE"
343-
# @valid: moments: "Rho", "Charge", "N", "Nppc", "T0i", "Tij"
343+
# @valid: moments: "Rho", "Charge", "N", "Nppc", "T0i", "Tij", "Vi"
344344
# @valid: for GR: "D", "H", "divD", "A"
345345
# @default: []
346346
# @note: For T, you can use unspecified indices, e.g., Tij, T0i, or specific ones, e.g., Ttt, T00, T02, T23

‎src/framework/domain/output.cpp‎

Lines changed: 63 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -333,6 +333,18 @@ namespace ntt {
333333
{},
334334
local_domain->fields.bckp,
335335
c);
336+
} else if (fld.id() == FldsID::V) {
337+
if constexpr (S != SimEngine::GRPIC) {
338+
ComputeMoments<S, M, FldsID::V>(params,
339+
local_domain->mesh,
340+
local_domain->species,
341+
fld.species,
342+
fld.comp[0],
343+
local_domain->fields.bckp,
344+
c);
345+
} else {
346+
raise::Error("Bulk velocity not supported for GRPIC", HERE);
347+
}
336348
} else {
337349
raise::Error("Wrong moment requested for output", HERE);
338350
}
@@ -362,16 +374,35 @@ namespace ntt {
362374
if (fld.is_moment()) {
363375
for (auto i = 0; i < 3; ++i) {
364376
const auto c = static_cast<unsigned short>(addresses[i]);
365-
raise::ErrorIf(fld.comp[i].size() != 2,
366-
"Wrong # of components requested for moment",
367-
HERE);
368-
ComputeMoments<S, M, FldsID::T>(params,
369-
local_domain->mesh,
370-
local_domain->species,
371-
fld.species,
372-
fld.comp[i],
373-
local_domain->fields.bckp,
374-
c);
377+
if (fld.id() == FldsID::T) {
378+
raise::ErrorIf(fld.comp[i].size() != 2,
379+
"Wrong # of components requested for moment",
380+
HERE);
381+
ComputeMoments<S, M, FldsID::T>(params,
382+
local_domain->mesh,
383+
local_domain->species,
384+
fld.species,
385+
fld.comp[i],
386+
local_domain->fields.bckp,
387+
c);
388+
} else if (fld.id() == FldsID::V) {
389+
raise::ErrorIf(fld.comp[i].size() != 1,
390+
"Wrong # of components requested for 3vel",
391+
HERE);
392+
if constexpr (S == SimEngine::SRPIC) {
393+
ComputeMoments<S, M, FldsID::V>(params,
394+
local_domain->mesh,
395+
local_domain->species,
396+
fld.species,
397+
fld.comp[i],
398+
local_domain->fields.bckp,
399+
c);
400+
} else {
401+
raise::Error("Bulk velocity not supported for GRPIC", HERE);
402+
}
403+
} else {
404+
raise::Error("Wrong moment requested for output", HERE);
405+
}
375406
}
376407
raise::ErrorIf(addresses[1] - addresses[0] !=
377408
addresses[2] - addresses[1],
@@ -380,6 +411,28 @@ namespace ntt {
380411
SynchronizeFields(*local_domain,
381412
Comm::Bckp,
382413
{ addresses[0], addresses[2] + 1 });
414+
if constexpr (S == SimEngine::SRPIC) {
415+
if (fld.id() == FldsID::V) {
416+
// normalize 3vel * rho (combuted above) by rho
417+
ComputeMoments<S, M, FldsID::Rho>(params,
418+
local_domain->mesh,
419+
local_domain->species,
420+
fld.species,
421+
{},
422+
local_domain->fields.bckp,
423+
0u);
424+
SynchronizeFields(*local_domain, Comm::Bckp, { 0, 1 });
425+
Kokkos::parallel_for("NormalizeVectorByRho",
426+
local_domain->mesh.rangeActiveCells(),
427+
kernel::NormalizeVectorByRho_kernel<M::Dim, 6>(
428+
local_domain->fields.bckp,
429+
local_domain->fields.bckp,
430+
0,
431+
addresses[0],
432+
addresses[1],
433+
addresses[2]));
434+
}
435+
}
383436
} else {
384437
// copy fields to bckp (:, 0, 1, 2)
385438
// if as-is specified ==> copy directly to 3, 4, 5

‎src/global/enums.h‎

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@
1313
* - enum ntt::PrtlPusher // boris, vay, photon, none
1414
* - enum ntt::Cooling // synchrotron, none
1515
* - enum ntt::FldsID // e, dive, d, divd, b, h, j,
16-
* a, t, rho, charge, n, nppc, custom
16+
* a, t, rho, charge, n, nppc, v, custom
1717
* @namespaces:
1818
* - ntt::
1919
* @note Enums of the same type can be compared with each other and with strings
@@ -287,17 +287,19 @@ namespace ntt {
287287
Charge = 11,
288288
N = 12,
289289
Nppc = 13,
290-
Custom = 14,
290+
V = 14,
291+
Custom = 15,
291292
};
292293

293294
constexpr FldsID(uint8_t c) : enums_hidden::BaseEnum<FldsID> { c } {}
294295

295-
static constexpr type variants[] = { E, divE, D, divD, B, H, J,
296-
A, T, Rho, Charge, N, Nppc, Custom };
297-
static constexpr const char* lookup[] = { "e", "dive", "d", "divd",
298-
"b", "h", "j", "a",
299-
"t", "rho", "charge", "n",
300-
"nppc", "custom" };
296+
static constexpr type variants[] = { E, divE, D, divD, B,
297+
H, J, A, T, Rho,
298+
Charge, N, Nppc, V, Custom };
299+
static constexpr const char* lookup[] = { "e", "dive", "d", "divd",
300+
"b", "h", "j", "a",
301+
"t", "rho", "charge", "n",
302+
"nppc", "v", "custom" };
301303
static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]);
302304
};
303305

‎src/global/tests/enums.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -66,9 +66,9 @@ auto main() -> int {
6666
enum_str_t all_particle_pushers = { "boris", "vay", "photon", "none" };
6767
enum_str_t all_coolings = { "synchrotron", "none" };
6868

69-
enum_str_t all_out_flds = { "e", "dive", "d", "divd", "b",
70-
"h", "j", "a", "t", "rho",
71-
"charge", "n", "nppc", "custom" };
69+
enum_str_t all_out_flds = { "e", "dive", "d", "divd", "b",
70+
"h", "j", "a", "t", "rho",
71+
"charge", "n", "nppc", "v", "custom" };
7272

7373
checkEnum<Coord>(all_coords);
7474
checkEnum<Metric>(all_metrics);

‎src/kernels/particle_moments.hpp‎

Lines changed: 166 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
#include "global.h"
1515

1616
#include "arch/kokkos_aliases.h"
17+
#include "utils/comparators.h"
1718
#include "utils/error.h"
1819
#include "utils/numeric.h"
1920

@@ -40,8 +41,10 @@ namespace kernel {
4041
static_assert(M::is_metric, "M must be a metric class");
4142
static constexpr auto D = M::Dim;
4243

43-
static_assert((F == FldsID::Rho) || (F == FldsID::Charge) ||
44-
(F == FldsID::N) || (F == FldsID::Nppc) || (F == FldsID::T),
44+
static_assert(!((S == SimEngine::GRPIC) && (F == FldsID::V)),
45+
"Bulk velocity not supported for GRPIC");
46+
static_assert((F == FldsID::Rho) || (F == FldsID::Charge) || (F == FldsID::N) ||
47+
(F == FldsID::Nppc) || (F == FldsID::T) || (F == FldsID::V),
4548
"Invalid field ID");
4649

4750
const unsigned short c1, c2;
@@ -89,8 +92,8 @@ namespace kernel {
8992
std::size_t ni2,
9093
real_t inv_n0,
9194
unsigned short window)
92-
: c1 { (components.size() == 2) ? components[0]
93-
: static_cast<unsigned short>(0) }
95+
: c1 { (components.size() > 0) ? components[0]
96+
: static_cast<unsigned short>(0) }
9497
, c2 { (components.size() == 2) ? components[1]
9598
: static_cast<unsigned short>(0) }
9699
, Buff { scatter_buff }
@@ -200,11 +203,97 @@ namespace kernel {
200203
coeff *= u_Phys[c - 1];
201204
}
202205
}
206+
} else if constexpr (F == FldsID::V) {
207+
real_t gamma { ZERO };
208+
// for bulk 3vel (tetrad basis)
209+
vec_t<Dim::_3D> u_Phys { ZERO };
210+
if constexpr (M::CoordType == Coord::Cart) {
211+
u_Phys[0] = ux1(p);
212+
u_Phys[1] = ux2(p);
213+
u_Phys[2] = ux3(p);
214+
} else {
215+
coord_t<M::PrtlDim> x_Code { ZERO };
216+
x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
217+
x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
218+
if constexpr (D == Dim::_3D) {
219+
x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
220+
} else {
221+
x_Code[2] = phi(p);
222+
}
223+
metric.template transform_xyz<Idx::XYZ, Idx::T>(x_Code,
224+
{ ux1(p), ux2(p), ux3(p) },
225+
u_Phys);
226+
}
227+
if (mass == ZERO) {
228+
gamma = NORM(u_Phys[0], u_Phys[1], u_Phys[2]);
229+
} else {
230+
gamma = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2]));
231+
}
232+
// compute the corresponding moment
233+
coeff = (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / gamma;
203234
} else {
204235
// for other cases, use the `contrib` defined above
205236
coeff = contrib;
206237
}
207238

239+
if constexpr (F == FldsID::V) {
240+
real_t gamma { ZERO };
241+
// for stress-energy tensor
242+
vec_t<Dim::_3D> u_Phys { ZERO };
243+
if constexpr (S == SimEngine::SRPIC) {
244+
// SR
245+
// stress-energy tensor for SR is computed in the tetrad (hatted) basis
246+
if constexpr (M::CoordType == Coord::Cart) {
247+
u_Phys[0] = ux1(p);
248+
u_Phys[1] = ux2(p);
249+
u_Phys[2] = ux3(p);
250+
} else {
251+
static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D");
252+
coord_t<M::PrtlDim> x_Code { ZERO };
253+
x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
254+
x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
255+
if constexpr (D == Dim::_3D) {
256+
x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
257+
} else {
258+
x_Code[2] = phi(p);
259+
}
260+
metric.template transform_xyz<Idx::XYZ, Idx::T>(
261+
x_Code,
262+
{ ux1(p), ux2(p), ux3(p) },
263+
u_Phys);
264+
}
265+
if (mass == ZERO) {
266+
gamma = NORM(u_Phys[0], u_Phys[1], u_Phys[2]);
267+
} else {
268+
gamma = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2]));
269+
}
270+
} else {
271+
// GR
272+
// stress-energy tensor for GR is computed in contravariant basis
273+
static_assert(D != Dim::_1D, "GRPIC 1D");
274+
coord_t<D> x_Code { ZERO };
275+
x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
276+
x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
277+
if constexpr (D == Dim::_3D) {
278+
x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
279+
}
280+
vec_t<Dim::_3D> u_Cntrv { ZERO };
281+
// compute u_i u^i for energy
282+
metric.template transform<Idx::D, Idx::U>(x_Code,
283+
{ ux1(p), ux2(p), ux3(p) },
284+
u_Cntrv);
285+
gamma = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) + u_Cntrv[2] * ux3(p);
286+
if (mass == ZERO) {
287+
gamma = math::sqrt(gamma);
288+
} else {
289+
gamma = math::sqrt(ONE + gamma);
290+
}
291+
metric.template transform<Idx::U, Idx::PU>(x_Code, u_Cntrv, u_Phys);
292+
}
293+
// compute the corresponding moment
294+
coeff = u_Phys[c1 - 1] / gamma;
295+
}
296+
208297
if constexpr (F != FldsID::Nppc) {
209298
// for nppc calculation ...
210299
// ... do not take volume, weights or smoothing into account
@@ -288,6 +377,79 @@ namespace kernel {
288377
}
289378
};
290379

380+
template <Dimension D, unsigned short N>
381+
class NormalizeVectorByRho_kernel {
382+
const ndfield_t<D, N> Rho;
383+
ndfield_t<D, N> Vector;
384+
const unsigned short c_rho, c_v1, c_v2, c_v3;
385+
386+
public:
387+
NormalizeVectorByRho_kernel(const ndfield_t<D, N>& rho,
388+
const ndfield_t<D, N>& vector,
389+
unsigned short crho,
390+
unsigned short cv1,
391+
unsigned short cv2,
392+
unsigned short cv3)
393+
: Rho { rho }
394+
, Vector { vector }
395+
, c_rho { crho }
396+
, c_v1 { cv1 }
397+
, c_v2 { cv2 }
398+
, c_v3 { cv3 } {
399+
raise::ErrorIf(c_rho >= N or c_v1 >= N or c_v2 >= N or c_v3 >= N,
400+
"Invalid component index",
401+
HERE);
402+
raise::ErrorIf(c_rho == c_v1 or c_rho == c_v2 or c_rho == c_v3,
403+
"Invalid component index",
404+
HERE);
405+
raise::ErrorIf(c_v1 == c_v2 or c_v1 == c_v3 or c_v2 == c_v3,
406+
"Invalid component index",
407+
HERE);
408+
}
409+
410+
Inline void operator()(index_t i1) const {
411+
if constexpr (D == Dim::_1D) {
412+
if (not cmp::AlmostZero(Rho(i1, c_rho))) {
413+
Vector(i1, c_v1) /= Rho(i1, c_rho);
414+
Vector(i1, c_v2) /= Rho(i1, c_rho);
415+
Vector(i1, c_v3) /= Rho(i1, c_rho);
416+
}
417+
} else {
418+
raise::KernelError(
419+
HERE,
420+
"1D implementation of NormalizeVectorByRho_kernel called for non-1D");
421+
}
422+
}
423+
424+
Inline void operator()(index_t i1, index_t i2) const {
425+
if constexpr (D == Dim::_2D) {
426+
if (not cmp::AlmostZero(Rho(i1, i2, c_rho))) {
427+
Vector(i1, i2, c_v1) /= Rho(i1, i2, c_rho);
428+
Vector(i1, i2, c_v2) /= Rho(i1, i2, c_rho);
429+
Vector(i1, i2, c_v3) /= Rho(i1, i2, c_rho);
430+
}
431+
} else {
432+
raise::KernelError(
433+
HERE,
434+
"2D implementation of NormalizeVectorByRho_kernel called for non-2D");
435+
}
436+
}
437+
438+
Inline void operator()(index_t i1, index_t i2, index_t i3) const {
439+
if constexpr (D == Dim::_3D) {
440+
if (not cmp::AlmostZero(Rho(i1, i2, i3, c_rho))) {
441+
Vector(i1, i2, i3, c_v1) /= Rho(i1, i2, i3, c_rho);
442+
Vector(i1, i2, i3, c_v2) /= Rho(i1, i2, i3, c_rho);
443+
Vector(i1, i2, i3, c_v3) /= Rho(i1, i2, i3, c_rho);
444+
}
445+
} else {
446+
raise::KernelError(
447+
HERE,
448+
"3D implementation of NormalizeVectorByRho_kernel called for non-3D");
449+
}
450+
}
451+
};
452+
291453
} // namespace kernel
292454

293455
#endif // KERNELS_PARTICLE_MOMENTS_HPP

0 commit comments

Comments
 (0)