From 3395bc006a7e8f438f721450f23a6069e06aa2b8 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Wed, 8 Nov 2023 15:19:16 -0800 Subject: [PATCH 1/2] Metgrid with SST and LANDMASK (#1280) * WIP initialization from met_em files * WIP metgrid initialization now timesteps in serial. * WIP metgrid initialization now checks met_em global attr * WIP started writing procedures for setting lateral boundary conditions from met_em files. * WIP progress on methods for setting lateral boundary conditions from met_em files. * WIP. Initialization & forcing from met_em files runs in serial but not yet in parallel. Currently limited to met_em generated from GFS because met_em variables differ slightly depending on the origin model processed. Runs with mositure hard-coded to 0 look great. Runs with moisture do not look great and there is either a bug or error in the methodology for calculating dry density and dry pressure. * WIP cleaned up errors from prior merge with development branch updates * WIP debugging initialization and forcing from met_em files. * WIP merging and fixing merge errors from recent commits from development. * WIP more debugging for WPS met_em*.nc initialization and forcing. * Add SST and LANDMASK read and copy. * This seems to be working but dies with MOST due to 0 eddy viscosity in ghost cell. * Amrex submodule. * Unify how we initialize time and do temporal interpolation with init from wrf and metgrid. * Cmake compile with metgrid. * Remove debug print. * Template build FABs from NETCDF to allow different data types. * Fix parallel issue on CPU. Need to implement GPU option. * Minor fixes for GPU. * Parallel run and GPU compile successful. * Cmake build with metgrid. --------- Co-authored-by: wiersema1 Co-authored-by: Ann Almgren Co-authored-by: Aaron Lattanzi --- CMake/BuildERFExe.cmake | 2 +- Exec/CMakeLists.txt | 1 + Exec/DevTests/MetGrid/CMakeLists.txt | 12 + Exec/DevTests/MetGrid/GNUmakefile | 7 +- Exec/DevTests/MetGrid/inputs | 27 +- Exec/DevTests/MetGrid/prob.H | 2 +- Exec/DevTests/MetGrid/prob.cpp | 2 +- Exec/LLJ/GNUmakefile | 2 +- Exec/RegTests/WPS_Test/GNUmakefile | 2 +- Source/BoundaryConditions/ABLMost.H | 156 ++- Source/BoundaryConditions/ABLMost.cpp | 48 +- .../BoundaryConditions_metgrid.cpp | 219 ++++ .../BoundaryConditions_wrfbdy.cpp | 4 +- Source/BoundaryConditions/ERF_FillPatch.cpp | 8 +- Source/BoundaryConditions/ERF_PhysBCFunct.cpp | 4 +- Source/BoundaryConditions/MOSTAverage.cpp | 11 +- Source/BoundaryConditions/Make.package | 1 + Source/EOS.H | 17 +- Source/ERF.H | 25 +- Source/ERF.cpp | 25 +- Source/ERF_make_new_level.cpp | 5 +- Source/IO/Make.package | 1 - Source/IO/NCBuildFABs.cpp | 152 --- Source/IO/NCInterface.cpp | 255 ++-- Source/IO/NCMultiFabFile.cpp | 8 +- Source/IO/NCPlotFile.cpp | 9 +- Source/IO/NCWpsFile.H | 192 ++- Source/IO/Plotfile.cpp | 4 +- Source/IO/ReadFromMetgrid.cpp | 129 +- Source/IO/ReadFromWRFBdy.cpp | 38 +- Source/IO/ReadFromWRFInput.cpp | 34 +- Source/IO/writeJobInfo.cpp | 4 +- Source/IndexDefines.H | 11 + .../Initialization/ERF_init_from_metgrid.cpp | 1079 +++++++++++++---- .../Initialization/ERF_init_from_wrfinput.cpp | 4 + Source/Initialization/Make.package | 1 + Source/Initialization/Metgrid_utils.H | 347 ++++++ Source/TimeIntegration/ERF_fast_rhs_T.cpp | 2 +- Source/TimeIntegration/TI_slow_rhs_fun.H | 29 +- Source/Utils/InteriorGhostCells.cpp | 63 +- Source/Utils/Utils.H | 3 +- 41 files changed, 2148 insertions(+), 797 deletions(-) create mode 100644 Exec/DevTests/MetGrid/CMakeLists.txt create mode 100644 Source/BoundaryConditions/BoundaryConditions_metgrid.cpp delete mode 100644 Source/IO/NCBuildFABs.cpp create mode 100644 Source/Initialization/Metgrid_utils.H diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 7f44ff5ff..215216e4c 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -45,7 +45,6 @@ function(build_erf_lib erf_lib_name) if(ERF_ENABLE_NETCDF) target_sources(${erf_lib_name} PRIVATE - ${SRC_DIR}/IO/NCBuildFABs.cpp ${SRC_DIR}/IO/NCInterface.cpp ${SRC_DIR}/IO/NCPlotFile.cpp ${SRC_DIR}/IO/NCCheckpoint.cpp @@ -113,6 +112,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/BoundaryConditions/BoundaryConditions_zvel.cpp ${SRC_DIR}/BoundaryConditions/BoundaryConditions_bndryreg.cpp ${SRC_DIR}/BoundaryConditions/BoundaryConditions_wrfbdy.cpp + ${SRC_DIR}/BoundaryConditions/BoundaryConditions_metgrid.cpp ${SRC_DIR}/BoundaryConditions/ERF_FillPatch.cpp ${SRC_DIR}/BoundaryConditions/ERF_FillPatcher.cpp ${SRC_DIR}/BoundaryConditions/ERF_PhysBCFunct.cpp diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index 6474d52a0..30ea921c9 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -24,4 +24,5 @@ else () add_subdirectory(DevTests/MovingTerrain) add_subdirectory(DevTests/ParticlesOverWoA) add_subdirectory(DevTests/MiguelDev) + add_subdirectory(DevTests/MetGrid) endif() diff --git a/Exec/DevTests/MetGrid/CMakeLists.txt b/Exec/DevTests/MetGrid/CMakeLists.txt new file mode 100644 index 000000000..621a7e0aa --- /dev/null +++ b/Exec/DevTests/MetGrid/CMakeLists.txt @@ -0,0 +1,12 @@ +set(erf_exe_name erf_metgrid_dev) + +add_executable(${erf_exe_name} "") +target_sources(${erf_exe_name} + PRIVATE + prob.cpp +) + +target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + +include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake) +build_erf_exe(${erf_exe_name}) diff --git a/Exec/DevTests/MetGrid/GNUmakefile b/Exec/DevTests/MetGrid/GNUmakefile index 838ce130e..0b412a6bb 100644 --- a/Exec/DevTests/MetGrid/GNUmakefile +++ b/Exec/DevTests/MetGrid/GNUmakefile @@ -19,7 +19,7 @@ USE_HIP = FALSE USE_SYCL = FALSE # Debugging -DEBUG = FALSE +#DEBUG = FALSE DEBUG = TRUE TEST = TRUE @@ -27,6 +27,11 @@ USE_ASSERTION = TRUE USE_NETCDF = TRUE +#USE_MOISTURE = FALSE +USE_MOISTURE = TRUE + +#USE_WARM_NO_PRECIP = TRUE + # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index cbe1f1f85..2be9c6f2e 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -1,11 +1,13 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 1 +max_step = 100 amrex.fpe_trap_invalid = 1 +amrex.fpe_trap_zero = 1 +amrex.fpe_trap_overflow = 1 # PROBLEM SIZE & GEOMETRY -geometry.prob_extent = 10600 5000 5000 -amr.n_cell = 140 80 60 +geometry.prob_extent = 28000 16000 8000 +amr.n_cell = 140 80 100 geometry.is_periodic = 0 0 0 @@ -17,9 +19,10 @@ zlo.type = "NoSlipWall" zhi.type = "SlipWall" # TIME STEP CONTROL -erf.fixed_dt = 0.1 # fixed time step depending on grid resolution +erf.fixed_dt = 1.0 # fixed time step depending on grid resolution erf.use_terrain = true +erf.terrain_smoothing = 2 # DIAGNOSTICS & VERBOSITY erf.sum_interval = -1 # timesteps between computing mass @@ -36,22 +39,26 @@ erf.check_int = 100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 10 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys +erf.plot_int_1 = 1 # number of timesteps between plotfiles +erf.plot_vars_1 = qv Rhoqv density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres KE QKE # SOLVER CHOICE -erf.alpha_T = 0.0 +erf.alpha_T = 1.0 erf.alpha_C = 1.0 erf.use_gravity = true erf.molec_diff_type = "None" erf.les_type = "Smagorinsky" +#erf.les_type = "Deardorff" erf.Cs = 0.1 +#erf.terrain_z_levels = 0 130 354 583 816 1054 1549 2068 2615 3193 3803 4450 5142 5892 6709 7603 8591 9702 10967 12442 14230 16610 18711 20752 22133 23960 26579 28493 31236 33699 36068 40000 + # INITIALIZATION WITH ATM DATA -erf.init_type = "metgrid" -erf.nc_init_file_0 = "met_em_d01.nc" -#erf.nc_bdy_file = "" +erf.metgrid_bdy_width = 5 +erf.metgrid_bdy_set_width = 1 +erf.init_type = "metgrid" +erf.nc_init_file_0 = "met_em.d01.2022-06-18_00:00:00.nc" "met_em.d01.2022-06-18_06:00:00.nc" "met_em.d01.2022-06-18_12:00:00.nc" "met_em.d01.2022-06-18_18:00:00.nc" #There will be no OpenMP tiling fabarray.mfiter_tile_size = 1024 1024 1024 diff --git a/Exec/DevTests/MetGrid/prob.H b/Exec/DevTests/MetGrid/prob.H index 18d9e2892..6ad4ec3fe 100644 --- a/Exec/DevTests/MetGrid/prob.H +++ b/Exec/DevTests/MetGrid/prob.H @@ -12,7 +12,7 @@ struct ProbParm : ProbParmDefaults { class Problem : public ProblemBase { public: - Problem(); + Problem(const amrex::Real* problo, const amrex::Real* probhi); protected: std::string name() override { return "Metgrid"; } diff --git a/Exec/DevTests/MetGrid/prob.cpp b/Exec/DevTests/MetGrid/prob.cpp index 4d9c43e75..bedd22dbd 100644 --- a/Exec/DevTests/MetGrid/prob.cpp +++ b/Exec/DevTests/MetGrid/prob.cpp @@ -8,5 +8,5 @@ amrex_probinit(const amrex_real* problo, const amrex_real* probhi) return std::make_unique(problo, probhi); } -Problem::Problem() +Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi) {} diff --git a/Exec/LLJ/GNUmakefile b/Exec/LLJ/GNUmakefile index 31c46b692..4a773e7b1 100644 --- a/Exec/LLJ/GNUmakefile +++ b/Exec/LLJ/GNUmakefile @@ -25,7 +25,7 @@ TEST = TRUE USE_ASSERTION = TRUE USE_NETCDF = TRUE -USE_HDF5 = TRUE +#USE_HDF5 = TRUE # GNU Make Bpack := ./Make.package diff --git a/Exec/RegTests/WPS_Test/GNUmakefile b/Exec/RegTests/WPS_Test/GNUmakefile index ad0a36000..f83bb0ff5 100644 --- a/Exec/RegTests/WPS_Test/GNUmakefile +++ b/Exec/RegTests/WPS_Test/GNUmakefile @@ -27,7 +27,7 @@ USE_ASSERTION = TRUE USE_MOISTURE = TRUE USE_NETCDF = TRUE -USE_HDF5 = TRUE +#USE_HDF5 = TRUE # GNU Make Bpack := ./Make.package diff --git a/Source/BoundaryConditions/ABLMost.H b/Source/BoundaryConditions/ABLMost.H index 9d6fd6bac..2b805a060 100644 --- a/Source/BoundaryConditions/ABLMost.H +++ b/Source/BoundaryConditions/ABLMost.H @@ -33,8 +33,14 @@ public: amrex::Vector>& vars_old, amrex::Vector>& Theta_prim, amrex::Vector>& z_phys_nd, - int ng_in) - : m_geom(geom), m_ma(geom,vars_old,Theta_prim,z_phys_nd) + amrex::Vector>>& sst_lev, + amrex::Vector>>& lmask_lev, + amrex::Real start_bdy_time = 0.0, + amrex::Real bdy_time_interval = 0.0) + : m_start_bdy_time(start_bdy_time), + m_bdy_time_interval(bdy_time_interval), + m_geom(geom), + m_ma(geom,vars_old,Theta_prim,z_phys_nd) { amrex::ParmParse pp("erf"); pp.query("most.z0", z0_const); @@ -86,6 +92,7 @@ public: amrex::Abort("Undefined MOST roughness type!"); } + // Size the MOST params for all levels int nlevs = m_geom.size(); z_0.resize(nlevs); u_star.resize(nlevs); @@ -93,11 +100,37 @@ public: t_surf.resize(nlevs); olen.resize(nlevs); - for (int lev = 0; lev < nlevs; lev++) - { - // Z0 heights + // Get pointers to SST and LANDMASK data + m_sst_lev.resize(nlevs); + m_lmask_lev.resize(nlevs); + for (int lev(0); lev(ba2d,dm,ncomp,ng); - u_star[lev]->setVal(1.E34); - - t_star[lev] = std::make_unique(ba2d,dm,ncomp,ng); - t_star[lev]->setVal(1.E34); - - olen[lev] = std::make_unique(ba2d,dm,ncomp,ng); - olen[lev]->setVal(1.E34); - - t_surf[lev] = std::make_unique(ba2d,dm,ncomp,ng); - if (erf_st) { - t_surf[lev]->setVal(surf_temp); - } else { - t_surf[lev]->setVal(0.0); - } + u_star[lev] = std::make_unique(ba2d,dm,ncomp,ng); + u_star[lev]->setVal(1.E34); + + t_star[lev] = std::make_unique(ba2d,dm,ncomp,ng); + t_star[lev]->setVal(1.E34); + + olen[lev] = std::make_unique(ba2d,dm,ncomp,ng); + olen[lev]->setVal(1.E34); + + t_surf[lev] = std::make_unique(ba2d,dm,ncomp,ng); + + if (m_sst_lev[lev][0]) { // Valid SST data at t==0 + theta_type = ThetaCalcType::SURFACE_TEMPERATURE; + amrex::MultiFab::Copy(*(t_surf[lev]), *(m_sst_lev[lev][0]), 0, 0, 1, ng); + } else if (erf_st) { // Constant temp + t_surf[lev]->setVal(surf_temp); + } else { + t_surf[lev]->setVal(0.0); } }// lev } void - update_surf_temp (amrex::Real cur_time) - { - if (surf_heating_rate != 0) { - int nlevs = m_geom.size(); - for (int lev = 0; lev < nlevs; lev++) - { - t_surf[lev]->setVal(surf_temp + surf_heating_rate*cur_time); - amrex::Print() << "Surface temp at t=" << cur_time - << ": " - << surf_temp + surf_heating_rate*cur_time - << std::endl; - } - } - } + update_fluxes (const int& lev, + const amrex::Real& time, + int max_iters = 25); + + template + void + compute_fluxes (const int& lev, + const int& max_iters, + const FluxIter& most_flux); void - impose_most_bcs (const int lev, + impose_most_bcs (const int& lev, const amrex::Vector& mfs, amrex::MultiFab* eddyDiffs, amrex::MultiFab* z_phys); template void - compute_most_bcs (const int lev, + compute_most_bcs (const int& lev, const amrex::Vector& mfs, amrex::MultiFab* eddyDiffs, amrex::MultiFab* z_phys, @@ -169,33 +187,42 @@ public: const FluxCalc& flux_comp); void - update_fluxes (int lev, - amrex::Real cur_time, - int max_iters = 25); + time_interp_tsurf(const int& lev, + const amrex::Real& time); - template void - compute_fluxes (const int& lev, - const int& max_iters, - const FluxIter& most_flux); + update_surf_temp (const amrex::Real& time) + { + if (surf_heating_rate != 0) { + int nlevs = m_geom.size(); + for (int lev = 0; lev < nlevs; lev++) + { + t_surf[lev]->setVal(surf_temp + surf_heating_rate*time); + amrex::Print() << "Surface temp at t=" << time + << ": " + << surf_temp + surf_heating_rate*time + << std::endl; + } + } + } void - update_mac_ptrs (int lev, + update_mac_ptrs (const int& lev, amrex::Vector>& vars_old, amrex::Vector>& Theta_prim) { m_ma.update_field_ptrs(lev,vars_old,Theta_prim); } const amrex::MultiFab* - get_u_star (int lev) { return u_star[lev].get(); } + get_u_star (const int& lev) { return u_star[lev].get(); } const amrex::MultiFab* - get_t_star (int lev) { return t_star[lev].get(); } + get_t_star (const int& lev) { return t_star[lev].get(); } const amrex::MultiFab* - get_olen (int lev) { return olen[lev].get(); } + get_olen (const int& lev) { return olen[lev].get(); } const amrex::MultiFab* - get_mac_avg (int lev, int comp) { return m_ma.get_average(lev,comp); } + get_mac_avg (const int& lev, int comp) { return m_ma.get_average(lev,comp); } enum struct FluxCalcType { MOENG = 0, ///< Moeng functional form @@ -225,6 +252,8 @@ private: amrex::Real surf_temp_flux{0}; amrex::Real cnk_a{0.0185}; amrex::Real depth{30.0}; + amrex::Real m_start_bdy_time; + amrex::Real m_bdy_time_interval; amrex::Vector m_geom; amrex::Vector z_0; @@ -233,6 +262,9 @@ private: amrex::Vector> t_star; amrex::Vector> olen; amrex::Vector> t_surf; + + amrex::Vector> m_sst_lev; + amrex::Vector> m_lmask_lev; }; #endif /* ABLMOST_H */ diff --git a/Source/BoundaryConditions/ABLMost.cpp b/Source/BoundaryConditions/ABLMost.cpp index c7daf5bbe..e088eeee3 100644 --- a/Source/BoundaryConditions/ABLMost.cpp +++ b/Source/BoundaryConditions/ABLMost.cpp @@ -10,10 +10,13 @@ using namespace amrex; * @param[in] max_iters maximum iterations to use */ void -ABLMost::update_fluxes (int lev, - Real cur_time, +ABLMost::update_fluxes (const int& lev, + const Real& time, int max_iters) { + // Update SST data if we have a valid pointer + if (m_sst_lev[lev][0]) time_interp_tsurf(lev, time); + // Compute plane averages for all vars (regardless of flux type) m_ma.compute_averages(lev); @@ -31,7 +34,7 @@ ABLMost::update_fluxes (int lev, compute_fluxes(lev, max_iters, most_flux); } } else if (theta_type == ThetaCalcType::SURFACE_TEMPERATURE) { - update_surf_temp(cur_time); + update_surf_temp(time); if (rough_type == RoughCalcType::CONSTANT) { surface_temp most_flux(m_ma.get_zref(), surf_temp_flux); compute_fluxes(lev, max_iters, most_flux); @@ -75,12 +78,9 @@ ABLMost::compute_fluxes (const int& lev, const auto *const tm_ptr = m_ma.get_average(lev,2); const auto *const umm_ptr = m_ma.get_average(lev,3); - // Ghost cells - amrex::IntVect ng = u_star[lev]->nGrowVect(); ng[2]=0; - for (MFIter mfi(*u_star[lev]); mfi.isValid(); ++mfi) { - amrex::Box gtbx = mfi.growntilebox(ng); + Box gtbx = mfi.growntilebox(); auto u_star_arr = u_star[lev]->array(mfi); auto t_star_arr = t_star[lev]->array(mfi); @@ -108,7 +108,7 @@ ABLMost::compute_fluxes (const int& lev, * @param[in] eddyDiffs Diffusion coefficients from turbulence model */ void -ABLMost::impose_most_bcs (const int lev, +ABLMost::impose_most_bcs (const int& lev, const Vector& mfs, MultiFab* eddyDiffs, MultiFab* z_phys) @@ -134,7 +134,7 @@ ABLMost::impose_most_bcs (const int lev, */ template void -ABLMost::compute_most_bcs (const int lev, +ABLMost::compute_most_bcs (const int& lev, const Vector& mfs, MultiFab* eddyDiffs, MultiFab* z_phys, @@ -218,3 +218,33 @@ ABLMost::compute_most_bcs (const int lev, } // var_idx } // mfiter } + +void +ABLMost::time_interp_tsurf(const int& lev, + const Real& time) +{ + // Time interpolation + Real dT = m_bdy_time_interval; + Real time_since_start = time - m_start_bdy_time; + int n_time = static_cast( time_since_start / dT); + amrex::Real alpha = (time_since_start - n_time * dT) / dT; + AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); + amrex::Real oma = 1.0 - alpha; + AMREX_ALWAYS_ASSERT( (n_time >= 0) && (n_time < (m_sst_lev[lev].size()-1))); + + // Populate t_surf + for (MFIter mfi(*t_surf[lev]); mfi.isValid(); ++mfi) + { + Box gtbx = mfi.growntilebox(); + + auto t_surf_arr = t_surf[lev]->array(mfi); + const auto sst_hi_arr = m_sst_lev[lev][n_time+1]->const_array(mfi); + const auto sst_lo_arr = m_sst_lev[lev][n_time ]->const_array(mfi); + + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + t_surf_arr(i,j,k) = oma * sst_lo_arr(i,j,k) + + alpha * sst_hi_arr(i,j,k); + }); + } +} diff --git a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp new file mode 100644 index 000000000..37f4dc04f --- /dev/null +++ b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp @@ -0,0 +1,219 @@ +#include "ERF.H" +#include "Utils.H" + +using namespace amrex; + +#ifdef ERF_USE_NETCDF +/* + * Impose boundary conditions using data read in from met_em files + * + * @param[out] mfs Vector of MultiFabs to be filled + * @param[in] time time at which the data should be filled + */ + +void +ERF::fill_from_metgrid (const Vector& mfs, + const Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons) +{ + int lev = 0; + + // Time interpolation + Real dT = bdy_time_interval; + Real time_since_start = time - start_bdy_time; + int n_time = static_cast( time_since_start / dT); + amrex::Real alpha = (time_since_start - n_time * dT) / dT; + AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); + amrex::Real oma = 1.0 - alpha; + + // Flags for read vars and index mapping + // See IndexDefines.H +#if defined(ERF_USE_MOISTURE) + // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQt RhoQp NumVars] + Vector cons_read = {1, 1, 0, 0, 0, 1, 0}; + Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0}; +#elif defined(ERF_USE_WARM_NO_PRECIP) + // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQv RhoQc NumVars] + Vector cons_read = {1, 1, 0, 0, 0, 1, 0}; + Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0}; +# else + // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar NumVars] + Vector cons_read = {1, 1, 0, 0, 0}; + Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0}; +#endif + + Vector> is_read; + is_read.push_back( cons_read ); + is_read.push_back( {1} ); // xvel + is_read.push_back( {1} ); // yvel + is_read.push_back( {0} ); // zvel + + Vector> ind_map; + ind_map.push_back( cons_map ); + ind_map.push_back( {MetGridBdyVars::U} ); // xvel + ind_map.push_back( {MetGridBdyVars::V} ); // yvel + ind_map.push_back( {0} ); // zvel + + // Nvars to loop over + Vector comp_var = {ncomp_cons, 1, 1, 1}; + + // End of vars loop + int var_idx_end = (cons_only) ? Vars::cons + 1 : Vars::NumTypes; + + // Loop over all variable types + for (int var_idx = Vars::cons; var_idx < var_idx_end; ++var_idx) + { + MultiFab& mf = *mfs[var_idx]; + + // + // Note that "domain" is mapped onto the type of box the data is in + // + Box domain = geom[lev].Domain(); + domain.convert(mf.boxArray().ixType()); + const auto& dom_lo = amrex::lbound(domain); + const auto& dom_hi = amrex::ubound(domain); + + // Offset only applys to cons (we may fill a subset of these vars) + int offset = (var_idx == Vars::cons) ? icomp_cons : 0; + + // Loop over each component + for (int comp_idx(0); comp_idx < comp_var[var_idx]; ++comp_idx) + { + int width = metgrid_bdy_set_width; + + // Variable can be read from met_em files + //------------------------------------ + if (is_read[var_idx][comp_idx]) + { + int ivar = ind_map[var_idx][comp_idx]; + IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; + +// if (ivar == MetGridBdyVars::U) { +// amrex::Print() << "fill_from_metgrid U var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; +// } else if (ivar == MetGridBdyVars::V) { +// amrex::Print() << "fill_from_metgrid V var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; +// } else if (ivar == MetGridBdyVars::R) { +// amrex::Print() << "fill_from_metgrid R var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; +// } else if (ivar == MetGridBdyVars::T) { +// amrex::Print() << "fill_from_metgrid T var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; +// } else if (ivar == MetGridBdyVars::QV) { +// amrex::Print() << "fill_from_metgrid QV var_idx=" << var_idx << " comp_idx=" << comp_idx << " ivar=" << ivar << std::endl; +// } else { +// amrex::Print() << "fill_from_metgrid UNKNOWN" << std::endl; +// } + + // We have data at fixed time intervals we will call dT + // Then to interpolate, given time, we can define n = (time/dT) + // and alpha = (time - n*dT) / dT, then we define the data at time + // as alpha * (data at time n+1) + (1 - alpha) * (data at time n) + const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar].const_array(); + const auto& bdatxlo_np1 = bdy_data_xlo[n_time+1][ivar].const_array(); + const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array(); + const auto& bdatxhi_np1 = bdy_data_xhi[n_time+1][ivar].const_array(); + const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array(); + const auto& bdatylo_np1 = bdy_data_ylo[n_time+1][ivar].const_array(); + const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array(); + const auto& bdatyhi_np1 = bdy_data_yhi[n_time+1][ivar].const_array(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // Grown tilebox so we fill exterior ghost cells as well + Box gbx = mfi.growntilebox(ng_vect); + const Array4& dest_arr = mf.array(mfi); + + // Call w/o interior ghost cells + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(gbx, domain, width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, ng_vect); + + // x-faces (includes exterior y ghost cells) + ParallelFor(bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + int ii = std::max(i , dom_lo.x); + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + dest_arr(i,j,k,comp_idx) = oma * bdatxlo_n (ii,jj,k,0) + + alpha * bdatxlo_np1(ii,jj,k,0); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + int ii = std::min(i , dom_hi.x); + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + dest_arr(i,j,k,comp_idx) = oma * bdatxhi_n (ii,jj,k,0) + + alpha * bdatxhi_np1(ii,jj,k,0); + }); + + // y-faces (do not include exterior x ghost cells) + ParallelFor(bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + int jj = std::max(j , dom_lo.y); + dest_arr(i,j,k,comp_idx) = oma * bdatylo_n (i,jj,k,0) + + alpha * bdatylo_np1(i,jj,k,0); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + int jj = std::min(j , dom_hi.y); + dest_arr(i,j,k,comp_idx) = oma * bdatyhi_n (i,jj,k,0) + + alpha * bdatyhi_np1(i,jj,k,0); + }); + } // mfi + + // Variable not read or computed from met_em files + //------------------------------------ + } else { + width = metgrid_bdy_width; + IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0; + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // Grown tilebox so we fill exterior ghost cells as well + Box gbx = mfi.growntilebox(ng_vect); + const Array4& dest_arr = mf.array(mfi); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + compute_interior_ghost_bxs_xy(gbx, domain, width, 0, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, ng_vect); + + // x-faces (includes y ghost cells) + ParallelFor(bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + dest_arr(i,j,k,comp_idx) = dest_arr(dom_lo.x+width,jj,k,comp_idx); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + dest_arr(i,j,k,comp_idx) = dest_arr(dom_hi.x-width,jj,k,comp_idx); + }); + + // y-faces (does not include x ghost cells) + ParallelFor(bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k,comp_idx) = dest_arr(i,dom_lo.y+width,k,comp_idx); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k,comp_idx) = dest_arr(i,dom_hi.y-width,k,comp_idx); + }); + } // mfi + } // is_read + } // comp + } // var +} +#endif diff --git a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp b/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp index 6abb2221d..0153430ab 100644 --- a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp @@ -22,8 +22,8 @@ ERF::fill_from_wrfbdy (const Vector& mfs, // Time interpolation Real dT = bdy_time_interval; - Real time_since_start = (time - start_bdy_time) / 1.e10; - int n_time = static_cast( time_since_start / dT); + Real time_since_start = time - start_bdy_time; + int n_time = static_cast( time_since_start / dT); amrex::Real alpha = (time_since_start - n_time * dT) / dT; AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); amrex::Real oma = 1.0 - alpha; diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index bf750a410..d4c863c1e 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -92,7 +92,8 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) #ifdef ERF_USE_NETCDF // We call this here because it is an ERF routine - if (init_type=="real" && lev==0) fill_from_wrfbdy(mfs, time); + if (init_type == "real" && lev==0) fill_from_wrfbdy(mfs,time); + if (init_type == "metgrid" && lev==0) fill_from_metgrid(mfs,time); #endif if (m_r2d) fill_from_bndryregs(mfs,time); @@ -128,7 +129,7 @@ ERF::FillPatchMoistVars (int lev, MultiFab& mf) IntVect ngvect_cons = mf.nGrowVect(); IntVect ngvect_vels = {0,0,0}; - if (init_type != "real") { + if ((init_type != "real") and (init_type != "metgrid")) { (*physbcs[lev])({&mf},icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,bccomp_cons); } @@ -248,7 +249,8 @@ ERF::FillIntermediatePatch (int lev, Real time, // --- i.e., cons_only and which cons indices (icomp_cons & ncomp_cons) // We call this here because it is an ERF routine - if (init_type=="real" && lev==0) fill_from_wrfbdy(mfs, time, cons_only, icomp_cons, ncomp_cons); + if (init_type == "real" && lev==0) fill_from_wrfbdy(mfs,time, cons_only, icomp_cons, ncomp_cons); + if (init_type == "metgrid" && lev==0) fill_from_metgrid(mfs,time, cons_only, icomp_cons, ncomp_cons); #endif if (m_r2d) fill_from_bndryregs(mfs,time); diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index ed3517ab2..eafea677a 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -82,7 +82,7 @@ void ERFPhysBCFunct::operator() (const Vector& mfs, int icomp_cons, i z_nd_arr = m_z_phys_nd->const_array(mfi); } - if (init_type != "real") + if ((init_type != "real") and (init_type != "metgrid")) { //! If there are cells not in the valid + periodic grown box @@ -109,7 +109,7 @@ void ERFPhysBCFunct::operator() (const Vector& mfs, int icomp_cons, i impose_lateral_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx,domain,z_nd_arr,dxInv,BCVars::zvel_bc); } // !cons_only - } // init_type != "real" + } // init_type != "real" and init_type != "metgrid" // Every grid touches the bottom and top domain boundary so we call the vertical bcs // for every box -- and we need to call these even if init_type == real diff --git a/Source/BoundaryConditions/MOSTAverage.cpp b/Source/BoundaryConditions/MOSTAverage.cpp index 26b812428..a70b7a0ec 100644 --- a/Source/BoundaryConditions/MOSTAverage.cpp +++ b/Source/BoundaryConditions/MOSTAverage.cpp @@ -304,12 +304,13 @@ MOSTAverage::set_k_indices_T() auto k_arr = m_k_indx[lev]->array(mfi); ParallelFor(npbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + k_arr(i,j,k) = 0; Real z_target = d_zref + z_phys_arr(i,j,k); for (int lk(0); lk<=kmax; ++lk) { Real z_lo = 0.25 * ( z_phys_arr(i,j ,lk ) + z_phys_arr(i+1,j ,lk ) - + z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) ); + + z_phys_arr(i,j+1,lk ) + z_phys_arr(i+1,j+1,lk ) ); Real z_hi = 0.25 * ( z_phys_arr(i,j ,lk+1) + z_phys_arr(i+1,j ,lk+1) - + z_phys_arr(i,j+1,lk+1) + z_phys_arr(i+1,j+1,lk+1) ); + + z_phys_arr(i,j+1,lk+1) + z_phys_arr(i+1,j+1,lk+1) ); if (z_target > z_lo && z_target < z_hi){ AMREX_ASSERT_WITH_MESSAGE(lk >= d_radius, "K index must be larger than averaging radius!"); @@ -376,9 +377,9 @@ MOSTAverage::set_norm_indices_T() Real z_target = delta_z + z_phys_arr(i,j,k); for (int lk(0); lk<=kmax; ++lk) { Real z_lo = 0.25 * ( z_phys_arr(i_new,j_new ,lk ) + z_phys_arr(i_new+1,j_new ,lk ) - + z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) ); + + z_phys_arr(i_new,j_new+1,lk ) + z_phys_arr(i_new+1,j_new+1,lk ) ); Real z_hi = 0.25 * ( z_phys_arr(i_new,j_new ,lk+1) + z_phys_arr(i_new+1,j_new ,lk+1) - + z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) ); + + z_phys_arr(i_new,j_new+1,lk+1) + z_phys_arr(i_new+1,j_new+1,lk+1) ); if (z_target > z_lo && z_target < z_hi){ AMREX_ASSERT_WITH_MESSAGE(lk >= d_radius, "K index must be larger than averaging radius!"); @@ -532,7 +533,7 @@ MOSTAverage::compute_plane_averages(int lev) // Peel back the level auto& fields = m_fields[lev]; auto& averages = m_averages[lev]; - const auto & geom = m_geom[lev]; + const auto & geom = m_geom[lev]; auto& z_phys = m_z_phys_nd[lev]; auto& x_pos = m_x_pos[lev]; diff --git a/Source/BoundaryConditions/Make.package b/Source/BoundaryConditions/Make.package index 56ffeafaf..b9bb5be57 100644 --- a/Source/BoundaryConditions/Make.package +++ b/Source/BoundaryConditions/Make.package @@ -4,6 +4,7 @@ CEXE_sources += BoundaryConditions_zvel.cpp CEXE_sources += BoundaryConditions_cons.cpp CEXE_sources += BoundaryConditions_bndryreg.cpp CEXE_sources += BoundaryConditions_wrfbdy.cpp +CEXE_sources += BoundaryConditions_metgrid.cpp CEXE_headers += ABLMost.H CEXE_sources += ABLMost.cpp diff --git a/Source/EOS.H b/Source/EOS.H index cc4604878..2ce071025 100644 --- a/Source/EOS.H +++ b/Source/EOS.H @@ -7,7 +7,20 @@ #include /** - * Function to return temperature given density and potential temperatue + * Function to return potential temperature given pressure and temperature + * + * @params[in] pressure + * @params[in] temperature + * @params[in] rd0cp ratio of R_d to c_p +*/ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real getThgivenPandT(const amrex::Real T, const amrex::Real P, const amrex::Real rdOcp) +{ + return T*std::pow(p_0/P, rdOcp); +} + +/** + * Function to return temperature given density and potential temperature * * @params[in] rho density * @params[in] rhotheta density times potential temperature theta @@ -24,7 +37,7 @@ amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, } /** - * Function to return potential temperature given density and temperatue + * Function to return potential temperature given density and temperature * * @params[in] rho density * @params[in] T temperature diff --git a/Source/ERF.H b/Source/ERF.H index e7e5d533d..60d649eab 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -301,6 +301,11 @@ public: bool cons_only = false, int icomp_cons = 0, int ncomp_cons = NVAR); + void fill_from_metgrid (const amrex::Vector& mfs, + amrex::Real time, + bool cons_only = false, + int icomp_cons = 0, + int ncomp_cons = NVAR); #endif #ifdef ERF_USE_PARTICLES @@ -308,6 +313,11 @@ public: static amrex::Vector tracer_particle_varnames; #endif +#ifdef ERF_USE_NETCDF + void init_from_wrfinput (int lev); + void init_from_metgrid (int lev); +#endif // ERF_USE_NETCDF + private: /////////////////////////// @@ -340,11 +350,6 @@ private: void initialize_integrator (int lev, amrex::MultiFab& cons_mf, amrex::MultiFab& vel_mf); -#ifdef ERF_USE_NETCDF - void init_from_wrfinput (int lev); - void init_from_metgrid (int lev); -#endif // ERF_USE_NETCDF - // more flexible version of AverageDown() that lets you average down across multiple levels void AverageDownTo (int crse_lev); // NOLINT @@ -546,6 +551,10 @@ private: amrex::Vector> eddyDiffs_lev; amrex::Vector> SmnSmn_lev; + // Sea Surface Temps and Land Masks (lev, ntimes) + amrex::Vector>> sst_lev; + amrex::Vector>> lmask_lev; + // Other SFS terms amrex::Vector> SFS_hfx1_lev, SFS_hfx2_lev, SFS_hfx3_lev; amrex::Vector> SFS_diss_lev; @@ -568,8 +577,6 @@ private: amrex::Vector> mapfac_u; amrex::Vector> mapfac_v; - amrex::Vector> sst; - amrex::Vector base_state; amrex::Vector base_state_new; @@ -698,6 +705,10 @@ private: int wrfbdy_width{0}; int wrfbdy_set_width{0}; + // NetCDF initialization (met_em) file + int metgrid_bdy_width{0}; + int metgrid_bdy_set_width{0}; + // Text input_sounding file static std::string input_sounding_file; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 6ecdc9384..9dff3a68b 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -178,6 +178,10 @@ ERF::ERF () eddyDiffs_lev.resize(nlevs_max); SmnSmn_lev.resize(nlevs_max); + // Sea surface temps + sst_lev.resize(nlevs_max); + lmask_lev.resize(nlevs_max); + // Metric terms z_phys_nd.resize(nlevs_max); z_phys_cc.resize(nlevs_max); @@ -554,8 +558,12 @@ ERF::InitData () // WritePlotFile calls FillPatch in order to compute gradients if (phys_bc_type[Orientation(Direction::z,Orientation::low)] == ERF_BC::MOST) { - int ng_for_most = ComputeGhostCells(solverChoice.advChoice,solverChoice.use_NumDiff)+1; - m_most = std::make_unique(geom,vars_old,Theta_prim,z_phys_nd,ng_for_most); + m_most = std::make_unique(geom, vars_old, Theta_prim, z_phys_nd, + sst_lev, lmask_lev +#ifdef ERF_USE_NETCDF + ,start_bdy_time, bdy_time_interval +#endif + ); // We now configure ABLMost params here so that we can print the averages at t=0 // Note we don't fill ghost cells here because this is just for diagnostics @@ -566,7 +574,7 @@ ERF::InitData () MultiFab::Copy( *Theta_prim[lev], S, Cons::RhoTheta, 0, 1, ng); MultiFab::Divide(*Theta_prim[lev], S, Cons::Rho , 0, 1, ng); m_most->update_mac_ptrs(lev, vars_new, Theta_prim); - m_most->update_fluxes(lev,t_new[lev]); + m_most->update_fluxes(lev, t_new[lev]); } } @@ -1033,6 +1041,13 @@ ERF::ReadParameters () AMREX_ALWAYS_ASSERT(wrfbdy_set_width >= 0); AMREX_ALWAYS_ASSERT(wrfbdy_width >= wrfbdy_set_width); + // Query the set and total widths for metgrid_bdy interior ghost cells + pp.query("metgrid_bdy_width", metgrid_bdy_width); + pp.query("metgrid_bdy_set_width", metgrid_bdy_set_width); + AMREX_ALWAYS_ASSERT(metgrid_bdy_width >= 0); + AMREX_ALWAYS_ASSERT(metgrid_bdy_set_width >= 0); + AMREX_ALWAYS_ASSERT(metgrid_bdy_width >= metgrid_bdy_set_width); + // Query the set and total widths for crse-fine interior ghost cells pp.query("cf_width", cf_width); pp.query("cf_set_width", cf_set_width); @@ -1380,6 +1395,10 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, eddyDiffs_lev.resize(nlevs_max); SmnSmn_lev.resize(nlevs_max); + // Sea surface temps + sst_lev.resize(nlevs_max); + lmask_lev.resize(nlevs_max); + // Metric terms z_phys_nd.resize(nlevs_max); z_phys_cc.resize(nlevs_max); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 21e635c46..c94b4e635 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -145,6 +145,9 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, z_t_rk[lev] = nullptr; } + sst_lev[lev].resize(1); sst_lev[lev][0] = nullptr; + lmask_lev[lev].resize(1); lmask_lev[lev][0] = nullptr; + // ******************************************************************************************** // Define Theta_prim storage if using MOST BC // ******************************************************************************************** @@ -167,7 +170,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, // since we may need to use the grid information before constructing // initial idealized data // ******************************************************************************************** - if (init_type == "real") { + if ((init_type == "real") || (init_type == "metgrid")) { // If called from restart, the data structures for terrain-related quantities // are built in the ReadCheckpoint routine. Otherwise we build them here. diff --git a/Source/IO/Make.package b/Source/IO/Make.package index e608b0756..079a2eb1c 100644 --- a/Source/IO/Make.package +++ b/Source/IO/Make.package @@ -15,7 +15,6 @@ ifeq ($(USE_NETCDF), TRUE) CEXE_sources += ReadFromWRFBdy.cpp CEXE_sources += ReadFromWRFInput.cpp CEXE_sources += ReadFromMetgrid.cpp - CEXE_sources += NCBuildFABs.cpp CEXE_sources += NCInterface.cpp CEXE_sources += NCPlotFile.cpp CEXE_sources += NCColumnFile.cpp diff --git a/Source/IO/NCBuildFABs.cpp b/Source/IO/NCBuildFABs.cpp deleted file mode 100644 index 048803254..000000000 --- a/Source/IO/NCBuildFABs.cpp +++ /dev/null @@ -1,152 +0,0 @@ -#include -#include -#include -#include - -#include "DataStruct.H" -#include "NCInterface.H" -#include "NCWpsFile.H" -#include "AMReX_FArrayBox.H" -#include "AMReX_IndexType.H" -#include "AMReX_Print.H" - -using namespace amrex; - -using RARRAY = NDArray; - -void -fill_fab_from_arrays(int iv, Vector& nc_arrays, - const std::string& var_name, - NC_Data_Dims_Type& NC_dim_type, - FArrayBox& temp); - -/** - * Function to read NetCDF variables and fill the corresponding Array4's - * - * @param fname Name of the NetCDF file to be read - * @param nc_var_names Variable names in the NetCDF file - * @param NC_dim_types NetCDF data dimension types - * @param fab_vars Fab data we are to fill - */ -void -BuildFABsFromNetCDFFile(const Box& domain, - const std::string &fname, - Vector nc_var_names, - Vector NC_dim_types, - Vector fab_vars) -{ - int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank - - amrex::Vector nc_arrays(nc_var_names.size()); - - if (amrex::ParallelDescriptor::IOProcessor()) - { - ReadNetCDFFile(fname, nc_var_names, nc_arrays); - } - - for (int iv = 0; iv < nc_var_names.size(); iv++) - { - amrex::FArrayBox tmp; - if (amrex::ParallelDescriptor::IOProcessor()) { - fill_fab_from_arrays(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp); - } - - Box box = tmp.box(); - int ncomp = tmp.nComp(); - ParallelDescriptor::Bcast(&box, 1, ioproc); - ParallelDescriptor::Bcast(&ncomp, 1, ioproc); - - if (!amrex::ParallelDescriptor::IOProcessor()) { -#ifdef AMREX_USE_GPU - tmp.resize(box,ncomp, The_Pinned_Arena()); -#else - tmp.resize(box,ncomp); -#endif - } - - ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc); - - // Shift box by the domain lower corner - Box fab_bx = tmp.box(); - Dim3 dom_lb = lbound(domain); - fab_bx += IntVect(dom_lb.x,dom_lb.y,dom_lb.z); - // fab_vars points to data on device - fab_vars[iv]->resize(fab_bx,1); -#ifdef AMREX_USE_GPU - Gpu::copy(Gpu::hostToDevice, tmp.dataPtr(), tmp.dataPtr() + tmp.size(), fab_vars[iv]->dataPtr()); -#else - // Provided by BaseFab inheritance through FArrayBox - fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1); -#endif - } -} - -/** - * Helper function for reading data from NetCDF file into a - * provided FAB. - * - * @param iv Index for which variable we are going to fill - * @param nc_arrays Arrays of data from NetCDF file - * @param var_name Variable name - * @param NC_dim_type Dimension type for the variable as stored in the NetCDF file - * @param temp FAB where we store the variable data from the NetCDF Arrays - */ -void -fill_fab_from_arrays(int iv, Vector& nc_arrays, - const std::string& var_name, - NC_Data_Dims_Type& NC_dim_type, - FArrayBox& temp) -{ - int ns1, ns2, ns3; - if (NC_dim_type == NC_Data_Dims_Type::Time_BT) { - ns1 = nc_arrays[iv].get_vshape()[1]; - ns2 = 1; - ns3 = 1; - // amrex::Print() << "TYPE BT " << ns1 << std::endl; - } else if (NC_dim_type == NC_Data_Dims_Type::Time_SN_WE) { - ns1 = 1; - ns2 = nc_arrays[iv].get_vshape()[1]; - ns3 = nc_arrays[iv].get_vshape()[2]; - // amrex::Print() << "TYPE SN WE " << ns2 << " " << ns3 << std::endl; - } else if (NC_dim_type == NC_Data_Dims_Type::Time_BT_SN_WE) { - ns1 = nc_arrays[iv].get_vshape()[1]; - ns2 = nc_arrays[iv].get_vshape()[2]; - ns3 = nc_arrays[iv].get_vshape()[3]; - // amrex::Print() << "TYPE BT SN WE " << ns1 << " " << ns2 << " " << ns3 << std::endl; - } else { - amrex::Abort("Dont know this NC_Data_Dims_Type"); - } - - // TODO: The box will only start at (0,0,0) at level 0 -- we need to generalize this - Box my_box(IntVect(0,0,0), IntVect(ns3-1,ns2-1,ns1-1)); - // amrex::Print() <<" MY BOX " << my_box << std::endl; - - if (var_name == "U" || var_name == "UU" || - var_name == "MACFAC_U" || var_name == "MAPFAC_UY") my_box.setType(amrex::IndexType(IntVect(1,0,0))); - if (var_name == "V" || var_name == "VV" || - var_name == "MACFAC_V" || var_name == "MAPFAC_VY") my_box.setType(amrex::IndexType(IntVect(0,1,0))); - if (var_name == "W" || var_name == "WW") my_box.setType(amrex::IndexType(IntVect(0,0,1))); - -#ifdef AMREX_USE_GPU - // Make sure temp lives on CPU since nc_arrays lives on CPU only - temp.resize(my_box,1,The_Pinned_Arena()); -#else - temp.resize(my_box,1); -#endif - Array4 fab_arr = temp.array(); - - int ioff = temp.box().smallEnd()[0]; - int joff = temp.box().smallEnd()[1]; - - auto num_pts = my_box.numPts(); - - // amrex::Print() <<" ns1 * ns2 * ns3 " << ns1 * ns2 * ns3 << std::endl; - // amrex::Print() <<" NUMPTS " << num_pts << std::endl; - - for (int n(0); n < num_pts; ++n) { - int k = n / (ns2*ns3); - int j = (n - k*(ns2*ns3)) / ns3 + joff; - int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff; - fab_arr(i,j,k,0) = static_cast(*(nc_arrays[iv].get_data()+n)); - } -} diff --git a/Source/IO/NCInterface.cpp b/Source/IO/NCInterface.cpp index fb61f1288..8a2000cb0 100644 --- a/Source/IO/NCInterface.cpp +++ b/Source/IO/NCInterface.cpp @@ -16,7 +16,7 @@ char recname[NC_MAX_NAME + 1]; * * @param ierr Error flag from NetCDF */ -void check_nc_error(int ierr) +void check_nc_error (int ierr) { if (ierr != NC_NOERR) { printf("\n%s\n\n", nc_strerror(ierr)); @@ -28,7 +28,7 @@ void check_nc_error(int ierr) /** * Error-checking wrapper for NetCDF function nc_inq_dimname */ -std::string NCDim::name() const +std::string NCDim::name () const { check_nc_error(nc_inq_dimname(ncid, dimid, recname)); return std::string(recname); @@ -37,7 +37,7 @@ std::string NCDim::name() const /** * Error-checking wrapper for NetCDF function nc_inq_dimlen */ -size_t NCDim::len() const +size_t NCDim::len () const { size_t dlen; check_nc_error(nc_inq_dimlen(ncid, dimid, &dlen)); @@ -47,7 +47,7 @@ size_t NCDim::len() const /** * Error-checking wrapper for NetCDF function nc_inq_varname */ -std::string NCVar::name() const +std::string NCVar::name () const { check_nc_error(nc_inq_varname(ncid, varid, recname)); return std::string(recname); @@ -56,7 +56,7 @@ std::string NCVar::name() const /** * Error-checking wrapper for NetCDF function nc_inq_varndims */ -int NCVar::ndim() const +int NCVar::ndim () const { int ndims; check_nc_error(nc_inq_varndims(ncid, varid, &ndims)); @@ -66,7 +66,7 @@ int NCVar::ndim() const /** * Error-checking function to get the length of each dimension from a NetCDF identity */ -std::vector NCVar::shape() const +std::vector NCVar::shape () const { int ndims = ndim(); std::vector dimids(ndims); @@ -86,7 +86,7 @@ std::vector NCVar::shape() const * * @param ptr Pointer to the data to put */ -void NCVar::put(const double* ptr) const +void NCVar::put (const double* ptr) const { check_nc_error(nc_put_var_double(ncid, varid, ptr)); } @@ -96,7 +96,7 @@ void NCVar::put(const double* ptr) const * * @param ptr Pointer to the data to put */ -void NCVar::put(const float* ptr) const +void NCVar::put (const float* ptr) const { check_nc_error(nc_put_var_float(ncid, varid, ptr)); } @@ -106,7 +106,7 @@ void NCVar::put(const float* ptr) const * * @param ptr Pointer to the data to put */ -void NCVar::put(const int* ptr) const +void NCVar::put (const int* ptr) const { check_nc_error(nc_put_var_int(ncid, varid, ptr)); } @@ -118,10 +118,9 @@ void NCVar::put(const int* ptr) const * @param start Starting indices * @param count Count sizes */ -void NCVar::put( - const double* dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::put (const double* dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_put_vara_double(ncid, varid, start.data(), count.data(), dptr)); @@ -135,11 +134,10 @@ void NCVar::put( * @param count Count sizes * @param stride Stride length for the data */ -void NCVar::put( - const double* dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::put (const double* dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_put_vars_double( ncid, varid, start.data(), count.data(), stride.data(), dptr)); @@ -152,10 +150,9 @@ void NCVar::put( * @param start Starting indices * @param count Count sizes */ -void NCVar::put( - const float* dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::put (const float* dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_put_vara_float(ncid, varid, start.data(), count.data(), dptr)); @@ -169,11 +166,10 @@ void NCVar::put( * @param count Count sizes * @param stride Stride length for the data */ -void NCVar::put( - const float* dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::put (const float* dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_put_vars_float( ncid, varid, start.data(), count.data(), stride.data(), dptr)); @@ -186,10 +182,9 @@ void NCVar::put( * @param start Starting indices * @param count Count sizes */ -void NCVar::put( - const int* dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::put (const int* dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_put_vara_int(ncid, varid, start.data(), count.data(), dptr)); @@ -203,11 +198,10 @@ void NCVar::put( * @param count Count sizes * @param stride Stride length for the data */ -void NCVar::put( - const int* dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::put (const int* dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_put_vars_int( ncid, varid, start.data(), count.data(), stride.data(), dptr)); @@ -220,10 +214,9 @@ void NCVar::put( * @param start Starting indices * @param count Count sizes */ -void NCVar::put( - const char** dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::put (const char** dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_put_vara_string(ncid, varid, start.data(), count.data(), dptr)); @@ -237,11 +230,10 @@ void NCVar::put( * @param count Count sizes * @param stride Stride length for the data */ -void NCVar::put( - const char** dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::put (const char** dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_put_vars_string( ncid, varid, start.data(), count.data(), stride.data(), dptr)); @@ -252,7 +244,7 @@ void NCVar::put( * * @param ptr Pointer to the data location we use to get */ -void NCVar::get(double* ptr) const +void NCVar::get (double* ptr) const { check_nc_error(nc_get_var_double(ncid, varid, ptr)); } @@ -262,7 +254,7 @@ void NCVar::get(double* ptr) const * * @param ptr Pointer to the data location we use to get */ -void NCVar::get(float* ptr) const +void NCVar::get (float* ptr) const { check_nc_error(nc_get_var_float(ncid, varid, ptr)); } @@ -272,7 +264,7 @@ void NCVar::get(float* ptr) const * * @param ptr Pointer to the data location we use to get */ -void NCVar::get(int* ptr) const +void NCVar::get (int* ptr) const { check_nc_error(nc_get_var_int(ncid, varid, ptr)); } @@ -284,10 +276,9 @@ void NCVar::get(int* ptr) const * @param start Starting indices * @param count Count sizes */ -void NCVar::get( - double* dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::get (double* dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_get_vara_double(ncid, varid, start.data(), count.data(), dptr)); @@ -301,11 +292,10 @@ void NCVar::get( * @param count Count sizes * @param stride Stride length for the data */ -void NCVar::get( - double* dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::get (double* dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_get_vars_double( ncid, varid, start.data(), count.data(), stride.data(), dptr)); @@ -318,10 +308,9 @@ void NCVar::get( * @param start Starting indices * @param count Count sizes */ -void NCVar::get( - float* dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::get (float* dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_get_vara_float(ncid, varid, start.data(), count.data(), dptr)); @@ -335,55 +324,50 @@ void NCVar::get( * @param count Count sizes * @param stride Stride length for the data */ -void NCVar::get( - float* dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::get (float* dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_get_vars_float( ncid, varid, start.data(), count.data(), stride.data(), dptr)); } -void NCVar::get( - int* dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::get (int* dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_get_vara_int(ncid, varid, start.data(), count.data(), dptr)); } -void NCVar::get( - int* dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::get (int* dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_get_vars_int( ncid, varid, start.data(), count.data(), stride.data(), dptr)); } -void NCVar::get( - char* dptr, - const std::vector& start, - const std::vector& count) const +void NCVar::get (char* dptr, + const std::vector& start, + const std::vector& count) const { check_nc_error( nc_get_vara_text(ncid, varid, start.data(), count.data(), dptr)); } -void NCVar::get( - char* dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const +void NCVar::get (char* dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const { check_nc_error(nc_get_vars_text( ncid, varid, start.data(), count.data(), stride.data(), dptr)); } -bool NCVar::has_attr(const std::string& name) const +bool NCVar::has_attr (const std::string& name) const { int ierr; size_t lenp; @@ -391,34 +375,31 @@ bool NCVar::has_attr(const std::string& name) const return (ierr == NC_NOERR); } -void NCVar::put_attr(const std::string& name, const std::string& value) const +void NCVar::put_attr (const std::string& name, const std::string& value) const { check_nc_error( nc_put_att_text(ncid, varid, name.data(), value.size(), value.data())); } -void NCVar::put_attr( - const std::string& name, const std::vector& value) const +void NCVar::put_attr (const std::string& name, const std::vector& value) const { check_nc_error(nc_put_att_double( ncid, varid, name.data(), NC_DOUBLE, value.size(), value.data())); } -void NCVar::put_attr( - const std::string& name, const std::vector& value) const +void NCVar::put_attr (const std::string& name, const std::vector& value) const { check_nc_error(nc_put_att_float( ncid, varid, name.data(), NC_FLOAT, value.size(), value.data())); } -void NCVar::put_attr( - const std::string& name, const std::vector& value) const +void NCVar::put_attr (const std::string& name, const std::vector& value) const { check_nc_error(nc_put_att_int( ncid, varid, name.data(), NC_INT, value.size(), value.data())); } -std::string NCVar::get_attr(const std::string& name) const +std::string NCVar::get_attr (const std::string& name) const { size_t lenp; std::vector aval; @@ -428,7 +409,7 @@ std::string NCVar::get_attr(const std::string& name) const return std::string{aval.begin(), aval.end()}; } -void NCVar::get_attr(const std::string& name, std::vector& values) const +void NCVar::get_attr (const std::string& name, std::vector& values) const { size_t lenp; check_nc_error(nc_inq_attlen(ncid, varid, name.data(), &lenp)); @@ -436,7 +417,7 @@ void NCVar::get_attr(const std::string& name, std::vector& values) const check_nc_error(nc_get_att_double(ncid, varid, name.data(), values.data())); } -void NCVar::get_attr(const std::string& name, std::vector& values) const +void NCVar::get_attr (const std::string& name, std::vector& values) const { size_t lenp; check_nc_error(nc_inq_attlen(ncid, varid, name.data(), &lenp)); @@ -444,7 +425,7 @@ void NCVar::get_attr(const std::string& name, std::vector& values) const check_nc_error(nc_get_att_float(ncid, varid, name.data(), values.data())); } -void NCVar::get_attr(const std::string& name, std::vector& values) const +void NCVar::get_attr (const std::string& name, std::vector& values) const { size_t lenp; check_nc_error(nc_inq_attlen(ncid, varid, name.data(), &lenp)); @@ -453,12 +434,12 @@ void NCVar::get_attr(const std::string& name, std::vector& values) const } //Uncomment for parallel NetCDF -void NCVar::par_access(const int cmode) const +void NCVar::par_access (const int cmode) const { check_nc_error(nc_var_par_access(ncid, varid, cmode)); } -std::string NCGroup::name() const +std::string NCGroup::name () const { size_t nlen; std::vector grpname; @@ -468,7 +449,7 @@ std::string NCGroup::name() const return std::string{grpname.begin(), grpname.end()}; } -std::string NCGroup::full_name() const +std::string NCGroup::full_name () const { size_t nlen; std::vector grpname; @@ -478,45 +459,44 @@ std::string NCGroup::full_name() const return std::string{grpname.begin(), grpname.end()}; } -NCGroup NCGroup::def_group(const std::string& name) const +NCGroup NCGroup::def_group (const std::string& name) const { int newid; check_nc_error(nc_def_grp(ncid, name.data(), &newid)); return NCGroup(newid, this); } -NCGroup NCGroup::group(const std::string& name) const +NCGroup NCGroup::group (const std::string& name) const { int newid; check_nc_error(nc_inq_ncid(ncid, name.data(), &newid)); return NCGroup(newid, this); } -NCDim NCGroup::dim(const std::string& name) const +NCDim NCGroup::dim (const std::string& name) const { int newid; check_nc_error(nc_inq_dimid(ncid, name.data(), &newid)); return NCDim{ncid, newid}; } -NCDim NCGroup::def_dim(const std::string& name, const size_t len) const +NCDim NCGroup::def_dim (const std::string& name, const size_t len) const { int newid; check_nc_error(nc_def_dim(ncid, name.data(), len, &newid)); return NCDim{ncid, newid}; } -NCVar NCGroup::def_scalar(const std::string& name, const nc_type dtype) const +NCVar NCGroup::def_scalar (const std::string& name, const nc_type dtype) const { int newid; check_nc_error(nc_def_var(ncid, name.data(), dtype, 0, nullptr, &newid)); return NCVar{ncid, newid}; } -NCVar NCGroup::def_array( - const std::string& name, - const nc_type dtype, - const std::vector& dnames) const +NCVar NCGroup::def_array (const std::string& name, + const nc_type dtype, + const std::vector& dnames) const { int newid; int ndims = dnames.size(); @@ -528,28 +508,28 @@ NCVar NCGroup::def_array( return NCVar{ncid, newid}; } -NCVar NCGroup::var(const std::string& name) const +NCVar NCGroup::var (const std::string& name) const { int varid; check_nc_error(nc_inq_varid(ncid, name.data(), &varid)); return NCVar{ncid, varid}; } -int NCGroup::num_groups() const +int NCGroup::num_groups () const { int ngrps; check_nc_error(nc_inq_grps(ncid, &ngrps, nullptr)); return ngrps; } -int NCGroup::num_dimensions() const +int NCGroup::num_dimensions () const { int ndims; check_nc_error(nc_inq(ncid, &ndims, nullptr, nullptr, nullptr)); return ndims; } -int NCGroup::num_attributes() const +int NCGroup::num_attributes () const { int nattrs; check_nc_error(nc_inq(ncid, nullptr, nullptr, &nattrs, nullptr)); @@ -563,25 +543,25 @@ int NCGroup::num_variables() const return nvars; } -bool NCGroup::has_group(const std::string& name) const +bool NCGroup::has_group (const std::string& name) const { int ierr = nc_inq_ncid(ncid, name.data(), nullptr); return (ierr == NC_NOERR); } -bool NCGroup::has_dim(const std::string& name) const +bool NCGroup::has_dim (const std::string& name) const { int ierr = nc_inq_dimid(ncid, name.data(), nullptr); return (ierr == NC_NOERR); } -bool NCGroup::has_var(const std::string& name) const +bool NCGroup::has_var (const std::string& name) const { int ierr = nc_inq_varid(ncid, name.data(), nullptr); return (ierr == NC_NOERR); } -bool NCGroup::has_attr(const std::string& name) const +bool NCGroup::has_attr (const std::string& name) const { int ierr; size_t lenp; @@ -589,34 +569,31 @@ bool NCGroup::has_attr(const std::string& name) const return (ierr == NC_NOERR); } -void NCGroup::put_attr(const std::string& name, const std::string& value) const +void NCGroup::put_attr (const std::string& name, const std::string& value) const { check_nc_error(nc_put_att_text( ncid, NC_GLOBAL, name.data(), value.size(), value.data())); } -void NCGroup::put_attr( - const std::string& name, const std::vector& value) const +void NCGroup::put_attr (const std::string& name, const std::vector& value) const { check_nc_error(nc_put_att_double( ncid, NC_GLOBAL, name.data(), NC_DOUBLE, value.size(), value.data())); } -void NCGroup::put_attr( - const std::string& name, const std::vector& value) const +void NCGroup::put_attr (const std::string& name, const std::vector& value) const { check_nc_error(nc_put_att_float( ncid, NC_GLOBAL, name.data(), NC_FLOAT, value.size(), value.data())); } -void NCGroup::put_attr( - const std::string& name, const std::vector& value) const +void NCGroup::put_attr (const std::string& name, const std::vector& value) const { check_nc_error(nc_put_att_int( ncid, NC_GLOBAL, name.data(), NC_INT, value.size(), value.data())); } -std::string NCGroup::get_attr(const std::string& name) const +std::string NCGroup::get_attr (const std::string& name) const { size_t lenp; std::vector aval; @@ -626,8 +603,7 @@ std::string NCGroup::get_attr(const std::string& name) const return std::string{aval.begin(), aval.end()}; } -void NCGroup::get_attr( - const std::string& name, std::vector& values) const +void NCGroup::get_attr (const std::string& name, std::vector& values) const { size_t lenp; check_nc_error(nc_inq_attlen(ncid, NC_GLOBAL, name.data(), &lenp)); @@ -636,8 +612,7 @@ void NCGroup::get_attr( nc_get_att_double(ncid, NC_GLOBAL, name.data(), values.data())); } -void NCGroup::get_attr( - const std::string& name, std::vector& values) const +void NCGroup::get_attr (const std::string& name, std::vector& values) const { size_t lenp; check_nc_error(nc_inq_attlen(ncid, NC_GLOBAL, name.data(), &lenp)); @@ -654,7 +629,7 @@ void NCGroup::get_attr(const std::string& name, std::vector& values) const check_nc_error(nc_get_att_int(ncid, NC_GLOBAL, name.data(), values.data())); } -std::vector NCGroup::all_groups() const +std::vector NCGroup::all_groups () const { std::vector grps; int ngrps = num_groups(); @@ -669,7 +644,7 @@ std::vector NCGroup::all_groups() const return grps; } -std::vector NCGroup::all_dims() const +std::vector NCGroup::all_dims () const { std::vector adims; int ndims = num_dimensions(); @@ -680,7 +655,7 @@ std::vector NCGroup::all_dims() const return adims; } -std::vector NCGroup::all_vars() const +std::vector NCGroup::all_vars () const { std::vector avars; int nvars = num_variables(); @@ -691,7 +666,7 @@ std::vector NCGroup::all_vars() const return avars; } -void NCGroup::enter_def_mode() const +void NCGroup::enter_def_mode () const { int ierr; ierr = nc_redef(ncid); @@ -702,24 +677,23 @@ void NCGroup::enter_def_mode() const check_nc_error(ierr); } -void NCGroup::exit_def_mode() const { check_nc_error(nc_enddef(ncid)); } +void NCGroup::exit_def_mode () const { check_nc_error(nc_enddef(ncid)); } -NCFile NCFile::create(const std::string& name, const int cmode) +NCFile NCFile::create (const std::string& name, const int cmode) { int ncid; check_nc_error(nc_create(name.data(), cmode, &ncid)); return NCFile(ncid); } -NCFile NCFile::open(const std::string& name, const int cmode) +NCFile NCFile::open (const std::string& name, const int cmode) { int ncid; check_nc_error(nc_open(name.data(), cmode, &ncid)); return NCFile(ncid); } //Uncomment for parallel NetCDF -NCFile NCFile::create_par( - const std::string& name, const int cmode, MPI_Comm comm, MPI_Info info) +NCFile NCFile::create_par (const std::string& name, const int cmode, MPI_Comm comm, MPI_Info info) { int ncid; check_nc_error(nc_create_par(name.data(), cmode, comm, info, &ncid)); @@ -727,20 +701,19 @@ NCFile NCFile::create_par( } //Uncomment for parallel NetCDF -NCFile NCFile::open_par( - const std::string& name, const int cmode, MPI_Comm comm, MPI_Info info) +NCFile NCFile::open_par (const std::string& name, const int cmode, MPI_Comm comm, MPI_Info info) { int ncid; check_nc_error(nc_open_par(name.data(), cmode, comm, info, &ncid)); return NCFile(ncid); } -NCFile::~NCFile() +NCFile::~NCFile () { if (is_open) check_nc_error(nc_close(ncid)); } -void NCFile::close() +void NCFile::close () { is_open = false; check_nc_error(nc_close(ncid)); diff --git a/Source/IO/NCMultiFabFile.cpp b/Source/IO/NCMultiFabFile.cpp index 61fce6227..3235380fe 100644 --- a/Source/IO/NCMultiFabFile.cpp +++ b/Source/IO/NCMultiFabFile.cpp @@ -20,10 +20,10 @@ using namespace amrex; void -ERF::ReadNCMultiFab(FabArray &mf, - const std::string &mf_name, - int /*coordinatorProc*/, - int /*allow_empty_mf*/) { +ERF::ReadNCMultiFab (FabArray &mf, + const std::string &mf_name, + int /*coordinatorProc*/, + int /*allow_empty_mf*/) { const std::string& FullPath = amrex::Concatenate(check_file,istep[0],5); diff --git a/Source/IO/NCPlotFile.cpp b/Source/IO/NCPlotFile.cpp index bbb094d71..807a1401f 100644 --- a/Source/IO/NCPlotFile.cpp +++ b/Source/IO/NCPlotFile.cpp @@ -19,10 +19,10 @@ using namespace amrex; void -ERF::writeNCPlotFile(int lev, int which_subdomain, const std::string& dir, - const Vector &plotMF, - const Vector &plot_var_names, - const Vector& level_steps, const Real time) const +ERF::writeNCPlotFile (int lev, int which_subdomain, const std::string& dir, + const Vector &plotMF, + const Vector &plot_var_names, + const Vector& level_steps, const Real time) const { // get the processor number int iproc = amrex::ParallelContext::MyProcAll(); @@ -132,6 +132,7 @@ ERF::writeNCPlotFile(int lev, int which_subdomain, const std::string& dir, ncf.put_attr("number_variables", std::vector{n_data_items}); ncf.put_attr("space_dimension", std::vector{AMREX_SPACEDIM}); ncf.put_attr("current_time", std::vector{time}); + ncf.put_attr("start_time", std::vector{start_bdy_time}); ncf.put_attr("CurrentLevel", std::vector{flev}); Real dx[AMREX_SPACEDIM]; diff --git a/Source/IO/NCWpsFile.H b/Source/IO/NCWpsFile.H index 097fa7502..7befc00cd 100644 --- a/Source/IO/NCWpsFile.H +++ b/Source/IO/NCWpsFile.H @@ -7,14 +7,19 @@ #include #include "AMReX_FArrayBox.H" +#include "AMReX_IArrayBox.H" #include "NCInterface.H" using PlaneVector = amrex::Vector; /* // Read from metgrid - NetCDF variables of dimensions Time_BT_SN_WE: "UU", "VV", "TT", "RH", "PRES" - NetCDF variables of dimensions Time_SN_WE : "HGT", "MAPFAC_U", "MAPFAC_V", "MAPFAC_M" + NetCDF variables of dimensions Time_BT_SN_WE: "UU", "VV", "TT", "RH", "PRES", "GHT" + NetCDF variables of dimensions Time_SN_WE : "HGT", "MAPFAC_U", "MAPFAC_V", "MAPFAC_M", "PSFC" + NetCDF global attributes of type int : "FLAG_PSFC", "FLAG_MAPFAC_U", "FLAG_MAPFAC_V", "FLAG_MAPFAC_M", + "FLAG_HGT_M", "WEST-EAST_GRID_DIMENSION", "SOUTH-NORTH_GRID_DIMENSION" + NetCDF global attributes of type string : "SIMULATION_START_DATE" + NetCDF global attributes of type real : "DX", "DY" // Read from wrfbdy NetCDF variables of dimensions Time_BdyWidth_BT_SN : "U_BXS", "U_BXE", "V_BXS", "V_BXE" etc. @@ -33,18 +38,6 @@ enum class NC_Data_Dims_Type { Time_BdyWidth_WE }; -void BuildFABsFromNetCDFFile(const amrex::Box& domain, - const std::string &fname, - amrex::Vector nc_var_names, - amrex::Vector NC_dim_types, - amrex::Vector fab_vars); - -int BuildFABsFromWRFBdyFile(const std::string &fname, - amrex::Vector>& bdy_data_xlo, - amrex::Vector>& bdy_data_xhi, - amrex::Vector>& bdy_data_ylo, - amrex::Vector>& bdy_data_yhi); - // // NDArray is the datatype designed to hold any data, including scalars, multidimensional // arrays, that read from the NetCDF file. @@ -65,10 +58,10 @@ struct NDArray } // default constructor - NDArray() : name{"null"}, ref_counted{1}, owned{false}, data{nullptr} {} + NDArray () : name{"null"}, ref_counted{1}, owned{false}, data{nullptr} {} // copy constructor - NDArray(const NDArray& array) { + NDArray (const NDArray& array) { name = array.name; shape = array.shape; data = array.data; @@ -87,29 +80,29 @@ struct NDArray } // destructor - ~NDArray() { + ~NDArray () { ref_counted.fetch_sub(1, std::memory_order_acq_rel); if (ref_counted == 1 && owned) delete[] data; } // get the data pointer - decltype(auto) get_data() { + decltype(auto) get_data () { ref_counted.fetch_add(1, std::memory_order_relaxed); return data; } // get the variable name - std::string get_vname() { + std::string get_vname () { return name; } // get the variable data shape - std::vector get_vshape() { + std::vector get_vshape () { return shape; } // return the total number of data - size_t ndim() { + size_t ndim () { size_t num = 1; int isize = static_cast(shape.size()); for (auto i=0; i vshape) { + void set_vshape (std::vector vshape) { shape = vshape; } @@ -129,9 +122,15 @@ struct NDArray DType* data; }; +int BuildFABsFromWRFBdyFile (const std::string &fname, + amrex::Vector>& bdy_data_xlo, + amrex::Vector>& bdy_data_xhi, + amrex::Vector>& bdy_data_ylo, + amrex::Vector>& bdy_data_yhi); + template -void ReadNetCDFFile(const std::string& fname, amrex::Vector names, - amrex::Vector >& arrays) +void ReadNetCDFFile (const std::string& fname, amrex::Vector names, + amrex::Vector >& arrays) { AMREX_ASSERT(arrays.size() == names.size()); @@ -160,7 +159,12 @@ void ReadNetCDFFile(const std::string& fname, amrex::Vector names, if (vname_to_read.substr(0,2) == "R_") { vname_to_read = names[n+4]; // This allows us to read "T" instead -- we will over-write this later } - // amrex::Print() << "About to read " << vname_to_read << " while filling the array for " << vname_to_write << std::endl; + + /* + amrex::AllPrint() << "About to read " << vname_to_read + << " while filling the array for " << vname_to_write << std::endl; + */ + std::vector shape = ncf.var(vname_to_read).shape(); arrays[n] = NDArray(vname_to_read,shape); DType* dataPtr = arrays[n].get_data(); @@ -180,4 +184,142 @@ void ReadNetCDFFile(const std::string& fname, amrex::Vector names, ncf.close(); } } + +/** + * Helper function for reading data from NetCDF file into a + * provided FAB. + * + * @param iv Index for which variable we are going to fill + * @param nc_arrays Arrays of data from NetCDF file + * @param var_name Variable name + * @param NC_dim_type Dimension type for the variable as stored in the NetCDF file + * @param temp FAB where we store the variable data from the NetCDF Arrays + */ +template +void +fill_fab_from_arrays (int iv, + amrex::Vector>& nc_arrays, + const std::string& var_name, + NC_Data_Dims_Type& NC_dim_type, + FAB& temp) +{ + int ns1, ns2, ns3; + if (NC_dim_type == NC_Data_Dims_Type::Time_BT) { + ns1 = nc_arrays[iv].get_vshape()[1]; + ns2 = 1; + ns3 = 1; + // amrex::Print() << "TYPE BT " << ns1 << std::endl; + } else if (NC_dim_type == NC_Data_Dims_Type::Time_SN_WE) { + ns1 = 1; + ns2 = nc_arrays[iv].get_vshape()[1]; + ns3 = nc_arrays[iv].get_vshape()[2]; + // amrex::Print() << "TYPE SN WE " << ns2 << " " << ns3 << std::endl; + } else if (NC_dim_type == NC_Data_Dims_Type::Time_BT_SN_WE) { + ns1 = nc_arrays[iv].get_vshape()[1]; + ns2 = nc_arrays[iv].get_vshape()[2]; + ns3 = nc_arrays[iv].get_vshape()[3]; + // amrex::Print() << "TYPE BT SN WE " << ns1 << " " << ns2 << " " << ns3 << std::endl; + } else { + amrex::Abort("Dont know this NC_Data_Dims_Type"); + } + + // TODO: The box will only start at (0,0,0) at level 0 -- we need to generalize this + amrex::Box my_box(amrex::IntVect(0,0,0), amrex::IntVect(ns3-1,ns2-1,ns1-1)); + // amrex::Print() <<" MY BOX " << my_box << std::endl; + + if (var_name == "U" || var_name == "UU" || + var_name == "MAPFAC_U" || var_name == "MAPFAC_UY") my_box.setType(amrex::IndexType(amrex::IntVect(1,0,0))); + if (var_name == "V" || var_name == "VV" || + var_name == "MAPFAC_V" || var_name == "MAPFAC_VY") my_box.setType(amrex::IndexType(amrex::IntVect(0,1,0))); + if (var_name == "W" || var_name == "WW") my_box.setType(amrex::IndexType(amrex::IntVect(0,0,1))); + +#ifdef AMREX_USE_GPU + // Make sure temp lives on CPU since nc_arrays lives on CPU only + temp.resize(my_box,1,amrex::The_Pinned_Arena()); +#else + temp.resize(my_box,1); +#endif + amrex::Array4 fab_arr = temp.array(); + + int ioff = temp.box().smallEnd()[0]; + int joff = temp.box().smallEnd()[1]; + + auto num_pts = my_box.numPts(); + + // amrex::Print() <<" ns1 * ns2 * ns3 " << ns1 * ns2 * ns3 << std::endl; + // amrex::Print() <<" NUMPTS " << num_pts << std::endl; + + for (int n(0); n < num_pts; ++n) { + int k = n / (ns2*ns3); + int j = (n - k*(ns2*ns3)) / ns3 + joff; + int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff; + fab_arr(i,j,k,0) = static_cast(*(nc_arrays[iv].get_data()+n)); + } +} + +/** + * Function to read NetCDF variables and fill the corresponding Array4's + * + * @param fname Name of the NetCDF file to be read + * @param nc_var_names Variable names in the NetCDF file + * @param NC_dim_types NetCDF data dimension types + * @param fab_vars Fab data we are to fill + */ +template +void +BuildFABsFromNetCDFFile (const amrex::Box& domain, + const std::string &fname, + amrex::Vector nc_var_names, + amrex::Vector NC_dim_types, + amrex::Vector fab_vars) +{ + int ioproc = amrex::ParallelDescriptor::IOProcessorNumber(); // I/O rank + + amrex::Vector> nc_arrays(nc_var_names.size()); + + if (amrex::ParallelDescriptor::IOProcessor()) + { + ReadNetCDFFile(fname, nc_var_names, nc_arrays); + } + + for (int iv = 0; iv < nc_var_names.size(); iv++) + { + FAB tmp; + if (amrex::ParallelDescriptor::IOProcessor()) { + fill_fab_from_arrays(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp); + } + + int ncomp = tmp.nComp(); + amrex::Box box = tmp.box(); + + amrex::ParallelDescriptor::Bcast(&box, 1, ioproc); + amrex::ParallelDescriptor::Bcast(&ncomp, 1, ioproc); + + if (!amrex::ParallelDescriptor::IOProcessor()) { +#ifdef AMREX_USE_GPU + tmp.resize(box,ncomp,amrex::The_Pinned_Arena()); +#else + tmp.resize(box,ncomp); +#endif + } + + amrex::ParallelDescriptor::Bcast(tmp.dataPtr(), tmp.size(), ioproc); + + // Shift box by the domain lower corner + amrex::Box fab_bx = tmp.box(); + amrex::Dim3 dom_lb = lbound(domain); + fab_bx += amrex::IntVect(dom_lb.x,dom_lb.y,dom_lb.z); + // fab_vars points to data on device + fab_vars[iv]->resize(fab_bx,1); +#ifdef AMREX_USE_GPU + amrex::Gpu::copy(amrex::Gpu::hostToDevice, + tmp.dataPtr(), tmp.dataPtr() + tmp.size(), + fab_vars[iv]->dataPtr()); +#else + // Provided by BaseFab inheritance through FArrayBox + fab_vars[iv]->copy(tmp,tmp.box(),0,fab_bx,0,1); +#endif + } +} + #endif diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index 5763a585a..41729e723 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -8,7 +8,7 @@ using namespace amrex; template -bool containerHasElement(const V& iterable, const T& query) { +bool containerHasElement (const V& iterable, const T& query) { return std::find(iterable.begin(), iterable.end(), query) != iterable.end(); } @@ -78,7 +78,7 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector // set plotfile variable names Vector -ERF::PlotFileVarNames ( Vector plot_var_names ) +ERF::PlotFileVarNames (Vector plot_var_names ) { Vector names; diff --git a/Source/IO/ReadFromMetgrid.cpp b/Source/IO/ReadFromMetgrid.cpp index 094469275..80748f083 100644 --- a/Source/IO/ReadFromMetgrid.cpp +++ b/Source/IO/ReadFromMetgrid.cpp @@ -1,42 +1,123 @@ -#include "NCWpsFile.H" -#include "AMReX_FArrayBox.H" +#include +#include +#include using namespace amrex; #ifdef ERF_USE_NETCDF + +// Converts UTC time string to a time_t value in UTC. +std::time_t +getEpochTime_metgrid (const std::string& dateTime, const std::string& dateTimeFormat) +{ + // Create a stream which we will use to parse the string, + // which we provide to constructor of stream to fill the buffer. + std::istringstream ss{ dateTime }; + + // Create a tm object to store the parsed date and time. + std::tm tmTime; + memset(&tmTime, 0, sizeof(tmTime)); + + // Now we read from buffer using get_time manipulator + // and formatting the input appropriately. + strptime(dateTime.c_str(), dateTimeFormat.c_str(), &tmTime); + + // Convert the tm structure to time_t value and return. + // Here we use timegm since the output should be relative to UTC. + auto epoch = timegm(&tmTime); + // Print() << "Time Stamp: "<< std::put_time(&tmTime, "%c") + // << " , Epoch: " << epoch << std::endl; + + return epoch; +} + void -read_from_metgrid(int lev, - const Box& domain, - const std::string& fname, - FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab, - FArrayBox& NC_temp_fab, FArrayBox& NC_rhum_fab, - FArrayBox& NC_pres_fab, FArrayBox& NC_hgt_fab, - FArrayBox& NC_msfu_fab, FArrayBox& NC_msfv_fab, - FArrayBox& NC_msfm_fab) +read_from_metgrid (int lev, const Box& domain, const std::string& fname, + std::string& NC_dateTime, Real& NC_epochTime, + int& flag_psfc, int& flag_msfu, int& flag_msfv, int& flag_msfm, + int& flag_hgt, int& flag_sst, int& flag_lmask, + int& NC_nx, int& NC_ny, + Real& NC_dx, Real& NC_dy, + FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab, + FArrayBox& NC_temp_fab, FArrayBox& NC_rhum_fab, + FArrayBox& NC_pres_fab, FArrayBox& NC_ght_fab, + FArrayBox& NC_hgt_fab, FArrayBox& NC_psfc_fab, + FArrayBox& NC_msfu_fab, FArrayBox& NC_msfv_fab, + FArrayBox& NC_msfm_fab, FArrayBox& NC_sst_fab, + IArrayBox& NC_lmask_iab) { amrex::Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; - int ncomp = 1; + if (amrex::ParallelDescriptor::IOProcessor()) { + auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); + { // Global Attributes (int) + std::vector attr; + ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0]; + ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0]; + ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0]; + ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0]; + ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0]; + ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; + ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0]; + ncf.get_attr("WEST-EAST_GRID_DIMENSION", attr); NC_nx = attr[0]; + ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0]; + } + { // Global Attributes (string) + NC_dateTime = ncf.get_attr("SIMULATION_START_DATE")+"UTC"; + const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S%Z"; + NC_epochTime = getEpochTime_metgrid(NC_dateTime, dateTimeFormat); + } + { // Global Attributes (Real) + std::vector attr; + ncf.get_attr("DX", attr); NC_dx = attr[0]; + ncf.get_attr("DY", attr); NC_dy = attr[0]; + } + ncf.close(); + } + int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank + ParallelDescriptor::Bcast(&flag_psfc, 1, ioproc); + ParallelDescriptor::Bcast(&flag_msfu, 1, ioproc); + ParallelDescriptor::Bcast(&flag_msfv, 1, ioproc); + ParallelDescriptor::Bcast(&flag_msfm, 1, ioproc); + ParallelDescriptor::Bcast(&flag_hgt, 1, ioproc); + ParallelDescriptor::Bcast(&flag_sst, 1, ioproc); + ParallelDescriptor::Bcast(&flag_lmask, 1, ioproc); + ParallelDescriptor::Bcast(&NC_nx, 1, ioproc); + ParallelDescriptor::Bcast(&NC_ny, 1, ioproc); + ParallelDescriptor::Bcast(&NC_epochTime, 1, ioproc); + ParallelDescriptor::Bcast(&NC_dx, 1, ioproc); + ParallelDescriptor::Bcast(&NC_dy, 1, ioproc); + + Vector NC_fabs; + Vector NC_iabs; + Vector NC_fnames; + Vector NC_inames; + Vector NC_fdim_types; + Vector NC_idim_types; - Vector NC_fabs; - Vector NC_names; - Vector NC_dim_types; + NC_fabs.push_back(&NC_xvel_fab); NC_fnames.push_back("UU"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); + NC_fabs.push_back(&NC_yvel_fab); NC_fnames.push_back("VV"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); + NC_fabs.push_back(&NC_temp_fab); NC_fnames.push_back("TT"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); + NC_fabs.push_back(&NC_rhum_fab); NC_fnames.push_back("RH"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); + NC_fabs.push_back(&NC_pres_fab); NC_fnames.push_back("PRES"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); + NC_fabs.push_back(&NC_ght_fab); NC_fnames.push_back("GHT"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); - NC_fabs.push_back(&NC_xvel_fab); NC_names.push_back("UU"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); - NC_fabs.push_back(&NC_yvel_fab); NC_names.push_back("VV"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); - NC_fabs.push_back(&NC_temp_fab); NC_names.push_back("TT"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); - NC_fabs.push_back(&NC_rhum_fab); NC_names.push_back("RH"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); - NC_fabs.push_back(&NC_pres_fab); NC_names.push_back("PRES"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); + if (flag_psfc) { NC_fabs.push_back(&NC_psfc_fab); NC_fnames.push_back("PSFC"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_msfu) { NC_fabs.push_back(&NC_msfu_fab); NC_fnames.push_back("MAPFAC_U"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_msfv) { NC_fabs.push_back(&NC_msfv_fab); NC_fnames.push_back("MAPFAC_V"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_msfm) { NC_fabs.push_back(&NC_msfm_fab); NC_fnames.push_back("MAPFAC_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_hgt) { NC_fabs.push_back(&NC_hgt_fab); NC_fnames.push_back("HGT_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_sst) { NC_fabs.push_back(&NC_sst_fab); NC_fnames.push_back("SST"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } - NC_fabs.push_back(&NC_hgt_fab); NC_names.push_back("HGT_M"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); - NC_fabs.push_back(&NC_msfu_fab); NC_names.push_back("MAPFAC_U"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); - NC_fabs.push_back(&NC_msfv_fab); NC_names.push_back("MAPFAC_V"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); - NC_fabs.push_back(&NC_msfm_fab); NC_names.push_back("MAPFAC_M"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); + if (flag_lmask) { NC_iabs.push_back(&NC_lmask_iab); NC_inames.push_back("LANDMASK"); NC_idim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } // Read the netcdf file and fill these FABs amrex::Print() << "Building initial FABS from file " << fname << std::endl; - BuildFABsFromNetCDFFile(domain, fname, NC_names, NC_dim_types, NC_fabs); + BuildFABsFromNetCDFFile(domain, fname, NC_fnames, NC_fdim_types, NC_fabs); + // Read the netcdf file and fill these IABs + amrex::Print() << "Building initial IABS from file " << fname << std::endl; + BuildFABsFromNetCDFFile(domain, fname, NC_inames, NC_idim_types, NC_iabs); // TODO: FIND OUT IF WE NEED TO DIVIDE VELS BY MAPFAC // diff --git a/Source/IO/ReadFromWRFBdy.cpp b/Source/IO/ReadFromWRFBdy.cpp index f111ea490..ad78f6a89 100644 --- a/Source/IO/ReadFromWRFBdy.cpp +++ b/Source/IO/ReadFromWRFBdy.cpp @@ -10,7 +10,6 @@ #include "DataStruct.H" #include "NCInterface.H" -#include "NCWpsFile.H" #include "AMReX_FArrayBox.H" #include "AMReX_Print.H" @@ -28,7 +27,7 @@ namespace WRFBdyTypes { } // Converts UTC time string to a time_t value. -std::time_t getEpochTime(const std::string& dateTime, const std::string& dateTimeFormat) +std::time_t getEpochTime (const std::string& dateTime, const std::string& dateTimeFormat) { // Create a stream which we will use to parse the string, // which we provide to constructor of stream to fill the buffer. @@ -51,12 +50,12 @@ std::time_t getEpochTime(const std::string& dateTime, const std::string& dateTim } Real -read_from_wrfbdy(const std::string& nc_bdy_file, const Box& domain, - Vector>& bdy_data_xlo, - Vector>& bdy_data_xhi, - Vector>& bdy_data_ylo, - Vector>& bdy_data_yhi, - int& width, Real& start_bdy_time) +read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, + Vector>& bdy_data_xlo, + Vector>& bdy_data_xhi, + Vector>& bdy_data_ylo, + Vector>& bdy_data_yhi, + int& width, Real& start_bdy_time) { amrex::Print() << "Loading boundary data from NetCDF file " << std::endl; @@ -96,10 +95,11 @@ read_from_wrfbdy(const std::string& nc_bdy_file, const Box& domain, auto epochTime = getEpochTime(date, dateTimeFormat); epochTimes.push_back(epochTime); - if (nt == 1) + if (nt == 1) { timeInterval = epochTimes[1] - epochTimes[0]; - else if (nt >= 1) + } else if (nt >= 1) { AMREX_ALWAYS_ASSERT(epochTimes[nt] - epochTimes[nt-1] == timeInterval); + } } start_bdy_time = epochTimes[0]; } @@ -547,15 +547,15 @@ read_from_wrfbdy(const std::string& nc_bdy_file, const Box& domain, } void -convert_wrfbdy_data(int which, const Box& domain, Vector>& bdy_data, - const FArrayBox& NC_MUB_fab, - const FArrayBox& NC_MSFU_fab, const FArrayBox& NC_MSFV_fab, - const FArrayBox& NC_MSFM_fab, - const FArrayBox& NC_PH_fab, const FArrayBox& NC_PHB_fab, - const FArrayBox& NC_C1H_fab, const FArrayBox& NC_C2H_fab, - const FArrayBox& NC_RDNW_fab, - const FArrayBox& NC_xvel_fab, const FArrayBox& NC_yvel_fab, - const FArrayBox& NC_rho_fab, const FArrayBox& NC_rhotheta_fab) +convert_wrfbdy_data (int which, const Box& domain, Vector>& bdy_data, + const FArrayBox& NC_MUB_fab, + const FArrayBox& NC_MSFU_fab, const FArrayBox& NC_MSFV_fab, + const FArrayBox& NC_MSFM_fab, + const FArrayBox& NC_PH_fab, const FArrayBox& NC_PHB_fab, + const FArrayBox& NC_C1H_fab, const FArrayBox& NC_C2H_fab, + const FArrayBox& NC_RDNW_fab, + const FArrayBox& NC_xvel_fab, const FArrayBox& NC_yvel_fab, + const FArrayBox& NC_rho_fab, const FArrayBox& NC_rhotheta_fab) { // These were filled from wrfinput Array4 c1h_arr = NC_C1H_fab.const_array(); diff --git a/Source/IO/ReadFromWRFInput.cpp b/Source/IO/ReadFromWRFInput.cpp index f2b69ea85..f930af6bb 100644 --- a/Source/IO/ReadFromWRFInput.cpp +++ b/Source/IO/ReadFromWRFInput.cpp @@ -5,25 +5,25 @@ using namespace amrex; #ifdef ERF_USE_NETCDF void -read_from_wrfinput(int lev, - const Box& domain, - const std::string& fname, - FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab, - FArrayBox& NC_zvel_fab, FArrayBox& NC_rho_fab, - FArrayBox& NC_rhop_fab, FArrayBox& NC_rhotheta_fab, - FArrayBox& NC_MUB_fab , - FArrayBox& NC_MSFU_fab, FArrayBox& NC_MSFV_fab, - FArrayBox& NC_MSFM_fab, FArrayBox& NC_SST_fab, - FArrayBox& NC_C1H_fab , FArrayBox& NC_C2H_fab, - FArrayBox& NC_RDNW_fab, +read_from_wrfinput (int lev, + const Box& domain, + const std::string& fname, + FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab, + FArrayBox& NC_zvel_fab, FArrayBox& NC_rho_fab, + FArrayBox& NC_rhop_fab, FArrayBox& NC_rhotheta_fab, + FArrayBox& NC_MUB_fab , + FArrayBox& NC_MSFU_fab, FArrayBox& NC_MSFV_fab, + FArrayBox& NC_MSFM_fab, FArrayBox& NC_SST_fab, + FArrayBox& NC_C1H_fab , FArrayBox& NC_C2H_fab, + FArrayBox& NC_RDNW_fab, #if defined(ERF_USE_MOISTURE) - FArrayBox& NC_QVAPOR_fab, - FArrayBox& NC_QCLOUD_fab, - FArrayBox& NC_QRAIN_fab, + FArrayBox& NC_QVAPOR_fab, + FArrayBox& NC_QCLOUD_fab, + FArrayBox& NC_QRAIN_fab, #elif defined(ERF_USE_WARM_NO_PRECIP) #endif - FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab, - FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab) + FArrayBox& NC_PH_fab , FArrayBox& NC_PHB_fab, + FArrayBox& NC_ALB_fab , FArrayBox& NC_PB_fab) { amrex::Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; @@ -60,7 +60,7 @@ read_from_wrfinput(int lev, // Read the netcdf file and fill these FABs amrex::Print() << "Building initial FABS from file " << fname << std::endl; - BuildFABsFromNetCDFFile(domain, fname, NC_names, NC_dim_types, NC_fabs); + BuildFABsFromNetCDFFile(domain, fname, NC_names, NC_dim_types, NC_fabs); // // Convert the velocities using the map factors diff --git a/Source/IO/writeJobInfo.cpp b/Source/IO/writeJobInfo.cpp index 0d461ed9d..6c02d1d75 100644 --- a/Source/IO/writeJobInfo.cpp +++ b/Source/IO/writeJobInfo.cpp @@ -4,7 +4,7 @@ extern std::string inputs_name; void -ERF::writeJobInfo(const std::string& dir) const +ERF::writeJobInfo (const std::string& dir) const { // job_info file with details about the run std::ofstream jobInfoFile; @@ -132,7 +132,7 @@ ERF::writeJobInfo(const std::string& dir) const } void -ERF::writeBuildInfo(std::ostream& os) +ERF::writeBuildInfo (std::ostream& os) { std::string PrettyLine = std::string(78, '=') + "\n"; std::string OtherLine = std::string(78, '-') + "\n"; diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index 306e7a8a6..9c90436a1 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -75,6 +75,17 @@ namespace WRFBdyVars { }; } +namespace MetGridBdyVars { + enum { + U = 0, + V = 1, + R = 2, + T = 3, + QV, + NumTypes + }; +} + namespace Vars { enum { cons = 0, diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 26453d2db..6037c0de7 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -2,52 +2,10 @@ * \file ERF_init_from_metgrid.cpp */ -#include -#include -#include -#include -#include +#include using namespace amrex; -void -read_from_metgrid (int lev, const Box& domain, const std::string& fname, - FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab, - FArrayBox& NC_temp_fab, FArrayBox& NC_rhum_fab, - FArrayBox& NC_pres_fab, FArrayBox& NC_hgt_fab, - FArrayBox& NC_msfu_fab, FArrayBox& NC_msfv_fab, - FArrayBox& NC_msfm_fab); -void -interpolate_column (int i, int j, int src_comp, int dest_comp, - const Array4& orig_z, const Array4& orig_data, - const Array4& new_z, const Array4& new_data); - -void -init_terrain_from_metgrid (int lev, FArrayBox& z_phys_nd_fab, - const Vector& NC_hgt_fab); - -void -init_state_from_metgrid (int lev, FArrayBox& state_fab, - FArrayBox& x_vel_fab, FArrayBox& y_vel_fab, - FArrayBox& z_vel_fab, FArrayBox& z_phys_nd_fab, - const Vector& NC_hgt_fab, - const Vector& NC_xvel_fab, - const Vector& NC_yvel_fab, - const Vector& NC_zvel_fab, - const Vector& NC_rho_fab, - const Vector& NC_rhotheta_fab); -void -init_msfs_from_metgrid (int lev, FArrayBox& msfu_fab, - FArrayBox& msfv_fab, FArrayBox& msfm_fab, - const Vector& NC_MSFU_fab, - const Vector& NC_MSFV_fab, - const Vector& NC_MSFM_fab); -void -init_base_state_from_metgrid (int lev, const Box& valid_bx, Real l_rdOcp, - FArrayBox& p_hse, FArrayBox& pi_hse, FArrayBox& r_hse, - const Vector& NC_ALB_fab, - const Vector& NC_PB_fab); - #ifdef ERF_USE_NETCDF /** * Initializes ERF data using metgrid data supplied by an external NetCDF file. @@ -57,20 +15,18 @@ init_base_state_from_metgrid (int lev, const Box& valid_bx, Real l_rdOcp, void ERF::init_from_metgrid (int lev) { - // *** FArrayBox's at this level for holding the INITIAL data - Vector NC_xvel_fab ; NC_xvel_fab.resize(num_boxes_at_level[lev]); - Vector NC_yvel_fab ; NC_yvel_fab.resize(num_boxes_at_level[lev]); - Vector NC_temp_fab ; NC_temp_fab.resize(num_boxes_at_level[lev]); - Vector NC_rhum_fab ; NC_rhum_fab.resize(num_boxes_at_level[lev]); - Vector NC_pres_fab ; NC_pres_fab.resize(num_boxes_at_level[lev]); - - Vector NC_hgt_fab ; NC_hgt_fab.resize(num_boxes_at_level[lev]); - - Vector NC_MSFU_fab ; NC_MSFU_fab.resize(num_boxes_at_level[lev]); - Vector NC_MSFV_fab ; NC_MSFV_fab.resize(num_boxes_at_level[lev]); - Vector NC_MSFM_fab ; NC_MSFM_fab.resize(num_boxes_at_level[lev]); +#ifndef AMREX_USE_GPU +#if defined(ERF_USE_MOISTURE) + amrex::Print() << "Init with met_em with ERF_USE_MOISTURE" << std::endl; +#elif defined(ERF_USE_WARM_NO_PRECIP) + amrex::Print() << "Init with met_em with ERF_USE_WARM_NO_PRECIP" << std::endl; +#else + amrex::Print() << "Init with met_em without moisture" << std::endl; +#endif +#endif int nboxes = num_boxes_at_level[lev]; + int ntimes = num_files_at_level[lev]; if (nc_init_file.empty()) amrex::Error("NetCDF initialization file name must be provided via input"); @@ -78,13 +34,109 @@ ERF::init_from_metgrid (int lev) if (nc_init_file[lev].empty()) amrex::Error("NetCDF initialization file name must be provided via input"); - for (int idx = 0; idx < nboxes; idx++) - { - read_from_metgrid(lev, boxes_at_level[lev][idx], nc_init_file[lev][idx], - NC_xvel_fab[idx], NC_yvel_fab[idx], - NC_temp_fab[idx], NC_rhum_fab[idx], - NC_pres_fab[idx], NC_hgt_fab[idx], - NC_MSFU_fab[idx], NC_MSFV_fab[idx], NC_MSFM_fab[idx] ); + // At least two met_em files are necessary to calculate tendency terms. + AMREX_ALWAYS_ASSERT(ntimes >= 2); + + // Size the SST and LANDMASK + sst_lev[lev].resize(ntimes); + lmask_lev[lev].resize(ntimes); + + // *** FArrayBox's at this level for holding the metgrid data + Vector NC_xvel_fab; NC_xvel_fab.resize(ntimes); + Vector NC_yvel_fab; NC_yvel_fab.resize(ntimes); + Vector NC_temp_fab; NC_temp_fab.resize(ntimes); + Vector NC_rhum_fab; NC_rhum_fab.resize(ntimes); + Vector NC_pres_fab; NC_pres_fab.resize(ntimes); + Vector NC_ght_fab; NC_ght_fab.resize( ntimes); + Vector NC_hgt_fab; NC_hgt_fab.resize( ntimes); + Vector NC_psfc_fab; NC_psfc_fab.resize(ntimes); + Vector NC_MSFU_fab; NC_MSFU_fab.resize(ntimes); + Vector NC_MSFV_fab; NC_MSFV_fab.resize(ntimes); + Vector NC_MSFM_fab; NC_MSFM_fab.resize(ntimes); + Vector NC_sst_fab; NC_sst_fab.resize (ntimes); + + // *** IArrayBox's at this level for holding mask data + Vector NC_lmask_iab; NC_lmask_iab.resize(ntimes); + + // *** Variables at this level for holding metgrid file global attributes + Vector flag_psfc; flag_psfc.resize( ntimes); + Vector flag_msfu; flag_msfu.resize( ntimes); + Vector flag_msfv; flag_msfv.resize( ntimes); + Vector flag_msfm; flag_msfm.resize( ntimes); + Vector flag_hgt; flag_hgt.resize( ntimes); + Vector flag_sst; flag_sst.resize( ntimes); + Vector flag_lmask; flag_lmask.resize( ntimes); + Vector NC_nx; NC_nx.resize( ntimes); + Vector NC_ny; NC_ny.resize( ntimes); + Vector NC_dateTime; NC_dateTime.resize( ntimes); + Vector NC_epochTime; NC_epochTime.resize(ntimes); + Vector NC_dx; NC_dx.resize( ntimes); + Vector NC_dy; NC_dy.resize( ntimes); + + for (int it = 0; it < ntimes; it++) { +#ifndef AMREX_USE_GPU + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << ": nc_init_file[" << lev << "][" << it << "]\t" << nc_init_file[lev][it] << std::endl; +#endif + read_from_metgrid(lev, boxes_at_level[lev][0], nc_init_file[lev][it], + NC_dateTime[it], NC_epochTime[it], + flag_psfc[it], flag_msfu[it], flag_msfv[it], flag_msfm[it], + flag_hgt[it], flag_sst[it], flag_lmask[it], + NC_nx[it], NC_ny[it], NC_dx[it], NC_dy[it], + NC_xvel_fab[it], NC_yvel_fab[it], + NC_temp_fab[it], NC_rhum_fab[it], NC_pres_fab[it], + NC_ght_fab[it], NC_hgt_fab[it], NC_psfc_fab[it], + NC_MSFU_fab[it], NC_MSFV_fab[it], NC_MSFM_fab[it], + NC_sst_fab[it], NC_lmask_iab[it]); +#ifndef AMREX_USE_GPU + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_psfc \t" << flag_psfc[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfu \t" << flag_msfu[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfv \t" << flag_msfv[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_msfm \t" << flag_msfm[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": flag_hgt \t" << flag_hgt[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_nx \t" << NC_nx[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_ny \t" << NC_ny[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_dx \t" << NC_dx[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_dy \t" << NC_dy[it] << std::endl; + amrex::AllPrint() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << " it-" << it << ": NC_epochTime\t" << NC_epochTime[it] << std::endl; + // NC_dateTime is only on the IOProcessor. + amrex::Print() << " DJW init_from_metgrid proc-" << ParallelDescriptor::MyProc() << ": NC_dateTime \t" << NC_dateTime[it] << std::endl; +#endif + } // it + + // Verify that files in nc_init_file[lev] are ordered from earliest to latest. + for (int it = 1; it < ntimes; it++) AMREX_ALWAYS_ASSERT(NC_epochTime[it] > NC_epochTime[it-1]); + + // Start at the earliest time in nc_init_file[lev]. + start_bdy_time = NC_epochTime[0]; + t_new[lev] = start_bdy_time; + t_old[lev] = start_bdy_time - 1.e200; + + // Determine the spacing between met_em files. + bdy_time_interval = NC_epochTime[1]-NC_epochTime[0]; + + // Verify that met_em files have even spacing in time. + for (int it = 1; it < ntimes; it++) { + Real NC_dt = NC_epochTime[it]-NC_epochTime[it-1]; +#ifndef AMREX_USE_GPU + amrex::Print() << " " << nc_init_file[lev][it-1] << " / " << nc_init_file[lev][it] << " are " << NC_dt << " seconds apart" << std::endl; +#endif + if (NC_dt != bdy_time_interval) amrex::Error("Time interval between consecutive met_em files must be consistent."); + } + + // Set up a FAB for mixing ratio and another for potential temperature. + // Necessary because the input data has relative humidity and temperature, not mixing ratio and potential temperature. + // TODO: add alternate pathways for other origin models where different combinations of variables may be present. + Vector mxrat_fab; mxrat_fab.resize(ntimes); + Vector theta_fab; theta_fab.resize(ntimes); + for (int it = 0; it < ntimes; it++) { + Box NC_box_unstag = NC_rhum_fab[it].box(); +#ifdef AMREX_USE_GPU + mxrat_fab[it].resize(NC_box_unstag, 1, The_Pinned_Arena()); + theta_fab[it].resize(NC_box_unstag, 1, The_Pinned_Arena()); +#else + mxrat_fab[it].resize(NC_box_unstag, 1); + theta_fab[it].resize(NC_box_unstag, 1); +#endif } auto& lev_new = vars_new[lev]; @@ -93,108 +145,418 @@ ERF::init_from_metgrid (int lev) AMREX_ALWAYS_ASSERT(solverChoice.use_terrain); - z_phys->setVal(0.); + // Verify that the terrain height (HGT_M) was in each met_em file. + for (int it = 0; it < ntimes; it++) AMREX_ALWAYS_ASSERT(flag_hgt[it] == 1); - for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { + z_phys->setVal(0.0); + + for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { // This defines only the z(i,j,0) values given the FAB filled from the NetCDF input FArrayBox& z_phys_nd_fab = (*z_phys)[mfi]; - init_terrain_from_metgrid(lev, z_phys_nd_fab, NC_hgt_fab); + init_terrain_from_metgrid(z_phys_nd_fab, NC_hgt_fab); } // mf // This defines all the z(i,j,k) values given z(i,j,0) from above. init_terrain_grid(geom[lev], *z_phys, zlevels_stag); + // Copy SST and LANDMASK data into MF and iMF data structures + auto& ba = lev_new[Vars::cons].boxArray(); + auto& dm = lev_new[Vars::cons].DistributionMap(); + auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0; + BoxList bl2d = ba.boxList(); + for (auto& b : bl2d) { + b.setRange(2,0); + } + BoxArray ba2d(std::move(bl2d)); + int i_lo = geom[lev].Domain().smallEnd(0); int i_hi = geom[lev].Domain().bigEnd(0); + int j_lo = geom[lev].Domain().smallEnd(1); int j_hi = geom[lev].Domain().bigEnd(1); + if (flag_sst[0]) { + for (int it = 0; it < ntimes; ++it) { + sst_lev[lev][it] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(sst_lev[lev][it]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Box gtbx = mfi.growntilebox(); + FArrayBox& dst = (*(sst_lev[lev][it]))[mfi]; + FArrayBox& src = NC_sst_fab[it]; + const Array4< Real>& dst_arr = dst.array(); + const Array4& src_arr = src.const_array(); + amrex::ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + dst_arr(i,j,0) = src_arr(li,lj,0); + }); + } + sst_lev[lev][it]->FillBoundary(geom[lev].periodicity()); + } + } else { + for (int it = 0; it < ntimes; ++it) sst_lev[lev][it] = nullptr; + } + if (flag_lmask[0]) { + for (int it = 0; it < ntimes; ++it) { + lmask_lev[lev][it] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(lmask_lev[lev][it]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Box gtbx = mfi.growntilebox(); + IArrayBox& dst = (*(lmask_lev[lev][it]))[mfi]; + IArrayBox& src = NC_lmask_iab[it]; + const Array4< int>& dst_arr = dst.array(); + const Array4& src_arr = src.const_array(); + amrex::ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + dst_arr(i,j,0) = src_arr(li,lj,0); + }); + } + lmask_lev[lev][it]->FillBoundary(geom[lev].periodicity()); + } + } else { + for (int it = 0; it < ntimes; ++it) lmask_lev[lev][it] = nullptr; + } + + for (int it = 0; it < ntimes; it++) { + // Verify that the grid size and resolution from met_em file matches that in geom (from ERF inputs file). + AMREX_ALWAYS_ASSERT(geom[lev].CellSizeArray()[0] == NC_dx[it]); + AMREX_ALWAYS_ASSERT(geom[lev].CellSizeArray()[1] == NC_dy[it]); + // NC_nx-2 because NC_nx is the number of staggered grid points indexed from 1. + AMREX_ALWAYS_ASSERT(geom[lev].Domain().bigEnd(0) == NC_nx[it]-2); + // NC_ny-2 because NC_ny is the number of staggered grid points indexed from 1. + AMREX_ALWAYS_ASSERT(geom[lev].Domain().bigEnd(1) == NC_ny[it]-2); + } // it + // This makes the Jacobian. - make_J (geom[lev],*z_phys, *detJ_cc[lev]); + make_J(geom[lev],*z_phys, *detJ_cc[lev]); // This defines z at w-cell faces. make_zcc(geom[lev],*z_phys,*z_phys_cc[lev]); + // Set up FABs to hold data that will be used to set lateral boundary conditions. +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + int MetGridBdyEnd = MetGridBdyVars::NumTypes; +#else + int MetGridBdyEnd = MetGridBdyVars::NumTypes-1; +#endif + //amrex::Vector > fabs_for_bcs; + amrex::Vector> fabs_for_bcs; + fabs_for_bcs.resize(ntimes); + for (int it(0); it < ntimes; it++) { + fabs_for_bcs[it].resize(MetGridBdyEnd); + + Box gdomain; + Box ldomain; + for (int nvar(0); nvar(0.0); + } + } + + + const Real l_rdOcp = solverChoice.rdOcp; + std::unique_ptr mask_c = OwnerMask(lev_new[Vars::cons], geom[lev].periodicity());//, lev_new[Vars::cons].nGrowVect()); + std::unique_ptr mask_u = OwnerMask(lev_new[Vars::xvel], geom[lev].periodicity());//, lev_new[Vars::xvel].nGrowVect()); + std::unique_ptr mask_v = OwnerMask(lev_new[Vars::yvel], geom[lev].periodicity());//, lev_new[Vars::yvel].nGrowVect()); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - // Define fabs for holding the initial data + for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + Box tbxc = mfi.tilebox(); + Box tbxu = mfi.tilebox(IntVect(1,0,0)); + Box tbxv = mfi.tilebox(IntVect(0,1,0)); + Box tbxw = mfi.tilebox(IntVect(0,0,1)); + + // Define FABs for hlding some of the initial data FArrayBox &cons_fab = lev_new[Vars::cons][mfi]; FArrayBox &xvel_fab = lev_new[Vars::xvel][mfi]; FArrayBox &yvel_fab = lev_new[Vars::yvel][mfi]; FArrayBox &zvel_fab = lev_new[Vars::zvel][mfi]; - - FArrayBox& z_phys_nd_fab = (*z_phys)[mfi]; - init_state_from_metgrid(lev, cons_fab, xvel_fab, yvel_fab, zvel_fab, + FArrayBox &z_phys_nd_fab = (*z_phys)[mfi]; + + const Array4& mask_c_arr = mask_c->const_array(mfi); + const Array4& mask_u_arr = mask_u->const_array(mfi); + const Array4& mask_v_arr = mask_v->const_array(mfi); + + // Fill state data using origin data (initialization and BC arrays) + // x_vel interpolated from origin levels + // y_vel interpolated from origin levels + // z_vel set to 0.0 + // theta calculate on origin levels then interpolate + // mxrat convert RH -> Q on origin levels then interpolate + init_state_from_metgrid(l_rdOcp, + tbxc, tbxu, tbxv, + cons_fab, xvel_fab, yvel_fab, zvel_fab, z_phys_nd_fab, - NC_hgt_fab, NC_xvel_fab, NC_yvel_fab, NC_temp_fab, - NC_rhum_fab, NC_pres_fab); + NC_hgt_fab, NC_ght_fab, NC_xvel_fab, + NC_yvel_fab, NC_temp_fab, NC_rhum_fab, + NC_pres_fab, theta_fab, mxrat_fab, + fabs_for_bcs, mask_c_arr, mask_u_arr, mask_v_arr); } // mf + #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - // Map scale factors common for "ideal" as well as "real" simulation - for ( MFIter mfi(*mapfac_u[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { + // Use map scale factors directly from the met_em files + for ( MFIter mfi(*mapfac_u[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { // Define fabs for holding the initial data FArrayBox &msfu_fab = (*mapfac_u[lev])[mfi]; FArrayBox &msfv_fab = (*mapfac_v[lev])[mfi]; FArrayBox &msfm_fab = (*mapfac_m[lev])[mfi]; - init_msfs_from_metgrid(lev, msfu_fab, msfv_fab, msfm_fab, + init_msfs_from_metgrid(msfu_fab, msfv_fab, msfm_fab, + flag_msfu[0], flag_msfv[0], flag_msfm[0], NC_MSFU_fab, NC_MSFV_fab, NC_MSFM_fab); } // mf + 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 MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component + for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + FArrayBox& p_hse_fab = p_hse[mfi]; + FArrayBox& pi_hse_fab = pi_hse[mfi]; + FArrayBox& r_hse_fab = r_hse[mfi]; + FArrayBox& cons_fab = lev_new[Vars::cons][mfi]; + FArrayBox& z_phys_nd_fab = (*z_phys)[mfi]; - const Real l_rdOcp = solverChoice.rdOcp; + const Array4& mask_c_arr = mask_c->const_array(mfi); + + // Fill base state data using origin data (initialization and BC arrays) + // p_hse calculate dry pressure + // r_hse calculate dry density + // pi_hse calculate Exner term given pressure + const Box valid_bx = mfi.validbox(); + init_base_state_from_metgrid(l_rdOcp, + valid_bx, + flag_psfc, + cons_fab, r_hse_fab, p_hse_fab, pi_hse_fab, + z_phys_nd_fab, NC_ght_fab, NC_psfc_fab, + fabs_for_bcs, mask_c_arr); + } // mf - if (init_type == "real") { - for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - FArrayBox& p_hse_fab = p_hse[mfi]; - FArrayBox& pi_hse_fab = pi_hse[mfi]; - FArrayBox& r_hse_fab = r_hse[mfi]; - const Box& bx = mfi.validbox(); - //init_base_state_from_metgrid(lev, bx, l_rdOcp, p_hse_fab, pi_hse_fab, r_hse_fab, - // NC_ALB_fab, NC_PB_fab); + // NOTE: fabs_for_bcs is defined over the whole domain on each rank. + // However, the operations needed to define the data on the ERF + // grid are done over MultiFab boxes that are local to the rank. + // So when we save the data in fabs_for_bc, only regions owned + // by the rank are populated. Use an allreduce sum to make the + // complete data set; initialized to 0 above. + for (int it(0); it < ntimes; it++) { + for (int nvar(0); nvar& R_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::R].const_array(); + + for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { + + auto xlo_arr = bdy_data_xlo[it][ivar].array(); + auto xhi_arr = bdy_data_xhi[it][ivar].array(); + auto ylo_arr = bdy_data_ylo[it][ivar].array(); + auto yhi_arr = bdy_data_yhi[it][ivar].array(); + const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); + + if (ivar == MetGridBdyVars::U) { + multiply_rho = false; + xlo_plane = xlo_plane_x_stag; xhi_plane = xhi_plane_x_stag; + ylo_plane = ylo_plane_x_stag; yhi_plane = yhi_plane_x_stag; + } else if (ivar == MetGridBdyVars::V) { + multiply_rho = false; + xlo_plane = xlo_plane_y_stag; xhi_plane = xhi_plane_y_stag; + ylo_plane = ylo_plane_y_stag; yhi_plane = yhi_plane_y_stag; + } else if (ivar == MetGridBdyVars::R) { + multiply_rho = false; + xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; + ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; + } else if (ivar == MetGridBdyVars::T) { + multiply_rho = false; + xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; + ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; + } else if (ivar == MetGridBdyVars::QV) { + multiply_rho = true; + xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; + ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; + } // MetGridBdyVars::QV + + // west boundary + amrex::ParallelFor(xlo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + }); + // xvel at east boundary + amrex::ParallelFor(xhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + }); + // xvel at south boundary + amrex::ParallelFor(ylo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + }); + // xvel at north boundary + amrex::ParallelFor(yhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + amrex::Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + }); + + } // ivar + } // it +} + /** - * Helper function to initialize terrain nodal z coordinates for a Fab - * given metgrid data. + * Helper function to initialize terrain nodal z coordinates given metgrid data. * - * @param lev Integer specifying the current level * @param z_phys_nd_fab FArrayBox (Fab) holding the nodal z coordinates for terrain data we want to fill * @param NC_hgt_fab Vector of FArrayBox objects holding height data read from NetCDF files for metgrid data */ void -init_terrain_from_metgrid (int lev, FArrayBox& z_phys_nd_fab, +init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, const Vector& NC_hgt_fab) { - int nboxes = NC_hgt_fab.size(); - - // NOTE NOTE NOTE -- this routine currently only fills the k=0 value - // TODO: we need to fill the rest of the values from the pressure (?) variable... + int ntimes = 1; // Use terrain from the first met_em file. - for (int idx = 0; idx < nboxes; idx++) - { + for (int it = 0; it < ntimes; it++) { #ifndef AMREX_USE_GPU - amrex::Print() << " SIZE OF HGT FAB " << NC_hgt_fab[idx].box() << std::endl; + amrex::Print() << " SIZE OF HGT FAB " << NC_hgt_fab[it].box() << std::endl; amrex::Print() << " SIZE OF ZP FAB " << z_phys_nd_fab.box() << std::endl; #endif // This copies from NC_zphys on z-faces to z_phys_nd on nodes const Array4& z_arr = z_phys_nd_fab.array(); - const Array4& nc_hgt_arr = NC_hgt_fab[idx].const_array(); + const Array4& nc_hgt_arr = NC_hgt_fab[it].const_array(); - const Box z_hgt_box = NC_hgt_fab[idx].box(); + const Box z_hgt_box = NC_hgt_fab[it].box(); int ilo = z_hgt_box.smallEnd()[0]; int ihi = z_hgt_box.bigEnd()[0]; @@ -202,272 +564,467 @@ init_terrain_from_metgrid (int lev, FArrayBox& z_phys_nd_fab, int jhi = z_hgt_box.bigEnd()[1]; Box z_phys_box = z_phys_nd_fab.box(); - Box from_box = surroundingNodes(NC_hgt_fab[idx].box()); from_box.growHi(2,-1); + Box from_box = surroundingNodes(NC_hgt_fab[it].box()); + from_box.growHi(2,-1); + Box bx = z_phys_box & from_box; + Box bxu = bx; bxu.growHi(0,1); + Box bxv = bx; bxv.growHi(1,1); #ifndef AMREX_USE_GPU amrex::Print() << "FROM BOX " << from_box << std::endl; amrex::Print() << "BX " << bx << std::endl; #endif - // - // We must be careful not to read out of bounds of the WPS data - // - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { int ii = std::max(std::min(i,ihi-1),ilo+1); int jj = std::max(std::min(j,jhi-1),jlo+1); z_arr(i,j,k) = 0.25 * ( nc_hgt_arr (ii,jj ,k) + nc_hgt_arr(ii-1,jj ,k) + nc_hgt_arr (ii,jj-1,k) + nc_hgt_arr(ii-1,jj-1,k) ); }); - } // idx + } // it } /** - * Helper function to initialize state and velocity data - * read from metgrid data. + * Helper function to initialize state and velocity data read from metgrid data. * - * @param lev Integer specifying the current level + * @param l_rdOcp Real constant specifying Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$) * @param state_fab FArrayBox holding the state data to initialize * @param x_vel_fab FArrayBox holding the x-velocity data to initialize * @param y_vel_fab FArrayBox holding the y-velocity data to initialize * @param z_vel_fab FArrayBox holding the z-velocity data to initialize * @param z_phys_nd_fab FArrayBox holding nodal z coordinate data for terrain - * @param NC_hgt_fab Vector of FArrayBox obects holding metgrid data for height + * @param NC_hgt_fab Vector of FArrayBox obects holding metgrid data for terrain height + * @param NC_ght_fab Vector of FArrayBox objects holding metgrid data for height of cell centers * @param NC_xvel_fab Vector of FArrayBox obects holding metgrid data for x-velocity * @param NC_yvel_fab Vector of FArrayBox obects holding metgrid data for y-velocity * @param NC_zvel_fab Vector of FArrayBox obects holding metgrid data for z-velocity - * @param NC_rho_fab Vector of FArrayBox obects holding metgrid data for density - * @param NC_rhotheta_fab Vector of FArrayBox obects holding metgrid data for (density * potential temperature) + * @param NC_temp_fab Vector of FArrayBox obects holding metgrid data for temperature + * @param NC_rhum_fab Vector of FArrayBox obects holding metgrid data for relative humidity + * @param NC_pres_fab Vector of FArrayBox obects holding metgrid data for pressure + * @param theta_fab Vector of FArrayBox obects holding potential temperature calculated from temperature and pressure + * @param mxrat_fab Vector of FArrayBox obects holding vapor mixing ratio calculated from relative humidity + * @param fabs_for_bcs Vector of Vector of FArrayBox objects holding MetGridBdyVars at each met_em time. */ void -init_state_from_metgrid (int lev, FArrayBox& state_fab, - FArrayBox& x_vel_fab, FArrayBox& y_vel_fab, - FArrayBox& z_vel_fab, FArrayBox& z_phys_nd_fab, +init_state_from_metgrid (const Real l_rdOcp, + Box& tbxc, + Box& tbxu, + Box& tbxv, + FArrayBox& state_fab, + FArrayBox& x_vel_fab, + FArrayBox& y_vel_fab, + FArrayBox& z_vel_fab, + FArrayBox& z_phys_nd_fab, const Vector& NC_hgt_fab, + const Vector& NC_ght_fab, const Vector& NC_xvel_fab, const Vector& NC_yvel_fab, - const Vector& NC_zvel_fab, - const Vector& NC_rho_fab, - const Vector& NC_rhotheta_fab) + const Vector& NC_temp_fab, + const Vector& NC_rhum_fab, + const Vector& NC_pres_fab, + Vector& theta_fab, + Vector& mxrat_fab, + amrex::Vector>& fabs_for_bcs, + const amrex::Array4& mask_c_arr, + const amrex::Array4& mask_u_arr, + const amrex::Array4& mask_v_arr) { - int nboxes = NC_hgt_fab.size(); - -#ifndef AMREX_USE_GPU - amrex::Print() << " U FROM NC " << NC_xvel_fab[0].box() << std::endl; - amrex::Print() << " U INTO FAB " << x_vel_fab.box() << std::endl; - exit(0); -#endif - for (int idx = 0; idx < nboxes; idx++) + int ntimes = NC_hgt_fab.size(); + for (int it = 0; it < ntimes; it++) { // ******************************************************** // U // ******************************************************** { - Box bx2d = NC_xvel_fab[idx].box() & x_vel_fab.box(); + Box bx2d = NC_xvel_fab[it].box() & tbxu; bx2d.setRange(2,0); - auto const orig_data = NC_xvel_fab[idx].const_array(); - auto const orig_z = NC_hgt_fab[idx].const_array(); - + auto const orig_data = NC_xvel_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); auto new_data = x_vel_fab.array(); + auto bc_data = fabs_for_bcs[it][MetGridBdyVars::U].array(); auto const new_z = z_phys_nd_fab.const_array(); - ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) + int kmax = amrex::ubound(tbxu).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - interpolate_column(i,j,0,0,orig_z,orig_data,new_z,new_data); + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid(i,j,k,'X',0,orig_z,orig_data,new_z); + if (mask_u_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,0) = Interp_Val; + } }); } + // ******************************************************** // V // ******************************************************** { - Box bx2d = NC_yvel_fab[idx].box() & y_vel_fab.box(); + Box bx2d = NC_yvel_fab[it].box() & tbxv; bx2d.setRange(2,0); - auto const orig_data = NC_yvel_fab[idx].const_array(); - auto const orig_z = NC_hgt_fab[idx].const_array(); - + auto const orig_data = NC_yvel_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); auto new_data = y_vel_fab.array(); + auto bc_data = fabs_for_bcs[it][MetGridBdyVars::V].array(); auto const new_z = z_phys_nd_fab.const_array(); - ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) + int kmax = amrex::ubound(tbxv).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - interpolate_column(i,j,0,0,orig_z,orig_data,new_z,new_data); + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid(i,j,k,'Y',0,orig_z,orig_data,new_z); + if (mask_v_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,0) = Interp_Val; + } }); } + // ******************************************************** // W // ******************************************************** - z_vel_fab.template setVal(0.); + if (it == 0) { // update at initialization + z_vel_fab.template setVal(0.0); + } + + // ******************************************************** + // Initialize all state_fab variables to zero + // ******************************************************** + if (it == 0) { // update at initialization + state_fab.template setVal(0.0); + } + // ******************************************************** - // rho + // theta // ******************************************************** + { // calculate potential temperature. + Box bx = NC_rhum_fab[it].box() & tbxc; + auto const temp = NC_temp_fab[it].const_array(); + auto const pres = NC_pres_fab[it].const_array(); + auto theta = theta_fab[it].array(); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + theta(i,j,k) = getThgivenPandT(temp(i,j,k),pres(i,j,k),l_rdOcp); + //theta(i,j,k) = 300.0; // TODO: Remove when not needed. Force an isothermal atmosphere for debugging. + }); + } + + // vertical interpolation of potential temperature. { - Box bx2d = NC_rho_fab[idx].box() & state_fab.box(); + Box bx2d = NC_temp_fab[it].box() & tbxc; bx2d.setRange(2,0); - - auto const orig_data = NC_rho_fab[idx].const_array(); - auto const orig_z = NC_hgt_fab[idx].const_array(); - auto new_data = state_fab.array(); + auto const orig_data = theta_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); + auto new_data = state_fab.array(); + auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); auto const new_z = z_phys_nd_fab.const_array(); - ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) + int kmax = amrex::ubound(tbxc).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - interpolate_column(i,j,0,Rho_comp,orig_z,orig_data,new_z,new_data); + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid(i,j,k,'M',0,orig_z,orig_data,new_z); + if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,RhoTheta_comp) = Interp_Val; + } }); } +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) // ******************************************************** - // rho_theta + // specific humidity / relative humidity / mixing ratio // ******************************************************** + // TODO: we will need to check what input data we have for moisture + // and then, if necessary, compute mixing ratio. For now, we will + // focus on the case where we have relative humidity. Alternate cases + // could be specific humidity or a mixing ratio. + // + { // calculate vapor mixing ratio from relative humidity. + Box bx = NC_temp_fab[it].box() & tbxc; + auto const rhum = NC_rhum_fab[it].const_array(); + auto const temp = NC_temp_fab[it].const_array(); + auto const pres = NC_pres_fab[it].const_array(); + auto mxrat = mxrat_fab[it].array(); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + rh_to_mxrat(i,j,k,rhum,temp,pres,mxrat); + }); + } + + // vertical interpolation of vapor mixing ratio. { - Box bx2d = NC_rhotheta_fab[idx].box() & state_fab.box(); + Box bx2d = NC_temp_fab[it].box() & tbxc; bx2d.setRange(2,0); - - auto const orig_data = NC_rhotheta_fab[idx].const_array(); - auto const orig_z = NC_hgt_fab[idx].const_array(); + auto const orig_data = mxrat_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); auto new_data = state_fab.array(); + auto bc_data = fabs_for_bcs[it][MetGridBdyVars::QV].array(); auto const new_z = z_phys_nd_fab.const_array(); - ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) + int kmax = amrex::ubound(tbxc).z; + +#if defined(ERF_USE_MOISTURE) + int state_indx = RhoQt_comp; +#elif defined(ERF_USE_WARM_NO_PRECIP) + int state_indx = RhoQv_comp; +#endif + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - interpolate_column(i,j,0,RhoTheta_comp,orig_z,orig_data,new_z,new_data); + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid(i,j,k,'M',0,orig_z,orig_data,new_z); + if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,state_indx) = Interp_Val; + } }); } - } // idx -} - -/** - * Helper function to initialize map factors from metgrid data - * - * @param lev Integer specifying current level - * @param msfu_fab FArrayBox specifying x-velocity map factors - * @param msfv_fab FArrayBox specifying y-velocity map factors - * @param msfm_fab FArrayBox specifying z-velocity map factors - * @param NC_MSFU_fab Vector of FArrayBox objects holding metgrid data for x-velocity map factors - * @param NC_MSFV_fab Vector of FArrayBox objects holding metgrid data for y-velocity map factors - * @param NC_MSFM_fab Vector of FArrayBox objects holding metgrid data for z-velocity map factors - */ -void -init_msfs_from_metgrid (int lev, FArrayBox& msfu_fab, - FArrayBox& msfv_fab, FArrayBox& msfm_fab, - const Vector& NC_MSFU_fab, - const Vector& NC_MSFV_fab, - const Vector& NC_MSFM_fab) -{ - int nboxes = NC_MSFU_fab.size(); - - for (int idx = 0; idx < nboxes; idx++) - { - // - // FArrayBox to FArrayBox copy does "copy on intersection" - // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks - // - // This copies mapfac_u - msfu_fab.template copy(NC_MSFU_fab[idx]); +#endif - // This copies mapfac_v - msfv_fab.template copy(NC_MSFV_fab[idx]); + // TODO: TEMPORARY CODE TO RUN QUIESCENT, REMOVE WHEN NOT NEEDED. +// if (it == 0) { +// x_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. +// y_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. +// } +// fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. +// fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. - // This copies mapfac_m - msfm_fab.template copy(NC_MSFM_fab[idx]); - } // idx + } // it } + /** * Helper function for initializing hydrostatic base state data from metgrid data * - * @param lev Integer specifying the current level - * @param valid_bx Box specifying the index space we are initializing * @param l_rdOcp Real constant specifying Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$) - * @param p_hse FArrayBox holding the hydrostatic base state pressure we are initializing - * @param pi_hse FArrayBox holding the hydrostatic base Exner pressure we are initializing - * @param r_hse FArrayBox holding the hydrostatic base state density we are initializing - * @param NC_ALB_fab Vector of FArrayBox objects holding metgrid data specifying 1/density - * @param NC_PB_fab Vector of FArrayBox objects holding metgrid data specifying pressure + * @param valid_bx Box specifying the index space we are to initialize + * @param flag_psfc Vector of Integer 1 if surface pressure is in metgrid data, 0 otherwise + * @param state_fab FArrayBox holding the state data to initialize + * @param r_hse_fab FArrayBox holding the hydrostatic base state density we are initializing + * @param p_hse_fab FArrayBox holding the hydrostatic base state pressure we are initializing + * @param pi_hse_fab FArrayBox holding the hydrostatic base Exner pressure we are initializing + * @param z_phys_nd_fab FArrayBox holding nodal z coordinate data for terrain + * @param NC_ght_fab Vector of FArrayBox objects holding metgrid data for height of cell centers + * @param NC_psfc_fab Vector of FArrayBox objects holding metgrid data for surface pressure + * @param fabs_for_bcs Vector of Vector of FArrayBox objects holding MetGridBdyVars at each met_em time. */ void -init_base_state_from_metgrid (int lev, const Box& valid_bx, const Real l_rdOcp, - FArrayBox& p_hse, FArrayBox& pi_hse, FArrayBox& r_hse, - const Vector& NC_ALB_fab, - const Vector& NC_PB_fab) +init_base_state_from_metgrid (const Real l_rdOcp, + const Box& valid_bx, + const Vector& flag_psfc, + FArrayBox& state_fab, + FArrayBox& r_hse_fab, + FArrayBox& p_hse_fab, + FArrayBox& pi_hse_fab, + FArrayBox& z_phys_nd_fab, + const Vector& NC_ght_fab, + const Vector& NC_psfc_fab, + Vector>& fabs_for_bcs, + const amrex::Array4& mask_c_arr) { - int nboxes = NC_ALB_fab.size(); +#if defined(ERF_USE_MOISTURE) + int RhoQ_comp = RhoQt_comp; +#elif defined(ERF_USE_WARM_NO_PRECIP) + int RhoQ_comp = RhoQv_comp; +#endif + int kmax = amrex::ubound(valid_bx).z; + + // Device vectors for columnwise operations + Gpu::DeviceVector z_vec_d(kmax+2,0); Real* z_vec = z_vec_d.data(); + Gpu::DeviceVector Thetad_vec_d(kmax+1,0); Real* Thetad_vec = Thetad_vec_d.data(); + Gpu::DeviceVector Thetam_vec_d(kmax+1,0); Real* Thetam_vec = Thetam_vec_d.data(); + Gpu::DeviceVector Rhod_vec_d(kmax+1,0); Real* Rhod_vec = Rhod_vec_d.data(); + Gpu::DeviceVector Rhom_vec_d(kmax+1,0); Real* Rhom_vec = Rhom_vec_d.data(); + Gpu::DeviceVector Pd_vec_d(kmax+1,0); Real* Pd_vec = Pd_vec_d.data(); + Gpu::DeviceVector Pm_vec_d(kmax+1,0); Real* Pm_vec = Pm_vec_d.data(); +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + Gpu::DeviceVector Q_vec_d(kmax+1,0); Real* Q_vec = Q_vec_d.data(); +#endif - for (int idx = 0; idx < nboxes; idx++) - { - // - // FArrayBox to FArrayBox copy does "copy on intersection" - // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks - // - const Array4& p_hse_arr = p_hse.array(); - const Array4& pi_hse_arr = pi_hse.array(); - const Array4& r_hse_arr = r_hse.array(); - const Array4& alpha_arr = NC_ALB_fab[idx].const_array(); - const Array4& nc_pb_arr = NC_PB_fab[idx].const_array(); - - amrex::ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - p_hse_arr(i,j,k) = nc_pb_arr(i,j,k); + // Device vectors for psfc flags + Gpu::DeviceVectorflag_psfc_d(flag_psfc.size()); + Gpu::copy(Gpu::hostToDevice, flag_psfc.begin(), flag_psfc.end(), flag_psfc_d.begin()); + int* flag_psfc_vec = flag_psfc_d.data(); + + { // set pressure and density at initialization. + const Array4& r_hse_arr = r_hse_fab.array(); + const Array4& p_hse_arr = p_hse_fab.array(); + const Array4& pi_hse_arr = pi_hse_fab.array(); + + // ******************************************************** + // calculate dry density and dry pressure + // ******************************************************** + // calculate density and dry pressure on the new grid. + Box valid_bx2d = valid_bx; + valid_bx2d.setRange(2,0); + auto const orig_psfc = NC_psfc_fab[0].const_array(); + auto new_data = state_fab.array(); + auto const new_z = z_phys_nd_fab.const_array(); + + amrex::ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + for (int k=0; k<=kmax; k++) { + z_vec[k] = new_z(i,j,k); + Thetad_vec[k] = new_data(i,j,k,RhoTheta_comp); +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + Q_vec[k] = new_data(i,j,k,RhoQ_comp); +#endif + } + z_vec[kmax+1] = new_z(i,j,kmax+1); + + calc_rho_p(kmax,flag_psfc_vec[0],orig_psfc(i,j,0),Thetad_vec,Thetam_vec, +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + Q_vec, +#endif + z_vec,Rhod_vec,Rhom_vec,Pd_vec,Pm_vec); + + for (int k=0; k<=kmax; k++) { + p_hse_arr(i,j,k) = Pd_vec[k]; + r_hse_arr(i,j,k) = Rhod_vec[k]; + } + }); + + amrex::ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k); + new_data(i,j,k,RhoScalar_comp) = 0.0; + // RhoTheta and RhoQt or RhoQv currently hold Theta and Qt or Qv. Multiply by Rho. + Real RhoTheta = r_hse_arr(i,j,k)*new_data(i,j,k,RhoTheta_comp); + new_data(i,j,k,RhoTheta_comp) = RhoTheta; +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + Real RhoQ = r_hse_arr(i,j,k)*new_data(i,j,k,RhoQ_comp); + new_data(i,j,k,RhoQ_comp) = RhoQ; +#endif pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp); - r_hse_arr(i,j,k) = 1.0 / alpha_arr(i,j,k); + //pi_hse_arr(i,j,k) = getExnergivenRTh(RhoTheta, l_rdOcp); + }); + } + int ntimes = NC_psfc_fab.size(); + for (int it=0; it& orig_z, const Array4& orig_data, - const Array4& new_z, const Array4& new_data) +init_msfs_from_metgrid (FArrayBox& msfu_fab, + FArrayBox& msfv_fab, + FArrayBox& msfm_fab, + const int& flag_msfu, + const int& flag_msfv, + const int& flag_msfm, + const Vector& NC_MSFU_fab, + const Vector& NC_MSFV_fab, + const Vector& NC_MSFM_fab) { - // CAVEAT: we only consider interpolation here - if we go past end of array this won't work for now +// int ntimes = NC_MSFU_fab.size(); + int ntimes = 1; + for (int it = 0; it < ntimes; it++) { + // + // FArrayBox to FArrayBox copy does "copy on intersection" + // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks + // - int kmax = amrex::ubound(Box(new_data)).z; - int kmax_old = amrex::ubound(Box(orig_data)).z; + // This copies or sets mapfac_m + if (flag_msfm == 1) { + msfm_fab.template copy(NC_MSFM_fab[it]); + } else { #ifndef AMREX_USE_GPU - amrex::Print() << "KMAX IN INTERP " << kmax << std::endl; - amrex::Print() << "KMAX_OLD IN INTERP " << kmax_old << std::endl; + amrex::Print() << " MAPFAC_M not present in met_em files. Setting to 1.0" << std::endl; #endif - int klast = 0; - for (int k = 0; k < kmax; k++) { - - Real z = new_z(i,j,k); + msfm_fab.template setVal(1.0); + } - bool kfound = false; + // This copies or sets mapfac_u + if (flag_msfu == 1) { + msfu_fab.template copy(NC_MSFU_fab[it]); + } else { +#ifndef AMREX_USE_GPU + amrex::Print() << " MAPFAC_U not present in met_em files. Setting to 1.0" << std::endl; +#endif + msfu_fab.template setVal(1.0); + } + // This copies or sets mapfac_v + if (flag_msfv == 1) { + msfv_fab.template copy(NC_MSFV_fab[it]); + } else { #ifndef AMREX_USE_GPU - amrex::Print() << "AT THIS K WE ARE USING KLAST " << k << " " << klast << std::endl; + amrex::Print() << " MAPFAC_V not present in met_em files. Setting to 1.0" << std::endl; #endif - for (int kk = klast; kk < kmax_old; kk++) { - if (z >= orig_z(i,j,k)) { - Real y0 = orig_data(i,j,kk ,src_comp); - Real y1 = orig_data(i,j,kk+1,src_comp); - Real z0 = orig_z(i,j,kk ); - Real z1 = orig_z(i,j,kk+1); - new_data(i,j,k,dest_comp) = y0 + (y1 - y0)*(z - z0) / (z1 - z0); - klast = kk; - return; - } + msfv_fab.template setVal(1.0); } - if (!kfound) amrex::Error("Wasnt able to interpolate"); - } -} + } // it +} #endif // ERF_USE_NETCDF diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index 0bfec76ea..d29fe922d 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -254,6 +254,10 @@ ERF::init_from_wrfinput (int lev) NC_C1H_fab[0], NC_C2H_fab[0], NC_RDNW_fab[0], NC_xvel_fab[0],NC_yvel_fab[0],NC_rho_fab[0],NC_rhoth_fab[0]); } + + // Start at the earliest time (read_from_wrfbdy) + t_new[lev] = start_bdy_time; + t_old[lev] = start_bdy_time - 1.e200; } /** diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 9115a0dac..d3bb99d4b 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -8,6 +8,7 @@ CEXE_sources += ERF_init_bcs.cpp CEXE_sources += ERF_init1d.cpp ifeq ($(USE_NETCDF),TRUE) +CEXE_headers += Metgrid_utils.H CEXE_sources += ERF_init_from_wrfinput.cpp CEXE_sources += ERF_init_from_metgrid.cpp endif diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H new file mode 100644 index 000000000..af06ecd5d --- /dev/null +++ b/Source/Initialization/Metgrid_utils.H @@ -0,0 +1,347 @@ +#ifndef _METGRIDUTIL_H_ +#define _METGRIDUTIL_H_ + +#include +#include +#include +#include +#include + +void +read_from_metgrid (int lev, + const amrex::Box& domain, + const std::string& fname, + std::string& NC_dateTime, + amrex::Real& NC_epochTime, + int& flag_psfc, + int& flag_msfu, + int& flag_msfv, + int& flag_msfm, + int& flag_hgt, + int& flag_sst, + int& flag_lmask, + int& NC_nx, + int& NC_ny, + amrex::Real& NC_dx, + amrex::Real& NC_dy, + amrex::FArrayBox& NC_xvel_fab, + amrex::FArrayBox& NC_yvel_fab, + amrex::FArrayBox& NC_temp_fab, + amrex::FArrayBox& NC_rhum_fab, + amrex::FArrayBox& NC_pres_fab, + amrex::FArrayBox& NC_ght_fab, + amrex::FArrayBox& NC_hgt_fab, + amrex::FArrayBox& NC_psfc_fab, + amrex::FArrayBox& NC_msfu_fab, + amrex::FArrayBox& NC_msfv_fab, + amrex::FArrayBox& NC_msfm_fab, + amrex::FArrayBox& NC_sst_fab, + amrex::IArrayBox& NC_lmask_iab); + + + +void +init_terrain_from_metgrid (amrex::FArrayBox& z_phys_nd_fab, + const amrex::Vector& NC_hgt_fab); + +void +init_state_from_metgrid (const amrex::Real l_rdOcp, + amrex::Box& tbxc, + amrex::Box& tbxu, + amrex::Box& tbxv, + amrex::FArrayBox& state_fab, + amrex::FArrayBox& x_vel_fab, + amrex::FArrayBox& y_vel_fab, + amrex::FArrayBox& z_vel_fab, + amrex::FArrayBox& z_phys_nd_fab, + const amrex::Vector& NC_hgt_fab, + const amrex::Vector& NC_ght_fab, + const amrex::Vector& NC_xvel_fab, + const amrex::Vector& NC_yvel_fab, + const amrex::Vector& NC_zvel_fab, + const amrex::Vector& NC_temp_fab, + const amrex::Vector& NC_rhum_fab, + amrex::Vector& theta_fab, + amrex::Vector& mxrat_fab, + amrex::Vector>& fabs_for_bcs, + const amrex::Array4& mask_c_arr, + const amrex::Array4& mask_u_arr, + const amrex::Array4& mask_v_arr); + +void +init_msfs_from_metgrid (amrex::FArrayBox& msfu_fab, + amrex::FArrayBox& msfv_fab, + amrex::FArrayBox& msfm_fab, + const int& flag_msfu, + const int& flag_msfv, + const int& flag_msfm, + const amrex::Vector& NC_MSFU_fab, + const amrex::Vector& NC_MSFV_fab, + const amrex::Vector& NC_MSFM_fab); + +void +init_base_state_from_metgrid (const amrex::Real l_rdOcp, + const amrex::Box& valid_bx, + const amrex::Vector& flag_psfc, + amrex::FArrayBox& state, + amrex::FArrayBox& r_hse_fab, + amrex::FArrayBox& p_hse_fab, + amrex::FArrayBox& pi_hse_fab, + amrex::FArrayBox& z_phys_nd_fab, + const amrex::Vector& NC_ght_fab, + const amrex::Vector& NC_psfc_fab, + amrex::Vector>& fabs_for_bcs, + const amrex::Array4& mask_c_arr); + +AMREX_FORCE_INLINE +AMREX_GPU_DEVICE +void +calc_rho_p (const int& kmax, + const int& flag_psfc, + const amrex::Real& psfc, + amrex::Real* Thetad_vec, + amrex::Real* Thetam_vec, +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + amrex::Real* Q_vec, +#endif + amrex::Real* z_vec, + amrex::Real* Rhod_vec, + amrex::Real* Rhom_vec, + amrex::Real* Pd_vec, + amrex::Real* Pm_vec) +{ + const int maxiter = 10; + + // Calculate or use moist pressure at the surface. + if (flag_psfc == 1) { + Pm_vec[0] = psfc; + } else { + amrex::Real t_0 = 290.0; // WRF's model_config_rec%base_temp + amrex::Real a = 50.0; // WRF's model_config_rec%base_lapse + Pm_vec[0] = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*CONST_GRAV*z_vec[0]/(a*R_d)), 0.5)); + } + + // calculate virtual potential temperature at the surface. + { +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[0]; +#else + amrex::Real qvf = 1.0; +#endif + Thetam_vec[0] = Thetad_vec[0]*qvf; + } + + // calculate moist density at the surface. + Rhom_vec[0] = 1.0/(R_d/p_0*Thetam_vec[0]*std::pow(Pm_vec[0]/p_0, -iGamma)); + + // integrate from the surface to the top boundary. + for (int k=1; k<=kmax; k++) { + amrex::Real dz = z_vec[k]-z_vec[k-1]; +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) + amrex::Real qvf = 1.0+(R_v/R_d+1.0)*Q_vec[k]; +#else + amrex::Real qvf = 1.0; +#endif + Thetam_vec[k] = Thetad_vec[k]*qvf; + Rhom_vec[k] = Rhom_vec[k-1]; // an initial guess. + for (int it=0; it=0; k--) { + amrex::Real dz = z_vec[k+1]-z_vec[k]; + Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess. + for (int it=0; it < maxiter; it++) { + Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV; + if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; + Rhod_vec[k] = 1.0/(R_d/p_0*Thetam_vec[k]*std::pow(Pd_vec[k]/p_0, -iGamma)); + } // it + } // k +} + +AMREX_FORCE_INLINE +AMREX_GPU_DEVICE +amrex::Real +interpolate_column_metgrid (const int& i, + const int& j, + const int& k, + char stag, + int src_comp, + const amrex::Array4& orig_z, + const amrex::Array4& orig_data, + const amrex::Array4& new_z) +{ + // This subroutine is a bit ham-handed and can be cleaned up later. + int imax_orig = amrex::ubound(amrex::Box(orig_data)).x; + int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y; + int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z; + + amrex::Real z; + if (stag == 'X') { + z = 0.25*(new_z(i,j,k)+new_z(i,j+1,k)+new_z(i,j,k+1)+new_z(i,j+1,k+1)); + } + else if (stag == 'Y') { + z = 0.25*(new_z(i,j,k)+new_z(i+1,j,k)+new_z(i,j,k+1)+new_z(i+1,j,k+1)); + } + else if (stag == 'M') { + z = 0.125*(new_z(i,j,k )+new_z(i,j+1,k )+new_z(i+1,j,k )+new_z(i+1,j+1,k )+ + new_z(i,j,k+1)+new_z(i,j+1,k+1)+new_z(i+1,j,k+1)+new_z(i+1,j+1,k+1)); + } + + amrex::Real z0, z1; + int klow = -1; + int khi0 = -1; + amrex::Real dzlow = 10000.0; + amrex::Real dzhi0 = -10000.0; + for (int kk = 0; kk < kmax_orig; kk++) { + amrex::Real orig_z_stag = 0.0; + if (stag == 'M') { + orig_z_stag = orig_z(i,j,kk); + } + if (stag == 'X') { + if (i == 0) { + orig_z_stag = orig_z(i,j,kk); + } + else if (i == imax_orig) { + orig_z_stag = orig_z(imax_orig-1,j,kk); + } + else { + orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk)); + } + } + else if (stag == 'Y') { + if (j == 0) { + orig_z_stag = orig_z(i,j,kk); + } + else if (j == jmax_orig) { + orig_z_stag = orig_z(i,jmax_orig-1,kk); + } + else { + orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk)); + } + } + amrex::Real dz = z - orig_z_stag; + if ((dz < 0.0) && (dz > dzhi0)) { + dzhi0 = dz; + khi0 = kk; + z1 = orig_z_stag; + } + if ((dz >= 0.0) && (dz < dzlow)) { + dzlow = dz; + klow = kk; + z0 = orig_z_stag; + } + } // kk + + if (klow == -1) { + // extrapolate + int khi1 = -1; + amrex::Real dzhi1 = -10000.0; + for (int kk = 0; kk < kmax_orig; kk++) { + amrex::Real orig_z_stag = 0.0; + if (stag == 'M') { + orig_z_stag = orig_z(i,j,kk); + } + else if (stag == 'X') { + if (i == 0) { + orig_z_stag = orig_z(i,j,kk); + } + else if (i == imax_orig) { + orig_z_stag = orig_z(imax_orig-1,j,kk); + } + else { + orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i-1,j,kk)); + } + } + else if (stag == 'Y') { + if (j == 0) { + orig_z_stag = orig_z(i,j,kk); + } + else if (j == jmax_orig) { + orig_z_stag = orig_z(i,jmax_orig-1,kk); + } + else { + orig_z_stag = 0.5*(orig_z(i,j,kk)+orig_z(i,j-1,kk)); + } + } + amrex::Real dz = z - orig_z_stag; + if ((dz < 0.0) && (dz > dzhi1) && (kk != khi0)) { + dzhi1 = dz; + khi1 = kk; + z1 = orig_z_stag; + } + } + amrex::Real y0 = orig_data(i,j,khi0,src_comp); + amrex::Real y1 = orig_data(i,j,khi1,src_comp); + return ( y0-(y1-y0)/(z1-z0)*(z0-z) ); + } else { + // interpolate + amrex::Real y0 = orig_data(i,j,klow,src_comp); + amrex::Real y1 = orig_data(i,j,khi0,src_comp); + return ( y0+(y1-y0)/(z1-z0)*(z-z0) ); + + } +} + +AMREX_FORCE_INLINE +AMREX_GPU_DEVICE +void +rh_to_mxrat (int i, + int j, + int k, + const amrex::Array4& rhum, + const amrex::Array4& temp, + const amrex::Array4& pres, + const amrex::Array4& mxrat) +{ + amrex::Real qv_max_p_safe = 10000.0; // WRF default value + amrex::Real qv_max_flag = 1.0e-5; // WRF default value + amrex::Real qv_max_value = 3.0e-6; // WRF default value + amrex::Real qv_min_p_safe = 110000.0; // WRF default value + amrex::Real qv_min_flag = 1.0e-6; // WRF default value + amrex::Real qv_min_value = 1.0e-6; // WRF default value + amrex::Real eps = 0.622; + amrex::Real svp1 = 0.6112; + amrex::Real svp2 = 17.67; + amrex::Real svp3 = 29.65; + amrex::Real svpt0 = 273.15; + // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior) + if (temp(i,j,k) != 0.0) { + amrex::Real es=0.01*rhum(i,j,k)*svp1*10.0*exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3)); + if (es >= pres(i,j,k)/100.0) { + // vapor pressure exceeds total pressure + mxrat(i,j,k) = std::pow(10.0, -6); + } + else { + mxrat(i,j,k) = amrex::max(eps*es/(pres(i,j,k)/100.0-es), 1.0e-6); + } + } + else { + // I don't know why there's a fringe case handled in WRF where T is absolute zero... + // Let's just deal with it here in case we also end up needing it. + mxrat(i,j,k) = 1.0e-6; + } + // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1. + // For pressures above a defined level, reasonable Qv values should be + // a certain value or smaller. If they are larger than this, the input data + // probably had "missing" RH, and we filled in some values. This is an + // attempt to catch those. Also, set the minimum value for the entire + // domain that is above the selected pressure level. + if (pres(i,j,k) < qv_max_p_safe) { + if (mxrat(i,j,k) > qv_max_flag) { + mxrat(i,j,k) = qv_max_value; + } + } + if (pres(i,j,k) < qv_min_p_safe) { + if (mxrat(i,j,k) < qv_min_flag) { + mxrat(i,j,k) = qv_min_value; + } + } +} +#endif diff --git a/Source/TimeIntegration/ERF_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index 51a3134f2..8cdd9be40 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -192,7 +192,7 @@ void erf_fast_rhs_T (int step, int /*level*/, const Array4 & stage_xmom = S_stage_data[IntVar::xmom].const_array(mfi); const Array4 & stage_ymom = S_stage_data[IntVar::ymom].const_array(mfi); -#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_MOISTURE) +#if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) const Array4 & prim = S_stage_prim.const_array(mfi); #endif diff --git a/Source/TimeIntegration/TI_slow_rhs_fun.H b/Source/TimeIntegration/TI_slow_rhs_fun.H index fd599247b..9be9eb6fc 100644 --- a/Source/TimeIntegration/TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/TI_slow_rhs_fun.H @@ -194,9 +194,17 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones - if (init_type=="real" && level==0) { - wrfbdy_compute_interior_ghost_rhs(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, - wrfbdy_width-1, wrfbdy_set_width, fine_geom, + if (((init_type == "real") || (init_type == "metgrid")) && level == 0) { + int width, set_width; + if (init_type == "real") { + width = wrfbdy_width; + set_width = wrfbdy_set_width; + } else if (init_type == "metgrid") { + width = metgrid_bdy_width; + set_width = metgrid_bdy_set_width; + } + wrfbdy_compute_interior_ghost_rhs(init_type, bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, + width-1, set_width, fine_geom, S_rhs, S_data, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi); } @@ -264,6 +272,7 @@ // Flag for moisture relaxation zone bool moist_zero = false; if (init_type=="real" && level==0 && wrfbdy_set_width>0) moist_zero = true; + if (init_type=="metgrid" && level==0 && metgrid_bdy_set_width > 0) moist_zero = true; #endif // Moving terrain @@ -341,9 +350,17 @@ #ifdef ERF_USE_NETCDF // Populate RHS for relaxation zones - if (init_type=="real" && level==0) { - wrfbdy_compute_interior_ghost_rhs(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, - wrfbdy_width-1, wrfbdy_set_width, fine_geom, + if (((init_type == "real") || (init_type == "metgrid")) && level == 0) { + int width, set_width; + if (init_type == "real") { + width = wrfbdy_width; + set_width = wrfbdy_set_width; + } else if (init_type == "metgrid") { + width = metgrid_bdy_width; + set_width = metgrid_bdy_set_width; + } + wrfbdy_compute_interior_ghost_rhs(init_type, bdy_time_interval, start_bdy_time, new_stage_time, slow_dt, + width-1, set_width, fine_geom, S_rhs, S_data, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi); } diff --git a/Source/Utils/InteriorGhostCells.cpp b/Source/Utils/InteriorGhostCells.cpp index 9ad271b00..4d9b5950e 100644 --- a/Source/Utils/InteriorGhostCells.cpp +++ b/Source/Utils/InteriorGhostCells.cpp @@ -89,6 +89,7 @@ compute_interior_ghost_bxs_xy (const Box& bx, /** * Compute the RHS in the relaxation zone * + * @param[in] init_type initialization method for this simulation * @param[in] bdy_time_interval time interval between boundary condition time stamps * @param[in] time current time * @param[in] delta_t timestep @@ -104,7 +105,8 @@ compute_interior_ghost_bxs_xy (const Box& bx, * @param[in] start_bdy_time time of the first boundary data read in */ void -wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, +wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, + const Real& bdy_time_interval, const Real& start_bdy_time, const Real& time, const Real& delta_t, @@ -126,8 +128,8 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, // Time interpolation Real dT = bdy_time_interval; - Real time_since_start = (time - start_bdy_time) / 1.e10; - int n_time = static_cast( time_since_start / dT); + Real time_since_start = time - start_bdy_time; + int n_time = static_cast( time_since_start / dT); amrex::Real alpha = (time_since_start - n_time * dT) / dT; AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); amrex::Real oma = 1.0 - alpha; @@ -145,13 +147,24 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, // Variable icomp map Vector comp_map = {0, 0, Rho_comp, RhoTheta_comp}; - // End of loop for WRFBdyVars - int WRFBdyEnd = WRFBdyVars::NumTypes-3; - + int BdyEnd, ivarU, ivarV, ivarR, ivarT; + if (init_type == "real") { + BdyEnd = WRFBdyVars::NumTypes-3; + ivarU = WRFBdyVars::U; + ivarV = WRFBdyVars::V; + ivarR = WRFBdyVars::R; + ivarT = WRFBdyVars::T; + } else if (init_type == "metgrid") { + BdyEnd = MetGridBdyVars::NumTypes-1; + ivarU = MetGridBdyVars::U; + ivarV = MetGridBdyVars::V; + ivarR = MetGridBdyVars::R; + ivarT = MetGridBdyVars::T; + } // Size the FABs //========================================================== - for (int ivar(WRFBdyVars::U); ivar < WRFBdyEnd; ivar++) + for (int ivar(ivarU); ivar < BdyEnd; ivar++) { int var_idx = var_map[ivar]; Box domain = geom.Domain(); @@ -167,16 +180,17 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, bx_ylo, bx_yhi, ng_vect, true); - if (ivar == WRFBdyVars::U) { + // Size the FABs + if (ivar == ivarU) { U_xlo.resize(bx_xlo,1); U_xhi.resize(bx_xhi,1); U_ylo.resize(bx_ylo,1); U_yhi.resize(bx_yhi,1); - } else if (ivar == WRFBdyVars::V) { + } else if (ivar == ivarV) { V_xlo.resize(bx_xlo,1); V_xhi.resize(bx_xhi,1); V_ylo.resize(bx_ylo,1); V_yhi.resize(bx_yhi,1); - } else if (ivar == WRFBdyVars::R) { + } else if (ivar == ivarR) { R_xlo.resize(bx_xlo,1); R_xhi.resize(bx_xhi,1); R_ylo.resize(bx_ylo,1); R_yhi.resize(bx_yhi,1); - } else if (ivar == WRFBdyVars::T){ + } else if (ivar == ivarT){ T_xlo.resize(bx_xlo,1); T_xhi.resize(bx_xhi,1); T_ylo.resize(bx_ylo,1); T_yhi.resize(bx_yhi,1); } else { @@ -197,10 +211,9 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, Elixir T_xlo_eli = T_xlo.elixir(); Elixir T_xhi_eli = T_xhi.elixir(); Elixir T_ylo_eli = T_ylo.elixir(); Elixir T_yhi_eli = T_yhi.elixir(); - // Populate FABs from boundary interpolation //========================================================== - for (int ivar(WRFBdyVars::U); ivar < WRFBdyEnd; ivar++) + for (int ivar(ivarU); ivar < BdyEnd; ivar++) { int var_idx = var_map[ivar]; Box domain = geom.Domain(); @@ -220,16 +233,16 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, Array4 arr_xlo; Array4 arr_xhi; Array4 arr_ylo; Array4 arr_yhi; - if (ivar == WRFBdyVars::U) { + if (ivar == ivarU) { arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); - } else if (ivar == WRFBdyVars::V) { + } else if (ivar == ivarV) { arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); - } else if (ivar == WRFBdyVars::R) { + } else if (ivar == ivarR) { arr_xlo = R_xlo.array(); arr_xhi = R_xhi.array(); arr_ylo = R_ylo.array(); arr_yhi = R_yhi.array(); - } else if (ivar == WRFBdyVars::T){ + } else if (ivar == ivarT){ arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array(); arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array(); } else { @@ -288,10 +301,9 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, }); } // ivar - // Velocity to momentum //========================================================== - for (int ivar(WRFBdyVars::U); ivar <= WRFBdyVars::V; ivar++) + for (int ivar(ivarU); ivar <= ivarV; ivar++) { int ivar_idx = ivar_map[ivar]; Box domain = geom.Domain(); @@ -312,7 +324,7 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, Array4 arr_xlo; Array4 arr_xhi; Array4 arr_ylo; Array4 arr_yhi; - if (ivar == WRFBdyVars::U) { + if (ivar == ivarU) { arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); @@ -425,10 +437,11 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, } } // mfi +// return; // DJW debugging return statement to shut off relaxation zone. // Compute RHS in relaxation region //========================================================== - for (int ivar(WRFBdyVars::U); ivar < WRFBdyEnd; ivar++) + for (int ivar(ivarU); ivar < BdyEnd; ivar++) { int ivar_idx = ivar_map[ivar]; int icomp = comp_map[ivar]; @@ -455,22 +468,22 @@ wrfbdy_compute_interior_ghost_rhs (const Real& bdy_time_interval, Array4 rhs_arr; Array4 data_arr; Array4 arr_xlo; Array4 arr_xhi; Array4 arr_ylo; Array4 arr_yhi; - if (ivar == WRFBdyVars::U) { + if (ivar == ivarU) { arr_xlo = U_xlo.array(); arr_xhi = U_xhi.array(); arr_ylo = U_ylo.array(); arr_yhi = U_yhi.array(); rhs_arr = S_rhs[IntVar::xmom].array(mfi); data_arr = S_data[IntVar::xmom].array(mfi); - } else if (ivar == WRFBdyVars::V) { + } else if (ivar == ivarV) { arr_xlo = V_xlo.array(); arr_xhi = V_xhi.array(); arr_ylo = V_ylo.array(); arr_yhi = V_yhi.array(); rhs_arr = S_rhs[IntVar::ymom].array(mfi); data_arr = S_data[IntVar::ymom].array(mfi); - } else if (ivar == WRFBdyVars::R) { + } else if (ivar == ivarR) { arr_xlo = R_xlo.array(); arr_xhi = R_xhi.array(); arr_ylo = R_ylo.array(); arr_yhi = R_yhi.array(); rhs_arr = S_rhs[IntVar::cons].array(mfi); data_arr = S_data[IntVar::cons].array(mfi); - } else if (ivar == WRFBdyVars::T){ + } else if (ivar == ivarT){ arr_xlo = T_xlo.array(); arr_xhi = T_xhi.array(); arr_ylo = T_ylo.array(); arr_yhi = T_yhi.array(); rhs_arr = S_rhs[IntVar::cons].array(mfi); diff --git a/Source/Utils/Utils.H b/Source/Utils/Utils.H index 6abc45a7d..e154a4946 100644 --- a/Source/Utils/Utils.H +++ b/Source/Utils/Utils.H @@ -68,7 +68,8 @@ void compute_interior_ghost_bxs_xy (const amrex::Box& bx, /* * Compute relaxation region RHS with wrfbdy */ -void wrfbdy_compute_interior_ghost_rhs (const amrex::Real& bdy_time_interval, +void wrfbdy_compute_interior_ghost_rhs (const std::string& init_type, + const amrex::Real& bdy_time_interval, const amrex::Real& start_bdy_time, const amrex::Real& time, const amrex::Real& delta_t, From ac2e9223b821fe23997d6f10be53afbb3704470e Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Thu, 9 Nov 2023 09:39:22 -0800 Subject: [PATCH 2/2] Unify epoch time and formatting. (#1287) Co-authored-by: Aaron Lattanzi --- Source/IO/NCInterface.H | 242 +++++++++++++++++----------------- Source/IO/NCWpsFile.H | 25 ++++ Source/IO/ReadFromMetgrid.cpp | 27 +--- Source/IO/ReadFromWRFBdy.cpp | 23 ---- 4 files changed, 149 insertions(+), 168 deletions(-) diff --git a/Source/IO/NCInterface.H b/Source/IO/NCInterface.H index 687d44a2a..386bd79e0 100644 --- a/Source/IO/NCInterface.H +++ b/Source/IO/NCInterface.H @@ -58,120 +58,128 @@ struct NCVar const int varid; //! Name of this variable - [[nodiscard]] std::string name() const; + [[nodiscard]] std::string name () const; //! Number of array dimensions for this variable - [[nodiscard]] int ndim() const; + [[nodiscard]] int ndim () const; //! Shape of the array (size in each array dimension) - [[nodiscard]] std::vector shape() const; + [[nodiscard]] std::vector shape () const; //! Write out the entire variable - void put(const double*) const; - void put(const float*) const; - void put(const int*) const; + void put (const double*) const; + void put (const float*) const; + void put (const int*) const; //! Write out a slice of data void - put(const double*, - const std::vector&, - const std::vector&) const; + put (const double*, + const std::vector&, + const std::vector&) const; //! Write out a slice of data with with strides (see hyperslab definition in //! NetCDF) void - put(const double*, - const std::vector&, - const std::vector&, - const std::vector&) const; + put (const double*, + const std::vector&, + const std::vector&, + const std::vector&) const; //! Write out a slice of data void - put(const float*, - const std::vector&, - const std::vector&) const; + put (const float*, + const std::vector&, + const std::vector&) const; //! Write out a slice of data with with strides (see hyperslab definition in //! NetCDF) void - put(const float*, - const std::vector&, - const std::vector&, - const std::vector&) const; + put (const float*, + const std::vector&, + const std::vector&, + const std::vector&) const; + + void put (const int*, const std::vector&, const std::vector&) const; - void put(const int*, const std::vector&, const std::vector&) const; void - put(const int*, - const std::vector&, - const std::vector&, - const std::vector&) const; + put (const int*, + const std::vector&, + const std::vector&, + const std::vector&) const; + + void put (const char**, const std::vector&, const std::vector&) const; - void put(const char**, const std::vector&, const std::vector&) const; void - put( - const char** dptr, - const std::vector& start, - const std::vector& count, - const std::vector& stride) const; + put (const char** dptr, + const std::vector& start, + const std::vector& count, + const std::vector& stride) const; //! Read the entire variable from file - void get(double*) const; - void get(float*) const; - void get(int*) const; + void get (double*) const; + void get (float*) const; + void get (int*) const; //! Read a chunk of data from the file void - get(double*, const std::vector&, const std::vector&) const; + get (double*, const std::vector&, const std::vector&) const; //! Read a chunk of data with strides void - get(double*, - const std::vector&, - const std::vector&, - const std::vector&) const; + get (double*, + const std::vector&, + const std::vector&, + const std::vector&) const; //! Read a chunk of data from the file void - get(float*, const std::vector&, const std::vector&) const; + get (float*, const std::vector&, const std::vector&) const; //! Read a chunk of data with strides void - get(float*, - const std::vector&, - const std::vector&, - const std::vector&) const; + get (float*, + const std::vector&, + const std::vector&, + const std::vector&) const; void - get(int*, const std::vector&, const std::vector&) const; + get (int*, const std::vector&, const std::vector&) const; + void - get(int*, - const std::vector&, - const std::vector&, - const std::vector&) const; + get (int*, + const std::vector&, + const std::vector&, + const std::vector&) const; //! Read a chunk of data from the file void - get(char*, const std::vector&, const std::vector&) const; + get (char*, const std::vector&, const std::vector&) const; //! Read a chunk of data with strides void - get(char*, - const std::vector&, - const std::vector&, - const std::vector&) const; + get (char*, + const std::vector&, + const std::vector&, + const std::vector&) const; + + [[nodiscard]] bool has_attr (const std::string& name) const; + + void put_attr (const std::string& name, const std::string& value) const; - [[nodiscard]] bool has_attr(const std::string& name) const; - void put_attr(const std::string& name, const std::string& value) const; void - put_attr(const std::string& name, const std::vector& value) const; + put_attr (const std::string& name, const std::vector& value) const; + void - put_attr(const std::string& name, const std::vector& value) const; - void put_attr(const std::string& name, const std::vector& value) const; - - [[nodiscard]] std::string get_attr(const std::string& name) const; - void get_attr(const std::string& name, std::vector& value) const; - void get_attr(const std::string& name, std::vector& value) const; - void get_attr(const std::string& name, std::vector& value) const; - void par_access(int cmode) const; //Uncomment for parallel NetCDF + put_attr (const std::string& name, const std::vector& value) const; + + void put_attr (const std::string& name, const std::vector& value) const; + + [[nodiscard]] std::string get_attr (const std::string& name) const; + + void get_attr (const std::string& name, std::vector& value) const; + void get_attr (const std::string& name, std::vector& value) const; + void get_attr (const std::string& name, std::vector& value) const; + + void par_access (int cmode) const; //Uncomment for parallel NetCDF }; //! Representation of a NetCDF group @@ -182,99 +190,97 @@ public: const int ncid; //! Name of this group - [[nodiscard]] std::string name() const; + [[nodiscard]] std::string name () const; //! Full name for this group - [[nodiscard]] std::string full_name() const; + [[nodiscard]] std::string full_name () const; //! Number of sub-groups within this group - [[nodiscard]] int num_groups() const; + [[nodiscard]] int num_groups () const; //! Number of dimensions in this group - [[nodiscard]] int num_dimensions() const; + [[nodiscard]] int num_dimensions () const; //! Number of variables within this group - [[nodiscard]] int num_variables() const; + [[nodiscard]] int num_variables () const; //! Number of attributes within this group - [[nodiscard]] int num_attributes() const; + [[nodiscard]] int num_attributes () const; //! Check if a group exists - [[nodiscard]] bool has_group(const std::string&) const; + [[nodiscard]] bool has_group (const std::string&) const; //! Check if a dimension exists by name - [[nodiscard]] bool has_dim(const std::string&) const; + [[nodiscard]] bool has_dim (const std::string&) const; //! Check if a variable exists by name - [[nodiscard]] bool has_var(const std::string&) const; + [[nodiscard]] bool has_var (const std::string&) const; //! Check if an attribute exists - [[nodiscard]] bool has_attr(const std::string&) const; + [[nodiscard]] bool has_attr (const std::string&) const; /** Get the group by name * * Throws error if the group doesn't exist, use `has_group` to check */ - [[nodiscard]] NCGroup group(const std::string&) const; + [[nodiscard]] NCGroup group (const std::string&) const; //! Get the dimension instance by name - [[nodiscard]] NCDim dim(const std::string&) const; + [[nodiscard]] NCDim dim (const std::string&) const; //! Get the variable instance by name - [[nodiscard]] NCVar var(const std::string&) const; + [[nodiscard]] NCVar var (const std::string&) const; //! Define new group - [[nodiscard]] NCGroup def_group(const std::string&) const; + [[nodiscard]] NCGroup def_group (const std::string&) const; //! Define new dimension - [[nodiscard]] NCDim def_dim(const std::string&, size_t len) const; + [[nodiscard]] NCDim def_dim (const std::string&, size_t len) const; //! Define a scalar variable, i.e., 0-dimensional array - [[nodiscard]] NCVar def_scalar(const std::string& name, nc_type dtype) const; + [[nodiscard]] NCVar def_scalar (const std::string& name, nc_type dtype) const; //! Define an array - [[nodiscard]] NCVar def_array( - const std::string& name, - nc_type dtype, - const std::vector&) const; + [[nodiscard]] NCVar def_array (const std::string& name, + nc_type dtype, + const std::vector&) const; //! Define a variable (wrapper for def_array) - [[nodiscard]] NCVar def_var( - const std::string& name, - const nc_type dtype, - const std::vector& dnames) const + [[nodiscard]] NCVar def_var (const std::string& name, + const nc_type dtype, + const std::vector& dnames) const { return def_array(name, dtype, dnames); } - void put_attr(const std::string& name, const std::string& value) const; - void put_attr(const std::string& name, const std::vector& value) const; - void put_attr(const std::string& name, const std::vector& value) const; - void put_attr(const std::string& name, const std::vector& value) const; + void put_attr (const std::string& name, const std::string& value) const; + void put_attr (const std::string& name, const std::vector& value) const; + void put_attr (const std::string& name, const std::vector& value) const; + void put_attr (const std::string& name, const std::vector& value) const; - [[nodiscard]] std::string get_attr(const std::string& name) const; - void get_attr(const std::string& name, std::vector& value) const; - void get_attr(const std::string& name, std::vector& value) const; - void get_attr(const std::string& name, std::vector& value) const; + [[nodiscard]] std::string get_attr (const std::string& name) const; + void get_attr (const std::string& name, std::vector& value) const; + void get_attr (const std::string& name, std::vector& value) const; + void get_attr (const std::string& name, std::vector& value) const; //! Return a list of all groups defined in this group - [[nodiscard]] std::vector all_groups() const; + [[nodiscard]] std::vector all_groups () const; //! Return a list of all dimensions defined in this group - [[nodiscard]] std::vector all_dims() const; + [[nodiscard]] std::vector all_dims () const; //! Return a list of all variables defined in this group - [[nodiscard]] std::vector all_vars() const; + [[nodiscard]] std::vector all_vars () const; //! Enter definition mode (not needed for NetCDF4 format) - void enter_def_mode() const; + void enter_def_mode () const; //! Exit definition mode - void exit_def_mode() const; + void exit_def_mode () const; protected: - NCGroup(const int id) : ncid(id) {} - NCGroup(const int id, const NCGroup* /*par*/) : ncid(id) {} + NCGroup (const int id) : ncid(id) {} + NCGroup (const int id, const NCGroup* /*par*/) : ncid(id) {} }; /** Representation of a NetCDF file @@ -285,28 +291,26 @@ class NCFile : public NCGroup { public: static NCFile - create(const std::string& name, const int cmode = NC_CLOBBER | NC_NETCDF4); + create (const std::string& name, const int cmode = NC_CLOBBER | NC_NETCDF4); - static NCFile open(const std::string& name, const int cmode = NC_NOWRITE); + static NCFile open (const std::string& name, const int cmode = NC_NOWRITE); - static NCFile create_par( - const std::string& name, - const int cmode = NC_CLOBBER | NC_NETCDF4 | NC_MPIIO, - MPI_Comm comm = MPI_COMM_WORLD, - MPI_Info info = MPI_INFO_NULL); + static NCFile create_par (const std::string& name, + const int cmode = NC_CLOBBER | NC_NETCDF4 | NC_MPIIO, + MPI_Comm comm = MPI_COMM_WORLD, + MPI_Info info = MPI_INFO_NULL); - static NCFile open_par( - const std::string& name, - const int cmode = NC_NOWRITE, - MPI_Comm comm = MPI_COMM_WORLD, - MPI_Info info = MPI_INFO_NULL); + static NCFile open_par (const std::string& name, + const int cmode = NC_NOWRITE, + MPI_Comm comm = MPI_COMM_WORLD, + MPI_Info info = MPI_INFO_NULL); - ~NCFile(); + ~NCFile (); - void close(); + void close (); protected: - NCFile(const int id) : NCGroup(id), is_open{true} {} + NCFile (const int id) : NCGroup(id), is_open{true} {} bool is_open{false}; }; diff --git a/Source/IO/NCWpsFile.H b/Source/IO/NCWpsFile.H index 7befc00cd..5e4994d64 100644 --- a/Source/IO/NCWpsFile.H +++ b/Source/IO/NCWpsFile.H @@ -128,6 +128,31 @@ int BuildFABsFromWRFBdyFile (const std::string &fname, amrex::Vector>& bdy_data_ylo, amrex::Vector>& bdy_data_yhi); +AMREX_FORCE_INLINE +std::time_t +getEpochTime (const std::string& dateTime, const std::string& dateTimeFormat) +{ + // Create a stream which we will use to parse the string, + // which we provide to constructor of stream to fill the buffer. + std::istringstream ss{ dateTime }; + + // Create a tm object to store the parsed date and time. + std::tm tmTime; + memset(&tmTime, 0, sizeof(tmTime)); + + // Now we read from buffer using get_time manipulator + // and formatting the input appropriately. + strptime(dateTime.c_str(), dateTimeFormat.c_str(), &tmTime); + + // Convert the tm structure to time_t value and return. + // Here we use timegm since the output should be relative to UTC. + auto epoch = timegm(&tmTime); + // Print() << "Time Stamp: "<< std::put_time(&tmTime, "%c") + // << " , Epoch: " << epoch << std::endl; + + return epoch; +} + template void ReadNetCDFFile (const std::string& fname, amrex::Vector names, amrex::Vector >& arrays) diff --git a/Source/IO/ReadFromMetgrid.cpp b/Source/IO/ReadFromMetgrid.cpp index 80748f083..43a82fb7d 100644 --- a/Source/IO/ReadFromMetgrid.cpp +++ b/Source/IO/ReadFromMetgrid.cpp @@ -6,31 +6,6 @@ using namespace amrex; #ifdef ERF_USE_NETCDF -// Converts UTC time string to a time_t value in UTC. -std::time_t -getEpochTime_metgrid (const std::string& dateTime, const std::string& dateTimeFormat) -{ - // Create a stream which we will use to parse the string, - // which we provide to constructor of stream to fill the buffer. - std::istringstream ss{ dateTime }; - - // Create a tm object to store the parsed date and time. - std::tm tmTime; - memset(&tmTime, 0, sizeof(tmTime)); - - // Now we read from buffer using get_time manipulator - // and formatting the input appropriately. - strptime(dateTime.c_str(), dateTimeFormat.c_str(), &tmTime); - - // Convert the tm structure to time_t value and return. - // Here we use timegm since the output should be relative to UTC. - auto epoch = timegm(&tmTime); - // Print() << "Time Stamp: "<< std::put_time(&tmTime, "%c") - // << " , Epoch: " << epoch << std::endl; - - return epoch; -} - void read_from_metgrid (int lev, const Box& domain, const std::string& fname, std::string& NC_dateTime, Real& NC_epochTime, @@ -65,7 +40,7 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, { // Global Attributes (string) NC_dateTime = ncf.get_attr("SIMULATION_START_DATE")+"UTC"; const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S%Z"; - NC_epochTime = getEpochTime_metgrid(NC_dateTime, dateTimeFormat); + NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); } { // Global Attributes (Real) std::vector attr; diff --git a/Source/IO/ReadFromWRFBdy.cpp b/Source/IO/ReadFromWRFBdy.cpp index ad78f6a89..434d190d1 100644 --- a/Source/IO/ReadFromWRFBdy.cpp +++ b/Source/IO/ReadFromWRFBdy.cpp @@ -26,29 +26,6 @@ namespace WRFBdyTypes { }; } -// Converts UTC time string to a time_t value. -std::time_t getEpochTime (const std::string& dateTime, const std::string& dateTimeFormat) -{ - // Create a stream which we will use to parse the string, - // which we provide to constructor of stream to fill the buffer. - std::istringstream ss{ dateTime }; - - // Create a tm object to store the parsed date and time. - std::tm tmTime; - memset(&tmTime, 0, sizeof(tmTime)); - - // Now we read from buffer using get_time manipulator - // and formatting the input appropriately. - strptime(dateTime.c_str(), dateTimeFormat.c_str(), &tmTime); - - // Convert the tm structure to time_t value and return. - auto epoch = std::mktime(&tmTime); - // Print() << "Time Stamp: "<< std::put_time(&tmTime, "%c") - // << " , Epoch: " << epoch << std::endl; - - return epoch; -} - Real read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, Vector>& bdy_data_xlo,