Issue459 oasis ocean coupling

CMakeLists.txt

@@ -9,6 +9,7 @@ project(nextsim_dg)
 option(WITH_THREADS "Build with support for openmp" OFF)
 option(ENABLE_MPI "Enable distributed-memory parallelization with MPI" OFF)
 option(ENABLE_XIOS "Enable XIOS library for IO" OFF)
+option(ENABLE_OASIS "Enable the OASIS interface" OFF)
 option(BUILD_TESTS "Build the tests" ON)
 set(DynamicsType
     "DG2"
@@ -145,13 +146,29 @@ if(ENABLE_XIOS)
     )
 endif()
 
+# Do we want to include the OASIS interface? It'll only work if we have MPI
+if (ENABLE_OASIS)
+    if (ENABLE_MPI)
+        message(STATUS "Building with OASIS support")
+        find_package(OASIS REQUIRED)
+        pkg_search_module(NETCDF REQUIRED netcdf-fortran)
+        target_compile_definitions(nextsimlib PUBLIC USE_OASIS)
+
+        target_link_directories(nextsimlib PUBLIC ${OASIS_LIBRARIES})
+        target_link_libraries(nextsimlib PUBLIC oasis.cbind mct mpeu psmile.MPI1 scrip netcdff)
+        target_include_directories(nextsimlib PRIVATE ${OASIS_INCLUDES})
+    else ()
+        message(FATAL_ERROR "Cannot build with OASIS support, because MPI is not enabled" .)
+    endif ()
+endif ()
+
 # OpenMP
 if(WITH_THREADS)
     find_package(OpenMP REQUIRED)
     target_link_libraries(nextsimlib PUBLIC OpenMP::OpenMP_CXX)
 endif()
 
-# Set an empty list of sources
+# Set an empty list of sources netcdf netcdfff
 set(NextsimSources "")
 # Set an empty list of include directories
 set(NextsimIncludeDirs "")

cmake/Findoasis.cmake

@@ -0,0 +1,21 @@
+# Find oasis
+#
+# Please pass the -DOASIS_DIR variable to cmake (location of the OASIS libraries).
+
+find_library (OASIS_LIBRARIES NAMES mct mpeu oasis psmile scrip HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH)
+
+get_filename_component (OASIS_LIBRARIES "${OASIS_LIBRARIES}" PATH)
+
+set (OASIS_DIR "${OASIS_LIBRARIES}/../")
+cmake_path(NORMAL_PATH OASIS_DIR)
+
+find_path (OASIS_INCLUDES NAMES oasis_c.h HINTS ${OASIS_DIR}/inc)
+
+# handle the QUIETLY and REQUIRED arguments and set OASIS_FOUND to TRUE if all listed variables are TRUE
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args(OASIS DEFAULT_MSG
+        OASIS_LIBRARIES OASIS_INCLUDES)
+
+# message (STATUS "OASIS_LIBRARIES: ${OASIS_LIBRARIES}")
+# message (STATUS "OASIS_INCLUDES: ${OASIS_INCLUDES}")
+# message (STATUS "OASIS_DIR: ${OASIS_DIR}")

core/src/Model.cpp

@@ -1,6 +1,6 @@
 /*!
  * @file Model.cpp
- * @date 12 Aug 2021
+ * @date 09 Sep 2024
  * @author Tim Spain <[email protected]>
  * @author Kacper Kornet <[email protected]>
  */
@@ -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,10 +139,16 @@ 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);
     pData.setData(initialState.data);
+    pData.setMetadata(m_etadata);
 }
 
 ConfigMap Model::getConfig() const

core/src/PrognosticData.cpp

@@ -1,7 +1,7 @@
 /*!
  * @file PrognosticData.cpp
  *
- * @date 1 Jul 2024
+ * @date 09 Sep 2024
  * @author Tim Spain <[email protected]>
  * @author Einar Ólason <[email protected]>
  */
@@ -12,6 +12,10 @@
 #include "include/Module.hpp"
 #include "include/gridNames.hpp"
 
+#ifdef USE_OASIS
+#include <oasis_c.h>
+#endif
+
 namespace Nextsim {
 
 PrognosticData::PrognosticData()
@@ -74,6 +78,17 @@ void PrognosticData::setData(const ModelState::DataMap& ms)
     iceGrowth.setData(ms);
 }
 
+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)
 {
     ModelArrayRef<Shared::T_ICE, RW> ticeUpd(getStore());

core/src/include/Model.hpp

@@ -1,6 +1,6 @@
 /*!
  * @file Model.hpp
- * @date 12 Aug 2021
+ * @date 09 Sep 2024
  * @author Tim Spain <[email protected]>
  * @author Kacper Kornet <[email protected]>
  */
@@ -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 <[email protected]>
  */
 
@@ -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/include/PrognosticData.hpp

@@ -38,6 +38,8 @@ class PrognosticData : public ModelComponent, public Configured<PrognosticData>
     ModelState getState(const OutputLevel& lvl) const override { return getState(); }
     ModelState getStateRecursive(const OutputSpec& os) const override;
 
+    void setMetadata(const ModelMetadata& metadata);
+
     static HelpMap& getHelpText(HelpMap& map, bool getAll);
     static HelpMap& getHelpRecursive(HelpMap& map, bool getAll);
 

core/src/main.cpp

@@ -1,6 +1,6 @@
 /*!
  * @file main.cpp
- * @date 11 Aug 2021
+ * @date 09 Sep 2024
  * @author Tim Spain <[email protected]>
  * @author Kacper Kornet <[email protected]>
  */
@@ -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

@@ -0,0 +1,38 @@
+/*!
+ * @file OASISCoupled.hpp
+ *
+ * @date 09 Sep 2024
+ * @author Einar Ólason <[email protected]>
+ */
+
+#ifndef OASISCOUPLED_HPP
+#define OASISCOUPLED_HPP
+
+#ifdef USE_OASIS
+#include <oasis_c.h>
+#endif
+
+namespace Nextsim {
+
+class OASISCoupled {
+public:
+    virtual std::string getName() const { return "OASISCoupled"; }
+
+#ifdef USE_OASIS
+    int OASISTime;
+
+    // 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
+    const std::string OASISError = std::string(": OASIS support not compiled in.\n");
+#endif
+};
+
+}
+
+#endif // OASISCOUPLED_HPP