Post-process dims in XIOS output files

core/src/ModelMetadata.cpp

@@ -1,12 +1,14 @@
 /*!
  * @author  Tim Spain <timothy.spain@nersc.no>
  * @author  Tom Meltzer <tdm39@cam.ac.uk>
+ * @author  Joe Wallwork <jw2423@cam.ac.uk>
  */
 
 #include "include/ModelMetadata.hpp"
 
 #include "include/Finalizer.hpp"
 #include "include/IStructure.hpp"
+#include "include/Logged.hpp"
 #include "include/ModelMPI.hpp"
 #include "include/NextsimModule.hpp"
 #ifdef USE_MPI
@@ -171,6 +173,10 @@ void ModelMetadata::setDimensionsFromFile(const std::string& filename)
         throw std::runtime_error(
             "ModelMetadata :: setDimensionsFromFile() called without input file.");
     }
+
+    // Set to record the names of dimensions found in the file
+    std::set<std::string> dimNames;
+
     try {
 #ifdef USE_MPI
         auto& modelMPI = ModelMPI::getInstance();
@@ -196,26 +202,13 @@ void ModelMetadata::setDimensionsFromFile(const std::string& filename)
             if (dim.isNull()) {
                 dim = ncFile.getDim(dimensionSpec.altName);
             }
-#ifdef USE_XIOS
-            // Account for the fact that XIOS writes dimensions differently if only one
-            // discretisation is written out. If the dimension is still null at this point then we
-            // assume that the only discretisation used is HDomain-based, i.e., HField, DGField, or
-            // DGSField.
             if (dim.isNull()) {
-                if (dimensionSpec.name == "x_dim") {
-                    dim = ncFile.getDim("x");
-                } else if (dimensionSpec.name == "y_dim") {
-                    dim = ncFile.getDim("y");
-                } else {
-                    continue;
-                }
-            }
-#endif
-            // If we didn't find a dimension with the dimensions name or altName, throw.
-            if (dim.isNull()) {
-                throw std::out_of_range(
+                Logged::warning(
                     "ModelMetadata: No netCDF dimension found corresponding to the dimension named "
                     + dimensionSpec.name + " or " + dimensionSpec.altName);
+                continue;
+            } else {
+                dimNames.insert(dimensionSpec.name);
             }
 #ifdef USE_MPI
             size_t localLength;
@@ -257,6 +250,12 @@ void ModelMetadata::setDimensionsFromFile(const std::string& filename)
         ncWhat += ": " + filename;
         throw std::runtime_error(ncWhat);
     }
+
+    // Throw an error if we didn't find any dimensions
+    if (dimNames.empty()) {
+        throw std::out_of_range(
+            "ModelMetadata: No netCDF dimensions in input file '" + filename + "'");
+    }
 }
 
 const ModelState& ModelMetadata::extractCoordinates(const ModelState& state)

core/src/Time.cpp

@@ -236,6 +236,8 @@ void Duration::setDurationSeconds(double secs)
 
 TimePoint Duration::operator+(const TimePoint& t) const { return t + *this; }
 
+TimePoint Duration::operator-(const TimePoint& t) const { return t - *this; }
+
 std::tm* TimePoint::gmtime() const
 {
     auto tt = Clock::to_time_t(m_t);

core/src/Xios.cpp

@@ -162,6 +162,7 @@ void Xios::close_context_definition()
 void Xios::context_finalize()
 {
     if (isEnabled) {
+        postprocessOutputFiles();
         cxios_context_finalize();
     }
 }
@@ -962,13 +963,10 @@ void Xios::setupFields()
 
         // Create map for field types
         const std::map<std::string, ModelArray::Type> dimensionKeys = {
-            { "yx", ModelArray::Type::H },
             { "ydimxdim", ModelArray::Type::H },
             { "y_dimx_dim", ModelArray::Type::H },
-            { "yxdg_comp", ModelArray::Type::DG },
             { "ydimxdimdg_comp", ModelArray::Type::DG },
             { "y_dimx_dimdg_comp", ModelArray::Type::DG },
-            { "yxdgstress_comp", ModelArray::Type::DGSTRESS },
             { "ydimxdimdgstress_comp", ModelArray::Type::DGSTRESS },
             { "y_dimx_dimdgstress_comp", ModelArray::Type::DGSTRESS },
             { "y_cgx_cg", ModelArray::Type::CG },
@@ -1385,6 +1383,77 @@ void Xios::setupFiles()
     }
 }
 
+/*!
+ * @brief   Postprocess output files after the simulation has completed.
+ *
+ * @details If only a single domain was written to file, rename the x and y dimensions and variables
+ *          to x_dim and y_dim, respectively, for compatibility with other model components.
+ */
+void Xios::postprocessOutputFiles()
+{
+    // Count how many domains were written out
+    int sum = 0;
+    for (const auto& [domainId, written] : domainWritten) {
+        if (written) {
+            sum++;
+        }
+    }
+
+    // If a single domain was written then modify x and y dimensions and variables in the output
+    // files
+    if (sum == 1 && mpi_rank == 0) {
+
+        // Consider both restart files and diagnostic files
+        for (std::string fileId : { outputFileId, diagnosticFileId }) {
+            bool exists;
+            cxios_file_valid_id(&exists, fileId.c_str(), fileId.length());
+            if (!exists) {
+                continue;
+            }
+            Duration step = getFileOutputFreq(fileId);
+
+            // Loop over the output window splits
+            ModelMetadata& metadata = ModelMetadata::getInstance();
+            TimePoint time = metadata.startTime();
+            TimePoint endTime = metadata.stopTime();
+            while (time < endTime) {
+
+                // Compute the end time of the window, subtracting 1 second to avoid overlap
+                TimePoint nextTime = time + step - Duration(1);
+
+                // Generate the filename used by XIOS
+                // TODO: Support file patterns (#898)
+                std::string filename = fileId + "_" + time.format("%Y%m%d%H%M%S") + "-"
+                    + nextTime.format("%Y%m%d%H%M%S") + ".nc";
+
+                // Increment the time then check if the file exists
+                time += step;
+                if (!std::filesystem::exists(filename)) {
+                    continue;
+                }
+
+                try {
+                    // Open the netCDF file for both reading and writing
+                    netCDF::NcFile ncFile(filename, netCDF::NcFile::write);
+
+                    // Rename the x and y dimensions with x_dim and y_dim, respectively
+                    ncFile.getDim("x").rename("x_dim");
+                    ncFile.getDim("y").rename("y_dim");
+
+                    // Rename the x and y variables with x_dim and y_dim, respectively
+                    ncFile.getVar("x").rename("x_dim");
+                    ncFile.getVar("y").rename("y_dim");
+
+                    // Ensure changes are flushed to disk before closing
+                    ncFile.sync();
+                } catch (const netCDF::exceptions::NcException& e) {
+                    std::cerr << "Error processing NetCDF file: " << e.what() << std::endl;
+                }
+            }
+        }
+    }
+}
+
 /*!
  * Send a field to the XIOS server to be written to file
  *
@@ -1402,6 +1471,7 @@ void Xios::write(const std::string& fieldId, const ModelArray& modelarray)
     }
     auto& dims = modelarray.dimensions();
     const ModelArray::Type& type = modelarray.getType();
+    domainWritten[domainIds[type]] = true;
     if ((type == ModelArray::Type::H) || (type == ModelArray::Type::CG)) {
         cxios_write_data_k82(
             fieldId.c_str(), fieldId.length(), modelarray.getData(), dims[0], dims[1], -1);

core/src/include/Time.hpp

@@ -47,6 +47,9 @@ class Duration {
     //! Add a Duration to a TimePoint to get a TimePoint (now + 7 days = next week).
     TimePoint operator+(const TimePoint& t) const;
 
+    //! Subtract a Duration from a TimePoint to get a TimePoint
+    TimePoint operator-(const TimePoint& t) const;
+
     //! Add-assign a Duration to this.
     Duration& operator+=(const Duration& a)
     {
@@ -181,6 +184,8 @@ class TimePoint {
 
     //! Calculate the Duration between two TimePoints.
     Duration operator-(const TimePoint& a) const { return Duration(m_t - a.m_t); }
+    //! Subtract a Duration from this TimePoint
+    TimePoint operator-(const Duration& d) const { return TimePoint(m_t - d.m_d); }
     //! Add-assign a Duration to this TimePoint (now + 7 days = next week).
     TimePoint& operator+=(const Duration& d)
     {

core/src/include/Xios.hpp

@@ -163,6 +163,11 @@ class Xios : public Configured<Xios> {
         // CG-based x- and y-dimensions (alt. names x_cg and y_cg)
         { ModelArray::Type::CG, "cg" },
     };
+    std::map<std::string, bool> domainWritten = {
+        { "dim", false },
+        { "vertex", false },
+        { "cg", false },
+    };
     xios::CDomainGroup* getDomainGroup();
     xios::CDomain* getDomain(const std::string& domainId);
     void setupDomains();
@@ -206,6 +211,7 @@ class Xios : public Configured<Xios> {
         { FORCING, forcingFileId },
     };
     void setupFiles();
+    void postprocessOutputFiles();
 
     /* I/O */
     void read(const std::string& fieldId, ModelArray& modelarray);