Skip to content

Commit

Permalink
Convert hbflx to C++
Browse files Browse the repository at this point in the history
  • Loading branch information
maxpkatz committed Sep 24, 2023
1 parent df7ce46 commit b208798
Show file tree
Hide file tree
Showing 8 changed files with 222 additions and 335 deletions.
203 changes: 203 additions & 0 deletions Source/radiation/HABEC.H
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
#ifndef _HABEC_H_
#define _HABEC_H_

#include <AMReX_Array4.H>
#include <AMReX_REAL.H>

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<Real> const flux,
Array4<Real> const er,
const Box& reg,
int cdir, int bct, int bho, Real bcl,
Array4<Real const> const bcval,
Array4<int const> const mask,
Array4<Real const> 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
69 changes: 0 additions & 69 deletions Source/radiation/HABEC_1D.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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), &
Expand Down
99 changes: 0 additions & 99 deletions Source/radiation/HABEC_2D.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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), &
Expand Down
Loading

0 comments on commit b208798

Please sign in to comment.