diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 78d438d47..a2ef237cb 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -203,10 +203,10 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Microphysics/Kessler/ERF_Init_Kessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_Kessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_Update_Kessler.cpp - ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp - ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp - ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp - ${SRC_DIR}/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp + ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp + ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp + ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp + ${SRC_DIR}/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp ${SRC_DIR}/LandSurfaceModel/SLM/ERF_SLM.cpp ${SRC_DIR}/LandSurfaceModel/MM5/ERF_MM5.cpp ) diff --git a/Exec/AWAKEN/CMakeLists.txt b/Exec/WindFarmTests/AWAKEN/CMakeLists.txt similarity index 100% rename from Exec/AWAKEN/CMakeLists.txt rename to Exec/WindFarmTests/AWAKEN/CMakeLists.txt diff --git a/Exec/AWAKEN/ERF_prob.H b/Exec/WindFarmTests/AWAKEN/ERF_prob.H similarity index 100% rename from Exec/AWAKEN/ERF_prob.H rename to Exec/WindFarmTests/AWAKEN/ERF_prob.H diff --git a/Exec/AWAKEN/ERF_prob.cpp b/Exec/WindFarmTests/AWAKEN/ERF_prob.cpp similarity index 100% rename from Exec/AWAKEN/ERF_prob.cpp rename to Exec/WindFarmTests/AWAKEN/ERF_prob.cpp diff --git a/Exec/AWAKEN/GNUmakefile b/Exec/WindFarmTests/AWAKEN/GNUmakefile similarity index 100% rename from Exec/AWAKEN/GNUmakefile rename to Exec/WindFarmTests/AWAKEN/GNUmakefile diff --git a/Exec/AWAKEN/Make.package b/Exec/WindFarmTests/AWAKEN/Make.package similarity index 100% rename from Exec/AWAKEN/Make.package rename to Exec/WindFarmTests/AWAKEN/Make.package diff --git a/Exec/AWAKEN/README b/Exec/WindFarmTests/AWAKEN/README similarity index 100% rename from Exec/AWAKEN/README rename to Exec/WindFarmTests/AWAKEN/README diff --git a/Exec/AWAKEN/inputs_All_EWP b/Exec/WindFarmTests/AWAKEN/inputs_All_EWP similarity index 100% rename from Exec/AWAKEN/inputs_All_EWP rename to Exec/WindFarmTests/AWAKEN/inputs_All_EWP diff --git a/Exec/AWAKEN/inputs_KingPlains_EWP b/Exec/WindFarmTests/AWAKEN/inputs_KingPlains_EWP similarity index 100% rename from Exec/AWAKEN/inputs_KingPlains_EWP rename to Exec/WindFarmTests/AWAKEN/inputs_KingPlains_EWP diff --git a/Exec/AWAKEN/inputs_KingPlains_Fitch b/Exec/WindFarmTests/AWAKEN/inputs_KingPlains_Fitch similarity index 100% rename from Exec/AWAKEN/inputs_KingPlains_Fitch rename to Exec/WindFarmTests/AWAKEN/inputs_KingPlains_Fitch diff --git a/Exec/AWAKEN/inputs_KingPlains_SimpleAD b/Exec/WindFarmTests/AWAKEN/inputs_KingPlains_SimpleAD similarity index 100% rename from Exec/AWAKEN/inputs_KingPlains_SimpleAD rename to Exec/WindFarmTests/AWAKEN/inputs_KingPlains_SimpleAD diff --git a/Exec/AWAKEN/windturbines_All.txt b/Exec/WindFarmTests/AWAKEN/windturbines_All.txt similarity index 100% rename from Exec/AWAKEN/windturbines_All.txt rename to Exec/WindFarmTests/AWAKEN/windturbines_All.txt diff --git a/Exec/AWAKEN/windturbines_ArmadilloFlats.txt b/Exec/WindFarmTests/AWAKEN/windturbines_ArmadilloFlats.txt similarity index 100% rename from Exec/AWAKEN/windturbines_ArmadilloFlats.txt rename to Exec/WindFarmTests/AWAKEN/windturbines_ArmadilloFlats.txt diff --git a/Exec/AWAKEN/windturbines_KingPlains.txt b/Exec/WindFarmTests/AWAKEN/windturbines_KingPlains.txt similarity index 100% rename from Exec/AWAKEN/windturbines_KingPlains.txt rename to Exec/WindFarmTests/AWAKEN/windturbines_KingPlains.txt diff --git a/Exec/AWAKEN/windturbines_spec_AWAKEN.tbl b/Exec/WindFarmTests/AWAKEN/windturbines_spec_AWAKEN.tbl similarity index 100% rename from Exec/AWAKEN/windturbines_spec_AWAKEN.tbl rename to Exec/WindFarmTests/AWAKEN/windturbines_spec_AWAKEN.tbl diff --git a/Exec/EWP/CMakeLists.txt b/Exec/WindFarmTests/EWP/CMakeLists.txt similarity index 100% rename from Exec/EWP/CMakeLists.txt rename to Exec/WindFarmTests/EWP/CMakeLists.txt diff --git a/Exec/EWP/ERF_prob.H b/Exec/WindFarmTests/EWP/ERF_prob.H similarity index 100% rename from Exec/EWP/ERF_prob.H rename to Exec/WindFarmTests/EWP/ERF_prob.H diff --git a/Exec/EWP/ERF_prob.cpp b/Exec/WindFarmTests/EWP/ERF_prob.cpp similarity index 100% rename from Exec/EWP/ERF_prob.cpp rename to Exec/WindFarmTests/EWP/ERF_prob.cpp diff --git a/Exec/EWP/GNUmakefile b/Exec/WindFarmTests/EWP/GNUmakefile similarity index 100% rename from Exec/EWP/GNUmakefile rename to Exec/WindFarmTests/EWP/GNUmakefile diff --git a/Exec/EWP/Make.package b/Exec/WindFarmTests/EWP/Make.package similarity index 100% rename from Exec/EWP/Make.package rename to Exec/WindFarmTests/EWP/Make.package diff --git a/Exec/EWP/README b/Exec/WindFarmTests/EWP/README similarity index 100% rename from Exec/EWP/README rename to Exec/WindFarmTests/EWP/README diff --git a/Exec/EWP/inputs_1WT_lat_lon b/Exec/WindFarmTests/EWP/inputs_1WT_lat_lon similarity index 100% rename from Exec/EWP/inputs_1WT_lat_lon rename to Exec/WindFarmTests/EWP/inputs_1WT_lat_lon diff --git a/Exec/EWP/inputs_1WT_x_y b/Exec/WindFarmTests/EWP/inputs_1WT_x_y similarity index 100% rename from Exec/EWP/inputs_1WT_x_y rename to Exec/WindFarmTests/EWP/inputs_1WT_x_y diff --git a/Exec/EWP/inputs_WindFarm_lat_lon b/Exec/WindFarmTests/EWP/inputs_WindFarm_lat_lon similarity index 100% rename from Exec/EWP/inputs_WindFarm_lat_lon rename to Exec/WindFarmTests/EWP/inputs_WindFarm_lat_lon diff --git a/Exec/EWP/inputs_WindFarm_x_y b/Exec/WindFarmTests/EWP/inputs_WindFarm_x_y similarity index 100% rename from Exec/EWP/inputs_WindFarm_x_y rename to Exec/WindFarmTests/EWP/inputs_WindFarm_x_y diff --git a/Exec/EWP/windturbines_loc_lat_lon_1WT.txt b/Exec/WindFarmTests/EWP/windturbines_loc_lat_lon_1WT.txt similarity index 100% rename from Exec/EWP/windturbines_loc_lat_lon_1WT.txt rename to Exec/WindFarmTests/EWP/windturbines_loc_lat_lon_1WT.txt diff --git a/Exec/EWP/windturbines_loc_lat_lon_WindFarm.txt b/Exec/WindFarmTests/EWP/windturbines_loc_lat_lon_WindFarm.txt similarity index 100% rename from Exec/EWP/windturbines_loc_lat_lon_WindFarm.txt rename to Exec/WindFarmTests/EWP/windturbines_loc_lat_lon_WindFarm.txt diff --git a/Exec/EWP/windturbines_loc_x_y_1WT.txt b/Exec/WindFarmTests/EWP/windturbines_loc_x_y_1WT.txt similarity index 100% rename from Exec/EWP/windturbines_loc_x_y_1WT.txt rename to Exec/WindFarmTests/EWP/windturbines_loc_x_y_1WT.txt diff --git a/Exec/EWP/windturbines_loc_x_y_WindFarm.txt b/Exec/WindFarmTests/EWP/windturbines_loc_x_y_WindFarm.txt similarity index 100% rename from Exec/EWP/windturbines_loc_x_y_WindFarm.txt rename to Exec/WindFarmTests/EWP/windturbines_loc_x_y_WindFarm.txt diff --git a/Exec/EWP/windturbines_spec_1WT.tbl b/Exec/WindFarmTests/EWP/windturbines_spec_1WT.tbl similarity index 100% rename from Exec/EWP/windturbines_spec_1WT.tbl rename to Exec/WindFarmTests/EWP/windturbines_spec_1WT.tbl diff --git a/Exec/EWP/windturbines_spec_WindFarm.tbl b/Exec/WindFarmTests/EWP/windturbines_spec_WindFarm.tbl similarity index 100% rename from Exec/EWP/windturbines_spec_WindFarm.tbl rename to Exec/WindFarmTests/EWP/windturbines_spec_WindFarm.tbl diff --git a/Exec/GeneralActuatorDisk/ERF_prob.cpp b/Exec/WindFarmTests/GeneralActuatorDisk/ERF_prob.cpp similarity index 100% rename from Exec/GeneralActuatorDisk/ERF_prob.cpp rename to Exec/WindFarmTests/GeneralActuatorDisk/ERF_prob.cpp diff --git a/Exec/SimpleActuatorDisk/CMakeLists.txt b/Exec/WindFarmTests/SimpleActuatorDisk/CMakeLists.txt similarity index 100% rename from Exec/SimpleActuatorDisk/CMakeLists.txt rename to Exec/WindFarmTests/SimpleActuatorDisk/CMakeLists.txt diff --git a/Exec/SimpleActuatorDisk/ERF_prob.H b/Exec/WindFarmTests/SimpleActuatorDisk/ERF_prob.H similarity index 100% rename from Exec/SimpleActuatorDisk/ERF_prob.H rename to Exec/WindFarmTests/SimpleActuatorDisk/ERF_prob.H diff --git a/Exec/SimpleActuatorDisk/ERF_prob.cpp b/Exec/WindFarmTests/SimpleActuatorDisk/ERF_prob.cpp similarity index 100% rename from Exec/SimpleActuatorDisk/ERF_prob.cpp rename to Exec/WindFarmTests/SimpleActuatorDisk/ERF_prob.cpp diff --git a/Exec/SimpleActuatorDisk/GNUmakefile b/Exec/WindFarmTests/SimpleActuatorDisk/GNUmakefile similarity index 100% rename from Exec/SimpleActuatorDisk/GNUmakefile rename to Exec/WindFarmTests/SimpleActuatorDisk/GNUmakefile diff --git a/Exec/SimpleActuatorDisk/Make.package b/Exec/WindFarmTests/SimpleActuatorDisk/Make.package similarity index 100% rename from Exec/SimpleActuatorDisk/Make.package rename to Exec/WindFarmTests/SimpleActuatorDisk/Make.package diff --git a/Exec/SimpleActuatorDisk/README b/Exec/WindFarmTests/SimpleActuatorDisk/README similarity index 100% rename from Exec/SimpleActuatorDisk/README rename to Exec/WindFarmTests/SimpleActuatorDisk/README diff --git a/Exec/SimpleActuatorDisk/inputs_1WT_lat_lon b/Exec/WindFarmTests/SimpleActuatorDisk/inputs_1WT_lat_lon similarity index 100% rename from Exec/SimpleActuatorDisk/inputs_1WT_lat_lon rename to Exec/WindFarmTests/SimpleActuatorDisk/inputs_1WT_lat_lon diff --git a/Exec/SimpleActuatorDisk/inputs_1WT_x_y b/Exec/WindFarmTests/SimpleActuatorDisk/inputs_1WT_x_y similarity index 100% rename from Exec/SimpleActuatorDisk/inputs_1WT_x_y rename to Exec/WindFarmTests/SimpleActuatorDisk/inputs_1WT_x_y diff --git a/Exec/SimpleActuatorDisk/inputs_WindFarm_lat_lon b/Exec/WindFarmTests/SimpleActuatorDisk/inputs_WindFarm_lat_lon similarity index 100% rename from Exec/SimpleActuatorDisk/inputs_WindFarm_lat_lon rename to Exec/WindFarmTests/SimpleActuatorDisk/inputs_WindFarm_lat_lon diff --git a/Exec/SimpleActuatorDisk/inputs_WindFarm_x_y b/Exec/WindFarmTests/SimpleActuatorDisk/inputs_WindFarm_x_y similarity index 100% rename from Exec/SimpleActuatorDisk/inputs_WindFarm_x_y rename to Exec/WindFarmTests/SimpleActuatorDisk/inputs_WindFarm_x_y diff --git a/Exec/SimpleActuatorDisk/windturbines_loc_lat_lon_1WT.txt b/Exec/WindFarmTests/SimpleActuatorDisk/windturbines_loc_lat_lon_1WT.txt similarity index 100% rename from Exec/SimpleActuatorDisk/windturbines_loc_lat_lon_1WT.txt rename to Exec/WindFarmTests/SimpleActuatorDisk/windturbines_loc_lat_lon_1WT.txt diff --git a/Exec/SimpleActuatorDisk/windturbines_loc_x_y_1WT.txt b/Exec/WindFarmTests/SimpleActuatorDisk/windturbines_loc_x_y_1WT.txt similarity index 100% rename from Exec/SimpleActuatorDisk/windturbines_loc_x_y_1WT.txt rename to Exec/WindFarmTests/SimpleActuatorDisk/windturbines_loc_x_y_1WT.txt diff --git a/Exec/SimpleActuatorDisk/windturbines_spec_1WT.tbl b/Exec/WindFarmTests/SimpleActuatorDisk/windturbines_spec_1WT.tbl similarity index 100% rename from Exec/SimpleActuatorDisk/windturbines_spec_1WT.tbl rename to Exec/WindFarmTests/SimpleActuatorDisk/windturbines_spec_1WT.tbl diff --git a/Exec/inputs.mesh_refinement b/Exec/inputs.mesh_refinement deleted file mode 100644 index 55da0226c..000000000 --- a/Exec/inputs.mesh_refinement +++ /dev/null @@ -1,15 +0,0 @@ - -# This file is meant to serve as an example of how to specify parameters -# that are only relevant when running with mesh refinement - -# REFINEMENT / REGRIDDING -amr.max_level = 1 # maximum level number allowed -amr.ref_ratio_vect = 2 2 1 3 3 1 # refinement ratio written as lev0_x lev0_y lev0_z lev1_x lev1_y lev1_z ... -amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est - -erf.regrid_int = 2 2 2 2 # how often to regrid - -# How to specify a region for static refinement -erf.refinement_indicators = box1 -erf.box1.in_box_lo = .15 .25 -erf.box1.in_box_hi = .55 .65 diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp index d2101ad4c..de97493ae 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp @@ -259,7 +259,7 @@ void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4& des { dest_arr(i,j,k,BaseState::r0_comp) = dest_arr(i,j,dom_lo.z,BaseState::r0_comp); dest_arr(i,j,k,BaseState::p0_comp) = p_0 - - dest_arr(i,j,k,BaseState::r0_comp) * hz * 9.81; + dest_arr(i,j,k,BaseState::r0_comp) * hz * CONST_GRAV; dest_arr(i,j,k,BaseState::pi0_comp) = dest_arr(i,j,dom_lo.z,BaseState::pi0_comp); dest_arr(i,j,k,BaseState::th0_comp) = dest_arr(i,j,dom_lo.z,BaseState::th0_comp); } diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp index 1a2edbd1c..5932ee22b 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp @@ -261,7 +261,7 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, const Array4& z_phys_nd, const GpuArray dxInv, - int icomp, int ncomp) + int icomp, int ncomp, bool do_terrain_adjustment) { BL_PROFILE_VAR("impose_vertical_cons_bcs()",impose_vertical_cons_bcs); const auto& dom_lo = lbound(domain); @@ -416,7 +416,7 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr ); } - if (m_z_phys_nd) { + if (do_terrain_adjustment && m_z_phys_nd) { const auto& bx_lo = lbound(bx); const auto& bx_hi = ubound(bx); const BCRec* bc_ptr_h = bcrs.data(); diff --git a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp index a494132f1..982040c74 100644 --- a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp +++ b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp @@ -61,11 +61,13 @@ ERF::FillIntermediatePatch (int lev, Real time, } } - // amrex::Print() << "CONS ONLY " << cons_only << std::endl; - // amrex::Print() << "ICOMP NCOMP " << icomp_cons << " " << ncomp_cons << " NGHOST " << ng_cons << std::endl; + // amrex::Print() << "LEVEL " << lev << " CONS ONLY " << cons_only << + // " ICOMP NCOMP " << icomp_cons << " " << ncomp_cons << " NGHOST " << ng_cons << std::endl; - AMREX_ALWAYS_ASSERT(mfs_mom.size() == IntVars::NumTypes); - AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes); + if (!cons_only) { + AMREX_ALWAYS_ASSERT(mfs_mom.size() == IntVars::NumTypes); + AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes); + } // Enforce no penetration for thin immersed body if (!cons_only) { @@ -122,7 +124,9 @@ ERF::FillIntermediatePatch (int lev, Real time, Vector ftime = {time,time}; // Impose physical bc's on fine data (note time and 0 are not used) - (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); + bool do_fb = true; bool do_terrain_adjustment = false; + (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc, + do_fb, do_terrain_adjustment); if ( (icomp_cons+ncomp_cons > 1) && (interpolation_type == StateInterpType::Perturbational) ) { diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index 756e4725d..6e201adb9 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -53,7 +53,7 @@ public: */ void operator() (amrex::MultiFab& mf, int icomp, int ncomp, amrex::IntVect const& nghost, const amrex::Real time, int bccomp_cons, - bool do_fb); + bool do_fb = true, bool do_terrain_adjustment = true); void impose_lateral_cons_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, @@ -62,7 +62,7 @@ public: const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& z_nd, const amrex::GpuArray dxInv, - int icomp, int ncomp); + int icomp, int ncomp, bool do_terrain_adjustment = true); private: int m_lev; diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index fb5a18f8f..49714b383 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -15,7 +15,7 @@ using namespace amrex; void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, IntVect const& nghost, const Real /*time*/, int /*bccomp*/, - bool do_fb) + bool do_fb, bool do_terrain_adjustment) { BL_PROFILE("ERFPhysBCFunct_cons::()"); @@ -32,21 +32,6 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, } } - MultiFab z_nd_mf_loc; - if (m_z_phys_nd) { - m_z_phys_nd->FillBoundary(m_geom.periodicity()); - BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList(); - for (auto& b : bl_z_phys) { - b.setSmall(2,0); - b.setBig(2,1); - } - BoxArray ba_z(std::move(bl_z_phys)); - - z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1],0)); - z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(), - z_nd_mf_loc.nGrowVect()); - } - // // We fill all of the interior and periodic ghost cells first, so we can fill // those directly inside the lateral and vertical calls. @@ -77,7 +62,7 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, if (m_z_phys_nd) { - z_nd_arr = z_nd_mf_loc.const_array(mfi); + z_nd_arr = m_z_phys_nd->const_array(mfi); } if (!gdomain.contains(cbx2)) @@ -91,7 +76,7 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, } // We send the full FAB box with ghost cells - impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp); + impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp,do_terrain_adjustment); } } // MFIter @@ -116,20 +101,6 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, } } - MultiFab z_nd_mf_loc; - if (m_z_phys_nd) { - m_z_phys_nd->FillBoundary(m_geom.periodicity()); - BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList(); - for (auto& b : bl_z_phys) { - b.setSmall(2,0); - b.setBig(2,1); - } - BoxArray ba_z(std::move(bl_z_phys)); - z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0]+1,nghost[1],0)); - z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(), - z_nd_mf_loc.nGrowVect()); - } - // // We fill all of the interior and periodic ghost cells first, so we can fill // those directly inside the lateral and vertical calls. @@ -163,7 +134,7 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, if (m_z_phys_nd) { - z_nd_arr = z_nd_mf_loc.const_array(mfi); + z_nd_arr = m_z_phys_nd->const_array(mfi); } if (!gdomainx.contains(xbx2)) @@ -202,20 +173,6 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, } } - MultiFab z_nd_mf_loc; - if (m_z_phys_nd) { - m_z_phys_nd->FillBoundary(m_geom.periodicity()); - BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList(); - for (auto& b : bl_z_phys) { - b.setSmall(2,0); - b.setBig(2,1); - } - BoxArray ba_z(std::move(bl_z_phys)); - z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1]+1,0)); - z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(), - z_nd_mf_loc.nGrowVect()); - } - // // We fill all of the interior and periodic ghost cells first, so we can fill // those directly inside the lateral and vertical calls. @@ -249,7 +206,7 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/, if (m_z_phys_nd) { - z_nd_arr = z_nd_mf_loc.const_array(mfi); + z_nd_arr = m_z_phys_nd->const_array(mfi); } if (!gdomainy.contains(ybx2)) @@ -293,21 +250,6 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, // if (gdomainz.smallEnd(2) == 0) gdomainz.setSmall(2,1); - Box ndomain = convert(domain,IntVect(1,1,1)); - - MultiFab z_nd_mf_loc; - if (m_z_phys_nd) { - BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList(); - for (auto& b : bl_z_phys) { - b &= ndomain; - } - BoxArray ba_z(std::move(bl_z_phys)); - z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1],0)); - z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(), - z_nd_mf_loc.nGrowVect()); - } - z_nd_mf_loc.FillBoundary(m_geom.periodicity()); - // // We fill all of the interior and periodic ghost cells first, so we can fill // those directly inside the lateral and vertical calls. @@ -339,7 +281,7 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, if (m_z_phys_nd) { - z_nd_arr = z_nd_mf_loc.const_array(mfi); + z_nd_arr = m_z_phys_nd->const_array(mfi); } // diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 377eb09be..f9d29ce2b 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -47,7 +47,9 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // Define dmap[lev] to be dm SetDistributionMap(lev, dm); - // amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl; + if (verbose) { + amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl; + } if (lev == 0) init_bcs(); @@ -179,7 +181,9 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, { AMREX_ALWAYS_ASSERT(lev > 0); - // amrex::Print() <<" NEW BA FROM COARSE AT LEVEL " << lev << " " << ba << std::endl; + if (verbose) { + amrex::Print() <<" NEW BA FROM COARSE AT LEVEL " << lev << " " << ba << std::endl; + } //******************************************************************************************** // This allocates all kinds of things, including but not limited to: solution arrays, @@ -283,7 +287,9 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, void ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapping& dm) { - // amrex::Print() <<" REMAKING WITH NEW BA AT LEVEL " << lev << " " << ba << std::endl; + if (verbose) { + amrex::Print() <<" REMAKING WITH NEW BA AT LEVEL " << lev << " " << ba << std::endl; + } AMREX_ALWAYS_ASSERT(lev > 0); AMREX_ALWAYS_ASSERT(solverChoice.terrain_type != TerrainType::Moving); @@ -291,7 +297,9 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp BoxArray ba_old(vars_new[lev][Vars::cons].boxArray()); DistributionMapping dm_old(vars_new[lev][Vars::cons].DistributionMap()); - // amrex::Print() <<" OLD BA AT LEVEL " << lev << " " << ba_old << std::endl; + if (verbose) { + amrex::Print() <<" OLD BA AT LEVEL " << lev << " " << ba_old << std::endl; + } int ncomp_cons = vars_new[lev][Vars::cons].nComp(); IntVect ngrow_state = vars_new[lev][Vars::cons].nGrowVect(); diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index c8b7294e9..100ed7458 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -1410,10 +1410,15 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p int lev0 = 0; int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]); - bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 0) || (ref_ratio[lev0][1] == 0) ) || - (ref_ratio[lev0][2] == 0) ); + bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 1) || (ref_ratio[lev0][1] == 1) ) || + (ref_ratio[lev0][2] == 1) ); + for (int lev = 1; lev < finest_level; lev++) { + any_ratio_one = any_ratio_one || + ( ( (ref_ratio[lev][0] == 1) || (ref_ratio[lev][1] == 1) ) || + (ref_ratio[lev][2] == 1) ); + } - if (any_ratio_one == 1 && m_expand_plotvars_to_unif_rr) + if (any_ratio_one && m_expand_plotvars_to_unif_rr) { Vector r2(finest_level); Vector g2(finest_level+1); @@ -1435,9 +1440,9 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p for (int lev = 1; lev <= finest_level; ++lev) { if (lev > 1) { - r2[lev-1][0] *= desired_ratio / ref_ratio[lev-1][0]; - r2[lev-1][1] *= desired_ratio / ref_ratio[lev-1][1]; - r2[lev-1][2] *= desired_ratio / ref_ratio[lev-1][2]; + r2[lev-1][0] = r2[lev-2][0] * desired_ratio / ref_ratio[lev-1][0]; + r2[lev-1][1] = r2[lev-2][1] * desired_ratio / ref_ratio[lev-1][1]; + r2[lev-1][2] = r2[lev-2][2] * desired_ratio / ref_ratio[lev-1][2]; } mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0); diff --git a/Source/TimeIntegration/ERF_MRI.H b/Source/TimeIntegration/ERF_MRI.H index 0a38cc63c..9e69858d7 100644 --- a/Source/TimeIntegration/ERF_MRI.H +++ b/Source/TimeIntegration/ERF_MRI.H @@ -56,11 +56,8 @@ private: int force_stage1_single_substep; /** - * \brief The pre_update function is called by the integrator on stage data before using it to evaluate a right-hand side. - * \brief The post_update function is called by the integrator on stage data at the end of the stage + * \brief The no_substep function is called when we have no acoustic substepping */ - std::function pre_update; - std::function post_update; std::function no_substep; @@ -164,16 +161,6 @@ public: return slow_fast_timestep_ratio; } - void set_pre_update (std::function F) - { - pre_update = F; - } - - void set_post_update (std::function F) - { - post_update = F; - } - void set_no_substep (std::function F) { no_substep = F; @@ -247,8 +234,6 @@ public: // RK3 for compressible integrator for (int nrk = 0; nrk < 3; nrk++) { - // amrex::Print() << "Starting RK3: Step " << nrk+1 << std::endl; - // Capture the time we got to in the previous RK step old_time_stage = time_stage; @@ -285,12 +270,6 @@ public: // step 2 starts with S_stage = S^* and we always start substepping at the old time // step 3 starts with S_stage = S^** and we always start substepping at the old time - // All pre_update does is call cons_to_prim, and we have done this with the old - // data already before starting the RK steps - if (nrk > 0) { - pre_update(S_new, S_new[IntVars::cons].nGrow()); - } - // S_scratch also holds the average momenta over the fast iterations -- // to be used to update the slow variables -- we will initialize with // the momenta used in the first call to the slow_rhs, then update @@ -327,10 +306,6 @@ public: // (because we did update the fast variables in the substepping) // **************************************************** slow_rhs_post(*F_slow, S_old, S_new, *S_sum, *S_scratch, time, old_time_stage, time_stage, nrk); - - // Call the post-update hook for S_new after all the fast steps completed - // This will update S_prim that is used in the slow RHS - post_update(S_new, time + nsubsteps*dtau, S_new[IntVars::cons].nGrow(), S_new[IntVars::xmom].nGrow()); } // nrk } else { @@ -343,12 +318,6 @@ public: if (nrk == 0) { nsubsteps = 1; dtau = timestep; time_stage = time + timestep; } if (nrk == 1) { nsubsteps = 1; dtau = timestep; time_stage = time + timestep; } - // All pre_update does is call cons_to_prim, and we have done this with the old - // data already before starting the RK steps - if (nrk > 0) { - pre_update(S_new, S_new[IntVars::cons].nGrow()); - } - // S_scratch also holds the average momenta over the fast iterations -- // to be used to update the slow variables -- we will initialize with // the momenta used in the first call to the slow_rhs, then update @@ -358,9 +327,14 @@ public: no_substep(*S_sum, S_old, *F_slow, time + nsubsteps*dtau, nsubsteps*dtau, nrk); + // **************************************************** + // Evaluate F_slow(S_stage) only for the slow variables + // Note that we are using the current stage versions (in S_new) of the slow variables + // (because we didn't update the slow variables in the substepping) + // but we are using the "new" versions (in S_sum) of the velocities + // (because we did update the fast variables in the substepping) + // **************************************************** slow_rhs_post(*F_slow, S_old, S_new, *S_sum, *S_scratch, time, old_time_stage, time_stage, nrk); - - post_update(S_new, time + nsubsteps*dtau, S_new[IntVars::cons].nGrow(), S_new[IntVars::xmom].nGrow()); } // nrk } diff --git a/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H b/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H index 1213b8cf6..bddb14930 100644 --- a/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H @@ -13,7 +13,9 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk, const Real new_substep_time) { BL_PROFILE("fast_rhs_fun"); - if (verbose) amrex::Print() << "Calling fast rhs at level " << level << " with dt = " << dtau << std::endl; + if (verbose) amrex::Print() << "Fast time integration at level " << level + << " from " << old_substep_time << " to " << new_substep_time + << " with dt = " << dtau << std::endl; // Define beta_s here so that it is consistent between where we make the fast coefficients // and where we use them diff --git a/Source/TimeIntegration/ERF_TI_no_substep_fun.H b/Source/TimeIntegration/ERF_TI_no_substep_fun.H index f49314a7c..1b7987799 100644 --- a/Source/TimeIntegration/ERF_TI_no_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_no_substep_fun.H @@ -16,6 +16,10 @@ const amrex::GpuArray scomp_fast = {0,0,0,0}; const amrex::GpuArray ncomp_fast = {2,1,1,1}; + if (verbose) amrex::Print() << " No-substepping time integration at level " << level + << " to " << time_for_fp + << " with dt = " << slow_dt << std::endl; + // Update S_sum = S_stage only for the fast variables #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -127,9 +131,8 @@ // Even if we update all the conserved variables we don't need // to fillpatch the slow ones every acoustic substep - int ng_cons = S_sum[IntVars::cons].nGrow(); - int ng_vel = S_sum[IntVars::xmom].nGrow(); - apply_bcs(S_sum, time_for_fp, ng_cons, ng_vel, fast_only=true, vel_and_mom_synced=false); + apply_bcs(S_sum, time_for_fp, S_sum[IntVars::cons].nGrow(), S_sum[IntVars::xmom].nGrow(), + fast_only=true, vel_and_mom_synced=false); #ifdef ERF_USE_POISSON_SOLVE if (solverChoice.anelastic[level]) { diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index b96c2e096..af489069c 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -12,6 +12,15 @@ const Real new_stage_time, const int nrk) { + // + // Define primitive variables for all later RK stages + // (We have already done this for the first RK step) + // + if (nrk > 0) { + int ng_cons = S_data[IntVars::cons].nGrow(); + cons_to_prim(S_data[IntVars::cons], ng_cons); + } + BL_PROFILE("slow_rhs_fun_pre"); if (verbose) Print() << "Making slow rhs at time " << old_stage_time << " for fast variables advancing from " << old_step_time << " to " << new_stage_time << std::endl; @@ -281,23 +290,6 @@ #endif }; // end slow_rhs_fun_pre - // ************************************************************* - // This called before RK stage - // ************************************************************* - auto pre_update_fun = [&](Vector& S_data, int ng_cons) - { - cons_to_prim(S_data[IntVars::cons], ng_cons); - }; - - // ************************************************************* - // This called after every RK stage -- from MRI or SRI - // ************************************************************* - auto post_update_fun = [&](Vector& S_data, - const Real time_for_fp, int ng_cons, int ng_vel) - { - apply_bcs(S_data, time_for_fp, ng_cons, ng_vel, fast_only=false, vel_and_mom_synced=false); - }; - // ************************************************************* // The "slow" integrator for MRI and the only integrator for SRI // ************************************************************* @@ -312,15 +304,17 @@ const int nrk) { amrex::ignore_unused(nrk); - if (verbose) Print() << "Making slow rhs at time " << old_stage_time << - " for slow variables advancing from " << - old_step_time << " to " << new_stage_time << std::endl; // Note that the "old" and "new" metric terms correspond to // t^n and the RK stage (either t^*, t^** or t^{n+1} that this source // will be used to advance to Real slow_dt = new_stage_time - old_step_time; + if (verbose) amrex::Print() << "Time integration of scalars at level " << level + << " from " << old_step_time << " to " << new_stage_time + << " with dt = " << slow_dt + << " using RHS created at " << old_stage_time << std::endl; + int n_qstate = micro->Get_Qstate_Size(); #if defined(ERF_USE_NETCDF) @@ -385,6 +379,12 @@ #endif fr_as_crse, fr_as_fine); } + + // Apply boundary conditions on all the state variables that have been updated + // in both the fast and slow integrators + apply_bcs(S_new, new_stage_time, S_new[IntVars::cons].nGrow(), S_new[IntVars::xmom].nGrow(), + fast_only=false, vel_and_mom_synced=false); + }; // end slow_rhs_fun_post auto slow_rhs_fun_inc = [&](Vector& S_rhs, @@ -399,6 +399,14 @@ BL_PROFILE("slow_rhs_fun_inc"); if (verbose) Print() << "Making slow rhs at time " << old_stage_time << " for fast variables advancing from " << old_step_time << " to " << new_stage_time << std::endl; + // + // Define primitive variables for all later RK stages + // (We have already done this for the first RK step) + // + if (nrk > 0) { + int ng_cons = S_data[IntVars::cons].nGrow(); + cons_to_prim(S_data[IntVars::cons], ng_cons); + } Real slow_dt = new_stage_time - old_step_time; diff --git a/Source/TimeIntegration/ERF_TI_utils.H b/Source/TimeIntegration/ERF_TI_utils.H index 7870508a4..d6e581d48 100644 --- a/Source/TimeIntegration/ERF_TI_utils.H +++ b/Source/TimeIntegration/ERF_TI_utils.H @@ -45,10 +45,7 @@ }; /** - * This routine is called before the first step of the time integration, *and* in the case - * of a multi-stage method like RK3, this is called from "pre_update_fun" which is called - * before every subsequent stage. Since we advance the variables in conservative form, - * we must convert momentum to velocity before imposing the bcs. + * This routine is called after the scalars are updated for each RK stage */ auto apply_bcs = [&](Vector& S_data, const Real time_for_fp, int ng_cons, int ng_vel, diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 182aea709..64313e133 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -288,8 +288,6 @@ void ERF::advance_dycore(int level, // any state data (e.g. at RK stages or at the end of a timestep) mri_integrator.set_slow_rhs_pre(slow_rhs_fun_pre); mri_integrator.set_slow_rhs_post(slow_rhs_fun_post); - mri_integrator.set_pre_update (pre_update_fun); - mri_integrator.set_post_update(post_update_fun); if (solverChoice.anelastic[level]) { mri_integrator.set_slow_rhs_inc(slow_rhs_fun_inc);