Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Convert hbflx to C++ #2583

Merged
merged 1 commit into from
Sep 25, 2023
Merged

Convert hbflx to C++ #2583

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
203 changes: 203 additions & 0 deletions Source/radiation/HABEC.H
View file
Open in desktop
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
View file
Open in desktop
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
View file
Open in desktop
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