Skip to content

Commit

Permalink
fix checkpointing to use old base_state size and allow multiple plotf… (
Browse files Browse the repository at this point in the history 
#1910)

* fix checkpointing to use old base_state size and allow multiple plotfile types

* remove unused

* allow us to write two different types of plotfiles in same run; fix restart with new base state ncomp/nghost

* remove unused and fix ncomp

* fix unused

* comment unused

* fix restart
  • Loading branch information
@asalmgren
asalmgren authored Oct 25, 2024
1 parent cee2934 commit 65276a8
Show file tree
Hide file tree
Showing 10 changed files with 213 additions and 96 deletions.
59 changes: 39 additions & 20 deletions Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,11 @@ using namespace amrex;
* @param[in] domain simulation domain
*/

void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain)
void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
int ncomp, const IntVect& nghost)
{
BL_PROFILE_VAR("impose_lateral_base_bcs()",impose_lateral_base_bcs);

int icomp = 0;
int ncomp = 3;

const int* bxlo = bx.loVect();
const int* bxhi = bx.hiVect();
const int* dlo = domain.loVect();
Expand Down Expand Up @@ -63,16 +61,16 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
if (!is_periodic_in_x)
{
// Populate ghost cells on lo-x and hi-x domain boundaries
Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,1);
if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,1);
if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,1);
if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,1);
Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-nghost[0]);
Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+nghost[0]);
if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,nghost[2]);
if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,nghost[2]);
if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,nghost[2]);
if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,nghost[2]);
ParallelFor(
bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = icomp+n;
int dest_comp = n;
int l_bc_type = bc_ptr[n].lo(0);
int iflip = dom_lo.x - 1 - i;
if (l_bc_type == ERFBCType::foextrap) {
Expand All @@ -85,7 +83,7 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
},
bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = icomp+n;
int dest_comp = n;
int h_bc_type = bc_ptr[n].hi(0);
int iflip = 2*dom_hi.x + 1 - i;
if (h_bc_type == ERFBCType::foextrap) {
Expand All @@ -102,16 +100,16 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
if (!is_periodic_in_y)
{
// Populate ghost cells on lo-y and hi-y domain boundaries
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,1);
if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,1);
if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,1);
if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,1);
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-nghost[1]);
Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+nghost[1]);
if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,nghost[2]);
if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,nghost[2]);
if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,nghost[2]);
if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,nghost[2]);
ParallelFor(
bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = icomp+n;
int dest_comp = n;
int l_bc_type = bc_ptr[n].lo(1);
int jflip = dom_lo.y - 1 - j;
if (l_bc_type == ERFBCType::foextrap) {
Expand All @@ -125,7 +123,7 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
},
bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
int dest_comp = icomp+n;
int dest_comp = n;
int h_bc_type = bc_ptr[n].hi(1);
int jflip = 2*dom_hi.y + 1 - j;
if (h_bc_type == ERFBCType::foextrap) {
Expand All @@ -141,3 +139,24 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4<Real>& dest
Gpu::streamSynchronize();
}

void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
int ncomp, const IntVect& /*nghost*/)
{
BL_PROFILE_VAR("impose_vertical_base_bcs()",impose_vertical_base_bcs);

const auto& dom_lo = lbound(domain);
const auto& dom_hi = ubound(domain);

Box bx_zlo(bx); bx_zlo.setBig(2,dom_lo.z-2);
Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+2);
ParallelFor(
bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
dest_arr(i,j,k,n) = dest_arr(i,j,dom_lo.z-1,n);
},
bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
dest_arr(i,j,k,n) = dest_arr(i,j,dom_hi.z+1,n);
}
);
}
9 changes: 6 additions & 3 deletions Source/BoundaryConditions/ERF_PhysBCFunct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -279,11 +279,14 @@ public:
* @param[in] nghost number of ghost cells to be filled for conserved variables
* @param[in] time time at which the data should be filled
*/
void operator() (amrex::MultiFab& mf, int icomp, int ncomp,
amrex::IntVect const& nghost, const amrex::Real time, int bccomp_cons);
void operator() (amrex::MultiFab& mf, int icomp, int ncomp, amrex::IntVect const& nghost);

void impose_lateral_basestate_bcs (const amrex::Array4<amrex::Real>& dest_arr,
const amrex::Box& bx, const amrex::Box& domain);
const amrex::Box& bx, const amrex::Box& domain,
int ncomp, const amrex::IntVect& nghost);
void impose_vertical_basestate_bcs (const amrex::Array4<amrex::Real>& dest_arr,
const amrex::Box& bx, const amrex::Box& domain,
int ncomp, const amrex::IntVect& nghost);

private:
int m_lev;
Expand Down
8 changes: 4 additions & 4 deletions Source/BoundaryConditions/ERF_PhysBCFunct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,8 +366,7 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,
} // OpenMP
} // operator()

void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
IntVect const& nghost, const Real /*time*/, int /*bccomp*/)
void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int ncomp, IntVect const& nghost)
{
BL_PROFILE("ERFPhysBCFunct_base::()");

Expand Down Expand Up @@ -403,9 +402,10 @@ void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/

if (!gdomain.contains(cbx2))
{
const Array4<Real> cons_arr = mf.array(mfi);
const Array4<Real> base_arr = mf.array(mfi);

impose_lateral_basestate_bcs(cons_arr,cbx1,domain);
impose_lateral_basestate_bcs(base_arr,cbx1,domain,ncomp,nghost);
impose_vertical_basestate_bcs(base_arr,cbx2,domain,ncomp,nghost);
}

} // MFIter
Expand Down
15 changes: 12 additions & 3 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,14 @@ AMREX_ENUM(InitType,
None, Input_Sounding, Ideal, Real, Metgrid, Uniform
);

/**
* Enum of possible plotfile types
*/
AMREX_ENUM(PlotFileType,
None, Amrex, Netcdf, HDF5
);


/**
* Enum of possible coarse/fine interpolation options
*/
Expand Down Expand Up @@ -318,7 +326,7 @@ public:
void derive_upwp (amrex::Vector<amrex::Real>& h_havg);

// write plotfile to disk
void WritePlotFile (int which, amrex::Vector<std::string> plot_var_names);
void WritePlotFile (int which, PlotFileType plotfile_type, amrex::Vector<std::string> plot_var_names);

void WriteMultiLevelPlotfileWithTerrain (const std::string &plotfilename,
int nlevels,
Expand Down Expand Up @@ -964,8 +972,9 @@ private:
static int pert_interval;
static amrex::Real sum_per;

// Native or NetCDF
static std::string plotfile_type;
// Write in native AMReX or NetCDF format for each plotfile
static PlotFileType plotfile_type_1;
static PlotFileType plotfile_type_2;

static InitType init_type;

Expand Down
65 changes: 47 additions & 18 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,8 @@ Real ERF::sum_per = -1.0;
int ERF::pert_interval = -1;

// Native AMReX vs NetCDF
std::string ERF::plotfile_type = "amrex";
PlotFileType ERF::plotfile_type_1 = PlotFileType::None;
PlotFileType ERF::plotfile_type_2 = PlotFileType::None;

InitType ERF::init_type;

Expand Down Expand Up @@ -369,11 +370,11 @@ ERF::Evolve ()

if (writeNow(cur_time, dt[0], step+1, m_plot_int_1, m_plot_per_1)) {
last_plot_file_step_1 = step+1;
WritePlotFile(1,plot_var_names_1);
WritePlotFile(1,plotfile_type_1,plot_var_names_1);
}
if (writeNow(cur_time, dt[0], step+1, m_plot_int_2, m_plot_per_2)) {
last_plot_file_step_2 = step+1;
WritePlotFile(2,plot_var_names_2);
WritePlotFile(2,plotfile_type_2,plot_var_names_2);
}

if (writeNow(cur_time, dt[0], step+1, m_check_int, m_check_per)) {
Expand Down Expand Up @@ -401,10 +402,10 @@ ERF::Evolve ()

// Write plotfiles at final time
if ( (m_plot_int_1 > 0 || m_plot_per_1 > 0.) && istep[0] > last_plot_file_step_1 ) {
WritePlotFile(1,plot_var_names_1);
WritePlotFile(1,plotfile_type_1,plot_var_names_1);
}
if ( (m_plot_int_2 > 0 || m_plot_per_2 > 0.) && istep[0] > last_plot_file_step_2) {
WritePlotFile(2,plot_var_names_2);
WritePlotFile(2,plotfile_type_1,plot_var_names_2);
}

if ( (m_check_int > 0 || m_check_per > 0.) && istep[0] > last_check_file_step) {
Expand Down Expand Up @@ -721,8 +722,10 @@ ERF::InitData_post ()
restart();

// Create the physbc objects for {cons, u, v, w, base state}
// We fill the additional base state ghost cells just in case we have read the old format
for (int lev(0); lev <= max_level; ++lev) {
make_physbcs(lev);
(*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect());
}
}

Expand Down Expand Up @@ -1082,12 +1085,12 @@ ERF::InitData_post ()
{
if (m_plot_int_1 > 0 || m_plot_per_1 > 0.)
{
WritePlotFile(1,plot_var_names_1);
WritePlotFile(1,plotfile_type_1,plot_var_names_1);
last_plot_file_step_1 = istep[0];
}
if (m_plot_int_2 > 0 || m_plot_per_2 > 0.)
{
WritePlotFile(2,plot_var_names_2);
WritePlotFile(2,plotfile_type_2,plot_var_names_2);
last_plot_file_step_2 = istep[0];
}
}
Expand Down Expand Up @@ -1475,8 +1478,42 @@ ERF::ReadParameters ()
// Flag to trigger initialization from input_sounding like WRF's ideal.exe
pp.query("init_sounding_ideal", init_sounding_ideal);

// Output format
pp.query("plotfile_type", plotfile_type);
PlotFileType plotfile_type_temp = PlotFileType::None;
pp.query_enum_case_insensitive("plotfile_type" ,plotfile_type_temp);
pp.query_enum_case_insensitive("plotfile_type_1",plotfile_type_1);
pp.query_enum_case_insensitive("plotfile_type_2",plotfile_type_2);
//
// This option is for backward consistency -- if only plotfile_type is set,
// then it will be used for both 1 and 2 if and only if they are not set
//
// Default is native amrex if no type is specified
//
if (plotfile_type_temp == PlotFileType::None) {
if (plotfile_type_1 == PlotFileType::None) {
plotfile_type_1 = PlotFileType::Amrex;
}
if (plotfile_type_2 == PlotFileType::None) {
plotfile_type_2 = PlotFileType::Amrex;
}
} else {
if (plotfile_type_1 == PlotFileType::None) {
plotfile_type_1 = plotfile_type_temp;
} else {
amrex::Abort("You must set either plotfile_type or plotfile_type_1, not both");
}
if (plotfile_type_2 == PlotFileType::None) {
plotfile_type_2 = plotfile_type_temp;
} else {
amrex::Abort("You must set either plotfile_type or plotfile_type_2, not both");
}
}
#ifndef ERF_USE_NETCDF
if (plotfile_type_1 == PlotFileType::Netcdf ||
plotfile_type_2 == PlotFileType::Netcdf) {
amrex::Abort("Plotfile type = Netcdf is not allowed without USE_NETCDF = TRUE");
}
#endif

pp.query("plot_file_1", plot_file_1);
pp.query("plot_file_2", plot_file_2);
pp.query("plot_int_1" , m_plot_int_1);
Expand Down Expand Up @@ -1582,7 +1619,7 @@ ERF::ParameterSanityChecks ()
{
AMREX_ALWAYS_ASSERT(cfl > 0. || fixed_dt[0] > 0.);

// We don't allow use_real_bcs to be true if init_type is not either InitType::Rreal or InitType::Metgrid
// We don't allow use_real_bcs to be true if init_type is not either InitType::Real or InitType::Metgrid
AMREX_ALWAYS_ASSERT(!use_real_bcs || ((init_type == InitType::Real) || (init_type == InitType::Metgrid)) );

AMREX_ALWAYS_ASSERT(real_width >= 0);
Expand All @@ -1602,14 +1639,6 @@ ERF::ParameterSanityChecks ()
}
}

if (plotfile_type != "amrex" &&
plotfile_type != "netcdf" && plotfile_type != "NetCDF" &&
plotfile_type != "hdf5" && plotfile_type != "HDF5" )
{
Print() << "User selected plotfile_type = " << plotfile_type << std::endl;
Abort("Dont know this plotfile_type");
}

// If fixed_mri_dt_ratio is set, it must be even
if (fixed_mri_dt_ratio > 0 && (fixed_mri_dt_ratio%2 != 0) )
{
Expand Down
5 changes: 1 addition & 4 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -346,9 +346,6 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
// from previous (pre-regrid) base_state array
// ********************************************************************************************
if (lev > 0) {
// Interp all three components: rho, p, pi
int icomp = 0; int bccomp = 0; int ncomp = 3;

Interpolater* mapper = &cell_cons_interp;

Vector<MultiFab*> fmf = {&base_state[lev ], &base_state[lev ]};
Expand All @@ -364,7 +361,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
BCVars::base_bc);

// Impose bc's outside the domain
(*physbcs_base[lev])(temp_base_state,icomp,ncomp,base_state[lev].nGrowVect(),time,bccomp);
(*physbcs_base[lev])(temp_base_state,0,temp_base_state.nComp(),base_state[lev].nGrowVect());

std::swap(temp_base_state, base_state[lev]);
}
Expand Down
46 changes: 39 additions & 7 deletions Source/IO/ERF_Checkpoint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,14 +131,23 @@ ERF::WriteCheckpointFile () const
VisMF::Write(zvel, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace"));

// Note that we write the ghost cells of the base state (unlike above)
IntVect ng = base_state[lev].nGrowVect();
MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ng);
MultiFab::Copy(base,base_state[lev],0,0,base.nComp(),ng);
// For backward compatibility we only write the first components and 1 ghost cell
IntVect ng; int ncomp;
bool write_old_base_state = true;
if (write_old_base_state) {
ng = IntVect{1};
ncomp = 3;
} else {
ng = base_state[lev].nGrowVect();
ncomp = base_state[lev].nComp();
}
MultiFab base(grids[lev],dmap[lev],ncomp,ng);
MultiFab::Copy(base,base_state[lev],0,0,ncomp,ng);
VisMF::Write(base, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "BaseState"));

if (solverChoice.use_terrain) {
// Note that we also write the ghost cells of z_phys_nd
ng = z_phys_nd[lev]->nGrowVect();
IntVect ng = z_phys_nd[lev]->nGrowVect();
MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng);
MultiFab::Copy(z_height,*z_phys_nd[lev],0,0,1,ng);
VisMF::Write(z_height, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Z_Phys_nd"));
Expand Down Expand Up @@ -379,10 +388,33 @@ ERF::ReadCheckpointFile ()
vars_new[lev][Vars::zvel].setBndry(1.0e34);

// Note that we read the ghost cells of the base state (unlike above)
IntVect ng = base_state[lev].nGrowVect();
MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ng);

// The original base state only had 3 components and 1 ghost cell -- we read this
// here to be consistent with the old style
IntVect ng; int ncomp;
bool read_old_base_state = true;
if (read_old_base_state) {
ng = IntVect{1};
ncomp = 3;
} else {
ng = base_state[lev].nGrowVect();
ncomp = base_state[lev].nComp();
}
MultiFab base(grids[lev],dmap[lev],ncomp,ng);
VisMF::Read(base, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "BaseState"));
MultiFab::Copy(base_state[lev],base,0,0,base.nComp(),ng);
MultiFab::Copy(base_state[lev],base,0,0,ncomp,ng);
if (read_old_base_state) {
for (MFIter mfi(base_state[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.growntilebox(1);
Array4<Real> const& fab = base_state[lev].array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
fab(i,j,k,BaseState::th0_comp) = getRhoThetagivenP(fab(i,j,k,BaseState::p0_comp))
/ fab(i,j,k,BaseState::r0_comp);
});
}
}
base_state[lev].FillBoundary(geom[lev].periodicity());

if (solverChoice.use_terrain) {
Expand Down
Loading

0 comments on commit 65276a8

Please sign in to comment.