halo exchange integration

.gitignore

@@ -24,6 +24,10 @@ __pycache__/
 *.sif
 *.smesh
 .*sw?
+/run/*.cdl
+.cache/
+.clangd
+compile_commands.json
 *.vtk
 
 # Specific files

core/src/CMakeLists.txt

@@ -13,6 +13,7 @@ else()
 endif()
 
 set(BaseSources
+    "ModelMPI.cpp"
     "Logged.cpp"
     "Model.cpp"
     "Iterator.cpp"

core/src/CommonRestartMetadata.cpp

@@ -3,6 +3,7 @@
  */
 
 #include "include/CommonRestartMetadata.hpp"
+#include "include/ModelMetadata.hpp"
 
 #include "include/MissingData.hpp"
 
@@ -15,17 +16,17 @@
 
 namespace Nextsim {
 
-netCDF::NcFile& CommonRestartMetadata::writeStructureType(
-    netCDF::NcFile& ncFile, const ModelMetadata& metadata)
+netCDF::NcFile& CommonRestartMetadata::writeStructureType(netCDF::NcFile& ncFile)
 {
+    auto& metadata = ModelMetadata::getInstance();
     ncFile.putAtt(IStructure::structureNodeName(), metadata.structureName());
     return ncFile;
 }
 
-netCDF::NcFile& CommonRestartMetadata::writeRestartMetadata(
-    netCDF::NcFile& ncFile, const ModelMetadata& metadata)
+netCDF::NcFile& CommonRestartMetadata::writeRestartMetadata(netCDF::NcFile& ncFile)
 {
     // Structure type
+    auto& metadata = ModelMetadata::getInstance();
     ncFile.putAtt(IStructure::structureNodeName(), metadata.structureName());
 
     // As Unix time

core/src/DevStep.cpp

@@ -12,7 +12,6 @@ namespace Nextsim {
 
 DevStep::DevStep()
     : pData(nullptr)
-    , mData(nullptr)
     , m_restartPeriod(0)
 {
 }
@@ -38,14 +37,15 @@ void DevStep::iterate(const TimestepTime& tst)
     pData->update(tst);
     // The state of the model has now advanced by one timestep, so update the
     // model metadata timestamp.
-    mData->incrementTime(tst.step);
-    if ((m_restartPeriod.seconds() > 0) && (mData->time() >= lastOutput + m_restartPeriod)) {
-        std::string currentFileName = mData->time().format(m_restartFileName);
-        pData->writeRestartFile(currentFileName, *mData);
-        lastOutput = mData->time();
+    auto& mData = ModelMetadata::getInstance();
+    mData.incrementTime(tst.step);
+    if ((m_restartPeriod.seconds() > 0) && (mData.time() >= lastOutput + m_restartPeriod)) {
+        std::string currentFileName = mData.time().format(m_restartFileName);
+        pData->writeRestartFile(currentFileName);
+        lastOutput = mData.time();
     }
     // XIOS wants all the fields, every timestep, so I guess that's what everyone gets
-    Module::getImplementation<IDiagnosticOutput>().outputState(pData->getStateDiagnostic(), *mData);
+    Module::getImplementation<IDiagnosticOutput>().outputState(pData->getStateDiagnostic());
 }
 
 void DevStep::setRestartDetails(const Duration& restartPeriod, const std::string& fileName)

core/src/Model.cpp

@@ -42,16 +42,9 @@ static const std::map<int, std::string> keyMap = {
     { Model::RESTARTOUTFILE_KEY, "model.restart_file" },
 };
 
-#ifdef USE_MPI
-Model::Model(MPI_Comm comm)
-#else
 Model::Model()
-#endif
     : iterator(modelStep)
 {
-#ifdef USE_MPI
-    m_etadata.setMpiMetadata(comm);
-#endif
     finalFileName = std::string("restart") + TimePoint::ymdhmsFormat + ".nc";
 }
 
@@ -62,25 +55,6 @@ void Model::configure()
     // Configure logging
     Logged::configure();
 
-    // Start/stop times. Run length will override stop time, if present.
-    std::string startTimeStr
-        = Configured::getConfiguration(keyMap.at(STARTTIME_KEY), std::string());
-    std::string stopTimeStr = Configured::getConfiguration(keyMap.at(STOPTIME_KEY), std::string());
-    std::string runLengthStr
-        = Configured::getConfiguration(keyMap.at(RUNLENGTH_KEY), std::string());
-    std::string stepStr = Configured::getConfiguration(keyMap.at(TIMESTEP_KEY), std::string());
-
-    if (runLengthStr.empty()) {
-        if (stopTimeStr.empty()) {
-            throw std::invalid_argument(std::string("At least one of ") + keyMap.at(STOPTIME_KEY)
-                + " or " + keyMap.at(RUNLENGTH_KEY) + " must be set");
-        } else {
-            m_etadata.setTimes(startTimeStr, TimePoint(stopTimeStr), stepStr);
-        }
-    } else {
-        m_etadata.setTimes(startTimeStr, Duration(runLengthStr), stepStr);
-    }
-    iterator.setStartStopStep(m_etadata.startTime(), m_etadata.stopTime(), m_etadata.stepLength());
     // Configure the missing data value
     MissingData::setValue(
         Configured::getConfiguration(keyMap.at(MISSINGVALUE_KEY), MissingData::defaultValue));
@@ -97,25 +71,42 @@ void Model::configure()
 #ifdef USE_MPI
     std::string partitionFile
         = Configured::getConfiguration(keyMap.at(PARTITIONFILE_KEY), std::string("partition.nc"));
-    m_etadata.getPartitionMetadata(partitionFile);
+    auto& metadata = ModelMetadata::getInstance(partitionFile);
+#else
+    auto& metadata = ModelMetadata::getInstance();
 #endif
 
-#ifdef USE_MPI
-    ModelState initialState(StructureFactory::stateFromFile(initialFileName, m_etadata));
-#else
+    // Start/stop times. Run length will override stop time, if present.
+    std::string startTimeStr
+        = Configured::getConfiguration(keyMap.at(STARTTIME_KEY), std::string());
+    std::string stopTimeStr = Configured::getConfiguration(keyMap.at(STOPTIME_KEY), std::string());
+    std::string runLengthStr
+        = Configured::getConfiguration(keyMap.at(RUNLENGTH_KEY), std::string());
+    std::string stepStr = Configured::getConfiguration(keyMap.at(TIMESTEP_KEY), std::string());
+
+    if (runLengthStr.empty()) {
+        if (stopTimeStr.empty()) {
+            throw std::invalid_argument(std::string("At least one of ") + keyMap.at(STOPTIME_KEY)
+                + " or " + keyMap.at(RUNLENGTH_KEY) + " must be set");
+        } else {
+            metadata.setTimes(startTimeStr, TimePoint(stopTimeStr), stepStr);
+        }
+    } else {
+        metadata.setTimes(startTimeStr, Duration(runLengthStr), stepStr);
+    }
+    iterator.setStartStopStep(metadata.startTime(), metadata.stopTime(), metadata.stepLength());
+
     ModelState initialState(StructureFactory::stateFromFile(initialFileName));
-#endif
 
     // The period with which to write restart files.
     std::string restartPeriodStr
         = Configured::getConfiguration(keyMap.at(RESTARTPERIOD_KEY), std::string("0"));
     restartPeriod = Duration(restartPeriodStr);
 
     // Get the coordinates from the ModelState for persistence
-    m_etadata.extractCoordinates(initialState);
+    metadata.extractCoordinates(initialState);
 
     modelStep.setData(pData);
-    modelStep.setMetadata(m_etadata);
     modelStep.setRestartDetails(restartPeriod, finalFileName);
     pData.setData(initialState.data);
 }
@@ -186,9 +177,9 @@ void Model::run()
 //! Write a restart file for the model.
 void Model::writeRestartFile()
 {
-    std::string formattedFileName = m_etadata.time().format(finalFileName);
-    pData.writeRestartFile(formattedFileName, m_etadata);
+    auto& metadata = ModelMetadata::getInstance();
+    std::string formattedFileName = metadata.time().format(finalFileName);
+    pData.writeRestartFile(formattedFileName);
 }
 
-ModelMetadata& Model::metadata() { return m_etadata; }
 } /* namespace Nextsim */

core/src/ModelMPI.cpp

@@ -0,0 +1,36 @@
+/*!
+ * @author  Tom Meltzer <tdm39@cam.ac.uk>
+ */
+
+#ifdef USE_MPI
+#include "include/ModelMPI.hpp"
+
+namespace Nextsim {
+
+ModelMPI::ModelMPI(MPI_Comm comm)
+    : m_comm(comm)
+{
+    MPI_Comm_size(m_comm, &m_size);
+    MPI_Comm_rank(m_comm, &m_rank);
+    static bool doneOnce = doOnce();
+    isInitialized = true;
+}
+
+MPI_Comm ModelMPI::getComm() const { return m_comm; }
+
+int ModelMPI::getSize() const { return m_size; }
+
+int ModelMPI::getRank() const { return m_rank; }
+
+bool ModelMPI::doOnce()
+{
+    // Register the finalization function here
+    Finalizer::registerUnique(finalize);
+    return true;
+}
+
+void ModelMPI::finalize() { }
+
+} // namespace Nextsim
+
+#endif // USE_MPI

core/src/ModelMetadata.cpp

@@ -1,17 +1,25 @@
 /*!
  * @author  Tim Spain <timothy.spain@nersc.no>
+ * @author  Tom Meltzer <tdm39@cam.ac.uk>
  */
 
 #include "include/ModelMetadata.hpp"
 
+#include "include/Finalizer.hpp"
 #include "include/IStructure.hpp"
+#include "include/ModelMPI.hpp"
 #include "include/NextsimModule.hpp"
 #ifdef USE_XIOS
 #include "include/Xios.hpp"
 #endif
 #include "include/gridNames.hpp"
+#include <array>
+#include <cstddef>
+#include <functional>
+#include <vector>
 
 #ifdef USE_MPI
+#include "mpi.h"
 #include <ncDim.h>
 #include <ncFile.h>
 #include <ncGroup.h>
@@ -26,25 +34,86 @@ const std::string& ModelMetadata::structureName() const
 }
 
 #ifdef USE_MPI
-ModelMetadata::ModelMetadata(std::string partitionFile, MPI_Comm comm)
+ModelMetadata::ModelMetadata(std::string partitionFile)
 {
-    setMpiMetadata(comm);
+    if (partitionFile.empty()) {
+        throw std::runtime_error(
+            "ModelMetadata :: getInstance() called without partition file in MPI build.");
+    }
     getPartitionMetadata(partitionFile);
+    static bool doneOnce = doOnce();
+    isInitialized = true;
 }
 
-void ModelMetadata::setMpiMetadata(MPI_Comm comm)
+void ModelMetadata::readNeighbourData(netCDF::NcFile& ncFile)
 {
-    mpiComm = comm;
-    MPI_Comm_size(mpiComm, &mpiSize);
-    MPI_Comm_rank(mpiComm, &mpiMyRank);
+    netCDF::NcGroup neighbourGroup(ncFile.getGroup(neighbourName));
+    std::string varName {};
+    auto& modelMPI = ModelMPI::getInstance();
+    auto mpiSize = modelMPI.getSize();
+    auto mpiMyRank = modelMPI.getRank();
+    enum BoundaryType { nonPeriodic, periodic };
+    for (BoundaryType btype : { nonPeriodic, periodic }) {
+        // Use btype as needed
+        std::array<std::string, 4> suffixes = { "_neighbour_ids", "_neighbour_halos",
+            "_neighbour_halo_send", "_neighbour_halo_recv" };
+        if (btype == periodic) {
+            for (auto& suffix : suffixes) {
+                suffix += "_periodic";
+            }
+        }
+        for (auto edge : edges) {
+            size_t nStart = 0; // start point in metadata arrays
+            size_t count = 0; // 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);
+            std::vector<std::reference_wrapper<std::vector<int>>> arrays;
+
+            if (btype == nonPeriodic) {
+                arrays = { neighbourRanks[edge], neighbourExtents[edge], neighbourHaloSend[edge],
+                    neighbourHaloRecv[edge] };
+            } else if (btype == periodic) {
+                arrays = { neighbourRanksPeriodic[edge], neighbourExtentsPeriodic[edge],
+                    neighbourHaloSendPeriodic[edge], neighbourHaloRecvPeriodic[edge] };
+            }
+
+            varName = edgeNames[edge] + "_neighbours";
+            if (btype == periodic) {
+                varName += "_periodic";
+            }
+            neighbourGroup.getVar(varName).getVar(
+                { 0 }, { static_cast<size_t>(mpiSize) }, numNeighbours.data());
+
+            // compute start index for each process
+            MPI_Exscan(&numNeighbours[mpiMyRank], &nStart, 1, MPI_INT, MPI_SUM, modelMPI.getComm());
+            if (mpiMyRank == 0) {
+                // MPI_Exscan is undefined on the first rank. So to be safe we manually set nStart
+                // to 0. (see e.g., https://www.open-mpi.org/doc/v4.1/man3/MPI_Exscan.3.php)
+                nStart = 0;
+            }
+            // how many elements to read for each process
+            count = numNeighbours[mpiMyRank];
+
+            if (count) {
+                // initialize neighbour info to zero
+                for (size_t i = 0; i < arrays.size(); ++i) {
+                    arrays[i].get().resize(count, 0);
+                    varName = edgeNames[edge] + suffixes[i];
+                    neighbourGroup.getVar(varName).getVar(
+                        { nStart }, { count }, arrays[i].get().data());
+                }
+            }
+        }
+    }
 }
 
 void ModelMetadata::getPartitionMetadata(std::string partitionFile)
 {
-    // TODO: Move the reading of the partition file to its own class
     netCDF::NcFile ncFile(partitionFile, netCDF::NcFile::read);
     int sizes = ncFile.getDim("L").getSize();
     int nBoxes = ncFile.getDim("P").getSize();
+    auto& modelMPI = ModelMPI::getInstance();
+    auto mpiSize = modelMPI.getSize();
     if (nBoxes != mpiSize) {
         std::string errorMsg = "Number of MPI ranks " + std::to_string(mpiSize) + " <> "
             + std::to_string(nBoxes) + "\n";
@@ -53,14 +122,32 @@ 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, modelMPI.getRank());
+    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();
 }
 
+int ModelMetadata::getLocalCornerX() const { return localCornerX; }
+int ModelMetadata::getLocalCornerY() const { return localCornerY; }
+int ModelMetadata::getLocalExtentX() const { return localExtentX; }
+int ModelMetadata::getLocalExtentY() const { return localExtentY; }
+int ModelMetadata::getGlobalExtentX() const { return globalExtentX; }
+int ModelMetadata::getGlobalExtentY() const { return globalExtentY; }
+#else
+
+ModelMetadata::ModelMetadata()
+{
+    isInitialized = true;
+    static bool doneOnce = doOnce();
+}
+
 #endif
 
 const ModelState& ModelMetadata::extractCoordinates(const ModelState& state)
@@ -146,4 +233,13 @@ void ModelMetadata::incrementTime(const Duration& step)
 #endif
 }
 
+void ModelMetadata::finalize() { }
+
+bool ModelMetadata::doOnce()
+{
+    // Register the finalization function here
+    Finalizer::registerUnique(finalize);
+    return true;
+}
+
 } /* namespace Nextsim */