diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp index e7a47ad52..867489d9a 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp @@ -11,13 +11,11 @@ using namespace amrex; * @param[in] domain simulation domain */ -void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest_arr, const Box& bx, const Box& domain) +void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& 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(); @@ -63,16 +61,16 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& 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) { @@ -85,7 +83,7 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& 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) { @@ -102,16 +100,16 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& 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) { @@ -125,7 +123,7 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& 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) { @@ -141,3 +139,27 @@ void ERFPhysBCFunct_base::impose_lateral_basestate_bcs (const Array4& dest Gpu::streamSynchronize(); } +void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4& 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 int* bxlo = bx.loVect(); + const int* bxhi = bx.hiVect(); + + const auto& dom_lo = lbound(domain); + const auto& dom_hi = ubound(domain); + + Box bx_zlo(bx); bx_zlo.setSmall(2,dom_lo.z-nghost[2]); bx_zlo.setBig(2,dom_lo.z-2); + Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+ 2); bx_zlo.setBig(2,dom_hi.z+nghost[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); + } + ); +} diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index d74a9834e..3666ea06c 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -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& 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& dest_arr, + const amrex::Box& bx, const amrex::Box& domain, + int ncomp, const amrex::IntVect& nghost); private: int m_lev; diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index 5ab47fb88..095c0c9cf 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -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::()"); @@ -403,9 +402,10 @@ void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/ if (!gdomain.contains(cbx2)) { - const Array4 cons_arr = mf.array(mfi); + const Array4 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,cbx1,domain,ncomp,nghost); } } // MFIter diff --git a/Source/ERF.H b/Source/ERF.H index 6c0532a9d..b2184c884 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -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 */ @@ -318,7 +326,7 @@ public: void derive_upwp (amrex::Vector& h_havg); // write plotfile to disk - void WritePlotFile (int which, amrex::Vector plot_var_names); + void WritePlotFile (int which, PlotFileType plotfile_type, amrex::Vector plot_var_names); void WriteMultiLevelPlotfileWithTerrain (const std::string &plotfilename, int nlevels, @@ -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; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index f28eaf737..b61ef8cba 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -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; @@ -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)) { @@ -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) { @@ -1082,12 +1083,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]; } } @@ -1475,8 +1476,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); @@ -1582,7 +1617,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); @@ -1602,14 +1637,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) ) { diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 3ed7ad42e..94cac0ee9 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -364,7 +364,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,icomp,ncomp,base_state[lev].nGrowVect()); std::swap(temp_base_state, base_state[lev]); } diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index 7a3452a51..55b1ab710 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -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")); @@ -379,10 +388,34 @@ 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 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); + }); + } + } + (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); base_state[lev].FillBoundary(geom[lev].periodicity()); if (solverChoice.use_terrain) { diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 7b6635b4d..00350a865 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -183,8 +183,14 @@ ERF::PlotFileVarNames (Vector plot_var_names ) // Write plotfile to disk void -ERF::WritePlotFile (int which, Vector plot_var_names) +ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector plot_var_names) { + if (plotfile_type == PlotFileType::Amrex) { + amrex::Print() << "WRITE AMREX " << std::endl; + } else if (plotfile_type == PlotFileType::Netcdf) { + amrex::Print() << "WRITE NETCDF " << std::endl; + } + const Vector varnames = PlotFileVarNames(plot_var_names); const int ncomp_mf = varnames.size(); @@ -349,15 +355,17 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp += 1; } - MultiFab r_hse(base_state[lev], make_alias, BaseState::r0_comp, 1); - MultiFab p_hse(base_state[lev], make_alias, BaseState::p0_comp, 1); + MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component + MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component if (containerHasElement(plot_var_names, "pres_hse")) { + // p_0 is second component of base_state MultiFab::Copy(mf[lev],p_hse,0,mf_comp,1,0); mf_comp += 1; } if (containerHasElement(plot_var_names, "dens_hse")) { + // r_0 is first component of base_state MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0); mf_comp += 1; } @@ -529,7 +537,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); if (solverChoice.anelastic[lev] == 1) { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - p_arr(i,j,k) = hse_arr(i,j,k,BaseState::p0_comp); + p_arr(i,j,k) = hse_arr(i,j,k,1); }); } else { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -623,7 +631,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); if (solverChoice.anelastic[lev] == 1) { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - p_arr(i,j,k) = hse_arr(i,j,k,BaseState::p0_comp); + p_arr(i,j,k) = hse_arr(i,j,k,1); }); } else { ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -1336,10 +1344,11 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } std::string plotfilename; - if (which == 1) + if (which == 1) { plotfilename = Concatenate(plot_file_1, istep[0], 5); - else if (which == 2) + } else if (which == 2) { plotfilename = Concatenate(plot_file_2, istep[0], 5); + } // LSM writes it's own data if (which==1 && plot_lsm) { @@ -1357,7 +1366,8 @@ ERF::WritePlotFile (int which, Vector plot_var_names) // Single level if (finest_level == 0) { - if (plotfile_type == "amrex") { + if (plotfile_type == PlotFileType::Amrex) + { Print() << "Writing native plotfile " << plotfilename << "\n"; if (solverChoice.use_terrain) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, @@ -1377,7 +1387,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) particleData.writePlotFile(plotfilename); #endif #ifdef ERF_USE_HDF5 - } else if (plotfile_type == "hdf5" || plotfile_type == "HDF5") { + } else if (plotfile_type == PlotFileType::HDF5) { Print() << "Writing plotfile " << plotfilename+"d01.h5" << "\n"; WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1, GetVecOfConstPtrs(mf), @@ -1385,19 +1395,19 @@ ERF::WritePlotFile (int which, Vector plot_var_names) Geom(), t_new[0], istep, refRatio()); #endif #ifdef ERF_USE_NETCDF - } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { + } else if (plotfile_type == PlotFileType::Netcdf) { int lev = 0; int l_which = 0; writeNCPlotFile(lev, l_which, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]); #endif } else { - Print() << "User specified plot_filetype = " << plotfile_type << std::endl; - Abort("Dont know this plot_filetype"); + // Here we assume the plotfile_type is PlotFileType::None + Print() << "Writing no plotfile since plotfile_type is none" << std::endl; } } else { // Multilevel - if (plotfile_type == "amrex") { + if (plotfile_type == PlotFileType::Amrex) { int lev0 = 0; int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]); @@ -1490,7 +1500,7 @@ ERF::WritePlotFile (int which, Vector plot_var_names) #endif #ifdef ERF_USE_NETCDF - } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { + } else if (plotfile_type == PlotFileType::Netcdf) { for (int lev = 0; lev <= finest_level; ++lev) { for (int which_box = 0; which_box < num_boxes_at_level[lev]; which_box++) { writeNCPlotFile(lev, which_box, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]); diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 745d8b469..9db03b244 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -30,7 +30,7 @@ ERF::initHSE (int lev) bool all_boxes_touch_bottom = true; Box domain(geom[lev].Domain()); - int icomp = 0; int bccomp = BCVars::base_bc; int ncomp = BaseState::num_comps; + int icomp = 0; int ncomp = BaseState::num_comps; Real time = 0.; @@ -64,7 +64,7 @@ ERF::initHSE (int lev) // 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. - (*physbcs_base[lev])(base_state[lev],icomp,ncomp,base_state[lev].nGrowVect(),time,bccomp); + (*physbcs_base[lev])(base_state[lev],icomp,ncomp,base_state[lev].nGrowVect()); } if (all_boxes_touch_bottom || lev > 0) { @@ -125,7 +125,7 @@ ERF::initHSE (int lev) // Impose physical bc's on the base state -- the values outside the fine region // but inside the domain have already been filled in the call above to InterpFromCoarseLevel // - (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect(),time,bccomp); + (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); base_state[lev].FillBoundary(geom[lev].periodicity()); } @@ -193,6 +193,9 @@ ERF::erf_enforce_hse (int lev, const Real rdOcp = solverChoice.rdOcp; + IntVect ngvect = base_state[lev].nGrowVect(); + int ng_z = ngvect[2]; + ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { // Set value at surface from Newton iteration for rho @@ -213,6 +216,11 @@ ERF::erf_enforce_hse (int lev, pres_arr(i,j,klo-1) = p_0 + hz * rho_arr(i,j,klo) * l_gravity; pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp); th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1); + for (int kk = 2; kk <= ng_z; kk++) { + pres_arr(i,j,klo-kk) = pres_arr(i,j,klo-1); + pi_arr(i,j,klo-kk) = pi_arr(i,j,klo-1); + th_arr(i,j,klo-kk) = th_arr(i,j,klo-1); + } } else { // If level > 0 and klo > 0, we need to use the value of pres_arr(i,j,klo-1) which was @@ -228,10 +236,12 @@ ERF::erf_enforce_hse (int lev, pres_arr(i,j,klo) = pres_arr(i,j,klo-1) - dz_loc * dens_interp * l_gravity; pi_arr(i,j,klo ) = getExnergivenP(pres_arr(i,j,klo ), rdOcp); - pi_arr(i,j,klo-1) = getExnergivenP(pres_arr(i,j,klo-1), rdOcp); - th_arr(i,j,klo ) = getRhoThetagivenP(pres_arr(i,j,klo )) / rho_arr(i,j,klo ); - th_arr(i,j,klo-1) = getRhoThetagivenP(pres_arr(i,j,klo-1)) / rho_arr(i,j,klo-1); + + for (int kk = 1; kk <= ng_z; kk++) { + pi_arr(i,j,klo-kk) = getExnergivenP(pres_arr(i,j,klo-kk), rdOcp); + th_arr(i,j,klo-kk) = getRhoThetagivenP(pres_arr(i,j,klo-kk)) / rho_arr(i,j,klo-kk); + } } @@ -256,12 +266,13 @@ ERF::erf_enforce_hse (int lev, int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0); int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1); - + if (pres[mfi].box().smallEnd(0) < domlo_x) { Box bx = mfi.nodaltilebox(2); + int ng = ngvect[0]; bx.setSmall(0,domlo_x-1); - bx.setBig(0,domlo_x-1); + bx.setBig(0,domlo_x-ng_z); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(domlo_x,j,k); pi_arr(i,j,k) = pi_arr(domlo_x,j,k); @@ -272,8 +283,9 @@ ERF::erf_enforce_hse (int lev, if (pres[mfi].box().bigEnd(0) > domhi_x) { Box bx = mfi.nodaltilebox(2); + int ng = ngvect[0]; bx.setSmall(0,domhi_x+1); - bx.setBig(0,domhi_x+1); + bx.setBig(0,domhi_x+ng); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(domhi_x,j,k); pi_arr(i,j,k) = pi_arr(domhi_x,j,k); @@ -284,7 +296,8 @@ ERF::erf_enforce_hse (int lev, if (pres[mfi].box().smallEnd(1) < domlo_y) { Box bx = mfi.nodaltilebox(2); - bx.setSmall(1,domlo_y-1); + int ng = ngvect[1]; + bx.setSmall(1,domlo_y-ng); bx.setBig(1,domlo_y-1); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(i,domlo_y,k); @@ -296,8 +309,9 @@ ERF::erf_enforce_hse (int lev, if (pres[mfi].box().bigEnd(1) > domhi_y) { Box bx = mfi.nodaltilebox(2); + int ng = ngvect[1]; bx.setSmall(1,domhi_y+1); - bx.setBig(1,domhi_y+1); + bx.setBig(1,domhi_y+ng); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { pres_arr(i,j,k) = pres_arr(i,domhi_y,k); pi_arr(i,j,k) = pi_arr(i,domhi_y,k); diff --git a/Source/Initialization/ERF_init_from_input_sounding.cpp b/Source/Initialization/ERF_init_from_input_sounding.cpp index ec31635e5..e05e97043 100644 --- a/Source/Initialization/ERF_init_from_input_sounding.cpp +++ b/Source/Initialization/ERF_init_from_input_sounding.cpp @@ -88,9 +88,7 @@ ERF::init_from_input_sounding (int lev) // 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; - Real dummy_time = 0.; - (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect(),dummy_time,bccomp); + (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); } auto& lev_new = vars_new[lev]; @@ -279,20 +277,28 @@ init_bx_scalars_from_input_sounding_hse (const Box &bx, pi_hse_arr(i,j,k) = getExnergivenRTh(rhoTh_k, l_rdOcp); th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k)) / r_hse_arr(i,j,k); + // TODO: we should be setting this to the number of ghost cells of base_state[lev] + // instead of hard-wiring it here! + int ng = 3; + // FOEXTRAP hse arrays 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); - pi_hse_arr(i, j, k-1) = pi_hse_arr(i,j,k); - th_hse_arr(i, j, k-1) = th_hse_arr(i,j,k); + for (int kk = 1; kk <= ng; kk++) { + r_hse_arr(i, j, k-kk) = r_hse_arr(i,j,k); + p_hse_arr(i, j, k-kk) = p_hse_arr(i,j,k); + pi_hse_arr(i, j, k-kk) = pi_hse_arr(i,j,k); + th_hse_arr(i, j, k-kk) = th_hse_arr(i,j,k); + } } else if (k==ktop) { - 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); - pi_hse_arr(i, j, k+1) = pi_hse_arr(i,j,k); - th_hse_arr(i, j, k+1) = th_hse_arr(i,j,k); + for (int kk = 1; kk <= ng; kk++) { + r_hse_arr(i, j, k+kk) = r_hse_arr(i,j,k); + p_hse_arr(i, j, k+kk) = p_hse_arr(i,j,k); + pi_hse_arr(i, j, k+kk) = pi_hse_arr(i,j,k); + th_hse_arr(i, j, k+kk) = th_hse_arr(i,j,k); + } } // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there