Skip to content

Commit

Permalink
move the fourth order center <-> average stuff to C++ (#884)
Browse files Browse the repository at this point in the history 
All of the averaging routines are now available in C++
  • Loading branch information
@zingale
zingale authored May 6, 2020
1 parent 324682d commit a85c95b
Show file tree
Hide file tree
Showing 11 changed files with 515 additions and 561 deletions.
73 changes: 44 additions & 29 deletions Source/driver/Castro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1109,16 +1109,20 @@ Castro::initData ()
AmrLevel::FillPatch(*this, Sborder, NUM_GROW, cur_time, State_Type, 0, NUM_STATE);

// note: this cannot be tiled
const int* domain_lo = geom.Domain().loVect();
const int* domain_hi = geom.Domain().hiVect();
auto domain_lo = geom.Domain().loVect3d();
auto domain_hi = geom.Domain().hiVect3d();

FArrayBox tmp;

for (MFIter mfi(S_new); mfi.isValid(); ++mfi)
{
const Box& box = mfi.validbox();
const Box& box = mfi.validbox();

tmp.resize(box, 1);
Elixir elix_tmp = tmp.elixir();
auto tmp_arr = tmp.array();

ca_make_fourth_in_place(BL_TO_FORTRAN_BOX(box),
BL_TO_FORTRAN_FAB(Sborder[mfi]),
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
make_fourth_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi);
}

// now copy back the averages
Expand All @@ -1132,16 +1136,20 @@ Castro::initData ()
AmrLevel::FillPatch(*this, Sborder, NUM_GROW, cur_time, State_Type, 0, NUM_STATE);

// convert to centers -- not tile safe
const int* domain_lo = geom.Domain().loVect();
const int* domain_hi = geom.Domain().hiVect();
auto domain_lo = geom.Domain().loVect3d();
auto domain_hi = geom.Domain().hiVect3d();

FArrayBox tmp;

for (MFIter mfi(S_new); mfi.isValid(); ++mfi)
{
const Box& box = mfi.growntilebox(2);

ca_make_cell_center_in_place(BL_TO_FORTRAN_BOX(box),
BL_TO_FORTRAN_FAB(Sborder[mfi]),
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
tmp.resize(box, 1);
Elixir elix_tmp = tmp.elixir();
auto tmp_arr = tmp.array();

make_cell_center_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi);
}

// reset the energy -- do this in one ghost cell so we can average in place below
Expand Down Expand Up @@ -1184,9 +1192,11 @@ Castro::initData ()
{
const Box& box = mfi.validbox();

ca_make_fourth_in_place(BL_TO_FORTRAN_BOX(box),
BL_TO_FORTRAN_FAB(Sborder[mfi]),
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
tmp.resize(box, 1);
Elixir elix_tmp = tmp.elixir();
auto tmp_arr = tmp.array();

make_fourth_in_place(box, Sborder.array(mfi), tmp_arr, domain_lo, domain_hi);
}

// now copy back the averages for UEINT and UTEMP only
Expand Down Expand Up @@ -3472,22 +3482,24 @@ Castro::computeTemp(

// convert to cell centers -- this will result in Stemp being
// cell centered only on 1 ghost cells
const int* domain_lo = geom.Domain().loVect();
const int* domain_hi = geom.Domain().hiVect();
auto domain_lo = geom.Domain().loVect3d();
auto domain_hi = geom.Domain().hiVect3d();

FArrayBox tmp;

for (MFIter mfi(Stemp); mfi.isValid(); ++mfi) {
const Box& bx = mfi.growntilebox(1);
const Box& bx0 = mfi.tilebox();
const int idx = mfi.tileIndex();

ca_compute_lap_term(BL_TO_FORTRAN_BOX(bx0),
BL_TO_FORTRAN_FAB(Stemp[mfi]),
BL_TO_FORTRAN_ANYD(Eint_lap[mfi]), UEINT,
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
compute_lap_term(bx0, Stemp.array(mfi), Eint_lap.array(mfi), UEINT,
domain_lo, domain_hi);

tmp.resize(bx, 1);
Elixir elix_tmp = tmp.elixir();
auto tmp_arr = tmp.array();

ca_make_cell_center_in_place(BL_TO_FORTRAN_BOX(bx),
BL_TO_FORTRAN_FAB(Stemp[mfi]),
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
make_cell_center_in_place(bx, Stemp.array(mfi), tmp_arr, domain_lo, domain_hi);

}

Expand Down Expand Up @@ -3578,19 +3590,22 @@ Castro::computeTemp(
// cell-averages -- this is 4th-order and will be a no-op for
// those zones where e wasn't changed.

const int* domain_lo = geom.Domain().loVect();
const int* domain_hi = geom.Domain().hiVect();
auto domain_lo = geom.Domain().loVect3d();
auto domain_hi = geom.Domain().hiVect3d();

FArrayBox tmp;

for (MFIter mfi(Stemp); mfi.isValid(); ++mfi) {

const Box& bx = mfi.tilebox();
const int idx = mfi.tileIndex();

// only temperature
ca_make_fourth_in_place_n(BL_TO_FORTRAN_BOX(bx),
BL_TO_FORTRAN_FAB(Stemp[mfi]), UTEMP,
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
tmp.resize(bx, 1);
Elixir elix_tmp = tmp.elixir();
auto tmp_arr = tmp.array();

// only temperature
make_fourth_in_place_n(bx, Stemp.array(mfi), UTEMP, tmp_arr, domain_lo, domain_hi);
}

// correct UEINT
Expand Down
26 changes: 0 additions & 26 deletions Source/driver/Castro_F.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,32 +185,6 @@ extern "C"
const amrex::Real* dx, amrex::Real* dt);
#endif

void ca_make_fourth_average(const int lo[], const int hi[],
BL_FORT_FAB_ARG_3D(q), const int* nc,
const BL_FORT_FAB_ARG_3D(q_bar), const int* nc_bar,
const int* domlo, const int* domhi);

void ca_make_fourth_in_place(const int lo[], const int hi[],
BL_FORT_FAB_ARG_3D(q), const int* nc,
const int* domlo, const int* domhi);

void ca_make_fourth_in_place_n(const int lo[], const int hi[],
BL_FORT_FAB_ARG_3D(q), const int* nc, const int ncomp,
const int* domlo, const int* domhi);

void ca_compute_lap_term(const int lo[], const int hi[],
const BL_FORT_FAB_ARG_3D(U), const int* nc,
BL_FORT_FAB_ARG_ANYD(lap), const int ncomp,
const int* domlo, const int* domhi);

void ca_make_cell_center(const int lo[], const int hi[],
const BL_FORT_FAB_ARG_3D(U), const int* nc,
BL_FORT_FAB_ARG_3D(U_cc), const int* nc_cc,
const int* domlo, const int* domhi);

void ca_make_cell_center_in_place(const int lo[], const int hi[],
BL_FORT_FAB_ARG_3D(q), const int* nc,
const int* domlo, const int* domhi);

#ifdef RADIATION
void ca_inelastic_sct
Expand Down
51 changes: 26 additions & 25 deletions Source/driver/Castro_advance_sdc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,8 @@ Castro::do_advance_sdc (Real time,
MultiFab& S_old = get_old_data(State_Type);
MultiFab& S_new = get_new_data(State_Type);

const int* domain_lo = geom.Domain().loVect();
const int* domain_hi = geom.Domain().hiVect();
auto domain_lo = geom.Domain().loVect3d();
auto domain_hi = geom.Domain().hiVect3d();

// Perform initialization steps.

Expand Down Expand Up @@ -104,10 +104,8 @@ Castro::do_advance_sdc (Real time,

for (MFIter mfi(S_new); mfi.isValid(); ++mfi) {
const Box& gbx = mfi.growntilebox(1);
ca_make_cell_center(BL_TO_FORTRAN_BOX(gbx),
BL_TO_FORTRAN_FAB(Sborder[mfi]),
BL_TO_FORTRAN_FAB(Sburn[mfi]),
AMREX_INT_ANYD(domain_lo), AMREX_INT_ANYD(domain_hi));

make_cell_center(gbx, Sborder.array(mfi), Sburn.array(mfi), domain_lo, domain_hi);

}

Expand All @@ -123,11 +121,16 @@ Castro::do_advance_sdc (Real time,
AmrLevel::FillPatch(*this, old_source, old_source.nGrow(), prev_time, Source_Type, 0, NSRC);

// Now convert to cell averages. This loop cannot be tiled.
FArrayBox tmp;

for (MFIter mfi(S_new); mfi.isValid(); ++mfi) {
const Box& bx = mfi.tilebox();
ca_make_fourth_in_place(BL_TO_FORTRAN_BOX(bx),
BL_TO_FORTRAN_FAB(old_source[mfi]),
AMREX_INT_ANYD(domain_lo), AMREX_INT_ANYD(domain_hi));

tmp.resize(bx, 1);
Elixir elix_tmp = tmp.elixir();
auto tmp_arr = tmp.array();

make_fourth_in_place(bx, old_source.array(mfi), tmp_arr, domain_lo, domain_hi);
}

} else {
Expand Down Expand Up @@ -306,6 +309,7 @@ Castro::do_advance_sdc (Real time,

FArrayBox U_center;
FArrayBox R_center;
FArrayBox tmp;

// this cannot be tiled
for (MFIter mfi(R_new); mfi.isValid(); ++mfi) {
Expand All @@ -314,33 +318,30 @@ Castro::do_advance_sdc (Real time,

if (sdc_order == 4) {

const int* domain_lo = geom.Domain().loVect();
const int* domain_hi = geom.Domain().hiVect();

// convert S_new to cell-centers
U_center.resize(obx, NUM_STATE);
ca_make_cell_center(BL_TO_FORTRAN_BOX(obx),
BL_TO_FORTRAN_FAB(Sborder[mfi]),
BL_TO_FORTRAN_FAB(U_center),
AMREX_INT_ANYD(domain_lo), AMREX_INT_ANYD(domain_hi));
Elixir elix_u_center = U_center.elixir();
auto const U_center_arr = U_center.array();

make_cell_center(obx, Sborder.array(mfi), U_center_arr, domain_lo, domain_hi);

// pass in the reaction source and state at centers, including one ghost cell
// and derive everything that is needed including 1 ghost cell
R_center.resize(obx, R_new.nComp());
Elixir elix_r_center = R_center.elixir();
auto const R_center_arr = R_center.array();

Array4<const Real> const Sburn_arr = Sburn.array(mfi);
Array4<const Real> const U_center_arr = U_center.array();
Array4<Real> const R_center_arr = R_center.array();

// we don't worry about the difference between centers and averages
ca_store_reaction_state(obx,
Sburn_arr,
U_center_arr,
R_center_arr);
ca_store_reaction_state(obx, Sburn_arr, U_center_arr, R_center_arr);

// convert R_new from centers to averages in place
ca_make_fourth_in_place(BL_TO_FORTRAN_BOX(bx),
BL_TO_FORTRAN_FAB(R_center),
AMREX_INT_ANYD(domain_lo), AMREX_INT_ANYD(domain_hi));
tmp.resize(bx, 1);
Elixir elix_tmp = tmp.elixir();
auto const tmp_arr = tmp.array();

make_fourth_in_place(bx, R_center_arr, tmp_arr, domain_lo, domain_hi);

// store
R_new[mfi].copy(R_center, bx, 0, bx, 0, R_new.nComp());
Expand Down
47 changes: 47 additions & 0 deletions Source/hydro/Castro_hydro.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -541,4 +541,51 @@
/// @param ng Number of ghost cells
///
void linear_to_hybrid_momentum(amrex::MultiFab& state, int ng = 0);

#endif

#ifdef TRUE_SDC
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Real compute_laplacian(int i, int j, int k, int n,
Array4<Real const> const a,
GpuArray<bool, AMREX_SPACEDIM> const& lo_periodic, GpuArray<bool, AMREX_SPACEDIM> const& hi_periodic,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Real trans_laplacian(int i, int j, int k, int n,
int idir,
Array4<Real const> const a,
GpuArray<bool, AMREX_SPACEDIM> const& lo_periodic, GpuArray<bool, AMREX_SPACEDIM> const& hi_periodic,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

void make_cell_center(const Box& bx,
Array4<Real const> const& U,
Array4<Real> const& U_cc,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

void make_cell_center_in_place(const Box& bx,
Array4<Real> const& U,
Array4<Real> const& tmp,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

void compute_lap_term(const Box& bx,
Array4<Real const> const& U,
Array4<Real> const& lap, const int ncomp,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

void make_fourth_average(const Box& bx,
Array4<Real> const& q,
Array4<Real const> const& q_bar,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

void make_fourth_in_place(const Box& bx,
Array4<Real> const& q,
Array4<Real> const& tmp,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

void make_fourth_in_place_n(const Box& bx,
Array4<Real> const& q, const int ncomp,
Array4<Real> const& tmp,
GpuArray<int, 3> const& domlo, GpuArray<int, 3> const& domhi);

#endif
26 changes: 10 additions & 16 deletions Source/hydro/Castro_hydro.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,11 +124,13 @@ Castro::cons_to_prim_fourth(const Real time)
// convert the conservative state cell averages to primitive cell
// averages with 4th order accuracy

const int* domain_lo = geom.Domain().loVect();
const int* domain_hi = geom.Domain().hiVect();
auto domain_lo = geom.Domain().loVect3d();
auto domain_hi = geom.Domain().hiVect3d();

MultiFab& S_new = get_new_data(State_Type);

FArrayBox U_cc;

// we don't support radiation here
#ifdef RADIATION
amrex::Abort("radiation not supported to fourth order");
Expand All @@ -147,13 +149,11 @@ Castro::cons_to_prim_fourth(const Real time)
// convert U_avg to U_cc -- this will use a Laplacian
// operation and will result in U_cc defined only on
// NUM_GROW-1 ghost cells at the end.
FArrayBox U_cc;
U_cc.resize(qbx, NUM_STATE);
Elixir elix_u_cc = U_cc.elixir();
auto const U_cc_arr = U_cc.array();

ca_make_cell_center(BL_TO_FORTRAN_BOX(qbxm1),
BL_TO_FORTRAN_FAB(Sborder[mfi]),
BL_TO_FORTRAN_FAB(U_cc),
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
make_cell_center(qbxm1, Sborder.array(mfi), U_cc_arr, domain_lo, domain_hi);

// enforce the minimum density on the new cell-centered state
ca_enforce_minimum_density
Expand Down Expand Up @@ -182,7 +182,6 @@ Castro::cons_to_prim_fourth(const Real time)
// convert U_cc to q_cc (we'll store this temporarily in q,
// qaux). This will remain valid only on the NUM_GROW-1 ghost
// cells.
auto U_cc_arr = U_cc.array();
auto q_arr = q.array(mfi);
auto qaux_arr = qaux.array(mfi);

Expand Down Expand Up @@ -218,17 +217,12 @@ Castro::cons_to_prim_fourth(const Real time)
// this will create q, qaux in NUM_GROW-1 ghost cells, but that's
// we need here

ca_make_fourth_average(BL_TO_FORTRAN_BOX(qbxm1),
BL_TO_FORTRAN_FAB(q[mfi]),
BL_TO_FORTRAN_FAB(q_bar[mfi]),
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));
make_fourth_average(qbxm1, q.array(mfi), q_bar.array(mfi), domain_lo, domain_hi);

// not sure if we need to convert qaux this way, or if we can
// just evaluate it (we may not need qaux at all actually)
ca_make_fourth_average(BL_TO_FORTRAN_BOX(qbxm1),
BL_TO_FORTRAN_FAB(qaux[mfi]),
BL_TO_FORTRAN_FAB(qaux_bar[mfi]),
AMREX_ARLIM_ANYD(domain_lo), AMREX_ARLIM_ANYD(domain_hi));

make_fourth_average(qbxm1, qaux.array(mfi), qaux_bar.array(mfi), domain_lo, domain_hi);

}

Expand Down
Loading

0 comments on commit a85c95b

Please sign in to comment.