From aa61e7f5e147f48b2da5c1bbcb176a9a3c8a9632 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 26 Jan 2024 10:35:24 -0800 Subject: [PATCH] Add edge sign computations --- components/omega/src/ocn/HorzMesh.cpp | 101 +++++++++++++++------ components/omega/src/ocn/HorzMesh.h | 10 ++ components/omega/test/ocn/HorzMeshTest.cpp | 74 +++++++++++++++ 3 files changed, 155 insertions(+), 30 deletions(-) diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index c1ea8986f2d6..d9a0a6e9f8c1 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -84,6 +84,8 @@ HorzMesh::HorzMesh(Decomp *MeshDecomp) { // TODO: add ability to compute (rather than read in) // dependent mesh quantities + computeEdgeSign(); + } // end constructor //------------------------------------------------------------------------------ @@ -117,6 +119,10 @@ void HorzMesh::clear() { VerticesOnEdge.deallocate(); CellsOnVertex.deallocate(); EdgesOnVertex.deallocate(); + EdgeSignOnCell.deallocate(); + EdgeSignOnCellH.deallocate(); + EdgeSignOnVertex.deallocate(); + EdgeSignOnVertexH.deallocate(); } // end clear @@ -169,7 +175,7 @@ void HorzMesh::initParallelIO(Decomp *MeshDecomp) { MaxEdges2 = 2 * MaxEdges; std::vector OnEdgeDims2{MeshDecomp->NEdgesGlobal, MaxEdges2}; I4 OnEdgeSize2 = NEdgesAll * MaxEdges2; - std::vector OnEdgeOffset2(OnEdgeSize2, 0); + std::vector OnEdgeOffset2(OnEdgeSize2, -1); for (int Edge = 0; Edge < NEdgesAll; Edge++) { for (int i = 0; i < MaxEdges2; i++) { I4 GlobalID = EdgeID[Edge] * MaxEdges2 + i; @@ -186,7 +192,7 @@ void HorzMesh::initParallelIO(Decomp *MeshDecomp) { // Create the IO decomp for arrays with (NVertices, VertexDegree) dimensions std::vector OnVertexDims{MeshDecomp->NVerticesGlobal, VertexDegree}; I4 OnVertexSize = NVerticesAll * VertexDegree; - std::vector OnVertexOffset(OnVertexSize, 0); + std::vector OnVertexOffset(OnVertexSize, -1); for (int Vertex = 0; Vertex < NVerticesAll; Vertex++) { for (int i = 0; i < VertexDegree; i++) { I4 GlobalID = VertexID[Vertex] * VertexDegree + i; @@ -209,35 +215,35 @@ 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) @@ -245,35 +251,35 @@ void HorzMesh::readCoordinates() { // 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) @@ -281,35 +287,35 @@ void HorzMesh::readCoordinates() { // 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) @@ -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) @@ -340,42 +346,42 @@ 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) @@ -383,7 +389,7 @@ void HorzMesh::readMeasurements() { int KiteAreasOnVertexID; KiteAreasOnVertexH = - ArrayHost2DR8("kiteAreasOnVertex", NVerticesAll, VertexDegree); + ArrayHost2DR8("KiteAreasOnVertex", NVerticesAll, VertexDegree); Err = IO::readArray(KiteAreasOnVertexH.data(), NVerticesAll * VertexDegree, "kiteAreasOnVertex", MeshFileID, OnVertexDecompR8, KiteAreasOnVertexID); @@ -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); @@ -415,21 +421,21 @@ 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) @@ -437,6 +443,42 @@ void HorzMesh::readCoriolis() { } // 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() { @@ -444,7 +486,6 @@ void HorzMesh::copyToDevice() { AreaCell = AreaCellH.createDeviceCopy(); AreaTriangle = AreaTriangleH.createDeviceCopy(); KiteAreasOnVertex = KiteAreasOnVertexH.createDeviceCopy(); - DvEdge = DvEdgeH.createDeviceCopy(); DcEdge = DcEdgeH.createDeviceCopy(); AngleEdge = AngleEdgeH.createDeviceCopy(); WeightsOnEdge = WeightsOnEdgeH.createDeviceCopy(); diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index 21b9be04da45..4c16b05e0b66 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -39,6 +39,8 @@ class HorzMesh { void readCoriolis(); + void computeEdgeSign(); + void copyToDevice(); // int computeMesh(); @@ -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 diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 9dabc4a5dcfa..c50cdf3f898c 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -554,6 +554,80 @@ int main(int argc, char *argv[]) { LOG_INFO("HorzMeshTest: weightsOnEdge test FAIL"); } + // Test edgeSignOnCell + // Check that the sign corresponds with convention + // Tests that the edge sign values were calculated correctly + count = 0; + for (int Edge = 0; Edge < LocEdges; Edge++) { + int Cell0 = Mesh.CellsOnEdgeH(Edge, 0); + int iEdge0; + for (int i = 0; i < Mesh.NEdgesOnCellH(Cell0); i++) { + if (Mesh.EdgesOnCellH(Cell0, i) == Edge) { + iEdge0 = i; + break; + } + } + if (abs(Mesh.EdgeSignOnCellH(Cell0, iEdge0) + 1.0) > tol) { + count++; + } + + int Cell1 = Mesh.CellsOnEdgeH(Edge, 1); + if (Cell1 < DefDecomp->NCellsAll) { + int iEdge1; + for (int i = 0; i < Mesh.NEdgesOnCellH(Cell1); i++) { + if (Mesh.EdgesOnCellH(Cell1, i) == Edge) { + iEdge1 = i; + break; + } + } + if (abs(Mesh.EdgeSignOnCellH(Cell1, iEdge1) - 1.0) > tol) { + count++; + } + } + } + + if (count == 0) { + LOG_INFO("HorzMeshTest: edgeSignOnCell test PASS"); + } else { + LOG_INFO("HorzMeshTest: edgeSignOnCell test FAIL"); + } + + // Test edgeSignOnVertex + // Check that the sign corresponds with convention + // Tests that the edge sign vlues were calculated correctly + count = 0; + for (int Edge = 0; Edge < LocEdges; Edge++) { + int Vertex0 = Mesh.VerticesOnEdgeH(Edge, 0); + int iEdge0; + for (int i = 0; i < Mesh.VertexDegree; i++) { + if (Mesh.EdgesOnVertexH(Vertex0, i) == Edge) { + iEdge0 = i; + break; + } + } + if (abs(Mesh.EdgeSignOnVertexH(Vertex0, iEdge0) + 1.0) > tol) { + count++; + } + + int Vertex1 = Mesh.VerticesOnEdgeH(Edge, 1); + int iEdge1; + for (int i = 0; i < Mesh.VertexDegree; i++) { + if (Mesh.EdgesOnVertexH(Vertex1, i) == Edge) { + iEdge1 = i; + break; + } + } + if (abs(Mesh.EdgeSignOnVertex(Vertex1, iEdge1) - 1.0) > tol) { + count++; + } + } + + if (count == 0) { + LOG_INFO("HorzMeshTest: edgeSignOnVertex test PASS"); + } else { + LOG_INFO("HorzMeshTest: edgeSignOnVertex test FAIL"); + } + // TODO: Test that device arrays are identical // Finalize Omega objects