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"
@@ -155,6 +156,20 @@ if(ENABLE_XIOS)
     target_link_libraries(nextsimlib PUBLIC "${FORTRAN_RUNTIME_LIB}")
 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)
+        target_compile_definitions(nextsimlib PUBLIC USE_OASIS)
+
+        target_link_libraries(nextsimlib PUBLIC ${OASIS_LIBRARIES})
+        target_include_directories(nextsimlib PUBLIC ${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)

Dockerfiles/Dockerfile.devenv

@@ -50,6 +50,8 @@ COPY --from=builder /opt/views /opt/views
 # Create entrypoint script
 RUN { \
       echo '#!/bin/sh'; \
+      echo 'export OMPI_ALLOW_RUN_AS_ROOT=1'; \
+      echo 'export OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1'; \
       echo '. /opt/spack-environment/activate.sh'; \
       echo 'exec "$@"'; \
     } > /entrypoint.sh \

Dockerfiles/Dockerfile.oasis

@@ -0,0 +1,30 @@
+FROM nextsimdg-dev-env:latest
+
+# A temproary Dockerfile with commands that should find their way into Dockerfile.devenv in the end
+
+# Clone the repository
+RUN git clone https://gitlab.com/cerfacs/oasis3-mct.git
+WORKDIR /oasis3-mct/util/make_dir
+
+# Set some ENV's and copy the make file include
+ENV OASIS_COUPLE=/oasis3-mct
+ENV OASIS_ENV=oasis_docker_spack_gcc
+ENV OASIS_DIR=$OASIS_COUPLE/build.$OASIS_ENV
+COPY make.oasis_docker make.$OASIS_ENV
+
+# make!
+RUN make -f TopMakefileOasis3 static-libs shared-libs pyoasis
+
+# Install python packages that pyoasis may need
+RUN apt-get install -y python3-mpi4py python3-isodate
+
+# Append the OASIS library directory to LD_LIBRARY_PATH and append it to PYTHONPATH
+# We assume the last line of the file is 'exec "$@"' and try to retain that
+RUN sed -i~ '$d' /entrypoint.sh
+RUN { \
+    echo 'export LD_LIBRARY_PATH=$LD_LIBRARY_PATH':$OASIS_DIR/lib; \
+    echo 'export PYTHONPATH=$PYTHONPATH':$OASIS_DIR/python; \
+    echo 'exec "$@"'; \
+    } >> /entrypoint.sh
+
+WORKDIR /

Dockerfiles/make.oasis_docker

@@ -0,0 +1,82 @@
+#
+# Include file for OASIS3 Makefile for a Linux docker image using spack and gcc
+#
+###############################################################################
+#
+# CHAN	: communication technique used in OASIS3 (MPI1/MPI2)
+CHAN            = MPI1
+#
+# Paths for libraries, object files and binaries
+#
+# COUPLE        : path for oasis3-mct main directory
+COUPLE          := $(OASIS_COUPLE)
+#
+# ARCHDIR       : directory created when compiling
+ARCHDIR         := $(COUPLE)/build.$(OASIS_ENV)
+#
+# MPI library
+MPIDIR      = /opt/view/
+MPIBIN      = $(MPIDIR)/bin
+MPI_INCLUDE = $(MPIDIR)/include
+MPILIB      = -L$(MPIDIR)/lib -lmpi
+MPIRUN      = $(MPIBIN)/mpirun
+#
+# NETCDF library of the system
+NETCDF_INCLUDE = /opt/view/include
+NETCDF_LIBRARY = -L/opt/view/lib -Wl,-rpath,/opt/view/lib -lnetcdff -lnetcdf -lm
+#
+# Compiling and other commands
+MAKE        = gmake
+F90         = $(MPIBIN)/mpif90 -I$(MPI_INCLUDE)
+F           = $(F90)
+f90         = $(F90)
+f           = $(F90)
+CC          = $(MPIBIN)/mpicc -I$(MPI_INCLUDE)
+LD          = $(MPIBIN)/mpif90 $(MPILIB)
+DYNOPT      = -fPIC
+LDDYNOPT    = -shared
+AR          = ar
+ARFLAGS     = -ruv
+# Fortran libraries for C linking
+F2C_LIBS    = -lmpi_mpifh -lgfortran -lm
+#
+# CPP keys and compiler options
+#  
+CPPDEF    = -Duse_comm_$(CHAN) -D__VERBOSE -DTREAT_OVERLAY
+#
+FCBASEFLAGS := -fallow-argument-mismatch -ffree-line-length-0 
+ifeq ($(OASIS_DEBUG), )
+	ifeq ($(OASIS_OPENMP), )
+		FCBASEFLAGS += -O2
+	else
+		FCBASEFLAGS += -O2 -fopenmp
+	endif
+else
+	ifeq ($(OASIS_OPENMP), )
+                FCBASEFLAGS += -g -fbounds-check
+        else
+                FCBASEFLAGS += -g -fbounds-check -fopenmp
+        endif
+endif
+CCBASEFLAGS  := 
+ifeq ($(OASIS_OPENMP), )
+	CCBASEFLAGS  +=
+else
+	CCBASEFLAGS  += -fopenmp
+endif
+#
+# INC_DIR : includes all *mod for each library
+ INC_DIR = -I$(ARCHDIR)/include
+# FLIBS : for toys when linking in local Makefile
+ FLIBS=${NETCDF_LIBRARY}
+###################
+#
+F90FLAGS    = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE)
+f90FLAGS    = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE)
+FFLAGS      = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE)
+fFLAGS      = $(FCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE)
+CCFLAGS     = $(CCBASEFLAGS) $(INC_DIR) $(CPPDEF) -I$(NETCDF_INCLUDE)
+LDFLAGS     = $(FCBASEFLAGS)
+F2C_LDFLAGS = $(F2C_LIBS)
+#
+#############################################################################

cmake/FindOASIS.cmake

@@ -0,0 +1,31 @@
+# Find oasis
+#
+# Please pass the -DOASIS_DIR variable to cmake (location of the OASIS libraries).
+
+# Search for all library files individually and merge to a single OASIS_LIBRARIES variable for the parent CMakeLists.txt
+find_library(LIB_OASIS NAMES oasis.cbind HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH)
+find_library(LIB_MCT NAMES mct HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH)
+find_library(LIB_MPEU NAMES mpeu HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH)
+find_library(LIB_PSMILE NAMES psmile.MPI1 psmile.MPI2 HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH)
+find_library(LIB_SCRIP NAMES scrip HINTS ${OASIS_DIR}/lib ENV LD_LIBRARY_PATH)
+
+# Only continue if all libraries are found
+if (LIB_OASIS AND LIB_MCT AND LIB_MPEU AND LIB_PSMILE AND LIB_SCRIP)
+    set(OASIS_LIBRARIES ${LIB_OASIS} ${LIB_MCT} ${LIB_MPEU} ${LIB_PSMILE} ${LIB_SCRIP})
+endif ()
+
+# Set OASIS_DIR (not strictly necessary, but nice to have)
+get_filename_component(OASIS_LIB_DIR "${LIB_OASIS}" PATH)
+set(OASIS_DIR "${OASIS_LIB_DIR}/../")
+cmake_path(NORMAL_PATH OASIS_DIR)
+
+# The oasis include file should be where oasis.cbind lives
+find_path(OASIS_INCLUDES NAMES oasis_c.h HINTS ${OASIS_DIR}/include)
+
+# 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_DIR OASIS_LIBRARIES OASIS_INCLUDES)
+
+# message (STATUS "OASIS_LIBRARIES: ${OASIS_LIBRARIES}")
+# message (STATUS "OASIS_INCLUDES: ${OASIS_INCLUDES}")
+# message (STATUS "OASIS_DIR: ${OASIS_DIR}")

core/src/DevStep.cpp

@@ -41,6 +41,8 @@ void DevStep::iterate(const TimestepTime& tst)
     mData.incrementTime(tst.step);
     if ((m_restartPeriod.seconds() > 0) && (mData.time() >= lastOutput + m_restartPeriod)) {
         std::string currentFileName = mData.time().format(m_restartFileName);
+        // Some MPI-IO implementations does not like colon in file names
+        std::replace(currentFileName.begin(), currentFileName.end(), ':', '_');
         pData->writeRestartFile(currentFileName);
         lastOutput = mData.time();
     }

core/src/Model.cpp

@@ -40,6 +40,9 @@ static const std::map<int, std::string> keyMap = {
     { 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
 };
 
 Model::Model()
@@ -122,9 +125,15 @@ void Model::configure()
     // Get the coordinates from the ModelState for persistence
     metadata.extractCoordinates(initialState);
 
+#ifdef USE_OASIS
+    const bool writeOasisGrid = Configured::getConfiguration(keyMap.at(WRITEOASISGRID_KEY), false);
+    metadata.initOasis(writeOasisGrid);
+#endif
+
     modelStep.setData(pData);
     modelStep.setRestartDetails(metadata.restartPeriod, metadata.finalFileName);
     pData.setData(initialState.data);
+    pData.setMetadata(metadata);
 }
 
 Model::HelpMap& Model::getHelpText(HelpMap& map, bool getAll)
@@ -195,6 +204,8 @@ void Model::writeRestartFile()
 {
     auto& metadata = ModelMetadata::getInstance();
     std::string formattedFileName = metadata.time().format(metadata.finalFileName);
+    // Some MPI-IO implementations does not like colon in file names
+    std::replace(formattedFileName.begin(), formattedFileName.end(), ':', '_');
     pData.writeRestartFile(formattedFileName);
 }
 

core/src/ModelArray.cpp

@@ -175,6 +175,27 @@ ModelArray& ModelArray::clampBelow(const ModelArray& minArr)
     return *this;
 }
 
+ModelArray ModelArray::sqrt()
+{
+    ModelArray sqrted = ModelArray(type);
+    sqrted.m_data.array() = m_data.array().sqrt();
+    return sqrted;
+}
+
+ModelArray ModelArray::sin()
+{
+    auto sined = ModelArray(type);
+    sined.m_data.array() = m_data.array().sin();
+    return sined;
+}
+
+ModelArray ModelArray::cos()
+{
+    auto cosed = ModelArray(type);
+    cosed.m_data.array() = m_data.array().cos();
+    return cosed;
+}
+
 void ModelArray::setData(double value)
 {
     resize();

core/src/ModelMetadata.cpp

@@ -25,6 +25,9 @@
 #include <ncGroup.h>
 #include <ncVar.h>
 #endif
+#ifdef USE_OASIS
+#include <oasis_c.h>
+#endif
 
 namespace Nextsim {
 
@@ -191,6 +194,71 @@ ModelState& ModelMetadata::affixCoordinates(ModelState& state) const
     return state;
 }
 
+#ifdef USE_OASIS
+void ModelMetadata::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) {
+        /* This needs to be figured out, but it's not a priority. Grid writing is
+         * not necessary for the type of coupling we'll start with.
+
+        const std::string gridName = "nxts";
+
+        int flag = 1;
+        OASIS_CHECK_ERR(oasis_c_start_grids_writing(&flag));
+
+        OASIS_CHECK_ERR(oasis_c_write_grid(
+            gridName.c_str(), nx, ny, nx_loc, ny_loc, lon, lat, OASISPartitionId));
+        OASIS_CHECK_ERR(oasis_c_write_corner(
+            gridName.c_str(), nx, ny, nx_loc, ny_loc, clo, cla, OASISPartitionId));
+        OASIS_CHECK_ERR(oasis_c_write_area(
+            gridName.c_str(), nx, ny, nx_loc, ny_loc, area, OASISPartitionId));
+        OASIS_CHECK_ERR(oasis_c_write_mask(
+            gridName.c_str(), nx, ny, nx_loc, ny_loc, angle, OASISPartitionId));
+
+        std::string companion = "land area fraction";
+        OASIS_CHECK_ERR(oasis_c_write_frac(
+                            gridName.c_str(), nx, ny, nx_loc, ny_loc, mask, OASISPartitionId),
+            companion.c_str());
+        companion = "land sea mask";
+        OASIS_CHECK_ERR(oasis_c_write_mask(
+                            gridName.c_str(), nx, ny, nx_loc, ny_loc, mask, OASISPartitionId),
+            companion.c_str());
+        */
+    }
+}
+#endif
+
 void ModelMetadata::setTimes(const TimePoint& start, const TimePoint& stop, const Duration& step)
 {
     this->start = start;

core/src/PrognosticData.cpp

@@ -11,6 +11,10 @@
 #include "include/NextsimModule.hpp"
 #include "include/gridNames.hpp"
 
+#ifdef USE_OASIS
+#include <oasis_c.h>
+#endif
+
 namespace Nextsim {
 
 static const std::string checkFieldsKey = "debug.check_fields";
@@ -111,6 +115,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)
 {
     // Prepare everything

core/src/include/Model.hpp

@@ -43,6 +43,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 { return ConfigMap(); };

core/src/include/ModelArray.hpp

@@ -297,6 +297,18 @@ class ModelArray {
      * @param minArr the array of clamp minimum target values.
      */
     ModelArray& clampBelow(const ModelArray& minArr);
+    /*!
+     * @brief Returns the per-element square root of the array
+     */
+    ModelArray sqrt();
+    /*!
+     * @brief Returns the per-element sin of the array
+     */
+    ModelArray sin();
+    /*!
+     * @brief Returns the per-element cos of the array
+     */
+    ModelArray cos();
 
     using MultiDim = std::vector<size_t>;