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