diff --git a/Source/radiation/HABEC.H b/Source/radiation/HABEC.H index 0b3d718099..2c9df6b2e3 100644 --- a/Source/radiation/HABEC.H +++ b/Source/radiation/HABEC.H @@ -4,6 +4,8 @@ #include #include +#include + using namespace amrex; namespace HABEC @@ -35,37 +37,37 @@ namespace HABEC #if AMREX_SPACEDIM == 1 if (cdir == 0) { h = dx[0]; - x_left= true; + x_left = true; } else if (cdir == 1) { h = dx[0]; - x_right= true; + x_right = true; } #elif AMREX_SPACEDIM == 2 if (cdir == 0) { h = dx[0]; - x_left= true; + x_left = true; } else if (cdir == 2) { h = dx[0]; - x_right= true; + x_right = true; } else if (cdir == 1) { h = dx[1]; - y_left= true; + y_left = true; } else if (cdir == 3) { h = dx[1]; - y_right= true; + y_right = true; } #else if (cdir == 0) { h = dx[0]; - x_left= true; + x_left = true; } else if (cdir == 3) { h = dx[0]; - x_right= true; + x_right = true; } else if (cdir == 1) { h = dx[1]; @@ -198,6 +200,173 @@ namespace HABEC } } + AMREX_INLINE + void hbflx3 (Array4 const flux, + Array4 const er, + const Box& reg, + int cdir, int bctype, + Array4 const tf, + int bho, Real bcl, + Array4 const bcval, + Array4 const mask, + Array4 const b, + Real beta, const Real* const dx, Real c, + const Orientation& ori, + const GeometryData& geomdata, int inhom, + Array4 const spa) + { + Real h; + + // Index shift for whichever edge we're checking. + // Negative means we're looking at the lo edge, + // positive means we're looking at the hi edge. + // The first group is for cell centers, the second + // group is for cell edges. + + int icp = 0; + int jcp = 0; + int kcp = 0; + + int iep = 0; + int jep = 0; + int kep = 0; + +#if AMREX_SPACEDIM == 1 + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 1) { + h = dx[0]; + icp = 1; + iep = 1; + } +#elif AMREX_SPACEDIM == 2 + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 2) { + h = dx[0]; + icp = 1; + iep = 1; + } + else if (cdir == 1) { + h = dx[1]; + jcp = -1; + } + else if (cdir == 3) { + h = dx[1]; + jcp = 1; + jep = 1; + } +#else + if (cdir == 0) { + h = dx[0]; + icp = -1; + } + else if (cdir == 3) { + h = dx[0]; + icp = 1; + iep = 1; + } + else if (cdir == 1) { + h = dx[1]; + jcp = -1; + } + else if (cdir == 4) { + h = dx[1]; + jcp = 1; + jep = 1; + } + else if (cdir == 2) { + h = dx[2]; + kcp = -1; + } + else if (cdir == 5) { + h = dx[2]; + kcp = 1; + kep = 1; + } +#endif + + amrex::LoopOnCpu(reg, [=] (int i, int j, int k) noexcept + { + Real r; + face_metric(i, j, k, reg.loVect()[0], reg.hiVect()[0], geomdata, ori.coordDir(), ori.isLow(), r); + + Real bfv, bfm, bfm2, h2, th2; + int bct; + + if (mask.contains(i+icp,j+jcp,k+kcp)) { + if (mask(i+icp,j+jcp,k+kcp) > 0) { + if (bctype == -1) { + bct = tf(i+icp,j+jcp,k+kcp); + } + else { + bct = bctype; + } + if (bct == LO_DIRICHLET) { + if (bho >= 1) { + h2 = 0.5e0_rt * h; + th2 = 3.e0_rt * h2; + bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j,k); + bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j,k); + bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j,k); + } + else { + bfv = beta / (0.5e0_rt * h + bcl) * b(i,j,k); + bfm = bfv; + } + } + else if (bct == LO_NEUMANN) { + bfv = beta * r; + bfm = 0.e0_rt; + bfm2 = 0.e0_rt; + } + else if (bct == LO_MARSHAK) { + bfv = 2.e0_rt * beta * r; + if (bho >= 1) { + bfm = 0.375e0_rt * c * bfv; + bfm2 = -0.125e0_rt * c * bfv; + } + else { + bfm = 0.25e0_rt * c * bfv; + } + } + else if (bct == LO_SANCHEZ_POMRANING) { + bfv = 2.e0_rt * beta * r; + if (bho >= 1) { + bfm = 1.5e0_rt * spa(i,j,k) * c * bfv; + bfm2 = -0.5e0_rt * spa(i,j,k) * c * bfv; + } + else { + bfm = spa(i,j,k) * c * bfv; + } + } + else { + amrex::Error("hbflx3: unsupported boundary type"); + } + if (inhom == 0) { + bfv = 0.e0_rt; + } + + Real flux_sign = 1.0_rt; + if (iep != 0 || jep != 0 || kep != 0) { + // right edge + flux_sign = -1.0_rt; + } + + flux(i+iep,j+jep,k+kep) = flux_sign * (bfv * bcval(i+icp,j+jcp,k+kcp) - bfm * er(i,j,k)); + + if (bho >= 1) { + Real df = -bfm2 * er(i-icp,j-jcp,k-kcp); + flux(i+iep,j+jep,k+kep) -= flux_sign * df; + } + } + } + }); + } } // namespace HABEC #endif diff --git a/Source/radiation/HABEC_1D.F90 b/Source/radiation/HABEC_1D.F90 index 3f2aafc52b..544297195e 100644 --- a/Source/radiation/HABEC_1D.F90 +++ b/Source/radiation/HABEC_1D.F90 @@ -13,150 +13,6 @@ module habec_module contains -subroutine hbflx3(flux, & - DIMS(fbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bho, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - b, DIMS(bbox), & - beta, dx, c, r, inhom, & - spa, DIMS(spabox)) bind(C, name="hbflx3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(fbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(bbox) - integer :: DIMDEC(spabox) - integer :: cdir, bctype, tf(DIMV(bcv)), bho, inhom - real(rt) :: bcl, beta, dx(1), c - real(rt) :: flux(DIMV(fbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: b(DIMV(bbox)) - real(rt) :: spa(DIMV(spabox)) - real(rt) :: r(1) - real(rt) :: h, bfm, bfv, r0 - real(rt) :: bfm2, h2, th2 - integer :: i, bct - h = dx(1) - ! r is passed as an array but actually has only one element, which is - ! appropriate for the face we are doing here. - r0 = r(1) - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - if (mask(i-1) > 0) then - if (bctype == -1) then - bct = tf(i-1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta * r0 - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 0.375e0_rt * c * bfv - bfm2 = -0.125e0_rt * c * bfv - else - bfm = 0.25e0_rt * c * bfv - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 1.5e0_rt * spa(i) * c * bfv - bfm2 = -0.5e0_rt * spa(i) * c * bfv - else - bfm = spa(i) * c * bfv - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i) = (bfv * bcval(i-1) - bfm * er(i)) - if (bho >= 1) then - flux(i) = flux(i) - bfm2 * er(i+1) - endif - endif - else if (cdir == 1) then - ! Right face of grid - i = reg_h1 - if (mask(i+1) > 0) then - if (bctype == -1) then - bct = tf(i+1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i+1) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i+1) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i+1) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i+1) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta * r0 - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 0.375e0_rt * c * bfv - bfm2 = -0.125e0_rt * c * bfv - else - bfm = 0.25e0_rt * c * bfv - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 1.5e0_rt * spa(i) * c * bfv - bfm2 = -0.5e0_rt * spa(i) * c * bfv - else - bfm = spa(i) * c * bfv - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i+1) = -(bfv * bcval(i+1) - bfm * er(i)) - if (bho >= 1) then - flux(i+1) = flux(i+1) + bfm2 * er(i-1) - endif - endif - else - print *, "hbflx3: impossible face orientation" - endif -end subroutine hbflx3 - subroutine hdterm(dterm, & DIMS(dtbox), & er, DIMS(ebox), & diff --git a/Source/radiation/HABEC_2D.F90 b/Source/radiation/HABEC_2D.F90 index 2119c4150c..837e29d219 100644 --- a/Source/radiation/HABEC_2D.F90 +++ b/Source/radiation/HABEC_2D.F90 @@ -14,267 +14,6 @@ module habec_module contains -subroutine hbflx3(flux, & - DIMS(fbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bho, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - b, DIMS(bbox), & - beta, dx, c, r, inhom, & - spa, DIMS(spabox)) bind(C, name="hbflx3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(fbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(bbox) - integer :: DIMDEC(spabox) - integer :: cdir, bctype, tf(DIMV(bcv)), bho, inhom - real(rt) :: bcl, beta, dx(2), c - real(rt) :: flux(DIMV(fbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: b(DIMV(bbox)) - real(rt) :: spa(DIMV(spabox)) - real(rt) :: r(reg_l1:reg_h1) - real(rt) :: h, bfm, bfv, r0 - real(rt) :: bfm2, h2, th2 - integer :: i, j, bct - if (cdir == 0 .OR. cdir == 2) then - h = dx(1) - ! For the left and right faces, r is constant and the array actually has - ! only one element. This following "r(reg_l1)" looks wrong, but what it - ! really means is that the array dimensions are meaningless for this case. - r0 = r(reg_l1) - else - h = dx(2) - endif - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - do j = reg_l2, reg_h2 - if (mask(i-1,j) > 0) then - if (bctype == -1) then - bct = tf(i-1,j) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta * r0 - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 0.375e0_rt * c * bfv - bfm2 = -0.125e0_rt * c * bfv - else - bfm = 0.25e0_rt * c * bfv - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 1.5e0_rt * spa(i,j) * c * bfv - bfm2 = -0.5e0_rt * spa(i,j) * c * bfv - else - bfm = spa(i,j) * c * bfv - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j) = (bfv * bcval(i-1,j) - bfm * er(i,j)) - if (bho >= 1) then - flux(i,j) = flux(i,j) - bfm2 * er(i+1,j) - endif - endif - enddo - else if (cdir == 2) then - ! Right face of grid - i = reg_h1 - do j = reg_l2, reg_h2 - if (mask(i+1,j) > 0) then - if (bctype == -1) then - bct = tf(i+1,j) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i+1,j) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i+1,j) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i+1,j) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i+1,j) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta * r0 - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 0.375e0_rt * c * bfv - bfm2 = -0.125e0_rt * c * bfv - else - bfm = 0.25e0_rt * c * bfv - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta * r0 - if (bho >= 1) then - bfm = 1.5e0_rt * spa(i,j) * c * bfv - bfm2 = -0.5e0_rt * spa(i,j) * c * bfv - else - bfm = spa(i,j) * c * bfv - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i+1,j) = -(bfv * bcval(i+1,j) - bfm * er(i,j)) - if (bho >= 1) then - flux(i+1,j) = flux(i+1,j) + bfm2 * er(i-1,j) - endif - endif - enddo - else if (cdir == 1) then - ! Bottom face of grid - j = reg_l2 - do i = reg_l1, reg_h1 - if (mask(i,j-1) > 0) then - if (bctype == -1) then - bct = tf(i,j-1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta * r(i) - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta * r(i) - if (bho >= 1) then - bfm = 0.375e0_rt * c * bfv - bfm2 = -0.125e0_rt * c * bfv - else - bfm = 0.25e0_rt * c * bfv - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta * r(i) - if (bho >= 1) then - bfm = 1.5e0_rt * spa(i,j) * c * bfv - bfm2 = -0.5e0_rt * spa(i,j) * c * bfv - else - bfm = spa(i,j) * c * bfv - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j) = (bfv * bcval(i,j-1) - bfm * er(i,j)) - if (bho >= 1) then - flux(i,j) = flux(i,j) - bfm2 * er(i,j+1) - endif - endif - enddo - else if (cdir == 3) then - ! Top face of grid - j = reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j+1) > 0) then - if (bctype == -1) then - bct = tf(i,j+1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j+1) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j+1) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j+1) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j+1) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta * r(i) - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta * r(i) - if (bho >= 1) then - bfm = 0.375e0_rt * c * bfv - bfm2 = -0.125e0_rt * c * bfv - else - bfm = 0.25e0_rt * c * bfv - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta * r(i) - if (bho >= 1) then - bfm = 1.5e0_rt * spa(i,j) * c * bfv - bfm2 = -0.5e0_rt * spa(i,j) * c * bfv - else - bfm = spa(i,j) * c * bfv - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j+1) = -(bfv * bcval(i,j+1) - bfm * er(i,j)) - if (bho >= 1) then - flux(i,j+1) = flux(i,j+1) + bfm2 * er(i,j-1) - endif - endif - enddo - else - print *, "hbflx3: impossible face orientation" - endif -end subroutine hbflx3 - subroutine hdterm(dterm, & DIMS(dtbox), & er, DIMS(ebox), & diff --git a/Source/radiation/HABEC_3D.F90 b/Source/radiation/HABEC_3D.F90 index dcd02ebf88..d59edec61a 100644 --- a/Source/radiation/HABEC_3D.F90 +++ b/Source/radiation/HABEC_3D.F90 @@ -14,379 +14,6 @@ module habec_module contains -subroutine hbflx3(flux, & - DIMS(fbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bho, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - b, DIMS(bbox), & - beta, dx, c, r, inhom, & - spa, DIMS(spabox)) bind(C, name="hbflx3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(fbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(bbox) - integer :: DIMDEC(spabox) - integer :: cdir, bctype, tf(DIMV(bcv)), bho, inhom - real(rt) :: bcl, beta, dx(3), c - real(rt) :: flux(DIMV(fbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: b(DIMV(bbox)) - real(rt) :: spa(DIMV(spabox)) - real(rt) :: r(1) - real(rt) :: h, bfm, bfv - real(rt) :: bfm2, h2, th2 - integer :: i, j, k, bct - if (cdir == 0 .OR. cdir == 3) then - h = dx(1) - elseif (cdir == 1 .OR. cdir == 4) then - h = dx(2) - else - h = dx(3) - endif - if (cdir == 0) then - i = reg_l1 - do k = reg_l3, reg_h3 - do j = reg_l2, reg_h2 - if (mask(i-1,j,k) > 0) then - if (bctype == -1) then - bct = tf(i-1,j,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j,k) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j,k) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j,k) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j,k) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 0.75e0_rt * beta * c - bfm2 = -0.25e0_rt * beta * c - else - bfm = 0.5e0_rt * beta * c - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 3.0e0_rt * spa(i,j,k) * beta * c - bfm2 = -1.0e0_rt * spa(i,j,k) * beta * c - else - bfm = 2.0e0_rt * spa(i,j,k) * beta * c - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j,k) = (bfv * bcval(i-1,j,k) - bfm * er(i,j,k)) - if (bho >= 1) then - flux(i,j,k) = flux(i,j,k) - bfm2 * er(i+1,j,k) - endif - endif - enddo - enddo - else if (cdir == 3) then - i = reg_h1 - do k = reg_l3, reg_h3 - do j = reg_l2, reg_h2 - if (mask(i+1,j,k) > 0) then - if (bctype == -1) then - bct = tf(i+1,j,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i+1,j,k) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i+1,j,k) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i+1,j,k) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i+1,j,k) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 0.75e0_rt * beta * c - bfm2 = -0.25e0_rt * beta * c - else - bfm = 0.5e0_rt * beta * c - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 3.0e0_rt * spa(i,j,k) * beta * c - bfm2 = -1.0e0_rt * spa(i,j,k) * beta * c - else - bfm = 2.0e0_rt * spa(i,j,k) * beta * c - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i+1,j,k) = -(bfv * bcval(i+1,j,k) - bfm * er(i,j,k)) - if (bho >= 1) then - flux(i+1,j,k) = flux(i+1,j,k) + bfm2 * er(i-1,j,k) - endif - endif - enddo - enddo - else if (cdir == 1) then - j = reg_l2 - do k = reg_l3, reg_h3 - do i = reg_l1, reg_h1 - if (mask(i,j-1,k) > 0) then - if (bctype == -1) then - bct = tf(i,j-1,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j,k) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j,k) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j,k) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j,k) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 0.75e0_rt * beta * c - bfm2 = -0.25e0_rt * beta * c - else - bfm = 0.5e0_rt * beta * c - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 3.0e0_rt * spa(i,j,k) * beta * c - bfm2 = -1.0e0_rt * spa(i,j,k) * beta * c - else - bfm = 2.0e0_rt * spa(i,j,k) * beta * c - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j,k) = (bfv * bcval(i,j-1,k) - bfm * er(i,j,k)) - if (bho >= 1) then - flux(i,j,k) = flux(i,j,k) - bfm2 * er(i,j+1,k) - endif - endif - enddo - enddo - else if (cdir == 4) then - j = reg_h2 - do k = reg_l3, reg_h3 - do i = reg_l1, reg_h1 - if (mask(i,j+1,k) > 0) then - if (bctype == -1) then - bct = tf(i,j+1,k) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j+1,k) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j+1,k) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j+1,k) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j+1,k) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 0.75e0_rt * beta * c - bfm2 = -0.25e0_rt * beta * c - else - bfm = 0.5e0_rt * beta * c - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 3.0e0_rt * spa(i,j,k) * beta * c - bfm2 = -1.0e0_rt * spa(i,j,k) * beta * c - else - bfm = 2.0e0_rt * spa(i,j,k) * beta * c - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j+1,k) = -(bfv * bcval(i,j+1,k) - bfm * er(i,j,k)) - if (bho >= 1) then - flux(i,j+1,k) = flux(i,j+1,k) + bfm2 * er(i,j-1,k) - endif - endif - enddo - enddo - else if (cdir == 2) then - k = reg_l3 - do j = reg_l2, reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j,k-1) > 0) then - if (bctype == -1) then - bct = tf(i,j,k-1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j,k) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j,k) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j,k) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j,k) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 0.75e0_rt * beta * c - bfm2 = -0.25e0_rt * beta * c - else - bfm = 0.5e0_rt * beta * c - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 3.0e0_rt * spa(i,j,k) * beta * c - bfm2 = -1.0e0_rt * spa(i,j,k) * beta * c - else - bfm = 2.0e0_rt * spa(i,j,k) * beta * c - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j,k) = (bfv * bcval(i,j,k-1) - bfm * er(i,j,k)) - if (bho >= 1) then - flux(i,j,k) = flux(i,j,k) - bfm2 * er(i,j,k+1) - endif - endif - enddo - enddo - else if (cdir == 5) then - k = reg_h3 - do j = reg_l2, reg_h2 - do i = reg_l1, reg_h1 - if (mask(i,j,k+1) > 0) then - if (bctype == -1) then - bct = tf(i,j,k+1) - else - bct = bctype - endif - if (bct == LO_DIRICHLET) then - if (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) * b(i,j,k+1) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) * b(i,j,k+1) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) * b(i,j,k+1) - else - bfv = beta / (0.5e0_rt * h + bcl) * b(i,j,k+1) - bfm = bfv - endif - else if (bct == LO_NEUMANN) then - bfv = beta - bfm = 0.e0_rt - bfm2 = 0.e0_rt - else if (bct == LO_MARSHAK) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 0.75e0_rt * beta * c - bfm2 = -0.25e0_rt * beta * c - else - bfm = 0.5e0_rt * beta * c - endif - else if (bct == LO_SANCHEZ_POMRANING) then - bfv = 2.e0_rt * beta - if (bho >= 1) then - bfm = 3.0e0_rt * spa(i,j,k) * beta * c - bfm2 = -1.0e0_rt * spa(i,j,k) * beta * c - else - bfm = 2.0e0_rt * spa(i,j,k) * beta * c - endif - else - print *, "hbflx3: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - flux(i,j,k+1) = -(bfv * bcval(i,j,k+1) - bfm * er(i,j,k)) - if (bho >= 1) then - flux(i,j,k+1) = flux(i,j,k+1) + bfm2 * er(i,j,k-1) - endif - endif - enddo - enddo - else - print *, "hbflx3: impossible face orientation" - endif -end subroutine hbflx3 - subroutine hdterm(dterm, & DIMS(dtbox), & er, DIMS(ebox), & diff --git a/Source/radiation/HABEC_F.H b/Source/radiation/HABEC_F.H index 1ad5995275..72a3f7be97 100644 --- a/Source/radiation/HABEC_F.H +++ b/Source/radiation/HABEC_F.H @@ -10,19 +10,6 @@ extern "C" { #define RadBoundCond int #endif - void hbflx3(BL_FORT_FAB_ARG(flux), - BL_FORT_FAB_ARG(soln), - ARLIM_P(reglo), ARLIM_P(reghi), - const int& cdir, const int& bctype, const int* tf, - const int& bho, const amrex::Real& bcl, - const BL_FORT_FAB_ARG(bcval), - const BL_FORT_IFAB_ARG(mask), - BL_FORT_FAB_ARG(bcoefs), - const amrex::Real& beta, const amrex::Real* dx, - const amrex::Real& flux_factor, const amrex::Real* r, - const int& inhom, - const amrex::Real* spa, ARLIM_P(splo), ARLIM_P(sphi)); - void hdterm(BL_FORT_FAB_ARG(dterm), BL_FORT_FAB_ARG(soln), ARLIM_P(reglo), ARLIM_P(reghi), diff --git a/Source/radiation/HypreABec.cpp b/Source/radiation/HypreABec.cpp index f1a8213986..e3ac5a4931 100644 --- a/Source/radiation/HypreABec.cpp +++ b/Source/radiation/HypreABec.cpp @@ -255,12 +255,11 @@ void HypreABec::boundaryFlux(MultiFab* Flux, MultiFab& Soln, int icomp, const NGBndry& bd = getBndry(); const Box& domain = bd.getDomain(); - + #ifdef _OPENMP #pragma omp parallel #endif { - Vector r; Real foo=1.e200; for (MFIter si(Soln); si.isValid(); ++si) { @@ -275,37 +274,33 @@ void HypreABec::boundaryFlux(MultiFab* Flux, MultiFab& Soln, int icomp, const Mask &msk = bd.bndryMasks(oitr(),i); if (reg[oitr()] == domain[oitr()]) { - const int *tfp = NULL; int bctype = bct; + Array4 tf_arr; if (bd.mixedBndry(oitr())) { const BaseFab &tf = *(bd.bndryTypes(oitr())[i]); - tfp = tf.dataPtr(); + tf_arr = tf.array(); bctype = -1; } // In normal code operation only the fluxes at internal // Dirichlet boundaries are used. Some diagnostics use the // fluxes computed at domain boundaries but these do not // influence the evolution of the interior solution. - Real* pSPa; - Box SPabox; + Array4 sp_arr; if (SPa != 0) { - pSPa = (*SPa)[si].dataPtr(); - SPabox = (*SPa)[si].box(); - } - else { - pSPa = &foo; - SPabox = Box(IntVect::TheZeroVector(),IntVect::TheZeroVector()); + sp_arr = (*SPa)[si].array(); } - getFaceMetric(r, reg, oitr(), geom); - hbflx3(BL_TO_FORTRAN(Flux[idim][si]), - BL_TO_FORTRAN_N(Soln[si], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bctype, tfp, bho, bcl, - BL_TO_FORTRAN_N(fs, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*bcoefs[idim])[si]), - beta, dx, flux_factor, r.dataPtr(), inhom, - pSPa, ARLIM(SPabox.loVect()), ARLIM(SPabox.hiVect())); + HABEC::hbflx3(Flux[idim][si].array(), + Soln[si].array(icomp), + reg, + cdir, bctype, + tf_arr, + bho, bcl, + fs.array(bdcomp), + msk.array(), + (*bcoefs[idim])[si].array(), + beta, dx, flux_factor, + oitr(), geom.data(), inhom, + sp_arr); } else { HABEC::hbflx(Flux[idim][si].array(), diff --git a/Source/radiation/HypreMultiABec.cpp b/Source/radiation/HypreMultiABec.cpp index bc98e86bbd..b548767211 100644 --- a/Source/radiation/HypreMultiABec.cpp +++ b/Source/radiation/HypreMultiABec.cpp @@ -3909,7 +3909,6 @@ void HypreMultiABec::boundaryFlux(int level, #pragma omp parallel #endif { - Vector r; Real foo=1.e200; for (MFIter mfi(Soln); mfi.isValid(); ++mfi) { @@ -3927,38 +3926,35 @@ void HypreMultiABec::boundaryFlux(int level, const Box &msb = msk.box(); const Box &bbox = (*bcoefs[level])[idim][mfi].box(); if (reg[oitr()] == domain[oitr()]) { - const int *tfp = NULL; int bctype = bct; + Array4 tf_arr; if (bd[level]->mixedBndry(oitr())) { const BaseFab &tf = *(bd[level]->bndryTypes(oitr())[i]); - tfp = tf.dataPtr(); + tf_arr = tf.array(); bctype = -1; } // In normal code operation only the fluxes at internal // Dirichlet boundaries are used. Some diagnostics use the // fluxes computed at domain boundaries but these do not // influence the evolution of the interior solution. - Real* pSPa; - Box SPabox; + Array4 sp_arr; if (SPa[level]) { - pSPa = (*SPa[level])[mfi].dataPtr(); - SPabox = (*SPa[level])[mfi].box(); + sp_arr = (*SPa[level])[mfi].array(); } - else { - pSPa = &foo; - SPabox = Box(IntVect::TheZeroVector(),IntVect::TheZeroVector()); - } - getFaceMetric(r, reg, oitr(), geom[level]); - hbflx3(BL_TO_FORTRAN(Flux[idim][mfi]), - BL_TO_FORTRAN_N(Soln[mfi], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bctype, tfp, bho, bcl, - BL_TO_FORTRAN_N(fs, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*bcoefs[level])[idim][mfi]), - beta, geom[level].CellSize(), - flux_factor, r.dataPtr(), inhom, - pSPa, ARLIM(SPabox.loVect()), ARLIM(SPabox.hiVect())); + + HABEC::hbflx3(Flux[idim][mfi].array(), + Soln[mfi].array(icomp), + reg, + cdir, bctype, + tf_arr, + bho, bcl, + fs.array(bdcomp), + msk.array(), + (*bcoefs[level])[idim][mfi].array(), + beta, geom[level].CellSize(), + flux_factor, oitr(), + geom[level].data(), inhom, + sp_arr); } else { HABEC::hbflx(Flux[idim][mfi].array(),