From 47af535dcc17d6f786f004b2d0194b6d12abb672 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Tue, 28 Apr 2026 20:45:06 -0700 Subject: [PATCH 01/27] WIP storing real boundary data in a "erfbdy" file of the AMReX format. This allows for the minimum number of boundary data times to be stored in memory. Stored boundary times are updated by reading from the boundary file as-needed during the run. This capability enables ERF to be run as a preprocessor by setting max_steps=0, which will initialize the simulation and produce checkpoint and boundary files. These features are working in serial for init_type="metgrid". --- Exec/GNUmakefile | 3 +- Source/ERF.H | 6 + Source/ERF.cpp | 8 +- Source/IO/ERF_Checkpoint.cpp | 35 +++++ Source/IO/ERF_ReadFromERFBdy.H | 28 ++++ Source/IO/ERF_ReadFromERFBdy.cpp | 137 +++++++++++++++++ Source/IO/ERF_WriteERFBdy.H | 25 +++ Source/IO/ERF_WriteERFBdy.cpp | 125 +++++++++++++++ Source/IO/Make.package | 5 + Source/Initialization/ERF_InitFromMetgrid.cpp | 71 +++++++++ .../Initialization/ERF_InitFromWRFInput.cpp | 143 +++++++++++++++--- Source/TimeIntegration/ERF_TimeStep.cpp | 89 +++++++---- 12 files changed, 618 insertions(+), 57 deletions(-) create mode 100644 Source/IO/ERF_ReadFromERFBdy.H create mode 100644 Source/IO/ERF_ReadFromERFBdy.cpp create mode 100644 Source/IO/ERF_WriteERFBdy.H create mode 100644 Source/IO/ERF_WriteERFBdy.cpp diff --git a/Exec/GNUmakefile b/Exec/GNUmakefile index 098a539412..100bf8d3c8 100644 --- a/Exec/GNUmakefile +++ b/Exec/GNUmakefile @@ -1,6 +1,7 @@ # AMReX COMP = gnu -PRECISION = DOUBLE +#PRECISION = DOUBLE +PRECISION = SINGLE # Profiling PROFILE = FALSE diff --git a/Source/ERF.H b/Source/ERF.H index 08bc5eb4fe..84b4e6fa51 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -1269,6 +1269,12 @@ private: int metgrid_order{2}; int metgrid_force_sfc_k{6}; + // Options for ERF boundary files (wrfinput/metgrid) + bool write_erfbdy{false}; + bool use_erfbdy{false}; + std::string erfbdy_file{"erfbdy"}; + int nvars_erfbdy{0}; + amrex::Vector ba1d; amrex::Vector ba2d; std::unique_ptr mf_C1H; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index a6ae009a07..45ace7dca1 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -52,6 +52,7 @@ Real ERF::final_low_time = -one; Real ERF::bdy_time_interval = std::numeric_limits::max(); Real ERF::low_time_interval = std::numeric_limits::max(); + #endif // Time step control @@ -1102,7 +1103,7 @@ ERF::InitData_post () // This follows init_from_wrfinput() // bool use_moist = (solverChoice.moisture_type != MoistureType::None); - if (solverChoice.use_real_bcs) { + if (solverChoice.use_real_bcs && solverChoice.init_type == InitType::WRFInput) { if ( geom[0].isPeriodic(0) || geom[0].isPeriodic(1) ) { amrex::Error("Cannot set periodic lateral boundary conditions when reading in real boundary values"); @@ -2500,6 +2501,11 @@ ERF::ReadParameters () pp.query("metgrid_order", metgrid_order); pp.query("metgrid_force_sfc_k", metgrid_force_sfc_k); + // Options for boundary file. + pp.query("write_erfbdy", write_erfbdy); + pp.query("use_erfbdy", use_erfbdy); + pp.query("erfbdy_file", erfbdy_file); + // Set default to FullState for now ... later we will try Perturbation interpolation_type = StateInterpType::FullState; pp.query_enum_case_insensitive("interpolation_type" ,interpolation_type); diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index fdcc72bceb..762ba64943 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -6,6 +6,7 @@ #include "ERF.H" #include "AMReX_PlotFileUtil.H" +#include "ERF_ReadFromERFBdy.H" using namespace amrex; @@ -1067,6 +1068,40 @@ ERF::ReadCheckpointFile () } // init_type == WRFInput or Metgrid #endif #endif + + // Load boundary data from erfbdy. + if (((solverChoice.init_type == InitType::WRFInput) || (solverChoice.init_type == InitType::Metgrid)) && + solverChoice.use_real_bcs && use_erfbdy) { + Print() << "Restart: Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; + + // Read metadata and times from erfbdy. + int ntimes_erfbdy; + Vector bdy_times; + bdy_time_interval = read_times_from_erfbdy(erfbdy_file, + ntimes_erfbdy, nvars_erfbdy, real_width, + bdy_times, start_bdy_time, final_bdy_time); + + Print() << "Restart: erfbdy file contains " << ntimes_erfbdy << " times" << std::endl; + + bdy_data_xlo.resize(ntimes_erfbdy); + bdy_data_xhi.resize(ntimes_erfbdy); + bdy_data_ylo.resize(ntimes_erfbdy); + bdy_data_yhi.resize(ntimes_erfbdy); + + // Determine which times we need based on current simulation time. + Real time_since_start_bdy = t_new[0] + start_time - start_bdy_time; + int n_time_old = std::min(static_cast(time_since_start_bdy / bdy_time_interval), ntimes_erfbdy-1); + int n_time_new = n_time_old + 1; + + // Load the necessary time slices. + for (int itime = n_time_old; itime <= std::min(n_time_new + 1, ntimes_erfbdy - 1); ++itime) { + read_from_erfbdy(itime, erfbdy_file, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + nvars_erfbdy, real_width); + Print() << "Restart: Loaded erfbdy time index " << itime << std::endl; + } + } } /** diff --git a/Source/IO/ERF_ReadFromERFBdy.H b/Source/IO/ERF_ReadFromERFBdy.H new file mode 100644 index 0000000000..2923d9c42a --- /dev/null +++ b/Source/IO/ERF_ReadFromERFBdy.H @@ -0,0 +1,28 @@ +#ifndef ERF_READ_FROM_ERFBDY_H_ +#define ERF_READ_FROM_ERFBDY_H_ + +#include +#include + +// Read metadata and time values from erfbdy Header +// Returns bdy_time_interval +amrex::Real +read_times_from_erfbdy(const std::string& bdy_file_name, + int& ntimes, + int& nvars, + int& real_width, + amrex::Vector& bdy_times, + amrex::Real& start_bdy_time, + amrex::Real& final_bdy_time); + +// Read a single time slice from erfbdy file +void +read_from_erfbdy(int itime, + const std::string& bdy_file_name, + amrex::Vector>& bdy_data_xlo, + amrex::Vector>& bdy_data_xhi, + amrex::Vector>& bdy_data_ylo, + amrex::Vector>& bdy_data_yhi, + int nvars, int real_width); + +#endif diff --git a/Source/IO/ERF_ReadFromERFBdy.cpp b/Source/IO/ERF_ReadFromERFBdy.cpp new file mode 100644 index 0000000000..7c53ce3e85 --- /dev/null +++ b/Source/IO/ERF_ReadFromERFBdy.cpp @@ -0,0 +1,137 @@ +#include "ERF_ReadFromERFBdy.H" +#include +#include +#include +#include + +using namespace amrex; + +Real +read_times_from_erfbdy(const std::string& bdy_file_name, + int& ntimes, + int& nvars, + int& real_width, + Vector& bdy_times, + Real& start_bdy_time, + Real& final_bdy_time) +{ + std::string HeaderFileName = bdy_file_name + "/Header"; + + // Read header file + std::ifstream HeaderFile; + HeaderFile.open(HeaderFileName.c_str(), std::ifstream::in); + + if (!HeaderFile.good()) { + amrex::FileOpenFailed(HeaderFileName); + } + + std::string line; + + // Read title line + std::getline(HeaderFile, line); + + // Read metadata + HeaderFile >> ntimes; + HeaderFile >> nvars; + HeaderFile >> real_width; + + // Read time values + bdy_times.resize(ntimes); + for (int i = 0; i < ntimes; ++i) { + HeaderFile >> bdy_times[i]; + } + + // Read domain box (stored but not used in this function) + int sml[3], big[3]; + HeaderFile >> sml[0] >> sml[1] >> sml[2]; + HeaderFile >> big[0] >> big[1] >> big[2]; + + HeaderFile.close(); + + // Set start and final times + start_bdy_time = bdy_times[0]; + final_bdy_time = bdy_times[ntimes - 1]; + + // Calculate time interval + Real bdy_time_interval = (ntimes > 1) ? (bdy_times[1] - bdy_times[0]) : 0.0; + + return bdy_time_interval; +} + +void +read_from_erfbdy(int itime, + const std::string& bdy_file_name, + Vector>& bdy_data_xlo, + Vector>& bdy_data_xhi, + Vector>& bdy_data_ylo, + Vector>& bdy_data_yhi, + int nvars, int real_width) +{ + // Construct time subdirectory path + std::string time_dir = bdy_file_name + "/Time_" + Concatenate("", itime, 6); + + // Ensure vectors are sized correctly + if (bdy_data_xlo[itime].empty()) { + bdy_data_xlo[itime].resize(nvars); + bdy_data_xhi[itime].resize(nvars); + bdy_data_ylo[itime].resize(nvars); + bdy_data_yhi[itime].resize(nvars); + } + + // Read each variable for each boundary direction using FArrayBox::readFrom + if (ParallelDescriptor::IOProcessor()) + { + for (int ivar = 0; ivar < nvars; ++ivar) + { + // X-low boundary + { + std::string filename = time_dir + "/BdyData_xlo_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_xlo[itime][ivar].readFrom(ifs); + ifs.close(); + } + + // X-high boundary + { + std::string filename = time_dir + "/BdyData_xhi_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_xhi[itime][ivar].readFrom(ifs); + ifs.close(); + } + + // Y-low boundary + { + std::string filename = time_dir + "/BdyData_ylo_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_ylo[itime][ivar].readFrom(ifs); + ifs.close(); + } + + // Y-high boundary + { + std::string filename = time_dir + "/BdyData_yhi_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_yhi[itime][ivar].readFrom(ifs); + ifs.close(); + } + } + } + + // Broadcast data to all processors + ParallelDescriptor::Barrier(); + for (int ivar = 0; ivar < nvars; ++ivar) + { + ParallelDescriptor::Bcast(bdy_data_xlo[itime][ivar].dataPtr(), + bdy_data_xlo[itime][ivar].size(), + ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Bcast(bdy_data_xhi[itime][ivar].dataPtr(), + bdy_data_xhi[itime][ivar].size(), + ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Bcast(bdy_data_ylo[itime][ivar].dataPtr(), + bdy_data_ylo[itime][ivar].size(), + ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::Bcast(bdy_data_yhi[itime][ivar].dataPtr(), + bdy_data_yhi[itime][ivar].size(), + ParallelDescriptor::IOProcessorNumber()); + } +} diff --git a/Source/IO/ERF_WriteERFBdy.H b/Source/IO/ERF_WriteERFBdy.H new file mode 100644 index 0000000000..d35d6b30b2 --- /dev/null +++ b/Source/IO/ERF_WriteERFBdy.H @@ -0,0 +1,25 @@ +#ifndef ERF_WRITE_ERFBDY_H_ +#define ERF_WRITE_ERFBDY_H_ + +#include +#include +#include + +// Initialize erfbdy file and write header (call once at start) +void InitERFBdyFile(const std::string& bdy_file_name, + int ntimes, + const amrex::Vector& bdy_times, + const amrex::Box& domain, + int nvars, + int real_width); + +// Write a single time slice (call once per time) +void WriteERFBdyTimeSlice(const std::string& bdy_file_name, + int itime, + const amrex::Vector& bdy_data_xlo, + const amrex::Vector& bdy_data_xhi, + const amrex::Vector& bdy_data_ylo, + const amrex::Vector& bdy_data_yhi, + int nvars); + +#endif diff --git a/Source/IO/ERF_WriteERFBdy.cpp b/Source/IO/ERF_WriteERFBdy.cpp new file mode 100644 index 0000000000..01e19ddccf --- /dev/null +++ b/Source/IO/ERF_WriteERFBdy.cpp @@ -0,0 +1,125 @@ +#include "ERF_WriteERFBdy.H" +#include +#include +#include +#include + +using namespace amrex; + +void InitERFBdyFile(const std::string& bdy_file_name, + int ntimes, + const Vector& bdy_times, + const Box& domain, + int nvars, + int real_width) +{ + // Only I/O processor creates the directory and writes header + if (ParallelDescriptor::IOProcessor()) + { + // Create erfbdy directory + if (!amrex::UtilCreateDirectory(bdy_file_name, 0755)) { + amrex::CreateDirectoryFailed(bdy_file_name); + } + + // Write Header file + std::string HeaderFileName = bdy_file_name + "/Header"; + std::ofstream HeaderFile; + HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out | std::ofstream::trunc); + + if (!HeaderFile.good()) { + amrex::FileOpenFailed(HeaderFileName); + } + + HeaderFile.precision(17); + + // Title + HeaderFile << "ERFBdy file for ERF\n"; + + // Metadata + HeaderFile << ntimes << "\n"; + HeaderFile << nvars << "\n"; + HeaderFile << real_width << "\n"; + + // Time values + for (int i = 0; i < ntimes; ++i) { + HeaderFile << bdy_times[i]; + if (i < ntimes - 1) { + HeaderFile << " "; + } + } + HeaderFile << "\n"; + + // Domain box + HeaderFile << domain.smallEnd(0) << " " << domain.smallEnd(1) << " " << domain.smallEnd(2) << "\n"; + HeaderFile << domain.bigEnd(0) << " " << domain.bigEnd(1) << " " << domain.bigEnd(2) << "\n"; + + HeaderFile.close(); + } + + // Barrier to ensure directory is created before any writes + ParallelDescriptor::Barrier(); +} + +void WriteERFBdyTimeSlice(const std::string& bdy_file_name, + int itime, + const Vector& bdy_data_xlo, + const Vector& bdy_data_xhi, + const Vector& bdy_data_ylo, + const Vector& bdy_data_yhi, + int nvars) +{ + // Create time subdirectory + std::string time_dir = bdy_file_name + "/Time_" + Concatenate("", itime, 6); + + if (ParallelDescriptor::IOProcessor()) + { + if (!amrex::UtilCreateDirectory(time_dir, 0755)) { + amrex::CreateDirectoryFailed(time_dir); + } + } + + // Barrier to ensure directory exists + ParallelDescriptor::Barrier(); + + // Write each variable for each boundary direction using FArrayBox::writeOn + if (ParallelDescriptor::IOProcessor()) + { + for (int ivar = 0; ivar < nvars; ++ivar) + { + // X-low boundary + if (bdy_data_xlo.size() > 0 && ivar < bdy_data_xlo.size()) { + std::string filename = time_dir + "/BdyData_xlo_var" + std::to_string(ivar); + std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); + bdy_data_xlo[ivar].writeOn(ofs); + ofs.close(); + } + + // X-high boundary + if (bdy_data_xhi.size() > 0 && ivar < bdy_data_xhi.size()) { + std::string filename = time_dir + "/BdyData_xhi_var" + std::to_string(ivar); + std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); + bdy_data_xhi[ivar].writeOn(ofs); + ofs.close(); + } + + // Y-low boundary + if (bdy_data_ylo.size() > 0 && ivar < bdy_data_ylo.size()) { + std::string filename = time_dir + "/BdyData_ylo_var" + std::to_string(ivar); + std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); + bdy_data_ylo[ivar].writeOn(ofs); + ofs.close(); + } + + // Y-high boundary + if (bdy_data_yhi.size() > 0 && ivar < bdy_data_yhi.size()) { + std::string filename = time_dir + "/BdyData_yhi_var" + std::to_string(ivar); + std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); + bdy_data_yhi[ivar].writeOn(ofs); + ofs.close(); + } + } + } + + // Barrier to ensure all writes complete + ParallelDescriptor::Barrier(); +} diff --git a/Source/IO/Make.package b/Source/IO/Make.package index f5fb45b12a..276d26573c 100644 --- a/Source/IO/Make.package +++ b/Source/IO/Make.package @@ -16,6 +16,11 @@ CEXE_sources += ERF_WriteSubvolume.cpp CEXE_sources += ERF_ConsoleIO.cpp +CEXE_headers += ERF_WriteERFBdy.H +CEXE_sources += ERF_WriteERFBdy.cpp +CEXE_headers += ERF_ReadFromERFBdy.H +CEXE_sources += ERF_ReadFromERFBdy.cpp + CEXE_headers += ERF_SampleData.H ifeq ($(USE_NETCDF), TRUE) diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 5d7d070e01..5373085c18 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -3,6 +3,8 @@ */ #include #include +#include +#include using namespace amrex; @@ -53,6 +55,45 @@ ERF::init_from_metgrid (int lev) Print() << "Init with met_em without moisture model." << std::endl; } + // If using erfbdy file, load boundary data and then skip the rest of met_em processing. + if (lev == 0 && use_erfbdy) { + Print() << "Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; + + int ntimes_erfbdy; + Vector bdy_times; + bdy_time_interval = read_times_from_erfbdy(erfbdy_file, + ntimes_erfbdy, nvars_erfbdy, real_width, + bdy_times, start_bdy_time, final_bdy_time); + + AMREX_ALWAYS_ASSERT(ntimes_erfbdy >= 2); + + Print() << "erfbdy file has " << ntimes_erfbdy << " times" << std::endl; + Print() << "start_bdy_time = " << start_bdy_time << std::endl; + Print() << "final_bdy_time = " << final_bdy_time << std::endl; + Print() << "bdy_time_interval = " << bdy_time_interval << std::endl; + + bdy_data_xlo.resize(ntimes_erfbdy); + bdy_data_xhi.resize(ntimes_erfbdy); + bdy_data_ylo.resize(ntimes_erfbdy); + bdy_data_yhi.resize(ntimes_erfbdy); + + // Load the first 2 times. + for (int itime = 0; itime < 2; ++itime) { + read_from_erfbdy(itime, erfbdy_file, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + nvars_erfbdy, real_width); + Print() << "Loaded erfbdy time slice " << itime << std::endl; + } + + // Set the simulation start time. + t_new[lev] = zero; + t_old[lev] = -Real(1.e200); + + Print() << "Loaded boundaries from erfbdy, skipping met_em processing" << std::endl; + return; // Skip the rest of met_em processing. + } + int ntimes = num_files_at_level[lev]; Print() << ntimes << " met_em.d0" << lev+1 << "*.nc files are listed" << std::endl; @@ -223,6 +264,8 @@ ERF::init_from_metgrid (int lev) MultiFab th_hse(base_state[lev], make_alias, BaseState::th0_comp, 1); MultiFab qv_hse(base_state[lev], make_alias, BaseState::qv0_comp, 1); + Vector bdy_times(ntimes); + for (int itime(0); itime < ntimes; itime++) { Print() << " init_from_metgrid: reading nc_init_file[" << lev << "][" << itime << "]\t" << nc_init_file[lev][itime] << std::endl; read_from_metgrid(lev, itime, @@ -239,6 +282,9 @@ ERF::init_from_metgrid (int lev) NC_lmask_iab, geom[lev]); if (lev == 0) { + // Collect time for erfbdy + bdy_times[itime] = NC_epochTime[itime]; + if (itime == 0) { // Start at the earliest time in nc_init_file[lev]. start_bdy_time = NC_epochTime[itime]; @@ -250,6 +296,13 @@ ERF::init_from_metgrid (int lev) << " from metgrid file but note that time variable in simulation is elapsed time" << std::endl; t_new[lev] = zero; t_old[lev] = -Real(1.e200); + + // Initialize erfbdy file header now that we know domain and real_width + if (write_erfbdy) { + InitERFBdyFile(erfbdy_file, ntimes, bdy_times, + geom[lev].Domain(), MetGridBdyVars::NumTypes, real_width); + Print() << "Initialized erfbdy file: " << erfbdy_file << std::endl; + } } else { // Verify that files in nc_init_file[lev] are ordered from earliest to latest. AMREX_ALWAYS_ASSERT(NC_epochTime[itime] > NC_epochTime[itime-1]); @@ -553,6 +606,24 @@ ERF::init_from_metgrid (int lev) bdy_data_yhi[itime][nvar].size(), ParallelContext::CommunicatorAll()); } // nvar + + // Write this time slice to erfbdy immediately to avoid memory accumulation + if (write_erfbdy) { + WriteERFBdyTimeSlice(erfbdy_file, itime, + bdy_data_xlo[itime], bdy_data_xhi[itime], + bdy_data_ylo[itime], bdy_data_yhi[itime], + MetGridBdyVars::NumTypes); + Print() << "Wrote erfbdy time slice " << itime << " of " << ntimes-1 << std::endl; + + // CRITICAL: Clear this time slice from memory after writing + // (unless it's one of the first two times needed for simulation start) + if (itime > 1) { + bdy_data_xlo[itime].clear(); + bdy_data_xhi[itime].clear(); + bdy_data_ylo[itime].clear(); + bdy_data_yhi[itime].clear(); + } + } } // itime } // lev==0 diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index 9011b07f1b..796b0b44e5 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -11,6 +11,8 @@ #include #include +#include +#include using namespace amrex; @@ -1034,36 +1036,129 @@ ERF::init_from_wrfinput (int lev, amrex::Error("Cannot set periodic lateral boundary conditions when reading in real boundary values"); } - if (nc_bdy_file.empty()) { - amrex::Error("NetCDF boundary file name must be provided via input"); + // Path 1: Load from existing erfbdy file + if (use_erfbdy) { + Print() << "Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; + + // Read metadata and times from erfbdy + int ntimes_erfbdy; + Vector bdy_times; + bdy_time_interval = read_times_from_erfbdy(erfbdy_file, + ntimes_erfbdy, nvars_erfbdy, real_width, + bdy_times, start_bdy_time, final_bdy_time); + + Print() << "erfbdy file contains " << ntimes_erfbdy << " time slices" << std::endl; + Print() << "start_bdy_time = " << start_bdy_time << std::endl; + Print() << "final_bdy_time = " << final_bdy_time << std::endl; + Print() << "bdy_time_interval = " << bdy_time_interval << std::endl; + + // Resize boundary data structures + bdy_data_xlo.resize(ntimes_erfbdy); + bdy_data_xhi.resize(ntimes_erfbdy); + bdy_data_ylo.resize(ntimes_erfbdy); + bdy_data_yhi.resize(ntimes_erfbdy); + + // Load first 2 time slices for simulation start + for (int itime = 0; itime < std::min(2, ntimes_erfbdy); ++itime) { + read_from_erfbdy(itime, erfbdy_file, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + nvars_erfbdy, real_width); + Print() << "Loaded erfbdy time slice " << itime << std::endl; + } + + Print() << "Read in boundary data with width " << real_width << std::endl; + Print() << "Running with relaxation width: " << real_width << std::endl; } + // Path 2: Load from wrfbdy and optionally write to erfbdy + else { + if (nc_bdy_file.empty()) { + amrex::Error("NetCDF boundary file name must be provided via input"); + } - bdy_time_interval = read_times_from_wrfbdy(nc_bdy_file, - bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, - start_bdy_time, final_bdy_time); + bdy_time_interval = read_times_from_wrfbdy(nc_bdy_file, + bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, + start_bdy_time, final_bdy_time); - // ******************************************************************************************* - // We intentionally only read in the first three slices here ... we will read the rest in - // as needed during the time stepping procedure - // ******************************************************************************************* - int ntimes = bdy_data_xlo.size(); ntimes = amrex::min(ntimes, 3); - for (int itime = 0; itime < ntimes; itime++) - { - read_from_wrfbdy(itime, nc_bdy_file, geom[0].Domain(), - bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, - real_width); + int ntimes_total = bdy_data_xlo.size(); + Vector bdy_times(ntimes_total); - if (itime == 0) { - Print() << "Read in boundary data with width " << real_width << std::endl; - Print() << "Running with relaxation width: " << real_width << std::endl; + // Initialize erfbdy file if writing is enabled + if (write_erfbdy) { + // Collect times (we'll need to read them from wrfbdy to fill this properly) + // For now, construct from start_bdy_time and bdy_time_interval + for (int itime = 0; itime < ntimes_total; ++itime) { + bdy_times[itime] = start_bdy_time + itime * bdy_time_interval; + } + + InitERFBdyFile(erfbdy_file, ntimes_total, bdy_times, + geom[lev].Domain(), WRFBdyVars::NumTypes, real_width); + Print() << "Initialized erfbdy file: " << erfbdy_file << std::endl; } - convert_all_wrfbdy_data(itime, geom[lev].Domain(), - bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, - mf_MUB_lev, mf_C1H_lev, mf_C2H_lev, - lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], - geom[lev], use_moist); - } // itime + // ******************************************************************************************* + // We intentionally only read in the first three slices here ... we will read the rest in + // as needed during the time stepping procedure + // ******************************************************************************************* + int ntimes = bdy_data_xlo.size(); ntimes = amrex::min(ntimes, 3); + for (int itime = 0; itime < ntimes; itime++) + { + read_from_wrfbdy(itime, nc_bdy_file, geom[0].Domain(), + bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, + real_width); + + if (itime == 0) { + Print() << "Read in boundary data with width " << real_width << std::endl; + Print() << "Running with relaxation width: " << real_width << std::endl; + } + + convert_all_wrfbdy_data(itime, geom[lev].Domain(), + bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, + mf_MUB_lev, mf_C1H_lev, mf_C2H_lev, + lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], + geom[lev], use_moist); + + // Write this time slice to erfbdy immediately after conversion + if (write_erfbdy) { + WriteERFBdyTimeSlice(erfbdy_file, itime, + bdy_data_xlo[itime], bdy_data_xhi[itime], + bdy_data_ylo[itime], bdy_data_yhi[itime], + WRFBdyVars::NumTypes); + Print() << "Wrote erfbdy time index " << itime << " of " << ntimes_total-1 << std::endl; + + // Note: We don't clear time slices here since we only load first 3 + } + } // itime + + // If writing erfbdy and we have more times, process them now + if (write_erfbdy && ntimes_total > 3) { + Print() << "Processing remaining " << ntimes_total - 3 << " boundary times..." << std::endl; + for (int itime = 3; itime < ntimes_total; ++itime) { + read_from_wrfbdy(itime, nc_bdy_file, geom[0].Domain(), + bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, + real_width); + + convert_all_wrfbdy_data(itime, geom[lev].Domain(), + bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, + mf_MUB_lev, mf_C1H_lev, mf_C2H_lev, + lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], + geom[lev], use_moist); + + WriteERFBdyTimeSlice(erfbdy_file, itime, + bdy_data_xlo[itime], bdy_data_xhi[itime], + bdy_data_ylo[itime], bdy_data_yhi[itime], + WRFBdyVars::NumTypes); + Print() << "Wrote erfbdy time index " << itime << " of " << ntimes_total-1 << std::endl; + + // Clear this time from memory after writing. + bdy_data_xlo[itime].clear(); + bdy_data_xhi[itime].clear(); + bdy_data_ylo[itime].clear(); + bdy_data_yhi[itime].clear(); + } + Print() << "Completed writing erfbdy times" << std::endl; + } + } // // Start at the earliest time (read_from_wrfbdy) diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 891f03dc4b..234f052f5a 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -1,6 +1,7 @@ #include #include #include +#include using namespace amrex; @@ -38,40 +39,66 @@ ERF::timeStep (int lev, Real time, int /*iteration*/) int n_time_old = std::min(static_cast( (time_since_start_bdy ) / bdy_time_interval), ntimes-1); int n_time_new = std::min(static_cast( (time_since_start_bdy+dt[lev]) / bdy_time_interval), ntimes-1); - for (int itime = 0; itime < ntimes; itime++) - { - /* - if (bdy_data_xlo[itime].size() > 0) { - amrex::Print() << "HAVE BDY DATA AT TIME " << itime << std::endl; - } else { - amrex::Print() << " NO BDY DATA AT TIME " << itime << std::endl; - } - */ + // Handle erfbdy files (AMReX native format) + if (use_erfbdy) { + for (int itime = 0; itime < ntimes; itime++) + { + bool clear_itime = (itime < n_time_old); - bool clear_itime = (itime < n_time_old); + if (clear_itime && bdy_data_xlo[itime].size() > 0) { + bdy_data_xlo[itime].clear(); + bdy_data_xhi[itime].clear(); + bdy_data_ylo[itime].clear(); + bdy_data_yhi[itime].clear(); + } - if (clear_itime && bdy_data_xlo[itime].size() > 0) { - bdy_data_xlo[itime].clear(); - bdy_data_xhi[itime].clear(); - bdy_data_ylo[itime].clear(); - bdy_data_yhi[itime].clear(); - //amrex::Print() << "CLEAR BDY DATA AT TIME " << itime << std::endl; - } + bool need_itime = (itime >= n_time_old && itime <= n_time_new+1); - bool need_itime = (itime >= n_time_old && itime <= n_time_new+1); - //if (need_itime) { amrex::Print() << "NEED BDY DATA AT TIME " << itime << std::endl; } - - if (bdy_data_xlo[itime].size() == 0 && need_itime) { - read_from_wrfbdy(itime,nc_bdy_file,geom[0].Domain(), - bdy_data_xlo,bdy_data_xhi,bdy_data_ylo,bdy_data_yhi, - real_width); - - convert_all_wrfbdy_data(itime, geom[0].Domain(), bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, - *mf_MUB, *mf_C1H, *mf_C2H, - vars_new[lev][Vars::xvel], vars_new[lev][Vars::yvel], vars_new[lev][Vars::cons], - geom[lev], use_moist); - } - } // itime + if (bdy_data_xlo[itime].size() == 0 && need_itime) { + read_from_erfbdy(itime, erfbdy_file, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + nvars_erfbdy, real_width); + } + } // itime + } + // Handle wrfbdy files (NetCDF format) + else { + for (int itime = 0; itime < ntimes; itime++) + { + /* + if (bdy_data_xlo[itime].size() > 0) { + amrex::Print() << "HAVE BDY DATA AT TIME " << itime << std::endl; + } else { + amrex::Print() << " NO BDY DATA AT TIME " << itime << std::endl; + } + */ + + bool clear_itime = (itime < n_time_old); + + if (clear_itime && bdy_data_xlo[itime].size() > 0) { + bdy_data_xlo[itime].clear(); + bdy_data_xhi[itime].clear(); + bdy_data_ylo[itime].clear(); + bdy_data_yhi[itime].clear(); + //amrex::Print() << "CLEAR BDY DATA AT TIME " << itime << std::endl; + } + + bool need_itime = (itime >= n_time_old && itime <= n_time_new+1); + //if (need_itime) { amrex::Print() << "NEED BDY DATA AT TIME " << itime << std::endl; } + + if (bdy_data_xlo[itime].size() == 0 && need_itime) { + read_from_wrfbdy(itime,nc_bdy_file,geom[0].Domain(), + bdy_data_xlo,bdy_data_xhi,bdy_data_ylo,bdy_data_yhi, + real_width); + + convert_all_wrfbdy_data(itime, geom[0].Domain(), bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, + *mf_MUB, *mf_C1H, *mf_C2H, + vars_new[lev][Vars::xvel], vars_new[lev][Vars::yvel], vars_new[lev][Vars::cons], + geom[lev], use_moist); + } + } // itime + } } // use_real_bcs && lev == 0 if (!nc_low_file.empty() && (lev==0)) From 462bd36a40314874c29f6a1c6b28fa62178adf12 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Tue, 28 Apr 2026 21:51:46 -0700 Subject: [PATCH 02/27] Resolved parallel issues with boundary file read. --- Source/IO/ERF_ReadFromERFBdy.cpp | 83 +++++++++++++------------------- 1 file changed, 33 insertions(+), 50 deletions(-) diff --git a/Source/IO/ERF_ReadFromERFBdy.cpp b/Source/IO/ERF_ReadFromERFBdy.cpp index 7c53ce3e85..95045f458e 100644 --- a/Source/IO/ERF_ReadFromERFBdy.cpp +++ b/Source/IO/ERF_ReadFromERFBdy.cpp @@ -78,60 +78,43 @@ read_from_erfbdy(int itime, bdy_data_yhi[itime].resize(nvars); } - // Read each variable for each boundary direction using FArrayBox::readFrom - if (ParallelDescriptor::IOProcessor()) + // All processors read the data from files + // Since erfbdy files are in a parallel filesystem, all processors can read simultaneously + for (int ivar = 0; ivar < nvars; ++ivar) { - for (int ivar = 0; ivar < nvars; ++ivar) + // X-low boundary + { + std::string filename = time_dir + "/BdyData_xlo_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_xlo[itime][ivar].readFrom(ifs); + ifs.close(); + } + + // X-high boundary + { + std::string filename = time_dir + "/BdyData_xhi_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_xhi[itime][ivar].readFrom(ifs); + ifs.close(); + } + + // Y-low boundary { - // X-low boundary - { - std::string filename = time_dir + "/BdyData_xlo_var" + std::to_string(ivar); - std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_xlo[itime][ivar].readFrom(ifs); - ifs.close(); - } - - // X-high boundary - { - std::string filename = time_dir + "/BdyData_xhi_var" + std::to_string(ivar); - std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_xhi[itime][ivar].readFrom(ifs); - ifs.close(); - } - - // Y-low boundary - { - std::string filename = time_dir + "/BdyData_ylo_var" + std::to_string(ivar); - std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_ylo[itime][ivar].readFrom(ifs); - ifs.close(); - } - - // Y-high boundary - { - std::string filename = time_dir + "/BdyData_yhi_var" + std::to_string(ivar); - std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_yhi[itime][ivar].readFrom(ifs); - ifs.close(); - } + std::string filename = time_dir + "/BdyData_ylo_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_ylo[itime][ivar].readFrom(ifs); + ifs.close(); + } + + // Y-high boundary + { + std::string filename = time_dir + "/BdyData_yhi_var" + std::to_string(ivar); + std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); + bdy_data_yhi[itime][ivar].readFrom(ifs); + ifs.close(); } } - // Broadcast data to all processors + // Barrier to ensure all reads complete ParallelDescriptor::Barrier(); - for (int ivar = 0; ivar < nvars; ++ivar) - { - ParallelDescriptor::Bcast(bdy_data_xlo[itime][ivar].dataPtr(), - bdy_data_xlo[itime][ivar].size(), - ParallelDescriptor::IOProcessorNumber()); - ParallelDescriptor::Bcast(bdy_data_xhi[itime][ivar].dataPtr(), - bdy_data_xhi[itime][ivar].size(), - ParallelDescriptor::IOProcessorNumber()); - ParallelDescriptor::Bcast(bdy_data_ylo[itime][ivar].dataPtr(), - bdy_data_ylo[itime][ivar].size(), - ParallelDescriptor::IOProcessorNumber()); - ParallelDescriptor::Bcast(bdy_data_yhi[itime][ivar].dataPtr(), - bdy_data_yhi[itime][ivar].size(), - ParallelDescriptor::IOProcessorNumber()); - } } From 394b24d8369a7a5b3db4dc5cea3db92fdbe678a4 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 30 Apr 2026 11:30:49 -0700 Subject: [PATCH 03/27] Update Docs to include erfbdy description and paramters. --- Docs/sphinx_doc/BoundaryConditions.rst | 31 +++++--- Docs/sphinx_doc/Inputs.rst | 13 ++++ Exec/GNUmakefile | 3 +- Source/ERF.H | 10 +-- Source/ERF.cpp | 20 ++++- Source/IO/ERF_Checkpoint.cpp | 74 +++++++++++-------- Source/Initialization/ERF_InitFromMetgrid.cpp | 8 +- .../Initialization/ERF_InitFromWRFInput.cpp | 4 + 8 files changed, 113 insertions(+), 50 deletions(-) diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst index 779758ee50..e82dcf0218 100644 --- a/Docs/sphinx_doc/BoundaryConditions.rst +++ b/Docs/sphinx_doc/BoundaryConditions.rst @@ -193,16 +193,17 @@ state data. Real Domain BCs ---------------------- -When using real lateral boundary conditions, time-dependent observation data is read -from a file. In ERF, the observation data is utilized to nudge the solution state towards -the observation data in the interior of the domain. We explicitly note that the wall normal -velocity on a domain boundary is **not** nudged toward the observation data but is computed -from a zero normal gradient (outflow) condition that uses a one-sided difference on the domain -interior. Therefore, a key difference between real domain BCs in ERF and WRF is that ERF allows -the normal velocity at a domain wall to adjust, given the nudged adjacent velocities, while -WRF imposes a Dirichlet condition; see the yellow region in Figure 8 below. -The user may specify (in the inputs file) the total width of the interior boundary region with -``erf.real_width = `` (yellow + blue). The real BCs are only imposed for +When using real lateral boundary conditions, ``init_type`` of ``metgrid`` or ``wrfinput``, +time-dependent observation data is read from one or more files. In ERF, the observation data +is utilized to nudge the solution state towards the observation data in the interior of the +domain. We explicitly note that the wall normal velocity on a domain boundary is **not** +nudged toward the observation data but is computed from a zero normal gradient (outflow) +condition that uses a one-sided difference on the domain interior. Therefore, a key +difference between real domain BCs in ERF and WRF is that ERF allows the normal velocity at +a domain wall to adjust, given the nudged adjacent velocities, while WRF imposes a Dirichlet +condition; see the yellow region in Figure 8 below. The user may specify (in the inputs +file) the total width of the interior boundary region with ``erf.real_width = `` (yellow ++ blue). The real BCs are only imposed for :math:`\psi = \left\{ \theta; \; q_v; \; u; \; v \right\}`. .. |wrfbdy| image:: figures/wrfbdy_BCs.png @@ -237,3 +238,13 @@ that is 1 at the domain boundary and 0 at the edge of the nudge region. For turbulent inflow perturbations, including the ``erf.perturbation_type`` options and their placement on level subdomains, see :doc:`InflowTurbulenceGeneration`. + +ERF Boundary Data Files +---------------------- + +For real lateral boundary conditions, ``init_type`` of ``metgrid`` or ``wrfinput``, ERF can +write lateral boundary data to an AMReX-native file during initialization. The file name is +configurable with ``erf.erfbdy_file`` (default ``"erfbdy"``). Use of the the erfbdy file is +mandatory for ``init_type`` of ``metgrid``, but can be optionally disabled for ``init_type`` +of ``wrfinput`` with ``erf.write_erfbdy = false`` so that the boundary data is processed +as-needed during time integration. diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 618f715c54..096c11be13 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -625,6 +625,19 @@ PlotFiles See :ref:`sec:Plotfiles` for how to control the types and frequency of plotfile generation. +Boundary Files +============== + ++----------------------+------------------------------+-------------------+----------------------+ +| Parameter | Definition | Acceptable Values | Default | ++======================+==============================+===================+======================+ +| **erf.write_erfbdy** | Write AMReX-native format | true / false | true for non-restart | +| | boundary file for real-data | | real data cases, | +| | cases only | | otherwise false | ++----------------------+------------------------------+-------------------+----------------------+ +| **erf.erfbdy_file** | Name of the boundary file | String | "erfbdy" | ++----------------------+------------------------------+-------------------+----------------------+ + Additional Plotfile Controls ---------------------------- diff --git a/Exec/GNUmakefile b/Exec/GNUmakefile index 100bf8d3c8..098a539412 100644 --- a/Exec/GNUmakefile +++ b/Exec/GNUmakefile @@ -1,7 +1,6 @@ # AMReX COMP = gnu -#PRECISION = DOUBLE -PRECISION = SINGLE +PRECISION = DOUBLE # Profiling PROFILE = FALSE diff --git a/Source/ERF.H b/Source/ERF.H index 84b4e6fa51..73c486b440 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -1269,11 +1269,11 @@ private: int metgrid_order{2}; int metgrid_force_sfc_k{6}; - // Options for ERF boundary files (wrfinput/metgrid) - bool write_erfbdy{false}; - bool use_erfbdy{false}; - std::string erfbdy_file{"erfbdy"}; - int nvars_erfbdy{0}; + // Options for ERF boundary files for real simulations. + bool write_erfbdy{false}; // User-settable: write erfbdy during initialization + bool use_erfbdy{false}; // Internal: set during init or restart + std::string erfbdy_file{"erfbdy"}; // User-settable: erfbdy file name + int nvars_erfbdy{0}; // Internal: number of variables in erfbdy amrex::Vector ba1d; amrex::Vector ba2d; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 45ace7dca1..cd1e1a8b91 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -52,7 +52,6 @@ Real ERF::final_low_time = -one; Real ERF::bdy_time_interval = std::numeric_limits::max(); Real ERF::low_time_interval = std::numeric_limits::max(); - #endif // Time step control @@ -2503,9 +2502,26 @@ ERF::ReadParameters () // Options for boundary file. pp.query("write_erfbdy", write_erfbdy); - pp.query("use_erfbdy", use_erfbdy); pp.query("erfbdy_file", erfbdy_file); + // Set a default value for write_erfbdy. + // write_erfbdy should be falst for restarts. + // write_erfbdy should be true for clean starts of the metgrid or wrfinput pathways. + bool is_restart = !restart_chkfile.empty(); + if (is_restart) { + // Restart run: write_erfbdy must be false + if (write_erfbdy) { + Abort("Cannot set erf.write_erfbdy = true during restart. erfbdy should only be written during initial runs."); + } + } else { + // Initial run: default write_erfbdy to true for metgrid/wrfinput + if (!pp.contains("write_erfbdy")) { + if ((solverChoice.init_type == InitType::Metgrid) || (solverChoice.init_type == InitType::WRFInput)) { + write_erfbdy = true; + } + } + } + // Set default to FullState for now ... later we will try Perturbation interpolation_type = StateInterpType::FullState; pp.query_enum_case_insensitive("interpolation_type" ,interpolation_type); diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index 762ba64943..3fe9ec2ba1 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -1069,37 +1069,51 @@ ERF::ReadCheckpointFile () #endif #endif - // Load boundary data from erfbdy. + // Load boundary data from erfbdy during restart for metgrid or wrfinput. if (((solverChoice.init_type == InitType::WRFInput) || (solverChoice.init_type == InitType::Metgrid)) && - solverChoice.use_real_bcs && use_erfbdy) { - Print() << "Restart: Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; - - // Read metadata and times from erfbdy. - int ntimes_erfbdy; - Vector bdy_times; - bdy_time_interval = read_times_from_erfbdy(erfbdy_file, - ntimes_erfbdy, nvars_erfbdy, real_width, - bdy_times, start_bdy_time, final_bdy_time); - - Print() << "Restart: erfbdy file contains " << ntimes_erfbdy << " times" << std::endl; - - bdy_data_xlo.resize(ntimes_erfbdy); - bdy_data_xhi.resize(ntimes_erfbdy); - bdy_data_ylo.resize(ntimes_erfbdy); - bdy_data_yhi.resize(ntimes_erfbdy); - - // Determine which times we need based on current simulation time. - Real time_since_start_bdy = t_new[0] + start_time - start_bdy_time; - int n_time_old = std::min(static_cast(time_since_start_bdy / bdy_time_interval), ntimes_erfbdy-1); - int n_time_new = n_time_old + 1; - - // Load the necessary time slices. - for (int itime = n_time_old; itime <= std::min(n_time_new + 1, ntimes_erfbdy - 1); ++itime) { - read_from_erfbdy(itime, erfbdy_file, - bdy_data_xlo, bdy_data_xhi, - bdy_data_ylo, bdy_data_yhi, - nvars_erfbdy, real_width); - Print() << "Restart: Loaded erfbdy time index " << itime << std::endl; + solverChoice.use_real_bcs) { + + // Check for erfbdy file. + std::string erfbdy_header = erfbdy_file + "/Header"; + use_erfbdy = FileSystem::Exists(erfbdy_header); + + // Handle metgrid case: erfbdy is required + if (solverChoice.init_type == InitType::Metgrid) { + if (!use_erfbdy) { + Abort("Restart with init_type=metgrid requires erfbdy file: " + erfbdy_file); + } + } + + // Load from erfbdy if it exists. + if (use_erfbdy) { + Print() << "Restart: Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; + + int ntimes_erfbdy; + Vector bdy_times; + bdy_time_interval = read_times_from_erfbdy(erfbdy_file, + ntimes_erfbdy, nvars_erfbdy, real_width, + bdy_times, start_bdy_time, final_bdy_time); + + Print() << "Restart: erfbdy file contains " << ntimes_erfbdy << " times" << std::endl; + + bdy_data_xlo.resize(ntimes_erfbdy); + bdy_data_xhi.resize(ntimes_erfbdy); + bdy_data_ylo.resize(ntimes_erfbdy); + bdy_data_yhi.resize(ntimes_erfbdy); + + // Determine which times we need based on current simulation time. + Real time_since_start_bdy = t_new[0] + start_time - start_bdy_time; + int n_time_old = std::min(static_cast(time_since_start_bdy / bdy_time_interval), ntimes_erfbdy-1); + int n_time_new = n_time_old + 1; + + // Read the necessary times into memory. + for (int itime = n_time_old; itime <= std::min(n_time_new + 1, ntimes_erfbdy - 1); ++itime) { + read_from_erfbdy(itime, erfbdy_file, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + nvars_erfbdy, real_width); + Print() << "Restart: Loaded erfbdy time index " << itime << std::endl; + } } } } diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 5373085c18..98399597e2 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -55,7 +55,13 @@ ERF::init_from_metgrid (int lev) Print() << "Init with met_em without moisture model." << std::endl; } - // If using erfbdy file, load boundary data and then skip the rest of met_em processing. + // Check for erfbdy file. + if (lev == 0) { + std::string erfbdy_header = erfbdy_file + "/Header"; + use_erfbdy = FileSystem::Exists(erfbdy_header); + } + + // If erfbdy file exists, load boundary data and skip met_em processing. if (lev == 0 && use_erfbdy) { Print() << "Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index 796b0b44e5..a42fce3db3 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -1036,6 +1036,10 @@ ERF::init_from_wrfinput (int lev, amrex::Error("Cannot set periodic lateral boundary conditions when reading in real boundary values"); } + // Check for erfbdy file. + std::string erfbdy_header = erfbdy_file + "/Header"; + use_erfbdy = FileSystem::Exists(erfbdy_header); + // Path 1: Load from existing erfbdy file if (use_erfbdy) { Print() << "Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; From 74744af1a137c989e2c113edeb1f0b3b791005f7 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Mon, 4 May 2026 13:09:18 -0700 Subject: [PATCH 04/27] Cleaned up comments and added assertion for catching empty erfbdy files. --- Source/ERF.cpp | 9 ++--- Source/IO/ERF_Checkpoint.cpp | 1 - Source/IO/ERF_ReadFromERFBdy.H | 6 +-- Source/IO/ERF_ReadFromERFBdy.cpp | 37 ++++++++----------- Source/IO/ERF_WriteERFBdy.H | 4 +- Source/IO/ERF_WriteERFBdy.cpp | 32 ++++++++-------- Source/Initialization/ERF_InitFromMetgrid.cpp | 9 ++--- .../Initialization/ERF_InitFromWRFInput.cpp | 20 ++++------ Source/TimeIntegration/ERF_TimeStep.cpp | 4 +- 9 files changed, 54 insertions(+), 68 deletions(-) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index c786a499eb..18e8748e55 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -2506,17 +2506,16 @@ ERF::ReadParameters () pp.query("write_erfbdy", write_erfbdy); pp.query("erfbdy_file", erfbdy_file); - // Set a default value for write_erfbdy. - // write_erfbdy should be falst for restarts. - // write_erfbdy should be true for clean starts of the metgrid or wrfinput pathways. + // Set a default value for write_erfbdy following these rules. + // Prioritize write_erfbdy provided by user. + // write_erfbdy must be false for restarts. + // write_erfbdy defaults to true for clean starts of the metgrid or wrfinput pathways. bool is_restart = !restart_chkfile.empty(); if (is_restart) { - // Restart run: write_erfbdy must be false if (write_erfbdy) { Abort("Cannot set erf.write_erfbdy = true during restart. erfbdy should only be written during initial runs."); } } else { - // Initial run: default write_erfbdy to true for metgrid/wrfinput if (!pp.contains("write_erfbdy")) { if ((solverChoice.init_type == InitType::Metgrid) || (solverChoice.init_type == InitType::WRFInput)) { write_erfbdy = true; diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index 3fe9ec2ba1..0a1ff97c29 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -1077,7 +1077,6 @@ ERF::ReadCheckpointFile () std::string erfbdy_header = erfbdy_file + "/Header"; use_erfbdy = FileSystem::Exists(erfbdy_header); - // Handle metgrid case: erfbdy is required if (solverChoice.init_type == InitType::Metgrid) { if (!use_erfbdy) { Abort("Restart with init_type=metgrid requires erfbdy file: " + erfbdy_file); diff --git a/Source/IO/ERF_ReadFromERFBdy.H b/Source/IO/ERF_ReadFromERFBdy.H index 2923d9c42a..834fa69395 100644 --- a/Source/IO/ERF_ReadFromERFBdy.H +++ b/Source/IO/ERF_ReadFromERFBdy.H @@ -4,8 +4,8 @@ #include #include -// Read metadata and time values from erfbdy Header -// Returns bdy_time_interval +// Read metadata and time values from erfbdy Header. +// Returns bdy_time_interval. amrex::Real read_times_from_erfbdy(const std::string& bdy_file_name, int& ntimes, @@ -15,7 +15,7 @@ read_times_from_erfbdy(const std::string& bdy_file_name, amrex::Real& start_bdy_time, amrex::Real& final_bdy_time); -// Read a single time slice from erfbdy file +// Read a single time index from erfbdy file. void read_from_erfbdy(int itime, const std::string& bdy_file_name, diff --git a/Source/IO/ERF_ReadFromERFBdy.cpp b/Source/IO/ERF_ReadFromERFBdy.cpp index 95045f458e..cc46365556 100644 --- a/Source/IO/ERF_ReadFromERFBdy.cpp +++ b/Source/IO/ERF_ReadFromERFBdy.cpp @@ -17,7 +17,7 @@ read_times_from_erfbdy(const std::string& bdy_file_name, { std::string HeaderFileName = bdy_file_name + "/Header"; - // Read header file + // Read header file. std::ifstream HeaderFile; HeaderFile.open(HeaderFileName.c_str(), std::ifstream::in); @@ -27,33 +27,36 @@ read_times_from_erfbdy(const std::string& bdy_file_name, std::string line; - // Read title line + // Read title line. std::getline(HeaderFile, line); - // Read metadata + // Read metadata. HeaderFile >> ntimes; HeaderFile >> nvars; HeaderFile >> real_width; - // Read time values + // Read time values. bdy_times.resize(ntimes); for (int i = 0; i < ntimes; ++i) { HeaderFile >> bdy_times[i]; } - // Read domain box (stored but not used in this function) + // Read domain box (stored but not used here). int sml[3], big[3]; HeaderFile >> sml[0] >> sml[1] >> sml[2]; HeaderFile >> big[0] >> big[1] >> big[2]; HeaderFile.close(); - // Set start and final times + // Ensure the file hodls at least two times. + AMREX_ALWAYS_ASSERT(ntimes >= 2); + + // Set start and final times. start_bdy_time = bdy_times[0]; final_bdy_time = bdy_times[ntimes - 1]; - // Calculate time interval - Real bdy_time_interval = (ntimes > 1) ? (bdy_times[1] - bdy_times[0]) : 0.0; + // Calculate time interval. + Real bdy_time_interval = bdy_times[1] - bdy_times[0]; return bdy_time_interval; } @@ -67,10 +70,8 @@ read_from_erfbdy(int itime, Vector>& bdy_data_yhi, int nvars, int real_width) { - // Construct time subdirectory path std::string time_dir = bdy_file_name + "/Time_" + Concatenate("", itime, 6); - // Ensure vectors are sized correctly if (bdy_data_xlo[itime].empty()) { bdy_data_xlo[itime].resize(nvars); bdy_data_xhi[itime].resize(nvars); @@ -78,36 +79,30 @@ read_from_erfbdy(int itime, bdy_data_yhi[itime].resize(nvars); } - // All processors read the data from files - // Since erfbdy files are in a parallel filesystem, all processors can read simultaneously for (int ivar = 0; ivar < nvars; ++ivar) { - // X-low boundary - { + { // X-low boundary std::string filename = time_dir + "/BdyData_xlo_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); bdy_data_xlo[itime][ivar].readFrom(ifs); ifs.close(); } - // X-high boundary - { + { // X-high boundary std::string filename = time_dir + "/BdyData_xhi_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); bdy_data_xhi[itime][ivar].readFrom(ifs); ifs.close(); } - // Y-low boundary - { + { // Y-low boundary std::string filename = time_dir + "/BdyData_ylo_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); bdy_data_ylo[itime][ivar].readFrom(ifs); ifs.close(); } - // Y-high boundary - { + { // Y-high boundary std::string filename = time_dir + "/BdyData_yhi_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); bdy_data_yhi[itime][ivar].readFrom(ifs); @@ -115,6 +110,6 @@ read_from_erfbdy(int itime, } } - // Barrier to ensure all reads complete + // Barrier to ensure all reads complete. ParallelDescriptor::Barrier(); } diff --git a/Source/IO/ERF_WriteERFBdy.H b/Source/IO/ERF_WriteERFBdy.H index d35d6b30b2..05b820a4a4 100644 --- a/Source/IO/ERF_WriteERFBdy.H +++ b/Source/IO/ERF_WriteERFBdy.H @@ -5,7 +5,7 @@ #include #include -// Initialize erfbdy file and write header (call once at start) +// Initialize erfbdy file and write header. void InitERFBdyFile(const std::string& bdy_file_name, int ntimes, const amrex::Vector& bdy_times, @@ -13,7 +13,7 @@ void InitERFBdyFile(const std::string& bdy_file_name, int nvars, int real_width); -// Write a single time slice (call once per time) +// Write a single time index. void WriteERFBdyTimeSlice(const std::string& bdy_file_name, int itime, const amrex::Vector& bdy_data_xlo, diff --git a/Source/IO/ERF_WriteERFBdy.cpp b/Source/IO/ERF_WriteERFBdy.cpp index 01e19ddccf..0f632ffdf4 100644 --- a/Source/IO/ERF_WriteERFBdy.cpp +++ b/Source/IO/ERF_WriteERFBdy.cpp @@ -13,15 +13,15 @@ void InitERFBdyFile(const std::string& bdy_file_name, int nvars, int real_width) { - // Only I/O processor creates the directory and writes header + // Only I/O processor creates the directory and writes header. if (ParallelDescriptor::IOProcessor()) { - // Create erfbdy directory + // Create erfbdy directory. if (!amrex::UtilCreateDirectory(bdy_file_name, 0755)) { amrex::CreateDirectoryFailed(bdy_file_name); } - // Write Header file + // Write Header file. std::string HeaderFileName = bdy_file_name + "/Header"; std::ofstream HeaderFile; HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out | std::ofstream::trunc); @@ -32,15 +32,15 @@ void InitERFBdyFile(const std::string& bdy_file_name, HeaderFile.precision(17); - // Title + // Title. HeaderFile << "ERFBdy file for ERF\n"; - // Metadata + // Metadata. HeaderFile << ntimes << "\n"; HeaderFile << nvars << "\n"; HeaderFile << real_width << "\n"; - // Time values + // Time values. for (int i = 0; i < ntimes; ++i) { HeaderFile << bdy_times[i]; if (i < ntimes - 1) { @@ -49,14 +49,14 @@ void InitERFBdyFile(const std::string& bdy_file_name, } HeaderFile << "\n"; - // Domain box + // Domain box. HeaderFile << domain.smallEnd(0) << " " << domain.smallEnd(1) << " " << domain.smallEnd(2) << "\n"; HeaderFile << domain.bigEnd(0) << " " << domain.bigEnd(1) << " " << domain.bigEnd(2) << "\n"; HeaderFile.close(); } - // Barrier to ensure directory is created before any writes + // Barrier to ensure directory is created before any writes. ParallelDescriptor::Barrier(); } @@ -68,7 +68,7 @@ void WriteERFBdyTimeSlice(const std::string& bdy_file_name, const Vector& bdy_data_yhi, int nvars) { - // Create time subdirectory + // Create time subdirectory. std::string time_dir = bdy_file_name + "/Time_" + Concatenate("", itime, 6); if (ParallelDescriptor::IOProcessor()) @@ -78,15 +78,15 @@ void WriteERFBdyTimeSlice(const std::string& bdy_file_name, } } - // Barrier to ensure directory exists + // Barrier to ensure directory exists. ParallelDescriptor::Barrier(); - // Write each variable for each boundary direction using FArrayBox::writeOn + // Write each variable for each boundary direction using FArrayBox::writeOn. if (ParallelDescriptor::IOProcessor()) { for (int ivar = 0; ivar < nvars; ++ivar) { - // X-low boundary + // X-low boundary. if (bdy_data_xlo.size() > 0 && ivar < bdy_data_xlo.size()) { std::string filename = time_dir + "/BdyData_xlo_var" + std::to_string(ivar); std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); @@ -94,7 +94,7 @@ void WriteERFBdyTimeSlice(const std::string& bdy_file_name, ofs.close(); } - // X-high boundary + // X-high boundary. if (bdy_data_xhi.size() > 0 && ivar < bdy_data_xhi.size()) { std::string filename = time_dir + "/BdyData_xhi_var" + std::to_string(ivar); std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); @@ -102,7 +102,7 @@ void WriteERFBdyTimeSlice(const std::string& bdy_file_name, ofs.close(); } - // Y-low boundary + // Y-low boundary. if (bdy_data_ylo.size() > 0 && ivar < bdy_data_ylo.size()) { std::string filename = time_dir + "/BdyData_ylo_var" + std::to_string(ivar); std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); @@ -110,7 +110,7 @@ void WriteERFBdyTimeSlice(const std::string& bdy_file_name, ofs.close(); } - // Y-high boundary + // Y-high boundary. if (bdy_data_yhi.size() > 0 && ivar < bdy_data_yhi.size()) { std::string filename = time_dir + "/BdyData_yhi_var" + std::to_string(ivar); std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary); @@ -120,6 +120,6 @@ void WriteERFBdyTimeSlice(const std::string& bdy_file_name, } } - // Barrier to ensure all writes complete + // Barrier to ensure all writes complete. ParallelDescriptor::Barrier(); } diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 98399597e2..b0f56c37b3 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -288,7 +288,6 @@ ERF::init_from_metgrid (int lev) NC_lmask_iab, geom[lev]); if (lev == 0) { - // Collect time for erfbdy bdy_times[itime] = NC_epochTime[itime]; if (itime == 0) { @@ -303,7 +302,7 @@ ERF::init_from_metgrid (int lev) t_new[lev] = zero; t_old[lev] = -Real(1.e200); - // Initialize erfbdy file header now that we know domain and real_width + // Initialize erfbdy file header now that we have the domain and real_width. if (write_erfbdy) { InitERFBdyFile(erfbdy_file, ntimes, bdy_times, geom[lev].Domain(), MetGridBdyVars::NumTypes, real_width); @@ -613,7 +612,6 @@ ERF::init_from_metgrid (int lev) ParallelContext::CommunicatorAll()); } // nvar - // Write this time slice to erfbdy immediately to avoid memory accumulation if (write_erfbdy) { WriteERFBdyTimeSlice(erfbdy_file, itime, bdy_data_xlo[itime], bdy_data_xhi[itime], @@ -621,8 +619,9 @@ ERF::init_from_metgrid (int lev) MetGridBdyVars::NumTypes); Print() << "Wrote erfbdy time slice " << itime << " of " << ntimes-1 << std::endl; - // CRITICAL: Clear this time slice from memory after writing - // (unless it's one of the first two times needed for simulation start) + // Clear this time from memory after writing unless + // it's one of the first two times, which are needed + // at initialization. if (itime > 1) { bdy_data_xlo[itime].clear(); bdy_data_xhi[itime].clear(); diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index a42fce3db3..b23cccf521 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -1040,11 +1040,11 @@ ERF::init_from_wrfinput (int lev, std::string erfbdy_header = erfbdy_file + "/Header"; use_erfbdy = FileSystem::Exists(erfbdy_header); - // Path 1: Load from existing erfbdy file + // Path 1: Load from existing erfbdy file. if (use_erfbdy) { Print() << "Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; - // Read metadata and times from erfbdy + // Read metadata and times from erfbdy. int ntimes_erfbdy; Vector bdy_times; bdy_time_interval = read_times_from_erfbdy(erfbdy_file, @@ -1056,13 +1056,12 @@ ERF::init_from_wrfinput (int lev, Print() << "final_bdy_time = " << final_bdy_time << std::endl; Print() << "bdy_time_interval = " << bdy_time_interval << std::endl; - // Resize boundary data structures bdy_data_xlo.resize(ntimes_erfbdy); bdy_data_xhi.resize(ntimes_erfbdy); bdy_data_ylo.resize(ntimes_erfbdy); bdy_data_yhi.resize(ntimes_erfbdy); - // Load first 2 time slices for simulation start + // Load the first 2 times for simulation initialization. for (int itime = 0; itime < std::min(2, ntimes_erfbdy); ++itime) { read_from_erfbdy(itime, erfbdy_file, bdy_data_xlo, bdy_data_xhi, @@ -1074,7 +1073,7 @@ ERF::init_from_wrfinput (int lev, Print() << "Read in boundary data with width " << real_width << std::endl; Print() << "Running with relaxation width: " << real_width << std::endl; } - // Path 2: Load from wrfbdy and optionally write to erfbdy + // Path 2: Load from wrfbdy and optionally write to erfbdy. else { if (nc_bdy_file.empty()) { amrex::Error("NetCDF boundary file name must be provided via input"); @@ -1087,10 +1086,8 @@ ERF::init_from_wrfinput (int lev, int ntimes_total = bdy_data_xlo.size(); Vector bdy_times(ntimes_total); - // Initialize erfbdy file if writing is enabled + // Initialize erfbdy file. if (write_erfbdy) { - // Collect times (we'll need to read them from wrfbdy to fill this properly) - // For now, construct from start_bdy_time and bdy_time_interval for (int itime = 0; itime < ntimes_total; ++itime) { bdy_times[itime] = start_bdy_time + itime * bdy_time_interval; } @@ -1122,19 +1119,17 @@ ERF::init_from_wrfinput (int lev, lev_new[Vars::xvel], lev_new[Vars::yvel], lev_new[Vars::cons], geom[lev], use_moist); - // Write this time slice to erfbdy immediately after conversion + // Write this time to erfbdy. if (write_erfbdy) { WriteERFBdyTimeSlice(erfbdy_file, itime, bdy_data_xlo[itime], bdy_data_xhi[itime], bdy_data_ylo[itime], bdy_data_yhi[itime], WRFBdyVars::NumTypes); Print() << "Wrote erfbdy time index " << itime << " of " << ntimes_total-1 << std::endl; - - // Note: We don't clear time slices here since we only load first 3 } } // itime - // If writing erfbdy and we have more times, process them now + // If writing erfbdy and we have more than 3 times, then process the remaining times. if (write_erfbdy && ntimes_total > 3) { Print() << "Processing remaining " << ntimes_total - 3 << " boundary times..." << std::endl; for (int itime = 3; itime < ntimes_total; ++itime) { @@ -1154,7 +1149,6 @@ ERF::init_from_wrfinput (int lev, WRFBdyVars::NumTypes); Print() << "Wrote erfbdy time index " << itime << " of " << ntimes_total-1 << std::endl; - // Clear this time from memory after writing. bdy_data_xlo[itime].clear(); bdy_data_xhi[itime].clear(); bdy_data_ylo[itime].clear(); diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 234f052f5a..ed00ce5058 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -39,7 +39,7 @@ ERF::timeStep (int lev, Real time, int /*iteration*/) int n_time_old = std::min(static_cast( (time_since_start_bdy ) / bdy_time_interval), ntimes-1); int n_time_new = std::min(static_cast( (time_since_start_bdy+dt[lev]) / bdy_time_interval), ntimes-1); - // Handle erfbdy files (AMReX native format) + // Handle erfbdy files (AMReX native format). if (use_erfbdy) { for (int itime = 0; itime < ntimes; itime++) { @@ -62,7 +62,7 @@ ERF::timeStep (int lev, Real time, int /*iteration*/) } } // itime } - // Handle wrfbdy files (NetCDF format) + // Handle wrfbdy files (NetCDF format). else { for (int itime = 0; itime < ntimes; itime++) { From e22f5cc81cd80016eef9dbe26751a9cb4d3a7a9f Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Mon, 4 May 2026 21:14:40 -0700 Subject: [PATCH 05/27] Fixed timing issue with parsing write_erfbdy. Switched from times in metgrid being parsed from the global attr SIMULATION_START_DATE to the variable Times. --- Source/ERF.cpp | 37 ++++----- Source/IO/ERF_ReadFromMetgrid.cpp | 16 +++- Source/Initialization/ERF_InitFromMetgrid.cpp | 75 +++++++++++-------- 3 files changed, 79 insertions(+), 49 deletions(-) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 18e8748e55..8f6e0f6b6a 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -2506,23 +2506,6 @@ ERF::ReadParameters () pp.query("write_erfbdy", write_erfbdy); pp.query("erfbdy_file", erfbdy_file); - // Set a default value for write_erfbdy following these rules. - // Prioritize write_erfbdy provided by user. - // write_erfbdy must be false for restarts. - // write_erfbdy defaults to true for clean starts of the metgrid or wrfinput pathways. - bool is_restart = !restart_chkfile.empty(); - if (is_restart) { - if (write_erfbdy) { - Abort("Cannot set erf.write_erfbdy = true during restart. erfbdy should only be written during initial runs."); - } - } else { - if (!pp.contains("write_erfbdy")) { - if ((solverChoice.init_type == InitType::Metgrid) || (solverChoice.init_type == InitType::WRFInput)) { - write_erfbdy = true; - } - } - } - // Set default to FullState for now ... later we will try Perturbation interpolation_type = StateInterpType::FullState; pp.query_enum_case_insensitive("interpolation_type" ,interpolation_type); @@ -2734,6 +2717,26 @@ ERF::ReadParameters () solverChoice.init_params(max_level,pp_prefix); + // Set a default value for write_erfbdy following these rules. + // Prioritize write_erfbdy provided by user. + // write_erfbdy must be false for restarts. + // write_erfbdy defaults to true for clean starts of the metgrid or wrfinput pathways. + { + ParmParse pp(pp_prefix); + bool is_restart = !restart_chkfile.empty(); + if (is_restart) { + if (write_erfbdy) { + Abort("Cannot set erf.write_erfbdy = true during restart. erfbdy should only be written during initial runs."); + } + } else { + if (!pp.contains("write_erfbdy")) { + if ((solverChoice.init_type == InitType::Metgrid) || (solverChoice.init_type == InitType::WRFInput)) { + write_erfbdy = true; + } + } + } + } + { ParmParse pp_no_prefix; // Traditionally, max_step and stop_time do not have prefix. pp_no_prefix.query("max_step", max_step); diff --git a/Source/IO/ERF_ReadFromMetgrid.cpp b/Source/IO/ERF_ReadFromMetgrid.cpp index e9adcba71c..7e80685763 100644 --- a/Source/IO/ERF_ReadFromMetgrid.cpp +++ b/Source/IO/ERF_ReadFromMetgrid.cpp @@ -76,8 +76,20 @@ read_from_metgrid (int lev, int itime, ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0]; } - { // Global Attributes (string) - NC_dateTime = ncf.get_attr("SIMULATION_START_DATE"); + { // Read Times variable (character array). + auto times_var = ncf.var("Times"); + std::vector start = {0, 0}; // Time index 0, character index 0. + std::vector count = times_var.shape(); + count[0] = 1; // Read the first time entry. + + // Allocate buffer for the time string. + std::vector time_chars(count[1]); + times_var.get(time_chars.data(), start, count); + + // Convert to std::string and trim trailing whitespace/null characters. + NC_dateTime = std::string(time_chars.begin(), time_chars.end()); + NC_dateTime.erase(NC_dateTime.find_last_not_of(" \0") + 1); + const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); } diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index b0f56c37b3..38c163d95f 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -24,7 +24,19 @@ read_start_time_from_metgrid(int lev, const std::string& fname) if (ParallelDescriptor::IOProcessor()) { auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); - NC_dateTime = ncf.get_attr("SIMULATION_START_DATE"); + // Read the Times variable (character array). + auto times_var = ncf.var("Times"); + std::vector start = {0, 0}; // Time index 0, character index 0. + std::vector count = times_var.shape(); + count[0] = 1; // Read only the first time entry. + + // Allocate buffer for the time string. + std::vector time_chars(count[1]); + times_var.get(time_chars.data(), start, count); + + // Convert to std::string and trim trailing whitespace/null characters. + NC_dateTime = std::string(time_chars.begin(), time_chars.end()); + NC_dateTime.erase(NC_dateTime.find_last_not_of(" \0") + 1); const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); @@ -272,6 +284,19 @@ ERF::init_from_metgrid (int lev) Vector bdy_times(ntimes); + // Read times from met_em files (necessary for erfbdy header initialization). + if (lev == 0 && write_erfbdy) { + Print() << "Reading times from " << ntimes << " met_em files for erfbdy initialization" << std::endl; + for (int itime(0); itime < ntimes; itime++) { + bdy_times[itime] = read_start_time_from_metgrid(lev, nc_init_file[lev][itime]); + } + + // Initialize erfbdy file header. + InitERFBdyFile(erfbdy_file, ntimes, bdy_times, + geom[lev].Domain(), MetGridBdyVars::NumTypes, real_width); + Print() << "Initialized erfbdy file: " << erfbdy_file << std::endl; + } + for (int itime(0); itime < ntimes; itime++) { Print() << " init_from_metgrid: reading nc_init_file[" << lev << "][" << itime << "]\t" << nc_init_file[lev][itime] << std::endl; read_from_metgrid(lev, itime, @@ -301,13 +326,6 @@ ERF::init_from_metgrid (int lev) << " from metgrid file but note that time variable in simulation is elapsed time" << std::endl; t_new[lev] = zero; t_old[lev] = -Real(1.e200); - - // Initialize erfbdy file header now that we have the domain and real_width. - if (write_erfbdy) { - InitERFBdyFile(erfbdy_file, ntimes, bdy_times, - geom[lev].Domain(), MetGridBdyVars::NumTypes, real_width); - Print() << "Initialized erfbdy file: " << erfbdy_file << std::endl; - } } else { // Verify that files in nc_init_file[lev] are ordered from earliest to latest. AMREX_ALWAYS_ASSERT(NC_epochTime[itime] > NC_epochTime[itime-1]); @@ -580,23 +598,13 @@ ERF::init_from_metgrid (int lev) } // itime==0 - } // itime - - // FillBoundary to populate the internal halo cells - r_hse.FillBoundary(geom[lev].periodicity()); - p_hse.FillBoundary(geom[lev].periodicity()); - pi_hse.FillBoundary(geom[lev].periodicity()); - th_hse.FillBoundary(geom[lev].periodicity()); - qv_hse.FillBoundary(geom[lev].periodicity()); - - // 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. - if (lev == 0) { - for (int itime(0); itime < ntimes; itime++) { + if (lev == 0) { + // 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 nvar(0); nvar 1) { bdy_data_xlo[itime].clear(); bdy_data_xhi[itime].clear(); @@ -629,8 +636,16 @@ ERF::init_from_metgrid (int lev) bdy_data_yhi[itime].clear(); } } - } // itime - } // lev==0 + } // lev==0 + + } // itime + + // FillBoundary to populate the internal halo cells + r_hse.FillBoundary(geom[lev].periodicity()); + p_hse.FillBoundary(geom[lev].periodicity()); + pi_hse.FillBoundary(geom[lev].periodicity()); + th_hse.FillBoundary(geom[lev].periodicity()); + qv_hse.FillBoundary(geom[lev].periodicity()); Print() << "Running with relaxation width: " << real_width << std::endl; } From 14787440ecb8aab622b796891ec737ee537c8f18 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Tue, 5 May 2026 17:03:21 -0700 Subject: [PATCH 06/27] Bugfix for MAPFAC vars being 3D in wrfinput pathway. --- .../Initialization/ERF_InitFromWRFInput.cpp | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index b23cccf521..d47cfa8c23 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -273,11 +273,17 @@ ERF::init_from_wrfinput (int lev, // to be filled int nx = var_fab_from_file.box().length(0); int ny = var_fab_from_file.box().length(1); + int nz = var_fab_from_file.box().length(2); Box subdomain_tmp(subdomain_to_fill); if (nx == 1 and ny == 1) { subdomain_tmp.setBig(0,subdomain_tmp.smallEnd(0)); subdomain_tmp.setBig(1,subdomain_tmp.smallEnd(1)); } + // Handle 2D horizontal data (single layer in z) + // For 2D variables like MAPFAC_U, MAPFAC_V, PSFC, MUB, etc. + if (nz == 1) { + subdomain_tmp.setBig(2,subdomain_tmp.smallEnd(2)); + } Box subdomain_to_fill_typed(convert(subdomain_tmp,var_fab_from_file.box().ixType())); Box subdomain_crse(subdomain_to_fill_typed); @@ -731,6 +737,9 @@ ERF::init_from_wrfinput (int lev, // Initialize MapFac U if ( var_name == "MAPFAC_U" ) { Real max_val = var_fab.template max(); +// Box bx2d = var_fab.box(); +// bx2d.setRange(2, 0, 1); +// Real max_val = var_fab.template max(bx2d, 0); if (std::fabs(max_val) < std::numeric_limits::epsilon()) { Print() << "MAPFAC_U cannot be 0, resetting to 1!\n"; var_fab.template setVal(1); @@ -763,6 +772,9 @@ ERF::init_from_wrfinput (int lev, // Initialize MapFac V if ( var_name == "MAPFAC_V" ) { Real max_val = var_fab.template max(); +// Box bx2d = var_fab.box(); +// bx2d.setRange(2, 0, 1); +// Real max_val = var_fab.template max(bx2d, 0); if (std::fabs(max_val) < std::numeric_limits::epsilon()) { Print() << "MAPFAC_V cannot be 0, resetting to 1!\n"; var_fab.template setVal(1); @@ -795,6 +807,9 @@ ERF::init_from_wrfinput (int lev, // Initialize MapFac M if ( var_name == "MAPFAC_M" ) { Real max_val = var_fab.template max(); +// Box bx2d = var_fab.box(); +// bx2d.setRange(2, 0, 1); +// Real max_val = var_fab.template max(bx2d, 0); if (std::fabs(max_val) < std::numeric_limits::epsilon()) { Print() << "MAPFAC_M cannot be 0, resetting to 1!\n"; var_fab.template setVal(1); @@ -1048,8 +1063,8 @@ ERF::init_from_wrfinput (int lev, int ntimes_erfbdy; Vector bdy_times; bdy_time_interval = read_times_from_erfbdy(erfbdy_file, - ntimes_erfbdy, nvars_erfbdy, real_width, - bdy_times, start_bdy_time, final_bdy_time); + ntimes_erfbdy, nvars_erfbdy, real_width, + bdy_times, start_bdy_time, final_bdy_time); Print() << "erfbdy file contains " << ntimes_erfbdy << " time slices" << std::endl; Print() << "start_bdy_time = " << start_bdy_time << std::endl; From 01d7198ca6f834d0676713b796c364aedfe2828f Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Wed, 13 May 2026 13:34:14 -0700 Subject: [PATCH 07/27] Modified CMake/BuildERFExe.cmake to include IO/ERF_ReadFromERFBdy.cpp and IO/ERF_WriteERFBdy.cpp --- CMake/BuildERFExe.cmake | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 2d48b2e666..a29afa842c 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -352,6 +352,8 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/IO/ERF_WriteSubvolume.cpp ${SRC_DIR}/IO/ERF_WriteJobInfo.cpp ${SRC_DIR}/IO/ERF_ConsoleIO.cpp + ${SRC_DIR}/IO/ERF_ReadFromERFBdy.cpp + ${SRC_DIR}/IO/ERF_WriteERFBdy.cpp ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve.cpp ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve_tb.cpp ${SRC_DIR}/LinearSolvers/ERF_PoissonWallDist.cpp From c737a956d4a2817a58dda75fc5af78d546700e32 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 14 May 2026 14:47:43 -0700 Subject: [PATCH 08/27] Clean up of erfbdy processing. --- Source/IO/ERF_ReadFromERFBdy.cpp | 1 + Source/Initialization/ERF_InitFromMetgrid.cpp | 9 ++++++--- Source/Initialization/ERF_InitFromWRFInput.cpp | 9 +++++---- Source/TimeIntegration/ERF_TimeStep.cpp | 6 +++--- 4 files changed, 15 insertions(+), 10 deletions(-) diff --git a/Source/IO/ERF_ReadFromERFBdy.cpp b/Source/IO/ERF_ReadFromERFBdy.cpp index cc46365556..c2b930a8d9 100644 --- a/Source/IO/ERF_ReadFromERFBdy.cpp +++ b/Source/IO/ERF_ReadFromERFBdy.cpp @@ -70,6 +70,7 @@ read_from_erfbdy(int itime, Vector>& bdy_data_yhi, int nvars, int real_width) { + Print() << "Reading ERF boundary data for time index " << itime << std::endl; std::string time_dir = bdy_file_name + "/Time_" + Concatenate("", itime, 6); if (bdy_data_xlo[itime].empty()) { diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 38c163d95f..72fa105d91 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -72,6 +72,7 @@ ERF::init_from_metgrid (int lev) std::string erfbdy_header = erfbdy_file + "/Header"; use_erfbdy = FileSystem::Exists(erfbdy_header); } + if (use_erfbdy || write_erfbdy) nvars_erfbdy = MetGridBdyVars::NumTypes; // If erfbdy file exists, load boundary data and skip met_em processing. if (lev == 0 && use_erfbdy) { @@ -112,6 +113,8 @@ ERF::init_from_metgrid (int lev) return; // Skip the rest of met_em processing. } + use_erfbdy = true; + int ntimes = num_files_at_level[lev]; Print() << ntimes << " met_em.d0" << lev+1 << "*.nc files are listed" << std::endl; @@ -272,7 +275,7 @@ ERF::init_from_metgrid (int lev) bdy_data_xhi[itime][nvar].template setVal(0); bdy_data_ylo[itime][nvar].template setVal(0); bdy_data_yhi[itime][nvar].template setVal(0); - } + } // nvar } // itime } // lev==0 @@ -293,7 +296,7 @@ ERF::init_from_metgrid (int lev) // Initialize erfbdy file header. InitERFBdyFile(erfbdy_file, ntimes, bdy_times, - geom[lev].Domain(), MetGridBdyVars::NumTypes, real_width); + geom[lev].Domain(), nvars_erfbdy, real_width); Print() << "Initialized erfbdy file: " << erfbdy_file << std::endl; } @@ -624,7 +627,7 @@ ERF::init_from_metgrid (int lev) WriteERFBdyTimeSlice(erfbdy_file, itime, bdy_data_xlo[itime], bdy_data_xhi[itime], bdy_data_ylo[itime], bdy_data_yhi[itime], - MetGridBdyVars::NumTypes); + nvars_erfbdy); Print() << "Wrote erfbdy time slice " << itime << " of " << ntimes-1 << std::endl; // Clear this time from memory after writing unless it's one diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index e4a79875a6..4a1ed99f81 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -1044,6 +1044,7 @@ ERF::init_from_wrfinput (int lev, // Check for erfbdy file. std::string erfbdy_header = erfbdy_file + "/Header"; use_erfbdy = FileSystem::Exists(erfbdy_header); + if (use_erfbdy || write_erfbdy) nvars_erfbdy = WRFBdyVars::NumTypes; // Path 1: Load from existing erfbdy file. if (use_erfbdy) { @@ -1069,9 +1070,9 @@ ERF::init_from_wrfinput (int lev, // Load the first 2 times for simulation initialization. for (int itime = 0; itime < std::min(2, ntimes_erfbdy); ++itime) { read_from_erfbdy(itime, erfbdy_file, - bdy_data_xlo, bdy_data_xhi, - bdy_data_ylo, bdy_data_yhi, - nvars_erfbdy, real_width); + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + nvars_erfbdy, real_width); Print() << "Loaded erfbdy time slice " << itime << std::endl; } @@ -1098,7 +1099,7 @@ ERF::init_from_wrfinput (int lev, } InitERFBdyFile(erfbdy_file, ntimes_total, bdy_times, - geom[lev].Domain(), WRFBdyVars::NumTypes, real_width); + geom[lev].Domain(), nvars_erfbdy, real_width); Print() << "Initialized erfbdy file: " << erfbdy_file << std::endl; } diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index ed00ce5058..bb672fda86 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -56,9 +56,9 @@ ERF::timeStep (int lev, Real time, int /*iteration*/) if (bdy_data_xlo[itime].size() == 0 && need_itime) { read_from_erfbdy(itime, erfbdy_file, - bdy_data_xlo, bdy_data_xhi, - bdy_data_ylo, bdy_data_yhi, - nvars_erfbdy, real_width); + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + nvars_erfbdy, real_width); } } // itime } From b782acdab298c3105c0e7c8da9c06ce3c63587e3 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 13:01:55 -0700 Subject: [PATCH 09/27] compiled in debug needs testing. --- CMake/BuildERFExe.cmake | 1 + .../ERF_BoundaryConditionsRealbdy.cpp | 534 ++++++++++++++++++ .../ERF_FillIntermediatePatch.cpp | 13 +- Source/BoundaryConditions/ERF_FillPatch.cpp | 13 +- Source/BoundaryConditions/ERF_PhysBCFunct.cpp | 46 +- Source/BoundaryConditions/Make.package | 1 + Source/ERF.H | 20 + Source/ERF.cpp | 17 + Source/Utils/ERF_InteriorGhostCells.cpp | 184 ++++-- Source/Utils/ERF_Utils.H | 13 + 10 files changed, 776 insertions(+), 66 deletions(-) create mode 100644 Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 2d48b2e666..eb71af3b03 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -300,6 +300,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/BoundaryConditions/ERF_BoundaryConditionsZvel.cpp ${SRC_DIR}/BoundaryConditions/ERF_BoundaryConditionsBaseState.cpp ${SRC_DIR}/BoundaryConditions/ERF_BoundaryConditionsBndryReg.cpp + ${SRC_DIR}/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp ${SRC_DIR}/BoundaryConditions/ERF_FillPatch.cpp ${SRC_DIR}/BoundaryConditions/ERF_FillCoarsePatch.cpp ${SRC_DIR}/BoundaryConditions/ERF_FillIntermediatePatch.cpp diff --git a/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp b/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp new file mode 100644 index 0000000000..579f8967b1 --- /dev/null +++ b/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp @@ -0,0 +1,534 @@ +#include "ERF.H" +#include "ERF_Utils.H" + +using namespace amrex; + +#ifdef ERF_USE_NETCDF +/* + * Impose boundary conditions using data read in from wrfbdy 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_realbdy (const Vector& mfs, + const Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons, + IntVect ngvect_cons, + IntVect ngvect_vels) +{ + int lev = 0; + + // We do not operate on the z ghost cells + ngvect_cons[2] = 0; + ngvect_vels[2] = 0; + + // Time interpolation + Real dT = bdy_time_interval; + + // Note this is because we define "time" to be time since start_bdy_time + Real time_since_start = time; + + int n_time = static_cast( time_since_start / dT); + Real alpha = (time_since_start - n_time * dT) / dT; + AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); + Real oma = 1.0 - alpha; + + int n_time_p1 = n_time + 1; + if ((time == stop_time-start_time) && (alpha==0)) { + // stop time coincides with final bdy snapshot -- don't try to read in + // another snapshot + n_time_p1 = n_time; + } + + // Flags for read vars and index mapping + Vector cons_read = {0, 1, 0, 0, + 1, 0, 0, + 0, 0, 0}; + + Vector cons_map = {Rho_comp, RealBdyVars::T, RhoKE_comp, RhoScalar_comp, + RealBdyVars::QV, RhoQ2_comp, RhoQ3_comp, + RhoQ4_comp, RhoQ5_comp, RhoQ6_comp}; + + 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( {RealBdyVars::U} ); // xvel + ind_map.push_back( {RealBdyVars::V} ); // yvel + ind_map.push_back( {0} ); // zvel NOTE: Loop is forward (need u/v filled first) + + // 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]; + + mf.FillBoundary(geom[lev].periodicity()); + + + // 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 = lbound(domain); + const auto& dom_hi = ubound(domain); + + // Offset only applies to cons (we may fill a subset of these vars) + int offset = (var_idx == Vars::cons) ? icomp_cons : 0; + + // Ghost cells to be filled + IntVect ng_vect = (var_idx == Vars::cons) ? ngvect_cons : ngvect_vels; + + // Set region width + int set_width = 1; + + // Loop over each component + for (int comp_idx(offset); comp_idx < (comp_var[var_idx]+offset); ++comp_idx) + { + + // Variable can be read from wrf bdy + //------------------------------------ + if (is_read[var_idx][comp_idx]) + { + int ivar = ind_map[var_idx][comp_idx]; + + // 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_p1][ivar].const_array(); + const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array(); + const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivar].const_array(); + const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array(); + const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivar].const_array(); + const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array(); + const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][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); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + realbdy_bc_bxs_xy(gbx, domain, set_width, + 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); + if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + }, + [=] 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); + if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + }); + + // 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); + if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + }, + [=] 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); + if (var_idx == Vars::cons) dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + }); + } // mfi + + // Variable not read from wrf bdy + //------------------------------------ + } else { +#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); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + realbdy_bc_bxs_xy(gbx, domain, set_width, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, + ng_vect); + + // Bounding + int i_xlo = bx_xlo.bigEnd(0) + set_width; + int i_xhi = bx_xhi.smallEnd(0) - set_width; + int j_ylo = bx_ylo.bigEnd(1) + set_width; + int j_yhi = bx_yhi.smallEnd(1) - set_width; + + // Destination array + const Array4& dest_arr = mf.array(mfi); + + // 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) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i_xlo,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) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i_xhi,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) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i,j_ylo,k,comp_idx); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k,comp_idx) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i,j_yhi,k,comp_idx); + }); + } // mfi + } // is_read + } // comp + } // var +} + +void +ERF::fill_from_realbdy_upwind (const Vector& mfs, + const Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons, + IntVect ngvect_cons, + IntVect ngvect_vels) +{ + int lev = 0; + + // We do not operate on the z ghost cells + ngvect_cons[2] = 0; + ngvect_vels[2] = 0; + + // Time interpolation + Real dT = bdy_time_interval; + + // Note this is because we define "time" to be time since start_bdy_time + Real time_since_start = time; + + int n_time = static_cast( time_since_start / dT); + Real alpha = (time_since_start - n_time * dT) / dT; + AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); + Real oma = 1.0 - alpha; + + int n_time_p1 = n_time + 1; + if ((time == stop_time-start_time) && (alpha==0)) { + // stop time coincides with final bdy snapshot -- don't try to read in + // another snapshot + n_time_p1 = n_time; + } + + // Map loop index to variable index + // NOTE: Loop backwards and do v/u first for WfromOmega + Vector var_idx_map = {Vars::cons, Vars::zvel, Vars::xvel, Vars::yvel}; + + // Variables read from bdy (indexed with var_idx) + Vector> is_read; + is_read.push_back( {0, 1, 0, 0, + 1, 0, 0, + 0, 0, 0} ); + is_read.push_back( {1} ); // xvel + is_read.push_back( {1} ); // yvel + is_read.push_back( {0} ); // zvel + + // Variable component (indexed with var_idx) + Vector> ind_map; + ind_map.push_back( {Rho_comp, RealBdyVars::T, RhoKE_comp, RhoScalar_comp, + RealBdyVars::QV, RhoQ2_comp, RhoQ3_comp, + RhoQ4_comp, RhoQ5_comp, RhoQ6_comp} ); + ind_map.push_back( {RealBdyVars::U} ); // xvel + ind_map.push_back( {RealBdyVars::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 -- note we intentionally do + // the velocities before the scalars since we may use + // the normal velocity on each face for upwinding + for (int var_lp_idx = var_idx_end-1; var_lp_idx >= Vars::cons; --var_lp_idx) + { + int var_idx = var_idx_map[var_lp_idx]; + MultiFab& mf = *mfs[var_idx]; + + mf.FillBoundary(geom[lev].periodicity()); + + // CC domain + Box domain = geom[lev].Domain(); + const auto& dom_cc_lo = lbound(domain); + const auto& dom_cc_hi = ubound(domain); + + // Map domain onto index type of the data + domain.convert(mf.boxArray().ixType()); + const auto& dom_lo = lbound(domain); + const auto& dom_hi = ubound(domain); + IntVect iv = domain.type(); + + // Offset only applies to cons (we may fill a subset of these vars) + int offset = (var_idx == Vars::cons) ? icomp_cons : 0; + + // Ghost cells to be filled + IntVect ng_vect = (var_idx == Vars::cons) ? ngvect_cons : ngvect_vels; + + // Set region width + int set_width = 1; + + // Loop over each component + for (int comp_idx(offset); comp_idx < (comp_var[var_idx]+offset); ++comp_idx) + { + // Variable can be read from wrf bdy + //------------------------------------ + if (is_read[var_idx][comp_idx]) + { + int ivar = ind_map[var_idx][comp_idx]; + + // 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_p1][ivar].const_array(); + const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar].const_array(); + const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivar].const_array(); + const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar].const_array(); + const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivar].const_array(); + const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar].const_array(); + const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][ivar].const_array(); + + // U/V for upwind masking + const auto& bdatxlo_n_m = bdy_data_xlo[n_time ][RealBdyVars::U].const_array(); + const auto& bdatxlo_np1_m = bdy_data_xlo[n_time_p1][RealBdyVars::U].const_array(); + const auto& bdatxhi_n_m = bdy_data_xhi[n_time ][RealBdyVars::U].const_array(); + const auto& bdatxhi_np1_m = bdy_data_xhi[n_time_p1][RealBdyVars::U].const_array(); + const auto& bdatylo_n_m = bdy_data_ylo[n_time ][RealBdyVars::V].const_array(); + const auto& bdatylo_np1_m = bdy_data_ylo[n_time_p1][RealBdyVars::V].const_array(); + const auto& bdatyhi_n_m = bdy_data_yhi[n_time ][RealBdyVars::V].const_array(); + const auto& bdatyhi_np1_m = bdy_data_yhi[n_time_p1][RealBdyVars::V].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); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + realbdy_bc_bxs_xy(gbx, domain, set_width, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, + ng_vect); + + // NOTE: Xlo/hi boxes own corner cells (Ylo/hi) + ParallelFor(bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Limit for BDY FAB data + int ii = std::max(i , dom_lo.x); + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + + // Compute bdy velocities + int jju = std::min(std::max(j, dom_cc_lo.y), dom_cc_hi.y); + int iiv = std::min(std::max(i, dom_cc_lo.x), dom_cc_hi.x); + Real u_bdy = oma * bdatxlo_n_m (dom_cc_lo.x,jju,k,0) + + alpha * bdatxlo_np1_m(dom_cc_lo.x,jju,k,0); + Real v_bdy_lo = oma * bdatylo_n_m (iiv,dom_cc_lo.y,k,0) + + alpha * bdatylo_np1_m(iiv,dom_cc_lo.y,k,0); + Real v_bdy_hi = oma * bdatyhi_n_m (iiv,dom_cc_hi.y+1,k,0) + + alpha * bdatyhi_np1_m(iiv,dom_cc_hi.y+1,k,0); + + if ( (u_bdy >= 0.0) || + ((jj == dom_cc_lo.y ) && (v_bdy_lo >= 0.0)) || + ((jj == dom_cc_hi.y+iv[1]) && (v_bdy_hi <= 0.0)) ) { + dest_arr(i,j,k,comp_idx) = oma * bdatxlo_n (ii,jj,k,0) + + alpha * bdatxlo_np1(ii,jj,k,0); + if (var_idx == Vars::cons) { + dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + } + } else { + dest_arr(i,j,k,comp_idx) = dest_arr(dom_lo.x,jj,k,comp_idx); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Limit for BDY FAB data + int ii = std::min(i , dom_hi.x); + int jj = std::max(j , dom_lo.y); + jj = std::min(jj, dom_hi.y); + + // Compute bdy velocities + int jju = std::min(std::max(j, dom_cc_lo.y), dom_cc_hi.y); + int iiv = std::min(std::max(i, dom_cc_lo.x), dom_cc_hi.x); + Real u_bdy = oma * bdatxhi_n_m (dom_cc_hi.x+1,jju,k,0) + + alpha * bdatxhi_np1_m(dom_cc_hi.x+1,jju,k,0); + Real v_bdy_lo = oma * bdatylo_n_m (iiv,dom_cc_lo.y,k,0) + + alpha * bdatylo_np1_m(iiv,dom_cc_lo.y,k,0); + Real v_bdy_hi = oma * bdatyhi_n_m (iiv,dom_cc_hi.y+1,k,0) + + alpha * bdatyhi_np1_m(iiv,dom_cc_hi.y+1,k,0); + + if ( (u_bdy <= 0.0) || + ((jj == dom_cc_lo.y ) && (v_bdy_lo >= 0.0)) || + ((jj == dom_cc_hi.y+iv[1]) && (v_bdy_hi <= 0.0)) ) { + dest_arr(i,j,k,comp_idx) = oma * bdatxhi_n (ii,jj,k,0) + + alpha * bdatxhi_np1(ii,jj,k,0); + if (var_idx == Vars::cons) { + dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + } + } else { + dest_arr(i,j,k,comp_idx) = dest_arr(dom_hi.x,jj,k,comp_idx); + } + }); + + // NOTE: Ylo/hi boxes do not own corner cells + ParallelFor(bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Limit for BDY FAB data + int jj = std::max(j, dom_lo.y); + + // Compute bdy velocities + int iiv = std::min(std::max(i, dom_cc_lo.x), dom_cc_hi.x); + Real v_bdy = oma * bdatylo_n_m (iiv,dom_cc_lo.y,k,0) + + alpha * bdatylo_np1_m(iiv,dom_cc_lo.y,k,0); + + if (v_bdy >= 0.0) { + dest_arr(i,j,k,comp_idx) = oma * bdatylo_n (i,jj,k,0) + + alpha * bdatylo_np1(i,jj,k,0); + if (var_idx == Vars::cons) { + dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + } + } else { + dest_arr(i,j,k,comp_idx) = dest_arr(i,dom_lo.y,k,comp_idx); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Limit for BDY FAB data + int jj = std::min(j, dom_hi.y); + + // Compute bdy velocities + int iiv = std::min(std::max(i, dom_cc_lo.x), dom_cc_hi.x); + Real v_bdy = oma * bdatyhi_n_m (iiv,dom_cc_hi.y+1,k,0) + + alpha * bdatyhi_np1_m(iiv,dom_cc_hi.y+1,k,0); + + if (v_bdy <= 0.0) { + dest_arr(i,j,k,comp_idx) = oma * bdatyhi_n (i,jj,k,0) + + alpha * bdatyhi_np1(i,jj,k,0); + if (var_idx == Vars::cons) { + dest_arr(i,j,k,comp_idx) *= dest_arr(i,j,k,Rho_comp); + } + } else { + dest_arr(i,j,k,comp_idx) = dest_arr(i,dom_hi.y,k,comp_idx); + } + }); + } // mfi + + // Variable not read from wrf bdy + //------------------------------------ + } else { +#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); + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + realbdy_bc_bxs_xy(gbx, domain, set_width, + bx_xlo, bx_xhi, + bx_ylo, bx_yhi, + ng_vect); + + // Bounding + int i_xlo = bx_xlo.bigEnd(0) + set_width; + int i_xhi = bx_xhi.smallEnd(0) - set_width; + int j_ylo = bx_ylo.bigEnd(1) + set_width; + int j_yhi = bx_yhi.smallEnd(1) - set_width; + + // Destination array + const Array4& dest_arr = mf.array(mfi); + + // 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) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i_xlo,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) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i_xhi,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) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i,j_ylo,k,comp_idx); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + dest_arr(i,j,k,comp_idx) = (comp_idx > RhoQ1_comp) ? Real(0.) : + dest_arr(i,j_yhi,k,comp_idx); + }); + } // mfi + } // is_read + } // comp + } // var +} +#endif diff --git a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp index 15e3031613..279f50e639 100644 --- a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp +++ b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp @@ -320,7 +320,18 @@ ERF::FillIntermediatePatch (int lev, Real time, bool do_fb = true; - if (m_r2d && !solverChoice.use_real_bcs) fill_from_bndryregs(mfs_vel,time); +#ifdef ERF_USE_NETCDF + if (solverChoice.use_real_bcs && (lev==0)) { + if (solverChoice.upwind_real_bcs) { + fill_from_realbdy_upwind(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels); + } else { + fill_from_realbdy(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels); + } + do_fb = false; + } +#endif + + if (m_r2d && !solverChoice.use_real_bcs) { fill_from_bndryregs(mfs_vel,time); } // We call this even if use_real_bcs is true because these will fill the vertical bcs (*physbcs_cons[lev])(*mfs_vel[Vars::cons],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel], diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 5b84c747e4..cb339ad3d1 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -352,7 +352,18 @@ ERF::FillPatchCrseLevel (int lev, Real time, bool do_fb = true; - if (m_r2d && !solverChoice.use_real_bcs) fill_from_bndryregs(mfs_vel,time); +#ifdef ERF_USE_NETCDF + if(solverChoice.use_real_bcs && (lev==0)) { + if (solverChoice.upwind_real_bcs) { + fill_from_realbdy_upwind(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels); + } else { + fill_from_realbdy(mfs_vel,time,cons_only,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels); + } + do_fb = false; + } +#endif + + if (m_r2d && !solverChoice.use_real_bcs) { fill_from_bndryregs(mfs_vel,time); } // We call this even if use_real_bcs is true because these will fill the vertical bcs // Note that we call FillBoundary inside the physbcs call diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index 976af98f24..236d99b39b 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -62,19 +62,19 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yv Array4 z_nd_arr; - if (m_z_phys_nd) - { + if (m_z_phys_nd) { z_nd_arr = m_z_phys_nd->const_array(mfi); } - if (!gdomain.contains(cbx2)) - { + if (!gdomain.contains(cbx2)) { Array4< Real> const& cons_arr = mf.array(mfi); Array4 const& velx_arr = xvel.const_array(mfi); Array4 const& vely_arr = yvel.const_array(mfi); - // We send a box with ghost cells in the lateral directions only - impose_lateral_cons_bcs(cons_arr,velx_arr,vely_arr,cbx1,domain,icomp,ncomp,nghost,time); + if (!m_use_real_bcs) { + // We send a box with ghost cells in the lateral directions only + impose_lateral_cons_bcs(cons_arr,velx_arr,vely_arr,cbx1,domain,icomp,ncomp,nghost,time); + } // We send the full FAB box with ghost cells impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp,time,do_terrain_adjustment); @@ -134,19 +134,16 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, Array4 z_nd_arr; - if (m_z_phys_nd) - { + if (m_z_phys_nd) { z_nd_arr = m_z_phys_nd->const_array(mfi); } - if (!gdomainx.contains(xbx2)) - { + if (!gdomainx.contains(xbx2)) { Array4< Real> const& dest_arr = mf.array(mfi); Array4 const& velx_arr = xvel.const_array(mfi); Array4 const& vely_arr = yvel.const_array(mfi); - if (!gdomainx.contains(xbx1)) - { + if (!gdomainx.contains(xbx1) && !m_use_real_bcs) { impose_lateral_xvel_bcs(dest_arr,velx_arr,vely_arr,xbx1,domain,bccomp,time); } @@ -206,18 +203,18 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, Array4 z_nd_arr; - if (m_z_phys_nd) - { + if (m_z_phys_nd) { z_nd_arr = m_z_phys_nd->const_array(mfi); } - if (!gdomainy.contains(ybx2)) - { + if (!gdomainy.contains(ybx2)) { Array4< Real> const& dest_arr = mf.array(mfi); Array4 const& velx_arr = xvel.const_array(mfi); Array4 const& vely_arr = yvel.const_array(mfi); - impose_lateral_yvel_bcs(dest_arr,velx_arr,vely_arr,ybx1,domain,bccomp,time); + if (!m_use_real_bcs) { + impose_lateral_yvel_bcs(dest_arr,velx_arr,vely_arr,ybx1,domain,bccomp,time); + } impose_vertical_yvel_bcs(dest_arr,ybx2,domain,z_nd_arr,dxInv,bccomp,time); } @@ -289,22 +286,19 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel, const Array4& mf_u = m_mapfac_u->const_array(mfi); const Array4& mf_v = m_mapfac_v->const_array(mfi); - if (m_z_phys_nd) - { + if (m_z_phys_nd) { z_nd_arr = m_z_phys_nd->const_array(mfi); } // // Recall that gdomainz.smallEnd(2) = 1 not 0! // - if (!gdomainz.contains(zbx)) - { + if (!gdomainz.contains(zbx)) { Array4 const& velx_arr = xvel.const_array(mfi); Array4 const& vely_arr = yvel.const_array(mfi); Array4< Real> const& velz_arr = mf.array(mfi); - if (!gdomainz.contains(zbx)) - { + if (!gdomainz.contains(zbx) && !m_use_real_bcs) { impose_lateral_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx,domain, mf_u,mf_v,z_nd_arr,dxInv,m_terrain_type,bccomp_w,time); } @@ -358,13 +352,11 @@ void ERFPhysBCFunct_base::operator() (MultiFab& mf, int /*icomp*/, int ncomp, In Array4 z_nd_arr; - if (m_z_phys_nd) - { + if (m_z_phys_nd) { z_nd_arr = m_z_phys_nd->const_array(mfi); } - if (!gdomain.contains(cbx2)) - { + if (!gdomain.contains(cbx2)) { const Array4 base_arr = mf.array(mfi); impose_lateral_basestate_bcs(base_arr,cbx1,domain,ncomp,nghost); diff --git a/Source/BoundaryConditions/Make.package b/Source/BoundaryConditions/Make.package index dc5491df72..7386f700cb 100644 --- a/Source/BoundaryConditions/Make.package +++ b/Source/BoundaryConditions/Make.package @@ -4,6 +4,7 @@ CEXE_sources += ERF_BoundaryConditionsZvel.cpp CEXE_sources += ERF_BoundaryConditionsCons.cpp CEXE_sources += ERF_BoundaryConditionsBaseState.cpp CEXE_sources += ERF_BoundaryConditionsBndryReg.cpp +CEXE_sources += ERF_BoundaryConditionsRealbdy.cpp CEXE_headers += ERF_SurfaceLayer.H CEXE_sources += ERF_SurfaceLayer.cpp diff --git a/Source/ERF.H b/Source/ERF.H index 61289e60ad..66e9ae7cd4 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -560,6 +560,26 @@ public: void fill_from_bndryregs (const amrex::Vector& mfs, amrex::Real time); +#ifdef ERF_USE_NETCDF + // Version in which we always set the velocity values but don't set theta + void fill_from_realbdy (const amrex::Vector& mfs, + amrex::Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons, + amrex::IntVect ngvect_cons, + amrex::IntVect ngvect_vels); + + // Version in which we only set the boundary value if inflowing, includes theta + void fill_from_realbdy_upwind (const amrex::Vector& mfs, + amrex::Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons, + amrex::IntVect ngvect_cons, + amrex::IntVect ngvect_vels); +#endif + #ifdef ERF_USE_NETCDF void init_from_wrfinput (int lev, amrex::MultiFab& mf_C1H, diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 1903db87fc..ac64e9c998 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1357,6 +1357,23 @@ ERF::InitData_post () int ncomp_cons = lev_new[Vars::cons].nComp(); bool do_fb = true; +#ifdef ERF_USE_NETCDF + if (solverChoice.use_real_bcs && (lev==0)) { + int icomp_cons = 0; + bool cons_only = false; + Vector mfs_vec = {&lev_new[Vars::cons],&lev_new[Vars::xvel], + &lev_new[Vars::yvel],&lev_new[Vars::zvel]}; + if (solverChoice.upwind_real_bcs) { + fill_from_realbdy_upwind(mfs_vec,t_new[lev],cons_only,icomp_cons, + ncomp_cons,ngvect_cons,ngvect_vels); + } else { + fill_from_realbdy(mfs_vec,t_new[lev],cons_only,icomp_cons, + ncomp_cons,ngvect_cons,ngvect_vels); + } + do_fb = false; + } +#endif + (*physbcs_cons[lev])(lev_new[Vars::cons],lev_new[Vars::xvel],lev_new[Vars::yvel],0,ncomp_cons, ngvect_cons,t_new[lev],BCVars::cons_bc,do_fb); ( *physbcs_u[lev])(lev_new[Vars::xvel],lev_new[Vars::xvel],lev_new[Vars::yvel], diff --git a/Source/Utils/ERF_InteriorGhostCells.cpp b/Source/Utils/ERF_InteriorGhostCells.cpp index 821e963279..5f93f79f2d 100644 --- a/Source/Utils/ERF_InteriorGhostCells.cpp +++ b/Source/Utils/ERF_InteriorGhostCells.cpp @@ -1,4 +1,4 @@ -#include +#include "ERF_Utils.H" using namespace amrex; @@ -74,6 +74,71 @@ realbdy_interior_bxs_xy (const Box& bx, } +/** + * Get the boxes for looping over interior/exterior ghost cells + * for use by fillpatch, erf_slow_rhs_pre, and erf_slow_rhs_post. + * + * @param[in] bx box to intersect with 4 halo regions + * @param[in] domain box of the whole domain + * @param[in] width number of cells in (relaxation+specified) zone + * @param[in] set_width number of cells in (specified) zone + * @param[out] bx_xlo halo box at x_lo boundary + * @param[out] bx_xhi halo box at x_hi boundary + * @param[out] bx_ylo halo box at y_lo boundary + * @param[out] bx_yhi halo box at y_hi boundary + * @param[in] ng_vect number of ghost cells in each direction + * @param[in] get_int_ng flag to get ghost cells inside the domain + */ +void +realbdy_bc_bxs_xy (const Box& bx, + const Box& domain, + const int& set_width, + Box& bx_xlo, + Box& bx_xhi, + Box& bx_ylo, + Box& bx_yhi, + const IntVect& ng_vect) +{ + AMREX_ALWAYS_ASSERT(bx.ixType() == domain.ixType()); + + // Domain bounds without ghost cells + const auto& dom_lo = lbound(domain); + const auto& dom_hi = ubound(domain); + + // Four boxes matching the domain + Box gdom_xlo(domain); Box gdom_xhi(domain); + Box gdom_ylo(domain); Box gdom_yhi(domain); + + // Get offsets from box index type + IntVect iv_type = bx.ixType().toIntVect(); + int offx = (iv_type[0]==1) ? 0 : -1; + int offy = (iv_type[1]==1) ? 0 : -1; + + // Stagger the boxes based upon index type + gdom_xlo += IntVect(offx,0,0); gdom_xhi += IntVect(-offx,0,0); + gdom_ylo += IntVect(0,offy,0); gdom_yhi += IntVect(0,-offy,0); + + // Trim the boxes to only include internal ghost cells + gdom_xlo.setBig(0,dom_lo.x+set_width+offx-1); gdom_xhi.setSmall(0,dom_hi.x-set_width-offx+1); + gdom_ylo.setBig(1,dom_lo.y+set_width+offy-1); gdom_yhi.setSmall(1,dom_hi.y-set_width-offy+1); + + // Remove overlapping corners from y-face boxes + gdom_ylo.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_ylo.setBig(0,gdom_xhi.smallEnd(0)-1); + gdom_yhi.setSmall(0,gdom_xlo.bigEnd(0)+1); gdom_yhi.setBig(0,gdom_xhi.smallEnd(0)-1); + + // Grow boxes to get external ghost cells only + gdom_xlo.growLo(0,ng_vect[0]+offx); gdom_xhi.growHi(0,ng_vect[0]+offx); + gdom_xlo.grow (1,ng_vect[1] ); gdom_xhi.grow (1,ng_vect[1] ); + gdom_ylo.growLo(1,ng_vect[1]+offy); gdom_yhi.growHi(1,ng_vect[1]+offy); + + // Populate everything + bx_xlo = (bx & gdom_xlo); + bx_xhi = (bx & gdom_xhi); + bx_ylo = (bx & gdom_ylo); + bx_yhi = (bx & gdom_yhi); +} + + /** * Compute the RHS in the relaxation zone * @@ -408,7 +473,7 @@ realbdy_compute_interior_ghost_rhs (const Real& time, tbx_ylo, tbx_yhi, ng_vect); - Array4 rhs_arr; Array4 data_arr; + Array4 rhs_arr; Array4 data_arr; Array4 arr_xlo; Array4 arr_xhi; Array4 arr_ylo; Array4 arr_yhi; if (ivar == ivarU) { @@ -440,49 +505,94 @@ realbdy_compute_interior_ghost_rhs (const Real& time, arr_xlo , arr_xhi , arr_ylo , arr_yhi , u_xlo, u_xhi, v_xlo, v_xhi, v_ylo, v_yhi, data_arr, rhs_arr, do_upwind); + } // mfi + } // ivar - /* - // UNIT TEST DEBUG - realbdy_interior_bxs_xy(tbx, domain, width, - tbx_xlo, tbx_xhi, - tbx_ylo, tbx_yhi); - ParallelFor(tbx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + // Set normal velocity RHS at the boundary + //========================================================== +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(S_cur_data[ivarT],TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Box tbx = mfi.nodaltilebox(0); + Box tby = mfi.nodaltilebox(1); + + Box domain = geom.Domain(); + Box domainx = convert(domain, IntVect(1,0,0)); + Box domainy = convert(domain, IntVect(0,1,0)); + + int ilo = domainx.smallEnd(0); + int ihi = domainx.bigEnd(0); + int jlo = domainy.smallEnd(1); + int jhi = domainy.bigEnd(1); + + Box tbx_lo, tbx_hi; + if (tbx.smallEnd(0) == ilo) { + tbx_lo = makeSlab(tbx,0,ilo); + } else if (tbx.bigEnd(0) == ihi) { + tbx_hi = makeSlab(tbx,0,ihi); + } + + Box tby_lo, tby_hi; + if (tby.smallEnd(1) == jlo) { + tby_lo = makeSlab(tby,1,jlo); + } else if (tby.bigEnd(1) == jhi) { + tby_hi = makeSlab(tby,1,jhi); + } + + Array4 rhs_xmom = S_rhs[IntVars::xmom].array(mfi); + Array4 rhs_ymom = S_rhs[IntVars::ymom].array(mfi); + + Array4 rhs_cons = S_rhs[IntVars::cons].const_array(mfi); + Array4 cons_arr = S_cur_data[IntVars::cons].const_array(mfi); + + const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivarU].const_array(); + const auto& bdatxlo_np1 = bdy_data_xlo[n_time_p1][ivarU].const_array(); + const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivarU].const_array(); + const auto& bdatxhi_np1 = bdy_data_xhi[n_time_p1][ivarU].const_array(); + + const auto& bdatylo_n = bdy_data_ylo[n_time ][ivarV].const_array(); + const auto& bdatylo_np1 = bdy_data_ylo[n_time_p1][ivarV].const_array(); + const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivarV].const_array(); + const auto& bdatyhi_np1 = bdy_data_yhi[n_time_p1][ivarV].const_array(); + + ParallelFor(tbx_lo, tbx_hi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - if (std::fabs(arr_xlo(i,j,k) - data_arr(i,j,k,icomp)) > Real(0.01)) { - Print() << "ERROR XLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; - Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_xlo(i,j,k) << "\n"; - exit(0); - } - }); - ParallelFor(tbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + Real rho_tend = rhs_cons(i,j,k); + Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k)); + Real u_tend = (bdatxlo_np1(i,j,k) - bdatxlo_n(i,j,k)) / bdy_time_interval; + Real u_val = oma * bdatxlo_n (i,j,k) + alpha * bdatxlo_np1(i,j,k); + rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend; + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - if (std::fabs(arr_xhi(i,j,k) - data_arr(i,j,k,icomp)) > Real(0.01)) { - Print() << "ERROR XHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; - Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_xhi(i,j,k) << "\n"; - exit(0); - } + Real rho_tend = rhs_cons(i,j,k); + Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k)); + Real u_tend = (bdatxhi_np1(i,j,k) - bdatxhi_n(i,j,k)) / bdy_time_interval; + Real u_val = oma * bdatxhi_n (i,j,k) + alpha * bdatxhi_np1(i,j,k); + rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend; }); - ParallelFor(tbx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + + ParallelFor(tby_lo, tby_hi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - if (std::fabs(arr_ylo(i,j,k) - data_arr(i,j,k,icomp)) > Real(0.01)) { - Print() << "ERROR YLO: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; - Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_ylo(i,j,k) << "\n"; - exit(0); - } - }); - ParallelFor(tbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + Real rho_tend = rhs_cons(i,j,k); + Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i,j-1,k)); + Real v_tend = (bdatylo_np1(i,j,k) - bdatylo_n(i,j,k)) / bdy_time_interval; + Real v_val = oma * bdatxlo_n (i,j,k) + alpha * bdatxlo_np1(i,j,k); + rhs_ymom(i,j,k) = rho_val * v_tend + v_val * rho_tend; + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - if (std::fabs(arr_yhi(i,j,k)-data_arr(i,j,k,icomp)) > Real(0.01)) { - Print() << "ERROR YHI: " << ivar << ' ' << icomp << ' ' << IntVect(i,j,k) << "\n"; - Print() << "DATA: " << data_arr(i,j,k,icomp) << ' ' << arr_yhi(i,j,k) << "\n"; - exit(0); - } + Real rho_tend = rhs_cons(i,j,k); + Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i,j-1,k)); + Real v_tend = (bdatyhi_np1(i,j,k) - bdatyhi_n(i,j,k)) / bdy_time_interval; + Real v_val = oma * bdatyhi_n(i,j,k) + alpha * bdatyhi_np1(i,j,k); + rhs_ymom(i,j,k) = rho_val * v_tend + v_val * rho_tend; }); - */ + } // mfi - } // ivar - //ParallelDescriptor::Barrier(); - //exit(0); } /** diff --git a/Source/Utils/ERF_Utils.H b/Source/Utils/ERF_Utils.H index 9e6dbad75f..5601774c26 100644 --- a/Source/Utils/ERF_Utils.H +++ b/Source/Utils/ERF_Utils.H @@ -106,6 +106,19 @@ void realbdy_interior_bxs_xy (const amrex::Box& bx, const amrex::IntVect& ng_vect=amrex::IntVect(0,0,0), const bool get_int_ng=false); +/* + * Compute boxes for looping over set region cells + * for use by fillpatch, erf_slow_rhs_pre, and erf_slow_rhs_post + */ +void realbdy_bc_bxs_xy (const amrex::Box& bx, + const amrex::Box& domain, + const int& set_width, + amrex::Box& bx_xlo, + amrex::Box& bx_xhi, + amrex::Box& bx_ylo, + amrex::Box& bx_yhi, + const amrex::IntVect& ng_vect=amrex::IntVect(0,0,0)); + /* * Compute relaxation region RHS with wrfbdy */ From f32008acfb2aef060134e0b7cd85ce8e1ada7ebb Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 13:14:17 -0700 Subject: [PATCH 10/27] Do microphysics everywhere. --- Source/Microphysics/ERF_EulerianMicrophysics.H | 7 ------- Source/Microphysics/ERF_LagrangianMicrophysics.H | 7 ------- Source/Microphysics/ERF_Microphysics.H | 4 ---- Source/Microphysics/Kessler/ERF_Kessler.H | 6 ------ Source/Microphysics/Kessler/ERF_Kessler.cpp | 8 -------- Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp | 8 -------- Source/Microphysics/Morrison/ERF_Morrison.H | 6 ------ Source/Microphysics/Null/ERF_NullMoist.H | 4 ---- Source/Microphysics/SAM/ERF_CloudSAM.cpp | 4 ---- Source/Microphysics/SAM/ERF_Precip.cpp | 4 ---- Source/Microphysics/SAM/ERF_SAM.H | 6 ------ Source/Microphysics/SatAdj/ERF_SatAdj.H | 6 ------ Source/Microphysics/SatAdj/ERF_SatAdj.cpp | 4 ---- Source/Microphysics/WSM6/ERF_WSM6.H | 3 --- Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp | 1 - 15 files changed, 78 deletions(-) diff --git a/Source/Microphysics/ERF_EulerianMicrophysics.H b/Source/Microphysics/ERF_EulerianMicrophysics.H index c6d364cd4d..ca381c2861 100644 --- a/Source/Microphysics/ERF_EulerianMicrophysics.H +++ b/Source/Microphysics/ERF_EulerianMicrophysics.H @@ -160,13 +160,6 @@ public: m_moist_model[a_lev]->Set_dzmin(dz_min); } - /*! \brief Populates real_width in micro model for box limiting */ - virtual void Set_RealWidth (const int a_lev, - const int real_width) const override - { - m_moist_model[a_lev]->Set_RealWidth(real_width); - } - protected: /*! \brief Create and set the specified moisture model */ diff --git a/Source/Microphysics/ERF_LagrangianMicrophysics.H b/Source/Microphysics/ERF_LagrangianMicrophysics.H index d59eaff851..ad7b0afbb1 100644 --- a/Source/Microphysics/ERF_LagrangianMicrophysics.H +++ b/Source/Microphysics/ERF_LagrangianMicrophysics.H @@ -212,13 +212,6 @@ public: m_moist_model->Set_dzmin(dz_min); } - /*! \brief Populates real_width in micro model for box limiting */ - virtual void Set_RealWidth (const int /*a_lev*/, - const int real_width) const override - { - m_moist_model->Set_RealWidth(real_width); - } - protected: /*! \brief Create and set the specified moisture model */ diff --git a/Source/Microphysics/ERF_Microphysics.H b/Source/Microphysics/ERF_Microphysics.H index 1c67258cc3..9b5db35699 100644 --- a/Source/Microphysics/ERF_Microphysics.H +++ b/Source/Microphysics/ERF_Microphysics.H @@ -95,10 +95,6 @@ public: virtual void Set_dzmin (const int lev, const amrex::Real dz_min) const = 0; - /*! \brief Define the nudging+set width for real BCs */ - virtual void Set_RealWidth (const int lev, - const int real_width) const = 0; - /*! \brief query if a specified moisture model is Eulerian or Lagrangian */ static MoistureModelType modelType (const MoistureType a_moisture_type) { diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index c6ace4afc2..883bb2debf 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -121,9 +121,6 @@ public: int Qstate_Moist_Size () override { return Kessler::n_qstate_moist_size; } - void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } - void Qmoist_Restart_Vars (const SolverChoice& /*a_sc*/, std::vector& a_idx, @@ -152,9 +149,6 @@ private: // geometry amrex::Geometry m_geom; - // Nudging + set width - int m_real_width{0}; - // timestep amrex::Real dt; diff --git a/Source/Microphysics/Kessler/ERF_Kessler.cpp b/Source/Microphysics/Kessler/ERF_Kessler.cpp index 3aeaf47d11..5c2fe13431 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.cpp +++ b/Source/Microphysics/Kessler/ERF_Kessler.cpp @@ -43,10 +43,6 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); auto tbx = mfi.tilebox(); - if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); } - if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); } - if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); } - if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); } Real d_fac_cond = m_fac_cond; @@ -220,10 +216,6 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) auto pres_array = mic_fab_vars[MicVar_Kess::pres]->array(mfi); auto tbx = mfi.tilebox(); - if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); } - if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); } - if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); } - if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); } Real d_fac_cond = m_fac_cond; diff --git a/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp b/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp index 102b5ba8e6..c6746cbb22 100644 --- a/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp +++ b/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp @@ -531,14 +531,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE { auto box = mfi.tilebox(); -#ifndef ERF_USE_MORR_FORT - // NOTE: Trimming the box with FORTRAN breaks the pointer dereferencing - if (box.smallEnd(0) == i_lo) { box.growLo(0,-m_real_width); } - if (box.bigEnd(0) == i_hi) { box.growHi(0,-m_real_width); } - if (box.smallEnd(1) == j_lo) { box.growLo(1,-m_real_width); } - if (box.bigEnd(1) == j_hi) { box.growHi(1,-m_real_width); } -#endif - if (!box.ok()) { // Avoid going farther if the box is inverted (i.e., ilo > ihi or jlo > jhi). continue; } diff --git a/Source/Microphysics/Morrison/ERF_Morrison.H b/Source/Microphysics/Morrison/ERF_Morrison.H index a4d22b3d0c..bc06dbd979 100644 --- a/Source/Microphysics/Morrison/ERF_Morrison.H +++ b/Source/Microphysics/Morrison/ERF_Morrison.H @@ -134,9 +134,6 @@ public: int Qstate_Moist_Size () override { return Morrison::n_qstate_moist_size; } - void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real @@ -300,9 +297,6 @@ private: // geometry amrex::Geometry m_geom; - // Nudging + set width - int m_real_width{0}; - // number of vertical levels int nlev, zlo, zhi; diff --git a/Source/Microphysics/Null/ERF_NullMoist.H b/Source/Microphysics/Null/ERF_NullMoist.H index a229a2cc56..f29d2a354d 100644 --- a/Source/Microphysics/Null/ERF_NullMoist.H +++ b/Source/Microphysics/Null/ERF_NullMoist.H @@ -110,10 +110,6 @@ public: void Set_dzmin (const amrex::Real /*dz_min*/) { } - virtual - void - Set_RealWidth (const int /*real_width*/) { } - private: int m_qmoist_size = 0; diff --git a/Source/Microphysics/SAM/ERF_CloudSAM.cpp b/Source/Microphysics/SAM/ERF_CloudSAM.cpp index 34b92576d3..ba48b3f9b3 100644 --- a/Source/Microphysics/SAM/ERF_CloudSAM.cpp +++ b/Source/Microphysics/SAM/ERF_CloudSAM.cpp @@ -45,10 +45,6 @@ SAM::Cloud (const SolverChoice& sc) auto pres_array = mic_fab_vars[MicVar::pres]->array(mfi); auto tbx = mfi.tilebox(); - if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); } - if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); } - if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); } - if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); } ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { diff --git a/Source/Microphysics/SAM/ERF_Precip.cpp b/Source/Microphysics/SAM/ERF_Precip.cpp index 28ab3c5d9b..31b4afb1c9 100644 --- a/Source/Microphysics/SAM/ERF_Precip.cpp +++ b/Source/Microphysics/SAM/ERF_Precip.cpp @@ -72,10 +72,6 @@ SAM::Precip (const SolverChoice& sc) auto qp_array = mic_fab_vars[MicVar::qp]->array(mfi); auto tbx = mfi.tilebox(); - if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); } - if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); } - if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); } - if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); } ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { diff --git a/Source/Microphysics/SAM/ERF_SAM.H b/Source/Microphysics/SAM/ERF_SAM.H index 637f14ab18..6f1fe2e87b 100644 --- a/Source/Microphysics/SAM/ERF_SAM.H +++ b/Source/Microphysics/SAM/ERF_SAM.H @@ -153,9 +153,6 @@ public: int Qstate_Moist_Size () override { return SAM::n_qstate_moist_size; } - void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real @@ -301,9 +298,6 @@ private: // geometry amrex::Geometry m_geom; - // Nudging + set width - int m_real_width{0}; - // timestep amrex::Real dt; diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.H b/Source/Microphysics/SatAdj/ERF_SatAdj.H index 707dbcc148..2e1ba8e4de 100644 --- a/Source/Microphysics/SatAdj/ERF_SatAdj.H +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.H @@ -125,9 +125,6 @@ public: a_names.clear(); } - void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real @@ -264,9 +261,6 @@ private: // geometry amrex::Geometry m_geom; - // Nudging + set width - int m_real_width{0}; - // timestep amrex::Real dt; diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.cpp b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp index 98f818b2ce..cda74e8e8a 100644 --- a/Source/Microphysics/SatAdj/ERF_SatAdj.cpp +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp @@ -25,10 +25,6 @@ void SatAdj::AdvanceSatAdj (const SolverChoice& /*solverChoice*/) for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto tbx = mfi.tilebox(); - if (tbx.smallEnd(0) == i_lo) { tbx.growLo(0,-m_real_width); } - if (tbx.bigEnd(0) == i_hi) { tbx.growHi(0,-m_real_width); } - if (tbx.smallEnd(1) == j_lo) { tbx.growLo(1,-m_real_width); } - if (tbx.bigEnd(1) == j_hi) { tbx.growHi(1,-m_real_width); } auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); diff --git a/Source/Microphysics/WSM6/ERF_WSM6.H b/Source/Microphysics/WSM6/ERF_WSM6.H index 26b0d9ce6d..cb957c0f15 100644 --- a/Source/Microphysics/WSM6/ERF_WSM6.H +++ b/Source/Microphysics/WSM6/ERF_WSM6.H @@ -82,8 +82,6 @@ public: int Qmoist_Size() override { return m_qmoist_size; } int Qstate_Moist_Size() override { return n_qstate_moist_size; } - void Set_RealWidth(const int real_width) override { m_real_width = real_width; } - void Qmoist_Restart_Vars(const SolverChoice&, std::vector& a_idx, std::vector& a_names) const override @@ -99,7 +97,6 @@ private: amrex::Vector MicVarMap; amrex::Geometry m_geom; - int m_real_width{0}; amrex::Real dt{0.0}; amrex::Real m_dzmin{0.0}; int nlev{0}, zlo{0}, zhi{0}; diff --git a/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp b/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp index dbb2cf62c4..2e7160c59f 100644 --- a/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp +++ b/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp @@ -9,7 +9,6 @@ void ERF::advance_microphysics (int lev, const Real& time ) { if (solverChoice.moisture_type != MoistureType::None) { - micro->Set_RealWidth(lev, real_width); if (lev > 0) { MultiFab& U_new = vars_new[lev][Vars::xvel]; MultiFab& V_new = vars_new[lev][Vars::yvel]; From 09edfcd4cc1a7bf0a0c251ad76f7e7968df75df3 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 13:20:52 -0700 Subject: [PATCH 11/27] remove unused. --- Source/Microphysics/SatAdj/ERF_SatAdj.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.cpp b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp index cda74e8e8a..772070f110 100644 --- a/Source/Microphysics/SatAdj/ERF_SatAdj.cpp +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp @@ -15,12 +15,6 @@ void SatAdj::AdvanceSatAdj (const SolverChoice& /*solverChoice*/) Real d_fac_cond = m_fac_cond; Real rdOcp = m_rdOcp; - auto domain = m_geom.Domain(); - int i_lo = domain.smallEnd(0); - int i_hi = domain.bigEnd(0); - int j_lo = domain.smallEnd(1); - int j_hi = domain.bigEnd(1); - // get the temperature, density, theta, qt and qc from input for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { From e84528bc9cd3ef21eb0848717a7cb34c53c0b65f Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 13:24:20 -0700 Subject: [PATCH 12/27] More unused. --- Source/Microphysics/Kessler/ERF_Kessler.cpp | 4 ---- Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp | 7 ------- Source/Microphysics/SAM/ERF_CloudSAM.cpp | 6 ------ Source/Microphysics/SAM/ERF_Precip.cpp | 6 ------ 4 files changed, 23 deletions(-) diff --git a/Source/Microphysics/Kessler/ERF_Kessler.cpp b/Source/Microphysics/Kessler/ERF_Kessler.cpp index 5c2fe13431..2ec19d95d8 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.cpp +++ b/Source/Microphysics/Kessler/ERF_Kessler.cpp @@ -16,10 +16,6 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) bool do_cond = m_do_cond; auto tabs = mic_fab_vars[MicVar_Kess::tabs]; auto domain = m_geom.Domain(); - int i_lo = domain.smallEnd(0); - int i_hi = domain.bigEnd(0); - int j_lo = domain.smallEnd(1); - int j_hi = domain.bigEnd(1); int k_lo = domain.smallEnd(2); int k_hi = domain.bigEnd(2); if (solverChoice.moisture_type == MoistureType::Kessler) diff --git a/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp b/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp index c6746cbb22..f7b2d05b93 100644 --- a/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp +++ b/Source/Microphysics/Morrison/ERF_AdvanceMorrison.cpp @@ -506,13 +506,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE // Store timestep Real dt = dt_advance; - // Bounds - auto domain = m_geom.Domain(); - int i_lo = domain.smallEnd(0); - int i_hi = domain.bigEnd(0); - int j_lo = domain.smallEnd(1); - int j_hi = domain.bigEnd(1); - // Check if CPP or FORT answer is used ParmParse pp("erf"); bool run_morr_cpp = true; diff --git a/Source/Microphysics/SAM/ERF_CloudSAM.cpp b/Source/Microphysics/SAM/ERF_CloudSAM.cpp index ba48b3f9b3..ba1f1fac8b 100644 --- a/Source/Microphysics/SAM/ERF_CloudSAM.cpp +++ b/Source/Microphysics/SAM/ERF_CloudSAM.cpp @@ -27,12 +27,6 @@ SAM::Cloud (const SolverChoice& sc) SAM_moisture_type = 2; } - auto domain = m_geom.Domain(); - int i_lo = domain.smallEnd(0); - int i_hi = domain.bigEnd(0); - int j_lo = domain.smallEnd(1); - int j_hi = domain.bigEnd(1); - for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto qt_array = mic_fab_vars[MicVar::qt]->array(mfi); auto qn_array = mic_fab_vars[MicVar::qn]->array(mfi); diff --git a/Source/Microphysics/SAM/ERF_Precip.cpp b/Source/Microphysics/SAM/ERF_Precip.cpp index 31b4afb1c9..4054093686 100644 --- a/Source/Microphysics/SAM/ERF_Precip.cpp +++ b/Source/Microphysics/SAM/ERF_Precip.cpp @@ -46,12 +46,6 @@ SAM::Precip (const SolverChoice& sc) SAM_moisture_type = 2; } - auto domain = m_geom.Domain(); - int i_lo = domain.smallEnd(0); - int i_hi = domain.bigEnd(0); - int j_lo = domain.smallEnd(1); - int j_hi = domain.bigEnd(1); - // get the temperature, density, theta, qt and qp from input for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]),TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi); From b6c64ddbfa683461c66546499dedf7e42e402117 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 13:33:02 -0700 Subject: [PATCH 13/27] fix copy and past name error. --- Source/Utils/ERF_InteriorGhostCells.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/Utils/ERF_InteriorGhostCells.cpp b/Source/Utils/ERF_InteriorGhostCells.cpp index 5f93f79f2d..0ed1079f35 100644 --- a/Source/Utils/ERF_InteriorGhostCells.cpp +++ b/Source/Utils/ERF_InteriorGhostCells.cpp @@ -562,7 +562,7 @@ realbdy_compute_interior_ghost_rhs (const Real& time, Real rho_tend = rhs_cons(i,j,k); Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k)); Real u_tend = (bdatxlo_np1(i,j,k) - bdatxlo_n(i,j,k)) / bdy_time_interval; - Real u_val = oma * bdatxlo_n (i,j,k) + alpha * bdatxlo_np1(i,j,k); + Real u_val = oma * bdatxlo_n(i,j,k) + alpha * bdatxlo_np1(i,j,k); rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -570,7 +570,7 @@ realbdy_compute_interior_ghost_rhs (const Real& time, Real rho_tend = rhs_cons(i,j,k); Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i-1,j,k)); Real u_tend = (bdatxhi_np1(i,j,k) - bdatxhi_n(i,j,k)) / bdy_time_interval; - Real u_val = oma * bdatxhi_n (i,j,k) + alpha * bdatxhi_np1(i,j,k); + Real u_val = oma * bdatxhi_n(i,j,k) + alpha * bdatxhi_np1(i,j,k); rhs_xmom(i,j,k) = rho_val * u_tend + u_val * rho_tend; }); @@ -580,7 +580,7 @@ realbdy_compute_interior_ghost_rhs (const Real& time, Real rho_tend = rhs_cons(i,j,k); Real rho_val = Real(0.5) * (cons_arr(i,j,k) + cons_arr(i,j-1,k)); Real v_tend = (bdatylo_np1(i,j,k) - bdatylo_n(i,j,k)) / bdy_time_interval; - Real v_val = oma * bdatxlo_n (i,j,k) + alpha * bdatxlo_np1(i,j,k); + Real v_val = oma * bdatylo_n(i,j,k) + alpha * bdatylo_np1(i,j,k); rhs_ymom(i,j,k) = rho_val * v_tend + v_val * rho_tend; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept From ccaf469251653168ddc4190ff63debf918562755 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 13:46:35 -0700 Subject: [PATCH 14/27] Fix mfiter. --- Source/Utils/ERF_InteriorGhostCells.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Utils/ERF_InteriorGhostCells.cpp b/Source/Utils/ERF_InteriorGhostCells.cpp index 0ed1079f35..1cb7c336a1 100644 --- a/Source/Utils/ERF_InteriorGhostCells.cpp +++ b/Source/Utils/ERF_InteriorGhostCells.cpp @@ -513,7 +513,7 @@ realbdy_compute_interior_ghost_rhs (const Real& time, #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(S_cur_data[ivarT],TilingIfNotGPU()); mfi.isValid(); ++mfi) { + for (MFIter mfi(S_cur_data[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box tbx = mfi.nodaltilebox(0); Box tby = mfi.nodaltilebox(1); From 3db3b5b0f387fe25b18c48b6bcc19427ea01b5f9 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 14:52:14 -0700 Subject: [PATCH 15/27] Fix g test. --- Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerCommon.H | 3 +-- .../Kessler/ERF_GTestKesslerPhysicalProperties.cpp | 3 +-- Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjMultiFab.cpp | 4 +--- 3 files changed, 3 insertions(+), 7 deletions(-) diff --git a/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerCommon.H b/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerCommon.H index 37b03f0575..b579f2d2c2 100644 --- a/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerCommon.H +++ b/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerCommon.H @@ -687,7 +687,6 @@ inline void run_kessler_public_flow (Kessler& kessler, { run_and_sync([&] { kessler.Define(sc); - kessler.Set_RealWidth(0); kessler.Set_dzmin(geom.CellSize(2)); kessler.Init(cons, cons.boxArray(), geom, dt, z_phys_nd, detJ_cc); kessler.Copy_State_to_Micro(cons); @@ -899,4 +898,4 @@ inline std::string case_label (const KernelCase& test_case) } // namespace kessler_test -#endif \ No newline at end of file +#endif diff --git a/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerPhysicalProperties.cpp b/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerPhysicalProperties.cpp index 5401e10f00..02ba30283d 100644 --- a/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerPhysicalProperties.cpp +++ b/Tests/Unit/Microphysics/Kessler/ERF_GTestKesslerPhysicalProperties.cpp @@ -23,7 +23,6 @@ void initialize_kessler (Kessler& kessler, std::unique_ptr& detJ_cc) { kessler.Define(sc); - kessler.Set_RealWidth(0); kessler.Set_dzmin(geom.CellSize(2)); kessler.Init(cons, cons.boxArray(), geom, kDefaultDt, z_phys_nd, detJ_cc); } @@ -579,4 +578,4 @@ TEST(KesslerPhysicalProperties, KesslerFullRain_DoCondFalseStillRunsRainPath) EXPECT_NEAR(qv_after, qv_before, property_accumulation_tol(4, qv_before)); EXPECT_LT(qc_after, qc_before); EXPECT_GT(rain_accum_after, amrex::Real(0.0)); -} \ No newline at end of file +} diff --git a/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjMultiFab.cpp b/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjMultiFab.cpp index 109256a64c..c9c1c3d4d6 100644 --- a/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjMultiFab.cpp +++ b/Tests/Unit/Microphysics/SatAdj/ERF_GTestSatAdjMultiFab.cpp @@ -154,7 +154,6 @@ TEST(SatAdjMultiFab, PublicFlowPreservesSatAdjInvariantsPortable) SatAdj satadj; SolverChoice sc = make_solver_choice(false); satadj.Define(sc); - satadj.Set_RealWidth(0); run_public_flow(satadj, sc, geom, cons); compute_satadj_invariant_errors(cons_initial, cons, err); @@ -203,7 +202,6 @@ TEST(SatAdjMultiFab, PublicFlowConservesWaterAndLatentEnergyPortable) SatAdj satadj; SolverChoice sc = make_solver_choice(false); satadj.Define(sc); - satadj.Set_RealWidth(0); run_public_flow(satadj, sc, geom, cons); compute_satadj_conservation_errors(cons_initial, cons, err); @@ -309,4 +307,4 @@ TEST(SatAdjKernel, AdjustCellMatchesHostReferencePortable) EXPECT_NEAR(qc_out[idx], host_reference[idx].qc, scaled_tol(host_reference[idx].qc, kernel_tol_factor)); } -} \ No newline at end of file +} From f4cef0bcd6b2ce031182ff0de8db361411f38e8a Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 17:51:06 -0700 Subject: [PATCH 16/27] don't over run the last bdy file. --- .../ERF_BoundaryConditionsRealbdy.cpp | 42 +++++++++---------- Source/Utils/ERF_InteriorGhostCells.cpp | 11 ++--- 2 files changed, 23 insertions(+), 30 deletions(-) diff --git a/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp b/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp index 579f8967b1..9fba8f94c9 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp @@ -29,21 +29,20 @@ ERF::fill_from_realbdy (const Vector& mfs, // Time interpolation Real dT = bdy_time_interval; - // Note this is because we define "time" to be time since start_bdy_time - Real time_since_start = time; - - int n_time = static_cast( time_since_start / dT); - Real alpha = (time_since_start - n_time * dT) / dT; - AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); - Real oma = 1.0 - alpha; - + int n_time = static_cast( (time-start_bdy_time) / dT); int n_time_p1 = n_time + 1; - if ((time == stop_time-start_time) && (alpha==0)) { - // stop time coincides with final bdy snapshot -- don't try to read in - // another snapshot + Real alpha = ((time-start_bdy_time) - n_time * dT) / dT; + + // Do not over run the last bdy file + if (time >= final_bdy_time) { + n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); n_time_p1 = n_time; + alpha = zero; } + AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one); + Real oma = one - alpha; + // Flags for read vars and index mapping Vector cons_read = {0, 1, 0, 0, 1, 0, 0, @@ -248,21 +247,20 @@ ERF::fill_from_realbdy_upwind (const Vector& mfs, // Time interpolation Real dT = bdy_time_interval; - // Note this is because we define "time" to be time since start_bdy_time - Real time_since_start = time; - - int n_time = static_cast( time_since_start / dT); - Real alpha = (time_since_start - n_time * dT) / dT; - AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0); - Real oma = 1.0 - alpha; - + int n_time = static_cast( (time-start_bdy_time) / dT); int n_time_p1 = n_time + 1; - if ((time == stop_time-start_time) && (alpha==0)) { - // stop time coincides with final bdy snapshot -- don't try to read in - // another snapshot + Real alpha = ((time-start_bdy_time) - n_time * dT) / dT; + + // Do not over run the last bdy file + if (time >= final_bdy_time) { + n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); n_time_p1 = n_time; + alpha = zero; } + AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one); + Real oma = one - alpha; + // Map loop index to variable index // NOTE: Loop backwards and do v/u first for WfromOmega Vector var_idx_map = {Vars::cons, Vars::zvel, Vars::xvel, Vars::yvel}; diff --git a/Source/Utils/ERF_InteriorGhostCells.cpp b/Source/Utils/ERF_InteriorGhostCells.cpp index 1cb7c336a1..5b0a93d7e2 100644 --- a/Source/Utils/ERF_InteriorGhostCells.cpp +++ b/Source/Utils/ERF_InteriorGhostCells.cpp @@ -206,19 +206,14 @@ realbdy_compute_interior_ghost_rhs (const Real& time, // Do not over run the last bdy file if (time >= final_bdy_time) { - n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); - n_time_p1 = n_time; - alpha = zero; + n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); + n_time_p1 = n_time; + alpha = zero; } AMREX_ALWAYS_ASSERT( alpha >= zero && alpha <= one); Real oma = one - alpha; - /* - // UNIT TEST DEBUG - oma = one; alpha = zero; - */ - // Temporary FABs for storage (owned/filled on all ranks) FArrayBox U_xlo, U_xhi, U_ylo, U_yhi; FArrayBox V_xlo, V_xhi, V_ylo, V_yhi; From 79ed23f26d338bbac15b3eed58b73022d59d79e4 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 26 May 2026 18:24:18 -0700 Subject: [PATCH 17/27] don't over run last bdy file. --- .../ERF_BoundaryConditionsRealbdy.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp b/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp index 9fba8f94c9..f2bd97f3a5 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditionsRealbdy.cpp @@ -29,12 +29,14 @@ ERF::fill_from_realbdy (const Vector& mfs, // Time interpolation Real dT = bdy_time_interval; - int n_time = static_cast( (time-start_bdy_time) / dT); + Real time_tot = time + start_time; + Real time_since_start_bdy = time_tot - start_bdy_time; + int n_time = static_cast( time_since_start_bdy / dT); int n_time_p1 = n_time + 1; - Real alpha = ((time-start_bdy_time) - n_time * dT) / dT; + Real alpha = (time_since_start_bdy - n_time * dT) / dT; // Do not over run the last bdy file - if (time >= final_bdy_time) { + if (time_tot >= final_bdy_time) { n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); n_time_p1 = n_time; alpha = zero; @@ -247,12 +249,14 @@ ERF::fill_from_realbdy_upwind (const Vector& mfs, // Time interpolation Real dT = bdy_time_interval; - int n_time = static_cast( (time-start_bdy_time) / dT); + Real time_tot = time + start_time; + Real time_since_start_bdy = time_tot - start_bdy_time; + int n_time = static_cast( time_since_start_bdy / dT); int n_time_p1 = n_time + 1; - Real alpha = ((time-start_bdy_time) - n_time * dT) / dT; + Real alpha = (time_since_start_bdy - n_time * dT) / dT; // Do not over run the last bdy file - if (time >= final_bdy_time) { + if (time_tot >= final_bdy_time) { n_time = static_cast( (final_bdy_time - start_bdy_time)/ dT); n_time_p1 = n_time; alpha = zero; From b874d05c2e153cd5343739ed77494884ee5dd693 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Wed, 27 May 2026 07:48:40 -0700 Subject: [PATCH 18/27] fix develop merge. --- Source/Microphysics/Morrison/ERF_Morrison.H | 3 --- 1 file changed, 3 deletions(-) diff --git a/Source/Microphysics/Morrison/ERF_Morrison.H b/Source/Microphysics/Morrison/ERF_Morrison.H index 98e0cbf64e..f54d651890 100644 --- a/Source/Microphysics/Morrison/ERF_Morrison.H +++ b/Source/Microphysics/Morrison/ERF_Morrison.H @@ -134,9 +134,6 @@ public: int Qstate_Moist_Size () override { return Morrison::n_qstate_moist_size; } - void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } - void Qmoist_Restart_Vars (const SolverChoice& /*a_sc*/, std::vector& a_idx, From d677778f75745107f2a85a0eeb68ffd825fcfeed Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Wed, 27 May 2026 08:43:50 -0700 Subject: [PATCH 19/27] explicitly zero fast RHS at walls with real bcs for safety. --- Source/TimeIntegration/ERF_Substep_MT.cpp | 1 + Source/TimeIntegration/ERF_Substep_NS.cpp | 16 ++++++++++++++-- Source/TimeIntegration/ERF_Substep_T.cpp | 14 ++++++++++++-- Source/TimeIntegration/ERF_TI_fast_headers.H | 6 +++--- Source/TimeIntegration/ERF_TI_substep_fun.H | 6 +++--- 5 files changed, 33 insertions(+), 10 deletions(-) diff --git a/Source/TimeIntegration/ERF_Substep_MT.cpp b/Source/TimeIntegration/ERF_Substep_MT.cpp index 940cb91316..1d7db7ffa7 100644 --- a/Source/TimeIntegration/ERF_Substep_MT.cpp +++ b/Source/TimeIntegration/ERF_Substep_MT.cpp @@ -83,6 +83,7 @@ void erf_substep_MT (int step, int /*nrk*/, YAFluxRegister* fr_as_fine, bool l_use_moisture, bool l_reflux, + bool /*l_real_bc*/, const Real* sinesq_stag_d, const Real l_damp_coef) { diff --git a/Source/TimeIntegration/ERF_Substep_NS.cpp b/Source/TimeIntegration/ERF_Substep_NS.cpp index 211bc8784f..6f9b24222c 100644 --- a/Source/TimeIntegration/ERF_Substep_NS.cpp +++ b/Source/TimeIntegration/ERF_Substep_NS.cpp @@ -67,6 +67,7 @@ void erf_substep_NS (int step, int nrk, YAFluxRegister* fr_as_fine, bool l_use_moisture, bool l_reflux, + bool l_real_bc, const amrex::Real* sinesq_stag_d, const Real l_damp_coef) { @@ -76,6 +77,15 @@ void erf_substep_NS (int step, int nrk, BL_PROFILE_REGION("erf_substep_S()"); + const Box& domain = geom.Domain(); + auto const domlo = lbound(domain); + auto const domhi = ubound(domain); + + int ilo = domlo.x; + int ihi = domhi.x + 1; + int jlo = domlo.y; + int jhi = domhi.y + 1; + Real beta_1 = myhalf * (one - beta_s); // multiplies explicit terms Real beta_2 = myhalf * (one + beta_s); // multiplies implicit terms @@ -237,7 +247,8 @@ void erf_substep_NS (int step, int nrk, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Add (negative) gradient of (rho theta) multiplied by lagged "pi" - Real gpx = (theta_extrap(i,j,k) - theta_extrap(i-1,j,k))*dxi; + Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ? Real(0.) : + (theta_extrap(i,j,k) - theta_extrap(i-1,j,k))*dxi; gpx *= mf_ux(i,j,0); Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) : zero; @@ -256,7 +267,8 @@ void erf_substep_NS (int step, int nrk, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Add (negative) gradient of (rho theta) multiplied by lagged "pi" - Real gpy = (theta_extrap(i,j,k) - theta_extrap(i,j-1,k))*dyi; + Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ? Real(0.) : + (theta_extrap(i,j,k) - theta_extrap(i,j-1,k))*dyi; gpy *= mf_vy(i,j,0); Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) : zero; diff --git a/Source/TimeIntegration/ERF_Substep_T.cpp b/Source/TimeIntegration/ERF_Substep_T.cpp index fd6b596766..94dd5e2fc4 100644 --- a/Source/TimeIntegration/ERF_Substep_T.cpp +++ b/Source/TimeIntegration/ERF_Substep_T.cpp @@ -69,6 +69,7 @@ void erf_substep_T (int step, int /*nrk*/, YAFluxRegister* fr_as_fine, bool l_use_moisture, bool l_reflux, + bool l_real_bc, const Real* sinesq_stag_d, const Real l_damp_coef) { @@ -78,6 +79,11 @@ void erf_substep_T (int step, int /*nrk*/, auto const domlo = lbound(domain); auto const domhi = ubound(domain); + int ilo = domlo.x; + int ihi = domhi.x + 1; + int jlo = domlo.y; + int jhi = domhi.y + 1; + Real beta_1 = myhalf * (one - beta_s); // multiplies explicit terms Real beta_2 = myhalf * (one + beta_s); // multiplies implicit terms @@ -291,7 +297,9 @@ void erf_substep_T (int step, int /*nrk*/, - theta_extrap(i-1,j,k ) - theta_extrap(i,j,k ) ) : fourth * dzi * ( theta_extrap(i-1,j,k+1) + theta_extrap(i,j,k+1) - theta_extrap(i-1,j,k-1) - theta_extrap(i,j,k-1) ); - Real gpx = gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface; + Real gpx = (l_real_bc && (level==0) && (i==ilo || i==ihi)) ? Real(0.) : + gp_xi - (met_h_xi / met_h_zeta) * gp_zeta_on_iface; + gpx *= mf_ux(i,j,0); Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i-1,j,k)) : zero; @@ -323,7 +331,9 @@ void erf_substep_T (int step, int /*nrk*/, - theta_extrap(i,j,k ) - theta_extrap(i,j-1,k ) ) : fourth * dzi * ( theta_extrap(i,j,k+1) + theta_extrap(i,j-1,k+1) - theta_extrap(i,j,k-1) - theta_extrap(i,j-1,k-1) ); - Real gpy = gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface; + Real gpy = (l_real_bc && (level==0) && (j==jlo || j==jhi)) ? Real(0.) : + gp_eta - (met_h_eta / met_h_zeta) * gp_zeta_on_jface; + gpy *= mf_vy(i,j,0); Real q = (l_use_moisture) ? myhalf * (qt_arr(i,j,k) + qt_arr(i,j-1,k)) : zero; diff --git a/Source/TimeIntegration/ERF_TI_fast_headers.H b/Source/TimeIntegration/ERF_TI_fast_headers.H index bda450af56..5b0776e183 100644 --- a/Source/TimeIntegration/ERF_TI_fast_headers.H +++ b/Source/TimeIntegration/ERF_TI_fast_headers.H @@ -41,7 +41,7 @@ void erf_substep_NS (int step, int nrk, int level, int finest_level, amrex::Vector>& mapfac, amrex::YAFluxRegister* fr_as_crse, amrex::YAFluxRegister* fr_as_fine, - bool l_use_moisture, bool l_reflux, + bool l_use_moisture, bool l_reflux, bool l_real_bc, const amrex::Real* sinesq_at_lev_d, const amrex::Real l_damp_coef); @@ -75,7 +75,7 @@ void erf_substep_T (int step, int nrk, int level, int finest_level, amrex::Vector>& mapfac, amrex::YAFluxRegister* fr_as_crse, amrex::YAFluxRegister* fr_as_fine, - bool l_use_moisture, bool l_reflux, + bool l_use_moisture, bool l_reflux, bool l_real_bc, const amrex::Real* sinesq_at_lev_d, const amrex::Real l_damp_coef); @@ -116,7 +116,7 @@ void erf_substep_MT (int step, int nrk, int level, int finest_level, amrex::Vector>& mapfac, amrex::YAFluxRegister* fr_as_crse, amrex::YAFluxRegister* fr_as_fine, - bool l_use_moisture, bool l_reflux, + bool l_use_moisture, bool l_reflux, bool l_real_bc, const amrex::Real* sinesq_at_lev_d, const amrex::Real l_damp_coef); diff --git a/Source/TimeIntegration/ERF_TI_substep_fun.H b/Source/TimeIntegration/ERF_TI_substep_fun.H index 4496dc6785..098a4a4ba4 100644 --- a/Source/TimeIntegration/ERF_TI_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_substep_fun.H @@ -140,7 +140,7 @@ auto acoustic_substepping_fun = [&](int fast_step, int n_sub, int nrk, z_phys_nd[level], z_phys_nd_new[level], z_phys_nd_src[level], detJ_cc[level], detJ_cc_new[level], detJ_cc_src[level], dtau, beta_s, inv_fac, mapfac[level], - fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, solverChoice.use_real_bcs, (l_rayleigh_implicit) ? d_sinesq_ptrs[level].data() : nullptr, l_damp_coef); } else if (solverChoice.mesh_type == MeshType::VariableDz) { @@ -152,7 +152,7 @@ auto acoustic_substepping_fun = [&](int fast_step, int n_sub, int nrk, fine_geom, solverChoice.gravity, z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac, mapfac[level], - fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, solverChoice.use_real_bcs, (l_rayleigh_implicit) ? d_sinesq_ptrs[level].data() : nullptr, l_damp_coef); } else { @@ -163,7 +163,7 @@ auto acoustic_substepping_fun = [&](int fast_step, int n_sub, int nrk, cc_src, xmom_src, ymom_src, zmom_src, fine_geom, solverChoice.gravity, stretched_dz_d[level], dtau, beta_s, inv_fac, mapfac[level], - fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, + fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, solverChoice.use_real_bcs, (l_rayleigh_implicit) ? d_sinesq_ptrs[level].data() : nullptr, l_damp_coef); } From 1bc08f8814d74d49250dd06709a0dd1a47486cb3 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Wed, 27 May 2026 09:35:40 -0700 Subject: [PATCH 20/27] WIP debugging wrfinput on HIP crash seemingly due to moisture treatment at the boundaries --- Source/IO/ERF_ReadFromERFBdy.cpp | 33 +++++-- Source/IO/ERF_ReadFromWRFBdy.cpp | 14 ++- Source/IO/ERF_WriteERFBdy.cpp | 2 +- Source/Initialization/ERF_InitFromMetgrid.cpp | 4 +- .../Initialization/ERF_InitFromWRFInput.cpp | 30 +++--- .../Microphysics/ERF_EulerianMicrophysics.H | 2 +- .../Microphysics/ERF_LagrangianMicrophysics.H | 4 +- Source/Microphysics/Kessler/ERF_Kessler.H | 4 +- Source/Microphysics/Kessler/ERF_Kessler.cpp | 5 + Source/Microphysics/Morrison/ERF_Morrison.H | 2 +- Source/Microphysics/Null/ERF_NullMoist.H | 2 +- Source/Microphysics/SAM/ERF_SAM.H | 2 +- Source/Microphysics/SatAdj/ERF_SatAdj.H | 2 +- Source/Microphysics/WSM6/ERF_WSM6.H | 2 +- Source/SourceTerms/ERF_MoistSetRhs.cpp | 39 +++++++- Source/SourceTerms/ERF_SrcHeaders.H | 3 +- .../ERF_AdvanceMicrophysics.cpp | 98 +++++++++++++++++++ Source/TimeIntegration/ERF_SlowRhsPost.cpp | 2 +- 18 files changed, 207 insertions(+), 43 deletions(-) diff --git a/Source/IO/ERF_ReadFromERFBdy.cpp b/Source/IO/ERF_ReadFromERFBdy.cpp index c2b930a8d9..be8aea4464 100644 --- a/Source/IO/ERF_ReadFromERFBdy.cpp +++ b/Source/IO/ERF_ReadFromERFBdy.cpp @@ -35,7 +35,7 @@ read_times_from_erfbdy(const std::string& bdy_file_name, HeaderFile >> nvars; HeaderFile >> real_width; - // Read time values. + // Read boundary times. bdy_times.resize(ntimes); for (int i = 0; i < ntimes; ++i) { HeaderFile >> bdy_times[i]; @@ -48,14 +48,14 @@ read_times_from_erfbdy(const std::string& bdy_file_name, HeaderFile.close(); - // Ensure the file hodls at least two times. + // Ensure the file holds at least two times. AMREX_ALWAYS_ASSERT(ntimes >= 2); - // Set start and final times. + // Set the first and last boundary times. start_bdy_time = bdy_times[0]; final_bdy_time = bdy_times[ntimes - 1]; - // Calculate time interval. + // Calculate the interval between boundary times. Real bdy_time_interval = bdy_times[1] - bdy_times[0]; return bdy_time_interval; @@ -73,6 +73,11 @@ read_from_erfbdy(int itime, Print() << "Reading ERF boundary data for time index " << itime << std::endl; std::string time_dir = bdy_file_name + "/Time_" + Concatenate("", itime, 6); + Arena* Arena_Used = The_Arena(); +#ifdef AMREX_USE_GPU + Arena_Used = The_Pinned_Arena(); +#endif + if (bdy_data_xlo[itime].empty()) { bdy_data_xlo[itime].resize(nvars); bdy_data_xhi[itime].resize(nvars); @@ -85,28 +90,40 @@ read_from_erfbdy(int itime, { // X-low boundary std::string filename = time_dir + "/BdyData_xlo_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_xlo[itime][ivar].readFrom(ifs); + FArrayBox tmp_fab; + tmp_fab.readFrom(ifs); + bdy_data_xlo[itime][ivar].resize(tmp_fab.box(), tmp_fab.nComp(), Arena_Used); + bdy_data_xlo[itime][ivar].template copy(tmp_fab, 0, 0, tmp_fab.nComp()); ifs.close(); } { // X-high boundary std::string filename = time_dir + "/BdyData_xhi_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_xhi[itime][ivar].readFrom(ifs); + FArrayBox tmp_fab; + tmp_fab.readFrom(ifs); + bdy_data_xhi[itime][ivar].resize(tmp_fab.box(), tmp_fab.nComp(), Arena_Used); + bdy_data_xhi[itime][ivar].template copy(tmp_fab, 0, 0, tmp_fab.nComp()); ifs.close(); } { // Y-low boundary std::string filename = time_dir + "/BdyData_ylo_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_ylo[itime][ivar].readFrom(ifs); + FArrayBox tmp_fab; + tmp_fab.readFrom(ifs); + bdy_data_ylo[itime][ivar].resize(tmp_fab.box(), tmp_fab.nComp(), Arena_Used); + bdy_data_ylo[itime][ivar].template copy(tmp_fab, 0, 0, tmp_fab.nComp()); ifs.close(); } { // Y-high boundary std::string filename = time_dir + "/BdyData_yhi_var" + std::to_string(ivar); std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary); - bdy_data_yhi[itime][ivar].readFrom(ifs); + FArrayBox tmp_fab; + tmp_fab.readFrom(ifs); + bdy_data_yhi[itime][ivar].resize(tmp_fab.box(), tmp_fab.nComp(), Arena_Used); + bdy_data_yhi[itime][ivar].template copy(tmp_fab, 0, 0, tmp_fab.nComp()); ifs.close(); } } diff --git a/Source/IO/ERF_ReadFromWRFBdy.cpp b/Source/IO/ERF_ReadFromWRFBdy.cpp index 1cc29fc386..d98fae4c32 100644 --- a/Source/IO/ERF_ReadFromWRFBdy.cpp +++ b/Source/IO/ERF_ReadFromWRFBdy.cpp @@ -124,7 +124,7 @@ read_from_wrfbdy (const int itime, const std::string& nc_bdy_file, const Box& do // ****************************************************************** // Read the netcdf file and fill these FABs // NOTE: the order and number of these must match the WRFBdyVars enum! - // WRFBdyVars: U, V, R, T, QV, MU, PC + // WRFBdyVars: U, V, T, QV, MU, PC // // These fields are at myhalf levels (unstaggered) // ****************************************************************** @@ -143,15 +143,18 @@ read_from_wrfbdy (const int itime, const std::string& nc_bdy_file, const Box& do Vector tslice(nc_var_names.size()); int width; // size of bdy_width from wrfbdy + Vector success(nc_var_names.size()); if (ParallelDescriptor::IOProcessor()) { - Vector success(nc_var_names.size()); ReadTimeSliceFromNetCDFFile(nc_bdy_file, itime, nc_var_names, tslice, success); - for (auto &istat:success) { - AMREX_ALWAYS_ASSERT(istat==1); + // Check that all required variables were read successfully + for (int i = 0; i < success.size(); ++i) { + if (success[i] != 1) { + Abort("Missing required wrfbdy variable " + nc_var_names[i]); + } } // Width of the boundary region @@ -328,6 +331,9 @@ read_from_wrfbdy (const int itime, const std::string& nc_bdy_file, const Box& do // Now fill the data if (ParallelDescriptor::IOProcessor()) { + // Skip filling data for variables that weren't successfully read from the file + if (success[iv] != 1) continue; + // Print() << "SHAPE0 " << tslice[iv].get_vshape()[0] << std::endl; // Print() << "SHAPE1 " << tslice[iv].get_vshape()[1] << std::endl; // Print() << "SHAPE2 " << tslice[iv].get_vshape()[2] << std::endl; diff --git a/Source/IO/ERF_WriteERFBdy.cpp b/Source/IO/ERF_WriteERFBdy.cpp index 0f632ffdf4..5b90080e26 100644 --- a/Source/IO/ERF_WriteERFBdy.cpp +++ b/Source/IO/ERF_WriteERFBdy.cpp @@ -81,7 +81,7 @@ void WriteERFBdyTimeSlice(const std::string& bdy_file_name, // Barrier to ensure directory exists. ParallelDescriptor::Barrier(); - // Write each variable for each boundary direction using FArrayBox::writeOn. + // Write each variable for each boundary direction. if (ParallelDescriptor::IOProcessor()) { for (int ivar = 0; ivar < nvars; ++ivar) diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 72fa105d91..ad64ef6d7c 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -67,14 +67,14 @@ ERF::init_from_metgrid (int lev) Print() << "Init with met_em without moisture model." << std::endl; } - // Check for erfbdy file. + // Check for an erfbdy file. if (lev == 0) { std::string erfbdy_header = erfbdy_file + "/Header"; use_erfbdy = FileSystem::Exists(erfbdy_header); } if (use_erfbdy || write_erfbdy) nvars_erfbdy = MetGridBdyVars::NumTypes; - // If erfbdy file exists, load boundary data and skip met_em processing. + // If the erfbdy file exists, load boundary data and skip met_em processing. if (lev == 0 && use_erfbdy) { Print() << "Loading boundary data from erfbdy file: " << erfbdy_file << std::endl; diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index 4a1ed99f81..63a7956e2b 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -158,22 +158,20 @@ ERF::init_from_wrfinput (int lev, NC_names.push_back("XLONG_U"); // 22 if (use_moist) { NC_names.push_back("QVAPOR"); // 23 - NC_names.push_back("QCLOUD"); // 24 - NC_names.push_back("QRAIN"); // 25 } - NC_names.push_back("IVGTYP"); // 26 - NC_names.push_back("ISLTYP"); // 27 + NC_names.push_back("IVGTYP"); // 24 + NC_names.push_back("ISLTYP"); // 25 if (use_lsm) { - NC_names.push_back("TSLB"); // 28 - NC_names.push_back("SMOIS"); // 29 - NC_names.push_back("SH2O"); // 30 - NC_names.push_back("LAI"); // 31 - NC_names.push_back("ZS"); // 32 - NC_names.push_back("DZS"); // 33 - NC_names.push_back("VEGFRA"); // 34 - NC_names.push_back("TMN"); // 35 - NC_names.push_back("SHDMIN"); // 36 - NC_names.push_back("SHDMAX"); // 37 + NC_names.push_back("TSLB"); // 26 + NC_names.push_back("SMOIS"); // 27 + NC_names.push_back("SH2O"); // 28 + NC_names.push_back("LAI"); // 29 + NC_names.push_back("ZS"); // 30 + NC_names.push_back("DZS"); // 31 + NC_names.push_back("VEGFRA"); // 32 + NC_names.push_back("TMN"); // 33 + NC_names.push_back("SHDMIN"); // 34 + NC_names.push_back("SHDMAX"); // 35 // --- debugging --- // print LSM varname->WRF input name map @@ -246,8 +244,8 @@ ERF::init_from_wrfinput (int lev, auto& var_fab_from_file = NC_fab_var[idx][ivar]; bool has_fallback_behavior = (var_name == "U") || (var_name == "V") || (var_name == "W") || - (var_name == "THM") || (var_name == "QVAPOR") || (var_name == "QCLOUD") || - (var_name == "QRAIN") || (var_name == "PH") || (var_name == "PHB"); + (var_name == "THM") || (var_name == "QVAPOR") || + (var_name == "PH") || (var_name == "PHB"); if (!success && !has_fallback_behavior) { amrex::Abort(std::string("ERF::init_from_wrfinput: failed to read required variable " + var_name).c_str()); } diff --git a/Source/Microphysics/ERF_EulerianMicrophysics.H b/Source/Microphysics/ERF_EulerianMicrophysics.H index c6d364cd4d..4614045251 100644 --- a/Source/Microphysics/ERF_EulerianMicrophysics.H +++ b/Source/Microphysics/ERF_EulerianMicrophysics.H @@ -164,7 +164,7 @@ public: virtual void Set_RealWidth (const int a_lev, const int real_width) const override { - m_moist_model[a_lev]->Set_RealWidth(real_width); + m_moist_model[a_lev]->Set_RealWidth(a_lev, real_width); } protected: diff --git a/Source/Microphysics/ERF_LagrangianMicrophysics.H b/Source/Microphysics/ERF_LagrangianMicrophysics.H index 71147b4b64..9f059c62a4 100644 --- a/Source/Microphysics/ERF_LagrangianMicrophysics.H +++ b/Source/Microphysics/ERF_LagrangianMicrophysics.H @@ -199,10 +199,10 @@ public: } /*! \brief Populates real_width in micro model for box limiting */ - virtual void Set_RealWidth (const int /*a_lev*/, + virtual void Set_RealWidth (const int a_lev, const int real_width) const override { - m_moist_model->Set_RealWidth(real_width); + m_moist_model->Set_RealWidth(a_lev, real_width); } protected: diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index 86703721f7..5e2cc5176d 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -121,7 +121,9 @@ public: Qstate_Moist_Size () override { return Kessler::n_qstate_moist_size; } void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } + Set_RealWidth (const int a_lev, const int real_width) override { + m_real_width = real_width; + } void Qmoist_Restart_Vars (const SolverChoice& /*a_sc*/, diff --git a/Source/Microphysics/Kessler/ERF_Kessler.cpp b/Source/Microphysics/Kessler/ERF_Kessler.cpp index 72d64bbe62..c200129356 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.cpp +++ b/Source/Microphysics/Kessler/ERF_Kessler.cpp @@ -1,5 +1,6 @@ #include #include +#include #include "ERF_Kessler.H" #include "ERF_DataStruct.H" @@ -25,6 +26,7 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) auto ba = tabs->boxArray(); auto dm = tabs->DistributionMap(); fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells + fz.setVal(zero); // Initialize to zero for boundary cells. Real dtn = dt; Real coef = dtn/m_dzmin; @@ -136,6 +138,9 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) }); } + // Initialize fz to zero (including ghost cells) before computing terminal velocity + fz.setVal(Real(0.0)); + // Precompute terminal velocity for substepping for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); diff --git a/Source/Microphysics/Morrison/ERF_Morrison.H b/Source/Microphysics/Morrison/ERF_Morrison.H index aacdf6b8a6..3cfaf80ffd 100644 --- a/Source/Microphysics/Morrison/ERF_Morrison.H +++ b/Source/Microphysics/Morrison/ERF_Morrison.H @@ -135,7 +135,7 @@ public: Qstate_Moist_Size () override { return Morrison::n_qstate_moist_size; } void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } + Set_RealWidth (const int /*a_lev*/, const int real_width) override { m_real_width = real_width; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Microphysics/Null/ERF_NullMoist.H b/Source/Microphysics/Null/ERF_NullMoist.H index 0f5e24ee99..e5816429fa 100644 --- a/Source/Microphysics/Null/ERF_NullMoist.H +++ b/Source/Microphysics/Null/ERF_NullMoist.H @@ -91,7 +91,7 @@ public: virtual void - Set_RealWidth (const int /*real_width*/) { } + Set_RealWidth (const int /*a_lev*/, const int /*real_width*/) { } private: diff --git a/Source/Microphysics/SAM/ERF_SAM.H b/Source/Microphysics/SAM/ERF_SAM.H index a2fdf314fb..41f59f8e17 100644 --- a/Source/Microphysics/SAM/ERF_SAM.H +++ b/Source/Microphysics/SAM/ERF_SAM.H @@ -153,7 +153,7 @@ public: Qstate_Moist_Size () override { return SAM::n_qstate_moist_size; } void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } + Set_RealWidth (const int /*a_lev*/, const int real_width) override { m_real_width = real_width; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.H b/Source/Microphysics/SatAdj/ERF_SatAdj.H index 17fb146914..fff1952d86 100644 --- a/Source/Microphysics/SatAdj/ERF_SatAdj.H +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.H @@ -126,7 +126,7 @@ public: } void - Set_RealWidth (const int real_width) override { m_real_width = real_width; } + Set_RealWidth (const int /*a_lev*/, const int real_width) override { m_real_width = real_width; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Microphysics/WSM6/ERF_WSM6.H b/Source/Microphysics/WSM6/ERF_WSM6.H index 26b0d9ce6d..32d11a3cba 100644 --- a/Source/Microphysics/WSM6/ERF_WSM6.H +++ b/Source/Microphysics/WSM6/ERF_WSM6.H @@ -82,7 +82,7 @@ public: int Qmoist_Size() override { return m_qmoist_size; } int Qstate_Moist_Size() override { return n_qstate_moist_size; } - void Set_RealWidth(const int real_width) override { m_real_width = real_width; } + void Set_RealWidth(const int /*a_lev*/, const int real_width) override { m_real_width = real_width; } void Qmoist_Restart_Vars(const SolverChoice&, std::vector& a_idx, diff --git a/Source/SourceTerms/ERF_MoistSetRhs.cpp b/Source/SourceTerms/ERF_MoistSetRhs.cpp index c967254c68..53990979b6 100644 --- a/Source/SourceTerms/ERF_MoistSetRhs.cpp +++ b/Source/SourceTerms/ERF_MoistSetRhs.cpp @@ -2,6 +2,7 @@ #include #include +#include using namespace amrex; @@ -28,7 +29,8 @@ moist_set_rhs (const Geometry& geom, Vector>& bdy_data_xhi, Vector>& bdy_data_ylo, Vector>& bdy_data_yhi, - std::unique_ptr& m_r2d) + std::unique_ptr& m_r2d, + int n_qstate) { // HACK HACK HACK // Get bndry data @@ -254,5 +256,40 @@ moist_set_rhs (const Geometry& geom, }); exit(0); */ + + // Zero out hydrometeors in boundary regions based on active microphysics scheme + ParallelFor(tbx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; + }); + ParallelFor(tbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; + }); + ParallelFor(tbx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; + }); + ParallelFor(tbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; + }); + } // moist_set_rhs #endif diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H index bbe5b29c0c..57adfa5bbd 100644 --- a/Source/SourceTerms/ERF_SrcHeaders.H +++ b/Source/SourceTerms/ERF_SrcHeaders.H @@ -140,7 +140,8 @@ moist_set_rhs (const amrex::Geometry& geom, amrex::Vector>& bdy_data_xhi, amrex::Vector>& bdy_data_ylo, amrex::Vector>& bdy_data_yhi, - std::unique_ptr& m_r2d); + std::unique_ptr& m_r2d, + int n_qstate); #endif void ApplySpongeZoneBCsForCC (const SpongeChoice& spongeChoice, diff --git a/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp b/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp index 9ddfbc02bd..6ad8f5f418 100644 --- a/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp +++ b/Source/TimeIntegration/ERF_AdvanceMicrophysics.cpp @@ -1,4 +1,5 @@ #include +#include using namespace amrex; @@ -11,8 +12,105 @@ void ERF::advance_microphysics (int lev, if (solverChoice.moisture_type != MoistureType::None) { micro->Set_RealWidth(lev, real_width); cons.FillBoundary(geom[lev].periodicity()); + + // Get number of moisture species for the microphysics scheme. + int n_qstate = micro->Get_Qstate_Moist_Size(); + + // Zero condensed moisture in boundary cells when using real BCs. + // Only zero if we have hydrometeors (n_qstate > 1 means QV + at least one hydrometeor) + if (real_width > 0 && lev == 0 && n_qstate > 1) { + amrex::Print() << "DEBUG: Zeroing boundary hydrometeors at time=" << time + << " n_qstate=" << n_qstate << " real_width=" << real_width << "\n"; + const auto& domain = geom[lev].Domain(); + int i_lo = domain.smallEnd(0); + int i_hi = domain.bigEnd(0); + int j_lo = domain.smallEnd(1); + int j_hi = domain.bigEnd(1); + + for (MFIter mfi(cons, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + auto cons_arr = cons.array(mfi); + Box tbx = mfi.tilebox(); + + // Zero boundary cells on domain faces + if (tbx.smallEnd(0) == i_lo) { + Box bx_xlo(tbx); bx_xlo.setBig(0, i_lo + real_width - 1); + amrex::Print() << "DEBUG: Zeroing xlo boundary, box=" << bx_xlo << "\n"; + ParallelFor(bx_xlo, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Zero only hydrometeors used by the current microphysics scheme + if (n_qstate >= 2) cons_arr(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cons_arr(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cons_arr(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cons_arr(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cons_arr(i,j,k,RhoQ6_comp) = zero; + }); + } + if (tbx.bigEnd(0) == i_hi) { + Box bx_xhi(tbx); bx_xhi.setSmall(0, i_hi - real_width + 1); + amrex::Print() << "DEBUG: Zeroing xhi boundary, box=" << bx_xhi << "\n"; + ParallelFor(bx_xhi, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (n_qstate >= 2) cons_arr(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cons_arr(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cons_arr(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cons_arr(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cons_arr(i,j,k,RhoQ6_comp) = zero; + }); + } + if (tbx.smallEnd(1) == j_lo) { + Box bx_ylo(tbx); bx_ylo.setBig(1, j_lo + real_width - 1); + amrex::Print() << "DEBUG: Zeroing ylo boundary, box=" << bx_ylo << "\n"; + ParallelFor(bx_ylo, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (n_qstate >= 2) cons_arr(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cons_arr(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cons_arr(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cons_arr(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cons_arr(i,j,k,RhoQ6_comp) = zero; + }); + } + if (tbx.bigEnd(1) == j_hi) { + Box bx_yhi(tbx); bx_yhi.setSmall(1, j_hi - real_width + 1); + amrex::Print() << "DEBUG: Zeroing yhi boundary, box=" << bx_yhi << "\n"; + ParallelFor(bx_yhi, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (n_qstate >= 2) cons_arr(i,j,k,RhoQ2_comp) = zero; + if (n_qstate >= 3) cons_arr(i,j,k,RhoQ3_comp) = zero; + if (n_qstate >= 4) cons_arr(i,j,k,RhoQ4_comp) = zero; + if (n_qstate >= 5) cons_arr(i,j,k,RhoQ5_comp) = zero; + if (n_qstate >= 6) cons_arr(i,j,k,RhoQ6_comp) = zero; + }); + } + } + } + + // Check boundary values before microphysics + if (real_width > 0 && lev == 0 && n_qstate > 1) { + const auto& domain = geom[lev].Domain(); + int i_lo = domain.smallEnd(0); + Real max_q2 = cons.max(RhoQ2_comp); + Real max_q3 = cons.max(RhoQ3_comp); + amrex::Print() << "DEBUG: Before micro->Advance at time=" << time + << " max(RhoQ2)=" << max_q2 << " max(RhoQ3)=" << max_q3 << "\n"; + // Sample boundary cell + for (MFIter mfi(cons); mfi.isValid(); ++mfi) { + auto cons_arr = cons.array(mfi); + Box tbx = mfi.tilebox(); + if (tbx.smallEnd(0) == i_lo) { + amrex::Print() << "DEBUG: Sample xlo boundary at i=" << i_lo + << " RhoQ2=" << cons_arr(i_lo,tbx.smallEnd(1),tbx.smallEnd(2),RhoQ2_comp) + << " RhoQ3=" << cons_arr(i_lo,tbx.smallEnd(1),tbx.smallEnd(2),RhoQ3_comp) << "\n"; + break; + } + } + } + micro->Update_Micro_Vars_Lev(lev, cons); micro->Advance(lev, dt_advance, iteration, time, solverChoice, vars_new, z_phys_nd, phys_bc_type); micro->Update_State_Vars_Lev(lev, cons, *z_phys_nd[lev]); + + // Check boundary values after microphysics + if (real_width > 0 && lev == 0 && n_qstate > 1) { + Real max_q2 = cons.max(RhoQ2_comp); + Real max_q3 = cons.max(RhoQ3_comp); + amrex::Print() << "DEBUG: After micro->Advance at time=" << time + << " max(RhoQ2)=" << max_q2 << " max(RhoQ3)=" << max_q3 << "\n"; + } } } diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp index 4d9b0bc67c..775b05ec83 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp @@ -489,7 +489,7 @@ void erf_slow_rhs_post (int level, int finest_level, old_stage_time_total, dt, start_bdy_time, final_bdy_time, bdy_time_interval, bdy_factor, width, do_upwind, domain, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, - m_r2d); + m_r2d, n_qstate); } #endif From 46d9b9f533b9f92642184ef914e47e455eda0e71 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Wed, 27 May 2026 21:19:27 -0700 Subject: [PATCH 21/27] Cleanup from debugging hydrometeors bcs. --- Source/SourceTerms/ERF_MoistSetRhs.cpp | 39 ++-------------------- Source/SourceTerms/ERF_SrcHeaders.H | 3 +- Source/TimeIntegration/ERF_SlowRhsPost.cpp | 2 +- 3 files changed, 4 insertions(+), 40 deletions(-) diff --git a/Source/SourceTerms/ERF_MoistSetRhs.cpp b/Source/SourceTerms/ERF_MoistSetRhs.cpp index 53990979b6..a67dbef514 100644 --- a/Source/SourceTerms/ERF_MoistSetRhs.cpp +++ b/Source/SourceTerms/ERF_MoistSetRhs.cpp @@ -2,7 +2,7 @@ #include #include -#include +//#include using namespace amrex; @@ -29,8 +29,7 @@ moist_set_rhs (const Geometry& geom, Vector>& bdy_data_xhi, Vector>& bdy_data_ylo, Vector>& bdy_data_yhi, - std::unique_ptr& m_r2d, - int n_qstate) + std::unique_ptr& m_r2d) { // HACK HACK HACK // Get bndry data @@ -257,39 +256,5 @@ moist_set_rhs (const Geometry& geom, exit(0); */ - // Zero out hydrometeors in boundary regions based on active microphysics scheme - ParallelFor(tbx_xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; - if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; - if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; - if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; - if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; - }); - ParallelFor(tbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; - if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; - if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; - if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; - if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; - }); - ParallelFor(tbx_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; - if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; - if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; - if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; - if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; - }); - ParallelFor(tbx_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (n_qstate >= 2) cell_rhs(i,j,k,RhoQ2_comp) = zero; - if (n_qstate >= 3) cell_rhs(i,j,k,RhoQ3_comp) = zero; - if (n_qstate >= 4) cell_rhs(i,j,k,RhoQ4_comp) = zero; - if (n_qstate >= 5) cell_rhs(i,j,k,RhoQ5_comp) = zero; - if (n_qstate >= 6) cell_rhs(i,j,k,RhoQ6_comp) = zero; - }); - } // moist_set_rhs #endif diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H index 57adfa5bbd..bbe5b29c0c 100644 --- a/Source/SourceTerms/ERF_SrcHeaders.H +++ b/Source/SourceTerms/ERF_SrcHeaders.H @@ -140,8 +140,7 @@ moist_set_rhs (const amrex::Geometry& geom, amrex::Vector>& bdy_data_xhi, amrex::Vector>& bdy_data_ylo, amrex::Vector>& bdy_data_yhi, - std::unique_ptr& m_r2d, - int n_qstate); + std::unique_ptr& m_r2d); #endif void ApplySpongeZoneBCsForCC (const SpongeChoice& spongeChoice, diff --git a/Source/TimeIntegration/ERF_SlowRhsPost.cpp b/Source/TimeIntegration/ERF_SlowRhsPost.cpp index 775b05ec83..4d9b0bc67c 100644 --- a/Source/TimeIntegration/ERF_SlowRhsPost.cpp +++ b/Source/TimeIntegration/ERF_SlowRhsPost.cpp @@ -489,7 +489,7 @@ void erf_slow_rhs_post (int level, int finest_level, old_stage_time_total, dt, start_bdy_time, final_bdy_time, bdy_time_interval, bdy_factor, width, do_upwind, domain, bdy_data_xlo, bdy_data_xhi, bdy_data_ylo, bdy_data_yhi, - m_r2d, n_qstate); + m_r2d); } #endif From df7a35b6cedaba1fcf484b64ec56fd9e1c429a06 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 28 May 2026 14:15:03 -0700 Subject: [PATCH 22/27] Switch from using times from wrfinput and metgrid global attribute SIMULATION_START_DATE to the variable Times. --- Source/IO/ERF_ReadFromMetgrid.cpp | 16 ++++++++++++++-- Source/Initialization/ERF_InitFromMetgrid.cpp | 14 +++++++++++++- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/Source/IO/ERF_ReadFromMetgrid.cpp b/Source/IO/ERF_ReadFromMetgrid.cpp index e9adcba71c..7e80685763 100644 --- a/Source/IO/ERF_ReadFromMetgrid.cpp +++ b/Source/IO/ERF_ReadFromMetgrid.cpp @@ -76,8 +76,20 @@ read_from_metgrid (int lev, int itime, ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0]; } - { // Global Attributes (string) - NC_dateTime = ncf.get_attr("SIMULATION_START_DATE"); + { // Read Times variable (character array). + auto times_var = ncf.var("Times"); + std::vector start = {0, 0}; // Time index 0, character index 0. + std::vector count = times_var.shape(); + count[0] = 1; // Read the first time entry. + + // Allocate buffer for the time string. + std::vector time_chars(count[1]); + times_var.get(time_chars.data(), start, count); + + // Convert to std::string and trim trailing whitespace/null characters. + NC_dateTime = std::string(time_chars.begin(), time_chars.end()); + NC_dateTime.erase(NC_dateTime.find_last_not_of(" \0") + 1); + const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); } diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index a017172a4e..91162f8ceb 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -22,7 +22,19 @@ read_start_time_from_metgrid(int lev, const std::string& fname) if (ParallelDescriptor::IOProcessor()) { auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); - NC_dateTime = ncf.get_attr("SIMULATION_START_DATE"); + // Read the Times variable (character array). + auto times_var = ncf.var("Times"); + std::vector start = {0, 0}; // Time index 0, character index 0. + std::vector count = times_var.shape(); + count[0] = 1; // Read only the first time entry. + + // Allocate buffer for the time string. + std::vector time_chars(count[1]); + times_var.get(time_chars.data(), start, count); + + // Convert to std::string and trim trailing whitespace/null characters. + NC_dateTime = std::string(time_chars.begin(), time_chars.end()); + NC_dateTime.erase(NC_dateTime.find_last_not_of(" \0") + 1); const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); From 562343f8683ff2dca0855cfa395dd9fa51070952 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 28 May 2026 16:01:32 -0700 Subject: [PATCH 23/27] Make syntax consistent with what was already done for wrfbdy processing of Times. --- Source/Initialization/ERF_InitFromMetgrid.cpp | 49 ++++++++++--------- .../Initialization/ERF_InitFromWRFInput.cpp | 37 ++++++++++---- 2 files changed, 54 insertions(+), 32 deletions(-) diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 91162f8ceb..152866601c 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -17,31 +17,36 @@ using namespace amrex; Real read_start_time_from_metgrid(int lev, const std::string& fname) { - std::string NC_dateTime; - Real NC_epochTime; - if (ParallelDescriptor::IOProcessor()) { - auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); - - // Read the Times variable (character array). - auto times_var = ncf.var("Times"); - std::vector start = {0, 0}; // Time index 0, character index 0. - std::vector count = times_var.shape(); - count[0] = 1; // Read only the first time entry. - - // Allocate buffer for the time string. - std::vector time_chars(count[1]); - times_var.get(time_chars.data(), start, count); - - // Convert to std::string and trim trailing whitespace/null characters. - NC_dateTime = std::string(time_chars.begin(), time_chars.end()); - NC_dateTime.erase(NC_dateTime.find_last_not_of(" \0") + 1); + Real NC_epochTime; + const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; - const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; - NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); + if (ParallelDescriptor::IOProcessor()) { + // Read the time stamps + using CharArray = NDArray; + Vector array_ts(1); + Vector success(1); + ReadNetCDFFile(fname, {"Times"}, array_ts, success); + + int ntimes = array_ts[0].get_vshape()[0]; + auto dateStrLen = array_ts[0].get_vshape()[1]; + char timeStamps[ntimes][dateStrLen]; + + // Fill up the characters read + int str_len = static_cast(dateStrLen); + for (int nt(0); nt < ntimes; nt++) { + for (int dateStrCt(0); dateStrCt < str_len; dateStrCt++) { + auto n = nt*dateStrLen + dateStrCt; + timeStamps[nt][dateStrCt] = *(array_ts[0].get_data() + n); + } + } - ncf.close(); + // Extract the first time entry + std::string date(&timeStamps[0][0], &timeStamps[0][dateStrLen-1]+1); + auto epochTime = getEpochTime(date, dateTimeFormat); + Print() << " metgrid datetime 0 : " << date << " " << epochTime << std::endl; + NC_epochTime = static_cast(epochTime); - amrex::Print() << "Have read start_time string at level "<< lev << " is " << NC_dateTime << std::endl; + amrex::Print() << "Have read start_time string at level "<< lev << " is " << date << std::endl; amrex::Print() << "Have read start_time number at level "<< lev << " is " << NC_epochTime << std::endl; } diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index 521fcdc84f..896ea61a96 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -91,19 +91,36 @@ init_base_state_from_wrfinput (const Box& subdomain, Real read_start_time_from_wrfinput (int lev, const std::string& fname) { - std::string NC_dateTime; - Real NC_epochTime; - if (ParallelDescriptor::IOProcessor()) { - auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); + Real NC_epochTime; + const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; - NC_dateTime = ncf.get_attr("SIMULATION_START_DATE"); - - const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; - NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); + if (ParallelDescriptor::IOProcessor()) { + // Read the time stamps + using CharArray = NDArray; + Vector array_ts(1); + Vector success(1); + ReadNetCDFFile(fname, {"Times"}, array_ts, success); + + int ntimes = array_ts[0].get_vshape()[0]; + auto dateStrLen = array_ts[0].get_vshape()[1]; + char timeStamps[ntimes][dateStrLen]; + + // Fill up the characters read + int str_len = static_cast(dateStrLen); + for (int nt(0); nt < ntimes; nt++) { + for (int dateStrCt(0); dateStrCt < str_len; dateStrCt++) { + auto n = nt*dateStrLen + dateStrCt; + timeStamps[nt][dateStrCt] = *(array_ts[0].get_data() + n); + } + } - ncf.close(); + // Extract the first time entry + std::string date(&timeStamps[0][0], &timeStamps[0][dateStrLen-1]+1); + auto epochTime = getEpochTime(date, dateTimeFormat); + Print() << " wrfinput datetime 0 : " << date << " " << epochTime << std::endl; + NC_epochTime = static_cast(epochTime); - Print() << "Have read start_time string at level "<< lev << " is " << NC_dateTime << std::endl; + Print() << "Have read start_time string at level "<< lev << " is " << date << std::endl; Print() << "Have read start_time number at level "<< lev << " is " << NC_epochTime << std::endl; } From 84f068173ec597ee26422c17e05a2963932ba642 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 28 May 2026 16:42:22 -0700 Subject: [PATCH 24/27] Updated syntax for Time variable parsing in metgrid reader. --- Source/IO/ERF_ReadFromMetgrid.cpp | 65 ++++++++++++++----------------- 1 file changed, 29 insertions(+), 36 deletions(-) diff --git a/Source/IO/ERF_ReadFromMetgrid.cpp b/Source/IO/ERF_ReadFromMetgrid.cpp index 7e80685763..8ff641d96f 100644 --- a/Source/IO/ERF_ReadFromMetgrid.cpp +++ b/Source/IO/ERF_ReadFromMetgrid.cpp @@ -19,19 +19,15 @@ read_subdomain_from_metgrid(int /*lev*/, const std::string& fname, int& ratio, i ncf.get_attr("WEST-EAST_GRID_DIMENSION" , attr); nx = attr[0]-1; ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); ny = attr[0]-1; ncf.get_attr("BOTTOM-TOP_GRID_DIMENSION" , attr); nz = attr[0]-1; -#ifndef AMREX_USE_GPU Print() << "Have read (nx,ny,nz) = " << nx << " " << ny << " " << nz << std::endl; -#endif ncf.get_attr("i_parent_start", attr) ; is = attr[0]-1; ncf.get_attr("j_parent_start", attr) ; js = attr[0]-1; ncf.get_attr("parent_grid_ratio", attr); ratio = attr[0]; } ncf.close(); -#ifndef AMREX_USE_GPU amrex::Print() << "Have read (parent_ilo,parent_jlo) = " << is << " " << js << std::endl; amrex::Print() << "Have read refinement ratio = " << ratio << std::endl; -#endif } ParallelDescriptor::Bcast(&is ,1,amrex::ParallelDescriptor::IOProcessorNumber()); @@ -63,9 +59,7 @@ read_from_metgrid (int lev, int itime, IArrayBox& NC_lmask_iab, Geometry& geom) { -#ifndef AMREX_USE_GPU Print() << "Loading header data from NetCDF file at level " << lev << std::endl; -#endif if (ParallelDescriptor::IOProcessor()) { auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); @@ -76,50 +70,57 @@ read_from_metgrid (int lev, int itime, ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0]; } - { // Read Times variable (character array). - auto times_var = ncf.var("Times"); - std::vector start = {0, 0}; // Time index 0, character index 0. - std::vector count = times_var.shape(); - count[0] = 1; // Read the first time entry. - - // Allocate buffer for the time string. - std::vector time_chars(count[1]); - times_var.get(time_chars.data(), start, count); - - // Convert to std::string and trim trailing whitespace/null characters. - NC_dateTime = std::string(time_chars.begin(), time_chars.end()); - NC_dateTime.erase(NC_dateTime.find_last_not_of(" \0") + 1); - - const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; - NC_epochTime = getEpochTime(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(); + // Read the time stamps + const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S"; + using CharArray = NDArray; + Vector array_ts(1); + Vector success(1); + ReadNetCDFFile(fname, {"Times"}, array_ts, success); + + int ntimes = array_ts[0].get_vshape()[0]; + auto dateStrLen = array_ts[0].get_vshape()[1]; + char timeStamps[ntimes][dateStrLen]; + + // Fill up the characters read + int str_len = static_cast(dateStrLen); + for (int nt(0); nt < ntimes; nt++) { + for (int dateStrCt(0); dateStrCt < str_len; dateStrCt++) { + auto n = nt*dateStrLen + dateStrCt; + timeStamps[nt][dateStrCt] = *(array_ts[0].get_data() + n); + } + } + + // Extract the first time entry + std::string date(&timeStamps[0][0], &timeStamps[0][dateStrLen-1]+1); + auto epochTime = getEpochTime(date, dateTimeFormat); + Print() << "parsed metgrid datetime " << date << " " << epochTime << std::endl; + + NC_dateTime = date; + NC_epochTime = static_cast(epochTime); + // Verify the inputs geometry matches what the NETCDF file has if (lev == 0) { Real tol = Real(1.0e-3); Real Len_x = NC_dx * Real(NC_nx-1); Real Len_y = NC_dy * Real(NC_ny-1); if (std::fabs(Len_x - (geom.ProbHi(0) - geom.ProbLo(0))) > tol) { -#ifndef AMREX_USE_GPU Print() << "X problem extent " << (geom.ProbHi(0) - geom.ProbLo(0)) << " does not match NETCDF file " << Len_x << "!\n"; Print() << "dx: " << NC_dx << ' ' << "Nx: " << NC_nx-1 << "\n"; -#endif Abort("Domain specification error"); } if (std::fabs(Len_y - (geom.ProbHi(1) - geom.ProbLo(1))) > tol) { -#ifndef AMREX_USE_GPU Print() << "Y problem extent " << (geom.ProbHi(1) - geom.ProbLo(1)) << " does not match NETCDF file " << Len_y << "!\n"; Print() << "dy: " << NC_dy << ' ' << "Ny: " << NC_ny-1 << "\n"; -#endif Abort("Domain specification error"); } } // lev == 0 @@ -132,9 +133,7 @@ read_from_metgrid (int lev, int itime, ParallelDescriptor::Bcast(&NC_dx, 1, ioproc); ParallelDescriptor::Bcast(&NC_dy, 1, ioproc); -#ifndef AMREX_USE_GPU Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; -#endif Vector NC_fabs; Vector NC_iabs; @@ -165,9 +164,7 @@ read_from_metgrid (int lev, int itime, 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 -#ifndef AMREX_USE_GPU Print() << "Building initial FABS from file " << fname << std::endl; -#endif Vector success; success.resize(NC_fabs.size()); BuildFABsFromNetCDFFile(domain, fname, NC_fnames, NC_fdim_types, NC_fabs, success); @@ -179,22 +176,18 @@ read_from_metgrid (int lev, int itime, } // Read the netcdf file and fill these IABs -#ifndef AMREX_USE_GPU Print() << "Building initial IABS from file " << fname << std::endl; -#endif Vector success_i; success_i.resize(NC_iabs.size()); BuildFABsFromNetCDFFile(domain, fname, NC_inames, NC_idim_types, NC_iabs, success_i); for (int i = 0; i < success_i.size(); i++) { flag_lmask = (NC_inames[i] == "LANDMASK" && success_i[i] == 1) ? 1 : 0; } -#ifndef AMREX_USE_GPU Print() << " flag_psfc: " << flag_psfc << std::endl; Print() << " flag_sst: " << flag_sst << std::endl; Print() << " flag_tsk: " << flag_tsk << std::endl; Print() << " flag_msf: " << flag_msf << std::endl; Print() << " flag_lmask: " << flag_lmask << std::endl; -#endif // TODO: FIND OUT IF WE NEED TO DIVIDE VELS BY MAPFAC // From 10c992598a423c57496280f4df0cf962389aa197 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 28 May 2026 20:37:31 -0700 Subject: [PATCH 25/27] Rollback format changes left over from debugging --- Source/SourceTerms/ERF_MoistSetRhs.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Source/SourceTerms/ERF_MoistSetRhs.cpp b/Source/SourceTerms/ERF_MoistSetRhs.cpp index a67dbef514..c967254c68 100644 --- a/Source/SourceTerms/ERF_MoistSetRhs.cpp +++ b/Source/SourceTerms/ERF_MoistSetRhs.cpp @@ -2,7 +2,6 @@ #include #include -//#include using namespace amrex; @@ -255,6 +254,5 @@ moist_set_rhs (const Geometry& geom, }); exit(0); */ - } // moist_set_rhs #endif From f7c1775c0714f5dc3acfffc00bd7520dbc1c900f Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Thu, 28 May 2026 22:13:30 -0700 Subject: [PATCH 26/27] Avoid writing QV to erfbdy if no moisture model is enabled. --- Source/Initialization/ERF_InitFromMetgrid.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 478e6e4950..fb6299fe45 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -77,7 +77,10 @@ ERF::init_from_metgrid (int lev) std::string erfbdy_header = erfbdy_file + "/Header"; use_erfbdy = FileSystem::Exists(erfbdy_header); } - if (use_erfbdy || write_erfbdy) nvars_erfbdy = MetGridBdyVars::NumTypes; + // Set nvars_erfbdy based on whether moisture is enabled + if (use_erfbdy || write_erfbdy) { + nvars_erfbdy = use_moisture ? MetGridBdyVars::NumTypes : (MetGridBdyVars::NumTypes - 1); + } // If the erfbdy file exists, load boundary data and skip met_em processing. if (lev == 0 && use_erfbdy) { From 344dfc9b6f2244d7e1683f4b0935eaee56954264 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Tue, 23 Jun 2026 16:56:21 -0700 Subject: [PATCH 27/27] Accidentally mucked up Noah-MP submodule. Reverted to match development --- Submodules/Noah-MP | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Submodules/Noah-MP b/Submodules/Noah-MP index c08dd40676..1900918c1f 160000 --- a/Submodules/Noah-MP +++ b/Submodules/Noah-MP @@ -1 +1 @@ -Subproject commit c08dd406761330ac823a3e268bd2e98eb6889b40 +Subproject commit 1900918c1fe2cb218bca97a3cdbcb8bcf855a935