Skip to content

Commit

Permalink
fix for multilevel base state with init_type = input_sounding (erf-mo…
Browse files Browse the repository at this point in the history 
…del#1868)
  • Loading branch information
@asalmgren
asalmgren authored Oct 9, 2024
1 parent 7957c46 commit 1c78d9e
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 28 deletions.
4 changes: 4 additions & 0 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,8 +315,10 @@ public:
const amrex::Vector<const amrex::MultiFab*> &mf,
const amrex::Vector<const amrex::MultiFab*> &mf_nd,
const amrex::Vector<std::string> &varnames,
const amrex::Vector<amrex::Geometry>& my_geom,
amrex::Real time,
const amrex::Vector<int> &level_steps,
const amrex::Vector<amrex::IntVect>& my_ref_ratio,
const std::string &versionName = "HyperCLaw-V1.1",
const std::string &levelPrefix = "Level_",
const std::string &mfPrefix = "Cell",
Expand All @@ -327,8 +329,10 @@ public:
int nlevels,
const amrex::Vector<amrex::BoxArray> &bArray,
const amrex::Vector<std::string> &varnames,
const amrex::Vector<amrex::Geometry>& my_geom,
amrex::Real time,
const amrex::Vector<int> &level_steps,
const amrex::Vector<amrex::IntVect>& my_ref_ratio,
const std::string &versionName,
const std::string &levelPrefix,
const std::string &mfPrefix) const;
Expand Down
6 changes: 5 additions & 1 deletion Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1301,11 +1301,15 @@ ERF::init_only (int lev, Real time)
// The base state is initialized by integrating vertically through the
// input sounding, if the init_sounding_ideal flag is set; otherwise
// it is set by initHSE()
init_from_input_sounding(lev);

// The physbc's need the terrain but are needed for initHSE
// We have already made the terrain in the call to init_zphys
// in MakeNewLevelFromScratch
make_physbcs(lev);

// Now init the base state and the data itself
init_from_input_sounding(lev);

if (init_sounding_ideal) {
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(solverChoice.use_gravity,
"Gravity should be on to be consistent with sounding initialization.");
Expand Down
58 changes: 33 additions & 25 deletions Source/IO/ERF_Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1368,7 +1368,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
GetVecOfConstPtrs(mf),
GetVecOfConstPtrs(mf_nd),
varnames,
t_new[0], istep);
Geom(), t_new[0], istep, refRatio());
} else {
WriteMultiLevelPlotfile(plotfilename, finest_level+1,
GetVecOfConstPtrs(mf),
Expand Down Expand Up @@ -1403,7 +1403,11 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)

if (plotfile_type == "amrex") {

if (ref_ratio[0][2] == 1) {
// NOTE: this assumes that -- at least -- there is refinement in the x-direction
// Otherwise the logic will be wrong
int lev0 = 0;
int desired_ratio = ref_ratio[lev0][0];
if (ref_ratio[lev0][2] == 1) {

Vector<IntVect> r2(finest_level);
Vector<Geometry> g2(finest_level+1);
Expand All @@ -1416,15 +1420,15 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)

// Define a new multi-level array of Geometry's so that we pass the new "domain" at lev > 0
Array<int,AMREX_SPACEDIM> periodicity =
{Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)};
g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data());
{Geom()[lev0].isPeriodic(0),Geom()[lev0].isPeriodic(1),Geom()[lev0].isPeriodic(2)};
g2[lev0].define(Geom()[lev0].Domain(),&(Geom()[lev0].ProbDomain()),0,periodicity.data());

r2[0] = IntVect(1,1,ref_ratio[0][0]);
r2[0] = IntVect(1,desired_ratio, desired_ratio);
for (int lev = 1; lev <= finest_level; ++lev) {
if (lev > 1) {
r2[lev-1][0] = 1;
r2[lev-1][1] = 1;
r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0];
r2[lev-1][1] *= desired_ratio / ref_ratio[lev0][1];
r2[lev-1][2] *= desired_ratio / ref_ratio[lev0][2];
}

mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
Expand All @@ -1449,16 +1453,16 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
// Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile
Vector<IntVect> rr(finest_level);
for (int lev = 0; lev < finest_level; ++lev) {
rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]);
rr[lev] = IntVect(desired_ratio);
}

Print() << "Writing plotfile " << plotfilename << "\n";
if (solverChoice.use_terrain) {
WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1,
GetVecOfConstPtrs(mf),
GetVecOfConstPtrs(mf2),
GetVecOfConstPtrs(mf_nd),
varnames,
t_new[0], istep);
g2, t_new[0], istep, rr);
} else {
WriteMultiLevelPlotfile(plotfilename, finest_level+1,
GetVecOfConstPtrs(mf2), varnames,
Expand All @@ -1471,7 +1475,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
GetVecOfConstPtrs(mf),
GetVecOfConstPtrs(mf_nd),
varnames,
t_new[0], istep);
geom, t_new[0], istep, ref_ratio);
} else {
WriteMultiLevelPlotfile(plotfilename, finest_level+1,
GetVecOfConstPtrs(mf), varnames,
Expand Down Expand Up @@ -1502,8 +1506,10 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl
const Vector<const MultiFab*>& mf,
const Vector<const MultiFab*>& mf_nd,
const Vector<std::string>& varnames,
const Vector<Geometry>& geom,
Real time,
const Vector<int>& level_steps,
const Vector<IntVect>& ref_ratio,
const std::string &versionName,
const std::string &levelPrefix,
const std::string &mfPrefix,
Expand Down Expand Up @@ -1542,7 +1548,7 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl
std::ofstream::binary);
if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName);
WriteGenericPlotfileHeaderWithTerrain(HeaderFile, nlevels, boxArrays, varnames,
time, level_steps, versionName,
geom, time, level_steps, ref_ratio, versionName,
levelPrefix, mfPrefix);
};

Expand Down Expand Up @@ -1587,14 +1593,16 @@ ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile,
int nlevels,
const Vector<BoxArray> &bArray,
const Vector<std::string> &varnames,
Real time,
const Vector<int> &level_steps,
const Vector<Geometry>& my_geom,
Real my_time,
const Vector<int>& level_steps,
const Vector<IntVect>& my_ref_ratio,
const std::string &versionName,
const std::string &levelPrefix,
const std::string &mfPrefix) const
{
AMREX_ALWAYS_ASSERT(nlevels <= bArray.size());
AMREX_ALWAYS_ASSERT(nlevels <= ref_ratio.size()+1);
AMREX_ALWAYS_ASSERT(nlevels <= my_ref_ratio.size()+1);
AMREX_ALWAYS_ASSERT(nlevels <= level_steps.size());

HeaderFile.precision(17);
Expand All @@ -1608,22 +1616,22 @@ ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile,
HeaderFile << varnames[ivar] << "\n";
}
HeaderFile << AMREX_SPACEDIM << '\n';
HeaderFile << time << '\n';
HeaderFile << my_time << '\n';
HeaderFile << finest_level << '\n';
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
HeaderFile << geom[0].ProbLo(i) << ' ';
HeaderFile << my_geom[0].ProbLo(i) << ' ';
}
HeaderFile << '\n';
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
HeaderFile << geom[0].ProbHi(i) << ' ';
HeaderFile << my_geom[0].ProbHi(i) << ' ';
}
HeaderFile << '\n';
for (int i = 0; i < finest_level; ++i) {
HeaderFile << ref_ratio[i][0] << ' ';
HeaderFile << my_ref_ratio[i][0] << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
HeaderFile << geom[i].Domain() << ' ';
HeaderFile << my_geom[i].Domain() << ' ';
}
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
Expand All @@ -1632,25 +1640,25 @@ ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile,
HeaderFile << '\n';
for (int i = 0; i <= finest_level; ++i) {
for (int k = 0; k < AMREX_SPACEDIM; ++k) {
HeaderFile << geom[i].CellSize()[k] << ' ';
HeaderFile << my_geom[i].CellSize()[k] << ' ';
}
HeaderFile << '\n';
}
HeaderFile << (int) geom[0].Coord() << '\n';
HeaderFile << (int) my_geom[0].Coord() << '\n';
HeaderFile << "0\n";

for (int level = 0; level <= finest_level; ++level) {
HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n';
HeaderFile << level << ' ' << bArray[level].size() << ' ' << my_time << '\n';
HeaderFile << level_steps[level] << '\n';

const IntVect& domain_lo = geom[level].Domain().smallEnd();
const IntVect& domain_lo = my_geom[level].Domain().smallEnd();
for (int i = 0; i < bArray[level].size(); ++i)
{
// Need to shift because the RealBox ctor we call takes the
// physical location of index (0,0,0). This does not affect
// the usual cases where the domain index starts with 0.
const Box& b = shift(bArray[level][i], -domain_lo);
RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo());
RealBox loc = RealBox(b, my_geom[level].CellSize(), my_geom[level].ProbLo());
for (int n = 0; n < AMREX_SPACEDIM; ++n) {
HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n';
}
Expand Down
31 changes: 29 additions & 2 deletions Source/Initialization/ERF_init_from_input_sounding.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,31 @@ ERF::init_from_input_sounding (int lev)
// this will calculate the hydrostatically balanced density and pressure
// profiles following WRF ideal.exe
if (init_sounding_ideal) input_sounding_data.calc_rho_p(0);

} else {
//
// We need to do this interp from coarse level in order to set the values of
// the base state inside the domain but outside of the fine region
//
base_state[lev-1].FillBoundary(geom[lev-1].periodicity());
//
// NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab
// have been pre-filled - this includes ghost cells both inside and outside
// the domain
//
InterpFromCoarseLevel(base_state[lev], base_state[lev].nGrowVect(),
IntVect(0,0,0), // do not fill ghost cells outside the domain
base_state[lev-1], 0, 0, 3,
geom[lev-1], geom[lev],
refRatio(lev-1), &cell_cons_interp,
domain_bcs_type, BCVars::base_bc);

// We need to do this here because the interpolation above may leave corners unfilled
// when the corners need to be filled by, for example, reflection of the fine ghost
// cell outside the fine region but inide the domain.
int bccomp = BCVars::base_bc;
int icomp = 0; int ncomp = 3; Real dummy_time = 0.;
(*physbcs_base[lev])(base_state[lev],icomp,ncomp,base_state[lev].nGrowVect(),dummy_time,bccomp);
}

auto& lev_new = vars_new[lev];
Expand Down Expand Up @@ -211,12 +236,14 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx,
const int inp_sound_size = inputSoundingData.size(0);

// Geometry
int ktop = bx.bigEnd(2);
const Real* prob_lo = geomdata.ProbLo();
const Real* dx = geomdata.CellSize();
const Real z_lo = prob_lo[2];
const Real dz = dx[2];

int kbot = geomdata.Domain().smallEnd(2);
int ktop = geomdata.Domain().bigEnd(2);

// We want to set the lateral BC values, too
Box gbx = bx; // Copy constructor
gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
Expand Down Expand Up @@ -247,7 +274,7 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx,
pi_hse_arr(i, j, k) = getExnergivenRTh(rhoTh_k, l_rdOcp);

// FOEXTRAP hse arrays
if (k==0)
if (k==kbot)
{
r_hse_arr (i, j, k-1) = r_hse_arr (i,j,k);
p_hse_arr (i, j, k-1) = p_hse_arr (i,j,k);
Expand Down

0 comments on commit 1c78d9e

Please sign in to comment.