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 34896490 : 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 34896490 : 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 138920567 : 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 128698198 : 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 392530511 : 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 138920567 : if (stencil[1] == StencilType::Lo)
112 29696247 : return (f(i, j, k, m) - f(i, j - 1, k, m)) / dx[1];
113 129021818 : 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 394472231 : 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 592639558 : if (stencil[0] == StencilType::Central)
198 2370278904 : 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 139664 : 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 592676584 : if (stencil[1] == StencilType::Central)
214 2358120608 : 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 6292864 : 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 56164970 : int ihi = 1, ilo = 1, jhi = 1, jlo = 1;
246 56164970 : Set::Scalar ifac = 0.5, jfac = 0.5;
247 56164970 : if (stencil[0] == StencilType::Hi) { ilo = 0; ifac = 1.0; }
248 56164970 : if (stencil[0] == StencilType::Lo) { ihi = 0; ifac = 1.0; }
249 56164970 : if (stencil[1] == StencilType::Hi) { jlo = 0; jfac = 1.0; }
250 56164970 : if (stencil[1] == StencilType::Lo) { jhi = 0; jfac = 1.0; }
251 :
252 280824850 : 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 536474588 : Set::Scalar ret = 0.0;
552 536474588 : ret += (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx, stencil));
553 : #if AMREX_SPACEDIM > 1
554 536474588 : 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 536474588 : 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 50960883 : Set::Vector ret;
678 50960883 : ret(0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, m, dx, stencil));
679 : #if AMREX_SPACEDIM > 1
680 50960883 : 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 50960883 : 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 6272000 : Set::Matrix ret;
741 12544000 : ret(0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
742 : #if AMREX_SPACEDIM > 1
743 12544000 : ret(0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
744 12544000 : ret(1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
745 12544000 : 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 6272000 : 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::Matrix
811 : NodeGradientOnCell(const amrex::Array4<const Set::Vector>& f,
812 : const int& i, const int& j, const int& k,
813 : const Set::Scalar dx[AMREX_SPACEDIM])
814 : {
815 0 : Set::Matrix ret;
816 : #if AMREX_SPACEDIM == 1
817 : ret.col(0) = (f(i + 1, j, k) - f(i, j, k)) / dx[0];
818 : #elif AMREX_SPACEDIM == 2
819 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];
820 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];
821 : #elif AMREX_SPACEDIM == 3
822 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];
823 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];
824 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];
825 : #endif
826 0 : return ret;
827 : }
828 :
829 : [[nodiscard]] AMREX_FORCE_INLINE
830 : Set::Vector
831 : NodeGradientOnCell(const amrex::Array4<const Set::Scalar>& f,
832 : const int& i, const int& j, const int& k, const int& m,
833 : const Set::Scalar dx[AMREX_SPACEDIM])
834 : {
835 : Set::Vector ret;
836 : #if AMREX_SPACEDIM == 1
837 : ret(0) = (f(i + 1, j, k, m) - f(i, j, k, m)) / dx[0];
838 : #elif AMREX_SPACEDIM == 2
839 : 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];
840 : 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];
841 : #elif AMREX_SPACEDIM == 3
842 : 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];
843 : 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];
844 : 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];
845 : #endif
846 : return ret;
847 : }
848 :
849 : [[nodiscard]] AMREX_FORCE_INLINE
850 : Set::Matrix3
851 : Gradient(const amrex::Array4<const Set::Matrix>& f,
852 : const int& i, const int& j, const int& k,
853 : const Set::Scalar dx[AMREX_SPACEDIM],
854 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
855 : {
856 : Set::Matrix3 ret;
857 :
858 : #if AMREX_SPACEDIM > 0
859 : ret[0] = Numeric::Stencil<Set::Matrix, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil);
860 : #endif
861 : #if AMREX_SPACEDIM > 1
862 : ret[1] = Numeric::Stencil<Set::Matrix, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil);
863 : #endif
864 : #if AMREX_SPACEDIM > 2
865 : ret[2] = Numeric::Stencil<Set::Matrix, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil);
866 : #endif
867 :
868 : return ret;
869 : }
870 : [[nodiscard]] AMREX_FORCE_INLINE
871 : Set::Matrix3
872 : NodeGradientOnCell(const amrex::Array4<const Set::Matrix>& f,
873 : const int& i, const int& j, const int& k,
874 : const Set::Scalar dx[AMREX_SPACEDIM])
875 : {
876 370688 : Set::Matrix3 ret;
877 :
878 : #if AMREX_SPACEDIM > 0
879 1482752 : ret[0] = (f(i+1,j,k) - f(i,j,k)) / dx[0];
880 : #endif
881 : #if AMREX_SPACEDIM > 1
882 1482752 : ret[1] = (f(i,j+1,k) - f(i,j,k)) / dx[1];
883 : #endif
884 : #if AMREX_SPACEDIM > 2
885 0 : ret[2] = (f(i,j,k+1) - f(i,j,k)) / dx[2];
886 : #endif
887 :
888 370688 : return ret;
889 : }
890 :
891 :
892 : [[nodiscard]] AMREX_FORCE_INLINE
893 : Set::Matrix3
894 : MatrixGradient(const amrex::Array4<const Set::Scalar>& f,
895 : const int& i, const int& j, const int& k,
896 : const Set::Scalar dx[AMREX_SPACEDIM],
897 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
898 : {
899 : Set::Matrix3 ret;
900 : #if AMREX_SPACEDIM == 1
901 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
902 : #elif AMREX_SPACEDIM == 2
903 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
904 : ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
905 : ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
906 : ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
907 : ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
908 : ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
909 : ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
910 : ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
911 : #elif AMREX_SPACEDIM == 3
912 : ret[0](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 0, dx, stencil));
913 : ret[0](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 1, dx, stencil));
914 : ret[0](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 2, dx, stencil));
915 : ret[1](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 3, dx, stencil));
916 : ret[1](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 4, dx, stencil));
917 : ret[1](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 5, dx, stencil));
918 : ret[2](0, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 6, dx, stencil));
919 : ret[2](1, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 7, dx, stencil));
920 : ret[2](2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(f, i, j, k, 8, dx, stencil));
921 : ret[0](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 0, dx, stencil));
922 : ret[0](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 1, dx, stencil));
923 : ret[0](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 2, dx, stencil));
924 : ret[1](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 3, dx, stencil));
925 : ret[1](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 4, dx, stencil));
926 : ret[1](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 5, dx, stencil));
927 : ret[2](0, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 6, dx, stencil));
928 : ret[2](1, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 7, dx, stencil));
929 : ret[2](2, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(f, i, j, k, 8, dx, stencil));
930 : ret[0](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 0, dx, stencil));
931 : ret[0](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 1, dx, stencil));
932 : ret[0](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 2, dx, stencil));
933 : ret[1](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 3, dx, stencil));
934 : ret[1](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 4, dx, stencil));
935 : ret[1](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 5, dx, stencil));
936 : ret[2](0, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 6, dx, stencil));
937 : ret[2](1, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 7, dx, stencil));
938 : ret[2](2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(f, i, j, k, 8, dx, stencil));
939 : #endif
940 : return ret;
941 : }
942 :
943 :
944 : [[nodiscard]] AMREX_FORCE_INLINE
945 : Set::Matrix
946 : Hessian(const amrex::Array4<const Set::Scalar>& f,
947 : const int& i, const int& j, const int& k, const int& m,
948 : const Set::Scalar dx[AMREX_SPACEDIM],
949 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType
950 : )
951 : {
952 14945721 : Set::Matrix ret;
953 14945721 : ret(0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, m, dx, stencil));
954 : #if AMREX_SPACEDIM > 1
955 14945721 : ret(1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, m, dx, stencil));
956 14945721 : ret(0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, m, dx, stencil));
957 14945721 : ret(1, 0) = ret(0, 1);
958 : #if AMREX_SPACEDIM > 2
959 0 : ret(2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, m, dx, stencil));
960 0 : ret(1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, m, dx, stencil));
961 0 : ret(2, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, m, dx, stencil));
962 0 : ret(2, 1) = ret(1, 2);
963 0 : ret(0, 2) = ret(2, 0);
964 : #endif
965 : #endif
966 14945721 : return ret;
967 : }
968 :
969 : [[nodiscard]] AMREX_FORCE_INLINE
970 : Set::Matrix3
971 : Hessian(const amrex::Array4<const Set::Scalar>& f,
972 : const int& i, const int& j, const int& k,
973 : const Set::Scalar DX[AMREX_SPACEDIM],
974 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
975 : {
976 6272000 : Set::Matrix3 ret;
977 : // 1D
978 : #if AMREX_SPACEDIM>0
979 12544000 : ret(0, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 0, DX, stencil));
980 : #endif
981 : // 2D
982 : #if AMREX_SPACEDIM>1
983 18816000 : ret(0, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 0, DX, stencil));
984 6272000 : ret(0, 1, 0) = ret(0, 0, 1);
985 12544000 : ret(0, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 0, DX, stencil));
986 12544000 : ret(1, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 1, DX, stencil));
987 18816000 : ret(1, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 1, DX, stencil));
988 6272000 : ret(1, 1, 0) = ret(1, 0, 1);
989 12544000 : ret(1, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 1, DX, stencil));
990 : #endif
991 : // 3D
992 : #if AMREX_SPACEDIM>2
993 : ret(0, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 0, DX, stencil));
994 : ret(0, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 0, DX, stencil));
995 : ret(0, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 0, DX, stencil));
996 : ret(0, 2, 0) = ret(0, 0, 2);
997 : ret(0, 2, 1) = ret(0, 1, 2);;
998 : ret(1, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 1, DX, stencil));
999 : ret(1, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 1, DX, stencil));
1000 : ret(1, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 1, DX, stencil));
1001 : ret(1, 2, 0) = ret(1, 0, 2);
1002 : ret(1, 2, 1) = ret(1, 1, 2);;
1003 : ret(2, 0, 0) = (Numeric::Stencil<Set::Scalar, 2, 0, 0>::D(f, i, j, k, 2, DX, stencil));
1004 : ret(2, 1, 1) = (Numeric::Stencil<Set::Scalar, 0, 2, 0>::D(f, i, j, k, 2, DX, stencil));
1005 : ret(2, 2, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 2>::D(f, i, j, k, 2, DX, stencil));
1006 : ret(2, 0, 1) = (Numeric::Stencil<Set::Scalar, 1, 1, 0>::D(f, i, j, k, 2, DX, stencil));
1007 : ret(2, 1, 0) = ret(2, 0, 1);
1008 : ret(2, 0, 2) = (Numeric::Stencil<Set::Scalar, 1, 0, 1>::D(f, i, j, k, 2, DX, stencil));
1009 : ret(2, 1, 2) = (Numeric::Stencil<Set::Scalar, 0, 1, 1>::D(f, i, j, k, 2, DX, stencil));
1010 : ret(2, 2, 0) = ret(2, 0, 2);
1011 : ret(2, 2, 1) = ret(2, 1, 2);;
1012 : #endif
1013 6272000 : return ret;
1014 : }
1015 :
1016 :
1017 : // Returns Hessian of a vector field.
1018 : // Return value: ret[i](j,k) = ret_{i,jk}
1019 : [[nodiscard]] AMREX_FORCE_INLINE
1020 : Set::Matrix3
1021 : Hessian(const amrex::Array4<const Set::Vector>& f,
1022 : const int& i, const int& j, const int& k,
1023 : const Set::Scalar dx[AMREX_SPACEDIM],
1024 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1025 : {
1026 : Set::Matrix3 ret;
1027 :
1028 : #if AMREX_SPACEDIM>0
1029 : Set::Vector f_11 = Numeric::Stencil<Set::Vector, 2, 0, 0>::D(f, i, j, k, 0, dx, stencil);
1030 : ret[0](0, 0) = f_11(0);
1031 : #endif
1032 : #if AMREX_SPACEDIM>1
1033 : ret[1](0, 0) = f_11(1);
1034 : Set::Vector f_12 = Numeric::Stencil<Set::Vector, 1, 1, 0>::D(f, i, j, k, 0, dx, stencil);
1035 : ret[0](1, 0) = ret[0](0, 1) = f_12(0);
1036 : ret[1](1, 0) = ret[1](0, 1) = f_12(1);
1037 : Set::Vector f_22 = Numeric::Stencil<Set::Vector, 0, 2, 0>::D(f, i, j, k, 0, dx, stencil);
1038 : ret[0](1, 1) = f_22(0);
1039 : ret[1](1, 1) = f_22(1);
1040 : #endif
1041 : #if AMREX_SPACEDIM>2
1042 : ret[2](0, 0) = f_11(2);
1043 : ret[2](1, 0) = ret[2](0, 1) = f_12(2);
1044 : Set::Vector f_13 = Numeric::Stencil<Set::Vector, 1, 0, 1>::D(f, i, j, k, 0, dx, stencil);
1045 : ret[0](2, 0) = ret[0](0, 2) = f_13(0);
1046 : ret[1](2, 0) = ret[1](0, 2) = f_13(1);
1047 : ret[2](2, 0) = ret[2](0, 2) = f_13(2);
1048 : ret[2](1, 1) = f_22(2);
1049 : Set::Vector f_23 = Numeric::Stencil<Set::Vector, 0, 1, 1>::D(f, i, j, k, 0, dx, stencil);
1050 : ret[0](1, 2) = ret[0](2, 1) = f_23(0);
1051 : ret[1](1, 2) = ret[1](2, 1) = f_23(1);
1052 : ret[2](1, 2) = ret[2](2, 1) = f_23(2);
1053 : Set::Vector f_33 = Numeric::Stencil<Set::Vector, 0, 0, 2>::D(f, i, j, k, 0, dx, stencil);
1054 : ret[0](2, 2) = f_33(0);
1055 : ret[1](2, 2) = f_33(1);
1056 : ret[2](2, 2) = f_33(2);
1057 : #endif
1058 : return ret;
1059 : }
1060 :
1061 :
1062 : [[nodiscard]] AMREX_FORCE_INLINE
1063 : Set::Matrix
1064 : FieldToMatrix(const amrex::Array4<const Set::Scalar>& f,
1065 : const int& i, const int& j, const int& k)
1066 : {
1067 : Set::Matrix ret;
1068 : #if AMREX_SPACEDIM == 1
1069 : ret(0, 0) = f(i, j, k, 0);
1070 :
1071 : #elif AMREX_SPACEDIM == 2
1072 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1073 : ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1074 :
1075 : #elif AMREX_SPACEDIM == 3
1076 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1); ret(0, 2) = f(i, j, k, 2);
1077 : ret(1, 0) = f(i, j, k, 3); ret(1, 1) = f(i, j, k, 4); ret(1, 2) = f(i, j, k, 5);
1078 : ret(2, 0) = f(i, j, k, 6); ret(2, 1) = f(i, j, k, 7); ret(2, 2) = f(i, j, k, 8);
1079 : #endif
1080 :
1081 : return ret;
1082 : }
1083 :
1084 : [[nodiscard]] AMREX_FORCE_INLINE
1085 : Set::Matrix
1086 : FieldToMatrix(const amrex::Array4<Set::Scalar>& f,
1087 : const int& i, const int& j, const int& k)
1088 : {
1089 : Set::Matrix ret;
1090 : #if AMREX_SPACEDIM == 1
1091 : ret(0, 0) = f(i, j, k, 0);
1092 :
1093 : #elif AMREX_SPACEDIM == 2
1094 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1);
1095 : ret(1, 0) = f(i, j, k, 2); ret(1, 1) = f(i, j, k, 3);
1096 :
1097 : #elif AMREX_SPACEDIM == 3
1098 : ret(0, 0) = f(i, j, k, 0); ret(0, 1) = f(i, j, k, 1); ret(0, 2) = f(i, j, k, 2);
1099 : ret(1, 0) = f(i, j, k, 3); ret(1, 1) = f(i, j, k, 4); ret(1, 2) = f(i, j, k, 5);
1100 : ret(2, 0) = f(i, j, k, 6); ret(2, 1) = f(i, j, k, 7); ret(2, 2) = f(i, j, k, 8);
1101 : #endif
1102 :
1103 : return ret;
1104 : }
1105 :
1106 : [[nodiscard]] AMREX_FORCE_INLINE
1107 : Set::Vector
1108 : FieldToVector(const amrex::Array4<const Set::Scalar>& f,
1109 : const int& i, const int& j, const int& k)
1110 : {
1111 : Set::Vector ret;
1112 : ret(0) = f(i, j, k, 0);
1113 : #if AMREX_SPACEDIM > 1
1114 : ret(1) = f(i, j, k, 1);
1115 : #if AMREX_SPACEDIM > 2
1116 : ret(2) = f(i, j, k, 2);
1117 : #endif
1118 : #endif
1119 : return ret;
1120 : }
1121 :
1122 : [[nodiscard]] AMREX_FORCE_INLINE
1123 : Set::Vector
1124 : FieldToVector(const amrex::Array4<Set::Scalar>& f,
1125 : const int& i, const int& j, const int& k)
1126 : {
1127 : Set::Vector ret;
1128 : ret(0) = f(i, j, k, 0);
1129 : #if AMREX_SPACEDIM > 1
1130 : ret(1) = f(i, j, k, 1);
1131 : #if AMREX_SPACEDIM > 2
1132 : ret(2) = f(i, j, k, 2);
1133 : #endif
1134 : #endif
1135 : return ret;
1136 : }
1137 :
1138 : AMREX_FORCE_INLINE
1139 : void
1140 : MatrixToField(const amrex::Array4<Set::Scalar>& f,
1141 : const int& i, const int& j, const int& k,
1142 : Set::Matrix matrix)
1143 : {
1144 : #if AMREX_SPACEDIM == 1
1145 : f(i, j, k, 0) = matrix(0, 0);
1146 : #elif AMREX_SPACEDIM == 2
1147 : f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1);
1148 : f(i, j, k, 2) = matrix(1, 0); f(i, j, k, 3) = matrix(1, 1);
1149 : #elif AMREX_SPACEDIM == 3
1150 : f(i, j, k, 0) = matrix(0, 0); f(i, j, k, 1) = matrix(0, 1); f(i, j, k, 2) = matrix(0, 2);
1151 : f(i, j, k, 3) = matrix(1, 0); f(i, j, k, 4) = matrix(1, 1); f(i, j, k, 5) = matrix(1, 2);
1152 : f(i, j, k, 6) = matrix(2, 0); f(i, j, k, 7) = matrix(2, 1); f(i, j, k, 8) = matrix(2, 2);
1153 : #endif
1154 : }
1155 :
1156 : AMREX_FORCE_INLINE
1157 : void
1158 : VectorToField(const amrex::Array4<Set::Scalar>& f,
1159 : const int& i, const int& j, const int& k,
1160 : Set::Vector vector)
1161 : {
1162 : f(i, j, k, 0) = vector(0);
1163 : #if AMREX_SPACEDIM > 1
1164 : f(i, j, k, 1) = vector(1);
1165 : #if AMREX_SPACEDIM > 2
1166 : f(i, j, k, 2) = vector(2);
1167 : #endif
1168 : #endif
1169 : }
1170 :
1171 : template<int index, int SYM>
1172 : [[nodiscard]]
1173 : Set::Matrix3
1174 : Divergence(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, SYM>>&,
1175 : const int, const int, const int, const Set::Scalar[AMREX_SPACEDIM],
1176 : std::array<StencilType, AMREX_SPACEDIM> /*stecil*/ = DefaultType)
1177 : {
1178 : Util::Abort(INFO, "Not implemented yet"); return Set::Matrix3::Zero();
1179 : }
1180 :
1181 : template<>
1182 : [[nodiscard]] AMREX_FORCE_INLINE
1183 : Set::Matrix3 Divergence<2, Set::Sym::Isotropic>(const amrex::Array4<const Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>>& C,
1184 : const int i, const int j, const int k, const Set::Scalar dx[AMREX_SPACEDIM],
1185 : std::array<StencilType, AMREX_SPACEDIM> stencil)
1186 : {
1187 : Set::Matrix3 ret;
1188 :
1189 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCx =
1190 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 1, 0, 0>::D(C, i, j, k, 0, dx, stencil);
1191 : #if AMREX_SPACEDIM>1
1192 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCy =
1193 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 1, 0>::D(C, i, j, k, 0, dx, stencil);
1194 : #endif
1195 : #if AMREX_SPACEDIM>2
1196 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic> gradCz =
1197 : Numeric::Stencil<Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Isotropic>, 0, 0, 1>::D(C, i, j, k, 0, dx, stencil);
1198 : #endif
1199 :
1200 : for (int i = 0; i < AMREX_SPACEDIM; i++)
1201 : for (int k = 0; k < AMREX_SPACEDIM; k++)
1202 : for (int l = 0; l < AMREX_SPACEDIM; l++)
1203 : {
1204 : ret(i, k, l) = 0.0;
1205 : ret(i, k, l) = gradCx(i, 0, k, l);
1206 : #if AMREX_SPACEDIM>1
1207 : ret(i, k, l) = gradCy(i, 1, k, l);
1208 : #endif
1209 : #if AMREX_SPACEDIM>2
1210 : ret(i, k, l) = gradCz(i, 2, k, l);
1211 : #endif
1212 : }
1213 : return ret;
1214 : }
1215 :
1216 :
1217 :
1218 :
1219 : template<int dim>
1220 : [[nodiscard]] AMREX_FORCE_INLINE
1221 : Set::Matrix4<dim, Set::Sym::Full>
1222 : DoubleHessian(const amrex::Array4<const Set::Scalar>& f,
1223 : const int& i, const int& j, const int& k, const int& m,
1224 : const Set::Scalar dx[AMREX_SPACEDIM]);
1225 :
1226 : template<>
1227 : [[nodiscard]] AMREX_FORCE_INLINE
1228 : Set::Matrix4<2, Set::Sym::Full>
1229 : DoubleHessian<2>(const amrex::Array4<const Set::Scalar>& f,
1230 : const int& i, const int& j, const int& k, const int& m,
1231 : const Set::Scalar dx[AMREX_SPACEDIM])
1232 : {
1233 1362712 : Set::Matrix4<2, Set::Sym::Full> ret;
1234 : // [0,0,0,0]
1235 2725424 : ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1236 : // [0, 0, 0, 1]
1237 2725424 : ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1238 : // [0, 0, 1, 1]
1239 2725424 : ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1240 : // [0, 1, 1, 1]
1241 2725424 : ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1242 : // [1, 1, 1, 1]
1243 1362712 : ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1244 1362712 : return ret;
1245 : }
1246 :
1247 : template<>
1248 : [[nodiscard]] AMREX_FORCE_INLINE
1249 : Set::Matrix4<3, Set::Sym::Full>
1250 : DoubleHessian<3>(const amrex::Array4<const Set::Scalar>& f,
1251 : const int& i, const int& j, const int& k, const int& m,
1252 : const Set::Scalar dx[AMREX_SPACEDIM])
1253 : {
1254 0 : Set::Matrix4<3, Set::Sym::Full> ret;
1255 : // [0,0,0,0]
1256 0 : ret(0, 0, 0, 0) = Stencil<Set::Scalar, 4, 0, 0>::D(f, i, j, k, m, dx);
1257 : // [0, 0, 0, 1]
1258 0 : ret(0, 0, 0, 1) = Stencil<Set::Scalar, 3, 1, 0>::D(f, i, j, k, m, dx);
1259 : // [0, 0, 0, 2]
1260 0 : ret(0, 0, 0, 2) = Stencil<Set::Scalar, 3, 0, 1>::D(f, i, j, k, m, dx);
1261 : // [0, 0, 1, 1]
1262 0 : ret(0, 0, 1, 1) = Stencil<Set::Scalar, 2, 2, 0>::D(f, i, j, k, m, dx);
1263 : // [0, 0, 1, 2]
1264 0 : ret(0, 0, 1, 2) = Stencil<Set::Scalar, 2, 1, 1>::D(f, i, j, k, m, dx);
1265 : // [0, 0, 2, 2]
1266 0 : ret(0, 0, 2, 2) = Stencil<Set::Scalar, 2, 0, 2>::D(f, i, j, k, m, dx);
1267 : // [0, 1, 1, 1]
1268 0 : ret(0, 1, 1, 1) = Stencil<Set::Scalar, 1, 3, 0>::D(f, i, j, k, m, dx);
1269 : // [0, 1, 1, 2]
1270 0 : ret(0, 1, 1, 2) = Stencil<Set::Scalar, 1, 2, 1>::D(f, i, j, k, m, dx);
1271 : // [0, 1, 2, 2]
1272 0 : ret(0, 1, 2, 2) = Stencil<Set::Scalar, 1, 1, 2>::D(f, i, j, k, m, dx);
1273 : // [0, 2, 2, 2]
1274 0 : ret(0, 2, 2, 2) = Stencil<Set::Scalar, 1, 0, 3>::D(f, i, j, k, m, dx);
1275 : // [1, 1, 1, 1]
1276 0 : ret(1, 1, 1, 1) = Stencil<Set::Scalar, 0, 4, 0>::D(f, i, j, k, m, dx);
1277 : // [1, 1, 1, 2]
1278 0 : ret(1, 1, 1, 2) = Stencil<Set::Scalar, 0, 3, 1>::D(f, i, j, k, m, dx);
1279 : // [1, 1, 2, 2]
1280 0 : ret(1, 1, 2, 2) = Stencil<Set::Scalar, 0, 2, 2>::D(f, i, j, k, m, dx);
1281 : // [1, 2, 2, 2]
1282 0 : ret(1, 2, 2, 2) = Stencil<Set::Scalar, 0, 1, 3>::D(f, i, j, k, m, dx);
1283 : // [2, 2, 2, 2]
1284 0 : ret(2, 2, 2, 2) = Stencil<Set::Scalar, 0, 0, 4>::D(f, i, j, k, m, dx);
1285 0 : return ret;
1286 : }
1287 :
1288 : struct Interpolate
1289 : {
1290 : public:
1291 : template<class T>
1292 : [[nodiscard]] AMREX_FORCE_INLINE
1293 : static T CellToNodeAverage(const amrex::Array4<const T>& f,
1294 : const int& i, const int& j, const int& k, const int& m,
1295 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1296 : {
1297 183397 : int AMREX_D_DECL(
1298 : ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1299 : jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1300 : klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1301 :
1302 733588 : return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1303 : ,
1304 : +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1305 : ,
1306 : +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1307 : + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1308 183397 : )) * fac;
1309 : }
1310 : template<class T>
1311 : [[nodiscard]] AMREX_FORCE_INLINE
1312 : static T CellToNodeAverage(const amrex::Array4<T>& f,
1313 : const int& i, const int& j, const int& k, const int& m,
1314 : std::array<StencilType, AMREX_SPACEDIM> stencil = DefaultType)
1315 : {
1316 4620094 : int AMREX_D_DECL(
1317 : ilo = (stencil[0] == Numeric::StencilType::Lo ? 0 : 1),
1318 : jlo = (stencil[1] == Numeric::StencilType::Lo ? 0 : 1),
1319 : klo = (stencil[2] == Numeric::StencilType::Lo ? 0 : 1));
1320 :
1321 18480376 : return (AMREX_D_TERM(f(i, j, k, m) + f(i - ilo, j, k, m)
1322 : ,
1323 : +f(i, j - jlo, k, m) + f(i - ilo, j - jlo, k, m)
1324 : ,
1325 : +f(i, j, k - klo, m) + f(i - ilo, j, k - klo, m)
1326 : + f(i, j - jlo, k - klo, m) + f(i - ilo, j - jlo, k - klo, m)
1327 4620094 : )) * fac;
1328 : }
1329 : template<class T>
1330 : [[nodiscard]] AMREX_FORCE_INLINE
1331 : static T NodeToCellAverage(const amrex::Array4<const T>& f,
1332 : const int& i, const int& j, const int& k, const int& m)
1333 : {
1334 21086176 : return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1335 : ,
1336 : +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1337 : ,
1338 : +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1339 : + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1340 5271544 : )) * fac;
1341 : }
1342 : template<class T>
1343 : [[nodiscard]] AMREX_FORCE_INLINE
1344 : static T NodeToCellAverage(const amrex::Array4<T>& f,
1345 : const int& i, const int& j, const int& k, const int& m)
1346 : {
1347 : return (AMREX_D_TERM(f(i, j, k, m) + f(i + 1, j, k, m)
1348 : ,
1349 : +f(i, j + 1, k, m) + f(i + 1, j + 1, k, m)
1350 : ,
1351 : +f(i, j, k + 1, m) + f(i + 1, j, k + 1, m)
1352 : + f(i, j + 1, k + 1, m) + f(i + 1, j + 1, k + 1, m)
1353 : )) * fac;
1354 : }
1355 : constexpr static Set::Scalar fac = AMREX_D_PICK(0.5, 0.25, 0.125);
1356 : };
1357 :
1358 : }
1359 : #endif
|
