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

fix checkpointing to use old base_state size and allow multiple plotf… #1910

Merged
merged 7 commits into from
Oct 25, 2024
Merged

fix checkpointing to use old base_state size and allow multiple plotf… #1910

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
89ab11b
fix checkpointing to use old base_state size and allow multiple plotf…
asalmgren Oct 25, 2024
a67b4ce
remove unused
asalmgren Oct 25, 2024
81fe88d
allow us to write two different types of plotfiles in same run; fix r…
asalmgren Oct 25, 2024
38f2fb2
remove unused and fix ncomp
asalmgren Oct 25, 2024
4517936
fix unused
asalmgren Oct 25, 2024
dce97f7
comment unused
asalmgren Oct 25, 2024
5bf6ed5
fix restart
asalmgren Oct 25, 2024
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
59 changes: 39 additions & 20 deletions Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading