1 #ifndef AMREX_MLABECLAP_3D_K_H_
2 #define AMREX_MLABECLAP_3D_K_H_
3 #include <AMReX_Config.H>
16 T alpha, T beta) noexcept
18 const T dhx = beta*dxinv[0]*dxinv[0];
19 const T dhy = beta*dxinv[1]*dxinv[1];
20 const T dhz = beta*dxinv[2]*dxinv[2];
21 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n)
22 - dhx * (bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i ,j,k,n))
23 - bX(i ,j,k,n)*(
x(i ,j,k,n) -
x(i-1,j,k,n)))
24 - dhy * (bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j ,k,n))
25 - bY(i,j ,k,n)*(
x(i,j ,k,n) -
x(i,j-1,k,n)))
26 - dhz * (bZ(i,j,k+1,n)*(
x(i,j,k+1,n) -
x(i,j,k ,n))
27 - bZ(i,j,k ,n)*(
x(i,j,k ,n) -
x(i,j,k-1,n)));
40 T alpha, T beta) noexcept
42 const T dhx = beta*dxinv[0]*dxinv[0];
43 const T dhy = beta*dxinv[1]*dxinv[1];
44 const T dhz = beta*dxinv[2]*dxinv[2];
45 if (osm(i,j,k) == 0) {
48 y(i,j,k,n) = alpha*a(i,j,k)*
x(i,j,k,n)
49 - dhx * (bX(i+1,j,k,n)*(
x(i+1,j,k,n) -
x(i ,j,k,n))
50 - bX(i ,j,k,n)*(
x(i ,j,k,n) -
x(i-1,j,k,n)))
51 - dhy * (bY(i,j+1,k,n)*(
x(i,j+1,k,n) -
x(i,j ,k,n))
52 - bY(i,j ,k,n)*(
x(i,j ,k,n) -
x(i,j-1,k,n)))
53 - dhz * (bZ(i,j,k+1,n)*(
x(i,j,k+1,n) -
x(i,j,k ,n))
54 - bZ(i,j,k ,n)*(
x(i,j,k ,n) -
x(i,j,k-1,n)));
66 T alpha, T beta) noexcept
68 const T dhx = beta*dxinv[0]*dxinv[0];
69 const T dhy = beta*dxinv[1]*dxinv[1];
70 const T dhz = beta*dxinv[2]*dxinv[2];
71 x(i,j,k,n) /= alpha*a(i,j,k)
72 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
73 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
74 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
80 Array4<T const>
const& bx, T fac,
int ncomp) noexcept
85 for (
int n = 0; n < ncomp; ++n) {
86 for (
int k = lo.z; k <= hi.z; ++k) {
87 for (
int j = lo.y; j <= hi.y; ++j) {
89 for (
int i = lo.x; i <= hi.x; ++i) {
90 fx(i,j,k,n) = -fac*bx(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
100 Array4<T const>
const& bx, T fac,
int xlen,
int ncomp) noexcept
105 for (
int n = 0; n < ncomp; ++n) {
106 for (
int k = lo.z; k <= hi.z; ++k) {
107 for (
int j = lo.y; j <= hi.y; ++j) {
109 fx(i,j,k,n) = -fac*bx(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
111 fx(i,j,k,n) = -fac*bx(i,j,k,n)*(sol(i,j,k,n)-sol(i-1,j,k,n));
117 template <
typename T>
120 Array4<T const>
const& by, T fac,
int ncomp) noexcept
125 for (
int n = 0; n < ncomp; ++n) {
126 for (
int k = lo.z; k <= hi.z; ++k) {
127 for (
int j = lo.y; j <= hi.y; ++j) {
129 for (
int i = lo.x; i <= hi.x; ++i) {
130 fy(i,j,k,n) = -fac*by(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
137 template <
typename T>
140 Array4<T const>
const& by, T fac,
int ylen,
int ncomp) noexcept
145 for (
int n = 0; n < ncomp; ++n) {
146 for (
int k = lo.z; k <= hi.z; ++k) {
149 for (
int i = lo.x; i <= hi.x; ++i) {
150 fy(i,j,k,n) = -fac*by(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
154 for (
int i = lo.x; i <= hi.x; ++i) {
155 fy(i,j,k,n) = -fac*by(i,j,k,n)*(sol(i,j,k,n)-sol(i,j-1,k,n));
161 template <
typename T>
169 for (
int n = 0; n < ncomp; ++n) {
170 for (
int k = lo.z; k <= hi.z; ++k) {
171 for (
int j = lo.y; j <= hi.y; ++j) {
173 for (
int i = lo.x; i <= hi.x; ++i) {
174 fz(i,j,k,n) = -fac*bz(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
181 template <
typename T>
189 for (
int n = 0; n < ncomp; ++n) {
191 for (
int j = lo.y; j <= hi.y; ++j) {
193 for (
int i = lo.x; i <= hi.x; ++i) {
194 fz(i,j,k,n) = -fac*bz(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
199 for (
int j = lo.y; j <= hi.y; ++j) {
201 for (
int i = lo.x; i <= hi.x; ++i) {
202 fz(i,j,k,n) = -fac*bz(i,j,k,n)*(sol(i,j,k,n)-sol(i,j,k-1,n));
208 template <
typename T>
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;
266 template <
typename T>
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;
330 template <
typename T>
346 Box const& vbox) noexcept
351 T cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
352 ? f0(vlo.x,j,k,n) : T(0.0);
353 T cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
354 ? f1(i,vlo.y,k,n) : T(0.0);
355 T cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
356 ? f2(i,j,vlo.z,n) : T(0.0);
357 T cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
358 ? f3(vhi.x,j,k,n) : T(0.0);
359 T cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
360 ? f4(i,vhi.y,k,n) : T(0.0);
361 T cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
362 ? f5(i,j,vhi.z,n) : T(0.0);
364 T gamma = alpha*a(i,j,k)
365 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
366 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
367 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
370 - (dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
371 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
372 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5));
374 phi(i,j,k,n) += T(2.0/3.0) * (rhs(i,j,k,n) - Ax(i,j,k,n)) / g_m_d;
377 template <
typename T>
395 Box const& vbox) noexcept
397 if (osm(i,j,k) == 0) {
398 phi(i,j,k,n) = T(0.0);
403 T cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
404 ? f0(vlo.x,j,k,n) : T(0.0);
405 T cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
406 ? f1(i,vlo.y,k,n) : T(0.0);
407 T cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
408 ? f2(i,j,vlo.z,n) : T(0.0);
409 T cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
410 ? f3(vhi.x,j,k,n) : T(0.0);
411 T cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
412 ? f4(i,vhi.y,k,n) : T(0.0);
413 T cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
414 ? f5(i,j,vhi.z,n) : T(0.0);
416 T gamma = alpha*a(i,j,k)
417 + dhx*(bX(i,j,k,n)+bX(i+1,j,k,n))
418 + dhy*(bY(i,j,k,n)+bY(i,j+1,k,n))
419 + dhz*(bZ(i,j,k,n)+bZ(i,j,k+1,n));
422 - (dhx*(bX(i,j,k,n)*cf0 + bX(i+1,j,k,n)*cf3)
423 + dhy*(bY(i,j,k,n)*cf1 + bY(i,j+1,k,n)*cf4)
424 + dhz*(bZ(i,j,k,n)*cf2 + bZ(i,j,k+1,n)*cf5));
426 phi(i,j,k,n) += T(2.0/3.0) * (rhs(i,j,k,n) - Ax(i,j,k,n)) / g_m_d;
430 template <
typename T>
437 u_ls(0) = r_ls(0) / bet;
439 for (
int i = 1; i <= ilen - 1; i++) {
440 gam(i) = c_ls(i-1) / bet;
441 bet = b_ls(i) - a_ls(i)*gam(i);
443 u_ls(i) = (r_ls(i)-a_ls(i)*u_ls(i-1)) / bet;
445 for (
int i = ilen-2; i >= 0; i--) {
446 u_ls(i) = u_ls(i) - gam(i+1)*u_ls(i+1);
450 template <
typename T>
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);
739 template <
typename T>
742 int ncomp, T osfac) noexcept
746 for (
int n = 0; n < ncomp; ++n) {
747 for (
int k = lo.z; k <= hi.z; ++k) {
748 for (
int j = lo.y; j <= hi.y; ++j) {
749 for (
int i = lo.x; i <= hi.x; ++i) {
750 if ((osm(i-1,j,k)+osm(i,j,k)) == 1) {
751 bX(i,j,k,n) *= osfac;
757 template <
typename T>
760 int ncomp, T osfac) noexcept
764 for (
int n = 0; n < ncomp; ++n) {
765 for (
int k = lo.z; k <= hi.z; ++k) {
766 for (
int j = lo.y; j <= hi.y; ++j) {
767 for (
int i = lo.x; i <= hi.x; ++i) {
768 if ((osm(i,j-1,k)+osm(i,j,k)) == 1) {
769 bY(i,j,k,n) *= osfac;
775 template <
typename T>
778 int ncomp, T osfac) noexcept
782 for (
int n = 0; n < ncomp; ++n) {
783 for (
int k = lo.z; k <= hi.z; ++k) {
784 for (
int j = lo.y; j <= hi.y; ++j) {
785 for (
int i = lo.x; i <= hi.x; ++i) {
786 if ((osm(i,j,k-1)+osm(i,j,k)) == 1) {
787 bZ(i,j,k,n) *= osfac;
#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
int idir
Definition: AMReX_HypreMLABecLap.cpp:1093
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void overset_rescale_bcoef_x(Box const &box, Array4< T > const &bX, Array4< int const > const &osm, int ncomp, T osfac) noexcept
Definition: AMReX_MLABecLap_1D_K.H:239
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_adotx_os(int i, int, int, int n, Array4< T > const &y, Array4< T const > const &x, Array4< T const > const &a, Array4< T const > const &bX, Array4< int const > const &osm, GpuArray< T, AMREX_SPACEDIM > const &dxinv, T alpha, T beta) noexcept
Definition: AMReX_MLABecLap_1D_K.H:24
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_x(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, Array4< T const > const &bx, T fac, int ncomp) noexcept
Definition: AMReX_MLABecLap_1D_K.H:56
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_z(Box const &box, Array4< T > const &fz, Array4< T const > const &sol, Array4< T const > const &bz, T fac, int ncomp) noexcept
Definition: AMReX_MLABecLap_3D_K.H:163
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_y(Box const &box, Array4< T > const &fy, Array4< T const > const &sol, Array4< T const > const &by, T fac, int ncomp) noexcept
Definition: AMReX_MLABecLap_2D_K.H:104
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_gsrb_os(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Array4< int const > const &osm, Box const &vbox, int redblack) noexcept
Definition: AMReX_MLABecLap_1D_K.H:121
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 void overset_rescale_bcoef_y(Box const &box, Array4< T > const &bY, Array4< int const > const &osm, int ncomp, T osfac) noexcept
Definition: AMReX_MLABecLap_2D_K.H:437
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_jacobi(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Box const &vbox) noexcept
Definition: AMReX_MLABecLap_1D_K.H:160
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_gsrb(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Box const &vbox, int redblack) noexcept
Definition: AMReX_MLABecLap_1D_K.H:87
AMREX_FORCE_INLINE void abec_gsrb_with_line_solve(Box const &, Array4< T > const &, Array4< T const > const &, T, Array4< T const > const &, T, Array4< T const > const &, Array4< int const > const &, Array4< int const > const &, Array4< T const > const &, Array4< T const > const &, Box const &, int, int) noexcept
Definition: AMReX_MLABecLap_1D_K.H:223
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_yface(Box const &box, Array4< T > const &fy, Array4< T const > const &sol, Array4< T const > const &by, T fac, int ylen, int ncomp) noexcept
Definition: AMReX_MLABecLap_2D_K.H:122
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_normalize(int i, int, int, int n, Array4< T > const &x, Array4< T const > const &a, Array4< T const > const &bX, GpuArray< T, AMREX_SPACEDIM > const &dxinv, T alpha, T beta) noexcept
Definition: AMReX_MLABecLap_1D_K.H:44
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition: AMReX.cpp:225
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_zface(Box const &box, Array4< T > const &fz, Array4< T const > const &sol, Array4< T const > const &bz, T fac, int zlen, int ncomp) noexcept
Definition: AMReX_MLABecLap_3D_K.H:183
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void abec_jacobi_os(int i, int, int, int n, Array4< T > const &phi, Array4< T const > const &rhs, Array4< T const > const &Ax, T alpha, Array4< T const > const &a, T dhx, Array4< T const > const &bX, Array4< int const > const &m0, Array4< int const > const &m1, Array4< T const > const &f0, Array4< T const > const &f1, Array4< int const > const &osm, Box const &vbox) noexcept
Definition: AMReX_MLABecLap_1D_K.H:189
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_flux_xface(Box const &box, Array4< T > const &fx, Array4< T const > const &sol, Array4< T const > const &bx, T fac, int xlen, int ncomp) noexcept
Definition: AMReX_MLABecLap_1D_K.H:72
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void overset_rescale_bcoef_z(Box const &box, Array4< T > const &bZ, Array4< int const > const &osm, int ncomp, T osfac) noexcept
Definition: AMReX_MLABecLap_3D_K.H:777
AMREX_FORCE_INLINE void tridiagonal_solve(Array1D< T, 0, 31 > &a_ls, Array1D< T, 0, 31 > &b_ls, Array1D< T, 0, 31 > &c_ls, Array1D< T, 0, 31 > &r_ls, Array1D< T, 0, 31 > &u_ls, Array1D< T, 0, 31 > &gam, int ilen) noexcept
Definition: AMReX_MLABecLap_3D_K.H:432
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlabeclap_adotx(int i, int, int, int n, Array4< T > const &y, Array4< T const > const &x, Array4< T const > const &a, Array4< T const > const &bX, GpuArray< T, AMREX_SPACEDIM > const &dxinv, T alpha, T beta) noexcept
Definition: AMReX_MLABecLap_1D_K.H:9
Definition: AMReX_Array.H:161
Definition: AMReX_Array4.H:61
Definition: AMReX_Array.H:34