From e5e028d6ee844135981082d731d9b113c9209d38 Mon Sep 17 00:00:00 2001 From: Philip W Jones Date: Tue, 7 Apr 2026 09:14:12 -0700 Subject: [PATCH 1/2] read omega HorzMesh from input IOStream Uses an input IOStream to read the input horizontal mesh variables. The same filename is used by Decomp. If the Decomp optional argument for a mesh file is supplied, that value overrides the mesh stream filename. All physical mesh fields are defined for IO and are added to a field group HorzMeshIn. Some cleanup of the HorzMesh test driver is included to reduce test output and conform to omega code conventions. User and developer docs have been updated with these changes. --- components/omega/configs/Default.yml | 10 + components/omega/doc/devGuide/Decomp.md | 13 +- components/omega/doc/devGuide/HorzMesh.md | 27 +- components/omega/doc/userGuide/Decomp.md | 7 +- components/omega/doc/userGuide/HorzMesh.md | 21 +- components/omega/src/base/Decomp.cpp | 24 +- components/omega/src/base/Decomp.h | 2 +- components/omega/src/infra/IOStream.cpp | 15 + components/omega/src/infra/IOStream.h | 7 + components/omega/src/ocn/HorzMesh.cpp | 1058 ++++++++++++----- components/omega/src/ocn/HorzMesh.h | 135 ++- components/omega/src/ocn/HorzOperators.cpp | 11 +- components/omega/src/ocn/OceanInit.cpp | 2 +- components/omega/test/infra/IOStreamTest.cpp | 2 +- .../omega/test/ocn/AuxiliaryStateTest.cpp | 8 +- .../omega/test/ocn/AuxiliaryVarsTest.cpp | 13 +- components/omega/test/ocn/EosTest.cpp | 15 +- components/omega/test/ocn/HorzMeshTest.cpp | 535 ++++----- .../omega/test/ocn/HorzOperatorsTest.cpp | 19 +- components/omega/test/ocn/StateTest.cpp | 8 +- components/omega/test/ocn/TendenciesTest.cpp | 10 +- .../omega/test/ocn/TendencyTermsTest.cpp | 8 +- components/omega/test/ocn/TracersTest.cpp | 11 +- components/omega/test/ocn/VertAdvTest.cpp | 11 +- components/omega/test/ocn/VertCoordTest.cpp | 11 +- components/omega/test/ocn/VertMixTest.cpp | 14 +- .../test/timeStepping/TimeStepperTest.cpp | 15 +- 27 files changed, 1239 insertions(+), 773 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index eb393b4d0982..a8cbf2417039 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -82,6 +82,16 @@ Omega: Alpha: 5.0 Exponent: 2.0 IOStreams: + HorzMeshIn: + UsePointerFile: false + Filename: OmegaMesh.nc + Mode: read + Precision: double + Freq: 1 + FreqUnits: OnStartup + UseStartEnd: false + Contents: + - HorzMeshIn InitialVertCoord: UsePointerFile: false Filename: OmegaMesh.nc diff --git a/components/omega/doc/devGuide/Decomp.md b/components/omega/doc/devGuide/Decomp.md index 3cce573d7012..3c4bfcf27351 100644 --- a/components/omega/doc/devGuide/Decomp.md +++ b/components/omega/doc/devGuide/Decomp.md @@ -20,11 +20,14 @@ with the call: Decomp::init(); ``` This must be called very early in the init process, just after initializing -the MachEnv, Config and IO. Mesh information is first read using parallel -IO into an equally-spaced linear decomposition, then partitioned by METIS or -ParMETIS into a more optimal decomposition. The parallel ParMETIS -implementation should be used to reduce the memory footprint for -high-resolution configurations. Although the algorithmic method is similar, +the MachEnv, Config and IO. By default, the mesh information is read from +the filename defined in the [HorzMeshIn IOStream](#omega-user-horz-mesh). +However, an optional filename argument to the init function above can be used +to override the input configuration (eg for unit testing). Mesh information is +first read using parallel IO into an equally-spaced linear decomposition, then +partitioned by METIS or ParMETIS into a more optimal decomposition. The +parallel ParMETIS implementation should be used to reduce the memory footprint +for high-resolution configurations. Although the algorithmic method is similar, ParMETIS leads to a somewhat different partitioning than the serial METIS version. The serial version should give the same partition as the serial offline gpmetis tool used by MPAS when the same version of the METIS library diff --git a/components/omega/doc/devGuide/HorzMesh.md b/components/omega/doc/devGuide/HorzMesh.md index e22565cc3e7a..a12b140e2231 100644 --- a/components/omega/doc/devGuide/HorzMesh.md +++ b/components/omega/doc/devGuide/HorzMesh.md @@ -7,17 +7,20 @@ Specification](https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf). A mesh object is created by the `init` method, which assumes that `Decomp` has already been initialized. ```c++ -Err = OMEGA::Decomp::init(); -Err = OMEGA::HorzMesh::init(); +OMEGA::Decomp::init(); +OMEGA::HorzMesh::init(ModelClock); ``` +The model clock should be the ocean model clock and is required by the IOStream +routines that read the HorzMesh variables from the input stream +[HorzMeshIn](#omega-user-horz-mesh). The constructor replicates the subdomain mesh cell/edge/vertex counts and connectivity information from Decomp so this information can be passed among the computational routines, alongside the other local mesh information. It then -creates several parallel I/O decompositions and reads in the remaining subdomain -mesh information. Finally, any mesh information needed on the device is copied -from the host to a device Kokkos array. Arrays such as the coordinate variables, -which are not involved in tendency calculations, are not transferred to the -device. These tasks are organized into several private methods. Eventually, +creates reads most of the remaining mesh variables (coordinates, lengths, areas) +from the input HorzMeshIn IOStream. A few variables, like masks and operator +scaling are computed internally. After reading, the init routine fills all +halos and synchronizes the host and device copies of the mesh variables. +These tasks are organized into several private methods. Eventually, dependent mesh variables will be computed from the minimum set of required mesh information. @@ -45,11 +48,9 @@ OMEGA::parallelFor({HMesh->NCellsOwned,HMesh->MaxEdges}, ``` For member variables that are host arrays, variable names are appended with an -`H`. Array variable names not ending in `H` are device arrays. The copy from -host to device array is performed in the constructor via: -```c++ -AreaCell = OMEGA::createDeviceMirrorCopy(AreaCellH); -``` +`H`. Array variable names not ending in `H` are device arrays. The device arrays are deallocated by the `HorzMesh::clear()` method, which is -necessary before calling `Kokkos::finalize`. +necessary before calling `Kokkos::finalize`. Because Field and IOStream also +have references to HorzMesh variables, the `IOStream::clear()` and +`Field::clear()` should also be called to ensure all device arrays are removed. diff --git a/components/omega/doc/userGuide/Decomp.md b/components/omega/doc/userGuide/Decomp.md index 93844628b4fa..3ef0b3636f80 100644 --- a/components/omega/doc/userGuide/Decomp.md +++ b/components/omega/doc/userGuide/Decomp.md @@ -41,9 +41,10 @@ An input mesh file must be provided that contains at a minimum VerticesOnCell, CellsOnEdge, EdgesOnEdge, CellsOnVertex, EdgesOnVertex. Again, a full description of the mesh is given in the [Developer's Guide](#omega-dev-decomp). -The mesh information is read via parallel IO into an initial linear domain -decomposition and then is partitioned by METIS and rearranged into the -final METIS parallel decomposition. +The file name for this input file is extracted from the HorzMeshIn input +mesh stream located in the IOStreams section of the input configuration. +This value can be overridden by the driver routine if needed (eg for unit +testing). Once the mesh is decomposed, all of the mesh index arrays are stored in a Decomp named Default which can be retrieved as described in the diff --git a/components/omega/doc/userGuide/HorzMesh.md b/components/omega/doc/userGuide/HorzMesh.md index 2160b51c82ae..8b68c4164c67 100644 --- a/components/omega/doc/userGuide/HorzMesh.md +++ b/components/omega/doc/userGuide/HorzMesh.md @@ -16,7 +16,26 @@ additional mesh variables, which are not required for Decomp. Currently, the Mesh class reads in all variables from the MPAS mesh specification except those read by [Decomp](#omega-dev-decomp). -This includes the following variables: +These are all read using the input [IOStream](#omega-user-iostreams) HorzMeshIn +```yaml + IOStreams: + HorzMeshIn: + UsePointerFile: false + Filename: OmegaMesh.nc + Mode: read + Precision: double + Freq: 1 + FreqUnits: OnStartup + UseStartEnd: false + Contents: + - HorzMeshIn +``` +Only the Filename should be changed by the user to point to the relevant input +mesh file. The mesh Filename is sometimes overridden by the driver routine in +the case of unit tests using an optional argument to the +[decomposition](#omega-dev-decomp). The input contents are defined by the +HorzMeshIn [FieldGroup](#omega-dev-field) that currently includes the following +variables: | Variable Name | Description | Units | | ------------- | ----------- | ----- | diff --git a/components/omega/src/base/Decomp.cpp b/components/omega/src/base/Decomp.cpp index 5bcc6326fb5a..6641eb935dc4 100644 --- a/components/omega/src/base/Decomp.cpp +++ b/components/omega/src/base/Decomp.cpp @@ -398,7 +398,7 @@ void readMesh(const int MeshFileID, // file ID for open mesh file // Initialize the decomposition and create the default decomposition with // (currently) one partition per MPI task using selected partition method -void Decomp::init(const std::string &MeshFileName) { +void Decomp::init(const std::string &InMeshFileName) { bool TimerFlag = Pacer::start("Decomp init", 0); Error Err; // default successful return code @@ -421,6 +421,26 @@ void Decomp::init(const std::string &MeshFileName) { PartMethod Method = getPartMethodFromStr(DecompMethodStr); + // If no mesh filename supplied, get the filename from the input + // HorzMeshIn stream configuration + std::string MeshFileNameTmp = InMeshFileName; + if (InMeshFileName == "") { + Config StreamConfigTmp("IOStreams"); + Err = OmegaConfig->get(StreamConfigTmp); + CHECK_ERROR_ABORT(Err, "Decomp: IOStreams not found in input Config"); + + Config HorzMeshInTmp("HorzMeshIn"); // retrieve HorzMeshIn stream info + Err = StreamConfigTmp.get(HorzMeshInTmp); + CHECK_ERROR_ABORT(Err, "Decomp: HorzMeshIn stream not found in Config"); + + Err = HorzMeshInTmp.get("Filename", MeshFileNameTmp); + CHECK_ERROR_ABORT(Err, "Decomp: Could not get mesh filename from Config"); + // Testing - re-retrieve a fresh copy of the stream and mesh filename + Err = OmegaConfig->get(StreamConfigTmp); + Err = StreamConfigTmp.get(HorzMeshInTmp); + Err = HorzMeshInTmp.get("Filename", MeshFileNameTmp); + } + // Retrieve the default machine environment MachEnv *DefEnv = MachEnv::getDefault(); @@ -429,7 +449,7 @@ void Decomp::init(const std::string &MeshFileName) { // Create the default decomposition and set pointer to it Decomp::DefaultDecomp = Decomp::create("Default", DefEnv, NParts, Method, - InHaloWidth, MeshFileName); + InHaloWidth, MeshFileNameTmp); TimerFlag = Pacer::stop("Decomp init", 0) && TimerFlag; if (!TimerFlag) diff --git a/components/omega/src/base/Decomp.h b/components/omega/src/base/Decomp.h index 1d4895d50158..1e98ef34cd9d 100644 --- a/components/omega/src/base/Decomp.h +++ b/components/omega/src/base/Decomp.h @@ -264,7 +264,7 @@ class Decomp { /// Initializes Omega decomposition info and creates the default /// decomposition based on the default MachEnv and configuration /// options. - static void init(const std::string &MeshFileName = "OmegaMesh.nc"); + static void init(const std::string &InMeshFileName = ""); // Creates a new decomposition using the constructor and puts it in the // AllDecomps map diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp index 3647c26b7a49..45df426a18d6 100644 --- a/components/omega/src/infra/IOStream.cpp +++ b/components/omega/src/infra/IOStream.cpp @@ -135,6 +135,16 @@ IOStream::get(const std::string &StreamName ///< [in] name of stream to retrieve } // End get stream +//------------------------------------------------------------------------------ +// Changes the filename associated with a stream (eg during unit testing) +void IOStream::changeFilename( + const std::string &StreamName, ///< [in] name of stream to modify + const std::string &NewFilename ///< [in] new filename for stream +) { + auto StreamPtr = get(StreamName); + StreamPtr->Filename = NewFilename; +} // End changeFilename + //------------------------------------------------------------------------------ // Adds a field to the contents of a stream. Because streams may be created // before all Fields have been defined, we only store the name. Validity @@ -2279,6 +2289,11 @@ Error IOStream::readStream( if (Mode != IO::ModeRead) ABORT_ERROR("IOStream read: cannot read stream defined as output stream"); + // Validate the stream if not already validated (this also expands group + // names into field names). + if (!validate()) + ABORT_ERROR("Unable to validate stream {}", Name); + // If it is not time to read, return if (!ForceRead) { if (!MyAlarm.isRinging() and !OnStartup) diff --git a/components/omega/src/infra/IOStream.h b/components/omega/src/infra/IOStream.h index e43fa63a4aad..14efa94cadcd 100644 --- a/components/omega/src/infra/IOStream.h +++ b/components/omega/src/infra/IOStream.h @@ -294,6 +294,13 @@ class IOStream { get(const std::string &StreamName ///< [in] name of stream to retrieve ); + //--------------------------------------------------------------------------- + /// Changes the filename in a stream (typically during unit testing) + static void changeFilename( + const std::string &StreamName, ///< [in] name of stream to modify + const std::string &NewFilename ///< [in] new file name + ); + //--------------------------------------------------------------------------- /// Adds a field to the contents of a stream. Because streams may be created /// before all Fields have been defined, we only store the name. Validity diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index ec3ff19e3336..4692e515716f 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -10,13 +10,13 @@ //===----------------------------------------------------------------------===// #include "HorzMesh.h" -#include "Config.h" #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" #include "Error.h" -#include "IO.h" -#include "Logging.h" +#include "Field.h" +#include "Halo.h" +#include "IOStream.h" #include "OmegaKokkos.h" namespace OMEGA { @@ -29,28 +29,28 @@ std::map> HorzMesh::AllHorzMeshes; // Initialize the mesh. Assumes that Decomp and VertCoord have already been // initialized. -void HorzMesh::init() { - - Error Err; // default successful error code +void HorzMesh::init(const Clock *ModelClock //< [in] Model clock for IO alarms +) { // Retrieve the default decomposition Decomp *DefDecomp = Decomp::getDefault(); // Create the default mesh and set pointer to it - HorzMesh::DefaultHorzMesh = create("Default", DefDecomp); + HorzMesh::DefaultHorzMesh = create("Default", DefDecomp, ModelClock); } //------------------------------------------------------------------------------ // Construct a new local mesh given a decomposition HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh - Decomp *MeshDecomp //< [in] Decomp for the new mesh - ) - : CellID(MeshDecomp->CellID) ///< global cell ID for each local cell -{ + Decomp *MeshDecomp, //< [in] Decomp for the new mesh + const Clock *ModelClock //< [in] Model clock for IO alarms +) { + + // Set mesh name based on input MeshName = Name; - // Retrieve mesh files name from Decomp + // Mesh filename should be the same as that used for the decomposition MeshFileName = MeshDecomp->MeshFileName; // Retrieve mesh cell/edge/vertex totals from Decomp @@ -68,6 +68,7 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh NEdgesSize = MeshDecomp->NEdgesSize; MaxCellsOnEdge = MeshDecomp->MaxCellsOnEdge; MaxEdges = MeshDecomp->MaxEdges; + MaxEdges2 = 2 * MaxEdges; NEdgesGlobal = MeshDecomp->NEdgesGlobal; NVerticesHalo = MeshDecomp->NVerticesHalo; @@ -78,6 +79,9 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh VertexDegree = MeshDecomp->VertexDegree; // Retrieve connectivity arrays from Decomp + + CellID = MeshDecomp->CellID; + CellsOnCellH = MeshDecomp->CellsOnCellH; EdgesOnCellH = MeshDecomp->EdgesOnCellH; NEdgesOnCellH = MeshDecomp->NEdgesOnCellH; @@ -100,62 +104,48 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh CellsOnVertex = MeshDecomp->CellsOnVertex; EdgesOnVertex = MeshDecomp->EdgesOnVertex; - // Open the mesh file for reading (assume IO has already been initialized) - IO::openFileRead(MeshFileID, MeshFileName); - // Create Omega Dimensions associated with this mesh createDimensions(MeshDecomp); - // Temporary - remove when IOStreams implemented - // Create the parallel IO decompositions required to read in mesh variables - initParallelIO(MeshDecomp); - - // Read x/y/z and lon/lat coordinates for cells, edges, and vertices - readCoordinates(); - - // Read the mesh areas, lengths, and angles - readMeasurements(); - - // Read the edge mesh weights - readWeights(); - - // Read the Coriolis parameter at the cells, edges, and vertices - readCoriolis(); - - // Destroy the parallel IO decompositions - finalizeParallelIO(); + // Allocate remaining device arrays and define all mesh fields for associated + // I/O and initialization. Equivalent host arrays will be allocated and + // initialized after the fields are filled on the device. + defineMeshFields(); - // Copy host data to device - copyToDevice(); + // Read the input mesh stream to fill most of the fields + Metadata ReqMetaData; // empty - we do not require any global metadata + Error Err = IOStream::read("HorzMeshIn", ModelClock, ReqMetaData); + CHECK_ERROR_ABORT(Err, "HorzMesh: error reading input mesh stream"); - // TODO: add ability to compute (rather than read in) - // dependent mesh quantities + // Complete read arrays (fill halos and duplicate on host) + completeReadArrays(); + // Compute additional mesh quantities // Compute EdgeSignOnCells and EdgeSignOnVertex computeEdgeSign(); - // set mesh scaling coefficients - setMeshScaling(); + // Compute mesh scaling coefficients + computeMeshScaling(); } // end horizontal mesh constructor -/// Creates a new mesh by calling the constructor and puts it in the -/// AllHorzMeshes map +//------------------------------------------------------------------------------ +// Creates a new mesh by calling the constructor and adding to the map of +// all horizontal meshes HorzMesh *HorzMesh::create(const std::string &Name, //< [in] Name for new mesh - Decomp *MeshDecomp //< [in] Decomp for the new mesh + Decomp *MeshDecomp, //< [in] Decomp for new mesh + const Clock *ModelClock //< [in] Model clock for IO ) { // Check to see if a mesh of the same name already exists and // if so, exit with an error - if (AllHorzMeshes.find(Name) != AllHorzMeshes.end()) { - LOG_ERROR("Attempted to create a HorzMesh with name {} but a HorzMesh of " - "that name already exists", - Name); - return nullptr; - } + if (AllHorzMeshes.find(Name) != AllHorzMeshes.end()) + ABORT_ERROR("Attempted to create a HorzMesh with name {} but a HorzMesh " + "of that name already exists", + Name); // create a new mesh on the heap and put it in a map of // unique_ptrs, which will manage its lifetime - auto *NewHorzMesh = new HorzMesh(Name, MeshDecomp); + auto *NewHorzMesh = new HorzMesh(Name, MeshDecomp, ModelClock); AllHorzMeshes.emplace(Name, NewHorzMesh); return NewHorzMesh; @@ -183,7 +173,7 @@ void HorzMesh::erase(std::string InName // [in] name of mesh to remove void HorzMesh::clear() { AllHorzMeshes.clear(); // removes all meshes from the list and in - // the porcess, calls the destructors for each + // the process, calls the destructors for each DefaultHorzMesh = nullptr; // prevent dangling pointer } // end clear @@ -193,12 +183,14 @@ void HorzMesh::clear() { void HorzMesh::createDimensions(Decomp *MeshDecomp) { // Create non-distributed dimensions - // If these have already been created (eg not the default mesh), skip + // If these have already been created, skip since the dims should not change if (!Dimension::exists("MaxCellsOnEdge")) auto MaxCellsOnEdgeDim = Dimension::create("MaxCellsOnEdge", MaxCellsOnEdge); if (!Dimension::exists("MaxEdges")) auto MaxEdgesDim = Dimension::create("MaxEdges", MaxEdges); + if (!Dimension::exists("MaxEdges2")) + auto MaxEdgesDim = Dimension::create("MaxEdges2", MaxEdges2); if (!Dimension::exists("VertexDegree")) auto VertexDegreeDim = Dimension::create("VertexDegree", VertexDegree); @@ -264,254 +256,10 @@ void HorzMesh::createDimensions(Decomp *MeshDecomp) { } // end createDimensions -//------------------------------------------------------------------------------ -// Initialize the parallel IO decompositions for the mesh variables -void HorzMesh::initParallelIO(Decomp *MeshDecomp) { - - I4 NDims = 1; - IO::Rearranger Rearr = IO::RearrBox; - - // Create the IO decomp for arrays with (NCells) dimensions - std::vector CellDims{MeshDecomp->NCellsGlobal}; - std::vector CellID(NCellsAll); - for (int Cell = 0; Cell < NCellsAll; ++Cell) { - CellID[Cell] = MeshDecomp->CellIDH(Cell) - 1; - } - - CellDecompR8 = IO::createDecomp(IO::IOTypeR8, NDims, CellDims, NCellsAll, - CellID, Rearr); - - // Create the IO decomp for arrays with (NEdges) dimensions - std::vector EdgeDims{MeshDecomp->NEdgesGlobal}; - std::vector EdgeID(NEdgesAll); - for (int Edge = 0; Edge < NEdgesAll; ++Edge) { - EdgeID[Edge] = MeshDecomp->EdgeIDH(Edge) - 1; - } - - EdgeDecompR8 = IO::createDecomp(IO::IOTypeR8, NDims, EdgeDims, NEdgesAll, - EdgeID, Rearr); - - // Create the IO decomp for arrays with (NVertices) dimensions - std::vector VertexDims{MeshDecomp->NVerticesGlobal}; - std::vector VertexID(NVerticesAll); - for (int Vertex = 0; Vertex < NVerticesAll; ++Vertex) { - VertexID[Vertex] = MeshDecomp->VertexIDH(Vertex) - 1; - } - - VertexDecompR8 = IO::createDecomp(IO::IOTypeR8, NDims, VertexDims, - NVerticesAll, VertexID, Rearr); - - // Create the IO decomp for arrays with (NEdges, 2*MaxEdges) dimensions - NDims = 2; - MaxEdges2 = 2 * MaxEdges; - std::vector OnEdgeDims2{MeshDecomp->NEdgesGlobal, MaxEdges2}; - I4 OnEdgeSize2 = NEdgesAll * MaxEdges2; - std::vector OnEdgeOffset2(OnEdgeSize2, -1); - for (int Edge = 0; Edge < NEdgesAll; Edge++) { - for (int i = 0; i < MaxEdges2; i++) { - I4 GlobalID = EdgeID[Edge] * MaxEdges2 + i; - - OnEdgeOffset2[Edge * MaxEdges2 + i] = GlobalID; - } - } - - OnEdgeDecompR8 = IO::createDecomp(IO::IOTypeR8, NDims, OnEdgeDims2, - OnEdgeSize2, OnEdgeOffset2, Rearr); - - // Create the IO decomp for arrays with (NVertices, VertexDegree) dimensions - std::vector OnVertexDims{MeshDecomp->NVerticesGlobal, VertexDegree}; - I4 OnVertexSize = NVerticesAll * VertexDegree; - std::vector OnVertexOffset(OnVertexSize, -1); - for (int Vertex = 0; Vertex < NVerticesAll; Vertex++) { - for (int i = 0; i < VertexDegree; i++) { - I4 GlobalID = VertexID[Vertex] * VertexDegree + i; - OnVertexOffset[Vertex * VertexDegree + i] = GlobalID; - } - } - - OnVertexDecompR8 = IO::createDecomp(IO::IOTypeR8, NDims, OnVertexDims, - OnVertexSize, OnVertexOffset, Rearr); - -} // end initParallelIO - -//------------------------------------------------------------------------------ -// Destroy parallel decompositions -void HorzMesh::finalizeParallelIO() { - - IO::destroyDecomp(CellDecompR8); - IO::destroyDecomp(EdgeDecompR8); - IO::destroyDecomp(VertexDecompR8); - IO::destroyDecomp(OnEdgeDecompR8); - IO::destroyDecomp(OnVertexDecompR8); - -} // end finalizeParallelIO - -// Read 1D vertex array -void HorzMesh::readVertexArray(HostArray1DReal &VertexArrayH, - const std::string &MPASName) { - Error Err; - - std::string OmegaName; - std::transform(MPASName.begin(), MPASName.end(), OmegaName.begin(), - [](unsigned char c) { return std::toupper(c); }); - - // Temporary double precision array for reading - HostArray1DR8 TmpArrayR8(OmegaName + "Tmp", NVerticesSize); - int ArrayID; - Err = IO::readArray(TmpArrayR8.data(), NVerticesAll, MPASName, MeshFileID, - VertexDecompR8, ArrayID); - CHECK_ERROR_ABORT(Err, "HorzMesh: error reading {}", MPASName); - - // Create host array of desired precision and copy the read data into it - VertexArrayH = HostArray1DReal(OmegaName + "H", NVerticesSize); - deepCopy(VertexArrayH, TmpArrayR8); -} - -// Read 1D edge array -void HorzMesh::readEdgeArray(HostArray1DReal &EdgeArrayH, - const std::string &MPASName) { - Error Err; - - std::string OmegaName; - std::transform(MPASName.begin(), MPASName.end(), OmegaName.begin(), - [](unsigned char c) { return std::toupper(c); }); - - // Temporary double precision array for reading - HostArray1DR8 TmpArrayR8(OmegaName + "Tmp", NEdgesSize); - int ArrayID; - Err = IO::readArray(TmpArrayR8.data(), NEdgesAll, MPASName, MeshFileID, - EdgeDecompR8, ArrayID); - CHECK_ERROR_ABORT(Err, "HorzMesh: error reading {}", MPASName); - - // Create host array of desired precision and copy the read data into it - EdgeArrayH = HostArray1DReal(OmegaName + "H", NEdgesSize); - deepCopy(EdgeArrayH, TmpArrayR8); -} - -// Read 1D cell array -void HorzMesh::readCellArray(HostArray1DReal &CellArrayH, - const std::string &MPASName) { - Error Err; - - std::string OmegaName; - std::transform(MPASName.begin(), MPASName.end(), OmegaName.begin(), - [](unsigned char c) { return std::toupper(c); }); - - // Temporary double precision array for reading - HostArray1DR8 TmpArrayR8(OmegaName + "Tmp", NCellsSize); - int ArrayID; - Err = IO::readArray(TmpArrayR8.data(), NCellsAll, MPASName, MeshFileID, - CellDecompR8, ArrayID); - CHECK_ERROR_ABORT(Err, "HorzMesh: error reading {}", MPASName); - - // Create host array of desired precision and copy the read data into it - CellArrayH = HostArray1DReal(OmegaName + "H", NCellsSize); - deepCopy(CellArrayH, TmpArrayR8); -} - -//------------------------------------------------------------------------------ -// Read x/y/z and lon/lat coordinates for cells, edges, and vertices -void HorzMesh::readCoordinates() { - - // 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 mesh areas (cell, triangle, and kite), -// lengths (between centers and vertices), and edge angles -void HorzMesh::readMeasurements() { - - readCellArray(AreaCellH, "areaCell"); - - readVertexArray(AreaTriangleH, "areaTriangle"); - - readEdgeArray(DvEdgeH, "dvEdge"); - - readEdgeArray(DcEdgeH, "dcEdge"); - - readEdgeArray(AngleEdgeH, "angleEdge"); - - readCellArray(MeshDensityH, "meshDensity"); - - // not using helper function since it kiteAreas is a 2d array - Error Err; - - // Read into a temporary double precision array - int KiteAreasOnVertexID; - - HostArray2DR8 TmpKiteAreasOnVertexR8("KiteAreasOnVertex", NVerticesSize, - VertexDegree); - Err = IO::readArray(TmpKiteAreasOnVertexR8.data(), - NVerticesAll * VertexDegree, "kiteAreasOnVertex", - MeshFileID, OnVertexDecompR8, KiteAreasOnVertexID); - CHECK_ERROR_ABORT(Err, "HorzMesh: error reading kiteAreasOnVertex"); - - // Create and fill array with Real precision - KiteAreasOnVertexH = - HostArray2DReal("KiteAreasOnVertex", NVerticesSize, VertexDegree); - deepCopy(KiteAreasOnVertexH, TmpKiteAreasOnVertexR8); - -} // end readMeasurements - -//------------------------------------------------------------------------------ -// Read the edge weights used in the discrete potential vorticity flux term -void HorzMesh::readWeights() { - - Error Err; - - int WeightsOnEdgeID; - HostArray2DR8 TmpWeightsOnEdgeR8("WeightsOnEdge", NEdgesSize, MaxEdges2); - Err = IO::readArray(TmpWeightsOnEdgeR8.data(), NEdgesAll * MaxEdges2, - "weightsOnEdge", MeshFileID, OnEdgeDecompR8, - WeightsOnEdgeID); - CHECK_ERROR_ABORT(Err, "HorzMesh: error reading weightsOnEdge"); - - WeightsOnEdgeH = HostArray2DReal("WeightsOnEdge", NEdgesSize, MaxEdges2); - deepCopy(WeightsOnEdgeH, TmpWeightsOnEdgeR8); - -} // end readWeights - -//------------------------------------------------------------------------------ -// Read the Coriolis parameter at the cells, edges, and vertices -void HorzMesh::readCoriolis() { - - readCellArray(FCellH, "fCell"); - - readVertexArray(FVertexH, "fVertex"); - - readEdgeArray(FEdgeH, "fEdge"); - -} // end readCoriolis - //------------------------------------------------------------------------------ // Compute the sign of edge contributions to a cell/vertex for each edge void HorzMesh::computeEdgeSign() { - EdgeSignOnCell = Array2DReal("EdgeSignOnCell", NCellsSize, MaxEdges); - OMEGA_SCOPE(o_NEdgesOnCell, NEdgesOnCell); OMEGA_SCOPE(o_EdgesOnCell, EdgesOnCell); OMEGA_SCOPE(o_CellsOnEdge, CellsOnEdge); @@ -561,10 +309,7 @@ void HorzMesh::computeEdgeSign() { //------------------------------------------------------------------------------ // Set mesh scaling coefficients for mixing terms in momentum and tracer // equations so viscosity and diffusion scale with mesh. -void HorzMesh::setMeshScaling() { - - MeshScalingDel2 = Array1DReal("MeshScalingDel2", NEdgesSize); - MeshScalingDel4 = Array1DReal("MeshScalingDel4", NEdgesSize); +void HorzMesh::computeMeshScaling() { OMEGA_SCOPE(o_MeshScalingDel2, MeshScalingDel2); OMEGA_SCOPE(o_MeshScalingDel4, MeshScalingDel4); @@ -580,27 +325,73 @@ void HorzMesh::setMeshScaling() { MeshScalingDel2H = createHostMirrorCopy(MeshScalingDel2); MeshScalingDel4H = createHostMirrorCopy(MeshScalingDel4); -} // end setMeshScaling +} // end computeMeshScaling //------------------------------------------------------------------------------ -// Perform copy to device for mesh variables -void HorzMesh::copyToDevice() { - - AreaCell = createDeviceMirrorCopy(AreaCellH); - AreaTriangle = createDeviceMirrorCopy(AreaTriangleH); - KiteAreasOnVertex = createDeviceMirrorCopy(KiteAreasOnVertexH); - DcEdge = createDeviceMirrorCopy(DcEdgeH); - DvEdge = createDeviceMirrorCopy(DvEdgeH); - AngleEdge = createDeviceMirrorCopy(AngleEdgeH); - WeightsOnEdge = createDeviceMirrorCopy(WeightsOnEdgeH); - FVertex = createDeviceMirrorCopy(FVertexH); - FEdge = createDeviceMirrorCopy(FEdgeH); - XCell = createDeviceMirrorCopy(XCellH); - YCell = createDeviceMirrorCopy(YCellH); - XEdge = createDeviceMirrorCopy(XEdgeH); - YEdge = createDeviceMirrorCopy(YEdgeH); - -} // end copyToDevice +// Finish filling arrays read from a file by filling halos and copying to host +void HorzMesh::completeReadArrays() { + + // Get precomputed halo information by name - mesh name is assumed to be + // same as halo name (including Default) + Halo *HorzMeshHalo = Halo::get(MeshName); + + // Fill halos for all arrays read from file + HorzMeshHalo->exchangeFullArrayHalo(XCell, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(YCell, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(ZCell, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(XEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(YEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(ZEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(XVertex, OnVertex); + HorzMeshHalo->exchangeFullArrayHalo(YVertex, OnVertex); + HorzMeshHalo->exchangeFullArrayHalo(ZVertex, OnVertex); + HorzMeshHalo->exchangeFullArrayHalo(LatCell, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(LonCell, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(LatEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(LonEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(LatVertex, OnVertex); + HorzMeshHalo->exchangeFullArrayHalo(LonVertex, OnVertex); + HorzMeshHalo->exchangeFullArrayHalo(AreaCell, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(MeshDensity, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(AreaTriangle, OnVertex); + HorzMeshHalo->exchangeFullArrayHalo(KiteAreasOnVertex, OnVertex); + HorzMeshHalo->exchangeFullArrayHalo(DcEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(DvEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(AngleEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(WeightsOnEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(FCell, OnCell); + HorzMeshHalo->exchangeFullArrayHalo(FEdge, OnEdge); + HorzMeshHalo->exchangeFullArrayHalo(FVertex, OnVertex); + + // Create host copies + XCellH = createHostMirrorCopy(XCell); + YCellH = createHostMirrorCopy(YCell); + ZCellH = createHostMirrorCopy(ZCell); + XEdgeH = createHostMirrorCopy(XEdge); + YEdgeH = createHostMirrorCopy(YEdge); + ZEdgeH = createHostMirrorCopy(ZEdge); + XVertexH = createHostMirrorCopy(XVertex); + YVertexH = createHostMirrorCopy(YVertex); + ZVertexH = createHostMirrorCopy(ZVertex); + LatCellH = createHostMirrorCopy(LatCell); + LonCellH = createHostMirrorCopy(LonCell); + LatEdgeH = createHostMirrorCopy(LatEdge); + LonEdgeH = createHostMirrorCopy(LonEdge); + LatVertexH = createHostMirrorCopy(LatVertex); + LonVertexH = createHostMirrorCopy(LonVertex); + AreaCellH = createHostMirrorCopy(AreaCell); + MeshDensityH = createHostMirrorCopy(MeshDensity); + AreaTriangleH = createHostMirrorCopy(AreaTriangle); + KiteAreasOnVertexH = createHostMirrorCopy(KiteAreasOnVertex); + DcEdgeH = createHostMirrorCopy(DcEdge); + DvEdgeH = createHostMirrorCopy(DvEdge); + AngleEdgeH = createHostMirrorCopy(AngleEdge); + WeightsOnEdgeH = createHostMirrorCopy(WeightsOnEdge); + FCellH = createHostMirrorCopy(FCell); + FEdgeH = createHostMirrorCopy(FEdge); + FVertexH = createHostMirrorCopy(FVertex); + +} // end completeReadArrays //------------------------------------------------------------------------------ // Get default mesh @@ -614,18 +405,631 @@ HorzMesh *HorzMesh::get(const std::string Name ///< [in] Name of mesh // look for an instance of this name auto it = AllHorzMeshes.find(Name); - // if found, return the mesh pointer - if (it != AllHorzMeshes.end()) { - return it->second.get(); + // if not found, abort + if (it == AllHorzMeshes.end()) + ABORT_ERROR("HorzMesh::get: Mesh {} not found. Mesh has not been defined " + "or has been removed", + Name); + + // found the mesh, return the pointer + return it->second.get(); - // otherwise print error and return null pointer - } else { - LOG_ERROR("HorzMesh::get: Attempt to retrieve non-existent mesh:"); - LOG_ERROR("{} has not been defined or has been removed", Name); - return nullptr; - } } // end get mesh +//------------------------------------------------------------------------------ +// Define all Horz mesh fields and associated mesh stream for any mesh I/O +void HorzMesh::defineMeshFields() { + + // First create a field group for mesh fields - simply HorzMesh for default, + // but add name of mesh for any other instances. We do not add the suffix to + // the field names because we expect any file containing the mesh fields + // to use the standard field names. + // We currently define a mesh group only for fields that are read in after + // the decomposition (HorzMeshIn). Other mesh groups may be defined later + // if fields are required in output, but connectivity arrays are not + // available after the mesh Decomposition has been defined. + + std::string MeshSuffix = ""; + if (MeshName != "Default") + MeshSuffix = MeshName; + + std::string MeshGroupName = "HorzMeshIn" + MeshSuffix; + auto MeshGroupIn = FieldGroup::create(MeshGroupName); + + // Check whether the file name in the input stream matches that of + // the Decomp. They are typically different only for unit testing where + // an internal variable overrides the default configuration. The stream + // name is assumed to be the same as the constructed group name above. + + IOStream::changeFilename(MeshGroupName, MeshFileName); + + // The connectivity arrays are computed in decomp and here contain only local + // connectivity and local addresses that are dependent on that decomposition. + // They should not be read or written and are therefore not included in + // mesh field groups. + + // CellID; < global cell ID for each local cell + // CellsOnCell, CellsOnCellH; < Indx of cells that neighbor each cell + // EdgesOnCell, EdgesOnCellH; < Indx of edges that border each cell + // NEdgesOnCell, NEdgesOnCellH; < Num of active edges around each cell + // VerticesOnCell, VerticesOnCellH; < Indx of vertices bordering each cell + // CellsOnEdge, CellsOnEdgeH; < Indx of cells straddling each edge + // EdgesOnEdge, EdgesOnEdgeH; < Indx of edges around cells across edge + // NEdgesOnEdge, NEdgesOnEdgeH; < Num of edges around cells across edge + // VerticesOnEdge, VerticesOnEdgeH; < Indx of vertices straddling each edge + // CellsOnVertex, CellsOnVertexH; < Indx of cells that share a vertex + // EdgesOnVertex, EdgesOnVertexH; < Indx of edges sharing vertex as endpt + + // These fields are all read from the mesh file so are added to both the + // full mesh group and the input mesh group + + // Coordinate arrays + // Cell coords + std::string FieldName = "XCell"; + XCell = Array1DReal("XCell", NCellsSize); // allocate space and init to zero + int NDims = 1; + std::vector DimNames(NDims); + DimNames[0] = "NCells"; + auto XCellField = + Field::create(FieldName, // field name + "X Coordinates of cell centers (m)", // long Name + "m", // units + "x", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, XCell); + + FieldName = "YCell"; + YCell = Array1DReal("YCell", NCellsSize); // allocate space and init to zero + auto YCellField = + Field::create(FieldName, // field name + "Y Coordinates of cell centers (m)", // long Name + "m", // units + "y", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, YCell); + + FieldName = "ZCell"; + ZCell = Array1DReal("ZCell", NCellsSize); // allocate space and init to zero + auto ZCellField = + Field::create(FieldName, // field name + "Z Coordinates of cell centers (m)", // long Name + "m", // units + "z", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, ZCell); + + FieldName = "LatCell"; + LatCell = Array1DReal("LatCell", NCellsSize); // allocate space init to zero + auto LatCellField = + Field::create(FieldName, // field name + "Latitude coordinates of cell centers", // long Name + "radians", // units + "latitude", // CF standard Name + -3.1415927, // min valid value + 3.1415927, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, LatCell); + + FieldName = "LonCell"; + LonCell = Array1DReal("LonCell", NCellsSize); // allocate space init to zero + auto LonCellField = + Field::create(FieldName, // field name + "Longitude coordinates of cell centers", // long Name + "radians", // units + "longitude", // CF standard Name + 0.0, // min valid value + 6.28319, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, LonCell); + + // Edge coords + DimNames[0] = "NEdges"; + FieldName = "XEdge"; + XEdge = Array1DReal("XEdge", NEdgesSize); // allocate space init to zero + auto XEdgeField = + Field::create(FieldName, // field name + "X Coordinates of cell edges (m)", // long Name + "m", // units + "x", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, XEdge); + + FieldName = "YEdge"; + YEdge = Array1DReal("YEdge", NEdgesSize); // allocate space init to zero + auto YEdgeField = + Field::create(FieldName, // field name + "Y Coordinates of cell edges (m)", // long Name + "m", // units + "y", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, YEdge); + + FieldName = "ZEdge"; + ZEdge = Array1DReal("ZEdge", NEdgesSize); // allocate space init to zero + auto ZEdgeField = + Field::create(FieldName, // field name + "Z Coordinates of cell edges (m)", // long Name + "m", // units + "z", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, ZEdge); + + FieldName = "LatEdge"; + LatEdge = Array1DReal("LatEdge", NEdgesSize); // allocate space init to zero + auto LatEdgeField = + Field::create(FieldName, // field name + "Latitude coordinates of cell edges", // long Name + "radians", // units + "latitude", // CF standard Name + -3.1415927, // min valid value + 3.1415927, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, LatEdge); + + FieldName = "LonEdge"; + LonEdge = Array1DReal("LonEdge", NEdgesSize); // allocate space init to zero + auto LonEdgeField = + Field::create(FieldName, // field name + "Longitude coordinates of cell edges", // long Name + "radians", // units + "longitude", // CF standard Name + 0.0, // min valid value + 6.28319, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, LonEdge); + + // Vertex coordinates + DimNames[0] = "NVertices"; + FieldName = "XVertex"; + XVertex = Array1DReal("XVertex", NVerticesSize); // allocate space + auto XVertexField = + Field::create(FieldName, // field name + "X Coordinates of cell vertices (m)", // long Name + "m", // units + "x", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, XVertex); + + FieldName = "YVertex"; + YVertex = Array1DReal("YVertex", NVerticesSize); // allocate space + auto YVertexField = + Field::create(FieldName, // field name + "Y Coordinates of cell vertices (m)", // long Name + "m", // units + "y", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, YVertex); + + FieldName = "ZVertex"; + ZVertex = Array1DReal("ZVertex", NVerticesSize); // allocate space + auto ZVertexField = + Field::create(FieldName, // field name + "Z Coordinates of cell vertices (m)", // long Name + "m", // units + "z", // CF standard Name + 0.0, // min valid value + 7.0E+6, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, ZVertex); + + FieldName = "LatVertex"; + LatVertex = Array1DReal("LatVertex", NVerticesSize); // allocate space + auto LatVertexField = + Field::create(FieldName, // field name + "Latitude coordinates of cell vertices", // long Name + "radians", // units + "latitude", // CF standard Name + -3.1415927, // min valid value + 3.1415927, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, LatVertex); + + FieldName = "LonVertex"; + LonVertex = Array1DReal("LonVertex", NVerticesSize); // allocate space + auto LonVertexField = + Field::create(FieldName, // field name + "Longitude coordinates of cell vertices", // long Name + "radians", // units + "longitude", // CF standard Name + 0.0, // min valid value + 6.28319, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, LonVertex); + + // Other mesh properties + // Mesh areas, lengths, and angles + DimNames[0] = "NCells"; + FieldName = "AreaCell"; + AreaCell = Array1DReal("AreaCell", NCellsSize); // allocate space + auto AreaCellField = Field::create(FieldName, // field name + "Area of each cell (m^2)", // long name + "m2", // units + "cell_area", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, AreaCell); + + FieldName = "MeshDensity"; + MeshDensity = Array1DReal("MeshDensity", NCellsSize); // allocate space + auto MeshDensityField = + Field::create(FieldName, // field name + "Value of density function used to generate mesh cell at" + " cell centers", // long name + "", // units + "", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, MeshDensity); + + DimNames[0] = "NVertices"; + FieldName = "AreaTriangle"; + AreaTriangle = Array1DReal("AreaTriangle", NVerticesSize); // allocate space + auto AreaTriangleField = + Field::create(FieldName, // field name + "Area of each triangle in the dual grid (m^2)", // lng name + "m2", // units + "cell_area", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, AreaTriangle); + + DimNames[0] = "NEdges"; + FieldName = "DvEdge"; + DvEdge = Array1DReal("DvEdge", NEdgesSize); // allocate space + auto DvEdgeField = + Field::create(FieldName, // field name + "Length of each edge, computed as the distance between" + " verticesOnEdge (m)", // long name + "m", // units + "", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, DvEdge); + + FieldName = "DcEdge"; + DcEdge = Array1DReal("DcEdge", NEdgesSize); // allocate space + auto DcEdgeField = + Field::create(FieldName, // field name + "Length of each edge, computed as the distance between" + " CellsOnEdge (m)", // long name + "m", // units + "", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, DcEdge); + + FieldName = "AngleEdge"; + AngleEdge = Array1DReal("AngleEdge", NEdgesSize); // allocate space + auto AngleEdgeField = + Field::create(FieldName, // field name + "Angle the edge normal makes with local eastward direction" + " (radians)", // long name + "radians", // units + "", // CF standard name + -3.1415927, // min valid value + 3.1415927, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, AngleEdge); + + NDims = 2; + DimNames.resize(2); + DimNames[0] = "NVertices"; + DimNames[1] = "VertexDegree"; + FieldName = "KiteAreasOnVertex"; + KiteAreasOnVertex = + Array2DReal("KiteAreasOnVertex", NVerticesSize, VertexDegree); + auto KiteAreasOnVertexField = + Field::create(FieldName, // field name + "Area of the portions of each dual cell that are part of " + "each cellsOnVertex (m^2)", // long name + "m2", // units + "", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, KiteAreasOnVertex); + + // Mesh weights + DimNames[0] = "NEdges"; + DimNames[1] = "MaxEdges2"; + FieldName = "WeightsOnEdge"; + WeightsOnEdge = Array2DReal("KiteAreasOnVertex", NEdgesSize, MaxEdges2); + auto WeightsOnEdgeField = + Field::create(FieldName, // field name + "Reconstruction weights associated with each of the" + " edgesOnEdge", // long name + "", // units + "", // CF standard name + -1.0, // min valid value + 1.0, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, WeightsOnEdge); + + // Coriolis parameter at the cells, edges, and vertices + NDims = 1; + DimNames.resize(1); + DimNames[0] = "NEdges"; + FieldName = "FEdge"; + FEdge = Array1DReal("FEdge", NEdgesSize); + auto FEdgeField = + Field::create(FieldName, // field name + "Coriolis parameter at edges (radians s^-1)", // long name + "radians s-1", // units + "coriolis_parameter", // CF standard name + -1.0E-3, // min valid value + 1.0E-3, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, FEdge); + + DimNames[0] = "NVertices"; + FieldName = "FVertex"; + FVertex = Array1DReal("FVertex", NVerticesSize); + auto FVertexField = + Field::create(FieldName, // field name + "Coriolis parameter at vertices (radians s^-1)", + "radians s-1", // units + "coriolis_parameter", // CF standard name + -1.0E-3, // min valid value + 1.0E-3, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, FVertex); + + DimNames[0] = "NCells"; + FieldName = "FCell"; + FCell = Array1DReal("FCell", NCellsSize); + auto FCellField = + Field::create(FieldName, // field name + "Coriolis parameter at cell centers (radians s^-1)", + "radians s-1", // units + "coriolis_parameter", // CF standard name + -1.0E-3, // min valid value + 1.0E-3, // max valid value + -9.99E+30, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); + MeshGroupIn->addField(FieldName); + Field::attachFieldData(FieldName, FCell); + + // The sign and scaling fields are computed internally so are allocated + // here. Because they are not read from a file or typically written to a + // file, we do not create field metadata, though the definitions are + // included here in comments in case they are needed later. + + EdgeSignOnCell = Array2DReal("EdgeSignOnCell", NCellsSize, MaxEdges); + EdgeSignOnVertex = + Array2DReal("EdgeSignOnVertex", NVerticesSize, VertexDegree); + MeshScalingDel2 = Array1DReal("MeshScalingDel2", NEdgesSize); + MeshScalingDel4 = Array1DReal("MeshScalingDel4", NEdgesSize); + + // Sign vectors + // + // NDims = 2; + // DimNames.resize(2); + // DimNames[0] = "NCells"; + // DimNames[1] = "MaxEdges"; + // FieldName = "EdgeSignOnCell"; + // auto EdgeSignOnCellField = + // Field::create(FieldName, // field name + // "Sign of vector connecting cells", // long name + // "", // units + // "", // CF standard name + // -1.0, // min valid value + // 1.0, // max valid value + // -9.99E+30, // scalar for undefined entries + // NDims, // num of dimensions + // DimNames, // dimension names + // false // not time dependent + // ); + // MeshGroupIn->addField(FieldName); + // Field::attachFieldData(FieldName, EdgeSignOnCell); + // + // DimNames[0] = "NVertices"; + // DimNames[1] = "VertexDegree"; + // FieldName = "EdgeSignOnVertex"; + // auto EdgeSignOnVertexField = + // Field::create(FieldName, // field name + // "Sign of vector connecting vertices", // long name + // "", // units + // "", // CF standard name + // -1.0, // min valid value + // 1.0, // max valid value + // -9.99E+30, // scalar for undefined entries + // NDims, // num of dimensions + // DimNames, // dimension names + // false // not time dependent + // ); + // MeshGroupIn->addField(FieldName); + // Field::attachFieldData(FieldName, EdgeSignOnVertex); + // + // Mesh scaling + // + // NDims = 1; + // DimNames.resize(1); + // DimNames[0] = "NEdges"; + // FieldName = "MeshScalingDel2"; + // auto MeshScalingDel2Field = + // Field::create(FieldName, // field name + // "Coefficient to Laplacian mixing terms", // long name + // "", // units + // "", // CF standard name + // -9.99E+30, // min valid value + // 9.99E+30, // max valid value + // -9.99E+30, // scalar for undefined entries + // NDims, // num of dimensions + // DimNames, // dimension names + // false // not time dependent + // ); + // MeshGroupIn->addField(FieldName); + // Field::attachFieldData(FieldName, MeshScalingDel2); + // + // FieldName = "MeshScalingDel4"; + // auto MeshScalingDel4Field = + // Field::create(FieldName, // field name + // "Coefficient to biharmonic mixing terms", // long name + // "", // units + // "", // CF standard name + // -9.99E+30, // min valid value + // 9.99E+30, // max valid value + // -9.99E+30, // scalar for undefined entries + // NDims, // num of dimensions + // DimNames, // dimension names + // false // not time dependent + // ); + // MeshGroupIn->addField(FieldName); + // Field::attachFieldData(FieldName, MeshScalingDel4); + +} // end defineMeshFields + +//------------------------------------------------------------------------------ + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index e36762730e57..cc5817cb98e7 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -14,6 +14,7 @@ #include "Decomp.h" #include "MachEnv.h" #include "OmegaKokkos.h" +#include "TimeMgr.h" #include #include @@ -29,41 +30,30 @@ namespace OMEGA { class HorzMesh { private: - void initParallelIO(Decomp *MeshDecomp); - - void finalizeParallelIO(); - - 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 readMeasurements(); - - void readWeights(); - - void readCoriolis(); + /// Creates dimensions associated with mesh if not already computed + void createDimensions(Decomp *MeshDecomp ///< [in] decomposition for mesh + ); - void copyToDevice(); + /// Creates metadata and allocates data for mesh variables not already + /// computed by Decomp + void defineMeshFields(); - // int computeMesh(); - I4 CellDecompR8; - I4 EdgeDecompR8; - I4 VertexDecompR8; - I4 OnEdgeDecompR8; - I4 OnVertexDecompR8; + /// Finishes filling arrays read from input file by filling halo regions + /// and creating host mirror copies + void completeReadArrays(); + /// The horizontal mesh defined for the default decomposition. Used to + /// provide a shortcut to this common mesh static HorzMesh *DefaultHorzMesh; + /// Map containing pointers to all defined meshes so that they can be + /// retrieved by name static std::map> AllHorzMeshes; /// Construct a new local mesh for a given decomposition HorzMesh(const std::string &Name, ///< [in] Name for mesh - Decomp *Decomp ///< [in] Decomposition for mesh + Decomp *Decomp, ///< [in] Decomposition for mesh + const Clock *ModelClock ///< [in] Model clock for IO ); // Forbid copy and move construction @@ -73,28 +63,31 @@ class HorzMesh { public: // KOKKOS_LAMBDA does not allow to have parallel_* functions inside of a // private function. + + /// Computes the sign of the edge normal at each edge void computeEdgeSign(); - void setMeshScaling(); + /// Computes the mesh scaling coefficients for any Del2 or Del4 mixing + /// operator + void computeMeshScaling(); // Variables // Since these are used frequently, we make them public to reduce the // number of retrievals required. - std::string MeshName; - std::string MeshFileName; - int MeshFileID; + std::string MeshName; ///< name of this mesh + std::string MeshFileName; ///< name (and full path) of input mesh file // Sizes and global IDs // Note that all sizes are actual counts (1-based) so that loop extents // should always use the 0:NCellsXX-1 form. - Array1DI4 NCellsHalo; ///< num cells owned+halo for halo layer - HostArray1DI4 NCellsHaloH; ///< num cells owned+halo for halo layer + Array1DI4 NCellsHalo; ///< num cells owned+halo for each halo layer + HostArray1DI4 NCellsHaloH; ///< num cells owned+halo for each halo layer I4 NCellsOwned; ///< Number of cells owned by this task - I4 NCellsAll; ///< Total number of local cells (owned + all halo) - I4 NCellsSize; ///< Array size (incl padding, bndy cell) for cell arrays - I4 NCellsGlobal; + I4 NCellsAll; ///< Total number of local cells (owned + all halo) + I4 NCellsSize; ///< Array size (incl padding, bndy cell) for cell arrays + I4 NCellsGlobal; ///< Total number of cells in global non-decomposed mesh Array1DI4 NEdgesHalo; ///< num cells owned+halo for halo layer HostArray1DI4 NEdgesHaloH; ///< num cells owned+halo for halo layer @@ -103,8 +96,8 @@ class HorzMesh { I4 NEdgesSize; ///< Array length (incl padding, bndy) for edge dim I4 MaxCellsOnEdge; ///< Max number of cells sharing an edge I4 MaxEdges; ///< Max number of edges around a cell - I4 NEdgesGlobal; - I4 MaxEdges2; ///< Max number of edges around a cell x2 + I4 MaxEdges2; ///< Max number of edges around a cell x2 + I4 NEdgesGlobal; ///< Total number of edges in global non-decomposed mesh Array1DI4 NVerticesHalo; ///< num cells owned+halo for halo layer HostArray1DI4 NVerticesHaloH; ///< num cells owned+halo for halo layer @@ -148,26 +141,36 @@ class HorzMesh { Array1DI4 CellID; ///< global cell ID for each local cell // Coordinates - Array1DReal XCell; ///< X Coordinates of cell centers (m) - HostArray1DReal XCellH; ///< X Coordinates of cell centers (m) - - Array1DReal YCell; ///< Y Coordinates of cell centers (m) - HostArray1DReal YCellH; ///< Y Coordinates of cell centers (m) + Array1DReal XCell; ///< X Coordinates of cell centers (m) + Array1DReal YCell; ///< Y Coordinates of cell centers (m) + Array1DReal ZCell; ///< Z Coordinates of cell centers (m) + Array1DReal LonCell; ///< Longitude location of cell centers (radians) + Array1DReal LatCell; ///< 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) - Array1DReal XEdge; ///< X Coordinate of edge midpoints (m) - HostArray1DReal XEdgeH; ///< X Coordinate of edge midpoints (m) - - Array1DReal YEdge; ///< Y Coordinate of edge midpoints (m) - HostArray1DReal YEdgeH; ///< Y Coordinate of edge midpoints (m) + Array1DReal XEdge; ///< X Coordinate of edge midpoints (m) + Array1DReal YEdge; ///< Y Coordinate of edge midpoints (m) + Array1DReal ZEdge; ///< Z Coordinate of edge midpoints (m) + Array1DReal LonEdge; ///< Longitude location of edge midpoints (radians) + Array1DReal LatEdge; ///< 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) + Array1DReal XVertex; ///< X Coordinate of vertices (m) + Array1DReal YVertex; ///< Y Coordinate of vertices (m) + Array1DReal ZVertex; ///< Z Coordinate of vertices (m) + Array1DReal LonVertex; ///< Longitude location of vertices (radians) + Array1DReal LatVertex; ///< 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) @@ -190,15 +193,13 @@ class HorzMesh { KiteAreasOnVertexH; ///< Area of the portions of each dual cell that /// are part of each cellsOnVertex (m^2) - Array1DReal - DvEdge; ///< 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) - Array1DReal - DcEdge; ///< 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) @@ -207,15 +208,15 @@ class HorzMesh { HostArray1DReal AngleEdgeH; ///< Angle the edge normal makes with local /// eastward direction (radians) - HostArray1DReal - MeshDensityH; ///< Value of density function used to generate a - /// particular mesh at cell centers + Array1DReal MeshDensity; ///< Value of density function used to generate + ///< mesh at cell centers + HostArray1DReal MeshDensityH; ///< Value of density function used to generate + ///< mesh at cell centers // Weights - Array2DReal - WeightsOnEdge; ///< 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 @@ -225,8 +226,7 @@ class HorzMesh { HostArray1DReal FEdgeH; ///< Coriolis parameter at edges (radians s^-1) Array1DReal FCell; ///< Coriolis parameter at cell centers (radians s^-1) - HostArray1DReal - FCellH; ///< Coriolis parameter at cell centers (radians s^-1) + HostArray1DReal FCellH; ///< Coriolis parameter cell centers (radians s^-1) Array1DReal FVertex; ///< Coriolis parameter at vertices (radians s^-1) HostArray1DReal FVertexH; ///< Coriolis parameter at vertices (radians s^-1) @@ -248,13 +248,15 @@ class HorzMesh { // Methods - /// Initialize Omega local mesh - static void init(); + /// Initialize Omega local mesh and create default mesh + static void init(const Clock *ModelClock ///< [in] Model clock for stream IO + ); /// Creates a new mesh by calling the constructor and puts it in the /// AllHorzMeshes map static HorzMesh *create(const std::string &Name, ///< [in] Name for mesh - Decomp *Decomp ///< [in] Decomposition for mesh + Decomp *Decomp, ///< [in] Mesh Decomposition + const Clock *ModelClock ///< [in] Clock for mesh IO ); /// Destructor - deallocates all memory and deletes a HorzMesh @@ -267,9 +269,12 @@ class HorzMesh { static void erase(std::string InName ///< [in] name of mesh to remove ); + /// Retrieves pointer to default mesh static HorzMesh *getDefault(); - static HorzMesh *get(std::string name); + /// Retrieves pointer to a mesh given the mesh name + static HorzMesh *get(std::string Name ///< [in] name of mesh to retrieve + ); }; // end class HorzMesh diff --git a/components/omega/src/ocn/HorzOperators.cpp b/components/omega/src/ocn/HorzOperators.cpp index e182dd4f77ed..f4832b9b4bb5 100644 --- a/components/omega/src/ocn/HorzOperators.cpp +++ b/components/omega/src/ocn/HorzOperators.cpp @@ -31,16 +31,11 @@ SecondDerivativeOnCell::SecondDerivativeOnCell(HorzMesh const *Mesh) : OnSphere(true), NCellsAll(Mesh->NCellsAll), MaxEdges(1 + Mesh->MaxEdges), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), CellsOnCell(Mesh->CellsOnCell), CellsOnEdge(Mesh->CellsOnEdge), - VerticesOnEdge(Mesh->VerticesOnEdge), - - XCell(Mesh->XCell), YCell(Mesh->YCell), - ZCell(createDeviceMirrorCopy(Mesh->ZCellH)), DvEdge(Mesh->DvEdge), + VerticesOnEdge(Mesh->VerticesOnEdge), XCell(Mesh->XCell), + YCell(Mesh->YCell), ZCell(Mesh->ZCell), DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AngleEdge(Mesh->AngleEdge), AreaCell(Mesh->AreaCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), - XVertex(createDeviceMirrorCopy(Mesh->XVertexH)), - YVertex(createDeviceMirrorCopy(Mesh->YVertexH)), - ZVertex(createDeviceMirrorCopy(Mesh->ZVertexH)), - + XVertex(Mesh->XVertex), YVertex(Mesh->YVertex), ZVertex(Mesh->ZVertex), ThetaAbs("SphereAngle", NCellsAll), XPCell("XP", NCellsAll, MaxEdges), YPCell("YP", NCellsAll, MaxEdges), Angle2DCell("Angle2D", NCellsAll, MaxEdges), diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index a22e3b3b0667..12fa3448b19a 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -131,7 +131,7 @@ int initOmegaModules(MPI_Comm Comm) { ABORT_ERROR("ocnInit: Error initializing default halo"); } - HorzMesh::init(); + HorzMesh::init(ModelClock); VertCoord::init(); Tracers::init(); VertAdv::init(); diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index 451ca4293edf..435e8674f768 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -94,7 +94,7 @@ void initIOStreamTest(Clock *&ModelClock // Model clock IOStream::init(ModelClock); // Initialize HorzMesh - this should read Mesh stream - HorzMesh::init(); + HorzMesh::init(ModelClock); HorzMesh *DefMesh = HorzMesh::getDefault(); // Initialize the vertical coordinate diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index a646f31dfce4..4a74669a3150 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -99,13 +99,17 @@ int initAuxStateTest(const std::string &mesh) { Config("Omega"); Config::readAll("omega.yml"); + // Initialize time stepper and get model clock TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); IO::init(DefComm); Decomp::init(mesh); // Initialize streams - IOStream::init(); + Field::init(ModelClock); + IOStream::init(ModelClock); int HaloErr = Halo::init(); if (HaloErr != 0) { @@ -113,7 +117,7 @@ int initAuxStateTest(const std::string &mesh) { LOG_ERROR("AuxStateTest: error initializing default halo"); } - HorzMesh::init(); + HorzMesh::init(ModelClock); VertCoord::init(); diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 44957762f953..28095626e595 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -7,6 +7,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OceanTestCommon.h" @@ -802,7 +803,17 @@ int initAuxVarsTest(const std::string &mesh) { LOG_ERROR("AuxVarsTest: error initializing default halo"); } - HorzMesh::init(); + // Create dummy clock for IO + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; + + // Read horizontal mesh + Field::init(ModelClock); + IOStream::init(ModelClock); + HorzMesh::init(ModelClock); // initialize vertical coordinate, do not read stream and use local // NVertLayers value diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index d5655033d82a..8b9c5f732ce7 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -14,7 +14,9 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Field.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OceanTestCommon.h" @@ -86,8 +88,17 @@ void initEosTest(const std::string &mesh) { /// Initialize Halo Halo::init(); - /// Initialize mesh - HorzMesh::init(); + /// Create dummy model clock + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; + + /// Read horizontal mesh + Field::init(ModelClock); + IOStream::init(ModelClock); + HorzMesh::init(ModelClock); /// Initialize vertical coordinate VertCoord::init(false); diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 0a0f96b710f7..a1849fb0cebd 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -14,25 +14,24 @@ #include "Decomp.h" #include "Dimension.h" #include "Error.h" +#include "Field.h" #include "Halo.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OmegaKokkos.h" #include "Pacer.h" +#include "TimeMgr.h" #include "mpi.h" -#include - using namespace OMEGA; //------------------------------------------------------------------------------ // The initialization routine for Mesh testing. It calls various // init routines, including the creation of the default decomposition. -int initHorzMeshTest() { - - int Err = 0; +void initHorzMeshTest() { // Initialize the Machine Environment class - this also creates // the default MachEnv. Then retrieve the default environment and @@ -43,6 +42,7 @@ int initHorzMeshTest() { // Initialize the Logging system initLogging(DefEnv); + LOG_INFO("------ Horz Mesh Unit Tests ------"); // Open config file Config("Omega"); @@ -55,78 +55,72 @@ int initHorzMeshTest() { Decomp::init(); // Initialize the default halo - Err = Halo::init(); + int Err = Halo::init(); if (Err != 0) - LOG_ERROR("HorzMeshTest: error initializing default halo"); + ABORT_ERROR("HorzMeshTest: error initializing default halo"); + + // Create a dummy model clock + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; // Initialize the default mesh - HorzMesh::init(); + Field::init(ModelClock); + IOStream::init(ModelClock); + HorzMesh::init(ModelClock); - return Err; + return; } //------------------------------------------------------------------------------ // Computes the distance of a x,y,z coordinate from the origin -R8 distance(R8 x, R8 y, R8 z) { +R8 distance(R8 X, R8 Y, R8 Z) { - R8 dist; - - dist = sqrt(x * x + y * y + z * z); - - return dist; + R8 Dist = sqrt(X * X + Y * Y + Z * Z); + return Dist; } //------------------------------------------------------------------------------ // Computes the distance between two lon/lat points on the sphere -R8 sphereDistance(R8 lon1, R8 lat1, R8 lon2, R8 lat2) { - - R8 arg; +R8 sphereDistance(R8 Lon1, R8 Lat1, R8 Lon2, R8 Lat2) { - arg = sqrt(pow(sin(0.5 * (lat1 - lat2)), 2) + - cos(lat2) * cos(lat1) * pow(sin(0.5 * (lon1 - lon2)), 2)); - return 2.0 * asin(arg); + R8 Arg = sqrt(pow(sin(0.5 * (Lat1 - Lat2)), 2) + + cos(Lat2) * cos(Lat1) * pow(sin(0.5 * (Lon1 - Lon2)), 2)); + return 2.0 * asin(Arg); } //------------------------------------------------------------------------------ // Computes the longitude of a point given its Cartesian coordinates -R8 computeLon(R8 x, R8 y, R8 z) { +R8 computeLon(R8 X, R8 Y, R8 Z) { - R8 lon; - lon = atan2(y, x); - - R8 pi; - pi = 4.0 * atan(1.0); - - if (lon < 0.0) { - lon = 2.0 * pi + lon; + R8 Lon = atan2(Y, X); + R8 Pi = 4.0 * atan(1.0); + if (Lon < 0.0) { + Lon = 2.0 * Pi + Lon; } - - return lon; + return Lon; } //------------------------------------------------------------------------------ // Computes the latitude of a point given its Cartesian coordinates -R8 computeLat(R8 x, R8 y, R8 z) { +R8 computeLat(R8 X, R8 Y, R8 Z) { - R8 dist; - dist = distance(x, y, z); + R8 Dist = distance(X, Y, Z); + R8 Lat = asin(Z / Dist); - R8 lat; - lat = asin(z / dist); - - return lat; + return Lat; } //------------------------------------------------------------------------------ // Computes coriolis parameter for a given latitude -R8 coriolis(R8 lat) { - - R8 f; - R8 omega = 7.29212e-5; +R8 coriolis(R8 Lat) { - f = 2.0 * omega * sin(lat); + R8 Omega = 7.29212e-5; + R8 F = 2.0 * Omega * sin(Lat); - return f; + return F; } //------------------------------------------------------------------------------ @@ -135,8 +129,6 @@ R8 coriolis(R8 lat) { // int main(int argc, char *argv[]) { - int RetVal = 0; - // Initialize the global MPI environment MPI_Init(&argc, &argv); Kokkos::initialize(); @@ -144,13 +136,12 @@ int main(int argc, char *argv[]) { Pacer::setPrefix("Omega:"); { - R8 tol = 1e-6; - R8 pi = 4.0 * atan(1.0); + int Err = 0; + R8 Tol = 1e-6; + R8 Pi = 4.0 * atan(1.0); // Call initialization routine to create the default decomposition - int Err = initHorzMeshTest(); - if (Err != 0) - LOG_CRITICAL("HorzMeshTest: Error initializing"); + initHorzMeshTest(); // Get MPI vars if needed MachEnv *DefEnv = MachEnv::getDefault(); @@ -161,12 +152,8 @@ int main(int argc, char *argv[]) { // Test retrieval of the default decomposition Decomp *DefDecomp = Decomp::getDefault(); - if (DefDecomp) { // true if non-null ptr - LOG_INFO("HorzMeshTest: Default decomp retrieval PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Default decomp retrieval FAIL"); - } + if (DefDecomp == nullptr) + ABORT_ERROR("HorzMeshTest: Default decomp retrieval FAIL"); // Retrieve default mesh HorzMesh *Mesh = HorzMesh::getDefault(); @@ -179,64 +166,54 @@ int main(int argc, char *argv[]) { I4 LocCells; LocCells = Mesh->NCellsOwned; Err = MPI_Allreduce(&LocCells, &SumCells, 1, MPI_INT32_T, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("HorzMeshTest: MPI error summing cells"); - if (SumCells == Mesh->NCellsGlobal) { - LOG_INFO("HorzMeshTest: Sum cell ID test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Sum cell ID test FAIL {} {}", SumCells, - Mesh->NCellsGlobal); - } + if (SumCells != Mesh->NCellsGlobal) + ABORT_ERROR("HorzMeshTest: Sum cell ID test FAIL {} {}", SumCells, + Mesh->NCellsGlobal); // Test that cell centers are on sphere // Check that all cell centers are a uniform distance from the origin // Tests that the Cartesian coordinates for cell centers have been read in // corectly - R8 sphere_radius = + R8 SphereRadius = distance(Mesh->XCellH(0), Mesh->YCellH(0), Mesh->ZCellH(0)); - R8 dist; - I4 count = 0; + R8 Dist; + I4 Count = 0; for (int Cell = 0; Cell < LocCells; Cell++) { - dist = distance(Mesh->XCellH(Cell), Mesh->YCellH(Cell), + Dist = distance(Mesh->XCellH(Cell), Mesh->YCellH(Cell), Mesh->ZCellH(Cell)); - if (abs(sphere_radius - dist) > tol) - count++; + if (abs(SphereRadius - Dist) > Tol) + Count++; } - if (count > 0) { - RetVal += 1; - LOG_INFO("HorzMeshTest: Cell sphere radius test FAIL"); - } else { - LOG_INFO("HorzMeshTest: Cell sphere radius test PASS"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: Cell sphere radius test FAIL"); // Test lon/lat coordinates of cell centers // Convert Cartesian coordinates to lon/lat and check these agree with the // values that have been read in // Tests that the lon/lat coordinates for cell // centers have been read in correctly - R8 lon; - R8 lat; - count = 0; + R8 Lon; + R8 Lat; + Count = 0; for (int Cell = 0; Cell < LocCells; Cell++) { - lon = computeLon(Mesh->XCellH(Cell), Mesh->YCellH(Cell), + Lon = computeLon(Mesh->XCellH(Cell), Mesh->YCellH(Cell), Mesh->ZCellH(Cell)); - lat = computeLat(Mesh->XCellH(Cell), Mesh->YCellH(Cell), + Lat = computeLat(Mesh->XCellH(Cell), Mesh->YCellH(Cell), Mesh->ZCellH(Cell)); - if (abs(lon - Mesh->LonCellH(Cell)) > tol) - count++; - if (abs(lat - Mesh->LatCellH(Cell)) > tol) - count++; + if (abs(Lon - Mesh->LonCellH(Cell)) > Tol) + Count++; + if (abs(Lat - Mesh->LatCellH(Cell)) > Tol) + Count++; } - if (count > 0) { - RetVal += 1; - LOG_INFO("HorzMeshTest: Cell lon/lat test FAIL"); - } else { - LOG_INFO("HorzMeshTest: Cell lon/lat test PASS"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: Cell lon/lat test FAIL"); // Test sum of local mesh edges // Get the global sum of all local edge counts @@ -247,60 +224,48 @@ int main(int argc, char *argv[]) { LocEdges = Mesh->NEdgesOwned; Err = MPI_Allreduce(&LocEdges, &SumEdges, 1, MPI_INT32_T, MPI_SUM, Comm); - if (SumEdges == Mesh->NEdgesGlobal) { - LOG_INFO("HorzMeshTest: Sum edge ID test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Sum edge ID test FAIL {} {}", SumEdges, - Mesh->NEdgesGlobal); - } + if (SumEdges != Mesh->NEdgesGlobal) + ABORT_ERROR("HorzMeshTest: Sum edge ID test FAIL {} {}", SumEdges, + Mesh->NEdgesGlobal); // Test that edge coordinates are on sphere // Check that all edge centers are a uniform distance from the origin // Tests that the Cartesian coordinates for edge centers have been read in // correctly - sphere_radius = + SphereRadius = distance(Mesh->XEdgeH(0), Mesh->YEdgeH(0), Mesh->ZEdgeH(0)); - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { - dist = distance(Mesh->XEdgeH(Edge), Mesh->YEdgeH(Edge), + Dist = distance(Mesh->XEdgeH(Edge), Mesh->YEdgeH(Edge), Mesh->ZEdgeH(Edge)); - if (abs(sphere_radius - dist) > tol) - count++; + if (abs(SphereRadius - Dist) > Tol) + Count++; } - if (count > 0) { - RetVal += 1; - LOG_INFO("HorzMeshTest: Edge sphere radius test FAIL"); - } else { - LOG_INFO("HorzMeshTest: Edge sphere radius test PASS"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: Edge sphere radius test FAIL"); // Test lon/lat coordinates of edge centers // Convert Cartesian coordinates to lon/lat and check these agree with the // values that have been read in // Tests that the lon/lat coordinates for edge centers have been read in // correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { - lon = computeLon(Mesh->XEdgeH(Edge), Mesh->YEdgeH(Edge), + Lon = computeLon(Mesh->XEdgeH(Edge), Mesh->YEdgeH(Edge), Mesh->ZEdgeH(Edge)); - lat = computeLat(Mesh->XEdgeH(Edge), Mesh->YEdgeH(Edge), + Lat = computeLat(Mesh->XEdgeH(Edge), Mesh->YEdgeH(Edge), Mesh->ZEdgeH(Edge)); - if (abs(lon - Mesh->LonEdgeH(Edge)) > tol) - count++; - if (abs(lat - Mesh->LatEdgeH(Edge)) > tol) - count++; + if (abs(Lon - Mesh->LonEdgeH(Edge)) > Tol) + Count++; + if (abs(Lat - Mesh->LatEdgeH(Edge)) > Tol) + Count++; } - if (count > 0) { - RetVal += 1; - LOG_INFO("HorzMeshTest: Edge lon/lat test FAIL"); - } else { - LOG_INFO("HorzMeshTest: Edge lon/lat test PASS"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: Edge lon/lat test FAIL"); // Test sum of local mesh vertices // Get the global sum of all local vertex counts @@ -311,62 +276,52 @@ int main(int argc, char *argv[]) { LocVertices = Mesh->NVerticesOwned; Err = MPI_Allreduce(&LocVertices, &SumVertices, 1, MPI_INT32_T, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("HorzMeshTest: MPI error summing vertices"); - if (SumVertices == DefDecomp->NVerticesGlobal) { - LOG_INFO("HorzMeshTest: Sum vertex ID test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Sum vertex ID test FAIL {} {}", SumVertices, - DefDecomp->NVerticesGlobal); - } + if (SumVertices != DefDecomp->NVerticesGlobal) + ABORT_ERROR("HorzMeshTest: Sum vertex ID test FAIL {} {}", SumVertices, + DefDecomp->NVerticesGlobal); // Test that vertex coordinates are on sphere // Check that all vertices are a uniform distance from the origin // Tests that the Cartesian coordinates for vertices have been read in // correctly - sphere_radius = + SphereRadius = distance(Mesh->XVertexH(0), Mesh->YVertexH(0), Mesh->ZVertexH(0)); - count = 0; + Count = 0; for (int Vertex = 0; Vertex < LocVertices; Vertex++) { - dist = distance(Mesh->XVertexH(Vertex), Mesh->YVertexH(Vertex), + Dist = distance(Mesh->XVertexH(Vertex), Mesh->YVertexH(Vertex), Mesh->ZVertexH(Vertex)); - if (abs(sphere_radius - dist) > tol) - count++; + if (abs(SphereRadius - Dist) > Tol) + Count++; } - if (count > 0) { - RetVal += 1; - LOG_INFO("HorzMeshTest: Vertex sphere radius test FAIL"); - } else { - LOG_INFO("HorzMeshTest: Vertex sphere radius test PASS"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: Vertex sphere radius test FAIL"); // Test lon/lat coordinates of vertices // Convert Cartesian coordinates to lon/lat and check these agree with the // values that have been read in // Tests that the lon/lat coordinates for vertices have been read in // correctly - count = 0; + Count = 0; for (int Vertex = 0; Vertex < LocVertices; Vertex++) { - lon = computeLon(Mesh->XVertexH(Vertex), Mesh->YVertexH(Vertex), + Lon = computeLon(Mesh->XVertexH(Vertex), Mesh->YVertexH(Vertex), Mesh->ZVertexH(Vertex)); - lat = computeLat(Mesh->XVertexH(Vertex), Mesh->YVertexH(Vertex), + Lat = computeLat(Mesh->XVertexH(Vertex), Mesh->YVertexH(Vertex), Mesh->ZVertexH(Vertex)); - if (abs(lon - Mesh->LonVertexH(Vertex)) > tol) - count++; + if (abs(Lon - Mesh->LonVertexH(Vertex)) > Tol) + Count++; - if (abs(lat - Mesh->LatVertexH(Vertex)) > tol) - count++; + if (abs(Lat - Mesh->LatVertexH(Vertex)) > Tol) + Count++; } - if (count > 0) { - RetVal += 1; - LOG_INFO("HorzMeshTest: Vertex lon/lat test FAIL"); - } else { - LOG_INFO("HorzMeshTest: Vertex lon/lat test PASS"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: Vertex lon/lat test FAIL"); // Test cell areas // Find the global sum of all the local cell areas @@ -379,14 +334,12 @@ int main(int argc, char *argv[]) { } Err = MPI_Allreduce(&LocSumArea, &SumCellArea, 1, MPI_DOUBLE, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("HorzMeshTest: MPI error summing cell area"); R8 OceanArea = 3.61e14; - if (abs(SumCellArea - OceanArea) / OceanArea < 0.05) { - LOG_INFO("HorzMeshTest: Cell area test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Cell area test FAIL"); - } + if (abs(SumCellArea - OceanArea) / OceanArea > 0.05) + ABORT_ERROR("HorzMeshTest: Cell area test FAIL"); // Test triangle areas // Find the global sum of all the local triangle areas @@ -399,13 +352,11 @@ int main(int argc, char *argv[]) { } Err = MPI_Allreduce(&LocSumArea, &SumTriangleArea, 1, MPI_DOUBLE, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("HorzMeshTest: MPI error summing triangle area"); - if (abs(SumTriangleArea - OceanArea) / OceanArea < 0.05) { - LOG_INFO("HorzMeshTest: Triangle area test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Triangle area test FAIL"); - } + if (abs(SumTriangleArea - OceanArea) / OceanArea > 0.05) + ABORT_ERROR("HorzMeshTest: Triangle area test FAIL"); // Test kite areas // Find the local sum of all the local kite areas @@ -420,48 +371,42 @@ int main(int argc, char *argv[]) { } Err = MPI_Allreduce(&LocSumArea, &SumKiteArea, 1, MPI_DOUBLE, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("HorzMeshTest: MPI error summing kite area"); - if (abs(SumKiteArea - OceanArea) / OceanArea < 0.05) { - LOG_INFO("HorzMeshTest: Kite area test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Kite area test FAIL"); - } + if (abs(SumKiteArea - OceanArea) / OceanArea > 0.05) + ABORT_ERROR("HorzMeshTest: Kite area test FAIL"); - // Test dcEdge + // Test DcEdge // Compute spherical distance between cell centers and compare to value // that was read in Tests that the distances between cell centers have // been read in correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { int Cell1 = Mesh->CellsOnEdgeH(Edge, 0); int Cell2 = Mesh->CellsOnEdgeH(Edge, 1); if ((Cell1 < DefDecomp->NCellsAll) && (Cell2 < DefDecomp->NCellsAll)) { - R8 dc = + R8 Dc = sphereDistance(Mesh->LonCellH(Cell1), Mesh->LatCellH(Cell1), Mesh->LonCellH(Cell2), Mesh->LatCellH(Cell2)); - dc = sphere_radius * dc; + Dc = SphereRadius * Dc; - if (abs((dc - Mesh->DcEdgeH(Edge)) / Mesh->DcEdgeH(Edge)) > tol) { - count++; + if (abs((Dc - Mesh->DcEdgeH(Edge)) / Mesh->DcEdgeH(Edge)) > Tol) { + Count++; } } } - if (count == 0) { - LOG_INFO("HorzMeshTest: dcEdge test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: dcEdge test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: DcEdge test FAIL"); - // Test dvEdge + // Test DvEdge // Compute spherical distance between vertices on edges and compare to // value that was read in Tests that the distances between vertices have // been read in correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { int Vertex1 = Mesh->VerticesOnEdgeH(Edge, 0); int Vertex2 = Mesh->VerticesOnEdgeH(Edge, 1); @@ -469,196 +414,165 @@ int main(int argc, char *argv[]) { if ((Vertex1 < DefDecomp->NVerticesAll) && (Vertex2 < DefDecomp->NVerticesAll)) { - R8 dv = sphereDistance( + R8 Dv = sphereDistance( Mesh->LonVertexH(Vertex1), Mesh->LatVertexH(Vertex1), Mesh->LonVertexH(Vertex2), Mesh->LatVertexH(Vertex2)); - dv = sphere_radius * dv; + Dv = SphereRadius * Dv; - if (abs((dv - Mesh->DvEdgeH(Edge)) / Mesh->DvEdgeH(Edge)) > tol) { - count++; + if (abs((Dv - Mesh->DvEdgeH(Edge)) / Mesh->DvEdgeH(Edge)) > Tol) { + Count++; } } } - if (count == 0) { - LOG_INFO("HorzMeshTest: dvEdge test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: dvEdge test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: DvEdge test FAIL"); - // Test angleEdge - // Check that the range of edge angles is between (-pi, pi) + // Test AngleEdge + // Check that the range of edge angles is between (-Pi, Pi) // Tests that the edge angles have been read in correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { - if (abs(Mesh->AngleEdgeH(Edge)) > pi) { - count++; + if (abs(Mesh->AngleEdgeH(Edge)) > Pi) { + Count++; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: angleEdge test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: angleEdge test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: angleEdge test FAIL"); - // Test fCell + // Test FCell // Compute the Coriolis parameter for cell centers and compare with values // that were read in // Tests that the cell Coriolis values were read in correctly - count = 0; + Count = 0; for (int Cell = 0; Cell < LocCells; Cell++) { - R8 f = coriolis(Mesh->LatCellH(Cell)); + R8 F = coriolis(Mesh->LatCellH(Cell)); - if (abs(f - Mesh->FCellH(Cell)) > tol) { - count++; + if (abs(F - Mesh->FCellH(Cell)) > Tol) { + Count++; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: fCell test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: fCell test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: fCell test FAIL"); - // Test fVertex + // Test FVertex // Compute the Coriolis parameter for vertices and compare with values // that were read in Tests that the vertex Coriolis values were read in // correctly - count = 0; + Count = 0; for (int Vertex = 0; Vertex < LocVertices; Vertex++) { - R8 f = coriolis(Mesh->LatVertexH(Vertex)); + R8 F = coriolis(Mesh->LatVertexH(Vertex)); - if (abs(f - Mesh->FVertexH(Vertex)) > tol) { - count++; + if (abs(F - Mesh->FVertexH(Vertex)) > Tol) { + Count++; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: fVertex test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: fVertex test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: fVertex test FAIL"); - // Test fEdge + // Test FEdge // Compute the Coriolis parameter for edges and compare with values that // were read in // Tests that the edge Coriolis values were read in correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { - R8 f = coriolis(Mesh->LatEdgeH(Edge)); + R8 F = coriolis(Mesh->LatEdgeH(Edge)); - if (abs(f - Mesh->FEdgeH(Edge)) > tol) { - count++; + if (abs(F - Mesh->FEdgeH(Edge)) > Tol) { + Count++; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: fEdge test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: fEdge test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: fEdge test FAIL"); // Test weightsOnEdge // Check the range of the edge weights // Tests that the edge weights were read in correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { for (int i = 0; i < Mesh->MaxEdges2; i++) { if (abs(Mesh->WeightsOnEdgeH(Edge, i)) > 1.0) { - count++; + Count++; } } } - if (count == 0) { - LOG_INFO("HorzMeshTest: weightsOnEdge test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: weightsOnEdge test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: weightsOnEdge test FAIL"); + // Test edgeSignOnCell // Check that the sign corresponds with convention // Tests that the edge sign values were calculated correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { int Cell0 = Mesh->CellsOnEdgeH(Edge, 0); - int iEdge0; + int IEdge0; for (int i = 0; i < Mesh->NEdgesOnCellH(Cell0); i++) { if (Mesh->EdgesOnCellH(Cell0, i) == Edge) { - iEdge0 = i; + IEdge0 = i; break; } } - if (abs(Mesh->EdgeSignOnCellH(Cell0, iEdge0) + 1.0) > tol) { - count++; + if (abs(Mesh->EdgeSignOnCellH(Cell0, IEdge0) + 1.0) > Tol) { + Count++; } int Cell1 = Mesh->CellsOnEdgeH(Edge, 1); if (Cell1 < DefDecomp->NCellsAll) { - int iEdge1; + int IEdge1; for (int i = 0; i < Mesh->NEdgesOnCellH(Cell1); i++) { if (Mesh->EdgesOnCellH(Cell1, i) == Edge) { - iEdge1 = i; + IEdge1 = i; break; } } - if (abs(Mesh->EdgeSignOnCellH(Cell1, iEdge1) - 1.0) > tol) { - count++; + if (abs(Mesh->EdgeSignOnCellH(Cell1, IEdge1) - 1.0) > Tol) { + Count++; } } } - if (count == 0) { - LOG_INFO("HorzMeshTest: edgeSignOnCell test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: edgeSignOnCell test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: edgeSignOnCell test FAIL"); // Test edgeSignOnVertex // Check that the sign corresponds with convention // Tests that the edge sign vlues were calculated correctly - count = 0; + Count = 0; for (int Edge = 0; Edge < LocEdges; Edge++) { int Vertex0 = Mesh->VerticesOnEdgeH(Edge, 0); - int iEdge0; + int IEdge0; for (int i = 0; i < Mesh->VertexDegree; i++) { if (Mesh->EdgesOnVertexH(Vertex0, i) == Edge) { - iEdge0 = i; + IEdge0 = i; break; } } - if (abs(Mesh->EdgeSignOnVertexH(Vertex0, iEdge0) + 1.0) > tol) { - count++; + if (abs(Mesh->EdgeSignOnVertexH(Vertex0, IEdge0) + 1.0) > Tol) { + Count++; } int Vertex1 = Mesh->VerticesOnEdgeH(Edge, 1); - int iEdge1; + int IEdge1; for (int i = 0; i < Mesh->VertexDegree; i++) { if (Mesh->EdgesOnVertexH(Vertex1, i) == Edge) { - iEdge1 = i; + IEdge1 = i; break; } } - if (abs(Mesh->EdgeSignOnVertexH(Vertex1, iEdge1) - 1.0) > tol) { - count++; + if (abs(Mesh->EdgeSignOnVertexH(Vertex1, IEdge1) - 1.0) > Tol) { + Count++; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: edgeSignOnVertex test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: edgeSignOnVertex test FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: edgeSignOnVertex test FAIL"); // Test cell halo values // Perform halo exhange on owned cell only array and compare @@ -666,7 +580,6 @@ int main(int argc, char *argv[]) { // Tests that halo values are read in correctly Halo *DefHalo = Halo::getDefault(); HostArray1DR8 XCellTest("XCellTest", Mesh->NCellsSize); - // Mesh->XCellH.deep_copy_to(XCellTest); deepCopy(XCellTest, Mesh->XCellH); for (int Cell = Mesh->NCellsOwned; Cell < Mesh->NCellsAll; Cell++) { @@ -674,27 +587,22 @@ int main(int argc, char *argv[]) { } DefHalo->exchangeFullArrayHalo(XCellTest, OnCell); - count = 0; + Count = 0; for (int Cell = 0; Cell < Mesh->NCellsAll; Cell++) { if (Mesh->XCellH(Cell) != XCellTest(Cell)) { - count++; + Count++; break; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: cell halo exhange PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: cell halo exhange FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: cell halo exhange FAIL"); // Test edge halo values // Perform halo exhange on owned edge only array and compare // read values // Tests that halo values are read in correctly HostArray1DR8 XEdgeTest("XEdgeTest", Mesh->NEdgesSize); - // Mesh->XEdgeH.deep_copy_to(XEdgeTest); deepCopy(XEdgeTest, Mesh->XEdgeH); for (int Edge = Mesh->NEdgesOwned; Edge < Mesh->NEdgesAll; Edge++) { @@ -702,27 +610,22 @@ int main(int argc, char *argv[]) { } DefHalo->exchangeFullArrayHalo(XEdgeTest, OnEdge); - count = 0; + Count = 0; for (int Edge = 0; Edge < Mesh->NEdgesAll; Edge++) { if (Mesh->XEdgeH(Edge) != XEdgeTest(Edge)) { - count++; + Count++; break; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: edge halo exhange PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: edge halo exhange FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: edge halo exhange FAIL"); // Test vertex halo values // Perform halo exhange on owned vertex only array and compare // read values // Tests that halo values are read in correctly HostArray1DR8 XVertexTest("XVertexTest", Mesh->NVerticesSize); - // Mesh->XVertexH.deep_copy_to(XVertexTest); deepCopy(XVertexTest, Mesh->XVertexH); for (int Vertex = Mesh->NVerticesOwned; Vertex < Mesh->NVerticesAll; @@ -731,20 +634,16 @@ int main(int argc, char *argv[]) { } DefHalo->exchangeFullArrayHalo(XVertexTest, OnVertex); - count = 0; + Count = 0; for (int Vertex = 0; Vertex < Mesh->NVerticesAll; Vertex++) { if (Mesh->XVertexH(Vertex) != XVertexTest(Vertex)) { - count++; + Count++; break; } } - if (count == 0) { - LOG_INFO("HorzMeshTest: vertex halo exhange PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: vertex halo exhange FAIL"); - } + if (Count > 0) + ABORT_ERROR("HorzMeshTest: vertex halo exhange FAIL"); // Test that all Cell IDs are accounted for by // summing the list of owned values by all tasks. The result should @@ -756,33 +655,31 @@ int main(int argc, char *argv[]) { for (int n = 0; n < Mesh->NCellsOwned; ++n) LocSumIDs += CellIDH(n); Err = MPI_Allreduce(&LocSumIDs, &SumIDs, 1, MPI_INT32_T, MPI_SUM, Comm); + if (Err != MPI_SUCCESS) + ABORT_ERROR("HorzMeshTest: MPI error summing cell IDs"); - if (SumIDs == RefSumIDs) { - LOG_INFO("DecompTest: Sum cell ID test PASS"); - } else { - RetVal += 1; - LOG_INFO("DecompTest: Sum cell ID test FAIL {} {}", SumIDs, RefSumIDs); - } + if (SumIDs != RefSumIDs) + ABORT_ERROR("DecompTest: Sum cell ID test FAIL {} {}", SumIDs, + RefSumIDs); // Finalize Omega objects + IOStream::finalize(); HorzMesh::clear(); + Field::clear(); Dimension::clear(); Halo::clear(); Decomp::clear(); MachEnv::removeAll(); - - // MPI_Status status; - if (Err == 0) - LOG_INFO("HorzMeshTest: Successful completion"); } Pacer::finalize(); Kokkos::finalize(); - MPI_Finalize(); - if (RetVal >= 256) - RetVal = 255; + // Success if we made it here + LOG_INFO("------ Horz Mesh unit tests successful ------"); + MPI_Finalize(); - return RetVal; + // Return success if we made it here + return 0; } // end of main //===-----------------------------------------------------------------------===/ diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index e4cfde75ac5e..dbefa58664e5 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -3,10 +3,12 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Field.h" #include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OceanTestCommon.h" @@ -926,6 +928,7 @@ int initOperatorsTest(const std::string &MeshFile) { // Initialize the Logging system initLogging(DefEnv); + LOG_INFO("------ Horz Operators unit tests ------"); // Open config file Config("Omega"); @@ -940,13 +943,25 @@ int initOperatorsTest(const std::string &MeshFile) { LOG_ERROR("OperatorsTest: error initializing default halo"); } - HorzMesh::init(); + // Create dummy model clock + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; + + // Read horizontal mesh + Field::init(ModelClock); + IOStream::init(ModelClock); + HorzMesh::init(ModelClock); return Err; } void finalizeOperatorsTest() { + IOStream::finalize(); HorzMesh::clear(); + Field::clear(); Dimension::clear(); Halo::clear(); Decomp::clear(); @@ -992,6 +1007,8 @@ int main(int argc, char *argv[]) { Pacer::finalize(); Kokkos::finalize(); + if (RetVal == 0) + LOG_INFO("------ Horz Operators unit tests successful ------"); MPI_Finalize(); if (RetVal >= 256) diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp index e93a32ad2030..96a36379804a 100644 --- a/components/omega/test/ocn/StateTest.cpp +++ b/components/omega/test/ocn/StateTest.cpp @@ -65,13 +65,13 @@ void initStateTest() { // Initialize the IO system IO::init(DefComm); + // Initialize Field infrastructure + Field::init(ModelClock); + // Initialize IOStreams - this does not yet validate the contents // of each file, only creates streams from Config IOStream::init(ModelClock); - // Initialize Field infrastructure - Field::init(ModelClock); - // Create the default decomposition (initializes the decomposition) Decomp::init(); @@ -81,7 +81,7 @@ void initStateTest() { ABORT_ERROR("State: error initializing default halo"); // Initialize the default mesh - HorzMesh::init(); + HorzMesh::init(ModelClock); // Initialize the vertical coordinate VertCoord::init(); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index c87828f6bccf..2b2c44438937 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -114,13 +114,17 @@ int initTendenciesTest(const std::string &mesh) { Config("Omega"); Config::readAll("omega.yml"); + // Initialize time stepping and model clock TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); IO::init(DefComm); Decomp::init(mesh); // Initialize streams - IOStream::init(); + Field::init(ModelClock); + IOStream::init(ModelClock); int HaloErr = Halo::init(); if (HaloErr != 0) { @@ -128,8 +132,10 @@ int initTendenciesTest(const std::string &mesh) { LOG_ERROR("TendenciesTest: error initializing default halo"); } - HorzMesh::init(); + // Read mesh + HorzMesh::init(ModelClock); VertCoord::init(); + Tracers::init(); VertAdv::init(); PressureGrad::init(); diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 980a7d759dcb..89ffd9248e0e 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -25,6 +25,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "OceanTestCommon.h" #include "OmegaKokkos.h" @@ -1026,7 +1027,10 @@ void initTendTest(const std::string &MeshFile, int NVertLayers) { Config::Initialize(); Config::readAll("omega.yml"); + // Initialize time stepping and get model clock TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); IO::init(DefComm); @@ -1037,7 +1041,9 @@ void initTendTest(const std::string &MeshFile, int NVertLayers) { ABORT_ERROR("TendencyTermsTest: error initializing default halo"); } - HorzMesh::init(); + Field::init(ModelClock); // Fields are used in streams + IOStream::init(ModelClock); // initialize streams for mesh reading + HorzMesh::init(ModelClock); // initialize vertical coordinate, do not read stream and use local // NVertLayers value diff --git a/components/omega/test/ocn/TracersTest.cpp b/components/omega/test/ocn/TracersTest.cpp index ec5a6889971f..738ee0f55505 100644 --- a/components/omega/test/ocn/TracersTest.cpp +++ b/components/omega/test/ocn/TracersTest.cpp @@ -17,6 +17,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OmegaKokkos.h" @@ -53,8 +54,10 @@ I4 initTracersTest() { Config("Omega"); Config::readAll("omega.yml"); - // Initialize the default time stepper + // Initialize the default time stepper and model clock TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); // Initialize the IO system IO::init(DefComm); @@ -69,8 +72,10 @@ I4 initTracersTest() { return Err; } - // Initialize the default mesh - HorzMesh::init(); + // Read in the default mesh + Field::init(ModelClock); + IOStream::init(ModelClock); + HorzMesh::init(ModelClock); // Initialize the vertical coordinate VertCoord::init(false); diff --git a/components/omega/test/ocn/VertAdvTest.cpp b/components/omega/test/ocn/VertAdvTest.cpp index e60a6586aed0..b1024f55a73c 100644 --- a/components/omega/test/ocn/VertAdvTest.cpp +++ b/components/omega/test/ocn/VertAdvTest.cpp @@ -114,8 +114,10 @@ void initVertAdvTest() { Config("Omega"); Config::readAll("omega.yml"); - // First step of time stepper initialization needed for IOstream + // First step of time stepper initialization needed for streams, model clock TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); // Initialize the IO system IO::init(DefComm); @@ -124,15 +126,16 @@ void initVertAdvTest() { Decomp::init(); // Initialize streams - IOStream::init(); + Field::init(ModelClock); + IOStream::init(ModelClock); // Initialize the default halo Err = Halo::init(); if (Err != 0) ABORT_ERROR("VertAdvTest: error initializing default halo"); - // Initialize the default mesh - HorzMesh::init(); + // Read and initialize the default mesh + HorzMesh::init(ModelClock); // Initialize the default vertical coordinate VertCoord::init(false); diff --git a/components/omega/test/ocn/VertCoordTest.cpp b/components/omega/test/ocn/VertCoordTest.cpp index e55dfa21f615..7cb37681c911 100644 --- a/components/omega/test/ocn/VertCoordTest.cpp +++ b/components/omega/test/ocn/VertCoordTest.cpp @@ -12,6 +12,7 @@ #include "Decomp.h" #include "Dimension.h" #include "Error.h" +#include "Field.h" #include "GlobalConstants.h" #include "Halo.h" #include "HorzMesh.h" @@ -45,8 +46,10 @@ int initVertCoordTest() { Config("Omega"); Config::readAll("omega.yml"); - // First step of time stepper initialization needed for IOstream + // First step of time stepper initialization needed for IOstream and clock TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); // Initialize the IO system IO::init(DefComm); @@ -55,7 +58,8 @@ int initVertCoordTest() { Decomp::init(); // Initialize streams - IOStream::init(); + Field::init(ModelClock); + IOStream::init(ModelClock); // Initialize the default halo Err = Halo::init(); @@ -63,7 +67,7 @@ int initVertCoordTest() { LOG_ERROR("VertCoordTest: error initializing default halo"); // Initialize the default mesh - HorzMesh::init(); + HorzMesh::init(ModelClock); // Initialize the default vertical coordinate VertCoord::init(); @@ -702,6 +706,7 @@ int main(int argc, char *argv[]) { TimeStepper::clear(); HorzMesh::clear(); VertCoord::clear(); + Field::clear(); Dimension::clear(); Halo::clear(); Decomp::clear(); diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index 96d1195e039f..54ee05769c7b 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -16,6 +16,7 @@ #include "Field.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OceanTestCommon.h" @@ -83,8 +84,19 @@ void initVertMixTest() { /// Initialize Halo Halo::init(); + /// Create dummy model clock for stream IO + Calendar::init("No Leap"); + TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); + TimeInterval TimeStep(1, TimeUnits::Hours); + Clock ModelClockTmp(StartTime, TimeStep); + Clock *ModelClock = &ModelClockTmp; + + /// Initialize IO streams for mesh IO + Field::init(ModelClock); + IOStream::init(ModelClock); + /// Initialize mesh - HorzMesh::init(); + HorzMesh::init(ModelClock); /// Initialize vertical coordinate VertCoord::init(false); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index aa2e2e89bbe2..c8eb24031b71 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -25,6 +25,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "IO.h" +#include "IOStream.h" #include "Logging.h" #include "MachEnv.h" #include "OceanState.h" @@ -149,9 +150,12 @@ int initTimeStepperTest(const std::string &mesh) { // Note that the default time stepper is not used in subsequent tests // but is initialized here because the number of time levels is needed - // to initialize the Tracers. If a later timestepper test uses more time - // levels than the default, this unit test may fail. + // to initialize the Tracers and a model clock is needed for stream IO. + // If a later timestepper test uses more time levels than the default, + // this unit test may fail. TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); IO::init(DefComm); Decomp::init(mesh); @@ -162,7 +166,12 @@ int initTimeStepperTest(const std::string &mesh) { LOG_ERROR("TimeStepperTest: error initializing default halo"); } - HorzMesh::init(); + // Initialize IO streams + Field::init(ModelClock); + IOStream::init(ModelClock); + + // Read and initialize horizontal mesh + HorzMesh::init(ModelClock); // initialize vertical coordinate, do not read stream and use local // NVertLayers value From 1eaa0c877c3a811a54c97d3559fe7da95efcc797 Mon Sep 17 00:00:00 2001 From: Philip W Jones Date: Tue, 7 Apr 2026 13:51:54 -0700 Subject: [PATCH 2/2] changes to PGradTest after rebase updated HorzMesh interfaces --- components/omega/test/ocn/PGradTest.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/components/omega/test/ocn/PGradTest.cpp b/components/omega/test/ocn/PGradTest.cpp index 09887a6d9ce6..f2fe555d5546 100644 --- a/components/omega/test/ocn/PGradTest.cpp +++ b/components/omega/test/ocn/PGradTest.cpp @@ -50,6 +50,8 @@ void initPGradTest() { // First step of time stepper initialization needed for IOstream TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); // Initialize the IO system IO::init(DefComm); @@ -57,8 +59,12 @@ void initPGradTest() { // Create the default decomposition (initializes the decomposition) Decomp::init(); - // Initialize streams - IOStream::init(); + // Initialize Field infrastructure + Field::init(ModelClock); + + // Initialize IOStreams - this does not yet validate the contents + // of each file, only creates streams from Config + IOStream::init(ModelClock); // Initialize the default halo Err1 = Halo::init(); @@ -68,7 +74,7 @@ void initPGradTest() { } // Initialize the default mesh - HorzMesh::init(); + HorzMesh::init(ModelClock); // Initialize the default vertical coordinate VertCoord::init();