Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLEBTensor_2D_K.H
Go to the documentation of this file.
1#ifndef AMREX_ML_EB_TENSOR_2D_K_H_
2#define AMREX_ML_EB_TENSOR_2D_K_H_
3#include <AMReX_Config.H>
4
6
7namespace amrex {
8
10void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
11 Array4<Real const> const& vel,
12 Array4<Real const> const& etax,
13 Array4<Real const> const& kapx,
14 Array4<Real const> const& apx,
15 Array4<EBCellFlag const> const& flag,
16 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
17{
18 const Real dyi = dxinv[1];
19 const auto lo = amrex::lbound(box);
20 const auto hi = amrex::ubound(box);
21 constexpr Real twoThirds = 2./3.;
22
23 int k = 0;
24 for (int j = lo.y; j <= hi.y; ++j) {
26 for (int i = lo.x; i <= hi.x; ++i) {
27 if (apx(i,j,0) == Real(0.0))
28 {
29 fx(i,j,0,0) = Real(0.0);
30 fx(i,j,0,1) = Real(0.0);
31 }
32 else
33 {
34 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
35 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
36 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
37 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
38 Real whi = mlebtensor_weight(jhip-jhim);
39 Real wlo = mlebtensor_weight(jlop-jlom);
40 Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
41 whi,wlo,jhip,jhim,jlop,jlom);
42 Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
43 whi,wlo,jhip,jhim,jlop,jlom);
44 Real divu = dvdy;
45 Real xif = kapx(i,j,0);
46 Real mun = Real(0.75)*(etax(i,j,0,0)-xif);// restore the original eta
47 Real mut = etax(i,j,0,1);
48 fx(i,j,0,0) = -mun*(-twoThirds*divu) - xif*divu;
49 fx(i,j,0,1) = -mut*dudy;
50 }
51 }
52 }
53}
54
56void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
57 Array4<Real const> const& vel,
58 Array4<Real const> const& etay,
59 Array4<Real const> const& kapy,
60 Array4<Real const> const& apy,
61 Array4<EBCellFlag const> const& flag,
62 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
63{
64 const Real dxi = dxinv[0];
65 const auto lo = amrex::lbound(box);
66 const auto hi = amrex::ubound(box);
67 constexpr Real twoThirds = 2./3.;
68
69 int k = 0;
70 for (int j = lo.y; j <= hi.y; ++j) {
72 for (int i = lo.x; i <= hi.x; ++i) {
73 if (apy(i,j,0) == Real(0.0))
74 {
75 fy(i,j,0,0) = Real(0.0);
76 fy(i,j,0,1) = Real(0.0);
77 }
78 else
79 {
80 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
81 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
82 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
83 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
84 Real whi = mlebtensor_weight(ihip-ihim);
85 Real wlo = mlebtensor_weight(ilop-ilom);
86 Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
87 whi,wlo,ihip,ihim,ilop,ilom);
88 Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
89 whi,wlo,ihip,ihim,ilop,ilom);
90 Real divu = dudx;
91 Real xif = kapy(i,j,0);
92 Real mun = Real(0.75)*(etay(i,j,0,1)-xif);// restore the original eta
93 Real mut = etay(i,j,0,0);
94 fy(i,j,0,0) = -mut*dvdx;
95 fy(i,j,0,1) = -mun*(-twoThirds*divu) - xif*divu;
96 }
97 }
98 }
99}
100
102void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
103 Array4<Real const> const& vel,
104 Array4<Real const> const& etax,
105 Array4<Real const> const& kapx,
106 Array4<Real const> const& apx,
107 Array4<EBCellFlag const> const& flag,
109 Array4<Real const> const& bvxlo,
110 Array4<Real const> const& bvxhi,
112 0,2*AMREX_SPACEDIM,
113 0,AMREX_SPACEDIM> const& bct,
114 Dim3 const& dlo, Dim3 const& dhi) noexcept
115{
116 const Real dyi = dxinv[1];
117 const auto lo = amrex::lbound(box);
118 const auto hi = amrex::ubound(box);
119 constexpr Real twoThirds = 2./3.;
120
121 int k = 0;
122 for (int j = lo.y; j <= hi.y; ++j) {
124 for (int i = lo.x; i <= hi.x; ++i) {
125 if (apx(i,j,0) == Real(0.0))
126 {
127 fx(i,j,0,0) = Real(0.0);
128 fx(i,j,0,1) = Real(0.0);
129 }
130 else
131 {
132 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
133 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
134 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
135 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
136 Real whi = mlebtensor_weight(jhip-jhim);
137 Real wlo = mlebtensor_weight(jlop-jlom);
138 Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
139 bvxlo,bvxhi,bct,dlo,dhi,
140 whi,wlo,jhip,jhim,jlop,jlom);
141 Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
142 bvxlo,bvxhi,bct,dlo,dhi,
143 whi,wlo,jhip,jhim,jlop,jlom);
144 Real divu = dvdy;
145 Real xif = kapx(i,j,0);
146 Real mun = Real(0.75)*(etax(i,j,0,0)-xif);// restore the original eta
147 Real mut = etax(i,j,0,1);
148 fx(i,j,0,0) = -mun*(-twoThirds*divu) - xif*divu;
149 fx(i,j,0,1) = -mut*dudy;
150 }
151 }
152 }
153}
154
156void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
157 Array4<Real const> const& vel,
158 Array4<Real const> const& etay,
159 Array4<Real const> const& kapy,
160 Array4<Real const> const& apy,
161 Array4<EBCellFlag const> const& flag,
163 Array4<Real const> const& bvylo,
164 Array4<Real const> const& bvyhi,
166 0,2*AMREX_SPACEDIM,
167 0,AMREX_SPACEDIM> const& bct,
168 Dim3 const& dlo, Dim3 const& dhi) noexcept
169{
170 const Real dxi = dxinv[0];
171 const auto lo = amrex::lbound(box);
172 const auto hi = amrex::ubound(box);
173 constexpr Real twoThirds = 2./3.;
174
175 int k = 0;
176 for (int j = lo.y; j <= hi.y; ++j) {
178 for (int i = lo.x; i <= hi.x; ++i) {
179 if (apy(i,j,0) == Real(0.0))
180 {
181 fy(i,j,0,0) = Real(0.0);
182 fy(i,j,0,1) = Real(0.0);
183 }
184 else
185 {
186 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
187 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
188 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
189 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
190 Real whi = mlebtensor_weight(ihip-ihim);
191 Real wlo = mlebtensor_weight(ilop-ilom);
192 Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
193 bvylo,bvyhi,bct,dlo,dhi,
194 whi,wlo,ihip,ihim,ilop,ilom);
195 Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
196 bvylo,bvyhi,bct,dlo,dhi,
197 whi,wlo,ihip,ihim,ilop,ilom);
198 Real divu = dudx;
199 Real xif = kapy(i,j,0);
200 Real mun = Real(0.75)*(etay(i,j,0,1)-xif);// restore the original eta
201 Real mut = etay(i,j,0,0);
202 fy(i,j,0,0) = -mut*dvdx;
203 fy(i,j,0,1) = -mun*(-twoThirds*divu) - xif*divu;
204 }
205 }
206 }
207}
208
210void mlebtensor_cross_terms (Box const& box, Array4<Real> const& Ax,
211 Array4<Real const> const& fx,
212 Array4<Real const> const& fy,
213 Array4<Real const> const& vel,
214 Array4<Real const> const& velb,
215 Array4<Real const> const& etab,
216 Array4<Real const> const& kapb,
217 Array4<int const> const& ccm,
218 Array4<EBCellFlag const> const& flag,
219 Array4<Real const> const& vol,
220 Array4<Real const> const& apx,
221 Array4<Real const> const& apy,
222 Array4<Real const> const& fcx,
223 Array4<Real const> const& fcy,
224 Array4<Real const> const& bc,
225 bool is_dirichlet, bool is_inhomog,
227 Real bscalar) noexcept
228{
229 const Real dxi = dxinv[0];
230 const Real dyi = dxinv[1];
231 const auto lo = amrex::lbound(box);
232 const auto hi = amrex::ubound(box);
233
234 for (int j = lo.y; j <= hi.y; ++j) {
236 for (int i = lo.x; i <= hi.x; ++i) {
237 if (flag(i,j,0).isRegular())
238 {
239 Ax(i,j,0,0) += bscalar*(dxi*(fx(i+1,j ,0,0) - fx(i,j,0,0))
240 + dyi*(fy(i ,j+1,0,0) - fy(i,j,0,0)));
241 Ax(i,j,0,1) += bscalar*(dxi*(fx(i+1,j ,0,1) - fx(i,j,0,1))
242 + dyi*(fy(i ,j+1,0,1) - fy(i,j,0,1)));
243 }
244 else if (flag(i,j,0).isSingleValued())
245 {
246 Real fxm_0 = fx(i,j,0,0);
247 Real fxm_1 = fx(i,j,0,1);
248 if (apx(i,j,0) > Real(0.0) && apx(i,j,0) < Real(1.0)) {
249 int jj = j + (fcx(i,j,0,0) >= Real(0.0) ? 1 : -1);
250 Real fracy = (ccm(i-1,jj,0) || ccm(i,jj,0)) ? std::abs(fcx(i,j,0,0)) : Real(0.0);
251 fxm_0 = (Real(1.0)-fracy)*fxm_0 + fracy*fx(i,jj,0,0);
252 fxm_1 = (Real(1.0)-fracy)*fxm_1 + fracy*fx(i,jj,0,1);
253 }
254
255 Real fxp_0 = fx(i+1,j,0,0);
256 Real fxp_1 = fx(i+1,j,0,1);
257 if (apx(i+1,j,0) > Real(0.0) && apx(i+1,j,0) < Real(1.0)) {
258 int jj = j + (fcx(i+1,j,0,0) >= Real(0.0) ? 1 : -1);
259 Real fracy = (ccm(i,jj,0) || ccm(i+1,jj,0)) ? std::abs(fcx(i+1,j,0,0)) : Real(0.0);
260 fxp_0 = (Real(1.0)-fracy)*fxp_0 + fracy*fx(i+1,jj,0,0);
261 fxp_1 = (Real(1.0)-fracy)*fxp_1 + fracy*fx(i+1,jj,0,1);
262 }
263
264 Real fym_0 = fy(i,j,0,0);
265 Real fym_1 = fy(i,j,0,1);
266 if (apy(i,j,0) > Real(0.0) && apy(i,j,0) < Real(1.0)) {
267 int ii = i + (fcy(i,j,0,0) >= Real(0.0) ? 1 : -1);
268 Real fracx = (ccm(ii,j-1,0) || ccm(ii,j,0)) ? std::abs(fcy(i,j,0,0)) : Real(0.0);
269 fym_0 = (Real(1.0)-fracx)*fym_0 + fracx*fy(ii,j,0,0);
270 fym_1 = (Real(1.0)-fracx)*fym_1 + fracx*fy(ii,j,0,1);
271 }
272
273 Real fyp_0 = fy(i,j+1,0,0);
274 Real fyp_1 = fy(i,j+1,0,1);
275 if (apy(i,j+1,0) > Real(0.0) && apy(i,j+1,0) < Real(1.0)) {
276 int ii = i + (fcy(i,j+1,0,0) >= Real(0.0) ? 1 : -1);
277 Real fracx = (ccm(ii,j,0) || ccm(ii,j+1,0)) ? std::abs(fcy(i,j+1,0,0)) : Real(0.0);
278 fyp_0 = (Real(1.0)-fracx)*fyp_0 + fracx*fy(ii,j+1,0,0);
279 fyp_1 = (Real(1.0)-fracx)*fyp_1 + fracx*fy(ii,j+1,0,1);
280 }
281
282 Real kappa = vol(i,j,0);
283 Real feb_0 = Real(0.0);
284 Real feb_1 = Real(0.0);
285 if (is_dirichlet)
286 {
287 Real dapx = apx(i+1,j,0)-apx(i,j,0);
288 Real dapy = apy(i,j+1,0)-apy(i,j,0);
289 Real anorminv = Real(1.0)/std::sqrt(dapx*dapx+dapy*dapy);
290 Real anrmx = -dapx * anorminv;
291 Real anrmy = -dapy * anorminv;
292
293 Real velb_0 = 0, velb_1 = 0;
294
295 if (is_inhomog) {
296 velb_0 = velb(i,j,0,0);
297 velb_1 = velb(i,j,0,1);
298 }
299
300 Real dx_eb = amrex::max(Real(0.3), (kappa*kappa-Real(0.25))/(Real(2.)*kappa));
301 Real dg = dx_eb / amrex::max(std::abs(anrmx),std::abs(anrmy));
302 Real dginv = Real(1.0)/dg;
303 Real gx = bc(i,j,0,0) - dg*anrmx;
304 Real gy = bc(i,j,0,1) - dg*anrmy;
305 int isx = (anrmx > Real(0.0)) ? 1 : -1;
306 int isy = (anrmy > Real(0.0)) ? 1 : -1;
307 int ii = i - isx;
308 int jj = j - isy;
309
310 gx *= isx;
311 gy *= isy;
312 Real gxy = gx*gy;
313 Real oneggg = Real(1.0)+gx+gy+gxy;
314
315 Real velg = oneggg * vel(i ,j ,0,0)
316 + (-gy - gxy) * vel(i ,jj,0,0)
317 + (-gx - gxy) * vel(ii,j ,0,0)
318 + gxy * vel(ii,jj,0,0);
319 Real dudn = (velb_0-velg) * dginv;
320
321 velg = oneggg * vel(i ,j ,0,1)
322 + (-gy - gxy) * vel(i ,jj,0,1)
323 + (-gx - gxy) * vel(ii,j ,0,1)
324 + gxy * vel(ii,jj,0,1);
325 Real dvdn = (velb_1-velg) * dginv;
326
327 // transverse derivatives are zero on the wall
328 Real dudx = dudn * anrmx;
329 Real dudy = dudn * anrmy;
330 Real dvdx = dvdn * anrmx;
331 Real dvdy = dvdn * anrmy;
332 Real divu = dudx + dvdy;
333 Real xi = kapb(i,j,0);
334 Real mu = etab(i,j,0);
335 Real tautmp = (xi-(2./3.)*mu)*divu;
336 // Note that mu*(grad v) has been included already in MLEBABecLap
337 Real tauxx = mu*dudx + tautmp;
338 Real tauyx = mu*dvdx;
339 Real tauyy = mu*dvdy + tautmp;
340 Real tauxy = mu*dudy;
341 // assuming dxi == dyi == dzi
342 feb_0 = dxi*(dapx*tauxx + dapy*tauyx);
343 feb_1 = dxi*(dapx*tauxy + dapy*tauyy);
344 }
345
346 Real volinv = bscalar / kappa;
347 Ax(i,j,0,0) += volinv * (dxi*(apx(i,j,0)*fxm_0-apx(i+1,j,0)*fxp_0)
348 +dyi*(apy(i,j,0)*fym_0-apy(i,j+1,0)*fyp_0)
349 -dxi* feb_0);
350 Ax(i,j,0,1) += volinv * (dxi*(apx(i,j,0)*fxm_1-apx(i+1,j,0)*fxp_1)
351 +dyi*(apy(i,j,0)*fym_1-apy(i,j+1,0)*fyp_1)
352 -dxi* feb_1);
353 }
354 }
355 }
356}
357
358
360void mlebtensor_flux_0 (Box const& box,
361 Array4<Real> const& Ax,
362 Array4<Real const> const& fx,
363 Array4<Real const> const& ap,
364 Real bscalar) noexcept
365{
366 const auto lo = amrex::lbound(box);
367 const auto hi = amrex::ubound(box);
368
369 int k = 0;
370 for (int j = lo.y; j <= hi.y; ++j) {
372 for (int i = lo.x; i <= hi.x; ++i) {
373 if (ap(i,j,k) != Real(0.0)) {
374 for (int n=0; n<AMREX_SPACEDIM; n++) {
375 Ax(i,j,k,n) += bscalar*fx(i,j,k,n);
376 }
377 }
378 //else MLEBABec::compFlux already set Ax = 0
379 }
380 }
381}
382
383
385void mlebtensor_flux_x (Box const& box, Array4<Real> const& Ax,
386 Array4<Real const> const& fx, Array4<Real const> const& apx,
387 Array4<Real const> const& fcx,
388 Real const bscalar, Array4<int const> const& ccm,
389 int face_only, Box const& xbox) noexcept
390{
391 int lof = xbox.smallEnd(0);
392 int hif = xbox.bigEnd(0);
393 amrex::LoopConcurrent(box, [=] (int i, int j, int k) noexcept
394 {
395 if (!face_only || lof == i || hif == i) {
396 if (apx(i,j,k) == Real(1.0)) {
397 for (int n=0; n<AMREX_SPACEDIM; n++) {
398 Ax(i,j,k,n) += bscalar*fx(i,j,k,n);
399 }
400 }
401 else if (apx(i,j,k) != 0.) {
402 Real fxm_0 = fx(i,j,0,0);
403 Real fxm_1 = fx(i,j,0,1);
404
405 int jj = j + (fcx(i,j,0,0) >= Real(0.0) ? 1 : -1);
406 Real fracy = (ccm(i-1,jj,0) || ccm(i,jj,0)) ? std::abs(fcx(i,j,0,0)) : Real(0.0);
407 fxm_0 = (Real(1.0)-fracy)*fxm_0 + fracy*fx(i,jj,0,0);
408 fxm_1 = (Real(1.0)-fracy)*fxm_1 + fracy*fx(i,jj,0,1);
409
410 Ax(i,j,k,0) += bscalar*fxm_0;
411 Ax(i,j,k,1) += bscalar*fxm_1;
412 }
413 //else MLEBABec::compFlux already set Ax = 0
414 }
415 });
416}
417
419void mlebtensor_flux_y (Box const& box, Array4<Real> const& Ay,
420 Array4<Real const> const& fy, Array4<Real const> const& apy,
421 Array4<Real const> const& fcy,
422 Real const bscalar, Array4<int const> const& ccm,
423 int face_only, Box const& ybox) noexcept
424{
425 int lof = ybox.smallEnd(1);
426 int hif = ybox.bigEnd(1);
427 amrex::LoopConcurrent(box, [=] (int i, int j, int k) noexcept
428 {
429 if (!face_only || lof == j || hif == j) {
430 if (apy(i,j,k) == Real(1.0)) {
431 for (int n=0; n<AMREX_SPACEDIM; n++) {
432 Ay(i,j,k,n) += bscalar*fy(i,j,k,n);
433 }
434 }
435 else if (apy(i,j,k) != 0.) {
436 Real fym_0 = fy(i,j,0,0);
437 Real fym_1 = fy(i,j,0,1);
438
439 int ii = i + (fcy(i,j,0,0) >= Real(0.0) ? 1 : -1);
440 Real fracx = (ccm(ii,j-1,0) || ccm(ii,j,0)) ? std::abs(fcy(i,j,0,0)) : Real(0.0);
441 fym_0 = (Real(1.0)-fracx)*fym_0 + fracx*fy(ii,j,0,0);
442 fym_1 = (Real(1.0)-fracx)*fym_1 + fracx*fy(ii,j,0,1);
443
444 Ay(i,j,k,0) += bscalar*fym_0;
445 Ay(i,j,k,1) += bscalar*fym_1;
446 }
447 //else MLEBABec::compFlux already set Ay = 0
448 }
449 });
450}
451
453void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
454 Array4<Real const> const& vel, Array4<Real const> const& apx,
455 Array4<EBCellFlag const> const& flag,
456 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
457{
458 const Real dxi = dxinv[0];
459 const Real dyi = dxinv[1];
460 const auto lo = amrex::lbound(box);
461 const auto hi = amrex::ubound(box);
462
463 int k = 0;
464 for (int j = lo.y; j <= hi.y; ++j) {
466 for (int i = lo.x; i <= hi.x; ++i) {
467 if (apx(i,j,0) == Real(0.0))
468 {
469 fx(i,j,0,0) = Real(0.0);
470 fx(i,j,0,1) = Real(0.0);
471 fx(i,j,0,2) = Real(0.0);
472 fx(i,j,0,3) = Real(0.0);
473 }
474 else
475 {
476 Real dudx = (vel(i,j,0,0) - vel(i-1,j,0,0))*dxi;
477 Real dvdx = (vel(i,j,0,1) - vel(i-1,j,0,1))*dxi;
478
479 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
480 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
481 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
482 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
483 Real whi = mlebtensor_weight(jhip-jhim);
484 Real wlo = mlebtensor_weight(jlop-jlom);
485 Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
486 whi,wlo,jhip,jhim,jlop,jlom);
487 Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
488 whi,wlo,jhip,jhim,jlop,jlom);
489 fx(i,j,0,0) = dudx;
490 fx(i,j,0,1) = dvdx;
491 fx(i,j,0,2) = dudy;
492 fx(i,j,0,3) = dvdy;
493 }
494 }
495 }
496}
497
499void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
500 Array4<Real const> const& vel, Array4<Real const> const& apy,
501 Array4<EBCellFlag const> const& flag,
502 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
503{
504 const Real dxi = dxinv[0];
505 const Real dyi = dxinv[1];
506 const auto lo = amrex::lbound(box);
507 const auto hi = amrex::ubound(box);
508
509 int k = 0;
510 for (int j = lo.y; j <= hi.y; ++j) {
512 for (int i = lo.x; i <= hi.x; ++i) {
513 if (apy(i,j,0) == 0.) {
514 fy(i,j,0,0) = Real(0.0);
515 fy(i,j,0,1) = Real(0.0);
516 fy(i,j,0,2) = Real(0.0);
517 fy(i,j,0,3) = Real(0.0);
518 }
519 else
520 {
521 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
522 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
523 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
524 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
525 Real whi = mlebtensor_weight(ihip-ihim);
526 Real wlo = mlebtensor_weight(ilop-ilom);
527 Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
528 whi,wlo,ihip,ihim,ilop,ilom);
529 Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
530 whi,wlo,ihip,ihim,ilop,ilom);
531
532 Real dudy = (vel(i,j,0,0) - vel(i,j-1,0,0))*dyi;
533 Real dvdy = (vel(i,j,0,1) - vel(i,j-1,0,1))*dyi;
534 fy(i,j,0,0) = dudx;
535 fy(i,j,0,1) = dvdx;
536 fy(i,j,0,2) = dudy;
537 fy(i,j,0,3) = dvdy;
538 }
539 }
540 }
541}
542
544void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
545 Array4<Real const> const& vel, Array4<Real const> const& apx,
546 Array4<EBCellFlag const> const& flag,
548 Array4<Real const> const& bvxlo,
549 Array4<Real const> const& bvxhi,
551 0,2*AMREX_SPACEDIM,
552 0,AMREX_SPACEDIM> const& bct,
553 Dim3 const& dlo, Dim3 const& dhi) noexcept
554{
555 const Real dxi = dxinv[0];
556 const Real dyi = dxinv[1];
557 const auto lo = amrex::lbound(box);
558 const auto hi = amrex::ubound(box);
559
560 int k = 0;
561 for (int j = lo.y; j <= hi.y; ++j) {
562 for (int i = lo.x; i <= hi.x; ++i) {
563 if (apx(i,j,0) == Real(0.0))
564 {
565 fx(i,j,0,0) = Real(0.0);
566 fx(i,j,0,1) = Real(0.0);
567 fx(i,j,0,2) = Real(0.0);
568 fx(i,j,0,3) = Real(0.0);
569 }
570 else
571 {
572 Real dudx = (vel(i,j,0,0) - vel(i-1,j,0,0))*dxi;
573 Real dvdx = (vel(i,j,0,1) - vel(i-1,j,0,1))*dxi;
574
575 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
576 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
577 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
578 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
579 Real whi = mlebtensor_weight(jhip-jhim);
580 Real wlo = mlebtensor_weight(jlop-jlom);
581 Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
582 bvxlo,bvxhi,bct,dlo,dhi,
583 whi,wlo,jhip,jhim,jlop,jlom);
584 Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
585 bvxlo,bvxhi,bct,dlo,dhi,
586 whi,wlo,jhip,jhim,jlop,jlom);
587 fx(i,j,0,0) = dudx;
588 fx(i,j,0,1) = dvdx;
589 fx(i,j,0,2) = dudy;
590 fx(i,j,0,3) = dvdy;
591 }
592 }
593 }
594}
595
597void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
598 Array4<Real const> const& vel, Array4<Real const> const& apy,
599 Array4<EBCellFlag const> const& flag,
601 Array4<Real const> const& bvylo,
602 Array4<Real const> const& bvyhi,
604 0,2*AMREX_SPACEDIM,
605 0,AMREX_SPACEDIM> const& bct,
606 Dim3 const& dlo, Dim3 const& dhi) noexcept
607{
608 const Real dxi = dxinv[0];
609 const Real dyi = dxinv[1];
610 const auto lo = amrex::lbound(box);
611 const auto hi = amrex::ubound(box);
612
613 int k = 0;
614 for (int j = lo.y; j <= hi.y; ++j) {
615 for (int i = lo.x; i <= hi.x; ++i) {
616 if (apy(i,j,0) == 0.) {
617 fy(i,j,0,0) = Real(0.0);
618 fy(i,j,0,1) = Real(0.0);
619 fy(i,j,0,2) = Real(0.0);
620 fy(i,j,0,3) = Real(0.0);
621 }
622 else
623 {
624 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
625 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
626 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
627 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
628 Real whi = mlebtensor_weight(ihip-ihim);
629 Real wlo = mlebtensor_weight(ilop-ilom);
630 Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
631 bvylo,bvyhi,bct,dlo,dhi,
632 whi,wlo,ihip,ihim,ilop,ilom);
633 Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
634 bvylo,bvyhi,bct,dlo,dhi,
635 whi,wlo,ihip,ihim,ilop,ilom);
636
637 Real dudy = (vel(i,j,0,0) - vel(i,j-1,0,0))*dyi;
638 Real dvdy = (vel(i,j,0,1) - vel(i,j-1,0,1))*dyi;
639 fy(i,j,0,0) = dudx;
640 fy(i,j,0,1) = dvdx;
641 fy(i,j,0,2) = dudy;
642 fy(i,j,0,3) = dvdy;
643 }
644 }
645 }
646}
647
649void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
650 Array4<Real const> const& vel, Array4<Real const> const& apx,
651 Array4<EBCellFlag const> const& flag,
652 Array4<int const> const& ccm,
653 Array4<Real const> const& fcx,
654 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
655{
656 const Real dxi = dxinv[0];
657 const Real dyi = dxinv[1];
658 const auto lo = amrex::lbound(box);
659 const auto hi = amrex::ubound(box);
660
661 int k = 0;
662 for (int j = lo.y; j <= hi.y; ++j) {
664 for (int i = lo.x; i <= hi.x; ++i) {
665 if (apx(i,j,0) == Real(0.0))
666 {
667 fx(i,j,0,0) = Real(0.0);
668 fx(i,j,0,1) = Real(0.0);
669 fx(i,j,0,2) = Real(0.0);
670 fx(i,j,0,3) = Real(0.0);
671 }
672 else
673 {
674 Real dudx = (vel(i,j,0,0) - vel(i-1,j,0,0))*dxi;
675 Real dvdx = (vel(i,j,0,1) - vel(i-1,j,0,1))*dxi;
676 if (apx(i,j,0) < Real(1.0)) {
677 int jj = j + (fcx(i,j,0,0) >= Real(0.0) ? 1 : -1);
678 Real fracy = (ccm(i-1,jj,0) || ccm(i,jj,0)) ? std::abs(fcx(i,j,0,0)) : Real(0.0);
679 dudx = (Real(1.0)-fracy)*dudx
680 + fracy *(vel(i,jj,0,0) - vel(i-1,jj,0,0))*dxi;
681 dvdx = (Real(1.0)-fracy)*dvdx
682 + fracy *(vel(i,jj,0,1) - vel(i-1,jj,0,1))*dxi;
683 }
684
685 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
686 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
687 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
688 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
689 Real whi = mlebtensor_weight(jhip-jhim);
690 Real wlo = mlebtensor_weight(jlop-jlom);
691 Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
692 whi,wlo,jhip,jhim,jlop,jlom);
693 Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
694 whi,wlo,jhip,jhim,jlop,jlom);
695 fx(i,j,0,0) = dudx;
696 fx(i,j,0,1) = dvdx;
697 fx(i,j,0,2) = dudy;
698 fx(i,j,0,3) = dvdy;
699 }
700 }
701 }
702}
703
705void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
706 Array4<Real const> const& vel, Array4<Real const> const& apy,
707 Array4<EBCellFlag const> const& flag,
708 Array4<int const> const& ccm,
709 Array4<Real const> const& fcy,
710 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
711{
712 const Real dxi = dxinv[0];
713 const Real dyi = dxinv[1];
714 const auto lo = amrex::lbound(box);
715 const auto hi = amrex::ubound(box);
716
717 int k = 0;
718 for (int j = lo.y; j <= hi.y; ++j) {
720 for (int i = lo.x; i <= hi.x; ++i) {
721 if (apy(i,j,0) == 0.) {
722 fy(i,j,0,0) = Real(0.0);
723 fy(i,j,0,1) = Real(0.0);
724 fy(i,j,0,2) = Real(0.0);
725 fy(i,j,0,3) = Real(0.0);
726 }
727 else
728 {
729 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
730 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
731 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
732 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
733 Real whi = mlebtensor_weight(ihip-ihim);
734 Real wlo = mlebtensor_weight(ilop-ilom);
735 Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
736 whi,wlo,ihip,ihim,ilop,ilom);
737 Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
738 whi,wlo,ihip,ihim,ilop,ilom);
739
740 Real dudy = (vel(i,j,0,0) - vel(i,j-1,0,0))*dyi;
741 Real dvdy = (vel(i,j,0,1) - vel(i,j-1,0,1))*dyi;
742 if (apy(i,j,0) < Real(1.0)) {
743 int ii = i + (fcy(i,j,0,0) >= Real(0.0) ? 1 : -1);
744 Real fracx = (ccm(ii,j-1,0) || ccm(ii,j,0)) ? std::abs(fcy(i,j,0,0)) : Real(0.0);
745 dudy = (Real(1.0)-fracx)*dudy
746 + fracx *(vel(ii,j,0,0) - vel(ii,j-1,0,0))*dyi;
747 dvdy = (Real(1.0)-fracx)*dvdy
748 + fracx *(vel(ii,j,0,1) - vel(ii,j-1,0,1))*dyi;
749 }
750
751 fy(i,j,0,0) = dudx;
752 fy(i,j,0,1) = dvdx;
753 fy(i,j,0,2) = dudy;
754 fy(i,j,0,3) = dvdy;
755 }
756 }
757 }
758}
759
761void mlebtensor_vel_grads_fx (Box const& box, Array4<Real> const& fx,
762 Array4<Real const> const& vel, Array4<Real const> const& apx,
763 Array4<EBCellFlag const> const& flag,
764 Array4<int const> const& ccm,
765 Array4<Real const> const& fcx,
767 Array4<Real const> const& bvxlo,
768 Array4<Real const> const& bvxhi,
770 0,2*AMREX_SPACEDIM,
771 0,AMREX_SPACEDIM> const& bct,
772 Dim3 const& dlo, Dim3 const& dhi) noexcept
773{
774 const Real dxi = dxinv[0];
775 const Real dyi = dxinv[1];
776 const auto lo = amrex::lbound(box);
777 const auto hi = amrex::ubound(box);
778
779 int k = 0;
780 for (int j = lo.y; j <= hi.y; ++j) {
781 for (int i = lo.x; i <= hi.x; ++i) {
782 if (apx(i,j,0) == Real(0.0))
783 {
784 fx(i,j,0,0) = Real(0.0);
785 fx(i,j,0,1) = Real(0.0);
786 fx(i,j,0,2) = Real(0.0);
787 fx(i,j,0,3) = Real(0.0);
788 }
789 else
790 {
791 Real dudx = (vel(i,j,0,0) - vel(i-1,j,0,0))*dxi;
792 Real dvdx = (vel(i,j,0,1) - vel(i-1,j,0,1))*dxi;
793 if (apx(i,j,0) < Real(1.0)) {
794 int jj = j + (fcx(i,j,0,0) >= Real(0.0) ? 1 : -1);
795 Real fracy = (ccm(i-1,jj,0) || ccm(i,jj,0)) ? std::abs(fcx(i,j,0,0)) : Real(0.0);
796 dudx = (Real(1.0)-fracy)*dudx
797 + fracy *(vel(i,jj,0,0) - vel(i-1,jj,0,0))*dxi;
798 dvdx = (Real(1.0)-fracy)*dvdx
799 + fracy *(vel(i,jj,0,1) - vel(i-1,jj,0,1))*dxi;
800 }
801
802 int jhip = j + flag(i ,j,0).isConnected(0, 1,0);
803 int jhim = j - flag(i ,j,0).isConnected(0,-1,0);
804 int jlop = j + flag(i-1,j,0).isConnected(0, 1,0);
805 int jlom = j - flag(i-1,j,0).isConnected(0,-1,0);
806 Real whi = mlebtensor_weight(jhip-jhim);
807 Real wlo = mlebtensor_weight(jlop-jlom);
808 Real dudy = mlebtensor_dy_on_xface(i,j,k,0,vel,dyi,
809 bvxlo,bvxhi,bct,dlo,dhi,
810 whi,wlo,jhip,jhim,jlop,jlom);
811 Real dvdy = mlebtensor_dy_on_xface(i,j,k,1,vel,dyi,
812 bvxlo,bvxhi,bct,dlo,dhi,
813 whi,wlo,jhip,jhim,jlop,jlom);
814 fx(i,j,0,0) = dudx;
815 fx(i,j,0,1) = dvdx;
816 fx(i,j,0,2) = dudy;
817 fx(i,j,0,3) = dvdy;
818 }
819 }
820 }
821}
822
824void mlebtensor_vel_grads_fy (Box const& box, Array4<Real> const& fy,
825 Array4<Real const> const& vel, Array4<Real const> const& apy,
826 Array4<EBCellFlag const> const& flag,
827 Array4<int const> const& ccm,
828 Array4<Real const> const& fcy,
830 Array4<Real const> const& bvylo,
831 Array4<Real const> const& bvyhi,
833 0,2*AMREX_SPACEDIM,
834 0,AMREX_SPACEDIM> const& bct,
835 Dim3 const& dlo, Dim3 const& dhi) noexcept
836{
837 const Real dxi = dxinv[0];
838 const Real dyi = dxinv[1];
839 const auto lo = amrex::lbound(box);
840 const auto hi = amrex::ubound(box);
841
842 int k = 0;
843 for (int j = lo.y; j <= hi.y; ++j) {
844 for (int i = lo.x; i <= hi.x; ++i) {
845 if (apy(i,j,0) == 0.) {
846 fy(i,j,0,0) = Real(0.0);
847 fy(i,j,0,1) = Real(0.0);
848 fy(i,j,0,2) = Real(0.0);
849 fy(i,j,0,3) = Real(0.0);
850 }
851 else
852 {
853 int ihip = i + flag(i,j ,0).isConnected( 1,0,0);
854 int ihim = i - flag(i,j ,0).isConnected(-1,0,0);
855 int ilop = i + flag(i,j-1,0).isConnected( 1,0,0);
856 int ilom = i - flag(i,j-1,0).isConnected(-1,0,0);
857 Real whi = mlebtensor_weight(ihip-ihim);
858 Real wlo = mlebtensor_weight(ilop-ilom);
859 Real dudx = mlebtensor_dx_on_yface(i,j,k,0,vel,dxi,
860 bvylo,bvyhi,bct,dlo,dhi,
861 whi,wlo,ihip,ihim,ilop,ilom);
862 Real dvdx = mlebtensor_dx_on_yface(i,j,k,1,vel,dxi,
863 bvylo,bvyhi,bct,dlo,dhi,
864 whi,wlo,ihip,ihim,ilop,ilom);
865
866 Real dudy = (vel(i,j,0,0) - vel(i,j-1,0,0))*dyi;
867 Real dvdy = (vel(i,j,0,1) - vel(i,j-1,0,1))*dyi;
868 if (apy(i,j,0) < Real(1.0)) {
869 int ii = i + (fcy(i,j,0,0) >= Real(0.0) ? 1 : -1);
870 Real fracx = (ccm(ii,j-1,0) || ccm(ii,j,0)) ? std::abs(fcy(i,j,0,0)) : Real(0.0);
871 dudy = (Real(1.0)-fracx)*dudy
872 + fracx *(vel(ii,j,0,0) - vel(ii,j-1,0,0))*dyi;
873 dvdy = (Real(1.0)-fracx)*dvdy
874 + fracx *(vel(ii,j,0,1) - vel(ii,j-1,0,1))*dyi;
875 }
876
877 fy(i,j,0,0) = dudx;
878 fy(i,j,0,1) = dvdx;
879 fy(i,j,0,2) = dudy;
880 fy(i,j,0,3) = dvdy;
881 }
882 }
883 }
884}
885
886}
887#endif
#define AMREX_PRAGMA_SIMD
Definition AMReX_Extension.H:80
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
Maintain an identifier for boundary condition types.
Definition AMReX_BoundCond.H:20
Definition AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlebtensor_weight(int d)
Definition AMReX_MLEBTensor_K.H:10
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_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_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_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_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_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
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_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_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_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_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void LoopConcurrent(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:150
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:35
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
Definition AMReX_Array.H:282
Definition AMReX_Array4.H:61
Definition AMReX_Dim3.H:12
Definition AMReX_Array.H:34