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