45template<
typename REAL>
inline
47 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
50 const std::array<Real, 3>& gridSpacing
53 const auto dx = gridSpacing[0];
54 const auto dy = gridSpacing[1];
55 const auto dz = gridSpacing[2];
56 return -(pC[
a_zz] * BGBZ) / dz + (pC[
a_z] * BGBZ) / dz - (pC[
a_yz] * BGBZ) / (2 * dz) -
58 (pC[
c_xy] * BGBZ) / (2 *
dx) - (pC[
c_x] * BGBZ) /
dx - (pC[
a_yz] * BGBY) / (2 * dy) - (pC[
a_yy] * BGBY) / dy +
59 (pC[
a_y] * BGBY) / dy + (pC[
b_xz] * BGBY) / (2 *
dx) - (pC[
b_xyz] * BGBY) / (4 *
dx) -
62 (pC[
a_zz] * pC[
c_z]) / (2 * dz) - (pC[
a_z] * pC[
c_z]) / (2 * dz) + (pC[
a_yz] * pC[
c_z]) / (4 * dz) -
64 (pC[
a_zz] * pC[
c_y]) / (2 * dz) - (pC[
a_z] * pC[
c_y]) / (2 * dz) + (pC[
a_yz] * pC[
c_y]) / (4 * dz) +
68 (pC[
a_yz] * pC[
b_z]) / (4 * dy) + (pC[
a_yy] * pC[
b_z]) / (2 * dy) - (pC[
a_y] * pC[
b_z]) / (2 * dy) -
71 (pC[
a_yz] * pC[
b_y]) / (4 * dy) + (pC[
a_yy] * pC[
b_y]) / (2 * dy) - (pC[
a_y] * pC[
b_y]) / (2 * dy) +
109template<
typename REAL>
inline
111 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
114 const std::array<Real, 3>& gridSpacing
117 const auto dx = gridSpacing[0];
118 const auto dy = gridSpacing[1];
119 const auto dz = gridSpacing[2];
120 return -(pC[
a_zz] * BGBZ) / dz + (pC[
a_z] * BGBZ) / dz + (pC[
a_yz] * BGBZ) / (2 * dz) -
122 (pC[
c_xy] * BGBZ) / (2 *
dx) - (pC[
c_x] * BGBZ) /
dx - (pC[
a_yz] * BGBY) / (2 * dy) + (pC[
a_yy] * BGBY) / dy +
123 (pC[
a_y] * BGBY) / dy + (pC[
b_xz] * BGBY) / (2 *
dx) + (pC[
b_xyz] * BGBY) / (4 *
dx) -
124 (pC[
b_xyy] * BGBY) / (6 *
dx) - (pC[
b_xy] * BGBY) / (2 *
dx) - (pC[
b_x] * BGBY) /
dx -
126 (pC[
a_zz] * pC[
c_z]) / (2 * dz) - (pC[
a_z] * pC[
c_z]) / (2 * dz) - (pC[
a_yz] * pC[
c_z]) / (4 * dz) +
128 (pC[
a_zz] * pC[
c_y]) / (2 * dz) + (pC[
a_z] * pC[
c_y]) / (2 * dz) + (pC[
a_yz] * pC[
c_y]) / (4 * dz) +
132 (pC[
a_yz] * pC[
b_z]) / (4 * dy) - (pC[
a_yy] * pC[
b_z]) / (2 * dy) - (pC[
a_y] * pC[
b_z]) / (2 * dy) +
135 (pC[
a_yz] * pC[
b_y]) / (4 * dy) + (pC[
a_yy] * pC[
b_y]) / (2 * dy) + (pC[
a_y] * pC[
b_y]) / (2 * dy) -
173template<
typename REAL>
inline
175 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
178 const std::array<Real, 3>& gridSpacing
181 const auto dx = gridSpacing[0];
182 const auto dy = gridSpacing[1];
183 const auto dz = gridSpacing[2];
184 return (pC[
a_zz] * BGBZ) / dz + (pC[
a_z] * BGBZ) / dz - (pC[
a_yz] * BGBZ) / (2 * dz) -
186 (pC[
c_xy] * BGBZ) / (2 *
dx) - (pC[
c_x] * BGBZ) /
dx + (pC[
a_yz] * BGBY) / (2 * dy) - (pC[
a_yy] * BGBY) / dy +
187 (pC[
a_y] * BGBY) / dy - (pC[
b_xz] * BGBY) / (2 *
dx) + (pC[
b_xyz] * BGBY) / (4 *
dx) -
188 (pC[
b_xyy] * BGBY) / (6 *
dx) + (pC[
b_xy] * BGBY) / (2 *
dx) - (pC[
b_x] * BGBY) /
dx +
190 (pC[
a_zz] * pC[
c_z]) / (2 * dz) + (pC[
a_z] * pC[
c_z]) / (2 * dz) - (pC[
a_yz] * pC[
c_z]) / (4 * dz) -
192 (pC[
a_zz] * pC[
c_y]) / (2 * dz) - (pC[
a_z] * pC[
c_y]) / (2 * dz) + (pC[
a_yz] * pC[
c_y]) / (4 * dz) +
196 (pC[
a_yz] * pC[
b_z]) / (4 * dy) - (pC[
a_yy] * pC[
b_z]) / (2 * dy) + (pC[
a_y] * pC[
b_z]) / (2 * dy) -
199 (pC[
a_yz] * pC[
b_y]) / (4 * dy) + (pC[
a_yy] * pC[
b_y]) / (2 * dy) - (pC[
a_y] * pC[
b_y]) / (2 * dy) -
237template<
typename REAL>
inline
239 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
242 const std::array<Real, 3>& gridSpacing
245 const auto dx = gridSpacing[0];
246 const auto dy = gridSpacing[1];
247 const auto dz = gridSpacing[2];
248 return (pC[
a_zz] * BGBZ) / dz + (pC[
a_z] * BGBZ) / dz + (pC[
a_yz] * BGBZ) / (2 * dz) -
250 (pC[
c_xy] * BGBZ) / (2 *
dx) - (pC[
c_x] * BGBZ) /
dx + (pC[
a_yz] * BGBY) / (2 * dy) + (pC[
a_yy] * BGBY) / dy +
251 (pC[
a_y] * BGBY) / dy - (pC[
b_xz] * BGBY) / (2 *
dx) - (pC[
b_xyz] * BGBY) / (4 *
dx) -
252 (pC[
b_xyy] * BGBY) / (6 *
dx) - (pC[
b_xy] * BGBY) / (2 *
dx) - (pC[
b_x] * BGBY) /
dx +
254 (pC[
a_zz] * pC[
c_z]) / (2 * dz) + (pC[
a_z] * pC[
c_z]) / (2 * dz) + (pC[
a_yz] * pC[
c_z]) / (4 * dz) +
256 (pC[
a_zz] * pC[
c_y]) / (2 * dz) + (pC[
a_z] * pC[
c_y]) / (2 * dz) + (pC[
a_yz] * pC[
c_y]) / (4 * dz) +
260 (pC[
a_yz] * pC[
b_z]) / (4 * dy) + (pC[
a_yy] * pC[
b_z]) / (2 * dy) + (pC[
a_y] * pC[
b_z]) / (2 * dy) +
263 (pC[
a_yz] * pC[
b_y]) / (4 * dy) + (pC[
a_yy] * pC[
b_y]) / (2 * dy) + (pC[
a_y] * pC[
b_y]) / (2 * dy) +
302template<
typename REAL>
inline
304 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
307 const std::array<Real, 3>& gridSpacing
310 const auto dx = gridSpacing[0];
311 const auto dy = gridSpacing[1];
312 const auto dz = gridSpacing[2];
313 return -(pC[
b_zz] * BGBZ) / dz + (pC[
b_z] * BGBZ) / dz - (pC[
b_xz] * BGBZ) / (2 * dz) -
314 (pC[
c_yzz] * BGBZ) / (6 * dy) + (pC[
c_yz] * BGBZ) / (2 * dy) - (pC[
c_y] * BGBZ) / dy -
315 (pC[
c_xyz] * BGBZ) / (4 * dy) + (pC[
c_xy] * BGBZ) / (2 * dy) + (pC[
a_yz] * BGBX) / (2 * dy) -
316 (pC[
a_y] * BGBX) / dy - (pC[
a_xyz] * BGBX) / (4 * dy) + (pC[
a_xy] * BGBX) / (2 * dy) -
317 (pC[
a_xxy] * BGBX) / (6 * dy) - (pC[
b_xz] * BGBX) / (2 *
dx) - (pC[
b_xx] * BGBX) /
dx +
319 (pC[
b_xz] * pC[
c_zz]) / (12 * dz) + (pC[
b_zz] * pC[
c_z]) / (2 * dz) - (pC[
b_z] * pC[
c_z]) / (2 * dz) +
323 (pC[
b_xz] * pC[
c_xz]) / (8 * dz) + (pC[
b_zz] * pC[
c_x]) / (2 * dz) - (pC[
b_z] * pC[
c_x]) / (2 * dz) +
327 (pC[
c_yzz] * pC[
c_z]) / (12 * dy) - (pC[
c_yz] * pC[
c_z]) / (4 * dy) + (pC[
c_y] * pC[
c_z]) / (2 * dy) +
332 (pC[
c_xz] * pC[
c_y]) / (4 * dy) + (pC[
c_x] * pC[
c_y]) / (2 * dy) - (pC[
c_0] * pC[
c_y]) / dy -
335 (pC[
a_yz] * pC[
a_z]) / (4 * dy) + (pC[
a_y] * pC[
a_z]) / (2 * dy) + (pC[
a_xyz] * pC[
a_z]) / (8 * dy) -
339 (pC[
a_xyy] * pC[
a_y]) / (12 * dy) - (pC[
a_xx] * pC[
a_y]) / (6 * dy) + (pC[
a_x] * pC[
a_y]) / (2 * dy) -
366template<
typename REAL>
inline
368 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
371 const std::array<Real, 3>& gridSpacing
374 const auto dx = gridSpacing[0];
375 const auto dy = gridSpacing[1];
376 const auto dz = gridSpacing[2];
377 return -(pC[
b_zz] * BGBZ) / dz + (pC[
b_z] * BGBZ) / dz + (pC[
b_xz] * BGBZ) / (2 * dz) -
378 (pC[
c_yzz] * BGBZ) / (6 * dy) + (pC[
c_yz] * BGBZ) / (2 * dy) - (pC[
c_y] * BGBZ) / dy +
379 (pC[
c_xyz] * BGBZ) / (4 * dy) - (pC[
c_xy] * BGBZ) / (2 * dy) + (pC[
a_yz] * BGBX) / (2 * dy) -
380 (pC[
a_y] * BGBX) / dy + (pC[
a_xyz] * BGBX) / (4 * dy) - (pC[
a_xy] * BGBX) / (2 * dy) -
381 (pC[
a_xxy] * BGBX) / (6 * dy) - (pC[
b_xz] * BGBX) / (2 *
dx) + (pC[
b_xx] * BGBX) /
dx +
383 (pC[
b_xz] * pC[
c_zz]) / (12 * dz) + (pC[
b_zz] * pC[
c_z]) / (2 * dz) - (pC[
b_z] * pC[
c_z]) / (2 * dz) -
387 (pC[
b_xz] * pC[
c_xz]) / (8 * dz) - (pC[
b_zz] * pC[
c_x]) / (2 * dz) + (pC[
b_z] * pC[
c_x]) / (2 * dz) +
391 (pC[
c_yzz] * pC[
c_z]) / (12 * dy) - (pC[
c_yz] * pC[
c_z]) / (4 * dy) + (pC[
c_y] * pC[
c_z]) / (2 * dy) -
396 (pC[
c_xz] * pC[
c_y]) / (4 * dy) - (pC[
c_x] * pC[
c_y]) / (2 * dy) - (pC[
c_0] * pC[
c_y]) / dy -
399 (pC[
a_yz] * pC[
a_z]) / (4 * dy) + (pC[
a_y] * pC[
a_z]) / (2 * dy) - (pC[
a_xyz] * pC[
a_z]) / (8 * dy) +
403 (pC[
a_xyy] * pC[
a_y]) / (12 * dy) - (pC[
a_xx] * pC[
a_y]) / (6 * dy) - (pC[
a_x] * pC[
a_y]) / (2 * dy) -
430template<
typename REAL>
inline
432 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
435 const std::array<Real, 3>& gridSpacing
438 const auto dx = gridSpacing[0];
439 const auto dy = gridSpacing[1];
440 const auto dz = gridSpacing[2];
441 return (pC[
b_zz] * BGBZ) / dz + (pC[
b_z] * BGBZ) / dz - (pC[
b_xz] * BGBZ) / (2 * dz) -
442 (pC[
c_yzz] * BGBZ) / (6 * dy) - (pC[
c_yz] * BGBZ) / (2 * dy) - (pC[
c_y] * BGBZ) / dy +
443 (pC[
c_xyz] * BGBZ) / (4 * dy) + (pC[
c_xy] * BGBZ) / (2 * dy) - (pC[
a_yz] * BGBX) / (2 * dy) -
444 (pC[
a_y] * BGBX) / dy + (pC[
a_xyz] * BGBX) / (4 * dy) + (pC[
a_xy] * BGBX) / (2 * dy) -
445 (pC[
a_xxy] * BGBX) / (6 * dy) + (pC[
b_xz] * BGBX) / (2 *
dx) - (pC[
b_xx] * BGBX) /
dx +
447 (pC[
b_xz] * pC[
c_zz]) / (12 * dz) + (pC[
b_zz] * pC[
c_z]) / (2 * dz) + (pC[
b_z] * pC[
c_z]) / (2 * dz) -
451 (pC[
b_xz] * pC[
c_xz]) / (8 * dz) - (pC[
b_zz] * pC[
c_x]) / (2 * dz) - (pC[
b_z] * pC[
c_x]) / (2 * dz) +
455 (pC[
c_yzz] * pC[
c_z]) / (12 * dy) - (pC[
c_yz] * pC[
c_z]) / (4 * dy) - (pC[
c_y] * pC[
c_z]) / (2 * dy) +
460 (pC[
c_xz] * pC[
c_y]) / (4 * dy) + (pC[
c_x] * pC[
c_y]) / (2 * dy) - (pC[
c_0] * pC[
c_y]) / dy -
463 (pC[
a_yz] * pC[
a_z]) / (4 * dy) - (pC[
a_y] * pC[
a_z]) / (2 * dy) + (pC[
a_xyz] * pC[
a_z]) / (8 * dy) +
467 (pC[
a_xyy] * pC[
a_y]) / (12 * dy) - (pC[
a_xx] * pC[
a_y]) / (6 * dy) + (pC[
a_x] * pC[
a_y]) / (2 * dy) -
494template<
typename REAL>
inline
496 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
499 const std::array<Real, 3>& gridSpacing
502 const auto dx = gridSpacing[0];
503 const auto dy = gridSpacing[1];
504 const auto dz = gridSpacing[2];
505 return (pC[
b_zz] * BGBZ) / dz + (pC[
b_z] * BGBZ) / dz + (pC[
b_xz] * BGBZ) / (2 * dz) -
506 (pC[
c_yzz] * BGBZ) / (6 * dy) - (pC[
c_yz] * BGBZ) / (2 * dy) - (pC[
c_y] * BGBZ) / dy -
507 (pC[
c_xyz] * BGBZ) / (4 * dy) - (pC[
c_xy] * BGBZ) / (2 * dy) - (pC[
a_yz] * BGBX) / (2 * dy) -
508 (pC[
a_y] * BGBX) / dy - (pC[
a_xyz] * BGBX) / (4 * dy) - (pC[
a_xy] * BGBX) / (2 * dy) -
509 (pC[
a_xxy] * BGBX) / (6 * dy) + (pC[
b_xz] * BGBX) / (2 *
dx) + (pC[
b_xx] * BGBX) /
dx +
511 (pC[
b_xz] * pC[
c_zz]) / (12 * dz) + (pC[
b_zz] * pC[
c_z]) / (2 * dz) + (pC[
b_z] * pC[
c_z]) / (2 * dz) +
515 (pC[
b_xz] * pC[
c_xz]) / (8 * dz) + (pC[
b_zz] * pC[
c_x]) / (2 * dz) + (pC[
b_z] * pC[
c_x]) / (2 * dz) +
519 (pC[
c_yzz] * pC[
c_z]) / (12 * dy) - (pC[
c_yz] * pC[
c_z]) / (4 * dy) - (pC[
c_y] * pC[
c_z]) / (2 * dy) -
524 (pC[
c_xz] * pC[
c_y]) / (4 * dy) - (pC[
c_x] * pC[
c_y]) / (2 * dy) - (pC[
c_0] * pC[
c_y]) / dy -
527 (pC[
a_yz] * pC[
a_z]) / (4 * dy) - (pC[
a_y] * pC[
a_z]) / (2 * dy) - (pC[
a_xyz] * pC[
a_z]) / (8 * dy) -
531 (pC[
a_xyy] * pC[
a_y]) / (12 * dy) - (pC[
a_xx] * pC[
a_y]) / (6 * dy) - (pC[
a_x] * pC[
a_y]) / (2 * dy) -
559template<
typename REAL>
inline
561 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
564 const std::array<Real, 3>& gridSpacing
567 const auto dx = gridSpacing[0];
568 const auto dy = gridSpacing[1];
569 const auto dz = gridSpacing[2];
570 return -(pC[
b_z] * BGBY) / dz + (pC[
b_yz] * BGBY) / (2 * dz) - (pC[
b_yyz] * BGBY) / (6 * dz) +
571 (pC[
b_xz] * BGBY) / (2 * dz) - (pC[
b_xyz] * BGBY) / (4 * dz) - (pC[
c_yy] * BGBY) / dy +
572 (pC[
c_y] * BGBY) / dy - (pC[
c_xy] * BGBY) / (2 * dy) - (pC[
a_z] * BGBX) / dz + (pC[
a_yz] * BGBX) / (2 * dz) +
573 (pC[
a_xz] * BGBX) / (2 * dz) - (pC[
a_xyz] * BGBX) / (4 * dz) - (pC[
a_xxz] * BGBX) / (6 * dz) -
576 (pC[
b_yy] * pC[
b_z]) / (6 * dz) + (pC[
b_y] * pC[
b_z]) / (2 * dz) - (pC[
b_xy] * pC[
b_z]) / (4 * dz) +
586 (pC[
a_xy] * pC[
a_z]) / (4 * dz) - (pC[
a_xx] * pC[
a_z]) / (6 * dz) + (pC[
a_x] * pC[
a_z]) / (2 * dz) -
598 (pC[
b_xy] * pC[
c_y]) / (4 * dy) - (pC[
b_x] * pC[
c_y]) / (2 * dy) + (pC[
b_0] * pC[
c_y]) / dy -
623template<
typename REAL>
inline
625 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
628 const std::array<Real, 3>& gridSpacing
631 const auto dx = gridSpacing[0];
632 const auto dy = gridSpacing[1];
633 const auto dz = gridSpacing[2];
634 return -(pC[
b_z] * BGBY) / dz + (pC[
b_yz] * BGBY) / (2 * dz) - (pC[
b_yyz] * BGBY) / (6 * dz) -
635 (pC[
b_xz] * BGBY) / (2 * dz) + (pC[
b_xyz] * BGBY) / (4 * dz) - (pC[
c_yy] * BGBY) / dy +
636 (pC[
c_y] * BGBY) / dy + (pC[
c_xy] * BGBY) / (2 * dy) - (pC[
a_z] * BGBX) / dz + (pC[
a_yz] * BGBX) / (2 * dz) -
637 (pC[
a_xz] * BGBX) / (2 * dz) + (pC[
a_xyz] * BGBX) / (4 * dz) - (pC[
a_xxz] * BGBX) / (6 * dz) -
640 (pC[
b_yy] * pC[
b_z]) / (6 * dz) + (pC[
b_y] * pC[
b_z]) / (2 * dz) + (pC[
b_xy] * pC[
b_z]) / (4 * dz) -
650 (pC[
a_xy] * pC[
a_z]) / (4 * dz) - (pC[
a_xx] * pC[
a_z]) / (6 * dz) - (pC[
a_x] * pC[
a_z]) / (2 * dz) -
662 (pC[
b_xy] * pC[
c_y]) / (4 * dy) + (pC[
b_x] * pC[
c_y]) / (2 * dy) + (pC[
b_0] * pC[
c_y]) / dy +
687template<
typename REAL>
inline
689 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
692 const std::array<Real, 3>& gridSpacing
695 const auto dx = gridSpacing[0];
696 const auto dy = gridSpacing[1];
697 const auto dz = gridSpacing[2];
698 return -(pC[
b_z] * BGBY) / dz - (pC[
b_yz] * BGBY) / (2 * dz) - (pC[
b_yyz] * BGBY) / (6 * dz) +
699 (pC[
b_xz] * BGBY) / (2 * dz) + (pC[
b_xyz] * BGBY) / (4 * dz) + (pC[
c_yy] * BGBY) / dy +
700 (pC[
c_y] * BGBY) / dy - (pC[
c_xy] * BGBY) / (2 * dy) - (pC[
a_z] * BGBX) / dz - (pC[
a_yz] * BGBX) / (2 * dz) +
701 (pC[
a_xz] * BGBX) / (2 * dz) + (pC[
a_xyz] * BGBX) / (4 * dz) - (pC[
a_xxz] * BGBX) / (6 * dz) +
704 (pC[
b_yy] * pC[
b_z]) / (6 * dz) - (pC[
b_y] * pC[
b_z]) / (2 * dz) + (pC[
b_xy] * pC[
b_z]) / (4 * dz) +
714 (pC[
a_xy] * pC[
a_z]) / (4 * dz) - (pC[
a_xx] * pC[
a_z]) / (6 * dz) + (pC[
a_x] * pC[
a_z]) / (2 * dz) -
726 (pC[
b_xy] * pC[
c_y]) / (4 * dy) - (pC[
b_x] * pC[
c_y]) / (2 * dy) + (pC[
b_0] * pC[
c_y]) / dy -
751template<
typename REAL>
inline
753 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
756 const std::array<Real, 3>& gridSpacing
759 const auto dx = gridSpacing[0];
760 const auto dy = gridSpacing[1];
761 const auto dz = gridSpacing[2];
762 return -(pC[
b_z] * BGBY) / dz - (pC[
b_yz] * BGBY) / (2 * dz) - (pC[
b_yyz] * BGBY) / (6 * dz) -
763 (pC[
b_xz] * BGBY) / (2 * dz) - (pC[
b_xyz] * BGBY) / (4 * dz) + (pC[
c_yy] * BGBY) / dy +
764 (pC[
c_y] * BGBY) / dy + (pC[
c_xy] * BGBY) / (2 * dy) - (pC[
a_z] * BGBX) / dz - (pC[
a_yz] * BGBX) / (2 * dz) -
765 (pC[
a_xz] * BGBX) / (2 * dz) - (pC[
a_xyz] * BGBX) / (4 * dz) - (pC[
a_xxz] * BGBX) / (6 * dz) +
768 (pC[
b_yy] * pC[
b_z]) / (6 * dz) - (pC[
b_y] * pC[
b_z]) / (2 * dz) - (pC[
b_xy] * pC[
b_z]) / (4 * dz) -
778 (pC[
a_xy] * pC[
a_z]) / (4 * dz) - (pC[
a_xx] * pC[
a_z]) / (6 * dz) - (pC[
a_x] * pC[
a_z]) / (2 * dz) -
790 (pC[
b_xy] * pC[
c_y]) / (4 * dy) + (pC[
b_x] * pC[
c_y]) / (2 * dy) + (pC[
b_0] * pC[
c_y]) / dy +
803template <
typename REAL>
805 Real BGBZ,
const std::array<Real, 3>& gridSpacing) {
871 const std::array<Real, 3>& gridSpacing,
872 const std::array<Real, Rec::N_REC_COEFFICIENTS>& perturbedCoefficients,
873 const fsgrid::FsStencil& stencil) {
874 const auto ooo = stencil.ooo();
875 const auto& bgb = bgbs[ooo];
876 const auto& perb = perbs[ooo];
877 const auto& dperb = dperbs[ooo];
878 const auto& moment = moments[ooo];
879 auto& ehall = ehalls[ooo];
885 auto computeHallRhoq = [&moments, &moment](
const std::array<size_t, 4>& indices) {
887 const auto max = std::numeric_limits<Real>::max();
899 cerr << __FILE__ << __LINE__ <<
"You shouldn't be in a Hall term function if Parameters::ohmHallTerm == 0."
917 const Real EXHall = (Bz * (xdz - zdx) - By * (ydx - xdy)) * invHallRhoqMU0;
923 const Real EYHall = (Bx * (ydx - xdy) - Bz * (zdy - ydz)) * invHallRhoqMU0;
929 const Real EZHall = (By * (zdy - ydz) - Bx * (xdz - zdx)) * invHallRhoqMU0;
938 auto computeEHall = [&perturbedCoefficients, &bgbx, &bgby, &bgbz, &gridSpacing](
fsgrids::ehall term,
943 const auto omo = stencil.omo();
944 const auto opo = stencil.opo();
945 const auto moo = stencil.moo();
946 const auto poo = stencil.poo();
947 const auto oom = stencil.oom();
948 const auto oop = stencil.oop();
950 const std::array<std::array<size_t, 4>, 12> indices = {
1026 const std::array<fsgrids::ehall, 12> terms = {
1034 ehall[terms[
index]] = computeEHall(terms[
index], computeHallRhoq(indices[
index]));
1041 cerr << __FILE__ <<
":" << __LINE__ <<
"You are welcome to code higher-order Hall term correction terms." << endl;
1066 SysBoundary& sysBoundaries,
const std::array<Real, 3>& gridSpacing) {
1068 if (!stencil.cellExists(0, 0, 0)) {
1069 cerr <<
"Out-of-bounds access in " << __FILE__ <<
":" << __LINE__ << endl;
1074 const auto& tech = technical[stencil.ooo()];
1075 cuint cellSysBoundaryFlag = tech.sysBoundaryFlag;
1076 cuint cellSysBoundaryLayer = tech.sysBoundaryLayer;
1085 const std::array<Real, Rec::N_REC_COEFFICIENTS> perturbedCoefficients =
1089 auto*
const sb = sysBoundaries.
getSysBoundary(cellSysBoundaryFlag);
1091 sb->fieldSolverBoundaryCondHallElectricField(ehall, stencil, 1);
1092 sb->fieldSolverBoundaryCondHallElectricField(ehall, stencil, 2);
1130 SysBoundary& sysBoundaries, int32_t RKCase,
const bool communicateMomentsDerivatives) {
1135 moments = momentsdt2;
1136 dmoments = dmomentsdt2;
1138 phiprof::Timer hallTimer{
"Calculate Hall term"};
1140 const size_t numCells =
fsgrid.getNumCells();
1142 phiprof::Timer mpiTimer{
"EHall ghost updates MPI", {
"MPI"}};
1143 fsgrid.updateGhostCells(dperb);
1145 fsgrid.updateGhostCells(dmoments);
1149 fsgrid.parallel_for([](
int timerId) -> phiprof::Timer {
return phiprof::Timer{timerId}; },
1150 phiprof::initializeTimer(
"EHall compute cells"), technical,
1151 [=, &sysBoundaries](
const fsgrid::Coordinates &coordinates,
const fsgrid::FsStencil& stencil,
cuint sysBoundaryFlag,
cuint sysBoundaryLayer) {
1152 calculateHallTerm(perb, ehall, moments, dperb, bgb, technical, stencil, sysBoundaries, coordinates.physicalGridSpacing);
1155 hallTimer.stop(numCells,
"Spatial Cells");
virtual void fieldSolverBoundaryCondHallElectricField(fsgrids::ehallspan ehall, const fsgrid::FsStencil &stencil, cuint component)=0
SysBoundary contains the SysBoundaryConditions used in the simulation.
SBC::SysBoundaryCondition * getSysBoundary(cuint sysBoundaryType) const
fsgrid::FsGrid< FS_STENCIL_WIDTH > FieldSolverGrid
std::array< Real, Rec::N_REC_COEFFICIENTS > reconstructionCoefficients(fsgrids::perbspan perb, fsgrids::constdperbspan dperb, const fsgrid::FsStencil &stencil, Real reconstructionOrder)
Low-level helper function.
REAL JXB(fsgrids::ehall term, const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, Real BGBX, Real BGBY, Real BGBZ, const std::array< Real, 3 > &gridSpacing)
REAL JXBY_100_110(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBX_011_111(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBY, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
void calculateHallTerm(fsgrids::perbspan perb, fsgrids::ehallspan ehall, fsgrids::constmomentsspan moments, fsgrids::constdperbspan dperb, fsgrids::constbgbspan bgb, fsgrids::consttechnicalspan technical, const fsgrid::FsStencil &stencil, SysBoundary &sysBoundaries, const std::array< Real, 3 > &gridSpacing)
Calculate the numerator of the Hall term on all given cells.
REAL JXBX_000_100(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBY, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBZ_000_001(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBY, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBY_001_011(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBY_000_010(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBZ_010_011(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBY, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBZ_100_101(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBY, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBX_001_101(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBY, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
REAL JXBX_010_110(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBY, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
void calculateHallTermSimple(fsgrids::perbspan perb, fsgrids::perbspan perbdt2, fsgrids::ehallspan ehall, fsgrids::momentsspan moments, fsgrids::momentsspan momentsdt2, fsgrids::dperbspan dperb, fsgrids::dmomentsspan dmoments, fsgrids::dmomentsspan dmomentsdt2, fsgrids::bgbspan bgb, fsgrids::technicalspan technical, FieldSolverGrid &fsgrid, SysBoundary &sysBoundaries, int32_t RKCase, const bool communicateMomentsDerivatives)
High-level function computing the Hall term.
REAL JXBZ_110_111(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBY, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
void calculateEdgeHallTermComponents(fsgrids::perbspan perbs, fsgrids::ehallspan ehalls, fsgrids::constmomentsspan moments, fsgrids::constdperbspan dperbs, fsgrids::constbgbspan bgbs, const std::array< Real, 3 > &gridSpacing, const std::array< Real, Rec::N_REC_COEFFICIENTS > &perturbedCoefficients, const fsgrid::FsStencil &stencil)
Low-level function computing the Hall term numerator x components.
REAL JXBY_101_111(const std::array< REAL, Rec::N_REC_COEFFICIENTS > &pC, creal BGBX, creal BGBZ, const std::array< Real, 3 > &gridSpacing)
Low-level Hall component computation.
std::span< std::array< Real, fsgrids::bfield::N_BFIELD > > perbspan
std::span< const std::array< Real, bgbfield::N_BGB > > constbgbspan
std::span< const std::array< Real, fsgrids::moments::N_MOMENTS > > constmomentsspan
std::span< std::array< Real, fsgrids::moments::N_MOMENTS > > momentsspan
std::span< std::array< Real, fsgrids::dmoments::N_DMOMENTS > > dmomentsspan
std::span< technical > technicalspan
std::span< std::array< Real, bgbfield::N_BGB > > bgbspan
std::span< std::array< Real, fsgrids::dperb::N_DPERB > > dperbspan
std::span< const technical > consttechnicalspan
std::span< std::array< Real, fsgrids::ehall::N_EHALL > > ehallspan
std::span< const std::array< Real, fsgrids::dperb::N_DPERB > > constdperbspan
static uint ohmGradPeTerm
static Real hallMinimumRhoq
static ARCH_HOSTDEV VecSimple< T > min(VecSimple< T > const &l, VecSimple< T > const &r)
static ARCH_HOSTDEV VecSimple< T > max(VecSimple< T > const &l, VecSimple< T > const &r)