Line data Source code
1 : #ifndef AMREX_MULTIFAB_UTIL_3D_C_H_
2 : #define AMREX_MULTIFAB_UTIL_3D_C_H_
3 : #include <AMReX_Config.H>
4 :
5 : #include <AMReX_Gpu.H>
6 : #include <AMReX_FArrayBox.H>
7 : #include <AMReX_IArrayBox.H>
8 : #include <cmath>
9 :
10 : namespace amrex {
11 :
12 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
13 : void amrex_avg_nd_to_cc (int i, int j, int k, int n,
14 : Array4<Real > const& cc,
15 : Array4<Real const> const& nd,
16 : int cccomp, int ndcomp) noexcept
17 : {
18 : cc(i,j,k,n+cccomp) = Real(0.125)*( nd(i,j ,k ,n+ndcomp) + nd(i+1,j ,k ,n+ndcomp)
19 : + nd(i,j+1,k ,n+ndcomp) + nd(i+1,j+1,k ,n+ndcomp)
20 : + nd(i,j ,k+1,n+ndcomp) + nd(i+1,j ,k+1,n+ndcomp)
21 : + nd(i,j+1,k+1,n+ndcomp) + nd(i+1,j+1,k+1,n+ndcomp));
22 : }
23 :
24 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
25 : void amrex_avg_eg_to_cc (int i, int j, int k,
26 : Array4<Real > const& cc,
27 : Array4<Real const> const& Ex,
28 : Array4<Real const> const& Ey,
29 : Array4<Real const> const& Ez,
30 : int cccomp) noexcept
31 : {
32 : cc(i,j,k,0+cccomp) = Real(0.25) * ( Ex(i,j,k) + Ex(i,j+1,k) + Ex(i,j,k+1) + Ex(i,j+1,k+1) );
33 : cc(i,j,k,1+cccomp) = Real(0.25) * ( Ey(i,j,k) + Ey(i+1,j,k) + Ey(i,j,k+1) + Ey(i+1,j,k+1) );
34 : cc(i,j,k,2+cccomp) = Real(0.25) * ( Ez(i,j,k) + Ez(i+1,j,k) + Ez(i,j+1,k) + Ez(i+1,j+1,k) );
35 : }
36 :
37 : template <typename CT, typename FT>
38 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
39 : void amrex_avg_fc_to_cc (int i, int j, int k,
40 : Array4<CT > const& cc,
41 : Array4<FT const> const& fx,
42 : Array4<FT const> const& fy,
43 : Array4<FT const> const& fz,
44 : int cccomp) noexcept
45 : {
46 : cc(i,j,k,0+cccomp) = CT(0.5) * CT( fx(i,j,k) + fx(i+1,j,k) );
47 : cc(i,j,k,1+cccomp) = CT(0.5) * CT( fy(i,j,k) + fy(i,j+1,k) );
48 : cc(i,j,k,2+cccomp) = CT(0.5) * CT( fz(i,j,k) + fz(i,j,k+1) );
49 : }
50 :
51 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
52 : void amrex_avg_cc_to_fc (int i, int j, int k, int n, Box const& xbx, Box const& ybx, Box const& zbx,
53 : Array4<Real> const& fx, Array4<Real> const& fy,
54 : Array4<Real> const& fz, Array4<Real const> const& cc,
55 : bool use_harmonic_averaging) noexcept
56 : {
57 : if (use_harmonic_averaging)
58 : {
59 : if (xbx.contains(i,j,k)) {
60 : if (cc(i-1,j,k,n) == Real(0.0) || cc(i,j,k,n) == Real(0.0)) {
61 : fx(i,j,k,n) = Real(0.0);
62 : } else {
63 : fx(i,j,k,n) = Real(2.0) / (Real(1.0)/cc(i-1,j,k,n) + Real(1.0)/cc(i,j,k,n));
64 : }
65 : }
66 : if (ybx.contains(i,j,k)) {
67 : if (cc(i,j-1,k,n) == Real(0.0) || cc(i,j,k,n) == Real(0.0)) {
68 : fy(i,j,k,n) = Real(0.0);
69 : } else {
70 : fy(i,j,k,n) = Real(2.0) / (Real(1.0)/cc(i,j-1,k,n) + Real(1.0)/cc(i,j,k,n));
71 : }
72 : }
73 : if (zbx.contains(i,j,k)) {
74 : if (cc(i,j,k-1,n) == Real(0.0) || cc(i,j,k,n) == Real(0.0)) {
75 : fz(i,j,k,n) = Real(0.0);
76 : } else {
77 : fz(i,j,k,n) = Real(2.0) / (Real(1.0)/cc(i,j,k-1,n) + Real(1.0)/cc(i,j,k,n));
78 : }
79 : }
80 : } else {
81 : if (xbx.contains(i,j,k)) {
82 : fx(i,j,k,n) = Real(0.5)*(cc(i-1,j,k,n) + cc(i,j,k,n));
83 : }
84 : if (ybx.contains(i,j,k)) {
85 : fy(i,j,k,n) = Real(0.5)*(cc(i,j-1,k,n) + cc(i,j,k,n));
86 : }
87 : if (zbx.contains(i,j,k)) {
88 : fz(i,j,k,n) = Real(0.5)*(cc(i,j,k-1,n) + cc(i,j,k,n));
89 : }
90 : }
91 : }
92 :
93 : template <typename T>
94 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
95 : void amrex_avgdown_faces (Box const& bx, Array4<T> const& crse,
96 : Array4<T const> const& fine,
97 : int ccomp, int fcomp, int ncomp,
98 : IntVect const& ratio, int idir) noexcept
99 : {
100 : const auto clo = lbound(bx);
101 : const auto chi = ubound(bx);
102 : const int facx = ratio[0];
103 : const int facy = ratio[1];
104 : const int facz = ratio[2];
105 :
106 : switch (idir) {
107 : case 0:
108 : {
109 : T facInv = T(1.0) / (facy*facz);
110 : for (int n = 0; n < ncomp; ++n) {
111 : for (int k = clo.z; k <= chi.z; ++k) {
112 : for (int j = clo.y; j <= chi.y; ++j) {
113 : for (int i = clo.x; i <= chi.x; ++i) {
114 : int ii = i*facx;
115 : int jj = j*facy;
116 : int kk = k*facz;
117 : T c = T(0.);
118 : for (int kref = 0; kref < facz; ++kref) {
119 : for (int jref = 0; jref < facy; ++jref) {
120 : c += fine(ii,jj+jref,kk+kref,n+fcomp);
121 : }}
122 : crse(i,j,k,n+ccomp) = c * facInv;
123 : }}}
124 : }
125 : break;
126 : }
127 : case 1:
128 : {
129 : T facInv = T(1.0) / (facx*facz);
130 : for (int n = 0; n < ncomp; ++n) {
131 : for (int k = clo.z; k <= chi.z; ++k) {
132 : for (int j = clo.y; j <= chi.y; ++j) {
133 : for (int i = clo.x; i <= chi.x; ++i) {
134 : int ii = i*facx;
135 : int jj = j*facy;
136 : int kk = k*facz;
137 : T c = T(0.);
138 : for (int kref = 0; kref < facz; ++kref) {
139 : for (int iref = 0; iref < facx; ++iref) {
140 : c += fine(ii+iref,jj,kk+kref,n+fcomp);
141 : }}
142 : crse(i,j,k,n+ccomp) = c * facInv;
143 : }}}
144 : }
145 : break;
146 : }
147 : case 2:
148 : {
149 : T facInv = T(1.0) / (facx*facy);
150 : for (int n = 0; n < ncomp; ++n) {
151 : for (int k = clo.z; k <= chi.z; ++k) {
152 : for (int j = clo.y; j <= chi.y; ++j) {
153 : for (int i = clo.x; i <= chi.x; ++i) {
154 : int ii = i*facx;
155 : int jj = j*facy;
156 : int kk = k*facz;
157 : T c = T(0.);
158 : for (int jref = 0; jref < facy; ++jref) {
159 : for (int iref = 0; iref < facx; ++iref) {
160 : c += fine(ii+iref,jj+jref,kk,n+fcomp);
161 : }}
162 : crse(i,j,k,n+ccomp) = c * facInv;
163 : }}}
164 : }
165 : break;
166 : }
167 : default: { break; }
168 : }
169 : }
170 :
171 : template <typename T>
172 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
173 : void amrex_avgdown_faces (int i, int j, int k, int n, Array4<T> const& crse,
174 : Array4<T const> const& fine,
175 : int ccomp, int fcomp, IntVect const& ratio, int idir) noexcept
176 : {
177 : const int facx = ratio[0];
178 : const int facy = ratio[1];
179 : const int facz = ratio[2];
180 : const int ii = i*facx;
181 : const int jj = j*facy;
182 : const int kk = k*facz;
183 :
184 : switch (idir) {
185 : case 0:
186 : {
187 : const T facInv = T(1.0) / (facy*facz);
188 : T c = T(0.);
189 : for (int kref = 0; kref < facz; ++kref) {
190 : for (int jref = 0; jref < facy; ++jref) {
191 : c += fine(ii,jj+jref,kk+kref,n+fcomp);
192 : }}
193 : crse(i,j,k,n+ccomp) = c * facInv;
194 : break;
195 : }
196 : case 1:
197 : {
198 : const T facInv = T(1.0) / (facx*facz);
199 : T c = T(0.);
200 : for (int kref = 0; kref < facz; ++kref) {
201 : for (int iref = 0; iref < facx; ++iref) {
202 : c += fine(ii+iref,jj,kk+kref,n+fcomp);
203 : }}
204 : crse(i,j,k,n+ccomp) = c * facInv;
205 : break;
206 : }
207 : case 2:
208 : {
209 : const T facInv = T(1.0) / (facx*facy);
210 : T c = T(0.);
211 : for (int jref = 0; jref < facy; ++jref) {
212 : for (int iref = 0; iref < facx; ++iref) {
213 : c += fine(ii+iref,jj+jref,kk,n+fcomp);
214 : }}
215 : crse(i,j,k,n+ccomp) = c * facInv;
216 : break;
217 : }
218 : default: { break; }
219 : }
220 : }
221 :
222 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
223 : void amrex_avgdown_edges (Box const& bx, Array4<Real> const& crse,
224 : Array4<Real const> const& fine,
225 : int ccomp, int fcomp, int ncomp,
226 : IntVect const& ratio, int idir) noexcept
227 : {
228 : const auto clo = lbound(bx);
229 : const auto chi = ubound(bx);
230 : const int facx = ratio[0];
231 : const int facy = ratio[1];
232 : const int facz = ratio[2];
233 :
234 : switch (idir) {
235 : case 0:
236 : {
237 : Real facInv = Real(1.0) / static_cast<Real>(facx);
238 : for (int n = 0; n < ncomp; ++n) {
239 : for (int k = clo.z; k <= chi.z; ++k) {
240 : for (int j = clo.y; j <= chi.y; ++j) {
241 : for (int i = clo.x; i <= chi.x; ++i) {
242 : int ii = i*facx;
243 : int jj = j*facy;
244 : int kk = k*facz;
245 : Real c = 0.;
246 : for (int iref = 0; iref < facx; ++iref) {
247 : c += fine(ii+iref,jj,kk,n+fcomp);
248 : }
249 : crse(i,j,k,n+ccomp) = c * facInv;
250 : }}}
251 : }
252 : break;
253 : }
254 : case 1:
255 : {
256 : Real facInv = Real(1.0) / static_cast<Real>(facy);
257 : for (int n = 0; n < ncomp; ++n) {
258 : for (int k = clo.z; k <= chi.z; ++k) {
259 : for (int j = clo.y; j <= chi.y; ++j) {
260 : for (int i = clo.x; i <= chi.x; ++i) {
261 : int ii = i*facx;
262 : int jj = j*facy;
263 : int kk = k*facz;
264 : Real c = 0.;
265 : for (int jref = 0; jref < facy; ++jref) {
266 : c += fine(ii,jj+jref,kk,n+fcomp);
267 : }
268 : crse(i,j,k,n+ccomp) = c * facInv;
269 : }}}
270 : }
271 : break;
272 : }
273 : case 2:
274 : {
275 : Real facInv = Real(1.0) / static_cast<Real>(facz);
276 : for (int n = 0; n < ncomp; ++n) {
277 : for (int k = clo.z; k <= chi.z; ++k) {
278 : for (int j = clo.y; j <= chi.y; ++j) {
279 : for (int i = clo.x; i <= chi.x; ++i) {
280 : int ii = i*facx;
281 : int jj = j*facy;
282 : int kk = k*facz;
283 : Real c = 0.;
284 : for (int kref = 0; kref < facz; ++kref) {
285 : c += fine(ii,jj,kk+kref,n+fcomp);
286 : }
287 : crse(i,j,k,n+ccomp) = c * facInv;
288 : }}}
289 : }
290 : break;
291 : }
292 : default: { break; }
293 : }
294 : }
295 :
296 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
297 : void amrex_avgdown_edges (int i, int j, int k, int n, Array4<Real> const& crse,
298 : Array4<Real const> const& fine,
299 : int ccomp, int fcomp, IntVect const& ratio, int idir) noexcept
300 : {
301 : const int facx = ratio[0];
302 : const int facy = ratio[1];
303 : const int facz = ratio[2];
304 : const int ii = i*facx;
305 : const int jj = j*facy;
306 : const int kk = k*facz;
307 :
308 : switch (idir) {
309 : case 0:
310 : {
311 : const Real facInv = Real(1.0) / static_cast<Real>(facx);
312 : Real c = 0.;
313 : for (int iref = 0; iref < facx; ++iref) {
314 : c += fine(ii+iref,jj,kk,n+fcomp);
315 : }
316 : crse(i,j,k,n+ccomp) = c * facInv;
317 : break;
318 : }
319 : case 1:
320 : {
321 : const Real facInv = Real(1.0) / static_cast<Real>(facy);
322 : Real c = 0.;
323 : for (int jref = 0; jref < facy; ++jref) {
324 : c += fine(ii,jj+jref,kk,n+fcomp);
325 : }
326 : crse(i,j,k,n+ccomp) = c * facInv;
327 : break;
328 : }
329 : case 2:
330 : {
331 : const Real facInv = Real(1.0) / static_cast<Real>(facz);
332 : Real c = 0.;
333 : for (int kref = 0; kref < facz; ++kref) {
334 : c += fine(ii,jj,kk+kref,n+fcomp);
335 : }
336 : crse(i,j,k,n+ccomp) = c * facInv;
337 : break;
338 : }
339 : default: { break; }
340 : }
341 : }
342 :
343 : template <typename T>
344 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
345 : void amrex_avgdown (Box const& bx, Array4<T> const& crse,
346 : Array4<T const> const& fine,
347 : int ccomp, int fcomp, int ncomp,
348 : IntVect const& ratio) noexcept
349 : {
350 : const auto clo = lbound(bx);
351 : const auto chi = ubound(bx);
352 : const int facx = ratio[0];
353 : const int facy = ratio[1];
354 : const int facz = ratio[2];
355 : const T volfrac = T(1.0)/T(facx*facy*facz);
356 :
357 : for (int n = 0; n < ncomp; ++n) {
358 : for (int k = clo.z; k <= chi.z; ++k) {
359 : for (int j = clo.y; j <= chi.y; ++j) {
360 : for (int i = clo.x; i <= chi.x; ++i) {
361 : int ii = i*facx;
362 : int jj = j*facy;
363 : int kk = k*facz;
364 : T c = 0;
365 : for (int kref = 0; kref < facz; ++kref) {
366 : for (int jref = 0; jref < facy; ++jref) {
367 : for (int iref = 0; iref < facx; ++iref) {
368 : c += fine(ii+iref,jj+jref,kk+kref,n+fcomp);
369 : }}}
370 : crse(i,j,k,n+ccomp) = volfrac * c;
371 : }}}
372 : }
373 : }
374 :
375 : template <typename T>
376 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
377 : void amrex_avgdown (int i, int j, int k, int n, Array4<T> const& crse,
378 : Array4<T const> const& fine,
379 : int ccomp, int fcomp, IntVect const& ratio) noexcept
380 : {
381 119120 : const int facx = ratio[0];
382 119120 : const int facy = ratio[1];
383 119120 : const int facz = ratio[2];
384 119120 : const T volfrac = T(1.0)/T(facx*facy*facz);
385 119120 : const int ii = i*facx;
386 119120 : const int jj = j*facy;
387 119120 : const int kk = k*facz;
388 119120 : T c = 0;
389 357360 : for (int kref = 0; kref < facz; ++kref) {
390 714720 : for (int jref = 0; jref < facy; ++jref) {
391 1429440 : for (int iref = 0; iref < facx; ++iref) {
392 1905920 : c += fine(ii+iref,jj+jref,kk+kref,n+fcomp);
393 : }}}
394 119120 : crse(i,j,k,n+ccomp) = volfrac * c;
395 119120 : }
396 :
397 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
398 : void amrex_avgdown_with_vol (int i, int j, int k, int n, Array4<Real> const& crse,
399 : Array4<Real const> const& fine,
400 : Array4<Real const> const& fv,
401 : int ccomp, int fcomp, IntVect const& ratio) noexcept
402 : {
403 : const int facx = ratio[0];
404 : const int facy = ratio[1];
405 : const int facz = ratio[2];
406 : const int ii = i*facx;
407 : const int jj = j*facy;
408 : const int kk = k*facz;
409 : Real cd = 0.0, cv = 0.0;
410 : for (int kref = 0; kref < facz; ++kref) {
411 : for (int jref = 0; jref < facy; ++jref) {
412 : for (int iref = 0; iref < facx; ++iref) {
413 : cv += fv(ii+iref,jj+jref,kk+kref);
414 : cd += fine(ii+iref,jj+jref,kk+kref,n+fcomp)*fv(ii+iref,jj+jref,kk+kref);
415 : }}}
416 : crse(i,j,k,n+ccomp) = cd/cv;
417 : }
418 :
419 : template <typename T>
420 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
421 : void amrex_avgdown_nodes (Box const& bx, Array4<T> const& crse,
422 : Array4<T const> const& fine,
423 : int ccomp, int fcomp, int ncomp,
424 : IntVect const& ratio) noexcept
425 : {
426 : const auto clo = lbound(bx);
427 : const auto chi = ubound(bx);
428 0 : const int facx = ratio[0];
429 0 : const int facy = ratio[1];
430 0 : const int facz = ratio[2];
431 :
432 0 : for (int n = 0; n < ncomp; ++n) {
433 0 : for (int k = clo.z; k <= chi.z; ++k) {
434 0 : int kk = k*facz;
435 0 : for (int j = clo.y; j <= chi.y; ++j) {
436 0 : int jj = j*facy;
437 : AMREX_PRAGMA_SIMD
438 0 : for (int i = clo.x; i <= chi.x; ++i) {
439 0 : crse(i,j,k,n+ccomp) = fine(i*facx,jj,kk,n+fcomp);
440 : }
441 : }
442 : }
443 : }
444 0 : }
445 :
446 : template<typename T>
447 : AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
448 : void amrex_avgdown_nodes (int i, int j, int k, int n, Array4<T> const& crse,
449 : Array4<T const> const& fine,
450 : int ccomp, int fcomp, IntVect const& ratio) noexcept
451 : {
452 1931850 : crse(i,j,k,n+ccomp) = fine(i*ratio[0],j*ratio[1],k*ratio[2],n+fcomp);
453 386370 : }
454 :
455 : AMREX_GPU_HOST_DEVICE
456 : inline
457 : void amrex_compute_divergence (Box const& bx, Array4<Real> const& divu,
458 : Array4<Real const> const& u,
459 : Array4<Real const> const& v,
460 : Array4<Real const> const& w,
461 : GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
462 : {
463 : const auto lo = lbound(bx);
464 : const auto hi = ubound(bx);
465 : const Real dxi = dxinv[0];
466 : const Real dyi = dxinv[1];
467 : const Real dzi = dxinv[2];
468 :
469 : for (int n = 0; n < divu.ncomp; ++n) {
470 : for (int k = lo.z; k <= hi.z; ++k) {
471 : for (int j = lo.y; j <= hi.y; ++j) {
472 : AMREX_PRAGMA_SIMD
473 : for (int i = lo.x; i <= hi.x; ++i) {
474 : divu(i,j,k,n) = dxi * (u(i+1,j,k,n)-u(i,j,k,n))
475 : + dyi * (v(i,j+1,k,n)-v(i,j,k,n))
476 : + dzi * (w(i,j,k+1,n)-w(i,j,k,n));
477 : }
478 : }
479 : }
480 : }
481 : }
482 :
483 : AMREX_GPU_HOST_DEVICE
484 : inline
485 : void amrex_compute_gradient (Box const& bx, Array4<Real> const& grad,
486 : Array4<Real const> const& u,
487 : Array4<Real const> const& v,
488 : Array4<Real const> const& w,
489 : GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
490 : {
491 : const auto lo = lbound(bx);
492 : const auto hi = ubound(bx);
493 : const Real dxi = dxinv[0];
494 : const Real dyi = dxinv[1];
495 : const Real dzi = dxinv[2];
496 :
497 : for (int k = lo.z; k <= hi.z; ++k) {
498 : for (int j = lo.y; j <= hi.y; ++j) {
499 : AMREX_PRAGMA_SIMD
500 : for (int i = lo.x; i <= hi.x; ++i) {
501 : grad(i,j,k,0) = dxi * (u(i+1,j,k,0)-u(i,j,k,0));
502 : grad(i,j,k,1) = dyi * (v(i,j+1,k,0)-v(i,j,k,0));
503 : grad(i,j,k,2) = dzi * (w(i,j,k+1,0)-w(i,j,k,0));
504 : }
505 : }
506 : }
507 : }
508 :
509 : AMREX_GPU_HOST_DEVICE
510 : inline
511 : void amrex_compute_convective_difference (Box const& bx, Array4<Real> const& diff,
512 : Array4<Real const> const& u_face,
513 : Array4<Real const> const& v_face,
514 : Array4<Real const> const& w_face,
515 : Array4<Real const> const& s_on_x_face,
516 : Array4<Real const> const& s_on_y_face,
517 : Array4<Real const> const& s_on_z_face,
518 : GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
519 : {
520 : const auto lo = lbound(bx);
521 : const auto hi = ubound(bx);
522 : const Real dxi = dxinv[0];
523 : const Real dyi = dxinv[1];
524 : const Real dzi = dxinv[2];
525 :
526 : for (int n = 0; n < diff.ncomp; ++n) {
527 : for (int k = lo.z; k <= hi.z; ++k) {
528 : for (int j = lo.y; j <= hi.y; ++j) {
529 : AMREX_PRAGMA_SIMD
530 : for (int i = lo.x; i <= hi.x; ++i) {
531 : diff(i,j,k,n) = Real(0.5)*dxi * (u_face(i+1,j,k)+u_face(i,j,k)) *
532 : (s_on_x_face(i+1,j,k,n)-s_on_x_face(i,j,k,n))
533 : + Real(0.5)*dyi * (v_face(i,j+1,k)+v_face(i,j,k)) *
534 : (s_on_y_face(i,j+1,k,n)-s_on_y_face(i,j,k,n))
535 : + Real(0.5)*dzi * (w_face(i,j,k+1)+w_face(i,j,k)) *
536 : (s_on_z_face(i,j,k+1,n)-s_on_z_face(i,j,k,n));
537 : }
538 : }
539 : }
540 : }
541 : }
542 :
543 : } // namespace amrex
544 :
545 : #endif
|
