Issue459 oasis ocean coupling

core/src/Model.cpp

@@ -1,6 +1,6 @@
 /*!
  * @file Model.cpp
- * @date 12 Aug 2021
+ * @date 09 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Kacper Kornet <kk562@cam.ac.uk>
  */
@@ -47,6 +47,9 @@ static std::map<int, std::string> modelConfigKeyMap = {
     { Model::MISSINGVALUE_KEY, "model.missing_value" },
     { Model::RESTARTPERIOD_KEY, "model.restart_period" },
     { Model::RESTARTOUTFILE_KEY, "model.restart_file" },
+#ifdef USE_OASIS
+    { Model::WRITEOASISGRID_KEY, "oasis.write_grid" },
+#endif
 };
 
 // Merge the configuration from ModelConfig into the Model keyMap.
@@ -136,6 +139,11 @@ void Model::configure()
     // Get the coordinates from the ModelState for persistence
     m_etadata.extractCoordinates(initialState);
 
+#ifdef USE_OASIS
+    const bool writeOasisGrid = Configured::getConfiguration(keyMap.at(WRITEOASISGRID_KEY), false);
+    m_etadata.initOasis(writeOasisGrid);
+#endif
+
     modelStep.setData(pData);
     modelStep.setMetadata(m_etadata);
     modelStep.setRestartDetails(restartPeriod, finalFileName);

core/src/PrognosticData.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file PrognosticData.cpp
  *
- * @date 1 Jul 2024
+ * @date 09 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -12,6 +12,10 @@
 #include "include/Module.hpp"
 #include "include/gridNames.hpp"
 
+#ifdef USE_OASIS
+#include <oasis_c.h>
+#endif
+
 namespace Nextsim {
 
 PrognosticData::PrognosticData()
@@ -78,6 +82,11 @@ void PrognosticData::setMetadata(const Nextsim::ModelMetadata& metadata)
 {
     pAtmBdy->setMetadata(metadata);
     pOcnBdy->setMetadata(metadata);
+
+#ifdef USE_OASIS
+    // OASIS finalising definition - can only be called once
+    OASIS_CHECK_ERR(oasis_c_enddef());
+#endif
 }
 
 void PrognosticData::update(const TimestepTime& tst)

core/src/include/Model.hpp

@@ -1,6 +1,6 @@
 /*!
  * @file Model.hpp
- * @date 12 Aug 2021
+ * @date 09 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Kacper Kornet <kk562@cam.ac.uk>
  */
@@ -49,6 +49,9 @@ class Model : public Configured<Model> {
         // Other Model configuration keys, not to be written to the restart file.
         RESTARTPERIOD_KEY,
         RESTARTOUTFILE_KEY,
+#ifdef USE_OASIS
+        WRITEOASISGRID_KEY,
+#endif
     };
 
     ConfigMap getConfig() const;

core/src/include/ModelMetadata.hpp

@@ -1,7 +1,7 @@
 /*!
  * @file ModelMetadata.hpp
  *
- * @date Jun 29, 2022
+ * @date 09 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -18,6 +18,9 @@
 #ifdef USE_MPI
 #include <mpi.h>
 #endif
+#ifdef USE_OASIS
+#include <oasis_c.h>
+#endif
 
 namespace Nextsim {
 
@@ -42,6 +45,10 @@ class ModelMetadata {
     ModelMetadata() = default;
 #endif
 
+#ifdef USE_OASIS
+    int OASISPartitionId;
+#endif
+
     /*!
      * @brief Sets the initial or current model time
      *
@@ -99,6 +106,44 @@ class ModelMetadata {
     int localCornerX, localCornerY, localExtentX, localExtentY, globalExtentX, globalExtentY;
 #endif
 
+#ifdef USE_OASIS
+    void initOasis(const bool writeOasisGrid)
+    {
+        // Set the partitioning
+        /* From the manual: "[ig_paral is a] vector of integers describing the local grid partition
+         * in the global index space; has a different expression depending on the type of the
+         * partition; in OASIS3-MCT, 5 types of partition are supported: Serial (no partition),
+         * Apple, Box, Orange, and Points" - it looks like we should use "Box", so partInfo[0] = 2
+         * (aka. ig_paral).
+         *
+         * #Box partition#
+         * Each partition is a rectangular region of the global domain, described by the global
+         * offset of its upper left corner, and its local extents in the X and Y dimensions. The
+         * global extent in the X dimension must also be given. In this case, we have ig_paral(1:5):
+         *  - ig_paral(1) = 2 (indicates a Box partition)
+         *  - ig_paral(2) = the upper left corner global offset
+         *  - ig paral(3) = the local extent in x
+         *  - ig_paral(4) = the local extent in y
+         *  - ig_paral(5) = the global extent in x.
+         *
+         * metdatata contains: localCornerX, localCornerY, localExtentX, localExtentY,
+         * globalExtentX, globalExtentY;
+         */
+        // TODO: The contents of metadata is not certain!
+        const int offset = globalExtentX * localCornerY + localCornerX;
+        const std::vector<int> partInfo
+            = { OASIS_Box, offset, localExtentX, localExtentY, globalExtentX };
+
+        const int globalSize = globalExtentX * globalExtentY;
+        const std::string compName = "nextsim"; // Not useful for any setups we have in mind
+        OASIS_CHECK_ERR(oasis_c_def_partition(
+            &OASISPartitionId, OASIS_Box_Params, &partInfo[0], globalSize, compName.c_str()));
+
+        // TODO: Writing out grid information should be possible, but optional
+        if (writeOasisGrid) { }
+    }
+#endif
+
 private:
     TimePoint m_time;
     ConfigMap m_config;

core/src/main.cpp

@@ -1,6 +1,6 @@
 /*!
  * @file main.cpp
- * @date 11 Aug 2021
+ * @date 09 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Kacper Kornet <kk562@cam.ac.uk>
  */
@@ -9,6 +9,9 @@
 #ifdef USE_MPI
 #include <mpi.h>
 #endif
+#ifdef USE_OASIS
+#include <oasis_c.h>
+#endif
 
 #include "include/CommandLineParser.hpp"
 #include "include/ConfigurationHelpPrinter.hpp"
@@ -20,7 +23,18 @@
 int main(int argc, char* argv[])
 {
 #ifdef USE_MPI
+    MPI_Comm modelCommunicator;
     MPI_Init(&argc, &argv);
+#ifdef USE_OASIS
+    /* We must call these oasis routines before any MPI communication takes place, to make sure we
+     * have the right communicator, i.e. modelCommunictor and not MPI_COMM_WORLD. */
+    int compID; // Not actually used. Only useful for debugging
+    const std::string compName = "nextsim"; // Not useful for any setups we have in mind
+    OASIS_CHECK_ERR(oasis_c_init_comp(&compID, compName.c_str(), OASIS_COUPLED));
+    OASIS_CHECK_ERR(oasis_c_get_localcomm(&modelCommunicator));
+#else
+    modelCommunicator = MPI_COMM_WORLD;
+#endif // USE_OASIS
 #endif // USE_MPI
 
     // Pass the command line to Configurator to handle
@@ -49,7 +63,7 @@ int main(int argc, char* argv[])
     } else {
         // Construct the Model
 #ifdef USE_MPI
-        Nextsim::Model model(MPI_COMM_WORLD);
+        Nextsim::Model model(modelCommunicator);
 #else
         Nextsim::Model model;
 #endif
@@ -59,6 +73,9 @@ int main(int argc, char* argv[])
         model.run();
     }
 #ifdef USE_MPI
+#ifdef USE_OASIS
+    OASIS_CHECK_ERR(oasis_c_terminate());
+#endif
     MPI_Finalize();
 #endif
 

physics/src/include/OASISCoupled.hpp

@@ -1,14 +1,13 @@
 /*!
  * @file OASISCoupled.hpp
  *
- * @date 29 Aug 2024
+ * @date 09 Sep 2024
  * @author Einar Ólason <einar.olason@nersc.no>
  */
 
 #ifndef OASISCOUPLED_HPP
 #define OASISCOUPLED_HPP
 
-#include "include/ModelMetadata.hpp"
 #ifdef USE_OASIS
 #include <oasis_c.h>
 #endif
@@ -17,82 +16,20 @@ namespace Nextsim {
 
 class OASISCoupled {
 public:
-#ifdef USE_OASIS
-    ~OASISCoupled() { OASIS_CHECK_ERR(!oasis_c_terminate()); }
-#else
-    ~OASISCoupled() { }
-    const std::string OASISError
-        = "Cannot access OASIS interface, as OASIS support was not activated at compile time.\n";
-#endif
-
     virtual std::string getName() const { return "OASISCoupled"; }
 
 #ifdef USE_OASIS
-    int partitionID, OASISTime;
-    void setMetadata(const ModelMetadata& metadata)
-    {
-        // Set the "OASIS time" (seconds since start) to zero
-        // TODO: We need to figure out if this is ok to do when restarting
-        OASISTime = 0;
-
-        // Set the communicators
-        // Some of these three may need to be global to the class
-        int compID;
-        MPI_Comm localComm;
-        MPI_Comm coupledComm;
-
-        // compName could be configurable, but I don't see the need
-        const std::string compName = "nextsim";
-        OASIS_CHECK_ERR(oasis_c_init_comp_with_comm(
-            &compID, compName.c_str(), OASIS_COUPLED, metadata.mpiComm));
-
-        // Eric didn't include those, but it seems they should be necessary
-        OASIS_CHECK_ERR(oasis_c_get_localcomm(&localComm));
-        OASIS_CHECK_ERR(oasis_c_create_couplcomm(OASIS_COUPLED, localComm, &coupledComm));
-
-        // Set the partitioning
-        /* From the manual: "[ig_paral is a] vector of integers describing the local grid partition
-         * in the global index space; has a different expression depending on the type of the
-         * partition; in OASIS3-MCT, 5 types of partition are supported: Serial (no partition),
-         * Apple, Box, Orange, and Points" - it looks like we should use "Box", so partInfo[0] = 2
-         * (aka. ig_paral).
-         *
-         * #Box partition#
-         * Each partition is a rectangular region of the global domain, described by the global
-         * offset of its upper left corner, and its local extents in the X and Y dimensions. The
-         * global extent in the X dimension must also be given. In this case, we have ig_paral(1:5):
-         *  - ig_paral(1) = 2 (indicates a Box partition)
-         *  - ig_paral(2) = the upper left corner global offset
-         *  - ig paral(3) = the local extent in x
-         *  - ig_paral(4) = the local extent in y
-         *  - ig_paral(5) = the global extent in x.
-         *
-         * metdatata contains: localCornerX, localCornerY, localExtentX, localExtentY,
-         * globalExtentX, globalExtentY;
-         */
-        // TODO: The contents of metadata is not certain!
-        const int offset = metadata.globalExtentX * metadata.localCornerY + metadata.localCornerX;
-        const std::vector<int> partInfo = { OASIS_Box, offset, metadata.localExtentX,
-            metadata.localExtentY, metadata.globalExtentX };
-
-        const int globalSize = metadata.globalExtentX * metadata.globalExtentY;
-        OASIS_CHECK_ERR(oasis_c_def_partition(
-            &partitionID, OASIS_Box_Params, &partInfo[0], globalSize, compName.c_str()));
-
-        // TODO: Writing out grid information should be possible, but optional
+    int OASISTime;
 
-        // def_var and end_def calls are called by the child class
-    }
+    // Set the "OASIS time" (seconds since start) to zero
+    OASISCoupled() { OASISTime = 0; }
 
     // Increment the "OASIS" time by the number of seconds in the time step
+    // Could be any time unit
     // Must be called at the end of the child class' update or updateAfter call.
     void updateOASISTime(const TimestepTime& tst) { OASISTime += tst.step.seconds(); }
 #else
-    virtual void setMetadata(const ModelMetadata& metadata)
-    {
-        std::string message = __func__ + std::string(": OASIS support not compiled in.\n");
-        throw std::runtime_error(message);
-    }
+    const std::string OASISError = std::string(": OASIS support not compiled in.\n");
 #endif
 };
 

physics/src/modules/OceanBoundaryModule/OASISCoupledOcean.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file OASISCoupledOcean.cpp
  *
- * @date 27 Aug 2024
+ * @date 09 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -15,30 +15,26 @@ namespace Nextsim {
 
 void OASISCoupledOcean::setMetadata(const ModelMetadata& metadata)
 {
-    // The parent class knows how to set the communicators and partitions
-    OASISCoupled::setMetadata(metadata);
-
 #ifdef USE_OASIS
     // OASIS defining variable
 
     /* Populate the couplingId map with the id string and number pair. We need to do this seperately
      * for the input (get) and output (put) variables. */
+    // TODO: coplingID should be a map of <std::string, std::pair> where the pair is idNumber and
+    // pointer to the ModelArray
     for (std::string idString : cplStringsIn) {
         int idNumber;
-        OASIS_CHECK_ERR(oasis_c_def_var(
-            &idNumber, idString.c_str(), partitionID, bundleSize, OASIS_IN, OASIS_DOUBLE));
+        OASIS_CHECK_ERR(oasis_c_def_var(&idNumber, idString.c_str(), metadata.OASISPartitionId,
+            bundleSize, OASIS_IN, OASIS_DOUBLE));
         couplingId[idString] = idNumber;
     }
 
     for (std::string idString : cplStringsOut) {
         int idNumber;
-        OASIS_CHECK_ERR(oasis_c_def_var(
-            &idNumber, idString.c_str(), partitionID, bundleSize, OASIS_OUT, OASIS_DOUBLE));
+        OASIS_CHECK_ERR(oasis_c_def_var(&idNumber, idString.c_str(), metadata.OASISPartitionId,
+            bundleSize, OASIS_OUT, OASIS_DOUBLE));
         couplingId[idString] = idNumber;
     }
-
-    // OASIS finalising definition
-    OASIS_CHECK_ERR(oasis_c_enddef());
 #else
     throw std::runtime_error(std::string(__func__) + ": " + OASISError);
 #endif
@@ -108,7 +104,7 @@ void OASISCoupledOcean::updateAfter(const TimestepTime& tst)
         OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo));
 
     OASIS_CHECK_ERR(oasis_c_put(couplingId.at(EMPKey), OASISTime, dimension0, dimension1, 1,
-        OASIS_DOUBLE, OASIS_COL_MAJOR, &emp[0], OASIS_No_Restart, &kinfo));
+        OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo));
 
     OASIS_CHECK_ERR(oasis_c_put(couplingId.at(QSWKey), OASISTime, dimension0, dimension1, 1,
         OASIS_DOUBLE, OASIS_COL_MAJOR, &dummy[0], OASIS_No_Restart, &kinfo));