Add mesh operators

components/omega/src/CMakeLists.txt

@@ -43,6 +43,8 @@ if(GKlib_FOUND)
     target_link_libraries(${OMEGA_LIB_NAME} gklib)
 endif()
 
+yakl_process_target(${OMEGA_LIB_NAME})
+
 # build Omega executable
 if(OMEGA_BUILD_EXECUTABLE)
 

components/omega/src/ocn/HorzMesh.cpp

@@ -512,6 +512,10 @@ void HorzMesh::readCoriolis() {
 void HorzMesh::computeEdgeSign() {
 
    EdgeSignOnCell = Array2DR8("EdgeSignOnCell", NCellsSize, MaxEdges);
+   YAKL_SCOPE(EdgeSignOnCell, this->EdgeSignOnCell);
+   YAKL_SCOPE(NEdgesOnCell, this->NEdgesOnCell);
+   YAKL_SCOPE(EdgesOnCell, this->EdgesOnCell);
+   YAKL_SCOPE(CellsOnEdge, this->CellsOnEdge);
    yakl::c::parallel_for(
        yakl::c::SimpleBounds<1>(NCellsAll), YAKL_LAMBDA(int Cell) {
           for (int i = 0; i < NEdgesOnCell(Cell); i++) {
@@ -529,6 +533,10 @@ void HorzMesh::computeEdgeSign() {
 
    EdgeSignOnVertex =
        Array2DR8("EdgeSignOnVertex", NVerticesSize, VertexDegree);
+   YAKL_SCOPE(EdgeSignOnVertex, this->EdgeSignOnVertex);
+   YAKL_SCOPE(EdgesOnVertex, this->EdgesOnVertex);
+   YAKL_SCOPE(VerticesOnEdge, this->VerticesOnEdge);
+   YAKL_SCOPE(VertexDegree, this->VertexDegree);
    yakl::c::parallel_for(
        yakl::c::SimpleBounds<1>(NVerticesAll), YAKL_LAMBDA(int Vertex) {
           for (int i = 0; i < VertexDegree; i++) {

components/omega/src/ocn/Operators.cpp

@@ -0,0 +1,24 @@
+#include "Operators.h"
+#include "DataTypes.h"
+#include "HorzMesh.h"
+
+namespace OMEGA {
+
+DivergenceOnCell::DivergenceOnCell(HorzMesh const *mesh)
+    : NEdgesOnCell(mesh->NEdgesOnCell), EdgesOnCell(mesh->EdgesOnCell),
+      DvEdge(mesh->DvEdge), AreaCell(mesh->AreaCell),
+      EdgeSignOnCell(mesh->EdgeSignOnCell) {}
+
+GradientOnEdge::GradientOnEdge(HorzMesh const *mesh)
+    : CellsOnEdge(mesh->CellsOnEdge), DcEdge(mesh->DcEdge) {}
+
+CurlOnVertex::CurlOnVertex(HorzMesh const *mesh)
+    : VertexDegree(mesh->VertexDegree), EdgesOnVertex(mesh->EdgesOnVertex),
+      DcEdge(mesh->DcEdge), AreaTriangle(mesh->AreaTriangle),
+      EdgeSignOnVertex(mesh->EdgeSignOnVertex) {}
+
+TangentialReconOnEdge::TangentialReconOnEdge(HorzMesh const *mesh)
+    : NEdgesOnEdge(mesh->NEdgesOnEdge), EdgesOnEdge(mesh->EdgesOnEdge),
+      WeightsOnEdge(mesh->WeightsOnEdge) {}
+
+} // namespace OMEGA

components/omega/src/ocn/Operators.h

@@ -0,0 +1,94 @@
+#ifndef OMEGA_OPERATORS_H
+#define OMEGA_OPERATORS_H
+
+#include "DataTypes.h"
+#include "HorzMesh.h"
+
+namespace OMEGA {
+
+class DivergenceOnCell {
+ public:
+   DivergenceOnCell(HorzMesh const *mesh);
+
+   YAKL_INLINE Real operator()(int ICell, const Array1DReal &VecEdge) const {
+      Real DivCell = 0;
+      for (int J = 0; J < NEdgesOnCell(ICell); ++J) {
+         const int JEdge = EdgesOnCell(ICell, J);
+         DivCell -= DvEdge(JEdge) * EdgeSignOnCell(ICell, J) * VecEdge(JEdge);
+      }
+      const Real InvAreaCell = 1. / AreaCell(ICell);
+      DivCell *= InvAreaCell;
+      return DivCell;
+   }
+
+ private:
+   Array1DI4 NEdgesOnCell;
+   Array2DI4 EdgesOnCell;
+   Array1DR8 DvEdge;
+   Array1DR8 AreaCell;
+   Array2DR8 EdgeSignOnCell;
+};
+
+class GradientOnEdge {
+ public:
+   GradientOnEdge(HorzMesh const *mesh);
+
+   YAKL_INLINE Real operator()(int IEdge, const Array1DReal &ScalarCell) const {
+      const auto JCell0    = CellsOnEdge(IEdge, 0);
+      const auto JCell1    = CellsOnEdge(IEdge, 1);
+      const Real InvDcEdge = 1. / DcEdge(IEdge);
+      const Real GradEdge =
+          InvDcEdge * (ScalarCell(JCell1) - ScalarCell(JCell0));
+      return GradEdge;
+   }
+
+ private:
+   Array2DI4 CellsOnEdge;
+   Array1DR8 DcEdge;
+};
+
+class CurlOnVertex {
+ public:
+   CurlOnVertex(HorzMesh const *mesh);
+
+   YAKL_INLINE Real operator()(int IVertex, const Array1DReal &VecEdge) const {
+      Real CurlVertex = 0;
+      for (int J = 0; J < VertexDegree; ++J) {
+         const int JEdge = EdgesOnVertex(IVertex, J);
+         CurlVertex +=
+             DcEdge(JEdge) * EdgeSignOnVertex(IVertex, J) * VecEdge(JEdge);
+      }
+      const Real InvAreaTriangle = 1. / AreaTriangle(IVertex);
+      CurlVertex *= InvAreaTriangle;
+      return CurlVertex;
+   }
+
+ private:
+   I4 VertexDegree;
+   Array2DI4 EdgesOnVertex;
+   Array1DR8 DcEdge;
+   Array1DR8 AreaTriangle;
+   Array2DR8 EdgeSignOnVertex;
+};
+
+class TangentialReconOnEdge {
+ public:
+   TangentialReconOnEdge(HorzMesh const *mesh);
+
+   YAKL_INLINE Real operator()(int IEdge, const Array1DReal &VecEdge) const {
+      Real ReconEdge = 0;
+      for (int J = 0; J < NEdgesOnEdge(IEdge); ++J) {
+         const int JEdge = EdgesOnEdge(IEdge, J);
+         ReconEdge += WeightsOnEdge(IEdge, J) * VecEdge(JEdge);
+      }
+      return ReconEdge;
+   }
+
+ private:
+   Array1DI4 NEdgesOnEdge;
+   Array2DI4 EdgesOnEdge;
+   Array2DR8 WeightsOnEdge;
+};
+
+} // namespace OMEGA
+#endif

components/omega/test/CMakeLists.txt

@@ -1,5 +1,51 @@
 # Omega Unit Tests
 
+##################
+# Test meshes
+##################
+
+include(FetchContent)
+set(FETCHCONTENT_UPDATES_DISCONNECTED ON)
+
+FetchContent_Declare(
+  QU240Mesh
+  URL      https://web.lcrc.anl.gov/public/e3sm/inputdata/ocn/mpas-o/oQU240/ocean.QU.240km.151209.nc
+  URL_HASH SHA256=a3758f88ceff3d91e86dba7922f6dd7d5672157b4793ef78214624ab8b2724ae
+  DOWNLOAD_DIR ${CMAKE_CURRENT_BINARY_DIR}/meshes
+  DOWNLOAD_NO_EXTRACT TRUE
+)
+
+FetchContent_Declare(
+  PlanarMesh
+  URL      https://gist.github.com/mwarusz/f8caf260398dbe140d2102ec46a41268/raw/e3c29afbadc835797604369114321d93fd69886d/PlanarPeriodic48x48.nc
+  URL_HASH SHA256=fbe73f796fcad0fcd0ea2df6f9deb40108755464c5fa3f802719372548833b44
+  DOWNLOAD_DIR ${CMAKE_CURRENT_BINARY_DIR}/meshes
+  DOWNLOAD_NO_EXTRACT TRUE
+)
+
+FetchContent_MakeAvailable(QU240Mesh PlanarMesh)
+
+add_test(NAME SetupQU240Mesh
+  COMMAND ${CMAKE_COMMAND} -E create_symlink
+    ${CMAKE_CURRENT_BINARY_DIR}/meshes/ocean.QU.240km.151209.nc
+    ${CMAKE_CURRENT_BINARY_DIR}/OmegaMesh.nc)
+set_tests_properties(SetupQU240Mesh PROPERTIES FIXTURES_SETUP QU240Fixture)
+
+add_test(NAME CleanupQU240Mesh
+  COMMAND ${CMAKE_COMMAND} -E rm ${CMAKE_CURRENT_BINARY_DIR}/OmegaMesh.nc)
+set_tests_properties(CleanupQU240Mesh PROPERTIES FIXTURES_CLEANUP QU240Fixture)
+
+add_test(NAME SetupPlanarMesh
+  COMMAND ${CMAKE_COMMAND} -E create_symlink
+    ${CMAKE_CURRENT_BINARY_DIR}/meshes/PlanarPeriodic48x48.nc
+    ${CMAKE_CURRENT_BINARY_DIR}/OmegaMesh.nc)
+set_tests_properties(SetupPlanarMesh PROPERTIES FIXTURES_SETUP PlanarFixture)
+set_tests_properties(SetupPlanarMesh PROPERTIES DEPENDS CleanupQU240Mesh)
+
+add_test(NAME CleanupPlanarMesh
+  COMMAND ${CMAKE_COMMAND} -E rm ${CMAKE_CURRENT_BINARY_DIR}/OmegaMesh.nc)
+set_tests_properties(CleanupPlanarMesh PROPERTIES FIXTURES_CLEANUP PlanarFixture)
+
 ##################
 # Data type test
 ##################
@@ -159,6 +205,7 @@ add_test(
   NAME DECOMP_TEST
   COMMAND ${MPI_EXEC} -n 8 -- ./${_TestDecompName}
 )
+set_tests_properties(DECOMP_TEST PROPERTIES FIXTURES_REQUIRED QU240Fixture)
 
 ##################
 # Halo test
@@ -189,6 +236,7 @@ add_test(
   NAME HALO_TEST
   COMMAND ${MPI_EXEC} -n 8 -- ./${_TestHaloName}
 )
+set_tests_properties(HALO_TEST PROPERTIES FIXTURES_REQUIRED QU240Fixture)
 
 ################
 # HorzMesh test
@@ -209,17 +257,64 @@ target_include_directories(
 )
 
 target_compile_options(
+  ${_TestHorzMeshName}
+  PRIVATE
+  ${OMEGA_CXX_FLAGS}
+)
+
+target_link_options(
   ${_TestHorzMeshName}
   PRIVATE
   ${OMEGA_LINK_OPTIONS}
 )
 
 target_link_libraries(${_TestHorzMeshName} ${OMEGA_LIB_NAME} spdlog yakl parmetis metis pioc)
+yakl_process_target(${_TestHorzMeshName})
 
 add_test(
   NAME HORZMESH_TEST
   COMMAND ${MPI_EXEC} -n 8 -- ./${_TestHorzMeshName}
 )
+set_tests_properties(HORZMESH_TEST PROPERTIES FIXTURES_REQUIRED QU240Fixture)
+
+################
+# Operators test
+################
+
+set(_TestOperatorsName testOperators.exe)
+
+# Add broadcast test
+add_executable(${_TestOperatorsName} ocn/OperatorsTest.cpp)
+
+target_include_directories(
+  ${_TestOperatorsName}
+  PRIVATE
+  ${OMEGA_SOURCE_DIR}/src/base
+  ${OMEGA_SOURCE_DIR}/src/infra
+  ${OMEGA_SOURCE_DIR}/src/ocn
+  ${Parmetis_INCLUDE_DIRS}
+)
+
+target_compile_options(
+  ${_TestOperatorsName}
+  PRIVATE
+  ${OMEGA_CXX_FLAGS}
+)
+
+target_link_options(
+  ${_TestOperatorsName}
+  PRIVATE
+  ${OMEGA_LINK_OPTIONS}
+)
+
+target_link_libraries(${_TestOperatorsName} ${OMEGA_LIB_NAME} spdlog yakl parmetis metis pioc)
+
+add_test(
+  NAME OPERATORS_TEST
+  COMMAND ${MPI_EXEC} -n 2 -- ./${_TestOperatorsName}
+)
+set_tests_properties(OPERATORS_TEST PROPERTIES FIXTURES_REQUIRED PlanarFixture)
+yakl_process_target(${_TestOperatorsName})
 
 #############
 # IO test
@@ -259,6 +354,7 @@ add_test(
   NAME IO_TEST
   COMMAND ${MPI_EXEC} -n 8 -- ./${_TestIOName}
 )
+set_tests_properties(IO_TEST PROPERTIES FIXTURES_REQUIRED QU240Fixture)
 
 ##################
 # Config test
@@ -357,6 +453,7 @@ set_tests_properties(
   DECOMP_TEST
   HALO_TEST
   HORZMESH_TEST
+  OPERATORS_TEST
   IO_TEST
   METADATA_TEST
   PROPERTIES FAIL_REGULAR_EXPRESSION "FAIL" PASS_REGULAR_EXPRESSION "PASS"

components/omega/test/ocn/HorzMeshTest.cpp

@@ -629,7 +629,7 @@ int main(int argc, char *argv[]) {
             break;
          }
       }
-      if (abs(Mesh->EdgeSignOnVertex(Vertex1, iEdge1) - 1.0) > tol) {
+      if (abs(Mesh->EdgeSignOnVertexH(Vertex1, iEdge1) - 1.0) > tol) {
          count++;
       }
    }