Read input file to determine dimensions

.github/workflows/test_suite.yml

@@ -47,7 +47,7 @@ jobs:
   # Serial tests on latest Ubuntu OS
   test-ubuntu-serial:
     runs-on: ubuntu-latest
-    timeout-minutes: 10
+    timeout-minutes: 15
     container:
       image: ghcr.io/nextsimhub/nextsimdg-dev-env:latest
     steps:
@@ -194,7 +194,7 @@ jobs:
   # Serial tests on latest Mac OS
   test-mac-serial:
     runs-on: macos-latest
-    timeout-minutes: 10
+    timeout-minutes: 15
     steps:
     - uses: actions/checkout@v2
       with:

core/src/Model.cpp

@@ -71,6 +71,9 @@ void Model::configureRestarts()
     std::string restartPeriodStr
         = Configured::getConfiguration(keyMap.at(RESTARTPERIOD_KEY), std::string("0"));
     metadata.restartPeriod = Duration(restartPeriodStr);
+
+    // Set global dimensions with an initial read of the input file
+    metadata.setDimensionsFromFile(metadata.initialFileName);
 }
 
 void Model::configureTime()

core/src/ModelMetadata.cpp

@@ -9,18 +9,23 @@
 #include "include/IStructure.hpp"
 #include "include/ModelMPI.hpp"
 #include "include/NextsimModule.hpp"
+#ifdef USE_MPI
+#include "include/ParallelNetcdfFile.hpp"
 #ifdef USE_XIOS
 #include "include/Xios.hpp"
+#else
+#include "include/Halo.hpp"
+#endif
 #endif
 #include "include/gridNames.hpp"
 #include <array>
 #include <cstddef>
 #include <functional>
 #include <vector>
 
+#include <ncDim.h>
 #ifdef USE_MPI
 #include "mpi.h"
-#include <ncDim.h>
 #include <ncFile.h>
 #include <ncGroup.h>
 #include <ncVar.h>
@@ -119,8 +124,18 @@ void ModelMetadata::getPartitionMetadata(std::string partitionFile)
             + std::to_string(nBoxes) + "\n";
         throw std::runtime_error(errorMsg);
     }
-    globalExtentX = ncFile.getDim("NX").getSize();
-    globalExtentY = ncFile.getDim("NY").getSize();
+    if (!globalExtentX) {
+        globalExtentX = ncFile.getDim("NX").getSize();
+    } else if (globalExtentX != ncFile.getDim("NX").getSize()) {
+        throw std::runtime_error("ModelMetadata: Inconsistent global x-extent between "
+                                 "partition and input files.");
+    }
+    if (!globalExtentY) {
+        globalExtentY = ncFile.getDim("NY").getSize();
+    } else if (globalExtentX != ncFile.getDim("NY").getSize()) {
+        throw std::runtime_error("ModelMetadata: Inconsistent global y-extent between "
+                                 "partition and input files.");
+    }
     netCDF::NcGroup bboxGroup(ncFile.getGroup(bboxName));
 
     std::vector<size_t> rank(1, modelMPI.getRank());
@@ -150,6 +165,85 @@ ModelMetadata::ModelMetadata()
 
 #endif
 
+void ModelMetadata::setDimensionsFromFile(const std::string& filename)
+{
+    if (filename.empty()) {
+        throw std::runtime_error(
+            "ModelMetadata :: setDimensionsFromFile() called without input file.");
+    }
+    try {
+#ifdef USE_MPI
+        auto& modelMPI = ModelMPI::getInstance();
+        netCDF::NcFilePar ncFile(filename, netCDF::NcFile::read, modelMPI.getComm());
+#else
+        netCDF::NcFile ncFile(filename, netCDF::NcFile::read);
+#endif
+
+        // Dimensions and DG components
+        std::multimap<std::string, netCDF::NcDim> dimMap = ncFile.getDims();
+        for (auto entry : ModelArray::definedDimensions) {
+            auto dimType = entry.first;
+            if (dimType == ModelArray::Dimension::DG || dimType == ModelArray::Dimension::DGSTRESS
+                || dimType == ModelArray::Dimension::NCOORDS) {
+                // TODO: Assert that DG in the file equals the compile time DG in the model (#205)
+                continue;
+            }
+
+            ModelArray::DimensionSpec& dimensionSpec = entry.second;
+            // Find dimensions in the netCDF file by their name in the ModelArray details
+            netCDF::NcDim dim = ncFile.getDim(dimensionSpec.name);
+            // Also check the old name
+            if (dim.isNull()) {
+                dim = ncFile.getDim(dimensionSpec.altName);
+            }
+            // If we didn't find a dimension with the dimensions name or altName, throw.
+            if (dim.isNull()) {
+                throw std::out_of_range(
+                    "ModelMetadata: No netCDF dimension found corresponding to the dimension named "
+                    + dimensionSpec.name + " or " + dimensionSpec.altName);
+            }
+#ifdef USE_MPI
+            size_t localLength;
+            size_t start;
+            if (dimType == ModelArray::Dimension::X) {
+                localLength = getLocalExtentX();
+                start = getLocalCornerX();
+            } else if (dimType == ModelArray::Dimension::Y) {
+                localLength = getLocalExtentY();
+                start = getLocalCornerY();
+            } else if (dimType == ModelArray::Dimension::XVERTEX) {
+                localLength = getLocalExtentX() + 1;
+                start = getLocalCornerX();
+            } else if (dimType == ModelArray::Dimension::YVERTEX) {
+                localLength = getLocalExtentY() + 1;
+                start = getLocalCornerY();
+            } else if (dimType == ModelArray::Dimension::XCG) {
+                localLength = CGDEGREE * getLocalExtentX() + 1;
+                start = CGDEGREE * getLocalCornerX();
+            } else if (dimType == ModelArray::Dimension::YCG) {
+                localLength = CGDEGREE * getLocalExtentY() + 1;
+                start = CGDEGREE * getLocalCornerY();
+            } else {
+                localLength = dim.getSize();
+                start = 0;
+            }
+#if USE_XIOS
+            ModelArray::setDimension(dimType, dim.getSize(), localLength, start);
+#else
+            ModelArray::setDimension(
+                dimType, dim.getSize(), localLength + 2 * Halo::haloWidth, start);
+#endif
+#else
+            ModelArray::setDimension(dimType, dim.getSize());
+#endif
+        }
+    } catch (const netCDF::exceptions::NcException& nce) {
+        std::string ncWhat(nce.what());
+        ncWhat += ": " + filename;
+        throw std::runtime_error(ncWhat);
+    }
+}
+
 const ModelState& ModelMetadata::extractCoordinates(const ModelState& state)
 {
     // More sophisticated grids include both vertex coordinates and grid azimuth values.

core/src/ParaGridIO.cpp

@@ -83,6 +83,10 @@ ModelState ParaGridIO::getModelState(const std::string& filePath)
 {
     ModelState state;
 
+    // Set global dimensions with an initial read of the input file
+    ModelMetadata& metadata = ModelMetadata::getInstance();
+    metadata.setDimensionsFromFile(filePath);
+
     try {
 #ifdef USE_MPI
         auto& modelMPI = ModelMPI::getInstance();
@@ -91,60 +95,7 @@ ModelState ParaGridIO::getModelState(const std::string& filePath)
         netCDF::NcFile ncFile(filePath, netCDF::NcFile::read);
 #endif
 
-        // Dimensions and DG components
-        std::multimap<std::string, netCDF::NcDim> dimMap = ncFile.getDims();
-        for (auto entry : ModelArray::definedDimensions) {
-            auto dimType = entry.first;
-            if (dimCompMap.count(dimType) > 0)
-                // TODO Assertions that DG in the file equals the compile time DG in the model. See
-                // #205
-                continue;
-
-            ModelArray::DimensionSpec& dimensionSpec = entry.second;
-            // Find dimensions in the netCDF file by their name in the ModelArray details
-            netCDF::NcDim dim = ncFile.getDim(dimensionSpec.name);
-            // Also check the old name
-            if (dim.isNull()) {
-                dim = ncFile.getDim(dimensionSpec.altName);
-            }
-            // If we didn't find a dimension with the dimensions name or altName, throw.
-            if (dim.isNull()) {
-                throw std::out_of_range(
-                    std::string("No netCDF dimension found corresponding to the dimension named ")
-                    + dimensionSpec.name + std::string(" or ") + dimensionSpec.altName);
-            }
-#ifdef USE_MPI
-            auto dimName = dim.getName();
-            size_t localLength = 0;
-            size_t start = 0;
-            auto& metadata = ModelMetadata::getInstance();
-            if (dimType == ModelArray::Dimension::X) {
-                localLength = metadata.getLocalExtentX();
-                start = metadata.getLocalCornerX();
-            } else if (dimType == ModelArray::Dimension::Y) {
-                localLength = metadata.getLocalExtentY();
-                start = metadata.getLocalCornerY();
-            } else if (dimType == ModelArray::Dimension::XVERTEX) {
-                localLength = metadata.getLocalExtentX() + 1;
-                start = metadata.getLocalCornerX();
-            } else if (dimType == ModelArray::Dimension::YVERTEX) {
-                localLength = metadata.getLocalExtentY() + 1;
-                start = metadata.getLocalCornerY();
-            } else {
-                localLength = dim.getSize();
-                start = 0;
-            }
-            // globalLength doesnt need to be padded with halo cells but localLength does
-            // setDimension(dim, globalLength, localLength, start)
-            ModelArray::setDimension(
-                dimType, dim.getSize(), localLength + 2 * Halo::haloWidth, start);
-#else
-            ModelArray::setDimension(dimType, dim.getSize());
-#endif
-        }
-
         // Get all valid variables and load them into a new ModelState
-
         for (auto entry : ncFile.getVars()) {
             const std::string& varName = entry.first;
             netCDF::NcVar& var = entry.second;

core/src/Xios.cpp

@@ -486,84 +486,35 @@ void Xios::parseInputFiles()
     auto& metadata = ModelMetadata::getInstance();
 
     // Initial read of the NetCDF file to deduce the dimensions
-    for (std::string filename : { inputFilename, forcingFilename }) {
-        if (filename.length() == 0) {
+    for (const std::string& filename : { inputFilename, forcingFilename }) {
+        if (filename.empty()) {
             break;
         }
-
+        metadata.setDimensionsFromFile(filename);
+
+        // Create map for field types
+        const std::map<std::string, ModelArray::Type> dimensionKeys = {
+            { "yx", ModelArray::Type::H },
+            { "ydimxdim", ModelArray::Type::H },
+            { "yxdg_comp", ModelArray::Type::DG },
+            { "ydimxdimdg_comp", ModelArray::Type::DG },
+            { "yxdgstress_comp", ModelArray::Type::DGSTRESS },
+            { "ydimxdimdgstress_comp", ModelArray::Type::DGSTRESS },
+            { "y_cgx_cg", ModelArray::Type::CG },
+            { "yvertexxvertexncoords", ModelArray::Type::VERTEX },
+        };
+
+        // Determine field types
+        std::set<std::string> configFieldIds;
+        if (filename == inputFilename) {
+            configFieldIds = configGetInputRestartFieldNames();
+        } else {
+            configFieldIds = configGetForcingFieldNames();
+        }
         try {
             auto& modelMPI = ModelMPI::getInstance();
             netCDF::NcFilePar ncFile(filename, netCDF::NcFile::read, modelMPI.getComm());
 
-            // Dimensions and DG components
-            std::multimap<std::string, netCDF::NcDim> dimMap = ncFile.getDims();
-            for (auto entry : ModelArray::definedDimensions) {
-                auto dimType = entry.first;
-                ModelArray::DimensionSpec& dimensionSpec = entry.second;
-                // TODO: Assert that DG in the file equals the compile time DG in the model (#205)
-
-                // Find dimensions in the netCDF file by their name in the ModelArray details
-                netCDF::NcDim dim = ncFile.getDim(dimensionSpec.name);
-
-                // Also check the old name
-                if (dim.isNull()) {
-                    dim = ncFile.getDim(dimensionSpec.altName);
-                }
-
-                // If we didn't find a dimension with the dimensions name or altName, throw.
-                if (dim.isNull()) {
-                    throw std::out_of_range(
-                        "Xios: No netCDF dimension found corresponding to the dimension named "
-                        + dimensionSpec.name + " or " + dimensionSpec.altName);
-                }
-
-                // Set corresponding dimensions
-                size_t localLength;
-                size_t start;
-                if (dimType == ModelArray::Dimension::X) {
-                    localLength = metadata.getLocalExtentX();
-                    start = metadata.getLocalCornerX();
-                } else if (dimType == ModelArray::Dimension::Y) {
-                    localLength = metadata.getLocalExtentY();
-                    start = metadata.getLocalCornerY();
-                } else if (dimType == ModelArray::Dimension::XVERTEX) {
-                    localLength = metadata.getLocalExtentX() + 1;
-                    start = metadata.getLocalCornerX();
-                } else if (dimType == ModelArray::Dimension::YVERTEX) {
-                    localLength = metadata.getLocalExtentY() + 1;
-                    start = metadata.getLocalCornerY();
-                } else if (dimType == ModelArray::Dimension::XCG) {
-                    localLength = CGDEGREE * metadata.getLocalExtentX() + 1;
-                    start = CGDEGREE * metadata.getLocalCornerX();
-                } else if (dimType == ModelArray::Dimension::YCG) {
-                    localLength = CGDEGREE * metadata.getLocalExtentY() + 1;
-                    start = CGDEGREE * metadata.getLocalCornerY();
-                } else {
-                    localLength = dim.getSize();
-                    start = 0;
-                }
-                ModelArray::setDimension(dimType, dim.getSize(), localLength, start);
-            }
-
-            // Create map for field types
-            const std::map<std::string, ModelArray::Type> dimensionKeys = {
-                { "yx", ModelArray::Type::H },
-                { "ydimxdim", ModelArray::Type::H },
-                { "yxdg_comp", ModelArray::Type::DG },
-                { "ydimxdimdg_comp", ModelArray::Type::DG },
-                { "yxdgstress_comp", ModelArray::Type::DGSTRESS },
-                { "ydimxdimdgstress_comp", ModelArray::Type::DGSTRESS },
-                { "y_cgx_cg", ModelArray::Type::CG },
-                { "yvertexxvertexncoords", ModelArray::Type::VERTEX },
-            };
-
-            // Determine field types
-            std::set<std::string> configFieldIds;
-            if (filename == inputFilename) {
-                configFieldIds = configGetInputRestartFieldNames();
-            } else {
-                configFieldIds = configGetForcingFieldNames();
-            }
             for (auto entry : ncFile.getVars()) {
                 const std::string& fieldId = entry.first;
                 // Only consider fields that appear in the config
@@ -600,8 +551,8 @@ void Xios::parseInputFiles()
 /*!
  * @brief   Create XIOS domains associated with each ModelArray type
  *
- * @details This function sets up the XIOS domains for each field type based on the configuration
- *          in the domainIds map and in the ModelMetadata class.
+ * @details This function sets up the XIOS domains for each field type based on the
+ *          configuration in the domainIds map and in the ModelMetadata class.
  */
 void Xios::setupDomains()
 {

core/src/include/ModelMetadata.hpp

@@ -71,6 +71,16 @@ class ModelMetadata {
     }
 #endif
 
+    /*!
+     * @brief   Set dimensions information based on the contents of an input file.
+     *
+     * @param   filename the name of the file
+     * @details If an input file hasn't been read yet, the dimensions are read from the file and
+     *          set. Otherwise, a consistency check is made against the dimensions read from file
+     *          and already set.
+     */
+    void setDimensionsFromFile(const std::string& filename);
+
     // finalize ModelMetadata
     static void finalize();
 

core/test/CMakeLists.txt

@@ -202,7 +202,11 @@ if(ENABLE_MPI)
         target_link_libraries(testXiosAxis_MPI3 PRIVATE nextsimlib doctest::doctest)
 
         add_executable(testXiosField_MPI3 "XiosField_test.cpp" "MainMPI.cpp")
-        target_compile_definitions(testXiosField_MPI3 PRIVATE USE_XIOS)
+        target_compile_definitions(testXiosField_MPI3
+            PRIVATE
+                USE_XIOS
+                TEST_FILES_DIR=\"${CMAKE_CURRENT_BINARY_DIR}\"
+        )
         target_include_directories(
             testXiosField_MPI3
             PRIVATE