Block-Structured AMR Software Framework
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 
6 namespace amrex {
7 
8 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
9 Real
10 EB_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 
50 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
51 void 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 
63 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
64 void 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 
76 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
77 void 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 
101 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
102 void 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 
126 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
127 void 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 
151 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
152 void 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 
176 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
177 void 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 
201 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
202 void 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 
226 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
227 void 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 
334 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
335 void 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 
372 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
373 void 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 
421 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
422 void 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 
488 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
489 void 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 
555 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
556 void 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 
622 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
623 void 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 
842 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
843 void 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 
1060 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1061 void 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 
1277 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1278 void 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 
1304 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1305 void 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 
1331 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1332 void 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 
1358 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1359 void 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::SundialsUserFun::f
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
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::abs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
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