Vlasiator ebf0dd394 on dev (v5.4.0 + 1054 commits)
Loading...
Searching...
No Matches
ldz_hall.cpp
Go to the documentation of this file.
1/*
2 * This file is part of Vlasiator.
3 * Copyright 2010-2016 Finnish Meteorological Institute
4 *
5 * For details of usage, see the COPYING file and read the "Rules of the Road"
6 * at http://www.physics.helsinki.fi/vlasiator/
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License along
19 * with this program; if not, write to the Free Software Foundation, Inc.,
20 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
21 */
22
23#include "fs_common.h"
24#include "ldz_hall.hpp"
25#include <limits>
26
27#ifdef DEBUG_VLASIATOR
28 #define DEBUG_FSOLVER
29#endif
30
31using namespace std;
32
45template<typename REAL> inline
47 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
48 creal BGBY,
49 creal BGBZ,
50 const std::array<Real, 3>& gridSpacing
51) {
52 using namespace Rec;
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) -
57 (pC[c_xzz] * BGBZ) / (6 * dx) + (pC[c_xz] * BGBZ) / (2 * dx) - (pC[c_xyz] * BGBZ) / (4 * dx) +
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) -
60 (pC[b_xyy] * BGBY) / (6 * dx) + (pC[b_xy] * BGBY) / (2 * dx) - (pC[b_x] * BGBY) / dx -
61 (pC[a_zz] * pC[c_zz]) / (6 * dz) + (pC[a_z] * pC[c_zz]) / (6 * dz) - (pC[a_yz] * pC[c_zz]) / (12 * dz) +
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) -
63 (pC[a_zz] * pC[c_yz]) / (4 * dz) + (pC[a_z] * pC[c_yz]) / (4 * dz) - (pC[a_yz] * pC[c_yz]) / (8 * 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) +
65 (pC[a_xzz] * pC[c_xz]) / (24 * dz) - (pC[a_xz] * pC[c_xz]) / (24 * dz) + (pC[a_xyz] * pC[c_xz]) / (48 * dz) -
66 (pC[a_xzz] * pC[c_x]) / (12 * dz) + (pC[a_xz] * pC[c_x]) / (12 * dz) - (pC[a_xyz] * pC[c_x]) / (24 * dz) -
67 (pC[a_zz] * pC[c_0]) / dz + (pC[a_z] * pC[c_0]) / dz - (pC[a_yz] * pC[c_0]) / (2 * 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) -
69 (pC[a_yz] * pC[b_yz]) / (8 * dy) - (pC[a_yy] * pC[b_yz]) / (4 * dy) + (pC[a_y] * pC[b_yz]) / (4 * dy) -
70 (pC[a_yz] * pC[b_yy]) / (12 * dy) - (pC[a_yy] * pC[b_yy]) / (6 * dy) + (pC[a_y] * pC[b_yy]) / (6 * 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) +
72 (pC[a_xyz] * pC[b_xy]) / (48 * dy) + (pC[a_xyy] * pC[b_xy]) / (24 * dy) - (pC[a_xy] * pC[b_xy]) / (24 * dy) -
73 (pC[a_xyz] * pC[b_x]) / (24 * dy) - (pC[a_xyy] * pC[b_x]) / (12 * dy) + (pC[a_xy] * pC[b_x]) / (12 * dy) -
74 (pC[a_yz] * pC[b_0]) / (2 * dy) - (pC[a_yy] * pC[b_0]) / dy + (pC[a_y] * pC[b_0]) / dy -
75 (pC[c_xzz] * pC[c_zz]) / (36 * dx) + (pC[c_xz] * pC[c_zz]) / (12 * dx) - (pC[c_xyz] * pC[c_zz]) / (24 * dx) +
76 (pC[c_xy] * pC[c_zz]) / (12 * dx) - (pC[c_x] * pC[c_zz]) / (6 * dx) + (pC[c_xzz] * pC[c_z]) / (12 * dx) -
77 (pC[c_xz] * pC[c_z]) / (4 * dx) + (pC[c_xyz] * pC[c_z]) / (8 * dx) - (pC[c_xy] * pC[c_z]) / (4 * dx) +
78 (pC[c_x] * pC[c_z]) / (2 * dx) - (pC[c_xzz] * pC[c_yz]) / (24 * dx) + (pC[c_xz] * pC[c_yz]) / (8 * dx) -
79 (pC[c_xyz] * pC[c_yz]) / (16 * dx) + (pC[c_xy] * pC[c_yz]) / (8 * dx) - (pC[c_x] * pC[c_yz]) / (4 * dx) +
80 (pC[c_xzz] * pC[c_y]) / (12 * dx) - (pC[c_xz] * pC[c_y]) / (4 * dx) + (pC[c_xyz] * pC[c_y]) / (8 * dx) -
81 (pC[c_xy] * pC[c_y]) / (4 * dx) + (pC[c_x] * pC[c_y]) / (2 * dx) - (pC[c_0] * pC[c_xzz]) / (6 * dx) -
82 (pC[c_xxz] * pC[c_xz]) / (24 * dx) + (pC[c_xx] * pC[c_xz]) / (12 * dx) + (pC[c_0] * pC[c_xz]) / (2 * dx) -
83 (pC[c_0] * pC[c_xyz]) / (4 * dx) + (pC[c_0] * pC[c_xy]) / (2 * dx) + (pC[c_x] * pC[c_xxz]) / (12 * dx) -
84 (pC[c_x] * pC[c_xx]) / (6 * dx) - (pC[c_0] * pC[c_x]) / dx - (pC[b_xz] * pC[b_z]) / (4 * dx) +
85 (pC[b_xyz] * pC[b_z]) / (8 * dx) + (pC[b_xyy] * pC[b_z]) / (12 * dx) - (pC[b_xy] * pC[b_z]) / (4 * dx) +
86 (pC[b_x] * pC[b_z]) / (2 * dx) + (pC[b_xz] * pC[b_yz]) / (8 * dx) - (pC[b_xyz] * pC[b_yz]) / (16 * dx) -
87 (pC[b_xyy] * pC[b_yz]) / (24 * dx) + (pC[b_xy] * pC[b_yz]) / (8 * dx) - (pC[b_x] * pC[b_yz]) / (4 * dx) +
88 (pC[b_xz] * pC[b_yy]) / (12 * dx) - (pC[b_xyz] * pC[b_yy]) / (24 * dx) - (pC[b_xyy] * pC[b_yy]) / (36 * dx) +
89 (pC[b_xy] * pC[b_yy]) / (12 * dx) - (pC[b_x] * pC[b_yy]) / (6 * dx) - (pC[b_xz] * pC[b_y]) / (4 * dx) +
90 (pC[b_xyz] * pC[b_y]) / (8 * dx) + (pC[b_xyy] * pC[b_y]) / (12 * dx) - (pC[b_xy] * pC[b_y]) / (4 * dx) +
91 (pC[b_x] * pC[b_y]) / (2 * dx) + (pC[b_0] * pC[b_xz]) / (2 * dx) - (pC[b_0] * pC[b_xyz]) / (4 * dx) -
92 (pC[b_0] * pC[b_xyy]) / (6 * dx) - (pC[b_xxy] * pC[b_xy]) / (24 * dx) + (pC[b_xx] * pC[b_xy]) / (12 * dx) +
93 (pC[b_0] * pC[b_xy]) / (2 * dx) + (pC[b_x] * pC[b_xxy]) / (12 * dx) - (pC[b_x] * pC[b_xx]) / (6 * dx) -
94 (pC[b_0] * pC[b_x]) / dx;
95}
96
109template<typename REAL> inline
111 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
112 creal BGBY,
113 creal BGBZ,
114 const std::array<Real, 3>& gridSpacing
115) {
116 using namespace Rec;
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) -
121 (pC[c_xzz] * BGBZ) / (6 * dx) + (pC[c_xz] * BGBZ) / (2 * dx) + (pC[c_xyz] * BGBZ) / (4 * dx) -
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 -
125 (pC[a_zz] * pC[c_zz]) / (6 * dz) + (pC[a_z] * pC[c_zz]) / (6 * dz) + (pC[a_yz] * pC[c_zz]) / (12 * dz) +
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) +
127 (pC[a_zz] * pC[c_yz]) / (4 * dz) - (pC[a_z] * pC[c_yz]) / (4 * dz) - (pC[a_yz] * pC[c_yz]) / (8 * 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) +
129 (pC[a_xzz] * pC[c_xz]) / (24 * dz) - (pC[a_xz] * pC[c_xz]) / (24 * dz) - (pC[a_xyz] * pC[c_xz]) / (48 * dz) -
130 (pC[a_xzz] * pC[c_x]) / (12 * dz) + (pC[a_xz] * pC[c_x]) / (12 * dz) + (pC[a_xyz] * pC[c_x]) / (24 * dz) -
131 (pC[a_zz] * pC[c_0]) / dz + (pC[a_z] * pC[c_0]) / dz + (pC[a_yz] * pC[c_0]) / (2 * 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) +
133 (pC[a_yz] * pC[b_yz]) / (8 * dy) - (pC[a_yy] * pC[b_yz]) / (4 * dy) - (pC[a_y] * pC[b_yz]) / (4 * dy) -
134 (pC[a_yz] * pC[b_yy]) / (12 * dy) + (pC[a_yy] * pC[b_yy]) / (6 * dy) + (pC[a_y] * pC[b_yy]) / (6 * 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) -
136 (pC[a_xyz] * pC[b_xy]) / (48 * dy) + (pC[a_xyy] * pC[b_xy]) / (24 * dy) + (pC[a_xy] * pC[b_xy]) / (24 * dy) -
137 (pC[a_xyz] * pC[b_x]) / (24 * dy) + (pC[a_xyy] * pC[b_x]) / (12 * dy) + (pC[a_xy] * pC[b_x]) / (12 * dy) -
138 (pC[a_yz] * pC[b_0]) / (2 * dy) + (pC[a_yy] * pC[b_0]) / dy + (pC[a_y] * pC[b_0]) / dy -
139 (pC[c_xzz] * pC[c_zz]) / (36 * dx) + (pC[c_xz] * pC[c_zz]) / (12 * dx) + (pC[c_xyz] * pC[c_zz]) / (24 * dx) -
140 (pC[c_xy] * pC[c_zz]) / (12 * dx) - (pC[c_x] * pC[c_zz]) / (6 * dx) + (pC[c_xzz] * pC[c_z]) / (12 * dx) -
141 (pC[c_xz] * pC[c_z]) / (4 * dx) - (pC[c_xyz] * pC[c_z]) / (8 * dx) + (pC[c_xy] * pC[c_z]) / (4 * dx) +
142 (pC[c_x] * pC[c_z]) / (2 * dx) + (pC[c_xzz] * pC[c_yz]) / (24 * dx) - (pC[c_xz] * pC[c_yz]) / (8 * dx) -
143 (pC[c_xyz] * pC[c_yz]) / (16 * dx) + (pC[c_xy] * pC[c_yz]) / (8 * dx) + (pC[c_x] * pC[c_yz]) / (4 * dx) -
144 (pC[c_xzz] * pC[c_y]) / (12 * dx) + (pC[c_xz] * pC[c_y]) / (4 * dx) + (pC[c_xyz] * pC[c_y]) / (8 * dx) -
145 (pC[c_xy] * pC[c_y]) / (4 * dx) - (pC[c_x] * pC[c_y]) / (2 * dx) - (pC[c_0] * pC[c_xzz]) / (6 * dx) -
146 (pC[c_xxz] * pC[c_xz]) / (24 * dx) + (pC[c_xx] * pC[c_xz]) / (12 * dx) + (pC[c_0] * pC[c_xz]) / (2 * dx) +
147 (pC[c_0] * pC[c_xyz]) / (4 * dx) - (pC[c_0] * pC[c_xy]) / (2 * dx) + (pC[c_x] * pC[c_xxz]) / (12 * dx) -
148 (pC[c_x] * pC[c_xx]) / (6 * dx) - (pC[c_0] * pC[c_x]) / dx - (pC[b_xz] * pC[b_z]) / (4 * dx) -
149 (pC[b_xyz] * pC[b_z]) / (8 * dx) + (pC[b_xyy] * pC[b_z]) / (12 * dx) + (pC[b_xy] * pC[b_z]) / (4 * dx) +
150 (pC[b_x] * pC[b_z]) / (2 * dx) - (pC[b_xz] * pC[b_yz]) / (8 * dx) - (pC[b_xyz] * pC[b_yz]) / (16 * dx) +
151 (pC[b_xyy] * pC[b_yz]) / (24 * dx) + (pC[b_xy] * pC[b_yz]) / (8 * dx) + (pC[b_x] * pC[b_yz]) / (4 * dx) +
152 (pC[b_xz] * pC[b_yy]) / (12 * dx) + (pC[b_xyz] * pC[b_yy]) / (24 * dx) - (pC[b_xyy] * pC[b_yy]) / (36 * dx) -
153 (pC[b_xy] * pC[b_yy]) / (12 * dx) - (pC[b_x] * pC[b_yy]) / (6 * dx) + (pC[b_xz] * pC[b_y]) / (4 * dx) +
154 (pC[b_xyz] * pC[b_y]) / (8 * dx) - (pC[b_xyy] * pC[b_y]) / (12 * dx) - (pC[b_xy] * pC[b_y]) / (4 * dx) -
155 (pC[b_x] * pC[b_y]) / (2 * dx) + (pC[b_0] * pC[b_xz]) / (2 * dx) + (pC[b_0] * pC[b_xyz]) / (4 * dx) -
156 (pC[b_0] * pC[b_xyy]) / (6 * dx) - (pC[b_xxy] * pC[b_xy]) / (24 * dx) - (pC[b_xx] * pC[b_xy]) / (12 * dx) -
157 (pC[b_0] * pC[b_xy]) / (2 * dx) - (pC[b_x] * pC[b_xxy]) / (12 * dx) - (pC[b_x] * pC[b_xx]) / (6 * dx) -
158 (pC[b_0] * pC[b_x]) / dx;
159}
160
173template<typename REAL> inline
175 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
176 creal BGBY,
177 creal BGBZ,
178 const std::array<Real, 3>& gridSpacing
179) {
180 using namespace Rec;
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) -
185 (pC[c_xzz] * BGBZ) / (6 * dx) - (pC[c_xz] * BGBZ) / (2 * dx) + (pC[c_xyz] * BGBZ) / (4 * dx) +
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 +
189 (pC[a_zz] * pC[c_zz]) / (6 * dz) + (pC[a_z] * pC[c_zz]) / (6 * dz) - (pC[a_yz] * pC[c_zz]) / (12 * dz) +
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) -
191 (pC[a_zz] * pC[c_yz]) / (4 * dz) - (pC[a_z] * pC[c_yz]) / (4 * dz) + (pC[a_yz] * pC[c_yz]) / (8 * 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) +
193 (pC[a_xzz] * pC[c_xz]) / (24 * dz) + (pC[a_xz] * pC[c_xz]) / (24 * dz) - (pC[a_xyz] * pC[c_xz]) / (48 * dz) +
194 (pC[a_xzz] * pC[c_x]) / (12 * dz) + (pC[a_xz] * pC[c_x]) / (12 * dz) - (pC[a_xyz] * pC[c_x]) / (24 * dz) +
195 (pC[a_zz] * pC[c_0]) / dz + (pC[a_z] * pC[c_0]) / dz - (pC[a_yz] * pC[c_0]) / (2 * 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) -
197 (pC[a_yz] * pC[b_yz]) / (8 * dy) + (pC[a_yy] * pC[b_yz]) / (4 * dy) - (pC[a_y] * pC[b_yz]) / (4 * dy) +
198 (pC[a_yz] * pC[b_yy]) / (12 * dy) - (pC[a_yy] * pC[b_yy]) / (6 * dy) + (pC[a_y] * pC[b_yy]) / (6 * 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) -
200 (pC[a_xyz] * pC[b_xy]) / (48 * dy) + (pC[a_xyy] * pC[b_xy]) / (24 * dy) - (pC[a_xy] * pC[b_xy]) / (24 * dy) +
201 (pC[a_xyz] * pC[b_x]) / (24 * dy) - (pC[a_xyy] * pC[b_x]) / (12 * dy) + (pC[a_xy] * pC[b_x]) / (12 * dy) +
202 (pC[a_yz] * pC[b_0]) / (2 * dy) - (pC[a_yy] * pC[b_0]) / dy + (pC[a_y] * pC[b_0]) / dy -
203 (pC[c_xzz] * pC[c_zz]) / (36 * dx) - (pC[c_xz] * pC[c_zz]) / (12 * dx) + (pC[c_xyz] * pC[c_zz]) / (24 * dx) +
204 (pC[c_xy] * pC[c_zz]) / (12 * dx) - (pC[c_x] * pC[c_zz]) / (6 * dx) - (pC[c_xzz] * pC[c_z]) / (12 * dx) -
205 (pC[c_xz] * pC[c_z]) / (4 * dx) + (pC[c_xyz] * pC[c_z]) / (8 * dx) + (pC[c_xy] * pC[c_z]) / (4 * dx) -
206 (pC[c_x] * pC[c_z]) / (2 * dx) + (pC[c_xzz] * pC[c_yz]) / (24 * dx) + (pC[c_xz] * pC[c_yz]) / (8 * dx) -
207 (pC[c_xyz] * pC[c_yz]) / (16 * dx) - (pC[c_xy] * pC[c_yz]) / (8 * dx) + (pC[c_x] * pC[c_yz]) / (4 * dx) +
208 (pC[c_xzz] * pC[c_y]) / (12 * dx) + (pC[c_xz] * pC[c_y]) / (4 * dx) - (pC[c_xyz] * pC[c_y]) / (8 * dx) -
209 (pC[c_xy] * pC[c_y]) / (4 * dx) + (pC[c_x] * pC[c_y]) / (2 * dx) - (pC[c_0] * pC[c_xzz]) / (6 * dx) -
210 (pC[c_xxz] * pC[c_xz]) / (24 * dx) - (pC[c_xx] * pC[c_xz]) / (12 * dx) - (pC[c_0] * pC[c_xz]) / (2 * dx) +
211 (pC[c_0] * pC[c_xyz]) / (4 * dx) + (pC[c_0] * pC[c_xy]) / (2 * dx) - (pC[c_x] * pC[c_xxz]) / (12 * dx) -
212 (pC[c_x] * pC[c_xx]) / (6 * dx) - (pC[c_0] * pC[c_x]) / dx - (pC[b_xz] * pC[b_z]) / (4 * dx) +
213 (pC[b_xyz] * pC[b_z]) / (8 * dx) - (pC[b_xyy] * pC[b_z]) / (12 * dx) + (pC[b_xy] * pC[b_z]) / (4 * dx) -
214 (pC[b_x] * pC[b_z]) / (2 * dx) + (pC[b_xz] * pC[b_yz]) / (8 * dx) - (pC[b_xyz] * pC[b_yz]) / (16 * dx) +
215 (pC[b_xyy] * pC[b_yz]) / (24 * dx) - (pC[b_xy] * pC[b_yz]) / (8 * dx) + (pC[b_x] * pC[b_yz]) / (4 * dx) -
216 (pC[b_xz] * pC[b_yy]) / (12 * dx) + (pC[b_xyz] * pC[b_yy]) / (24 * dx) - (pC[b_xyy] * pC[b_yy]) / (36 * dx) +
217 (pC[b_xy] * pC[b_yy]) / (12 * dx) - (pC[b_x] * pC[b_yy]) / (6 * dx) + (pC[b_xz] * pC[b_y]) / (4 * dx) -
218 (pC[b_xyz] * pC[b_y]) / (8 * dx) + (pC[b_xyy] * pC[b_y]) / (12 * dx) - (pC[b_xy] * pC[b_y]) / (4 * dx) +
219 (pC[b_x] * pC[b_y]) / (2 * dx) - (pC[b_0] * pC[b_xz]) / (2 * dx) + (pC[b_0] * pC[b_xyz]) / (4 * dx) -
220 (pC[b_0] * pC[b_xyy]) / (6 * dx) - (pC[b_xxy] * pC[b_xy]) / (24 * dx) + (pC[b_xx] * pC[b_xy]) / (12 * dx) +
221 (pC[b_0] * pC[b_xy]) / (2 * dx) + (pC[b_x] * pC[b_xxy]) / (12 * dx) - (pC[b_x] * pC[b_xx]) / (6 * dx) -
222 (pC[b_0] * pC[b_x]) / dx;
223}
224
237template<typename REAL> inline
239 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
240 creal BGBY,
241 creal BGBZ,
242 const std::array<Real, 3>& gridSpacing
243) {
244 using namespace Rec;
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) -
249 (pC[c_xzz] * BGBZ) / (6 * dx) - (pC[c_xz] * BGBZ) / (2 * dx) - (pC[c_xyz] * BGBZ) / (4 * dx) -
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 +
253 (pC[a_zz] * pC[c_zz]) / (6 * dz) + (pC[a_z] * pC[c_zz]) / (6 * dz) + (pC[a_yz] * pC[c_zz]) / (12 * dz) +
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) +
255 (pC[a_zz] * pC[c_yz]) / (4 * dz) + (pC[a_z] * pC[c_yz]) / (4 * dz) + (pC[a_yz] * pC[c_yz]) / (8 * 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) +
257 (pC[a_xzz] * pC[c_xz]) / (24 * dz) + (pC[a_xz] * pC[c_xz]) / (24 * dz) + (pC[a_xyz] * pC[c_xz]) / (48 * dz) +
258 (pC[a_xzz] * pC[c_x]) / (12 * dz) + (pC[a_xz] * pC[c_x]) / (12 * dz) + (pC[a_xyz] * pC[c_x]) / (24 * dz) +
259 (pC[a_zz] * pC[c_0]) / dz + (pC[a_z] * pC[c_0]) / dz + (pC[a_yz] * pC[c_0]) / (2 * 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) +
261 (pC[a_yz] * pC[b_yz]) / (8 * dy) + (pC[a_yy] * pC[b_yz]) / (4 * dy) + (pC[a_y] * pC[b_yz]) / (4 * dy) +
262 (pC[a_yz] * pC[b_yy]) / (12 * dy) + (pC[a_yy] * pC[b_yy]) / (6 * dy) + (pC[a_y] * pC[b_yy]) / (6 * 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) +
264 (pC[a_xyz] * pC[b_xy]) / (48 * dy) + (pC[a_xyy] * pC[b_xy]) / (24 * dy) + (pC[a_xy] * pC[b_xy]) / (24 * dy) +
265 (pC[a_xyz] * pC[b_x]) / (24 * dy) + (pC[a_xyy] * pC[b_x]) / (12 * dy) + (pC[a_xy] * pC[b_x]) / (12 * dy) +
266 (pC[a_yz] * pC[b_0]) / (2 * dy) + (pC[a_yy] * pC[b_0]) / dy + (pC[a_y] * pC[b_0]) / dy -
267 (pC[c_xzz] * pC[c_zz]) / (36 * dx) - (pC[c_xz] * pC[c_zz]) / (12 * dx) - (pC[c_xyz] * pC[c_zz]) / (24 * dx) -
268 (pC[c_xy] * pC[c_zz]) / (12 * dx) - (pC[c_x] * pC[c_zz]) / (6 * dx) - (pC[c_xzz] * pC[c_z]) / (12 * dx) -
269 (pC[c_xz] * pC[c_z]) / (4 * dx) - (pC[c_xyz] * pC[c_z]) / (8 * dx) - (pC[c_xy] * pC[c_z]) / (4 * dx) -
270 (pC[c_x] * pC[c_z]) / (2 * dx) - (pC[c_xzz] * pC[c_yz]) / (24 * dx) - (pC[c_xz] * pC[c_yz]) / (8 * dx) -
271 (pC[c_xyz] * pC[c_yz]) / (16 * dx) - (pC[c_xy] * pC[c_yz]) / (8 * dx) - (pC[c_x] * pC[c_yz]) / (4 * dx) -
272 (pC[c_xzz] * pC[c_y]) / (12 * dx) - (pC[c_xz] * pC[c_y]) / (4 * dx) - (pC[c_xyz] * pC[c_y]) / (8 * dx) -
273 (pC[c_xy] * pC[c_y]) / (4 * dx) - (pC[c_x] * pC[c_y]) / (2 * dx) - (pC[c_0] * pC[c_xzz]) / (6 * dx) -
274 (pC[c_xxz] * pC[c_xz]) / (24 * dx) - (pC[c_xx] * pC[c_xz]) / (12 * dx) - (pC[c_0] * pC[c_xz]) / (2 * dx) -
275 (pC[c_0] * pC[c_xyz]) / (4 * dx) - (pC[c_0] * pC[c_xy]) / (2 * dx) - (pC[c_x] * pC[c_xxz]) / (12 * dx) -
276 (pC[c_x] * pC[c_xx]) / (6 * dx) - (pC[c_0] * pC[c_x]) / dx - (pC[b_xz] * pC[b_z]) / (4 * dx) -
277 (pC[b_xyz] * pC[b_z]) / (8 * dx) - (pC[b_xyy] * pC[b_z]) / (12 * dx) - (pC[b_xy] * pC[b_z]) / (4 * dx) -
278 (pC[b_x] * pC[b_z]) / (2 * dx) - (pC[b_xz] * pC[b_yz]) / (8 * dx) - (pC[b_xyz] * pC[b_yz]) / (16 * dx) -
279 (pC[b_xyy] * pC[b_yz]) / (24 * dx) - (pC[b_xy] * pC[b_yz]) / (8 * dx) - (pC[b_x] * pC[b_yz]) / (4 * dx) -
280 (pC[b_xz] * pC[b_yy]) / (12 * dx) - (pC[b_xyz] * pC[b_yy]) / (24 * dx) - (pC[b_xyy] * pC[b_yy]) / (36 * dx) -
281 (pC[b_xy] * pC[b_yy]) / (12 * dx) - (pC[b_x] * pC[b_yy]) / (6 * dx) - (pC[b_xz] * pC[b_y]) / (4 * dx) -
282 (pC[b_xyz] * pC[b_y]) / (8 * dx) - (pC[b_xyy] * pC[b_y]) / (12 * dx) - (pC[b_xy] * pC[b_y]) / (4 * dx) -
283 (pC[b_x] * pC[b_y]) / (2 * dx) - (pC[b_0] * pC[b_xz]) / (2 * dx) - (pC[b_0] * pC[b_xyz]) / (4 * dx) -
284 (pC[b_0] * pC[b_xyy]) / (6 * dx) - (pC[b_xxy] * pC[b_xy]) / (24 * dx) - (pC[b_xx] * pC[b_xy]) / (12 * dx) -
285 (pC[b_0] * pC[b_xy]) / (2 * dx) - (pC[b_x] * pC[b_xxy]) / (12 * dx) - (pC[b_x] * pC[b_xx]) / (6 * dx) -
286 (pC[b_0] * pC[b_x]) / dx;
287}
288
289// Y
302template<typename REAL> inline
304 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
305 creal BGBX,
306 creal BGBZ,
307 const std::array<Real, 3>& gridSpacing
308) {
309 using namespace Rec;
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 +
318 (pC[b_x] * BGBX) / dx - (pC[b_zz] * pC[c_zz]) / (6 * dz) + (pC[b_z] * pC[c_zz]) / (6 * dz) -
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) +
320 (pC[b_xz] * pC[c_z]) / (4 * dz) + (pC[b_yzz] * pC[c_yz]) / (24 * dz) - (pC[b_yz] * pC[c_yz]) / (24 * dz) +
321 (pC[b_xyz] * pC[c_yz]) / (48 * dz) - (pC[b_yzz] * pC[c_y]) / (12 * dz) + (pC[b_yz] * pC[c_y]) / (12 * dz) -
322 (pC[b_xyz] * pC[c_y]) / (24 * dz) - (pC[b_zz] * pC[c_xz]) / (4 * dz) + (pC[b_z] * pC[c_xz]) / (4 * 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) +
324 (pC[b_xz] * pC[c_x]) / (4 * dz) - (pC[b_zz] * pC[c_0]) / dz + (pC[b_z] * pC[c_0]) / dz -
325 (pC[b_xz] * pC[c_0]) / (2 * dz) - (pC[c_yzz] * pC[c_zz]) / (36 * dy) + (pC[c_yz] * pC[c_zz]) / (12 * dy) -
326 (pC[c_y] * pC[c_zz]) / (6 * dy) - (pC[c_xyz] * pC[c_zz]) / (24 * dy) + (pC[c_xy] * pC[c_zz]) / (12 * dy) +
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) +
328 (pC[c_xyz] * pC[c_z]) / (8 * dy) - (pC[c_xy] * pC[c_z]) / (4 * dy) - (pC[c_xz] * pC[c_yzz]) / (24 * dy) +
329 (pC[c_x] * pC[c_yzz]) / (12 * dy) - (pC[c_0] * pC[c_yzz]) / (6 * dy) - (pC[c_yyz] * pC[c_yz]) / (24 * dy) +
330 (pC[c_yy] * pC[c_yz]) / (12 * dy) + (pC[c_xz] * pC[c_yz]) / (8 * dy) - (pC[c_x] * pC[c_yz]) / (4 * dy) +
331 (pC[c_0] * pC[c_yz]) / (2 * dy) + (pC[c_y] * pC[c_yyz]) / (12 * dy) - (pC[c_y] * pC[c_yy]) / (6 * 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 -
333 (pC[c_xyz] * pC[c_xz]) / (16 * dy) + (pC[c_xy] * pC[c_xz]) / (8 * dy) + (pC[c_x] * pC[c_xyz]) / (8 * dy) -
334 (pC[c_0] * pC[c_xyz]) / (4 * dy) - (pC[c_x] * pC[c_xy]) / (4 * dy) + (pC[c_0] * pC[c_xy]) / (2 * 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) -
336 (pC[a_xy] * pC[a_z]) / (4 * dy) + (pC[a_xxy] * pC[a_z]) / (12 * dy) + (pC[a_xz] * pC[a_yz]) / (8 * dy) +
337 (pC[a_xx] * pC[a_yz]) / (12 * dy) - (pC[a_x] * pC[a_yz]) / (4 * dy) + (pC[a_0] * pC[a_yz]) / (2 * dy) -
338 (pC[a_y] * pC[a_yy]) / (6 * dy) + (pC[a_xy] * pC[a_yy]) / (12 * dy) - (pC[a_xz] * pC[a_y]) / (4 * 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) -
340 (pC[a_0] * pC[a_y]) / dy - (pC[a_xyz] * pC[a_xz]) / (16 * dy) + (pC[a_xy] * pC[a_xz]) / (8 * dy) -
341 (pC[a_xxy] * pC[a_xz]) / (24 * dy) - (pC[a_xx] * pC[a_xyz]) / (24 * dy) + (pC[a_x] * pC[a_xyz]) / (8 * dy) -
342 (pC[a_0] * pC[a_xyz]) / (4 * dy) - (pC[a_xy] * pC[a_xyy]) / (24 * dy) + (pC[a_xx] * pC[a_xy]) / (12 * dy) -
343 (pC[a_x] * pC[a_xy]) / (4 * dy) + (pC[a_0] * pC[a_xy]) / (2 * dy) - (pC[a_xx] * pC[a_xxy]) / (36 * dy) +
344 (pC[a_x] * pC[a_xxy]) / (12 * dy) - (pC[a_0] * pC[a_xxy]) / (6 * dy) + (pC[a_z] * pC[b_xz]) / (4 * dx) -
345 (pC[a_xz] * pC[b_xz]) / (8 * dx) - (pC[a_xx] * pC[b_xz]) / (12 * dx) + (pC[a_x] * pC[b_xz]) / (4 * dx) -
346 (pC[a_0] * pC[b_xz]) / (2 * dx) - (pC[a_y] * pC[b_xyz]) / (24 * dx) + (pC[a_xy] * pC[b_xyz]) / (48 * dx) +
347 (pC[a_y] * pC[b_xy]) / (12 * dx) - (pC[a_xy] * pC[b_xy]) / (24 * dx) - (pC[a_y] * pC[b_xxy]) / (12 * dx) +
348 (pC[a_xy] * pC[b_xxy]) / (24 * dx) + (pC[a_z] * pC[b_xx]) / (2 * dx) - (pC[a_xz] * pC[b_xx]) / (4 * dx) -
349 (pC[a_xx] * pC[b_xx]) / (6 * dx) + (pC[a_x] * pC[b_xx]) / (2 * dx) - (pC[a_0] * pC[b_xx]) / dx -
350 (pC[a_z] * pC[b_x]) / (2 * dx) + (pC[a_xz] * pC[b_x]) / (4 * dx) + (pC[a_xx] * pC[b_x]) / (6 * dx) -
351 (pC[a_x] * pC[b_x]) / (2 * dx) + (pC[a_0] * pC[b_x]) / dx;
352}
353
366template<typename REAL> inline
368 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
369 creal BGBX,
370 creal BGBZ,
371 const std::array<Real, 3>& gridSpacing
372) {
373 using namespace Rec;
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 +
382 (pC[b_x] * BGBX) / dx - (pC[b_zz] * pC[c_zz]) / (6 * dz) + (pC[b_z] * pC[c_zz]) / (6 * dz) +
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) -
384 (pC[b_xz] * pC[c_z]) / (4 * dz) + (pC[b_yzz] * pC[c_yz]) / (24 * dz) - (pC[b_yz] * pC[c_yz]) / (24 * dz) -
385 (pC[b_xyz] * pC[c_yz]) / (48 * dz) - (pC[b_yzz] * pC[c_y]) / (12 * dz) + (pC[b_yz] * pC[c_y]) / (12 * dz) +
386 (pC[b_xyz] * pC[c_y]) / (24 * dz) + (pC[b_zz] * pC[c_xz]) / (4 * dz) - (pC[b_z] * pC[c_xz]) / (4 * 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) +
388 (pC[b_xz] * pC[c_x]) / (4 * dz) - (pC[b_zz] * pC[c_0]) / dz + (pC[b_z] * pC[c_0]) / dz +
389 (pC[b_xz] * pC[c_0]) / (2 * dz) - (pC[c_yzz] * pC[c_zz]) / (36 * dy) + (pC[c_yz] * pC[c_zz]) / (12 * dy) -
390 (pC[c_y] * pC[c_zz]) / (6 * dy) + (pC[c_xyz] * pC[c_zz]) / (24 * dy) - (pC[c_xy] * pC[c_zz]) / (12 * dy) +
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) -
392 (pC[c_xyz] * pC[c_z]) / (8 * dy) + (pC[c_xy] * pC[c_z]) / (4 * dy) + (pC[c_xz] * pC[c_yzz]) / (24 * dy) -
393 (pC[c_x] * pC[c_yzz]) / (12 * dy) - (pC[c_0] * pC[c_yzz]) / (6 * dy) - (pC[c_yyz] * pC[c_yz]) / (24 * dy) +
394 (pC[c_yy] * pC[c_yz]) / (12 * dy) - (pC[c_xz] * pC[c_yz]) / (8 * dy) + (pC[c_x] * pC[c_yz]) / (4 * dy) +
395 (pC[c_0] * pC[c_yz]) / (2 * dy) + (pC[c_y] * pC[c_yyz]) / (12 * dy) - (pC[c_y] * pC[c_yy]) / (6 * 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 -
397 (pC[c_xyz] * pC[c_xz]) / (16 * dy) + (pC[c_xy] * pC[c_xz]) / (8 * dy) + (pC[c_x] * pC[c_xyz]) / (8 * dy) +
398 (pC[c_0] * pC[c_xyz]) / (4 * dy) - (pC[c_x] * pC[c_xy]) / (4 * dy) - (pC[c_0] * pC[c_xy]) / (2 * 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) +
400 (pC[a_xy] * pC[a_z]) / (4 * dy) + (pC[a_xxy] * pC[a_z]) / (12 * dy) - (pC[a_xz] * pC[a_yz]) / (8 * dy) +
401 (pC[a_xx] * pC[a_yz]) / (12 * dy) + (pC[a_x] * pC[a_yz]) / (4 * dy) + (pC[a_0] * pC[a_yz]) / (2 * dy) -
402 (pC[a_y] * pC[a_yy]) / (6 * dy) - (pC[a_xy] * pC[a_yy]) / (12 * dy) + (pC[a_xz] * pC[a_y]) / (4 * 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) -
404 (pC[a_0] * pC[a_y]) / dy - (pC[a_xyz] * pC[a_xz]) / (16 * dy) + (pC[a_xy] * pC[a_xz]) / (8 * dy) +
405 (pC[a_xxy] * pC[a_xz]) / (24 * dy) + (pC[a_xx] * pC[a_xyz]) / (24 * dy) + (pC[a_x] * pC[a_xyz]) / (8 * dy) +
406 (pC[a_0] * pC[a_xyz]) / (4 * dy) - (pC[a_xy] * pC[a_xyy]) / (24 * dy) - (pC[a_xx] * pC[a_xy]) / (12 * dy) -
407 (pC[a_x] * pC[a_xy]) / (4 * dy) - (pC[a_0] * pC[a_xy]) / (2 * dy) - (pC[a_xx] * pC[a_xxy]) / (36 * dy) -
408 (pC[a_x] * pC[a_xxy]) / (12 * dy) - (pC[a_0] * pC[a_xxy]) / (6 * dy) + (pC[a_z] * pC[b_xz]) / (4 * dx) +
409 (pC[a_xz] * pC[b_xz]) / (8 * dx) - (pC[a_xx] * pC[b_xz]) / (12 * dx) - (pC[a_x] * pC[b_xz]) / (4 * dx) -
410 (pC[a_0] * pC[b_xz]) / (2 * dx) - (pC[a_y] * pC[b_xyz]) / (24 * dx) - (pC[a_xy] * pC[b_xyz]) / (48 * dx) +
411 (pC[a_y] * pC[b_xy]) / (12 * dx) + (pC[a_xy] * pC[b_xy]) / (24 * dx) + (pC[a_y] * pC[b_xxy]) / (12 * dx) +
412 (pC[a_xy] * pC[b_xxy]) / (24 * dx) - (pC[a_z] * pC[b_xx]) / (2 * dx) - (pC[a_xz] * pC[b_xx]) / (4 * dx) +
413 (pC[a_xx] * pC[b_xx]) / (6 * dx) + (pC[a_x] * pC[b_xx]) / (2 * dx) + (pC[a_0] * pC[b_xx]) / dx -
414 (pC[a_z] * pC[b_x]) / (2 * dx) - (pC[a_xz] * pC[b_x]) / (4 * dx) + (pC[a_xx] * pC[b_x]) / (6 * dx) +
415 (pC[a_x] * pC[b_x]) / (2 * dx) + (pC[a_0] * pC[b_x]) / dx;
416}
417
430template<typename REAL> inline
432 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
433 creal BGBX,
434 creal BGBZ,
435 const std::array<Real, 3>& gridSpacing
436) {
437 using namespace Rec;
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 +
446 (pC[b_x] * BGBX) / dx + (pC[b_zz] * pC[c_zz]) / (6 * dz) + (pC[b_z] * pC[c_zz]) / (6 * dz) -
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) -
448 (pC[b_xz] * pC[c_z]) / (4 * dz) + (pC[b_yzz] * pC[c_yz]) / (24 * dz) + (pC[b_yz] * pC[c_yz]) / (24 * dz) -
449 (pC[b_xyz] * pC[c_yz]) / (48 * dz) + (pC[b_yzz] * pC[c_y]) / (12 * dz) + (pC[b_yz] * pC[c_y]) / (12 * dz) -
450 (pC[b_xyz] * pC[c_y]) / (24 * dz) - (pC[b_zz] * pC[c_xz]) / (4 * dz) - (pC[b_z] * pC[c_xz]) / (4 * 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) +
452 (pC[b_xz] * pC[c_x]) / (4 * dz) + (pC[b_zz] * pC[c_0]) / dz + (pC[b_z] * pC[c_0]) / dz -
453 (pC[b_xz] * pC[c_0]) / (2 * dz) - (pC[c_yzz] * pC[c_zz]) / (36 * dy) - (pC[c_yz] * pC[c_zz]) / (12 * dy) -
454 (pC[c_y] * pC[c_zz]) / (6 * dy) + (pC[c_xyz] * pC[c_zz]) / (24 * dy) + (pC[c_xy] * pC[c_zz]) / (12 * dy) -
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) +
456 (pC[c_xyz] * pC[c_z]) / (8 * dy) + (pC[c_xy] * pC[c_z]) / (4 * dy) + (pC[c_xz] * pC[c_yzz]) / (24 * dy) +
457 (pC[c_x] * pC[c_yzz]) / (12 * dy) - (pC[c_0] * pC[c_yzz]) / (6 * dy) - (pC[c_yyz] * pC[c_yz]) / (24 * dy) -
458 (pC[c_yy] * pC[c_yz]) / (12 * dy) + (pC[c_xz] * pC[c_yz]) / (8 * dy) + (pC[c_x] * pC[c_yz]) / (4 * dy) -
459 (pC[c_0] * pC[c_yz]) / (2 * dy) - (pC[c_y] * pC[c_yyz]) / (12 * dy) - (pC[c_y] * pC[c_yy]) / (6 * 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 -
461 (pC[c_xyz] * pC[c_xz]) / (16 * dy) - (pC[c_xy] * pC[c_xz]) / (8 * dy) - (pC[c_x] * pC[c_xyz]) / (8 * dy) +
462 (pC[c_0] * pC[c_xyz]) / (4 * dy) - (pC[c_x] * pC[c_xy]) / (4 * dy) + (pC[c_0] * pC[c_xy]) / (2 * 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) +
464 (pC[a_xy] * pC[a_z]) / (4 * dy) - (pC[a_xxy] * pC[a_z]) / (12 * dy) + (pC[a_xz] * pC[a_yz]) / (8 * dy) -
465 (pC[a_xx] * pC[a_yz]) / (12 * dy) + (pC[a_x] * pC[a_yz]) / (4 * dy) - (pC[a_0] * pC[a_yz]) / (2 * dy) -
466 (pC[a_y] * pC[a_yy]) / (6 * dy) + (pC[a_xy] * pC[a_yy]) / (12 * dy) + (pC[a_xz] * pC[a_y]) / (4 * 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) -
468 (pC[a_0] * pC[a_y]) / dy - (pC[a_xyz] * pC[a_xz]) / (16 * dy) - (pC[a_xy] * pC[a_xz]) / (8 * dy) +
469 (pC[a_xxy] * pC[a_xz]) / (24 * dy) + (pC[a_xx] * pC[a_xyz]) / (24 * dy) - (pC[a_x] * pC[a_xyz]) / (8 * dy) +
470 (pC[a_0] * pC[a_xyz]) / (4 * dy) - (pC[a_xy] * pC[a_xyy]) / (24 * dy) + (pC[a_xx] * pC[a_xy]) / (12 * dy) -
471 (pC[a_x] * pC[a_xy]) / (4 * dy) + (pC[a_0] * pC[a_xy]) / (2 * dy) - (pC[a_xx] * pC[a_xxy]) / (36 * dy) +
472 (pC[a_x] * pC[a_xxy]) / (12 * dy) - (pC[a_0] * pC[a_xxy]) / (6 * dy) + (pC[a_z] * pC[b_xz]) / (4 * dx) -
473 (pC[a_xz] * pC[b_xz]) / (8 * dx) + (pC[a_xx] * pC[b_xz]) / (12 * dx) - (pC[a_x] * pC[b_xz]) / (4 * dx) +
474 (pC[a_0] * pC[b_xz]) / (2 * dx) + (pC[a_y] * pC[b_xyz]) / (24 * dx) - (pC[a_xy] * pC[b_xyz]) / (48 * dx) +
475 (pC[a_y] * pC[b_xy]) / (12 * dx) - (pC[a_xy] * pC[b_xy]) / (24 * dx) - (pC[a_y] * pC[b_xxy]) / (12 * dx) +
476 (pC[a_xy] * pC[b_xxy]) / (24 * dx) - (pC[a_z] * pC[b_xx]) / (2 * dx) + (pC[a_xz] * pC[b_xx]) / (4 * dx) -
477 (pC[a_xx] * pC[b_xx]) / (6 * dx) + (pC[a_x] * pC[b_xx]) / (2 * dx) - (pC[a_0] * pC[b_xx]) / dx +
478 (pC[a_z] * pC[b_x]) / (2 * dx) - (pC[a_xz] * pC[b_x]) / (4 * dx) + (pC[a_xx] * pC[b_x]) / (6 * dx) -
479 (pC[a_x] * pC[b_x]) / (2 * dx) + (pC[a_0] * pC[b_x]) / dx;
480}
481
494template<typename REAL> inline
496 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
497 creal BGBX,
498 creal BGBZ,
499 const std::array<Real, 3>& gridSpacing
500) {
501 using namespace Rec;
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 +
510 (pC[b_x] * BGBX) / dx + (pC[b_zz] * pC[c_zz]) / (6 * dz) + (pC[b_z] * pC[c_zz]) / (6 * dz) +
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) +
512 (pC[b_xz] * pC[c_z]) / (4 * dz) + (pC[b_yzz] * pC[c_yz]) / (24 * dz) + (pC[b_yz] * pC[c_yz]) / (24 * dz) +
513 (pC[b_xyz] * pC[c_yz]) / (48 * dz) + (pC[b_yzz] * pC[c_y]) / (12 * dz) + (pC[b_yz] * pC[c_y]) / (12 * dz) +
514 (pC[b_xyz] * pC[c_y]) / (24 * dz) + (pC[b_zz] * pC[c_xz]) / (4 * dz) + (pC[b_z] * pC[c_xz]) / (4 * 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) +
516 (pC[b_xz] * pC[c_x]) / (4 * dz) + (pC[b_zz] * pC[c_0]) / dz + (pC[b_z] * pC[c_0]) / dz +
517 (pC[b_xz] * pC[c_0]) / (2 * dz) - (pC[c_yzz] * pC[c_zz]) / (36 * dy) - (pC[c_yz] * pC[c_zz]) / (12 * dy) -
518 (pC[c_y] * pC[c_zz]) / (6 * dy) - (pC[c_xyz] * pC[c_zz]) / (24 * dy) - (pC[c_xy] * pC[c_zz]) / (12 * dy) -
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) -
520 (pC[c_xyz] * pC[c_z]) / (8 * dy) - (pC[c_xy] * pC[c_z]) / (4 * dy) - (pC[c_xz] * pC[c_yzz]) / (24 * dy) -
521 (pC[c_x] * pC[c_yzz]) / (12 * dy) - (pC[c_0] * pC[c_yzz]) / (6 * dy) - (pC[c_yyz] * pC[c_yz]) / (24 * dy) -
522 (pC[c_yy] * pC[c_yz]) / (12 * dy) - (pC[c_xz] * pC[c_yz]) / (8 * dy) - (pC[c_x] * pC[c_yz]) / (4 * dy) -
523 (pC[c_0] * pC[c_yz]) / (2 * dy) - (pC[c_y] * pC[c_yyz]) / (12 * dy) - (pC[c_y] * pC[c_yy]) / (6 * 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 -
525 (pC[c_xyz] * pC[c_xz]) / (16 * dy) - (pC[c_xy] * pC[c_xz]) / (8 * dy) - (pC[c_x] * pC[c_xyz]) / (8 * dy) -
526 (pC[c_0] * pC[c_xyz]) / (4 * dy) - (pC[c_x] * pC[c_xy]) / (4 * dy) - (pC[c_0] * pC[c_xy]) / (2 * 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) -
528 (pC[a_xy] * pC[a_z]) / (4 * dy) - (pC[a_xxy] * pC[a_z]) / (12 * dy) - (pC[a_xz] * pC[a_yz]) / (8 * dy) -
529 (pC[a_xx] * pC[a_yz]) / (12 * dy) - (pC[a_x] * pC[a_yz]) / (4 * dy) - (pC[a_0] * pC[a_yz]) / (2 * dy) -
530 (pC[a_y] * pC[a_yy]) / (6 * dy) - (pC[a_xy] * pC[a_yy]) / (12 * dy) - (pC[a_xz] * pC[a_y]) / (4 * 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) -
532 (pC[a_0] * pC[a_y]) / dy - (pC[a_xyz] * pC[a_xz]) / (16 * dy) - (pC[a_xy] * pC[a_xz]) / (8 * dy) -
533 (pC[a_xxy] * pC[a_xz]) / (24 * dy) - (pC[a_xx] * pC[a_xyz]) / (24 * dy) - (pC[a_x] * pC[a_xyz]) / (8 * dy) -
534 (pC[a_0] * pC[a_xyz]) / (4 * dy) - (pC[a_xy] * pC[a_xyy]) / (24 * dy) - (pC[a_xx] * pC[a_xy]) / (12 * dy) -
535 (pC[a_x] * pC[a_xy]) / (4 * dy) - (pC[a_0] * pC[a_xy]) / (2 * dy) - (pC[a_xx] * pC[a_xxy]) / (36 * dy) -
536 (pC[a_x] * pC[a_xxy]) / (12 * dy) - (pC[a_0] * pC[a_xxy]) / (6 * dy) + (pC[a_z] * pC[b_xz]) / (4 * dx) +
537 (pC[a_xz] * pC[b_xz]) / (8 * dx) + (pC[a_xx] * pC[b_xz]) / (12 * dx) + (pC[a_x] * pC[b_xz]) / (4 * dx) +
538 (pC[a_0] * pC[b_xz]) / (2 * dx) + (pC[a_y] * pC[b_xyz]) / (24 * dx) + (pC[a_xy] * pC[b_xyz]) / (48 * dx) +
539 (pC[a_y] * pC[b_xy]) / (12 * dx) + (pC[a_xy] * pC[b_xy]) / (24 * dx) + (pC[a_y] * pC[b_xxy]) / (12 * dx) +
540 (pC[a_xy] * pC[b_xxy]) / (24 * dx) + (pC[a_z] * pC[b_xx]) / (2 * dx) + (pC[a_xz] * pC[b_xx]) / (4 * dx) +
541 (pC[a_xx] * pC[b_xx]) / (6 * dx) + (pC[a_x] * pC[b_xx]) / (2 * dx) + (pC[a_0] * pC[b_xx]) / dx +
542 (pC[a_z] * pC[b_x]) / (2 * dx) + (pC[a_xz] * pC[b_x]) / (4 * dx) + (pC[a_xx] * pC[b_x]) / (6 * dx) +
543 (pC[a_x] * pC[b_x]) / (2 * dx) + (pC[a_0] * pC[b_x]) / dx;
544}
545
546// Z
559template<typename REAL> inline
561 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
562 creal BGBX,
563 creal BGBY,
564 const std::array<Real, 3>& gridSpacing
565) {
566 using namespace Rec;
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) -
574 (pC[c_xy] * BGBX) / (2 * dx) - (pC[c_xx] * BGBX) / dx + (pC[c_x] * BGBX) / dx -
575 (pC[b_z] * pC[b_zz]) / (6 * dz) + (pC[b_yz] * pC[b_zz]) / (12 * dz) + (pC[b_yzz] * pC[b_z]) / (12 * 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) +
577 (pC[b_x] * pC[b_z]) / (2 * dz) - (pC[b_0] * pC[b_z]) / dz - (pC[b_yz] * pC[b_yzz]) / (24 * dz) +
578 (pC[b_yy] * pC[b_yz]) / (12 * dz) - (pC[b_y] * pC[b_yz]) / (4 * dz) + (pC[b_xy] * pC[b_yz]) / (8 * dz) -
579 (pC[b_x] * pC[b_yz]) / (4 * dz) + (pC[b_0] * pC[b_yz]) / (2 * dz) - (pC[b_yy] * pC[b_yyz]) / (36 * dz) +
580 (pC[b_y] * pC[b_yyz]) / (12 * dz) - (pC[b_xy] * pC[b_yyz]) / (24 * dz) + (pC[b_x] * pC[b_yyz]) / (12 * dz) -
581 (pC[b_0] * pC[b_yyz]) / (6 * dz) + (pC[b_xz] * pC[b_yy]) / (12 * dz) - (pC[b_xyz] * pC[b_yy]) / (24 * dz) -
582 (pC[b_xz] * pC[b_y]) / (4 * dz) + (pC[b_xyz] * pC[b_y]) / (8 * dz) + (pC[b_xy] * pC[b_xz]) / (8 * dz) -
583 (pC[b_x] * pC[b_xz]) / (4 * dz) + (pC[b_0] * pC[b_xz]) / (2 * dz) - (pC[b_xy] * pC[b_xyz]) / (16 * dz) +
584 (pC[b_x] * pC[b_xyz]) / (8 * dz) - (pC[b_0] * pC[b_xyz]) / (4 * dz) - (pC[a_z] * pC[a_zz]) / (6 * dz) +
585 (pC[a_xz] * pC[a_zz]) / (12 * dz) + (pC[a_y] * pC[a_z]) / (2 * dz) + (pC[a_xzz] * pC[a_z]) / (12 * 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) -
587 (pC[a_0] * pC[a_z]) / dz - (pC[a_y] * pC[a_yz]) / (4 * dz) + (pC[a_xy] * pC[a_yz]) / (8 * dz) +
588 (pC[a_xx] * pC[a_yz]) / (12 * dz) - (pC[a_x] * pC[a_yz]) / (4 * dz) + (pC[a_0] * pC[a_yz]) / (2 * dz) -
589 (pC[a_xz] * pC[a_y]) / (4 * dz) + (pC[a_xyz] * pC[a_y]) / (8 * dz) + (pC[a_xxz] * pC[a_y]) / (12 * dz) -
590 (pC[a_xz] * pC[a_xzz]) / (24 * dz) + (pC[a_xy] * pC[a_xz]) / (8 * dz) + (pC[a_xx] * pC[a_xz]) / (12 * dz) -
591 (pC[a_x] * pC[a_xz]) / (4 * dz) + (pC[a_0] * pC[a_xz]) / (2 * dz) - (pC[a_xy] * pC[a_xyz]) / (16 * dz) -
592 (pC[a_xx] * pC[a_xyz]) / (24 * dz) + (pC[a_x] * pC[a_xyz]) / (8 * dz) - (pC[a_0] * pC[a_xyz]) / (4 * dz) -
593 (pC[a_xxz] * pC[a_xy]) / (24 * dz) - (pC[a_xx] * pC[a_xxz]) / (36 * dz) + (pC[a_x] * pC[a_xxz]) / (12 * dz) -
594 (pC[a_0] * pC[a_xxz]) / (6 * dz) + (pC[b_z] * pC[c_yz]) / (12 * dy) - (pC[b_yz] * pC[c_yz]) / (24 * dy) -
595 (pC[b_z] * pC[c_yyz]) / (12 * dy) + (pC[b_yz] * pC[c_yyz]) / (24 * dy) - (pC[b_yy] * pC[c_yy]) / (6 * dy) +
596 (pC[b_y] * pC[c_yy]) / (2 * dy) - (pC[b_xy] * pC[c_yy]) / (4 * dy) + (pC[b_x] * pC[c_yy]) / (2 * dy) -
597 (pC[b_0] * pC[c_yy]) / dy + (pC[b_yy] * pC[c_y]) / (6 * dy) - (pC[b_y] * pC[c_y]) / (2 * dy) +
598 (pC[b_xy] * pC[c_y]) / (4 * dy) - (pC[b_x] * pC[c_y]) / (2 * dy) + (pC[b_0] * pC[c_y]) / dy -
599 (pC[b_z] * pC[c_xyz]) / (24 * dy) + (pC[b_yz] * pC[c_xyz]) / (48 * dy) - (pC[b_yy] * pC[c_xy]) / (12 * dy) +
600 (pC[b_y] * pC[c_xy]) / (4 * dy) - (pC[b_xy] * pC[c_xy]) / (8 * dy) + (pC[b_x] * pC[c_xy]) / (4 * dy) -
601 (pC[b_0] * pC[c_xy]) / (2 * dy) + (pC[a_z] * pC[c_xz]) / (12 * dx) - (pC[a_xz] * pC[c_xz]) / (24 * dx) -
602 (pC[a_z] * pC[c_xyz]) / (24 * dx) + (pC[a_xz] * pC[c_xyz]) / (48 * dx) + (pC[a_y] * pC[c_xy]) / (4 * dx) -
603 (pC[a_xy] * pC[c_xy]) / (8 * dx) - (pC[a_xx] * pC[c_xy]) / (12 * dx) + (pC[a_x] * pC[c_xy]) / (4 * dx) -
604 (pC[a_0] * pC[c_xy]) / (2 * dx) - (pC[a_z] * pC[c_xxz]) / (12 * dx) + (pC[a_xz] * pC[c_xxz]) / (24 * dx) +
605 (pC[a_y] * pC[c_xx]) / (2 * dx) - (pC[a_xy] * pC[c_xx]) / (4 * dx) - (pC[a_xx] * pC[c_xx]) / (6 * dx) +
606 (pC[a_x] * pC[c_xx]) / (2 * dx) - (pC[a_0] * pC[c_xx]) / dx - (pC[a_y] * pC[c_x]) / (2 * dx) +
607 (pC[a_xy] * pC[c_x]) / (4 * dx) + (pC[a_xx] * pC[c_x]) / (6 * dx) - (pC[a_x] * pC[c_x]) / (2 * dx) +
608 (pC[a_0] * pC[c_x]) / dx;
609}
610
623template<typename REAL> inline
625 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
626 creal BGBX,
627 creal BGBY,
628 const std::array<Real, 3>& gridSpacing
629) {
630 using namespace Rec;
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) -
638 (pC[c_xy] * BGBX) / (2 * dx) + (pC[c_xx] * BGBX) / dx + (pC[c_x] * BGBX) / dx -
639 (pC[b_z] * pC[b_zz]) / (6 * dz) + (pC[b_yz] * pC[b_zz]) / (12 * dz) + (pC[b_yzz] * pC[b_z]) / (12 * 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) -
641 (pC[b_x] * pC[b_z]) / (2 * dz) - (pC[b_0] * pC[b_z]) / dz - (pC[b_yz] * pC[b_yzz]) / (24 * dz) +
642 (pC[b_yy] * pC[b_yz]) / (12 * dz) - (pC[b_y] * pC[b_yz]) / (4 * dz) - (pC[b_xy] * pC[b_yz]) / (8 * dz) +
643 (pC[b_x] * pC[b_yz]) / (4 * dz) + (pC[b_0] * pC[b_yz]) / (2 * dz) - (pC[b_yy] * pC[b_yyz]) / (36 * dz) +
644 (pC[b_y] * pC[b_yyz]) / (12 * dz) + (pC[b_xy] * pC[b_yyz]) / (24 * dz) - (pC[b_x] * pC[b_yyz]) / (12 * dz) -
645 (pC[b_0] * pC[b_yyz]) / (6 * dz) - (pC[b_xz] * pC[b_yy]) / (12 * dz) + (pC[b_xyz] * pC[b_yy]) / (24 * dz) +
646 (pC[b_xz] * pC[b_y]) / (4 * dz) - (pC[b_xyz] * pC[b_y]) / (8 * dz) + (pC[b_xy] * pC[b_xz]) / (8 * dz) -
647 (pC[b_x] * pC[b_xz]) / (4 * dz) - (pC[b_0] * pC[b_xz]) / (2 * dz) - (pC[b_xy] * pC[b_xyz]) / (16 * dz) +
648 (pC[b_x] * pC[b_xyz]) / (8 * dz) + (pC[b_0] * pC[b_xyz]) / (4 * dz) - (pC[a_z] * pC[a_zz]) / (6 * dz) -
649 (pC[a_xz] * pC[a_zz]) / (12 * dz) + (pC[a_y] * pC[a_z]) / (2 * dz) - (pC[a_xzz] * pC[a_z]) / (12 * 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) -
651 (pC[a_0] * pC[a_z]) / dz - (pC[a_y] * pC[a_yz]) / (4 * dz) - (pC[a_xy] * pC[a_yz]) / (8 * dz) +
652 (pC[a_xx] * pC[a_yz]) / (12 * dz) + (pC[a_x] * pC[a_yz]) / (4 * dz) + (pC[a_0] * pC[a_yz]) / (2 * dz) +
653 (pC[a_xz] * pC[a_y]) / (4 * dz) - (pC[a_xyz] * pC[a_y]) / (8 * dz) + (pC[a_xxz] * pC[a_y]) / (12 * dz) -
654 (pC[a_xz] * pC[a_xzz]) / (24 * dz) + (pC[a_xy] * pC[a_xz]) / (8 * dz) - (pC[a_xx] * pC[a_xz]) / (12 * dz) -
655 (pC[a_x] * pC[a_xz]) / (4 * dz) - (pC[a_0] * pC[a_xz]) / (2 * dz) - (pC[a_xy] * pC[a_xyz]) / (16 * dz) +
656 (pC[a_xx] * pC[a_xyz]) / (24 * dz) + (pC[a_x] * pC[a_xyz]) / (8 * dz) + (pC[a_0] * pC[a_xyz]) / (4 * dz) +
657 (pC[a_xxz] * pC[a_xy]) / (24 * dz) - (pC[a_xx] * pC[a_xxz]) / (36 * dz) - (pC[a_x] * pC[a_xxz]) / (12 * dz) -
658 (pC[a_0] * pC[a_xxz]) / (6 * dz) + (pC[b_z] * pC[c_yz]) / (12 * dy) - (pC[b_yz] * pC[c_yz]) / (24 * dy) -
659 (pC[b_z] * pC[c_yyz]) / (12 * dy) + (pC[b_yz] * pC[c_yyz]) / (24 * dy) - (pC[b_yy] * pC[c_yy]) / (6 * dy) +
660 (pC[b_y] * pC[c_yy]) / (2 * dy) + (pC[b_xy] * pC[c_yy]) / (4 * dy) - (pC[b_x] * pC[c_yy]) / (2 * dy) -
661 (pC[b_0] * pC[c_yy]) / dy + (pC[b_yy] * pC[c_y]) / (6 * dy) - (pC[b_y] * pC[c_y]) / (2 * dy) -
662 (pC[b_xy] * pC[c_y]) / (4 * dy) + (pC[b_x] * pC[c_y]) / (2 * dy) + (pC[b_0] * pC[c_y]) / dy +
663 (pC[b_z] * pC[c_xyz]) / (24 * dy) - (pC[b_yz] * pC[c_xyz]) / (48 * dy) + (pC[b_yy] * pC[c_xy]) / (12 * dy) -
664 (pC[b_y] * pC[c_xy]) / (4 * dy) - (pC[b_xy] * pC[c_xy]) / (8 * dy) + (pC[b_x] * pC[c_xy]) / (4 * dy) +
665 (pC[b_0] * pC[c_xy]) / (2 * dy) + (pC[a_z] * pC[c_xz]) / (12 * dx) + (pC[a_xz] * pC[c_xz]) / (24 * dx) -
666 (pC[a_z] * pC[c_xyz]) / (24 * dx) - (pC[a_xz] * pC[c_xyz]) / (48 * dx) + (pC[a_y] * pC[c_xy]) / (4 * dx) +
667 (pC[a_xy] * pC[c_xy]) / (8 * dx) - (pC[a_xx] * pC[c_xy]) / (12 * dx) - (pC[a_x] * pC[c_xy]) / (4 * dx) -
668 (pC[a_0] * pC[c_xy]) / (2 * dx) + (pC[a_z] * pC[c_xxz]) / (12 * dx) + (pC[a_xz] * pC[c_xxz]) / (24 * dx) -
669 (pC[a_y] * pC[c_xx]) / (2 * dx) - (pC[a_xy] * pC[c_xx]) / (4 * dx) + (pC[a_xx] * pC[c_xx]) / (6 * dx) +
670 (pC[a_x] * pC[c_xx]) / (2 * dx) + (pC[a_0] * pC[c_xx]) / dx - (pC[a_y] * pC[c_x]) / (2 * dx) -
671 (pC[a_xy] * pC[c_x]) / (4 * dx) + (pC[a_xx] * pC[c_x]) / (6 * dx) + (pC[a_x] * pC[c_x]) / (2 * dx) +
672 (pC[a_0] * pC[c_x]) / dx;
673}
674
687template<typename REAL> inline
689 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
690 creal BGBX,
691 creal BGBY,
692 const std::array<Real, 3>& gridSpacing
693) {
694 using namespace Rec;
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) +
702 (pC[c_xy] * BGBX) / (2 * dx) - (pC[c_xx] * BGBX) / dx + (pC[c_x] * BGBX) / dx -
703 (pC[b_z] * pC[b_zz]) / (6 * dz) - (pC[b_yz] * pC[b_zz]) / (12 * dz) - (pC[b_yzz] * pC[b_z]) / (12 * 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) +
705 (pC[b_x] * pC[b_z]) / (2 * dz) - (pC[b_0] * pC[b_z]) / dz - (pC[b_yz] * pC[b_yzz]) / (24 * dz) -
706 (pC[b_yy] * pC[b_yz]) / (12 * dz) - (pC[b_y] * pC[b_yz]) / (4 * dz) + (pC[b_xy] * pC[b_yz]) / (8 * dz) +
707 (pC[b_x] * pC[b_yz]) / (4 * dz) - (pC[b_0] * pC[b_yz]) / (2 * dz) - (pC[b_yy] * pC[b_yyz]) / (36 * dz) -
708 (pC[b_y] * pC[b_yyz]) / (12 * dz) + (pC[b_xy] * pC[b_yyz]) / (24 * dz) + (pC[b_x] * pC[b_yyz]) / (12 * dz) -
709 (pC[b_0] * pC[b_yyz]) / (6 * dz) + (pC[b_xz] * pC[b_yy]) / (12 * dz) + (pC[b_xyz] * pC[b_yy]) / (24 * dz) +
710 (pC[b_xz] * pC[b_y]) / (4 * dz) + (pC[b_xyz] * pC[b_y]) / (8 * dz) - (pC[b_xy] * pC[b_xz]) / (8 * dz) -
711 (pC[b_x] * pC[b_xz]) / (4 * dz) + (pC[b_0] * pC[b_xz]) / (2 * dz) - (pC[b_xy] * pC[b_xyz]) / (16 * dz) -
712 (pC[b_x] * pC[b_xyz]) / (8 * dz) + (pC[b_0] * pC[b_xyz]) / (4 * dz) - (pC[a_z] * pC[a_zz]) / (6 * dz) +
713 (pC[a_xz] * pC[a_zz]) / (12 * dz) - (pC[a_y] * pC[a_z]) / (2 * dz) + (pC[a_xzz] * pC[a_z]) / (12 * 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) -
715 (pC[a_0] * pC[a_z]) / dz - (pC[a_y] * pC[a_yz]) / (4 * dz) + (pC[a_xy] * pC[a_yz]) / (8 * dz) -
716 (pC[a_xx] * pC[a_yz]) / (12 * dz) + (pC[a_x] * pC[a_yz]) / (4 * dz) - (pC[a_0] * pC[a_yz]) / (2 * dz) +
717 (pC[a_xz] * pC[a_y]) / (4 * dz) + (pC[a_xyz] * pC[a_y]) / (8 * dz) - (pC[a_xxz] * pC[a_y]) / (12 * dz) -
718 (pC[a_xz] * pC[a_xzz]) / (24 * dz) - (pC[a_xy] * pC[a_xz]) / (8 * dz) + (pC[a_xx] * pC[a_xz]) / (12 * dz) -
719 (pC[a_x] * pC[a_xz]) / (4 * dz) + (pC[a_0] * pC[a_xz]) / (2 * dz) - (pC[a_xy] * pC[a_xyz]) / (16 * dz) +
720 (pC[a_xx] * pC[a_xyz]) / (24 * dz) - (pC[a_x] * pC[a_xyz]) / (8 * dz) + (pC[a_0] * pC[a_xyz]) / (4 * dz) +
721 (pC[a_xxz] * pC[a_xy]) / (24 * dz) - (pC[a_xx] * pC[a_xxz]) / (36 * dz) + (pC[a_x] * pC[a_xxz]) / (12 * dz) -
722 (pC[a_0] * pC[a_xxz]) / (6 * dz) + (pC[b_z] * pC[c_yz]) / (12 * dy) + (pC[b_yz] * pC[c_yz]) / (24 * dy) +
723 (pC[b_z] * pC[c_yyz]) / (12 * dy) + (pC[b_yz] * pC[c_yyz]) / (24 * dy) + (pC[b_yy] * pC[c_yy]) / (6 * dy) +
724 (pC[b_y] * pC[c_yy]) / (2 * dy) - (pC[b_xy] * pC[c_yy]) / (4 * dy) - (pC[b_x] * pC[c_yy]) / (2 * dy) +
725 (pC[b_0] * pC[c_yy]) / dy + (pC[b_yy] * pC[c_y]) / (6 * dy) + (pC[b_y] * pC[c_y]) / (2 * dy) -
726 (pC[b_xy] * pC[c_y]) / (4 * dy) - (pC[b_x] * pC[c_y]) / (2 * dy) + (pC[b_0] * pC[c_y]) / dy -
727 (pC[b_z] * pC[c_xyz]) / (24 * dy) - (pC[b_yz] * pC[c_xyz]) / (48 * dy) - (pC[b_yy] * pC[c_xy]) / (12 * dy) -
728 (pC[b_y] * pC[c_xy]) / (4 * dy) + (pC[b_xy] * pC[c_xy]) / (8 * dy) + (pC[b_x] * pC[c_xy]) / (4 * dy) -
729 (pC[b_0] * pC[c_xy]) / (2 * dy) + (pC[a_z] * pC[c_xz]) / (12 * dx) - (pC[a_xz] * pC[c_xz]) / (24 * dx) +
730 (pC[a_z] * pC[c_xyz]) / (24 * dx) - (pC[a_xz] * pC[c_xyz]) / (48 * dx) + (pC[a_y] * pC[c_xy]) / (4 * dx) -
731 (pC[a_xy] * pC[c_xy]) / (8 * dx) + (pC[a_xx] * pC[c_xy]) / (12 * dx) - (pC[a_x] * pC[c_xy]) / (4 * dx) +
732 (pC[a_0] * pC[c_xy]) / (2 * dx) - (pC[a_z] * pC[c_xxz]) / (12 * dx) + (pC[a_xz] * pC[c_xxz]) / (24 * dx) -
733 (pC[a_y] * pC[c_xx]) / (2 * dx) + (pC[a_xy] * pC[c_xx]) / (4 * dx) - (pC[a_xx] * pC[c_xx]) / (6 * dx) +
734 (pC[a_x] * pC[c_xx]) / (2 * dx) - (pC[a_0] * pC[c_xx]) / dx + (pC[a_y] * pC[c_x]) / (2 * dx) -
735 (pC[a_xy] * pC[c_x]) / (4 * dx) + (pC[a_xx] * pC[c_x]) / (6 * dx) - (pC[a_x] * pC[c_x]) / (2 * dx) +
736 (pC[a_0] * pC[c_x]) / dx;
737}
738
751template<typename REAL> inline
753 const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC,
754 creal BGBX,
755 creal BGBY,
756 const std::array<Real, 3>& gridSpacing
757) {
758 using namespace Rec;
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) +
766 (pC[c_xy] * BGBX) / (2 * dx) + (pC[c_xx] * BGBX) / dx + (pC[c_x] * BGBX) / dx -
767 (pC[b_z] * pC[b_zz]) / (6 * dz) - (pC[b_yz] * pC[b_zz]) / (12 * dz) - (pC[b_yzz] * pC[b_z]) / (12 * 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) -
769 (pC[b_x] * pC[b_z]) / (2 * dz) - (pC[b_0] * pC[b_z]) / dz - (pC[b_yz] * pC[b_yzz]) / (24 * dz) -
770 (pC[b_yy] * pC[b_yz]) / (12 * dz) - (pC[b_y] * pC[b_yz]) / (4 * dz) - (pC[b_xy] * pC[b_yz]) / (8 * dz) -
771 (pC[b_x] * pC[b_yz]) / (4 * dz) - (pC[b_0] * pC[b_yz]) / (2 * dz) - (pC[b_yy] * pC[b_yyz]) / (36 * dz) -
772 (pC[b_y] * pC[b_yyz]) / (12 * dz) - (pC[b_xy] * pC[b_yyz]) / (24 * dz) - (pC[b_x] * pC[b_yyz]) / (12 * dz) -
773 (pC[b_0] * pC[b_yyz]) / (6 * dz) - (pC[b_xz] * pC[b_yy]) / (12 * dz) - (pC[b_xyz] * pC[b_yy]) / (24 * dz) -
774 (pC[b_xz] * pC[b_y]) / (4 * dz) - (pC[b_xyz] * pC[b_y]) / (8 * dz) - (pC[b_xy] * pC[b_xz]) / (8 * dz) -
775 (pC[b_x] * pC[b_xz]) / (4 * dz) - (pC[b_0] * pC[b_xz]) / (2 * dz) - (pC[b_xy] * pC[b_xyz]) / (16 * dz) -
776 (pC[b_x] * pC[b_xyz]) / (8 * dz) - (pC[b_0] * pC[b_xyz]) / (4 * dz) - (pC[a_z] * pC[a_zz]) / (6 * dz) -
777 (pC[a_xz] * pC[a_zz]) / (12 * dz) - (pC[a_y] * pC[a_z]) / (2 * dz) - (pC[a_xzz] * pC[a_z]) / (12 * 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) -
779 (pC[a_0] * pC[a_z]) / dz - (pC[a_y] * pC[a_yz]) / (4 * dz) - (pC[a_xy] * pC[a_yz]) / (8 * dz) -
780 (pC[a_xx] * pC[a_yz]) / (12 * dz) - (pC[a_x] * pC[a_yz]) / (4 * dz) - (pC[a_0] * pC[a_yz]) / (2 * dz) -
781 (pC[a_xz] * pC[a_y]) / (4 * dz) - (pC[a_xyz] * pC[a_y]) / (8 * dz) - (pC[a_xxz] * pC[a_y]) / (12 * dz) -
782 (pC[a_xz] * pC[a_xzz]) / (24 * dz) - (pC[a_xy] * pC[a_xz]) / (8 * dz) - (pC[a_xx] * pC[a_xz]) / (12 * dz) -
783 (pC[a_x] * pC[a_xz]) / (4 * dz) - (pC[a_0] * pC[a_xz]) / (2 * dz) - (pC[a_xy] * pC[a_xyz]) / (16 * dz) -
784 (pC[a_xx] * pC[a_xyz]) / (24 * dz) - (pC[a_x] * pC[a_xyz]) / (8 * dz) - (pC[a_0] * pC[a_xyz]) / (4 * dz) -
785 (pC[a_xxz] * pC[a_xy]) / (24 * dz) - (pC[a_xx] * pC[a_xxz]) / (36 * dz) - (pC[a_x] * pC[a_xxz]) / (12 * dz) -
786 (pC[a_0] * pC[a_xxz]) / (6 * dz) + (pC[b_z] * pC[c_yz]) / (12 * dy) + (pC[b_yz] * pC[c_yz]) / (24 * dy) +
787 (pC[b_z] * pC[c_yyz]) / (12 * dy) + (pC[b_yz] * pC[c_yyz]) / (24 * dy) + (pC[b_yy] * pC[c_yy]) / (6 * dy) +
788 (pC[b_y] * pC[c_yy]) / (2 * dy) + (pC[b_xy] * pC[c_yy]) / (4 * dy) + (pC[b_x] * pC[c_yy]) / (2 * dy) +
789 (pC[b_0] * pC[c_yy]) / dy + (pC[b_yy] * pC[c_y]) / (6 * dy) + (pC[b_y] * pC[c_y]) / (2 * dy) +
790 (pC[b_xy] * pC[c_y]) / (4 * dy) + (pC[b_x] * pC[c_y]) / (2 * dy) + (pC[b_0] * pC[c_y]) / dy +
791 (pC[b_z] * pC[c_xyz]) / (24 * dy) + (pC[b_yz] * pC[c_xyz]) / (48 * dy) + (pC[b_yy] * pC[c_xy]) / (12 * dy) +
792 (pC[b_y] * pC[c_xy]) / (4 * dy) + (pC[b_xy] * pC[c_xy]) / (8 * dy) + (pC[b_x] * pC[c_xy]) / (4 * dy) +
793 (pC[b_0] * pC[c_xy]) / (2 * dy) + (pC[a_z] * pC[c_xz]) / (12 * dx) + (pC[a_xz] * pC[c_xz]) / (24 * dx) +
794 (pC[a_z] * pC[c_xyz]) / (24 * dx) + (pC[a_xz] * pC[c_xyz]) / (48 * dx) + (pC[a_y] * pC[c_xy]) / (4 * dx) +
795 (pC[a_xy] * pC[c_xy]) / (8 * dx) + (pC[a_xx] * pC[c_xy]) / (12 * dx) + (pC[a_x] * pC[c_xy]) / (4 * dx) +
796 (pC[a_0] * pC[c_xy]) / (2 * dx) + (pC[a_z] * pC[c_xxz]) / (12 * dx) + (pC[a_xz] * pC[c_xxz]) / (24 * dx) +
797 (pC[a_y] * pC[c_xx]) / (2 * dx) + (pC[a_xy] * pC[c_xx]) / (4 * dx) + (pC[a_xx] * pC[c_xx]) / (6 * dx) +
798 (pC[a_x] * pC[c_xx]) / (2 * dx) + (pC[a_0] * pC[c_xx]) / dx + (pC[a_y] * pC[c_x]) / (2 * dx) +
799 (pC[a_xy] * pC[c_x]) / (4 * dx) + (pC[a_xx] * pC[c_x]) / (6 * dx) + (pC[a_x] * pC[c_x]) / (2 * dx) +
800 (pC[a_0] * pC[c_x]) / dx;
801}
802
803template <typename REAL>
804inline REAL JXB(fsgrids::ehall term, const std::array<REAL, Rec::N_REC_COEFFICIENTS>& pC, Real BGBX, Real BGBY,
805 Real BGBZ, const std::array<Real, 3>& gridSpacing) {
806 switch (term) {
808 return JXBX_000_100(pC, BGBY, BGBZ, gridSpacing);
809 }
811 return JXBX_010_110(pC, BGBY, BGBZ, gridSpacing);
812 }
814 return JXBX_001_101(pC, BGBY, BGBZ, gridSpacing);
815 }
817 return JXBX_011_111(pC, BGBY, BGBZ, gridSpacing);
818 }
820 return JXBY_000_010(pC, BGBX, BGBZ, gridSpacing);
821 }
823 return JXBY_100_110(pC, BGBX, BGBZ, gridSpacing);
824 }
826 return JXBY_001_011(pC, BGBX, BGBZ, gridSpacing);
827 }
829 return JXBY_101_111(pC, BGBX, BGBZ, gridSpacing);
830 }
832 return JXBZ_000_001(pC, BGBX, BGBY, gridSpacing);
833 }
835 return JXBZ_100_101(pC, BGBX, BGBY, gridSpacing);
836 }
838 return JXBZ_010_011(pC, BGBX, BGBY, gridSpacing);
839 }
841 return JXBZ_110_111(pC, BGBX, BGBY, gridSpacing);
842 }
843 case fsgrids::N_EHALL: {
844 break;
845 }
846 }
847 return 0.0;
848}
849
867 fsgrids::ehallspan ehalls,
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];
880
881 const Real bgbx = bgb[fsgrids::bgbfield::BGBX];
882 const Real bgby = bgb[fsgrids::bgbfield::BGBY];
883 const Real bgbz = bgb[fsgrids::bgbfield::BGBZ];
884
885 auto computeHallRhoq = [&moments, &moment](const std::array<size_t, 4>& indices) {
886 const auto min = Parameters::hallMinimumRhoq;
887 const auto max = std::numeric_limits<Real>::max();
888
889 return std::clamp(
891 ? moment[fsgrids::moments::RHOQ]
892 : FOURTH * (moments[indices[0]][fsgrids::moments::RHOQ] + moments[indices[1]][fsgrids::moments::RHOQ] +
893 moments[indices[2]][fsgrids::moments::RHOQ] + moments[indices[3]][fsgrids::moments::RHOQ]),
894 min, max);
895 };
896
897 switch (Parameters::ohmHallTerm) {
898 case 0:
899 cerr << __FILE__ << __LINE__ << "You shouldn't be in a Hall term function if Parameters::ohmHallTerm == 0."
900 << endl;
901 break;
902
903 case 1: {
904 const Real Bx = perb[fsgrids::bfield::PERBX] + bgbx;
905 const Real By = perb[fsgrids::bfield::PERBY] + bgby;
906 const Real Bz = perb[fsgrids::bfield::PERBZ] + bgbz;
907
908 const Real invHallRhoqMU0 = 1.0 / (physicalconstants::MU_0 * computeHallRhoq({}));
909
910 const Real ydx = (bgb[fsgrids::bgbfield::dBGBydx] + dperb[fsgrids::dperb::dPERBydx]) / gridSpacing[0];
911 const Real zdx = (bgb[fsgrids::bgbfield::dBGBzdx] + dperb[fsgrids::dperb::dPERBzdx]) / gridSpacing[0];
912 const Real xdy = (bgb[fsgrids::bgbfield::dBGBxdy] + dperb[fsgrids::dperb::dPERBxdy]) / gridSpacing[1];
913 const Real zdy = (bgb[fsgrids::bgbfield::dBGBzdy] + dperb[fsgrids::dperb::dPERBzdy]) / gridSpacing[1];
914 const Real xdz = (bgb[fsgrids::bgbfield::dBGBxdz] + dperb[fsgrids::dperb::dPERBxdz]) / gridSpacing[2];
915 const Real ydz = (bgb[fsgrids::bgbfield::dBGBydz] + dperb[fsgrids::dperb::dPERBydz]) / gridSpacing[2];
916
917 const Real EXHall = (Bz * (xdz - zdx) - By * (ydx - xdy)) * invHallRhoqMU0;
918 ehall[fsgrids::ehall::EXHALL_000_100] = EXHall;
919 ehall[fsgrids::ehall::EXHALL_010_110] = EXHall;
920 ehall[fsgrids::ehall::EXHALL_001_101] = EXHall;
921 ehall[fsgrids::ehall::EXHALL_011_111] = EXHall;
922
923 const Real EYHall = (Bx * (ydx - xdy) - Bz * (zdy - ydz)) * invHallRhoqMU0;
924 ehall[fsgrids::ehall::EYHALL_000_010] = EYHall;
925 ehall[fsgrids::ehall::EYHALL_100_110] = EYHall;
926 ehall[fsgrids::ehall::EYHALL_101_111] = EYHall;
927 ehall[fsgrids::ehall::EYHALL_001_011] = EYHall;
928
929 const Real EZHall = (By * (zdy - ydz) - Bx * (xdz - zdx)) * invHallRhoqMU0;
930 ehall[fsgrids::ehall::EZHALL_000_001] = EZHall;
931 ehall[fsgrids::ehall::EZHALL_100_101] = EZHall;
932 ehall[fsgrids::ehall::EZHALL_110_111] = EZHall;
933 ehall[fsgrids::ehall::EZHALL_010_011] = EZHall;
934
935 break;
936 }
937 case 2: {
938 auto computeEHall = [&perturbedCoefficients, &bgbx, &bgby, &bgbz, &gridSpacing](fsgrids::ehall term,
939 Real hallRhoq) {
940 return JXB(term, perturbedCoefficients, bgbx, bgby, bgbz, gridSpacing) / (physicalconstants::MU_0 * hallRhoq);
941 };
942
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();
949 // clang-format off
950 const std::array<std::array<size_t, 4>, 12> indices = {
951 std::array{
952 ooo,
953 omo,
954 oom,
955 stencil.omm(),
956 },
957 std::array{
958 ooo,
959 moo,
960 oom,
961 stencil.mom(),
962 },
963 std::array{
964 ooo,
965 moo,
966 omo,
967 stencil.mmo(),
968 },
969 std::array{
970 ooo,
971 poo,
972 oom,
973 stencil.pom(),
974 },
975 std::array{
976 ooo,
977 poo,
978 omo,
979 stencil.pmo(),
980 },
981 std::array{
982 ooo,
983 opo,
984 oom,
985 stencil.opm(),
986 },
987 std::array{
988 ooo,
989 moo,
990 opo,
991 stencil.mpo(),
992 },
993 std::array{
994 ooo,
995 poo,
996 opo,
997 stencil.ppo(),
998 },
999 std::array{
1000 ooo,
1001 omo,
1002 oop,
1003 stencil.omp(),
1004 },
1005 std::array{
1006 ooo,
1007 moo,
1008 oop,
1009 stencil.mop(),
1010 },
1011 std::array{
1012 ooo,
1013 poo,
1014 oop,
1015 stencil.pop(),
1016 },
1017 std::array{
1018 ooo,
1019 opo,
1020 oop,
1021 stencil.opp(),
1022 },
1023 };
1024 // clang-format on
1025
1026 const std::array<fsgrids::ehall, 12> terms = {
1031 };
1032
1033 for (size_t index = 0; index < terms.size(); index++) {
1034 ehall[terms[index]] = computeEHall(terms[index], computeHallRhoq(indices[index]));
1035 }
1036
1037 break;
1038 }
1039
1040 default:
1041 cerr << __FILE__ << ":" << __LINE__ << "You are welcome to code higher-order Hall term correction terms." << endl;
1042 break;
1043 }
1044}
1045
1061 fsgrids::ehallspan ehall,
1065 fsgrids::consttechnicalspan technical, const fsgrid::FsStencil& stencil,
1066 SysBoundary& sysBoundaries, const std::array<Real, 3>& gridSpacing) {
1067#ifdef DEBUG_FSOLVER
1068 if (!stencil.cellExists(0, 0, 0)) {
1069 cerr << "Out-of-bounds access in " << __FILE__ << ":" << __LINE__ << endl;
1070 exit(1);
1071 }
1072#endif
1073
1074 const auto& tech = technical[stencil.ooo()];
1075 cuint cellSysBoundaryFlag = tech.sysBoundaryFlag;
1076 cuint cellSysBoundaryLayer = tech.sysBoundaryLayer;
1077
1078 if (cellSysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE ||
1079 cellSysBoundaryFlag == sysboundarytype::OUTER_BOUNDARY_PADDING) {
1080 return;
1081 }
1082
1083 // Reconstruction order of the fields after Balsara 2009, 2 used for general B, 3 used
1084 // here for 2nd-order Hall term
1085 const std::array<Real, Rec::N_REC_COEFFICIENTS> perturbedCoefficients =
1086 reconstructionCoefficients(perb, dperb, stencil, 3);
1087
1088 if ((cellSysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY) && (cellSysBoundaryLayer != 1)) {
1089 auto* const sb = sysBoundaries.getSysBoundary(cellSysBoundaryFlag);
1090 sb->fieldSolverBoundaryCondHallElectricField(ehall, stencil, 0);
1091 sb->fieldSolverBoundaryCondHallElectricField(ehall, stencil, 1);
1092 sb->fieldSolverBoundaryCondHallElectricField(ehall, stencil, 2);
1093 } else {
1094 calculateEdgeHallTermComponents(perb, ehall, moments, dperb, bgb, gridSpacing, perturbedCoefficients, stencil);
1095 }
1096}
1097
1121 fsgrids::perbspan perbdt2,
1122 fsgrids::ehallspan ehall,
1123 fsgrids::momentsspan moments,
1124 fsgrids::momentsspan momentsdt2,
1125 fsgrids::dperbspan dperb,
1126 fsgrids::dmomentsspan dmoments,
1127 fsgrids::dmomentsspan dmomentsdt2,
1128 fsgrids::bgbspan bgb,
1130 SysBoundary& sysBoundaries, int32_t RKCase, const bool communicateMomentsDerivatives) {
1131
1132
1133 if (not(RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2)) {
1134 perb = perbdt2;
1135 moments = momentsdt2;
1136 dmoments = dmomentsdt2;
1137 }
1138 phiprof::Timer hallTimer{"Calculate Hall term"};
1139
1140 const size_t numCells = fsgrid.getNumCells();
1141
1142 phiprof::Timer mpiTimer{"EHall ghost updates MPI", {"MPI"}};
1143 fsgrid.updateGhostCells(dperb);
1144 if (P::ohmGradPeTerm == 0 && communicateMomentsDerivatives) {
1145 fsgrid.updateGhostCells(dmoments);
1146 }
1147 mpiTimer.stop();
1148
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);
1153 });
1154
1155 hallTimer.stop(numCells, "Spatial Cells");
1156}
dx
Definition Dispersion.m:38
virtual void fieldSolverBoundaryCondHallElectricField(fsgrids::ehallspan ehall, const fsgrid::FsStencil &stencil, cuint component)=0
SysBoundary contains the SysBoundaryConditions used in the simulation.
Definition sysboundary.h:54
SBC::SysBoundaryCondition * getSysBoundary(cuint sysBoundaryType) const
@ RK_ORDER1
Definition common.h:508
@ RK_ORDER2_STEP2
Definition common.h:510
const uint32_t cuint
Definition definitions.h:50
float Real
Definition definitions.h:41
fsgrid::FsGrid< FS_STENCIL_WIDTH > FieldSolverGrid
Definition definitions.h:78
const float creal
Definition definitions.h:42
std::array< Real, Rec::N_REC_COEFFICIENTS > reconstructionCoefficients(fsgrids::perbspan perb, fsgrids::constdperbspan dperb, const fsgrid::FsStencil &stencil, Real reconstructionOrder)
Low-level helper function.
Definition fs_common.cpp:53
const Real FOURTH
Definition fs_common.h:53
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)
Definition ldz_hall.cpp:804
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.
Definition ldz_hall.cpp:367
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.
Definition ldz_hall.cpp:238
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.
Definition ldz_hall.cpp:46
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.
Definition ldz_hall.cpp:560
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.
Definition ldz_hall.cpp:431
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.
Definition ldz_hall.cpp:303
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.
Definition ldz_hall.cpp:688
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.
Definition ldz_hall.cpp:624
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.
Definition ldz_hall.cpp:174
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.
Definition ldz_hall.cpp:110
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.
Definition ldz_hall.cpp:752
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.
Definition ldz_hall.cpp:866
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.
Definition ldz_hall.cpp:495
#define index(i, j, k)
@ b_xx
Definition fs_common.h:89
@ c_xxz
Definition fs_common.h:90
@ b_xxy
Definition fs_common.h:89
@ b_xyy
Definition fs_common.h:89
@ b_yz
Definition fs_common.h:89
@ a_x
Definition fs_common.h:88
@ a_z
Definition fs_common.h:88
@ a_xx
Definition fs_common.h:88
@ a_xy
Definition fs_common.h:88
@ a_y
Definition fs_common.h:88
@ c_xyz
Definition fs_common.h:90
@ b_yzz
Definition fs_common.h:89
@ c_yz
Definition fs_common.h:90
@ c_xzz
Definition fs_common.h:90
@ a_xxy
Definition fs_common.h:88
@ b_xyz
Definition fs_common.h:89
@ a_0
Definition fs_common.h:88
@ a_xxz
Definition fs_common.h:88
@ b_x
Definition fs_common.h:89
@ c_yy
Definition fs_common.h:90
@ b_yyz
Definition fs_common.h:89
@ a_yz
Definition fs_common.h:88
@ b_yy
Definition fs_common.h:89
@ b_y
Definition fs_common.h:89
@ c_xy
Definition fs_common.h:90
@ c_zz
Definition fs_common.h:90
@ c_y
Definition fs_common.h:90
@ b_0
Definition fs_common.h:89
@ c_yzz
Definition fs_common.h:90
@ c_xz
Definition fs_common.h:90
@ c_yyz
Definition fs_common.h:90
@ a_xyz
Definition fs_common.h:88
@ a_yy
Definition fs_common.h:88
@ a_xz
Definition fs_common.h:88
@ c_x
Definition fs_common.h:90
@ a_xyy
Definition fs_common.h:88
@ c_z
Definition fs_common.h:90
@ b_xy
Definition fs_common.h:89
@ b_z
Definition fs_common.h:89
@ a_xzz
Definition fs_common.h:88
@ b_xz
Definition fs_common.h:89
@ a_zz
Definition fs_common.h:88
@ c_xx
Definition fs_common.h:90
@ b_zz
Definition fs_common.h:89
@ c_0
Definition fs_common.h:90
@ BGBY
Definition common.h:376
@ dBGBxdz
Definition common.h:385
@ dBGBydx
Definition common.h:386
@ BGBZ
Definition common.h:377
@ dBGBydz
Definition common.h:387
@ BGBX
Definition common.h:375
@ dBGBzdx
Definition common.h:388
@ dBGBxdy
Definition common.h:384
@ dBGBzdy
Definition common.h:389
std::span< std::array< Real, fsgrids::bfield::N_BFIELD > > perbspan
Definition common.h:434
std::span< const std::array< Real, bgbfield::N_BGB > > constbgbspan
Definition common.h:445
std::span< const std::array< Real, fsgrids::moments::N_MOMENTS > > constmomentsspan
Definition common.h:447
std::span< std::array< Real, fsgrids::moments::N_MOMENTS > > momentsspan
Definition common.h:446
std::span< std::array< Real, fsgrids::dmoments::N_DMOMENTS > > dmomentsspan
Definition common.h:448
@ RHOQ
Definition common.h:313
std::span< technical > technicalspan
Definition common.h:452
@ EZHALL_010_011
Definition common.h:295
@ N_EHALL
Definition common.h:301
@ EYHALL_101_111
Definition common.h:299
@ EYHALL_100_110
Definition common.h:292
@ EXHALL_010_110
Definition common.h:294
@ EZHALL_110_111
Definition common.h:296
@ EZHALL_000_001
Definition common.h:291
@ EYHALL_001_011
Definition common.h:298
@ EXHALL_001_101
Definition common.h:297
@ EYHALL_000_010
Definition common.h:290
@ EXHALL_000_100
Definition common.h:289
@ EXHALL_011_111
Definition common.h:300
@ EZHALL_100_101
Definition common.h:293
std::span< std::array< Real, bgbfield::N_BGB > > bgbspan
Definition common.h:444
std::span< std::array< Real, fsgrids::dperb::N_DPERB > > dperbspan
Definition common.h:442
std::span< const technical > consttechnicalspan
Definition common.h:453
@ PERBY
Definition common.h:276
@ PERBZ
Definition common.h:277
@ PERBX
Definition common.h:275
std::span< std::array< Real, fsgrids::ehall::N_EHALL > > ehallspan
Definition common.h:438
@ dPERBzdy
Definition common.h:329
@ dPERBydx
Definition common.h:326
@ dPERBzdx
Definition common.h:328
@ dPERBxdy
Definition common.h:324
@ dPERBxdz
Definition common.h:325
@ dPERBydz
Definition common.h:327
std::span< const std::array< Real, fsgrids::dperb::N_DPERB > > constdperbspan
Definition common.h:443
const Real MU_0
Definition common.h:570
@ OUTER_BOUNDARY_PADDING
Definition common.h:493
static uint ohmGradPeTerm
Definition parameters.h:144
static uint ohmHallTerm
Definition parameters.h:142
static Real hallMinimumRhoq
Definition parameters.h:161
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)