1 #ifndef AMREX_MLEBABECLAP_2D_K_H_
2 #define AMREX_MLEBABECLAP_2D_K_H_
3 #include <AMReX_Config.H>
21 const int& domlo_x,
const int& domlo_y,
22 const int& domhi_x,
const int& domhi_y,
23 const bool& on_x_face,
const bool& on_y_face,
24 bool is_eb_dirichlet,
bool is_eb_inhomog,
26 Real alpha, Real beta,
int ncomp) noexcept
28 Real dhx = beta*dxinv[0]*dxinv[0];
29 Real dhy = beta*dxinv[1]*dxinv[1];
30 Real dh = beta*dxinv[0]*dxinv[1];
32 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n) noexcept
34 if (flag(i,j,k).isCovered())
36 y(i,j,k,n) = Real(0.0);
38 else if (flag(i,j,k).isRegular() &&
39 ((flag(i-1,j ,k).isRegular() && flag(i+1,j ,k).isRegular() &&
40 flag(i ,j-1,k).isRegular() && flag(i ,j+1,k).isRegular()) ))
42 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n)
43 - dhx * (bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i ,j,k,n))
44 - bX(i ,j,k,n)*(
x(i ,j,k,n) -
x(i-1,j,k,n)))
45 - dhy * (bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j ,k,n))
46 - bY(i,j ,k,n)*(
x(i,j ,k,n) -
x(i,j-1,k,n)));
52 Real kappa = vfrc(i,j,k);
53 Real apxm = apx(i,j,k);
54 Real apxp = apx(i+1,j,k);
55 Real apym = apy(i,j,k);
56 Real apyp = apy(i,j+1,k);
59 bool needs_bdry_stencil = (i <= domlo_x) || (i >= domhi_x) ||
60 (j <= domlo_y) || (j >= domhi_y);
69 Real fxm = bX(i,j,k,n) * (
x(i,j,k,n)-
x(i-1,j,k,n));
70 if ( (apxm != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0)) )
72 Real yloc_on_xface = fcx(i,j,k);
74 if(needs_bdry_stencil) {
76 fxm =
grad_x_of_phi_on_centroids_extdir(i,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
77 yloc_on_xface,is_eb_dirichlet,is_eb_inhomog,
78 on_x_face, domlo_x, domhi_x,
79 on_y_face, domlo_y, domhi_y);
83 yloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
88 Real fxp = bX(i+1,j,k,n)*(
x(i+1,j,k,n)-
x(i,j,k,n));
89 if ( (apxp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i+1,j,k) != Real(1.0) || vfrc(i-1,j,k) != Real(1.0)) ) {
90 Real yloc_on_xface = fcx(i+1,j,k,0);
91 if(needs_bdry_stencil) {
92 fxp =
grad_x_of_phi_on_centroids_extdir(i+1,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
93 yloc_on_xface,is_eb_dirichlet,is_eb_inhomog,
94 on_x_face, domlo_x, domhi_x,
95 on_y_face, domlo_y, domhi_y);
99 yloc_on_xface,is_eb_dirichlet,is_eb_inhomog);
101 fxp *= bX(i+1,j,k,n);
105 Real fym = bY(i,j,k,n)*(
x(i,j,k,n)-
x(i,j-1,k,n));
106 if ( (apym != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0)) ) {
107 Real xloc_on_yface = fcy(i,j,k,0);
109 if(needs_bdry_stencil) {
111 fym =
grad_y_of_phi_on_centroids_extdir(i,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
112 xloc_on_yface,is_eb_dirichlet,is_eb_inhomog,
113 on_x_face, domlo_x, domhi_x,
114 on_y_face, domlo_y, domhi_y);
118 xloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
123 Real fyp = bY(i,j+1,k,n)*(
x(i,j+1,k,n)-
x(i,j,k,n));
124 if ( (apyp != Real(0.0)) && (vfrc(i,j,k) != Real(1.0) || vfrc(i,j+1,k) != Real(1.0) || vfrc(i,j-1,k) != Real(1.0)) ) {
125 Real xloc_on_yface = fcy(i,j+1,k,0);
126 if(needs_bdry_stencil) {
127 fyp =
grad_y_of_phi_on_centroids_extdir(i,j+1,k,n,
x,phieb,flag,ccent,bcent,vfrc,
128 xloc_on_yface,is_eb_dirichlet,is_eb_inhomog,
129 on_x_face, domlo_x, domhi_x,
130 on_y_face, domlo_y, domhi_y);
134 xloc_on_yface,is_eb_dirichlet,is_eb_inhomog);
136 fyp *= bY(i,j+1,k,n);
139 Real feb = Real(0.0);
140 if (is_eb_dirichlet && flag(i,j,k).isSingleValued())
142 Real dapx = (apxm-apxp)/dxinv[1];
143 Real dapy = (apym-apyp)/dxinv[0];
144 Real anorm = std::hypot(dapx,dapy);
145 Real anorminv = Real(1.0)/anorm;
146 Real anrmx = dapx * anorminv;
147 Real anrmy = dapy * anorminv;
150 feb =
grad_eb_of_phi_on_centroids_extdir(i,j,k,n,
x,phieb,flag,ccent,bcent,vfrc,
151 anrmx,anrmy,is_eb_inhomog,
152 on_x_face, domlo_x, domhi_x,
153 on_y_face, domlo_y, domhi_y);
155 feb *= ba(i,j,k) * beb(i,j,k,n);
159 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n) + (Real(1.0)/kappa) *
160 (dhx*(apxm*fxm-apxp*fxp) + dhy*(apym*fym-apyp*fyp) - dh*feb);
176 Real alpha, Real beta,
int ncomp,
177 bool beta_on_centroid,
bool phi_on_centroid) noexcept
179 Real dhx = beta*dxinv[0]*dxinv[0];
180 Real dhy = beta*dxinv[1]*dxinv[1];
181 Real dh = beta*dxinv[0]*dxinv[1];
184 bool beta_on_center = !(beta_on_centroid);
185 bool phi_on_center = !( phi_on_centroid);
187 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n) noexcept
189 if (flag(i,j,k).isCovered())
191 y(i,j,k,n) = Real(0.0);
193 else if (flag(i,j,k).isRegular())
195 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n)
196 - dhx * (bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i ,j,k,n))
197 - bX(i ,j,k,n)*(
x(i ,j,k,n) -
x(i-1,j,k,n)))
198 - dhy * (bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j ,k,n))
199 - bY(i,j ,k,n)*(
x(i,j ,k,n) -
x(i,j-1,k,n)));
203 Real kappa = vfrc(i,j,k);
204 Real apxm = apx(i,j,k);
205 Real apxp = apx(i+1,j,k);
206 Real apym = apy(i,j,k);
207 Real apyp = apy(i,j+1,k);
209 Real fxm = bX(i,j,k,n) * (
x(i,j,k,n)-
x(i-1,j,k,n));
210 if (apxm != Real(0.0) && apxm != Real(1.0)) {
211 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i,j,k)));
212 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ?
std::abs(fcx(i,j,k)) : Real(0.0);
213 if (beta_on_center && phi_on_center) {
214 fxm = (Real(1.0)-fracy)*fxm + fracy*bX(i,jj,k,n)*(
x(i,jj,k,n)-
x(i-1,jj,k,n));
215 }
else if (beta_on_centroid && phi_on_center) {
216 fxm = bX(i,j,k,n) * ( (Real(1.0)-fracy)*(
x(i, j,k,n)-
x(i-1, j,k,n))
217 + fracy *(
x(i,jj,k,n)-
x(i-1,jj,k,n)) );
221 Real fxp = bX(i+1,j,k,n)*(
x(i+1,j,k,n)-
x(i,j,k,n));
222 if (apxp != Real(0.0) && apxp != Real(1.0)) {
223 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k)));
224 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ?
std::abs(fcx(i+1,j,k)) : Real(0.0);
225 if (beta_on_center && phi_on_center) {
226 fxp = (Real(1.0)-fracy)*fxp + fracy*bX(i+1,jj,k,n)*(
x(i+1,jj,k,n)-
x(i,jj,k,n));
227 }
else if (beta_on_centroid && phi_on_center) {
228 fxp = bX(i+1,j,k,n) * ( (Real(1.0)-fracy)*(
x(i+1, j,k,n)-
x(i, j,k,n))
229 + fracy *(
x(i+1,jj,k,n)-
x(i,jj,k,n)) );
233 Real fym = bY(i,j,k,n)*(
x(i,j,k,n)-
x(i,j-1,k,n));
234 if (apym != Real(0.0) && apym != Real(1.0)) {
235 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k)));
236 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ?
std::abs(fcy(i,j,k)) : Real(0.0);
237 if (beta_on_center && phi_on_center) {
238 fym = (Real(1.0)-fracx)*fym + fracx*bY(ii,j,k,n)*(
x(ii,j,k,n)-
x(ii,j-1,k,n));
239 }
else if (beta_on_centroid && phi_on_center) {
240 fym = bY(i,j,k,n) * ( (Real(1.0)-fracx)*(
x( i,j,k,n)-
x( i,j-1,k,n))
241 + fracx *(
x(ii,j,k,n)-
x(ii,j-1,k,n)) );
245 Real fyp = bY(i,j+1,k,n)*(
x(i,j+1,k,n)-
x(i,j,k,n));
246 if (apyp != Real(0.0) && apyp != Real(1.0)) {
247 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k)));
248 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ?
std::abs(fcy(i,j+1,k)) : Real(0.0);
249 if (beta_on_center && phi_on_center) {
250 fyp = (Real(1.0)-fracx)*fyp + fracx*bY(ii,j+1,k,n)*(
x(ii,j+1,k,n)-
x(ii,j,k,n));
251 }
else if (beta_on_centroid && phi_on_center) {
252 fyp = bY(i,j+1,k,n) * ( (Real(1.0)-fracx)*(
x( i,j+1,k,n)-
x( i,j,k,n))
253 + fracx *(
x(ii,j+1,k,n)-
x(ii,j,k,n)) );
257 Real feb = Real(0.0);
259 Real dapx = (apxm-apxp)/dxinv[1];
260 Real dapy = (apym-apyp)/dxinv[0];
261 Real anorm = std::hypot(dapx,dapy);
262 Real anorminv = Real(1.0)/anorm;
263 Real anrmx = dapx * anorminv;
264 Real anrmy = dapy * anorminv;
266 Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
268 Real bctx = bc(i,j,k,0);
269 Real bcty = bc(i,j,k,1);
270 Real dx_eb = get_dx_eb(kappa);
272 Real dg, gx, gy, sx, sy;
278 gx = (bctx - dg*anrmx);
279 gy = (bcty - dg*anrmy);
280 sx = std::copysign(Real(1.0),anrmx);
281 sy = std::copysign(Real(1.0),anrmy);
283 int ii = i -
static_cast<int>(sx);
284 int jj = j -
static_cast<int>(sy);
286 Real phig = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy) *
x(i ,j ,k,n)
287 + ( - gx*sx - gx*gy*sx*sy) *
x(ii,j ,k,n)
288 + ( - gy*sy - gx*gy*sx*sy) *
x(i ,jj,k,n)
289 + ( + gx*gy*sx*sy) *
x(ii,jj,k,n) ;
291 Real dphidn = (phib-phig) / dg;
293 feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
296 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n) + (Real(1.0)/kappa) *
297 (dhx*(apxm*fxm-apxp*fxp) +
298 dhy*(apym*fym-apyp*fyp) - dh*feb);
318 if (!flag(i,j,k).isSingleValued())
320 feb(i,j,k,n) = Real(0.0);
324 Real kappa = vfrc(i,j,k);
325 Real apxm = apx(i,j,k);
326 Real apxp = apx(i+1,j,k);
327 Real apym = apy(i,j,k);
328 Real apyp = apy(i,j+1,k);
330 Real dapx = (apxm-apxp)/dxinv[1];
331 Real dapy = (apym-apyp)/dxinv[0];
332 Real anorm = std::hypot(dapx,dapy);
333 Real anorminv = Real(1.0)/anorm;
334 Real anrmx = dapx * anorminv;
335 Real anrmy = dapy * anorminv;
336 const Real bareascaling =
std::sqrt( (anrmx/dxinv[0])*(anrmx/dxinv[0]) +
337 (anrmy/dxinv[1])*(anrmy/dxinv[1]) );
340 Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
342 Real bctx = bc(i,j,k,0);
343 Real bcty = bc(i,j,k,1);
344 Real dx_eb = get_dx_eb(kappa);
346 Real dg, gx, gy, sx, sy;
352 gx = bctx - dg*anrmx;
353 gy = bcty - dg*anrmy;
354 sx = std::copysign(Real(1.0),anrmx);
355 sy = std::copysign(Real(1.0),anrmy);
357 int ii = i -
static_cast<int>(sx);
358 int jj = j -
static_cast<int>(sy);
360 Real phig = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy) *
x(i ,j ,k,n)
361 + ( - gx*sx - gx*gy*sx*sy) *
x(ii,j ,k,n)
362 + ( - gy*sy - gx*gy*sx*sy) *
x(i ,jj,k,n)
363 + ( + gx*gy*sx*sy) *
x(ii,jj,k,n) ;
365 Real dphidn = (phib-phig)/(dg * bareascaling);
366 feb(i,j,k,n) = -beb(i,j,k,n) * dphidn;
374 Real dhx, Real dhy, Real dh,
387 bool is_dirichlet,
bool beta_on_centroid,
bool phi_on_centroid,
388 Box const& vbox,
int redblack,
int ncomp) noexcept
393 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n) noexcept
395 if ((i+j+k+redblack) % 2 == 0)
397 if (flag(i,j,k).isCovered())
399 phi(i,j,k,n) = Real(0.0);
403 Real cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
404 ? f0(vlo.x,j,k,n) : Real(0.0);
405 Real cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
406 ? f1(i,vlo.y,k,n) : Real(0.0);
407 Real cf2 = (i == vhi.x && m2(vhi.x+1,j,k) > 0)
408 ? f2(vhi.x,j,k,n) : Real(0.0);
409 Real cf3 = (j == vhi.y && m3(i,vhi.y+1,k) > 0)
410 ? f3(i,vhi.y,k,n) : Real(0.0);
412 if (flag(i,j,k).isRegular())
414 Real gamma = alpha*a(i,j,k)
415 + dhx * (bX(i+1,j,k,n) + bX(i,j,k,n))
416 + dhy * (bY(i,j+1,k,n) + bY(i,j,k,n));
418 Real rho = dhx * (bX(i+1,j,k,n)*phi(i+1,j,k,n)
419 + bX(i ,j,k,n)*phi(i-1,j,k,n))
420 + dhy * (bY(i,j+1,k,n)*phi(i,j+1,k,n)
421 + bY(i,j ,k,n)*phi(i,j-1,k,n));
423 Real delta = dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf2)
424 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf3);
426 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
427 phi(i,j,k,n) += res/(gamma-delta);
431 Real kappa = vfrc(i,j,k);
432 Real apxm = apx(i,j,k);
433 Real apxp = apx(i+1,j,k);
434 Real apym = apy(i,j,k);
435 Real apyp = apy(i,j+1,k);
437 Real fxm = -bX(i,j,k,n)*phi(i-1,j,k,n);
438 Real oxm = -bX(i,j,k,n)*cf0;
439 Real sxm = bX(i,j,k,n);
440 if (apxm != Real(0.0) && apxm != Real(1.0)) {
441 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i,j,k)));
442 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
444 if (!beta_on_centroid && !phi_on_centroid)
446 fxm = (Real(1.0)-fracy)*fxm +
447 fracy *bX(i,jj,k,n)*(phi(i,jj,k,n)-phi(i-1,jj,k,n));
449 else if (beta_on_centroid && !phi_on_centroid)
451 fxm = (Real(1.0)-fracy)*( -phi(i-1,j,k,n)) +
452 fracy *(phi(i,jj,k,n)-phi(i-1,jj,k,n));
456 sxm = (Real(1.0)-fracy)*sxm;
459 Real fxp = bX(i+1,j,k,n)*phi(i+1,j,k,n);
460 Real oxp = bX(i+1,j,k,n)*cf2;
461 Real sxp = -bX(i+1,j,k,n);
462 if (apxp != Real(0.0) && apxp != Real(1.0)) {
463 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k)));
464 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
465 ?
std::abs(fcx(i+1,j,k)) : Real(0.0);
466 if (!beta_on_centroid && !phi_on_centroid)
468 fxp = (Real(1.0)-fracy)*fxp +
469 fracy *bX(i+1,jj,k,n)*(phi(i+1,jj,k,n)-phi(i,jj,k,n));
471 else if (beta_on_centroid && !phi_on_centroid)
473 fxp = (Real(1.0)-fracy)*(phi(i+1,j,k,n) ) +
474 fracy *(phi(i+1,jj,k,n)-phi(i,jj,k,n));
475 fxp *= bX(i+1,j,k,n);
478 sxp = (Real(1.0)-fracy)*sxp;
481 Real fym = -bY(i,j,k,n)*phi(i,j-1,k,n);
482 Real oym = -bY(i,j,k,n)*cf1;
483 Real sym = bY(i,j,k,n);
484 if (apym != Real(0.0) && apym != Real(1.0)) {
485 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k)));
486 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
488 if (!beta_on_centroid && !phi_on_centroid)
490 fym = (Real(1.0)-fracx)*fym +
491 fracx *bY(ii,j,k,n)*(phi(ii,j,k,n)-phi(ii,j-1,k,n));
493 else if (beta_on_centroid && !phi_on_centroid)
495 fym = (Real(1.0)-fracx)*( -phi( i,j-1,k,n)) +
496 fracx *(phi(ii,j,k,n)-phi(ii,j-1,k,n));
501 sym = (Real(1.0)-fracx)*sym;
504 Real fyp = bY(i,j+1,k,n)*phi(i,j+1,k,n);
505 Real oyp = bY(i,j+1,k,n)*cf3;
506 Real syp = -bY(i,j+1,k,n);
507 if (apyp != Real(0.0) && apyp != Real(1.0)) {
508 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k)));
509 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
510 ?
std::abs(fcy(i,j+1,k)) : Real(0.0);
511 if (!beta_on_centroid && !phi_on_centroid)
513 fyp = (Real(1.0)-fracx)*fyp +
514 fracx*bY(ii,j+1,k,n)*(phi(ii,j+1,k,n)-phi(ii,j,k,n));
516 else if (beta_on_centroid && !phi_on_centroid)
518 fyp = (Real(1.0)-fracx)*(phi(i,j+1,k,n) )+
519 fracx *(phi(ii,j+1,k,n)-phi(ii,j,k,n));
520 fyp *= bY(i,j+1,k,n);
523 syp = (Real(1.0)-fracx)*syp;
526 Real vfrcinv = (Real(1.0)/kappa);
527 Real gamma = alpha*a(i,j,k) + vfrcinv *
528 (dhx*(apxm*sxm-apxp*sxp) +
529 dhy*(apym*sym-apyp*syp));
530 Real rho = -vfrcinv *
531 (dhx*(apxm*fxm-apxp*fxp) +
532 dhy*(apym*fym-apyp*fyp));
534 Real delta = -vfrcinv *
535 (dhx*(apxm*oxm-apxp*oxp) +
536 dhy*(apym*oym-apyp*oyp));
539 Real dapx = (apxm-apxp)*dx[1];
540 Real dapy = (apym-apyp)*dx[0];
541 Real anorm = std::hypot(dapx,dapy);
542 Real anorminv = Real(1.0)/anorm;
543 Real anrmx = dapx * anorminv;
544 Real anrmy = dapy * anorminv;
546 Real bctx = bc(i,j,k,0);
547 Real bcty = bc(i,j,k,1);
548 Real dx_eb = get_dx_eb(vfrc(i,j,k));
550 Real dg, gx, gy, sx, sy;
556 gx = bctx - dg*anrmx;
557 gy = bcty - dg*anrmy;
558 sx = std::copysign(Real(1.0),anrmx);
559 sy = std::copysign(Real(1.0),anrmy);
561 int ii = i -
static_cast<int>(sx);
562 int jj = j -
static_cast<int>(sy);
564 Real phig_gamma = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy);
565 Real phig = ( - gx*sx - gx*gy*sx*sy) * phi(ii,j ,k,n)
566 + ( - gy*sy - gx*gy*sx*sy) * phi(i ,jj,k,n)
567 + ( + gx*gy*sx*sy) * phi(ii,jj,k,n);
570 Real dphidn = ( -phig)/dg;
572 Real feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
573 rho += -vfrcinv*(-dh)*feb;
575 Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
576 gamma += vfrcinv*(-dh)*feb_gamma;
579 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
580 phi(i,j,k,n) += res/(gamma-delta);
591 Real dhx,
int face_only,
int ncomp,
Box const& xbox,
592 bool beta_on_centroid,
bool phi_on_centroid) noexcept
594 int lof = xbox.smallEnd(0);
595 int hif = xbox.bigEnd(0);
598 if (!face_only || lof == i || hif == i) {
599 if (apx(i,j,k) == Real(0.0)) {
600 fx(i,j,k,n) = Real(0.0);
601 }
else if (apx(i,j,k) == Real(1.0)) {
602 fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
604 Real fxm = bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
605 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i,j,k)));
606 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ?
std::abs(fcx(i,j,k)) : Real(0.0);
607 if (!beta_on_centroid && !phi_on_centroid)
609 fxm = (Real(1.0)-fracy)*fxm + fracy*bX(i,jj,k,n)*(sol(i,jj,k,n)-sol(i-1,jj,k,n));
610 }
else if (beta_on_centroid && !phi_on_centroid) {
611 fxm = bX(i,j,k,n) * ( (Real(1.0)-fracy)*(sol(i, j,k,n)-sol(i-1, j,k,n)) +
612 fracy *(sol(i,jj,k,n)-sol(i-1,jj,k,n)) );
614 fx(i,j,k,n) = -fxm*dhx;
624 Real dhy,
int face_only,
int ncomp,
Box const& ybox,
625 bool beta_on_centroid,
bool phi_on_centroid) noexcept
627 int lof = ybox.smallEnd(1);
628 int hif = ybox.bigEnd(1);
631 if (!face_only || lof == j || hif == j) {
632 if (apy(i,j,k) == Real(0.0)) {
633 fy(i,j,k,n) = Real(0.0);
634 }
else if (apy(i,j,k) == Real(1.0)) {
635 fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
637 Real fym = bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
638 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k)));
639 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ?
std::abs(fcy(i,j,k)) : Real(0.0);
640 if (!beta_on_centroid && !phi_on_centroid)
642 fym = (Real(1.0)-fracx)*fym + fracx*bY(ii,j,k,n)*(sol(ii,j,k,n)-sol(ii,j-1,k,n));
643 }
else if (beta_on_centroid && !phi_on_centroid) {
644 fym = bY(i,j,k,n) * ( (Real(1.0)-fracx)*(sol( i,j,k,n)-sol( i,j-1,k,n)) +
645 fracx *(sol(ii,j,k,n)-sol(ii,j-1,k,n)) );
647 fy(i,j,k,n) = -fym*dhy;
656 Real dhx,
int face_only,
int ncomp,
Box const& xbox) noexcept
658 int lof = xbox.smallEnd(0);
659 int hif = xbox.bigEnd(0);
662 if (!face_only || lof == i || hif == i) {
663 if (apx(i,j,k) == Real(0.0)) {
664 fx(i,j,k,n) = Real(0.0);
666 fx(i,j,k,n) = -dhx*bX(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
675 Real dhy,
int face_only,
int ncomp,
Box const& ybox) noexcept
677 int lof = ybox.smallEnd(1);
678 int hif = ybox.bigEnd(1);
681 if (!face_only || lof == j || hif == j) {
682 if (apy(i,j,k) == Real(0.0)) {
683 fy(i,j,k,n) = Real(0.0);
685 fy(i,j,k,n) = -dhy*bY(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
695 Real dxi,
int ncomp,
bool phi_on_centroid) noexcept
699 if (apx(i,j,k) == Real(0.0)) {
700 gx(i,j,k,n) = Real(0.0);
701 }
else if (apx(i,j,k) == Real(1.0)) {
702 gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
704 Real gxm = (sol(i,j,k,n)-sol(i-1,j,k,n));
705 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i,j,k)));
706 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ?
std::abs(fcx(i,j,k)) : Real(0.0);
707 if (!phi_on_centroid) {
708 gxm = (Real(1.0)-fracy)*gxm + fracy*(sol(i,jj,k,n)-sol(i-1,jj,k,n));
710 gx(i,j,k,n) = gxm*dxi;
719 Real dyi,
int ncomp,
bool phi_on_centroid) noexcept
723 if (apy(i,j,k) == Real(0.0)) {
724 gy(i,j,k,n) = Real(0.0);
725 }
else if (apy(i,j,k) == Real(1.0)) {
726 gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
728 Real gym = (sol(i,j,k,n)-sol(i,j-1,k,n));
729 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k)));
730 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ?
std::abs(fcy(i,j,k)) : Real(0.0);
731 if (!phi_on_centroid) {
732 gym = (Real(1.0)-fracx)*gym + fracx*(sol(ii,j,k,n)-sol(ii,j-1,k,n));
734 gy(i,j,k,n) = gym*dyi;
745 if (apx(i,j,k) == Real(0.0)) {
746 gx(i,j,k,n) = Real(0.0);
748 gx(i,j,k,n) = dxi*(sol(i,j,k,n)-sol(i-1,j,k,n));
759 if (apy(i,j,k) == Real(0.0)) {
760 gy(i,j,k,n) = Real(0.0);
762 gy(i,j,k,n) = dyi*(sol(i,j,k,n)-sol(i,j-1,k,n));
770 Real dhx, Real dhy, Real dh,
779 bool is_dirichlet,
bool beta_on_centroid,
int ncomp) noexcept
781 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n) noexcept
783 if (flag(i,j,k).isRegular())
785 phi(i,j,k,n) /= alpha*a(i,j,k) + dhx*(bX(i,j,k,n) + bX(i+1,j,k,n))
786 + dhy*(bY(i,j,k,n) + bY(i,j+1,k,n));
788 else if (flag(i,j,k).isSingleValued())
790 Real kappa = vfrc(i,j,k);
791 Real apxm = apx(i,j,k);
792 Real apxp = apx(i+1,j,k);
793 Real apym = apy(i,j,k);
794 Real apyp = apy(i,j+1,k);
796 Real sxm = bX(i,j,k,n);
797 if (apxm != Real(0.0) && apxm != Real(1.0) && !beta_on_centroid) {
798 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i,j,k)));
799 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
801 sxm = (Real(1.0)-fracy)*sxm;
804 Real sxp = -bX(i+1,j,k,n);
805 if (apxp != Real(0.0) && apxp != Real(1.0) && !beta_on_centroid) {
806 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k)));
807 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
808 ?
std::abs(fcx(i+1,j,k)) : Real(0.0);
809 sxp = (Real(1.0)-fracy)*sxp;
812 Real sym = bY(i,j,k,n);
813 if (apym != Real(0.0) && apym != Real(1.0) && !beta_on_centroid) {
814 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k)));
815 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
817 sym = (Real(1.0)-fracx)*sym;
820 Real syp = -bY(i,j+1,k,n);
821 if (apyp != Real(0.0) && apyp != Real(1.0) && !beta_on_centroid) {
822 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k)));
823 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
824 ?
std::abs(fcy(i,j+1,k)) : Real(0.0);
825 syp = (Real(1.0)-fracx)*syp;
828 Real vfrcinv = (Real(1.0)/kappa);
829 Real gamma = alpha*a(i,j,k) + vfrcinv *
830 (dhx*(apxm*sxm-apxp*sxp) +
831 dhy*(apym*sym-apyp*syp));
834 Real dapx = (apxm-apxp)*dx[1];
835 Real dapy = (apym-apyp)*dx[0];
836 Real anorm = std::hypot(dapx,dapy);
837 Real anorminv = Real(1.0)/anorm;
838 Real anrmx = dapx * anorminv;
839 Real anrmy = dapy * anorminv;
841 Real bctx = bc(i,j,k,0);
842 Real bcty = bc(i,j,k,1);
843 Real dx_eb = get_dx_eb(vfrc(i,j,k));
845 Real dg, gx, gy, sx, sy;
851 gx = bctx - dg*anrmx;
852 gy = bcty - dg*anrmy;
853 sx = std::copysign(Real(1.0),anrmx);
854 sy = std::copysign(Real(1.0),anrmy);
856 Real phig_gamma = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy);
857 Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
858 gamma += vfrcinv*(-dh)*feb_gamma;
861 phi(i,j,k,n) /= gamma;
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_x_0(Box const &box, Array4< Real > const &gx, Array4< Real const > const &sol, Array4< Real const > const &apx, Real dxi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:740
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:331
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_y_0(Box const &box, Array4< Real > const &gy, Array4< Real const > const &sol, Array4< Real const > const &apy, Real dyi, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:754
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_adotx_centroid(Box const &box, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &a, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ccent, Array4< Real const > const &ba, Array4< Real const > const &bcent, Array4< Real const > const &beb, Array4< Real const > const &phieb, const int &domlo_x, const int &domlo_y, const int &domhi_x, const int &domhi_y, const bool &on_x_face, const bool &on_y_face, bool is_eb_dirichlet, bool is_eb_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Real alpha, Real beta, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:11
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_adotx(Box const &box, Array4< Real > const &y, Array4< Real const > const &x, Array4< Real const > const &a, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< const int > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, Array4< Real const > const &phieb, bool is_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv, Real alpha, Real beta, int ncomp, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:166
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_x_0(Box const &box, Array4< Real > const &fx, Array4< Real const > const &apx, Array4< Real const > const &sol, Array4< Real const > const &bX, Real dhx, int face_only, int ncomp, Box const &xbox) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:654
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_x(Box const &box, Array4< Real > const &gx, Array4< Real const > const &sol, Array4< Real const > const &apx, Array4< Real const > const &fcx, Array4< int const > const &ccm, Real dxi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:692
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:308
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_x_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:495
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_y_0(Box const &box, Array4< Real > const &fy, Array4< Real const > const &apy, Array4< Real const > const &sol, Array4< Real const > const &bY, Real dhy, int face_only, int ncomp, Box const &ybox) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:673
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_eb_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &nrmx, Real &nrmy, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:765
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_x(Box const &box, Array4< Real > const &fx, Array4< Real const > const &apx, Array4< Real const > const &fcx, Array4< Real const > const &sol, Array4< Real const > const &bX, Array4< int const > const &ccm, Real dhx, int face_only, int ncomp, Box const &xbox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:588
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_gsrb(Box const &box, Array4< Real > const &phi, Array4< Real const > const &rhs, Real alpha, Array4< Real const > const &a, Real dhx, Real dhy, Real dh, GpuArray< Real, AMREX_SPACEDIM > const &dx, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< int const > const &m0, Array4< int const > const &m2, Array4< int const > const &m1, Array4< int const > const &m3, Array4< Real const > const &f0, Array4< Real const > const &f2, Array4< Real const > const &f1, Array4< Real const > const &f3, Array4< const int > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, bool beta_on_centroid, bool phi_on_centroid, Box const &vbox, int redblack, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:371
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:632
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void LoopConcurrent(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:150
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_ebflux(int i, int j, int k, int n, Array4< Real > const &feb, Array4< Real const > const &x, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &bc, Array4< Real const > const &beb, Array4< Real const > const &phieb, bool is_inhomog, GpuArray< Real, AMREX_SPACEDIM > const &dxinv) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:304
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_normalize(Box const &box, Array4< Real > const &phi, Real alpha, Array4< Real const > const &a, Real dhx, Real dhy, Real dh, const amrex::GpuArray< Real, AMREX_SPACEDIM > &dx, Array4< Real const > const &bX, Array4< Real const > const &bY, Array4< const int > const &ccm, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vfrc, Array4< Real const > const &apx, Array4< Real const > const &apy, Array4< Real const > const &fcx, Array4< Real const > const &fcy, Array4< Real const > const &ba, Array4< Real const > const &bc, Array4< Real const > const &beb, bool is_dirichlet, bool beta_on_centroid, int ncomp) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:768
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_grad_y(Box const &box, Array4< Real > const &gy, Array4< Real const > const &sol, Array4< Real const > const &apy, Array4< Real const > const &fcy, Array4< int const > const &ccm, Real dyi, int ncomp, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:716
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_x_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:246
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlebabeclap_flux_y(Box const &box, Array4< Real > const &fy, Array4< Real const > const &apy, Array4< Real const > const &fcy, Array4< Real const > const &sol, Array4< Real const > const &bY, Array4< int const > const &ccm, Real dhy, int face_only, int ncomp, Box const &ybox, bool beta_on_centroid, bool phi_on_centroid) noexcept
Definition: AMReX_MLEBABecLap_2D_K.H:621
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373
AMREX_GPU_HOST_DEVICE AMREX_ATTRIBUTE_FLATTEN_FOR void Loop(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:125
Definition: AMReX_Array4.H:61