Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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
7namespace amrex {
8
9AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
10Real 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
19AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
20Real 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
29AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
30Real 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
39AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
40Real 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
49AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
50void 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
109AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
110void 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
169AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
170void 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
230AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
231Real 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
285AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
286Real 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
340AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
341Real 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
396AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
397Real 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
451AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
452void 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
522AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
523void 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
592AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
593void 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
663AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
664void 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
968AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
969void 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
994AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
995void 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
1041AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1042void 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
1088AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1089void 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
1136AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1137void 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
1211AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1212void 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
1286AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1287void 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
1361AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1362void 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
1448AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1449void 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
1535AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1536void 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
1622AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1623void 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
1717AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1718void 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
1812AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1813void 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
1907AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1908void 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
2014AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2015void 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
2121AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
2122void 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::yhi
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr int yhi() noexcept
Int value of the y-hi-face.
Definition AMReX_Orientation.H:110
amrex::Orientation::zlo
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr int zlo() noexcept
Int value of the z-lo-face.
Definition AMReX_Orientation.H:114
amrex::Orientation::xhi
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr int xhi() noexcept
Int value of the x-hi-face.
Definition AMReX_Orientation.H:102
amrex::Orientation::xlo
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr int xlo() noexcept
Int value of the x-lo-face.
Definition AMReX_Orientation.H:98
amrex::Orientation::zhi
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr int zhi() noexcept
Int value of the z-hi-face.
Definition AMReX_Orientation.H:118
amrex::Orientation::ylo
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr int ylo() noexcept
Int value of the y-lo-face.
Definition AMReX_Orientation.H:106
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::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::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::max
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:35
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::Array2D
Definition AMReX_Array.H:282
amrex::Array4
Definition AMReX_Array4.H:61
amrex::Dim3
Definition AMReX_Dim3.H:12
amrex::GpuArray
Definition AMReX_Array.H:34