Block-Structured AMR Software Framework
AMReX_MLEBTensor_3D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_ML_EB_TENSOR_3D_K_H_
2 #define AMREX_ML_EB_TENSOR_3D_K_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_MLEBABecLap_K.H>
6 
7 namespace amrex {
8 
9 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
10 Real mlebtensor_dz_on_xface (int i, int j, int, int n,
11  Array4<Real const> const& vel, Real dzi,
12  Real whi, Real wlo,
13  int khip, int khim, int klop, int klom) noexcept
14 {
15  return Real(0.5)*dzi * ((vel(i ,j,khip,n)-vel(i ,j,khim,n))*whi +
16  (vel(i-1,j,klop,n)-vel(i-1,j,klom,n))*wlo);
17 }
18 
19 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
20 Real mlebtensor_dz_on_yface (int i, int j, int, int n,
21  Array4<Real const> const& vel, Real dzi,
22  Real whi, Real wlo,
23  int khip, int khim, int klop, int klom) noexcept
24 {
25  return Real(0.5)*dzi * ((vel(i,j ,khip,n)-vel(i,j ,khim,n))*whi +
26  (vel(i,j-1,klop,n)-vel(i,j-1,klom,n))*wlo);
27 }
28 
29 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
30 Real mlebtensor_dx_on_zface (int, int j, int k, int n,
31  Array4<Real const> const& vel, Real dxi,
32  Real whi, Real wlo,
33  int ihip, int ihim, int ilop, int ilom) noexcept
34 {
35  return Real(0.5)*dxi * ((vel(ihip,j,k ,n)-vel(ihim,j,k ,n))*whi +
36  (vel(ilop,j,k-1,n)-vel(ilom,j,k-1,n))*wlo);
37 }
38 
39 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
40 Real mlebtensor_dy_on_zface (int i, int, int k, int n,
41  Array4<Real const> const& vel, Real dyi,
42  Real whi, Real wlo,
43  int jhip, int jhim, int jlop, int jlom) noexcept
44 {
45  return Real(0.5)*dyi * ((vel(i,jhip,k ,n)-vel(i,jhim,k ,n))*whi +
46  (vel(i,jlop,k-1,n)-vel(i,jlom,k-1,n))*wlo);
47 }
48 
49 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
50 void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
51  Array4<Real const> const& vel,
52  Array4<Real const> const& etax,
53  Array4<Real const> const& kapx,
54  Array4<Real const> const& apx,
55  Array4<EBCellFlag const> const& flag,
56  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
57 {
58  const Real dyi = dxinv[1];
59  const Real dzi = dxinv[2];
60  const auto lo = amrex::lbound(box);
61  const auto hi = amrex::ubound(box);
62  constexpr Real twoThirds = 2./3.;
63 
64  for (int k = lo.z; k <= hi.z; ++k) {
65  for (int j = lo.y; j <= hi.y; ++j) {
66  AMREX_PRAGMA_SIMD
67  for (int i = lo.x; i <= hi.x; ++i) {
68  if (apx(i,j,k) == Real(0.0))
69  {
70  fx(i,j,k,0) = Real(0.0);
71  fx(i,j,k,1) = Real(0.0);
72  fx(i,j,k,2) = Real(0.0);
73  }
74  else
75  {
76  int jhip = j + flag(i ,j,k).isConnected(0, 1,0);
77  int jhim = j - flag(i ,j,k).isConnected(0,-1,0);
78  int jlop = j + flag(i-1,j,k).isConnected(0, 1,0);
79  int jlom = j - flag(i-1,j,k).isConnected(0,-1,0);
80  Real whi = mlebtensor_weight(jhip-jhim);
81  Real wlo = mlebtensor_weight(jlop-jlom);
82  Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
83  whi,wlo,jhip,jhim,jlop,jlom);
84  Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
85  whi,wlo,jhip,jhim,jlop,jlom);
86  int khip = k + flag(i ,j,k).isConnected(0,0, 1);
87  int khim = k - flag(i ,j,k).isConnected(0,0,-1);
88  int klop = k + flag(i-1,j,k).isConnected(0,0, 1);
89  int klom = k - flag(i-1,j,k).isConnected(0,0,-1);
90  whi = mlebtensor_weight(khip-khim);
91  wlo = mlebtensor_weight(klop-klom);
92  Real dudz = mlebtensor_dz_on_xface(i,j,k,0,vel,dzi,
93  whi,wlo,khip,khim,klop,klom);
94  Real dwdz = mlebtensor_dz_on_xface(i,j,k,2,vel,dzi,
95  whi,wlo,khip,khim,klop,klom);
96  Real divu = dvdy + dwdz;
97  Real xif = kapx(i,j,k);
98  Real mun = Real(0.75)*(etax(i,j,k,0)-xif);// restore the original eta
99  Real mut = etax(i,j,k,1);
100  fx(i,j,k,0) = -mun*(-twoThirds*divu) - xif*divu;
101  fx(i,j,k,1) = -mut*dudy;
102  fx(i,j,k,2) = -mut*dudz;
103  }
104  }
105  }
106  }
107 }
108 
109 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
110 void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
111  Array4<Real const> const& vel,
112  Array4<Real const> const& etay,
113  Array4<Real const> const& kapy,
114  Array4<Real const> const& apy,
115  Array4<EBCellFlag const> const& flag,
116  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
117 {
118  const Real dxi = dxinv[0];
119  const Real dzi = dxinv[2];
120  const auto lo = amrex::lbound(box);
121  const auto hi = amrex::ubound(box);
122  constexpr Real twoThirds = 2./3.;
123 
124  for (int k = lo.z; k <= hi.z; ++k) {
125  for (int j = lo.y; j <= hi.y; ++j) {
126  AMREX_PRAGMA_SIMD
127  for (int i = lo.x; i <= hi.x; ++i) {
128  if (apy(i,j,k) == Real(0.0))
129  {
130  fy(i,j,k,0) = Real(0.0);
131  fy(i,j,k,1) = Real(0.0);
132  fy(i,j,k,2) = Real(0.0);
133  }
134  else
135  {
136  int ihip = i + flag(i,j ,k).isConnected( 1,0,0);
137  int ihim = i - flag(i,j ,k).isConnected(-1,0,0);
138  int ilop = i + flag(i,j-1,k).isConnected( 1,0,0);
139  int ilom = i - flag(i,j-1,k).isConnected(-1,0,0);
140  Real whi = mlebtensor_weight(ihip-ihim);
141  Real wlo = mlebtensor_weight(ilop-ilom);
142  Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
143  whi,wlo,ihip,ihim,ilop,ilom);
144  Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
145  whi,wlo,ihip,ihim,ilop,ilom);
146  int khip = k + flag(i,j ,k).isConnected(0,0, 1);
147  int khim = k - flag(i,j ,k).isConnected(0,0,-1);
148  int klop = k + flag(i,j-1,k).isConnected(0,0, 1);
149  int klom = k - flag(i,j-1,k).isConnected(0,0,-1);
150  whi = mlebtensor_weight(khip-khim);
151  wlo = mlebtensor_weight(klop-klom);
152  Real dvdz = mlebtensor_dz_on_yface(i,j,k,1,vel,dzi,
153  whi,wlo,khip,khim,klop,klom);
154  Real dwdz = mlebtensor_dz_on_yface(i,j,k,2,vel,dzi,
155  whi,wlo,khip,khim,klop,klom);
156  Real divu = dudx + dwdz;
157  Real xif = kapy(i,j,k);
158  Real mun = Real(0.75)*(etay(i,j,k,1)-xif);// restore the original eta
159  Real mut = etay(i,j,k,0);
160  fy(i,j,k,0) = -mut*dvdx;
161  fy(i,j,k,1) = -mun*(-twoThirds*divu) - xif*divu;
162  fy(i,j,k,2) = -mut*dvdz;
163  }
164  }
165  }
166  }
167 }
168 
169 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
170 void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
171  Array4<Real const> const& vel,
172  Array4<Real const> const& etaz,
173  Array4<Real const> const& kapz,
174  Array4<Real const> const& apz,
175  Array4<EBCellFlag const> const& flag,
176  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
177 {
178  const Real dxi = dxinv[0];
179  const Real dyi = dxinv[1];
180  const auto lo = amrex::lbound(box);
181  const auto hi = amrex::ubound(box);
182  constexpr Real twoThirds = 2./3.;
183 
184  for (int k = lo.z; k <= hi.z; ++k) {
185  for (int j = lo.y; j <= hi.y; ++j) {
186  AMREX_PRAGMA_SIMD
187  for (int i = lo.x; i <= hi.x; ++i) {
188  if (apz(i,j,k) == Real(0.0))
189  {
190  fz(i,j,k,0) = Real(0.0);
191  fz(i,j,k,1) = Real(0.0);
192  fz(i,j,k,2) = Real(0.0);
193  }
194  else
195  {
196  int ihip = i + flag(i,j,k ).isConnected( 1,0,0);
197  int ihim = i - flag(i,j,k ).isConnected(-1,0,0);
198  int ilop = i + flag(i,j,k-1).isConnected( 1,0,0);
199  int ilom = i - flag(i,j,k-1).isConnected(-1,0,0);
200  Real whi = mlebtensor_weight(ihip-ihim);
201  Real wlo = mlebtensor_weight(ilop-ilom);
202  Real dudx = mlebtensor_dx_on_zface(i,j,k,0,vel,dxi,
203  whi,wlo,ihip,ihim,ilop,ilom);
204  Real dwdx = mlebtensor_dx_on_zface(i,j,k,2,vel,dxi,
205  whi,wlo,ihip,ihim,ilop,ilom);
206  int jhip = j + flag(i,j,k ).isConnected(0, 1,0);
207  int jhim = j - flag(i,j,k ).isConnected(0,-1,0);
208  int jlop = j + flag(i,j,k-1).isConnected(0, 1,0);
209  int jlom = j - flag(i,j,k-1).isConnected(0,-1,0);
210  whi = mlebtensor_weight(jhip-jhim);
211  wlo = mlebtensor_weight(jlop-jlom);
212  Real dvdy = mlebtensor_dy_on_zface(i,j,k,1,vel,dyi,
213  whi,wlo,jhip,jhim,jlop,jlom);
214  Real dwdy = mlebtensor_dy_on_zface(i,j,k,2,vel,dyi,
215  whi,wlo,jhip,jhim,jlop,jlom);
216  Real divu = dudx + dvdy;
217  Real xif = kapz(i,j,k);
218  Real mun = Real(0.75)*(etaz(i,j,k,2)-xif);// restore the original eta
219  Real mut = etaz(i,j,k,0);
220 
221  fz(i,j,k,0) = -mut*dwdx;
222  fz(i,j,k,1) = -mut*dwdy;
223  fz(i,j,k,2) = -mun*(-twoThirds*divu) - xif*divu;
224  }
225  }
226  }
227  }
228 }
229 
230 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
231 Real mlebtensor_dz_on_xface (int i, int j, int k, int n,
232  Array4<Real const> const& vel, Real dzi,
233  Array4<Real const> const& bvxlo,
234  Array4<Real const> const& bvxhi,
235  Array2D<BoundCond,
236  0,2*AMREX_SPACEDIM,
237  0,AMREX_SPACEDIM> const& bct,
238  Dim3 const& dlo, Dim3 const& dhi,
239  Real whi, Real wlo,
240  int khip, int khim, int klop, int klom) noexcept
241 {
242  Real ddz;
243  if (i == dlo.x) {
244  if (bct(Orientation::xlo(),n) == AMREX_LO_DIRICHLET && bvxlo) {
245  if (k == dlo.z) {
246  ddz = (bvxlo(i-1,j,k ,n) * Real(-1.5) +
247  bvxlo(i-1,j,k+1,n) * Real(2.) +
248  bvxlo(i-1,j,k+2,n) * Real(-0.5)) * dzi;
249  } else if (k == dhi.z) {
250  ddz = -(bvxlo(i-1,j,k ,n) * Real(-1.5) +
251  bvxlo(i-1,j,k-1,n) * Real(2.) +
252  bvxlo(i-1,j,k-2,n) * Real(-0.5)) * dzi;
253  } else {
254  ddz = wlo*dzi*(bvxlo(i-1,j,klop,n)-bvxlo(i-1,j,klom,n));
255  }
256  } else if (bct(Orientation::xlo(),n) == AMREX_LO_NEUMANN) {
257  ddz = whi*dzi*(vel(i,j,khip,n)-vel(i,j,khim,n));
258  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
259  ddz = Real(0.);
260  }
261  } else if (i == dhi.x+1) {
262  if (bct(Orientation::xhi(),n) == AMREX_LO_DIRICHLET && bvxhi) {
263  if (k == dlo.z) {
264  ddz = (bvxhi(i,j,k ,n) * Real(-1.5) +
265  bvxhi(i,j,k+1,n) * Real(2.) +
266  bvxhi(i,j,k+2,n) * Real(-0.5)) * dzi;
267  } else if (k == dhi.z) {
268  ddz = -(bvxhi(i,j,k ,n) * Real(-1.5) +
269  bvxhi(i,j,k-1,n) * Real(2.) +
270  bvxhi(i,j,k-2,n) * Real(-0.5)) * dzi;
271  } else {
272  ddz = whi*dzi*(bvxhi(i,j,khip,n)-bvxhi(i,j,khim,n));
273  }
274  } else if (bct(Orientation::xhi(),n) == AMREX_LO_NEUMANN) {
275  ddz = wlo*dzi*(vel(i-1,j,klop,n)-vel(i-1,j,klom,n));
276  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
277  ddz = Real(0.);
278  }
279  } else {
280  ddz = mlebtensor_dz_on_xface(i,j,k,n,vel,dzi,whi,wlo,khip,khim,klop,klom);
281  }
282  return ddz;
283 }
284 
285 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
286 Real mlebtensor_dz_on_yface (int i, int j, int k, int n,
287  Array4<Real const> const& vel, Real dzi,
288  Array4<Real const> const& bvylo,
289  Array4<Real const> const& bvyhi,
290  Array2D<BoundCond,
291  0,2*AMREX_SPACEDIM,
292  0,AMREX_SPACEDIM> const& bct,
293  Dim3 const& dlo, Dim3 const& dhi,
294  Real whi, Real wlo,
295  int khip, int khim, int klop, int klom) noexcept
296 {
297  Real ddz;
298  if (j == dlo.y) {
299  if (bct(Orientation::ylo(),n) == AMREX_LO_DIRICHLET && bvylo) {
300  if (k == dlo.z) {
301  ddz = (bvylo(i,j-1,k ,n) * Real(-1.5) +
302  bvylo(i,j-1,k+1,n) * Real(2.) +
303  bvylo(i,j-1,k+2,n) * Real(-0.5)) * dzi;
304  } else if (k == dhi.z) {
305  ddz = -(bvylo(i,j-1,k ,n) * Real(-1.5) +
306  bvylo(i,j-1,k-1,n) * Real(2.) +
307  bvylo(i,j-1,k-2,n) * Real(-0.5)) * dzi;
308  } else {
309  ddz = wlo*dzi*(bvylo(i,j-1,klop,n)-bvylo(i,j-1,klom,n));
310  }
311  } else if (bct(Orientation::ylo(),n) == AMREX_LO_NEUMANN) {
312  ddz = whi*dzi*(vel(i,j,khip,n)-vel(i,j,khim,n));
313  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
314  ddz = Real(0.);
315  }
316  } else if (j == dhi.y+1) {
317  if (bct(Orientation::yhi(),n) == AMREX_LO_DIRICHLET && bvyhi) {
318  if (k == dlo.z) {
319  ddz = (bvyhi(i,j,k ,n) * Real(-1.5) +
320  bvyhi(i,j,k+1,n) * Real(2.) +
321  bvyhi(i,j,k+2,n) * Real(-0.5)) * dzi;
322  } else if (k == dhi.z) {
323  ddz = -(bvyhi(i,j,k ,n) * Real(-1.5) +
324  bvyhi(i,j,k-1,n) * Real(2.) +
325  bvyhi(i,j,k-2,n) * Real(-0.5)) * dzi;
326  } else {
327  ddz = whi*dzi*(bvyhi(i,j,khip,n)-bvyhi(i,j,khim,n));
328  }
329  } else if (bct(Orientation::yhi(),n) == AMREX_LO_NEUMANN) {
330  ddz = wlo*dzi*(vel(i,j-1,klop,n)-vel(i,j-1,klom,n));
331  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
332  ddz = Real(0.);
333  }
334  } else {
335  ddz = mlebtensor_dz_on_yface(i,j,k,n,vel,dzi,whi,wlo,khip,khim,klop,klom);
336  }
337  return ddz;
338 }
339 
340 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
341 Real mlebtensor_dx_on_zface (int i, int j, int k, int n,
342  Array4<Real const> const& vel, Real dxi,
343  Array4<Real const> const& bvzlo,
344  Array4<Real const> const& bvzhi,
345  Array2D<BoundCond,
346  0,2*AMREX_SPACEDIM,
347  0,AMREX_SPACEDIM> const& bct,
348  Dim3 const& dlo, Dim3 const& dhi,
349  Real whi, Real wlo,
350  int ihip, int ihim, int ilop, int ilom) noexcept
351 {
352  Real ddx;
353  if (k == dlo.z) {
354  if (bct(Orientation::zlo(),n) == AMREX_LO_DIRICHLET && bvzlo) {
355  if (i == dlo.x) {
356  ddx = (bvzlo(i ,j,k-1,n) * Real(-1.5) +
357  bvzlo(i+1,j,k-1,n) * Real(2.) +
358  bvzlo(i+2,j,k-1,n) * Real(-0.5)) * dxi;
359  } else if (i == dhi.x) {
360  ddx = -(bvzlo(i ,j,k-1,n) * Real(-1.5) +
361  bvzlo(i-1,j,k-1,n) * Real(2.) +
362  bvzlo(i-2,j,k-1,n) * Real(-0.5)) * dxi;
363  } else {
364  ddx = wlo*dxi*(bvzlo(ilop,j,k-1,n)-bvzlo(ilom,j,k-1,n));
365  }
366  } else if (bct(Orientation::zlo(),n) == AMREX_LO_NEUMANN) {
367  ddx = whi*dxi*(vel(ihip,j,k,n)-vel(ihim,j,k,n));
368  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
369  ddx = Real(0.);
370  }
371  } else if (k == dhi.z+1) {
372  if (bct(Orientation::zhi(),n) == AMREX_LO_DIRICHLET && bvzhi) {
373  if (i == dlo.x) {
374  ddx = (bvzhi(i ,j,k,n) * Real(-1.5) +
375  bvzhi(i+1,j,k,n) * Real(2.) +
376  bvzhi(i+2,j,k,n) * Real(-0.5)) * dxi;
377  } else if (i == dhi.x) {
378  ddx = -(bvzhi(i ,j,k,n) * Real(-1.5) +
379  bvzhi(i-1,j,k,n) * Real(2.) +
380  bvzhi(i-2,j,k,n) * Real(-0.5)) * dxi;
381  } else {
382  ddx = whi*dxi*(bvzhi(ihip,j,k,n)-bvzhi(ihim,j,k,n));
383  }
384  } else if (bct(Orientation::zhi(),n) == AMREX_LO_NEUMANN) {
385  ddx = wlo*dxi*(vel(ilop,j,k-1,n)-vel(ilom,j,k-1,n));
386  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
387  ddx = Real(0.);
388  }
389  } else {
390  ddx = mlebtensor_dx_on_zface(i,j,k,n,vel,dxi,whi,wlo,ihip,ihim,ilop,ilom);
391 
392  }
393  return ddx;
394 }
395 
396 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
397 Real mlebtensor_dy_on_zface (int i, int j, int k, int n,
398  Array4<Real const> const& vel, Real dyi,
399  Array4<Real const> const& bvzlo,
400  Array4<Real const> const& bvzhi,
401  Array2D<BoundCond,
402  0,2*AMREX_SPACEDIM,
403  0,AMREX_SPACEDIM> const& bct,
404  Dim3 const& dlo, Dim3 const& dhi,
405  Real whi, Real wlo,
406  int jhip, int jhim, int jlop, int jlom) noexcept
407 {
408  Real ddy;
409  if (k == dlo.z) {
410  if (bct(Orientation::zlo(),n) == AMREX_LO_DIRICHLET && bvzlo) {
411  if (j == dlo.y) {
412  ddy = (bvzlo(i,j ,k-1,n) * Real(-1.5) +
413  bvzlo(i,j+1,k-1,n) * Real(2.) +
414  bvzlo(i,j+2,k-1,n) * Real(-0.5)) * dyi;
415  } else if (j == dhi.y) {
416  ddy = -(bvzlo(i,j ,k-1,n) * Real(-1.5) +
417  bvzlo(i,j-1,k-1,n) * Real(2.) +
418  bvzlo(i,j-2,k-1,n) * Real(-0.5)) * dyi;
419  } else {
420  ddy = wlo*dyi*(bvzlo(i,jlop,k-1,n)-bvzlo(i,jlom,k-1,n));
421  }
422  } else if (bct(Orientation::zlo(),n) == AMREX_LO_NEUMANN) {
423  ddy = whi*dyi*(vel(i,jhip,k,n)-vel(i,jhim,k,n));
424  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
425  ddy = Real(0.);
426  }
427  } else if (k == dhi.z+1) {
428  if (bct(Orientation::zhi(),n) == AMREX_LO_DIRICHLET && bvzhi) {
429  if (j == dlo.y) {
430  ddy = (bvzhi(i,j ,k,n) * Real(-1.5) +
431  bvzhi(i,j+1,k,n) * Real(2.) +
432  bvzhi(i,j+2,k,n) * Real(-0.5)) * dyi;
433  } else if (j == dhi.y) {
434  ddy = -(bvzhi(i,j ,k,n) * Real(-1.5) +
435  bvzhi(i,j-1,k,n) * Real(2.) +
436  bvzhi(i,j-2,k,n) * Real(-0.5)) * dyi;
437  } else {
438  ddy = whi*dyi*(bvzhi(i,jhip,k,n)-bvzhi(i,jhim,k,n));
439  }
440  } else if (bct(Orientation::zhi(),n) == AMREX_LO_NEUMANN) {
441  ddy = wlo*dyi*(vel(i,jlop,k-1,n)-vel(i,jlom,k-1,n));
442  } else { // AMREX_LO_REFLECT_ODD or homogeneous Dirichlet
443  ddy = Real(0.);
444  }
445  } else {
446  ddy = mlebtensor_dy_on_zface(i,j,k,n,vel,dyi,whi,wlo,jhip,jhim,jlop,jlom);
447  }
448  return ddy;
449 }
450 
451 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
452 void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
453  Array4<Real const> const& vel,
454  Array4<Real const> const& etax,
455  Array4<Real const> const& kapx,
456  Array4<Real const> const& apx,
457  Array4<EBCellFlag const> const& flag,
458  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
459  Array4<Real const> const& bvxlo,
460  Array4<Real const> const& bvxhi,
461  Array2D<BoundCond,
462  0,2*AMREX_SPACEDIM,
463  0,AMREX_SPACEDIM> const& bct,
464  Dim3 const& dlo, Dim3 const& dhi) noexcept
465 
466 {
467  const Real dyi = dxinv[1];
468  const Real dzi = dxinv[2];
469  const auto lo = amrex::lbound(box);
470  const auto hi = amrex::ubound(box);
471  constexpr Real twoThirds = 2./3.;
472 
473  for (int k = lo.z; k <= hi.z; ++k) {
474  for (int j = lo.y; j <= hi.y; ++j) {
475  AMREX_PRAGMA_SIMD
476  for (int i = lo.x; i <= hi.x; ++i) {
477  if (apx(i,j,k) == Real(0.0))
478  {
479  fx(i,j,k,0) = Real(0.0);
480  fx(i,j,k,1) = Real(0.0);
481  fx(i,j,k,2) = Real(0.0);
482  }
483  else
484  {
485  int jhip = j + flag(i ,j,k).isConnected(0, 1,0);
486  int jhim = j - flag(i ,j,k).isConnected(0,-1,0);
487  int jlop = j + flag(i-1,j,k).isConnected(0, 1,0);
488  int jlom = j - flag(i-1,j,k).isConnected(0,-1,0);
489  Real whi = mlebtensor_weight(jhip-jhim);
490  Real wlo = mlebtensor_weight(jlop-jlom);
491  Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
492  bvxlo,bvxhi,bct,dlo,dhi,
493  whi,wlo,jhip,jhim,jlop,jlom);
494  Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
495  bvxlo,bvxhi,bct,dlo,dhi,
496  whi,wlo,jhip,jhim,jlop,jlom);
497  int khip = k + flag(i ,j,k).isConnected(0,0, 1);
498  int khim = k - flag(i ,j,k).isConnected(0,0,-1);
499  int klop = k + flag(i-1,j,k).isConnected(0,0, 1);
500  int klom = k - flag(i-1,j,k).isConnected(0,0,-1);
501  whi = mlebtensor_weight(khip-khim);
502  wlo = mlebtensor_weight(klop-klom);
503  Real dudz = mlebtensor_dz_on_xface(i,j,k,0,vel,dzi,
504  bvxlo,bvxhi,bct,dlo,dhi,
505  whi,wlo,khip,khim,klop,klom);
506  Real dwdz = mlebtensor_dz_on_xface(i,j,k,2,vel,dzi,
507  bvxlo,bvxhi,bct,dlo,dhi,
508  whi,wlo,khip,khim,klop,klom);
509  Real divu = dvdy + dwdz;
510  Real xif = kapx(i,j,k);
511  Real mun = Real(0.75)*(etax(i,j,k,0)-xif);// restore the original eta
512  Real mut = etax(i,j,k,1);
513  fx(i,j,k,0) = -mun*(-twoThirds*divu) - xif*divu;
514  fx(i,j,k,1) = -mut*dudy;
515  fx(i,j,k,2) = -mut*dudz;
516  }
517  }
518  }
519  }
520 }
521 
522 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
523 void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
524  Array4<Real const> const& vel,
525  Array4<Real const> const& etay,
526  Array4<Real const> const& kapy,
527  Array4<Real const> const& apy,
528  Array4<EBCellFlag const> const& flag,
529  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
530  Array4<Real const> const& bvylo,
531  Array4<Real const> const& bvyhi,
532  Array2D<BoundCond,
533  0,2*AMREX_SPACEDIM,
534  0,AMREX_SPACEDIM> const& bct,
535  Dim3 const& dlo, Dim3 const& dhi) noexcept
536 {
537  const Real dxi = dxinv[0];
538  const Real dzi = dxinv[2];
539  const auto lo = amrex::lbound(box);
540  const auto hi = amrex::ubound(box);
541  constexpr Real twoThirds = 2./3.;
542 
543  for (int k = lo.z; k <= hi.z; ++k) {
544  for (int j = lo.y; j <= hi.y; ++j) {
545  AMREX_PRAGMA_SIMD
546  for (int i = lo.x; i <= hi.x; ++i) {
547  if (apy(i,j,k) == Real(0.0))
548  {
549  fy(i,j,k,0) = Real(0.0);
550  fy(i,j,k,1) = Real(0.0);
551  fy(i,j,k,2) = Real(0.0);
552  }
553  else
554  {
555  int ihip = i + flag(i,j ,k).isConnected( 1,0,0);
556  int ihim = i - flag(i,j ,k).isConnected(-1,0,0);
557  int ilop = i + flag(i,j-1,k).isConnected( 1,0,0);
558  int ilom = i - flag(i,j-1,k).isConnected(-1,0,0);
559  Real whi = mlebtensor_weight(ihip-ihim);
560  Real wlo = mlebtensor_weight(ilop-ilom);
561  Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
562  bvylo,bvyhi,bct,dlo,dhi,
563  whi,wlo,ihip,ihim,ilop,ilom);
564  Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
565  bvylo,bvyhi,bct,dlo,dhi,
566  whi,wlo,ihip,ihim,ilop,ilom);
567  int khip = k + flag(i,j ,k).isConnected(0,0, 1);
568  int khim = k - flag(i,j ,k).isConnected(0,0,-1);
569  int klop = k + flag(i,j-1,k).isConnected(0,0, 1);
570  int klom = k - flag(i,j-1,k).isConnected(0,0,-1);
571  whi = mlebtensor_weight(khip-khim);
572  wlo = mlebtensor_weight(klop-klom);
573  Real dvdz = mlebtensor_dz_on_yface(i,j,k,1,vel,dzi,
574  bvylo,bvyhi,bct,dlo,dhi,
575  whi,wlo,khip,khim,klop,klom);
576  Real dwdz = mlebtensor_dz_on_yface(i,j,k,2,vel,dzi,
577  bvylo,bvyhi,bct,dlo,dhi,
578  whi,wlo,khip,khim,klop,klom);
579  Real divu = dudx + dwdz;
580  Real xif = kapy(i,j,k);
581  Real mun = Real(0.75)*(etay(i,j,k,1)-xif);// restore the original eta
582  Real mut = etay(i,j,k,0);
583  fy(i,j,k,0) = -mut*dvdx;
584  fy(i,j,k,1) = -mun*(-twoThirds*divu) - xif*divu;
585  fy(i,j,k,2) = -mut*dvdz;
586  }
587  }
588  }
589  }
590 }
591 
592 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
593 void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
594  Array4<Real const> const& vel,
595  Array4<Real const> const& etaz,
596  Array4<Real const> const& kapz,
597  Array4<Real const> const& apz,
598  Array4<EBCellFlag const> const& flag,
599  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
600  Array4<Real const> const& bvzlo,
601  Array4<Real const> const& bvzhi,
602  Array2D<BoundCond,
603  0,2*AMREX_SPACEDIM,
604  0,AMREX_SPACEDIM> const& bct,
605  Dim3 const& dlo, Dim3 const& dhi) noexcept
606 {
607  const Real dxi = dxinv[0];
608  const Real dyi = dxinv[1];
609  const auto lo = amrex::lbound(box);
610  const auto hi = amrex::ubound(box);
611  constexpr Real twoThirds = 2./3.;
612 
613  for (int k = lo.z; k <= hi.z; ++k) {
614  for (int j = lo.y; j <= hi.y; ++j) {
615  AMREX_PRAGMA_SIMD
616  for (int i = lo.x; i <= hi.x; ++i) {
617  if (apz(i,j,k) == Real(0.0))
618  {
619  fz(i,j,k,0) = Real(0.0);
620  fz(i,j,k,1) = Real(0.0);
621  fz(i,j,k,2) = Real(0.0);
622  }
623  else
624  {
625  int ihip = i + flag(i,j,k ).isConnected( 1,0,0);
626  int ihim = i - flag(i,j,k ).isConnected(-1,0,0);
627  int ilop = i + flag(i,j,k-1).isConnected( 1,0,0);
628  int ilom = i - flag(i,j,k-1).isConnected(-1,0,0);
629  Real whi = mlebtensor_weight(ihip-ihim);
630  Real wlo = mlebtensor_weight(ilop-ilom);
631  Real dudx = mlebtensor_dx_on_zface(i,j,k,0,vel,dxi,
632  bvzlo,bvzhi,bct,dlo,dhi,
633  whi,wlo,ihip,ihim,ilop,ilom);
634  Real dwdx = mlebtensor_dx_on_zface(i,j,k,2,vel,dxi,
635  bvzlo,bvzhi,bct,dlo,dhi,
636  whi,wlo,ihip,ihim,ilop,ilom);
637  int jhip = j + flag(i,j,k ).isConnected(0, 1,0);
638  int jhim = j - flag(i,j,k ).isConnected(0,-1,0);
639  int jlop = j + flag(i,j,k-1).isConnected(0, 1,0);
640  int jlom = j - flag(i,j,k-1).isConnected(0,-1,0);
641  whi = mlebtensor_weight(jhip-jhim);
642  wlo = mlebtensor_weight(jlop-jlom);
643  Real dvdy = mlebtensor_dy_on_zface(i,j,k,1,vel,dyi,
644  bvzlo,bvzhi,bct,dlo,dhi,
645  whi,wlo,jhip,jhim,jlop,jlom);
646  Real dwdy = mlebtensor_dy_on_zface(i,j,k,2,vel,dyi,
647  bvzlo,bvzhi,bct,dlo,dhi,
648  whi,wlo,jhip,jhim,jlop,jlom);
649  Real divu = dudx + dvdy;
650  Real xif = kapz(i,j,k);
651  Real mun = Real(0.75)*(etaz(i,j,k,2)-xif);// restore the original eta
652  Real mut = etaz(i,j,k,0);
653 
654  fz(i,j,k,0) = -mut*dwdx;
655  fz(i,j,k,1) = -mut*dwdy;
656  fz(i,j,k,2) = -mun*(-twoThirds*divu) - xif*divu;
657  }
658  }
659  }
660  }
661 }
662 
663 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
664 void mlebtensor_cross_terms (Box const& box, Array4<Real> const& Ax,
665  Array4<Real const> const& fx,
666  Array4<Real const> const& fy,
667  Array4<Real const> const& fz,
668  Array4<Real const> const& vel,
669  Array4<Real const> const& velb,
670  Array4<Real const> const& etab,
671  Array4<Real const> const& kapb,
672  Array4<int const> const& ccm,
673  Array4<EBCellFlag const> const& flag,
674  Array4<Real const> const& vol,
675  Array4<Real const> const& apx,
676  Array4<Real const> const& apy,
677  Array4<Real const> const& apz,
678  Array4<Real const> const& fcx,
679  Array4<Real const> const& fcy,
680  Array4<Real const> const& fcz,
681  Array4<Real const> const& bc,
682  bool is_dirichlet, bool is_inhomog,
683  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
684  Real bscalar) noexcept
685 {
686  const Real dxi = dxinv[0];
687  const Real dyi = dxinv[1];
688  const Real dzi = dxinv[2];
689  const auto lo = amrex::lbound(box);
690  const auto hi = amrex::ubound(box);
691 
692  for (int k = lo.z; k <= hi.z; ++k) {
693  for (int j = lo.y; j <= hi.y; ++j) {
694  AMREX_PRAGMA_SIMD
695  for (int i = lo.x; i <= hi.x; ++i) {
696  if (flag(i,j,k).isRegular())
697  {
698  Ax(i,j,k,0) += bscalar*(dxi*(fx(i+1,j ,k ,0) - fx(i,j,k,0))
699  + dyi*(fy(i ,j+1,k ,0) - fy(i,j,k,0))
700  + dzi*(fz(i ,j ,k+1,0) - fz(i,j,k,0)));
701  Ax(i,j,k,1) += bscalar*(dxi*(fx(i+1,j ,k ,1) - fx(i,j,k,1))
702  + dyi*(fy(i ,j+1,k ,1) - fy(i,j,k,1))
703  + dzi*(fz(i ,j ,k+1,1) - fz(i,j,k,1)));
704  Ax(i,j,k,2) += bscalar*(dxi*(fx(i+1,j ,k ,2) - fx(i,j,k,2))
705  + dyi*(fy(i ,j+1,k ,2) - fy(i,j,k,2))
706  + dzi*(fz(i ,j ,k+1,2) - fz(i,j,k,2)));
707  }
708  else if (flag(i,j,k).isSingleValued())
709  {
710  Real fxm_0 = fx(i,j,k,0);
711  Real fxm_1 = fx(i,j,k,1);
712  Real fxm_2 = fx(i,j,k,2);
713  if (apx(i,j,k) > Real(0.0) && apx(i,j,k) < Real(1.0)) {
714  int jj = j + (fcx(i,j,k,0) >= Real(0.0) ? 1 : -1);
715  int kk = k + (fcx(i,j,k,1) >= Real(0.0) ? 1 : -1);
716  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
717  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
718  fxm_0 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxm_0
719  + fracy*(Real(1.0)-fracz) * fx(i,jj,k ,0)
720  + fracz*(Real(1.0)-fracy) * fx(i,j ,kk,0)
721  + fracy* fracz * fx(i,jj,kk,0);
722  fxm_1 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxm_1
723  + fracy*(Real(1.0)-fracz) * fx(i,jj,k ,1)
724  + fracz*(Real(1.0)-fracy) * fx(i,j ,kk,1)
725  + fracy* fracz * fx(i,jj,kk,1);
726  fxm_2 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxm_2
727  + fracy*(Real(1.0)-fracz) * fx(i,jj,k ,2)
728  + fracz*(Real(1.0)-fracy) * fx(i,j ,kk,2)
729  + fracy* fracz * fx(i,jj,kk,2);
730  }
731 
732  Real fxp_0 = fx(i+1,j,k,0);
733  Real fxp_1 = fx(i+1,j,k,1);
734  Real fxp_2 = fx(i+1,j,k,2);
735  if (apx(i+1,j,k) > Real(0.0) && apx(i+1,j,k) < Real(1.0)) {
736  int jj = j + (fcx(i+1,j,k,0) >= Real(0.0) ? 1 : -1);
737  int kk = k + (fcx(i+1,j,k,1) >= Real(0.0) ? 1 : -1);
738  Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx(i+1,j,k,0)) : Real(0.0);
739  Real fracz = (ccm(i,j,kk) || ccm(i+1,j,kk)) ? std::abs(fcx(i+1,j,k,1)) : Real(0.0);
740  fxp_0 = (Real(1.0)-fracy)*(Real(1.0)-fracz) * fxp_0
741  + fracy*(Real(1.0)-fracz) * fx(i+1,jj,k ,0)
742  + fracz*(Real(1.0)-fracy) * fx(i+1,j ,kk,0)
743  + fracy* fracz * fx(i+1,jj,kk,0);
744  fxp_1 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxp_1
745  + fracy*(Real(1.0)-fracz) * fx(i+1,jj,k ,1)
746  + fracz*(Real(1.0)-fracy) * fx(i+1,j ,kk,1)
747  + fracy* fracz * fx(i+1,jj,kk,1);
748  fxp_2 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxp_2
749  + fracy*(Real(1.0)-fracz) * fx(i+1,jj,k ,2)
750  + fracz*(Real(1.0)-fracy) * fx(i+1,j ,kk,2)
751  + fracy* fracz * fx(i+1,jj,kk,2);
752  }
753 
754  Real fym_0 = fy(i,j,k,0);
755  Real fym_1 = fy(i,j,k,1);
756  Real fym_2 = fy(i,j,k,2);
757  if (apy(i,j,k) > Real(0.0) && apy(i,j,k) < Real(1.0)) {
758  int ii = i + (fcy(i,j,k,0) >= Real(0.0) ? 1 : -1);
759  int kk = k + (fcy(i,j,k,1) >= Real(0.0) ? 1 : -1);
760  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
761  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
762  fym_0 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fym_0
763  + fracx*(Real(1.0)-fracz) * fy(ii,j,k ,0)
764  + fracz*(Real(1.0)-fracx) * fy(i ,j,kk,0)
765  + fracx* fracz * fy(ii,j,kk,0);
766  fym_1 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fym_1
767  + fracx*(Real(1.0)-fracz) * fy(ii,j,k ,1)
768  + fracz*(Real(1.0)-fracx) * fy(i ,j,kk,1)
769  + fracx* fracz * fy(ii,j,kk,1);
770  fym_2 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fym_2
771  + fracx*(Real(1.0)-fracz) * fy(ii,j,k ,2)
772  + fracz*(Real(1.0)-fracx) * fy(i ,j,kk,2)
773  + fracx* fracz * fy(ii,j,kk,2);
774  }
775 
776  Real fyp_0 = fy(i,j+1,k,0);
777  Real fyp_1 = fy(i,j+1,k,1);
778  Real fyp_2 = fy(i,j+1,k,2);
779  if (apy(i,j+1,k) > Real(0.0) && apy(i,j+1,k) < Real(1.0)) {
780  int ii = i + (fcy(i,j+1,k,0) >= Real(0.0) ? 1 : -1);
781  int kk = k + (fcy(i,j+1,k,1) >= Real(0.0) ? 1 : -1);
782  Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy(i,j+1,k,0)) : Real(0.0);
783  Real fracz = (ccm(i,j,kk) || ccm(i,j+1,kk)) ? std::abs(fcy(i,j+1,k,1)) : Real(0.0);
784  fyp_0 = (Real(1.0)-fracx)*(Real(1.0)-fracz) * fyp_0
785  + fracx*(Real(1.0)-fracz) * fy(ii,j+1,k ,0)
786  + fracz*(Real(1.0)-fracx) * fy(i ,j+1,kk,0)
787  + fracx* fracz * fy(ii,j+1,kk,0);
788  fyp_1 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fyp_1
789  + fracx*(Real(1.0)-fracz) * fy(ii,j+1,k ,1)
790  + fracz*(Real(1.0)-fracx) * fy(i ,j+1,kk,1)
791  + fracx* fracz * fy(ii,j+1,kk,1);
792  fyp_2 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fyp_2
793  + fracx*(Real(1.0)-fracz) * fy(ii,j+1,k ,2)
794  + fracz*(Real(1.0)-fracx) * fy(i ,j+1,kk,2)
795  + fracx* fracz * fy(ii,j+1,kk,2);
796  }
797 
798  Real fzm_0 = fz(i,j,k,0);
799  Real fzm_1 = fz(i,j,k,1);
800  Real fzm_2 = fz(i,j,k,2);
801  if (apz(i,j,k) > Real(0.0) && apz(i,j,k) < Real(1.0)) {
802  int ii = i + (fcz(i,j,k,0) >= Real(0.0) ? 1 : -1);
803  int jj = j + (fcz(i,j,k,1) >= Real(0.0) ? 1 : -1);
804  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
805  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
806  fzm_0 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzm_0
807  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k,0)
808  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k,0)
809  + fracx* fracy * fz(ii,jj,k,0);
810  fzm_1 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzm_1
811  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k,1)
812  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k,1)
813  + fracx* fracy * fz(ii,jj,k,1);
814  fzm_2 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzm_2
815  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k,2)
816  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k,2)
817  + fracx* fracy * fz(ii,jj,k,2);
818  }
819 
820  Real fzp_0 = fz(i,j,k+1,0);
821  Real fzp_1 = fz(i,j,k+1,1);
822  Real fzp_2 = fz(i,j,k+1,2);
823  if (apz(i,j,k+1) > Real(0.0) && apz(i,j,k+1) < Real(1.0)) {
824  int ii = i + (fcz(i,j,k+1,0) >= Real(0.0) ? 1 : -1);
825  int jj = j + (fcz(i,j,k+1,1) >= Real(0.0) ? 1 : -1);
826  Real fracx = (ccm(ii,j,k) || ccm(ii,j,k+1)) ? std::abs(fcz(i,j,k+1,0)) : Real(0.0);
827  Real fracy = (ccm(i,jj,k) || ccm(i,jj,k+1)) ? std::abs(fcz(i,j,k+1,1)) : Real(0.0);
828  fzp_0 = (Real(1.0)-fracx)*(Real(1.0)-fracy) * fzp_0
829  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k+1,0)
830  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k+1,0)
831  + fracx* fracy * fz(ii,jj,k+1,0);
832  fzp_1 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzp_1
833  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k+1,1)
834  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k+1,1)
835  + fracx* fracy * fz(ii,jj,k+1,1);
836  fzp_2 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzp_2
837  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k+1,2)
838  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k+1,2)
839  + fracx* fracy * fz(ii,jj,k+1,2);
840  }
841 
842  Real kappa = vol(i,j,k);
843  Real feb_0 = Real(0.0);
844  Real feb_1 = Real(0.0);
845  Real feb_2 = Real(0.0);
846 
847  if (is_dirichlet)
848  {
849  Real dapx = apx(i+1,j,k)-apx(i,j,k);
850  Real dapy = apy(i,j+1,k)-apy(i,j,k);
851  Real dapz = apz(i,j,k+1)-apz(i,j,k);
852  Real anorminv = Real(1.0)/std::sqrt(dapx*dapx+dapy*dapy+dapz*dapz);
853  Real anrmx = -dapx * anorminv;
854  Real anrmy = -dapy * anorminv;
855  Real anrmz = -dapz * anorminv;
856 
857  Real velb_0 = 0, velb_1 = 0, velb_2 = 0;
858 
859  if (is_inhomog) {
860  velb_0 = velb(i,j,k,0);
861  velb_1 = velb(i,j,k,1);
862  velb_2 = velb(i,j,k,2);
863  }
864 
865  Real dx_eb = amrex::max(Real(0.3), (kappa*kappa-Real(0.25))/(Real(2.)*kappa));
866  Real dg = dx_eb / amrex::max(std::abs(anrmx),
867  std::abs(anrmy),
868  std::abs(anrmz));
869  Real dginv = Real(1.0)/dg;
870  Real gx = bc(i,j,k,0) - dg*anrmx;
871  Real gy = bc(i,j,k,1) - dg*anrmy;
872  Real gz = bc(i,j,k,2) - dg*anrmz;
873  int isx = (anrmx > Real(0.0)) ? 1 : -1;
874  int isy = (anrmy > Real(0.0)) ? 1 : -1;
875  int isz = (anrmz > Real(0.0)) ? 1 : -1;
876  int ii = i - isx;
877  int jj = j - isy;
878  int kk = k - isz;
879 
880  gx *= static_cast<Real>(isx);
881  gy *= static_cast<Real>(isy);
882  gz *= static_cast<Real>(isz);
883  Real gxy = gx*gy;
884  Real gxz = gx*gz;
885  Real gyz = gy*gz;
886  Real gxyz = gx*gy*gz;
887  Real oneggg = Real(1.0)+gx+gy+gz+gxy+gxz+gyz+gxyz;
888 
889  Real velg = oneggg * vel(i ,j ,k ,0)
890  + (-gz - gxz - gyz - gxyz) * vel(i ,j ,kk,0)
891  + (-gy - gxy - gyz - gxyz) * vel(i ,jj,k ,0)
892  + (gyz + gxyz) * vel(i ,jj,kk,0)
893  + (-gx - gxy - gxz - gxyz) * vel(ii,j ,k ,0)
894  + (gxz + gxyz) * vel(ii,j ,kk,0)
895  + (gxy + gxyz) * vel(ii,jj,k ,0)
896  + (-gxyz) * vel(ii,jj,kk,0);
897  Real dudn = (velb_0-velg) * dginv;
898 
899  velg = oneggg * vel(i ,j ,k ,1)
900  + (-gz - gxz - gyz - gxyz) * vel(i ,j ,kk,1)
901  + (-gy - gxy - gyz - gxyz) * vel(i ,jj,k ,1)
902  + (gyz + gxyz) * vel(i ,jj,kk,1)
903  + (-gx - gxy - gxz - gxyz) * vel(ii,j ,k ,1)
904  + (gxz + gxyz) * vel(ii,j ,kk,1)
905  + (gxy + gxyz) * vel(ii,jj,k ,1)
906  + (-gxyz) * vel(ii,jj,kk,1);
907  Real dvdn = (velb_1-velg) * dginv;
908 
909  velg = oneggg * vel(i ,j ,k ,2)
910  + (-gz - gxz - gyz - gxyz) * vel(i ,j ,kk,2)
911  + (-gy - gxy - gyz - gxyz) * vel(i ,jj,k ,2)
912  + (gyz + gxyz) * vel(i ,jj,kk,2)
913  + (-gx - gxy - gxz - gxyz) * vel(ii,j ,k ,2)
914  + (gxz + gxyz) * vel(ii,j ,kk,2)
915  + (gxy + gxyz) * vel(ii,jj,k ,2)
916  + (-gxyz) * vel(ii,jj,kk,2);
917  Real dwdn = (velb_2-velg) * dginv;
918 
919  // transverse derivatives are zero on the wall
920  Real dudx = dudn * anrmx;
921  Real dudy = dudn * anrmy;
922  Real dudz = dudn * anrmz;
923  Real dvdx = dvdn * anrmx;
924  Real dvdy = dvdn * anrmy;
925  Real dvdz = dvdn * anrmz;
926  Real dwdx = dwdn * anrmx;
927  Real dwdy = dwdn * anrmy;
928  Real dwdz = dwdn * anrmz;
929  Real divu = dudx + dvdy + dwdz;
930  Real xi = kapb(i,j,k);
931  Real mu = etab(i,j,k);
932  Real tautmp = (xi-Real(2./3.)*mu)*divu;
933  // Note that mu*(grad v) has been included already in MLEBABecLap
934  Real tauxx = mu*dudx + tautmp;
935  Real tauyx = mu*dvdx;
936  Real tauzx = mu*dwdx;
937  Real tauyy = mu*dvdy + tautmp;
938  Real tauxy = mu*dudy;
939  Real tauzy = mu*dwdy;
940  Real tauzz = mu*dwdz + tautmp;
941  Real tauxz = mu*dudz;
942  Real tauyz = mu*dvdz;
943  // assuming dxi == dyi == dzi
944  feb_0 = dxi*(dapx*tauxx + dapy*tauyx + dapz*tauzx);
945  feb_1 = dxi*(dapx*tauxy + dapy*tauyy + dapz*tauzy);
946  feb_2 = dxi*(dapx*tauxz + dapy*tauyz + dapz*tauzz);
947  }
948 
949  Real volinv = bscalar / kappa;
950  Ax(i,j,k,0) += volinv * (dxi*(apx(i+1,j,k)*fxp_0-apx(i,j,k)*fxm_0)
951  +dyi*(apy(i,j+1,k)*fyp_0-apy(i,j,k)*fym_0)
952  +dzi*(apz(i,j,k+1)*fzp_0-apz(i,j,k)*fzm_0)
953  +dxi* feb_0);
954  Ax(i,j,k,1) += volinv * (dxi*(apx(i+1,j,k)*fxp_1-apx(i,j,k)*fxm_1)
955  +dyi*(apy(i,j+1,k)*fyp_1-apy(i,j,k)*fym_1)
956  +dzi*(apz(i,j,k+1)*fzp_1-apz(i,j,k)*fzm_1)
957  +dxi* feb_1);
958  Ax(i,j,k,2) += volinv * (dxi*(apx(i+1,j,k)*fxp_2-apx(i,j,k)*fxm_2)
959  +dyi*(apy(i,j+1,k)*fyp_2-apy(i,j,k)*fym_2)
960  +dzi*(apz(i,j,k+1)*fzp_2-apz(i,j,k)*fzm_2)
961  +dxi* feb_2);
962  }
963  }
964  }
965  }
966 }
967 
968 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
969 void mlebtensor_flux_0 (Box const& box,
970  Array4<Real> const& Ax,
971  Array4<Real const> const& fx,
972  Array4<Real const> const& ap,
973  Real bscalar) noexcept
974  {
975  const auto lo = amrex::lbound(box);
976  const auto hi = amrex::ubound(box);
977 
978  for (int k = lo.z; k <= hi.z; ++k) {
979  for (int j = lo.y; j <= hi.y; ++j) {
980  AMREX_PRAGMA_SIMD
981  for (int i = lo.x; i <= hi.x; ++i) {
982  if (ap(i,j,k) != Real(0.0)) {
983  for (int n=0; n<AMREX_SPACEDIM; n++) {
984  Ax(i,j,k,n) += bscalar*fx(i,j,k,n);
985  }
986  }
987  //else MLEBABec::compFlux already set Ax = 0
988  }
989  }
990  }
991 }
992 
993 
994 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
995 void mlebtensor_flux_x (Box const& box, Array4<Real> const& Ax,
996  Array4<Real const> const& fx, Array4<Real const> const& apx,
997  Array4<Real const> const& fcx,
998  Real const bscalar, Array4<int const> const& ccm,
999  int face_only, Box const& xbox) noexcept
1000 {
1001  int lof = xbox.smallEnd(0);
1002  int hif = xbox.bigEnd(0);
1003  amrex::LoopConcurrent(box, [=] (int i, int j, int k) noexcept
1004  {
1005  if (!face_only || lof == i || hif == i) {
1006  if (apx(i,j,k) == Real(1.0)) {
1007  for (int n=0; n<AMREX_SPACEDIM; n++) {
1008  Ax(i,j,k,n) += bscalar*fx(i,j,k,n);
1009  }
1010  }
1011  else if (apx(i,j,k) != Real(0.0)) {
1012  Real fxm_0 = fx(i,j,k,0);
1013  Real fxm_1 = fx(i,j,k,1);
1014  Real fxm_2 = fx(i,j,k,2);
1015 
1016  int jj = j + (fcx(i,j,k,0) >= Real(0.0) ? 1 : -1);
1017  int kk = k + (fcx(i,j,k,1) >= Real(0.0) ? 1 : -1);
1018  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
1019  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
1020  fxm_0 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxm_0
1021  + fracy*(Real(1.0)-fracz) * fx(i,jj,k ,0)
1022  + fracz*(Real(1.0)-fracy) * fx(i,j ,kk,0)
1023  + fracy* fracz * fx(i,jj,kk,0);
1024  fxm_1 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxm_1
1025  + fracy*(Real(1.0)-fracz) * fx(i,jj,k ,1)
1026  + fracz*(Real(1.0)-fracy) * fx(i,j ,kk,1)
1027  + fracy* fracz * fx(i,jj,kk,1);
1028  fxm_2 = (Real(1.0)-fracy)*(Real(1.0)-fracz) *fxm_2
1029  + fracy*(Real(1.0)-fracz) * fx(i,jj,k ,2)
1030  + fracz*(Real(1.0)-fracy) * fx(i,j ,kk,2)
1031  + fracy* fracz * fx(i,jj,kk,2);
1032 
1033  Ax(i,j,k,0) += bscalar*fxm_0;
1034  Ax(i,j,k,1) += bscalar*fxm_1;
1035  Ax(i,j,k,2) += bscalar*fxm_2;
1036  }
1037  }
1038  });
1039 }
1040 
1041 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1042 void mlebtensor_flux_y (Box const& box, Array4<Real> const& Ay,
1043  Array4<Real const> const& fy, Array4<Real const> const& apy,
1044  Array4<Real const> const& fcy,
1045  Real const bscalar, Array4<int const> const& ccm,
1046  int face_only, Box const& ybox) noexcept
1047 {
1048  int lof = ybox.smallEnd(1);
1049  int hif = ybox.bigEnd(1);
1050  amrex::LoopConcurrent(box, [=] (int i, int j, int k) noexcept
1051  {
1052  if (!face_only || lof == j || hif == j) {
1053  if (apy(i,j,k) == Real(1.0)) {
1054  for (int n=0; n<AMREX_SPACEDIM; n++) {
1055  Ay(i,j,k,n) += bscalar*fy(i,j,k,n);
1056  }
1057  } else if (apy(i,j,k) != 0) {
1058  Real fym_0 = fy(i,j,k,0);
1059  Real fym_1 = fy(i,j,k,1);
1060  Real fym_2 = fy(i,j,k,2);
1061 
1062  int ii = i + (fcy(i,j,k,0) >= Real(0.0) ? 1 : -1);
1063  int kk = k + (fcy(i,j,k,1) >= Real(0.0) ? 1 : -1);
1064  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
1065  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
1066  fym_0 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fym_0
1067  + fracx*(Real(1.0)-fracz) * fy(ii,j,k ,0)
1068  + fracz*(Real(1.0)-fracx) * fy(i ,j,kk,0)
1069  + fracx* fracz * fy(ii,j,kk,0);
1070  fym_1 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fym_1
1071  + fracx*(Real(1.0)-fracz) * fy(ii,j,k ,1)
1072  + fracz*(Real(1.0)-fracx) * fy(i ,j,kk,1)
1073  + fracx* fracz * fy(ii,j,kk,1);
1074  fym_2 = (Real(1.0)-fracx)*(Real(1.0)-fracz) *fym_2
1075  + fracx*(Real(1.0)-fracz) * fy(ii,j,k ,2)
1076  + fracz*(Real(1.0)-fracx) * fy(i ,j,kk,2)
1077  + fracx* fracz * fy(ii,j,kk,2);
1078 
1079  Ay(i,j,k,0) += bscalar*fym_0;
1080  Ay(i,j,k,1) += bscalar*fym_1;
1081  Ay(i,j,k,2) += bscalar*fym_2;
1082  }
1083  //else MLEBABec::compFlux already set Ay = 0
1084  }
1085  });
1086 }
1087 
1088 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1089 void mlebtensor_flux_z (Box const& box, Array4<Real> const& Az,
1090  Array4<Real const> const& fz, Array4<Real const> const& apz,
1091  Array4<Real const> const& fcz,
1092  Real const bscalar, Array4<int const> const& ccm,
1093  int face_only, Box const& zbox) noexcept
1094 {
1095  int lof = zbox.smallEnd(2);
1096  int hif = zbox.bigEnd(2);
1097  amrex::LoopConcurrent(box, [=] (int i, int j, int k) noexcept
1098  {
1099  if (!face_only || lof == k || hif == k) {
1100  if (apz(i,j,k) == Real(1.0)) {
1101  for (int n=0; n<AMREX_SPACEDIM; n++) {
1102  Az(i,j,k,n) += bscalar*fz(i,j,k,n);
1103  }
1104  }
1105  else if (apz(i,j,k) != 0.) {
1106  Real fzm_0 = fz(i,j,k,0);
1107  Real fzm_1 = fz(i,j,k,1);
1108  Real fzm_2 = fz(i,j,k,2);
1109 
1110  int ii = i + (fcz(i,j,k,0) >= Real(0.0) ? 1 : -1);
1111  int jj = j + (fcz(i,j,k,1) >= Real(0.0) ? 1 : -1);
1112  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1113  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1114  fzm_0 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzm_0
1115  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k,0)
1116  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k,0)
1117  + fracx* fracy * fz(ii,jj,k,0);
1118  fzm_1 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzm_1
1119  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k,1)
1120  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k,1)
1121  + fracx* fracy * fz(ii,jj,k,1);
1122  fzm_2 = (Real(1.0)-fracx)*(Real(1.0)-fracy) *fzm_2
1123  + fracx*(Real(1.0)-fracy) * fz(ii,j ,k,2)
1124  + fracy*(Real(1.0)-fracx) * fz(i ,jj,k,2)
1125  + fracx* fracy * fz(ii,jj,k,2);
1126 
1127  Az(i,j,k,0) += bscalar*fzm_0;
1128  Az(i,j,k,1) += bscalar*fzm_1;
1129  Az(i,j,k,2) += bscalar*fzm_2;
1130  }
1131  //else MLEBABec::compFlux already set Az = 0
1132  }
1133  });
1134 }
1135 
1136 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1137 void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
1138  Array4<Real const> const& vel, Array4<Real const> const& apx,
1139  Array4<EBCellFlag const> const& flag,
1140  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1141 {
1142  const Real dxi = dxinv[0];
1143  const Real dyi = dxinv[1];
1144  const Real dzi = dxinv[2];
1145  const auto lo = amrex::lbound(box);
1146  const auto hi = amrex::ubound(box);
1147 
1148  for (int k = lo.z; k <= hi.z; ++k) {
1149  for (int j = lo.y; j <= hi.y; ++j) {
1150  AMREX_PRAGMA_SIMD
1151  for (int i = lo.x; i <= hi.x; ++i) {
1152  if (apx(i,j,k) == Real(0.0))
1153  {
1154  fx(i,j,k,0) = Real(0.0);
1155  fx(i,j,k,1) = Real(0.0);
1156  fx(i,j,k,2) = Real(0.0);
1157  fx(i,j,k,3) = Real(0.0);
1158  fx(i,j,k,4) = Real(0.0);
1159  fx(i,j,k,5) = Real(0.0);
1160  fx(i,j,k,6) = Real(0.0);
1161  fx(i,j,k,7) = Real(0.0);
1162  fx(i,j,k,8) = Real(0.0);
1163  }
1164  else
1165  {
1166  Real dudx = (vel(i,j,k,0) - vel(i-1,j,k,0))*dxi;
1167  Real dvdx = (vel(i,j,k,1) - vel(i-1,j,k,1))*dxi;
1168  Real dwdx = (vel(i,j,k,2) - vel(i-1,j,k,2))*dxi;
1169 
1170  int jhip = j + flag(i ,j,k).isConnected(0, 1,0);
1171  int jhim = j - flag(i ,j,k).isConnected(0,-1,0);
1172  int jlop = j + flag(i-1,j,k).isConnected(0, 1,0);
1173  int jlom = j - flag(i-1,j,k).isConnected(0,-1,0);
1174  Real whi = mlebtensor_weight(jhip-jhim);
1175  Real wlo = mlebtensor_weight(jlop-jlom);
1176  Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
1177  whi,wlo,jhip,jhim,jlop,jlom);
1178  Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
1179  whi,wlo,jhip,jhim,jlop,jlom);
1180  Real dwdy = mlebtensor_dy_on_xface(i,j,k,2,vel,dyi,
1181  whi,wlo,jhip,jhim,jlop,jlom);
1182 
1183  int khip = k + flag(i ,j,k).isConnected(0,0, 1);
1184  int khim = k - flag(i ,j,k).isConnected(0,0,-1);
1185  int klop = k + flag(i-1,j,k).isConnected(0,0, 1);
1186  int klom = k - flag(i-1,j,k).isConnected(0,0,-1);
1187  whi = mlebtensor_weight(khip-khim);
1188  wlo = mlebtensor_weight(klop-klom);
1189  Real dudz = mlebtensor_dz_on_xface(i,j,k,0,vel,dzi,
1190  whi,wlo,khip,khim,klop,klom);
1191  Real dvdz = mlebtensor_dz_on_xface(i,j,k,1,vel,dzi,
1192  whi,wlo,khip,khim,klop,klom);
1193  Real dwdz = mlebtensor_dz_on_xface(i,j,k,2,vel,dzi,
1194  whi,wlo,khip,khim,klop,klom);
1195 
1196  fx(i,j,k,0) = dudx;
1197  fx(i,j,k,1) = dvdx;
1198  fx(i,j,k,2) = dwdx;
1199  fx(i,j,k,3) = dudy;
1200  fx(i,j,k,4) = dvdy;
1201  fx(i,j,k,5) = dwdy;
1202  fx(i,j,k,6) = dudz;
1203  fx(i,j,k,7) = dvdz;
1204  fx(i,j,k,8) = dwdz;
1205  }
1206  }
1207  }
1208  }
1209 }
1210 
1211 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1212 void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
1213  Array4<Real const> const& vel, Array4<Real const> const& apy,
1214  Array4<EBCellFlag const> const& flag,
1215  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1216 {
1217  const Real dxi = dxinv[0];
1218  const Real dyi = dxinv[1];
1219  const Real dzi = dxinv[2];
1220  const auto lo = amrex::lbound(box);
1221  const auto hi = amrex::ubound(box);
1222 
1223  for (int k = lo.z; k <= hi.z; ++k) {
1224  for (int j = lo.y; j <= hi.y; ++j) {
1225  AMREX_PRAGMA_SIMD
1226  for (int i = lo.x; i <= hi.x; ++i) {
1227  if (apy(i,j,k) == Real(0.0))
1228  {
1229  fy(i,j,k,0) = Real(0.0);
1230  fy(i,j,k,1) = Real(0.0);
1231  fy(i,j,k,2) = Real(0.0);
1232  fy(i,j,k,3) = Real(0.0);
1233  fy(i,j,k,4) = Real(0.0);
1234  fy(i,j,k,5) = Real(0.0);
1235  fy(i,j,k,6) = Real(0.0);
1236  fy(i,j,k,7) = Real(0.0);
1237  fy(i,j,k,8) = Real(0.0);
1238  }
1239  else
1240  {
1241  int ihip = i + flag(i,j ,k).isConnected( 1,0,0);
1242  int ihim = i - flag(i,j ,k).isConnected(-1,0,0);
1243  int ilop = i + flag(i,j-1,k).isConnected( 1,0,0);
1244  int ilom = i - flag(i,j-1,k).isConnected(-1,0,0);
1245  Real whi = mlebtensor_weight(ihip-ihim);
1246  Real wlo = mlebtensor_weight(ilop-ilom);
1247  Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
1248  whi,wlo,ihip,ihim,ilop,ilom);
1249  Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
1250  whi,wlo,ihip,ihim,ilop,ilom);
1251  Real dwdx = mlebtensor_dx_on_yface(i,j,k,2,vel,dxi,
1252  whi,wlo,ihip,ihim,ilop,ilom);
1253 
1254  Real dudy = (vel(i,j,k,0) - vel(i,j-1,k,0))*dyi;
1255  Real dvdy = (vel(i,j,k,1) - vel(i,j-1,k,1))*dyi;
1256  Real dwdy = (vel(i,j,k,2) - vel(i,j-1,k,2))*dyi;
1257 
1258  int khip = k + flag(i,j ,k).isConnected(0,0, 1);
1259  int khim = k - flag(i,j ,k).isConnected(0,0,-1);
1260  int klop = k + flag(i,j-1,k).isConnected(0,0, 1);
1261  int klom = k - flag(i,j-1,k).isConnected(0,0,-1);
1262  whi = mlebtensor_weight(khip-khim);
1263  wlo = mlebtensor_weight(klop-klom);
1264  Real dudz = mlebtensor_dz_on_yface(i,j,k,0,vel,dzi,
1265  whi,wlo,khip,khim,klop,klom);
1266  Real dvdz = mlebtensor_dz_on_yface(i,j,k,1,vel,dzi,
1267  whi,wlo,khip,khim,klop,klom);
1268  Real dwdz = mlebtensor_dz_on_yface(i,j,k,2,vel,dzi,
1269  whi,wlo,khip,khim,klop,klom);
1270 
1271  fy(i,j,k,0) = dudx;
1272  fy(i,j,k,1) = dvdx;
1273  fy(i,j,k,2) = dwdx;
1274  fy(i,j,k,3) = dudy;
1275  fy(i,j,k,4) = dvdy;
1276  fy(i,j,k,5) = dwdy;
1277  fy(i,j,k,6) = dudz;
1278  fy(i,j,k,7) = dvdz;
1279  fy(i,j,k,8) = dwdz;
1280  }
1281  }
1282  }
1283  }
1284 }
1285 
1286 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1287 void mlebtensor_vel_grads_fz (Box const& box, Array4<Real> const& fz,
1288  Array4<Real const> const& vel, Array4<Real const> const& apz,
1289  Array4<EBCellFlag const> const& flag,
1290  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1291 {
1292  const Real dxi = dxinv[0];
1293  const Real dyi = dxinv[1];
1294  const Real dzi = dxinv[2];
1295  const auto lo = amrex::lbound(box);
1296  const auto hi = amrex::ubound(box);
1297 
1298  for (int k = lo.z; k <= hi.z; ++k) {
1299  for (int j = lo.y; j <= hi.y; ++j) {
1300  AMREX_PRAGMA_SIMD
1301  for (int i = lo.x; i <= hi.x; ++i) {
1302  if (apz(i,j,k) == Real(0.0))
1303  {
1304  fz(i,j,k,0) = Real(0.0);
1305  fz(i,j,k,1) = Real(0.0);
1306  fz(i,j,k,2) = Real(0.0);
1307  fz(i,j,k,3) = Real(0.0);
1308  fz(i,j,k,4) = Real(0.0);
1309  fz(i,j,k,5) = Real(0.0);
1310  fz(i,j,k,6) = Real(0.0);
1311  fz(i,j,k,7) = Real(0.0);
1312  fz(i,j,k,8) = Real(0.0);
1313  }
1314  else
1315  {
1316  int ihip = i + flag(i,j,k ).isConnected( 1,0,0);
1317  int ihim = i - flag(i,j,k ).isConnected(-1,0,0);
1318  int ilop = i + flag(i,j,k-1).isConnected( 1,0,0);
1319  int ilom = i - flag(i,j,k-1).isConnected(-1,0,0);
1320  Real whi = mlebtensor_weight(ihip-ihim);
1321  Real wlo = mlebtensor_weight(ilop-ilom);
1322  Real dudx = mlebtensor_dx_on_zface(i,j,k,0,vel,dxi,
1323  whi,wlo,ihip,ihim,ilop,ilom);
1324  Real dvdx = mlebtensor_dx_on_zface(i,j,k,1,vel,dxi,
1325  whi,wlo,ihip,ihim,ilop,ilom);
1326  Real dwdx = mlebtensor_dx_on_zface(i,j,k,2,vel,dxi,
1327  whi,wlo,ihip,ihim,ilop,ilom);
1328 
1329  int jhip = j + flag(i,j,k ).isConnected(0, 1,0);
1330  int jhim = j - flag(i,j,k ).isConnected(0,-1,0);
1331  int jlop = j + flag(i,j,k-1).isConnected(0, 1,0);
1332  int jlom = j - flag(i,j,k-1).isConnected(0,-1,0);
1333  whi = mlebtensor_weight(jhip-jhim);
1334  wlo = mlebtensor_weight(jlop-jlom);
1335  Real dudy = mlebtensor_dy_on_zface(i,j,k,0,vel,dyi,
1336  whi,wlo,jhip,jhim,jlop,jlom);
1337  Real dvdy = mlebtensor_dy_on_zface(i,j,k,1,vel,dyi,
1338  whi,wlo,jhip,jhim,jlop,jlom);
1339  Real dwdy = mlebtensor_dy_on_zface(i,j,k,2,vel,dyi,
1340  whi,wlo,jhip,jhim,jlop,jlom);
1341 
1342  Real dudz = (vel(i,j,k,0) - vel(i,j,k-1,0))*dzi;
1343  Real dvdz = (vel(i,j,k,1) - vel(i,j,k-1,1))*dzi;
1344  Real dwdz = (vel(i,j,k,2) - vel(i,j,k-1,2))*dzi;
1345 
1346  fz(i,j,k,0) = dudx;
1347  fz(i,j,k,1) = dvdx;
1348  fz(i,j,k,2) = dwdx;
1349  fz(i,j,k,3) = dudy;
1350  fz(i,j,k,4) = dvdy;
1351  fz(i,j,k,5) = dwdy;
1352  fz(i,j,k,6) = dudz;
1353  fz(i,j,k,7) = dvdz;
1354  fz(i,j,k,8) = dwdz;
1355  }
1356  }
1357  }
1358  }
1359 }
1360 
1361 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1362 void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
1363  Array4<Real const> const& vel,
1364  Array4<Real const> const& apx,
1365  Array4<EBCellFlag const> const& flag,
1366  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1367  Array4<Real const> const& bvxlo,
1368  Array4<Real const> const& bvxhi,
1369  Array2D<BoundCond,
1370  0,2*AMREX_SPACEDIM,
1371  0,AMREX_SPACEDIM> const& bct,
1372  Dim3 const& dlo, Dim3 const& dhi) noexcept
1373 {
1374  const Real dxi = dxinv[0];
1375  const Real dyi = dxinv[1];
1376  const Real dzi = dxinv[2];
1377  const auto lo = amrex::lbound(box);
1378  const auto hi = amrex::ubound(box);
1379 
1380  for (int k = lo.z; k <= hi.z; ++k) {
1381  for (int j = lo.y; j <= hi.y; ++j) {
1382  for (int i = lo.x; i <= hi.x; ++i) {
1383  if (apx(i,j,k) == Real(0.0))
1384  {
1385  fx(i,j,k,0) = Real(0.0);
1386  fx(i,j,k,1) = Real(0.0);
1387  fx(i,j,k,2) = Real(0.0);
1388  fx(i,j,k,3) = Real(0.0);
1389  fx(i,j,k,4) = Real(0.0);
1390  fx(i,j,k,5) = Real(0.0);
1391  fx(i,j,k,6) = Real(0.0);
1392  fx(i,j,k,7) = Real(0.0);
1393  fx(i,j,k,8) = Real(0.0);
1394  }
1395  else
1396  {
1397  Real dudx = (vel(i,j,k,0) - vel(i-1,j,k,0))*dxi;
1398  Real dvdx = (vel(i,j,k,1) - vel(i-1,j,k,1))*dxi;
1399  Real dwdx = (vel(i,j,k,2) - vel(i-1,j,k,2))*dxi;
1400 
1401  int jhip = j + flag(i ,j,k).isConnected(0, 1,0);
1402  int jhim = j - flag(i ,j,k).isConnected(0,-1,0);
1403  int jlop = j + flag(i-1,j,k).isConnected(0, 1,0);
1404  int jlom = j - flag(i-1,j,k).isConnected(0,-1,0);
1405  Real whi = mlebtensor_weight(jhip-jhim);
1406  Real wlo = mlebtensor_weight(jlop-jlom);
1407  Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
1408  bvxlo,bvxhi,bct,dlo,dhi,
1409  whi,wlo,jhip,jhim,jlop,jlom);
1410  Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
1411  bvxlo,bvxhi,bct,dlo,dhi,
1412  whi,wlo,jhip,jhim,jlop,jlom);
1413  Real dwdy = mlebtensor_dy_on_xface(i,j,k,2,vel,dyi,
1414  bvxlo,bvxhi,bct,dlo,dhi,
1415  whi,wlo,jhip,jhim,jlop,jlom);
1416 
1417  int khip = k + flag(i ,j,k).isConnected(0,0, 1);
1418  int khim = k - flag(i ,j,k).isConnected(0,0,-1);
1419  int klop = k + flag(i-1,j,k).isConnected(0,0, 1);
1420  int klom = k - flag(i-1,j,k).isConnected(0,0,-1);
1421  whi = mlebtensor_weight(khip-khim);
1422  wlo = mlebtensor_weight(klop-klom);
1423  Real dudz = mlebtensor_dz_on_xface(i,j,k,0,vel,dzi,
1424  bvxlo,bvxhi,bct,dlo,dhi,
1425  whi,wlo,khip,khim,klop,klom);
1426  Real dvdz = mlebtensor_dz_on_xface(i,j,k,1,vel,dzi,
1427  bvxlo,bvxhi,bct,dlo,dhi,
1428  whi,wlo,khip,khim,klop,klom);
1429  Real dwdz = mlebtensor_dz_on_xface(i,j,k,2,vel,dzi,
1430  bvxlo,bvxhi,bct,dlo,dhi,
1431  whi,wlo,khip,khim,klop,klom);
1432 
1433  fx(i,j,k,0) = dudx;
1434  fx(i,j,k,1) = dvdx;
1435  fx(i,j,k,2) = dwdx;
1436  fx(i,j,k,3) = dudy;
1437  fx(i,j,k,4) = dvdy;
1438  fx(i,j,k,5) = dwdy;
1439  fx(i,j,k,6) = dudz;
1440  fx(i,j,k,7) = dvdz;
1441  fx(i,j,k,8) = dwdz;
1442  }
1443  }
1444  }
1445  }
1446 }
1447 
1448 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1449 void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
1450  Array4<Real const> const& vel,
1451  Array4<Real const> const& apy,
1452  Array4<EBCellFlag const> const& flag,
1453  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1454  Array4<Real const> const& bvylo,
1455  Array4<Real const> const& bvyhi,
1456  Array2D<BoundCond,
1457  0,2*AMREX_SPACEDIM,
1458  0,AMREX_SPACEDIM> const& bct,
1459  Dim3 const& dlo, Dim3 const& dhi) noexcept
1460 {
1461  const Real dxi = dxinv[0];
1462  const Real dyi = dxinv[1];
1463  const Real dzi = dxinv[2];
1464  const auto lo = amrex::lbound(box);
1465  const auto hi = amrex::ubound(box);
1466 
1467  for (int k = lo.z; k <= hi.z; ++k) {
1468  for (int j = lo.y; j <= hi.y; ++j) {
1469  for (int i = lo.x; i <= hi.x; ++i) {
1470  if (apy(i,j,k) == Real(0.0))
1471  {
1472  fy(i,j,k,0) = Real(0.0);
1473  fy(i,j,k,1) = Real(0.0);
1474  fy(i,j,k,2) = Real(0.0);
1475  fy(i,j,k,3) = Real(0.0);
1476  fy(i,j,k,4) = Real(0.0);
1477  fy(i,j,k,5) = Real(0.0);
1478  fy(i,j,k,6) = Real(0.0);
1479  fy(i,j,k,7) = Real(0.0);
1480  fy(i,j,k,8) = Real(0.0);
1481  }
1482  else
1483  {
1484  int ihip = i + flag(i,j ,k).isConnected( 1,0,0);
1485  int ihim = i - flag(i,j ,k).isConnected(-1,0,0);
1486  int ilop = i + flag(i,j-1,k).isConnected( 1,0,0);
1487  int ilom = i - flag(i,j-1,k).isConnected(-1,0,0);
1488  Real whi = mlebtensor_weight(ihip-ihim);
1489  Real wlo = mlebtensor_weight(ilop-ilom);
1490  Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
1491  bvylo,bvyhi,bct,dlo,dhi,
1492  whi,wlo,ihip,ihim,ilop,ilom);
1493  Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
1494  bvylo,bvyhi,bct,dlo,dhi,
1495  whi,wlo,ihip,ihim,ilop,ilom);
1496  Real dwdx = mlebtensor_dx_on_yface(i,j,k,2,vel,dxi,
1497  bvylo,bvyhi,bct,dlo,dhi,
1498  whi,wlo,ihip,ihim,ilop,ilom);
1499 
1500  Real dudy = (vel(i,j,k,0) - vel(i,j-1,k,0))*dyi;
1501  Real dvdy = (vel(i,j,k,1) - vel(i,j-1,k,1))*dyi;
1502  Real dwdy = (vel(i,j,k,2) - vel(i,j-1,k,2))*dyi;
1503 
1504  int khip = k + flag(i,j ,k).isConnected(0,0, 1);
1505  int khim = k - flag(i,j ,k).isConnected(0,0,-1);
1506  int klop = k + flag(i,j-1,k).isConnected(0,0, 1);
1507  int klom = k - flag(i,j-1,k).isConnected(0,0,-1);
1508  whi = mlebtensor_weight(khip-khim);
1509  wlo = mlebtensor_weight(klop-klom);
1510  Real dudz = mlebtensor_dz_on_yface(i,j,k,0,vel,dzi,
1511  bvylo,bvyhi,bct,dlo,dhi,
1512  whi,wlo,khip,khim,klop,klom);
1513  Real dvdz = mlebtensor_dz_on_yface(i,j,k,1,vel,dzi,
1514  bvylo,bvyhi,bct,dlo,dhi,
1515  whi,wlo,khip,khim,klop,klom);
1516  Real dwdz = mlebtensor_dz_on_yface(i,j,k,2,vel,dzi,
1517  bvylo,bvyhi,bct,dlo,dhi,
1518  whi,wlo,khip,khim,klop,klom);
1519 
1520  fy(i,j,k,0) = dudx;
1521  fy(i,j,k,1) = dvdx;
1522  fy(i,j,k,2) = dwdx;
1523  fy(i,j,k,3) = dudy;
1524  fy(i,j,k,4) = dvdy;
1525  fy(i,j,k,5) = dwdy;
1526  fy(i,j,k,6) = dudz;
1527  fy(i,j,k,7) = dvdz;
1528  fy(i,j,k,8) = dwdz;
1529  }
1530  }
1531  }
1532  }
1533 }
1534 
1535 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1536 void mlebtensor_vel_grads_fz (Box const& box, Array4<Real> const& fz,
1537  Array4<Real const> const& vel,
1538  Array4<Real const> const& apz,
1539  Array4<EBCellFlag const> const& flag,
1540  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1541  Array4<Real const> const& bvzlo,
1542  Array4<Real const> const& bvzhi,
1543  Array2D<BoundCond,
1544  0,2*AMREX_SPACEDIM,
1545  0,AMREX_SPACEDIM> const& bct,
1546  Dim3 const& dlo, Dim3 const& dhi) noexcept
1547 {
1548  const Real dxi = dxinv[0];
1549  const Real dyi = dxinv[1];
1550  const Real dzi = dxinv[2];
1551  const auto lo = amrex::lbound(box);
1552  const auto hi = amrex::ubound(box);
1553 
1554  for (int k = lo.z; k <= hi.z; ++k) {
1555  for (int j = lo.y; j <= hi.y; ++j) {
1556  for (int i = lo.x; i <= hi.x; ++i) {
1557  if (apz(i,j,k) == Real(0.0))
1558  {
1559  fz(i,j,k,0) = Real(0.0);
1560  fz(i,j,k,1) = Real(0.0);
1561  fz(i,j,k,2) = Real(0.0);
1562  fz(i,j,k,3) = Real(0.0);
1563  fz(i,j,k,4) = Real(0.0);
1564  fz(i,j,k,5) = Real(0.0);
1565  fz(i,j,k,6) = Real(0.0);
1566  fz(i,j,k,7) = Real(0.0);
1567  fz(i,j,k,8) = Real(0.0);
1568  }
1569  else
1570  {
1571  int ihip = i + flag(i,j,k ).isConnected( 1,0,0);
1572  int ihim = i - flag(i,j,k ).isConnected(-1,0,0);
1573  int ilop = i + flag(i,j,k-1).isConnected( 1,0,0);
1574  int ilom = i - flag(i,j,k-1).isConnected(-1,0,0);
1575  Real whi = mlebtensor_weight(ihip-ihim);
1576  Real wlo = mlebtensor_weight(ilop-ilom);
1577  Real dudx = mlebtensor_dx_on_zface(i,j,k,0,vel,dxi,
1578  bvzlo,bvzhi,bct,dlo,dhi,
1579  whi,wlo,ihip,ihim,ilop,ilom);
1580  Real dvdx = mlebtensor_dx_on_zface(i,j,k,1,vel,dxi,
1581  bvzlo,bvzhi,bct,dlo,dhi,
1582  whi,wlo,ihip,ihim,ilop,ilom);
1583  Real dwdx = mlebtensor_dx_on_zface(i,j,k,2,vel,dxi,
1584  bvzlo,bvzhi,bct,dlo,dhi,
1585  whi,wlo,ihip,ihim,ilop,ilom);
1586 
1587  int jhip = j + flag(i,j,k ).isConnected(0, 1,0);
1588  int jhim = j - flag(i,j,k ).isConnected(0,-1,0);
1589  int jlop = j + flag(i,j,k-1).isConnected(0, 1,0);
1590  int jlom = j - flag(i,j,k-1).isConnected(0,-1,0);
1591  whi = mlebtensor_weight(jhip-jhim);
1592  wlo = mlebtensor_weight(jlop-jlom);
1593  Real dudy = mlebtensor_dy_on_zface(i,j,k,0,vel,dyi,
1594  bvzlo,bvzhi,bct,dlo,dhi,
1595  whi,wlo,jhip,jhim,jlop,jlom);
1596  Real dvdy = mlebtensor_dy_on_zface(i,j,k,1,vel,dyi,
1597  bvzlo,bvzhi,bct,dlo,dhi,
1598  whi,wlo,jhip,jhim,jlop,jlom);
1599  Real dwdy = mlebtensor_dy_on_zface(i,j,k,2,vel,dyi,
1600  bvzlo,bvzhi,bct,dlo,dhi,
1601  whi,wlo,jhip,jhim,jlop,jlom);
1602 
1603  Real dudz = (vel(i,j,k,0) - vel(i,j,k-1,0))*dzi;
1604  Real dvdz = (vel(i,j,k,1) - vel(i,j,k-1,1))*dzi;
1605  Real dwdz = (vel(i,j,k,2) - vel(i,j,k-1,2))*dzi;
1606 
1607  fz(i,j,k,0) = dudx;
1608  fz(i,j,k,1) = dvdx;
1609  fz(i,j,k,2) = dwdx;
1610  fz(i,j,k,3) = dudy;
1611  fz(i,j,k,4) = dvdy;
1612  fz(i,j,k,5) = dwdy;
1613  fz(i,j,k,6) = dudz;
1614  fz(i,j,k,7) = dvdz;
1615  fz(i,j,k,8) = dwdz;
1616  }
1617  }
1618  }
1619  }
1620 }
1621 
1622 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1623 void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
1624  Array4<Real const> const& vel, Array4<Real const> const& apx,
1625  Array4<EBCellFlag const> const& flag,
1626  Array4<int const> const& ccm,
1627  Array4<Real const> const& fcx,
1628  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1629 {
1630  const Real dxi = dxinv[0];
1631  const Real dyi = dxinv[1];
1632  const Real dzi = dxinv[2];
1633  const auto lo = amrex::lbound(box);
1634  const auto hi = amrex::ubound(box);
1635 
1636  for (int k = lo.z; k <= hi.z; ++k) {
1637  for (int j = lo.y; j <= hi.y; ++j) {
1638  AMREX_PRAGMA_SIMD
1639  for (int i = lo.x; i <= hi.x; ++i) {
1640  if (apx(i,j,k) == Real(0.0))
1641  {
1642  fx(i,j,k,0) = Real(0.0);
1643  fx(i,j,k,1) = Real(0.0);
1644  fx(i,j,k,2) = Real(0.0);
1645  fx(i,j,k,3) = Real(0.0);
1646  fx(i,j,k,4) = Real(0.0);
1647  fx(i,j,k,5) = Real(0.0);
1648  fx(i,j,k,6) = Real(0.0);
1649  fx(i,j,k,7) = Real(0.0);
1650  fx(i,j,k,8) = Real(0.0);
1651  }
1652  else
1653  {
1654  Real dudx = (vel(i,j,k,0) - vel(i-1,j,k,0))*dxi;
1655  Real dvdx = (vel(i,j,k,1) - vel(i-1,j,k,1))*dxi;
1656  Real dwdx = (vel(i,j,k,2) - vel(i-1,j,k,2))*dxi;
1657  if (apx(i,j,k) < Real(1.0)) {
1658  int jj = j + (fcx(i,j,k,0) >= Real(0.0) ? 1 : -1);
1659  int kk = k + (fcx(i,j,k,1) >= Real(0.0) ? 1 : -1);
1660  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
1661  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
1662  dudx = (Real(1.0)-fracy)*(Real(1.0)-fracz) * dudx
1663  + fracy*(Real(1.0)-fracz) * (vel(i,jj,k ,0) - vel(i-1,jj,k ,0))*dxi
1664  + fracz*(Real(1.0)-fracy) * (vel(i,j ,kk,0) - vel(i-1,j ,kk,0))*dxi
1665  + fracy* fracz * (vel(i,jj,kk,0) - vel(i-1,jj,kk,0))*dxi;
1666  dvdx = (Real(1.0)-fracy)*(Real(1.0)-fracz) * dvdx
1667  + fracy*(Real(1.0)-fracz) * (vel(i,jj,k ,1) - vel(i-1,jj,k ,1))*dxi
1668  + fracz*(Real(1.0)-fracy) * (vel(i,j ,kk,1) - vel(i-1,j ,kk,1))*dxi
1669  + fracy* fracz * (vel(i,jj,kk,1) - vel(i-1,jj,kk,1))*dxi;
1670  dwdx = (Real(1.0)-fracy)*(Real(1.0)-fracz) * dwdx
1671  + fracy*(Real(1.0)-fracz) * (vel(i,jj,k ,2) - vel(i-1,jj,k ,2))*dxi
1672  + fracz*(Real(1.0)-fracy) * (vel(i,j ,kk,2) - vel(i-1,j ,kk,2))*dxi
1673  + fracy* fracz * (vel(i,jj,kk,2) - vel(i-1,jj,kk,2))*dxi;
1674  }
1675 
1676  int jhip = j + flag(i ,j,k).isConnected(0, 1,0);
1677  int jhim = j - flag(i ,j,k).isConnected(0,-1,0);
1678  int jlop = j + flag(i-1,j,k).isConnected(0, 1,0);
1679  int jlom = j - flag(i-1,j,k).isConnected(0,-1,0);
1680  Real whi = mlebtensor_weight(jhip-jhim);
1681  Real wlo = mlebtensor_weight(jlop-jlom);
1682  Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
1683  whi,wlo,jhip,jhim,jlop,jlom);
1684  Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
1685  whi,wlo,jhip,jhim,jlop,jlom);
1686  Real dwdy = mlebtensor_dy_on_xface(i,j,k,2,vel,dyi,
1687  whi,wlo,jhip,jhim,jlop,jlom);
1688 
1689  int khip = k + flag(i ,j,k).isConnected(0,0, 1);
1690  int khim = k - flag(i ,j,k).isConnected(0,0,-1);
1691  int klop = k + flag(i-1,j,k).isConnected(0,0, 1);
1692  int klom = k - flag(i-1,j,k).isConnected(0,0,-1);
1693  whi = mlebtensor_weight(khip-khim);
1694  wlo = mlebtensor_weight(klop-klom);
1695  Real dudz = mlebtensor_dz_on_xface(i,j,k,0,vel,dzi,
1696  whi,wlo,khip,khim,klop,klom);
1697  Real dvdz = mlebtensor_dz_on_xface(i,j,k,1,vel,dzi,
1698  whi,wlo,khip,khim,klop,klom);
1699  Real dwdz = mlebtensor_dz_on_xface(i,j,k,2,vel,dzi,
1700  whi,wlo,khip,khim,klop,klom);
1701 
1702  fx(i,j,k,0) = dudx;
1703  fx(i,j,k,1) = dvdx;
1704  fx(i,j,k,2) = dwdx;
1705  fx(i,j,k,3) = dudy;
1706  fx(i,j,k,4) = dvdy;
1707  fx(i,j,k,5) = dwdy;
1708  fx(i,j,k,6) = dudz;
1709  fx(i,j,k,7) = dvdz;
1710  fx(i,j,k,8) = dwdz;
1711  }
1712  }
1713  }
1714  }
1715 }
1716 
1717 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1718 void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
1719  Array4<Real const> const& vel, Array4<Real const> const& apy,
1720  Array4<EBCellFlag const> const& flag,
1721  Array4<int const> const& ccm,
1722  Array4<Real const> const& fcy,
1723  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1724 {
1725  const Real dxi = dxinv[0];
1726  const Real dyi = dxinv[1];
1727  const Real dzi = dxinv[2];
1728  const auto lo = amrex::lbound(box);
1729  const auto hi = amrex::ubound(box);
1730 
1731  for (int k = lo.z; k <= hi.z; ++k) {
1732  for (int j = lo.y; j <= hi.y; ++j) {
1733  AMREX_PRAGMA_SIMD
1734  for (int i = lo.x; i <= hi.x; ++i) {
1735  if (apy(i,j,k) == Real(0.0))
1736  {
1737  fy(i,j,k,0) = Real(0.0);
1738  fy(i,j,k,1) = Real(0.0);
1739  fy(i,j,k,2) = Real(0.0);
1740  fy(i,j,k,3) = Real(0.0);
1741  fy(i,j,k,4) = Real(0.0);
1742  fy(i,j,k,5) = Real(0.0);
1743  fy(i,j,k,6) = Real(0.0);
1744  fy(i,j,k,7) = Real(0.0);
1745  fy(i,j,k,8) = Real(0.0);
1746  }
1747  else
1748  {
1749  int ihip = i + flag(i,j ,k).isConnected( 1,0,0);
1750  int ihim = i - flag(i,j ,k).isConnected(-1,0,0);
1751  int ilop = i + flag(i,j-1,k).isConnected( 1,0,0);
1752  int ilom = i - flag(i,j-1,k).isConnected(-1,0,0);
1753  Real whi = mlebtensor_weight(ihip-ihim);
1754  Real wlo = mlebtensor_weight(ilop-ilom);
1755  Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
1756  whi,wlo,ihip,ihim,ilop,ilom);
1757  Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
1758  whi,wlo,ihip,ihim,ilop,ilom);
1759  Real dwdx = mlebtensor_dx_on_yface(i,j,k,2,vel,dxi,
1760  whi,wlo,ihip,ihim,ilop,ilom);
1761 
1762  Real dudy = (vel(i,j,k,0) - vel(i,j-1,k,0))*dyi;
1763  Real dvdy = (vel(i,j,k,1) - vel(i,j-1,k,1))*dyi;
1764  Real dwdy = (vel(i,j,k,2) - vel(i,j-1,k,2))*dyi;
1765  if (apy(i,j,k) < Real(1.0)) {
1766  int ii = i + (fcy(i,j,k,0) >= Real(0.0) ? 1 : -1);
1767  int kk = k + (fcy(i,j,k,1) >= Real(0.0) ? 1 : -1);
1768  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
1769  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
1770  dudy = (Real(1.0)-fracx)*(Real(1.0)-fracz) * dudy
1771  + fracx*(Real(1.0)-fracz) * (vel(ii,j,k ,0) - vel(ii,j-1,k ,0))*dyi
1772  + fracz*(Real(1.0)-fracx) * (vel(i ,j,kk,0) - vel(i ,j-1,kk,0))*dyi
1773  + fracx* fracz * (vel(ii,j,kk,0) - vel(ii,j-1,kk,0))*dyi;
1774  dvdy = (Real(1.0)-fracx)*(Real(1.0)-fracz) * dvdy
1775  + fracx*(Real(1.0)-fracz) * (vel(ii,j,k ,1) - vel(ii,j-1,k ,1))*dyi
1776  + fracz*(Real(1.0)-fracx) * (vel(i ,j,kk,1) - vel(i ,j-1,kk,1))*dyi
1777  + fracx* fracz * (vel(ii,j,kk,1) - vel(ii,j-1,kk,1))*dyi;
1778  dwdy = (Real(1.0)-fracx)*(Real(1.0)-fracz) * dwdy
1779  + fracx*(Real(1.0)-fracz) * (vel(ii,j,k ,2) - vel(ii,j-1,k ,2))*dyi
1780  + fracz*(Real(1.0)-fracx) * (vel(i ,j,kk,2) - vel(i ,j-1,kk,2))*dyi
1781  + fracx* fracz * (vel(ii,j,kk,2) - vel(ii,j-1,kk,2))*dyi;
1782  }
1783 
1784  int khip = k + flag(i,j ,k).isConnected(0,0, 1);
1785  int khim = k - flag(i,j ,k).isConnected(0,0,-1);
1786  int klop = k + flag(i,j-1,k).isConnected(0,0, 1);
1787  int klom = k - flag(i,j-1,k).isConnected(0,0,-1);
1788  whi = mlebtensor_weight(khip-khim);
1789  wlo = mlebtensor_weight(klop-klom);
1790  Real dudz = mlebtensor_dz_on_yface(i,j,k,0,vel,dzi,
1791  whi,wlo,khip,khim,klop,klom);
1792  Real dvdz = mlebtensor_dz_on_yface(i,j,k,1,vel,dzi,
1793  whi,wlo,khip,khim,klop,klom);
1794  Real dwdz = mlebtensor_dz_on_yface(i,j,k,2,vel,dzi,
1795  whi,wlo,khip,khim,klop,klom);
1796 
1797  fy(i,j,k,0) = dudx;
1798  fy(i,j,k,1) = dvdx;
1799  fy(i,j,k,2) = dwdx;
1800  fy(i,j,k,3) = dudy;
1801  fy(i,j,k,4) = dvdy;
1802  fy(i,j,k,5) = dwdy;
1803  fy(i,j,k,6) = dudz;
1804  fy(i,j,k,7) = dvdz;
1805  fy(i,j,k,8) = dwdz;
1806  }
1807  }
1808  }
1809  }
1810 }
1811 
1812 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1813 void mlebtensor_vel_grads_fz (Box const& box, Array4<Real> const& fz,
1814  Array4<Real const> const& vel, Array4<Real const> const& apz,
1815  Array4<EBCellFlag const> const& flag,
1816  Array4<int const> const& ccm,
1817  Array4<Real const> const& fcz,
1818  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1819 {
1820  const Real dxi = dxinv[0];
1821  const Real dyi = dxinv[1];
1822  const Real dzi = dxinv[2];
1823  const auto lo = amrex::lbound(box);
1824  const auto hi = amrex::ubound(box);
1825 
1826  for (int k = lo.z; k <= hi.z; ++k) {
1827  for (int j = lo.y; j <= hi.y; ++j) {
1828  AMREX_PRAGMA_SIMD
1829  for (int i = lo.x; i <= hi.x; ++i) {
1830  if (apz(i,j,k) == Real(0.0))
1831  {
1832  fz(i,j,k,0) = Real(0.0);
1833  fz(i,j,k,1) = Real(0.0);
1834  fz(i,j,k,2) = Real(0.0);
1835  fz(i,j,k,3) = Real(0.0);
1836  fz(i,j,k,4) = Real(0.0);
1837  fz(i,j,k,5) = Real(0.0);
1838  fz(i,j,k,6) = Real(0.0);
1839  fz(i,j,k,7) = Real(0.0);
1840  fz(i,j,k,8) = Real(0.0);
1841  }
1842  else
1843  {
1844  int ihip = i + flag(i,j,k ).isConnected( 1,0,0);
1845  int ihim = i - flag(i,j,k ).isConnected(-1,0,0);
1846  int ilop = i + flag(i,j,k-1).isConnected( 1,0,0);
1847  int ilom = i - flag(i,j,k-1).isConnected(-1,0,0);
1848  Real whi = mlebtensor_weight(ihip-ihim);
1849  Real wlo = mlebtensor_weight(ilop-ilom);
1850  Real dudx = mlebtensor_dx_on_zface(i,j,k,0,vel,dxi,
1851  whi,wlo,ihip,ihim,ilop,ilom);
1852  Real dvdx = mlebtensor_dx_on_zface(i,j,k,1,vel,dxi,
1853  whi,wlo,ihip,ihim,ilop,ilom);
1854  Real dwdx = mlebtensor_dx_on_zface(i,j,k,2,vel,dxi,
1855  whi,wlo,ihip,ihim,ilop,ilom);
1856 
1857  int jhip = j + flag(i,j,k ).isConnected(0, 1,0);
1858  int jhim = j - flag(i,j,k ).isConnected(0,-1,0);
1859  int jlop = j + flag(i,j,k-1).isConnected(0, 1,0);
1860  int jlom = j - flag(i,j,k-1).isConnected(0,-1,0);
1861  whi = mlebtensor_weight(jhip-jhim);
1862  wlo = mlebtensor_weight(jlop-jlom);
1863  Real dudy = mlebtensor_dy_on_zface(i,j,k,0,vel,dyi,
1864  whi,wlo,jhip,jhim,jlop,jlom);
1865  Real dvdy = mlebtensor_dy_on_zface(i,j,k,1,vel,dyi,
1866  whi,wlo,jhip,jhim,jlop,jlom);
1867  Real dwdy = mlebtensor_dy_on_zface(i,j,k,2,vel,dyi,
1868  whi,wlo,jhip,jhim,jlop,jlom);
1869 
1870  Real dudz = (vel(i,j,k,0) - vel(i,j,k-1,0))*dzi;
1871  Real dvdz = (vel(i,j,k,1) - vel(i,j,k-1,1))*dzi;
1872  Real dwdz = (vel(i,j,k,2) - vel(i,j,k-1,2))*dzi;
1873  if (apz(i,j,k) < Real(1.0)) {
1874  int ii = i + (fcz(i,j,k,0) >= Real(0.0) ? 1 : -1);
1875  int jj = j + (fcz(i,j,k,1) >= Real(0.0) ? 1 : -1);
1876  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
1877  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
1878  dudz = (Real(1.0)-fracx)*(Real(1.0)-fracy) * dudz
1879  + fracx*(Real(1.0)-fracy) * (vel(ii,j ,k,0) - vel(ii,j ,k-1,0))*dzi
1880  + fracy*(Real(1.0)-fracx) * (vel(i ,jj,k,0) - vel(i ,jj,k-1,0))*dzi
1881  + fracx* fracy * (vel(ii,jj,k,0) - vel(ii,jj,k-1,0))*dzi;
1882  dvdz = (Real(1.0)-fracx)*(Real(1.0)-fracy) * dvdz
1883  + fracx*(Real(1.0)-fracy) * (vel(ii,j ,k,1) - vel(ii,j ,k-1,1))*dzi
1884  + fracy*(Real(1.0)-fracx) * (vel(i ,jj,k,1) - vel(i ,jj,k-1,1))*dzi
1885  + fracx* fracy * (vel(ii,jj,k,1) - vel(ii,jj,k-1,1))*dzi;
1886  dwdz = (Real(1.0)-fracx)*(Real(1.0)-fracy) * dwdz
1887  + fracx*(Real(1.0)-fracy) * (vel(ii,j ,k,2) - vel(ii,j ,k-1,2))*dzi
1888  + fracy*(Real(1.0)-fracx) * (vel(i ,jj,k,2) - vel(i ,jj,k-1,2))*dzi
1889  + fracx* fracy * (vel(ii,jj,k,2) - vel(ii,jj,k-1,2))*dzi;
1890  }
1891 
1892  fz(i,j,k,0) = dudx;
1893  fz(i,j,k,1) = dvdx;
1894  fz(i,j,k,2) = dwdx;
1895  fz(i,j,k,3) = dudy;
1896  fz(i,j,k,4) = dvdy;
1897  fz(i,j,k,5) = dwdy;
1898  fz(i,j,k,6) = dudz;
1899  fz(i,j,k,7) = dvdz;
1900  fz(i,j,k,8) = dwdz;
1901  }
1902  }
1903  }
1904  }
1905 }
1906 
1907 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1908 void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
1909  Array4<Real const> const& vel,
1910  Array4<Real const> const& apx,
1911  Array4<EBCellFlag const> const& flag,
1912  Array4<int const> const& ccm,
1913  Array4<Real const> const& fcx,
1914  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1915  Array4<Real const> const& bvxlo,
1916  Array4<Real const> const& bvxhi,
1917  Array2D<BoundCond,
1918  0,2*AMREX_SPACEDIM,
1919  0,AMREX_SPACEDIM> const& bct,
1920  Dim3 const& dlo, Dim3 const& dhi) noexcept
1921 {
1922  const Real dxi = dxinv[0];
1923  const Real dyi = dxinv[1];
1924  const Real dzi = dxinv[2];
1925  const auto lo = amrex::lbound(box);
1926  const auto hi = amrex::ubound(box);
1927 
1928  for (int k = lo.z; k <= hi.z; ++k) {
1929  for (int j = lo.y; j <= hi.y; ++j) {
1930  for (int i = lo.x; i <= hi.x; ++i) {
1931  if (apx(i,j,k) == Real(0.0))
1932  {
1933  fx(i,j,k,0) = Real(0.0);
1934  fx(i,j,k,1) = Real(0.0);
1935  fx(i,j,k,2) = Real(0.0);
1936  fx(i,j,k,3) = Real(0.0);
1937  fx(i,j,k,4) = Real(0.0);
1938  fx(i,j,k,5) = Real(0.0);
1939  fx(i,j,k,6) = Real(0.0);
1940  fx(i,j,k,7) = Real(0.0);
1941  fx(i,j,k,8) = Real(0.0);
1942  }
1943  else
1944  {
1945  Real dudx = (vel(i,j,k,0) - vel(i-1,j,k,0))*dxi;
1946  Real dvdx = (vel(i,j,k,1) - vel(i-1,j,k,1))*dxi;
1947  Real dwdx = (vel(i,j,k,2) - vel(i-1,j,k,2))*dxi;
1948  if (apx(i,j,k) < Real(1.0)) {
1949  int jj = j + (fcx(i,j,k,0) >= Real(0.0) ? 1 : -1);
1950  int kk = k + (fcx(i,j,k,1) >= Real(0.0) ? 1 : -1);
1951  Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx(i,j,k,0)) : Real(0.0);
1952  Real fracz = (ccm(i-1,j,kk) || ccm(i,j,kk)) ? std::abs(fcx(i,j,k,1)) : Real(0.0);
1953  dudx = (Real(1.0)-fracy)*(Real(1.0)-fracz) * dudx
1954  + fracy*(Real(1.0)-fracz) * (vel(i,jj,k ,0) - vel(i-1,jj,k ,0))*dxi
1955  + fracz*(Real(1.0)-fracy) * (vel(i,j ,kk,0) - vel(i-1,j ,kk,0))*dxi
1956  + fracy* fracz * (vel(i,jj,kk,0) - vel(i-1,jj,kk,0))*dxi;
1957  dvdx = (Real(1.0)-fracy)*(Real(1.0)-fracz) * dvdx
1958  + fracy*(Real(1.0)-fracz) * (vel(i,jj,k ,1) - vel(i-1,jj,k ,1))*dxi
1959  + fracz*(Real(1.0)-fracy) * (vel(i,j ,kk,1) - vel(i-1,j ,kk,1))*dxi
1960  + fracy* fracz * (vel(i,jj,kk,1) - vel(i-1,jj,kk,1))*dxi;
1961  dwdx = (Real(1.0)-fracy)*(Real(1.0)-fracz) * dwdx
1962  + fracy*(Real(1.0)-fracz) * (vel(i,jj,k ,2) - vel(i-1,jj,k ,2))*dxi
1963  + fracz*(Real(1.0)-fracy) * (vel(i,j ,kk,2) - vel(i-1,j ,kk,2))*dxi
1964  + fracy* fracz * (vel(i,jj,kk,2) - vel(i-1,jj,kk,2))*dxi;
1965  }
1966 
1967  int jhip = j + flag(i ,j,k).isConnected(0, 1,0);
1968  int jhim = j - flag(i ,j,k).isConnected(0,-1,0);
1969  int jlop = j + flag(i-1,j,k).isConnected(0, 1,0);
1970  int jlom = j - flag(i-1,j,k).isConnected(0,-1,0);
1971  Real whi = mlebtensor_weight(jhip-jhim);
1972  Real wlo = mlebtensor_weight(jlop-jlom);
1973  Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
1974  bvxlo,bvxhi,bct,dlo,dhi,
1975  whi,wlo,jhip,jhim,jlop,jlom);
1976  Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
1977  bvxlo,bvxhi,bct,dlo,dhi,
1978  whi,wlo,jhip,jhim,jlop,jlom);
1979  Real dwdy = mlebtensor_dy_on_xface(i,j,k,2,vel,dyi,
1980  bvxlo,bvxhi,bct,dlo,dhi,
1981  whi,wlo,jhip,jhim,jlop,jlom);
1982 
1983  int khip = k + flag(i ,j,k).isConnected(0,0, 1);
1984  int khim = k - flag(i ,j,k).isConnected(0,0,-1);
1985  int klop = k + flag(i-1,j,k).isConnected(0,0, 1);
1986  int klom = k - flag(i-1,j,k).isConnected(0,0,-1);
1987  whi = mlebtensor_weight(khip-khim);
1988  wlo = mlebtensor_weight(klop-klom);
1989  Real dudz = mlebtensor_dz_on_xface(i,j,k,0,vel,dzi,
1990  bvxlo,bvxhi,bct,dlo,dhi,
1991  whi,wlo,khip,khim,klop,klom);
1992  Real dvdz = mlebtensor_dz_on_xface(i,j,k,1,vel,dzi,
1993  bvxlo,bvxhi,bct,dlo,dhi,
1994  whi,wlo,khip,khim,klop,klom);
1995  Real dwdz = mlebtensor_dz_on_xface(i,j,k,2,vel,dzi,
1996  bvxlo,bvxhi,bct,dlo,dhi,
1997  whi,wlo,khip,khim,klop,klom);
1998 
1999  fx(i,j,k,0) = dudx;
2000  fx(i,j,k,1) = dvdx;
2001  fx(i,j,k,2) = dwdx;
2002  fx(i,j,k,3) = dudy;
2003  fx(i,j,k,4) = dvdy;
2004  fx(i,j,k,5) = dwdy;
2005  fx(i,j,k,6) = dudz;
2006  fx(i,j,k,7) = dvdz;
2007  fx(i,j,k,8) = dwdz;
2008  }
2009  }
2010  }
2011  }
2012 }
2013 
2014 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2015 void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
2016  Array4<Real const> const& vel,
2017  Array4<Real const> const& apy,
2018  Array4<EBCellFlag const> const& flag,
2019  Array4<int const> const& ccm,
2020  Array4<Real const> const& fcy,
2021  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2022  Array4<Real const> const& bvylo,
2023  Array4<Real const> const& bvyhi,
2024  Array2D<BoundCond,
2025  0,2*AMREX_SPACEDIM,
2026  0,AMREX_SPACEDIM> const& bct,
2027  Dim3 const& dlo, Dim3 const& dhi) noexcept
2028 {
2029  const Real dxi = dxinv[0];
2030  const Real dyi = dxinv[1];
2031  const Real dzi = dxinv[2];
2032  const auto lo = amrex::lbound(box);
2033  const auto hi = amrex::ubound(box);
2034 
2035  for (int k = lo.z; k <= hi.z; ++k) {
2036  for (int j = lo.y; j <= hi.y; ++j) {
2037  for (int i = lo.x; i <= hi.x; ++i) {
2038  if (apy(i,j,k) == Real(0.0))
2039  {
2040  fy(i,j,k,0) = Real(0.0);
2041  fy(i,j,k,1) = Real(0.0);
2042  fy(i,j,k,2) = Real(0.0);
2043  fy(i,j,k,3) = Real(0.0);
2044  fy(i,j,k,4) = Real(0.0);
2045  fy(i,j,k,5) = Real(0.0);
2046  fy(i,j,k,6) = Real(0.0);
2047  fy(i,j,k,7) = Real(0.0);
2048  fy(i,j,k,8) = Real(0.0);
2049  }
2050  else
2051  {
2052  int ihip = i + flag(i,j ,k).isConnected( 1,0,0);
2053  int ihim = i - flag(i,j ,k).isConnected(-1,0,0);
2054  int ilop = i + flag(i,j-1,k).isConnected( 1,0,0);
2055  int ilom = i - flag(i,j-1,k).isConnected(-1,0,0);
2056  Real whi = mlebtensor_weight(ihip-ihim);
2057  Real wlo = mlebtensor_weight(ilop-ilom);
2058  Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
2059  bvylo,bvyhi,bct,dlo,dhi,
2060  whi,wlo,ihip,ihim,ilop,ilom);
2061  Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
2062  bvylo,bvyhi,bct,dlo,dhi,
2063  whi,wlo,ihip,ihim,ilop,ilom);
2064  Real dwdx = mlebtensor_dx_on_yface(i,j,k,2,vel,dxi,
2065  bvylo,bvyhi,bct,dlo,dhi,
2066  whi,wlo,ihip,ihim,ilop,ilom);
2067 
2068  Real dudy = (vel(i,j,k,0) - vel(i,j-1,k,0))*dyi;
2069  Real dvdy = (vel(i,j,k,1) - vel(i,j-1,k,1))*dyi;
2070  Real dwdy = (vel(i,j,k,2) - vel(i,j-1,k,2))*dyi;
2071  if (apy(i,j,k) < Real(1.0)) {
2072  int ii = i + (fcy(i,j,k,0) >= Real(0.0) ? 1 : -1);
2073  int kk = k + (fcy(i,j,k,1) >= Real(0.0) ? 1 : -1);
2074  Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy(i,j,k,0)) : Real(0.0);
2075  Real fracz = (ccm(i,j-1,kk) || ccm(i,j,kk)) ? std::abs(fcy(i,j,k,1)) : Real(0.0);
2076  dudy = (Real(1.0)-fracx)*(Real(1.0)-fracz) * dudy
2077  + fracx*(Real(1.0)-fracz) * (vel(ii,j,k ,0) - vel(ii,j-1,k ,0))*dyi
2078  + fracz*(Real(1.0)-fracx) * (vel(i ,j,kk,0) - vel(i ,j-1,kk,0))*dyi
2079  + fracx* fracz * (vel(ii,j,kk,0) - vel(ii,j-1,kk,0))*dyi;
2080  dvdy = (Real(1.0)-fracx)*(Real(1.0)-fracz) * dvdy
2081  + fracx*(Real(1.0)-fracz) * (vel(ii,j,k ,1) - vel(ii,j-1,k ,1))*dyi
2082  + fracz*(Real(1.0)-fracx) * (vel(i ,j,kk,1) - vel(i ,j-1,kk,1))*dyi
2083  + fracx* fracz * (vel(ii,j,kk,1) - vel(ii,j-1,kk,1))*dyi;
2084  dwdy = (Real(1.0)-fracx)*(Real(1.0)-fracz) * dwdy
2085  + fracx*(Real(1.0)-fracz) * (vel(ii,j,k ,2) - vel(ii,j-1,k ,2))*dyi
2086  + fracz*(Real(1.0)-fracx) * (vel(i ,j,kk,2) - vel(i ,j-1,kk,2))*dyi
2087  + fracx* fracz * (vel(ii,j,kk,2) - vel(ii,j-1,kk,2))*dyi;
2088  }
2089 
2090  int khip = k + flag(i,j ,k).isConnected(0,0, 1);
2091  int khim = k - flag(i,j ,k).isConnected(0,0,-1);
2092  int klop = k + flag(i,j-1,k).isConnected(0,0, 1);
2093  int klom = k - flag(i,j-1,k).isConnected(0,0,-1);
2094  whi = mlebtensor_weight(khip-khim);
2095  wlo = mlebtensor_weight(klop-klom);
2096  Real dudz = mlebtensor_dz_on_yface(i,j,k,0,vel,dzi,
2097  bvylo,bvyhi,bct,dlo,dhi,
2098  whi,wlo,khip,khim,klop,klom);
2099  Real dvdz = mlebtensor_dz_on_yface(i,j,k,1,vel,dzi,
2100  bvylo,bvyhi,bct,dlo,dhi,
2101  whi,wlo,khip,khim,klop,klom);
2102  Real dwdz = mlebtensor_dz_on_yface(i,j,k,2,vel,dzi,
2103  bvylo,bvyhi,bct,dlo,dhi,
2104  whi,wlo,khip,khim,klop,klom);
2105 
2106  fy(i,j,k,0) = dudx;
2107  fy(i,j,k,1) = dvdx;
2108  fy(i,j,k,2) = dwdx;
2109  fy(i,j,k,3) = dudy;
2110  fy(i,j,k,4) = dvdy;
2111  fy(i,j,k,5) = dwdy;
2112  fy(i,j,k,6) = dudz;
2113  fy(i,j,k,7) = dvdz;
2114  fy(i,j,k,8) = dwdz;
2115  }
2116  }
2117  }
2118  }
2119 }
2120 
2121 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2122 void mlebtensor_vel_grads_fz (Box const& box, Array4<Real> const& fz,
2123  Array4<Real const> const& vel,
2124  Array4<Real const> const& apz,
2125  Array4<EBCellFlag const> const& flag,
2126  Array4<int const> const& ccm,
2127  Array4<Real const> const& fcz,
2128  GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2129  Array4<Real const> const& bvzlo,
2130  Array4<Real const> const& bvzhi,
2131  Array2D<BoundCond,
2132  0,2*AMREX_SPACEDIM,
2133  0,AMREX_SPACEDIM> const& bct,
2134  Dim3 const& dlo, Dim3 const& dhi) noexcept
2135 {
2136  const Real dxi = dxinv[0];
2137  const Real dyi = dxinv[1];
2138  const Real dzi = dxinv[2];
2139  const auto lo = amrex::lbound(box);
2140  const auto hi = amrex::ubound(box);
2141 
2142  for (int k = lo.z; k <= hi.z; ++k) {
2143  for (int j = lo.y; j <= hi.y; ++j) {
2144  for (int i = lo.x; i <= hi.x; ++i) {
2145  if (apz(i,j,k) == Real(0.0))
2146  {
2147  fz(i,j,k,0) = Real(0.0);
2148  fz(i,j,k,1) = Real(0.0);
2149  fz(i,j,k,2) = Real(0.0);
2150  fz(i,j,k,3) = Real(0.0);
2151  fz(i,j,k,4) = Real(0.0);
2152  fz(i,j,k,5) = Real(0.0);
2153  fz(i,j,k,6) = Real(0.0);
2154  fz(i,j,k,7) = Real(0.0);
2155  fz(i,j,k,8) = Real(0.0);
2156  }
2157  else
2158  {
2159  int ihip = i + flag(i,j,k ).isConnected( 1,0,0);
2160  int ihim = i - flag(i,j,k ).isConnected(-1,0,0);
2161  int ilop = i + flag(i,j,k-1).isConnected( 1,0,0);
2162  int ilom = i - flag(i,j,k-1).isConnected(-1,0,0);
2163  Real whi = mlebtensor_weight(ihip-ihim);
2164  Real wlo = mlebtensor_weight(ilop-ilom);
2165  Real dudx = mlebtensor_dx_on_zface(i,j,k,0,vel,dxi,
2166  bvzlo,bvzhi,bct,dlo,dhi,
2167  whi,wlo,ihip,ihim,ilop,ilom);
2168  Real dvdx = mlebtensor_dx_on_zface(i,j,k,1,vel,dxi,
2169  bvzlo,bvzhi,bct,dlo,dhi,
2170  whi,wlo,ihip,ihim,ilop,ilom);
2171  Real dwdx = mlebtensor_dx_on_zface(i,j,k,2,vel,dxi,
2172  bvzlo,bvzhi,bct,dlo,dhi,
2173  whi,wlo,ihip,ihim,ilop,ilom);
2174 
2175  int jhip = j + flag(i,j,k ).isConnected(0, 1,0);
2176  int jhim = j - flag(i,j,k ).isConnected(0,-1,0);
2177  int jlop = j + flag(i,j,k-1).isConnected(0, 1,0);
2178  int jlom = j - flag(i,j,k-1).isConnected(0,-1,0);
2179  whi = mlebtensor_weight(jhip-jhim);
2180  wlo = mlebtensor_weight(jlop-jlom);
2181  Real dudy = mlebtensor_dy_on_zface(i,j,k,0,vel,dyi,
2182  bvzlo,bvzhi,bct,dlo,dhi,
2183  whi,wlo,jhip,jhim,jlop,jlom);
2184  Real dvdy = mlebtensor_dy_on_zface(i,j,k,1,vel,dyi,
2185  bvzlo,bvzhi,bct,dlo,dhi,
2186  whi,wlo,jhip,jhim,jlop,jlom);
2187  Real dwdy = mlebtensor_dy_on_zface(i,j,k,2,vel,dyi,
2188  bvzlo,bvzhi,bct,dlo,dhi,
2189  whi,wlo,jhip,jhim,jlop,jlom);
2190 
2191  Real dudz = (vel(i,j,k,0) - vel(i,j,k-1,0))*dzi;
2192  Real dvdz = (vel(i,j,k,1) - vel(i,j,k-1,1))*dzi;
2193  Real dwdz = (vel(i,j,k,2) - vel(i,j,k-1,2))*dzi;
2194  if (apz(i,j,k) < Real(1.0)) {
2195  int ii = i + (fcz(i,j,k,0) >= Real(0.0) ? 1 : -1);
2196  int jj = j + (fcz(i,j,k,1) >= Real(0.0) ? 1 : -1);
2197  Real fracx = (ccm(ii,j,k-1) || ccm(ii,j,k)) ? std::abs(fcz(i,j,k,0)) : Real(0.0);
2198  Real fracy = (ccm(i,jj,k-1) || ccm(i,jj,k)) ? std::abs(fcz(i,j,k,1)) : Real(0.0);
2199  dudz = (Real(1.0)-fracx)*(Real(1.0)-fracy) * dudz
2200  + fracx*(Real(1.0)-fracy) * (vel(ii,j ,k,0) - vel(ii,j ,k-1,0))*dzi
2201  + fracy*(Real(1.0)-fracx) * (vel(i ,jj,k,0) - vel(i ,jj,k-1,0))*dzi
2202  + fracx* fracy * (vel(ii,jj,k,0) - vel(ii,jj,k-1,0))*dzi;
2203  dvdz = (Real(1.0)-fracx)*(Real(1.0)-fracy) * dvdz
2204  + fracx*(Real(1.0)-fracy) * (vel(ii,j ,k,1) - vel(ii,j ,k-1,1))*dzi
2205  + fracy*(Real(1.0)-fracx) * (vel(i ,jj,k,1) - vel(i ,jj,k-1,1))*dzi
2206  + fracx* fracy * (vel(ii,jj,k,1) - vel(ii,jj,k-1,1))*dzi;
2207  dwdz = (Real(1.0)-fracx)*(Real(1.0)-fracy) * dwdz
2208  + fracx*(Real(1.0)-fracy) * (vel(ii,j ,k,2) - vel(ii,j ,k-1,2))*dzi
2209  + fracy*(Real(1.0)-fracx) * (vel(i ,jj,k,2) - vel(i ,jj,k-1,2))*dzi
2210  + fracx* fracy * (vel(ii,jj,k,2) - vel(ii,jj,k-1,2))*dzi;
2211  }
2212 
2213  fz(i,j,k,0) = dudx;
2214  fz(i,j,k,1) = dvdx;
2215  fz(i,j,k,2) = dwdx;
2216  fz(i,j,k,3) = dudy;
2217  fz(i,j,k,4) = dvdy;
2218  fz(i,j,k,5) = dwdy;
2219  fz(i,j,k,6) = dudz;
2220  fz(i,j,k,7) = dvdz;
2221  fz(i,j,k,8) = dwdz;
2222  }
2223  }
2224  }
2225  }
2226 }
2227 
2228 }
2229 
2230 #endif
AMREX_PRAGMA_SIMD
#define AMREX_PRAGMA_SIMD
Definition: AMReX_Extension.H:80
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
AMREX_LO_NEUMANN
#define AMREX_LO_NEUMANN
Definition: AMReX_LO_BCTYPES.H:6
AMREX_LO_DIRICHLET
#define AMREX_LO_DIRICHLET
Definition: AMReX_LO_BCTYPES.H:5
AMReX_MLEBABecLap_K.H
amrex::BoundCond
Maintain an identifier for boundary condition types.
Definition: AMReX_BoundCond.H:20
amrex::Orientation::zhi
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE int zhi() noexcept
Int value of the z-hi-face.
Definition: AMReX_Orientation.H:118
amrex::Orientation::yhi
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE int yhi() noexcept
Int value of the y-hi-face.
Definition: AMReX_Orientation.H:110
amrex::Orientation::ylo
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE int ylo() noexcept
Int value of the y-lo-face.
Definition: AMReX_Orientation.H:106
amrex::Orientation::xlo
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE int xlo() noexcept
Int value of the x-lo-face.
Definition: AMReX_Orientation.H:98
amrex::Orientation::xhi
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE int xhi() noexcept
Int value of the x-hi-face.
Definition: AMReX_Orientation.H:102
amrex::Orientation::zlo
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE int zlo() noexcept
Int value of the z-lo-face.
Definition: AMReX_Orientation.H:114
amrex
Definition: AMReX_Amr.cpp:49
amrex::mlebtensor_dx_on_zface
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_dx_on_zface(int, int j, int k, int n, Array4< Real const > const &vel, Real dxi, Real whi, Real wlo, int ihip, int ihim, int ilop, int ilom) noexcept
Definition: AMReX_MLEBTensor_3D_K.H:30
amrex::mlebtensor_weight
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_weight(int d)
Definition: AMReX_MLEBTensor_K.H:10
amrex::Box
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
amrex::max
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
amrex::mlebtensor_flux_z
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_flux_z(Box const &box, Array4< Real > const &Az, Array4< Real const > const &fz, Array4< Real const > const &apz, Array4< Real const > const &fcz, Real const bscalar, Array4< int const > const &ccm, int face_only, Box const &zbox) noexcept
Definition: AMReX_MLEBTensor_3D_K.H:1089
amrex::mlebtensor_dz_on_yface
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_dz_on_yface(int i, int j, int, int n, Array4< Real const > const &vel, Real dzi, Real whi, Real wlo, int khip, int khim, int klop, int klom) noexcept
Definition: AMReX_MLEBTensor_3D_K.H:20
amrex::mlebtensor_vel_grads_fx
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_vel_grads_fx(Box const &box, Array4< Real > const &fx, Array4< Real const > const &vel, Array4< Real const > const &apx, Array4< EBCellFlag const > const &flag, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:453
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::mlebtensor_flux_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_flux_y(Box const &box, Array4< Real > const &Ay, Array4< Real const > const &fy, Array4< Real const > const &apy, Array4< Real const > const &fcy, Real const bscalar, Array4< int const > const &ccm, int face_only, Box const &ybox) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:419
amrex::mlebtensor_cross_terms_fy
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_cross_terms_fy(Box const &box, Array4< Real > const &fy, Array4< Real const > const &vel, Array4< Real const > const &etay, Array4< Real const > const &kapy, Array4< Real const > const &apy, Array4< EBCellFlag const > const &flag, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:56
amrex::mlebtensor_cross_terms
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_cross_terms(Box const &box, Array4< Real > const &Ax, Array4< Real const > const &fx, Array4< Real const > const &fy, Array4< Real const > const &vel, Array4< Real const > const &velb, Array4< Real const > const &etab, Array4< Real const > const &kapb, Array4< int const > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vol, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &bc, bool is_dirichlet, bool is_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Real bscalar) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:210
amrex::mlebtensor_vel_grads_fy
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_vel_grads_fy(Box const &box, Array4< Real > const &fy, Array4< Real const > const &vel, Array4< Real const > const &apy, Array4< EBCellFlag const > const &flag, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:499
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::mlebtensor_dx_on_yface
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_dx_on_yface(int, int j, int k, int n, Array4< Real const > const &vel, Real dxi, Real whi, Real wlo, int ihip, int ihim, int ilop, int ilom) noexcept
Definition: AMReX_MLEBTensor_K.H:25
amrex::mlebtensor_dy_on_xface
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_dy_on_xface(int i, int, int k, int n, Array4< Real const > const &vel, Real dyi, Real whi, Real wlo, int jhip, int jhim, int jlop, int jlom) noexcept
Definition: AMReX_MLEBTensor_K.H:15
amrex::mlebtensor_cross_terms_fx
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_cross_terms_fx(Box const &box, Array4< Real > const &fx, Array4< Real const > const &vel, Array4< Real const > const &etax, Array4< Real const > const &kapx, Array4< Real const > const &apx, Array4< EBCellFlag const > const &flag, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:10
amrex::mlebtensor_cross_terms_fz
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_cross_terms_fz(Box const &box, Array4< Real > const &fz, Array4< Real const > const &vel, Array4< Real const > const &etaz, Array4< Real const > const &kapz, Array4< Real const > const &apz, Array4< EBCellFlag const > const &flag, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBTensor_3D_K.H:170
amrex::mlebtensor_flux_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_flux_x(Box const &box, Array4< Real > const &Ax, Array4< Real const > const &fx, Array4< Real const > const &apx, Array4< Real const > const &fcx, Real const bscalar, Array4< int const > const &ccm, int face_only, Box const &xbox) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:385
amrex::LoopConcurrent
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void LoopConcurrent(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:150
amrex::mlebtensor_dy_on_zface
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_dy_on_zface(int i, int, int k, int n, Array4< Real const > const &vel, Real dyi, Real whi, Real wlo, int jhip, int jhim, int jlop, int jlom) noexcept
Definition: AMReX_MLEBTensor_3D_K.H:40
amrex::mlebtensor_dz_on_xface
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_dz_on_xface(int i, int j, int, int n, Array4< Real const > const &vel, Real dzi, Real whi, Real wlo, int khip, int khim, int klop, int klom) noexcept
Definition: AMReX_MLEBTensor_3D_K.H:10
amrex::mlebtensor_flux_0
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_flux_0(Box const &box, Array4< Real > const &Ax, Array4< Real const > const &fx, Array4< Real const > const &ap, Real bscalar) noexcept
Definition: AMReX_MLEBTensor_2D_K.H:360
amrex::mlebtensor_vel_grads_fz
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebtensor_vel_grads_fz(Box const &box, Array4< Real > const &fz, Array4< Real const > const &vel, Array4< Real const > const &apz, Array4< EBCellFlag const > const &flag, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBTensor_3D_K.H:1287
amrex::sqrt
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::Array2D
Definition: AMReX_Array.H:282
amrex::Array4
Definition: AMReX_Array4.H:61
amrex::Dim3
Definition: AMReX_Dim3.H:12