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 28624438 : 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 28624438 : 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 113832460 : if (stencil[0] == StencilType::Lo)
55 30667068 : return (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0]; // 1st order stencil
56 103610122 : else if (stencil[0] == StencilType::Hi)
57 30666978 : return (f(i + 1, j, k, m) - f(i, j, k, m)) / dx[0]; // 1st order stencil
58 : else
59 317266422 : 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 113832460 : if (stencil[1] == StencilType::Lo)
96 29696210 : return (f(i, j, k, m) - f(i, j - 1, k, m)) / dx[1];
97 103933776 : else if (stencil[1] == StencilType::Hi)
98 29696120 : return (f(i, j + 1, k, m) - f(i, j, k, m)) / dx[1];
99 : else
100 319207776 : 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 363407030 : if (stencil[0] == StencilType::Central)
182 1453461348 : 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 363444056 : if (stencil[1] == StencilType::Central)
198 1453607052 : 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 : std::array<StencilType, AMREX_SPACEDIM>& stencil = DefaultType)
534 : {
535 332329844 : Set::Scalar ret = 0.0;
536 332329844 : ret += (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx, stencil));
537 : #if AMREX_SPACEDIM > 1
538 332329844 : ret += (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, m, dx, stencil));
539 : #if AMREX_SPACEDIM > 2
540 835328 : ret += (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, m, dx, stencil));
541 : #endif
542 : #endif
543 332329844 : return ret;
544 : }
545 :
546 : [[nodiscard]] AMREX_FORCE_INLINE
547 : Set::Vector
548 : Laplacian(const amrex::Array4<const Set::Vector>& f,
549 : const int& i, const int& j, const int& k,
550 : const Set::Scalar dx[AMREX_SPACEDIM])
551 : {
552 0 : Set::Vector ret = Set::Vector::Zero();
553 0 : ret += (Numeric::Stencil<Set::Vector, 2, 0, 0>::D(f, i, j, k, 0, dx));
554 : #if AMREX_SPACEDIM > 1
555 0 : ret += (Numeric::Stencil<Set::Vector, 0, 2, 0>::D(f, i, j, k, 0, dx));
556 : #if AMREX_SPACEDIM > 2
557 0 : ret += (Numeric::Stencil<Set::Vector, 0, 0, 2>::D(f, i, j, k, 0, dx));
558 : #endif
559 : #endif
560 0 : return ret;
561 : }
562 :
563 : [[nodiscard]] AMREX_FORCE_INLINE
564 : Set::Vector
565 : Divergence(const amrex::Array4<const Set::Matrix>& dw,
566 : const int& i, const int& j, const int& k,
567 : const Set::Scalar DX[AMREX_SPACEDIM],
568 : std::array<StencilType, AMREX_SPACEDIM>& stencil = DefaultType)
569 : {
570 160392 : Set::Vector ret = Set::Vector::Zero();
571 :
572 160392 : if (stencil[0] == StencilType::Central)
573 : {
574 782610 : AMREX_D_TERM(ret(0) += (dw(i + 1, j, k)(0, 0) - dw(i - 1, j, k)(0, 0)) / 2. / DX[0];,
575 : ret(1) += (dw(i + 1, j, k)(1, 0) - dw(i - 1, j, k)(1, 0)) / 2. / DX[0];,
576 : ret(2) += (dw(i + 1, j, k)(2, 0) - dw(i - 1, j, k)(2, 0)) / 2. / DX[0];)
577 : }
578 7110 : else if (stencil[0] == StencilType::Lo)
579 : {
580 17700 : AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 0) - dw(i - 1, j, k)(0, 0)) / DX[0];,
581 : ret(1) += (dw(i, j, k)(1, 0) - dw(i - 1, j, k)(1, 0)) / DX[0];,
582 : ret(2) += (dw(i, j, k)(2, 0) - dw(i - 1, j, k)(2, 0)) / DX[0];)
583 : }
584 3570 : else if (stencil[0] == StencilType::Hi)
585 : {
586 17850 : AMREX_D_TERM(ret(0) += (dw(i + 1, j, k)(0, 0) - dw(i, j, k)(0, 0)) / DX[0];,
587 : ret(1) += (dw(i + 1, j, k)(1, 0) - dw(i, j, k)(1, 0)) / DX[0];,
588 : ret(2) += (dw(i + 1, j, k)(2, 0) - dw(i, j, k)(2, 0)) / DX[0];)
589 : }
590 :
591 : #if AMREX_SPACEDIM > 1
592 160392 : if (stencil[1] == StencilType::Central)
593 : {
594 770010 : AMREX_D_TERM(ret(0) += (dw(i, j + 1, k)(0, 1) - dw(i, j - 1, k)(0, 1)) / 2. / DX[1];,
595 : ret(1) += (dw(i, j + 1, k)(1, 1) - dw(i, j - 1, k)(1, 1)) / 2. / DX[1];,
596 : ret(2) += (dw(i, j + 1, k)(2, 1) - dw(i, j - 1, k)(2, 1)) / 2. / DX[1];)
597 : }
598 9630 : else if (stencil[1] == StencilType::Lo)
599 : {
600 24000 : AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 1) - dw(i, j - 1, k)(0, 1)) / DX[1];,
601 : ret(1) += (dw(i, j, k)(1, 1) - dw(i, j - 1, k)(1, 1)) / DX[1];,
602 : ret(2) += (dw(i, j, k)(2, 1) - dw(i, j - 1, k)(2, 1)) / DX[1];)
603 : }
604 4830 : else if (stencil[1] == StencilType::Hi)
605 : {
606 24150 : AMREX_D_TERM(ret(0) += (dw(i, j + 1, k)(0, 1) - dw(i, j, k)(0, 1)) / DX[1];,
607 : ret(1) += (dw(i, j + 1, k)(1, 1) - dw(i, j, k)(1, 1)) / DX[1];,
608 : ret(2) += (dw(i, j + 1, k)(2, 1) - dw(i, j, k)(2, 1)) / DX[1];)
609 : }
610 : #endif
611 : #if AMREX_SPACEDIM > 2
612 8100 : if (stencil[2] == StencilType::Central)
613 : {
614 56700 : AMREX_D_TERM(ret(0) += (dw(i, j, k + 1)(0, 2) - dw(i, j, k - 1)(0, 2)) / 2. / DX[2];,
615 : ret(1) += (dw(i, j, k + 1)(1, 2) - dw(i, j, k - 1)(1, 2)) / 2. / DX[2];,
616 : ret(2) += (dw(i, j, k + 1)(2, 2) - dw(i, j, k - 1)(2, 2)) / 2. / DX[2];)
617 : }
618 0 : else if (stencil[2] == StencilType::Lo)
619 : {
620 0 : AMREX_D_TERM(ret(0) += (dw(i, j, k)(0, 2) - dw(i, j, k - 1)(0, 2)) / DX[2];,
621 : ret(1) += (dw(i, j, k)(1, 2) - dw(i, j, k - 1)(1, 2)) / DX[2];,
622 : ret(2) += (dw(i, j, k)(2, 2) - dw(i, j, k - 1)(2, 2)) / DX[2];)
623 : }
624 0 : else if (stencil[2] == StencilType::Hi)
625 : {
626 0 : AMREX_D_TERM(ret(0) += (dw(i, j, k + 1)(0, 2) - dw(i, j, k)(0, 2)) / DX[2];,
627 : ret(1) += (dw(i, j, k + 1)(1, 2) - dw(i, j, k)(1, 2)) / DX[2];,
628 : ret(2) += (dw(i, j, k + 1)(2, 2) - dw(i, j, k)(2, 2)) / DX[2];)
629 : }
630 : #endif
631 160392 : return ret;
632 : }
633 :
634 :
635 : [[nodiscard]] AMREX_FORCE_INLINE
636 : Set::Scalar
637 : Divergence(const amrex::Array4<const Set::Scalar>& f,
638 : const int& i, const int& j, const int& k, const int& m,
639 : const Set::Scalar dx[AMREX_SPACEDIM],
640 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
641 : {
642 : Set::Scalar ret;
643 : ret = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, m, dx, stencil));
644 : #if AMREX_SPACEDIM > 1
645 : ret += (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, m, dx, stencil));
646 : #endif
647 : #if AMREX_SPACEDIM > 2
648 : ret += (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, m, dx, stencil));
649 : #endif
650 : return ret;
651 : }
652 :
653 :
654 : [[nodiscard]] AMREX_FORCE_INLINE
655 : Set::Vector
656 : Gradient(const amrex::Array4<const Set::Scalar>& f,
657 : const int& i, const int& j, const int& k, const int& m,
658 : const Set::Scalar dx[AMREX_SPACEDIM],
659 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
660 : {
661 38416853 : Set::Vector ret;
662 38416853 : ret(0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, m, dx, stencil));
663 : #if AMREX_SPACEDIM > 1
664 38416853 : ret(1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, m, dx, stencil));
665 : #if AMREX_SPACEDIM > 2
666 78784 : ret(2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, m, dx, stencil));
667 : #endif
668 : #endif
669 38416853 : return ret;
670 : }
671 :
672 : [[nodiscard]] AMREX_FORCE_INLINE
673 : Set::Vector
674 : CellGradientOnNode(const amrex::Array4<const Set::Scalar>& f,
675 : const int& i, const int& j, const int& k, const int& m,
676 : const Set::Scalar dx[AMREX_SPACEDIM])
677 : {
678 4358180 : Set::Vector ret;
679 : #if AMREX_SPACEDIM == 1
680 : ret(0) = (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0];
681 : #elif AMREX_SPACEDIM == 2
682 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];
683 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];
684 : #elif AMREX_SPACEDIM == 3
685 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];
686 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];
687 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];
688 : #endif
689 4358180 : return ret;
690 : }
691 :
692 :
693 : template<class T>
694 : [[nodiscard]] AMREX_FORCE_INLINE
695 : std::array<T, AMREX_SPACEDIM>
696 : CellGradientOnNode(const amrex::Array4<const T>& f,
697 : const int& i, const int& j, const int& k, const int& m,
698 : const Set::Scalar dx[AMREX_SPACEDIM])
699 : {
700 : std::array<T, AMREX_SPACEDIM> ret;
701 : Util::Abort(INFO);
702 : #if AMREX_SPACEDIM == 1
703 : ret[0] = (f(i, j, k, m) - f(i - 1, j, k, m)) / dx[0];
704 : #elif AMREX_SPACEDIM == 2
705 : 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];
706 : 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];
707 : #elif AMREX_SPACEDIM == 3
708 : 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];
709 : 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];
710 : 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];
711 : #endif
712 : return ret;
713 : }
714 :
715 :
716 :
717 : [[nodiscard]] AMREX_FORCE_INLINE
718 : Set::Matrix
719 : Gradient(const amrex::Array4<const Set::Scalar>& f,
720 : const int& i, const int& j, const int& k,
721 : const Set::Scalar dx[AMREX_SPACEDIM],
722 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
723 : {
724 0 : Set::Matrix ret;
725 0 : ret(0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
726 : #if AMREX_SPACEDIM > 1
727 0 : ret(0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
728 0 : ret(1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
729 0 : ret(1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
730 : #if AMREX_SPACEDIM > 2
731 0 : ret(0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil));
732 0 : ret(2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
733 0 : ret(1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 1, dx, stencil));
734 0 : ret(2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
735 0 : ret(2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 2, dx, stencil));
736 : #endif
737 : #endif
738 0 : return ret;
739 : }
740 :
741 : [[nodiscard]] AMREX_FORCE_INLINE
742 : Set::Matrix
743 : Gradient(const amrex::Array4<const Set::Vector>& f,
744 : const int& i, const int& j, const int& k,
745 : const Set::Scalar dx[AMREX_SPACEDIM],
746 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
747 : {
748 446379 : Set::Matrix ret;
749 :
750 : #if AMREX_SPACEDIM > 0
751 892758 : ret.col(0) = Numeric::Stencil<Set::Vector, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil);
752 : #endif
753 : #if AMREX_SPACEDIM > 1
754 892758 : ret.col(1) = Numeric::Stencil<Set::Vector, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil);
755 : #endif
756 : #if AMREX_SPACEDIM > 2
757 208800 : ret.col(2) = Numeric::Stencil<Set::Vector, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil);
758 : #endif
759 :
760 446379 : return ret;
761 : }
762 :
763 : [[nodiscard]] AMREX_FORCE_INLINE
764 : std::pair<Set::Vector,Set::Matrix>
765 : GradientSplit( const amrex::Array4<const Set::Vector>& f,
766 : const int& i, const int& j, const int& k,
767 : const Set::Scalar dx[AMREX_SPACEDIM],
768 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
769 : {
770 0 : Set::Vector diag;
771 0 : Set::Matrix offdiag;
772 :
773 0 : std::pair<Set::Scalar,Set::Vector> ret;
774 : #if AMREX_SPACEDIM > 0
775 0 : ret = Numeric::Stencil<Set::Vector, 1, 0, 0>::Dsplit(f, i, j, k, 0, dx, stencil);
776 0 : diag(0) = ret.first;
777 0 : offdiag.col(0) = ret.second;
778 : #endif
779 : #if AMREX_SPACEDIM > 1
780 0 : ret = Numeric::Stencil<Set::Vector, 0, 1, 0>::Dsplit(f, i, j, k, 0, dx, stencil);
781 0 : diag(1) = ret.first;
782 0 : offdiag.col(1) = ret.second;
783 : #endif
784 : #if AMREX_SPACEDIM > 2
785 0 : ret = Numeric::Stencil<Set::Vector, 0, 0, 1>::Dsplit(f, i, j, k, 0, dx, stencil);
786 0 : diag(2) = ret.first;
787 0 : offdiag.col(2) = ret.second;
788 : #endif
789 0 : return std::make_pair(diag,offdiag);
790 : }
791 :
792 :
793 : [[nodiscard]] AMREX_FORCE_INLINE
794 : Set::Matrix
795 : NodeGradientOnCell(const amrex::Array4<const Set::Vector>& f,
796 : const int& i, const int& j, const int& k,
797 : const Set::Scalar dx[AMREX_SPACEDIM])
798 : {
799 : Set::Matrix ret;
800 : #if AMREX_SPACEDIM == 1
801 : ret.col(0) = (f(i + 1, j, k) - f(i, j, k)) / dx[0];
802 : #elif AMREX_SPACEDIM == 2
803 : 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];
804 : 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];
805 : #elif AMREX_SPACEDIM == 3
806 : 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];
807 : 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];
808 : 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];
809 : #endif
810 : return ret;
811 : }
812 :
813 : [[nodiscard]] AMREX_FORCE_INLINE
814 : Set::Vector
815 : NodeGradientOnCell(const amrex::Array4<const Set::Scalar>& f,
816 : const int& i, const int& j, const int& k, const int& m,
817 : const Set::Scalar dx[AMREX_SPACEDIM])
818 : {
819 : Set::Vector ret;
820 : #if AMREX_SPACEDIM == 1
821 : ret(0) = (f(i + 1, j, k, m) - f(i, j, k, m)) / dx[0];
822 : #elif AMREX_SPACEDIM == 2
823 : ret(0) = 0.5 * (f(i + 1, j + 1, k, m) - f(i, j + 1, k, m) + f(i + 1, j, k, m) - f(i, j, k, m)) / dx[0];
824 : ret(1) = 0.5 * (f(i + 1, j + 1, k, m) - f(i + 1, j, k, m) + f(i, j + 1, k, m) - f(i, j, k, m)) / dx[1];
825 : #elif AMREX_SPACEDIM == 3
826 : ret(0) = 0.25 * (f(i + 1, j + 1, k + 1, m) - f(i, j + 1, k + 1, m) + f(i + 1, j, k + 1, m) - f(i, j, k + 1, m) + f(i + 1, j + 1, k, m) - f(i, j + 1, k, m) + f(i + 1, j, k, m) - f(i, j, k, m)) / dx[0];
827 : ret(1) = 0.25 * (f(i + 1, j + 1, k + 1, m) - f(i + 1, j, k + 1, m) + f(i, j + 1, k + 1, m) - f(i, j, k + 1, m) + f(i + 1, j + 1, k, m) - f(i + 1, j, k, m) + f(i, j + 1, k, m) - f(i, j, k, m)) / dx[1];
828 : ret(2) = 0.25 * (f(i + 1, j + 1, k + 1, m) - f(i + 1, j + 1, k, m) + f(i, j + 1, k + 1, m) - f(i, j + 1, k, m) + f(i + 1, j, k + 1, m) - f(i + 1, j, k, m) + f(i, j, k + 1, m) - f(i, j, k, m)) / dx[2];
829 : #endif
830 : return ret;
831 : }
832 :
833 : [[nodiscard]] AMREX_FORCE_INLINE
834 : Set::Matrix3
835 : Gradient(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 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
839 : {
840 : Set::Matrix3 ret;
841 :
842 : #if AMREX_SPACEDIM > 0
843 : ret[0] = Numeric::Stencil<Set::Matrix, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil);
844 : #endif
845 : #if AMREX_SPACEDIM > 1
846 : ret[1] = Numeric::Stencil<Set::Matrix, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil);
847 : #endif
848 : #if AMREX_SPACEDIM > 2
849 : ret[2] = Numeric::Stencil<Set::Matrix, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil);
850 : #endif
851 :
852 : return ret;
853 : }
854 : [[nodiscard]] AMREX_FORCE_INLINE
855 : Set::Matrix3
856 : NodeGradientOnCell(const amrex::Array4<const Set::Matrix>& f,
857 : const int& i, const int& j, const int& k,
858 : const Set::Scalar dx[AMREX_SPACEDIM])
859 : {
860 370688 : Set::Matrix3 ret;
861 :
862 : #if AMREX_SPACEDIM > 0
863 1482748 : ret[0] = (f(i+1,j,k) - f(i,j,k)) / dx[0];
864 : #endif
865 : #if AMREX_SPACEDIM > 1
866 1482748 : ret[1] = (f(i,j+1,k) - f(i,j,k)) / dx[1];
867 : #endif
868 : #if AMREX_SPACEDIM > 2
869 0 : ret[2] = (f(i,j,k+1) - f(i,j,k)) / dx[2];
870 : #endif
871 :
872 370688 : return ret;
873 : }
874 :
875 :
876 : [[nodiscard]] AMREX_FORCE_INLINE
877 : Set::Matrix3
878 : MatrixGradient(const amrex::Array4<const Set::Scalar>& f,
879 : const int& i, const int& j, const int& k,
880 : const Set::Scalar dx[AMREX_SPACEDIM],
881 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
882 : {
883 : Set::Matrix3 ret;
884 : #if AMREX_SPACEDIM == 1
885 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
886 : #elif AMREX_SPACEDIM == 2
887 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
888 : ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
889 : ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
890 : ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
891 : ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
892 : ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
893 : ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
894 : ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
895 : #elif AMREX_SPACEDIM == 3
896 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
897 : ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
898 : ret[0](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
899 : ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
900 : ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 4, dx, stencil));
901 : ret[1](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 5, dx, stencil));
902 : ret[2](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 6, dx, stencil));
903 : ret[2](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 7, dx, stencil));
904 : ret[2](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 8, dx, stencil));
905 : ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
906 : ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
907 : ret[0](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
908 : ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
909 : ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 4, dx, stencil));
910 : ret[1](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 5, dx, stencil));
911 : ret[2](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 6, dx, stencil));
912 : ret[2](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 7, dx, stencil));
913 : ret[2](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 8, dx, stencil));
914 : ret[0](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil));
915 : ret[0](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 1, dx, stencil));
916 : ret[0](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 2, dx, stencil));
917 : ret[1](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 3, dx, stencil));
918 : ret[1](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 4, dx, stencil));
919 : ret[1](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 5, dx, stencil));
920 : ret[2](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 6, dx, stencil));
921 : ret[2](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 7, dx, stencil));
922 : ret[2](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 8, dx, stencil));
923 : #endif
924 : return ret;
925 : }
926 :
927 :
928 : [[nodiscard]] AMREX_FORCE_INLINE
929 : Set::Matrix
930 : Hessian(const amrex::Array4<const Set::Scalar>& f,
931 : const int& i, const int& j, const int& k, const int& m,
932 : const Set::Scalar dx[AMREX_SPACEDIM],
933 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType
934 : )
935 : {
936 2401960 : Set::Matrix ret;
937 2401960 : ret(0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx, stencil));
938 : #if AMREX_SPACEDIM > 1
939 2401960 : ret(1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, m, dx, stencil));
940 2401960 : ret(0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, m, dx, stencil));
941 2401960 : ret(1, 0) = ret(0, 1);
942 : #if AMREX_SPACEDIM > 2
943 0 : ret(2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, m, dx, stencil));
944 0 : ret(1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, m, dx, stencil));
945 0 : ret(2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, m, dx, stencil));
946 0 : ret(2, 1) = ret(1, 2);
947 0 : ret(0, 2) = ret(2, 0);
948 : #endif
949 : #endif
950 2401960 : return ret;
951 : }
952 :
953 : [[nodiscard]] AMREX_FORCE_INLINE
954 : Set::Matrix3
955 : Hessian(const amrex::Array4<const Set::Scalar>& f,
956 : const int& i, const int& j, const int& k,
957 : const Set::Scalar DX[AMREX_SPACEDIM],
958 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
959 : {
960 : Set::Matrix3 ret;
961 : // 1D
962 : #if AMREX_SPACEDIM>0
963 : ret(0, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 0, DX, stencil));
964 : #endif
965 : // 2D
966 : #if AMREX_SPACEDIM>1
967 : ret(0, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 0, DX, stencil));
968 : ret(0, 1, 0) = ret(0, 0, 1);
969 : ret(0, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 0, DX, stencil));
970 : ret(1, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 1, DX, stencil));
971 : ret(1, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 1, DX, stencil));
972 : ret(1, 1, 0) = ret(1, 0, 1);
973 : ret(1, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 1, DX, stencil));
974 : #endif
975 : // 3D
976 : #if AMREX_SPACEDIM>2
977 : ret(0, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 0, DX, stencil));
978 : ret(0, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 0, DX, stencil));
979 : ret(0, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 0, DX, stencil));
980 : ret(0, 2, 0) = ret(0, 0, 2);
981 : ret(0, 2, 1) = ret(0, 1, 2);;
982 : ret(1, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 1, DX, stencil));
983 : ret(1, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 1, DX, stencil));
984 : ret(1, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 1, DX, stencil));
985 : ret(1, 2, 0) = ret(1, 0, 2);
986 : ret(1, 2, 1) = ret(1, 1, 2);;
987 : ret(2, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 2, DX, stencil));
988 : ret(2, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 2, DX, stencil));
989 : ret(2, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 2, DX, stencil));
990 : ret(2, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 2, DX, stencil));
991 : ret(2, 1, 0) = ret(2, 0, 1);
992 : ret(2, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 2, DX, stencil));
993 : ret(2, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 2, DX, stencil));
994 : ret(2, 2, 0) = ret(2, 0, 2);
995 : ret(2, 2, 1) = ret(2, 1, 2);;
996 : #endif
997 : return ret;
998 : }
999 :
1000 :
1001 : // Returns Hessian of a vector field.
1002 : // Return value: ret[i](j,k) = ret_{i,jk}
1003 : [[nodiscard]] AMREX_FORCE_INLINE
1004 : Set::Matrix3
1005 : Hessian(const amrex::Array4<const Set::Vector>& f,
1006 : const int& i, const int& j, const int& k,
1007 : const Set::Scalar dx[AMREX_SPACEDIM],
1008 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1009 : {
1010 : Set::Matrix3 ret;
1011 :
1012 : #if AMREX_SPACEDIM>0
1013 : Set::Vector f_11 = Numeric::Stencil<Set::Vector, 2, 0, 0>::D(f, i, j, k, 0, dx, stencil);
1014 : ret[0](0, 0) = f_11(0);
1015 : #endif
1016 : #if AMREX_SPACEDIM>1
1017 : ret[1](0, 0) = f_11(1);
1018 : Set::Vector f_12 = Numeric::Stencil<Set::Vector, 1, 1, 0>::D(f, i, j, k, 0, dx, stencil);
1019 : ret[0](1, 0) = ret[0](0, 1) = f_12(0);
1020 : ret[1](1, 0) = ret[1](0, 1) = f_12(1);
1021 : Set::Vector f_22 = Numeric::Stencil<Set::Vector, 0, 2, 0>::D(f, i, j, k, 0, dx, stencil);
1022 : ret[0](1, 1) = f_22(0);
1023 : ret[1](1, 1) = f_22(1);
1024 : #endif
1025 : #if AMREX_SPACEDIM>2
1026 : ret[2](0, 0) = f_11(2);
1027 : ret[2](1, 0) = ret[2](0, 1) = f_12(2);
1028 : Set::Vector f_13 = Numeric::Stencil<Set::Vector, 1, 0, 1>::D(f, i, j, k, 0, dx, stencil);
1029 : ret[0](2, 0) = ret[0](0, 2) = f_13(0);
1030 : ret[1](2, 0) = ret[1](0, 2) = f_13(1);
1031 : ret[2](2, 0) = ret[2](0, 2) = f_13(2);
1032 : ret[2](1, 1) = f_22(2);
1033 : Set::Vector f_23 = Numeric::Stencil<Set::Vector, 0, 1, 1>::D(f, i, j, k, 0, dx, stencil);
1034 : ret[0](1, 2) = ret[0](2, 1) = f_23(0);
1035 : ret[1](1, 2) = ret[1](2, 1) = f_23(1);
1036 : ret[2](1, 2) = ret[2](2, 1) = f_23(2);
1037 : Set::Vector f_33 = Numeric::Stencil<Set::Vector, 0, 0, 2>::D(f, i, j, k, 0, dx, stencil);
1038 : ret[0](2, 2) = f_33(0);
1039 : ret[1](2, 2) = f_33(1);
1040 : ret[2](2, 2) = f_33(2);
1041 : #endif
1042 : return ret;
1043 : }
1044 :
1045 :
1046 : [[nodiscard]] AMREX_FORCE_INLINE
1047 : Set::Matrix
1048 : FieldToMatrix(const amrex::Array4<const Set::Scalar>& f,
1049 : const int& i, const int& j, const int& k)
1050 : {
1051 0 : Set::Matrix ret;
1052 : #if AMREX_SPACEDIM == 1
1053 : ret(0, 0) = f(i, j, k, 0);
1054 :
1055 : #elif AMREX_SPACEDIM == 2
1056 0 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1057 0 : ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1058 :
1059 : #elif AMREX_SPACEDIM == 3
1060 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);
1061 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);
1062 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);
1063 : #endif
1064 :
1065 0 : return ret;
1066 : }
1067 :
1068 : [[nodiscard]] AMREX_FORCE_INLINE
1069 : Set::Matrix
1070 : FieldToMatrix(const amrex::Array4<Set::Scalar>& f,
1071 : const int& i, const int& j, const int& k)
1072 : {
1073 : Set::Matrix ret;
1074 : #if AMREX_SPACEDIM == 1
1075 : ret(0, 0) = f(i, j, k, 0);
1076 :
1077 : #elif AMREX_SPACEDIM == 2
1078 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1079 : ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1080 :
1081 : #elif AMREX_SPACEDIM == 3
1082 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1); ret(0, 2) = f(i, j, k, 2);
1083 : ret(1, 0) = f(i, j, k, 3); ret(1, 1) = f(i, j, k, 4); ret(1, 2) = f(i, j, k, 5);
1084 : ret(2, 0) = f(i, j, k, 6); ret(2, 1) = f(i, j, k, 7); ret(2, 2) = f(i, j, k, 8);
1085 : #endif
1086 :
1087 : return ret;
1088 : }
1089 :
1090 : [[nodiscard]] AMREX_FORCE_INLINE
1091 : Set::Vector
1092 : FieldToVector(const amrex::Array4<const Set::Scalar>& f,
1093 : const int& i, const int& j, const int& k)
1094 : {
1095 : Set::Vector ret;
1096 : ret(0) = f(i, j, k, 0);
1097 : #if AMREX_SPACEDIM > 1
1098 : ret(1) = f(i, j, k, 1);
1099 : #if AMREX_SPACEDIM > 2
1100 : ret(2) = f(i, j, k, 2);
1101 : #endif
1102 : #endif
1103 : return ret;
1104 : }
1105 :
1106 : [[nodiscard]] AMREX_FORCE_INLINE
1107 : Set::Vector
1108 : FieldToVector(const amrex::Array4<Set::Scalar>& f,
1109 : const int& i, const int& j, const int& k)
1110 : {
1111 : Set::Vector ret;
1112 : ret(0) = f(i, j, k, 0);
1113 : #if AMREX_SPACEDIM > 1
1114 : ret(1) = f(i, j, k, 1);
1115 : #if AMREX_SPACEDIM > 2
1116 : ret(2) = f(i, j, k, 2);
1117 : #endif
1118 : #endif
1119 : return ret;
1120 : }
1121 :
1122 : AMREX_FORCE_INLINE
1123 : void
1124 : MatrixToField(const amrex::Array4<Set::Scalar>& f,
1125 : const int& i, const int& j, const int& k,
1126 : Set::Matrix matrix)
1127 : {
1128 : #if AMREX_SPACEDIM == 1
1129 : f(i, j, k, 0) = matrix(0, 0);
1130 : #elif AMREX_SPACEDIM == 2
1131 : f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1);
1132 : f(i, j, k, 2) = matrix(1, 0); f(i, j, k, 3) = matrix(1, 1);
1133 : #elif AMREX_SPACEDIM == 3
1134 : f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1); f(i, j, k, 2) = matrix(0, 2);
1135 : f(i, j, k, 3) = matrix(1, 0); f(i, j, k, 4) = matrix(1, 1); f(i, j, k, 5) = matrix(1, 2);
1136 : f(i, j, k, 6) = matrix(2, 0); f(i, j, k, 7) = matrix(2, 1); f(i, j, k, 8) = matrix(2, 2);
1137 : #endif
1138 : }
1139 :
1140 : AMREX_FORCE_INLINE
1141 : void
1142 : VectorToField(const amrex::Array4<Set::Scalar>& f,
1143 : const int& i, const int& j, const int& k,
1144 : Set::Vector vector)
1145 : {
1146 : f(i, j, k, 0) = vector(0);
1147 : #if AMREX_SPACEDIM > 1
1148 : f(i, j, k, 1) = vector(1);
1149 : #if AMREX_SPACEDIM > 2
1150 : f(i, j, k, 2) = vector(2);
1151 : #endif
1152 : #endif
1153 : }
1154 :
1155 : template<int index, int SYM>
1156 : [[nodiscard]]
1157 : Set::Matrix3
1158 : Divergence(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, SYM>>&,
1159 : const int, const int, const int, const Set::Scalar[AMREX_SPACEDIM],
1160 : std::array<StencilType, AMREX_SPACEDIM> /*stecil*/ = DefaultType)
1161 : {
1162 : Util::Abort(INFO, "Not implemented yet"); return Set::Matrix3::Zero();
1163 : }
1164 :
1165 : template<>
1166 : [[nodiscard]] AMREX_FORCE_INLINE
1167 : Set::Matrix3 Divergence<2, Set::Sym::Isotropic>(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>>& C,
1168 : const int i, const int j, const int k, const Set::Scalar dx[AMREX_SPACEDIM],
1169 : std::array<StencilType, AMREX_SPACEDIM> stencil)
1170 : {
1171 : Set::Matrix3 ret;
1172 :
1173 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCx =
1174 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 1, 0, 0>::D(C, i, j, k, 0, dx, stencil);
1175 : #if AMREX_SPACEDIM>1
1176 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCy =
1177 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 1, 0>::D(C, i, j, k, 0, dx, stencil);
1178 : #endif
1179 : #if AMREX_SPACEDIM>2
1180 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCz =
1181 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 0, 1>::D(C, i, j, k, 0, dx, stencil);
1182 : #endif
1183 :
1184 : for (int i = 0; i < AMREX_SPACEDIM; i++)
1185 : for (int k = 0; k < AMREX_SPACEDIM; k++)
1186 : for (int l = 0; l < AMREX_SPACEDIM; l++)
1187 : {
1188 : ret(i, k, l) = 0.0;
1189 : ret(i, k, l) = gradCx(i, 0, k, l);
1190 : #if AMREX_SPACEDIM>1
1191 : ret(i, k, l) = gradCy(i, 1, k, l);
1192 : #endif
1193 : #if AMREX_SPACEDIM>2
1194 : ret(i, k, l) = gradCz(i, 2, k, l);
1195 : #endif
1196 : }
1197 : return ret;
1198 : }
1199 :
1200 :
1201 :
1202 :
1203 : template<int dim>
1204 : [[nodiscard]] AMREX_FORCE_INLINE
1205 : Set::Matrix4<dim, Set::Sym::Full>
1206 : DoubleHessian(const amrex::Array4<const Set::Scalar>& f,
1207 : const int& i, const int& j, const int& k, const int& m,
1208 : const Set::Scalar dx[AMREX_SPACEDIM]);
1209 :
1210 : template<>
1211 : [[nodiscard]] AMREX_FORCE_INLINE
1212 : Set::Matrix4<2, Set::Sym::Full>
1213 : DoubleHessian<2>(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 1362860 : Set::Matrix4<2, Set::Sym::Full> ret;
1218 : // [0,0,0,0]
1219 2725720 : ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1220 : // [0, 0, 0, 1]
1221 2725720 : ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1222 : // [0, 0, 1, 1]
1223 2725720 : ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1224 : // [0, 1, 1, 1]
1225 2725720 : ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1226 : // [1, 1, 1, 1]
1227 1362860 : ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1228 1362860 : return ret;
1229 : }
1230 :
1231 : template<>
1232 : [[nodiscard]] AMREX_FORCE_INLINE
1233 : Set::Matrix4<3, Set::Sym::Full>
1234 : DoubleHessian<3>(const amrex::Array4<const Set::Scalar>& f,
1235 : const int& i, const int& j, const int& k, const int& m,
1236 : const Set::Scalar dx[AMREX_SPACEDIM])
1237 : {
1238 0 : Set::Matrix4<3, Set::Sym::Full> ret;
1239 : // [0,0,0,0]
1240 0 : ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1241 : // [0, 0, 0, 1]
1242 0 : ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1243 : // [0, 0, 0, 2]
1244 0 : ret(0, 0, 0, 2) = Stencil<Set::Scalar, 3, 0, 1>::D(f, i, j, k, m, dx);
1245 : // [0, 0, 1, 1]
1246 0 : ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1247 : // [0, 0, 1, 2]
1248 0 : ret(0, 0, 1, 2) = Stencil<Set::Scalar, 2, 1, 1>::D(f, i, j, k, m, dx);
1249 : // [0, 0, 2, 2]
1250 0 : ret(0, 0, 2, 2) = Stencil<Set::Scalar, 2, 0, 2>::D(f, i, j, k, m, dx);
1251 : // [0, 1, 1, 1]
1252 0 : ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1253 : // [0, 1, 1, 2]
1254 0 : ret(0, 1, 1, 2) = Stencil<Set::Scalar, 1, 2, 1>::D(f, i, j, k, m, dx);
1255 : // [0, 1, 2, 2]
1256 0 : ret(0, 1, 2, 2) = Stencil<Set::Scalar, 1, 1, 2>::D(f, i, j, k, m, dx);
1257 : // [0, 2, 2, 2]
1258 0 : ret(0, 2, 2, 2) = Stencil<Set::Scalar, 1, 0, 3>::D(f, i, j, k, m, dx);
1259 : // [1, 1, 1, 1]
1260 0 : ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1261 : // [1, 1, 1, 2]
1262 0 : ret(1, 1, 1, 2) = Stencil<Set::Scalar, 0, 3, 1>::D(f, i, j, k, m, dx);
1263 : // [1, 1, 2, 2]
1264 0 : ret(1, 1, 2, 2) = Stencil<Set::Scalar, 0, 2, 2>::D(f, i, j, k, m, dx);
1265 : // [1, 2, 2, 2]
1266 0 : ret(1, 2, 2, 2) = Stencil<Set::Scalar, 0, 1, 3>::D(f, i, j, k, m, dx);
1267 : // [2, 2, 2, 2]
1268 0 : ret(2, 2, 2, 2) = Stencil<Set::Scalar, 0, 0, 4>::D(f, i, j, k, m, dx);
1269 0 : return ret;
1270 : }
1271 :
1272 : struct Interpolate
1273 : {
1274 : public:
1275 : template<class T>
1276 : [[nodiscard]] AMREX_FORCE_INLINE
1277 : static T CellToNodeAverage(const amrex::Array4<const T>& f,
1278 : const int& i, const int& j, const int& k, const int& m,
1279 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1280 : {
1281 183397 : int AMREX_D_DECL(
1282 : ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1283 : jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1284 : klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1285 :
1286 733588 : return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1287 : ,
1288 : +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1289 : ,
1290 : +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1291 : + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1292 183397 : )) * fac;
1293 : }
1294 : template<class T>
1295 : [[nodiscard]] AMREX_FORCE_INLINE
1296 : static T CellToNodeAverage(const amrex::Array4<T>& f,
1297 : const int& i, const int& j, const int& k, const int& m,
1298 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1299 : {
1300 4620090 : int AMREX_D_DECL(
1301 : ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1302 : jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1303 : klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1304 :
1305 18480400 : return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1306 : ,
1307 : +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1308 : ,
1309 : +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1310 : + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1311 4620090 : )) * fac;
1312 : }
1313 : template<class T>
1314 : [[nodiscard]] AMREX_FORCE_INLINE
1315 : static T NodeToCellAverage(const amrex::Array4<const T>& f,
1316 : const int& i, const int& j, const int& k, const int& m)
1317 : {
1318 21086196 : return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1319 : ,
1320 : +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1321 : ,
1322 : +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1323 : + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1324 5271544 : )) * fac;
1325 : }
1326 : template<class T>
1327 : [[nodiscard]] AMREX_FORCE_INLINE
1328 : static T NodeToCellAverage(const amrex::Array4<T>& f,
1329 : const int& i, const int& j, const int& k, const int& m)
1330 : {
1331 : return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1332 : ,
1333 : +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1334 : ,
1335 : +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1336 : + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1337 : )) * fac;
1338 : }
1339 : constexpr static Set::Scalar fac = AMREX_D_PICK(0.5, 0.25, 0.125);
1340 : };
1341 :
1342 : }
1343 : #endif
|
