Block-Structured AMR Software Framework
AMReX_MLEBABecLap_3D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_MLEBABECLAP_3D_K_H_
2 #define AMREX_MLEBABECLAP_3D_K_H_
3 #include <AMReX_Config.H>
4 #include <AMReX_REAL.H>
5 
7 
8 namespace amrex {
9 
11 void mlebabeclap_adotx_centroid (Box const& box, Array4<Real> const& y,
12  Array4<Real const> const& x, Array4<Real const> const& a,
13  Array4<Real const> const& bX, Array4<Real const> const& bY,
14  Array4<Real const> const& bZ,
15  Array4<EBCellFlag const> const& flag,
16  Array4<Real const> const& vfrc, Array4<Real const> const& apx,
17  Array4<Real const> const& apy, Array4<Real const> const& apz,
18  Array4<Real const> const& fcx, Array4<Real const> const& fcy,
19  Array4<Real const> const& fcz,
20  Array4<Real const> const& ccent, Array4<Real const> const& ba,
21  Array4<Real const> const& bcent, Array4<Real const> const& beb,
22  Array4<Real const> const& phieb,
23  const int& domlo_x, const int& domlo_y, const int& domlo_z,
24  const int& domhi_x, const int& domhi_y, const int& domhi_z,
25  const bool& on_x_face, const bool& on_y_face, const bool& on_z_face,
26  bool is_eb_dirichlet, bool is_eb_inhomog,
27  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
28  Real alpha, Real beta, int ncomp) noexcept
29 {
30  Real dhx = beta*dxinv[0]*dxinv[0];
31  Real dhy = beta*dxinv[1]*dxinv[1];
32  Real dhz = beta*dxinv[2]*dxinv[2];
33 
34  amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
35  {
36  if (flag(i,j,k).isCovered())
37  {
38  y(i,j,k,n) = Real(0.0);
39  }
40  else if (flag(i,j,k).isRegular() &&
41  ((flag(i-1,j ,k ).isRegular() && flag(i+1,j ,k ).isRegular() &&
42  flag(i ,j-1,k ).isRegular() && flag(i ,j+1,k ).isRegular() &&
43  flag(i ,j ,k-1).isRegular() && flag(i ,j ,k+1).isRegular()) ))
44  {
45  y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n)
46  - dhx * (bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i ,j,k,n))
47  -bX(i ,j,k,n)*(x(i ,j,k,n) - x(i-1,j,k,n)))
48  - dhy * (bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j ,k,n))
49  -bY(i,j ,k,n)*(x(i,j ,k,n) - x(i,j-1,k,n)))
50  - dhz * (bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k ,n))
51  -bZ(i,j,k ,n)*(x(i,j,k ,n) - x(i,j,k-1,n)));
52  }
53  else
54  {
55  Real kappa = vfrc(i,j,k);
56  Real apxm = apx(i,j,k);
57  Real apxp = apx(i+1,j,k);
58  Real apym = apy(i,j,k);
59  Real apyp = apy(i,j+1,k);
60  Real apzm = apz(i,j,k);
61  Real apzp = apz(i,j,k+1);
62 
63  // First get EB-aware slope that doesn't know about extdir
64  bool needs_bdry_stencil = (i <= domlo_x) || (i >= domhi_x) ||
65  (j <= domlo_y) || (j >= domhi_y) ||
66  (k <= domlo_z) || (k >= domhi_z);
67 
68  Real fxm = bX(i,j,k,n)*(x(i,j,k,n) - x(i-1,j,k,n));
69  if ( (apxm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0)) ) {
70  Real yloc_on_xface = fcx(i,j,k,0);
71  Real zloc_on_xface = fcx(i,j,k,1);
72 
73  if(needs_bdry_stencil) {
74  fxm = grad_x_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc,
75  yloc_on_xface,zloc_on_xface,
76  is_eb_dirichlet,is_eb_inhomog,
77  on_x_face,domlo_x,domhi_x,
78  on_y_face,domlo_y,domhi_y,
79  on_z_face,domlo_z,domhi_z);
80  } else {
81  fxm = grad_x_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent,
82  yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
83  }
84 
85  fxm *= bX(i,j,k,n);
86  }
87 
88  Real fxp = bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i,j,k,n));
89  if ( (apxp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0)) ) {
90  Real yloc_on_xface = fcx(i+1,j,k,0);
91  Real zloc_on_xface = fcx(i+1,j,k,1);
92 
93  if(needs_bdry_stencil) {
94  fxp = grad_x_of_phi_on_centroids_extdir(i+1,j,k,n,x,phieb,flag,ccent,bcent,vfrc,
95  yloc_on_xface,zloc_on_xface,
96  is_eb_dirichlet,is_eb_inhomog,
97  on_x_face,domlo_x,domhi_x,
98  on_y_face,domlo_y,domhi_y,
99  on_z_face,domlo_z,domhi_z);
100  } else {
101  fxp = grad_x_of_phi_on_centroids(i+1,j,k,n,x,phieb,flag,ccent,bcent,
102  yloc_on_xface,zloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
103  }
104 
105  fxp *= bX(i+1,j,k,n);
106  }
107 
108  Real fym = bY(i,j,k,n)*(x(i,j,k,n) - x(i,j-1,k,n));
109  if ( (apym != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0)) ) {
110  Real xloc_on_yface = fcy(i,j,k,0);
111  Real zloc_on_yface = fcy(i,j,k,1);
112 
113  if(needs_bdry_stencil) {
114  fym = grad_y_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc,
115  xloc_on_yface,zloc_on_yface,
116  is_eb_dirichlet,is_eb_inhomog,
117  on_x_face,domlo_x,domhi_x,
118  on_y_face,domlo_y,domhi_y,
119  on_z_face,domlo_z,domhi_z);
120  } else {
121  fym = grad_y_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent,
122  xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
123  }
124 
125  fym *= bY(i,j,k,n);
126  }
127 
128  Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j,k,n));
129  if ( (apyp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0)) ) {
130  Real xloc_on_yface = fcy(i,j+1,k,0);
131  Real zloc_on_yface = fcy(i,j+1,k,1);
132 
133  if(needs_bdry_stencil) {
134  fyp = grad_y_of_phi_on_centroids_extdir(i,j+1,k,n,x,phieb,flag,ccent,bcent,vfrc,
135  xloc_on_yface,zloc_on_yface,
136  is_eb_dirichlet,is_eb_inhomog,
137  on_x_face,domlo_x,domhi_x,
138  on_y_face,domlo_y,domhi_y,
139  on_z_face,domlo_z,domhi_z);
140  } else {
141  fyp = grad_y_of_phi_on_centroids(i,j+1,k,n,x,phieb,flag,ccent,bcent,
142  xloc_on_yface,zloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
143  }
144 
145  fyp *= bY(i,j+1,k,n);
146  }
147 
148  Real fzm = bZ(i,j,k,n)*(x(i,j,k,n) - x(i,j,k-1,n));
149  if ( (apzm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0)) ) {
150  Real xloc_on_zface = fcz(i,j,k,0);
151  Real yloc_on_zface = fcz(i,j,k,1);
152 
153  if(needs_bdry_stencil) {
154  fzm = grad_z_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc,
155  xloc_on_zface,yloc_on_zface,
156  is_eb_dirichlet,is_eb_inhomog,
157  on_x_face,domlo_x,domhi_x,
158  on_y_face,domlo_y,domhi_y,
159  on_z_face,domlo_z,domhi_z);
160  } else {
161  fzm = grad_z_of_phi_on_centroids(i,j,k,n,x,phieb,flag,ccent,bcent,
162  xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog);
163  }
164 
165  fzm *= bZ(i,j,k,n);
166  }
167 
168  Real fzp = bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k,n));
169  if ( (apzp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j,k+1) != Real(1.0) || vfrc(i,j,k-1) != Real(1.0)) ) {
170  Real xloc_on_zface = fcz(i,j,k+1,0);
171  Real yloc_on_zface = fcz(i,j,k+1,1);
172 
173  if(needs_bdry_stencil) {
174  fzp = grad_z_of_phi_on_centroids_extdir(i,j,k+1,n,x,phieb,flag,ccent,bcent,vfrc,
175  xloc_on_zface,yloc_on_zface,
176  is_eb_dirichlet,is_eb_inhomog,
177  on_x_face,domlo_x,domhi_x,
178  on_y_face,domlo_y,domhi_y,
179  on_z_face,domlo_z,domhi_z);
180  } else {
181  fzp = grad_z_of_phi_on_centroids(i,j,k+1,n,x,phieb,flag,ccent,bcent,
182  xloc_on_zface,yloc_on_zface,is_eb_dirichlet,is_eb_inhomog);
183  }
184 
185  fzp *= bZ(i,j,k+1,n);
186  }
187 
188  Real feb = Real(0.0);
189  if (is_eb_dirichlet && flag(i,j,k).isSingleValued()) {
190  Real dapx = apxm-apxp;
191  Real dapy = apym-apyp;
192  Real dapz = apzm-apzp;
193  Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
194  Real anorminv = Real(1.0)/anorm;
195  Real anrmx = dapx * anorminv;
196  Real anrmy = dapy * anorminv;
197  Real anrmz = dapz * anorminv;
198 
199  feb = grad_eb_of_phi_on_centroids_extdir(i,j,k,n,x,phieb,flag,ccent,bcent,vfrc,
200  anrmx,anrmy,anrmz,is_eb_inhomog,
201  on_x_face,domlo_x,domhi_x,
202  on_y_face,domlo_y,domhi_y,
203  on_z_face,domlo_z,domhi_z);
204  feb *= ba(i,j,k) * beb(i,j,k,n);
205  }
206 
207  y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) *
208  (dhx*(apxm*fxm - apxp*fxp) +
209  dhy*(apym*fym - apyp*fyp) +
210  dhz*(apzm*fzm - apzp*fzp) - dhx*feb);
211  }
212  });
213 }
214 
216 void mlebabeclap_adotx (Box const& box, Array4<Real> const& y,
217  Array4<Real const> const& x, Array4<Real const> const& a,
218  Array4<Real const> const& bX, Array4<Real const> const& bY,
219  Array4<Real const> const& bZ, Array4<const int> const& ccm,
220  Array4<EBCellFlag const> const& flag,
221  Array4<Real const> const& vfrc, Array4<Real const> const& apx,
222  Array4<Real const> const& apy, Array4<Real const> const& apz,
223  Array4<Real const> const& fcx, Array4<Real const> const& fcy,
224  Array4<Real const> const& fcz, Array4<Real const> const& ba,
225  Array4<Real const> const& bc, Array4<Real const> const& beb,
226  bool is_dirichlet, Array4<Real const> const& phieb,
227  bool is_inhomog, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
228  Real alpha, Real beta, int ncomp,
229  bool beta_on_centroid, bool phi_on_centroid) noexcept
230 {
231  Real dhx = beta*dxinv[0]*dxinv[0];
232  Real dhy = beta*dxinv[1]*dxinv[1];
233  Real dhz = beta*dxinv[2]*dxinv[2];
234 
235  bool beta_on_center = !(beta_on_centroid);
236  bool phi_on_center = !( phi_on_centroid);
237 
238  amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
239  {
240  if (flag(i,j,k).isCovered())
241  {
242  y(i,j,k,n) = Real(0.0);
243  }
244  else if (flag(i,j,k).isRegular())
245  {
246  y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n)
247  - dhx * (bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i ,j,k,n))
248  -bX(i ,j,k,n)*(x(i ,j,k,n) - x(i-1,j,k,n)))
249  - dhy * (bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j ,k,n))
250  -bY(i,j ,k,n)*(x(i,j ,k,n) - x(i,j-1,k,n)))
251  - dhz * (bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k ,n))
252  -bZ(i,j,k ,n)*(x(i,j,k ,n) - x(i,j,k-1,n)));
253  }
254  else
255  {
256  Real kappa = vfrc(i,j,k);
257  Real apxm = apx(i,j,k);
258  Real apxp = apx(i+1,j,k);
259  Real apym = apy(i,j,k);
260  Real apyp = apy(i,j+1,k);
261  Real apzm = apz(i,j,k);
262  Real apzp = apz(i,j,k+1);
263 
264  Real fxm = bX(i,j,k,n)*(x(i,j,k,n) - x(i-1,j,k,n));
265  if (apxm != Real(0.0) && apxm != Real(1.0)) {
266  int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
267  int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
268  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
269  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
270  if (beta_on_center && phi_on_center)
271  {
272  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm +
273  fracy*(Real(1.0)-fracz)*bX(i,jj,k ,n)*(x(i,jj,k ,n)-x(i-1,jj,k ,n)) +
274  fracz*(Real(1.0)-fracy)*bX(i,j ,kk,n)*(x(i,j ,kk,n)-x(i-1,j ,kk,n)) +
275  fracy* fracz *bX(i,jj,kk,n)*(x(i,jj,kk,n)-x(i-1,jj,kk,n));
276  }
277  else if (beta_on_centroid && phi_on_center)
278  {
279  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(x(i, j, k,n)-x(i-1, j, k,n)) +
280  fracy *(Real(1.0)-fracz)*(x(i,jj, k,n)-x(i-1,jj, k,n)) +
281  fracz *(Real(1.0)-fracy)*(x(i, j,kk,n)-x(i-1, j,kk,n)) +
282  fracy * fracz *(x(i,jj,kk,n)-x(i-1,jj,kk,n));
283  fxm *= bX(i,j,k,n);
284  }
285  }
286 
287  Real fxp = bX(i+1,j,k,n)*(x(i+1,j,k,n) - x(i,j,k,n));
288  if (apxp != Real(0.0) && apxp != Real(1.0)) {
289  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,0)));
290  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,1)));
291  Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k,0)) : Real(0.0);
292  Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk)) ? std::abs(fcx(i+1,j,k,1)) : Real(0.0);
293  if (beta_on_center && phi_on_center)
294  {
295  fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp +
296  fracy*(Real(1.0)-fracz)*bX(i+1,jj,k ,n)*(x(i+1,jj,k ,n)-x(i,jj,k ,n)) +
297  fracz*(Real(1.0)-fracy)*bX(i+1,j ,kk,n)*(x(i+1,j ,kk,n)-x(i,j ,kk,n)) +
298  fracy* fracz *bX(i+1,jj,kk,n)*(x(i+1,jj,kk,n)-x(i,jj,kk,n));
299  }
300  else if (beta_on_centroid && phi_on_center)
301  {
302  fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(x(i+1, j, k,n)-x(i, j, k,n)) +
303  fracy *(Real(1.0)-fracz)*(x(i+1,jj, k,n)-x(i,jj, k,n)) +
304  fracz *(Real(1.0)-fracy)*(x(i+1, j,kk,n)-x(i, j,kk,n)) +
305  fracy * fracz *(x(i+1,jj,kk,n)-x(i,jj,kk,n));
306  fxp *= bX(i+1,j,k,n);
307 
308  }
309  }
310 
311  Real fym = bY(i,j,k,n)*(x(i,j,k,n) - x(i,j-1,k,n));
312  if (apym != Real(0.0) && apym != Real(1.0)) {
313  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
314  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
315  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
316  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
317  if (beta_on_center && phi_on_center)
318  {
319  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym +
320  fracx*(Real(1.0)-fracz)*bY(ii,j,k ,n)*(x(ii,j,k ,n)-x(ii,j-1,k ,n)) +
321  fracz*(Real(1.0)-fracx)*bY(i ,j,kk,n)*(x(i ,j,kk,n)-x(i ,j-1,kk,n)) +
322  fracx* fracz *bY(ii,j,kk,n)*(x(ii,j,kk,n)-x(ii,j-1,kk,n));
323  }
324  else if (beta_on_centroid && phi_on_center)
325  {
326  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(x( i,j, k,n)-x( i,j-1, k,n)) +
327  fracx *(Real(1.0)-fracz)*(x(ii,j, k,n)-x(ii,j-1, k,n)) +
328  fracz *(Real(1.0)-fracx)*(x(i ,j,kk,n)-x( i,j-1,kk,n)) +
329  fracx * fracz *(x(ii,j,kk,n)-x(ii,j-1,kk,n));
330  fym *= bY(i,j,k,n);
331 
332  }
333  }
334 
335  Real fyp = bY(i,j+1,k,n)*(x(i,j+1,k,n) - x(i,j,k,n));
336  if (apyp != Real(0.0) && apyp != Real(1.0)) {
337  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,0)));
338  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,1)));
339  Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k,0)) : Real(0.0);
340  Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk)) ? std::abs(fcy(i,j+1,k,1)) : Real(0.0);
341  if (beta_on_center && phi_on_center)
342  {
343  fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp +
344  fracx*(Real(1.0)-fracz)*bY(ii,j+1,k ,n)*(x(ii,j+1,k ,n)-x(ii,j,k ,n)) +
345  fracz*(Real(1.0)-fracx)*bY(i ,j+1,kk,n)*(x(i ,j+1,kk,n)-x(i ,j,kk,n)) +
346  fracx* fracz *bY(ii,j+1,kk,n)*(x(ii,j+1,kk,n)-x(ii,j,kk,n));
347  }
348  else if (beta_on_centroid && phi_on_center)
349  {
350  fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(x( i,j+1, k,n)-x( i,j, k,n)) +
351  fracx *(Real(1.0)-fracz)*(x(ii,j+1, k,n)-x(ii,j, k,n)) +
352  fracz *(Real(1.0)-fracx)*(x( i,j+1,kk,n)-x( i,j,kk,n)) +
353  fracx * fracz *(x(ii,j+1,kk,n)-x(ii,j,kk,n));
354  fyp *= bY(i,j+1,k,n);
355 
356  }
357  }
358 
359  Real fzm = bZ(i,j,k,n)*(x(i,j,k,n) - x(i,j,k-1,n));
360  if (apzm != Real(0.0) && apzm != Real(1.0)) {
361  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
362  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
363  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
364  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
365  if (beta_on_center && phi_on_center)
366  {
367  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm +
368  fracx*(Real(1.0)-fracy)*bZ(ii,j ,k,n)*(x(ii,j ,k,n)-x(ii,j ,k-1,n)) +
369  fracy*(Real(1.0)-fracx)*bZ(i ,jj,k,n)*(x(i ,jj,k,n)-x(i ,jj,k-1,n)) +
370  fracx* fracy *bZ(ii,jj,k,n)*(x(ii,jj,k,n)-x(ii,jj,k-1,n));
371  }
372  else if (beta_on_centroid && phi_on_center)
373  {
374  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(x( i, j,k,n)-x( i, j,k-1,n)) +
375  fracx *(Real(1.0)-fracy)*(x(ii, j,k,n)-x(ii, j,k-1,n)) +
376  fracy *(Real(1.0)-fracx)*(x( i,jj,k,n)-x( i,jj,k-1,n)) +
377  fracx * fracy *(x(ii,jj,k,n)-x(ii,jj,k-1,n));
378  fzm *= bZ(i,j,k,n);
379 
380  }
381  }
382 
383  Real fzp = bZ(i,j,k+1,n)*(x(i,j,k+1,n) - x(i,j,k,n));
384  if (apzp != Real(0.0) && apzp != Real(1.0)) {
385  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,0)));
386  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,1)));
387  Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1)) ? std::abs(fcz(i,j,k+1,0)) : Real(0.0);
388  Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1)) ? std::abs(fcz(i,j,k+1,1)) : Real(0.0);
389  if (beta_on_center && phi_on_center)
390  {
391  fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp +
392  fracx*(Real(1.0)-fracy)*bZ(ii,j ,k+1,n)*(x(ii,j ,k+1,n)-x(ii,j ,k,n)) +
393  fracy*(Real(1.0)-fracx)*bZ(i ,jj,k+1,n)*(x(i ,jj,k+1,n)-x(i ,jj,k,n)) +
394  fracx* fracy *bZ(ii,jj,k+1,n)*(x(ii,jj,k+1,n)-x(ii,jj,k,n));
395  }
396  else if (beta_on_centroid && phi_on_center)
397  {
398  fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(x( i, j,k+1,n)-x( i, j,k,n)) +
399  fracx *(Real(1.0)-fracy)*(x(ii, j,k+1,n)-x(ii, j,k,n)) +
400  fracy *(Real(1.0)-fracx)*(x( i,jj,k+1,n)-x( i,jj,k,n)) +
401  fracx * fracy *(x(ii,jj,k+1,n)-x(ii,jj,k,n));
402  fzp *= bZ(i,j,k+1,n);
403 
404  }
405  }
406 
407  Real feb = Real(0.0);
408  if (is_dirichlet) {
409  Real dapx = apxm-apxp;
410  Real dapy = apym-apyp;
411  Real dapz = apzm-apzp;
412  Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
413  Real anorminv = Real(1.0)/anorm;
414  Real anrmx = dapx * anorminv;
415  Real anrmy = dapy * anorminv;
416  Real anrmz = dapz * anorminv;
417 
418  Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
419 
420  Real bctx = bc(i,j,k,0);
421  Real bcty = bc(i,j,k,1);
422  Real bctz = bc(i,j,k,2);
423  Real dx_eb = get_dx_eb(kappa);
424 
425  Real dg = dx_eb / amrex::max(std::abs(anrmx), std::abs(anrmy),
426  std::abs(anrmz));
427  Real gx = bctx - dg*anrmx;
428  Real gy = bcty - dg*anrmy;
429  Real gz = bctz - dg*anrmz;
430  Real sx = std::copysign(Real(1.0),anrmx);
431  Real sy = std::copysign(Real(1.0),anrmy);
432  Real sz = std::copysign(Real(1.0),anrmz);
433  int ii = i - static_cast<int>(sx);
434  int jj = j - static_cast<int>(sy);
435  int kk = k - static_cast<int>(sz);
436 
437  gx = sx*gx;
438  gy = sy*gy;
439  gz = sz*gz;
440  Real gxy = gx*gy;
441  Real gxz = gx*gz;
442  Real gyz = gy*gz;
443  Real gxyz = gx*gy*gz;
444  Real phig = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz) * x(i ,j ,k ,n)
445  + (-gz - gxz - gyz - gxyz) * x(i ,j ,kk,n)
446  + (-gy - gxy - gyz - gxyz) * x(i ,jj,k ,n)
447  + (gyz + gxyz) * x(i ,jj,kk,n)
448  + (-gx - gxy - gxz - gxyz) * x(ii,j ,k ,n)
449  + (gxz + gxyz) * x(ii,j ,kk,n)
450  + (gxy + gxyz) * x(ii,jj,k ,n)
451  + (-gxyz) * x(ii,jj,kk,n);
452 
453  Real dphidn = (phib-phig)/dg;
454 
455  feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
456  }
457 
458  y(i,j,k,n) = alpha*a(i,j,k)*x(i,j,k,n) + (Real(1.0)/kappa) *
459  (dhx*(apxm*fxm - apxp*fxp) +
460  dhy*(apym*fym - apyp*fyp) +
461  dhz*(apzm*fzm - apzp*fzp) - dhx*feb);
462  }
463  });
464 }
465 
467 void mlebabeclap_ebflux (int i, int j, int k, int n,
468  Array4<Real> const& feb,
469  Array4<Real const> const& x,
470  Array4<EBCellFlag const> const& flag,
471  Array4<Real const> const& vfrc,
472  Array4<Real const> const& apx,
473  Array4<Real const> const& apy,
474  Array4<Real const> const& apz,
475  Array4<Real const> const& bc,
476  Array4<Real const> const& beb,
477  Array4<Real const> const& phieb,
478  bool is_inhomog,
479  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
480 {
481  Real dhx = dxinv[0];
482 
483  if (!flag(i,j,k).isSingleValued())
484  {
485  feb(i,j,k,n) = Real(0.0);
486  }
487  else
488  {
489  Real kappa = vfrc(i,j,k);
490  Real apxm = apx(i,j,k);
491  Real apxp = apx(i+1,j,k);
492  Real apym = apy(i,j,k);
493  Real apyp = apy(i,j+1,k);
494  Real apzm = apz(i,j,k);
495  Real apzp = apz(i,j,k+1);
496 
497  Real dapx = apxm-apxp;
498  Real dapy = apym-apyp;
499  Real dapz = apzm-apzp;
500  Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
501  Real anorminv = Real(1.0)/anorm;
502  Real anrmx = dapx * anorminv;
503  Real anrmy = dapy * anorminv;
504  Real anrmz = dapz * anorminv;
505 
506  Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
507 
508  Real bctx = bc(i,j,k,0);
509  Real bcty = bc(i,j,k,1);
510  Real bctz = bc(i,j,k,2);
511  Real dx_eb = get_dx_eb(kappa);
512 
513  Real dg = dx_eb / amrex::max(std::abs(anrmx), std::abs(anrmy), std::abs(anrmz));
514  Real gx = bctx - dg*anrmx;
515  Real gy = bcty - dg*anrmy;
516  Real gz = bctz - dg*anrmz;
517  Real sx = std::copysign(Real(1.0),anrmx);
518  Real sy = std::copysign(Real(1.0),anrmy);
519  Real sz = std::copysign(Real(1.0),anrmz);
520  int ii = i - static_cast<int>(sx);
521  int jj = j - static_cast<int>(sy);
522  int kk = k - static_cast<int>(sz);
523 
524  gx = sx*gx;
525  gy = sy*gy;
526  gz = sz*gz;
527  Real gxy = gx*gy;
528  Real gxz = gx*gz;
529  Real gyz = gy*gz;
530  Real gxyz = gx*gy*gz;
531  Real phig = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz) * x(i ,j ,k ,n)
532  + (-gz - gxz - gyz - gxyz) * x(i ,j ,kk,n)
533  + (-gy - gxy - gyz - gxyz) * x(i ,jj,k ,n)
534  + (gyz + gxyz) * x(i ,jj,kk,n)
535  + (-gx - gxy - gxz - gxyz) * x(ii,j ,k ,n)
536  + (gxz + gxyz) * x(ii,j ,kk,n)
537  + (gxy + gxyz) * x(ii,jj,k ,n)
538  + (-gxyz) * x(ii,jj,kk,n);
539 
540  Real dphidn = dhx*(phib-phig)/dg;
541  feb(i,j,k,n) = -beb(i,j,k,n) * dphidn;
542  }
543 }
544 
546 void mlebabeclap_gsrb (Box const& box,
547  Array4<Real> const& phi, Array4<Real const> const& rhs,
548  Real alpha, Array4<Real const> const& a,
549  Real dhx, Real dhy, Real dhz,
550  Array4<Real const> const& bX, Array4<Real const> const& bY,
551  Array4<Real const> const& bZ,
552  Array4<int const> const& m0, Array4<int const> const& m2,
553  Array4<int const> const& m4,
554  Array4<int const> const& m1, Array4<int const> const& m3,
555  Array4<int const> const& m5,
556  Array4<Real const> const& f0, Array4<Real const> const& f2,
557  Array4<Real const> const& f4,
558  Array4<Real const> const& f1, Array4<Real const> const& f3,
559  Array4<Real const> const& f5,
560  Array4<const int> const& ccm, Array4<Real const> const& beb,
561  EBData const& ebdata,
562  bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid,
563  Box const& vbox, int redblack, int ncomp) noexcept
564 {
565  constexpr Real omega = 1.15;
566 
567  const auto vlo = amrex::lbound(vbox);
568  const auto vhi = amrex::ubound(vbox);
569 
570 // amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
571  // amrex::Loop here causes gcc 8 to crash.
572  const auto lo = amrex::lbound(box);
573  const auto hi = amrex::ubound(box);
574  for (int n = 0; n < ncomp; ++n) {
575  for (int k = lo.z; k <= hi.z; ++k) {
576  for (int j = lo.y; j <= hi.y; ++j) {
577  for (int i = lo.x; i <= hi.x; ++i)
578  {
579  if ((i+j+k+redblack) % 2 == 0)
580  {
581  auto const flag = ebdata.get<EBData_t::cellflag>(i,j,k);
582  if (flag.isCovered())
583  {
584  phi(i,j,k,n) = Real(0.0);
585  }
586  else
587  {
588  Real cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
589  ? f0(vlo.x,j,k,n) : Real(0.0);
590  Real cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
591  ? f1(i,vlo.y,k,n) : Real(0.0);
592  Real cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
593  ? f2(i,j,vlo.z,n) : Real(0.0);
594  Real cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
595  ? f3(vhi.x,j,k,n) : Real(0.0);
596  Real cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
597  ? f4(i,vhi.y,k,n) : Real(0.0);
598  Real cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
599  ? f5(i,j,vhi.z,n) : Real(0.0);
600 
601  if (flag.isRegular())
602  {
603  Real gamma = alpha*a(i,j,k)
604  + dhx*(bX(i+1,j,k,n) + bX(i,j,k,n))
605  + dhy*(bY(i,j+1,k,n) + bY(i,j,k,n))
606  + dhz*(bZ(i,j,k+1,n) + bZ(i,j,k,n));
607 
608  Real rho = dhx*(bX(i+1,j ,k ,n)*phi(i+1,j ,k ,n) +
609  bX(i ,j ,k ,n)*phi(i-1,j ,k ,n))
610  + dhy*(bY(i ,j+1,k ,n)*phi(i ,j+1,k ,n) +
611  bY(i ,j ,k ,n)*phi(i ,j-1,k ,n))
612  + dhz*(bZ(i ,j ,k+1,n)*phi(i ,j ,k+1,n) +
613  bZ(i ,j ,k ,n)*phi(i ,j ,k-1,n));
614 
615  Real delta = dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
616  + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
617  + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5);
618 
619  Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
620  phi(i,j,k,n) += omega*res/(gamma-delta);
621  }
622  else
623  {
624  Real kappa = ebdata.get<EBData_t::volfrac>(i,j,k);
625  Real apxm = ebdata.get<EBData_t::apx>(i ,j ,k );
626  Real apxp = ebdata.get<EBData_t::apx>(i+1,j ,k );
627  Real apym = ebdata.get<EBData_t::apy>(i ,j ,k );
628  Real apyp = ebdata.get<EBData_t::apy>(i ,j+1,k );
629  Real apzm = ebdata.get<EBData_t::apz>(i ,j ,k );
630  Real apzp = ebdata.get<EBData_t::apz>(i ,j ,k+1);
631 
632  Real fxm = -bX(i,j,k,n)*phi(i-1,j,k,n);
633  Real oxm = -bX(i,j,k,n)*cf0;
634  Real sxm = bX(i,j,k,n);
635  if (apxm != Real(0.0) && apxm != Real(1.0)) {
636  auto fcx0 = ebdata.get<EBData_t::fcx>(i,j,k,0);
637  auto fcx1 = ebdata.get<EBData_t::fcx>(i,j,k,1);
638  int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx0));
639  int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx1));
640  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
641  ? std::abs(fcx0) : Real(0.0);
642  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk))
643  ? std::abs(fcx1) : Real(0.0);
644  if (!beta_on_centroid && !phi_on_centroid)
645  {
646  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm
647  + fracy *(Real(1.0)-fracz)*bX(i,jj,k ,n)*(phi(i,jj,k ,n)-phi(i-1,jj,k ,n))
648  +(Real(1.0)-fracy)* fracz *bX(i,j ,kk,n)*(phi(i,j ,kk,n)-phi(i-1,j ,kk,n))
649  + fracy * fracz *bX(i,jj,kk,n)*(phi(i,jj,kk,n)-phi(i-1,jj,kk,n));
650  }
651  else if (beta_on_centroid && !phi_on_centroid)
652  {
653  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*( -phi(i-1, j, k,n))
654  + fracy *(Real(1.0)-fracz)*(phi(i,jj,k ,n)-phi(i-1,jj, k,n))
655  +(Real(1.0)-fracy)* fracz *(phi(i,j ,kk,n)-phi(i-1, j,kk,n))
656  + fracy * fracz *(phi(i,jj,kk,n)-phi(i-1,jj,kk,n));
657  fxm *= bX(i,j,k,n);
658 
659  }
660  oxm = Real(0.0);
661  sxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxm;
662  }
663 
664  Real fxp = bX(i+1,j,k,n)*phi(i+1,j,k,n);
665  Real oxp = bX(i+1,j,k,n)*cf3;
666  Real sxp = -bX(i+1,j,k,n);
667  if (apxp != Real(0.0) && apxp != Real(1.0)) {
668  auto fcx0 = ebdata.get<EBData_t::fcx>(i+1,j,k,0);
669  auto fcx1 = ebdata.get<EBData_t::fcx>(i+1,j,k,1);
670  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcx0));
671  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcx1));
672  Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
673  ? std::abs(fcx0) : Real(0.0);
674  Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk))
675  ? std::abs(fcx1) : Real(0.0);
676  if (!beta_on_centroid && !phi_on_centroid)
677  {
678  fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxp
679  + fracy *(Real(1.0)-fracz)*bX(i+1,jj,k ,n)*(phi(i+1,jj,k ,n)-phi(i,jj,k ,n))
680  +(Real(1.0)-fracy)* fracz *bX(i+1,j ,kk,n)*(phi(i+1,j ,kk,n)-phi(i,j ,kk,n))
681  + fracy * fracz *bX(i+1,jj,kk,n)*(phi(i+1,jj,kk,n)-phi(i,jj,kk,n));
682  }
683  else if (beta_on_centroid && !phi_on_centroid)
684  {
685  fxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(phi(i+1, j, k,n) ) +
686  fracy *(Real(1.0)-fracz)*(phi(i+1,jj, k,n)-phi(i,jj, k,n)) +
687  fracz *(Real(1.0)-fracy)*(phi(i+1, j,kk,n)-phi(i, j,kk,n)) +
688  fracy * fracz *(phi(i+1,jj,kk,n)-phi(i,jj,kk,n));
689  fxp *= bX(i+1,j,k,n);
690 
691  }
692 
693  oxp = Real(0.0);
694  sxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxp;
695  }
696 
697  Real fym = -bY(i,j,k,n)*phi(i,j-1,k,n);
698  Real oym = -bY(i,j,k,n)*cf1;
699  Real sym = bY(i,j,k,n);
700  if (apym != Real(0.0) && apym != Real(1.0)) {
701  auto fcy0 = ebdata.get<EBData_t::fcy>(i,j,k,0);
702  auto fcy1 = ebdata.get<EBData_t::fcy>(i,j,k,1);
703  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy0));
704  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy1));
705  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
706  ? std::abs(fcy0) : Real(0.0);
707  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk))
708  ? std::abs(fcy1) : Real(0.0);
709  if (!beta_on_centroid && !phi_on_centroid)
710  {
711  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym
712  + fracx *(Real(1.0)-fracz)*bY(ii,j,k ,n)*(phi(ii,j,k ,n)-phi(ii,j-1,k ,n))
713  + (Real(1.0)-fracx)* fracz *bY(i ,j,kk,n)*(phi(i ,j,kk,n)-phi(i ,j-1,kk,n))
714  + fracx * fracz *bY(ii,j,kk,n)*(phi(ii,j,kk,n)-phi(ii,j-1,kk,n));
715  }
716  else if (beta_on_centroid && !phi_on_centroid)
717  {
718  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*( -phi( i,j-1, k,n))
719  + fracx *(Real(1.0)-fracz)*(phi(ii,j,k ,n)-phi(ii,j-1, k,n))
720  + (Real(1.0)-fracx)* fracz *(phi(i ,j,kk,n)-phi( i,j-1,kk,n))
721  + fracx * fracz *(phi(ii,j,kk,n)-phi(ii,j-1,kk,n));
722  fym *= bY(i,j,k,n);
723 
724  }
725  oym = Real(0.0);
726  sym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*sym;
727  }
728 
729  Real fyp = bY(i,j+1,k,n)*phi(i,j+1,k,n);
730  Real oyp = bY(i,j+1,k,n)*cf4;
731  Real syp = -bY(i,j+1,k,n);
732  if (apyp != Real(0.0) && apyp != Real(1.0)) {
733  auto fcy0 = ebdata.get<EBData_t::fcy>(i,j+1,k,0);
734  auto fcy1 = ebdata.get<EBData_t::fcy>(i,j+1,k,1);
735  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy0));
736  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy1));
737  Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
738  ? std::abs(fcy0) : Real(0.0);
739  Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk))
740  ? std::abs(fcy1) : Real(0.0);
741  if (!beta_on_centroid && !phi_on_centroid)
742  {
743  fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fyp
744  + fracx *(Real(1.0)-fracz)*bY(ii,j+1,k ,n)*(phi(ii,j+1,k ,n)-phi(ii,j,k ,n))
745  + (Real(1.0)-fracx)* fracz *bY(i ,j+1,kk,n)*(phi(i ,j+1,kk,n)-phi(i ,j,kk,n))
746  + fracx * fracz *bY(ii,j+1,kk,n)*(phi(ii,j+1,kk,n)-phi(ii,j,kk,n));
747  }
748  else if (beta_on_centroid && !phi_on_centroid)
749  {
750  fyp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(phi( i,j+1, k,n) )
751  + fracx *(Real(1.0)-fracz)*(phi(ii,j+1, k,n)-phi(ii,j, k,n))
752  + (Real(1.0)-fracx)* fracz *(phi( i,j+1,kk,n)-phi( i,j,kk,n))
753  + fracx * fracz *(phi(ii,j+1,kk,n)-phi(ii,j,kk,n));
754  fyp *= bY(i,j+1,k,n);
755 
756  }
757  oyp = Real(0.0);
758  syp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*syp;
759  }
760 
761  Real fzm = -bZ(i,j,k,n)*phi(i,j,k-1,n);
762  Real ozm = -bZ(i,j,k,n)*cf2;
763  Real szm = bZ(i,j,k,n);
764  if (apzm != Real(0.0) && apzm != Real(1.0)) {
765  auto fcz0 = ebdata.get<EBData_t::fcz>(i,j,k,0);
766  auto fcz1 = ebdata.get<EBData_t::fcz>(i,j,k,1);
767  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz0));
768  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz1));
769  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k))
770  ? std::abs(fcz0) : Real(0.0);
771  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k))
772  ? std::abs(fcz1) : Real(0.0);
773  if (!beta_on_centroid && !phi_on_centroid)
774  {
775  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm
776  + fracx *(Real(1.0)-fracy)*bZ(ii, j,k,n)*(phi(ii, j,k,n)-phi(ii, j,k-1,n))
777  +(Real(1.0)-fracx)* fracy *bZ( i,jj,k,n)*(phi( i,jj,k,n)-phi( i,jj,k-1,n))
778  + fracx * fracy *bZ(ii,jj,k,n)*(phi(ii,jj,k,n)-phi(ii,jj,k-1,n));
779  }
780  else if (beta_on_centroid && !phi_on_centroid)
781  {
782  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*( -phi( i, j,k-1,n))
783  + fracx *(Real(1.0)-fracy)*(phi(ii, j,k,n)-phi(ii, j,k-1,n))
784  + (Real(1.0)-fracx)* fracy *(phi( i,jj,k,n)-phi(i ,jj,k-1,n))
785  + fracx * fracy *(phi(ii,jj,k,n)-phi(ii,jj,k-1,n));
786  fzm *= bZ(i,j,k,n);
787 
788  }
789  ozm = Real(0.0);
790  szm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szm;
791  }
792 
793  Real fzp = bZ(i,j,k+1,n)*phi(i,j,k+1,n);
794  Real ozp = bZ(i,j,k+1,n)*cf5;
795  Real szp = -bZ(i,j,k+1,n);
796  if (apzp != Real(0.0) && apzp != Real(1.0)) {
797  auto fcz0 = ebdata.get<EBData_t::fcz>(i,j,k+1,0);
798  auto fcz1 = ebdata.get<EBData_t::fcz>(i,j,k+1,1);
799  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz0));
800  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz1));
801  Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1))
802  ? std::abs(fcz0) : Real(0.0);
803  Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1))
804  ? std::abs(fcz1) : Real(0.0);
805  if (!beta_on_centroid && !phi_on_centroid)
806  {
807  fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzp
808  + fracx *(Real(1.0)-fracy)*bZ(ii,j ,k+1,n)*(phi(ii,j ,k+1,n)-phi(ii,j ,k,n))
809  + (Real(1.0)-fracx)* fracy *bZ(i ,jj,k+1,n)*(phi(i ,jj,k+1,n)-phi(i ,jj,k,n))
810  + fracx * fracy *bZ(ii,jj,k+1,n)*(phi(ii,jj,k+1,n)-phi(ii,jj,k,n));
811  }
812  else if (beta_on_centroid && !phi_on_centroid)
813  {
814  fzp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(phi( i, j,k+1,n) )
815  + fracx *(Real(1.0)-fracy)*(phi(ii, j,k+1,n)-phi(ii, j,k,n))
816  + (Real(1.0)-fracx)* fracy *(phi( i,jj,k+1,n)-phi( i,jj,k,n))
817  + fracx * fracy *(phi(ii,jj,k+1,n)-phi(ii,jj,k,n));
818  fzp *= bZ(i,j,k+1,n);
819 
820  }
821  ozp = Real(0.0);
822  szp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szp;
823  }
824 
825  Real vfrcinv = Real(1.0)/kappa;
826  Real gamma = alpha*a(i,j,k) + vfrcinv *
827  (dhx*(apxm*sxm-apxp*sxp) +
828  dhy*(apym*sym-apyp*syp) +
829  dhz*(apzm*szm-apzp*szp));
830 
831  Real rho = -vfrcinv *
832  (dhx*(apxm*fxm-apxp*fxp) +
833  dhy*(apym*fym-apyp*fyp) +
834  dhz*(apzm*fzm-apzp*fzp));
835 
836  Real delta = -vfrcinv *
837  (dhx*(apxm*oxm-apxp*oxp) +
838  dhy*(apym*oym-apyp*oyp) +
839  dhz*(apzm*ozm-apzp*ozp));
840 
841  if (is_dirichlet) {
842  Real dapx = apxm-apxp;
843  Real dapy = apym-apyp;
844  Real dapz = apzm-apzp;
845  Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
846  Real anorminv = Real(1.0)/anorm;
847  Real anrmx = dapx * anorminv;
848  Real anrmy = dapy * anorminv;
849  Real anrmz = dapz * anorminv;
850  Real bctx = ebdata.get<EBData_t::bndrycent>(i,j,k,0);
851  Real bcty = ebdata.get<EBData_t::bndrycent>(i,j,k,1);
852  Real bctz = ebdata.get<EBData_t::bndrycent>(i,j,k,2);
853  Real dx_eb = get_dx_eb(kappa);
854 
855  Real dg = dx_eb / amrex::max(std::abs(anrmx),std::abs(anrmy),
856  std::abs(anrmz));
857 
858  Real gx = bctx - dg*anrmx;
859  Real gy = bcty - dg*anrmy;
860  Real gz = bctz - dg*anrmz;
861  Real sx = std::copysign(Real(1.0),anrmx);
862  Real sy = std::copysign(Real(1.0),anrmy);
863  Real sz = std::copysign(Real(1.0),anrmz);
864  int ii = i - static_cast<int>(sx);
865  int jj = j - static_cast<int>(sy);
866  int kk = k - static_cast<int>(sz);
867 
868  gx *= sx;
869  gy *= sy;
870  gz *= sz;
871  Real gxy = gx*gy;
872  Real gxz = gx*gz;
873  Real gyz = gy*gz;
874  Real gxyz = gx*gy*gz;
875  Real phig_gamma = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz);
876  Real phig = (-gz - gxz - gyz - gxyz) * phi(i,j,kk,n)
877  + (-gy - gxy - gyz - gxyz) * phi(i,jj,k,n)
878  + (gyz + gxyz) * phi(i,jj,kk,n)
879  + (-gx - gxy - gxz - gxyz) * phi(ii,j,k,n)
880  + (gxz + gxyz) * phi(ii,j,kk,n)
881  + (gxy + gxyz) * phi(ii,jj,k,n)
882  + (-gxyz) * phi(ii,jj,kk,n);
883 
884  Real ba = ebdata.get<EBData_t::bndryarea>(i,j,k);
885 
886  Real dphidn = ( -phig)/dg;
887  Real feb_gamma = -phig_gamma/dg * ba * beb(i,j,k,n);
888  gamma += vfrcinv*(-dhx)*feb_gamma;
889  Real feb = dphidn * ba * beb(i,j,k,n);
890  rho += -vfrcinv*(-dhx)*feb;
891  }
892 
893  Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
894  phi(i,j,k,n) += omega*res/(gamma-delta);
895  }
896  }
897  }
898  }}}}
899 // });
900 }
901 
903 void mlebabeclap_flux_x (Box const& box, Array4<Real> const& fx, Array4<Real const> const& apx,
904  Array4<Real const> const& fcx, Array4<Real const> const& sol,
905  Array4<Real const> const& bX, Array4<int const> const& ccm,
906  Real dhx, int face_only, int ncomp, Box const& xbox,
907  bool beta_on_centroid, bool phi_on_centroid) noexcept
908 {
909  int lof = xbox.smallEnd(0);
910  int hif = xbox.bigEnd(0);
911  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
912  {
913  if (!face_only || lof == i || hif == i) {
914  if (apx(i,j,k) == Real(0.0)) {
915  fx(i,j,k,n) = Real(0.0);
916  } else if (apx(i,j,k) == Real(1.0)) {
917  fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
918  } else {
919  Real fxm = bX(i,j,k,n)*(sol(i,j,k,n) - sol(i-1,j,k,n));
920  int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
921  int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
922  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
923  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
924  if (!beta_on_centroid && !phi_on_centroid)
925  {
926  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*fxm +
927  fracy *(Real(1.0)-fracz)*bX(i,jj,k ,n)*(sol(i,jj,k ,n)-sol(i-1,jj,k ,n)) +
928  fracz *(Real(1.0)-fracy)*bX(i,j ,kk,n)*(sol(i,j ,kk,n)-sol(i-1,j ,kk,n)) +
929  fracy* fracz *bX(i,jj,kk,n)*(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
930  }
931  else if (beta_on_centroid && !phi_on_centroid)
932  {
933  fxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*(sol(i, j, k,n)-sol(i-1, j, k,n)) +
934  fracy *(Real(1.0)-fracz)*(sol(i,jj, k,n)-sol(i-1,jj, k,n)) +
935  fracz *(Real(1.0)-fracy)*(sol(i, j,kk,n)-sol(i-1, j,kk,n)) +
936  fracy* fracz *(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
937  fxm *= bX(i,j,k,n);
938 
939  }
940  fx(i,j,k,n) = -dhx*fxm;
941  }
942  }
943  });
944 }
945 
947 void mlebabeclap_flux_y (Box const& box, Array4<Real> const& fy, Array4<Real const> const& apy,
948  Array4<Real const> const& fcy, Array4<Real const> const& sol,
949  Array4<Real const> const& bY, Array4<int const> const& ccm,
950  Real dhy, int face_only, int ncomp, Box const& ybox,
951  bool beta_on_centroid, bool phi_on_centroid) noexcept
952 {
953  int lof = ybox.smallEnd(1);
954  int hif = ybox.bigEnd(1);
955  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
956  {
957  if (!face_only || lof == j || hif == j) {
958  if (apy(i,j,k) == Real(0.0)) {
959  fy(i,j,k,n) = Real(0.0);
960  } else if (apy(i,j,k) == Real(1.0)) {
961  fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
962  } else {
963  Real fym = bY(i,j,k,n)*(sol(i,j,k,n) - sol(i,j-1,k,n));
964  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
965  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
966  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
967  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
968  if (!beta_on_centroid && !phi_on_centroid)
969  {
970  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*fym +
971  fracx *(Real(1.0)-fracz)*bY(ii,j,k ,n)*(sol(ii,j,k ,n)-sol(ii,j-1,k ,n)) +
972  fracz *(Real(1.0)-fracx)*bY(i ,j,kk,n)*(sol(i ,j,kk,n)-sol(i ,j-1,kk,n)) +
973  fracx * fracz *bY(ii,j,kk,n)*(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
974  }
975  else if (beta_on_centroid && !phi_on_centroid)
976  {
977  fym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*(sol( i,j, k,n)-sol( i,j-1, k,n)) +
978  fracx *(Real(1.0)-fracz)*(sol(ii,j, k,n)-sol(ii,j-1, k,n)) +
979  fracz *(Real(1.0)-fracx)*(sol( i,j,kk,n)-sol( i,j-1,kk,n)) +
980  fracx * fracz *(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
981  fym *= bY(i,j,k,n);
982  }
983 
984  fy(i,j,k,n) = -dhy*fym;
985  }
986  }
987  });
988 }
989 
991 void mlebabeclap_flux_z (Box const& box, Array4<Real> const& fz, Array4<Real const> const& apz,
992  Array4<Real const> const& fcz, Array4<Real const> const& sol,
993  Array4<Real const> const& bZ, Array4<int const> const& ccm,
994  Real dhz, int face_only, int ncomp, Box const& zbox,
995  bool beta_on_centroid, bool phi_on_centroid) noexcept
996 {
997  int lof = zbox.smallEnd(2);
998  int hif = zbox.bigEnd(2);
999  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1000  {
1001  if (!face_only || lof == k || hif == k) {
1002  if (apz(i,j,k) == Real(0.0)) {
1003  fz(i,j,k,n) = Real(0.0);
1004  } else if (apz(i,j,k) == Real(1.0)) {
1005  fz(i,j,k,n) = -dhz*bZ(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
1006  } else {
1007  Real fzm = bZ(i,j,k,n)*(sol(i,j,k,n) - sol(i,j,k-1,n));
1008  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1009  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1010  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1011  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1012  if (!beta_on_centroid && !phi_on_centroid)
1013  {
1014  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*fzm +
1015  fracx*(Real(1.0)-fracy)*bZ(ii,j ,k,n)*(sol(ii,j ,k,n)-sol(ii,j ,k-1,n)) +
1016  fracy*(Real(1.0)-fracx)*bZ(i ,jj,k,n)*(sol(i ,jj,k,n)-sol(i ,jj,k-1,n)) +
1017  fracx* fracy *bZ(ii,jj,k,n)*(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1018  }
1019  else if (beta_on_centroid && !phi_on_centroid)
1020  {
1021  fzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*(sol( i, j,k,n)-sol( i, j,k-1,n)) +
1022  fracx *(Real(1.0)-fracy)*(sol(ii, j,k,n)-sol(ii, j,k-1,n)) +
1023  fracy *(Real(1.0)-fracx)*(sol( i,jj,k,n)-sol( i,jj,k-1,n)) +
1024  fracx * fracy *(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1025  fzm *= bZ(i,j,k,n);
1026 
1027  }
1028 
1029  fz(i,j,k,n) = -dhz*fzm;
1030  }
1031  }
1032  });
1033 }
1034 
1036 void mlebabeclap_flux_x_0 (Box const& box, Array4<Real> const& fx, Array4<Real const> const& apx,
1037  Array4<Real const> const& sol, Array4<Real const> const& bX,
1038  Real dhx, int face_only, int ncomp, Box const& xbox) noexcept
1039 {
1040  int lof = xbox.smallEnd(0);
1041  int hif = xbox.bigEnd(0);
1042  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1043  {
1044  if (!face_only || lof == i || hif == i) {
1045  if (apx(i,j,k) == Real(0.0)) {
1046  fx(i,j,k,n) = Real(0.0);
1047  } else {
1048  fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
1049  }
1050  }
1051  });
1052 }
1053 
1055 void mlebabeclap_flux_y_0 (Box const& box, Array4<Real> const& fy, Array4<Real const> const& apy,
1056  Array4<Real const> const& sol, Array4<Real const> const& bY,
1057  Real dhy, int face_only, int ncomp, Box const& ybox) noexcept
1058 {
1059  int lof = ybox.smallEnd(1);
1060  int hif = ybox.bigEnd(1);
1061  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1062  {
1063  if (!face_only || lof == j || hif == j) {
1064  if (apy(i,j,k) == Real(0.0)) {
1065  fy(i,j,k,n) = Real(0.0);
1066  } else {
1067  fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
1068  }
1069  }
1070  });
1071 }
1072 
1074 void mlebabeclap_flux_z_0 (Box const& box, Array4<Real> const& fz, Array4<Real const> const& apz,
1075  Array4<Real const> const& sol, Array4<Real const> const& bZ,
1076  Real dhz, int face_only, int ncomp, Box const& zbox) noexcept
1077 {
1078  int lof = zbox.smallEnd(2);
1079  int hif = zbox.bigEnd(2);
1080  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1081  {
1082  if (!face_only || lof == k || hif == k) {
1083  if (apz(i,j,k) == Real(0.0)) {
1084  fz(i,j,k,n) = Real(0.0);
1085  } else {
1086  fz(i,j,k,n) = -dhz*bZ(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
1087  }
1088  }
1089  });
1090 }
1091 
1093 void mlebabeclap_grad_x (Box const& box, Array4<Real> const& gx, Array4<Real const> const& sol,
1094  Array4<Real const> const& apx, Array4<Real const> const& fcx,
1095  Array4<int const> const& ccm,
1096  Real dxi, int ncomp, bool phi_on_centroid) noexcept
1097 {
1098  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1099  {
1100  if (apx(i,j,k) == Real(0.0)) {
1101  gx(i,j,k,n) = Real(0.0);
1102  } else if (apx(i,j,k) == Real(1.0)) {
1103  gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
1104  } else {
1105  Real gxm = (sol(i,j,k,n) - sol(i-1,j,k,n));
1106  int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
1107  int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
1108  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
1109  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
1110  if (!phi_on_centroid)
1111  {
1112  gxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*gxm +
1113  fracy*(Real(1.0)-fracz)*(sol(i,jj,k ,n)-sol(i-1,jj,k ,n)) +
1114  fracz*(Real(1.0)-fracy)*(sol(i,j ,kk,n)-sol(i-1,j ,kk,n)) +
1115  fracy* fracz *(sol(i,jj,kk,n)-sol(i-1,jj,kk,n));
1116  }
1117  gx(i,j,k,n) = dxi*gxm;
1118  }
1119  });
1120 }
1121 
1123 void mlebabeclap_grad_y (Box const& box, Array4<Real> const& gy, Array4<Real const> const& sol,
1124  Array4<Real const> const& apy, Array4<Real const> const& fcy,
1125  Array4<int const> const& ccm,
1126  Real dyi, int ncomp, bool phi_on_centroid) noexcept
1127 {
1128  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1129  {
1130  if (apy(i,j,k) == Real(0.0)) {
1131  gy(i,j,k,n) = Real(0.0);
1132  } else if (apy(i,j,k) == Real(1.0)) {
1133  gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
1134  } else {
1135  Real gym = (sol(i,j,k,n) - sol(i,j-1,k,n));
1136  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
1137  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
1138  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
1139  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
1140  if (!phi_on_centroid)
1141  {
1142  gym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*gym +
1143  fracx*(Real(1.0)-fracz)*(sol(ii,j,k ,n)-sol(ii,j-1,k ,n)) +
1144  fracz*(Real(1.0)-fracx)*(sol(i ,j,kk,n)-sol(i ,j-1,kk,n)) +
1145  fracx* fracz *(sol(ii,j,kk,n)-sol(ii,j-1,kk,n));
1146  }
1147  gy(i,j,k,n) = dyi*gym;
1148  }
1149  });
1150 }
1151 
1153 void mlebabeclap_grad_z (Box const& box, Array4<Real> const& gz, Array4<Real const> const& sol,
1154  Array4<Real const> const& apz, Array4<Real const> const& fcz,
1155  Array4<int const> const& ccm,
1156  Real dzi, int ncomp, bool phi_on_centroid) noexcept
1157 {
1158  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1159  {
1160  if (apz(i,j,k) == Real(0.0)) {
1161  gz(i,j,k,n) = Real(0.0);
1162  } else if (apz(i,j,k) == Real(1.0)) {
1163  gz(i,j,k,n) = dzi*(sol(i,j,k,n)-sol(i,j,k-1,n));
1164  } else {
1165  Real gzm = (sol(i,j,k,n) - sol(i,j,k-1,n));
1166  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1167  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1168  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1169  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1170  if (!phi_on_centroid)
1171  {
1172  gzm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*gzm +
1173  fracx*(Real(1.0)-fracy)*(sol(ii,j ,k,n)-sol(ii,j ,k-1,n)) +
1174  fracy*(Real(1.0)-fracx)*(sol(i ,jj,k,n)-sol(i ,jj,k-1,n)) +
1175  fracx* fracy *(sol(ii,jj,k,n)-sol(ii,jj,k-1,n));
1176  }
1177  gz(i,j,k,n) = dzi*gzm;
1178  }
1179  });
1180 }
1181 
1183 void mlebabeclap_grad_x_0 (Box const& box, Array4<Real> const& gx, Array4<Real const> const& sol,
1184  Array4<Real const> const& apx, Real dxi, int ncomp) noexcept
1185 {
1186  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1187  {
1188  if (apx(i,j,k) == Real(0.0)) {
1189  gx(i,j,k,n) = Real(0.0);
1190  } else {
1191  gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
1192  }
1193  });
1194 }
1195 
1197 void mlebabeclap_grad_y_0 (Box const& box, Array4<Real> const& gy, Array4<Real const> const& sol,
1198  Array4<Real const> const& apy, Real dyi, int ncomp) noexcept
1199 {
1200  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1201  {
1202  if (apy(i,j,k) == Real(0.0)) {
1203  gy(i,j,k,n) = Real(0.0);
1204  } else {
1205  gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
1206  }
1207  });
1208 }
1209 
1211 void mlebabeclap_grad_z_0 (Box const& box, Array4<Real> const& gz, Array4<Real const> const& sol,
1212  Array4<Real const> const& apz, Real dzi, int ncomp) noexcept
1213 {
1214  amrex::LoopConcurrent(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1215  {
1216  if (apz(i,j,k) == Real(0.0)) {
1217  gz(i,j,k,n) = Real(0.0);
1218  } else {
1219  gz(i,j,k,n) = dzi*(sol(i,j,k,n)-sol(i,j,k-1,n));
1220  }
1221  });
1222 }
1223 
1225 void mlebabeclap_normalize (Box const& box, Array4<Real> const& phi,
1226  Real alpha, Array4<Real const> const& a,
1227  Real dhx, Real dhy, Real dhz,
1228  Array4<Real const> const& bX, Array4<Real const> const& bY,
1229  Array4<Real const> const& bZ,
1230  Array4<const int> const& ccm, Array4<EBCellFlag const> const& flag,
1231  Array4<Real const> const& vfrc,
1232  Array4<Real const> const& apx, Array4<Real const> const& apy,
1233  Array4<Real const> const& apz,
1234  Array4<Real const> const& fcx, Array4<Real const> const& fcy,
1235  Array4<Real const> const& fcz,
1236  Array4<Real const> const& ba, Array4<Real const> const& bc,
1237  Array4<Real const> const& beb,
1238  bool is_dirichlet, bool beta_on_centroid, int ncomp) noexcept
1239 {
1240  amrex::Loop(box, ncomp, [=] (int i, int j, int k, int n) noexcept
1241  {
1242  if (flag(i,j,k).isRegular())
1243  {
1244  phi(i,j,k,n) /= alpha*a(i,j,k) + dhx*(bX(i,j,k,n) + bX(i+1,j,k,n))
1245  + dhy*(bY(i,j,k,n) + bY(i,j+1,k,n))
1246  + dhz*(bZ(i,j,k,n) + bZ(i,j,k+1,n));
1247  }
1248  else if (flag(i,j,k).isSingleValued())
1249  {
1250  Real kappa = vfrc(i,j,k);
1251  Real apxm = apx(i,j,k);
1252  Real apxp = apx(i+1,j,k);
1253  Real apym = apy(i,j,k);
1254  Real apyp = apy(i,j+1,k);
1255  Real apzm = apz(i,j,k);
1256  Real apzp = apz(i,j,k+1);
1257 
1258  Real sxm = bX(i,j,k,n);
1259  if (apxm != Real(0.0) && apxm != Real(1.0) && !beta_on_centroid) {
1260  int jj = j + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
1261  int kk = k + static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
1262  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
1263  ? std::abs(fcx(i,j,k,0)) : Real(0.0);
1264  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk))
1265  ? std::abs(fcx(i,j,k,1)) : Real(0.0);
1266  sxm = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxm;
1267  }
1268 
1269  Real sxp = -bX(i+1,j,k,n);
1270  if (apxp != Real(0.0) && apxp != Real(1.0) && !beta_on_centroid) {
1271  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,0)));
1272  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k,1)));
1273  Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
1274  ? std::abs(fcx(i+1,j,k,0)) : Real(0.0);
1275  Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk))
1276  ? std::abs(fcx(i+1,j,k,1)) : Real(0.0);
1277  sxp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*sxp;
1278  }
1279 
1280  Real sym = bY(i,j,k,n);
1281  if (apym != Real(0.0) && apym != Real(1.0) && !beta_on_centroid) {
1282  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,0)));
1283  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k,1)));
1284  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
1285  ? std::abs(fcy(i,j,k,0)) : Real(0.0);
1286  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk))
1287  ? std::abs(fcy(i,j,k,1)) : Real(0.0);
1288  sym = (Real(1.0)-fracx)*(Real(1.0)-fracz)*sym;
1289  }
1290 
1291  Real syp = -bY(i,j+1,k,n);
1292  if (apyp != Real(0.0) && apyp != Real(1.0) && !beta_on_centroid) {
1293  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,0)));
1294  int kk = k + static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k,1)));
1295  Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
1296  ? std::abs(fcy(i,j+1,k,0)) : Real(0.0);
1297  Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk))
1298  ? std::abs(fcy(i,j+1,k,1)) : Real(0.0);
1299  syp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*syp;
1300  }
1301 
1302  Real szm = bZ(i,j,k,n);
1303  if (apzm != Real(0.0) && apzm != Real(1.0) && !beta_on_centroid) {
1304  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,0)));
1305  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k,1)));
1306  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k))
1307  ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1308  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k))
1309  ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1310  szm = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szm;
1311  }
1312 
1313  Real szp = -bZ(i,j,k+1,n);
1314  if (apzp != Real(0.0) && apzp != Real(1.0) && !beta_on_centroid) {
1315  int ii = i + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,0)));
1316  int jj = j + static_cast<int>(std::copysign(Real(1.0),fcz(i,j,k+1,1)));
1317  Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1))
1318  ? std::abs(fcz(i,j,k+1,0)) : Real(0.0);
1319  Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1))
1320  ? std::abs(fcz(i,j,k+1,1)) : Real(0.0);
1321  szp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*szp;
1322  }
1323 
1324  Real vfrcinv = Real(1.0)/kappa;
1325  Real gamma = alpha*a(i,j,k) + vfrcinv *
1326  (dhx*(apxm*sxm-apxp*sxp) +
1327  dhy*(apym*sym-apyp*syp) +
1328  dhz*(apzm*szm-apzp*szp));
1329 
1330  if (is_dirichlet) {
1331  Real dapx = apxm-apxp;
1332  Real dapy = apym-apyp;
1333  Real dapz = apzm-apzp;
1334  Real anorm = std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
1335  Real anorminv = Real(1.0)/anorm;
1336  Real anrmx = dapx * anorminv;
1337  Real anrmy = dapy * anorminv;
1338  Real anrmz = dapz * anorminv;
1339  Real bctx = bc(i,j,k,0);
1340  Real bcty = bc(i,j,k,1);
1341  Real bctz = bc(i,j,k,2);
1342  Real dx_eb = get_dx_eb(vfrc(i,j,k));
1343 
1344  Real dg = dx_eb / amrex::max(std::abs(anrmx),std::abs(anrmy),
1345  std::abs(anrmz));
1346 
1347  Real gx = bctx - dg*anrmx;
1348  Real gy = bcty - dg*anrmy;
1349  Real gz = bctz - dg*anrmz;
1350  Real sx = std::copysign(Real(1.0),anrmx);
1351  Real sy = std::copysign(Real(1.0),anrmy);
1352  Real sz = std::copysign(Real(1.0),anrmz);
1353 
1354  gx *= sx;
1355  gy *= sy;
1356  gz *= sz;
1357  Real gxy = gx*gy;
1358  Real gxz = gx*gz;
1359  Real gyz = gy*gz;
1360  Real gxyz = gx*gy*gz;
1361  Real phig_gamma = (Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz);
1362  Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
1363  gamma += vfrcinv*(-dhx)*feb_gamma;
1364  }
1365 
1366  phi(i,j,k,n) /= gamma;
1367  }
1368  });
1369 }
1370 
1371 }
1372 
1373 #endif
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_x_0(Box const &box, Array4< Real > const &gx, Array4< Real const > const &sol, Array4< Real const > const &apx, Real dxi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:739
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:331
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_y_0(Box const &box, Array4< Real > const &gy, Array4< Real const > const &sol, Array4< Real const > const &apy, Real dyi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:753
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_adotx_centroid(Box const &box, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &a, Array4< Real const > const &bX, Array4< Real const > const &bY, 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, Array4< Real const > const &ccent, Array4< Real const > const &ba, Array4< Real const > const &bcent, Array4< Real const > const &beb, Array4< Real const > const &phieb, const int &domlo_x, const int &domlo_y, const int &domhi_x, const int &domhi_y, const bool &on_x_face, const bool &on_y_face, bool is_eb_dirichlet, bool is_eb_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Real alpha, Real beta, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:11
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_adotx(Box const &box, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &a, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< const int > 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, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, Array4< Real const > const &phieb, bool is_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Real alpha, Real beta, int ncomp, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:166
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_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_x_0(Box const &box, Array4< Real > const &fx, Array4< Real const > const &apx, Array4< Real const > const &sol, Array4< Real const > const &bX, Real dhx, int face_only, int ncomp, Box const &xbox) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:653
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_gsrb(Box const &box, Array4< Real > const &phi, Array4< Real const > const &rhs, Real alpha, Array4< Real const > const &a, Real dhx, Real dhy, Real dh, GpuArray< Real, AMREX_SPACEDIM > const &dx, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< int const > const &m0, Array4< int const > const &m2, Array4< int const > const &m1, Array4< int const > const &m3, Array4< Real const > const &f0, Array4< Real const > const &f2, Array4< Real const > const &f1, Array4< Real const > const &f3, Array4< const int > const &ccm, Array4< Real const > const &beb, EBData const &ebdata, bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid, Box const &vbox, int redblack, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:371
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_z(Box const &box, Array4< Real > const &gz, Array4< Real const > const &sol, Array4< Real const > const &apz, Array4< Real const > const &fcz, Array4< int const > const &ccm, Real dzi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:1153
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_x(Box const &box, Array4< Real > const &gx, Array4< Real const > const &sol, Array4< Real const > const &apx, Array4< Real const > const &fcx, Array4< int const > const &ccm, Real dxi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:691
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:308
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_x_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:495
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_z_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &xloc_on_zface, Real &yloc_on_zface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_3D_K.H:513
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_z_0(Box const &box, Array4< Real > const &fz, Array4< Real const > const &apz, Array4< Real const > const &sol, Array4< Real const > const &bZ, Real dhz, int face_only, int ncomp, Box const &zbox) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:1074
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_y_0(Box const &box, Array4< Real > const &fy, Array4< Real const > const &apy, Array4< Real const > const &sol, Array4< Real const > const &bY, Real dhy, int face_only, int ncomp, Box const &ybox) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:672
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_z_0(Box const &box, Array4< Real > const &gz, Array4< Real const > const &sol, Array4< Real const > const &apz, Real dzi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:1211
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_eb_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &nrmx, Real &nrmy, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:765
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_z_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &xloc_on_zface, Real &yloc_on_zface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y, const bool on_z_face, const int domlo_z, const int domhi_z)
Definition: AMReX_EB_LeastSquares_3D_K.H:1010
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_x(Box const &box, Array4< Real > const &fx, Array4< Real const > const &apx, Array4< Real const > const &fcx, Array4< Real const > const &sol, Array4< Real const > const &bX, Array4< int const > const &ccm, Real dhx, int face_only, int ncomp, Box const &xbox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:587
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:632
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void LoopConcurrent(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:150
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_ebflux(int i, int j, int k, int n, Array4< Real > const &feb, Array4< Real const > const &x, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &bc, Array4< Real const > const &beb, Array4< Real const > const &phieb, bool is_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:304
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_normalize(Box const &box, Array4< Real > const &phi, Real alpha, Array4< Real const > const &a, Real dhx, Real dhy, Real dh, const amrex::GpuArray< Real, AMREX_SPACEDIM > &dx, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< const int > 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, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, bool beta_on_centroid, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:767
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_z(Box const &box, Array4< Real > const &fz, Array4< Real const > const &apz, Array4< Real const > const &fcz, Array4< Real const > const &sol, Array4< Real const > const &bZ, Array4< int const > const &ccm, Real dhz, int face_only, int ncomp, Box const &zbox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_3D_K.H:991
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_y(Box const &box, Array4< Real > const &gy, Array4< Real const > const &sol, Array4< Real const > const &apy, Array4< Real const > const &fcy, Array4< int const > const &ccm, Real dyi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:715
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_x_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:246
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_y(Box const &box, Array4< Real > const &fy, Array4< Real const > const &apy, Array4< Real const > const &fcy, Array4< Real const > const &sol, Array4< Real const > const &bY, Array4< int const > const &ccm, Real dhy, int face_only, int ncomp, Box const &ybox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:620
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void Loop(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:125
Definition: AMReX_Array4.H:61
Definition: AMReX_EBData.H:26