From 9b81621ad4668811be60d992bcf0db3aee699c2e Mon Sep 17 00:00:00 2001 From: Max Katz Date: Tue, 10 Oct 2023 12:31:53 -0400 Subject: [PATCH] Convert hdterm and hdterm3 to C++ (#2625) --- Source/radiation/HABEC.H | 216 +++++++++++++++++ Source/radiation/HABEC_1D.F90 | 133 ----------- Source/radiation/HABEC_2D.F90 | 206 ----------------- Source/radiation/HABEC_3D.F90 | 308 ------------------------- Source/radiation/HABEC_F.H | 37 --- Source/radiation/HypreABec.cpp | 1 - Source/radiation/HypreExtMultiABec.cpp | 38 +-- Source/radiation/HypreMultiABec.cpp | 1 - Source/radiation/Make.package | 2 - Source/radiation/RadSolve.cpp | 1 - 10 files changed, 235 insertions(+), 708 deletions(-) delete mode 100644 Source/radiation/HABEC_1D.F90 delete mode 100644 Source/radiation/HABEC_2D.F90 delete mode 100644 Source/radiation/HABEC_3D.F90 delete mode 100644 Source/radiation/HABEC_F.H diff --git a/Source/radiation/HABEC.H b/Source/radiation/HABEC.H index 2c9df6b2e3..3039e759d3 100644 --- a/Source/radiation/HABEC.H +++ b/Source/radiation/HABEC.H @@ -367,6 +367,222 @@ namespace HABEC } }); } + + AMREX_INLINE + void hdterm (Array4 const dterm, + Array4 const er, + const Box& reg, + int cdir, int bct, Real bcl, + Array4 const bcval, + Array4 const mask, + Array4 const d, + const Real* dx) + { + 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 + { + if (mask.contains(i+icp,j+jcp,k+kcp)) { + if (mask(i+icp,j+jcp,k+kcp) > 0) { + if (bct == LO_DIRICHLET) { + Real d_sign = 1.0_rt; + if (iep != 0 || jep != 0 || kep != 0) { + // right edge + d_sign = -1.0_rt; + } + dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl); + } + } + } + }); + } + + AMREX_INLINE + void hdterm3 (Array4 const dterm, + Array4 const er, + const Box& reg, + int cdir, int bctype, + Array4 const tf, + Real bcl, + Array4 const bcval, + Array4 const mask, + Array4 const d, + const Real* const dx) + { + 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 + { + if (mask.contains(i+icp,j+jcp,k+kcp)) { + if (mask(i+icp,j+jcp,k+kcp) > 0) { + int bct; + if (bctype == -1) { + bct = tf(i+icp,j+jcp,k+kcp); + } + else { + bct = bctype; + } + if (bct == LO_DIRICHLET) { + Real d_sign = 1.0_rt; + if (iep != 0 || jep != 0 || kep != 0) { + // right edge + d_sign = -1.0_rt; + } + dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl); + } + else if (bct == LO_NEUMANN && bcval(i+icp,j+jcp,k+kcp) == 0.0_rt) { + dterm(i+iep,j+jep,k+kep) = 0.0_rt; + } + } + } + }); + } } // namespace HABEC #endif diff --git a/Source/radiation/HABEC_1D.F90 b/Source/radiation/HABEC_1D.F90 deleted file mode 100644 index 544297195e..0000000000 --- a/Source/radiation/HABEC_1D.F90 +++ /dev/null @@ -1,133 +0,0 @@ -#include -#include - -module habec_module - - ! habec is Hypre abec, where abec is the form of the linear equation - ! we are solving: - ! - ! alpha*phi - div(beta*grad phi) + div(\vec{c}*phi) - - use amrex_fort_module, only : rt => amrex_real - implicit none - -contains - -subroutine hdterm(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bct - real(rt) :: bcl, dx(1) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i - h = dx(1) - if (bct == LO_DIRICHLET) then - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - if (mask(i-1) > 0) then - dterm(i) = d(i)*(er(i) - bcval(i-1))/(0.5e0_rt*h+bcl) - endif - else if (cdir == 1) then - ! Right face of grid - i = reg_h1 - if (mask(i+1) > 0) then - dterm(i+1) = d(i+1)*(bcval(i+1)-er(i))/(0.5e0_rt*h+bcl) - endif - else - print *, "hdterm: impossible face orientation" - endif - else - print *, "hdterm: unsupported boundary type" - stop - endif -end subroutine hdterm - -subroutine hdterm3(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bctype, tf(DIMV(bcv)) - real(rt) :: bcl, dx(1) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i, bct - h = dx(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 - dterm(i) = d(i)*(er(i) - bcval(i-1))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i-1) == 0.e0_rt) then - dterm(i) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i+1) = d(i+1)*(bcval(i+1)-er(i))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i+1) == 0.e0_rt) then - dterm(i+1) = 0.e0_rt - else - print *, "hbterm3: unsupported boundary type" - stop - endif - endif - else - print *, "hdterm3: impossible face orientation" - endif -end subroutine hdterm3 - -end module habec_module diff --git a/Source/radiation/HABEC_2D.F90 b/Source/radiation/HABEC_2D.F90 deleted file mode 100644 index 837e29d219..0000000000 --- a/Source/radiation/HABEC_2D.F90 +++ /dev/null @@ -1,206 +0,0 @@ -#include -#include - - -module habec_module - - ! habec is Hypre abec, where abec is the form of the linear equation - ! we are solving: - ! - ! alpha*phi - div(beta*grad phi) + div(\vec{c}*phi) - - use amrex_fort_module, only : rt => amrex_real - implicit none - -contains - -subroutine hdterm(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bct - real(rt) :: bcl, dx(2) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i, j - if (cdir == 0 .OR. cdir == 2) then - h = dx(1) - else - h = dx(2) - endif - if (bct == LO_DIRICHLET) then - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - do j = reg_l2, reg_h2 - if (mask(i-1,j) > 0) then - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i-1,j))/(0.5e0_rt*h+bcl) - 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 - dterm(i+1,j) = d(i+1,j)*(bcval(i+1,j)-er(i,j))/(0.5e0_rt*h+bcl) - 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 - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i,j-1))/(0.5e0_rt*h+bcl) - 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 - dterm(i,j+1) = d(i,j+1)*(bcval(i,j+1)-er(i,j))/(0.5e0_rt*h+bcl) - endif - enddo - else - print *, "hdterm: impossible face orientation" - endif - else - print *, "hdterm: unsupported boundary type" - stop - endif -end subroutine hdterm - -subroutine hdterm3(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bctype, tf(DIMV(bcv)) - real(rt) :: bcl, dx(2) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h - integer :: i, j, bct - if (cdir == 0 .OR. cdir == 2) then - h = dx(1) - 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 - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i-1,j))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i-1,j) == 0.e0_rt) then - dterm(i,j) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i+1,j) = d(i+1,j)*(bcval(i+1,j)-er(i,j))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i+1,j) == 0.e0_rt) then - dterm(i+1,j) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i,j) = d(i,j)*(er(i,j)-bcval(i,j-1))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i,j-1) == 0.e0_rt) then - dterm(i,j) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i,j+1) = d(i,j+1)*(bcval(i,j+1)-er(i,j))/(0.5e0_rt*h+bcl) - else if (bct == LO_NEUMANN .AND. bcval(i,j+1) == 0.e0_rt) then - dterm(i,j+1) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - else - print *, "hdterm3: impossible face orientation" - endif -end subroutine hdterm3 - -end module habec_module diff --git a/Source/radiation/HABEC_3D.F90 b/Source/radiation/HABEC_3D.F90 deleted file mode 100644 index d59edec61a..0000000000 --- a/Source/radiation/HABEC_3D.F90 +++ /dev/null @@ -1,308 +0,0 @@ -#include -#include - - -module habec_module - - ! habec is Hypre abec, where abec is the form of the linear equation - ! we are solving: - ! - ! alpha*phi - div(beta*grad phi) + div(\vec{c}*phi) - - use amrex_fort_module, only : rt => amrex_real - implicit none - -contains - -subroutine hdterm(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bct - real(rt) :: bcl, dx(3) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h, bfm, bfv - integer :: i, j, k - 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 (bct == LO_DIRICHLET) then - 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 - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i-1,j,k)) & - / (0.5e0_rt*h + bcl) - 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 - dterm(i+1,j,k) = d(i+1,j,k) * & - (bcval(i+1,j,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - 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 - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j-1,k)) & - / (0.5e0_rt*h + bcl) - 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 - dterm(i,j+1,k) = d(i,j+1,k) * & - (bcval(i,j+1,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - 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 - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j,k-1)) & - / (0.5e0_rt*h + bcl) - 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 - dterm(i,j,k+1) = d(i,j,k+1) * & - (bcval(i,j,k+1) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - endif - enddo - enddo - else - print *, "hdterm: impossible face orientation" - endif - else - print *, "hdterm: unsupported boundary type" - stop - endif -end subroutine hdterm - -subroutine hdterm3(dterm, & - DIMS(dtbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bctype, tf, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - d, DIMS(dbox), & - dx) bind(C, name="hdterm3") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(dtbox) - integer :: DIMDEC(ebox) - integer :: DIMDEC(reg) - integer :: DIMDEC(bcv) - integer :: DIMDEC(msk) - integer :: DIMDEC(dbox) - integer :: cdir, bctype, tf(DIMV(bcv)) - real(rt) :: bcl, dx(3) - real(rt) :: dterm(DIMV(dtbox)) - real(rt) :: er(DIMV(ebox)) - real(rt) :: bcval(DIMV(bcv)) - integer :: mask(DIMV(msk)) - real(rt) :: d(DIMV(dbox)) - real(rt) :: h, bfm, bfv - 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 - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i-1,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i-1,j,k) == 0.e0_rt) then - dterm(i,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i+1,j,k) = d(i+1,j,k) * & - (bcval(i+1,j,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i+1,j,k) == 0.e0_rt) then - dterm(i+1,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j-1,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j-1,k) == 0.e0_rt) then - dterm(i,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i,j+1,k) = d(i,j+1,k) * & - (bcval(i,j+1,k) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j+1,k) == 0.e0_rt) then - dterm(i,j+1,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i,j,k) = d(i,j,k) * & - (er(i,j,k) - bcval(i,j,k-1)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j,k-1) == 0.e0_rt) then - dterm(i,j,k) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - 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 - dterm(i,j,k+1) = d(i,j,k+1) * & - (bcval(i,j,k+1) - er(i,j,k)) & - / (0.5e0_rt*h + bcl) - else if (bct == LO_NEUMANN & - .AND. bcval(i,j,k+1) == 0.e0_rt) then - dterm(i,j,k+1) = 0.e0_rt - else - print *, "hdterm3: unsupported boundary type" - stop - endif - endif - enddo - enddo - else - print *, "hdterm3: impossible face orientation" - endif -end subroutine hdterm3 - -end module habec_module diff --git a/Source/radiation/HABEC_F.H b/Source/radiation/HABEC_F.H deleted file mode 100644 index 72a3f7be97..0000000000 --- a/Source/radiation/HABEC_F.H +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef CASTRO_HABEC_F_H -#define CASTRO_HABEC_F_H - -#include -#include - -#ifdef __cplusplus -extern "C" { -#else -#define RadBoundCond int -#endif - - void hdterm(BL_FORT_FAB_ARG(dterm), - BL_FORT_FAB_ARG(soln), - ARLIM_P(reglo), ARLIM_P(reghi), - const int& cdir, const RadBoundCond& bct, - const amrex::Real& bcl, - const BL_FORT_FAB_ARG(bcval), - const BL_FORT_IFAB_ARG(mask), - BL_FORT_FAB_ARG(dcoefs), - const amrex::Real* dx); - - void hdterm3(BL_FORT_FAB_ARG(dterm), - BL_FORT_FAB_ARG(soln), - ARLIM_P(reglo), ARLIM_P(reghi), - const int& cdir, const int& bctype, const int* tf, - const amrex::Real& bcl, - const BL_FORT_FAB_ARG(bcval), - const BL_FORT_IFAB_ARG(mask), - BL_FORT_FAB_ARG(dcoefs), - const amrex::Real* dx); - -#ifdef __cplusplus -}; -#endif - -#endif diff --git a/Source/radiation/HypreABec.cpp b/Source/radiation/HypreABec.cpp index e3ac5a4931..638d4715bb 100644 --- a/Source/radiation/HypreABec.cpp +++ b/Source/radiation/HypreABec.cpp @@ -4,7 +4,6 @@ #include #include -#include #include #include diff --git a/Source/radiation/HypreExtMultiABec.cpp b/Source/radiation/HypreExtMultiABec.cpp index 9f958fae89..a65292751a 100644 --- a/Source/radiation/HypreExtMultiABec.cpp +++ b/Source/radiation/HypreExtMultiABec.cpp @@ -1,6 +1,6 @@ #include -#include +#include #include #include <_hypre_sstruct_mv.h> @@ -1040,31 +1040,31 @@ void HypreExtMultiABec::boundaryDterm(int level, const Mask &msk = bd[level]->bndryMasks(oitr(), i); if (reg[oitr()] == domain[oitr()]) { - const int *tfp = NULL; + Array4 tf_arr; int bctype = bct; if (bd[level]->mixedBndry(oitr())) { const BaseFab &tf = *(bd[level]->bndryTypes(oitr())[i]); - tfp = tf.dataPtr(); + tf_arr = tf.array(); bctype = -1; } - hdterm3(BL_TO_FORTRAN(Dterm[idim][mfi]), - BL_TO_FORTRAN_N(Soln[mfi], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bctype, tfp, bcl, - BL_TO_FORTRAN_N(bcv, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*d2coefs[level])[idim][mfi]), - geom[level].CellSize()); + HABEC::hdterm3(Dterm[idim][mfi].array(), + Soln[mfi].array(icomp), + reg, + cdir, bctype, tf_arr, bcl, + bcv.array(bdcomp), + msk.array(), + (*d2coefs[level])[idim][mfi].array(), + geom[level].CellSize()); } else { - hdterm(BL_TO_FORTRAN(Dterm[idim][mfi]), - BL_TO_FORTRAN_N(Soln[mfi], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bct, bcl, - BL_TO_FORTRAN_N(bcv, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*d2coefs[level])[idim][mfi]), - geom[level].CellSize()); + HABEC::hdterm(Dterm[idim][mfi].array(), + Soln[mfi].array(icomp), + reg, + cdir, bct, bcl, + bcv.array(bdcomp), + msk.array(), + (*d2coefs[level])[idim][mfi].array(), + geom[level].CellSize()); } } } diff --git a/Source/radiation/HypreMultiABec.cpp b/Source/radiation/HypreMultiABec.cpp index b548767211..311c78c0a4 100644 --- a/Source/radiation/HypreMultiABec.cpp +++ b/Source/radiation/HypreMultiABec.cpp @@ -4,7 +4,6 @@ #include #include -#include #include #include diff --git a/Source/radiation/Make.package b/Source/radiation/Make.package index 1def860164..eb9eb54a20 100644 --- a/Source/radiation/Make.package +++ b/Source/radiation/Make.package @@ -30,10 +30,8 @@ CEXE_headers += filter.H CEXE_headers += filt_prim.H FEXE_headers += RAD_F.H -FEXE_headers += HABEC_F.H ca_F90EXE_sources += RAD_$(DIM)D.F90 -ca_F90EXE_sources += HABEC_$(DIM)D.F90 CEXE_sources += trace_ppm_rad.cpp diff --git a/Source/radiation/RadSolve.cpp b/Source/radiation/RadSolve.cpp index 12c8fd7b5f..ea8e78d0df 100644 --- a/Source/radiation/RadSolve.cpp +++ b/Source/radiation/RadSolve.cpp @@ -8,7 +8,6 @@ #include #include #include -#include // only for nonsymmetric flux; may be changed? #include