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

move the fourth order center <-> average stuff to C++ #884

Merged
merged 11 commits into from
May 6, 2020
Merged

move the fourth order center <-> average stuff to C++ #884

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
85423e5
move the fourth order center <-> average stuff to C++
zingale Apr 27, 2020
3ab2be0
more work on porting the average stuff to C++
zingale Apr 27, 2020
557114e
this compiles now
zingale Apr 27, 2020
969b046
fix assert
zingale Apr 27, 2020
afce97b
fix the fix
zingale Apr 27, 2020
4322ff6
fix a copy for SDC -- wrong number of ghost cells were assumed
zingale Apr 27, 2020
d8c6daa
Merge branch 'sdc_source_fix' into cxx_fourth
zingale Apr 27, 2020
2415b8c
fix a box
zingale Apr 27, 2020
fd60515
fix compilation
zingale Apr 27, 2020
4de317c
Merge branch 'development' into cxx_fourth
zingale May 2, 2020
859a796
address some PR reviews
zingale May 6, 2020
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
73 changes: 44 additions & 29 deletions Source/driver/Castro.cpp
View file
Open in desktop
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 @@ -3499,22 +3509,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 @@ -3605,19 +3617,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,32 +195,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
Open in desktop
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
Open in desktop
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,
maxpkatz marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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