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