diff --git a/core/src/ModelMetadata.cpp b/core/src/ModelMetadata.cpp index 79c981c0b..538b9e529 100644 --- a/core/src/ModelMetadata.cpp +++ b/core/src/ModelMetadata.cpp @@ -131,6 +131,12 @@ void ModelMetadata::getPartitionMetadata(std::string partitionFile) readNeighbourData(ncFile); + // gather rank extents in X & Y direction for all processes + rankExtentsX.resize(modelMPI.getSize(), 0); + rankExtentsY.resize(modelMPI.getSize(), 0); + MPI_Allgather(&localExtentX, 1, MPI_INT, rankExtentsX.data(), 1, MPI_INT, modelMPI.getComm()); + MPI_Allgather(&localExtentY, 1, MPI_INT, rankExtentsY.data(), 1, MPI_INT, modelMPI.getComm()); + ncFile.close(); } @@ -140,6 +146,8 @@ int ModelMetadata::getLocalExtentX() const { return localExtentX; } int ModelMetadata::getLocalExtentY() const { return localExtentY; } int ModelMetadata::getGlobalExtentX() const { return globalExtentX; } int ModelMetadata::getGlobalExtentY() const { return globalExtentY; } +std::vector ModelMetadata::getRankExtentsX() const { return rankExtentsX; } +std::vector ModelMetadata::getRankExtentsY() const { return rankExtentsY; } #else ModelMetadata::ModelMetadata() diff --git a/core/src/include/Halo.hpp b/core/src/include/Halo.hpp index 327ec000b..5cae1c562 100644 --- a/core/src/include/Halo.hpp +++ b/core/src/include/Halo.hpp @@ -7,7 +7,7 @@ * * All functionality for halo exchange between MPI ranks is contained in this class. * - * Halo supports the main data structures of NextSim e.g., ModelArray and DGVector. + * Halo supports the main data structures of NextSim e.g., ModelArray, DGVector and CGVector. * * The halos are exchange via one-sided MPI communication using RMA. */ @@ -18,14 +18,17 @@ #include #include #include +#include #include #include "Slice.hpp" +#include "cgVector.hpp" #include "dgVector.hpp" #include "include/ModelArray.hpp" #include "include/ModelArraySlice.hpp" #include "include/ModelMPI.hpp" #include "include/ModelMetadata.hpp" +#include "include/dgVectorHolder.hpp" #include "mpi.h" #ifndef HALOWIDTH @@ -50,13 +53,35 @@ class Halo { intializeHaloMetadata(); } + /*! + * @brief Constructs a halo object from DGVector + */ + template Halo(DGVectorHolder& dgvh) + { + m_numComps = N; + setSpatialDims(); + intializeHaloMetadata(); + } + /*! * @brief Constructs a halo object from DGVector */ template Halo(DGVector& dgv) { m_numComps = N; - isVertex = false; + setSpatialDims(); + intializeHaloMetadata(); + } + + /*! + * @brief Constructs a halo object from CGVector + */ + template Halo(CGVector& cgv) + { + m_numComps = 1; + isCG = true; + CGdegree = N; + nCells = CGdegree; setSpatialDims(); intializeHaloMetadata(); } @@ -81,10 +106,21 @@ class Halo { m_innerNx += 1; m_innerNy += 1; } + // extend dimensions for CGVectors + if (isCG) { + m_innerNy = m_innerNy * CGdegree + 1; + m_innerNx = m_innerNx * CGdegree + 1; + } // inner dimension of domain excluding the halo cells m_Nx = m_innerNx + 2 * haloWidth; m_Ny = m_innerNy + 2 * haloWidth; + + // each additional halo cell add CGdegree more points + if (isCG) { + m_Nx = m_innerNx + 2 * haloWidth * CGdegree; + m_Ny = m_innerNy + 2 * haloWidth * CGdegree; + } } /** @@ -101,7 +137,12 @@ class Halo { void intializeHaloMetadata() { // number of halo cells (should be general for any halo width) - m_numHaloCells = 2 * haloWidth * (m_innerNx + m_innerNy + 2 * haloWidth); + if (not isCG) { + m_numHaloCells = 2 * haloWidth * (m_innerNx + m_innerNy + 2 * haloWidth); + } else { + m_numHaloCells + = 2 * haloWidth * CGdegree * (m_innerNx + m_innerNy + 2 * haloWidth * CGdegree); + } // need send / recv buffers for each component (e.g., each DGCOMP) send.resize(m_numComps); @@ -114,28 +155,6 @@ class Halo { // order is Bottom, Right, Top, Left m_edgeLengths = { m_innerNx, m_innerNy, m_innerNx, m_innerNy }; - - m_outerSlices = { - { Edge::BOTTOM, VBounds({ { 1, m_innerNx + haloWidth }, { 0 } }) }, - { Edge::RIGHT, VBounds({ { -1 }, { 1, m_innerNy + haloWidth } }) }, - { Edge::TOP, VBounds({ { 1, m_innerNx + haloWidth }, { -1 } }) }, - { Edge::LEFT, VBounds({ { 0 }, { 1, m_innerNy + haloWidth } }) }, - }; - if (isVertex) { - m_innerSlices = { - { Edge::BOTTOM, VBounds({ { 1, m_innerNx + haloWidth }, { 2 } }) }, - { Edge::RIGHT, VBounds({ { -3 }, { 1, m_innerNy + haloWidth } }) }, - { Edge::TOP, VBounds({ { 1, m_innerNx + haloWidth }, { -3 } }) }, - { Edge::LEFT, VBounds({ { 2 }, { 1, m_innerNy + haloWidth } }) }, - }; - } else { - m_innerSlices = { - { Edge::BOTTOM, VBounds({ { 1, m_innerNx + haloWidth }, { 1 } }) }, - { Edge::RIGHT, VBounds({ { -2 }, { 1, m_innerNy + haloWidth } }) }, - { Edge::TOP, VBounds({ { 1, m_innerNx + haloWidth }, { -2 } }) }, - { Edge::LEFT, VBounds({ { 1 }, { 1, m_innerNy + haloWidth } }) }, - }; - } } static const size_t haloWidth = HALOWIDTH; // how many cells wide is the halo region @@ -155,11 +174,11 @@ class Halo { size_t m_innerNy; // local extent in y-direction size_t m_numHaloCells; // number of halo cells size_t m_numComps; // number of DG components + int nCells = 1; // this is the number of points to grab per cell, per direction (1 for + // everything except CG Fields) std::array m_edgeLengths; // array containing length of each edge std::array edges = ModelMetadata::edges; // array of edge enums - std::map m_outerSlices; - std::map m_innerSlices; std::map oppositeEdge = { { Edge::LEFT, Edge::RIGHT }, { Edge::RIGHT, Edge::LEFT }, @@ -191,7 +210,9 @@ class Halo { MPI_Win m_win; // RMA memory window object (used for sharing send buffers between ranks) - bool isVertex; // some ModelArrays can be of type VERTEX + bool isVertex = false; // some ModelArrays can be of type VERTEX + bool isCG = false; // is layout CGVector + size_t CGdegree = 0; std::vector> send; // buffer to store halo region that will be read by other ranks @@ -234,35 +255,64 @@ class Halo { */ template void populateSendBuffers(S& source) { - tempBuffer.resize(m_numHaloCells, m_numComps); - tempBuffer = 0.; - - for (auto edge : edges) { - size_t beg = std::accumulate(m_edgeLengths.begin(), m_edgeLengths.begin() + edge, 0); - size_t num = m_edgeLengths.at(edge); - - Slice sourceSlice = m_innerSlices.at(edge); - SliceIter sourceIter(sourceSlice, { m_Nx, m_Ny }); - - if (isVertical(edge)) { - Slice tempBufferSlice = { { {}, { beg, beg + num } } }; - // spatial dims are flattened into 1-D. - SliceIter tempBufferIter(tempBufferSlice, { 1, m_numHaloCells }); - ModelArraySlice::copySliceWithIters(source, sourceIter, tempBuffer, tempBufferIter); - } else { - Slice tempBufferSlice = { { { beg, beg + num }, {} } }; - // spatial dims are flattened into 1-D. - SliceIter tempBufferIter(tempBufferSlice, { m_numHaloCells, 1 }); - ModelArraySlice::copySliceWithIters(source, sourceIter, tempBuffer, tempBufferIter); + size_t buffer_len = m_numHaloCells; + if (isCG) { + buffer_len = m_numHaloCells / CGdegree; + } + for (size_t comp = 0; comp < m_numComps; ++comp) { + Eigen::Map> source_map( + source.col(comp).data(), m_Nx, m_Ny, Eigen::InnerStride<>(m_numComps)); + Eigen::Map buffer_map(send[comp].data(), buffer_len, nCells); + for (auto edge : edges) { + sourceToSendBuffer(edge, buffer_map, source_map); } } - // we need to copy into std vector send buffer because MPI doesn't work with Eigen - // Arrays - for (size_t i = 0; i < m_numComps; i++) { - typedef Eigen::Map MapType; - MapType map(send[i].data(), m_numHaloCells, 1); - // map is connected with the send buffer so the following line sets the data in send - map = tempBuffer.col(i); + } + + /** + * @brief Calculate recv buffer positions and offsets for halo transfer + * + * @param[in,out] count The count of elements to be communicated, increased by haloWidth + * @param[in,out] disp The displacement value, adjusted based on opposite edge + * @param[in,out] recvOffset The receive offset, adjusted based on edge + * @param[in] edge The edge along which communication occurs + */ + void recvPositions(int& fromRank, size_t& count, size_t& disp, size_t& recvOffset, Edge edge, + const size_t neighbourIndex, const size_t cell, bool isPeriodic) + { + auto& metadata = ModelMetadata::getInstance(); + if (isPeriodic) { + fromRank = metadata.neighbourRanksPeriodic[edge][neighbourIndex]; + count = metadata.neighbourExtentsPeriodic[edge][neighbourIndex]; + disp = metadata.neighbourHaloSendPeriodic[edge][neighbourIndex]; + recvOffset = metadata.neighbourHaloRecvPeriodic[edge][neighbourIndex]; + } else { + fromRank = metadata.neighbourRanks[edge][neighbourIndex]; + count = metadata.neighbourExtents[edge][neighbourIndex]; + disp = metadata.neighbourHaloSend[edge][neighbourIndex]; + recvOffset = metadata.neighbourHaloRecv[edge][neighbourIndex]; + } + if (isVertex) { + count = count + 1; + disp = disp + oppositeEdge.at(edge); + recvOffset = recvOffset + edge; + } + if (isCG) { + count = CGdegree * count + 1; + disp = (disp > 0) ? CGdegree * disp + oppositeEdge.at(edge) : 0; + recvOffset = (recvOffset > 0) ? CGdegree * recvOffset + edge : 0; + + // recvOffset is the offset in the recv buffer and this belongs to the current rank + recvOffset = recvOffset + m_numHaloCells / nCells * cell; + + // disp is the offset in the "sending" buffer which belongs to rank "fromRank" + // Therefore we need to compute how many halo cells that rank has to work out the offset + // for each cell + auto extentX = CGdegree * metadata.getRankExtentsX()[fromRank] + 1; + auto extentY = CGdegree * metadata.getRankExtentsY()[fromRank] + 1; + auto fromRank_numHaloCells + = 2 * haloWidth * CGdegree * (extentX + extentY + 2 * haloWidth * CGdegree); + disp = disp + fromRank_numHaloCells / nCells * cell; } } @@ -273,48 +323,39 @@ class Halo { void populateRecvBuffers() { + int fromRank; + size_t count, disp, recvOffset; // do halo exchange for each component for (size_t comp = 0; comp < m_numComps; comp++) { // open memory window to send buffer on other ranks openMemoryWindow(comp); auto& metadata = ModelMetadata::getInstance(); - // get non-periodic neighbours and populate recv buffer (if the exist) for (auto edge : edges) { + + // get neighbours (if they exist) auto numNeighbours = metadata.neighbourRanks[edge].size(); if (numNeighbours) { // get data for each neighbour that exists along a given edge for (size_t i = 0; i < numNeighbours; ++i) { - int fromRank = metadata.neighbourRanks[edge][i]; - size_t count = metadata.neighbourExtents[edge][i]; - size_t disp = metadata.neighbourHaloSend[edge][i]; - size_t recvOffset = metadata.neighbourHaloRecv[edge][i]; - if (isVertex) { - vertexAdjustedPositions( - count = count, disp = disp, recvOffset = recvOffset, edge = edge); + for (size_t cell = 0; cell < nCells; ++cell) { + recvPositions(fromRank, count, disp, recvOffset, edge, i, cell, false); + MPI_Get(&recv[comp][recvOffset], count, MPI_DOUBLE, fromRank, disp, + count, MPI_DOUBLE, m_win); } - MPI_Get(&recv[comp][recvOffset], count, MPI_DOUBLE, fromRank, disp, count, - MPI_DOUBLE, m_win); } } - } - // get periodic neighbours and populate recv buffer (if they exist) - for (auto edge : edges) { - auto numNeighbours = metadata.neighbourRanksPeriodic[edge].size(); + // get periodic neighbours (if they exist) + numNeighbours = metadata.neighbourRanksPeriodic[edge].size(); if (numNeighbours) { // get data for each neighbour that exists along a given edge for (size_t i = 0; i < numNeighbours; ++i) { - int fromRank = metadata.neighbourRanksPeriodic[edge][i]; - size_t count = metadata.neighbourExtentsPeriodic[edge][i]; - size_t disp = metadata.neighbourHaloSendPeriodic[edge][i]; - size_t recvOffset = metadata.neighbourHaloRecvPeriodic[edge][i]; - if (isVertex) { - vertexAdjustedPositions( - count = count, disp = disp, recvOffset = recvOffset, edge = edge); + for (size_t cell = 0; cell < nCells; ++cell) { + recvPositions(fromRank, count, disp, recvOffset, edge, i, cell, true); + MPI_Get(&recv[comp][recvOffset], count, MPI_DOUBLE, fromRank, disp, + count, MPI_DOUBLE, m_win); } - MPI_Get(&recv[comp][recvOffset], count, MPI_DOUBLE, fromRank, disp, count, - MPI_DOUBLE, m_win); } } } @@ -323,33 +364,6 @@ class Halo { // moving on) closeMemoryWindow(); } - - // copy from the recv std vector into an eigen array - for (size_t i = 0; i < m_numComps; i++) { - typedef Eigen::Map MapType; - MapType map(recv[i].data(), m_numHaloCells, 1); - // map is connected with the recv buffer so the following line copies the data into - // tempBuffer - tempBuffer.col(i) = map; - } - } - - /** - * @brief Adjusts positions and offsets for vertex communications in the halo region - * - * Updates the count, displacement and receive offset values for vertex communications - * across halo boundaries according to the halo width and specified edge. - * - * @param[in,out] count The count of elements to be communicated, increased by haloWidth - * @param[in,out] disp The displacement value, adjusted based on opposite edge - * @param[in,out] recvOffset The receive offset, adjusted based on edge - * @param[in] edge The edge along which communication occurs - */ - void vertexAdjustedPositions(size_t& count, size_t& disp, size_t& recvOffset, const Edge& edge) - { - count = count + haloWidth; - disp = disp + oppositeEdge.at(edge) * haloWidth; - recvOffset = recvOffset + edge * haloWidth; } public: @@ -373,6 +387,106 @@ class Halo { */ size_t getInnerSize() { return m_innerNx * m_innerNy; } + void sendBufferPositions(int& idx_a, int& idx_b, int& offset, const Edge edge) + { + int vertexOffset = 0; + if (isVertex || isCG) { + vertexOffset = 1; + } + offset = std::accumulate(m_edgeLengths.begin(), m_edgeLengths.begin() + edge, 0); + switch (edge) { + case Edge::LEFT: + idx_a = nCells + vertexOffset; + idx_b = 2 * nCells - 1 + vertexOffset; + break; + case Edge::RIGHT: + idx_a = m_Nx - 2 * nCells - vertexOffset; + idx_b = m_Nx - nCells - 1 - vertexOffset; + break; + case Edge::BOTTOM: + idx_a = nCells + vertexOffset; + idx_b = 2 * nCells - 1 + vertexOffset; + break; + case Edge::TOP: + idx_a = m_Ny - 2 * nCells - vertexOffset; + idx_b = m_Ny - nCells - 1 - vertexOffset; + break; + default: + throw std::runtime_error("Unrecognised edge type"); + } + } + + template + void sourceToSendBuffer(const Edge edge, Eigen::Map& buffer_map, + Eigen::Map>& source_map) + { + int idx_a, idx_b, offset; + sendBufferPositions(idx_a, idx_b, offset, edge); + if (isVertical(edge)) { + buffer_map.transpose()(Eigen::all, Eigen::seq(offset, offset + m_innerNy - 1)) + = source_map(Eigen::seq(idx_a, idx_b), Eigen::seq(nCells, Eigen::last - nCells)); + } else { + buffer_map(Eigen::seq(offset, offset + m_innerNx - 1), Eigen::all) + = source_map(Eigen::seq(nCells, Eigen::last - nCells), Eigen::seq(idx_a, idx_b)); + } + } + + void recvBufferPositions(int& idx_a, int& idx_b, int& offset, const Edge edge) + { + offset = std::accumulate(m_edgeLengths.begin(), m_edgeLengths.begin() + edge, 0); + switch (edge) { + case Edge::LEFT: + idx_a = 0; + idx_b = nCells - 1; + break; + case Edge::RIGHT: + idx_a = m_Nx - nCells; + idx_b = m_Nx - 1; + break; + case Edge::BOTTOM: + idx_a = 0; + idx_b = nCells - 1; + break; + case Edge::TOP: + idx_a = m_Ny - nCells; + idx_b = m_Ny - 1; + break; + default: + throw std::runtime_error("Unrecognised edge type"); + } + } + + template + void recvBufferToTarget(const Edge edge, Eigen::Map& buffer_map, + Eigen::Map>& target_map) + { + int idx_a, idx_b, offset; + recvBufferPositions(idx_a, idx_b, offset, edge); + if (isVertical(edge)) { + target_map(Eigen::seq(idx_a, idx_b), Eigen::seq(nCells, Eigen::last - nCells)) + = buffer_map.transpose()(Eigen::all, Eigen::seq(offset, offset + m_innerNy - 1)); + } else { + target_map(Eigen::seq(nCells, Eigen::last - nCells), Eigen::seq(idx_a, idx_b)) + = buffer_map(Eigen::seq(offset, offset + m_innerNx - 1), Eigen::all); + } + } + + template void populateTarget(T& target) + { + size_t buffer_len = m_numHaloCells; + if (isCG) { + buffer_len = m_numHaloCells / CGdegree; + } + for (size_t comp = 0; comp < m_numComps; ++comp) { + Eigen::Map> target_map( + target.col(comp).data(), m_Nx, m_Ny, Eigen::InnerStride<>(m_numComps)); + Eigen::Map buffer_map(recv[comp].data(), buffer_len, nCells); + for (auto edge : edges) { + recvBufferToTarget(edge, buffer_map, target_map); + } + } + } + /*! * @brief Get inner block of ModelArray from tempData * @@ -402,47 +516,11 @@ class Halo { ModelArraySlice::copySliceWithIters(source, sourceIter, target, targetIter); } - /*! - * @brief Update a ModelArray with data from - * - * @params dgvec DGVector which we intend to update across MPI ranks based on halo cells - * note that the start index is offset by 1 and the loop limit is size() - 2 because - * the edge of each domain is 2 less than the length of the expanded halo region - * (see diagram below - the empty cells are skipped by going from i+1 to size()-2) - * ┌─┬─┬─┬─┐ - * │ │x│x│ │ - * ├─┼─┼─┼─┤ - * │x│o│o│x│ - * ├─┼─┼─┼─┤ - * │x│o│o│x│ o = original data - * ├─┼─┼─┼─┤ x = mpi halo data (from recv) - * │ │x│x│ │ (empty) = unused data in DGVector - * └─┴─┴─┴─┘ - */ template void exchangeHalos(T& target) { populateSendBuffers(target); populateRecvBuffers(); - - for (auto edge : edges) { - size_t beg = std::accumulate(m_edgeLengths.begin(), m_edgeLengths.begin() + edge, 0); - size_t num = m_edgeLengths.at(edge); - - Slice targetSlice = m_outerSlices.at(edge); - SliceIter targetIter(targetSlice, { m_Nx, m_Ny }); - - if (isVertical(edge)) { - Slice tempBufferSlice = { { {}, { beg, beg + num } } }; - // spatial dims are flattened into 1-D. - SliceIter tempBufferIter(tempBufferSlice, { 1, m_numHaloCells }); - ModelArraySlice::copySliceWithIters(tempBuffer, tempBufferIter, target, targetIter); - } else { - Slice tempBufferSlice = { { { beg, beg + num }, {} } }; - // spatial dims are flattened into 1-D. - SliceIter tempBufferIter(tempBufferSlice, { m_numHaloCells, 1 }); - ModelArraySlice::copySliceWithIters(tempBuffer, tempBufferIter, target, targetIter); - } - } + populateTarget(target); } }; } // end of nextsim namespace diff --git a/core/src/include/ModelMetadata.hpp b/core/src/include/ModelMetadata.hpp index 4b4855e35..6813fe30c 100644 --- a/core/src/include/ModelMetadata.hpp +++ b/core/src/include/ModelMetadata.hpp @@ -181,6 +181,16 @@ class ModelMetadata { * @return The total number of grid points in Y for the global domain. */ int getGlobalExtentY() const; + /*! + * @brief Gets the extents of the grid in the X direction for all ranks. + * @return A vector containing the number of grid points in X for each rank. + */ + std::vector getRankExtentsX() const; + /*! + * @brief Gets the extents of the grid in the Y direction for all ranks. + * @return A vector containing the number of grid points in Y for each rank. + */ + std::vector getRankExtentsY() const; enum Edge { BOTTOM, RIGHT, TOP, LEFT, N_EDGE }; // An array to allow the edges to be accessed in the correct order. @@ -235,6 +245,8 @@ class ModelMetadata { int localCornerX, localCornerY; int localExtentX, localExtentY; int globalExtentX, globalExtentY; + std::vector rankExtentsX; // vector containing domain extents for each rank x-direction + std::vector rankExtentsY; // vector containing domain extents for each rank y-direction const std::string bboxName = "bounding_boxes"; const std::string neighbourName = "connectivity"; #endif diff --git a/core/src/modules/DynamicsModule/BBMDynamics.cpp b/core/src/modules/DynamicsModule/BBMDynamics.cpp index bffeb7f4a..396256280 100644 --- a/core/src/modules/DynamicsModule/BBMDynamics.cpp +++ b/core/src/modules/DynamicsModule/BBMDynamics.cpp @@ -4,6 +4,9 @@ */ #include "include/BBMDynamics.hpp" +#ifdef USE_MPI +#include "include/ModelMPI.hpp" +#endif #include "include/constants.hpp" #include "include/gridNames.hpp" @@ -134,7 +137,14 @@ void BBMDynamics::setData(const ModelState::DataMap& ms) void BBMDynamics::update(const TimestepTime& tst) { +#ifdef USE_MPI + auto& modelMPI = ModelMPI::getInstance(); + if (modelMPI.getRank() == 0) { + std::cout << tst.start << std::endl; + } +#else std::cout << tst.start << std::endl; +#endif // set the forcing velocities kernel.setData(uWindName, uwind); diff --git a/core/src/modules/DynamicsModule/MEVPDynamics.cpp b/core/src/modules/DynamicsModule/MEVPDynamics.cpp index 0d1455b00..e2135a818 100644 --- a/core/src/modules/DynamicsModule/MEVPDynamics.cpp +++ b/core/src/modules/DynamicsModule/MEVPDynamics.cpp @@ -5,6 +5,9 @@ */ #include "include/MEVPDynamics.hpp" +#ifdef USE_MPI +#include "include/ModelMPI.hpp" +#endif #include "include/constants.hpp" #include "include/gridNames.hpp" @@ -102,7 +105,14 @@ void MEVPDynamics::setData(const ModelState::DataMap& ms) void MEVPDynamics::update(const TimestepTime& tst) { +#ifdef USE_MPI + auto& modelMPI = ModelMPI::getInstance(); + if (modelMPI.getRank() == 0) { + std::cout << tst.start << std::endl; + } +#else std::cout << tst.start << std::endl; +#endif // set the forcing velocities kernel.setData(uWindName, uwind); diff --git a/core/src/modules/DynamicsModule/include/FreeDriftDynamics.hpp b/core/src/modules/DynamicsModule/include/FreeDriftDynamics.hpp index c1de4848a..a13601fc5 100644 --- a/core/src/modules/DynamicsModule/include/FreeDriftDynamics.hpp +++ b/core/src/modules/DynamicsModule/include/FreeDriftDynamics.hpp @@ -7,6 +7,9 @@ #define FREEDRIFTDYNAMICS_HPP #include "include/FreeDriftDynamicsKernel.hpp" +#ifdef USE_MPI +#include "include/ModelMPI.hpp" +#endif #include "include/IDynamics.hpp" #include "include/dgVectorHolder.hpp" @@ -35,7 +38,14 @@ class FreeDriftDynamics : public IDynamics { std::string getName() const override { return "FreeDriftDynamics"; } void update(const TimestepTime& tst) override { +#ifdef USE_MPI + auto& modelMPI = ModelMPI::getInstance(); + if (modelMPI.getRank() == 0) { + std::cout << tst.start << std::endl; + } +#else std::cout << tst.start << std::endl; +#endif // set the forcing velocities kernel.setData(uOceanName, uocean); diff --git a/core/test/HaloExchangeCB_test.cpp b/core/test/HaloExchangeCB_test.cpp index 3d5ffc7c5..7c20933f3 100644 --- a/core/test/HaloExchangeCB_test.cpp +++ b/core/test/HaloExchangeCB_test.cpp @@ -14,8 +14,10 @@ #include "ModelMPI.hpp" #include "ModelMetadata.hpp" +#include "cgVector.hpp" #include "include/DGModelArray.hpp" #include "include/Halo.hpp" +#include "include/Interpolations.hpp" namespace Nextsim { @@ -23,6 +25,7 @@ const std::string testFilesDir = TEST_FILES_DIR; const std::string file = testFilesDir + "/partition_metadata_3_cb.nc"; static const int DG = 3; +static const int CGDEGREE = 2; static const bool debug = false; // reference data for each process @@ -258,13 +261,13 @@ MPI_TEST_CASE("DGVector", 3) ParametricMesh smesh(CARTESIAN); smesh.nx = localNx; smesh.ny = localNy; - smesh.nnodes = localNx * localNy; smesh.nelements = localNx * localNy; - smesh.vertices.resize(smesh.nelements, Eigen::NoChange); - for (size_t i = 0; i < localNx; ++i) { - for (size_t j = 0; j < localNy; ++j) { - smesh.vertices(i * localNy + j, 0) = i; - smesh.vertices(i * localNy + j, 1) = j; + smesh.nnodes = (localNx + 1) * (localNy + 1); + smesh.vertices.resize(smesh.nnodes, Eigen::NoChange); + for (size_t i = 0; i < localNx + 1; ++i) { + for (size_t j = 0; j < localNy + 1; ++j) { + smesh.vertices(i * (localNy + 1) + j, 0) = i; + smesh.vertices(i * (localNy + 1) + j, 1) = j; } } @@ -304,5 +307,124 @@ MPI_TEST_CASE("DGVector", 3) } } } + + // initialize CGVector (after halo exchange on DGVector) + CGVector cgVector; + cgVector.resize_by_mesh(smesh); + Interpolations::DG2CG(smesh, cgVector, testData); + + // create halo for testData model array + Halo haloCG(testData); + // haloCG.exchangeHalos(cgVector); +} + +MPI_TEST_CASE("CGVector", 3) +{ + auto& modelMPI = ModelMPI::getInstance(); + auto& metadata = ModelMetadata::getInstance(); + + const size_t nx = metadata.getGlobalExtentX(); + const size_t ny = metadata.getGlobalExtentY(); + const size_t localNx = metadata.getLocalExtentX() + 2 * Halo::haloWidth; + const size_t localNy = metadata.getLocalExtentY() + 2 * Halo::haloWidth; + const size_t offsetX = metadata.getLocalCornerX(); + const size_t offsetY = metadata.getLocalCornerY(); + + ModelArray::setDimension(ModelArray::Dimension::X, nx, localNx, offsetX); + ModelArray::setDimension(ModelArray::Dimension::Y, ny, localNy, offsetY); + ModelArray::setNComponents(ModelArray::Type::DG, DG); + + ParametricMesh smesh(CARTESIAN); + smesh.nx = localNx; + smesh.ny = localNy; + smesh.nelements = localNx * localNy; + smesh.nnodes = (localNx + 1) * (localNy + 1); + smesh.vertices.resize(smesh.nnodes, Eigen::NoChange); + for (size_t i = 0; i < localNx + 1; ++i) { + for (size_t j = 0; j < localNy + 1; ++j) { + smesh.vertices(i * (localNy + 1) + j, 0) = i; + smesh.vertices(i * (localNy + 1) + j, 1) = j; + } + } + + auto CGNx = localNx * CGDEGREE + 1; + auto CGNy = localNy * CGDEGREE + 1; + auto CGoffsetX = offsetX * CGDEGREE + 1; + auto CGoffsetY = offsetY * CGDEGREE + 1; + + // create CGVector + CGVector cgVector; + cgVector.resize_by_mesh(smesh); + cgVector.zero(); + + // Interpolations::DG2CG(smesh, cgVector, testData); + + // initialize data + for (size_t j = CGDEGREE; j < CGNy - CGDEGREE; ++j) { + for (size_t i = CGDEGREE; i < CGNx - CGDEGREE; ++i) { + cgVector(i + j * CGNx) + = ((i - CGDEGREE) + CGoffsetX) * 100 + ((j - CGDEGREE) + CGoffsetY); + } + } + + // create halo for testData model array + Halo halo(cgVector); + halo.exchangeHalos(cgVector); + + std::array, 3> cgRefDataAllProcs = { + std::vector({ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 101, 201, 301, 401, 501, 601, 701, + 801, 901, 1001, 1101, 1201, 1301, 1401, 1501, 1601, 1701, 0, 0, 102, 202, 302, 402, 502, + 602, 702, 802, 902, 1002, 1102, 1202, 1302, 1402, 1502, 1602, 1702, 0, 0, 103, 203, 303, + 403, 503, 603, 703, 803, 903, 1003, 1103, 1203, 1303, 1403, 1503, 1603, 1703, 0, 0, 104, + 204, 304, 404, 504, 604, 704, 804, 904, 1004, 1104, 1204, 1304, 1404, 1504, 1604, 1704, + 0, 0, 105, 205, 305, 405, 505, 605, 705, 805, 905, 1005, 1105, 1205, 1305, 1405, 1505, + 1605, 1705, 0, 0, 106, 206, 306, 406, 506, 606, 706, 806, 906, 1006, 1106, 1206, 1306, + 1406, 1506, 1606, 1706, 0, 0, 107, 207, 307, 407, 507, 607, 707, 807, 907, 1007, 1107, + 1207, 1307, 1407, 1507, 1607, 1707, 0, 0, 108, 208, 308, 408, 508, 608, 708, 808, 908, + 1008, 1108, 1208, 1308, 1408, 1508, 1608, 1708, 0, 0, 109, 209, 309, 409, 509, 609, 709, + 809, 909, 1009, 1109, 1209, 1309, 1409, 1509, 1609, 1709, 0, 0, 110, 210, 310, 410, 510, + 610, 710, 810, 910, 1010, 1110, 1210, 1310, 1410, 1510, 0, 0, 0, 0, 111, 211, 311, 411, + 511, 611, 711, 811, 911, 1011, 1111, 1211, 1311, 1411, 1511, 0, 0 }), + std::vector({ 0, 0, 107, 207, 307, 407, 507, 607, 707, 807, 907, 1007, 1107, 1207, + 1307, 1407, 1507, 0, 0, 0, 0, 108, 208, 308, 408, 508, 608, 708, 808, 908, 1008, 1108, + 1208, 1308, 1408, 1508, 0, 0, 0, 0, 109, 209, 309, 409, 509, 609, 709, 809, 909, 1009, + 1109, 1209, 1309, 1409, 1509, 1609, 1709, 0, 0, 110, 210, 310, 410, 510, 610, 710, 810, + 910, 1010, 1110, 1210, 1310, 1410, 1510, 1610, 1710, 0, 0, 111, 211, 311, 411, 511, 611, + 711, 811, 911, 1011, 1111, 1211, 1311, 1411, 1511, 1611, 1711, 0, 0, 112, 212, 312, 412, + 512, 612, 712, 812, 912, 1012, 1112, 1212, 1312, 1412, 1512, 1612, 1712, 0, 0, 113, 213, + 313, 413, 513, 613, 713, 813, 913, 1013, 1113, 1213, 1313, 1413, 1513, 1613, 1713, 0, 0, + 114, 214, 314, 414, 514, 614, 714, 814, 914, 1014, 1114, 1214, 1314, 1414, 1514, 1614, + 1714, 0, 0, 115, 215, 315, 415, 515, 615, 715, 815, 915, 1015, 1115, 1215, 1315, 1415, + 1515, 1615, 1715, 0, 0, 116, 216, 316, 416, 516, 616, 716, 816, 916, 1016, 1116, 1216, + 1316, 1416, 1516, 1616, 1716, 0, 0, 117, 217, 317, 417, 517, 617, 717, 817, 917, 1017, + 1117, 1217, 1317, 1417, 1517, 1617, 1717, 0, 0, 118, 218, 318, 418, 518, 618, 718, 818, + 918, 1018, 1118, 1218, 1318, 1418, 1518, 1618, 1718, 0, 0, 119, 219, 319, 419, 519, 619, + 719, 819, 919, 1019, 1119, 1219, 1319, 1419, 1519, 1619, 1719, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0 }), + std::vector({ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 1301, 1401, 1501, 1601, 1701, 1801, 1901, 2001, 2101, 0, 0, 1302, 1402, 1502, 1602, + 1702, 1802, 1902, 2002, 2102, 0, 0, 1303, 1403, 1503, 1603, 1703, 1803, 1903, 2003, + 2103, 0, 0, 1304, 1404, 1504, 1604, 1704, 1804, 1904, 2004, 2104, 0, 0, 1305, 1405, + 1505, 1605, 1705, 1805, 1905, 2005, 2105, 0, 0, 1306, 1406, 1506, 1606, 1706, 1806, + 1906, 2006, 2106, 0, 0, 1307, 1407, 1507, 1607, 1707, 1807, 1907, 2007, 2107, 0, 0, + 1308, 1408, 1508, 1608, 1708, 1808, 1908, 2008, 2108, 0, 0, 1309, 1409, 1509, 1609, + 1709, 1809, 1909, 2009, 2109, 0, 0, 1310, 1410, 1510, 1610, 1710, 1810, 1910, 2010, + 2110, 0, 0, 1311, 1411, 1511, 1611, 1711, 1811, 1911, 2011, 2111, 0, 0, 1312, 1412, + 1512, 1612, 1712, 1812, 1912, 2012, 2112, 0, 0, 1313, 1413, 1513, 1613, 1713, 1813, + 1913, 2013, 2113, 0, 0, 1314, 1414, 1514, 1614, 1714, 1814, 1914, 2014, 2114, 0, 0, + 1315, 1415, 1515, 1615, 1715, 1815, 1915, 2015, 2115, 0, 0, 1316, 1416, 1516, 1616, + 1716, 1816, 1916, 2016, 2116, 0, 0, 1317, 1417, 1517, 1617, 1717, 1817, 1917, 2017, + 2117, 0, 0, 1318, 1418, 1518, 1618, 1718, 1818, 1918, 2018, 2118, 0, 0, 1319, 1419, + 1519, 1619, 1719, 1819, 1919, 2019, 2119, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0 }), + }; + + for (size_t j = 0; j < CGNy; j++) { + for (size_t i = 0; i < CGNx; i++) { + REQUIRE(cgVector(i + j * CGNx) == cgRefDataAllProcs[test_rank][i + j * CGNx]); + } + } } } diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 4aa2e8a6f..e4a1f499e 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -19,6 +19,10 @@ #include +#ifdef USE_MPI +#include "include/Halo.hpp" +#endif + namespace Nextsim { template @@ -258,6 +262,14 @@ template void CGDynamicsKernel::prepareIteration( Interpolations::DG2CG(*smesh, cgA, data.at(ciceName)); VectorManipulations::CGAveragePeriodic(*smesh, cgA); +#ifdef USE_MPI + // Halo object only depends on shape of the array, so we can share it for all DGVectors of the + // same shape + Halo halo(cgH); + halo.exchangeHalos(cgH); + halo.exchangeHalos(cgA); +#endif + // Reinit the gradient of the sea surface height. Not done by // DataMap as seaSurfaceHeight is always dG(0) computeGradientOfSeaSurfaceHeight(DynamicsKernel::seaSurfaceHeight); diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index 9dedaca46..fbd70456e 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -15,6 +15,10 @@ #include "include/constants.hpp" #include +#ifdef USE_MPI +#include "include/Halo.hpp" +#endif + namespace Nextsim { // The brittle momentum solver for CG velocity fields @@ -116,6 +120,20 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern { advectDynamicsFields(tst.step.seconds()); + // halo exchange + // - damage, hice, cice + // - s11, s12, s22 + // - e11, e12, e22 + // - dStressY, dStressX + // - u, v + // Note: only need to create one halo object per array type, dimensionality and size. +#ifdef USE_MPI + Halo halo(hice); + halo.exchangeHalos(hice.getDataRef()); + halo.exchangeHalos(cice.getDataRef()); + halo.exchangeHalos(damage.getDataRef()); +#endif + prepareIteration({ { hiceName, hice }, { ciceName, cice } }); // The timestep for the brittle solver is the solver subtimestep @@ -128,18 +146,42 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern projectVelocityToStrain(); +#ifdef USE_MPI + Halo haloDGVector(e11); + haloDGVector.exchangeHalos(e11); + haloDGVector.exchangeHalos(e12); + haloDGVector.exchangeHalos(e22); +#endif + std::array>, N_TENSOR_ELEMENTS> stress = { s11, s12, s22 }; // Call the step function on the StressUpdateStep class // Call the step function on the StressUpdateStep class stressStep.stressUpdateHighOrder( params, *smesh, stress, { e11, e12, e22 }, hice, cice, deltaT); +#ifdef USE_MPI + haloDGVector.exchangeHalos(s11); + haloDGVector.exchangeHalos(s12); + haloDGVector.exchangeHalos(s22); +#endif + stressDivergence(); // Compute divergence of stress tensor +#ifdef USE_MPI + Halo haloCGVector(dStressX); + haloCGVector.exchangeHalos(dStressX); + haloCGVector.exchangeHalos(dStressY); +#endif + updateMomentum(tst); applyBoundaries(); +#ifdef USE_MPI + haloCGVector.exchangeHalos(u); + haloCGVector.exchangeHalos(v); +#endif + // Land mask } @@ -190,6 +232,7 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern #pragma omp parallel for for (size_t i = 0; i < u.rows(); ++i) { + // TOM look here if (pmap->cglandmask(i) == 0) continue; diff --git a/dynamics/src/include/dgVectorHolder.hpp b/dynamics/src/include/dgVectorHolder.hpp index 8b87a3160..2b4fa6e2f 100644 --- a/dynamics/src/include/dgVectorHolder.hpp +++ b/dynamics/src/include/dgVectorHolder.hpp @@ -41,6 +41,8 @@ template class DGVectorHolder { void zero() { ref->setZero(); } + EigenDGVector& getDataRef() { return *ref; } + private: EigenDGVector* ref; };