Line data Source code
1 : #ifndef NUMERIC_STENCIL_H_
2 : #define NUMERIC_STENCIL_H_
3 :
4 : #include <AMReX.H>
5 : #include <AMReX_MultiFab.H>
6 : #include "Set/Set.H"
7 :
8 : /// \brief This namespace contains some numerical tools
9 : namespace Numeric
10 : {
11 :
12 : enum StencilType { Lo, Hi, Central };
13 : static std::array<StencilType, AMREX_SPACEDIM>
14 : DefaultType = { AMREX_D_DECL(StencilType::Central, StencilType::Central, StencilType::Central) };
15 :
16 : [[nodiscard]]
17 : static
18 : AMREX_FORCE_INLINE
19 : std::array<StencilType, AMREX_SPACEDIM>
20 : GetStencil(const int i, const int j, const int k, const amrex::Box domain)
21 : {
22 : (void)i; (void)j; (void)k; // Suppress "unused variable" warnings.
23 : std::array<StencilType, AMREX_SPACEDIM> sten;
24 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
25 28625233 : AMREX_D_TERM(sten[0] = (i == lo.x ? Numeric::StencilType::Hi :
26 : i == hi.x ? Numeric::StencilType::Lo :
27 : Numeric::StencilType::Central);,
28 : sten[1] = (j == lo.y ? Numeric::StencilType::Hi :
29 : j == hi.y ? Numeric::StencilType::Lo :
30 : Numeric::StencilType::Central);,
31 : sten[2] = (k == lo.z ? Numeric::StencilType::Hi :
32 : k == hi.z ? Numeric::StencilType::Lo :
33 : Numeric::StencilType::Central););
34 28625233 : return sten;
35 : }
36 :
37 : template<class T, int x, int y, int z>
38 : struct Stencil
39 : {};
40 :
41 : //
42 : // FIRST order derivatives
43 : //
44 :
45 : template<class T>
46 : struct Stencil<T, 1, 0, 0>
47 : {
48 : [[nodiscard]] AMREX_FORCE_INLINE
49 : static T D(const amrex::Array4<const T>& f,
50 : const int& i, const int& j, const int& k, const int& m,
51 : const Set::Scalar dx[AMREX_SPACEDIM],
52 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
53 : {
54 113833255 : if (stencil[0] == StencilType::Lo)
55 30667257 : return (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0]; // 1st order stencil
56 103610854 : else if (stencil[0] == StencilType::Hi)
57 30667167 : return (f(i + 1, j, k, m) - f(i, j, k, m)) / dx[0]; // 1st order stencil
58 : else
59 317268429 : return (f(i + 1, j, k, m) - f(i - 1, j, k, m)) * 0.5 / dx[0];
60 : };
61 : [[nodiscard]] AMREX_FORCE_INLINE
62 : static std::pair<Set::Scalar,T>
63 : Dsplit( const amrex::Array4<const T>& f,
64 : const int& i, const int& j, const int& k, const int& m,
65 : const Set::Scalar dx[AMREX_SPACEDIM],
66 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
67 : {
68 0 : if (stencil[0] == StencilType::Lo)
69 : return std::make_pair(
70 0 : 1.0 / dx[0],
71 0 : -f(i - 1, j, k, m) / dx[0]
72 0 : ); // 1st order
73 0 : else if (stencil[0] == StencilType::Hi)
74 : return std::make_pair(
75 0 : -1.0 / dx[0],
76 0 : f(i + 1, j, k, m) / dx[0]
77 0 : ); // 1st order
78 : else
79 : return std::make_pair(
80 0 : 0.0,
81 0 : (f(i + 1, j, k, m) - f(i - 1, j, k, m)) * 0.5 / dx[0]
82 0 : );
83 : };
84 : };
85 :
86 : template<class T>
87 : struct Stencil<T, 0, 1, 0>
88 : {
89 : [[nodiscard]] AMREX_FORCE_INLINE
90 : static T D(const amrex::Array4<const T>& f,
91 : const int& i, const int& j, const int& k, const int& m,
92 : const Set::Scalar dx[AMREX_SPACEDIM],
93 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
94 : {
95 113833255 : if (stencil[1] == StencilType::Lo)
96 29696267 : return (f(i, j, k, m) - f(i, j - 1, k, m)) / dx[1];
97 103934552 : else if (stencil[1] == StencilType::Hi)
98 29696177 : return (f(i, j + 1, k, m) - f(i, j, k, m)) / dx[1];
99 : else
100 319210047 : return (f(i, j + 1, k, m) - f(i, j - 1, k, m)) * 0.5 / dx[1];
101 : };
102 : [[nodiscard]] AMREX_FORCE_INLINE
103 : static std::pair<Set::Scalar,T>
104 : Dsplit(const amrex::Array4<const T>& f,
105 : const int& i, const int& j, const int& k, const int& m,
106 : const Set::Scalar dx[AMREX_SPACEDIM],
107 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
108 : {
109 0 : if (stencil[1] == StencilType::Lo)
110 : return std::make_pair(
111 0 : 1.0 / dx[1],
112 0 : - f(i, j - 1, k, m) / dx[1]
113 0 : );
114 0 : else if (stencil[1] == StencilType::Hi)
115 : return std::make_pair(
116 0 : - 1.0 / dx[1],
117 0 : f(i, j + 1, k, m) / dx[1]
118 0 : );
119 : else
120 : return std::make_pair(
121 0 : 0.0,
122 0 : (f(i, j + 1, k, m) - f(i, j - 1, k, m)) * 0.5 / dx[1]
123 0 : );
124 : };
125 : };
126 :
127 : template<class T>
128 : struct Stencil<T, 0, 0, 1>
129 : {
130 : [[nodiscard]] AMREX_FORCE_INLINE
131 : static T D(const amrex::Array4<const T>& f,
132 : const int& i, const int& j, const int& k, const int& m,
133 : const Set::Scalar dx[AMREX_SPACEDIM],
134 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
135 : {
136 35200721 : if (stencil[2] == StencilType::Lo)
137 24720000 : return (f(i, j, k, m) - f(i, j, k - 1, m)) / dx[2];
138 26960721 : else if (stencil[2] == StencilType::Hi)
139 24720000 : return (f(i, j, k + 1, m) - f(i, j, k, m)) / dx[2];
140 : else
141 61256463 : return (f(i, j, k + 1, m) - f(i, j, k - 1, m)) * 0.5 / dx[2];
142 : };
143 : [[nodiscard]] AMREX_FORCE_INLINE
144 : static std::pair<Set::Scalar, T>
145 : Dsplit(const amrex::Array4<const T>& f,
146 : const int& i, const int& j, const int& k, const int& m,
147 : const Set::Scalar dx[AMREX_SPACEDIM],
148 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
149 : {
150 0 : if (stencil[2] == StencilType::Lo)
151 : return std::make_pair(
152 0 : 1.0 / dx[2],
153 0 : -f(i, j, k - 1, m) / dx[2]
154 0 : );
155 0 : else if (stencil[2] == StencilType::Hi)
156 : return std::make_pair(
157 0 : -1.0 / dx[2],
158 0 : f(i, j, k + 1, m) / dx[2]
159 0 : );
160 : else
161 : return std::make_pair(
162 0 : 0.0,
163 0 : (f(i, j, k + 1, m) - f(i, j, k - 1, m)) * 0.5 / dx[2]
164 0 : );
165 : };
166 : };
167 :
168 : //
169 : // SECOND order derivatives
170 : //
171 :
172 : template<class T>
173 : struct Stencil<T, 2, 0, 0>
174 : {
175 : [[nodiscard]] AMREX_FORCE_INLINE
176 : static T D(const amrex::Array4<const T>& f,
177 : const int& i, const int& j, const int& k, const int& m,
178 : const Set::Scalar dx[AMREX_SPACEDIM],
179 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
180 : {
181 363410026 : if (stencil[0] == StencilType::Central)
182 1453473332 : return (f(i + 1, j, k, m) - 2.0 * f(i, j, k, m) + f(i - 1, j, k, m)) / dx[0] / dx[0];
183 : else
184 83728 : return 0.0 * f(i, j, k, m); // TODO this is not a great way to do this
185 : };
186 : };
187 :
188 : template<class T>
189 : struct Stencil<T, 0, 2, 0>
190 : {
191 : [[nodiscard]] AMREX_FORCE_INLINE
192 : static T D(const amrex::Array4<const T>& f,
193 : const int& i, const int& j, const int& k, const int& m,
194 : const Set::Scalar dx[AMREX_SPACEDIM],
195 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
196 : {
197 363447052 : if (stencil[1] == StencilType::Central)
198 1453619036 : return (f(i, j + 1, k, m) - 2.0 * f(i, j, k, m) + f(i, j - 1, k, m)) / dx[1] / dx[1];
199 : else
200 84928 : return 0.0 * f(i, j, k, m); // TODO this is not a great way to do this
201 : };
202 : };
203 :
204 : template<class T>
205 : struct Stencil<T, 0, 0, 2>
206 : {
207 : [[nodiscard]] AMREX_FORCE_INLINE
208 : static T D(const amrex::Array4<const T>& f,
209 : const int& i, const int& j, const int& k, const int& m,
210 : const Set::Scalar dx[AMREX_SPACEDIM],
211 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
212 : {
213 5965565 : if (stencil[2] == StencilType::Central)
214 23862258 : return (f(i, j, k + 1, m) - 2.0 * f(i, j, k, m) + f(i, j, k - 1, m)) / dx[2] / dx[2];
215 : else
216 0 : return 0.0 * f(i, j, k, m); // TODO this is not a great way to do this
217 : };
218 : };
219 :
220 : template<class T>
221 : struct Stencil<T, 1, 1, 0>
222 : {
223 : [[nodiscard]] AMREX_FORCE_INLINE
224 : static T D(const amrex::Array4<const T>& f,
225 : const int& i, const int& j, const int& k, const int& m,
226 : const Set::Scalar dx[AMREX_SPACEDIM],
227 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
228 : {
229 31077186 : int ihi = 1, ilo = 1, jhi = 1, jlo = 1;
230 31077186 : Set::Scalar ifac = 0.5, jfac = 0.5;
231 31077186 : if (stencil[0] == StencilType::Hi) { ilo = 0; ifac = 1.0; }
232 31077186 : if (stencil[0] == StencilType::Lo) { ihi = 0; ifac = 1.0; }
233 31077186 : if (stencil[1] == StencilType::Hi) { jlo = 0; jfac = 1.0; }
234 31077186 : if (stencil[1] == StencilType::Lo) { jhi = 0; jfac = 1.0; }
235 :
236 155386430 : return ifac * jfac * (f(i + ihi, j + jhi, k, m) + f(i - ilo, j - jlo, k, m) - f(i + ihi, j - jlo, k, m) - f(i - ilo, j + jhi, k, m)) / (dx[0] * dx[1]);
237 : };
238 : };
239 : template<class T>
240 : struct Stencil<T, 1, 0, 1>
241 : {
242 : [[nodiscard]] AMREX_FORCE_INLINE
243 : static T D(const amrex::Array4<const T>& f,
244 : const int& i, const int& j, const int& k, const int& m,
245 : const Set::Scalar dx[AMREX_SPACEDIM],
246 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
247 : {
248 5130237 : int khi = 1, klo = 1, ihi = 1, ilo = 1;
249 5130237 : Set::Scalar kfac = 0.5, ifac = 0.5;
250 5130237 : if (stencil[0] == StencilType::Hi) { ilo = 0; ifac = 1.0; }
251 5130237 : if (stencil[0] == StencilType::Lo) { ihi = 0; ifac = 1.0; }
252 5130237 : if (stencil[2] == StencilType::Hi) { klo = 0; kfac = 1.0; }
253 5130237 : if (stencil[2] == StencilType::Lo) { khi = 0; kfac = 1.0; }
254 :
255 25651185 : return kfac * ifac * (f(i + ihi, j, k + khi, m) + f(i - ilo, j, k - klo, m) - f(i + ihi, j, k - klo, m) - f(i - ilo, j, k + khi, m)) / (dx[0] * dx[2]);
256 : };
257 : };
258 : template<class T>
259 : struct Stencil<T, 0, 1, 1>
260 : {
261 : [[nodiscard]] AMREX_FORCE_INLINE
262 : static T D(const amrex::Array4<const T>& f,
263 : const int& i, const int& j, const int& k, const int& m,
264 : const Set::Scalar dx[AMREX_SPACEDIM],
265 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
266 : {
267 5130237 : int jhi = 1, jlo = 1, khi = 1, klo = 1;
268 5130237 : Set::Scalar jfac = 0.5, kfac = 0.5;
269 5130237 : if (stencil[1] == StencilType::Hi) { jlo = 0; jfac = 1.0; }
270 5130237 : if (stencil[1] == StencilType::Lo) { jhi = 0; jfac = 1.0; }
271 5130237 : if (stencil[2] == StencilType::Hi) { klo = 0; kfac = 1.0; }
272 5130237 : if (stencil[2] == StencilType::Lo) { khi = 0; kfac = 1.0; }
273 :
274 25651185 : return jfac * kfac * (f(i, j + jhi, k + khi, m) + f(i, j - jlo, k - klo, m) - f(i, j + jhi, k - klo, m) - f(i, j - jlo, k + khi, m)) / (dx[1] * dx[2]);
275 : };
276 : };
277 :
278 : //
279 : // FOURTH order derivatives
280 : //
281 :
282 : template<class T>
283 : struct Stencil<T, 4, 0, 0>
284 : {
285 : [[nodiscard]] AMREX_FORCE_INLINE
286 : static T D(const amrex::Array4<const T>& f,
287 : const int& i, const int& j, const int& k, const int& m,
288 : const Set::Scalar dx[AMREX_SPACEDIM])
289 : {
290 6999430 : return ((f(i + 2, j, k, m)) - 4. * (f(i + 1, j, k, m)) + 6. * (f(i, j, k, m)) - 4. * (f(i - 1, j, k, m)) + (f(i - 2, j, k, m))) /
291 1399886 : (dx[0] * dx[0] * dx[0] * dx[0]);
292 : };
293 : };
294 : template<class T>
295 : struct Stencil<T, 0, 4, 0>
296 : {
297 : [[nodiscard]] AMREX_FORCE_INLINE
298 : static T D(const amrex::Array4<const T>& f,
299 : const int& i, const int& j, const int& k, const int& m,
300 : const Set::Scalar dx[AMREX_SPACEDIM])
301 : {
302 5636570 : return ((f(i, j + 2, k, m)) - 4. * (f(i, j + 1, k, m)) + 6. * (f(i, j, k, m)) - 4. * (f(i, j - 1, k, m)) + (f(i, j - 2, k, m))) /
303 1399886 : (dx[1] * dx[1] * dx[1] * dx[1]);
304 : };
305 : };
306 : template<class T>
307 : struct Stencil<T, 0, 0, 4>
308 : {
309 : [[nodiscard]] AMREX_FORCE_INLINE
310 : static T D(const amrex::Array4<const T>& f,
311 : const int& i, const int& j, const int& k, const int& m,
312 : const Set::Scalar dx[AMREX_SPACEDIM])
313 : {
314 179685 : return ((f(i, j, k + 2, m)) - 4. * (f(i, j, k + 1, m)) + 6. * (f(i, j, k, m)) - 4. * (f(i, j, k - 1, m)) + (f(i, j, k - 2, m))) /
315 35937 : (dx[2] * dx[2] * dx[2] * dx[2]);
316 : };
317 : };
318 :
319 : /// Compute
320 : /// \f[ \frac{\partial^4}{\partial^3 x_1\partial x_2\f]
321 : template<class T>
322 : struct Stencil<T, 3, 1, 0>
323 : {
324 : [[nodiscard]] AMREX_FORCE_INLINE
325 : static T D(const amrex::Array4<const T>& f,
326 : const int& i, const int& j, const int& k, const int& m,
327 : const Set::Scalar dx[AMREX_SPACEDIM])
328 : {
329 4236684 : return ((-f(i + 2, j + 2, k, m) + 8.0 * f(i + 2, j + 1, k, m) - 8.0 * f(i + 2, j - 1, k, m) + f(i + 2, j - 2, k, m))
330 5599544 : - 2 * (-f(i + 1, j + 2, k, m) + 8.0 * f(i + 1, j + 1, k, m) - 8.0 * f(i + 1, j - 1, k, m) + f(i + 1, j - 2, k, m))
331 5599544 : + 2 * (-f(i - 1, j + 2, k, m) + 8.0 * f(i - 1, j + 1, k, m) - 8.0 * f(i - 1, j - 1, k, m) + f(i - 1, j - 2, k, m))
332 5599544 : - (-f(i - 2, j + 2, k, m) + 8.0 * f(i - 2, j + 1, k, m) - 8.0 * f(i - 2, j - 1, k, m) + f(i - 2, j - 2, k, m))) /
333 1399886 : (24.0 * dx[0] * dx[0] * dx[0] * dx[1]);
334 : };
335 : };
336 :
337 : /// Compute
338 : /// \f[ \frac{\partial^4}{\partial^3 x_2\partial x_1\f]
339 : template<class T>
340 : struct Stencil<T, 1, 3, 0>
341 : {
342 : [[nodiscard]] AMREX_FORCE_INLINE
343 : static T D(const amrex::Array4<const T>& f,
344 : const int& i, const int& j, const int& k, const int& m,
345 : const Set::Scalar dx[AMREX_SPACEDIM])
346 : {
347 4236684 : return ((-f(i + 2, j + 2, k, m) + 8.0 * f(i + 1, j + 2, k, m) - 8.0 * f(i - 1, j + 2, k, m) + f(i - 2, j + 2, k, m))
348 5599544 : - 2 * (-f(i + 2, j + 1, k, m) + 8.0 * f(i + 1, j + 1, k, m) - 8.0 * f(i - 1, j + 1, k, m) + f(i - 2, j + 1, k, m))
349 5599544 : + 2 * (-f(i + 2, j - 1, k, m) + 8.0 * f(i + 1, j - 1, k, m) - 8.0 * f(i - 1, j - 1, k, m) + f(i - 2, j - 1, k, m))
350 5599544 : - (-f(i + 2, j - 2, k, m) + 8.0 * f(i + 1, j - 2, k, m) - 8.0 * f(i - 1, j - 2, k, m) + f(i - 2, j - 2, k, m))) /
351 1399886 : (24.0 * dx[0] * dx[1] * dx[1] * dx[1]);
352 : };
353 : };
354 :
355 : /// Compute
356 : /// \f[ \frac{\partial^4}{\partial^3 x_2\partial x_3\f]
357 : template<class T>
358 : struct Stencil<T, 0, 3, 1>
359 : {
360 : [[nodiscard]] AMREX_FORCE_INLINE
361 : static T D(const amrex::Array4<const T>& f,
362 : const int& i, const int& j, const int& k, const int& m,
363 : const Set::Scalar dx[AMREX_SPACEDIM])
364 : {
365 143748 : return ((-f(i, j + 2, k + 2, m) + 8.0 * f(i, j + 2, k + 1, m) - 8.0 * f(i, j + 2, k - 1, m) + f(i, j + 2, k - 2, m))
366 143748 : - 2 * (-f(i, j + 1, k + 2, m) + 8.0 * f(i, j + 1, k + 1, m) - 8.0 * f(i, j + 1, k - 1, m) + f(i, j + 1, k - 2, m))
367 143748 : + 2 * (-f(i, j - 1, k + 2, m) + 8.0 * f(i, j - 1, k + 1, m) - 8.0 * f(i, j - 1, k - 1, m) + f(i, j - 1, k - 2, m))
368 143748 : - (-f(i, j - 2, k + 2, m) + 8.0 * f(i, j - 2, k + 1, m) - 8.0 * f(i, j - 2, k - 1, m) + f(i, j - 2, k - 2, m))) /
369 35937 : (24.0 * dx[1] * dx[1] * dx[1] * dx[2]);
370 : };
371 : };
372 : /// \brief Compute \f[ \frac{\partial^4}{\partial^3 x_3\partial x_2}\f]
373 : template<class T>
374 : struct Stencil<T, 0, 1, 3>
375 : {
376 : [[nodiscard]] AMREX_FORCE_INLINE
377 : static T D(const amrex::Array4<const T>& f,
378 : const int& i, const int& j, const int& k, const int& m,
379 : const Set::Scalar dx[AMREX_SPACEDIM])
380 : {
381 143748 : return ((-f(i, j + 2, k + 2, m) + 8.0 * f(i, j + 1, k + 2, m) - 8.0 * f(i, j - 1, k + 2, m) + f(i, j - 2, k + 2, m))
382 143748 : - 2 * (-f(i, j + 2, k + 1, m) + 8.0 * f(i, j + 1, k + 1, m) - 8.0 * f(i, j - 1, k + 1, m) + f(i, j - 2, k + 1, m))
383 143748 : + 2 * (-f(i, j + 2, k - 1, m) + 8.0 * f(i, j + 1, k - 1, m) - 8.0 * f(i, j - 1, k - 1, m) + f(i, j - 2, k - 1, m))
384 143748 : - (-f(i, j + 2, k - 2, m) + 8.0 * f(i, j + 1, k - 2, m) - 8.0 * f(i, j - 1, k - 2, m) + f(i, j - 2, k - 2, m))) /
385 35937 : (24.0 * dx[1] * dx[2] * dx[2] * dx[2]);
386 : };
387 : };
388 : /// \brief Compute \f[ \frac{\partial^4}{\partial^3 x_3\partial x_1}\f]
389 : template<class T>
390 : struct Stencil<T, 1, 0, 3>
391 : {
392 : [[nodiscard]] AMREX_FORCE_INLINE
393 : static T D(const amrex::Array4<const T>& f,
394 : const int& i, const int& j, const int& k, const int& m,
395 : const Set::Scalar dx[AMREX_SPACEDIM])
396 : {
397 143748 : return ((-f(i + 2, j, k + 2, m) + 8.0 * f(i + 1, j, k + 2, m) - 8.0 * f(i - 1, j, k + 2, m) + f(i - 2, j, k + 2, m))
398 143748 : - 2 * (-f(i + 2, j, k + 1, m) + 8.0 * f(i + 1, j, k + 1, m) - 8.0 * f(i - 1, j, k + 1, m) + f(i - 2, j, k + 1, m))
399 143748 : + 2 * (-f(i + 2, j, k - 1, m) + 8.0 * f(i + 1, j, k - 1, m) - 8.0 * f(i - 1, j, k - 1, m) + f(i - 2, j, k - 1, m))
400 143748 : - (-f(i + 2, j, k - 2, m) + 8.0 * f(i + 1, j, k - 2, m) - 8.0 * f(i - 1, j, k - 2, m) + f(i - 2, j, k - 2, m))) /
401 35937 : (24.0 * dx[0] * dx[2] * dx[2] * dx[2]);
402 :
403 : };
404 : };
405 : /// \brief Compute \f[ \frac{\partial^4}{\partial^3 x_1\partial x_3}\f]
406 : template<class T>
407 : struct Stencil<T, 3, 0, 1>
408 : {
409 : [[nodiscard]] AMREX_FORCE_INLINE
410 : static T D(const amrex::Array4<const T>& f,
411 : const int& i, const int& j, const int& k, const int& m,
412 : const Set::Scalar dx[AMREX_SPACEDIM])
413 : {
414 143748 : return ((-f(i + 2, j, k + 2, m) + 8.0 * f(i + 2, j, k + 1, m) - 8.0 * f(i + 2, j, k - 1, m) + f(i + 2, j, k - 2, m))
415 143748 : - 2 * (-f(i + 1, j, k + 2, m) + 8.0 * f(i + 1, j, k + 1, m) - 8.0 * f(i + 1, j, k - 1, m) + f(i + 1, j, k - 2, m))
416 143748 : + 2 * (-f(i - 1, j, k + 2, m) + 8.0 * f(i - 1, j, k + 1, m) - 8.0 * f(i - 1, j, k - 1, m) + f(i - 1, j, k - 2, m))
417 143748 : - (-f(i - 2, j, k + 2, m) + 8.0 * f(i - 2, j, k + 1, m) - 8.0 * f(i - 2, j, k - 1, m) + f(i - 2, j, k - 2, m))) /
418 35937 : (24.0 * dx[0] * dx[0] * dx[0] * dx[2]);
419 :
420 : };
421 : };
422 : /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_1\partial^2 x_2}\f]
423 : template<class T>
424 : struct Stencil<T, 2, 2, 0>
425 : {
426 : [[nodiscard]] AMREX_FORCE_INLINE
427 : static T D(const amrex::Array4<const T>& f,
428 : const int& i, const int& j, const int& k, const int& m,
429 : const Set::Scalar dx[AMREX_SPACEDIM])
430 : {
431 5599544 : return (-(-f(i + 2, j + 2, k, m) + 16.0 * f(i + 1, j + 2, k, m) - 30.0 * f(i, j + 2, k, m) + 16.0 * f(i - 1, j + 2, k, m) - f(i - 2, j + 2, k, m))
432 7036456 : + 16 * (-f(i + 2, j + 1, k, m) + 16.0 * f(i + 1, j + 1, k, m) - 30.0 * f(i, j + 1, k, m) + 16.0 * f(i - 1, j + 1, k, m) - f(i - 2, j + 1, k, m))
433 6999430 : - 30 * (-f(i + 2, j, k, m) + 16.0 * f(i + 1, j, k, m) - 30.0 * f(i, j, k, m) + 16.0 * f(i - 1, j, k, m) - f(i - 2, j, k, m))
434 6999430 : + 16 * (-f(i + 2, j - 1, k, m) + 16.0 * f(i + 1, j - 1, k, m) - 30.0 * f(i, j - 1, k, m) + 16.0 * f(i - 1, j - 1, k, m) - f(i - 2, j - 1, k, m))
435 6999430 : - (-f(i + 2, j - 2, k, m) + 16.0 * f(i + 1, j - 2, k, m) - 30.0 * f(i, j - 2, k, m) + 16.0 * f(i - 1, j - 2, k, m) - f(i - 2, j - 2, k, m))) /
436 1399886 : (144.0 * dx[0] * dx[0] * dx[1] * dx[1]);
437 : };
438 : };
439 :
440 : /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_2\partial^2 x_3}\f]
441 : template<class T>
442 : struct Stencil<T, 0, 2, 2>
443 : {
444 : [[nodiscard]] AMREX_FORCE_INLINE
445 : static T D(const amrex::Array4<const T>& f,
446 : const int& i, const int& j, const int& k, const int& m,
447 : const Set::Scalar dx[AMREX_SPACEDIM])
448 : {
449 143748 : return (-(-f(i, j + 2, k + 2, m) + 16.0 * f(i, j + 2, k + 1, m) - 30.0 * f(i, j + 2, k, m) + 16.0 * f(i, j + 2, k - 1, m) - f(i, j + 2, k - 2, m))
450 215622 : + 16 * (-f(i, j + 1, k + 2, m) + 16.0 * f(i, j + 1, k + 1, m) - 30.0 * f(i, j + 1, k, m) + 16.0 * f(i, j + 1, k - 1, m) - f(i, j + 1, k - 2, m))
451 179685 : - 30 * (-f(i, j, k + 2, m) + 16.0 * f(i, j, k + 1, m) - 30.0 * f(i, j, k, m) + 16.0 * f(i, j, k - 1, m) - f(i, j, k - 2, m))
452 179685 : + 16 * (-f(i, j - 1, k + 2, m) + 16.0 * f(i, j - 1, k + 1, m) - 30.0 * f(i, j - 1, k, m) + 16.0 * f(i, j - 1, k - 1, m) - f(i, j - 1, k - 2, m))
453 179685 : - (-f(i, j - 2, k + 2, m) + 16.0 * f(i, j - 2, k + 1, m) - 30.0 * f(i, j - 2, k, m) + 16.0 * f(i, j - 2, k - 1, m) - f(i, j - 2, k - 2, m))) /
454 35937 : (144.0 * dx[0] * dx[0] * dx[2] * dx[2]);
455 : };
456 : };
457 : /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_1\partial^2 x_3}\f]
458 : template<class T>
459 : struct Stencil<T, 2, 0, 2>
460 : {
461 : [[nodiscard]] AMREX_FORCE_INLINE
462 : static T D(const amrex::Array4<const T>& f,
463 : const int& i, const int& j, const int& k, const int& m,
464 : const Set::Scalar dx[AMREX_SPACEDIM])
465 : {
466 143748 : return (-(-f(i + 2, j, k + 2, m) + 16.0 * f(i + 2, j, k + 1, m) - 30.0 * f(i + 2, j, k, m) + 16.0 * f(i + 2, j, k - 1, m) - f(i + 2, j, k - 2, m))
467 215622 : + 16 * (-f(i + 1, j, k + 2, m) + 16.0 * f(i + 1, j, k + 1, m) - 30.0 * f(i + 1, j, k, m) + 16.0 * f(i + 1, j, k - 1, m) - f(i + 1, j, k - 2, m))
468 179685 : - 30 * (-f(i, j, k + 2, m) + 16.0 * f(i, j, k + 1, m) - 30.0 * f(i, j, k, m) + 16.0 * f(i, j, k - 1, m) - f(i, j, k - 2, m))
469 179685 : + 16 * (-f(i - 1, j, k + 2, m) + 16.0 * f(i - 1, j, k + 1, m) - 30.0 * f(i - 1, j, k, m) + 16.0 * f(i - 1, j, k - 1, m) - f(i - 1, j, k - 2, m))
470 179685 : - (-f(i - 2, j, k + 2, m) + 16.0 * f(i - 2, j, k + 1, m) - 30.0 * f(i - 2, j, k, m) + 16.0 * f(i - 2, j, k - 1, m) - f(i - 2, j, k - 2, m))) /
471 35937 : (144.0 * dx[0] * dx[0] * dx[2] * dx[2]);
472 : };
473 : };
474 :
475 :
476 : /// \brief Compute \f[ \frac{\partial^4}{\partial^2 x_1\partial x_2\partial x_3}\f]
477 : template<class T>
478 : struct Stencil<T, 2, 1, 1>
479 : {
480 : [[nodiscard]] AMREX_FORCE_INLINE
481 : static T D(const amrex::Array4<const T>& f,
482 : const int& i, const int& j, const int& k, const int& m,
483 : const Set::Scalar dx[AMREX_SPACEDIM])
484 : {
485 107811 : return (+(f(i + 1, j + 1, k + 1, m) - 2.0 * f(i, j + 1, k + 1, m) + f(i - 1, j + 1, k + 1, m))
486 107811 : - (f(i + 1, j - 1, k + 1, m) - 2.0 * f(i, j - 1, k + 1, m) + f(i - 1, j - 1, k + 1, m))
487 107811 : - (f(i + 1, j + 1, k - 1, m) - 2.0 * f(i, j + 1, k - 1, m) + f(i - 1, j + 1, k - 1, m))
488 107811 : + (f(i + 1, j - 1, k - 1, m) - 2.0 * f(i, j - 1, k - 1, m) + f(i - 1, j - 1, k - 1, m)))
489 35937 : / (4.0 * dx[0] * dx[0] * dx[1] * dx[2]);
490 :
491 : };
492 : };
493 : /// \brief Compute \f[ \frac{\partial^4}{\partial x_1\partial^2 x_2\partial x_3}\f]
494 : template<class T>
495 : struct Stencil<T, 1, 2, 1>
496 : {
497 : [[nodiscard]] AMREX_FORCE_INLINE
498 : static T D(const amrex::Array4<const T>& f,
499 : const int& i, const int& j, const int& k, const int& m,
500 : const Set::Scalar dx[AMREX_SPACEDIM])
501 : {
502 107811 : return (+(f(i + 1, j + 1, k + 1, m) - 2.0 * f(i + 1, j, k + 1, m) + f(i + 1, j - 1, k + 1, m))
503 107811 : - (f(i - 1, j + 1, k + 1, m) - 2.0 * f(i - 1, j, k + 1, m) + f(i - 1, j - 1, k + 1, m))
504 107811 : - (f(i + 1, j + 1, k - 1, m) - 2.0 * f(i + 1, j, k - 1, m) + f(i + 1, j - 1, k - 1, m))
505 107811 : + (f(i - 1, j + 1, k - 1, m) - 2.0 * f(i - 1, j, k - 1, m) + f(i - 1, j - 1, k - 1, m)))
506 35937 : / (4.0 * dx[0] * dx[1] * dx[1] * dx[2]);
507 :
508 : };
509 : };
510 : /// \brief Compute \f[ \frac{\partial^4}{\partial x_1\partial x_2\partial^2 x_3}\f]
511 : template<class T>
512 : struct Stencil<T, 1, 1, 2>
513 : {
514 : [[nodiscard]] AMREX_FORCE_INLINE
515 : static T D(const amrex::Array4<const T>& f,
516 : const int& i, const int& j, const int& k, const int& m,
517 : const Set::Scalar dx[AMREX_SPACEDIM])
518 : {
519 107811 : return (+(f(i + 1, j + 1, k + 1, m) - 2.0 * f(i + 1, j + 1, k, m) + f(i + 1, j + 1, k - 1, m))
520 107811 : - (f(i - 1, j + 1, k + 1, m) - 2.0 * f(i - 1, j + 1, k, m) + f(i - 1, j + 1, k - 1, m))
521 107811 : - (f(i + 1, j - 1, k + 1, m) - 2.0 * f(i + 1, j - 1, k, m) + f(i + 1, j - 1, k - 1, m))
522 107811 : + (f(i - 1, j - 1, k + 1, m) - 2.0 * f(i - 1, j - 1, k, m) + f(i - 1, j - 1, k - 1, m)))
523 35937 : / (4.0 * dx[0] * dx[1] * dx[1] * dx[2]);
524 :
525 : };
526 : };
527 :
528 : [[nodiscard]] AMREX_FORCE_INLINE
529 : Set::Scalar
530 : Laplacian(const amrex::Array4<const Set::Scalar>& f,
531 : const int& i, const int& j, const int& k, const int& m,
532 : const Set::Scalar dx[AMREX_SPACEDIM])
533 : {
534 332332840 : Set::Scalar ret = 0.0;
535 332332840 : ret += (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx));
536 : #if AMREX_SPACEDIM > 1
537 332332840 : ret += (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, m, dx));
538 : #if AMREX_SPACEDIM > 2
539 835328 : ret += (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, m, dx));
540 : #endif
541 : #endif
542 332332840 : return ret;
543 : }
544 :
545 : [[nodiscard]] AMREX_FORCE_INLINE
546 : Set::Vector
547 : Laplacian(const amrex::Array4<const Set::Vector>& f,
548 : const int& i, const int& j, const int& k,
549 : const Set::Scalar dx[AMREX_SPACEDIM])
550 : {
551 0 : Set::Vector ret = Set::Vector::Zero();
552 0 : ret += (Numeric::Stencil<Set::Vector, 2, 0, 0>::D(f, i, j, k, 0, dx));
553 : #if AMREX_SPACEDIM > 1
554 0 : ret += (Numeric::Stencil<Set::Vector, 0, 2, 0>::D(f, i, j, k, 0, dx));
555 : #if AMREX_SPACEDIM > 2
556 0 : ret += (Numeric::Stencil<Set::Vector, 0, 0, 2>::D(f, i, j, k, 0, dx));
557 : #endif
558 : #endif
559 0 : return ret;
560 : }
561 :
562 : [[nodiscard]] AMREX_FORCE_INLINE
563 : Set::Vector
564 : Divergence(const amrex::Array4<const Set::Matrix>& dw,
565 : const int& i, const int& j, const int& k,
566 : const Set::Scalar DX[AMREX_SPACEDIM],
567 : std::array<StencilType, AMREX_SPACEDIM>& stencil = DefaultType)
568 : {
569 160392 : Set::Vector ret = Set::Vector::Zero();
570 :
571 160392 : if (stencil[0] == StencilType::Central)
572 : {
573 782610 : AMREX_D_TERM(ret(0) += (dw(i + 1, j, k)(0, 0) - dw(i - 1, j, k)(0, 0)) / 2. / DX[0];,
574 : ret(1) += (dw(i + 1, j, k)(1, 0) - dw(i - 1, j, k)(1, 0)) / 2. / DX[0];,
575 : ret(2) += (dw(i + 1, j, k)(2, 0) - dw(i - 1, j, k)(2, 0)) / 2. / DX[0];)
576 : }
577 7110 : else if (stencil[0] == StencilType::Lo)
578 : {
579 17700 : AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 0) - dw(i - 1, j, k)(0, 0)) / DX[0];,
580 : ret(1) += (dw(i, j, k)(1, 0) - dw(i - 1, j, k)(1, 0)) / DX[0];,
581 : ret(2) += (dw(i, j, k)(2, 0) - dw(i - 1, j, k)(2, 0)) / DX[0];)
582 : }
583 3570 : else if (stencil[0] == StencilType::Hi)
584 : {
585 17850 : AMREX_D_TERM(ret(0) += (dw(i + 1, j, k)(0, 0) - dw(i, j, k)(0, 0)) / DX[0];,
586 : ret(1) += (dw(i + 1, j, k)(1, 0) - dw(i, j, k)(1, 0)) / DX[0];,
587 : ret(2) += (dw(i + 1, j, k)(2, 0) - dw(i, j, k)(2, 0)) / DX[0];)
588 : }
589 :
590 : #if AMREX_SPACEDIM > 1
591 160392 : if (stencil[1] == StencilType::Central)
592 : {
593 770010 : AMREX_D_TERM(ret(0) += (dw(i, j + 1, k)(0, 1) - dw(i, j - 1, k)(0, 1)) / 2. / DX[1];,
594 : ret(1) += (dw(i, j + 1, k)(1, 1) - dw(i, j - 1, k)(1, 1)) / 2. / DX[1];,
595 : ret(2) += (dw(i, j + 1, k)(2, 1) - dw(i, j - 1, k)(2, 1)) / 2. / DX[1];)
596 : }
597 9630 : else if (stencil[1] == StencilType::Lo)
598 : {
599 24000 : AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 1) - dw(i, j - 1, k)(0, 1)) / DX[1];,
600 : ret(1) += (dw(i, j, k)(1, 1) - dw(i, j - 1, k)(1, 1)) / DX[1];,
601 : ret(2) += (dw(i, j, k)(2, 1) - dw(i, j - 1, k)(2, 1)) / DX[1];)
602 : }
603 4830 : else if (stencil[1] == StencilType::Hi)
604 : {
605 24150 : AMREX_D_TERM(ret(0) += (dw(i, j + 1, k)(0, 1) - dw(i, j, k)(0, 1)) / DX[1];,
606 : ret(1) += (dw(i, j + 1, k)(1, 1) - dw(i, j, k)(1, 1)) / DX[1];,
607 : ret(2) += (dw(i, j + 1, k)(2, 1) - dw(i, j, k)(2, 1)) / DX[1];)
608 : }
609 : #endif
610 : #if AMREX_SPACEDIM > 2
611 8100 : if (stencil[2] == StencilType::Central)
612 : {
613 56700 : AMREX_D_TERM(ret(0) += (dw(i, j, k + 1)(0, 2) - dw(i, j, k - 1)(0, 2)) / 2. / DX[2];,
614 : ret(1) += (dw(i, j, k + 1)(1, 2) - dw(i, j, k - 1)(1, 2)) / 2. / DX[2];,
615 : ret(2) += (dw(i, j, k + 1)(2, 2) - dw(i, j, k - 1)(2, 2)) / 2. / DX[2];)
616 : }
617 0 : else if (stencil[2] == StencilType::Lo)
618 : {
619 0 : AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 2) - dw(i, j, k - 1)(0, 2)) / DX[2];,
620 : ret(1) += (dw(i, j, k)(1, 2) - dw(i, j, k - 1)(1, 2)) / DX[2];,
621 : ret(2) += (dw(i, j, k)(2, 2) - dw(i, j, k - 1)(2, 2)) / DX[2];)
622 : }
623 0 : else if (stencil[2] == StencilType::Hi)
624 : {
625 0 : AMREX_D_TERM(ret(0) += (dw(i, j, k + 1)(0, 2) - dw(i, j, k)(0, 2)) / DX[2];,
626 : ret(1) += (dw(i, j, k + 1)(1, 2) - dw(i, j, k)(1, 2)) / DX[2];,
627 : ret(2) += (dw(i, j, k + 1)(2, 2) - dw(i, j, k)(2, 2)) / DX[2];)
628 : }
629 : #endif
630 160392 : return ret;
631 : }
632 :
633 :
634 : [[nodiscard]] AMREX_FORCE_INLINE
635 : Set::Scalar
636 : Divergence(const amrex::Array4<const Set::Scalar>& f,
637 : const int& i, const int& j, const int& k, const int& m,
638 : const Set::Scalar dx[AMREX_SPACEDIM],
639 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
640 : {
641 : Set::Scalar ret;
642 : ret = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, m, dx, stencil));
643 : #if AMREX_SPACEDIM > 1
644 : ret += (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, m, dx, stencil));
645 : #endif
646 : #if AMREX_SPACEDIM > 2
647 : ret += (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, m, dx, stencil));
648 : #endif
649 : return ret;
650 : }
651 :
652 :
653 : [[nodiscard]] AMREX_FORCE_INLINE
654 : Set::Vector
655 : Gradient(const amrex::Array4<const Set::Scalar>& f,
656 : const int& i, const int& j, const int& k, const int& m,
657 : const Set::Scalar dx[AMREX_SPACEDIM],
658 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
659 : {
660 38417648 : Set::Vector ret;
661 38417648 : ret(0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, m, dx, stencil));
662 : #if AMREX_SPACEDIM > 1
663 38417648 : ret(1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, m, dx, stencil));
664 : #if AMREX_SPACEDIM > 2
665 78784 : ret(2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, m, dx, stencil));
666 : #endif
667 : #endif
668 38417648 : return ret;
669 : }
670 :
671 : [[nodiscard]] AMREX_FORCE_INLINE
672 : Set::Vector
673 : CellGradientOnNode(const amrex::Array4<const Set::Scalar>& f,
674 : const int& i, const int& j, const int& k, const int& m,
675 : const Set::Scalar dx[AMREX_SPACEDIM])
676 : {
677 4358180 : Set::Vector ret;
678 : #if AMREX_SPACEDIM == 1
679 : ret(0) = (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0];
680 : #elif AMREX_SPACEDIM == 2
681 21790850 : ret(0) = 0.5 * (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m)) / dx[0];
682 21790850 : ret(1) = 0.5 * (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m)) / dx[1];
683 : #elif AMREX_SPACEDIM == 3
684 0 : ret(0) = 0.25 * (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k - 1, m) - f(i - 1, j - 1, k - 1, m)) / dx[0];
685 0 : ret(1) = 0.25 * (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i, j - 1, k - 1, m) + f(i - 1, j, k - 1, m) - f(i - 1, j - 1, k - 1, m)) / dx[1];
686 0 : ret(2) = 0.25 * (f(i, j, k, m) - f(i, j, k - 1, m) + f(i - 1, j, k, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k, m) - f(i, j - 1, k - 1, m) + f(i - 1, j - 1, k, m) - f(i - 1, j - 1, k - 1, m)) / dx[2];
687 : #endif
688 4358180 : return ret;
689 : }
690 :
691 :
692 : template<class T>
693 : [[nodiscard]] AMREX_FORCE_INLINE
694 : std::array<T, AMREX_SPACEDIM>
695 : CellGradientOnNode(const amrex::Array4<const T>& f,
696 : const int& i, const int& j, const int& k, const int& m,
697 : const Set::Scalar dx[AMREX_SPACEDIM])
698 : {
699 : std::array<T, AMREX_SPACEDIM> ret;
700 : Util::Abort(INFO);
701 : #if AMREX_SPACEDIM == 1
702 : ret[0] = (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0];
703 : #elif AMREX_SPACEDIM == 2
704 : ret[0] = (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m)) * 0.5 / dx[0];
705 : ret[1] = (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m)) * 0.5 / dx[1];
706 : #elif AMREX_SPACEDIM == 3
707 : ret[0] = (f(i, j, k, m) - f(i - 1, j, k, m) + f(i, j - 1, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k - 1, m) - f(i - 1, j - 1, k - 1, m)) * 0.25 / dx[0];
708 : ret[1] = (f(i, j, k, m) - f(i, j - 1, k, m) + f(i - 1, j, k, m) - f(i - 1, j - 1, k, m) + f(i, j, k - 1, m) - f(i, j - 1, k - 1, m) + f(i - 1, j, k - 1, m) - f(i - 1, j - 1, k - 1, m)) * 0.25 / dx[1];
709 : ret[2] = (f(i, j, k, m) - f(i, j, k - 1, m) + f(i - 1, j, k, m) - f(i - 1, j, k - 1, m) + f(i, j - 1, k, m) - f(i, j - 1, k - 1, m) + f(i - 1, j - 1, k, m) - f(i - 1, j - 1, k - 1, m)) * 0.25 / dx[2];
710 : #endif
711 : return ret;
712 : }
713 :
714 :
715 :
716 : [[nodiscard]] AMREX_FORCE_INLINE
717 : Set::Matrix
718 : Gradient(const amrex::Array4<const Set::Scalar>& f,
719 : const int& i, const int& j, const int& k,
720 : const Set::Scalar dx[AMREX_SPACEDIM],
721 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
722 : {
723 0 : Set::Matrix ret;
724 0 : ret(0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
725 : #if AMREX_SPACEDIM > 1
726 0 : ret(0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
727 0 : ret(1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
728 0 : ret(1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
729 : #if AMREX_SPACEDIM > 2
730 0 : ret(0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil));
731 0 : ret(2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
732 0 : ret(1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 1, dx, stencil));
733 0 : ret(2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
734 0 : ret(2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 2, dx, stencil));
735 : #endif
736 : #endif
737 0 : return ret;
738 : }
739 :
740 : [[nodiscard]] AMREX_FORCE_INLINE
741 : Set::Matrix
742 : Gradient(const amrex::Array4<const Set::Vector>& f,
743 : const int& i, const int& j, const int& k,
744 : const Set::Scalar dx[AMREX_SPACEDIM],
745 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
746 : {
747 446379 : Set::Matrix ret;
748 :
749 : #if AMREX_SPACEDIM > 0
750 892758 : ret.col(0) = Numeric::Stencil<Set::Vector, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil);
751 : #endif
752 : #if AMREX_SPACEDIM > 1
753 892758 : ret.col(1) = Numeric::Stencil<Set::Vector, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil);
754 : #endif
755 : #if AMREX_SPACEDIM > 2
756 208800 : ret.col(2) = Numeric::Stencil<Set::Vector, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil);
757 : #endif
758 :
759 446379 : return ret;
760 : }
761 :
762 : [[nodiscard]] AMREX_FORCE_INLINE
763 : std::pair<Set::Vector,Set::Matrix>
764 : GradientSplit( const amrex::Array4<const Set::Vector>& f,
765 : const int& i, const int& j, const int& k,
766 : const Set::Scalar dx[AMREX_SPACEDIM],
767 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
768 : {
769 0 : Set::Vector diag;
770 0 : Set::Matrix offdiag;
771 :
772 0 : std::pair<Set::Scalar,Set::Vector> ret;
773 : #if AMREX_SPACEDIM > 0
774 0 : ret = Numeric::Stencil<Set::Vector, 1, 0, 0>::Dsplit(f, i, j, k, 0, dx, stencil);
775 0 : diag(0) = ret.first;
776 0 : offdiag.col(0) = ret.second;
777 : #endif
778 : #if AMREX_SPACEDIM > 1
779 0 : ret = Numeric::Stencil<Set::Vector, 0, 1, 0>::Dsplit(f, i, j, k, 0, dx, stencil);
780 0 : diag(1) = ret.first;
781 0 : offdiag.col(1) = ret.second;
782 : #endif
783 : #if AMREX_SPACEDIM > 2
784 0 : ret = Numeric::Stencil<Set::Vector, 0, 0, 1>::Dsplit(f, i, j, k, 0, dx, stencil);
785 0 : diag(2) = ret.first;
786 0 : offdiag.col(2) = ret.second;
787 : #endif
788 0 : return std::make_pair(diag,offdiag);
789 : }
790 :
791 :
792 : [[nodiscard]] AMREX_FORCE_INLINE
793 : Set::Matrix
794 : NodeGradientOnCell(const amrex::Array4<const Set::Vector>& f,
795 : const int& i, const int& j, const int& k,
796 : const Set::Scalar dx[AMREX_SPACEDIM])
797 : {
798 : Set::Matrix ret;
799 : #if AMREX_SPACEDIM == 1
800 : ret.col(0) = (f(i + 1, j, k) - f(i, j, k)) / dx[0];
801 : #elif AMREX_SPACEDIM == 2
802 : ret.col(0) = (f(i + 1, j, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i, j + 1, k)) * 0.5 / dx[0];
803 : ret.col(1) = (f(i, j + 1, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i + 1, j, k)) * 0.5 / dx[1];
804 : #elif AMREX_SPACEDIM == 3
805 : ret.col(0) = (f(i + 1, j, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i, j + 1, k) + f(i + 1, j, k + 1) - f(i, j, k + 1) + f(i + 1, j + 1, k + 1) - f(i, j + 1, k + 1)) * 0.25 / dx[0];
806 : ret.col(1) = (f(i, j + 1, k) - f(i, j, k) + f(i + 1, j + 1, k) - f(i + 1, j, k) + f(i, j + 1, k + 1) - f(i, j, k + 1) + f(i + 1, j + 1, k + 1) - f(i + 1, j, k + 1)) * 0.25 / dx[1];
807 : ret.col(2) = (f(i, j, k + 1) - f(i, j, k) + f(i, j + 1, k + 1) - f(i, j + 1, k) + f(i + 1, j, k + 1) - f(i + 1, j, k) + f(i + 1, j + 1, k + 1) - f(i + 1, j + 1, k)) * 0.25 / dx[2];
808 : #endif
809 : return ret;
810 : }
811 :
812 : [[nodiscard]] AMREX_FORCE_INLINE
813 : Set::Matrix3
814 : Gradient(const amrex::Array4<const Set::Matrix>& f,
815 : const int& i, const int& j, const int& k,
816 : const Set::Scalar dx[AMREX_SPACEDIM],
817 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
818 : {
819 : Set::Matrix3 ret;
820 :
821 : #if AMREX_SPACEDIM > 0
822 : ret[0] = Numeric::Stencil<Set::Matrix, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil);
823 : #endif
824 : #if AMREX_SPACEDIM > 1
825 : ret[1] = Numeric::Stencil<Set::Matrix, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil);
826 : #endif
827 : #if AMREX_SPACEDIM > 2
828 : ret[2] = Numeric::Stencil<Set::Matrix, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil);
829 : #endif
830 :
831 : return ret;
832 : }
833 : [[nodiscard]] AMREX_FORCE_INLINE
834 : Set::Matrix3
835 : NodeGradientOnCell(const amrex::Array4<const Set::Matrix>& f,
836 : const int& i, const int& j, const int& k,
837 : const Set::Scalar dx[AMREX_SPACEDIM])
838 : {
839 371408 : Set::Matrix3 ret;
840 :
841 : #if AMREX_SPACEDIM > 0
842 1485628 : ret[0] = (f(i+1,j,k) - f(i,j,k)) / dx[0];
843 : #endif
844 : #if AMREX_SPACEDIM > 1
845 1485628 : ret[1] = (f(i,j+1,k) - f(i,j,k)) / dx[1];
846 : #endif
847 : #if AMREX_SPACEDIM > 2
848 0 : ret[2] = (f(i,j,k+1) - f(i,j,k)) / dx[2];
849 : #endif
850 :
851 371408 : return ret;
852 : }
853 :
854 :
855 : [[nodiscard]] AMREX_FORCE_INLINE
856 : Set::Matrix3
857 : MatrixGradient(const amrex::Array4<const Set::Scalar>& f,
858 : const int& i, const int& j, const int& k,
859 : const Set::Scalar dx[AMREX_SPACEDIM],
860 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
861 : {
862 : Set::Matrix3 ret;
863 : #if AMREX_SPACEDIM == 1
864 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
865 : #elif AMREX_SPACEDIM == 2
866 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
867 : ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
868 : ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
869 : ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
870 : ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
871 : ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
872 : ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
873 : ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
874 : #elif AMREX_SPACEDIM == 3
875 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
876 : ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
877 : ret[0](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
878 : ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
879 : ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 4, dx, stencil));
880 : ret[1](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 5, dx, stencil));
881 : ret[2](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 6, dx, stencil));
882 : ret[2](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 7, dx, stencil));
883 : ret[2](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 8, dx, stencil));
884 : ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
885 : ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
886 : ret[0](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
887 : ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
888 : ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 4, dx, stencil));
889 : ret[1](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 5, dx, stencil));
890 : ret[2](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 6, dx, stencil));
891 : ret[2](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 7, dx, stencil));
892 : ret[2](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 8, dx, stencil));
893 : ret[0](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil));
894 : ret[0](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 1, dx, stencil));
895 : ret[0](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 2, dx, stencil));
896 : ret[1](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 3, dx, stencil));
897 : ret[1](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 4, dx, stencil));
898 : ret[1](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 5, dx, stencil));
899 : ret[2](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 6, dx, stencil));
900 : ret[2](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 7, dx, stencil));
901 : ret[2](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 8, dx, stencil));
902 : #endif
903 : return ret;
904 : }
905 :
906 :
907 : [[nodiscard]] AMREX_FORCE_INLINE
908 : Set::Matrix
909 : Hessian(const amrex::Array4<const Set::Scalar>& f,
910 : const int& i, const int& j, const int& k, const int& m,
911 : const Set::Scalar dx[AMREX_SPACEDIM],
912 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType
913 : )
914 : {
915 2401960 : Set::Matrix ret;
916 2401960 : ret(0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx, stencil));
917 : #if AMREX_SPACEDIM > 1
918 2401960 : ret(1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, m, dx, stencil));
919 2401960 : ret(0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, m, dx, stencil));
920 2401960 : ret(1, 0) = ret(0, 1);
921 : #if AMREX_SPACEDIM > 2
922 0 : ret(2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, m, dx, stencil));
923 0 : ret(1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, m, dx, stencil));
924 0 : ret(2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, m, dx, stencil));
925 0 : ret(2, 1) = ret(1, 2);
926 0 : ret(0, 2) = ret(2, 0);
927 : #endif
928 : #endif
929 2401960 : return ret;
930 : }
931 :
932 : [[nodiscard]] AMREX_FORCE_INLINE
933 : Set::Matrix3
934 : Hessian(const amrex::Array4<const Set::Scalar>& f,
935 : const int& i, const int& j, const int& k,
936 : const Set::Scalar DX[AMREX_SPACEDIM],
937 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
938 : {
939 : Set::Matrix3 ret;
940 : // 1D
941 : #if AMREX_SPACEDIM>0
942 : ret(0, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 0, DX, stencil));
943 : #endif
944 : // 2D
945 : #if AMREX_SPACEDIM>1
946 : ret(0, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 0, DX, stencil));
947 : ret(0, 1, 0) = ret(0, 0, 1);
948 : ret(0, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 0, DX, stencil));
949 : ret(1, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 1, DX, stencil));
950 : ret(1, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 1, DX, stencil));
951 : ret(1, 1, 0) = ret(1, 0, 1);
952 : ret(1, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 1, DX, stencil));
953 : #endif
954 : // 3D
955 : #if AMREX_SPACEDIM>2
956 : ret(0, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 0, DX, stencil));
957 : ret(0, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 0, DX, stencil));
958 : ret(0, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 0, DX, stencil));
959 : ret(0, 2, 0) = ret(0, 0, 2);
960 : ret(0, 2, 1) = ret(0, 1, 2);;
961 : ret(1, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 1, DX, stencil));
962 : ret(1, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 1, DX, stencil));
963 : ret(1, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 1, DX, stencil));
964 : ret(1, 2, 0) = ret(1, 0, 2);
965 : ret(1, 2, 1) = ret(1, 1, 2);;
966 : ret(2, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 2, DX, stencil));
967 : ret(2, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 2, DX, stencil));
968 : ret(2, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 2, DX, stencil));
969 : ret(2, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 2, DX, stencil));
970 : ret(2, 1, 0) = ret(2, 0, 1);
971 : ret(2, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 2, DX, stencil));
972 : ret(2, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 2, DX, stencil));
973 : ret(2, 2, 0) = ret(2, 0, 2);
974 : ret(2, 2, 1) = ret(2, 1, 2);;
975 : #endif
976 : return ret;
977 : }
978 :
979 :
980 : // Returns Hessian of a vector field.
981 : // Return value: ret[i](j,k) = ret_{i,jk}
982 : [[nodiscard]] AMREX_FORCE_INLINE
983 : Set::Matrix3
984 : Hessian(const amrex::Array4<const Set::Vector>& f,
985 : const int& i, const int& j, const int& k,
986 : const Set::Scalar dx[AMREX_SPACEDIM],
987 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
988 : {
989 : Set::Matrix3 ret;
990 :
991 : #if AMREX_SPACEDIM>0
992 : Set::Vector f_11 = Numeric::Stencil<Set::Vector, 2, 0, 0>::D(f, i, j, k, 0, dx, stencil);
993 : ret[0](0, 0) = f_11(0);
994 : #endif
995 : #if AMREX_SPACEDIM>1
996 : ret[1](0, 0) = f_11(1);
997 : Set::Vector f_12 = Numeric::Stencil<Set::Vector, 1, 1, 0>::D(f, i, j, k, 0, dx, stencil);
998 : ret[0](1, 0) = ret[0](0, 1) = f_12(0);
999 : ret[1](1, 0) = ret[1](0, 1) = f_12(1);
1000 : Set::Vector f_22 = Numeric::Stencil<Set::Vector, 0, 2, 0>::D(f, i, j, k, 0, dx, stencil);
1001 : ret[0](1, 1) = f_22(0);
1002 : ret[1](1, 1) = f_22(1);
1003 : #endif
1004 : #if AMREX_SPACEDIM>2
1005 : ret[2](0, 0) = f_11(2);
1006 : ret[2](1, 0) = ret[2](0, 1) = f_12(2);
1007 : Set::Vector f_13 = Numeric::Stencil<Set::Vector, 1, 0, 1>::D(f, i, j, k, 0, dx, stencil);
1008 : ret[0](2, 0) = ret[0](0, 2) = f_13(0);
1009 : ret[1](2, 0) = ret[1](0, 2) = f_13(1);
1010 : ret[2](2, 0) = ret[2](0, 2) = f_13(2);
1011 : ret[2](1, 1) = f_22(2);
1012 : Set::Vector f_23 = Numeric::Stencil<Set::Vector, 0, 1, 1>::D(f, i, j, k, 0, dx, stencil);
1013 : ret[0](1, 2) = ret[0](2, 1) = f_23(0);
1014 : ret[1](1, 2) = ret[1](2, 1) = f_23(1);
1015 : ret[2](1, 2) = ret[2](2, 1) = f_23(2);
1016 : Set::Vector f_33 = Numeric::Stencil<Set::Vector, 0, 0, 2>::D(f, i, j, k, 0, dx, stencil);
1017 : ret[0](2, 2) = f_33(0);
1018 : ret[1](2, 2) = f_33(1);
1019 : ret[2](2, 2) = f_33(2);
1020 : #endif
1021 : return ret;
1022 : }
1023 :
1024 :
1025 : [[nodiscard]] AMREX_FORCE_INLINE
1026 : Set::Matrix
1027 : FieldToMatrix(const amrex::Array4<const Set::Scalar>& f,
1028 : const int& i, const int& j, const int& k)
1029 : {
1030 0 : Set::Matrix ret;
1031 : #if AMREX_SPACEDIM == 1
1032 : ret(0, 0) = f(i, j, k, 0);
1033 :
1034 : #elif AMREX_SPACEDIM == 2
1035 0 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1036 0 : ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1037 :
1038 : #elif AMREX_SPACEDIM == 3
1039 0 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1); ret(0, 2) = f(i, j, k, 2);
1040 0 : ret(1, 0) = f(i, j, k, 3); ret(1, 1) = f(i, j, k, 4); ret(1, 2) = f(i, j, k, 5);
1041 0 : ret(2, 0) = f(i, j, k, 6); ret(2, 1) = f(i, j, k, 7); ret(2, 2) = f(i, j, k, 8);
1042 : #endif
1043 :
1044 0 : return ret;
1045 : }
1046 :
1047 : [[nodiscard]] AMREX_FORCE_INLINE
1048 : Set::Matrix
1049 : FieldToMatrix(const amrex::Array4<Set::Scalar>& f,
1050 : const int& i, const int& j, const int& k)
1051 : {
1052 : Set::Matrix ret;
1053 : #if AMREX_SPACEDIM == 1
1054 : ret(0, 0) = f(i, j, k, 0);
1055 :
1056 : #elif AMREX_SPACEDIM == 2
1057 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1058 : ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1059 :
1060 : #elif AMREX_SPACEDIM == 3
1061 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1); ret(0, 2) = f(i, j, k, 2);
1062 : ret(1, 0) = f(i, j, k, 3); ret(1, 1) = f(i, j, k, 4); ret(1, 2) = f(i, j, k, 5);
1063 : ret(2, 0) = f(i, j, k, 6); ret(2, 1) = f(i, j, k, 7); ret(2, 2) = f(i, j, k, 8);
1064 : #endif
1065 :
1066 : return ret;
1067 : }
1068 :
1069 : [[nodiscard]] AMREX_FORCE_INLINE
1070 : Set::Vector
1071 : FieldToVector(const amrex::Array4<const Set::Scalar>& f,
1072 : const int& i, const int& j, const int& k)
1073 : {
1074 : Set::Vector ret;
1075 : ret(0) = f(i, j, k, 0);
1076 : #if AMREX_SPACEDIM > 1
1077 : ret(1) = f(i, j, k, 1);
1078 : #if AMREX_SPACEDIM > 2
1079 : ret(2) = f(i, j, k, 2);
1080 : #endif
1081 : #endif
1082 : return ret;
1083 : }
1084 :
1085 : [[nodiscard]] AMREX_FORCE_INLINE
1086 : Set::Vector
1087 : FieldToVector(const amrex::Array4<Set::Scalar>& f,
1088 : const int& i, const int& j, const int& k)
1089 : {
1090 : Set::Vector ret;
1091 : ret(0) = f(i, j, k, 0);
1092 : #if AMREX_SPACEDIM > 1
1093 : ret(1) = f(i, j, k, 1);
1094 : #if AMREX_SPACEDIM > 2
1095 : ret(2) = f(i, j, k, 2);
1096 : #endif
1097 : #endif
1098 : return ret;
1099 : }
1100 :
1101 : AMREX_FORCE_INLINE
1102 : void
1103 : MatrixToField(const amrex::Array4<Set::Scalar>& f,
1104 : const int& i, const int& j, const int& k,
1105 : Set::Matrix matrix)
1106 : {
1107 : #if AMREX_SPACEDIM == 1
1108 : f(i, j, k, 0) = matrix(0, 0);
1109 : #elif AMREX_SPACEDIM == 2
1110 : f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1);
1111 : f(i, j, k, 2) = matrix(1, 0); f(i, j, k, 3) = matrix(1, 1);
1112 : #elif AMREX_SPACEDIM == 3
1113 : f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1); f(i, j, k, 2) = matrix(0, 2);
1114 : f(i, j, k, 3) = matrix(1, 0); f(i, j, k, 4) = matrix(1, 1); f(i, j, k, 5) = matrix(1, 2);
1115 : f(i, j, k, 6) = matrix(2, 0); f(i, j, k, 7) = matrix(2, 1); f(i, j, k, 8) = matrix(2, 2);
1116 : #endif
1117 : }
1118 :
1119 : AMREX_FORCE_INLINE
1120 : void
1121 : VectorToField(const amrex::Array4<Set::Scalar>& f,
1122 : const int& i, const int& j, const int& k,
1123 : Set::Vector vector)
1124 : {
1125 : f(i, j, k, 0) = vector(0);
1126 : #if AMREX_SPACEDIM > 1
1127 : f(i, j, k, 1) = vector(1);
1128 : #if AMREX_SPACEDIM > 2
1129 : f(i, j, k, 2) = vector(2);
1130 : #endif
1131 : #endif
1132 : }
1133 :
1134 : template<int index, int SYM>
1135 : [[nodiscard]]
1136 : Set::Matrix3
1137 : Divergence(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, SYM>>&,
1138 : const int, const int, const int, const Set::Scalar[AMREX_SPACEDIM],
1139 : std::array<StencilType, AMREX_SPACEDIM> /*stecil*/ = DefaultType)
1140 : {
1141 : Util::Abort(INFO, "Not implemented yet"); return Set::Matrix3::Zero();
1142 : }
1143 :
1144 : template<>
1145 : [[nodiscard]] AMREX_FORCE_INLINE
1146 : Set::Matrix3 Divergence<2, Set::Sym::Isotropic>(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>>& C,
1147 : const int i, const int j, const int k, const Set::Scalar dx[AMREX_SPACEDIM],
1148 : std::array<StencilType, AMREX_SPACEDIM> stencil)
1149 : {
1150 : Set::Matrix3 ret;
1151 :
1152 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCx =
1153 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 1, 0, 0>::D(C, i, j, k, 0, dx, stencil);
1154 : #if AMREX_SPACEDIM>1
1155 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCy =
1156 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 1, 0>::D(C, i, j, k, 0, dx, stencil);
1157 : #endif
1158 : #if AMREX_SPACEDIM>2
1159 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCz =
1160 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 0, 1>::D(C, i, j, k, 0, dx, stencil);
1161 : #endif
1162 :
1163 : for (int i = 0; i < AMREX_SPACEDIM; i++)
1164 : for (int k = 0; k < AMREX_SPACEDIM; k++)
1165 : for (int l = 0; l < AMREX_SPACEDIM; l++)
1166 : {
1167 : ret(i, k, l) = 0.0;
1168 : ret(i, k, l) = gradCx(i, 0, k, l);
1169 : #if AMREX_SPACEDIM>1
1170 : ret(i, k, l) = gradCy(i, 1, k, l);
1171 : #endif
1172 : #if AMREX_SPACEDIM>2
1173 : ret(i, k, l) = gradCz(i, 2, k, l);
1174 : #endif
1175 : }
1176 : return ret;
1177 : }
1178 :
1179 :
1180 :
1181 :
1182 : template<int dim>
1183 : [[nodiscard]] AMREX_FORCE_INLINE
1184 : Set::Matrix4<dim, Set::Sym::Full>
1185 : DoubleHessian(const amrex::Array4<const Set::Scalar>& f,
1186 : const int& i, const int& j, const int& k, const int& m,
1187 : const Set::Scalar dx[AMREX_SPACEDIM]);
1188 :
1189 : template<>
1190 : [[nodiscard]] AMREX_FORCE_INLINE
1191 : Set::Matrix4<2, Set::Sym::Full>
1192 : DoubleHessian<2>(const amrex::Array4<const Set::Scalar>& f,
1193 : const int& i, const int& j, const int& k, const int& m,
1194 : const Set::Scalar dx[AMREX_SPACEDIM])
1195 : {
1196 1362860 : Set::Matrix4<2, Set::Sym::Full> ret;
1197 : // [0,0,0,0]
1198 2725720 : ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1199 : // [0, 0, 0, 1]
1200 2725720 : ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1201 : // [0, 0, 1, 1]
1202 2725720 : ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1203 : // [0, 1, 1, 1]
1204 2725720 : ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1205 : // [1, 1, 1, 1]
1206 1362860 : ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1207 1362860 : return ret;
1208 : }
1209 :
1210 : template<>
1211 : [[nodiscard]] AMREX_FORCE_INLINE
1212 : Set::Matrix4<3, Set::Sym::Full>
1213 : DoubleHessian<3>(const amrex::Array4<const Set::Scalar>& f,
1214 : const int& i, const int& j, const int& k, const int& m,
1215 : const Set::Scalar dx[AMREX_SPACEDIM])
1216 : {
1217 0 : Set::Matrix4<3, Set::Sym::Full> ret;
1218 : // [0,0,0,0]
1219 0 : ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1220 : // [0, 0, 0, 1]
1221 0 : ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1222 : // [0, 0, 0, 2]
1223 0 : ret(0, 0, 0, 2) = Stencil<Set::Scalar, 3, 0, 1>::D(f, i, j, k, m, dx);
1224 : // [0, 0, 1, 1]
1225 0 : ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1226 : // [0, 0, 1, 2]
1227 0 : ret(0, 0, 1, 2) = Stencil<Set::Scalar, 2, 1, 1>::D(f, i, j, k, m, dx);
1228 : // [0, 0, 2, 2]
1229 0 : ret(0, 0, 2, 2) = Stencil<Set::Scalar, 2, 0, 2>::D(f, i, j, k, m, dx);
1230 : // [0, 1, 1, 1]
1231 0 : ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1232 : // [0, 1, 1, 2]
1233 0 : ret(0, 1, 1, 2) = Stencil<Set::Scalar, 1, 2, 1>::D(f, i, j, k, m, dx);
1234 : // [0, 1, 2, 2]
1235 0 : ret(0, 1, 2, 2) = Stencil<Set::Scalar, 1, 1, 2>::D(f, i, j, k, m, dx);
1236 : // [0, 2, 2, 2]
1237 0 : ret(0, 2, 2, 2) = Stencil<Set::Scalar, 1, 0, 3>::D(f, i, j, k, m, dx);
1238 : // [1, 1, 1, 1]
1239 0 : ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1240 : // [1, 1, 1, 2]
1241 0 : ret(1, 1, 1, 2) = Stencil<Set::Scalar, 0, 3, 1>::D(f, i, j, k, m, dx);
1242 : // [1, 1, 2, 2]
1243 0 : ret(1, 1, 2, 2) = Stencil<Set::Scalar, 0, 2, 2>::D(f, i, j, k, m, dx);
1244 : // [1, 2, 2, 2]
1245 0 : ret(1, 2, 2, 2) = Stencil<Set::Scalar, 0, 1, 3>::D(f, i, j, k, m, dx);
1246 : // [2, 2, 2, 2]
1247 0 : ret(2, 2, 2, 2) = Stencil<Set::Scalar, 0, 0, 4>::D(f, i, j, k, m, dx);
1248 0 : return ret;
1249 : }
1250 :
1251 : struct Interpolate
1252 : {
1253 : public:
1254 : template<class T>
1255 : [[nodiscard]] AMREX_FORCE_INLINE
1256 : static T CellToNodeAverage(const amrex::Array4<const T>& f,
1257 : const int& i, const int& j, const int& k, const int& m,
1258 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1259 : {
1260 184516 : int AMREX_D_DECL(
1261 : ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1262 : jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1263 : klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1264 :
1265 738064 : return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1266 : ,
1267 : +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1268 : ,
1269 : +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1270 : + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1271 184516 : )) * fac;
1272 : }
1273 : template<class T>
1274 : [[nodiscard]] AMREX_FORCE_INLINE
1275 : static T CellToNodeAverage(const amrex::Array4<T>& f,
1276 : const int& i, const int& j, const int& k, const int& m,
1277 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1278 : {
1279 4620090 : int AMREX_D_DECL(
1280 : ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1281 : jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1282 : klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1283 :
1284 18480400 : return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1285 : ,
1286 : +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1287 : ,
1288 : +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1289 : + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1290 4620090 : )) * fac;
1291 : }
1292 : template<class T>
1293 : [[nodiscard]] AMREX_FORCE_INLINE
1294 : static T NodeToCellAverage(const amrex::Array4<const T>& f,
1295 : const int& i, const int& j, const int& k, const int& m)
1296 : {
1297 21110164 : return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1298 : ,
1299 : +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1300 : ,
1301 : +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1302 : + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1303 5277536 : )) * fac;
1304 : }
1305 : template<class T>
1306 : [[nodiscard]] AMREX_FORCE_INLINE
1307 : static T NodeToCellAverage(const amrex::Array4<T>& f,
1308 : const int& i, const int& j, const int& k, const int& m)
1309 : {
1310 : return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1311 : ,
1312 : +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1313 : ,
1314 : +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1315 : + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1316 : )) * fac;
1317 : }
1318 : constexpr static Set::Scalar fac = AMREX_D_PICK(0.5, 0.25, 0.125);
1319 : };
1320 :
1321 : }
1322 : #endif
|
