Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_EBMultiFabUtil_3D_C.H
Go to the documentation of this file.
1#ifndef AMREX_EB_MULTIFAB_UTIL_3D_C_H_
2#define AMREX_EB_MULTIFAB_UTIL_3D_C_H_
3#include <AMReX_Config.H>
4#include <AMReX_REAL.H>
5
6namespace amrex {
7
8AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
9Real
10EB_interp_in_quad(Real xint,Real yint,Real v0,Real v1,Real v2,Real v3,
11 Real x0,Real y0, Real x1,Real y1,
12 Real x2,Real y2, Real x3,Real y3)
13{
14//
15// 2D Stencil
16//
17// v3 @ (x3,y3) v2 @ (x2,y2)
18//
19// (xint,yint)
20//
21// v0 @ (x0,y0) v1 @ (x1,y1)
22
23 Real va = v3 - v0;
24 Real vb = v2 - v0;
25 Real vc = v1 - v0;
26
27 Real xa = x3 - x0;
28 Real xb = x2 - x0;
29 Real xc = x1 - x0;
30
31 Real ya = y3 - y0;
32 Real yb = y2 - y0;
33 Real yc = y1 - y0;
34
35 Real xder = vc*(-xa + xb)*ya*yb + vb*(xa - xc)*ya*yc + va*(-xb + xc)*yb*yc;
36 Real yder = vc*xa*xb*(ya - yb) + va*xb*xc*(yb - yc) + vb*xa*xc*(-ya + yc);
37 Real xyder = -(vc*xb*ya) + vb*xc*ya + vc*xa*yb - va*xc*yb - vb*xa*yc + va*xb*yc ;
38
39 Real det = -(xa*xc*ya*yb) + xb*xc*ya*yb + xa*xb*ya*yc - xb*xc*ya*yc - xa*xb*yb*yc + xa*xc*yb*yc ;
40
41 xder = xder / det;
42 yder = yder / det;
43 xyder = xyder / det;
44
45 Real value = v0 + xder*(xint-x0)+yder*(yint-y0)+ xyder*(xint-x0)*(yint-y0);
46
47 return value;
48}
49
50AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
51void eb_set_covered_nodes (int i, int j, int k, int n, int icomp, Array4<Real> const& d,
52 Array4<EBCellFlag const> const& f, Real v)
53{
54 if (f(i-1,j-1,k-1).isCovered() && f(i ,j-1,k-1).isCovered() &&
55 f(i-1,j ,k-1).isCovered() && f(i ,j ,k-1).isCovered() &&
56 f(i-1,j-1,k ).isCovered() && f(i ,j-1,k ).isCovered() &&
57 f(i-1,j ,k ).isCovered() && f(i ,j ,k ).isCovered())
58 {
59 d(i,j,k,n+icomp) = v;
60 }
61}
62
63AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
64void eb_set_covered_nodes (int i, int j, int k, int n, int icomp, Array4<Real> const& d,
65 Array4<EBCellFlag const> const& f, Real const * AMREX_RESTRICT v)
66{
67 if (f(i-1,j-1,k-1).isCovered() && f(i ,j-1,k-1).isCovered() &&
68 f(i-1,j ,k-1).isCovered() && f(i ,j ,k-1).isCovered() &&
69 f(i-1,j-1,k ).isCovered() && f(i ,j-1,k ).isCovered() &&
70 f(i-1,j ,k ).isCovered() && f(i ,j ,k ).isCovered())
71 {
72 d(i,j,k,n+icomp) = v[n];
73 }
74}
75
76AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
77void eb_avgdown_with_vol (int i, int j, int k,
78 Array4<Real const> const& fine, int fcomp,
79 Array4<Real> const& crse, int ccomp,
80 Array4<Real const> const& fv, Array4<Real const> const& vfrc,
81 Dim3 const& ratio, int ncomp)
82{
83 for (int n = 0; n < ncomp; ++n) {
84 Real c = Real(0.0);
85 Real cv = Real(0.0);
86 for (int kk = k*ratio.z; kk < (k+1)*ratio.z; ++kk) {
87 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
88 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
89 Real tmp = fv(ii,jj,kk)*vfrc(ii,jj,kk);
90 c += fine(ii,jj,kk,n+fcomp)*tmp;
91 cv += tmp;
92 }}}
93 if (cv > Real(1.e-30)) {
94 crse(i,j,k,n+ccomp) = c/cv;
95 } else {
96 crse(i,j,k,n+ccomp) = fine(i*ratio.x,j*ratio.y,k*ratio.z,n+fcomp);
97 }
98 }
99}
100
101AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
102void eb_avgdown (int i, int j, int k,
103 Array4<Real const> const& fine, int fcomp,
104 Array4<Real> const& crse, int ccomp,
105 Array4<Real const> const& vfrc,
106 Dim3 const& ratio, int ncomp)
107{
108 for (int n = 0; n < ncomp; ++n) {
109 Real c = Real(0.0);
110 Real cv = Real(0.0);
111 for (int kk = k*ratio.z; kk < (k+1)*ratio.z; ++kk) {
112 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
113 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
114 Real tmp = vfrc(ii,jj,kk);
115 c += fine(ii,jj,kk,n+fcomp)*tmp;
116 cv += tmp;
117 }}}
118 if (cv > Real(1.e-30)) {
119 crse(i,j,k,n+ccomp) = c/cv;
120 } else {
121 crse(i,j,k,n+ccomp) = fine(i*ratio.x,j*ratio.y,k*ratio.z,n+fcomp);
122 }
123 }
124}
125
126AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
127void eb_avgdown_face_x (int i, int j, int k,
128 Array4<Real const> const& fine, int fcomp,
129 Array4<Real> const& crse, int ccomp,
130 Array4<Real const> const& area,
131 Dim3 const& ratio, int ncomp)
132{
133 int ii = i*ratio.x;
134 for (int n = 0; n < ncomp; ++n) {
135 Real c = Real(0.0);
136 Real a = Real(0.0);
137 for (int kk = k*ratio.z; kk < (k+1)*ratio.z; ++kk) {
138 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
139 Real tmp = area(ii,jj,kk);
140 c += tmp*fine(ii,jj,kk,n+fcomp);
141 a += tmp;
142 }}
143 if (a > Real(1.e-30)) {
144 crse(i,j,k,n+ccomp) = c/a;
145 } else {
146 crse(i,j,k,n+ccomp) = fine(ii,j*ratio.y,k*ratio.z,n+fcomp);
147 }
148 }
149}
150
151AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
152void eb_avgdown_face_y (int i, int j, int k,
153 Array4<Real const> const& fine, int fcomp,
154 Array4<Real> const& crse, int ccomp,
155 Array4<Real const> const& area,
156 Dim3 const& ratio, int ncomp)
157{
158 int jj = j*ratio.y;
159 for (int n = 0; n < ncomp; ++n) {
160 Real c = Real(0.0);
161 Real a = Real(0.0);
162 for (int kk = k*ratio.z; kk < (k+1)*ratio.z; ++kk) {
163 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
164 Real tmp = area(ii,jj,kk);
165 c += tmp*fine(ii,jj,kk,n+fcomp);
166 a += tmp;
167 }}
168 if (a > Real(1.e-30)) {
169 crse(i,j,k,n+ccomp) = c/a;
170 } else {
171 crse(i,j,k,n+ccomp) = fine(i*ratio.x,jj,k*ratio.z,n+fcomp);
172 }
173 }
174}
175
176AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
177void eb_avgdown_face_z (int i, int j, int k,
178 Array4<Real const> const& fine, int fcomp,
179 Array4<Real> const& crse, int ccomp,
180 Array4<Real const> const& area,
181 Dim3 const& ratio, int ncomp)
182{
183 int kk = k*ratio.z;
184 for (int n = 0; n < ncomp; ++n) {
185 Real c = Real(0.0);
186 Real a = Real(0.0);
187 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
188 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
189 Real tmp = area(ii,jj,kk);
190 c += tmp*fine(ii,jj,kk,n+fcomp);
191 a += tmp;
192 }}
193 if (a > Real(1.e-30)) {
194 crse(i,j,k,n+ccomp) = c/a;
195 } else {
196 crse(i,j,k,n+ccomp) = fine(i*ratio.x,j*ratio.y,kk,n+fcomp);
197 }
198 }
199}
200
201AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
202void eb_avgdown_boundaries (int i, int j, int k,
203 Array4<Real const> const& fine, int fcomp,
204 Array4<Real> const& crse, int ccomp,
205 Array4<Real const> const& ba,
206 Dim3 const& ratio, int ncomp)
207{
208 for (int n = 0; n < ncomp; ++n) {
209 Real c = Real(0.0);
210 Real cv = Real(0.0);
211 for (int kk = k*ratio.z; kk < (k+1)*ratio.z; ++kk) {
212 for (int jj = j*ratio.y; jj < (j+1)*ratio.y; ++jj) {
213 for (int ii = i*ratio.x; ii < (i+1)*ratio.x; ++ii) {
214 Real tmp = ba(ii,jj,kk);
215 c += fine(ii,jj,kk,n+fcomp)*tmp;
216 cv += tmp;
217 }}}
218 if (cv > Real(1.e-30)) {
219 crse(i,j,k,n+ccomp) = c/cv;
220 } else {
221 crse(i,j,k,n+ccomp) = Real(0.0);
222 }
223 }
224}
225
226AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
227void eb_compute_divergence (int i, int j, int k, int n, Array4<Real> const& divu,
228 Array4<Real const> const& u, Array4<Real const> const& v,
229 Array4<Real const> const& w, Array4<int const> const& ccm,
230 Array4<EBCellFlag const> const& flag, Array4<Real const> const& vfrc,
231 Array4<Real const> const& apx, Array4<Real const> const& apy,
232 Array4<Real const> const& apz, Array4<Real const> const& fcx,
233 Array4<Real const> const& fcy, Array4<Real const> const& fcz,
234 GpuArray<Real,3> const& dxinv, bool already_on_centroids)
235{
236 if (flag(i,j,k).isCovered())
237 {
238 divu(i,j,k,n) = Real(0.0);
239 }
240 else if (flag(i,j,k).isRegular())
241 {
242 divu(i,j,k,n) = dxinv[0] * (u(i+1,j,k,n)-u(i,j,k,n))
243 + dxinv[1] * (v(i,j+1,k,n)-v(i,j,k,n))
244 + dxinv[2] * (w(i,j,k+1,n)-w(i,j,k,n));
245 }
246 else if (already_on_centroids)
247 {
248 divu(i,j,k,n) = (Real(1.0)/vfrc(i,j,k)) * (
249 dxinv[0] * (apx(i+1,j,k)*u(i+1,j,k,n)-apx(i,j,k)*u(i,j,k,n))
250 + dxinv[1] * (apy(i,j+1,k)*v(i,j+1,k,n)-apy(i,j,k)*v(i,j,k,n))
251 + dxinv[2] * (apz(i,j,k+1)*w(i,j,k+1,n)-apz(i,j,k)*w(i,j,k,n)) );
252 }
253 else
254 {
255 Real fxm = u(i,j,k,n);
256 if (apx(i,j,k) != Real(0.0) && apx(i,j,k) != Real(1.0)) {
257 int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
258 int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
259 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
260 Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
261 fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm
262 + fracy *(Real(1.0)-fracz)*u(i,jj,k ,n)
263 + fracz *(Real(1.0)-fracy)*u(i,j ,kk,n)
264 + fracy * fracz *u(i,jj,kk,n);
265 }
266
267 Real fxp = u(i+1,j,k,n);
268 if (apx(i+1,j,k) != Real(0.0) && apx(i+1,j,k) != Real(1.0)) {
269 int jj = j + static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,0)));
270 int kk = k + static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,1)));
271 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k,0)) : Real(0.0);
272 Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk)) ? std::abs(fcx(i+1,j,k,1)) : Real(0.0);
273 fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp
274 + fracy *(Real(1.0)-fracz)*u(i+1,jj,k ,n)
275 + fracz *(Real(1.0)-fracy)*u(i+1,j ,kk,n)
276 + fracy * fracz *u(i+1,jj,kk,n);
277 }
278
279 Real fym = v(i,j,k,n);
280 if (apy(i,j,k) != Real(0.0) && apy(i,j,k) != Real(1.0)) {
281 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
282 int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
283 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
284 Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
285 fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym
286 + fracx *(Real(1.0)-fracz)*v(ii,j,k ,n)
287 + fracz *(Real(1.0)-fracx)*v(i ,j,kk,n)
288 + fracx * fracz *v(ii,j,kk,n);
289 }
290
291 Real fyp = v(i,j+1,k,n);
292 if (apy(i,j+1,k) != Real(0.0) && apy(i,j+1,k) != Real(1.0)) {
293 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,0)));
294 int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,1)));
295 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k,0)) : Real(0.0);
296 Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk)) ? std::abs(fcy(i,j+1,k,1)) : Real(0.0);
297 fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp
298 + fracx *(Real(1.0)-fracz)*v(ii,j+1,k ,n)
299 + fracz *(Real(1.0)-fracx)*v(i ,j+1,kk,n)
300 + fracx * fracz *v(ii,j+1,kk,n);
301 }
302
303 Real fzm = w(i,j,k,n);
304 if (apz(i,j,k) != Real(0.0) && apz(i,j,k) != Real(1.0)) {
305 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
306 int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
307 Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
308 Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
309 fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm
310 + fracx *(Real(1.0)-fracy)*w(ii,j ,k,n)
311 + fracy *(Real(1.0)-fracx)*w(i ,jj,k,n)
312 + fracx * fracy *w(ii,jj,k,n);
313 }
314
315 Real fzp = w(i,j,k+1,n);
316 if (apz(i,j,k+1) != Real(0.0) && apz(i,j,k+1) != Real(1.0)) {
317 int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,0)));
318 int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,1)));
319 Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1)) ? std::abs(fcz(i,j,k+1,0)) : Real(0.0);
320 Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1)) ? std::abs(fcz(i,j,k+1,1)) : Real(0.0);
321 fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp
322 + fracx *(Real(1.0)-fracy)*w(ii,j ,k+1,n)
323 + fracy *(Real(1.0)-fracx)*w(i ,jj,k+1,n)
324 + fracx * fracy *w(ii,jj,k+1,n);
325 }
326
327 divu(i,j,k,n) = (Real(1.0)/vfrc(i,j,k)) *
328 ( dxinv[0] * (apx(i+1,j,k)*fxp-apx(i,j,k)*fxm)
329 + dxinv[1] * (apy(i,j+1,k)*fyp-apy(i,j,k)*fym)
330 + dxinv[2] * (apz(i,j,k+1)*fzp-apz(i,j,k)*fzm) );
331 }
332}
333
334AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
335void eb_avg_fc_to_cc (int i, int j, int k, int n, Array4<Real> const& cc,
336 Array4<Real const> const& fx, Array4<Real const> const& fy,
337 Array4<Real const> const& fz, Array4<Real const> const& ax,
338 Array4<Real const> const& ay, Array4<Real const> const& az,
339 Array4<EBCellFlag const> const& flag)
340{
341 if (flag(i,j,k).isCovered()) {
342 cc(i,j,k,n+0) = Real(0.0);
343 cc(i,j,k,n+1) = Real(0.0);
344 cc(i,j,k,n+2) = Real(0.0);
345 } else {
346 if (ax(i,j,k) == Real(0.0)) {
347 cc(i,j,k,n+0) = fx(i+1,j,k);
348 } else if (ax(i+1,j,k) == Real(0.0)) {
349 cc(i,j,k,n+0) = fx(i,j,k);
350 } else {
351 cc(i,j,k,n+0) = Real(0.5) * (fx(i,j,k) + fx(i+1,j,k));
352 }
353
354 if (ay(i,j,k) == Real(0.0)) {
355 cc(i,j,k,n+1) = fy(i,j+1,k);
356 } else if (ay(i,j+1,k) == Real(0.0)) {
357 cc(i,j,k,n+1) = fy(i,j,k);
358 } else {
359 cc(i,j,k,n+1) = Real(0.5) * (fy(i,j,k) + fy(i,j+1,k));
360 }
361
362 if (az(i,j,k) == Real(0.0)) {
363 cc(i,j,k,n+2) = fz(i,j,k+1);
364 } else if (az(i,j,k+1) == Real(0.0)) {
365 cc(i,j,k,n+2) = fz(i,j,k);
366 } else {
367 cc(i,j,k,n+2) = Real(0.5) * (fz(i,j,k) + fz(i,j,k+1));
368 }
369 }
370}
371
372AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
373void eb_interp_cc2cent (Box const& box,
374 const Array4<Real>& phicent,
375 Array4<Real const > const& phicc,
376 Array4<EBCellFlag const> const& flag,
377 Array4<Real const> const& cent,
378 int ncomp) noexcept
379{
380 amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
381 {
382 if (flag(i,j,k).isCovered())
383 {
384 phicent(i,j,k,n) = phicc(i,j,k,n); //do nothing
385 }
386 else
387 {
388 if (flag(i,j,k).isRegular())
389 {
390 phicent(i,j,k,n) = phicc(i,j,k,n);
391 }
392 else
393 {
394 Real gx = cent(i,j,k,0);
395 Real gy = cent(i,j,k,1);
396 Real gz = cent(i,j,k,2);
397 int ii = (gx < Real(0.0)) ? i - 1 : i + 1;
398 int jj = (gy < Real(0.0)) ? j - 1 : j + 1;
399 int kk = (gz < Real(0.0)) ? k - 1 : k + 1;
400 gx = std::abs(gx);
401 gy = std::abs(gy);
402 gz = std::abs(gz);
403 Real gxy = gx*gy;
404 Real gxz = gx*gz;
405 Real gyz = gy*gz;
406 Real gxyz = gx*gy*gz;
407 phicent(i,j,k,n)
408 = ( Real(1.0) - gx - gy - gz + gxy + gxz + gyz - gxyz) * phicc(i ,j ,k ,n)
409 + ( gz - gxz - gyz + gxyz) * phicc(i ,j ,kk,n)
410 + ( gy - gxy - gyz + gxyz) * phicc(i ,jj,k ,n)
411 + ( gyz - gxyz) * phicc(i ,jj,kk,n)
412 + ( gx - gxy - gxz + gxyz) * phicc(ii,j ,k ,n)
413 + ( gxz - gxyz) * phicc(ii,j ,kk,n)
414 + ( gxy - gxyz) * phicc(ii,jj,k ,n)
415 + ( gxyz) * phicc(ii,jj,kk,n);
416 }
417 }
418 });
419}
420
421AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
422void eb_interp_cc2facecent_x (Box const& ubx,
423 Array4<Real const> const& phi,
424 Array4<Real const> const& apx,
425 Array4<Real const> const& fcx,
426 Array4<Real> const& edg_x,
427 int ncomp,
428 const Box& domain,
429 const BCRec* bc) noexcept
430{
431 const Dim3 domlo = amrex::lbound(domain);
432 const Dim3 domhi = amrex::ubound(domain);
433 amrex::Loop(ubx, ncomp, [=] (int i, int j, int k, int n) noexcept
434 {
435 if (apx(i,j,k) == 0)
436 {
437#ifdef AMREX_USE_FLOAT
438 edg_x(i,j,k,n) = Real(1e20);
439#else
440 edg_x(i,j,k,n) = 1e40;
441#endif
442 }
443 else
444 {
445 if (apx(i,j,k) == Real(1.))
446 {
447 if ( (i == domlo.x) && (bc[n].lo(0) == BCType::ext_dir) )
448 {
449 edg_x(i,j,k,n) = phi(domlo.x-1,j,k,n);
450 }
451 else if ( (i == domhi.x+1) && (bc[n].hi(0) == BCType::ext_dir) )
452 {
453 edg_x(i,j,k,n) = phi(domhi.x+1,j,k,n);
454 }
455 else
456 {
457 edg_x(i,j,k,n) = Real(0.5) * ( phi(i,j,k,n) + phi(i-1,j,k,n) );
458 }
459 }
460 else
461 {
462 Real gx = Real(0.5);
463 Real gy = fcx(i,j,k,0);
464 Real gz = fcx(i,j,k,1);
465 int ii = i - 1;
466 int jj = (gy < Real(0.0)) ? j - 1 : j + 1;
467 int kk = (gz < Real(0.0)) ? k - 1 : k + 1;
468 gy = std::abs(gy);
469 gz = std::abs(gz);
470 Real gxy = gx*gy;
471 Real gxz = gx*gz;
472 Real gyz = gy*gz;
473 Real gxyz = gx*gy*gz;
474 edg_x(i,j,k,n) =
475 ( Real(1.0) - gx - gy - gz + gxy + gxz + gyz - gxyz) * phi(i ,j ,k ,n)
476 + ( gz - gxz - gyz + gxyz) * phi(i ,j ,kk,n)
477 + ( gy - gxy - gyz + gxyz) * phi(i ,jj,k ,n)
478 + ( gyz - gxyz) * phi(i ,jj,kk,n)
479 + ( gx - gxy - gxz + gxyz) * phi(ii,j ,k ,n)
480 + ( gxz - gxyz) * phi(ii,j ,kk,n)
481 + ( gxy - gxyz) * phi(ii,jj,k ,n)
482 + ( gxyz) * phi(ii,jj,kk,n);
483 }
484 }
485 });
486}
487
488AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
489void eb_interp_cc2facecent_y (Box const& vbx,
490 Array4<Real const> const& phi,
491 Array4<Real const> const& apy,
492 Array4<Real const> const& fcy,
493 Array4<Real> const& edg_y,
494 int ncomp,
495 const Box& domain,
496 const BCRec* bc) noexcept
497{
498 const Dim3 domlo = amrex::lbound(domain);
499 const Dim3 domhi = amrex::ubound(domain);
500 amrex::Loop(vbx, ncomp, [=] (int i, int j, int k, int n) noexcept
501 {
502 if (apy(i,j,k) == 0)
503 {
504#ifdef AMREX_USE_FLOAT
505 edg_y(i,j,k,n) = Real(1e20);
506#else
507 edg_y(i,j,k,n) = 1e40;
508#endif
509 }
510 else
511 {
512 if (apy(i,j,k) == Real(1.))
513 {
514 if ( (j == domlo.y) && (bc[n].lo(1) == BCType::ext_dir) )
515 {
516 edg_y(i,j,k,n) = phi(i,domlo.y-1,k,n);
517 }
518 else if ( (j == domhi.y+1) && (bc[n].hi(1) == BCType::ext_dir) )
519 {
520 edg_y(i,j,k,n) = phi(i,domhi.y+1,k,n);
521 }
522 else
523 {
524 edg_y(i,j,k,n) = Real(0.5) * (phi(i,j ,k,n) + phi(i,j-1,k,n));
525 }
526 }
527 else
528 {
529 Real gx = fcy(i,j,k,0);
530 Real gy = Real(0.5);
531 Real gz = fcy(i,j,k,1);
532 int ii = (gx < Real(0.0)) ? i - 1 : i + 1;
533 int jj = j - 1;
534 int kk = (gz < Real(0.0)) ? k - 1 : k + 1;
535 gx = std::abs(gx);
536 gz = std::abs(gz);
537 Real gxy = gx*gy;
538 Real gxz = gx*gz;
539 Real gyz = gy*gz;
540 Real gxyz = gx*gy*gz;
541 edg_y(i,j,k,n) =
542 ( Real(1.0) - gx - gy - gz + gxy + gxz + gyz - gxyz) * phi(i ,j ,k ,n)
543 + ( gz - gxz - gyz + gxyz) * phi(i ,j ,kk,n)
544 + ( gy - gxy - gyz + gxyz) * phi(i ,jj,k ,n)
545 + ( gyz - gxyz) * phi(i ,jj,kk,n)
546 + ( gx - gxy - gxz + gxyz) * phi(ii,j ,k ,n)
547 + ( gxz - gxyz) * phi(ii,j ,kk,n)
548 + ( gxy - gxyz) * phi(ii,jj,k ,n)
549 + ( gxyz) * phi(ii,jj,kk,n);
550 }
551 }
552 });
553}
554
555AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
556void eb_interp_cc2facecent_z (Box const& wbx,
557 Array4<Real const> const& phi,
558 Array4<Real const> const& apz,
559 Array4<Real const> const& fcz,
560 Array4<Real> const& edg_z,
561 int ncomp,
562 const Box& domain,
563 const BCRec* bc) noexcept
564{
565 const Dim3 domlo = amrex::lbound(domain);
566 const Dim3 domhi = amrex::ubound(domain);
567 amrex::Loop(wbx, ncomp, [=] (int i, int j, int k, int n) noexcept
568 {
569 if (apz(i,j,k) == 0)
570 {
571#ifdef AMREX_USE_FLOAT
572 edg_z(i,j,k,n) = Real(1e20);
573#else
574 edg_z(i,j,k,n) = 1e40;
575#endif
576 }
577 else
578 {
579 if (apz(i,j,k) == Real(1.))
580 {
581 if ( (k == domlo.z) && (bc[n].lo(2) == BCType::ext_dir) )
582 {
583 edg_z(i,j,k,n) = phi(i,j,domlo.z-1,n);
584 }
585 else if ( (k == domhi.z+1) && (bc[n].hi(2) == BCType::ext_dir) )
586 {
587 edg_z(i,j,k,n) = phi(i,j,domhi.z+1,n);
588 }
589 else
590 {
591 edg_z(i,j,k,n) = Real(0.5) * (phi(i,j ,k,n) + phi(i,j,k-1,n));
592 }
593 }
594 else
595 {
596 Real gx = fcz(i,j,k,0);
597 Real gy = fcz(i,j,k,1);
598 Real gz = Real(0.5);
599 int ii = (gx < Real(0.0)) ? i - 1 : i + 1;
600 int jj = (gy < Real(0.0)) ? j - 1 : j + 1;
601 int kk = k -1;
602 gx = std::abs(gx);
603 gy = std::abs(gy);
604 Real gxy = gx*gy;
605 Real gxz = gx*gz;
606 Real gyz = gy*gz;
607 Real gxyz = gx*gy*gz;
608 edg_z(i,j,k,n) =
609 ( Real(1.0) - gx - gy - gz + gxy + gxz + gyz - gxyz) * phi(i ,j ,k ,n)
610 + ( gz - gxz - gyz + gxyz) * phi(i ,j ,kk,n)
611 + ( gy - gxy - gyz + gxyz) * phi(i ,jj,k ,n)
612 + ( gyz - gxyz) * phi(i ,jj,kk,n)
613 + ( gx - gxy - gxz + gxyz) * phi(ii,j ,k ,n)
614 + ( gxz - gxyz) * phi(ii,j ,kk,n)
615 + ( gxy - gxyz) * phi(ii,jj,k ,n)
616 + ( gxyz) * phi(ii,jj,kk,n);
617 }
618 }
619 });
620}
621
622AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
623void eb_interp_centroid2facecent_x (Box const& ubx,
624 Array4<Real const> const& phi,
625 Array4<Real const> const& apx,
626 Array4<Real const> const& cvol,
627 Array4<Real const> const& ccent,
628 Array4<Real const> const& fcx,
629 Array4<Real> const& phi_x,
630 int ncomp,
631 const Box& domain,
632 const BCRec* bc) noexcept
633{
634 const Dim3 domlo = amrex::lbound(domain);
635 const Dim3 domhi = amrex::ubound(domain);
636 // Note that ccent holds (x,y,z) of the cell centroids as components (0/1/2)
637 // fcx holds ( y,z) of the x-face centroid as components ( /0/1)
638 amrex::Loop(ubx, ncomp, [=] (int i, int j, int k, int n) noexcept
639 {
640 if (apx(i,j,k) == 0)
641 {
642#ifdef AMREX_USE_FLOAT
643 phi_x(i,j,k,n) = Real(1e20);
644#else
645 phi_x(i,j,k,n) = 1e40;
646#endif
647 }
648 else if ( (i == domlo.x) && (bc[n].lo(0) == BCType::ext_dir) )
649 {
650 phi_x(i,j,k,n) = phi(domlo.x-1,j,k,n);
651 }
652 else if ( (i == domhi.x+1) && (bc[n].hi(0) == BCType::ext_dir) )
653 {
654 phi_x(i,j,k,n) = phi(domhi.x+1,j,k,n);
655 }
656 else if (apx(i,j,k) == Real(1.) && cvol(i,j,k) == Real(1.) && cvol(i-1,j,k) == Real(1.))
657 {
658 phi_x(i,j,k,n) = Real(0.5) * ( phi(i,j,k,n) + phi(i-1,j,k,n) );
659 }
660 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)) &&
661 (std::abs(ccent(i,j,k,2)-ccent(i-1,j,k,2)) < Real(1.e-8)) )
662 {
663 Real d0 = Real(0.5) + ccent(i ,j,k,0); // distance in x-dir from centroid(i ,j) to face(i,j)
664 Real d1 = Real(0.5) - ccent(i-1,j,k,0); // distance in x-dir from centroid(i-1,j) to face(i,j)
665 Real a0 = d1 / (d0 + d1); // coefficient of phi(i ,j,k,n)
666 Real a1 = d0 / (d0 + d1); // coefficient of phi(i-1,j,k,n)
667 phi_x(i,j,k,n) = a0*phi(i,j,k,n) + a1*phi(i-1,j,k,n); // interpolation of phi in x-direction
668 }
669 else
670 {
671 // We must add additional tests to avoid the case where fcx is very close to zero, but j-1/j+1 or k-1/k+1
672 // might be covered cells -- this can happen when the EB is exactly aligned with the horizontal grid planes
673 int jj,kk;
674 if (std::abs(fcx(i,j,k,0)) > Real(1.e-8)) {
675 jj = (fcx(i,j,k,0) < Real(0.0)) ? j - 1 : j + 1;
676 } else if (apx(i,j-1,k) > Real(0.)) {
677 jj = j-1;
678 } else {
679 jj = j+1;
680 }
681
682 if (std::abs(fcx(i,j,k,1)) > Real(1.e-8)) {
683 kk = (fcx(i,j,k,1) < Real(0.0)) ? k - 1 : k + 1;
684 } else if (apx(i,j,k-1) > Real(0.)) {
685 kk = k-1;
686 } else {
687 kk = k+1;
688 }
689
690 // If any of these cells has zero volume we don't want to use this stencil
691 Real test_zero = cvol(i-1,jj,k ) * cvol(i-1,j ,kk) * cvol(i-1,jj,kk) *
692 cvol(i ,jj,k ) * cvol(i ,j ,kk) * cvol(i ,jj,kk);
693
694 if (test_zero == Real(0.))
695 {
696 int jalt = 2*j - jj;
697 int kalt = 2*k - kk;
698 Real test_zero_jalt = cvol(i-1,jalt,k) * cvol(i-1,j,kk ) * cvol(i-1,jalt,kk ) *
699 cvol(i ,jalt,k) * cvol(i ,j,kk ) * cvol(i ,jalt,kk );
700
701 Real test_zero_kalt = cvol(i-1,jj ,k) * cvol(i-1,j,kalt) * cvol(i-1,jj ,kalt) *
702 cvol(i ,jj ,k) * cvol(i ,j,kalt) * cvol(i ,jj ,kalt);
703
704 Real test_zero_jkalt = cvol(i-1,jalt,k) * cvol(i-1,j,kalt) * cvol(i-1,jalt,kalt) *
705 cvol(i ,jalt,k) * cvol(i ,j,kalt) * cvol(i ,jalt,kalt);
706
707 if (test_zero_jalt > Real(0.)) {
708 jj = jalt;
709 } else if (test_zero_kalt > Real(0.)) {
710 kk = kalt;
711 } else if (test_zero_jkalt > Real(0.)) {
712 jj = jalt;
713 kk = kalt;
714 }
715 }
716
717 Real d0 = Real(0.5) + ccent(i , j, k,0); // distance in x-dir from centroid(i , j, k) to face(i,j,k)
718 Real d1 = Real(0.5) - ccent(i-1, j, k,0); // distance in x-dir from centroid(i-1, j, k) to face(i,j,k)
719
720 Real d0_jj = Real(0.5) + ccent(i ,jj, k,0); // distance in x-dir from centroid(i ,jj, k) to face(i,jj,k)
721 Real d1_jj = Real(0.5) - ccent(i-1,jj, k,0); // distance in x-dir from centroid(i-1,jj, k) to face(i,jj,k)
722
723 Real d0_kk = Real(0.5) + ccent(i , j,kk,0); // distance in x-dir from centroid(i , j,kk) to face(i,j,kk)
724 Real d1_kk = Real(0.5) - ccent(i-1, j,kk,0); // distance in x-dir from centroid(i-1, j,kk) to face(i,j,kk)
725
726 Real d0_jj_kk = Real(0.5) + ccent(i ,jj,kk,0); // distance in x-dir from centroid(i ,jj,kk) to face(i,jj,kk)
727 Real d1_jj_kk = Real(0.5) - ccent(i-1,jj,kk,0); // distance in x-dir from centroid(i-1,jj,kk) to face(i,jj,kk)
728
729 Real a0 = d1 / (d0 + d1); // coefficient of phi(i , j, k,n)
730 Real a1 = Real(1.0) - a0; // coefficient of phi(i-1, j, k,n)
731
732 Real a0_jj = d1_jj / (d0_jj + d1_jj); // coefficient of phi(i ,jj, k,n)
733 Real a1_jj = Real(1.0) - a0_jj; // coefficient of phi(i-1,jj, k,n)
734
735 Real a0_kk = d1_kk / (d0_kk + d1_kk); // coefficient of phi(i , j,kk,n)
736 Real a1_kk = Real(1.0) - a0_kk; // coefficient of phi(i-1, j,kk,n)
737
738 Real a0_jj_kk = d1_jj_kk / (d0_jj_kk + d1_jj_kk); // coefficient of phi(i ,jj,kk,n)
739 Real a1_jj_kk = Real(1.0) - a0_jj_kk; // coefficient of phi(i-1,jj,kk,n)
740
741 Real phi01 = a0 * phi(i, j, k,n) + a1 * phi(i-1, j, k,n); // interpolation in x-dir of phi
742 Real y01 = a0 * ccent(i, j, k,1) + a1 * ccent(i-1, j, k,1); // interpolation in x-dir of y-loc
743 Real z01 = a0 * ccent(i, j, k,2) + a1 * ccent(i-1, j, k,2); // interpolation in x-dir of z-loc
744
745 Real phi01_jj = a0_jj * phi(i,jj, k,n) + a1_jj * phi(i-1,jj, k,n); // interpolation in x-dir of phi
746 Real y01_jj = a0_jj * ccent(i,jj, k,1) + a1_jj * ccent(i-1,jj, k,1); // interpolation in x-dir of y-loc
747 Real z01_jj = a0_jj * ccent(i,jj, k,2) + a1_jj * ccent(i-1,jj, k,2); // interpolation in x-dir of z-loc
748
749 Real phi01_kk = a0_kk * phi(i, j,kk,n) + a1_kk * phi(i-1, j,kk,n); // interpolation in x-dir of phi
750 Real y01_kk = a0_kk * ccent(i, j,kk,1) + a1_kk * ccent(i-1, j,kk,1); // interpolation in x-dir of y-loc
751 Real z01_kk = a0_kk * ccent(i, j,kk,2) + a1_kk * ccent(i-1, j,kk,2); // interpolation in x-dir of z-loc
752
753 Real phi01_jj_kk = a0_jj_kk * phi(i,jj,kk,n) + a1_jj_kk * phi(i-1,jj,kk,n); // interpolation in x-dir of phi
754 Real y01_jj_kk = a0_jj_kk * ccent(i,jj,kk,1) + a1_jj_kk * ccent(i-1,jj,kk,1); // interpolation in x-dir of y-loc
755 Real z01_jj_kk = a0_jj_kk * ccent(i,jj,kk,2) + a1_jj_kk * ccent(i-1,jj,kk,2); // interpolation in x-dir of z-loc
756
757 // 2D interpolation on x-face from interpolated points at (y01,z01), (y01_jj,z01_jj), (y01_kk,z01_kk), (y01_jj_kk,z01_jj_kk)
758 // to centroid of x-face(i,j,k) at (fcx(i,j,k,0), fcx(i,j,k,1))
759
760 // We order these in order to form a set of four points on the x-face ...
761 // (x0,y0) : lower left
762 // (x1,y1) : lower right
763 // (x2,y2) : upper right
764 // (x3,y3) : upper left
765
766 Real y,z;
767
768 // This is the location of the face centroid relative to the central node
769 // Recall fcx holds (y,z) of the x-face centroid as components ( /0/1)
770 if (j < jj) {
771 y = Real(-0.5) + fcx(i,j,k,0); // (j,k) is in lower half of stencil so y < 0
772 } else {
773 y = Real(0.5) + fcx(i,j,k,0); // (j,k) is in upper half of stencil so y > 0
774 }
775
776 if (k < kk) {
777 z = Real(-0.5) + fcx(i,j,k,1); // (j,k) is in lower half of stencil so z < 0
778 } else {
779 z = Real(0.5) + fcx(i,j,k,1); // (j,k) is in upper half of stencil so z > 0
780 }
781
782 if (j < jj && k < kk) // (j,k) is lower left, (j+1,k+1) is upper right
783 {
784 Real y_ll = Real(-0.5)+y01;
785 Real z_ll = Real(-0.5)+z01;
786 Real y_hl = Real(0.5)+y01_jj;
787 Real z_hl = Real(-0.5)+z01_jj;
788 Real y_hh = Real(0.5)+y01_jj_kk;
789 Real z_hh = Real(0.5)+z01_jj_kk;
790 Real y_lh = Real(-0.5)+y01_kk;
791 Real z_lh = Real(0.5)+z01_kk;
792 phi_x(i,j,k,n) = EB_interp_in_quad(y,z, phi01, phi01_jj, phi01_jj_kk, phi01_kk,
793 y_ll, z_ll, y_hl, z_hl, y_hh, z_hh, y_lh, z_lh);
794 }
795 else if (jj < j && k < kk) // (j-1,k) is lower left, (j,k+1) is upper right
796 {
797 Real y_ll = Real(-0.5)+y01_jj;
798 Real z_ll = Real(-0.5)+z01_jj;
799 Real y_hl = Real(0.5)+y01;
800 Real z_hl = Real(-0.5)+z01;
801 Real y_hh = Real(0.5)+y01_kk;
802 Real z_hh = Real(0.5)+z01_kk;
803 Real y_lh = Real(-0.5)+y01_jj_kk;
804 Real z_lh = Real(0.5)+z01_jj_kk;
805 phi_x(i,j,k,n) = EB_interp_in_quad(y,z,phi01_jj,phi01,phi01_kk, phi01_jj_kk,
806 y_ll, z_ll, y_hl, z_hl, y_hh, z_hh, y_lh, z_lh);
807 }
808 else if (j < jj && kk < k) // (j,k-1) is lower left, (j+1,k) is upper right
809 {
810 Real y_ll = Real(-0.5)+y01_kk;
811 Real z_ll = Real(-0.5)+z01_kk;
812 Real y_hl = Real(0.5)+y01_jj_kk;
813 Real z_hl = Real(-0.5)+z01_jj_kk;
814 Real y_hh = Real(0.5)+y01_jj;
815 Real z_hh = Real(0.5)+z01_jj;
816 Real y_lh = Real(-0.5)+y01;
817 Real z_lh = Real(0.5)+z01;
818 phi_x(i,j,k,n) = EB_interp_in_quad(y,z,phi01_kk,phi01_jj_kk,phi01_jj, phi01,
819 y_ll, z_ll, y_hl, z_hl, y_hh, z_hh, y_lh, z_lh);
820 }
821 else if (jj < j && kk < k) // (j-1,k-1) is lower left, (j,k) is upper right
822 {
823 Real y_ll = Real(-0.5)+y01_jj_kk;
824 Real z_ll = Real(-0.5)+z01_jj_kk;
825 Real y_hl = Real(0.5)+y01_kk;
826 Real z_hl = Real(-0.5)+z01_kk;
827 Real y_hh = Real(0.5)+y01;
828 Real z_hh = Real(0.5)+z01;
829 Real y_lh = Real(-0.5)+y01_jj;
830 Real z_lh = Real(0.5)+z01_jj;
831 phi_x(i,j,k,n) = EB_interp_in_quad(y,z,phi01_jj_kk,phi01_kk,phi01, phi01_jj,
832 y_ll, z_ll, y_hl, z_hl, y_hh, z_hh, y_lh, z_lh);
833 }
834 else
835 {
836 amrex::Abort("Bad option in interpolation from cell centroid to x-face centroid!");
837 }
838 }
839 });
840}
841
842AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
843void eb_interp_centroid2facecent_y (Box const& vbx,
844 Array4<Real const> const& phi,
845 Array4<Real const> const& apy,
846 Array4<Real const> const& cvol,
847 Array4<Real const> const& ccent,
848 Array4<Real const> const& fcy,
849 Array4<Real> const& phi_y,
850 int ncomp,
851 const Box& domain,
852 const BCRec* bc) noexcept
853{
854 const Dim3 domlo = amrex::lbound(domain);
855 const Dim3 domhi = amrex::ubound(domain);
856 // Note that ccent holds (x,y,z) of the cell centroids as components (0/1/2)
857 // fcy holds (x, z) of the y-face centroid as components (0/ /1)
858 amrex::Loop(vbx, ncomp, [=] (int i, int j, int k, int n) noexcept
859 {
860 if (apy(i,j,k) == 0)
861 {
862#ifdef AMREX_USE_FLOAT
863 phi_y(i,j,k,n) = Real(1e20);
864#else
865 phi_y(i,j,k,n) = 1e40;
866#endif
867 }
868 else if ( (j == domlo.y) && (bc[n].lo(1) == BCType::ext_dir) )
869 {
870 phi_y(i,j,k,n) = phi(i,domlo.y-1,k,n);
871 }
872 else if ( (j == domhi.y+1) && (bc[n].hi(1) == BCType::ext_dir) )
873 {
874 phi_y(i,j,k,n) = phi(i,domhi.y+1,k,n);
875 }
876 else if (apy(i,j,k) == Real(1.) && cvol(i,j,k) == Real(1.) && cvol(i,j-1,k) == Real(1.))
877 {
878 phi_y(i,j,k,n) = Real(0.5) * (phi(i,j ,k,n) + phi(i,j-1,k,n));
879 }
880 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)) &&
881 (std::abs(ccent(i,j,k,2)-ccent(i,j-1,k,2)) < Real(1.e-8)) )
882 {
883 Real d0 = Real(0.5) + ccent(i,j ,k,1); // distance in y-dir from centroid(i,j ,k) to face(i,j,k)
884 Real d1 = Real(0.5) - ccent(i,j-1,k,1); // distance in y-dir from centroid(i,j-1,k) to face(i,j,k)
885 Real a0 = d1 / (d0 + d1); // coefficient of phi(i,j ,k,n)
886 Real a1 = d0 / (d0 + d1); // coefficient of phi(i,j-1,k,n)
887 phi_y(i,j,k,n) = a0*phi(i,j,k,n) + a1*phi(i,j-1,k,n); // interpolation of phi in y-direction
888 }
889 else
890 {
891
892 // We must add additional tests to avoid the case where fcy is very close to zero, but i-1/i+1 or k-1/k+1
893 // might be covered cells -- this can happen when the EB is exactly aligned with the grid planes
894 int ii,kk;
895 if (std::abs(fcy(i,j,k,0)) > Real(1.e-8)) {
896 ii = (fcy(i,j,k,0) < Real(0.0)) ? i - 1 : i + 1;
897 } else if (apy(i-1,j,k) > Real(0.)) {
898 ii = i-1;
899 } else {
900 ii = i+1;
901 }
902
903 if (std::abs(fcy(i,j,k,1)) > Real(1.e-8)) {
904 kk = (fcy(i,j,k,1) < Real(0.0)) ? k - 1 : k + 1;
905 } else if (apy(i,j,k-1) > Real(0.)) {
906 kk = k-1;
907 } else {
908 kk = k+1;
909 }
910
911 // If any of these cells has zero volume we don't want to use this stencil
912 Real test_zero = cvol(ii,j-1,k) * cvol(i,j-1,kk) * cvol(ii,j-1,kk) *
913 cvol(ii,j ,k) * cvol(i,j ,kk) * cvol(ii,j ,kk);
914
915 if (test_zero == Real(0.))
916 {
917 int ialt = 2*i - ii;
918 int kalt = 2*k - kk;
919 Real test_zero_ialt = cvol(ialt,j-1,k) * cvol(i,j-1,kk ) * cvol(ialt,j-1,kk ) *
920 cvol(ialt,j ,k) * cvol(i,j ,kk ) * cvol(ialt,j ,kk );
921
922 Real test_zero_kalt = cvol(ii ,j-1,k) * cvol(i,j-1,kalt) * cvol(ii ,j-1,kalt) *
923 cvol(ii ,j ,k) * cvol(i,j ,kalt) * cvol(ii ,j ,kalt);
924
925 Real test_zero_ikalt = cvol(ialt,j-1,k) * cvol(i,j-1,kalt) * cvol(ialt,j-1,kalt) *
926 cvol(ialt,j ,k) * cvol(i,j ,kalt) * cvol(ialt,j ,kalt);
927
928 if (test_zero_ialt > Real(0.)) {
929 ii = ialt;
930 } else if (test_zero_kalt > Real(0.)) {
931 kk = kalt;
932 } else if (test_zero_ikalt > Real(0.)) {
933 ii = ialt;
934 kk = kalt;
935 }
936 }
937
938 Real d0 = Real(0.5) + ccent( i,j , k,1); // distance in y-dir from centroid( i,j , k) to face(i,j,k)
939 Real d1 = Real(0.5) - ccent( i,j-1, k,1); // distance in y-dir from centroid( i,j-1, k) to face(i,j,k)
940
941 Real d0_ii = Real(0.5) + ccent(ii,j , k,1); // distance in y-dir from centroid(ii,j , k) to face(ii,j,k)
942 Real d1_ii = Real(0.5) - ccent(ii,j-1, k,1); // distance in y-dir from centroid(ii,j-1, k) to face(ii,j,k)
943
944 Real d0_kk = Real(0.5) + ccent( i,j ,kk,1); // distance in y-dir from centroid( i,j ,kk) to face( i,j,kk)
945 Real d1_kk = Real(0.5) - ccent( i,j-1,kk,1); // distance in y-dir from centroid( i,j-1,kk) to face( i,j,kk)
946
947 Real d0_ii_kk = Real(0.5) + ccent(ii,j ,kk,1); // distance in y-dir from centroid(ii,j ,kk) to face(ii,j,kk)
948 Real d1_ii_kk = Real(0.5) - ccent(ii,j-1,kk,1); // distance in y-dir from centroid(ii,j-1,kk) to face(ii,j,kk)
949
950 Real a0 = d1 / (d0 + d1); // coefficient of phi( i,j ,k,n)
951 Real a1 = d0 / (d0 + d1); // coefficient of phi( i,j-1,k,n)
952
953 Real a0_ii = d1_ii / (d0_ii + d1_ii); // coefficient of phi(ii,j ,k,n)
954 Real a1_ii = d0_ii / (d0_ii + d1_ii); // coefficient of phi(ii,j-1,k,n)
955
956 Real a0_kk = d1_kk / (d0_kk + d1_kk); // coefficient of phi( i,j ,kk,n)
957 Real a1_kk = d0_kk / (d0_kk + d1_kk); // coefficient of phi( i,j-1,kk,n)
958
959 Real a0_ii_kk = d1_ii_kk / (d0_ii_kk + d1_ii_kk); // coefficient of phi(ii,j ,kk,n)
960 Real a1_ii_kk = d0_ii_kk / (d0_ii_kk + d1_ii_kk); // coefficient of phi(ii,j-1,kk,n)
961
962 Real phi01 = a0 * phi( i, j, k,n) + a1 * phi( i, j-1, k,n); // interpolation in y-dir of phi
963 Real x01 = a0 * ccent( i, j, k,0) + a1 * ccent( i, j-1, k,0); // interpolation in y-dir of x-loc
964 Real z01 = a0 * ccent( i, j, k,2) + a1 * ccent( i, j-1, k,2); // interpolation in y-dir of z-loc
965
966 Real phi01_ii = a0_ii * phi(ii, j, k,n) + a1_ii * phi(ii, j-1, k,n); // interpolation in y-dir of phi
967 Real x01_ii = a0_ii * ccent(ii, j, k,0) + a1_ii * ccent(ii, j-1, k,0); // interpolation in y-dir of x-loc
968 Real z01_ii = a0_ii * ccent(ii, j, k,2) + a1_ii * ccent(ii, j-1, k,2); // interpolation in y-dir of z-loc
969
970 Real phi01_kk = a0_kk * phi( i, j,kk,n) + a1_kk * phi( i, j-1,kk,n); // interpolation in y-dir of phi
971 Real x01_kk = a0_kk * ccent( i, j,kk,0) + a1_kk * ccent( i, j-1,kk,0); // interpolation in y-dir of x-loc
972 Real z01_kk = a0_kk * ccent( i, j,kk,2) + a1_kk * ccent( i, j-1,kk,2); // interpolation in y-dir of z-loc
973
974 Real phi01_ii_kk = a0_ii_kk * phi(ii, j,kk,n) + a1_ii_kk * phi(ii, j-1,kk,n); // interpolation in y-dir of phi
975 Real x01_ii_kk = a0_ii_kk * ccent(ii, j,kk,0) + a1_ii_kk * ccent(ii, j-1,kk,0); // interpolation in y-dir of x-loc
976 Real z01_ii_kk = a0_ii_kk * ccent(ii, j,kk,2) + a1_ii_kk * ccent(ii, j-1,kk,2); // interpolation in y-dir of z-loc
977
978 // We order these in order to form a set of four points on the x-face ...
979 // (x0,y0) : lower left
980 // (x1,y1) : lower right
981 // (x2,y2) : upper right
982 // (x3,y3) : upper left
983
984 Real x,z;
985
986 // This is the location of the face centroid relative to the central node
987 // Recall fcy holds (x,z) of the x-face centroid as components (0/ /1)
988 if (i < ii) {
989 x = Real(-0.5) + fcy(i,j,k,0); // (i,k) is in lower half of stencil so x < 0
990 } else {
991 x = Real(0.5) + fcy(i,j,k,0); // (i,k) is in upper half of stencil so x > 0
992 }
993
994 if (k < kk) {
995 z = Real(-0.5) + fcy(i,j,k,1); // (i,k) is in lower half of stencil so z < 0
996 } else {
997 z = Real(0.5) + fcy(i,j,k,1); // (i,k) is in upper half of stencil so z > 0
998 }
999
1000 if (i < ii && k < kk) // (i,k) is lower left, (i+1,k+1) is upper right
1001 {
1002 Real x_ll = Real(-0.5)+x01;
1003 Real z_ll = Real(-0.5)+z01;
1004 Real x_hl = Real(0.5)+x01_ii;
1005 Real z_hl = Real(-0.5)+z01_ii;
1006 Real x_hh = Real(0.5)+x01_ii_kk;
1007 Real z_hh = Real(0.5)+z01_ii_kk;
1008 Real x_lh = Real(-0.5)+x01_kk;
1009 Real z_lh = Real(0.5)+z01_kk;
1010 phi_y(i,j,k,n) = EB_interp_in_quad(x,z, phi01, phi01_ii, phi01_ii_kk, phi01_kk,
1011 x_ll, z_ll, x_hl, z_hl, x_hh, z_hh, x_lh, z_lh);
1012 }
1013 else if (ii < i && k < kk) // (i-1,k) is lower left, (i,k+1) is upper right
1014 {
1015 Real x_ll = Real(-0.5)+x01_ii;
1016 Real z_ll = Real(-0.5)+z01_ii;
1017 Real x_hl = Real(0.5)+x01;
1018 Real z_hl = Real(-0.5)+z01;
1019 Real x_hh = Real(0.5)+x01_kk;
1020 Real z_hh = Real(0.5)+z01_kk;
1021 Real x_lh = Real(-0.5)+x01_ii_kk;
1022 Real z_lh = Real(0.5)+z01_ii_kk;
1023 phi_y(i,j,k,n) = EB_interp_in_quad(x,z,phi01_ii,phi01,phi01_kk, phi01_ii_kk,
1024 x_ll, z_ll, x_hl, z_hl, x_hh, z_hh, x_lh, z_lh);
1025 }
1026 else if (i < ii && kk < k) // (i,k-1) is lower left, (i+1,k) is upper right
1027 {
1028 Real x_ll = Real(-0.5)+x01_kk;
1029 Real z_ll = Real(-0.5)+z01_kk;
1030 Real x_hl = Real(0.5)+x01_ii_kk;
1031 Real z_hl = Real(-0.5)+z01_ii_kk;
1032 Real x_hh = Real(0.5)+x01_ii;
1033 Real z_hh = Real(0.5)+z01_ii;
1034 Real x_lh = Real(-0.5)+x01;
1035 Real z_lh = Real(0.5)+z01;
1036 phi_y(i,j,k,n) = EB_interp_in_quad(x,z,phi01_kk,phi01_ii_kk,phi01_ii, phi01,
1037 x_ll, z_ll, x_hl, z_hl, x_hh, z_hh, x_lh, z_lh);
1038 }
1039 else if (ii < i && kk < k) // (i-1,k-1) is lower left, (i,k) is upper right
1040 {
1041 Real x_ll = Real(-0.5)+x01_ii_kk;
1042 Real z_ll = Real(-0.5)+z01_ii_kk;
1043 Real x_hl = Real(0.5)+x01_kk;
1044 Real z_hl = Real(-0.5)+z01_kk;
1045 Real x_hh = Real(0.5)+x01;
1046 Real z_hh = Real(0.5)+z01;
1047 Real x_lh = Real(-0.5)+x01_ii;
1048 Real z_lh = Real(0.5)+z01_ii;
1049 phi_y(i,j,k,n) = EB_interp_in_quad(x,z,phi01_ii_kk,phi01_kk,phi01, phi01_ii,
1050 x_ll, z_ll, x_hl, z_hl, x_hh, z_hh, x_lh, z_lh);
1051 }
1052 else
1053 {
1054 amrex::Abort("Bad option in interpolation from cell centroid to y-face centroid!");
1055 }
1056 }
1057 });
1058}
1059
1060AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1061void eb_interp_centroid2facecent_z (Box const& wbx,
1062 Array4<Real const> const& phi,
1063 Array4<Real const> const& apz,
1064 Array4<Real const> const& cvol,
1065 Array4<Real const> const& ccent,
1066 Array4<Real const> const& fcz,
1067 Array4<Real> const& phi_z,
1068 int ncomp,
1069 const Box& domain,
1070 const BCRec* bc) noexcept
1071{
1072 const Dim3 domlo = amrex::lbound(domain);
1073 const Dim3 domhi = amrex::ubound(domain);
1074 // Note that ccent holds (x,y,z) of the cell centroids as components (0/1/2)
1075 // fcz holds (x,y ) of the z-face centroid as components (0/1 )
1076 amrex::Loop(wbx, ncomp, [=] (int i, int j, int k, int n) noexcept
1077 {
1078 if (apz(i,j,k) == 0)
1079 {
1080#ifdef AMREX_USE_FLOAT
1081 phi_z(i,j,k,n) = Real(1e20);
1082#else
1083 phi_z(i,j,k,n) = 1e40;
1084#endif
1085 }
1086 else if ( (k == domlo.z) && (bc[n].lo(2) == BCType::ext_dir) )
1087 {
1088 phi_z(i,j,k,n) = phi(i,j,domlo.z-1,n);
1089 }
1090 else if ( (k == domhi.z+1) && (bc[n].hi(2) == BCType::ext_dir) )
1091 {
1092 phi_z(i,j,k,n) = phi(i,j,domhi.z+1,n);
1093 }
1094 else if (apz(i,j,k) == Real(1.) && cvol(i,j,k) == Real(1.) && cvol(i,j,k-1) == Real(1.))
1095 {
1096 phi_z(i,j,k,n) = Real(0.5) * (phi(i,j,k,n) + phi(i,j,k-1,n));
1097 }
1098 else if ( apz(i,j,k) == Real(1.) && (std::abs(ccent(i,j,k,0)-ccent(i,j,k-1,0)) < Real(1.e-8)) &&
1099 (std::abs(ccent(i,j,k,1)-ccent(i,j,k-1,1)) < Real(1.e-8)) )
1100 {
1101 Real d0 = Real(0.5) + ccent(i,j,k ,2); // distance in z-dir from centroid(i,j,k ) to face(i,j,k)
1102 Real d1 = Real(0.5) - ccent(i,j,k-1,2); // distance in z-dir from centroid(i,j,k-1) to face(i,j,k)
1103 Real a0 = d1 / (d0 + d1); // coefficient of phi(i,j,k ,n)
1104 Real a1 = d0 / (d0 + d1); // coefficient of phi(i,j,k-1,n)
1105 phi_z(i,j,k,n) = a0*phi(i,j,k,n) + a1*phi(i,j,k-1,n); // interpolation of phi in z-direction
1106 }
1107 else
1108 {
1109 // We must add additional tests to avoid the case where fcz is very close to zero, but i-1/i+1 or j-1/j+1
1110 // might be covered cells -- this can happen when the EB is exactly aligned with the grid planes
1111 int ii,jj;
1112 if (std::abs(fcz(i,j,k,0)) > Real(1.e-8)) {
1113 ii = (fcz(i,j,k,0) < Real(0.0)) ? i - 1 : i + 1;
1114 } else if (apz(i-1,j,k) > Real(0.)) {
1115 ii = i-1;
1116 } else {
1117 ii = i+1;
1118 }
1119
1120 if (std::abs(fcz(i,j,k,1)) > Real(1.e-8)) {
1121 jj = (fcz(i,j,k,1) < Real(0.0)) ? j - 1 : j + 1;
1122 } else if (apz(i,j-1,k) > Real(0.)) {
1123 jj = j-1;
1124 } else {
1125 jj = j+1;
1126 }
1127
1128 // If any of these cells has zero volume we don't want to use this stencil
1129 Real test_zero = cvol(ii,j,k-1) * cvol(i,jj,k-1) * cvol(ii,jj,k-1) *
1130 cvol(ii,j,k ) * cvol(i,jj,k ) * cvol(ii,jj,k );
1131
1132 if (test_zero == Real(0.))
1133 {
1134 int ialt = 2*i - ii;
1135 int jalt = 2*j - jj;
1136 Real test_zero_ialt = cvol(ialt,j,k-1) * cvol(i,jj ,k-1) * cvol(ialt,jj ,k-1) *
1137 cvol(ialt,j,k ) * cvol(i,jj ,k ) * cvol(ialt,jj ,k );
1138
1139 Real test_zero_jalt = cvol(ii ,j,k-1) * cvol(i,jalt,k-1) * cvol(ii ,jalt,k-1) *
1140 cvol(ii ,j,k ) * cvol(i,jalt,k ) * cvol(ii ,jalt,k );
1141
1142 Real test_zero_ijalt = cvol(ialt,j,k-1) * cvol(i,jalt,k-1) * cvol(ialt,jalt,k-1) *
1143 cvol(ialt,j,k ) * cvol(i,jalt,k ) * cvol(ialt,jalt,k );
1144
1145 if (test_zero_ialt > Real(0.)) {
1146 ii = ialt;
1147 } else if (test_zero_jalt > Real(0.)) {
1148 jj = jalt;
1149 } else if (test_zero_ijalt > Real(0.)) {
1150 ii = ialt;
1151 jj = jalt;
1152 }
1153 }
1154
1155 Real d0 = Real(0.5) + ccent( i, j,k ,2); // distance in z-dir from centroid(i,j,k ) to face(i,j,k)
1156 Real d1 = Real(0.5) - ccent( i, j,k-1,2); // distance in z-dir from centroid(i,j,k-1) to face(i,j,k)
1157
1158 Real d0_ii = Real(0.5) + ccent(ii, j,k ,2); // distance in z-dir from centroid(ii,j,k ) to face(ii,j,k)
1159 Real d1_ii = Real(0.5) - ccent(ii, j,k-1,2); // distance in z-dir from centroid(ii,j,k-1) to face(ii,j,k)
1160
1161 Real d0_jj = Real(0.5) + ccent( i,jj,k ,2); // distance in z-dir from centroid(ii,j,k ) to face(ii,j,k)
1162 Real d1_jj = Real(0.5) - ccent( i,jj,k-1,2); // distance in z-dir from centroid(ii,j,k-1) to face(ii,j,k)
1163
1164 Real d0_ii_jj = Real(0.5) + ccent(ii,jj,k ,2); // distance in z-dir from centroid(ii,j,k ) to face(ii,j,k)
1165 Real d1_ii_jj = Real(0.5) - ccent(ii,jj,k-1,2); // distance in z-dir from centroid(ii,j,k-1) to face(ii,j,k)
1166
1167 Real a0 = d1 / (d0 + d1); // coefficient of phi( i, j,k ,n)
1168 Real a1 = d0 / (d0 + d1); // coefficient of phi( i, j,k-1,n)
1169
1170 Real a0_ii = d1_ii / (d0_ii + d1_ii); // coefficient of phi(ii, j,k ,n)
1171 Real a1_ii = d0_ii / (d0_ii + d1_ii); // coefficient of phi(ii, j,k-1,n)
1172
1173 Real a0_jj = d1_jj / (d0_jj + d1_jj); // coefficient of phi( i,jj,k ,n)
1174 Real a1_jj = d0_jj / (d0_jj + d1_jj); // coefficient of phi( i,jj,k-1,n)
1175
1176 Real a0_ii_jj = d1_ii_jj / (d0_ii_jj + d1_ii_jj); // coefficient of phi(ii,jj,k ,n)
1177 Real a1_ii_jj = d0_ii_jj / (d0_ii_jj + d1_ii_jj); // coefficient of phi(ii,jj,k-1,n)
1178
1179 Real phi01 = a0 * phi( i, j,k,n) + a1 * phi( i, j,k-1,n); // interpolation of phi in z-direction
1180 Real x01 = a0 * ccent( i, j,k,0) + a1 * ccent( i, j,k-1,0); // interpolation of x-loc in z-direction
1181 Real y01 = a0 * ccent( i, j,k,1) + a1 * ccent( i, j,k-1,1); // interpolation of y-loc in z-direction
1182
1183 Real phi01_ii = a0_ii * phi(ii, j,k,n) + a1_ii * phi(ii, j,k-1,n); // interpolation of phi in z-direction
1184 Real x01_ii = a0_ii * ccent(ii, j,k,0) + a1_ii * ccent(ii, j,k-1,0); // interpolation of x-loc in z-direction
1185 Real y01_ii = a0_ii * ccent(ii, j,k,1) + a1_ii * ccent(ii, j,k-1,1); // interpolation of y-loc in z-direction
1186
1187 Real phi01_jj = a0_jj * phi( i,jj,k,n) + a1_jj * phi( i,jj,k-1,n); // interpolation of phi in z-direction
1188 Real x01_jj = a0_jj * ccent( i,jj,k,0) + a1_jj * ccent( i,jj,k-1,0); // interpolation of x-loc in z-direction
1189 Real y01_jj = a0_jj * ccent( i,jj,k,1) + a1_jj * ccent( i,jj,k-1,1); // interpolation of y-loc in z-direction
1190
1191 Real phi01_ii_jj = a0_ii_jj * phi(ii,jj,k,n) + a1_ii_jj * phi(ii,jj,k-1,n); // interpolation of phi in z-direction
1192 Real x01_ii_jj = a0_ii_jj * ccent(ii,jj,k,0) + a1_ii_jj * ccent(ii,jj,k-1,0); // interpolation of x-loc in z-direction
1193 Real y01_ii_jj = a0_ii_jj * ccent(ii,jj,k,1) + a1_ii_jj * ccent(ii,jj,k-1,1); // interpolation of y-loc in z-direction
1194
1195 // We order these in order to form a set of four points on the x-face ...
1196 // (x0,y0) : lower left
1197 // (x1,y1) : lower right
1198 // (x2,y2) : upper right
1199 // (x3,y3) : upper left
1200
1201 Real x,y;
1202
1203 // This is the location of the face centroid relative to the central node
1204 // Recall fcz holds (x,y) of the x-face centroid as components (0/1/ )
1205 if (i < ii) {
1206 x = Real(-0.5) + fcz(i,j,k,0); // (i,k) is in lower half of stencil so x < 0
1207 } else {
1208 x = Real(0.5) + fcz(i,j,k,0); // (i,k) is in upper half of stencil so x > 0
1209 }
1210
1211 if (j < jj) {
1212 y = Real(-0.5) + fcz(i,j,k,1); // (i,k) is in lower half of stencil so z < 0
1213 } else {
1214 y = Real(0.5) + fcz(i,j,k,1); // (i,k) is in upper half of stencil so z > 0
1215 }
1216
1217 if (i < ii && j < jj) // (i,j) is lower left, (i+1,j+1) is upper right
1218 {
1219 Real x_ll = Real(-0.5)+x01;
1220 Real y_ll = Real(-0.5)+y01;
1221 Real x_hl = Real(0.5)+x01_ii;
1222 Real y_hl = Real(-0.5)+y01_ii;
1223 Real x_hh = Real(0.5)+x01_ii_jj;
1224 Real y_hh = Real(0.5)+y01_ii_jj;
1225 Real x_lh = Real(-0.5)+x01_jj;
1226 Real y_lh = Real(0.5)+y01_jj;
1227 phi_z(i,j,k,n) = EB_interp_in_quad(x,y, phi01, phi01_ii, phi01_ii_jj, phi01_jj,
1228 x_ll, y_ll, x_hl, y_hl, x_hh, y_hh, x_lh, y_lh);
1229 }
1230 else if (ii < i && j < jj) // (i-1,j) is lower left, (i,j+1) is upper right
1231 {
1232 Real x_ll = Real(-0.5)+x01_ii;
1233 Real y_ll = Real(-0.5)+y01_ii;
1234 Real x_hl = Real(0.5)+x01;
1235 Real y_hl = Real(-0.5)+y01;
1236 Real x_hh = Real(0.5)+x01_jj;
1237 Real y_hh = Real(0.5)+y01_jj;
1238 Real x_lh = Real(-0.5)+x01_ii_jj;
1239 Real y_lh = Real(0.5)+y01_ii_jj;
1240 phi_z(i,j,k,n) = EB_interp_in_quad(x,y,phi01_ii,phi01,phi01_jj, phi01_ii_jj,
1241 x_ll, y_ll, x_hl, y_hl, x_hh, y_hh, x_lh, y_lh);
1242 }
1243 else if (i < ii && jj < j) // (i,j-1) is lower left, (i+1,j) is upper right
1244 {
1245 Real x_ll = Real(-0.5)+x01_jj;
1246 Real y_ll = Real(-0.5)+y01_jj;
1247 Real x_hl = Real(0.5)+x01_ii_jj;
1248 Real y_hl = Real(-0.5)+y01_ii_jj;
1249 Real x_hh = Real(0.5)+x01_ii;
1250 Real y_hh = Real(0.5)+y01_ii;
1251 Real x_lh = Real(-0.5)+x01;
1252 Real y_lh = Real(0.5)+y01;
1253 phi_z(i,j,k,n) = EB_interp_in_quad(x,y,phi01_jj,phi01_ii_jj,phi01_ii, phi01,
1254 x_ll, y_ll, x_hl, y_hl, x_hh, y_hh, x_lh, y_lh);
1255 }
1256 else if (ii < i && jj < j) // (i-1,j-1) is lower left, (i,j) is upper right
1257 {
1258 Real x_ll = Real(-0.5)+x01_ii_jj;
1259 Real y_ll = Real(-0.5)+y01_ii_jj;
1260 Real x_hl = Real(0.5)+x01_jj;
1261 Real y_hl = Real(-0.5)+y01_jj;
1262 Real x_hh = Real(0.5)+x01;
1263 Real y_hh = Real(0.5)+y01;
1264 Real x_lh = Real(-0.5)+x01_ii;
1265 Real y_lh = Real(0.5)+y01_ii;
1266 phi_z(i,j,k,n) = EB_interp_in_quad(x,y,phi01_ii_jj,phi01_jj,phi01, phi01_ii,
1267 x_ll, y_ll, x_hl, y_hl, x_hh, y_hh, x_lh, y_lh);
1268 }
1269 else
1270 {
1271 amrex::Abort("Bad option in interpolation from cell centroid to z-face centroid!");
1272 }
1273 }
1274 });
1275}
1276
1277AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1278void eb_interp_cc2face_x (Box const& ubx,
1279 Array4<Real const> const& phi,
1280 Array4<Real> const& edg_x,
1281 int ncomp,
1282 const Box& domain,
1283 const BCRec* bc) noexcept
1284{
1285 const Dim3 domlo = amrex::lbound(domain);
1286 const Dim3 domhi = amrex::ubound(domain);
1287 amrex::Loop(ubx, ncomp, [=] (int i, int j, int k, int n) noexcept
1288 {
1289 if ( (i == domlo.x) && (bc[n].lo(0) == BCType::ext_dir) )
1290 {
1291 edg_x(i,j,k,n) = phi(domlo.x-1,j,k,n);
1292 }
1293 else if ( (i == domhi.x+1) && (bc[n].hi(0) == BCType::ext_dir) )
1294 {
1295 edg_x(i,j,k,n) = phi(domhi.x+1,j,k,n);
1296 }
1297 else
1298 {
1299 edg_x(i,j,k,n) = Real(0.5) * ( phi(i,j,k,n) + phi(i-1,j,k,n) );
1300 }
1301 });
1302}
1303
1304AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1305void eb_interp_cc2face_y (Box const& vbx,
1306 Array4<Real const> const& phi,
1307 Array4<Real> const& edg_y,
1308 int ncomp,
1309 const Box& domain,
1310 const BCRec* bc) noexcept
1311{
1312 const Dim3 domlo = amrex::lbound(domain);
1313 const Dim3 domhi = amrex::ubound(domain);
1314 amrex::Loop(vbx, ncomp, [=] (int i, int j, int k, int n) noexcept
1315 {
1316 if ( (j == domlo.y) && (bc[n].lo(1) == BCType::ext_dir) )
1317 {
1318 edg_y(i,j,k,n) = phi(i,domlo.y-1,k,n);
1319 }
1320 else if ( (j == domhi.y+1) && (bc[n].hi(1) == BCType::ext_dir) )
1321 {
1322 edg_y(i,j,k,n) = phi(i,domhi.y+1,k,n);
1323 }
1324 else
1325 {
1326 edg_y(i,j,k,n) = Real(0.5) * (phi(i,j ,k,n) + phi(i,j-1,k,n));
1327 }
1328 });
1329}
1330
1331AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1332void eb_interp_cc2face_z (Box const& wbx,
1333 Array4<Real const> const& phi,
1334 Array4<Real> const& edg_z,
1335 int ncomp,
1336 const Box& domain,
1337 const BCRec* bc) noexcept
1338{
1339 const Dim3 domlo = amrex::lbound(domain);
1340 const Dim3 domhi = amrex::ubound(domain);
1341 amrex::Loop(wbx, ncomp, [=] (int i, int j, int k, int n) noexcept
1342 {
1343 if ( (k == domlo.z) && (bc[n].lo(2) == BCType::ext_dir) )
1344 {
1345 edg_z(i,j,k,n) = phi(i,j,domlo.z-1,n);
1346 }
1347 else if ( (k == domhi.z+1) && (bc[n].hi(2) == BCType::ext_dir) )
1348 {
1349 edg_z(i,j,k,n) = phi(i,j,domhi.z+1,n);
1350 }
1351 else
1352 {
1353 edg_z(i,j,k,n) = Real(0.5) * (phi(i,j ,k,n) + phi(i,j,k-1,n));
1354 }
1355 });
1356}
1357
1358AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1359void eb_add_divergence_from_flow (int i, int j, int k, int n, Array4<Real> const& divu,
1360 Array4<Real const> const& vel_eb, Array4<EBCellFlag const> const& flag,
1361 Array4<Real const> const& vfrc, Array4<Real const> const& bnorm,
1362 Array4<Real const> const& barea, GpuArray<Real,3> const& dxinv)
1363{
1364 if (flag(i,j,k).isSingleValued())
1365 {
1366 Real Ueb_dot_n = vel_eb(i,j,k,0)*bnorm(i,j,k,0)
1367 + vel_eb(i,j,k,1)*bnorm(i,j,k,1)
1368 + vel_eb(i,j,k,2)*bnorm(i,j,k,2);
1369
1370 divu(i,j,k,n) += Ueb_dot_n * barea(i,j,k) * dxinv[0] / vfrc(i,j,k);
1371 }
1372}
1373
1374
1375}
1376#endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMREX_RESTRICT
#define AMREX_RESTRICT
Definition AMReX_Extension.H:37
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
fine
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
crse
Array4< Real const > crse
Definition AMReX_InterpFaceRegister.cpp:92
AMReX_REAL.H
if
if(!(yy_init))
Definition amrex_iparser.lex.nolint.H:935
amrex::BCRec
Boundary Condition Records. Necessary information and functions for computing boundary conditions.
Definition AMReX_BCRec.H:17
amrex
Definition AMReX_Amr.cpp:49
amrex::eb_avgdown_face_y
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::eb_interp_centroid2facecent_y
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::Box
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
amrex::eb_add_divergence_from_flow
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::eb_avgdown_face_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_avgdown_face_z(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_3D_C.H:177
amrex::eb_compute_divergence
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::eb_avgdown_face_x
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::eb_avg_fc_to_cc
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::eb_avgdown_boundaries
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::eb_avgdown
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::eb_avgdown_with_vol
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::ubound
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:315
amrex::lbound
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
amrex::eb_set_covered_nodes
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::eb_interp_cc2facecent_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_cc2facecent_z(Box const &wbx, Array4< Real const > const &phi, Array4< Real const > const &apz, Array4< Real const > const &fcz, Array4< Real > const &edg_z, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_3D_C.H:556
amrex::eb_interp_cc2facecent_y
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::eb_interp_cc2face_x
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::eb_interp_centroid2facecent_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_centroid2facecent_z(Box const &wbx, Array4< Real const > const &phi, Array4< Real const > const &apz, Array4< Real const > const &cvol, Array4< Real const > const &ccent, Array4< Real const > const &fcz, Array4< Real > const &phi_z, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_3D_C.H:1061
amrex::eb_interp_cc2face_y
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::Abort
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:225
amrex::eb_interp_cc2cent
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::EB_interp_in_quad
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real EB_interp_in_quad(Real xint, Real yint, Real v0, Real v1, Real v2, Real v3, Real x0, Real y0, Real x1, Real y1, Real x2, Real y2, Real x3, Real y3)
Definition AMReX_EBMultiFabUtil_3D_C.H:10
amrex::eb_interp_cc2facecent_x
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::eb_interp_cc2face_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void eb_interp_cc2face_z(Box const &wbx, Array4< Real const > const &phi, Array4< Real > const &edg_z, int ncomp, const Box &domain, const BCRec *bc) noexcept
Definition AMReX_EBMultiFabUtil_3D_C.H:1332
amrex::eb_interp_centroid2facecent_x
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::Loop
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void Loop(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:125
amrex::Array4
Definition AMReX_Array4.H:61
amrex::Dim3
Definition AMReX_Dim3.H:12
amrex::Dim3::x
int x
Definition AMReX_Dim3.H:12
amrex::Dim3::z
int z
Definition AMReX_Dim3.H:12
amrex::Dim3::y
int y
Definition AMReX_Dim3.H:12
amrex::GpuArray
Definition AMReX_Array.H:34