XIOS implementation of readForcingTimeStatic

core/src/ParaGridIO_Xios.cpp

@@ -100,7 +100,7 @@ ModelState ParaGridIO::getModelState(const std::string& filePath)
 
     // Get all variables in the file and load them into a new ModelState
     const bool readAccess = true;
-    for (std::string fieldId : xiosHandler.configGetFieldNames(readAccess)) {
+    for (std::string fieldId : xiosHandler.configGetInputRestartFieldNames()) {
         ModelArray::Type type = xiosHandler.getFieldType(fieldId);
         if (type == ModelArray::Type::H) {
             HField field(ModelArray::Type::H);
@@ -135,56 +135,40 @@ ModelState ParaGridIO::getModelState(const std::string& filePath)
 ModelState ParaGridIO::readForcingTimeStatic(
     const std::set<std::string>& forcings, const TimePoint& time, const std::string& filePath)
 {
-    // TODO: XIOS implementation
-
     ModelState state;
+    Xios& xiosHandler = Xios::getInstance();
 
-    try {
-        netCDF::NcFile ncFile(filePath, netCDF::NcFile::read);
-
-        // Read the time axis
-        netCDF::NcDim timeDim = ncFile.getDim(timeName);
-        // Read the time variable
-        netCDF::NcVar timeVar = ncFile.getVar(timeName);
-        // Calculate the index of the largest time value on the axis below our target
-        size_t targetTIndex;
-        // Get the time axis as a vector
-        std::vector<double> timeVec(timeDim.getSize());
-        timeVar.getVar(timeVec.data());
-        // Get the index of the largest TimePoint less than the target.
-        targetTIndex = std::find_if(std::begin(timeVec), std::end(timeVec), [time](double t) {
-            return (TimePoint() + Duration(t)) > time;
-        }) - timeVec.begin();
-        // Rather than the first that is greater than, get the last that is less
-        // than or equal to. But don't go out of bounds.
-        if (targetTIndex > 0)
-            --targetTIndex;
+    // Increment the XIOS calendar until it reaches the requested time
+    while (xiosHandler.getCurrentDate() < time) {
+        xiosHandler.incrementCalendar();
+    }
+    TimePoint xiosTime = xiosHandler.getCurrentDate();
+    if (xiosTime > time) {
+        throw std::runtime_error("ParaGridIO::readForcingTimeStatic: requested time point does"
+                                 " not align with the calendar and timestep used by XIOS.");
+    }
+
+    // Get all forcings and load them into a new ModelState
+    const bool readAccess = true;
+    std::set<std::string> forcingFieldIds = xiosHandler.configGetForcingFieldNames();
+    for (const std::string& fieldId : forcings) {
+        if (forcingFieldIds.count(fieldId) == 0 || !xiosHandler.getFieldReadAccess(fieldId)) {
+            throw std::runtime_error("ParaGridIO::readForcingTimeStatic: forcing " + fieldId
+                + " is not configured for reading, but is being read from file.");
+        }
         // ASSUME all forcings are HFields: finite volume fields on the same
         // grid as ice thickness
-        std::vector<size_t> indexArray = { targetTIndex };
-        std::vector<size_t> extentArray = { 1 };
-
-        // Loop over the dimensions of H
-        const std::vector<ModelArray::Dimension>& dimensions
-            = ModelArray::typeDimensions.at(ModelArray::Type::H);
-        for (auto riter = dimensions.rbegin(); riter != dimensions.rend(); ++riter) {
-            indexArray.push_back(0);
-            extentArray.push_back(ModelArray::definedDimensions.at(*riter).localLength);
-        }
-
-        for (const std::string& varName : forcings) {
-            netCDF::NcVar var = ncFile.getVar(varName);
-            state.data[varName] = ModelArray(ModelArray::Type::H);
-            ModelArray& data = state.data.at(varName);
-            data.resize();
+        HField field(ModelArray::Type::H);
+        field.resize();
+        state.merge(ModelState { { { fieldId, field } }, {} });
+    }
 
-            var.getVar(indexArray, extentArray, &data[0]);
+    // Read all forcings from file
+    for (auto& entry : state.data) {
+        const std::string fieldId = entry.first;
+        if (forcings.count(fieldId)) {
+            xiosHandler.read(fieldId, entry.second);
         }
-        ncFile.close();
-    } catch (const netCDF::exceptions::NcException& nce) {
-        std::string ncWhat(nce.what());
-        ncWhat += ": " + filePath;
-        throw std::runtime_error(ncWhat);
     }
     return state;
 }

core/src/Xios.cpp

@@ -44,6 +44,7 @@ namespace Nextsim {
 
 static const std::string xOutputPfx = "XiosOutput";
 static const std::string xInputPfx = "XiosInput";
+static const std::string xForcingPfx = "XiosForcing";
 static const std::map<int, std::string> keyMap
     = { { Xios::ENABLED_KEY, "xios.enable" }, { Xios::STARTTIME_KEY, "model.start" },
           { Xios::STOPTIME_KEY, "model.stop" }, { Xios::TIME_STEP_KEY, "model.time_step" },
@@ -52,7 +53,10 @@ static const std::map<int, std::string> keyMap
           { Xios::OUTPUT_FIELD_NAMES_KEY, xOutputPfx + ".field_names" },
           { Xios::INPUT_RESTARTPERIOD_KEY, xInputPfx + ".period" },
           { Xios::INPUT_RESTARTFILE_KEY, xInputPfx + ".filename" },
-          { Xios::INPUT_FIELD_NAMES_KEY, xInputPfx + ".field_names" } };
+          { Xios::INPUT_FIELD_NAMES_KEY, xInputPfx + ".field_names" },
+          { Xios::FORCING_PERIOD_KEY, xForcingPfx + ".period" },
+          { Xios::FORCING_FILE_KEY, xForcingPfx + ".filename" },
+          { Xios::FORCING_FIELD_NAMES_KEY, xForcingPfx + ".field_names" } };
 
 //! Enable XIOS in the 'config'
 void enableXios()
@@ -79,14 +83,19 @@ Xios::Xios(const std::string contextid, const std::string calendartype)
     if (firstTime) {
         istringstream(Configured::getConfiguration(keyMap.at(INPUT_RESTARTFILE_KEY), std::string()))
             >> inputFilename;
-        if (inputFilename.length() > 0) {
-            inputFileId = ((std::filesystem::path)inputFilename).replace_extension();
-        }
+        inputFileId = ((std::filesystem::path)inputFilename).replace_extension();
         istringstream(
             Configured::getConfiguration(keyMap.at(OUTPUT_RESTARTFILE_KEY), std::string()))
             >> outputFilename;
-        if (outputFilename.length() > 0) {
-            outputFileId = ((std::filesystem::path)outputFilename).replace_extension();
+        outputFileId = ((std::filesystem::path)outputFilename).replace_extension();
+        istringstream(Configured::getConfiguration(keyMap.at(FORCING_FILE_KEY), std::string()))
+            >> forcingFilename;
+        forcingFileId = ((std::filesystem::path)forcingFilename).replace_extension();
+
+        // TODO: Account for separate restart and forcing files (#929)
+        if (forcingFilename.length() > 0 && inputFilename.length() > 0
+            && inputFilename != forcingFilename) {
+            throw std::runtime_error("Xios: Separate forcing and restart files not yet supported");
         }
 
         for (std::string fileId : { inputFileId, outputFileId }) {
@@ -130,7 +139,7 @@ Xios::HelpMap& Xios::getHelpText(HelpMap& map, bool getAll)
             "an ISO8601 duration (P prefix) or number of seconds. A value of zero assumes no "
             "intermediate restart files." },
         { keyMap.at(INPUT_RESTARTFILE_KEY), ConfigType::STRING, {}, "", "",
-            // TODO: Support format "restart%Y-%m-%dT%H:%M:%SZ.nc"
+            // TODO: Support format "restart%Y-%m-%dT%H:%M:%SZ.nc" (#898)
             "The file name to be used for input." },
         { keyMap.at(INPUT_FIELD_NAMES_KEY), ConfigType::STRING, {}, "", "",
             "Comma-separated list of field names to be read from the input file." },
@@ -141,11 +150,23 @@ Xios::HelpMap& Xios::getHelpText(HelpMap& map, bool getAll)
             "duration (P prefix) or number of seconds. A value of zero "
             "ensures no intermediate restart files are written." },
         { keyMap.at(OUTPUT_RESTARTFILE_KEY), ConfigType::STRING, {}, "", "",
-            // TODO: Support format "restart%Y-%m-%dT%H:%M:%SZ.nc"
+            // TODO: Support format "restart%Y-%m-%dT%H:%M:%SZ.nc" (#898)
             "The file name to be used for output." },
         { keyMap.at(OUTPUT_FIELD_NAMES_KEY), ConfigType::STRING, {}, "", "",
             "Comma-separated list of field names to be written to the output file." },
     };
+    map["XiosForcing"] = {
+        { keyMap.at(FORCING_PERIOD_KEY), ConfigType::STRING, {}, "0", "",
+            "The period between forcing file outputs expected in a file to be "
+            "read, formatted as an ISO8601 duration (P prefix) or number of "
+            "seconds. A value of zero assumes no intermediate restart files." },
+        { keyMap.at(FORCING_FILE_KEY), ConfigType::STRING, {}, "", "",
+            // TODO: Support format "restart%Y-%m-%dT%H:%M:%SZ.nc" (#898)
+            "The file name to be used for forcings." },
+        { keyMap.at(FORCING_FIELD_NAMES_KEY), ConfigType::STRING, {}, "", "",
+            "Comma-separated list of field names to be read from the forcings "
+            "file." },
+    };
 
     return map;
 }
@@ -484,11 +505,11 @@ int Xios::getCalendarStep() { return clientCalendar->getCalendar()->getStep(); }
  *
  * @return current calendar date
  */
-std::string Xios::getCurrentDate(const bool isoFormat)
+TimePoint Xios::getCurrentDate()
 {
     cxios_date xiosDate;
     cxios_get_current_date(&xiosDate);
-    return convertXiosDatetimeToString(xiosDate, isoFormat);
+    return TimePoint(convertXiosDatetimeToString(xiosDate, true));
 }
 
 /*!
@@ -1028,35 +1049,79 @@ xios::CField* Xios::getField(const std::string fieldId)
     return field;
 }
 
-// Extract the field_names entry from the XiosInput or XiosOutput section of the config
-std::set<std::string> Xios::configGetFieldNames(const bool reading)
+// Split a string into a set by some delimiter.
+std::set<std::string> str2set(std::string asStr, const char delim = ',')
 {
-    std::string fieldsStr;
-    if (reading) {
-        istringstream(Configured::getConfiguration(keyMap.at(INPUT_FIELD_NAMES_KEY), std::string()))
-            >> fieldsStr;
-    } else {
-        istringstream(
-            Configured::getConfiguration(keyMap.at(OUTPUT_FIELD_NAMES_KEY), std::string()))
-            >> fieldsStr;
-    }
-    std::set<std::string> fieldNames;
-    if (fieldsStr.length() > 0) {
+    std::set<std::string> asSet;
+    if (asStr.length() > 0) {
         const char delim = ',';
-        std::istringstream iss(fieldsStr);
+        std::istringstream iss(asStr);
         std::string item;
         while (std::getline(iss, item, delim)) {
-            fieldNames.insert(item);
+            asSet.insert(item);
         }
     }
-    return fieldNames;
+    return asSet;
+}
+
+// Extract the field_names entry from the XiosInput section of the config.
+std::set<std::string> Xios::configGetInputRestartFieldNames()
+{
+    std::string fieldsStr;
+    istringstream(Configured::getConfiguration(keyMap.at(INPUT_FIELD_NAMES_KEY), std::string()))
+        >> fieldsStr;
+    return str2set(fieldsStr);
+}
+
+/*!
+ * Extract the field_names entry from the XiosForcing section of the config.
+ */
+std::set<std::string> Xios::configGetForcingFieldNames()
+{
+    std::string fieldsStr;
+    istringstream(Configured::getConfiguration(keyMap.at(FORCING_FIELD_NAMES_KEY), std::string()))
+        >> fieldsStr;
+    return str2set(fieldsStr);
+}
+
+// Extract the field_names entry from the XiosInput and XiosForcing sections of the config.
+std::set<std::string> Xios::configGetInputFieldNames()
+{
+    std::set<std::string> restartFieldNames = configGetInputRestartFieldNames();
+    std::set<std::string> forcingFieldNames = configGetForcingFieldNames();
+    restartFieldNames.insert(forcingFieldNames.begin(), forcingFieldNames.end());
+    return restartFieldNames;
+}
+
+// Extract the field_names entry from the XiosOutput section of the config.
+std::set<std::string> Xios::configGetOutputFieldNames()
+{
+    // TODO: Account for diagnostics (#917)
+    std::string fieldsStr;
+    istringstream(Configured::getConfiguration(keyMap.at(OUTPUT_FIELD_NAMES_KEY), std::string()))
+        >> fieldsStr;
+    return str2set(fieldsStr);
+}
+
+/*!
+ * Extract the field_names entry from the XIOS config.
+ *
+ * @param readAccess true if the fields are to be read, false if written
+ */
+std::set<std::string> Xios::configGetFieldNames(const bool readAccess)
+{
+    if (readAccess) {
+        return configGetInputFieldNames();
+    } else {
+        return configGetOutputFieldNames();
+    }
 }
 
 // Check whether a fieldId exists in a string of field names separated by commas, as determined by
 // the map key
-bool Xios::configCheckField(const std::string fieldId, const bool reading)
+bool Xios::configCheckField(const std::string fieldId, const bool readAccess)
 {
-    std::set<std::string> fieldNames = configGetFieldNames(reading);
+    std::set<std::string> fieldNames = configGetFieldNames(readAccess);
     return fieldNames.find(fieldId) != fieldNames.end();
 }
 
@@ -1300,7 +1365,9 @@ void Xios::createFile(const std::string fileId)
 
     // Determine whether the file is configured for reading or writing
     bool readAccess = ((inputFileId.length() > 0) && (inputFileId == fileId));
+    // TODO: Account for separate restart and forcing files (#929)
     bool writeAccess = ((outputFileId.length() > 0) && (outputFileId == fileId));
+    // TODO: Account for diagnostics (#917)
 
     // Check that the filename is not in both the XiosOutput and XiosInput config sections
     if (readAccess && writeAccess) {
@@ -1337,10 +1404,12 @@ void Xios::createFile(const std::string fileId)
         istringstream(
             Configured::getConfiguration(keyMap.at(INPUT_RESTARTPERIOD_KEY), std::string()))
             >> periodStr;
+        // TODO: Account for forcing period being different to restart period (#929)
     } else {
         istringstream(
             Configured::getConfiguration(keyMap.at(OUTPUT_RESTARTPERIOD_KEY), std::string()))
             >> periodStr;
+        // TODO: Account for diagnostics (#917)
     }
     if (periodStr.length() == 0 || periodStr == "0") {
         setFileOutputFreq(fileId, stopTime - startTime);
@@ -1349,7 +1418,8 @@ void Xios::createFile(const std::string fileId)
     }
 
     // Create all fields found in the config based off the field names found in the
-    // XiosInput.field_names or XiosOutput.field_names entries in the config.
+    // XiosInput.field_names, XiosOutput.field_names, or XiosForcing.field_names entries in the
+    // config.
     for (std::string fieldId : configGetFieldNames(readAccess)) {
         createField(fieldId);
         fileAddField(fileId, fieldId);

core/src/include/Xios.hpp

@@ -86,7 +86,7 @@ class Xios : public Configured<Xios> {
     TimePoint getCalendarStart();
     Duration getCalendarTimestep();
     int getCalendarStep();
-    std::string getCurrentDate(const bool isoFormat = true);
+    TimePoint getCurrentDate();
 
     /* Axis */
     void createAxis(const std::string axisId);
@@ -112,6 +112,7 @@ class Xios : public Configured<Xios> {
     std::string getFieldGridRef(const std::string fieldId);
     bool getFieldReadAccess(const std::string fieldId);
     Duration getFieldFreqOffset(const std::string fieldId);
+    std::set<std::string> configGetForcingFieldNames();
     ModelArray::Type getFieldType(const std::string fieldId);
     void setFieldType(const std::string fieldId, ModelArray::Type type);
 
@@ -144,6 +145,9 @@ class Xios : public Configured<Xios> {
         INPUT_RESTARTPERIOD_KEY,
         INPUT_RESTARTFILE_KEY,
         INPUT_FIELD_NAMES_KEY,
+        FORCING_PERIOD_KEY,
+        FORCING_FILE_KEY,
+        FORCING_FIELD_NAMES_KEY,
     };
 
     /* Length of C-strings passed to XIOS */
@@ -199,8 +203,11 @@ class Xios : public Configured<Xios> {
     xios::CFieldGroup* getFieldGroup();
     xios::CField* getField(const std::string fieldId);
     void setFieldReadAccess(const std::string fieldId, const bool readAccess);
-    std::set<std::string> configGetFieldNames(const bool reading);
-    bool configCheckField(const std::string fieldId, const bool reading);
+    std::set<std::string> configGetInputRestartFieldNames();
+    std::set<std::string> configGetInputFieldNames();
+    std::set<std::string> configGetOutputFieldNames();
+    std::set<std::string> configGetFieldNames(const bool readAccess);
+    bool configCheckField(const std::string fieldId, const bool readAccess);
     std::map<std::string, ModelArray::Type> fieldTypes;
 
     /* Grid */
@@ -220,6 +227,8 @@ class Xios : public Configured<Xios> {
     std::string inputFileId;
     std::string outputFilename;
     std::string outputFileId;
+    std::string forcingFilename;
+    std::string forcingFileId;
 
     /* I/O */
     void write(const std::string fieldId, ModelArray& modelarray);

core/test/XiosCalendar_test.cpp

@@ -67,12 +67,12 @@ MPI_TEST_CASE("TestXiosCalendar", 1)
 
     // --- Tests for getCurrentDate method
     REQUIRE(xiosHandler.getCalendarStep() == 0);
-    REQUIRE(xiosHandler.getCurrentDate() == "2023-03-17T17:11:00Z");
-    REQUIRE(xiosHandler.getCurrentDate(false) == "2023-03-17 17:11:00");
+    REQUIRE(xiosHandler.getCurrentDate().format() == "2023-03-17T17:11:00Z");
 
     // -- Tests that the timestep is set up correctly
     xiosHandler.setCalendarStep(1);
-    REQUIRE(xiosHandler.getCurrentDate() == "2023-03-17T18:41:00Z");
+    REQUIRE(xiosHandler.getCalendarStep() == 1);
+    REQUIRE(xiosHandler.getCurrentDate().format() == "2023-03-17T18:41:00Z");
 
     xiosHandler.context_finalize();
     Finalizer::finalize();