223 Box const& vbox,
int redblack)
noexcept
225 constexpr T omega = T(1.15);
227 if ((i+j+k+redblack)%2 == 0) {
231 T cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
232 ? f0(vlo.x,j,k,n) : T(0.0);
233 T cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
234 ? f1(i,vlo.y,k,n) : T(0.0);
235 T cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
236 ? f2(i,j,vlo.z,n) : T(0.0);
237 T cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
238 ? f3(vhi.x,j,k,n) : T(0.0);
239 T cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
240 ? f4(i,vhi.y,k,n) : T(0.0);
241 T cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
242 ? f5(i,j,vhi.z,n) : T(0.0);
244 T gamma = alpha*a(i,j,k)
245 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
246 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
247 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
250 - (dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
251 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
252 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5));
254 T rho = dhx*( bX(i ,j,k,n)*phi(i-1,j,k,n)
255 + bX(i+1,j,k,n)*phi(i+1,j,k,n) )
256 + dhy*( bY(i,j ,k,n)*phi(i,j-1,k,n)
257 + bY(i,j+1,k,n)*phi(i,j+1,k,n) )
258 + dhz*( bZ(i,j,k ,n)*phi(i,j,k-1,n)
259 + bZ(i,j,k+1,n)*phi(i,j,k+1,n) );
261 T res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
262 phi(i,j,k,n) = phi(i,j,k,n) + omega/g_m_d * res;
283 Box const& vbox,
int redblack)
noexcept
285 constexpr T omega = T(1.15);
287 if ((i+j+k+redblack)%2 == 0) {
288 if (osm(i,j,k) == 0) {
289 phi(i,j,k,n) = T(0.0);
294 T cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
295 ? f0(vlo.x,j,k,n) : T(0.0);
296 T cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
297 ? f1(i,vlo.y,k,n) : T(0.0);
298 T cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
299 ? f2(i,j,vlo.z,n) : T(0.0);
300 T cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
301 ? f3(vhi.x,j,k,n) : T(0.0);
302 T cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
303 ? f4(i,vhi.y,k,n) : T(0.0);
304 T cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
305 ? f5(i,j,vhi.z,n) : T(0.0);
307 T gamma = alpha*a(i,j,k)
308 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
309 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
310 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
313 - (dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
314 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
315 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5));
317 T rho = dhx*( bX(i ,j,k,n)*phi(i-1,j,k,n)
318 + bX(i+1,j,k,n)*phi(i+1,j,k,n) )
319 + dhy*( bY(i,j ,k,n)*phi(i,j-1,k,n)
320 + bY(i,j+1,k,n)*phi(i,j+1,k,n) )
321 + dhz*( bZ(i,j,k ,n)*phi(i,j,k-1,n)
322 + bZ(i,j,k+1,n)*phi(i,j,k+1,n) );
324 T res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
325 phi(i,j,k,n) = phi(i,j,k,n) + omega/g_m_d * res;
466 Box const& vbox,
int redblack,
int nc)
noexcept
477 int ilen = hi.z - lo.z + 1;
479 if ( (dhx <= dhy) && (dhz <= dhy) ) {
481 ilen = hi.y - lo.y + 1;
483 if ( (dhy <= dhx) && (dhz <= dhx) ) {
485 ilen = hi.x - lo.x + 1;
489 if (ilen > 32) {
amrex::Abort(
"abec_gsrb_with_line_solve is hard-wired to be no longer than 32"); }
499 for (
int n = 0; n < nc; ++n) {
500 for (
int j = lo.y; j <= hi.y; ++j) {
502 for (
int i = lo.x; i <= hi.x; ++i) {
503 if ((i+j+redblack)%2 == 0) {
505 for (
int k = lo.z; k <= hi.z; ++k)
507 T gamma = alpha*a(i,j,k)
508 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
509 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
510 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
512 T cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
513 ? f0(vlo.x,j,k,n) : T(0.0);
514 T cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
515 ? f1(i,vlo.y,k,n) : T(0.0);
516 T cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
517 ? f2(i,j,vlo.z,n) : T(0.0);
518 T cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
519 ? f3(vhi.x,j,k,n) : T(0.0);
520 T cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
521 ? f4(i,vhi.y,k,n) : T(0.0);
522 T cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
523 ? f5(i,j,vhi.z,n) : T(0.0);
526 - (dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
527 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
528 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5));
530 T rho = dhx*( bX(i ,j,k,n)*phi(i-1,j,k,n)
531 + bX(i+1,j,k,n)*phi(i+1,j,k,n) )
532 + dhy*( bY(i,j ,k,n)*phi(i,j-1,k,n)
533 + bY(i,j+1,k,n)*phi(i,j+1,k,n) );
536 if (i == vlo.x && m0(vlo.x-1,j,k) > 0) {
537 rho -= dhx*bX(i ,j,k,n)*phi(i-1,j,k,n);
539 if (i == vhi.x && m3(vhi.x+1,j,k) > 0) {
540 rho -= dhx*bX(i+1,j,k,n)*phi(i+1,j,k,n);
542 if (j == vlo.y && m1(i,vlo.y-1,k) > 0) {
543 rho -= dhy*bY(i,j ,k,n)*phi(i,j-1,k,n);
545 if (j == vhi.y && m4(i,vhi.y+1,k) > 0) {
546 rho -= dhy*bY(i,j+1,k,n)*phi(i,j+1,k,n);
549 a_ls(k-lo.z) = -dhz*bZ(i,j,k,n);
550 b_ls(k-lo.z) = g_m_d;
551 c_ls(k-lo.z) = -dhz*bZ(i,j,k+1,n);
552 u_ls(k-lo.z) = T(0.);
553 r_ls(k-lo.z) = rhs(i,j,k,n) + rho;
558 a_ls(k-lo.z) = T(0.);
559 if (!(m2(i,j,vlo.z-1) > 0)) { r_ls(k-lo.z) += dhz*bZ(i,j,k,n)*phi(i,j,k-1,n); }
563 c_ls(k-lo.z) = T(0.);
564 if (!(m5(i,j,vhi.z+1) > 0)) { r_ls(k-lo.z) += dhz*bZ(i,j,k+1,n)*phi(i,j,k+1,n); }
570 for (
int k = lo.z; k <= hi.z; ++k)
572 phi(i,j,k,n) = u_ls(k-lo.z);
578 }
else if (
idir == 1) {
579 for (
int n = 0; n < nc; ++n) {
580 for (
int i = lo.x; i <= hi.x; ++i) {
582 for (
int k = lo.z; k <= hi.z; ++k) {
583 if ((i+k+redblack)%2 == 0) {
585 for (
int j = lo.y; j <= hi.y; ++j)
587 T gamma = alpha*a(i,j,k)
588 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
589 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
590 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
592 T cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
593 ? f0(vlo.x,j,k,n) : T(0.0);
594 T cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
595 ? f1(i,vlo.y,k,n) : T(0.0);
596 T cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
597 ? f2(i,j,vlo.z,n) : T(0.0);
598 T cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
599 ? f3(vhi.x,j,k,n) : T(0.0);
600 T cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
601 ? f4(i,vhi.y,k,n) : T(0.0);
602 T cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
603 ? f5(i,j,vhi.z,n) : T(0.0);
606 - (dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
607 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
608 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5));
610 T rho = dhx*( bX(i ,j,k,n)*phi(i-1,j,k,n)
611 + bX(i+1,j,k,n)*phi(i+1,j,k,n) )
612 + dhz*( bZ(i,j ,k,n)*phi(i,j,k-1,n)
613 + bZ(i,j,k+1,n)*phi(i,j,k+1,n) );
616 if (i == vlo.x && m0(vlo.x-1,j,k) > 0) {
617 rho -= dhx*bX(i ,j,k,n)*phi(i-1,j,k,n);
619 if (i == vhi.x && m3(vhi.x+1,j,k) > 0) {
620 rho -= dhx*bX(i+1,j,k,n)*phi(i+1,j,k,n);
622 if (k == vlo.z && m2(i,j,vlo.z-1) > 0) {
623 rho -= dhz*bZ(i,j ,k,n)*phi(i,j,k-1,n);
625 if (k == vhi.z && m5(i,j,vhi.z+1) > 0) {
626 rho -= dhz*bZ(i,j,k+1,n)*phi(i,j,k+1,n);
629 a_ls(j-lo.y) = -dhy*bY(i,j,k,n);
630 b_ls(j-lo.y) = g_m_d;
631 c_ls(j-lo.y) = -dhy*bY(i,j+1,k,n);
632 u_ls(j-lo.y) = T(0.);
633 r_ls(j-lo.y) = rhs(i,j,k,n) + rho;
637 a_ls(j-lo.y) = T(0.);
638 if (!(m1(i,vlo.y-1,k) > 0)) { r_ls(j-lo.y) += dhy*bY(i,j,k,n)*phi(i,j-1,k,n); }
642 c_ls(j-lo.y) = T(0.);
643 if (!(m4(i,vhi.y+1,k) > 0)) { r_ls(j-lo.y) += dhy*bY(i,j+1,k,n)*phi(i,j+1,k,n); }
649 for (
int j = lo.y; j <= hi.y; ++j)
651 phi(i,j,k,n) = u_ls(j-lo.y);
657 }
else if (
idir == 0) {
658 for (
int n = 0; n < nc; ++n) {
659 for (
int j = lo.y; j <= hi.y; ++j) {
661 for (
int k = lo.z; k <= hi.z; ++k) {
662 if ((j+k+redblack)%2 == 0) {
664 for (
int i = lo.x; i <= hi.x; ++i)
666 T gamma = alpha*a(i,j,k)
667 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
668 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
669 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
671 T cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
672 ? f0(vlo.x,j,k,n) : T(0.0);
673 T cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
674 ? f1(i,vlo.y,k,n) : T(0.0);
675 T cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
676 ? f2(i,j,vlo.z,n) : T(0.0);
677 T cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
678 ? f3(vhi.x,j,k,n) : T(0.0);
679 T cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
680 ? f4(i,vhi.y,k,n) : T(0.0);
681 T cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
682 ? f5(i,j,vhi.z,n) : T(0.0);
685 - (dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
686 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
687 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5));
689 T rho = dhy*( bY(i,j ,k,n)*phi(i,j-1,k,n)
690 + bY(i,j+1,k,n)*phi(i,j+1,k,n) )
691 + dhz*( bZ(i,j ,k,n)*phi(i,j,k-1,n)
692 + bZ(i,j,k+1,n)*phi(i,j,k+1,n) );
695 if (j == vlo.y && m1(i,vlo.y-1,k) > 0) {
696 rho -= dhy*bY(i,j ,k,n)*phi(i,j-1,k,n);
698 if (j == vhi.y && m4(i,vhi.y+1,k) > 0) {
699 rho -= dhy*bY(i,j+1,k,n)*phi(i,j+1,k,n);
701 if (k == vlo.z && m2(i,j,vlo.z-1) > 0) {
702 rho -= dhz*bZ(i,j ,k,n)*phi(i,j,k-1,n);
704 if (k == vhi.z && m5(i,j,vhi.z+1) > 0) {
705 rho -= dhz*bZ(i,j,k+1,n)*phi(i,j,k+1,n);
708 a_ls(i-lo.x) = -dhx*bX(i,j,k,n);
709 b_ls(i-lo.x) = g_m_d;
710 c_ls(i-lo.x) = -dhx*bX(i+1,j,k,n);
711 u_ls(i-lo.x) = T(0.);
712 r_ls(i-lo.x) = rhs(i,j,k,n) + rho;
716 a_ls(i-lo.x) = T(0.);
717 if (!(m0(vlo.x-1,j,k) > 0)) { r_ls(i-lo.x) += dhx*bX(i,j,k,n)*phi(i-1,j,k,n); }
721 c_ls(i-lo.x) = T(0.);
722 if (!(m3(vhi.x+1,j,k) > 0)) { r_ls(i-lo.x) += dhx*bX(i+1,j,k,n)*phi(i+1,j,k,n); }
728 for (
int i = lo.x; i <= hi.x; ++i)
730 phi(i,j,k,n) = u_ls(i-lo.x);