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