From df7ce4657f0718f5fe826554b2fb8c5d36c21706 Mon Sep 17 00:00:00 2001 From: Max Katz Date: Sun, 24 Sep 2023 10:16:39 -0400 Subject: [PATCH] Convert sphe to C++ (#2579) --- Source/radiation/RAD_1D.F90 | 15 ------------- Source/radiation/RAD_2D.F90 | 32 -------------------------- Source/radiation/RAD_3D.F90 | 32 -------------------------- Source/radiation/RAD_F.H | 3 --- Source/radiation/RadSolve.cpp | 42 ++++++++++++++++++++++++++++++----- 5 files changed, 37 insertions(+), 87 deletions(-) diff --git a/Source/radiation/RAD_1D.F90 b/Source/radiation/RAD_1D.F90 index f89fe8359f..bef1b8620e 100644 --- a/Source/radiation/RAD_1D.F90 +++ b/Source/radiation/RAD_1D.F90 @@ -11,21 +11,6 @@ module rad_module contains -subroutine sphe(r, s, n, & - DIMS(reg), dx) bind(C, name="sphe") - use amrex_fort_module, only : rt => amrex_real - implicit none - integer :: DIMDEC(reg) - real(rt) :: r(reg_l1:reg_h1) - real(rt) :: s(1) - integer :: n - real(rt) :: dx(1) - integer :: i - do i = reg_l1, reg_h1 - r(i) = r(i)**2 - enddo -end subroutine sphe - subroutine rfface(fine, & DIMS(fbox), & crse, & diff --git a/Source/radiation/RAD_2D.F90 b/Source/radiation/RAD_2D.F90 index 34f3e07d00..f32f5a8ea5 100644 --- a/Source/radiation/RAD_2D.F90 +++ b/Source/radiation/RAD_2D.F90 @@ -11,38 +11,6 @@ module rad_module contains -subroutine sphe(r, s, n, & - DIMS(reg), dx) bind(C, name="sphe") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(reg) - real(rt) :: r(reg_l1:reg_h1) - real(rt) :: s(reg_l2:reg_h2) - integer :: n - real(rt) :: dx(2) - real(rt) :: h1, h2, d1, d2 - integer :: i, j - if (n == 0) then - do i = reg_l1, reg_h1 - r(i) = r(i)**2 - enddo - h2 = 0.5e0_rt * dx(2) - d2 = 1.e0_rt / dx(2) - do j = reg_l2, reg_h2 - s(j) = d2 * (cos(s(j) - h2) - cos(s(j) + h2)) - enddo - else - h1 = 0.5e0_rt * dx(1) - d1 = 1.e0_rt / (3.e0_rt * dx(1)) - do i = reg_l1, reg_h1 - r(i) = d1 * ((r(i) + h1)**3 - (r(i) - h1)**3) - enddo - do j = reg_l2, reg_h2 - s(j) = sin(s(j)) - enddo - endif -end subroutine sphe - subroutine rfface(fine, & DIMS(fbox), & crse, & diff --git a/Source/radiation/RAD_3D.F90 b/Source/radiation/RAD_3D.F90 index f3245fe54b..ebdbdee21a 100644 --- a/Source/radiation/RAD_3D.F90 +++ b/Source/radiation/RAD_3D.F90 @@ -14,38 +14,6 @@ module rad_module ! The following routines implement metric terms in 2D and are included ! in the 3D source only to enable the code to link. -subroutine sphe(r, s, n, & - DIMS(reg), dx) bind(C, name="sphe") - - use amrex_fort_module, only : rt => amrex_real - integer :: DIMDEC(reg) - real(rt) :: r(reg_l1:reg_h1) - real(rt) :: s(reg_l2:reg_h2) - integer :: n - real(rt) :: dx(2) - real(rt) :: h1, h2, d1, d2 - integer :: i, j - if (n == 0) then - do i = reg_l1, reg_h1 - r(i) = r(i)**2 - enddo - h2 = 0.5e0_rt * dx(2) - d2 = 1.e0_rt / dx(2) - do j = reg_l2, reg_h2 - s(j) = d2 * (cos(s(j) - h2) - cos(s(j) + h2)) - enddo - else - h1 = 0.5e0_rt * dx(1) - d1 = 1.e0_rt / (3.e0_rt * dx(1)) - do i = reg_l1, reg_h1 - r(i) = d1 * ((r(i) + h1)**3 - (r(i) - h1)**3) - enddo - do j = reg_l2, reg_h2 - s(j) = sin(s(j)) - enddo - endif -end subroutine sphe - subroutine rfface(fine, & DIMS(fbox), & crse, & diff --git a/Source/radiation/RAD_F.H b/Source/radiation/RAD_F.H index be3fe79d23..bb8e1f5d86 100644 --- a/Source/radiation/RAD_F.H +++ b/Source/radiation/RAD_F.H @@ -126,9 +126,6 @@ extern "C" #ifdef __cplusplus extern "C" { #endif - void sphe(amrex::Real* r, amrex::Real* s, const int& idim, - ARLIM_P(edgeboxlo), ARLIM_P(edgeboxhi), const amrex::Real* dx); - void fkpn(const int* lo, const int* hi, BL_FORT_FAB_ARG_3D(fkp), amrex::Real con, amrex::Real em, amrex::Real en, diff --git a/Source/radiation/RadSolve.cpp b/Source/radiation/RadSolve.cpp index fb8926f5b5..2789476e07 100644 --- a/Source/radiation/RadSolve.cpp +++ b/Source/radiation/RadSolve.cpp @@ -877,9 +877,18 @@ void RadSolve::levelDterm(int level, MultiFab& Dterm, MultiFab& Er, int igroup) parent->Geom(level).GetCellLoc(rc, reg, 0); parent->Geom(level).GetCellLoc(s, reg, 0); const Box &dbox = Dterm_face[0][fi].box(); - sphe(re.dataPtr(), s.dataPtr(), 0, - ARLIM(dbox.loVect()), ARLIM(dbox.hiVect()), dx.data()); - + + for (int i = dbox.loVect()[0]; i <= dbox.hiVect()[0]; ++i) { + re[i] *= re[i]; + } +#if AMREX_SPACEDIM >= 2 + Real h2 = 0.5e0_rt * dx[1]; + Real d2 = 1.e0_rt / dx[1]; + for (int j = dbox.loVect()[1]; j <= dbox.hiVect()[1]; ++j) { + s[j] = d2 * (std::cos(s[j] - h2) - std::cos(s[j] + h2)); + } +#endif + ca_correct_dterm(D_DECL(BL_TO_FORTRAN(Dterm_face[0][fi]), BL_TO_FORTRAN(Dterm_face[1][fi]), BL_TO_FORTRAN(Dterm_face[2][fi])), @@ -1160,8 +1169,31 @@ void RadSolve::getEdgeMetric(int idim, const Geometry& geom, const Box& edgebox, geom.GetEdgeLoc(s, reg, I); } const Real *dx = geom.CellSize(); - sphe(r.dataPtr(), s.dataPtr(), idim, - ARLIM(edgebox.loVect()), ARLIM(edgebox.hiVect()), dx); + + if (idim == 0) { + for (int i = edgebox.loVect()[0]; i <= edgebox.hiVect()[0]; ++i) { + r[i] *= r[i]; + } +#if AMREX_SPACEDIM >= 2 + Real h2 = 0.5e0_rt * dx[1]; + Real d2 = 1.e0_rt / dx[1]; + for (int j = edgebox.loVect()[1]; j <= edgebox.hiVect()[1]; ++j) { + s[j] = d2 * (std::cos(s[j] - h2) - std::cos(s[j] + h2)); + } +#endif + } + else { + Real h1 = 0.5e0_rt * dx[0]; + Real d1 = 1.e0_rt / (3.e0_rt * dx[0]); + for (int i = edgebox.loVect()[0]; i <= edgebox.hiVect()[0]; ++i) { + r[i] = d1 * (std::pow(r[i] + h1, 3) - std::pow(r[i] - h1, 3)); + } +#if AMREX_SPACEDIM >= 2 + for (int j = edgebox.loVect()[1]; j <= edgebox.hiVect()[1]; ++j) { + s[j] = std::sin(s[j]); + } +#endif + } } }