Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions core/src/ModelMetadata.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();
}

Expand All @@ -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<int> ModelMetadata::getRankExtentsX() const { return rankExtentsX; }
std::vector<int> ModelMetadata::getRankExtentsY() const { return rankExtentsY; }
#else

ModelMetadata::ModelMetadata()
Expand Down
362 changes: 220 additions & 142 deletions core/src/include/Halo.hpp

Large diffs are not rendered by default.

12 changes: 12 additions & 0 deletions core/src/include/ModelMetadata.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<int> 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<int> getRankExtentsY() const;

enum Edge { BOTTOM, RIGHT, TOP, LEFT, N_EDGE };
// An array to allow the edges to be accessed in the correct order.
Expand Down Expand Up @@ -235,6 +245,8 @@ class ModelMetadata {
int localCornerX, localCornerY;
int localExtentX, localExtentY;
int globalExtentX, globalExtentY;
std::vector<int> rankExtentsX; // vector containing domain extents for each rank x-direction
std::vector<int> rankExtentsY; // vector containing domain extents for each rank y-direction
const std::string bboxName = "bounding_boxes";
const std::string neighbourName = "connectivity";
#endif
Expand Down
10 changes: 10 additions & 0 deletions core/src/modules/DynamicsModule/BBMDynamics.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,9 @@
*/

#include "include/BBMDynamics.hpp"
#ifdef USE_MPI
#include "include/ModelMPI.hpp"
#endif
#include "include/constants.hpp"
#include "include/gridNames.hpp"

Expand Down Expand Up @@ -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;
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we should probably use a logger here instead. (Can we see if boost log can handle MPI? or can we make it handle MPI?)

}
#else
std::cout << tst.start << std::endl;
#endif

// set the forcing velocities
kernel.setData(uWindName, uwind);
Expand Down
10 changes: 10 additions & 0 deletions core/src/modules/DynamicsModule/MEVPDynamics.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,9 @@
*/

#include "include/MEVPDynamics.hpp"
#ifdef USE_MPI
#include "include/ModelMPI.hpp"
#endif
#include "include/constants.hpp"
#include "include/gridNames.hpp"

Expand Down Expand Up @@ -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);
Expand Down
10 changes: 10 additions & 0 deletions core/src/modules/DynamicsModule/include/FreeDriftDynamics.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -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);
Expand Down
134 changes: 128 additions & 6 deletions core/test/HaloExchangeCB_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,18 @@

#include "ModelMPI.hpp"
#include "ModelMetadata.hpp"
#include "cgVector.hpp"
#include "include/DGModelArray.hpp"
#include "include/Halo.hpp"
#include "include/Interpolations.hpp"

namespace Nextsim {

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
Expand Down Expand Up @@ -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;
}
}

Expand Down Expand Up @@ -304,5 +307,124 @@ MPI_TEST_CASE("DGVector", 3)
}
}
}

// initialize CGVector (after halo exchange on DGVector)
CGVector<CGDEGREE> 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<CGDEGREE> 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<std::vector<double>, 3> cgRefDataAllProcs = {
std::vector<double>({ 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<double>({ 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<double>({ 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]);
}
}
}
}
12 changes: 12 additions & 0 deletions dynamics/src/CGDynamicsKernel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,10 @@

#include <limits>

#ifdef USE_MPI
#include "include/Halo.hpp"
#endif

namespace Nextsim {

template <int DGadvection>
Expand Down Expand Up @@ -258,6 +262,14 @@ template <int DGadvection> void CGDynamicsKernel<DGadvection>::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<DGadvection, DGstressComp>::seaSurfaceHeight);
Expand Down
43 changes: 43 additions & 0 deletions dynamics/src/include/BrittleCGDynamicsKernel.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,10 @@
#include "include/constants.hpp"
#include <cmath>

#ifdef USE_MPI
#include "include/Halo.hpp"
#endif

namespace Nextsim {

// The brittle momentum solver for CG velocity fields
Expand Down Expand Up @@ -116,6 +120,20 @@ template <int DGadvection> 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
Expand All @@ -128,18 +146,42 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern

projectVelocityToStrain();

#ifdef USE_MPI
Halo haloDGVector(e11);
haloDGVector.exchangeHalos(e11);
haloDGVector.exchangeHalos(e12);
haloDGVector.exchangeHalos(e22);
#endif

std::array<std::reference_wrapper<DGVector<DGstressComp>>, 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
}

Expand Down Expand Up @@ -190,6 +232,7 @@ template <int DGadvection> 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;

Expand Down
Loading
Loading