diff --git a/Source/radiation/HABEC.H b/Source/radiation/HABEC.H new file mode 100644 index 0000000000..0b3d718099 --- /dev/null +++ b/Source/radiation/HABEC.H @@ -0,0 +1,203 @@ +#ifndef _HABEC_H_ +#define _HABEC_H_ + +#include +#include + +using namespace amrex; + +namespace HABEC +{ + // 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) + + AMREX_INLINE + void hbflx (Array4 const flux, + Array4 const er, + const Box& reg, + int cdir, int bct, int bho, Real bcl, + Array4 const bcval, + Array4 const mask, + Array4 const b, + Real beta, const Real* dx, int inhom) + { + Real h; + + bool x_left = false; + bool x_right = false; + bool y_left = false; + bool y_right = false; + bool z_left = false; + bool z_right = false; + +#if AMREX_SPACEDIM == 1 + if (cdir == 0) { + h = dx[0]; + x_left= true; + } + else if (cdir == 1) { + h = dx[0]; + x_right= true; + } +#elif AMREX_SPACEDIM == 2 + if (cdir == 0) { + h = dx[0]; + x_left= true; + } + else if (cdir == 2) { + h = dx[0]; + x_right= true; + } + else if (cdir == 1) { + h = dx[1]; + y_left= true; + } + else if (cdir == 3) { + h = dx[1]; + y_right= true; + } +#else + if (cdir == 0) { + h = dx[0]; + x_left= true; + } + else if (cdir == 3) { + h = dx[0]; + x_right= true; + } + else if (cdir == 1) { + h = dx[1]; + y_left = true; + } + else if (cdir == 4) { + h = dx[1]; + y_right = true; + } + else if (cdir == 2) { + h = dx[2]; + z_left = true; + } + else if (cdir == 5) { + h = dx[2]; + z_right = true; + } +#endif + + Real bfv, bfm, bfm2; + + if (bct == LO_DIRICHLET) { + if (bho >= 1) { + Real h2 = 0.5e0_rt * h; + Real th2 = 3.e0_rt * h2; + bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)); + bfm = (beta / h) * (th2 - bcl) / (bcl + h2); + bfm2 = (beta / h) * (bcl - h2) / (bcl + th2); + } + else { + bfv = beta / (0.5e0_rt * h + bcl); + bfm = bfv; + } + } + else { + amrex::Error("hbflx: unsupported boundary type"); + } + + if (inhom == 0) { + bfv = 0.e0_rt; + } + + int reg_ilo = reg.loVect3d()[0]; + int reg_ihi = reg.hiVect3d()[0]; + int reg_jlo = reg.loVect3d()[1]; + int reg_jhi = reg.hiVect3d()[1]; + int reg_klo = reg.loVect3d()[2]; + int reg_khi = reg.hiVect3d()[2]; + + if (x_left) { + int i = reg_ilo; + for (int k = reg_klo; k <= reg_khi; ++k) { + for (int j = reg_jlo; j <= reg_jhi; ++j) { + if (mask(i-1,j,k) > 0) { + flux(i,j,k) = b(i,j,k) * (bfv * bcval(i-1,j,k) - bfm * er(i,j,k)); + if (bho >= 1) { + flux(i,j,k) = flux(i,j,k) - b(i,j,k) * bfm2 * er(i+1,j,k); + } + } + } + } + } + else if (x_right) { + int i = reg_ihi; + for (int k = reg_klo; k <= reg_khi; ++k) { + for (int j = reg_jlo; j <= reg_jhi; ++j) { + if (mask(i+1,j,k) > 0) { + flux(i+1,j,k) = -b(i+1,j,k) * (bfv * bcval(i+1,j,k) - bfm * er(i,j,k)); + if (bho >= 1) { + flux(i+1,j,k) = flux(i+1,j,k) + b(i+1,j,k) * bfm2 * er(i-1,j,k); + } + } + } + } + } + else if (y_left) { + int j = reg_jlo; + for (int k = reg_klo; k <= reg_khi; ++k) { + for (int i = reg_ilo; i <= reg_ihi; ++i) { + if (mask(i,j-1,k) > 0) { + flux(i,j,k) = b(i,j,k) * (bfv * bcval(i,j-1,k) - bfm * er(i,j,k)); + if (bho >= 1) { + flux(i,j,k) = flux(i,j,k) - b(i,j,k) * bfm2 * er(i,j+1,k); + } + } + } + } + } + else if (y_right) { + int j = reg_jhi; + for (int k = reg_klo; k <= reg_khi; ++k) { + for (int i = reg_ilo; i <= reg_ihi; ++i) { + if (mask(i,j+1,k) > 0) { + flux(i,j+1,k) = -b(i,j+1,k) * (bfv * bcval(i,j+1,k) - bfm * er(i,j,k)); + if (bho >= 1) { + flux(i,j+1,k) = flux(i,j+1,k) + b(i,j+1,k) * bfm2 * er(i,j-1,k); + } + } + } + } + } + else if (z_left) { + int k = reg_klo; + for (int j = reg_jlo; j <= reg_jhi; ++j) { + for (int i = reg_ilo; i <= reg_ihi; ++i) { + if (mask(i,j,k-1) > 0) { + flux(i,j,k) = b(i,j,k) * (bfv * bcval(i,j,k-1) - bfm * er(i,j,k)); + if (bho >= 1) { + flux(i,j,k) = flux(i,j,k) - b(i,j,k) * bfm2 * er(i,j,k+1); + } + } + } + } + } + else if (z_right) { + int k = reg_khi; + for (int j = reg_jlo; j <= reg_jhi; ++j) { + for (int i = reg_ilo; i <= reg_ihi; ++i) { + if (mask(i,j,k+1) > 0) { + flux(i,j,k+1) = -b(i,j,k+1) * (bfv * bcval(i,j,k+1) - bfm * er(i,j,k)); + if (bho >= 1) { + flux(i,j,k+1) = flux(i,j,k+1) + b(i,j,k+1) * bfm2 * er(i,j,k-1); + } + } + } + } + } + else { + std::cout << "hbflx: impossible face orientation" << std::endl; + } + } + +} // namespace HABEC + +#endif diff --git a/Source/radiation/HABEC_1D.F90 b/Source/radiation/HABEC_1D.F90 index a8d2f446b4..3f2aafc52b 100644 --- a/Source/radiation/HABEC_1D.F90 +++ b/Source/radiation/HABEC_1D.F90 @@ -13,75 +13,6 @@ module habec_module contains -subroutine hbflx(flux, & - DIMS(fbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bho, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - b, DIMS(bbox), & - beta, dx, inhom) bind(C, name="hbflx") - - 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 :: cdir, bct, bho, inhom - real(rt) :: bcl, beta, dx(1) - 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) :: h, bfm, bfv - real(rt) :: bfm2, h2, th2 - integer :: i - h = dx(1) - 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)) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) - else - bfv = beta / (0.5e0_rt * h + bcl) - bfm = bfv - endif - else - print *, "hbflx: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - endif - if (cdir == 0) then - ! Left face of grid - i = reg_l1 - if (mask(i-1) > 0) then - flux(i) = b(i) * (bfv * bcval(i-1) - bfm * er(i)) - if (bho >= 1) then - flux(i) = flux(i) - b(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 - flux(i+1) = -b(i+1) * (bfv * bcval(i+1) - bfm * er(i)) - if (bho >= 1) then - flux(i+1) = flux(i+1) + b(i+1) * bfm2 * er(i-1) - endif - endif - else - print *, "hbflx: impossible face orientation" - endif -end subroutine hbflx - subroutine hbflx3(flux, & DIMS(fbox), & er, DIMS(ebox), & diff --git a/Source/radiation/HABEC_2D.F90 b/Source/radiation/HABEC_2D.F90 index c1f37450b8..2119c4150c 100644 --- a/Source/radiation/HABEC_2D.F90 +++ b/Source/radiation/HABEC_2D.F90 @@ -14,105 +14,6 @@ module habec_module contains -subroutine hbflx(flux, & - DIMS(fbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bho, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - b, DIMS(bbox), & - beta, dx, inhom) bind(C, name="hbflx") - - 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 :: cdir, bct, bho, inhom - real(rt) :: bcl, beta, dx(2) - 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) :: h, bfm, bfv - real(rt) :: bfm2, h2, th2 - integer :: i, j - if (cdir == 0 .OR. cdir == 2) then - h = dx(1) - else - h = dx(2) - 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)) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) - else - bfv = beta / (0.5e0_rt * h + bcl) - bfm = bfv - endif - else - print *, "hbflx: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - 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 - flux(i,j) = b(i,j) * (bfv * bcval(i-1,j) - bfm * er(i,j)) - if (bho >= 1) then - flux(i,j) = flux(i,j) - b(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 - flux(i+1,j) = -b(i+1,j) * (bfv * bcval(i+1,j) - bfm * er(i,j)) - if (bho >= 1) then - flux(i+1,j) = flux(i+1,j) + b(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 - flux(i,j) = b(i,j) * (bfv * bcval(i,j-1) - bfm * er(i,j)) - if (bho >= 1) then - flux(i,j) = flux(i,j) - b(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 - flux(i,j+1) = -b(i,j+1) * (bfv * bcval(i,j+1) - bfm * er(i,j)) - if (bho >= 1) then - flux(i,j+1) = flux(i,j+1) + b(i,j+1) * bfm2 * er(i,j-1) - endif - endif - enddo - else - print *, "hbflx: impossible face orientation" - endif -end subroutine hbflx - subroutine hbflx3(flux, & DIMS(fbox), & er, DIMS(ebox), & diff --git a/Source/radiation/HABEC_3D.F90 b/Source/radiation/HABEC_3D.F90 index 17aa922c66..dcd02ebf88 100644 --- a/Source/radiation/HABEC_3D.F90 +++ b/Source/radiation/HABEC_3D.F90 @@ -14,147 +14,6 @@ module habec_module contains -subroutine hbflx(flux, & - DIMS(fbox), & - er, DIMS(ebox), & - DIMS(reg), & - cdir, bct, bho, bcl, & - bcval, DIMS(bcv), & - mask, DIMS(msk), & - b, DIMS(bbox), & - beta, dx, inhom) bind(C, name="hbflx") - - 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 :: cdir, bct, bho, inhom - real(rt) :: bcl, beta, dx(3) - 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) :: h, bfm, bfv - real(rt) :: bfm2, h2, th2 - 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 (bho >= 1) then - h2 = 0.5e0_rt * h - th2 = 3.e0_rt * h2 - bfv = 2.e0_rt * beta * h / ((bcl + h2) * (bcl + th2)) - bfm = (beta / h) * (th2 - bcl) / (bcl + h2) - bfm2 = (beta / h) * (bcl - h2) / (bcl + th2) - else - bfv = beta / (0.5e0_rt * h + bcl) - bfm = bfv - endif - else - print *, "hbflx: unsupported boundary type" - stop - endif - if (inhom == 0) then - bfv = 0.e0_rt - 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 - flux(i,j,k) = b(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) - & - b(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 - flux(i+1,j,k) = -b(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) + & - b(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 - flux(i,j,k) = b(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) - & - b(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 - flux(i,j+1,k) = -b(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) + & - b(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 - flux(i,j,k) = b(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) - & - b(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 - flux(i,j,k+1) = -b(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) + & - b(i,j,k+1) * bfm2 * er(i,j,k-1) - endif - endif - enddo - enddo - else - print *, "hbflx: impossible face orientation" - endif -end subroutine hbflx - subroutine hbflx3(flux, & DIMS(fbox), & er, DIMS(ebox), & diff --git a/Source/radiation/HABEC_F.H b/Source/radiation/HABEC_F.H index 4df2a1318c..1ad5995275 100644 --- a/Source/radiation/HABEC_F.H +++ b/Source/radiation/HABEC_F.H @@ -10,16 +10,6 @@ extern "C" { #define RadBoundCond int #endif - void hbflx(BL_FORT_FAB_ARG(flux), - BL_FORT_FAB_ARG(soln), - ARLIM_P(reglo), ARLIM_P(reghi), - const int& cdir, const RadBoundCond& bct, - const int& bho, const amrex::Real& bcl, - const amrex::Real* bcval, ARLIM_P(fslo), ARLIM_P(fshi), - const BL_FORT_IFAB_ARG(mask), - BL_FORT_FAB_ARG(bcoefs), - const amrex::Real& beta, const amrex::Real* dx, const int& inhom); - void hbflx3(BL_FORT_FAB_ARG(flux), BL_FORT_FAB_ARG(soln), ARLIM_P(reglo), ARLIM_P(reghi), diff --git a/Source/radiation/HypreABec.cpp b/Source/radiation/HypreABec.cpp index 9d24aef815..f1a8213986 100644 --- a/Source/radiation/HypreABec.cpp +++ b/Source/radiation/HypreABec.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -307,14 +308,14 @@ void HypreABec::boundaryFlux(MultiFab* Flux, MultiFab& Soln, int icomp, pSPa, ARLIM(SPabox.loVect()), ARLIM(SPabox.hiVect())); } else { - hbflx(BL_TO_FORTRAN(Flux[idim][si]), - BL_TO_FORTRAN_N(Soln[si], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bct, bho, bcl, - BL_TO_FORTRAN_N(fs, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*bcoefs[idim])[si]), - beta, dx, inhom); + HABEC::hbflx(Flux[idim][si].array(), + Soln[si].array(icomp), + reg, + cdir, bct, bho, bcl, + fs.array(bdcomp), + msk.array(), + (*bcoefs[idim])[si].array(), + beta, dx, inhom); } } } diff --git a/Source/radiation/HypreMultiABec.cpp b/Source/radiation/HypreMultiABec.cpp index a151205ce7..bc98e86bbd 100644 --- a/Source/radiation/HypreMultiABec.cpp +++ b/Source/radiation/HypreMultiABec.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -3960,14 +3961,14 @@ void HypreMultiABec::boundaryFlux(int level, pSPa, ARLIM(SPabox.loVect()), ARLIM(SPabox.hiVect())); } else { - hbflx(BL_TO_FORTRAN(Flux[idim][mfi]), - BL_TO_FORTRAN_N(Soln[mfi], icomp), - ARLIM(reg.loVect()), ARLIM(reg.hiVect()), - cdir, bct, bho, bcl, - BL_TO_FORTRAN_N(fs, bdcomp), - BL_TO_FORTRAN(msk), - BL_TO_FORTRAN((*bcoefs[level])[idim][mfi]), - beta, geom[level].CellSize(), inhom); + HABEC::hbflx(Flux[idim][mfi].array(), + Soln[mfi].array(icomp), + reg, + cdir, bct, bho, bcl, + fs.array(bdcomp), + msk.array(), + (*bcoefs[level])[idim][mfi].array(), + beta, geom[level].CellSize(), inhom); } } } diff --git a/Source/radiation/Make.package b/Source/radiation/Make.package index 4b0479540b..f564af29cf 100644 --- a/Source/radiation/Make.package +++ b/Source/radiation/Make.package @@ -25,6 +25,7 @@ CEXE_headers += RadSolve.H CEXE_headers += RadBndry.H CEXE_headers += RadTypes.H CEXE_headers += MGRadBndry.H +CEXE_headers += HABEC.H CEXE_headers += filter.H CEXE_headers += filt_prim.H