diff --git a/Source/ERF.H b/Source/ERF.H index 5b526f7b8..2533b0ffe 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -315,8 +315,10 @@ public: const amrex::Vector &mf, const amrex::Vector &mf_nd, const amrex::Vector &varnames, + const amrex::Vector& my_geom, amrex::Real time, const amrex::Vector &level_steps, + const amrex::Vector& my_ref_ratio, const std::string &versionName = "HyperCLaw-V1.1", const std::string &levelPrefix = "Level_", const std::string &mfPrefix = "Cell", @@ -327,8 +329,10 @@ public: int nlevels, const amrex::Vector &bArray, const amrex::Vector &varnames, + const amrex::Vector& my_geom, amrex::Real time, const amrex::Vector &level_steps, + const amrex::Vector& my_ref_ratio, const std::string &versionName, const std::string &levelPrefix, const std::string &mfPrefix) const; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index ab1c85133..b433ab704 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -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."); diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index bc02bedfe..f3bd61633 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -1368,7 +1368,7 @@ ERF::WritePlotFile (int which, Vector 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), @@ -1403,7 +1403,11 @@ ERF::WritePlotFile (int which, Vector 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 r2(finest_level); Vector g2(finest_level+1); @@ -1416,15 +1420,15 @@ ERF::WritePlotFile (int which, Vector plot_var_names) // Define a new multi-level array of Geometry's so that we pass the new "domain" at lev > 0 Array 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); @@ -1449,16 +1453,16 @@ ERF::WritePlotFile (int which, Vector plot_var_names) // Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile Vector 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, @@ -1471,7 +1475,7 @@ ERF::WritePlotFile (int which, Vector 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, @@ -1502,8 +1506,10 @@ ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nl const Vector& mf, const Vector& mf_nd, const Vector& varnames, + const Vector& geom, Real time, const Vector& level_steps, + const Vector& ref_ratio, const std::string &versionName, const std::string &levelPrefix, const std::string &mfPrefix, @@ -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); }; @@ -1587,14 +1593,16 @@ ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile, int nlevels, const Vector &bArray, const Vector &varnames, - Real time, - const Vector &level_steps, + const Vector& my_geom, + Real my_time, + const Vector& level_steps, + const Vector& 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); @@ -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) { @@ -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'; } diff --git a/Source/Initialization/ERF_init_from_input_sounding.cpp b/Source/Initialization/ERF_init_from_input_sounding.cpp index e0a07954f..abbf40565 100644 --- a/Source/Initialization/ERF_init_from_input_sounding.cpp +++ b/Source/Initialization/ERF_init_from_input_sounding.cpp @@ -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]; @@ -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 @@ -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);