Add halo exchange logic

CMakeLists.txt

@@ -37,6 +37,7 @@ else()
         set(NSDG_NetCDF_Library "netcdf-cxx4")
     else()
         set(NSDG_NetCDF_Library "netcdf_c++4")
+        set(NSDG_NetCDF_Library "netcdf-cxx4")
     endif()
 endif()
 target_include_directories(nextsimlib PUBLIC "${netCDF_INCLUDE_DIR}")

core/src/ModelArraySlice.cpp

@@ -39,7 +39,7 @@ ModelArraySlice& ModelArraySlice::operator=(ModelArraySlice& other)
             throw std::out_of_range("ModelArraySlice shape mismatch");
     }
 
-    copySliceWithIters(other.data, otherIter, data, thisIter);
+    copySliceWithIters(other.data.m_data, otherIter, data.m_data, thisIter);
     return *this;
 }
 
@@ -55,24 +55,6 @@ ModelArray& ModelArraySlice::copyToModelArray(ModelArray& target) const
     return target;
 }
 
-ModelArray::DataType& ModelArraySlice::copyToDataSlice(
-    ModelArray::DataType& target, SliceIter& targetIter) const
-{
-    SliceIter iter(slice, data.dimensions());
-
-    copySliceWithItersData(data.m_data, iter, target, targetIter);
-    return target;
-}
-
-ModelArraySlice& ModelArraySlice::copyFromDataSlice(
-    const ModelArray::DataType& source, SliceIter& sourceIter)
-{
-    SliceIter iter(slice, data.dimensions());
-
-    copySliceWithItersData(source, sourceIter, data.m_data, iter);
-    return *this;
-}
-
 ModelArraySlice::iterator ModelArraySlice::begin()
 {
     iterator iter(*this);
@@ -127,31 +109,9 @@ void ModelArraySlice::copyBetweenMAandMASlice(
 
     // Copy the data in the correct direction
     if (toSlice) {
-        copySliceWithIters(ma, maIter, mas.data, masIter);
+        copySliceWithIters(ma.m_data, maIter, mas.data.m_data, masIter);
     } else {
-        copySliceWithIters(mas.data, masIter, ma, maIter);
-    }
-}
-
-void ModelArraySlice::copySliceWithIters(
-    ModelArray& source, SliceIter& sourceIter, ModelArray& target, SliceIter targetIter)
-{
-    copySliceWithItersData(source.m_data, sourceIter, target.m_data, targetIter);
-}
-
-void ModelArraySlice::copySliceWithItersData(const ModelArray::DataType& source,
-    SliceIter& sourceIter, ModelArray::DataType& target, SliceIter targetIter)
-{
-    const size_t targetNEl = targetIter.nElements(0);
-    const size_t sourceNEl = sourceIter.nElements(0);
-
-    while (!targetIter.isEnd() && !sourceIter.isEnd()) {
-        const size_t targetIndex = targetIter.index();
-        const size_t sourceIndex = sourceIter.index();
-        target(Eigen::seqN(targetIndex, targetNEl, targetIter.step(0)), Eigen::all)
-            = source(Eigen::seqN(sourceIndex, sourceNEl, sourceIter.step(0)), Eigen::all);
-        targetIter.incrementDim(1);
-        sourceIter.incrementDim(1);
+        copySliceWithIters(mas.data.m_data, masIter, ma.m_data, maIter);
     }
 }
 

core/src/ModelMetadata.cpp

@@ -10,6 +10,9 @@
 #include "include/IStructure.hpp"
 #include "include/NextsimModule.hpp"
 #include "include/gridNames.hpp"
+#include "mpi.h"
+#include <cstddef>
+#include <vector>
 
 #ifdef USE_MPI
 #include <ncDim.h>
@@ -39,6 +42,75 @@ void ModelMetadata::setMpiMetadata(MPI_Comm comm)
     MPI_Comm_rank(mpiComm, &mpiMyRank);
 }
 
+void ModelMetadata::readNeighbourData(netCDF::NcFile& ncFile)
+{
+    netCDF::NcGroup neighbourGroup(ncFile.getGroup(neighbourName));
+    std::string varName {};
+    for (auto edge : edges) {
+        size_t nStart {}; // start point in metadata arrays
+        size_t count {}; // number of elements to read from metadata arrays
+        std::vector<int> numNeighbours = std::vector<int>(mpiSize, 0);
+        std::vector<int> offsets = std::vector<int>(mpiSize, 0);
+
+        // non-periodic neighbours
+        varName = edgeNames[edge] + "_neighbours";
+        neighbourGroup.getVar(varName).getVar(
+            { 0 }, { static_cast<size_t>(mpiSize) }, &numNeighbours[0]);
+
+        // compute start index for each process
+        MPI_Exscan(&numNeighbours[mpiMyRank], &nStart, 1, MPI_INT, MPI_SUM, mpiComm);
+        // how many elements to read for each process
+        count = numNeighbours[mpiMyRank];
+
+        if (count) {
+            // initialize neighbour info to zero and correct size
+            neighbourRanks[edge].resize(count, 0);
+            neighbourExtents[edge].resize(count, 0);
+            neighbourHaloStarts[edge].resize(count, 0);
+
+            varName = edgeNames[edge] + "_neighbour_ids";
+            neighbourGroup.getVar(varName).getVar({ nStart }, { count }, &neighbourRanks[edge][0]);
+
+            varName = edgeNames[edge] + "_neighbour_halos";
+            neighbourGroup.getVar(varName).getVar(
+                { nStart }, { count }, &neighbourExtents[edge][0]);
+
+            varName = edgeNames[edge] + "_neighbour_halo_starts";
+            neighbourGroup.getVar(varName).getVar(
+                { nStart }, { count }, &neighbourHaloStarts[edge][0]);
+        }
+
+        // periodic neighbours
+        varName = edgeNames[edge] + "_neighbours_periodic";
+        neighbourGroup.getVar(varName).getVar(
+            { 0 }, { static_cast<size_t>(mpiSize) }, &numNeighbours[0]);
+
+        // compute start index for each process
+        MPI_Exscan(&numNeighbours[mpiMyRank], &nStart, 1, MPI_INT, MPI_SUM, mpiComm);
+        // how many elements to read for each process
+        count = numNeighbours[mpiMyRank];
+
+        if (count) {
+            // initialize neighbour info to zero and correct size
+            neighbourRanksPeriodic[edge].resize(count, 0);
+            neighbourExtentsPeriodic[edge].resize(count, 0);
+            neighbourHaloStartsPeriodic[edge].resize(count, 0);
+
+            varName = edgeNames[edge] + "_neighbour_ids_periodic";
+            neighbourGroup.getVar(varName).getVar(
+                { nStart }, { count }, &neighbourRanksPeriodic[edge][0]);
+
+            varName = edgeNames[edge] + "_neighbour_halos_periodic";
+            neighbourGroup.getVar(varName).getVar(
+                { nStart }, { count }, &neighbourExtentsPeriodic[edge][0]);
+
+            varName = edgeNames[edge] + "_neighbour_halo_starts_periodic";
+            neighbourGroup.getVar(varName).getVar(
+                { nStart }, { count }, &neighbourHaloStartsPeriodic[edge][0]);
+        }
+    }
+}
+
 void ModelMetadata::getPartitionMetadata(std::string partitionFile)
 {
     // TODO: Move the reading of the partition file to its own class
@@ -53,11 +125,15 @@ void ModelMetadata::getPartitionMetadata(std::string partitionFile)
     globalExtentX = ncFile.getDim("NX").getSize();
     globalExtentY = ncFile.getDim("NY").getSize();
     netCDF::NcGroup bboxGroup(ncFile.getGroup(bboxName));
-    std::vector<size_t> index(1, mpiMyRank);
-    bboxGroup.getVar("domain_x").getVar(index, &localCornerX);
-    bboxGroup.getVar("domain_y").getVar(index, &localCornerY);
-    bboxGroup.getVar("domain_extent_x").getVar(index, &localExtentX);
-    bboxGroup.getVar("domain_extent_y").getVar(index, &localExtentY);
+
+    std::vector<size_t> rank(1, mpiMyRank);
+    bboxGroup.getVar("domain_x").getVar(rank, &localCornerX);
+    bboxGroup.getVar("domain_y").getVar(rank, &localCornerY);
+    bboxGroup.getVar("domain_extent_x").getVar(rank, &localExtentX);
+    bboxGroup.getVar("domain_extent_y").getVar(rank, &localExtentY);
+
+    readNeighbourData(ncFile);
+
     ncFile.close();
 }
 

core/src/include/ModelArray.hpp

@@ -391,8 +391,7 @@ class ModelArray {
     template <typename T = size_t, typename C, typename I, typename... Args>
     static inline T indexr(C dims, I first, Args... args)
     {
-        std::initializer_list<I> loc { first, args... };
-        return indexer(dims, loc);
+        return indexer(dims, { static_cast<size_t>(first), static_cast<size_t>(args)... });
     }
 
     // Indices as a Dimensions object

core/src/include/ModelArraySlice.hpp

@@ -118,9 +118,44 @@ class ModelArraySlice {
         return buffer;
     }
 
-    ModelArray::DataType& copyToDataSlice(
-        ModelArray::DataType& target, SliceIter& targetIter) const;
-    ModelArraySlice& copyFromDataSlice(const ModelArray::DataType& source, SliceIter& sourceIter);
+private:
+    template <typename S, typename T = S>
+    static void copySliceWithIters(
+        const S& source, SliceIter& sourceIter, T& target, SliceIter targetIter)
+    {
+        const size_t targetNEl = targetIter.nElements(0);
+        const size_t sourceNEl = sourceIter.nElements(0);
+
+        while (!targetIter.isEnd() && !sourceIter.isEnd()) {
+            const size_t targetIndex = targetIter.index();
+            const size_t sourceIndex = sourceIter.index();
+            target(Eigen::seqN(targetIndex, targetNEl, targetIter.step(0)), Eigen::all)
+                = source(Eigen::seqN(sourceIndex, sourceNEl, sourceIter.step(0)), Eigen::all);
+            targetIter.incrementDim(1);
+            sourceIter.incrementDim(1);
+        }
+    }
+
+public:
+    template <typename T>
+    const ModelArraySlice& copyToSlicedBuffer(T& target, SliceIter& targetIter) const
+    {
+        SliceIter iter(slice, data.dimensions());
+
+        copySliceWithIters(data.m_data, iter, target, targetIter);
+        // Return this, even though it is unchanged. The code looks weird if
+        // the return value is the target buffer.
+        return *this;
+    }
+
+    template <typename S>
+    ModelArraySlice& copyFromSlicedBuffer(const S& source, SliceIter& sourceIter)
+    {
+        SliceIter iter(slice, data.dimensions());
+
+        copySliceWithIters(source, sourceIter, data.m_data, iter);
+        return *this;
+    }
 
     iterator begin();
     iterator end();
@@ -132,10 +167,6 @@ class ModelArraySlice {
 private:
     static void copyBetweenMAandMASlice(
         ModelArray& ma, const ModelArraySlice& mas, bool toSlice, const std::string& functionName);
-    static void copySliceWithIters(
-        ModelArray& source, SliceIter& sourceIter, ModelArray& target, SliceIter targetIter);
-    static void copySliceWithItersData(const ModelArray::DataType& source, SliceIter& sourceIter,
-        ModelArray::DataType& target, SliceIter targetIter);
 
     ModelArray& data;
 };

core/src/include/ModelMetadata.hpp

@@ -13,7 +13,9 @@
 #include "include/ModelState.hpp"
 #include "include/Time.hpp"
 
+#include <ncFile.h>
 #include <string>
+#include <vector>
 
 #ifdef USE_MPI
 #include <mpi.h>
@@ -93,13 +95,34 @@ class ModelMetadata {
      */
     void getPartitionMetadata(std::string partitionFile);
 
+    enum Edge { BOTTOM, RIGHT, TOP, LEFT, N_EDGE };
+    // An array to allow the edges to be accessed in the correct order.
+    static constexpr std::array<Edge, N_EDGE> edges = { BOTTOM, RIGHT, TOP, LEFT };
+    std::array<std::string, N_EDGE> edgeNames = { "bottom", "right", "top", "left" };
+
     MPI_Comm mpiComm;
     int mpiSize = 0;
     int mpiMyRank = -1;
-    int localCornerX, localCornerY, localExtentX, localExtentY, globalExtentX, globalExtentY;
+    int localCornerX, localCornerY;
+    int localExtentX, localExtentY;
+    int globalExtentX, globalExtentY;
+    // mpi rank ID and extent for each edge direction
+    std::array<std::vector<int>, N_EDGE> neighbourRanks;
+    std::array<std::vector<int>, N_EDGE> neighbourExtents;
+    std::array<std::vector<int>, N_EDGE> neighbourHaloStarts;
+    std::array<std::vector<int>, N_EDGE> neighbourRanksPeriodic;
+    std::array<std::vector<int>, N_EDGE> neighbourExtentsPeriodic;
+    std::array<std::vector<int>, N_EDGE> neighbourHaloStartsPeriodic;
 #endif
 
 private:
+    /*!
+     * @brief Read neighbour metadata from netCDF file
+     *
+     * @param netCDF file with partition metadata
+     */
+    void readNeighbourData(netCDF::NcFile& ncFile);
+
     TimePoint m_time;
     ConfigMap m_config;
 
@@ -116,6 +139,7 @@ class ModelMetadata {
     bool hasParameters;
 #ifdef USE_MPI
     const std::string bboxName = "bounding_boxes";
+    const std::string neighbourName = "connectivity";
 #endif
 };