Skip to content

Commit

Permalink
Give FlashFluxRegisters ways to accumulate data in registers
Browse files Browse the repository at this point in the history 
I introduce variant methods that allow adding to previously
stored data in Flash(-X) flux registers, rather than just
copying over the older data.

For Flash-X'ers: Something like this is needed to get flux
correction working in the nontelescoping variant of the current
Spark Hydro implementation.
  • Loading branch information
@kweide
kweide committed Oct 18, 2023
1 parent 7ee2912 commit d77d93b
Show file tree
Hide file tree
Showing 4 changed files with 399 additions and 10 deletions.
12 changes: 12 additions & 0 deletions Src/F_Interfaces/AmrCore/AMReX_FlashFluxRegister.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,18 @@ public:
void store (int fine_global_index, int dir, FArrayBox const& fine_flux, FArrayBox const& area,
const int* isFluxDensity, Real sf);

// flux_in_register += scaling_factor * \sum{fine_flux} / (ref_ratio[0]*ref_ratio[1]*ref_ratio[2])
void add (int fine_global_index, int dir, FArrayBox const& fine_flux, Real sf);

// flux_in_register += scaling_factor * \sum{fine_flux * area}
void add (int fine_global_index, int dir, FArrayBox const& fine_flux, FArrayBox const& area,
Real sf);

// flux_in_register += scaling_factor * \sum{fine_flux * area}, if the component is flux density
// scaling_factor * \sum{fine_flux} , otherwise
void add (int fine_global_index, int dir, FArrayBox const& fine_flux, FArrayBox const& area,
const int* isFluxDensity, Real sf);

void communicate ();

// crse_flux = flux_in_register * scaling_factor
Expand Down
277 changes: 277 additions & 0 deletions Src/F_Interfaces/AmrCore/AMReX_FlashFluxRegister.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -453,6 +453,283 @@ void FlashFluxRegister::store (int fine_global_index, int dir, FArrayBox const&
}
}

void FlashFluxRegister::add (int fine_global_index, int dir, FArrayBox const& fine_flux, Real sf)
{
AMREX_ASSERT(dir < AMREX_SPACEDIM);
auto found = m_fine_map.find(fine_global_index);
if (found != m_fine_map.end()) {
const int ncomp = m_ncomp;
Array<FArrayBox*,AMREX_SPACEDIM> const& fab_a = found->second;
if (fab_a[dir]) {
Box const& b = fab_a[dir]->box();
Array4<Real> const& dest = fab_a[dir]->array();
Array4<Real const> const& src = fine_flux.const_array();
if (dir == 0) {
#if (AMREX_SPACEDIM == 1)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(j,k);
dest(i,0,0,n) += src(2*i,0,0,n)*sf;
});
#endif
#if (AMREX_SPACEDIM == 2)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(k);
dest(i,j,0,n) += (src(2*i,2*j ,0,n) +
src(2*i,2*j+1,0,n)) * (Real(0.5)*sf);
});
#endif
#if (AMREX_SPACEDIM == 3)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
dest(i,j,k,n) += (src(2*i,2*j ,2*k ,n) +
src(2*i,2*j+1,2*k ,n) +
src(2*i,2*j ,2*k+1,n) +
src(2*i,2*j+1,2*k+1,n)) * (Real(0.25)*sf);
});
#endif
}
#if (AMREX_SPACEDIM >= 2)
else if (dir == 1) {
#if (AMREX_SPACEDIM == 2)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(k);
dest(i,j,0,n) += (src(2*i ,2*j,0,n) +
src(2*i+1,2*j,0,n)) * (Real(0.5)*sf);
});
#endif
#if (AMREX_SPACEDIM == 3)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
dest(i,j,k,n) += (src(2*i ,2*j,2*k ,n) +
src(2*i+1,2*j,2*k ,n) +
src(2*i ,2*j,2*k+1,n) +
src(2*i+1,2*j,2*k+1,n)) * (Real(0.25)*sf);
});
#endif
}
#if (AMREX_SPACEDIM == 3)
else {
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
dest(i,j,k,n) += (src(2*i ,2*j ,2*k,n) +
src(2*i+1,2*j ,2*k,n) +
src(2*i ,2*j+1,2*k,n) +
src(2*i+1,2*j+1,2*k,n)) * (Real(0.25)*sf);
});
}
#endif
#endif
}
}
}

void FlashFluxRegister::add (int fine_global_index, int dir, FArrayBox const& fine_flux,
FArrayBox const& fine_area, Real sf)
{
AMREX_ASSERT(dir < AMREX_SPACEDIM);
auto found = m_fine_map.find(fine_global_index);
if (found != m_fine_map.end()) {
const int ncomp = m_ncomp;
Array<FArrayBox*,AMREX_SPACEDIM> const& fab_a = found->second;
if (fab_a[dir]) {
Box const& b = fab_a[dir]->box();
Array4<Real> const& dest = fab_a[dir]->array();
Array4<Real const> const& src = fine_flux.const_array();
Array4<Real const> const& area = fine_area.const_array();
if (dir == 0) {
#if (AMREX_SPACEDIM == 1)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(j,k);
dest(i,0,0,n) += src(2*i,0,0,n)*area(2*i,0,0)*sf;
});
#endif
#if (AMREX_SPACEDIM == 2)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(k);
dest(i,j,0,n) += (src(2*i,2*j ,0,n)*area(2*i,2*j ,0) +
src(2*i,2*j+1,0,n)*area(2*i,2*j+1,0)) * sf;
});
#endif
#if (AMREX_SPACEDIM == 3)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
dest(i,j,k,n) += (src(2*i,2*j ,2*k ,n)*area(2*i,2*j ,2*k ) +
src(2*i,2*j+1,2*k ,n)*area(2*i,2*j+1,2*k ) +
src(2*i,2*j ,2*k+1,n)*area(2*i,2*j ,2*k+1) +
src(2*i,2*j+1,2*k+1,n)*area(2*i,2*j+1,2*k+1)) * sf;
});
#endif
}
#if (AMREX_SPACEDIM >= 2)
else if (dir == 1) {
#if (AMREX_SPACEDIM == 2)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(k);
dest(i,j,0,n) += (src(2*i ,2*j,0,n)*area(2*i ,2*j,0) +
src(2*i+1,2*j,0,n)*area(2*i+1,2*j,0)) * sf;
});
#endif
#if (AMREX_SPACEDIM == 3)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
dest(i,j,k,n) += (src(2*i ,2*j,2*k ,n)*area(2*i ,2*j,2*k ) +
src(2*i+1,2*j,2*k ,n)*area(2*i+1,2*j,2*k ) +
src(2*i ,2*j,2*k+1,n)*area(2*i ,2*j,2*k+1) +
src(2*i+1,2*j,2*k+1,n)*area(2*i+1,2*j,2*k+1)) * sf;
});
#endif
}
#if (AMREX_SPACEDIM == 3)
else {
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
dest(i,j,k,n) += (src(2*i ,2*j ,2*k,n)*area(2*i ,2*j ,2*k) +
src(2*i+1,2*j ,2*k,n)*area(2*i+1,2*j ,2*k) +
src(2*i ,2*j+1,2*k,n)*area(2*i ,2*j+1,2*k) +
src(2*i+1,2*j+1,2*k,n)*area(2*i+1,2*j+1,2*k)) * sf;
});
}
#endif
#endif
}
}
}

void FlashFluxRegister::add (int fine_global_index, int dir, FArrayBox const& fine_flux,
FArrayBox const& fine_area, const int* isFluxDensity, Real sf)
{
auto& h_ifd = m_h_ifd[OpenMP::get_thread_num()];
auto& d_ifd = m_d_ifd[OpenMP::get_thread_num()];

AMREX_ASSERT(dir < AMREX_SPACEDIM);
auto found = m_fine_map.find(fine_global_index);
if (found != m_fine_map.end()) {
const int ncomp = m_ncomp;
Array<FArrayBox*,AMREX_SPACEDIM> const& fab_a = found->second;
if (fab_a[dir]) {
bool allsame = true;
for (int n = 0; n < m_ncomp; ++n) {
if (h_ifd[n] != isFluxDensity[n]) {
allsame = false;
h_ifd[n] = isFluxDensity[n];
}
}
if (d_ifd.empty()) {
allsame = false;
d_ifd.resize(m_ncomp);
}
if (! allsame) {
Gpu::copyAsync(Gpu::HostToDevice(), h_ifd.begin(), h_ifd.end(), d_ifd.begin());
}

Box const& b = fab_a[dir]->box();
Array4<Real> const& dest = fab_a[dir]->array();
Array4<Real const> const& src = fine_flux.const_array();
Array4<Real const> const& area = fine_area.const_array();
const int* ifd = d_ifd.data();
if (dir == 0) {
#if (AMREX_SPACEDIM == 1)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(j,k);
if (ifd[n]) {
dest(i,0,0,n) += src(2*i,0,0,n)*area(2*i,0,0)*sf;
} else {
dest(i,0,0,n) += src(2*i,0,0,n)*sf;
}
});
#endif
#if (AMREX_SPACEDIM == 2)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(k);
if (ifd[n]) {
dest(i,j,0,n) += (src(2*i,2*j ,0,n)*area(2*i,2*j ,0) +
src(2*i,2*j+1,0,n)*area(2*i,2*j+1,0)) * sf;
} else {
dest(i,j,0,n) += (src(2*i,2*j ,0,n) +
src(2*i,2*j+1,0,n)) * sf;
}
});
#endif
#if (AMREX_SPACEDIM == 3)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
if (ifd[n]) {
dest(i,j,k,n) += (src(2*i,2*j ,2*k ,n)*area(2*i,2*j ,2*k ) +
src(2*i,2*j+1,2*k ,n)*area(2*i,2*j+1,2*k ) +
src(2*i,2*j ,2*k+1,n)*area(2*i,2*j ,2*k+1) +
src(2*i,2*j+1,2*k+1,n)*area(2*i,2*j+1,2*k+1)) * sf;
} else {
dest(i,j,k,n) += (src(2*i,2*j ,2*k ,n) +
src(2*i,2*j+1,2*k ,n) +
src(2*i,2*j ,2*k+1,n) +
src(2*i,2*j+1,2*k+1,n)) * sf;
}
});
#endif
}
#if (AMREX_SPACEDIM >= 2)
else if (dir == 1) {
#if (AMREX_SPACEDIM == 2)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
amrex::ignore_unused(k);
if (ifd[n]) {
dest(i,j,0,n) += (src(2*i ,2*j,0,n)*area(2*i ,2*j,0) +
src(2*i+1,2*j,0,n)*area(2*i+1,2*j,0)) * sf;
} else {
dest(i,j,0,n) += (src(2*i ,2*j,0,n) +
src(2*i+1,2*j,0,n)) * sf;
}
});
#endif
#if (AMREX_SPACEDIM == 3)
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
if (ifd[n]) {
dest(i,j,k,n) += (src(2*i ,2*j,2*k ,n)*area(2*i ,2*j,2*k ) +
src(2*i+1,2*j,2*k ,n)*area(2*i+1,2*j,2*k ) +
src(2*i ,2*j,2*k+1,n)*area(2*i ,2*j,2*k+1) +
src(2*i+1,2*j,2*k+1,n)*area(2*i+1,2*j,2*k+1)) * sf;
} else {
dest(i,j,k,n) += (src(2*i ,2*j,2*k ,n) +
src(2*i+1,2*j,2*k ,n) +
src(2*i ,2*j,2*k+1,n) +
src(2*i+1,2*j,2*k+1,n)) * sf;
}
});
#endif
}
#if (AMREX_SPACEDIM == 3)
else {
AMREX_HOST_DEVICE_PARALLEL_FOR_4D (b, ncomp, i, j, k, n,
{
if (ifd[n]) {
dest(i,j,k,n) += (src(2*i ,2*j ,2*k,n)*area(2*i ,2*j ,2*k) +
src(2*i+1,2*j ,2*k,n)*area(2*i+1,2*j ,2*k) +
src(2*i ,2*j+1,2*k,n)*area(2*i ,2*j+1,2*k) +
src(2*i+1,2*j+1,2*k,n)*area(2*i+1,2*j+1,2*k)) * sf;
} else {
dest(i,j,k,n) += (src(2*i ,2*j ,2*k,n) +
src(2*i+1,2*j ,2*k,n) +
src(2*i ,2*j+1,2*k,n) +
src(2*i+1,2*j+1,2*k,n)) * sf;
}
});
}
#endif
#endif
}
}
}

void FlashFluxRegister::communicate ()
{
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
Expand Down
35 changes: 35 additions & 0 deletions Src/F_Interfaces/AmrCore/AMReX_flash_fluxregister_fi.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,41 @@ extern "C" {
flux_reg->store(fgid, dir, fab, areafab, ifd, scaling_factor);
}

void amrex_fi_flash_fluxregister_add (FlashFluxRegister* flux_reg, int fgid, int dir,
Real const* flux, const int* flo, const int* fhi, int nc,
Real scaling_factor)
{
Box bx;
bx = Box(IntVect(flo), IntVect(fhi));
bx.shiftHalf(dir,-1);
const FArrayBox fab(bx,nc,const_cast<Real*>(flux));
flux_reg->add(fgid, dir, fab, scaling_factor);
}

void amrex_fi_flash_fluxregister_add_area (FlashFluxRegister* flux_reg, int fgid, int dir,
Real const* flux, const int* flo, const int* fhi, int nc,
Real const* area, Real scaling_factor)
{
Box bx;
bx = Box(IntVect(flo), IntVect(fhi));
bx.shiftHalf(dir,-1);
const FArrayBox fab(bx,nc,const_cast<Real*>(flux));
const FArrayBox areafab(bx,1,const_cast<Real*>(area));
flux_reg->add(fgid, dir, fab, areafab, scaling_factor);
}

void amrex_fi_flash_fluxregister_add_area_ifd (FlashFluxRegister* flux_reg, int fgid, int dir,
Real const* flux, const int* flo, const int* fhi, int nc,
Real const* area, const int* ifd, Real scaling_factor)
{
Box bx;
bx = Box(IntVect(flo), IntVect(fhi));
bx.shiftHalf(dir,-1);
const FArrayBox fab(bx,nc,const_cast<Real*>(flux));
const FArrayBox areafab(bx,1,const_cast<Real*>(area));
flux_reg->add(fgid, dir, fab, areafab, ifd, scaling_factor);
}

void amrex_fi_flash_fluxregister_communicate (FlashFluxRegister* flux_reg)
{
flux_reg->communicate();
Expand Down
Loading

0 comments on commit d77d93b

Please sign in to comment.