From 22e28c51677c5391ba20d0393a88929e37a3f16b Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Wed, 7 Feb 2024 10:29:08 -0700 Subject: [PATCH 1/4] Add operators --- components/omega/src/CMakeLists.txt | 2 + components/omega/src/ocn/Operators.cpp | 24 +++++++ components/omega/src/ocn/Operators.h | 94 ++++++++++++++++++++++++++ 3 files changed, 120 insertions(+) create mode 100644 components/omega/src/ocn/Operators.cpp create mode 100644 components/omega/src/ocn/Operators.h diff --git a/components/omega/src/CMakeLists.txt b/components/omega/src/CMakeLists.txt index b72215858d45..01fce4b3b0de 100644 --- a/components/omega/src/CMakeLists.txt +++ b/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) diff --git a/components/omega/src/ocn/Operators.cpp b/components/omega/src/ocn/Operators.cpp new file mode 100644 index 000000000000..4e8629d02195 --- /dev/null +++ b/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 diff --git a/components/omega/src/ocn/Operators.h b/components/omega/src/ocn/Operators.h new file mode 100644 index 000000000000..166e847e2fbf --- /dev/null +++ b/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 From b142ce58850a77b23d769e62a53dc752997af71b Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Wed, 7 Feb 2024 11:12:50 -0700 Subject: [PATCH 2/4] Add operators test --- components/omega/test/CMakeLists.txt | 39 ++ components/omega/test/ocn/OperatorsTest.cpp | 458 ++++++++++++++++++++ 2 files changed, 497 insertions(+) create mode 100644 components/omega/test/ocn/OperatorsTest.cpp diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 518f76dc2f06..b69b487d527a 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -221,6 +221,44 @@ add_test( COMMAND ${MPI_EXEC} -n 8 -- ./${_TestHorzMeshName} ) +################ +# 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} +) +yakl_process_target(${_TestOperatorsName}) + ############# # IO test ############# @@ -357,6 +395,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" diff --git a/components/omega/test/ocn/OperatorsTest.cpp b/components/omega/test/ocn/OperatorsTest.cpp new file mode 100644 index 000000000000..a00a1c2343ef --- /dev/null +++ b/components/omega/test/ocn/OperatorsTest.cpp @@ -0,0 +1,458 @@ +#include "Operators.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "IO.h" +#include "Logging.h" +#include "MachEnv.h" +#include "mpi.h" + +#include +#include + +using namespace OMEGA; +using yakl::c::parallel_for; + +bool isApprox(Real X, Real Y, Real RTol) { + return std::abs(X - Y) <= RTol * std::max(std::abs(X), std::abs(Y)); +} + +struct ExactFunctions { + Real PI = M_PI; + // TODO: get this from the mesh once we support periodic planar meshes + Real Lx = 1; + Real Ly = std::sqrt(3) / 2; + + YAKL_INLINE Real exactScalar(Real X, Real Y) const { + return std::sin(2 * PI * X / Lx) * std::sin(2 * PI * Y / Ly); + } + + YAKL_INLINE Real exactGradScalarX(Real X, Real Y) const { + return 2 * PI / Lx * std::cos(2 * PI * X / Lx) * + std::sin(2 * PI * Y / Ly); + } + + YAKL_INLINE Real exactGradScalarY(Real X, Real Y) const { + return 2 * PI / Ly * std::sin(2 * PI * X / Lx) * + std::cos(2 * PI * Y / Ly); + } + + YAKL_INLINE Real exactVecX(Real X, Real Y) const { + return std::sin(2 * PI * X / Lx) * std::cos(2 * PI * Y / Ly); + } + + YAKL_INLINE Real exactVecY(Real X, Real Y) const { + return std::cos(2 * PI * X / Lx) * std::sin(2 * PI * Y / Ly); + } + + YAKL_INLINE Real exactDivVec(Real X, Real Y) const { + return 2 * PI * (1. / Lx + 1. / Ly) * std::cos(2 * PI * X / Lx) * + std::cos(2 * PI * Y / Ly); + } + + YAKL_INLINE Real exactCurlVec(Real X, Real Y) const { + return 2 * PI * (-1. / Lx + 1. / Ly) * std::sin(2 * PI * X / Lx) * + std::sin(2 * PI * Y / Ly); + } +}; + +int testDivergence(Real RTol) { + int Err; + ExactFunctions EF; + + const auto &mesh = HorzMesh::getDefault(); + auto XEdge = mesh->XEdgeH.createDeviceCopy(); + auto YEdge = mesh->YEdgeH.createDeviceCopy(); + auto &AngleEdge = mesh->AngleEdge; + + // Prepare operator input + Array1DReal VecEdge("VecEdge", mesh->NEdgesSize); + parallel_for( + mesh->NEdgesOwned, YAKL_LAMBDA(int IEdge) { + const Real X = XEdge(IEdge); + const Real Y = YEdge(IEdge); + + const Real VecX = EF.exactVecX(X, Y); + const Real VecY = EF.exactVecY(X, Y); + + const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); + const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); + + VecEdge(IEdge) = EdgeNormalX * VecX + EdgeNormalY * VecY; + }); + + // Perform halo exchange + Halo MyHalo(MachEnv::getDefaultEnv(), Decomp::getDefault()); + auto VecEdgeH = VecEdge.createHostCopy(); + MyHalo.exchangeFullArrayHalo(VecEdgeH, OnEdge); + VecEdgeH.deep_copy_to(VecEdge); + + auto XCell = mesh->XCellH.createDeviceCopy(); + auto YCell = mesh->YCellH.createDeviceCopy(); + auto &AreaCell = mesh->AreaCell; + + // Compute element-wise errors + Array1DReal LInfCell("LInfCell", mesh->NCellsOwned); + Array1DReal L2Cell("L2Cell", mesh->NCellsOwned); + DivergenceOnCell DivergenceCell(mesh); + parallel_for( + mesh->NCellsOwned, YAKL_LAMBDA(int ICell) { + // Numerical result + const Real DivCellNum = DivergenceCell(ICell, VecEdge); + + // Exact result + const Real X = XCell(ICell); + const Real Y = YCell(ICell); + const Real DivCellExact = EF.exactDivVec(X, Y); + + // Errors + LInfCell(ICell) = std::abs(DivCellNum - DivCellExact); + L2Cell(ICell) = AreaCell(ICell) * LInfCell(ICell) * LInfCell(ICell); + }); + + // Compute global error norms + const Real LInfErrorLoc = yakl::intrinsics::maxval(LInfCell); + const Real L2ErrorLoc = yakl::intrinsics::sum(L2Cell); + + MPI_Comm Comm = MachEnv::getDefaultEnv()->getComm(); + Real LInfError; + Err = + MPI_Allreduce(&LInfErrorLoc, &LInfError, 1, MPI_RealKind, MPI_MAX, Comm); + + Real L2Error; + Err = MPI_Allreduce(&L2ErrorLoc, &L2Error, 1, MPI_RealKind, MPI_SUM, Comm); + L2Error = std::sqrt(L2Error); + + // Check error values + const Real ExpectedLInfError = 0.016907664729341576; + const Real ExpectedL2Error = 0.00786717747637947704; + + if (Err == 0 && isApprox(LInfError, ExpectedLInfError, RTol) && + isApprox(L2Error, ExpectedL2Error, RTol)) { + return 0; + } else { + return 1; + } +} + +int testGradient(Real RTol) { + int Err; + ExactFunctions EF; + + const auto &mesh = HorzMesh::getDefault(); + const auto XCell = mesh->XCellH.createDeviceCopy(); + const auto YCell = mesh->YCellH.createDeviceCopy(); + + // Prepare operator input + Array1DReal ScalarCell("ScalarCell", mesh->NCellsSize); + parallel_for( + mesh->NCellsOwned, YAKL_LAMBDA(int ICell) { + const Real X = XCell(ICell); + const Real Y = YCell(ICell); + ScalarCell(ICell) = EF.exactScalar(X, Y); + }); + + // Perform halo exchange + Halo MyHalo(MachEnv::getDefaultEnv(), Decomp::getDefault()); + auto ScalarCellH = ScalarCell.createHostCopy(); + MyHalo.exchangeFullArrayHalo(ScalarCellH, OnCell); + ScalarCellH.deep_copy_to(ScalarCell); + + const auto XEdge = mesh->XEdgeH.createDeviceCopy(); + const auto YEdge = mesh->YEdgeH.createDeviceCopy(); + const auto &AngleEdge = mesh->AngleEdge; + const auto &DcEdge = mesh->DcEdge; + const auto &DvEdge = mesh->DvEdge; + + // Compute element-wise errors + Array1DReal LInfEdge("LInfEdge", mesh->NEdgesOwned); + Array1DReal L2Edge("L2Edge", mesh->NEdgesOwned); + GradientOnEdge GradientEdge(mesh); + parallel_for( + mesh->NEdgesOwned, YAKL_LAMBDA(int IEdge) { + // Numerical result + const Real GradScalarNum = GradientEdge(IEdge, ScalarCell); + + // Exact result + const Real X = XEdge(IEdge); + const Real Y = YEdge(IEdge); + const Real GradScalarExactX = EF.exactGradScalarX(X, Y); + const Real GradScalarExactY = EF.exactGradScalarY(X, Y); + const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); + const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); + const Real GradScalarExact = + EdgeNormalX * GradScalarExactX + EdgeNormalY * GradScalarExactY; + + // Errors + LInfEdge(IEdge) = std::abs(GradScalarNum - GradScalarExact); + const Real AreaEdge = DcEdge(IEdge) * DvEdge(IEdge) / 2; + L2Edge(IEdge) = AreaEdge * LInfEdge(IEdge) * LInfEdge(IEdge); + }); + + // Compute global error norms + const Real LInfErrorLoc = yakl::intrinsics::maxval(LInfEdge); + const Real L2ErrorLoc = yakl::intrinsics::sum(L2Edge); + + MPI_Comm Comm = MachEnv::getDefaultEnv()->getComm(); + Real LInfError; + Err = + MPI_Allreduce(&LInfErrorLoc, &LInfError, 1, MPI_RealKind, MPI_MAX, Comm); + + Real L2Error; + Err = MPI_Allreduce(&L2ErrorLoc, &L2Error, 1, MPI_RealKind, MPI_SUM, Comm); + L2Error = std::sqrt(L2Error); + + // Check error values + const Real ExpectedLInfError = 0.0078388002934621781; + const Real ExpectedL2Error = 0.00424268862440643248; + + if (Err == 0 && isApprox(LInfError, ExpectedLInfError, RTol) && + isApprox(L2Error, ExpectedL2Error, RTol)) { + return 0; + } else { + return 1; + } +} + +int testCurl(Real RTol) { + int Err; + ExactFunctions EF; + + const auto &mesh = HorzMesh::getDefault(); + const auto XEdge = mesh->XEdgeH.createDeviceCopy(); + const auto YEdge = mesh->YEdgeH.createDeviceCopy(); + const auto &AngleEdge = mesh->AngleEdge; + + // Prepare operator input + Array1DReal VecEdge("VecEdge", mesh->NEdgesSize); + parallel_for( + mesh->NEdgesOwned, YAKL_LAMBDA(int IEdge) { + const Real X = XEdge(IEdge); + const Real Y = YEdge(IEdge); + + const Real VecExactX = EF.exactVecX(X, Y); + const Real VecExactY = EF.exactVecY(X, Y); + const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); + const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); + VecEdge(IEdge) = EdgeNormalX * VecExactX + EdgeNormalY * VecExactY; + }); + + // Perform halo exchange + Halo MyHalo(MachEnv::getDefaultEnv(), Decomp::getDefault()); + auto VecEdgeH = VecEdge.createHostCopy(); + MyHalo.exchangeFullArrayHalo(VecEdgeH, OnEdge); + VecEdgeH.deep_copy_to(VecEdge); + + const auto XVertex = mesh->XVertexH.createDeviceCopy(); + const auto YVertex = mesh->YVertexH.createDeviceCopy(); + const auto &AreaTriangle = mesh->AreaTriangle; + + // Compute element-wise errors + Array1DReal LInfVertex("LInfVertex", mesh->NVerticesOwned); + Array1DReal L2Vertex("L2Vertex", mesh->NVerticesOwned); + CurlOnVertex CurlVertex(mesh); + parallel_for( + mesh->NVerticesOwned, YAKL_LAMBDA(int IVertex) { + // Numerical result + const Real CurlNum = CurlVertex(IVertex, VecEdge); + + // Exact result + const Real X = XVertex(IVertex); + const Real Y = YVertex(IVertex); + const Real CurlExact = EF.exactCurlVec(X, Y); + + // Errors + LInfVertex(IVertex) = std::abs(CurlNum - CurlExact); + L2Vertex(IVertex) = + AreaTriangle(IVertex) * LInfVertex(IVertex) * LInfVertex(IVertex); + }); + + // Compute global error norms + const Real LInfErrorLoc = yakl::intrinsics::maxval(LInfVertex); + const Real L2ErrorLoc = yakl::intrinsics::sum(L2Vertex); + + MPI_Comm Comm = MachEnv::getDefaultEnv()->getComm(); + Real LInfError; + Err = + MPI_Allreduce(&LInfErrorLoc, &LInfError, 1, MPI_RealKind, MPI_MAX, Comm); + + Real L2Error; + Err = MPI_Allreduce(&L2ErrorLoc, &L2Error, 1, MPI_RealKind, MPI_SUM, Comm); + L2Error = std::sqrt(L2Error); + + // Check error values + const Real ExpectedLInfError = 0.156364592741396718; + const Real ExpectedL2Error = 0.0729744189366629548; + + if (Err == 0 && isApprox(LInfError, ExpectedLInfError, RTol) && + isApprox(L2Error, ExpectedL2Error, RTol)) { + return 0; + } else { + return 1; + } +} + +int testRecon(Real RTol) { + int Err; + ExactFunctions EF; + + const auto &mesh = HorzMesh::getDefault(); + const auto XEdge = mesh->XEdgeH.createDeviceCopy(); + const auto YEdge = mesh->YEdgeH.createDeviceCopy(); + const auto &AngleEdge = mesh->AngleEdge; + + // Prepare operator input + Array1DReal VecEdge("VecEdge", mesh->NEdgesSize); + parallel_for( + mesh->NEdgesOwned, YAKL_LAMBDA(int IEdge) { + const Real X = XEdge(IEdge); + const Real Y = YEdge(IEdge); + + const Real VecExactX = EF.exactVecX(X, Y); + const Real VecExactY = EF.exactVecY(X, Y); + const Real EdgeNormalX = std::cos(AngleEdge(IEdge)); + const Real EdgeNormalY = std::sin(AngleEdge(IEdge)); + VecEdge(IEdge) = EdgeNormalX * VecExactX + EdgeNormalY * VecExactY; + }); + + // Perform halo exchange + Halo MyHalo(MachEnv::getDefaultEnv(), Decomp::getDefault()); + auto VecEdgeH = VecEdge.createHostCopy(); + MyHalo.exchangeFullArrayHalo(VecEdgeH, OnEdge); + VecEdgeH.deep_copy_to(VecEdge); + + const auto &DcEdge = mesh->DcEdge; + const auto &DvEdge = mesh->DvEdge; + + // Compute element-wise errors + Array1DReal LInfEdge("LInfEdge", mesh->NEdgesOwned); + Array1DReal L2Edge("L2Edge", mesh->NEdgesOwned); + TangentialReconOnEdge TanReconEdge(mesh); + parallel_for( + mesh->NEdgesOwned, YAKL_LAMBDA(int IEdge) { + // Numerical result + const Real VecReconNum = TanReconEdge(IEdge, VecEdge); + + // Exact result + const Real X = XEdge(IEdge); + const Real Y = YEdge(IEdge); + const Real VecX = EF.exactVecX(X, Y); + const Real VecY = EF.exactVecY(X, Y); + const Real EdgeTangentX = -std::sin(AngleEdge(IEdge)); + const Real EdgeTangentY = std::cos(AngleEdge(IEdge)); + const Real VecReconExact = EdgeTangentX * VecX + EdgeTangentY * VecY; + + // Errors + LInfEdge(IEdge) = std::abs(VecReconNum - VecReconExact); + const Real AreaEdge = DcEdge(IEdge) * DvEdge(IEdge) / 2; + L2Edge(IEdge) = AreaEdge * LInfEdge(IEdge) * LInfEdge(IEdge); + }); + + // Compute global error norms + const Real LInfErrorLoc = yakl::intrinsics::maxval(LInfEdge); + const Real L2ErrorLoc = yakl::intrinsics::sum(L2Edge); + + MPI_Comm Comm = MachEnv::getDefaultEnv()->getComm(); + Real LInfError; + Err = + MPI_Allreduce(&LInfErrorLoc, &LInfError, 1, MPI_RealKind, MPI_MAX, Comm); + + Real L2Error; + Err = MPI_Allreduce(&L2ErrorLoc, &L2Error, 1, MPI_RealKind, MPI_SUM, Comm); + L2Error = std::sqrt(L2Error); + + // Check error values + const Real ExpectedLInfError = 0.00449932090822358077; + const Real ExpectedL2Error = 0.00194202022746063673; + + if (Err == 0 && isApprox(LInfError, ExpectedLInfError, RTol) && + isApprox(L2Error, ExpectedL2Error, RTol)) { + return 0; + } else { + return 1; + } +} + +//------------------------------------------------------------------------------ +// The initialization routine for Operators testing +int initOperatorsTest(int argc, char *argv[]) { + + MPI_Init(&argc, &argv); + yakl::init(); + + int Err = 0; + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefaultEnv(); + MPI_Comm DefComm = DefEnv->getComm(); + + Err = IO::init(DefComm); + if (Err != 0) + LOG_ERROR("HorzMeshTest: error initializing parallel IO"); + + Err = Decomp::init(); + if (Err != 0) + LOG_ERROR("HorzMeshTest: error initializing default decomposition"); + + Err = HorzMesh::init(); + if (Err != 0) + LOG_ERROR("HorzMeshTest: error initializing default mesh"); + + return Err; +} + +void finalizeOperatorsTest() { + HorzMesh::clear(); + Decomp::clear(); + MachEnv::removeAll(); + yakl::finalize(); + MPI_Finalize(); +} + +int main(int argc, char *argv[]) { + int Err = initOperatorsTest(argc, argv); + if (Err != 0) + LOG_CRITICAL("OperatorsTest: Error initializing"); + + const Real RTol = sizeof(Real) == 4 ? 1e-2 : 1e-10; + + int DivErr = testDivergence(RTol); + if (DivErr == 0) { + LOG_INFO("OperatorsTest: Divergence PASS"); + } else { + Err = DivErr; + LOG_INFO("OperatorsTest: Divergence FAIL"); + } + + int GradErr = testGradient(RTol); + if (GradErr == 0) { + LOG_INFO("OperatorsTest: Gradient PASS"); + } else { + Err = GradErr; + LOG_INFO("OperatorsTest: Gradient FAIL"); + } + + int CurlErr = testCurl(RTol); + if (CurlErr == 0) { + LOG_INFO("OperatorsTest: Curl PASS"); + } else { + Err = CurlErr; + LOG_INFO("OperatorsTest: Curl FAIL"); + } + + int ReconErr = testRecon(RTol); + if (Err == 0) { + LOG_INFO("OperatorsTest: Recon PASS"); + } else { + Err = ReconErr; + LOG_INFO("OperatorsTest: Recon FAIL"); + } + + if (Err == 0) { + LOG_INFO("OperatorsTest: Successful completion"); + } + + finalizeOperatorsTest(); +} // end of main +//===-----------------------------------------------------------------------===/ From 024b03e65395ed0a3e3879b17c9f79e93f02562a Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Wed, 7 Feb 2024 16:22:54 -0700 Subject: [PATCH 3/4] Setup test meshes using CMake --- components/omega/test/CMakeLists.txt | 51 ++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index b69b487d527a..18872a03deaa 100644 --- a/components/omega/test/CMakeLists.txt +++ b/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 @@ -220,6 +268,7 @@ add_test( NAME HORZMESH_TEST COMMAND ${MPI_EXEC} -n 8 -- ./${_TestHorzMeshName} ) +set_tests_properties(HORZMESH_TEST PROPERTIES FIXTURES_REQUIRED QU240Fixture) ################ # Operators test @@ -257,6 +306,7 @@ add_test( NAME OPERATORS_TEST COMMAND ${MPI_EXEC} -n 2 -- ./${_TestOperatorsName} ) +set_tests_properties(OPERATORS_TEST PROPERTIES FIXTURES_REQUIRED PlanarFixture) yakl_process_target(${_TestOperatorsName}) ############# @@ -297,6 +347,7 @@ add_test( NAME IO_TEST COMMAND ${MPI_EXEC} -n 8 -- ./${_TestIOName} ) +set_tests_properties(IO_TEST PROPERTIES FIXTURES_REQUIRED QU240Fixture) ################## # Config test From 8258752e02254ec748bafd5d6512d584f6d2a82a Mon Sep 17 00:00:00 2001 From: Maciej Waruszewski Date: Fri, 9 Feb 2024 16:55:24 -0700 Subject: [PATCH 4/4] GPU fixes in HorzMesh --- components/omega/src/ocn/HorzMesh.cpp | 8 ++++++++ components/omega/test/CMakeLists.txt | 7 +++++++ components/omega/test/ocn/HorzMeshTest.cpp | 2 +- 3 files changed, 16 insertions(+), 1 deletion(-) diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 1f6a354725c0..39946122659b 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/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++) { diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 18872a03deaa..c24da46e26ee 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -257,12 +257,19 @@ 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 diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index f603f48fa513..b69d6054758e 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/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++; } }