ERF boundary file (erfbdy)

CMake/BuildERFExe.cmake

@@ -353,6 +353,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

Docs/sphinx_doc/BoundaryConditions.rst

@@ -233,3 +233,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.

Docs/sphinx_doc/Inputs.rst

@@ -635,6 +635,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
 ----------------------------
 

Source/ERF.H

@@ -1301,6 +1301,12 @@ private:
     int         metgrid_order{2};
     int         metgrid_force_sfc_k{6};
 
+    // 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<amrex::BoxArray> ba1d;
     amrex::Vector<amrex::BoxArray> ba2d;
 

Source/ERF.cpp

@@ -2551,6 +2551,10 @@ 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("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);
@@ -2762,6 +2766,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);

Source/IO/ERF_Checkpoint.cpp

@@ -6,6 +6,7 @@
 
 #include "ERF.H"
 #include "AMReX_PlotFileUtil.H"
+#include "ERF_ReadFromERFBdy.H"
 
 using namespace amrex;
 
@@ -1143,6 +1144,55 @@ ERF::ReadCheckpointFile ()
     } // init_type == WRFInput or Metgrid
 #endif
 #endif
+
+#ifdef ERF_USE_NETCDF
+    // 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) {
+
+        // Check for erfbdy file.
+        std::string erfbdy_header = erfbdy_file + "/Header";
+        use_erfbdy = FileSystem::Exists(erfbdy_header);
+
+        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<Real> 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<int>(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;
+            }
+        }
+    }
+#endif
 }
 
 /**

Source/IO/ERF_ReadFromERFBdy.H

@@ -0,0 +1,28 @@
+#ifndef ERF_READ_FROM_ERFBDY_H_
+#define ERF_READ_FROM_ERFBDY_H_
+
+#include <AMReX_FArrayBox.H>
+#include <AMReX_Vector.H>
+
+// 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<amrex::Real>& bdy_times,
+                       amrex::Real& start_bdy_time,
+                       amrex::Real& final_bdy_time);
+
+// Read a single time index from erfbdy file.
+void
+read_from_erfbdy(int itime,
+                 const std::string& bdy_file_name,
+                 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_xlo,
+                 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_xhi,
+                 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_ylo,
+                 amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_yhi,
+                 int nvars, int real_width);
+
+#endif

Source/IO/ERF_ReadFromERFBdy.cpp

@@ -0,0 +1,133 @@
+#include "ERF_ReadFromERFBdy.H"
+#include <AMReX_VisMF.H>
+#include <AMReX_ParallelDescriptor.H>
+#include <AMReX_Print.H>
+#include <fstream>
+
+using namespace amrex;
+
+Real
+read_times_from_erfbdy(const std::string& bdy_file_name,
+                       int& ntimes,
+                       int& nvars,
+                       int& real_width,
+                       Vector<Real>& 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 boundary times.
+    bdy_times.resize(ntimes);
+    for (int i = 0; i < ntimes; ++i) {
+        HeaderFile >> bdy_times[i];
+    }
+
+    // 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();
+
+    // Ensure the file holds at least two times.
+    AMREX_ALWAYS_ASSERT(ntimes >= 2);
+
+    // Set the first and last boundary times.
+    start_bdy_time = bdy_times[0];
+    final_bdy_time = bdy_times[ntimes - 1];
+
+    // Calculate the interval between boundary times.
+    Real bdy_time_interval = bdy_times[1] - bdy_times[0];
+
+    return bdy_time_interval;
+}
+
+void
+read_from_erfbdy(int itime,
+                 const std::string& bdy_file_name,
+                 Vector<Vector<FArrayBox>>& bdy_data_xlo,
+                 Vector<Vector<FArrayBox>>& bdy_data_xhi,
+                 Vector<Vector<FArrayBox>>& bdy_data_ylo,
+                 Vector<Vector<FArrayBox>>& 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);
+
+    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);
+        bdy_data_ylo[itime].resize(nvars);
+        bdy_data_yhi[itime].resize(nvars);
+    }
+
+    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);
+            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<RunOn::Host>(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);
+            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<RunOn::Host>(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);
+            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<RunOn::Host>(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);
+            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<RunOn::Host>(tmp_fab, 0, 0, tmp_fab.nComp());
+            ifs.close();
+        }
+    }
+
+    // Barrier to ensure all reads complete.
+    ParallelDescriptor::Barrier();
+}

Source/IO/ERF_WriteERFBdy.H

@@ -0,0 +1,25 @@
+#ifndef ERF_WRITE_ERFBDY_H_
+#define ERF_WRITE_ERFBDY_H_
+
+#include <AMReX_Box.H>
+#include <AMReX_FArrayBox.H>
+#include <AMReX_Vector.H>
+
+// Initialize erfbdy file and write header.
+void InitERFBdyFile(const std::string& bdy_file_name,
+                    int ntimes,
+                    const amrex::Vector<amrex::Real>& bdy_times,
+                    const amrex::Box& domain,
+                    int nvars,
+                    int real_width);
+
+// Write a single time index.
+void WriteERFBdyTimeSlice(const std::string& bdy_file_name,
+                          int itime,
+                          const amrex::Vector<amrex::FArrayBox>& bdy_data_xlo,
+                          const amrex::Vector<amrex::FArrayBox>& bdy_data_xhi,
+                          const amrex::Vector<amrex::FArrayBox>& bdy_data_ylo,
+                          const amrex::Vector<amrex::FArrayBox>& bdy_data_yhi,
+                          int nvars);
+
+#endif