Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_EBMultiFabUtil_2D_C.H
Go to the documentation of this file.
1#ifndef AMREX_EB_MULTIFAB_UTIL_2D_C_H_
2#define AMREX_EB_MULTIFAB_UTIL_2D_C_H_
3#include <AMReX_Config.H>
4#include <AMReX_REAL.H>
5
6namespace amrex {
7
9void eb_set_covered_nodes (int i, int j, int k, int n, int icomp, Array4<Real> const& d,
10 Array4<EBCellFlag const> const& f, Real v)
11{
12 if (f(i-1,j-1,k).isCovered() && f(i ,j-1,k).isCovered() &&
13 f(i-1,j ,k).isCovered() && f(i ,j ,k).isCovered())
14 {
15 d(i,j,k,n+icomp) = v;
16 }
17}
18
20void eb_set_covered_nodes (int i, int j, int k, int n, int icomp, Array4<Real> const& d,
21 Array4<EBCellFlag const> const& f, Real const * AMREX_RESTRICT v)
22{
23 if (f(i-1,j-1,k).isCovered() && f(i ,j-1,k).isCovered() &&
24 f(i-1,j ,k).isCovered() && f(i ,j ,k).isCovered())
25 {
26 d(i,j,k,n+icomp) = v[n];
27 }
28}
29
31void eb_avgdown_with_vol (int i, int j, int k,
32 Array4<Real const> const& fine, int fcomp,
33 Array4<Real> const& crse, int ccomp,
34 Array4<Real const> const& fv, Array4<Real const> const& vfrc,
35 Dim3 const& ratio, int ncomp)
36{
37 for (int n = 0; n < ncomp; ++n) {
38 Real c = 0.0_rt;
39 Real cv = 0.0_rt;
40 constexpr int kk = 0;
41 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
42 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
43 Real tmp = fv(ii,jj,kk)*vfrc(ii,jj,kk);
44 c += fine(ii,jj,kk,n+fcomp)*tmp;
45 cv += tmp;
46 }}
47 if (cv > Real(1.e-30)) {
48 crse(i,j,k,n+ccomp) = c/cv;
49 } else {
50 crse(i,j,k,n+ccomp) = fine(i*ratio.x,j*ratio.y,kk,n+fcomp);
51 }
52 }
53}
54
56void eb_avgdown (int i, int j, int k,
57 Array4<Real const> const& fine, int fcomp,
58 Array4<Real> const& crse, int ccomp,
59 Array4<Real const> const& vfrc,
60 Dim3 const& ratio, int ncomp)
61{
62 for (int n = 0; n < ncomp; ++n) {
63 Real c = 0.0_rt;
64 Real cv = 0.0_rt;
65 constexpr int kk = 0;
66 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
67 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
68 Real tmp = vfrc(ii,jj,kk);
69 c += fine(ii,jj,kk,n+fcomp)*tmp;
70 cv += tmp;
71 }}
72 if (cv > Real(1.e-30)) {
73 crse(i,j,k,n+ccomp) = c/cv;
74 } else {
75 crse(i,j,k,n+ccomp) = fine(i*ratio.x,j*ratio.y,kk,n+fcomp);
76 }
77 }
78}
79
81void eb_avgdown_face_x (int i, int j, int k,
82 Array4<Real const> const& fine, int fcomp,
83 Array4<Real> const& crse, int ccomp,
84 Array4<Real const> const& area,
85 Dim3 const& ratio, int ncomp)
86{
87 int ii = i*ratio.x;
88 constexpr int kk = 0;
89 for (int n = 0; n < ncomp; ++n) {
90 Real c = 0.0_rt;
91 Real a = 0.0_rt;
92 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
93 Real tmp = area(ii,jj,kk);
94 c += tmp*fine(ii,jj,kk,n+fcomp);
95 a += tmp;
96 }
97 if (a > Real(1.e-30)) {
98 crse(i,j,k,n+ccomp) = c/a;
99 } else {
100 crse(i,j,k,n+ccomp) = fine(ii,j*ratio.y,kk,n+fcomp);
101 }
102 }
103}
104
106void eb_avgdown_face_y (int i, int j, int k,
107 Array4<Real const> const& fine, int fcomp,
108 Array4<Real> const& crse, int ccomp,
109 Array4<Real const> const& area,
110 Dim3 const& ratio, int ncomp)
111{
112 int jj = j*ratio.y;
113 constexpr int kk = 0;
114 for (int n = 0; n < ncomp; ++n) {
115 Real c = 0.0_rt;
116 Real a = 0.0_rt;
117 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
118 Real tmp = area(ii,jj,kk);
119 c += tmp*fine(ii,jj,kk,n+fcomp);
120 a += tmp;
121 }
122 if (a > Real(1.e-30)) {
123 crse(i,j,k,n+ccomp) = c/a;
124 } else {
125 crse(i,j,k,n+ccomp) = fine(i*ratio.x,jj,kk,n+fcomp);
126 }
127 }
128}
129
131void eb_avgdown_boundaries (int i, int j, int k,
132 Array4<Real const> const& fine, int fcomp,
133 Array4<Real> const& crse, int ccomp,
134 Array4<Real const> const& ba,
135 Dim3 const& ratio, int ncomp)
136{
137 for (int n = 0; n < ncomp; ++n) {
138 Real c = 0.0_rt;
139 Real cv = 0.0_rt;
140 constexpr int kk = 0;
141 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
142 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
143 Real tmp = ba(ii,jj,kk);
144 c += fine(ii,jj,kk,n+fcomp)*tmp;
145 cv += tmp;
146 }}
147 if (cv > Real(1.e-30)) {
148 crse(i,j,k,n+ccomp) = c/cv;
149 } else {
150 crse(i,j,k,n+ccomp) = 0.0_rt;
151 }
152 }
153}
154
156void eb_compute_divergence (int i, int j, int k, int n, Array4<Real> const& divu,
157 Array4<Real const> const& u, Array4<Real const> const& v,
158 Array4<int const> const& ccm, Array4<EBCellFlag const> const& flag,
159 Array4<Real const> const& vfrc, Array4<Real const> const& apx,
160 Array4<Real const> const& apy, Array4<Real const> const& fcx,
161 Array4<Real const> const& fcy, GpuArray<Real,2> const& dxinv,
162 bool already_on_centroids)
163{
164 if (flag(i,j,k).isCovered())
165 {
166 divu(i,j,k,n) = 0.0_rt;
167 }
168 else if (flag(i,j,k).isRegular())
169 {
170 divu(i,j,k,n) = dxinv[0] * (u(i+1,j,k,n)-u(i,j,k,n))
171 + dxinv[1] * (v(i,j+1,k,n)-v(i,j,k,n));
172 }
173 else if (already_on_centroids)
174 {
175 divu(i,j,k,n) = (1.0_rt/vfrc(i,j,k)) *
176 ( dxinv[0] * (apx(i+1,j,k)*u(i+1,j,k,n)-apx(i,j,k)*u(i,j,k,n))
177 + dxinv[1] * (apy(i,j+1,k)*v(i,j+1,k,n)-apy(i,j,k)*v(i,j,k,n)) );
178 }
179 else
180 {
181 Real fxm = u(i,j,k,n);
182 if (apx(i,j,k) != 0.0_rt && apx(i,j,k) != 1.0_rt) {
183 int jj = j + static_cast<int>(std::copysign(1.0_rt,fcx(i,j,k)));
184 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k)) : 0.0_rt;
185 fxm = (1.0_rt-fracy)*fxm + fracy*u(i,jj,k,n);
186 }
187
188 Real fxp = u(i+1,j,k,n);
189 if (apx(i+1,j,k) != 0.0_rt && apx(i+1,j,k) != 1.0_rt) {
190 int jj = j + static_cast<int>(std::copysign(1.0_rt,fcx(i+1,j,k)));
191 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k)) : 0.0_rt;
192 fxp = (1.0_rt-fracy)*fxp + fracy*u(i+1,jj,k,n);
193 }
194
195 Real fym = v(i,j,k,n);
196 if (apy(i,j,k) != 0.0_rt && apy(i,j,k) != 1.0_rt) {
197 int ii = i + static_cast<int>(std::copysign(1.0_rt,fcy(i,j,k)));
198 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k)) : 0.0_rt;
199 fym = (1.0_rt-fracx)*fym + fracx*v(ii,j,k,n);
200 }
201
202 Real fyp = v(i,j+1,k,n);
203 if (apy(i,j+1,k) != 0.0_rt && apy(i,j+1,k) != 1.0_rt) {
204 int ii = i + static_cast<int>(std::copysign(1.0_rt,fcy(i,j+1,k)));
205 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k)) : 0.0_rt;
206 fyp = (1.0_rt-fracx)*fyp + fracx*v(ii,j+1,k,n);
207 }
208
209 divu(i,j,k,n) = (1.0_rt/vfrc(i,j,k)) *
210 ( dxinv[0] * (apx(i+1,j,k)*fxp-apx(i,j,k)*fxm)
211 + dxinv[1] * (apy(i,j+1,k)*fyp-apy(i,j,k)*fym) );
212 }
213}
214
216void eb_avg_fc_to_cc (int i, int j, int k, int n, Array4<Real> const& cc,
217 Array4<Real const> const& fx, Array4<Real const> const& fy,
218 Array4<Real const> const& ax, Array4<Real const> const& ay,
219 Array4<EBCellFlag const> const& flag)
220{
221 if (flag(i,j,k).isCovered()) {
222 cc(i,j,k,n+0) = 0.0_rt;
223 cc(i,j,k,n+1) = 0.0_rt;
224 } else {
225 if (ax(i,j,k) == 0.0_rt) {
226 cc(i,j,k,n+0) = fx(i+1,j,k);
227 } else if (ax(i+1,j,k) == 0.0_rt) {
228 cc(i,j,k,n+0) = fx(i,j,k);
229 } else {
230 cc(i,j,k,n+0) = 0.5_rt * (fx(i,j,k) + fx(i+1,j,k));
231 }
232
233 if (ay(i,j,k) == 0.0_rt) {
234 cc(i,j,k,n+1) = fy(i,j+1,k);
235 } else if (ay(i,j+1,k) == 0.0_rt) {
236 cc(i,j,k,n+1) = fy(i,j,k);
237 } else {
238 cc(i,j,k,n+1) = 0.5_rt * (fy(i,j,k) + fy(i,j+1,k));
239 }
240 }
241}
242
244void eb_interp_cc2cent (Box const& box,
245 const Array4<Real>& phicent,
246 Array4<Real const> const& phicc,
247 Array4<EBCellFlag const> const& flag,
248 Array4<Real const> const& cent,
249 int ncomp) noexcept
250{
251 amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
252 {
253 if (flag(i,j,k).isCovered())
254 {
255 phicent(i,j,k,n) = phicc(i,j,k,n); //do nothing
256 }
257 else
258 {
259 if (flag(i,j,k).isRegular())
260 {
261 phicent(i,j,k,n) = phicc(i,j,k,n);
262 }
263 else
264 {
265 Real gx = cent(i,j,k,0);
266 Real gy = cent(i,j,k,1);
267 int ii = (gx < 0.0_rt) ? i - 1 : i + 1;
268 int jj = (gy < 0.0_rt) ? j - 1 : j + 1;
269 gx = std::abs(gx);
270 gy = std::abs(gy);
271 Real gxy = gx*gy;
272
273 phicent(i,j,k,n) = ( 1.0_rt - gx - gy + gxy ) * phicc(i ,j ,k ,n)
274 + ( gy - gxy ) * phicc(i ,jj,k ,n)
275 + ( gx - gxy ) * phicc(ii,j ,k ,n)
276 + ( gxy ) * phicc(ii,jj,k ,n);
277 }
278 }
279 });
280}
281
284 Array4<Real const> const& phi,
285 Array4<Real const> const& apx,
286 Array4<Real const> const& fcx,
287 Array4<Real> const& edg_x,
288 int ncomp,
289 const Box& domain,
290 const BCRec* bc) noexcept
291{
292 const Dim3 domlo = amrex::lbound(domain);
293 const Dim3 domhi = amrex::ubound(domain);
294 amrex::Loop(ubx, ncomp, [=] (int i, int j, int k, int n) noexcept
295 {
296 if (apx(i,j,k) == 0)
297 {
298 edg_x(i,j,k,n) = 1e40;
299 }
300 else
301 {
302 if (apx(i,j,k) == Real(1.))
303 {
304 if ( (i == domlo.x) && (bc[n].lo(0) == BCType::ext_dir) )
305 {
306 edg_x(i,j,k,n) = phi(domlo.x-1,j,k,n);
307 }
308 else if ( (i == domhi.x+1) && (bc[n].hi(0) == BCType::ext_dir) )
309 {
310 edg_x(i,j,k,n) = phi(domhi.x+1,j,k,n);
311 }
312 else
313 {
314 edg_x(i,j,k,n) = 0.5_rt * ( phi(i,j,k,n) + phi(i-1,j,k,n) );
315 }
316 }
317 else
318 {
319 int ii = i - 1;
320 int jj = (fcx(i,j,k) < 0.0_rt) ? j - 1 : j + 1;
321 Real gx = 0.5_rt;
322 Real gy = std::abs(fcx(i,j,k));
323 Real gxy = gx*gy;
324 edg_x(i,j,k,n) = ( 1.0_rt - gx - gy + gxy ) * phi(i ,j ,k ,n)
325 + ( gy - gxy ) * phi(i ,jj,k ,n)
326 + ( gx - gxy ) * phi(ii,j ,k ,n)
327 + ( gxy ) * phi(ii,jj,k ,n);
328 }
329 }
330 });
331}
332
335 Array4<Real const> const& phi,
336 Array4<Real const> const& apy,
337 Array4<Real const> const& fcy,
338 Array4<Real> const& edg_y,
339 int ncomp,
340 const Box& domain,
341 const BCRec* bc) noexcept
342{
343 const Dim3 domlo = amrex::lbound(domain);
344 const Dim3 domhi = amrex::ubound(domain);
345 amrex::Loop(vbx, ncomp, [=] (int i, int j, int k, int n) noexcept
346 {
347 if (apy(i,j,k) == 0)
348 {
349 edg_y(i,j,k,n) = 1e40;
350 }
351 else
352 {
353 if (apy(i,j,k) == Real(1.))
354 {
355 if ( (j == domlo.y) && (bc[n].lo(1) == BCType::ext_dir) )
356 {
357 edg_y(i,j,k,n) = phi(i,domlo.y-1,k,n);
358 }
359 else if ( (j == domhi.y+1) && (bc[n].hi(1) == BCType::ext_dir) )
360 {
361 edg_y(i,j,k,n) = phi(i,domhi.y+1,k,n);
362 }
363 else
364 {
365 edg_y(i,j,k,n) = 0.5_rt * (phi(i,j ,k,n) + phi(i,j-1,k,n));
366 }
367 }
368 else
369 {
370 int ii = (fcy(i,j,k) < 0.0_rt) ? i - 1 : i + 1;
371 int jj = j - 1;
372 Real gx = std::abs(fcy(i,j,k));
373 Real gy = 0.5_rt;
374 Real gxy = gx*gy;
375 edg_y(i,j,k,n) = ( 1.0_rt - gx - gy + gxy ) * phi(i ,j ,k ,n)
376 + ( gy - gxy ) * phi(i ,jj,k ,n)
377 + ( gx - gxy ) * phi(ii,j ,k ,n)
378 + ( gxy ) * phi(ii,jj,k ,n);
379 }
380 }
381 });
382}
383
386 Array4<Real const> const& phi,
387 Array4<Real const> const& apx,
388 Array4<Real const> const& cvol,
389 Array4<Real const> const& ccent,
390 Array4<Real const> const& fcx,
391 Array4<Real> const& edg_x,
392 int ncomp,
393 const Box& domain,
394 const BCRec* bc) noexcept
395{
396 const Dim3 domlo = amrex::lbound(domain);
397 const Dim3 domhi = amrex::ubound(domain);
398 amrex::Loop(ubx, ncomp, [=] (int i, int j, int k, int n) noexcept
399 {
400 if (apx(i,j,k) == 0)
401 {
402 edg_x(i,j,k,n) = 1e40;
403 }
404 else if ( (i == domlo.x) && (bc[n].lo(0) == BCType::ext_dir) )
405 {
406 edg_x(i,j,k,n) = phi(domlo.x-1,j,k,n);
407 }
408 else if ( (i == domhi.x+1) && (bc[n].hi(0) == BCType::ext_dir) )
409 {
410 edg_x(i,j,k,n) = phi(domhi.x+1,j,k,n);
411 }
412 else if ( apx(i,j,k) == Real(1.) && (cvol(i,j,k) == Real(1.) && cvol(i-1,j,k) == Real(1.)) )
413 {
414 edg_x(i,j,k,n) = 0.5_rt * ( phi(i,j,k,n) + phi(i-1,j,k,n) );
415 }
416 else if (apx(i,j,k) == Real(1.) && (std::abs(ccent(i,j,k,1)-ccent(i-1,j,k,1)) < Real(1.e-8)) )
417 {
418 Real d0 = 0.5_rt + ccent(i ,j,k,0); // distance in x-dir from centroid(i ,j) to face(i,j)
419 Real d1 = 0.5_rt - ccent(i-1,j,k,0); // distance in x-dir from centroid(i-1,j) to face(i,j)
420 Real a0 = d1 / (d0 + d1); // coefficient of phi(i ,j,k,n)
421 Real a1 = d0 / (d0 + d1); // coefficient of phi(i-1,j,k,n)
422 edg_x(i,j,k,n) = a0*phi(i,j,k,n) + a1*phi(i-1,j,k,n); // interpolation of phi in x-direction
423 }
424 else
425 {
426 int ii = i - 1;
427 int jj;
428 if (std::abs(fcx(i,j,k)) > Real(1.e-8)) {
429 jj = (fcx(i,j,k) < 0.0_rt) ? j - 1 : j + 1;
430 } else {
431
432 // We must add an additional test to avoid the case where fcx is very close to zero,
433 // but the cells at (j-1) or (j+1) might be covered or very small cells
434 Real min_vol_lo = std::min(cvol(i,j-1,k),cvol(ii,j-1,k));
435 Real min_vol_hi = std::min(cvol(i,j+1,k),cvol(ii,j+1,k));
436 jj = (min_vol_lo > min_vol_hi) ? j - 1 : j + 1;
437 }
438
439 Real d0 = 0.5_rt + ccent( i,j,k,0); // distance in x-dir from centroid(i ,j) to face(i,j)
440 Real d1 = 0.5_rt - ccent(ii,j,k,0); // distance in x-dir from centroid(i-1,j) to face(i,j)
441
442 Real d0_jj = 0.5_rt + ccent( i,jj,k,0); // distance in x-dir from centroid(i ,jj) to face(i,jj)
443 Real d1_jj = 0.5_rt - ccent(ii,jj,k,0); // distance in x-dir from centroid(i-1,jj) to face(i,jj)
444
445 Real a0 = d1 / (d0 + d1); // coefficient of phi(i ,j,k,n)
446 Real a1 = d0 / (d0 + d1); // coefficient of phi(i-1,j,k,n)
447
448 Real a0_jj = d1_jj / (d0_jj + d1_jj); // coefficient of phi(i ,jj,k,n)
449 Real a1_jj = d0_jj / (d0_jj + d1_jj); // coefficient of phi(i-1,jj,k,n)
450
451 Real phi01 = a0 * phi(i, j,k,n) + a1 * phi(ii, j,k,n); // interpolation of phi in x-direction
452 Real y01 = a0 * ccent(i, j,k,1) + a1 * ccent(ii, j,k,1); // interpolation of y-loc in x-direction
453
454 Real phi01_jj = a0_jj * phi(i,jj,k,n) + a1_jj * phi(ii,jj,k,n); // interpolation of phi in x-direction
455 Real y01_jj = a0_jj * ccent(i,jj,k,1) + a1_jj * ccent(ii,jj,k,1); // interpolation of y-loc in x-direction
456
457 edg_x(i,j,k,n) = ( ( fcx(i,j,k) - y01 ) * phi01_jj +
458 (Real(1.) - fcx(i,j,k) + y01_jj) * phi01 ) / (Real(1.) - y01 + y01_jj);
459 }
460 });
461}
462
465 Array4<Real const> const& phi,
466 Array4<Real const> const& apy,
467 Array4<Real const> const& cvol,
468 Array4<Real const> const& ccent,
469 Array4<Real const> const& fcy,
470 Array4<Real> const& edg_y,
471 int ncomp,
472 const Box& domain,
473 const BCRec* bc) noexcept
474{
475 const Dim3 domlo = amrex::lbound(domain);
476 const Dim3 domhi = amrex::ubound(domain);
477 amrex::Loop(vbx, ncomp, [=] (int i, int j, int k, int n) noexcept
478 {
479 if (apy(i,j,k) == 0)
480 {
481 edg_y(i,j,k,n) = 1e40;
482 }
483 else if ( (j == domlo.y) && (bc[n].lo(1) == BCType::ext_dir) )
484 {
485 edg_y(i,j,k,n) = phi(i,domlo.y-1,k,n);
486 }
487 else if ( (j == domhi.y+1) && (bc[n].hi(1) == BCType::ext_dir) )
488 {
489 edg_y(i,j,k,n) = phi(i,domhi.y+1,k,n);
490 }
491 else if (apy(i,j,k) == Real(1.) && cvol(i,j,k) == Real(1.) && cvol(i,j-1,k) == Real(1.))
492 {
493 edg_y(i,j,k,n) = 0.5_rt * (phi(i,j ,k,n) + phi(i,j-1,k,n));
494 }
495 else if (apy(i,j,k) == Real(1.) && (std::abs(ccent(i,j,k,0)-ccent(i,j-1,k,0)) < Real(1.e-8)) )
496 {
497 Real d0 = 0.5_rt + ccent(i,j ,k,1); // distance in y-dir from centroid(i,j ) to face(i,j)
498 Real d1 = 0.5_rt - ccent(i,j-1,k,1); // distance in y-dir from centroid(i,j-1) to face(i,j)
499 Real a0 = d1 / (d0 + d1); // coefficient of phi(i,j ,k,n)
500 Real a1 = d0 / (d0 + d1); // coefficient of phi(i,j-1,k,n)
501 edg_y(i,j,k,n) = a0*phi(i,j,k,n) + a1*phi(i,j-1,k,n); // interpolation of phi in y-direction
502 }
503 else
504 {
505 int jj = j - 1;
506 int ii;
507
508 // We must add an additional test to avoid the case where fcy is very close to zero, but (i-1) or (i+1)
509 // might be covered cells -- this can happen when the EB is exactly aligned with the horizontal grid lines
510 if (std::abs(fcy(i,j,k)) > Real(1.e-8)) {
511 ii = (fcy(i,j,k) < 0.0_rt) ? i - 1 : i + 1;
512 } else {
513 // We must add an additional test to avoid the case where fcx is very close to zero,
514 // but the cells at (i-1) or (i+1) might be covered or very small cells
515 Real min_vol_lo = std::min(cvol(i-1,j,k),cvol(i-1,jj,k));
516 Real min_vol_hi = std::min(cvol(i+1,j,k),cvol(i+1,jj,k));
517 ii = (min_vol_lo > min_vol_hi) ? i - 1 : i + 1;
518 }
519
520 Real d0 = 0.5_rt + ccent(i, j,k,1); // distance in y-dir from centroid(i,j ) to face(i,j)
521 Real d1 = 0.5_rt - ccent(i,jj,k,1); // distance in y-dir from centroid(i,j-1) to face(i,j)
522
523 Real d0_ii = 0.5_rt + ccent(ii, j,k,0); // distance from centroid(ii,j ) to face(ii,j)
524 Real d1_ii = 0.5_rt - ccent(ii,jj,k,0); // distance from centroid(ii,j-1) to face(ii,j)
525
526 Real a0 = d1 / (d0 + d1); // coefficient of phi(i,j ,k,n)
527 Real a1 = d0 / (d0 + d1); // coefficient of phi(i,j-1,k,n)
528
529 Real a0_ii = d1_ii / (d0_ii + d1_ii); // coefficient of phi(i ,jj,k,n)
530 Real a1_ii = d0_ii / (d0_ii + d1_ii); // coefficient of phi(i-1,jj,k,n)
531
532 Real phi01 = a0 * phi( i,j,k,n) + a1 * phi( i,jj,k,n); // interpolation of phi in y-direction
533 Real x01 = a0 * ccent( i,j,k,0) + a1 * ccent( i,jj,k,0); // interpolation of x-loc in y-direction
534
535 Real phi01_ii = a0_ii * phi(ii,j,k,n) + a1_ii * phi(ii,jj,k,n); // interpolation of phi in y-direction
536 Real x01_ii = a0_ii * ccent(ii,j,k,0) + a1_ii * ccent(ii,jj,k,0); // interpolation of x-loc in y-direction
537
538 edg_y(i,j,k,n) = ( ( fcy(i,j,k) - x01 ) * phi01_ii +
539 (Real(1.) - fcy(i,j,k) + x01_ii) * phi01 ) / (Real(1.) - x01 + x01_ii);
540 }
541 });
542}
543
545void eb_interp_cc2face_x (Box const& ubx,
546 Array4<Real const> const& phi,
547 Array4<Real> const& edg_x,
548 int ncomp,
549 const Box& domain,
550 const BCRec* bc) noexcept
551{
552 const Dim3 domlo = amrex::lbound(domain);
553 const Dim3 domhi = amrex::ubound(domain);
554 amrex::Loop(ubx, ncomp, [=] (int i, int j, int k, int n) noexcept
555 {
556 if ( (i == domlo.x) && (bc[n].lo(0) == BCType::ext_dir) )
557 {
558 edg_x(i,j,k,n) = phi(domlo.x-1,j,k,n);
559 }
560 else if ( (i == domhi.x+1) && (bc[n].hi(0) == BCType::ext_dir) )
561 {
562 edg_x(i,j,k,n) = phi(domhi.x+1,j,k,n);
563 }
564 else
565 {
566 edg_x(i,j,k,n) = 0.5_rt * ( phi(i,j,k,n) + phi(i-1,j,k,n) );
567 }
568 });
569}
570
572void eb_interp_cc2face_y (Box const& vbx,
573 Array4<Real const> const& phi,
574 Array4<Real> const& edg_y,
575 int ncomp,
576 const Box& domain,
577 const BCRec* bc) noexcept
578{
579 const Dim3 domlo = amrex::lbound(domain);
580 const Dim3 domhi = amrex::ubound(domain);
581 amrex::Loop(vbx, ncomp, [=] (int i, int j, int k, int n) noexcept
582 {
583 if ( (j == domlo.y) && (bc[n].lo(1) == BCType::ext_dir) )
584 {
585 edg_y(i,j,k,n) = phi(i,domlo.y-1,k,n);
586 }
587 else if ( (j == domhi.y+1) && (bc[n].hi(1) == BCType::ext_dir) )
588 {
589 edg_y(i,j,k,n) = phi(i,domhi.y+1,k,n);
590 }
591 else
592 {
593 edg_y(i,j,k,n) = 0.5_rt * (phi(i,j ,k,n) + phi(i,j-1,k,n));
594 }
595 });
596}
597
599void eb_add_divergence_from_flow (int i, int j, int k, int n, Array4<Real> const& divu,
600 Array4<Real const> const& vel_eb, Array4<EBCellFlag const> const& flag,
601 Array4<Real const> const& vfrc, Array4<Real const> const& bnorm,
602 Array4<Real const> const& barea, GpuArray<Real,2> const& dxinv)
603{
604 if (flag(i,j,k).isSingleValued())
605 {
606 Real Ueb_dot_n = vel_eb(i,j,k,0)*bnorm(i,j,k,0)
607 + vel_eb(i,j,k,1)*bnorm(i,j,k,1);
608
609 divu(i,j,k,n) += Ueb_dot_n * barea(i,j,k) * dxinv[0] / vfrc(i,j,k);
610 }
611}
612
613}
614
615#endif
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_RESTRICT
Definition AMReX_Extension.H:37
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
Array4< Real const > crse
Definition AMReX_InterpFaceRegister.cpp:92
if(!(yy_init))
Definition amrex_iparser.lex.nolint.H:935
Boundary Condition Records. Necessary information and functions for computing boundary conditions.
Definition AMReX_BCRec.H:17
Definition AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_avgdown_face_y(int i, int j, int k, Array4< Real const > const &fine, int fcomp, Array4< Real > const &crse, int ccomp, Array4< Real const > const &area, Dim3 const &ratio, int ncomp)
Definition AMReX_EBMultiFabUtil_2D_C.H:106
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_centroid2facecent_y(Box const &vbx, Array4< Real const > const &phi, Array4< Real const > const &apy, Array4< Real const > const &cvol, Array4< Real const > const &ccent, Array4< Real const > const &fcy, Array4< Real > const &edg_y, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_2D_C.H:464
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_add_divergence_from_flow(int i, int j, int k, int n, Array4< Real > const &divu, Array4< Real const > const &vel_eb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &bnorm, Array4< Real const > const &barea, GpuArray< Real, 2 > const &dxinv)
Definition AMReX_EBMultiFabUtil_2D_C.H:599
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_compute_divergence(int i, int j, int k, int n, Array4< Real > const &divu, Array4< Real const > const &u, Array4< Real const > const &v, Array4< int const > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, GpuArray< Real, 2 > const &dxinv, bool already_on_centroids)
Definition AMReX_EBMultiFabUtil_2D_C.H:156
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_avgdown_face_x(int i, int j, int k, Array4< Real const > const &fine, int fcomp, Array4< Real > const &crse, int ccomp, Array4< Real const > const &area, Dim3 const &ratio, int ncomp)
Definition AMReX_EBMultiFabUtil_2D_C.H:81
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_avg_fc_to_cc(int i, int j, int k, int n, Array4< Real > const &cc, Array4< Real const > const &fx, Array4< Real const > const &fy, Array4< Real const > const &ax, Array4< Real const > const &ay, Array4< EBCellFlag const > const &flag)
Definition AMReX_EBMultiFabUtil_2D_C.H:216
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_avgdown_boundaries(int i, int j, int k, Array4< Real const > const &fine, int fcomp, Array4< Real > const &crse, int ccomp, Array4< Real const > const &ba, Dim3 const &ratio, int ncomp)
Definition AMReX_EBMultiFabUtil_2D_C.H:131
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_avgdown(int i, int j, int k, Array4< Real const > const &fine, int fcomp, Array4< Real > const &crse, int ccomp, Array4< Real const > const &vfrc, Dim3 const &ratio, int ncomp)
Definition AMReX_EBMultiFabUtil_2D_C.H:56
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_avgdown_with_vol(int i, int j, int k, Array4< Real const > const &fine, int fcomp, Array4< Real > const &crse, int ccomp, Array4< Real const > const &fv, Array4< Real const > const &vfrc, Dim3 const &ratio, int ncomp)
Definition AMReX_EBMultiFabUtil_2D_C.H:31
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_set_covered_nodes(int i, int j, int k, int n, int icomp, Array4< Real > const &d, Array4< EBCellFlag const > const &f, Real v)
Definition AMReX_EBMultiFabUtil_2D_C.H:9
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_cc2facecent_y(Box const &vbx, Array4< Real const > const &phi, Array4< Real const > const &apy, Array4< Real const > const &fcy, Array4< Real > const &edg_y, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_2D_C.H:334
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_cc2face_x(Box const &ubx, Array4< Real const > const &phi, Array4< Real > const &edg_x, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_2D_C.H:545
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_cc2face_y(Box const &vbx, Array4< Real const > const &phi, Array4< Real > const &edg_y, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_2D_C.H:572
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_cc2cent(Box const &box, const Array4< Real > &phicent, Array4< Real const > const &phicc, Array4< EBCellFlag const > const &flag, Array4< Real const > const &cent, int ncomp) noexcept
Definition AMReX_EBMultiFabUtil_2D_C.H:244
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_cc2facecent_x(Box const &ubx, Array4< Real const > const &phi, Array4< Real const > const &apx, Array4< Real const > const &fcx, Array4< Real > const &edg_x, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_2D_C.H:283
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_centroid2facecent_x(Box const &ubx, Array4< Real const > const &phi, Array4< Real const > const &apx, Array4< Real const > const &cvol, Array4< Real const > const &ccent, Array4< Real const > const &fcx, Array4< Real > const &edg_x, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_2D_C.H:385
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void Loop(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:125
Definition AMReX_Array4.H:61
Definition AMReX_Dim3.H:12
int x
Definition AMReX_Dim3.H:12
int y
Definition AMReX_Dim3.H:12
Definition AMReX_Array.H:34