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;
273 if (std::abs(anrmx) > std::abs(anrmy)) {
274 dg = dx_eb / std::abs(anrmx);
276 dg = dx_eb / std::abs(anrmy);
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);
374 Real dhx, Real dhy, Real dh,
383 bool is_dirichlet,
bool beta_on_centroid,
bool phi_on_centroid,
384 Box const& vbox,
int redblack,
int ncomp)
noexcept
389 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n)
noexcept
391 if ((i+j+k+redblack) % 2 == 0)
394 if (flag.isCovered())
396 phi(i,j,k,n) = Real(0.0);
400 Real cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
401 ? f0(vlo.x,j,k,n) : Real(0.0);
402 Real cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
403 ? f1(i,vlo.y,k,n) : Real(0.0);
404 Real cf2 = (i == vhi.x && m2(vhi.x+1,j,k) > 0)
405 ? f2(vhi.x,j,k,n) : Real(0.0);
406 Real cf3 = (j == vhi.y && m3(i,vhi.y+1,k) > 0)
407 ? f3(i,vhi.y,k,n) : Real(0.0);
409 if (flag.isRegular())
411 Real gamma = alpha*a(i,j,k)
412 + dhx * (bX(i+1,j,k,n) + bX(i,j,k,n))
413 + dhy * (bY(i,j+1,k,n) + bY(i,j,k,n));
415 Real rho = dhx * (bX(i+1,j,k,n)*phi(i+1,j,k,n)
416 + bX(i ,j,k,n)*phi(i-1,j,k,n))
417 + dhy * (bY(i,j+1,k,n)*phi(i,j+1,k,n)
418 + bY(i,j ,k,n)*phi(i,j-1,k,n));
420 Real delta = dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf2)
421 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf3);
423 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
424 phi(i,j,k,n) += res/(gamma-delta);
429 Real apxm = ebdata.get<EBData_t::apx>(i ,j ,k);
430 Real apxp = ebdata.get<EBData_t::apx>(i+1,j ,k);
431 Real apym = ebdata.get<EBData_t::apy>(i ,j ,k);
432 Real apyp = ebdata.get<EBData_t::apy>(i ,j+1,k);
434 Real fxm = -bX(i,j,k,n)*phi(i-1,j,k,n);
435 Real oxm = -bX(i,j,k,n)*cf0;
436 Real sxm = bX(i,j,k,n);
437 if (apxm != Real(0.0) && apxm != Real(1.0)) {
438 Real fcx = ebdata.get<EBData_t::fcx>(i,j,k);
439 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx));
440 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k)) ? std::abs(fcx) : Real(0.0);
441 if (!beta_on_centroid && !phi_on_centroid)
443 fxm = (Real(1.0)-fracy)*fxm +
444 fracy *bX(i,jj,k,n)*(phi(i,jj,k,n)-phi(i-1,jj,k,n));
446 else if (beta_on_centroid && !phi_on_centroid)
448 fxm = (Real(1.0)-fracy)*( -phi(i-1,j,k,n)) +
449 fracy *(phi(i,jj,k,n)-phi(i-1,jj,k,n));
453 sxm = (Real(1.0)-fracy)*sxm;
456 Real fxp = bX(i+1,j,k,n)*phi(i+1,j,k,n);
457 Real oxp = bX(i+1,j,k,n)*cf2;
458 Real sxp = -bX(i+1,j,k,n);
459 if (apxp != Real(0.0) && apxp != Real(1.0)) {
460 Real fcx = ebdata.get<EBData_t::fcx>(i+1,j,k);
461 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx));
462 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k)) ? std::abs(fcx) : Real(0.0);
463 if (!beta_on_centroid && !phi_on_centroid)
465 fxp = (Real(1.0)-fracy)*fxp +
466 fracy *bX(i+1,jj,k,n)*(phi(i+1,jj,k,n)-phi(i,jj,k,n));
468 else if (beta_on_centroid && !phi_on_centroid)
470 fxp = (Real(1.0)-fracy)*(phi(i+1,j,k,n) ) +
471 fracy *(phi(i+1,jj,k,n)-phi(i,jj,k,n));
472 fxp *= bX(i+1,j,k,n);
475 sxp = (Real(1.0)-fracy)*sxp;
478 Real fym = -bY(i,j,k,n)*phi(i,j-1,k,n);
479 Real oym = -bY(i,j,k,n)*cf1;
480 Real sym = bY(i,j,k,n);
481 if (apym != Real(0.0) && apym != Real(1.0)) {
482 Real fcy = ebdata.get<EBData_t::fcy>(i,j,k);
483 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy));
484 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k)) ? std::abs(fcy) : Real(0.0);
485 if (!beta_on_centroid && !phi_on_centroid)
487 fym = (Real(1.0)-fracx)*fym +
488 fracx *bY(ii,j,k,n)*(phi(ii,j,k,n)-phi(ii,j-1,k,n));
490 else if (beta_on_centroid && !phi_on_centroid)
492 fym = (Real(1.0)-fracx)*( -phi( i,j-1,k,n)) +
493 fracx *(phi(ii,j,k,n)-phi(ii,j-1,k,n));
498 sym = (Real(1.0)-fracx)*sym;
501 Real fyp = bY(i,j+1,k,n)*phi(i,j+1,k,n);
502 Real oyp = bY(i,j+1,k,n)*cf3;
503 Real syp = -bY(i,j+1,k,n);
504 if (apyp != Real(0.0) && apyp != Real(1.0)) {
505 Real fcy = ebdata.get<EBData_t::fcy>(i,j+1,k);
506 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy));
507 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k)) ? std::abs(fcy) : Real(0.0);
508 if (!beta_on_centroid && !phi_on_centroid)
510 fyp = (Real(1.0)-fracx)*fyp +
511 fracx*bY(ii,j+1,k,n)*(phi(ii,j+1,k,n)-phi(ii,j,k,n));
513 else if (beta_on_centroid && !phi_on_centroid)
515 fyp = (Real(1.0)-fracx)*(phi(i,j+1,k,n) )+
516 fracx *(phi(ii,j+1,k,n)-phi(ii,j,k,n));
517 fyp *= bY(i,j+1,k,n);
520 syp = (Real(1.0)-fracx)*syp;
523 Real vfrcinv = (Real(1.0)/kappa);
524 Real gamma = alpha*a(i,j,k) + vfrcinv *
525 (dhx*(apxm*sxm-apxp*sxp) +
526 dhy*(apym*sym-apyp*syp));
527 Real rho = -vfrcinv *
528 (dhx*(apxm*fxm-apxp*fxp) +
529 dhy*(apym*fym-apyp*fyp));
531 Real delta = -vfrcinv *
532 (dhx*(apxm*oxm-apxp*oxp) +
533 dhy*(apym*oym-apyp*oyp));
536 Real dapx = (apxm-apxp)*dx[1];
537 Real dapy = (apym-apyp)*dx[0];
538 Real anorm = std::hypot(dapx,dapy);
539 Real anorminv = Real(1.0)/anorm;
540 Real anrmx = dapx * anorminv;
541 Real anrmy = dapy * anorminv;
545 Real dx_eb = get_dx_eb(kappa);
547 Real dg, gx, gy, sx, sy;
548 if (std::abs(anrmx) > std::abs(anrmy)) {
549 dg = dx_eb / std::abs(anrmx);
551 dg = dx_eb / std::abs(anrmy);
553 gx = bctx - dg*anrmx;
554 gy = bcty - dg*anrmy;
555 sx = std::copysign(Real(1.0),anrmx);
556 sy = std::copysign(Real(1.0),anrmy);
558 int ii = i -
static_cast<int>(sx);
559 int jj = j -
static_cast<int>(sy);
561 Real phig_gamma = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy);
562 Real phig = ( - gx*sx - gx*gy*sx*sy) * phi(ii,j ,k,n)
563 + ( - gy*sy - gx*gy*sx*sy) * phi(i ,jj,k,n)
564 + ( + gx*gy*sx*sy) * phi(ii,jj,k,n);
567 Real dphidn = ( -phig)/dg;
571 Real feb = dphidn * ba * beb(i,j,k,n);
572 rho += -vfrcinv*(-dh)*feb;
574 Real feb_gamma = -phig_gamma/dg * ba * beb(i,j,k,n);
575 gamma += vfrcinv*(-dh)*feb_gamma;
578 Real res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
579 phi(i,j,k,n) += res/(gamma-delta);
769 Real dhx, Real dhy, Real dh,
778 bool is_dirichlet,
bool beta_on_centroid,
int ncomp)
noexcept
780 amrex::Loop(box, ncomp, [=] (
int i,
int j,
int k,
int n)
noexcept
782 if (flag(i,j,k).isRegular())
784 phi(i,j,k,n) /= alpha*a(i,j,k) + dhx*(bX(i,j,k,n) + bX(i+1,j,k,n))
785 + dhy*(bY(i,j,k,n) + bY(i,j+1,k,n));
787 else if (flag(i,j,k).isSingleValued())
789 Real kappa = vfrc(i,j,k);
790 Real apxm = apx(i,j,k);
791 Real apxp = apx(i+1,j,k);
792 Real apym = apy(i,j,k);
793 Real apyp = apy(i,j+1,k);
795 Real sxm = bX(i,j,k,n);
796 if (apxm != Real(0.0) && apxm != Real(1.0) && !beta_on_centroid) {
797 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i,j,k)));
798 Real fracy = (ccm(i-1,jj,k) || ccm(i,jj,k))
799 ? std::abs(fcx(i,j,k)) : Real(0.0);
800 sxm = (Real(1.0)-fracy)*sxm;
803 Real sxp = -bX(i+1,j,k,n);
804 if (apxp != Real(0.0) && apxp != Real(1.0) && !beta_on_centroid) {
805 int jj = j +
static_cast<int>(std::copysign(Real(1.0),fcx(i+1,j,k)));
806 Real fracy = (ccm(i,jj,k) || ccm(i+1,jj,k))
807 ? std::abs(fcx(i+1,j,k)) : Real(0.0);
808 sxp = (Real(1.0)-fracy)*sxp;
811 Real sym = bY(i,j,k,n);
812 if (apym != Real(0.0) && apym != Real(1.0) && !beta_on_centroid) {
813 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j,k)));
814 Real fracx = (ccm(ii,j-1,k) || ccm(ii,j,k))
815 ? std::abs(fcy(i,j,k)) : Real(0.0);
816 sym = (Real(1.0)-fracx)*sym;
819 Real syp = -bY(i,j+1,k,n);
820 if (apyp != Real(0.0) && apyp != Real(1.0) && !beta_on_centroid) {
821 int ii = i +
static_cast<int>(std::copysign(Real(1.0),fcy(i,j+1,k)));
822 Real fracx = (ccm(ii,j,k) || ccm(ii,j+1,k))
823 ? std::abs(fcy(i,j+1,k)) : Real(0.0);
824 syp = (Real(1.0)-fracx)*syp;
827 Real vfrcinv = (Real(1.0)/kappa);
828 Real gamma = alpha*a(i,j,k) + vfrcinv *
829 (dhx*(apxm*sxm-apxp*sxp) +
830 dhy*(apym*sym-apyp*syp));
833 Real dapx = (apxm-apxp)*dx[1];
834 Real dapy = (apym-apyp)*dx[0];
835 Real anorm = std::hypot(dapx,dapy);
836 Real anorminv = Real(1.0)/anorm;
837 Real anrmx = dapx * anorminv;
838 Real anrmy = dapy * anorminv;
840 Real bctx = bc(i,j,k,0);
841 Real bcty = bc(i,j,k,1);
842 Real dx_eb = get_dx_eb(vfrc(i,j,k));
844 Real dg, gx, gy, sx, sy;
845 if (std::abs(anrmx) > std::abs(anrmy)) {
846 dg = dx_eb / std::abs(anrmx);
848 dg = dx_eb / std::abs(anrmy);
850 gx = bctx - dg*anrmx;
851 gy = bcty - dg*anrmy;
852 sx = std::copysign(Real(1.0),anrmx);
853 sy = std::copysign(Real(1.0),anrmy);
855 Real phig_gamma = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy);
856 Real feb_gamma = -phig_gamma/dg * ba(i,j,k) * beb(i,j,k,n);
857 gamma += vfrcinv*(-dh)*feb_gamma;
860 phi(i,j,k,n) /= gamma;