Skip to content

Commit

Permalink
Add edge sign computations
Browse files Browse the repository at this point in the history
  • Loading branch information
@sbrus89
sbrus89 committed Jan 26, 2024
1 parent 5633b40 commit 2f5bff7
Show file tree
Hide file tree
Showing 3 changed files with 155 additions and 30 deletions.
101 changes: 71 additions & 30 deletions components/omega/src/ocn/HorzMesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,8 @@ HorzMesh::HorzMesh(Decomp *MeshDecomp) {
// TODO: add ability to compute (rather than read in)
// dependent mesh quantities

computeEdgeSign();

} // end constructor

//------------------------------------------------------------------------------
Expand Down Expand Up @@ -117,6 +119,10 @@ void HorzMesh::clear() {
VerticesOnEdge.deallocate();
CellsOnVertex.deallocate();
EdgesOnVertex.deallocate();
EdgeSignOnCell.deallocate();
EdgeSignOnCellH.deallocate();
EdgeSignOnVertex.deallocate();
EdgeSignOnVertexH.deallocate();

} // end clear

Expand Down Expand Up @@ -169,7 +175,7 @@ void HorzMesh::initParallelIO(Decomp *MeshDecomp) {
MaxEdges2 = 2 * MaxEdges;
std::vector<I4> OnEdgeDims2{MeshDecomp->NEdgesGlobal, MaxEdges2};
I4 OnEdgeSize2 = NEdgesAll * MaxEdges2;
std::vector<I4> OnEdgeOffset2(OnEdgeSize2, 0);
std::vector<I4> OnEdgeOffset2(OnEdgeSize2, -1);
for (int Edge = 0; Edge < NEdgesAll; Edge++) {
for (int i = 0; i < MaxEdges2; i++) {
I4 GlobalID = EdgeID[Edge] * MaxEdges2 + i;
Expand All @@ -186,7 +192,7 @@ void HorzMesh::initParallelIO(Decomp *MeshDecomp) {
// Create the IO decomp for arrays with (NVertices, VertexDegree) dimensions
std::vector<I4> OnVertexDims{MeshDecomp->NVerticesGlobal, VertexDegree};
I4 OnVertexSize = NVerticesAll * VertexDegree;
std::vector<I4> OnVertexOffset(OnVertexSize, 0);
std::vector<I4> OnVertexOffset(OnVertexSize, -1);
for (int Vertex = 0; Vertex < NVerticesAll; Vertex++) {
for (int i = 0; i < VertexDegree; i++) {
I4 GlobalID = VertexID[Vertex] * VertexDegree + i;
Expand All @@ -209,107 +215,107 @@ void HorzMesh::readCoordinates() {

// Read mesh cell coordinates
int XCellID;
XCellH = ArrayHost1DR8("xCell", NCellsAll);
XCellH = ArrayHost1DR8("XCell", NCellsAll);
Err = IO::readArray(XCellH.data(), NCellsAll, "xCell", MeshFileID,
CellDecompR8, XCellID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading xCell");

int YCellID;
YCellH = ArrayHost1DR8("yCell", NCellsAll);
YCellH = ArrayHost1DR8("YCell", NCellsAll);
Err = IO::readArray(YCellH.data(), NCellsAll, "yCell", MeshFileID,
CellDecompR8, YCellID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading yCell");

int ZCellID;
ZCellH = ArrayHost1DR8("zCell", NCellsAll);
ZCellH = ArrayHost1DR8("ZCell", NCellsAll);
Err = IO::readArray(ZCellH.data(), NCellsAll, "zCell", MeshFileID,
CellDecompR8, ZCellID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading zCell");

int LonCellID;
LonCellH = ArrayHost1DR8("lonCell", NCellsAll);
LonCellH = ArrayHost1DR8("LonCell", NCellsAll);
Err = IO::readArray(LonCellH.data(), NCellsAll, "lonCell", MeshFileID,
CellDecompR8, LonCellID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading lonCell");

int LatCellID;
LatCellH = ArrayHost1DR8("latCell", NCellsAll);
LatCellH = ArrayHost1DR8("LatCell", NCellsAll);
Err = IO::readArray(LatCellH.data(), NCellsAll, "latCell", MeshFileID,
CellDecompR8, LatCellID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading latCell");

// Read mesh edge coordinateID
int XEdgeID;
XEdgeH = ArrayHost1DR8("xEdge", NEdgesAll);
XEdgeH = ArrayHost1DR8("XEdge", NEdgesAll);
Err = IO::readArray(XEdgeH.data(), NEdgesAll, "xEdge", MeshFileID,
EdgeDecompR8, XEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading xEdge");

int YEdgeID;
YEdgeH = ArrayHost1DR8("yEdge", NEdgesAll);
YEdgeH = ArrayHost1DR8("YEdge", NEdgesAll);
Err = IO::readArray(YEdgeH.data(), NEdgesAll, "yEdge", MeshFileID,
EdgeDecompR8, YEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading yEdge");

int ZEdgeID;
ZEdgeH = ArrayHost1DR8("zEdge", NEdgesAll);
ZEdgeH = ArrayHost1DR8("ZEdge", NEdgesAll);
Err = IO::readArray(ZEdgeH.data(), NEdgesAll, "zEdge", MeshFileID,
EdgeDecompR8, ZEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading zEdge");

int LonEdgeID;
LonEdgeH = ArrayHost1DR8("lonEdge", NEdgesAll);
LonEdgeH = ArrayHost1DR8("LonEdge", NEdgesAll);
Err = IO::readArray(LonEdgeH.data(), NEdgesAll, "lonEdge", MeshFileID,
EdgeDecompR8, LonEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading lonEdge");

int LatEdgeID;
LatEdgeH = ArrayHost1DR8("latEdge", NEdgesAll);
LatEdgeH = ArrayHost1DR8("LatEdge", NEdgesAll);
Err = IO::readArray(LatEdgeH.data(), NEdgesAll, "latEdge", MeshFileID,
EdgeDecompR8, LatEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading latEdge");

// Read mesh vertex coordinates
int XVertexID;
XVertexH = ArrayHost1DR8("xVertex", NVerticesAll);
XVertexH = ArrayHost1DR8("XVertex", NVerticesAll);
Err = IO::readArray(XVertexH.data(), NVerticesAll, "xVertex", MeshFileID,
VertexDecompR8, XVertexID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading xVertex");

int YVertexID;
YVertexH = ArrayHost1DR8("yVertex", NVerticesAll);
YVertexH = ArrayHost1DR8("YVertex", NVerticesAll);
Err = IO::readArray(YVertexH.data(), NVerticesAll, "yVertex", MeshFileID,
VertexDecompR8, YVertexID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading yVertex");

int ZVertexID;
ZVertexH = ArrayHost1DR8("zVertex", NVerticesAll);
ZVertexH = ArrayHost1DR8("ZVertex", NVerticesAll);
Err = IO::readArray(ZVertexH.data(), NVerticesAll, "zVertex", MeshFileID,
VertexDecompR8, ZVertexID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading zVertex");

int LonVertexID;
LonVertexH = ArrayHost1DR8("lonVertex", NVerticesAll);
LonVertexH = ArrayHost1DR8("LonVertex", NVerticesAll);
Err = IO::readArray(LonVertexH.data(), NVerticesAll, "lonVertex", MeshFileID,
VertexDecompR8, LonVertexID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading lonVertex");

int LatVertexID;
LatVertexH = ArrayHost1DR8("latVertex", NVerticesAll);
LatVertexH = ArrayHost1DR8("LatVertex", NVerticesAll);
Err = IO::readArray(LatVertexH.data(), NVerticesAll, "latVertex", MeshFileID,
VertexDecompR8, LatVertexID);
if (Err != 0)
Expand All @@ -324,7 +330,7 @@ void HorzMesh::readBottomDepth() {
I4 Err;

int BottomDepthID;
BottomDepthH = ArrayHost1DR8("bottomDepth", NCellsAll);
BottomDepthH = ArrayHost1DR8("BottomDepth", NCellsAll);
Err = IO::readArray(BottomDepthH.data(), NCellsAll, "bottomDepth",
MeshFileID, CellDecompR8, BottomDepthID);
if (Err != 0)
Expand All @@ -340,50 +346,50 @@ void HorzMesh::readMeasurements() {
I4 Err;

int AreaCellID;
AreaCellH = ArrayHost1DR8("areaCell", NCellsAll);
AreaCellH = ArrayHost1DR8("AreaCell", NCellsAll);
Err = IO::readArray(AreaCellH.data(), NCellsAll, "areaCell", MeshFileID,
CellDecompR8, AreaCellID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading areaCell");

int AreaTriangleID;
AreaTriangleH = ArrayHost1DR8("areaTriangle", NVerticesAll);
AreaTriangleH = ArrayHost1DR8("AreaTriangle", NVerticesAll);
Err = IO::readArray(AreaTriangleH.data(), NVerticesAll, "areaTriangle",
MeshFileID, VertexDecompR8, AreaTriangleID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading areaTriangle");

int DvEdgeID;
DvEdgeH = ArrayHost1DR8("dvEdge", NEdgesAll);
DvEdgeH = ArrayHost1DR8("DvEdge", NEdgesAll);
Err = IO::readArray(DvEdgeH.data(), NEdgesAll, "dvEdge", MeshFileID,
EdgeDecompR8, DvEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading dvEdge");

int DcEdgeID;
DcEdgeH = ArrayHost1DR8("dcEdge", NEdgesAll);
DcEdgeH = ArrayHost1DR8("DcEdge", NEdgesAll);
Err = IO::readArray(DcEdgeH.data(), NEdgesAll, "dcEdge", MeshFileID,
EdgeDecompR8, DcEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading dcEdge");

int AngleEdgeID;
AngleEdgeH = ArrayHost1DR8("angleEdge", NEdgesAll);
AngleEdgeH = ArrayHost1DR8("AngleEdge", NEdgesAll);
Err = IO::readArray(AngleEdgeH.data(), NEdgesAll, "angleEdge", MeshFileID,
EdgeDecompR8, AngleEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading angleEdge");

int MeshDensityID;
MeshDensityH = ArrayHost1DR8("meshDensity", NCellsAll);
MeshDensityH = ArrayHost1DR8("MeshDensity", NCellsAll);
Err = IO::readArray(MeshDensityH.data(), NCellsAll, "meshDensity",
MeshFileID, CellDecompR8, MeshDensityID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading meshDensity");

int KiteAreasOnVertexID;
KiteAreasOnVertexH =
ArrayHost2DR8("kiteAreasOnVertex", NVerticesAll, VertexDegree);
ArrayHost2DR8("KiteAreasOnVertex", NVerticesAll, VertexDegree);
Err = IO::readArray(KiteAreasOnVertexH.data(), NVerticesAll * VertexDegree,
"kiteAreasOnVertex", MeshFileID, OnVertexDecompR8,
KiteAreasOnVertexID);
Expand All @@ -399,7 +405,7 @@ void HorzMesh::readWeights() {
I4 Err;

int WeightsOnEdgeID;
WeightsOnEdgeH = ArrayHost2DR8("weightsOnEdge", NEdgesAll, MaxEdges2);
WeightsOnEdgeH = ArrayHost2DR8("WeightsOnEdge", NEdgesAll, MaxEdges2);
Err = IO::readArray(WeightsOnEdgeH.data(), NEdgesAll * MaxEdges2,
"weightsOnEdge", MeshFileID, OnEdgeDecompR8,
WeightsOnEdgeID);
Expand All @@ -415,36 +421,71 @@ void HorzMesh::readCoriolis() {
int Err;

int FCellID;
FCellH = ArrayHost1DR8("fCell", NCellsAll);
FCellH = ArrayHost1DR8("FCell", NCellsAll);
Err = IO::readArray(FCellH.data(), NCellsAll, "fCell", MeshFileID,
CellDecompR8, FCellID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading fCell");

int FVertexID;
FVertexH = ArrayHost1DR8("fVertex", NVerticesAll);
FVertexH = ArrayHost1DR8("FVertex", NVerticesAll);
Err = IO::readArray(FVertexH.data(), NVerticesAll, "fVertex", MeshFileID,
VertexDecompR8, FVertexID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading fVertex");

int FEdgeID;
FEdgeH = ArrayHost1DR8("fEdge", NEdgesAll);
FEdgeH = ArrayHost1DR8("FEdge", NEdgesAll);
Err = IO::readArray(FEdgeH.data(), NEdgesAll, "fEdge", MeshFileID,
EdgeDecompR8, FEdgeID);
if (Err != 0)
LOG_CRITICAL("HorzMesh: error reading fEdge");

} // end readCoriolis

//------------------------------------------------------------------------------
// Compute the sign of edge contributions to a cell/vertex for each ecge
void HorzMesh::computeEdgeSign() {

EdgeSignOnCell = Array2DR8("EdgeSignOnCell", NCellsAll, MaxEdges);
yakl::c::parallel_for(yakl::c::SimpleBounds<1>(NCellsAll), YAKL_LAMBDA(int Cell) {
for (int i = 0; i < NEdgesOnCell(Cell); i++) {
int Edge = EdgesOnCell(Cell, i);

// Vector points from cell 0 to cell 1
if (Cell == CellsOnEdge(Edge, 0)) {
EdgeSignOnCell(Cell, i) = -1.0;
} else {
EdgeSignOnCell(Cell, i) = 1.0;
}
}
});
EdgeSignOnCellH = EdgeSignOnCell.createHostCopy();

EdgeSignOnVertex = Array2DR8("EdgeSignOnVertex", NVerticesAll, VertexDegree);
yakl::c::parallel_for(yakl::c::SimpleBounds<1>(NVerticesAll), YAKL_LAMBDA(int Vertex) {
for (int i = 0; i < VertexDegree; i++) {
int Edge = EdgesOnVertex(Vertex, i);

// Vector points from vertex 0 to vertex 1
if (Vertex == VerticesOnEdge(Edge, 0)) {
EdgeSignOnVertex(Vertex, i) = -1.0;
} else {
EdgeSignOnVertex(Vertex, i) = 1.0;
}
}
});
EdgeSignOnVertexH = EdgeSignOnVertex.createHostCopy();

} //end computeEdgeSign

//------------------------------------------------------------------------------
// Perform copy to device for mesh variables
void HorzMesh::copyToDevice() {

AreaCell = AreaCellH.createDeviceCopy();
AreaTriangle = AreaTriangleH.createDeviceCopy();
KiteAreasOnVertex = KiteAreasOnVertexH.createDeviceCopy();
DvEdge = DvEdgeH.createDeviceCopy();
DcEdge = DcEdgeH.createDeviceCopy();
AngleEdge = AngleEdgeH.createDeviceCopy();
WeightsOnEdge = WeightsOnEdgeH.createDeviceCopy();
Expand Down
10 changes: 10 additions & 0 deletions components/omega/src/ocn/HorzMesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ class HorzMesh {

void readCoriolis();

void computeEdgeSign();

void copyToDevice();

// int computeMesh();
Expand Down Expand Up @@ -184,6 +186,14 @@ class HorzMesh {
Array1DR8 BottomDepth; ///< Depth of the bottom of the ocean (m)
ArrayHost1DR8 BottomDepthH; ///< Depth of the bottom of the ocean (m)

// Edge sign

Array2DR8 EdgeSignOnCell; ///< Sign of vector connecting cells
ArrayHost2DR8 EdgeSignOnCellH; ///< Sign of vector connecting cells

Array2DR8 EdgeSignOnVertex; ///< Sign of vector connecting vertices
ArrayHost2DR8 EdgeSignOnVertexH; ///< Sign of vector connecting vertices

// Methods

/// Construct a new local mesh for a given decomposition
Expand Down
Loading

0 comments on commit 2f5bff7

Please sign in to comment.