diff --git a/core/src/ModelMetadata.cpp b/core/src/ModelMetadata.cpp index 79c981c0b..c971cbc3f 100644 --- a/core/src/ModelMetadata.cpp +++ b/core/src/ModelMetadata.cpp @@ -104,34 +104,105 @@ void ModelMetadata::readNeighbourData(netCDF::NcFile& ncFile) } } } + + std::array cornerSuffixes = { "_neighbour_ids", "_neighbour_send" }; + if (btype == periodic) { + continue; + for (auto& suffix : cornerSuffixes) { + suffix += "_periodic"; + } + } + + for (auto corner : corners) { + size_t nStart = 0; // start index in the netCDF variable + size_t count = 0; // number of elements to read for this rank + std::vector numCorners = std::vector(mpiSize, 0); + std::vector offsets = std::vector(mpiSize, 0); + std::vector>> arrays; + + // pick the correct set of corner arrays (periodic / non‑periodic) + if (btype == nonPeriodic) { + arrays = { cornerRanks[corner], cornerHaloSend[corner] }; + } else { // periodic + arrays = { cornerRanksPeriodic[corner], cornerHaloSendPeriodic[corner] }; + } + + // variable that stores *how many* corner entries each rank has + varName = cornerNames[corner] + "_neighbours"; + if (btype == periodic) { + varName += "_periodic"; + } + neighbourGroup.getVar(varName).getVar( + { 0 }, { static_cast(mpiSize) }, numCorners.data()); + + // compute the global offset for this rank + MPI_Exscan(&numCorners[mpiMyRank], &nStart, 1, MPI_INT, MPI_SUM, modelMPI.getComm()); + if (mpiMyRank == 0) { + nStart = 0; // MPI_Exscan undefined on rank 0 → set manually + } + count = numCorners[mpiMyRank]; + + if (count) { + // allocate and read each corner‑related array + for (size_t i = 0; i < arrays.size(); ++i) { + arrays[i].get().resize(count, 0); + varName = cornerNames[corner] + cornerSuffixes[i]; + neighbourGroup.getVar(varName).getVar( + { nStart }, { count }, arrays[i].get().data()); + } + } + } } } void ModelMetadata::getPartitionMetadata(std::string partitionFile) { - 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"; - throw std::runtime_error(errorMsg); - } - globalExtentX = ncFile.getDim("NX").getSize(); - globalExtentY = ncFile.getDim("NY").getSize(); - netCDF::NcGroup bboxGroup(ncFile.getGroup(bboxName)); + try { + 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"; + throw std::runtime_error(errorMsg); + } + globalExtentX = ncFile.getDim("NX").getSize(); + globalExtentY = ncFile.getDim("NY").getSize(); + netCDF::NcGroup bboxGroup(ncFile.getGroup(bboxName)); - std::vector 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); + std::vector 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); + readNeighbourData(ncFile); - ncFile.close(); + // cornerHaloRecv doesn't need to be read because it can be easily calculated. + for (auto corner : corners) { + if (cornerRanks[corner].size()) { + cornerHaloRecv[corner].resize(1); + cornerHaloRecv[corner][0] = 2 * (localExtentX + localExtentY) + corner; + } + } + + // 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(); + } catch (netCDF::exceptions::NcException& e) { + std::cerr << "Failed to open partition file [" << partitionFile << "] :: " << e.what() + << std::endl; + auto& modelMPI = ModelMPI::getInstance(); + MPI_Abort(modelMPI.getComm(), 1); + } } int ModelMetadata::getLocalCornerX() const { return localCornerX; } @@ -140,6 +211,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..4e8ea2c49 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,41 +137,26 @@ 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) + sendBufferSize = m_numHaloCells - 4 * nCells * nCells; + recvBufferSize = m_numHaloCells; send.resize(m_numComps); recv.resize(m_numComps); for (size_t i = 0; i < m_numComps; i++) { // allocate size and initialize to zero - send[i].resize(m_numHaloCells, 0.0); - recv[i].resize(m_numHaloCells, 0.0); + send[i].resize(sendBufferSize, 0.0); + recv[i].resize(recvBufferSize, 0.0); } // 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 @@ -144,6 +165,7 @@ class Halo { using Slice = ArraySlicer::Slice; using SliceIter = ArraySlicer::SliceIter; using Edge = ModelMetadata::Edge; + using Corner = ModelMetadata::Corner; using VBounds = ArraySlicer::Slice::VBounds; const typedef ArraySlicer::SliceIter::MultiDim MultiDim; @@ -155,17 +177,12 @@ 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 }, - { Edge::TOP, Edge::BOTTOM }, - { Edge::BOTTOM, Edge::TOP }, - }; // map to opposite edge + std::array corners = ModelMetadata::corners; // array of edge enums /** * @brief Return true if the provided edge is vertical (LEFT or RIGHT). @@ -191,13 +208,16 @@ 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; + size_t sendBufferSize = 0; + size_t recvBufferSize = 0; std::vector> send; // buffer to store halo region that will be read by other ranks std::vector> recv; // buffer to store halo region which is read from other ranks - ModelArray::DataType tempBuffer; /*! * @brief Open memory window to exchange send buffer between MPI ranks. @@ -209,7 +229,7 @@ class Halo { { // create a RMA memory window which all ranks will be able to access auto& modelMPI = ModelMPI::getInstance(); - MPI_Win_create(&send[idx][0], m_numHaloCells * sizeof(double), sizeof(double), + MPI_Win_create(&send[idx][0], sendBufferSize * sizeof(double), sizeof(double), MPI_INFO_NULL, modelMPI.getComm(), &m_win); // remove fence and check that no proceding RMA calls have been made MPI_Win_fence(MPI_MODE_NOPRECEDE, m_win); @@ -234,35 +254,145 @@ 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 = sendBufferSize; + if (isCG) { + buffer_len = sendBufferSize / 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) { + // populate edge data into send buffer + 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 which edge the send position sendPos would fall under in the send buffer + * + */ + Edge edgeFromSendPos(int sendPos, int fromRank) + { + auto& metadata = ModelMetadata::getInstance(); + + // extents of sending domain + auto extentX = metadata.getRankExtentsX()[fromRank]; + auto extentY = metadata.getRankExtentsY()[fromRank]; + + if (sendPos - (2 * extentX + extentY) >= 0) { + return Edge::LEFT; + } else if (sendPos - (extentX + extentY) >= 0) { + return Edge::TOP; + } else if (sendPos - extentX >= 0) { + return Edge::RIGHT; + } else { + return Edge::BOTTOM; + } + } + + /** + * @brief Calculate recv buffer positions and offsets for halo transfer + * + */ + 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]; + } + auto sendEdge = edgeFromSendPos(disp, fromRank); + if (isVertex) { + count = count + 1; + disp = disp + sendEdge; + recvOffset = recvOffset + edge; + } + if (isCG) { + count = CGdegree * count + 1; + disp = (disp > 0) ? CGdegree * disp + sendEdge : 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 + recvBufferSize / 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; + // TODO this is too big. Should only be 2 x CGdegree * haloWidth (extentX + extentY) + auto fromRankSendBufferSize = 2 * haloWidth * CGdegree * (extentX + extentY); + disp = disp + fromRankSendBufferSize / nCells * cell; + } + } + + /** + * @brief Calculate recv buffer positions and offsets for halo transfer of corner cells + */ + void recvPositions(int& fromRank, size_t& count, size_t& disp, size_t& recvOffset, + Corner corner, const size_t cell, bool isPeriodic) + { + count = 1; // we only have maximum of 1 corner neighbour for each corner + auto& metadata = ModelMetadata::getInstance(); + if (isPeriodic) { + // fromRank = metadata.neighbourRanksPeriodic[corner][neighbourIndex]; + // disp = metadata.neighbourHaloSendPeriodic[corner][neighbourIndex]; + // recvOffset = metadata.neighbourHaloRecvPeriodic[corner][neighbourIndex]; + } else { + fromRank = metadata.cornerRanks[corner][0]; + disp = metadata.cornerHaloSend[corner][0]; + recvOffset = metadata.cornerHaloRecv[corner][0]; + } + auto sendEdge = edgeFromSendPos(disp, fromRank); + if (isVertex) { + count = 1; + disp = disp + sendEdge; + recvOffset = recvOffset + 4; + // this is messy + // Essentially account for the fact that the vertex field is split differently to the + // face centered fields + if ((sendEdge == Edge::TOP or sendEdge == Edge::BOTTOM) + and (corner == Corner::TOP_RIGHT or corner == Corner::BOTTOM_RIGHT)) { + disp = disp + 1; + } + if ((sendEdge == Edge::LEFT or sendEdge == Edge::RIGHT) + and (corner == Corner::TOP_RIGHT or corner == Corner::TOP_LEFT)) { + disp = disp + 1; + } + } + if (isCG) { + count = CGdegree * count; + disp = (disp > 0) ? CGdegree * disp + sendEdge : 0; + recvOffset = CGdegree * recvOffset + 4; + + // recvOffset is the offset in the recv buffer and this belongs to the current rank + recvOffset = recvOffset + recvBufferSize / 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 fromRankSendBufferSize = 2 * haloWidth * CGdegree * (extentX + extentY); + disp = disp + fromRankSendBufferSize / nCells * cell; + if ((sendEdge == Edge::TOP or sendEdge == Edge::BOTTOM) + and (corner == Corner::TOP_RIGHT or corner == Corner::BOTTOM_RIGHT)) { + disp = disp + 1; + } + if ((sendEdge == Edge::LEFT or sendEdge == Edge::RIGHT) + and (corner == Corner::TOP_RIGHT or corner == Corner::TOP_LEFT)) { + disp = disp + 1; + } } } @@ -273,83 +403,69 @@ 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); } + } + } + } + + // get non-periodic corner neighbours and populate recv buffer (if the exist) + for (auto corner : corners) { + + // get neighbours (if they exist) + auto hasCorner = metadata.cornerRanks[corner].size(); + // hasCorner will either be 0 or 1 + if (hasCorner) { + for (size_t cell = 0; cell < nCells; ++cell) { + recvPositions(fromRank, count, disp, recvOffset, corner, cell, false); MPI_Get(&recv[comp][recvOffset], count, MPI_DOUBLE, fromRank, disp, count, MPI_DOUBLE, m_win); } } + + // get periodic neighbours (if they exist) + hasCorner = metadata.cornerRanksPeriodic[corner].size(); + if (hasCorner) { + // TODO implement periodic corner neighbours + throw std::runtime_error("Periodic Corner Neighbours not implemented yet."); + } } // close memory window (essentially make sure all communications are done before // 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 +489,175 @@ 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 recvBufferToTarget(const Corner corner, Eigen::Map& buffer_map, + Eigen::Map>& target_map) + { + int perimeter = 2 * (m_innerNx + m_innerNy); + int offset = perimeter + corner * nCells; + + switch (corner) { + case Corner::TOP_LEFT: + target_map(Eigen::seq(0, nCells - 1), Eigen::seq(Eigen::last - nCells + 1, Eigen::last)) + = buffer_map(Eigen::seq(offset, offset + nCells - 1), Eigen::all); + break; + case Corner::TOP_RIGHT: + target_map(Eigen::seq(Eigen::last - nCells + 1, Eigen::last), + Eigen::seq(Eigen::last - nCells + 1, Eigen::last)) + = buffer_map(Eigen::seq(offset, offset + nCells - 1), Eigen::all); + break; + case Corner::BOTTOM_RIGHT: + target_map(Eigen::seq(Eigen::last - nCells + 1, Eigen::last), Eigen::seq(0, nCells - 1)) + = buffer_map(Eigen::seq(offset, offset + nCells - 1), Eigen::all); + break; + case Corner::BOTTOM_LEFT: + target_map(Eigen::seq(0, nCells - 1), Eigen::seq(0, nCells - 1)) + = buffer_map(Eigen::seq(offset, offset + nCells - 1), Eigen::all); + break; + default: + throw std::runtime_error("Unrecognised corner type"); + } + } + + /* + * @brief transpose corners received from vertical edges + * + * If the sending edge was from a vertical side (e.g., LEFT or RIGHT) then it has been + * transposed when flattened into the 1D send buffer. We need to account for that. + * + * */ + void transposeCorners() + { + auto& metadata = ModelMetadata::getInstance(); + int fromRank; + size_t disp; + for (auto corner : corners) { + auto hasCorner = metadata.cornerRanks[corner].size(); + if (hasCorner) { + // if (isPeriodic) { + // // fromRank = metadata.neighbourRanksPeriodic[corner][neighbourIndex]; + // // disp = metadata.neighbourHaloSendPeriodic[corner][neighbourIndex]; + // } else { + fromRank = metadata.cornerRanks[corner][0]; + disp = metadata.cornerHaloSend[corner][0]; + + auto sendEdge = edgeFromSendPos(disp, fromRank); + + if (isVertical(sendEdge)) { + auto buffer_len = recvBufferSize / CGdegree; + for (size_t comp = 0; comp < m_numComps; ++comp) { + Eigen::Map rmap(recv[comp].data(), buffer_len, nCells); + auto offset = buffer_len - (4 - corner) * nCells; + rmap.block(offset, 0, nCells, nCells).transposeInPlace(); + } + } + } + } + } + + template void populateTarget(T& target) + { + size_t buffer_len = recvBufferSize; + if (isCG) { + buffer_len = recvBufferSize / 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); + } + for (auto corner : corners) { + recvBufferToTarget(corner, buffer_map, target_map); + } + } + } + /*! * @brief Get inner block of ModelArray from tempData * @@ -402,47 +687,14 @@ 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); - } + if (isCG) { + transposeCorners(); } + populateTarget(target); } }; } // end of nextsim namespace diff --git a/core/src/include/ModelMetadata.hpp b/core/src/include/ModelMetadata.hpp index 4b4855e35..c50c8a273 100644 --- a/core/src/include/ModelMetadata.hpp +++ b/core/src/include/ModelMetadata.hpp @@ -181,12 +181,28 @@ 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. static constexpr std::array edges = { BOTTOM, RIGHT, TOP, LEFT }; std::array edgeNames = { "bottom", "right", "top", "left" }; + enum Corner { TOP_LEFT, TOP_RIGHT, BOTTOM_RIGHT, BOTTOM_LEFT, N_CORNER }; + static constexpr std::array corners + = { TOP_LEFT, TOP_RIGHT, BOTTOM_RIGHT, BOTTOM_LEFT }; + std::array cornerNames + = { "top_left", "top_right", "bottom_right", "bottom_left" }; + typedef std::array, N_EDGE> neighbourArray; neighbourArray neighbourRanks; neighbourArray neighbourExtents; @@ -196,6 +212,13 @@ class ModelMetadata { neighbourArray neighbourExtentsPeriodic; neighbourArray neighbourHaloSendPeriodic; neighbourArray neighbourHaloRecvPeriodic; + + typedef std::array, N_CORNER> cornerArray; + cornerArray cornerRanks; + cornerArray cornerHaloSend; + cornerArray cornerHaloRecv; + cornerArray cornerRanksPeriodic; + cornerArray cornerHaloSendPeriodic; #endif private: @@ -235,6 +258,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/test/HaloExchangeCB_test.cpp b/core/test/HaloExchangeCB_test.cpp index 3d5ffc7c5..87adbf2e5 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,17 +25,18 @@ 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 static std::array, 3> refDataAllProcs = { std::vector({ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 100, 110, 120, 130, 140, 150, 160, 370, 0, 101, 111, 121, 131, 141, 151, 161, 371, 0, 102, 112, 122, 132, 142, 152, 162, 372, 0, 103, - 113, 123, 133, 143, 153, 163, 373, 0, 204, 214, 224, 234, 244, 254, 264, 0 }), - std::vector({ 0, 103, 113, 123, 133, 143, 153, 163, 0, 0, 204, 214, 224, 234, 244, 254, - 264, 374, 0, 205, 215, 225, 235, 245, 255, 265, 375, 0, 206, 216, 226, 236, 246, 256, 266, - 376, 0, 207, 217, 227, 237, 247, 257, 267, 377, 0, 208, 218, 228, 238, 248, 258, 268, 378, - 0, 0, 0, 0, 0, 0, 0, 0, 0 }), + 113, 123, 133, 143, 153, 163, 373, 0, 204, 214, 224, 234, 244, 254, 264, 374 }), + std::vector({ 0, 103, 113, 123, 133, 143, 153, 163, 373, 0, 204, 214, 224, 234, 244, + 254, 264, 374, 0, 205, 215, 225, 235, 245, 255, 265, 375, 0, 206, 216, 226, 236, 246, 256, + 266, 376, 0, 207, 217, 227, 237, 247, 257, 267, 377, 0, 208, 218, 228, 238, 248, 258, 268, + 378, 0, 0, 0, 0, 0, 0, 0, 0, 0 }), std::vector({ 0, 0, 0, 0, 0, 160, 370, 380, 390, 0, 161, 371, 381, 391, 0, 162, 372, 382, 392, 0, 163, 373, 383, 393, 0, 264, 374, 384, 394, 0, 265, 375, 385, 395, 0, 266, 376, 386, 396, 0, 267, 377, 387, 397, 0, 268, 378, 388, 398, 0, 0, 0, 0, 0, 0 }), @@ -145,20 +148,22 @@ MPI_TEST_CASE("test halo exchange on 3 proc grid", 3) -200., -450., -200., -350., -200., -250., -200., -150., -200., -50., -200., 50., -200., 150., -200., 250., -200., 0., 0., -550., -100., -450., -100., -350., -100., -250., -100., -150., -100., -50., -100., 50., -100., 150., -100., 250., -100., 0., 0., -550., - 0., -450., 0., -350., 0., -250., 0., -150., 0., -50., 0., 50., 0., 150., 0., 0., 0. }); + 0., -450., 0., -350., 0., -250., 0., -150., 0., -50., 0., 50., 0., 150., 0., 250., + 0. }); } if (test_rank == 1) { refDataVertex = std::vector({ 0., 0., -550., -200., -450., -200., -350., -200., - -250., -200., -150., -200., -50., -200., 50., -200., 150., -200., 0., 0., 0., 0., -550., - -100., -450., -100., -350., -100., -250., -100., -150., -100., -50., -100., 50., -100., - 150., -100., 250., -100., 0., 0., -550., 0., -450., 0., -350., 0., -250., 0., -150., 0., - -50., 0., 50., 0., 150., 0., 250., 0., 0., 0., -550., 100., -450., 100., -350., 100., - -250., 100., -150., 100., -50., 100., 50., 100., 150., 100., 250., 100., 0., 0., -550., - 200., -450., 200., -350., 200., -250., 200., -150., 200., -50., 200., 50., 200., 150., - 200., 250., 200., 0., 0., -550., 300., -450., 300., -350., 300., -250., 300., -150., - 300., -50., 300., 50., 300., 150., 300., 250., 300., 0., 0., -550., 400., -450., 400., - -350., 400., -250., 400., -150., 400., -50., 400., 50., 400., 150., 400., 250., 400., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. }); + -250., -200., -150., -200., -50., -200., 50., -200., 150., -200., 250., -200., 0., 0., + -550., -100., -450., -100., -350., -100., -250., -100., -150., -100., -50., -100., 50., + -100., 150., -100., 250., -100., 0., 0., -550., 0., -450., 0., -350., 0., -250., 0., + -150., 0., -50., 0., 50., 0., 150., 0., 250., 0., 0., 0., -550., 100., -450., 100., + -350., 100., -250., 100., -150., 100., -50., 100., 50., 100., 150., 100., 250., 100., + 0., 0., -550., 200., -450., 200., -350., 200., -250., 200., -150., 200., -50., 200., + 50., 200., 150., 200., 250., 200., 0., 0., -550., 300., -450., 300., -350., 300., -250., + 300., -150., 300., -50., 300., 50., 300., 150., 300., 250., 300., 0., 0., -550., 400., + -450., 400., -350., 400., -250., 400., -150., 400., -50., 400., 50., 400., 150., 400., + 250., 400., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0. }); } if (test_rank == 2) { refDataVertex = std::vector({ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 50., @@ -258,13 +263,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 +309,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, 1610, 1710, 0, 0, 111, 211, 311, + 411, 511, 611, 711, 811, 911, 1011, 1111, 1211, 1311, 1411, 1511, 1611, 1711 }), + std::vector({ 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, 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/core/test/ModelMetadataCB_test.cpp b/core/test/ModelMetadataCB_test.cpp index 1e1493ed2..fb2234fb0 100644 --- a/core/test/ModelMetadataCB_test.cpp +++ b/core/test/ModelMetadataCB_test.cpp @@ -23,6 +23,11 @@ constexpr ModelMetadata::Edge RIGHT = ModelMetadata::Edge::RIGHT; constexpr ModelMetadata::Edge TOP = ModelMetadata::Edge::TOP; constexpr ModelMetadata::Edge LEFT = ModelMetadata::Edge::LEFT; +constexpr ModelMetadata::Corner TOP_LEFT = ModelMetadata::Corner::TOP_LEFT; +constexpr ModelMetadata::Corner TOP_RIGHT = ModelMetadata::Corner::TOP_RIGHT; +constexpr ModelMetadata::Corner BOTTOM_RIGHT = ModelMetadata::Corner::BOTTOM_RIGHT; +constexpr ModelMetadata::Corner BOTTOM_LEFT = ModelMetadata::Corner::BOTTOM_LEFT; + typedef std::vector vec; // these tests are the same for closed boundary conditions (BC) and peridic BC @@ -30,6 +35,7 @@ static void testNonPeriodicBC(int test_rank) { auto& meta = ModelMetadata::getInstance(); if (test_rank == 0) { + // edges REQUIRE(meta.neighbourRanks[LEFT].size() == 0); REQUIRE(meta.neighbourRanks[RIGHT] == vec { 2 }); REQUIRE(meta.neighbourExtents[RIGHT] == vec { 4 }); @@ -40,7 +46,15 @@ static void testNonPeriodicBC(int test_rank) REQUIRE(meta.neighbourExtents[TOP] == vec { 7 }); REQUIRE(meta.neighbourHaloSend[TOP] == vec { 0 }); REQUIRE(meta.neighbourHaloRecv[TOP] == vec { 11 }); + + // corners + REQUIRE(meta.cornerRanks[BOTTOM_LEFT].size() == 0); + REQUIRE(meta.cornerRanks[BOTTOM_RIGHT].size() == 0); + REQUIRE(meta.cornerRanks[TOP_LEFT].size() == 0); + REQUIRE(meta.cornerRanks[TOP_RIGHT] == vec { 2 }); + REQUIRE(meta.cornerHaloSend[TOP_RIGHT] == vec { 12 }); } else if (test_rank == 1) { + // edges REQUIRE(meta.neighbourRanks[LEFT].size() == 0); REQUIRE(meta.neighbourRanks[RIGHT] == vec { 2 }); REQUIRE(meta.neighbourExtents[RIGHT] == vec { 5 }); @@ -51,7 +65,15 @@ static void testNonPeriodicBC(int test_rank) REQUIRE(meta.neighbourHaloSend[BOTTOM] == vec { 11 }); REQUIRE(meta.neighbourHaloRecv[BOTTOM] == vec { 0 }); REQUIRE(meta.neighbourRanks[TOP].size() == 0); + + // corners + REQUIRE(meta.cornerRanks[BOTTOM_LEFT].size() == 0); + REQUIRE(meta.cornerRanks[BOTTOM_RIGHT] == vec { 2 }); + REQUIRE(meta.cornerHaloSend[BOTTOM_RIGHT] == vec { 9 }); + REQUIRE(meta.cornerRanks[TOP_LEFT].size() == 0); + REQUIRE(meta.cornerRanks[TOP_RIGHT].size() == 0); } else if (test_rank == 2) { + // edges REQUIRE(meta.neighbourRanks[LEFT] == vec { 0, 1 }); REQUIRE(meta.neighbourExtents[LEFT] == vec { 4, 5 }); REQUIRE(meta.neighbourHaloSend[LEFT] == vec { 7, 7 }); @@ -59,6 +81,12 @@ static void testNonPeriodicBC(int test_rank) REQUIRE(meta.neighbourRanks[RIGHT].size() == 0); REQUIRE(meta.neighbourRanks[BOTTOM].size() == 0); REQUIRE(meta.neighbourRanks[TOP].size() == 0); + + // corners + REQUIRE(meta.cornerRanks[BOTTOM_LEFT].size() == 0); + REQUIRE(meta.cornerRanks[BOTTOM_RIGHT].size() == 0); + REQUIRE(meta.cornerRanks[TOP_LEFT].size() == 0); + REQUIRE(meta.cornerRanks[TOP_RIGHT].size() == 0); } } @@ -78,6 +106,10 @@ MPI_TEST_CASE("Test getPartitionMetadata closed boundary", 3) REQUIRE(meta.neighbourRanksPeriodic[RIGHT].size() == 0); REQUIRE(meta.neighbourRanksPeriodic[BOTTOM].size() == 0); REQUIRE(meta.neighbourRanksPeriodic[TOP].size() == 0); + REQUIRE(meta.cornerRanksPeriodic[LEFT].size() == 0); + REQUIRE(meta.cornerRanksPeriodic[RIGHT].size() == 0); + REQUIRE(meta.cornerRanksPeriodic[BOTTOM].size() == 0); + REQUIRE(meta.cornerRanksPeriodic[TOP].size() == 0); } } diff --git a/core/test/partition_metadata_3_cb.cdl b/core/test/partition_metadata_3_cb.cdl index 789136a45..e6e2d8aaf 100644 --- a/core/test/partition_metadata_3_cb.cdl +++ b/core/test/partition_metadata_3_cb.cdl @@ -1,4 +1,4 @@ -netcdf partition_metadata_3_cb { +netcdf partition_metadata_3 { dimensions: NX = 10 ; NY = 9 ; @@ -7,10 +7,18 @@ dimensions: R = 2 ; B = 1 ; T = 1 ; + TL = UNLIMITED ; // (0 currently) + TR = 1 ; + BR = 1 ; + BL = UNLIMITED ; // (0 currently) L_periodic = UNLIMITED ; // (0 currently) R_periodic = UNLIMITED ; // (0 currently) B_periodic = UNLIMITED ; // (0 currently) T_periodic = UNLIMITED ; // (0 currently) + TL_periodic = UNLIMITED ; // (0 currently) + TR_periodic = UNLIMITED ; // (0 currently) + BR_periodic = UNLIMITED ; // (0 currently) + BL_periodic = UNLIMITED ; // (0 currently) group: bounding_boxes { variables: @@ -51,6 +59,18 @@ group: connectivity { int top_neighbour_halos(T) ; int top_neighbour_halo_send(T) ; int top_neighbour_halo_recv(T) ; + int top_left_neighbours(P) ; + int top_left_neighbour_ids(TL) ; + int top_left_neighbour_send(TL) ; + int top_right_neighbours(P) ; + int top_right_neighbour_ids(TR) ; + int top_right_neighbour_send(TR) ; + int bottom_right_neighbours(P) ; + int bottom_right_neighbour_ids(BR) ; + int bottom_right_neighbour_send(BR) ; + int bottom_left_neighbours(P) ; + int bottom_left_neighbour_ids(BL) ; + int bottom_left_neighbour_send(BL) ; int left_neighbours_periodic(P) ; int left_neighbour_ids_periodic(L_periodic) ; int left_neighbour_halos_periodic(L_periodic) ; @@ -71,6 +91,18 @@ group: connectivity { int top_neighbour_halos_periodic(T_periodic) ; int top_neighbour_halo_send_periodic(T_periodic) ; int top_neighbour_halo_recv_periodic(T_periodic) ; + int top_left_neighbours_periodic(P) ; + int top_left_neighbour_ids_periodic(TL_periodic) ; + int top_left_neighbour_send_periodic(TL_periodic) ; + int top_right_neighbours_periodic(P) ; + int top_right_neighbour_ids_periodic(TR_periodic) ; + int top_right_neighbour_send_periodic(TR_periodic) ; + int bottom_right_neighbours_periodic(P) ; + int bottom_right_neighbour_ids_periodic(BR_periodic) ; + int bottom_right_neighbour_send_periodic(BR_periodic) ; + int bottom_left_neighbours_periodic(P) ; + int bottom_left_neighbour_ids_periodic(BL_periodic) ; + int bottom_left_neighbour_send_periodic(BL_periodic) ; data: left_neighbours = 0, 0, 2 ; @@ -113,6 +145,22 @@ group: connectivity { top_neighbour_halo_recv = 11 ; + top_left_neighbours = 0, 0, 0 ; + + top_right_neighbours = 1, 0, 0 ; + + top_right_neighbour_ids = 2 ; + + top_right_neighbour_send = 19 ; + + bottom_right_neighbours = 0, 1, 0 ; + + bottom_right_neighbour_ids = 2 ; + + bottom_right_neighbour_send = 18 ; + + bottom_left_neighbours = 0, 0, 0 ; + left_neighbours_periodic = 0, 0, 0 ; right_neighbours_periodic = 0, 0, 0 ; @@ -120,5 +168,13 @@ group: connectivity { bottom_neighbours_periodic = 0, 0, 0 ; top_neighbours_periodic = 0, 0, 0 ; + + top_left_neighbours_periodic = 0, 0, 0 ; + + top_right_neighbours_periodic = 0, 0, 0 ; + + bottom_right_neighbours_periodic = 0, 0, 0 ; + + bottom_left_neighbours_periodic = 0, 0, 0 ; } // group connectivity } 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..e3ec32a55 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(static_cast&>(hice)); + halo.exchangeHalos(static_cast&>(cice)); + halo.exchangeHalos(static_cast&>(damage)); +#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;