diff --git a/components/omega/doc/devGuide/DataTypes.md b/components/omega/doc/devGuide/DataTypes.md index 79e8bd157e28..ee95c69269ec 100644 --- a/components/omega/doc/devGuide/DataTypes.md +++ b/components/omega/doc/devGuide/DataTypes.md @@ -9,7 +9,7 @@ integer and floating point variables. Note that these exist in the Omega namespace so use the scoped form OMEGA::I4, etc. In addition, we define a Real data type that is, by default, double precision (8 bytes/64-bit) but if the code is built with a -`-DSINGLE_PRECISION` (see insert link to build system) preprocessor flag, +`-DOMEGA_SINGLE_PRECISION` (see insert link to build system) preprocessor flag, the default Real becomes single precision (4-byte/32-bit). For floating point variables, developers should use this Real type instead of the specific R4 or R8 forms unless the specific form is required diff --git a/components/omega/doc/userGuide/DataTypes.md b/components/omega/doc/userGuide/DataTypes.md index 03fc8bc586fe..2b6fdb28ecdc 100644 --- a/components/omega/doc/userGuide/DataTypes.md +++ b/components/omega/doc/userGuide/DataTypes.md @@ -7,7 +7,7 @@ types to guarantee a specific level of precision. There is only one user-configurable option for precision. When a specific floating point precision is not required, we use a Real data type that is, by default, double precision (8 bytes/64-bit) but if the code is built with a -`-DSINGLE_PRECISION` (see insert link to build system) preprocessor flag, +`-DOMEGA_SINGLE_PRECISION` (see insert link to build system) preprocessor flag, the default Real becomes single precision (4-byte/32-bit). Users are encouraged to use the default double precision unless exploring the performance or accuracy characteristics of single precision. diff --git a/components/omega/src/CMakeLists.txt b/components/omega/src/CMakeLists.txt index 71a3a306b5cd..9e0d0575e6ff 100644 --- a/components/omega/src/CMakeLists.txt +++ b/components/omega/src/CMakeLists.txt @@ -21,6 +21,14 @@ target_compile_definitions( OMEGA_ARCH=${OMEGA_ARCH} ) +if(OMEGA_SINGLE_PRECISION) + target_compile_definitions( + OmegaLibFlags + INTERFACE + OMEGA_SINGLE_PRECISION=1 + ) +endif() + target_link_options( OmegaLibFlags INTERFACE @@ -57,6 +65,26 @@ target_link_libraries( OmegaLibFlags ) +# some tests require single precision version of the library +# if OmegaLib is configured to use double precision we need to +# build another version +if(OMEGA_BUILD_TEST AND NOT OMEGA_SINGLE_PRECISION) + add_library(${OMEGA_LIB_NAME}_SP ${_LIBSRC_FILES}) + + target_link_libraries( + ${OMEGA_LIB_NAME}_SP + PRIVATE + OmegaLibFlags + ) + + target_compile_definitions( + ${OMEGA_LIB_NAME}_SP + PUBLIC + OMEGA_SINGLE_PRECISION=1 + ) +endif() + + # build Omega executable if(OMEGA_BUILD_EXECUTABLE) diff --git a/components/omega/src/base/DataTypes.h b/components/omega/src/base/DataTypes.h index 4e571c936048..11e35aa404e8 100644 --- a/components/omega/src/base/DataTypes.h +++ b/components/omega/src/base/DataTypes.h @@ -27,7 +27,7 @@ using R4 = float; ///< alias for 32-bit (single prec) real using R8 = double; ///< alias for 64-bit (double prec) real /// generic real 64-bit (default) or 32-bit (if -DSINGLE_PRECISION used) -#ifdef SINGLE_PRECISION +#ifdef OMEGA_SINGLE_PRECISION using Real = float; #else using Real = double; diff --git a/components/omega/src/base/Halo.cpp b/components/omega/src/base/Halo.cpp index a5f87071b74d..12bdd3499bff 100644 --- a/components/omega/src/base/Halo.cpp +++ b/components/omega/src/base/Halo.cpp @@ -647,7 +647,7 @@ int Halo::startReceives() { I4 BufferSize = TotSize * MyNeighbor->RecvLists[MyElem].NTot; MyNeighbor->RecvBuffer.resize(BufferSize); IErr[INghbr] = MPI_Irecv(&MyNeighbor->RecvBuffer[0], BufferSize, - MPI_RealKind, MyNeighbor->TaskID, MPI_ANY_TAG, + MPI_DOUBLE, MyNeighbor->TaskID, MPI_ANY_TAG, MyComm, &MyNeighbor->RReq); if (IErr[INghbr] != 0) { LOG_ERROR("MPI error {} on task {} receive from task {}", @@ -676,7 +676,7 @@ int Halo::startSends() { MyNeighbor = &Neighbors[INghbr]; I4 BufferSize = TotSize * MyNeighbor->SendLists[MyElem].NTot; IErr[INghbr] = - MPI_Isend(&MyNeighbor->SendBuffer[0], BufferSize, MPI_RealKind, + MPI_Isend(&MyNeighbor->SendBuffer[0], BufferSize, MPI_DOUBLE, MyNeighbor->TaskID, 0, MyComm, &MyNeighbor->SReq); if (IErr[INghbr] != 0) { LOG_ERROR("MPI error {} on task {} send to task {}", IErr[INghbr], @@ -709,7 +709,7 @@ int Halo::packBuffer(const HostArray1DI4 Array) { for (int IExch = 0; IExch < MyList->NList[ILayer]; ++IExch) { I4 IBuff = MyList->Offsets[ILayer] + IExch; MyNeighbor->SendBuffer[IBuff] = - reinterpret_cast(Array(MyList->Ind[ILayer][IExch])); + reinterpret_cast(Array(MyList->Ind[ILayer][IExch])); } } @@ -726,7 +726,7 @@ int Halo::packBuffer(const HostArray1DI8 Array) { for (int IExch = 0; IExch < MyList->NList[ILayer]; ++IExch) { I4 IBuff = MyList->Offsets[ILayer] + IExch; MyNeighbor->SendBuffer[IBuff] = - reinterpret_cast(Array(MyList->Ind[ILayer][IExch])); + reinterpret_cast(Array(MyList->Ind[ILayer][IExch])); } } @@ -777,7 +777,7 @@ int Halo::packBuffer(const HostArray2DI4 Array) { for (int J = 0; J < NJ; ++J) { I4 IBuff = (MyList->Offsets[ILayer] + IExch) * NJ + J; MyNeighbor->SendBuffer[IBuff] = - reinterpret_cast(Array(MyList->Ind[ILayer][IExch], J)); + reinterpret_cast(Array(MyList->Ind[ILayer][IExch], J)); } } } @@ -797,7 +797,7 @@ int Halo::packBuffer(const HostArray2DI8 Array) { for (int J = 0; J < NJ; ++J) { I4 IBuff = (MyList->Offsets[ILayer] + IExch) * NJ + J; MyNeighbor->SendBuffer[IBuff] = - reinterpret_cast(Array(MyList->Ind[ILayer][IExch], J)); + reinterpret_cast(Array(MyList->Ind[ILayer][IExch], J)); } } } @@ -860,7 +860,7 @@ int Halo::packBuffer(const HostArray3DI4 Array) { I4 IBuff = (K * MyList->NTot + MyList->Offsets[ILayer] + IExch) * NJ + J; - MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( + MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( Array(K, MyList->Ind[ILayer][IExch], J)); } } @@ -885,7 +885,7 @@ int Halo::packBuffer(const HostArray3DI8 Array) { I4 IBuff = (K * MyList->NTot + MyList->Offsets[ILayer] + IExch) * NJ + J; - MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( + MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( Array(K, MyList->Ind[ILayer][IExch], J)); } } @@ -963,7 +963,7 @@ int Halo::packBuffer(const HostArray4DI4 Array) { MyList->Offsets[ILayer] + IExch) * NJ + J; - MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( + MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( Array(L, K, MyList->Ind[ILayer][IExch], J)); } } @@ -992,7 +992,7 @@ int Halo::packBuffer(const HostArray4DI8 Array) { MyList->Offsets[ILayer] + IExch) * NJ + J; - MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( + MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( Array(L, K, MyList->Ind[ILayer][IExch], J)); } } @@ -1081,7 +1081,7 @@ int Halo::packBuffer(const HostArray5DI4 Array) { MyList->Offsets[ILayer] + IExch) * NJ + J; - MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( + MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( Array(M, L, K, MyList->Ind[ILayer][IExch], J)); } } @@ -1113,7 +1113,7 @@ int Halo::packBuffer(const HostArray5DI8 Array) { MyList->Offsets[ILayer] + IExch) * NJ + J; - MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( + MyNeighbor->SendBuffer[IBuff] = reinterpret_cast( Array(M, L, K, MyList->Ind[ILayer][IExch], J)); } } diff --git a/components/omega/src/base/Halo.h b/components/omega/src/base/Halo.h index 4c4920c8ee56..93c4ea47e696 100644 --- a/components/omega/src/base/Halo.h +++ b/components/omega/src/base/Halo.h @@ -30,7 +30,7 @@ namespace OMEGA { // Set the default MPI real data type as single or double precision based on // the default real data type used by OMEGA. -#ifdef SINGLE_PRECISION +#ifdef OMEGA_SINGLE_PRECISION static const MPI_Datatype MPI_RealKind = MPI_FLOAT; #else static const MPI_Datatype MPI_RealKind = MPI_DOUBLE; @@ -132,7 +132,7 @@ class Halo { /// index space. 0 = OnCell, 1 = OnEdge, 2 = OnVertex ExchList SendLists[3], RecvLists[3]; /// Buffers for MPI communication - std::vector SendBuffer, RecvBuffer; + std::vector SendBuffer, RecvBuffer; /// MPI request handles for non-blocking MPI communication MPI_Request RReq, SReq; diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 4f2c6b7c5297..a186a1f88ebc 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -404,135 +404,109 @@ void HorzMesh::finalizeParallelIO() { } // end finalizeParallelIO -//------------------------------------------------------------------------------ -// Read x/y/z and lon/lat coordinates for cells, edges, and vertices -void HorzMesh::readCoordinates() { +// Read 1D vertex array +void HorzMesh::readVertexArray(HostArray1DReal &VertexArrayH, + const std::string &MPASName) { + int Err; - I4 Err; + std::string OmegaName; + std::transform(MPASName.begin(), MPASName.end(), OmegaName.begin(), + [](unsigned char c) { return std::toupper(c); }); - // Read mesh cell coordinates - int XCellID; - XCellH = HostArray1DR8("XCell", NCellsSize); - Err = IO::readArray(XCellH.data(), NCellsAll, "xCell", MeshFileID, - CellDecompR8, XCellID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading xCell"); + // Temporary double precision array for reading + HostArray1DR8 TmpArrayR8(OmegaName + "Tmp", NVerticesSize); + int ArrayID; + Err = IO::readArray(TmpArrayR8.data(), NVerticesAll, MPASName, MeshFileID, + VertexDecompR8, ArrayID); - int YCellID; - YCellH = HostArray1DR8("YCell", NCellsSize); - Err = IO::readArray(YCellH.data(), NCellsAll, "yCell", MeshFileID, - CellDecompR8, YCellID); if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading yCell"); + LOG_CRITICAL("HorzMesh: error reading {}", MPASName); - int ZCellID; - ZCellH = HostArray1DR8("ZCell", NCellsSize); - Err = IO::readArray(ZCellH.data(), NCellsAll, "zCell", MeshFileID, - CellDecompR8, ZCellID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading zCell"); + // Create host array of desired precision and copy the read data into it + VertexArrayH = HostArray1DReal(OmegaName + "H", NVerticesSize); + deepCopy(VertexArrayH, TmpArrayR8); +} - int LonCellID; - LonCellH = HostArray1DR8("LonCell", NCellsSize); - Err = IO::readArray(LonCellH.data(), NCellsAll, "lonCell", MeshFileID, - CellDecompR8, LonCellID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading lonCell"); +// Read 1D edge array +void HorzMesh::readEdgeArray(HostArray1DReal &EdgeArrayH, + const std::string &MPASName) { + int Err; - int LatCellID; - LatCellH = HostArray1DR8("LatCell", NCellsSize); - Err = IO::readArray(LatCellH.data(), NCellsAll, "latCell", MeshFileID, - CellDecompR8, LatCellID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading latCell"); + std::string OmegaName; + std::transform(MPASName.begin(), MPASName.end(), OmegaName.begin(), + [](unsigned char c) { return std::toupper(c); }); - // Read mesh edge coordinateID - int XEdgeID; - XEdgeH = HostArray1DR8("XEdge", NEdgesSize); - Err = IO::readArray(XEdgeH.data(), NEdgesAll, "xEdge", MeshFileID, - EdgeDecompR8, XEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading xEdge"); + // Temporary double precision array for reading + HostArray1DR8 TmpArrayR8(OmegaName + "Tmp", NEdgesSize); + int ArrayID; + Err = IO::readArray(TmpArrayR8.data(), NEdgesAll, MPASName, MeshFileID, + EdgeDecompR8, ArrayID); - int YEdgeID; - YEdgeH = HostArray1DR8("YEdge", NEdgesSize); - Err = IO::readArray(YEdgeH.data(), NEdgesAll, "yEdge", MeshFileID, - EdgeDecompR8, YEdgeID); if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading yEdge"); + LOG_CRITICAL("HorzMesh: error reading {}", MPASName); - int ZEdgeID; - ZEdgeH = HostArray1DR8("ZEdge", NEdgesSize); - Err = IO::readArray(ZEdgeH.data(), NEdgesAll, "zEdge", MeshFileID, - EdgeDecompR8, ZEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading zEdge"); + // Create host array of desired precision and copy the read data into it + EdgeArrayH = HostArray1DReal(OmegaName + "H", NEdgesSize); + deepCopy(EdgeArrayH, TmpArrayR8); +} - int LonEdgeID; - LonEdgeH = HostArray1DR8("LonEdge", NEdgesSize); - Err = IO::readArray(LonEdgeH.data(), NEdgesAll, "lonEdge", MeshFileID, - EdgeDecompR8, LonEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading lonEdge"); +// Read 1D cell array +void HorzMesh::readCellArray(HostArray1DReal &CellArrayH, + const std::string &MPASName) { + int Err; - int LatEdgeID; - LatEdgeH = HostArray1DR8("LatEdge", NEdgesSize); - Err = IO::readArray(LatEdgeH.data(), NEdgesAll, "latEdge", MeshFileID, - EdgeDecompR8, LatEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading latEdge"); + std::string OmegaName; + std::transform(MPASName.begin(), MPASName.end(), OmegaName.begin(), + [](unsigned char c) { return std::toupper(c); }); - // Read mesh vertex coordinates - int XVertexID; - XVertexH = HostArray1DR8("XVertex", NVerticesSize); - Err = IO::readArray(XVertexH.data(), NVerticesAll, "xVertex", MeshFileID, - VertexDecompR8, XVertexID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading xVertex"); + // Temporary double precision array for reading + HostArray1DR8 TmpArrayR8(OmegaName + "Tmp", NCellsSize); + int ArrayID; + Err = IO::readArray(TmpArrayR8.data(), NCellsAll, MPASName, MeshFileID, + CellDecompR8, ArrayID); - int YVertexID; - YVertexH = HostArray1DR8("YVertex", NVerticesSize); - Err = IO::readArray(YVertexH.data(), NVerticesAll, "yVertex", MeshFileID, - VertexDecompR8, YVertexID); if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading yVertex"); + LOG_CRITICAL("HorzMesh: error reading {}", MPASName); - int ZVertexID; - ZVertexH = HostArray1DR8("ZVertex", NVerticesSize); - Err = IO::readArray(ZVertexH.data(), NVerticesAll, "zVertex", MeshFileID, - VertexDecompR8, ZVertexID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading zVertex"); + // Create host array of desired precision and copy the read data into it + CellArrayH = HostArray1DReal(OmegaName + "H", NCellsSize); + deepCopy(CellArrayH, TmpArrayR8); +} - int LonVertexID; - LonVertexH = HostArray1DR8("LonVertex", NVerticesSize); - Err = IO::readArray(LonVertexH.data(), NVerticesAll, "lonVertex", MeshFileID, - VertexDecompR8, LonVertexID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading lonVertex"); +//------------------------------------------------------------------------------ +// Read x/y/z and lon/lat coordinates for cells, edges, and vertices +void HorzMesh::readCoordinates() { - int LatVertexID; - LatVertexH = HostArray1DR8("LatVertex", NVerticesSize); - Err = IO::readArray(LatVertexH.data(), NVerticesAll, "latVertex", MeshFileID, - VertexDecompR8, LatVertexID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading latVertex"); + // Read mesh cell coordinates + readCellArray(XCellH, "xCell"); + readCellArray(YCellH, "yCell"); + readCellArray(ZCellH, "zCell"); + + readCellArray(LonCellH, "lonCell"); + readCellArray(LatCellH, "latCell"); + + // Read mesh edge coordinate + readEdgeArray(XEdgeH, "xEdge"); + readEdgeArray(YEdgeH, "yEdge"); + readEdgeArray(ZEdgeH, "zEdge"); + + readEdgeArray(LonEdgeH, "lonEdge"); + readEdgeArray(LatEdgeH, "latEdge"); + + // Read mesh vertex coordinates + readVertexArray(XVertexH, "xVertex"); + readVertexArray(YVertexH, "yVertex"); + readVertexArray(ZVertexH, "zVertex"); + + readVertexArray(LonVertexH, "lonVertex"); + readVertexArray(LatVertexH, "latVertex"); } // end readCoordinates //------------------------------------------------------------------------------ // Read the cell-centered bottom depth void HorzMesh::readBottomDepth() { - - I4 Err; - - int BottomDepthID; - BottomDepthH = HostArray1DR8("BottomDepth", NCellsSize); - Err = IO::readArray(BottomDepthH.data(), NCellsAll, "bottomDepth", - MeshFileID, CellDecompR8, BottomDepthID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading bottomDepth"); - + readCellArray(BottomDepthH, "bottomDepth"); } // end readDepth //------------------------------------------------------------------------------ @@ -540,59 +514,37 @@ void HorzMesh::readBottomDepth() { // lengths (between centers and vertices), and edge angles void HorzMesh::readMeasurements() { - I4 Err; + readCellArray(AreaCellH, "areaCell"); - int AreaCellID; - AreaCellH = HostArray1DR8("AreaCell", NCellsSize); - Err = IO::readArray(AreaCellH.data(), NCellsAll, "areaCell", MeshFileID, - CellDecompR8, AreaCellID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading areaCell"); + readVertexArray(AreaTriangleH, "areaTriangle"); - int AreaTriangleID; - AreaTriangleH = HostArray1DR8("AreaTriangle", NVerticesSize); - Err = IO::readArray(AreaTriangleH.data(), NVerticesAll, "areaTriangle", - MeshFileID, VertexDecompR8, AreaTriangleID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading areaTriangle"); + readEdgeArray(DvEdgeH, "dvEdge"); - int DvEdgeID; - DvEdgeH = HostArray1DR8("DvEdge", NEdgesSize); - Err = IO::readArray(DvEdgeH.data(), NEdgesAll, "dvEdge", MeshFileID, - EdgeDecompR8, DvEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading dvEdge"); + readEdgeArray(DcEdgeH, "dcEdge"); - int DcEdgeID; - DcEdgeH = HostArray1DR8("DcEdge", NEdgesSize); - Err = IO::readArray(DcEdgeH.data(), NEdgesAll, "dcEdge", MeshFileID, - EdgeDecompR8, DcEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading dcEdge"); + readEdgeArray(AngleEdgeH, "angleEdge"); - int AngleEdgeID; - AngleEdgeH = HostArray1DR8("AngleEdge", NEdgesSize); - Err = IO::readArray(AngleEdgeH.data(), NEdgesAll, "angleEdge", MeshFileID, - EdgeDecompR8, AngleEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading angleEdge"); + readCellArray(MeshDensityH, "meshDensity"); - int MeshDensityID; - MeshDensityH = HostArray1DR8("MeshDensity", NCellsSize); - Err = IO::readArray(MeshDensityH.data(), NCellsAll, "meshDensity", - MeshFileID, CellDecompR8, MeshDensityID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading meshDensity"); + // not using helper function since it kiteAreas is a 2d array + I4 Err; + // Read into a temporary double precision array int KiteAreasOnVertexID; - KiteAreasOnVertexH = - HostArray2DR8("KiteAreasOnVertex", NVerticesSize, VertexDegree); - Err = IO::readArray(KiteAreasOnVertexH.data(), NVerticesAll * VertexDegree, - "kiteAreasOnVertex", MeshFileID, OnVertexDecompR8, - KiteAreasOnVertexID); + + HostArray2DR8 TmpKiteAreasOnVertexR8("KiteAreasOnVertex", NVerticesSize, + VertexDegree); + Err = IO::readArray(TmpKiteAreasOnVertexR8.data(), + NVerticesAll * VertexDegree, "kiteAreasOnVertex", + MeshFileID, OnVertexDecompR8, KiteAreasOnVertexID); if (Err != 0) LOG_CRITICAL("HorzMesh: error reading kiteAreasOnVertex"); + // Create and fill array with Real precision + KiteAreasOnVertexH = + HostArray2DReal("KiteAreasOnVertex", NVerticesSize, VertexDegree); + deepCopy(KiteAreasOnVertexH, TmpKiteAreasOnVertexR8); + } // end readMeasurements //------------------------------------------------------------------------------ @@ -602,13 +554,16 @@ void HorzMesh::readWeights() { I4 Err; int WeightsOnEdgeID; - WeightsOnEdgeH = HostArray2DR8("WeightsOnEdge", NEdgesSize, MaxEdges2); - Err = IO::readArray(WeightsOnEdgeH.data(), NEdgesAll * MaxEdges2, - "weightsOnEdge", MeshFileID, OnEdgeDecompR8, - WeightsOnEdgeID); + HostArray2DR8 TmpWeightsOnEdgeR8("WeightsOnEdge", NEdgesSize, MaxEdges2); + Err = IO::readArray(TmpWeightsOnEdgeR8.data(), NEdgesAll * MaxEdges2, + "weightsOnEdge", MeshFileID, OnEdgeDecompR8, + WeightsOnEdgeID); if (Err != 0) LOG_CRITICAL("HorzMesh: error reading weightsOnEdge"); + WeightsOnEdgeH = HostArray2DReal("WeightsOnEdge", NEdgesSize, MaxEdges2); + deepCopy(WeightsOnEdgeH, TmpWeightsOnEdgeR8); + } // end readWeights //------------------------------------------------------------------------------ @@ -617,26 +572,11 @@ void HorzMesh::readCoriolis() { int Err; - int FCellID; - FCellH = HostArray1DR8("FCell", NCellsSize); - Err = IO::readArray(FCellH.data(), NCellsAll, "fCell", MeshFileID, - CellDecompR8, FCellID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading fCell"); + readCellArray(FCellH, "fCell"); - int FVertexID; - FVertexH = HostArray1DR8("FVertex", NVerticesSize); - Err = IO::readArray(FVertexH.data(), NVerticesAll, "fVertex", MeshFileID, - VertexDecompR8, FVertexID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading fVertex"); + readVertexArray(FVertexH, "fVertex"); - int FEdgeID; - FEdgeH = HostArray1DR8("FEdge", NEdgesSize); - Err = IO::readArray(FEdgeH.data(), NEdgesAll, "fEdge", MeshFileID, - EdgeDecompR8, FEdgeID); - if (Err != 0) - LOG_CRITICAL("HorzMesh: error reading fEdge"); + readEdgeArray(FEdgeH, "fEdge"); } // end readCoriolis @@ -644,7 +584,7 @@ void HorzMesh::readCoriolis() { // Compute the sign of edge contributions to a cell/vertex for each edge void HorzMesh::computeEdgeSign() { - EdgeSignOnCell = Array2DR8("EdgeSignOnCell", NCellsSize, MaxEdges); + EdgeSignOnCell = Array2DReal("EdgeSignOnCell", NCellsSize, MaxEdges); OMEGA_SCOPE(o_NEdgesOnCell, NEdgesOnCell); OMEGA_SCOPE(o_EdgesOnCell, EdgesOnCell); @@ -668,7 +608,7 @@ void HorzMesh::computeEdgeSign() { EdgeSignOnCellH = createHostMirrorCopy(EdgeSignOnCell); EdgeSignOnVertex = - Array2DR8("EdgeSignOnVertex", NVerticesSize, VertexDegree); + Array2DReal("EdgeSignOnVertex", NVerticesSize, VertexDegree); OMEGA_SCOPE(o_VertexDegree, VertexDegree); OMEGA_SCOPE(o_EdgesOnVertex, EdgesOnVertex); @@ -698,7 +638,7 @@ void HorzMesh::computeEdgeSign() { // and vertices void HorzMesh::setMasks(int NVertLevels) { - EdgeMask = Array2DR8("EdgeMask", NEdgesSize, NVertLevels); + EdgeMask = Array2DReal("EdgeMask", NEdgesSize, NVertLevels); OMEGA_SCOPE(O_EdgeMask, EdgeMask); @@ -718,8 +658,8 @@ void HorzMesh::setMasks(int NVertLevels) { // equations so viscosity and diffusion scale with mesh. void HorzMesh::setMeshScaling() { - MeshScalingDel2 = Array1DR8("MeshScalingDel2", NEdgesSize); - MeshScalingDel4 = Array1DR8("MeshScalingDel4", NEdgesSize); + MeshScalingDel2 = Array1DReal("MeshScalingDel2", NEdgesSize); + MeshScalingDel4 = Array1DReal("MeshScalingDel4", NEdgesSize); OMEGA_SCOPE(o_MeshScalingDel2, MeshScalingDel2); OMEGA_SCOPE(o_MeshScalingDel4, MeshScalingDel4); diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index 0b5f0c90588c..25be940e3466 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -35,6 +35,11 @@ class HorzMesh { void createDimensions(Decomp *MeshDecomp); + void readVertexArray(HostArray1DReal &VertexArrayH, + const std::string &MPASName); + void readEdgeArray(HostArray1DReal &EdgeArrayH, const std::string &MPASName); + void readCellArray(HostArray1DReal &CellArrayH, const std::string &MPASName); + void readCoordinates(); void readBottomDepth(); @@ -149,100 +154,106 @@ class HorzMesh { // Coordinates - HostArray1DR8 XCellH; ///< X Coordinates of cell centers (m) - HostArray1DR8 YCellH; ///< Y Coordinates of cell centers (m) - HostArray1DR8 ZCellH; ///< Z Coordinates of cell centers (m) - HostArray1DR8 LonCellH; ///< Longitude location of cell centers (radians) - HostArray1DR8 LatCellH; ///< Latitude location of cell centers (radians) + HostArray1DReal XCellH; ///< X Coordinates of cell centers (m) + HostArray1DReal YCellH; ///< Y Coordinates of cell centers (m) + HostArray1DReal ZCellH; ///< Z Coordinates of cell centers (m) + HostArray1DReal LonCellH; ///< Longitude location of cell centers (radians) + HostArray1DReal LatCellH; ///< Latitude location of cell centers (radians) - HostArray1DR8 XEdgeH; ///< X Coordinate of edge midpoints (m) - HostArray1DR8 YEdgeH; ///< Y Coordinate of edge midpoints (m) - HostArray1DR8 ZEdgeH; ///< Z Coordinate of edge midpoints (m) - HostArray1DR8 LonEdgeH; ///< Longitude location of edge midpoints (radians) - HostArray1DR8 LatEdgeH; ///< Latitude location of edge midpoints (radians) + HostArray1DReal XEdgeH; ///< X Coordinate of edge midpoints (m) + HostArray1DReal YEdgeH; ///< Y Coordinate of edge midpoints (m) + HostArray1DReal ZEdgeH; ///< Z Coordinate of edge midpoints (m) + HostArray1DReal LonEdgeH; ///< Longitude location of edge midpoints (radians) + HostArray1DReal LatEdgeH; ///< Latitude location of edge midpoints (radians) - HostArray1DR8 XVertexH; ///< X Coordinate of vertices (m) - HostArray1DR8 YVertexH; ///< Y Coordinate of vertices (m) - HostArray1DR8 ZVertexH; ///< Z Coordinate of vertices (m) - HostArray1DR8 LonVertexH; ///< Longitude location of vertices (radians) - HostArray1DR8 LatVertexH; ///< Latitude location of vertices (radians) + HostArray1DReal XVertexH; ///< X Coordinate of vertices (m) + HostArray1DReal YVertexH; ///< Y Coordinate of vertices (m) + HostArray1DReal ZVertexH; ///< Z Coordinate of vertices (m) + HostArray1DReal LonVertexH; ///< Longitude location of vertices (radians) + HostArray1DReal LatVertexH; ///< Latitude location of vertices (radians) // Mesh measurements - Array1DR8 AreaCell; ///< Area of each cell (m^2) - HostArray1DR8 AreaCellH; ///< Area of each cell (m^2) + Array1DReal AreaCell; ///< Area of each cell (m^2) + HostArray1DReal AreaCellH; ///< Area of each cell (m^2) - Array1DR8 AreaTriangle; ///< Area of each triangle in the dual grid (m^2) - HostArray1DR8 + Array1DReal AreaTriangle; ///< Area of each triangle in the dual grid (m^2) + HostArray1DReal AreaTriangleH; ///< Area of each triangle in the dual grid (m^2) - Array2DR8 KiteAreasOnVertex; ///< Area of the portions of each dual cell that - /// are part of each cellsOnVertex (m^2) - HostArray2DR8 + Array2DReal + KiteAreasOnVertex; ///< Area of the portions of each dual cell that + /// are part of each cellsOnVertex (m^2) + HostArray2DReal KiteAreasOnVertexH; ///< Area of the portions of each dual cell that /// are part of each cellsOnVertex (m^2) - Array1DR8 DvEdge; ///< Length of each edge, computed as the distance between - /// verticesOnEdge (m) - HostArray1DR8 DvEdgeH; ///< Length of each edge, computed as the distance - /// between verticesOnEdge (m) + Array1DReal + DvEdge; ///< Length of each edge, computed as the distance between + /// verticesOnEdge (m) + HostArray1DReal DvEdgeH; ///< Length of each edge, computed as the distance + /// between verticesOnEdge (m) - Array1DR8 DcEdge; ///< Length of each edge, computed as the distance between - /// CellsOnEdge (m) - HostArray1DR8 DcEdgeH; ///< Length of each edge, computed as the distance - /// between CellsOnEdge (m) + Array1DReal + DcEdge; ///< Length of each edge, computed as the distance between + /// CellsOnEdge (m) + HostArray1DReal DcEdgeH; ///< Length of each edge, computed as the distance + /// between CellsOnEdge (m) - Array1DR8 AngleEdge; ///< Angle the edge normal makes with local eastward - /// direction (radians) - HostArray1DR8 AngleEdgeH; ///< Angle the edge normal makes with local - /// eastward direction (radians) + Array1DReal AngleEdge; ///< Angle the edge normal makes with local eastward + /// direction (radians) + HostArray1DReal AngleEdgeH; ///< Angle the edge normal makes with local + /// eastward direction (radians) - HostArray1DR8 MeshDensityH; ///< Value of density function used to generate a - /// particular mesh at cell centers + HostArray1DReal + MeshDensityH; ///< Value of density function used to generate a + /// particular mesh at cell centers // Weights - Array2DR8 WeightsOnEdge; ///< Reconstruction weights associated with each of - /// the edgesOnEdge - HostArray2DR8 WeightsOnEdgeH; ///< Reconstruction weights associated with - /// each of the edgesOnEdge + Array2DReal + WeightsOnEdge; ///< Reconstruction weights associated with each of + /// the edgesOnEdge + HostArray2DReal WeightsOnEdgeH; ///< Reconstruction weights associated with + /// each of the edgesOnEdge // Coriolis - Array1DR8 FEdge; ///< Coriolis parameter at edges (radians s^-1) - HostArray1DR8 FEdgeH; ///< Coriolis parameter at edges (radians s^-1) + Array1DReal FEdge; ///< Coriolis parameter at edges (radians s^-1) + HostArray1DReal FEdgeH; ///< Coriolis parameter at edges (radians s^-1) - Array1DR8 FCell; ///< Coriolis parameter at cell centers (radians s^-1) - HostArray1DR8 FCellH; ///< Coriolis parameter at cell centers (radians s^-1) + Array1DReal FCell; ///< Coriolis parameter at cell centers (radians s^-1) + HostArray1DReal + FCellH; ///< Coriolis parameter at cell centers (radians s^-1) - Array1DR8 FVertex; ///< Coriolis parameter at vertices (radians s^-1) - HostArray1DR8 FVertexH; ///< Coriolis parameter at vertices (radians s^-1) + Array1DReal FVertex; ///< Coriolis parameter at vertices (radians s^-1) + HostArray1DReal FVertexH; ///< Coriolis parameter at vertices (radians s^-1) // Depth - Array1DR8 BottomDepth; ///< Depth of the bottom of the ocean (m) - HostArray1DR8 BottomDepthH; ///< Depth of the bottom of the ocean (m) + Array1DReal BottomDepth; ///< Depth of the bottom of the ocean (m) + HostArray1DReal BottomDepthH; ///< Depth of the bottom of the ocean (m) // Edge sign - Array2DR8 EdgeSignOnCell; ///< Sign of vector connecting cells - HostArray2DR8 EdgeSignOnCellH; ///< Sign of vector connecting cells + Array2DReal EdgeSignOnCell; ///< Sign of vector connecting cells + HostArray2DReal EdgeSignOnCellH; ///< Sign of vector connecting cells - Array2DR8 EdgeSignOnVertex; ///< Sign of vector connecting vertices - HostArray2DR8 EdgeSignOnVertexH; ///< Sign of vector connecting vertices + Array2DReal EdgeSignOnVertex; ///< Sign of vector connecting vertices + HostArray2DReal EdgeSignOnVertexH; ///< Sign of vector connecting vertices // Masks - Array2DR8 EdgeMask; ///< Mask to determine if computations should be - /// done on edge - HostArray2DR8 EdgeMaskH; ///< Mask to determine if computations should be - /// done on edge + Array2DReal EdgeMask; ///< Mask to determine if computations should be + /// done on edge + HostArray2DReal EdgeMaskH; ///< Mask to determine if computations should be + /// done on edge // Mesh scaling - Array1DR8 MeshScalingDel2; /// Coef to Laplacian mixing terms - HostArray1DR8 MeshScalingDel2H; /// Coef to Laplacian mixing terms + Array1DReal MeshScalingDel2; /// Coef to Laplacian mixing terms + HostArray1DReal MeshScalingDel2H; /// Coef to Laplacian mixing terms - Array1DR8 MeshScalingDel4; /// Coef to biharmonic mixing terms - HostArray1DR8 MeshScalingDel4H; /// Coef to biharmonic mixing terms + Array1DReal MeshScalingDel4; /// Coef to biharmonic mixing terms + HostArray1DReal MeshScalingDel4H; /// Coef to biharmonic mixing terms // Methods diff --git a/components/omega/src/ocn/HorzOperators.h b/components/omega/src/ocn/HorzOperators.h index bad06ba86072..9d26f255d5cd 100644 --- a/components/omega/src/ocn/HorzOperators.h +++ b/components/omega/src/ocn/HorzOperators.h @@ -36,9 +36,9 @@ class DivergenceOnCell { private: Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; - Array1DR8 DvEdge; - Array1DR8 AreaCell; - Array2DR8 EdgeSignOnCell; + Array1DReal DvEdge; + Array1DReal AreaCell; + Array2DReal EdgeSignOnCell; }; class GradientOnEdge { @@ -62,7 +62,7 @@ class GradientOnEdge { private: Array2DI4 CellsOnEdge; - Array1DR8 DcEdge; + Array1DReal DcEdge; }; class CurlOnVertex { @@ -96,9 +96,9 @@ class CurlOnVertex { private: I4 VertexDegree; Array2DI4 EdgesOnVertex; - Array1DR8 DcEdge; - Array1DR8 AreaTriangle; - Array2DR8 EdgeSignOnVertex; + Array1DReal DcEdge; + Array1DReal AreaTriangle; + Array2DReal EdgeSignOnVertex; }; class TangentialReconOnEdge { @@ -129,7 +129,7 @@ class TangentialReconOnEdge { private: Array1DI4 NEdgesOnEdge; Array2DI4 EdgesOnEdge; - Array2DR8 WeightsOnEdge; + Array2DReal WeightsOnEdge; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index ee4aaa46da94..a128a8258f59 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/components/omega/src/ocn/OceanState.cpp @@ -84,10 +84,10 @@ OceanState::OceanState( // Allocate state host arrays for (int I = 0; I < NTimeLevels; I++) { - LayerThicknessH[I] = HostArray2DR8("LayerThickness" + std::to_string(I), - NCellsSize, NVertLevels); - NormalVelocityH[I] = HostArray2DR8("NormalVelocity" + std::to_string(I), - NEdgesSize, NVertLevels); + LayerThicknessH[I] = HostArray2DReal("LayerThickness" + std::to_string(I), + NCellsSize, NVertLevels); + NormalVelocityH[I] = HostArray2DReal("NormalVelocity" + std::to_string(I), + NEdgesSize, NVertLevels); } // Create device arrays and copy host data @@ -296,11 +296,11 @@ void OceanState::defineFields() { StateGroupName); // Associate Field with data - Err = NormalVelocityField->attachData(NormalVelocity[CurLevel]); + Err = NormalVelocityField->attachData(NormalVelocity[CurLevel]); if (Err != 0) LOG_ERROR("Error attaching data array to field {}", NormalVelocityFldName); - Err = LayerThicknessField->attachData(LayerThickness[CurLevel]); + Err = LayerThicknessField->attachData(LayerThickness[CurLevel]); if (Err != 0) LOG_ERROR("Error attaching data array to field {}", LayerThicknessFldName); @@ -331,22 +331,32 @@ void OceanState::read(int StateFileID, I4 CellDecompR8, I4 EdgeDecompR8) { I4 Err; - // Read LayerThickness + // Read LayerThickness into a temporary double-precision array int LayerThicknessID; - Err = IO::readArray(LayerThicknessH[CurLevel].data(), NCellsAll, - "layerThickness", StateFileID, CellDecompR8, - LayerThicknessID); + HostArray2DR8 TmpLayerThicknessR8("TmpLayerThicknessR8", NCellsSize, + NVertLevels); + Err = IO::readArray(TmpLayerThicknessR8.data(), NCellsAll, "layerThickness", + StateFileID, CellDecompR8, LayerThicknessID); if (Err != 0) LOG_CRITICAL("OceanState: error reading layerThickness"); - // Read NormalVelocity + // Copy the thickness data into the final state array of user-specified + // precision + deepCopy(LayerThicknessH[CurLevel], TmpLayerThicknessR8); + + // Read NormalVelocity into a temporary double-precision array int NormalVelocityID; - Err = IO::readArray(NormalVelocityH[CurLevel].data(), NEdgesAll, - "normalVelocity", StateFileID, EdgeDecompR8, - NormalVelocityID); + HostArray2DR8 TmpNormalVelocityR8("TmpNormalVelocityR8", NEdgesSize, + NVertLevels); + Err = IO::readArray(TmpNormalVelocityR8.data(), NEdgesAll, "normalVelocity", + StateFileID, EdgeDecompR8, NormalVelocityID); if (Err != 0) LOG_CRITICAL("OceanState: error reading normalVelocity"); + // Copy the velocity data into the final state array of user-specified + // precision + deepCopy(NormalVelocityH[CurLevel], TmpNormalVelocityR8); + } // end read //------------------------------------------------------------------------------ @@ -397,10 +407,10 @@ void OceanState::updateTimeLevels() { // Update IOField data associations int Err = 0; - Err = Field::attachFieldData(NormalVelocityFldName, - NormalVelocity[CurLevel]); - Err = Field::attachFieldData(LayerThicknessFldName, - LayerThickness[CurLevel]); + Err = Field::attachFieldData(NormalVelocityFldName, + NormalVelocity[CurLevel]); + Err = Field::attachFieldData(LayerThicknessFldName, + LayerThickness[CurLevel]); } // end updateTimeLevels diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 7089defe327c..b695c7759737 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -35,7 +35,7 @@ class ThicknessFluxDivOnCell { /// The functor takes cell index, vertical chunk index, and thickness flux /// array as inputs, outputs the tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 ICell, I4 KChunk, - const Array2DR8 &ThicknessFlux, + const Array2DReal &ThicknessFlux, const Array2DReal &NormalVelEdge) const { const I4 KStart = KChunk * VecLength; @@ -62,9 +62,9 @@ class ThicknessFluxDivOnCell { private: Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; - Array1DR8 DvEdge; - Array1DR8 AreaCell; - Array2DR8 EdgeSignOnCell; + Array1DReal DvEdge; + Array1DReal AreaCell; + Array2DReal EdgeSignOnCell; }; /// Horizontal advection of potential vorticity defined on edges, for @@ -81,10 +81,10 @@ class PotentialVortHAdvOnEdge { /// thickness on edges, and normal velocity on edges as inputs, /// outputs the tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, - const Array2DR8 &NormRVortEdge, - const Array2DR8 &NormFEdge, - const Array2DR8 &FluxLayerThickEdge, - const Array2DR8 &NormVelEdge) const { + const Array2DReal &NormRVortEdge, + const Array2DReal &NormFEdge, + const Array2DReal &FluxLayerThickEdge, + const Array2DReal &NormVelEdge) const { const I4 KStart = KChunk * VecLength; Real VortTmp[VecLength] = {0}; @@ -112,7 +112,7 @@ class PotentialVortHAdvOnEdge { private: Array1DI4 NEdgesOnEdge; Array2DI4 EdgesOnEdge; - Array2DR8 WeightsOnEdge; + Array2DReal WeightsOnEdge; }; /// Gradient of kinetic energy defined on edges, for momentum equation @@ -126,7 +126,7 @@ class KEGradOnEdge { /// The functor takes edge index, vertical chunk index, and kinetic energy /// array as inputs, outputs the tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, - const Array2DR8 &KECell) const { + const Array2DReal &KECell) const { const I4 KStart = KChunk * VecLength; const I4 JCell0 = CellsOnEdge(IEdge, 0); @@ -141,7 +141,7 @@ class KEGradOnEdge { private: Array2DI4 CellsOnEdge; - Array1DR8 DcEdge; + Array1DReal DcEdge; }; /// Gradient of sea surface height defined on edges multipled by gravitational @@ -171,9 +171,9 @@ class SSHGradOnEdge { } private: - R8 Grav = 9.80665_Real; + Real Grav = 9.80665_Real; Array2DI4 CellsOnEdge; - Array1DR8 DcEdge; + Array1DReal DcEdge; }; /// Laplacian horizontal mixing, for momentum equation @@ -181,7 +181,7 @@ class VelocityDiffusionOnEdge { public: bool Enabled; - R8 ViscDel2; + Real ViscDel2; /// constructor declaration VelocityDiffusionOnEdge(const HorzMesh *Mesh); @@ -190,8 +190,8 @@ class VelocityDiffusionOnEdge { /// divergence of horizontal velocity (defined at cell centers) and relative /// vorticity (defined at vertices), outputs tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, - const Array2DR8 &DivCell, - const Array2DR8 &RVortVertex) const { + const Array2DReal &DivCell, + const Array2DReal &RVortVertex) const { const I4 KStart = KChunk * VecLength; const I4 ICell0 = CellsOnEdge(IEdge, 0); @@ -218,10 +218,10 @@ class VelocityDiffusionOnEdge { private: Array2DI4 CellsOnEdge; Array2DI4 VerticesOnEdge; - Array1DR8 DcEdge; - Array1DR8 DvEdge; - Array1DR8 MeshScalingDel2; - Array2DR8 EdgeMask; + Array1DReal DcEdge; + Array1DReal DvEdge; + Array1DReal MeshScalingDel2; + Array2DReal EdgeMask; }; /// Biharmonic horizontal mixing, for momentum equation @@ -238,8 +238,8 @@ class VelocityHyperDiffOnEdge { /// the laplacian of divergence of horizontal velocity and the laplacian of /// the relative vorticity, outputs tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, - const Array2DR8 &Del2DivCell, - const Array2DR8 &Del2RVortVertex) const { + const Array2DReal &Del2DivCell, + const Array2DReal &Del2RVortVertex) const { const I4 KStart = KChunk * VecLength; const I4 ICell0 = CellsOnEdge(IEdge, 0); @@ -266,10 +266,10 @@ class VelocityHyperDiffOnEdge { private: Array2DI4 CellsOnEdge; Array2DI4 VerticesOnEdge; - Array1DR8 DcEdge; - Array1DR8 DvEdge; - Array1DR8 MeshScalingDel4; - Array2DR8 EdgeMask; + Array1DReal DcEdge; + Array1DReal DvEdge; + Array1DReal MeshScalingDel4; + Array2DReal EdgeMask; }; // Tracer horizontal advection term @@ -280,8 +280,8 @@ class TracerHorzAdvOnCell { TracerHorzAdvOnCell(const HorzMesh *Mesh); KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, - I4 KChunk, const Array2DR8 &NormVelEdge, - const Array3DR8 &HTracersOnEdge) const { + I4 KChunk, const Array2DReal &NormVelEdge, + const Array3DReal &HTracersOnEdge) const { const I4 KStart = KChunk * VecLength; const Real InvAreaCell = 1._Real / AreaCell(ICell); @@ -308,9 +308,9 @@ class TracerHorzAdvOnCell { Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; Array2DI4 CellsOnEdge; - Array2DR8 EdgeSignOnCell; - Array1DR8 DvEdge; - Array1DR8 AreaCell; + Array2DReal EdgeSignOnCell; + Array1DReal DvEdge; + Array1DReal AreaCell; }; // Tracer horizontal diffusion term @@ -322,9 +322,10 @@ class TracerDiffOnCell { TracerDiffOnCell(const HorzMesh *Mesh); - KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, - I4 KChunk, const Array3DR8 &TracerCell, - const Array2DR8 &MeanLayerThickEdge) const { + KOKKOS_FUNCTION void + operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, + const Array3DReal &TracerCell, + const Array2DReal &MeanLayerThickEdge) const { const I4 KStart = KChunk * VecLength; const Real InvAreaCell = 1._Real / AreaCell(ICell); @@ -359,11 +360,11 @@ class TracerDiffOnCell { Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; Array2DI4 CellsOnEdge; - Array2DR8 EdgeSignOnCell; - Array1DR8 DvEdge; - Array1DR8 DcEdge; - Array1DR8 AreaCell; - Array1DR8 MeshScalingDel2; + Array2DReal EdgeSignOnCell; + Array1DReal DvEdge; + Array1DReal DcEdge; + Array1DReal AreaCell; + Array1DReal MeshScalingDel2; }; // Tracer biharmonic horizontal mixing term @@ -377,7 +378,7 @@ class TracerHyperDiffOnCell { KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, - const Array3DR8 &TrDel2Cell) const { + const Array3DReal &TrDel2Cell) const { const I4 KStart = KChunk * VecLength; const Real InvAreaCell = 1._Real / AreaCell(ICell); @@ -411,11 +412,11 @@ class TracerHyperDiffOnCell { Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; Array2DI4 CellsOnEdge; - Array2DR8 EdgeSignOnCell; - Array1DR8 DvEdge; - Array1DR8 DcEdge; - Array1DR8 AreaCell; - Array1DR8 MeshScalingDel4; + Array2DReal EdgeSignOnCell; + Array1DReal DvEdge; + Array1DReal DcEdge; + Array1DReal AreaCell; + Array1DReal MeshScalingDel4; }; /// A class that can be used to calculate the thickness and diff --git a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h index e7854010e087..572bae9ba2d7 100644 --- a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.h @@ -53,10 +53,10 @@ class KineticAuxVars { private: Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; - Array2DR8 EdgeSignOnCell; - Array1DR8 DcEdge; - Array1DR8 DvEdge; - Array1DR8 AreaCell; + Array2DReal EdgeSignOnCell; + Array1DReal DcEdge; + Array1DReal DvEdge; + Array1DReal AreaCell; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h index b34e42c41caa..1436dc773af2 100644 --- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h @@ -97,10 +97,10 @@ class TracerAuxVars { Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; Array2DI4 CellsOnEdge; - Array2DR8 EdgeSignOnCell; - Array1DR8 DcEdge; - Array1DR8 DvEdge; - Array1DR8 AreaCell; + Array2DReal EdgeSignOnCell; + Array1DReal DcEdge; + Array1DReal DvEdge; + Array1DReal AreaCell; }; } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h index 4b922c0100bf..e359942aa8f1 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h @@ -95,15 +95,15 @@ class VelocityDel2AuxVars { private: Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; - Array2DR8 EdgeSignOnCell; - Array1DR8 DcEdge; - Array1DR8 DvEdge; - Array1DR8 AreaCell; + Array2DReal EdgeSignOnCell; + Array1DReal DcEdge; + Array1DReal DvEdge; + Array1DReal AreaCell; Array2DI4 EdgesOnVertex; Array2DI4 CellsOnEdge; Array2DI4 VerticesOnEdge; - Array2DR8 EdgeSignOnVertex; - Array1DR8 AreaTriangle; + Array2DReal EdgeSignOnVertex; + Array1DReal AreaTriangle; I4 VertexDegree; }; diff --git a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h index d6cba92ca2b3..82dc8c489c02 100644 --- a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.h @@ -83,12 +83,12 @@ class VorticityAuxVars { I4 VertexDegree; Array2DI4 CellsOnVertex; Array2DI4 EdgesOnVertex; - Array2DR8 EdgeSignOnVertex; - Array1DR8 DcEdge; - Array2DR8 KiteAreasOnVertex; - Array1DR8 AreaTriangle; + Array2DReal EdgeSignOnVertex; + Array1DReal DcEdge; + Array2DReal KiteAreasOnVertex; + Array1DReal AreaTriangle; Array2DI4 VerticesOnEdge; - Array1DR8 FVertex; + Array1DReal FVertex; }; } // namespace OMEGA diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index eda0b500aefa..00770bc41567 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -191,7 +191,7 @@ void TimeStepper::updateThicknessByTend(OceanState *State1, int TimeLevel1, const auto &LayerThickTend = Tend->LayerThicknessTend; const int NVertLevels = LayerThickTend.extent_int(1); - Real CoeffSeconds; + R8 CoeffSeconds; Coeff.get(CoeffSeconds, TimeUnits::Seconds); parallelFor( @@ -213,7 +213,7 @@ void TimeStepper::updateVelocityByTend(OceanState *State1, int TimeLevel1, const auto &NormalVelTend = Tend->NormalVelocityTend; const int NVertLevels = NormalVelTend.extent_int(1); - Real CoeffSeconds; + R8 CoeffSeconds; Coeff.get(CoeffSeconds, TimeUnits::Seconds); parallelFor( @@ -247,7 +247,7 @@ void TimeStepper::updateTracersByTend(const Array3DReal &NextTracers, const int NTracers = TracerTend.extent(0); const int NVertLevels = TracerTend.extent(2); - Real CoeffSeconds; + R8 CoeffSeconds; Err = Coeff.get(CoeffSeconds, TimeUnits::Seconds); parallelFor( @@ -286,7 +286,7 @@ void TimeStepper::accumulateTracersUpdate(const Array3DReal &AccumTracer, const int NTracers = TracerTend.extent(0); const int NVertLevels = TracerTend.extent(2); - Real CoeffSeconds; + R8 CoeffSeconds; int Err = Coeff.get(CoeffSeconds, TimeUnits::Seconds); parallelFor( diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index a11b3a65c3e6..91be32856b4a 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -6,6 +6,9 @@ ##################### function(add_omega_test test_name exe_name source_files mpi_args) + # Copy extra arguments to a local variable so they can be treated + # as a list + set(extra_args ${ARGN}) # Create the executable add_executable(${exe_name} ${source_files}) @@ -14,10 +17,35 @@ function(add_omega_test test_name exe_name source_files mpi_args) target_link_libraries( ${exe_name} PRIVATE - ${OMEGA_LIB_NAME} OmegaLibFlags ) + list(FIND extra_args single_precision use_sp) + + # Link the single precision version if requested and not already configured to + # use single precision + if (NOT use_sp EQUAL -1) + if (OMEGA_SINGLE_PRECISION) + target_link_libraries( + ${exe_name} + PRIVATE + ${OMEGA_LIB_NAME} + ) + else() + target_link_libraries( + ${exe_name} + PRIVATE + ${OMEGA_LIB_NAME}_SP + ) + endif() + else() + target_link_libraries( + ${exe_name} + PRIVATE + ${OMEGA_LIB_NAME} + ) + endif() + set_target_properties(${exe_name} PROPERTIES LINKER_LANGUAGE C) # Add the test command @@ -261,6 +289,14 @@ target_compile_definitions( TENDENCYTERMS_TEST_PLANE ) +add_omega_test( + TEND_PLANE_SINGLE_PRECISION_TEST + testTendencyTermsPlaneSinglePrecision.exe + ocn/TendencyTermsTest.cpp + "-n 8;--cpu-bind=cores" + single_precision +) + add_omega_test( TEND_SPHERE_TEST testTendencyTermsSphere.exe diff --git a/components/omega/test/base/ReductionsTest.cpp b/components/omega/test/base/ReductionsTest.cpp index 9f8935292fc0..a9d24154bb10 100644 --- a/components/omega/test/base/ReductionsTest.cpp +++ b/components/omega/test/base/ReductionsTest.cpp @@ -41,8 +41,8 @@ int main(int argc, char *argv[]) { I4 MyInt4 = 1, MyResI4 = 0; I8 MyInt8 = 2, MyResI8 = 0; R4 MyR4 = 3.000001, MyResR4 = 0.0; - R8 MyR8 = 4.0000000000001, MyResR8 = 0.0, MyResReal = 0.0; - Real MyReal = 5.000001; + R8 MyR8 = 4.0000000000001, MyResR8 = 0.0; + Real MyReal = 5.000001, MyResReal = 0.0; // test SUM for scalars err = globalSum(&MyInt4, Comm, &MyResI4); @@ -212,10 +212,10 @@ int main(int argc, char *argv[]) { LocalSum2D = complex(t1 + t2, t2 - ((t1 + t2) - t1)); } } - Sum1DR8 = real(LocalSum1D); - Sum2DR8 = real(LocalSum2D); - MyResReal = 0.0; - err = globalSum(HostArr1DR8, Comm, &MyResReal); + Sum1DR8 = real(LocalSum1D); + Sum2DR8 = real(LocalSum2D); + MyResR8 = 0.0; + err = globalSum(HostArr1DR8, Comm, &MyResR8); // perform serial sum across all MPI tasks complex SerialSum(0.0, 0.0); for (i = 0; i < MySize; i++) { @@ -227,22 +227,22 @@ int main(int argc, char *argv[]) { } expR8 = real(SerialSum); res = "FAIL"; - if (err == 0 && MyResReal == expR8) + if (err == 0 && MyResR8 == expR8) res = "PASS"; else RetVal += 1; printf("Global sum A1DR8: %s (exp,act=%.13lf,%.13lf)\n", res, expR8, - MyResReal); + MyResR8); - err = globalSum(HostArr2DR8, Comm, &MyResReal); + err = globalSum(HostArr2DR8, Comm, &MyResR8); expR8 = Sum2DR8 * MySize; res = "FAIL"; - if (err == 0 && MyResReal == expR8) + if (err == 0 && MyResR8 == expR8) res = "PASS"; else RetVal += 1; printf("Global sum A2DR8: %s (exp,act=%.13lf,%.13lf)\n", res, expR8, - MyResReal); + MyResR8); //========================================================================== // test MIN, MAX of scalars diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 1b8bae8d7a35..7e5f9e1d8f79 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -183,42 +183,43 @@ struct TestSetupSphere { } KOKKOS_FUNCTION Real velocityX(Real Lon, Real Lat) const { - return -Radius * std::pow(std::sin(Lon), 2) * std::pow(std::cos(Lat), 3); + return -std::pow(std::sin(Lon), 2) * std::pow(std::cos(Lat), 3); } KOKKOS_FUNCTION Real velocityY(Real Lon, Real Lat) const { - return -4 * Radius * std::sin(Lon) * std::cos(Lon) * - std::pow(std::cos(Lat), 3) * std::sin(Lat); + return -4 * std::sin(Lon) * std::cos(Lon) * std::pow(std::cos(Lat), 3) * + std::sin(Lat); } KOKKOS_FUNCTION Real relativeVorticity(Real Lon, Real Lat) const { return -4 * std::pow(std::cos(Lon), 2) * std::pow(std::cos(Lat), 2) * - std::sin(Lat); + std::sin(Lat) / Radius; } KOKKOS_FUNCTION Real divergence(Real Lon, Real Lat) const { return std::sin(Lon) * std::cos(Lon) * std::pow(std::cos(Lat), 2) * - (20 * std::pow(std::sin(Lat), 2) - 6); + (20 * std::pow(std::sin(Lat), 2) - 6) / Radius; } KOKKOS_FUNCTION Real velocityDel2X(Real Lon, Real Lat) const { - return 1 / Radius * + return 1 / (Radius * Radius) * (std::pow(std::cos(Lon), 2) - std::pow(std::sin(Lon), 2)) * std::cos(Lat) * (20 * std::pow(std::sin(Lat), 2) - 6) + - 4 / Radius * std::pow(cos(Lon), 2) * + 4 / (Radius * Radius) * std::pow(cos(Lon), 2) * (std::pow(cos(Lat), 3) - 2 * std::cos(Lat) * std::pow(sin(Lat), 2)); } KOKKOS_FUNCTION Real velocityDel2Y(Real Lon, Real Lat) const { - return 1 / Radius * std::sin(Lon) * std::cos(Lon) * std::sin(Lat) * - std::cos(Lat) * (80 * std::pow(std::cos(Lat), 2) - 28) + - 8 / Radius * std::sin(Lon) * std::cos(Lon) * std::sin(Lat) * - std::cos(Lat); + return 1 / (Radius * Radius) * std::sin(Lon) * std::cos(Lon) * + std::sin(Lat) * std::cos(Lat) * + (80 * std::pow(std::cos(Lat), 2) - 28) + + 8 / (Radius * Radius) * std::sin(Lon) * std::cos(Lon) * + std::sin(Lat) * std::cos(Lat); } KOKKOS_FUNCTION Real velocityDel2Div(Real Lon, Real Lat) const { - return 1 / (Radius * Radius) * + return 1 / (Radius * Radius * Radius) * (-2 * std::sin(Lon) * std::cos(Lon) * (28 * std::pow(sin(Lat), 2) - 8) + std::sin(Lon) * std::cos(Lon) * @@ -228,7 +229,7 @@ struct TestSetupSphere { } KOKKOS_FUNCTION Real velocityDel2Curl(Real Lon, Real Lat) const { - return 1 / (Radius * Radius) * + return 1 / (Radius * Radius * Radius) * (-std::sin(Lat) * (std::pow(std::cos(Lat), 2) * (56 * std::pow(cos(Lon), 2) - 40) - 2 * (std::pow(cos(Lon), 2) * diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 24c4b9083e4f..74be0bec52ca 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -36,20 +36,20 @@ struct TestSetupPlane { Real Lx = 1; Real Ly = std::sqrt(3) / 2; - Real ExpectedDivErrorLInf = 0.00124886886594453264; - Real ExpectedDivErrorL2 = 0.00124886886590977139; - Real ExpectedPVErrorLInf = 0.00807347170900282914; - Real ExpectedPVErrorL2 = 0.00794755105765788429; - Real ExpectedGradErrorLInf = 0.00125026071878537952; - Real ExpectedGradErrorL2 = 0.00134354611117262161; - Real ExpectedLaplaceErrorLInf = 0.00113090174765822192; - Real ExpectedLaplaceErrorL2 = 0.00134324628763667899; - Real ExpectedTrHAdvErrorLInf = 0.00205864372747571571; - Real ExpectedTrHAdvErrorL2 = 0.00172418025417940784; - Real ExpectedTrDel2ErrorLInf = 0.00334357193650093847; - Real ExpectedTrDel2ErrorL2 = 0.00290978146207349032; - Real ExpectedTrDel4ErrorLInf = 0.00508833446725232875; - Real ExpectedTrDel4ErrorL2 = 0.00523080740758275625; + ErrorMeasures ExpectedDivErrors = {0.00124886886594453264, + 0.00124886886590977139}; + ErrorMeasures ExpectedPVErrors = {0.00807347170900282914, + 0.00794755105765788429}; + ErrorMeasures ExpectedGradErrors = {0.00125026071878537952, + 0.00134354611117262161}; + ErrorMeasures ExpectedLaplaceErrors = {0.00113090174765822192, + 0.00134324628763667899}; + ErrorMeasures ExpectedTrHAdvErrors = {0.00205864372747571571, + 0.00172418025417940784}; + ErrorMeasures ExpectedTrDel2Errors = {0.00334357193650093847, + 0.00290978146207349032}; + ErrorMeasures ExpectedTrDel4Errors = {0.00508833446725232875, + 0.00523080740758275625}; KOKKOS_FUNCTION Real vectorX(Real X, Real Y) const { return std::sin(2 * Pi * X / Lx) * std::cos(2 * Pi * Y / Ly); @@ -154,20 +154,20 @@ struct TestSetupSphere { Real Pi = M_PI; - Real ExpectedDivErrorLInf = 0.0136595773989796766; - Real ExpectedDivErrorL2 = 0.00367052484586384131; - Real ExpectedPVErrorLInf = 0.0219217796608757037; - Real ExpectedPVErrorL2 = 0.0122537418367830303; - Real ExpectedGradErrorLInf = 0.00187912292540623471; - Real ExpectedGradErrorL2 = 0.00149841802817334935; - Real ExpectedLaplaceErrorLInf = 0.281930203304510130; - Real ExpectedLaplaceErrorL2 = 0.270530313560271740; - Real ExpectedTrHAdvErrorLInf = 0.0132310202299444034; - Real ExpectedTrHAdvErrorL2 = 0.0038523368564029538; - Real ExpectedTrDel2ErrorLInf = 0.0486107109846934185; - Real ExpectedTrDel2ErrorL2 = 0.00507514214194892694; - Real ExpectedTrDel4ErrorLInf = 0.000819552466009620408; - Real ExpectedTrDel4ErrorL2 = 0.00064700084412871962; + ErrorMeasures ExpectedDivErrors = {0.0136595773989796766, + 0.00367052484586384131}; + ErrorMeasures ExpectedPVErrors = {0.0219217796608757037, + 0.0122537418367830303}; + ErrorMeasures ExpectedGradErrors = {0.00187912292540623471, + 0.00149841802817334935}; + ErrorMeasures ExpectedLaplaceErrors = {0.281930203304510130, + 0.270530313560271740}; + ErrorMeasures ExpectedTrHAdvErrors = {0.0132310202299444034, + 0.0038523368564029538}; + ErrorMeasures ExpectedTrDel2Errors = {0.0486107109846934185, + 0.00507514214194892694}; + ErrorMeasures ExpectedTrDel4Errors = {0.000819552466009620408, + 0.00064700084412871962}; KOKKOS_FUNCTION Real vectorX(Real Lon, Real Lat) const { return -Radius * std::pow(std::sin(Lon), 2) * std::pow(std::cos(Lat), 3); @@ -291,10 +291,10 @@ int testThickFluxDiv(int NVertLevels, Real RTol) { ExactThickFluxDiv, Geom, Mesh, OnCell, NVertLevels, ExchangeHalos::No); // Set input array - Array2DR8 ThickFluxEdge("ThickFluxEdge", Mesh->NEdgesSize, NVertLevels); + Array2DReal ThickFluxEdge("ThickFluxEdge", Mesh->NEdgesSize, NVertLevels); // TODO(mwarusz) temporary fix for this test - Array2DR8 OnesEdge("OnesEdge", Mesh->NEdgesSize, NVertLevels); + Array2DReal OnesEdge("OnesEdge", Mesh->NEdgesSize, NVertLevels); deepCopy(OnesEdge, 1); Err += setVectorEdge( @@ -320,15 +320,8 @@ int testThickFluxDiv(int NVertLevels, Real RTol) { OnCell, NVertLevels); // Check error values - if (!isApprox(TFDivErrors.LInf, Setup.ExpectedDivErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: ThickFluxDiv LInf FAIL"); - } - - if (!isApprox(TFDivErrors.L2, Setup.ExpectedDivErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: ThickFluxDiv L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "ThickFluxDiv", TFDivErrors, + Setup.ExpectedDivErrors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: ThickFluxDiv PASS"); @@ -359,25 +352,26 @@ int testPotVortHAdv(int NVertLevels, Real RTol) { ExchangeHalos::No); // Set input arrays - Array2DR8 NormRelVortEdge("NormRelVortEdge", Mesh->NEdgesSize, NVertLevels); + Array2DReal NormRelVortEdge("NormRelVortEdge", Mesh->NEdgesSize, + NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normRelVort(X, Y); }, NormRelVortEdge, Geom, Mesh, OnEdge, NVertLevels); - Array2DR8 NormPlanetVortEdge("NormPlanetVortEdge", Mesh->NEdgesSize, - NVertLevels); + Array2DReal NormPlanetVortEdge("NormPlanetVortEdge", Mesh->NEdgesSize, + NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.normPlanetVort(X, Y); }, NormPlanetVortEdge, Geom, Mesh, OnEdge, NVertLevels); - Array2DR8 LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); + Array2DReal LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.layerThick(X, Y); }, LayerThickEdge, Geom, Mesh, OnEdge, NVertLevels); - Array2DR8 NormVelEdge("NormVelEdge", Mesh->NEdgesSize, NVertLevels); + Array2DReal NormVelEdge("NormVelEdge", Mesh->NEdgesSize, NVertLevels); Err += setVectorEdge( KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { @@ -402,15 +396,8 @@ int testPotVortHAdv(int NVertLevels, Real RTol) { Mesh, OnEdge, NVertLevels); // Check error values - if (!isApprox(PotVortHAdvErrors.LInf, Setup.ExpectedPVErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: PotVortHAdv LInf FAIL"); - } - - if (!isApprox(PotVortHAdvErrors.L2, Setup.ExpectedPVErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: PotVortHAdv L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "PotVortHAdv", PotVortHAdvErrors, + Setup.ExpectedPVErrors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: PotVortHAdv PASS"); @@ -459,15 +446,8 @@ int testKEGrad(int NVertLevels, Real RTol) { NVertLevels); // Check error values - if (!isApprox(KEGradErrors.LInf, Setup.ExpectedGradErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: KEGrad LInf FAIL"); - } - - if (!isApprox(KEGradErrors.L2, Setup.ExpectedGradErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: KEGrad L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "KEGrad", KEGradErrors, + Setup.ExpectedGradErrors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: KEGrad PASS"); @@ -516,19 +496,8 @@ int testSSHGrad(int NVertLevels, Real RTol) { NVertLevels); // Check error values - if (!isApprox(SSHGradErrors.LInf, Setup.ExpectedGradErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: SSHGrad LInf FAIL"); - } - - if (!isApprox(SSHGradErrors.L2, Setup.ExpectedGradErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: SSHGrad L2 FAIL"); - } - - if (Err == 0) { - LOG_INFO("TendencyTermsTest: SSHGrad PASS"); - } + Err += checkErrors("TendencyTermsTest", "SSHGrad", SSHGradErrors, + Setup.ExpectedGradErrors, RTol); return Err; } // end testSSHGrad @@ -592,15 +561,8 @@ int testVelDiff(int NVertLevels, Real RTol) { NVertLevels); // Check error values - if (!isApprox(VelDiffErrors.LInf, Setup.ExpectedLaplaceErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: VelDiff LInf FAIL"); - } - - if (!isApprox(VelDiffErrors.L2, Setup.ExpectedLaplaceErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: VelDiff L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "VelDiff", VelDiffErrors, + Setup.ExpectedLaplaceErrors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: VelDiff PASS"); @@ -670,16 +632,8 @@ int testVelHyperDiff(int NVertLevels, Real RTol) { Mesh, OnEdge, NVertLevels); // Check error values - if (!isApprox(VelHyperDiffErrors.LInf, Setup.ExpectedLaplaceErrorLInf, - RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: VelHyperDiff LInf FAIL"); - } - - if (!isApprox(VelHyperDiffErrors.L2, Setup.ExpectedLaplaceErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: VelHyperDiff L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "VelHyperDiff", VelHyperDiffErrors, + Setup.ExpectedLaplaceErrors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: VelHyperDiff PASS"); @@ -705,7 +659,7 @@ int testTracerHorzAdvOnCell(int NVertLevels, int NTracers, Real RTol) { ExchangeHalos::No); // Set input arrays - Array2DR8 NormalVelocity("NormalVelocity", Mesh->NEdgesSize, NVertLevels); + Array2DReal NormalVelocity("NormalVelocity", Mesh->NEdgesSize, NVertLevels); Err += setVectorEdge( KOKKOS_LAMBDA(Real(&VecField)[2], Real X, Real Y) { @@ -714,7 +668,7 @@ int testTracerHorzAdvOnCell(int NVertLevels, int NTracers, Real RTol) { }, NormalVelocity, EdgeComponent::Normal, Geom, Mesh, NVertLevels); - Array3DR8 HTrOnEdge("HTrOnEdge", NTracers, Mesh->NEdgesSize, NVertLevels); + Array3DReal HTrOnEdge("HTrOnEdge", NTracers, Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return -Setup.layerThick(X, Y); }, @@ -735,15 +689,8 @@ int testTracerHorzAdvOnCell(int NVertLevels, int NTracers, Real RTol) { Err += computeErrors(TrHAdvErrors, NumTrFluxDiv, ExactTrFluxDiv, Mesh, OnCell, NVertLevels, NTracers); - if (!isApprox(TrHAdvErrors.LInf, Setup.ExpectedTrHAdvErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: TracerHorzAdv LInf FAIL"); - } - - if (!isApprox(TrHAdvErrors.L2, Setup.ExpectedTrHAdvErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: TracerHorzAdv L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "TracerHorzAdv", TrHAdvErrors, + Setup.ExpectedTrHAdvErrors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: TracerHorzAdv PASS"); @@ -769,13 +716,14 @@ int testTracerDiffOnCell(int NVertLevels, int NTracers, Real RTol) { ExchangeHalos::No); // Set input arrays - Array3DR8 TracerCell("TracerCell", NTracers, Mesh->NCellsSize, NVertLevels); + Array3DReal TracerCell("TracerCell", NTracers, Mesh->NCellsSize, + NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, TracerCell, Geom, Mesh, OnCell, NVertLevels, NTracers); - Array2DR8 LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); + Array2DReal LayerThickEdge("LayerThickEdge", Mesh->NEdgesSize, NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, @@ -798,15 +746,8 @@ int testTracerDiffOnCell(int NVertLevels, int NTracers, Real RTol) { Err += computeErrors(TrDiffErrors, NumTracerDiff, ExactTracerDiff, Mesh, OnCell, NVertLevels, NTracers); - if (!isApprox(TrDiffErrors.LInf, Setup.ExpectedTrDel2ErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: TracerDiff LInf FAIL"); - } - - if (!isApprox(TrDiffErrors.L2, Setup.ExpectedTrDel2ErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: TracerDiff L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "TracerDiff", TrDiffErrors, + Setup.ExpectedTrDel2Errors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: TracerDiff PASS"); @@ -832,7 +773,8 @@ int testTracerHyperDiffOnCell(int NVertLevels, int NTracers, Real RTol) { ExchangeHalos::No); // Set input arrays - Array3DR8 TrDel2Cell("TracerCell", NTracers, Mesh->NCellsSize, NVertLevels); + Array3DReal TrDel2Cell("TracerCell", NTracers, Mesh->NCellsSize, + NVertLevels); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarC(X, Y); }, @@ -854,15 +796,8 @@ int testTracerHyperDiffOnCell(int NVertLevels, int NTracers, Real RTol) { computeErrors(TrHyperDiffErrors, NumTracerHyperDiff, ExactTracerHyperDiff, Mesh, OnCell, NVertLevels, NTracers); - if (!isApprox(TrHyperDiffErrors.LInf, Setup.ExpectedTrDel4ErrorLInf, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: TracerHyperDiff LInf FAIL"); - } - - if (!isApprox(TrHyperDiffErrors.L2, Setup.ExpectedTrDel4ErrorL2, RTol)) { - Err++; - LOG_ERROR("TendencyTermsTest: TracerHyperDiff L2 FAIL"); - } + Err += checkErrors("TendencyTermsTest", "TracerHyperDiff", TrHyperDiffErrors, + Setup.ExpectedTrDel4Errors, RTol); if (Err == 0) { LOG_INFO("TendencyTermsTest: TracerHyperDiff PASS"); @@ -938,7 +873,7 @@ int tendencyTermsTest(const std::string &mesh = DefaultMeshFile) { int NVertLevels = 16; int NTracers = 3; - const Real RTol = sizeof(Real) == 4 ? 1e-2 : 1e-10; + const Real RTol = sizeof(Real) == 4 ? 2e-2 : 1e-10; Err += testThickFluxDiv(NVertLevels, RTol); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 68629f86b046..bdb4d806e4a0 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -242,7 +242,7 @@ int initTimeStepperTest(const std::string &mesh) { // of steps int adjustTimeStep(TimeStepper *Stepper, Real TimeEnd) { TimeInterval TimeStep = Stepper->getTimeStep(); - Real TimeStepSeconds; + R8 TimeStepSeconds; TimeStep.get(TimeStepSeconds, TimeUnits::Seconds); const int NSteps = std::ceil(TimeEnd / TimeStepSeconds); @@ -310,7 +310,7 @@ int testTimeStepper(const std::string &Name, TimeStepperType Type, const Real TimeEnd = 1; - const Real BaseTimeStepSeconds = 0.2; + const R8 BaseTimeStepSeconds = 0.2; // This creates global exact solution and needs to be done only once const static bool CallOnlyOnce = [=]() { @@ -318,7 +318,7 @@ int testTimeStepper(const std::string &Name, TimeStepperType Type, return true; }(); - Real TimeStepSeconds = BaseTimeStepSeconds; + R8 TimeStepSeconds = BaseTimeStepSeconds; // Convergence loop for (int RefLevel = 0; RefLevel < NRefinements; ++RefLevel) {