From 548c9fe08f03a651e8ae052b17a520c182cbd205 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 25 Aug 2026 13:39:49 -0500 Subject: [PATCH 1/7] Add multi-layer support to VectorReconOnCell Add an overload that reconstructs one layer of a multi-layer field. The reconstruction is independent in each layer, so this is the existing single-layer form with a vertical index added. Factor the Cartesian to local geographic rotation into a shared helper so the spherical and planar paths stay identical between the two overloads. Extend the operator unit test to exercise the multi-layer form. Each layer holds the exact field scaled by a layer-dependent factor so that a reconstruction which mixed layers would show up as an error. The error measures are normalized by the exact field, so the same expected errors apply as for the single-layer case. Co-Authored-By: Claude Opus 5 --- components/omega/src/ocn/HorzOperators.h | 47 ++++++++++++--- .../omega/test/ocn/HorzOperatorsTest.cpp | 59 +++++++++++++++++-- 2 files changed, 93 insertions(+), 13 deletions(-) diff --git a/components/omega/src/ocn/HorzOperators.h b/components/omega/src/ocn/HorzOperators.h index bcb96b12dcab..3469ac65e464 100644 --- a/components/omega/src/ocn/HorzOperators.h +++ b/components/omega/src/ocn/HorzOperators.h @@ -657,8 +657,10 @@ class MasksAndCoefficients { class VectorReconOnCell { public: VectorReconOnCell(HorzMesh const *Mesh); - // Currently only support computing Zonal/Meridional (X/Y) for - // spherical (planar) meshes on a single vertical layer + + // Reconstruct a single vertical layer. The reconstructed components are + // zonal/meridional on spherical meshes and the Cartesian X/Y components + // on planar meshes. KOKKOS_FUNCTION void operator()(const Array1DReal &UReconX, const Array1DReal &UReconY, int ICell, const Array1DReal &VecEdge) const { @@ -674,6 +676,38 @@ class VectorReconOnCell { Uz += ReconWeightsCell(ICell, 2, J) * Field; } + cartesianToLocal(UReconX(ICell), UReconY(ICell), ICell, Ux, Uy, Uz); + } + + // Reconstruct one layer of a multi-layer field. The reconstruction is + // independent in each layer, so this is the single-layer form above + // with a vertical index added. + KOKKOS_FUNCTION void operator()(const Array2DReal &UReconX, + const Array2DReal &UReconY, int ICell, int K, + const Array2DReal &VecEdge) const { + + Real Ux = 0._Real, Uy = 0._Real, Uz = 0._Real; + + for (int J = 0; J < NEdgesReconOnCell(ICell); ++J) { + const I4 JEdge = ReconStencilCell(ICell, J); + const Real Field = VecEdge(JEdge, K); + + Ux += ReconWeightsCell(ICell, 0, J) * Field; + Uy += ReconWeightsCell(ICell, 1, J) * Field; + Uz += ReconWeightsCell(ICell, 2, J) * Field; + } + + cartesianToLocal(UReconX(ICell, K), UReconY(ICell, K), ICell, Ux, Uy, Uz); + } + + private: + // Convert a reconstructed Cartesian vector at a cell center into the + // components stored in the output arrays: local geographic + // (zonal/meridional) on spherical meshes, Cartesian X/Y on planar meshes + // where the vector already lies in the plane of the mesh. + KOKKOS_FUNCTION void cartesianToLocal(Real &UReconX, Real &UReconY, + int ICell, Real Ux, Real Uy, + Real Uz) const { if (OnSphere) { const Real CLat = Kokkos::cos(LatCell(ICell)); const Real SLat = Kokkos::sin(LatCell(ICell)); @@ -681,15 +715,14 @@ class VectorReconOnCell { const Real SLon = Kokkos::sin(LonCell(ICell)); // cartesian to local geographic - UReconX(ICell) = -SLon * Ux + CLon * Uy; - UReconY(ICell) = -(CLon * Ux + SLon * Uy) * SLat + Uz * CLat; + UReconX = -SLon * Ux + CLon * Uy; + UReconY = -(CLon * Ux + SLon * Uy) * SLat + Uz * CLat; } else { - UReconX(ICell) = Ux; - UReconY(ICell) = Uy; + UReconX = Ux; + UReconY = Uy; } } - private: bool OnSphere; Array1DI4 NEdgesReconOnCell; Array2DI4 ReconStencilCell; diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp index f4af88fd369c..6abfd657c6de 100644 --- a/components/omega/test/ocn/HorzOperatorsTest.cpp +++ b/components/omega/test/ocn/HorzOperatorsTest.cpp @@ -378,8 +378,9 @@ int testTangentRecon(Real RTol) { // Reconstructs the Cartesian vector field at cell centers from edge-normal // values using the least-squares weights/stencil in the mesh and compares // the magnitude of the reconstructed vector against the exact magnitude at -// cell centers. This currently only supports spherical meshes, so it is a -// no-op for planar meshes. +// cell centers. Both the single-layer and the multi-layer form of the +// operator are exercised. This currently only supports spherical meshes, +// so it is a no-op for planar meshes. int testVectorRecon(Real RTol) { int Err = 0; @@ -391,11 +392,10 @@ int testVectorRecon(Real RTol) { TestSetup Setup; - const auto &Mesh = HorzMesh::getDefault(); + const auto &Mesh = HorzMesh::getDefault(); + const int NVertLayers = 16; // Prepare operator input: edge-normal component of the exact vector field - // (VectorReconOnCell currently only supports a single vertical - // layer, so we use rank-1 arrays here) Array1DReal VecEdge("VecEdge", Mesh->NEdgesSize); Err += setVectorEdge( KOKKOS_LAMBDA(Real(&VecField)[2], Real Lon, Real Lat) { @@ -413,7 +413,7 @@ int testVectorRecon(Real RTol) { }, ExactMagCell, Geom, Mesh, OnCell, ExchangeHalos::No); - // Compute numerical reconstruction at cell centers + // Compute numerical reconstruction at cell centers, single layer Array1DReal UReconX("UReconX", Mesh->NCellsOwned); Array1DReal UReconY("UReconY", Mesh->NCellsOwned); VectorReconOnCell ReconCell(Mesh); @@ -436,6 +436,53 @@ int testVectorRecon(Real RTol) { Err += checkErrors("OperatorsTest", "VectorRecon", ReconErrors, Setup.ExpectedVectorReconErrors, RTol); + // Repeat with the multi-layer form of the operator. Each layer holds the + // exact field scaled by a layer-dependent factor, so a reconstruction + // that mixed layers or reused a single layer would show up as an error. + // Since the error measures are normalized by the exact field, a per-layer + // constant factor leaves them unchanged and the same expected errors + // apply as in the single-layer case above. + Array2DReal VecEdgeMulti("VecEdgeMulti", Mesh->NEdgesSize, NVertLayers); + Err += setVectorEdge( + KOKKOS_LAMBDA(Real(&VecField)[2], Real Lon, Real Lat) { + VecField[0] = Setup.exactVecX(Lon, Lat); + VecField[1] = Setup.exactVecY(Lon, Lat); + }, + VecEdgeMulti, EdgeComponent::Normal, Geom, Mesh); + + parallelFor( + {Mesh->NEdgesSize, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int K) { + VecEdgeMulti(IEdge, K) *= (1._Real + K); + }); + + Array2DReal ExactMagCellMulti("ExactMagCellMulti", Mesh->NCellsOwned, + NVertLayers); + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int ICell, int K) { + ExactMagCellMulti(ICell, K) = (1._Real + K) * ExactMagCell(ICell); + }); + + Array2DReal UReconXMulti("UReconXMulti", Mesh->NCellsOwned, NVertLayers); + Array2DReal UReconYMulti("UReconYMulti", Mesh->NCellsOwned, NVertLayers); + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int ICell, int K) { + ReconCell(UReconXMulti, UReconYMulti, ICell, K, VecEdgeMulti); + }); + + Array2DReal NumMagCellMulti("NumMagCellMulti", Mesh->NCellsOwned, + NVertLayers); + parallelFor( + {Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int ICell, int K) { + NumMagCellMulti(ICell, K) = + vecMagnitude(UReconXMulti(ICell, K), UReconYMulti(ICell, K)); + }); + + ErrorMeasures ReconErrorsMulti; + Err += computeErrors(ReconErrorsMulti, NumMagCellMulti, ExactMagCellMulti, + Mesh, OnCell); + Err += checkErrors("OperatorsTest", "VectorReconMultiLayer", + ReconErrorsMulti, Setup.ExpectedVectorReconErrors, RTol); + if (Err == 0) { LOG_INFO("OperatorsTest: VectorRecon PASS"); } From 46b925359092662d01c0d3478ef5c31bbf5eb8bc Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 25 Aug 2026 13:44:31 -0500 Subject: [PATCH 2/7] Gate vector reconstruction on mesh contents, not geometry The reconstruction stencil and weights are precomputed as a mesh preprocessing step rather than by Omega, so whether they are available is a property of the mesh file, not of whether the mesh is spherical. Detect their presence by attempting the read of NEdgesReconOnCell and treating a failure as "not present", and carry the answer as Decomp::HasVectorRecon and HorzMesh::HasVectorRecon. This lets planar meshes carry reconstruction data while leaving planar meshes that predate it readable, and makes VectorReconOnCell abort with a message naming the three arrays it needs rather than refusing every planar mesh. The planar branch already present in the operator becomes reachable as a result. Decomp no longer parses the OnSphere attribute, which existed only to gate these reads; HorzMesh still determines it from the full mesh stream as before. Co-Authored-By: Claude Opus 5 --- components/omega/doc/devGuide/Decomp.md | 16 ++- .../omega/doc/devGuide/HorzOperators.md | 11 ++ components/omega/doc/userGuide/Decomp.md | 8 +- components/omega/src/base/Decomp.cpp | 102 ++++++++---------- components/omega/src/base/Decomp.h | 6 +- components/omega/src/ocn/HorzMesh.cpp | 19 ++-- components/omega/src/ocn/HorzMesh.h | 14 +-- components/omega/src/ocn/HorzOperators.cpp | 8 +- components/omega/test/base/DecompTest.cpp | 12 +-- 9 files changed, 107 insertions(+), 89 deletions(-) diff --git a/components/omega/doc/devGuide/Decomp.md b/components/omega/doc/devGuide/Decomp.md index 9c02ada2f207..9ca3fd023fdb 100644 --- a/components/omega/doc/devGuide/Decomp.md +++ b/components/omega/doc/devGuide/Decomp.md @@ -88,11 +88,19 @@ described in the mesh specification above. In particular, it contains - NEdgesOnCell(NCellsSize): the number of actual edges on each cell - NEdgesOnEdge(NEdgesSize): the number of actual edges on each edge - NEdgesReconOnCell(NCellsSize): number of edges in the vector - reconstruction stencil for each cell (spherical meshes only) + reconstruction stencil for each cell (only for meshes that supply + the reconstruction arrays) - ReconStencilCell(NCellsSize,MaxEdges2): edge indices in the - vector reconstruction stencil for each cell (spherical meshes only) - - OnSphere: whether the mesh is spherical, read from the mesh file to - gate the reconstruction stencil arrays above + vector reconstruction stencil for each cell (only for meshes that + supply the reconstruction arrays) + - HasVectorRecon: whether the mesh file supplied the reconstruction + stencil arrays above. These are precomputed as a mesh preprocessing + step rather than by Omega, so they are absent from mesh files that + have not been through that step. Decomp detects this by attempting + the read of NEdgesReconOnCell and treating a failure as "not + present" rather than as an error, so the error logged by that read + is expected for such meshes. The paired ReconWeightsCell is read by + HorzMesh under the same flag. For each of the arrays above, there is a copy of the array on the host and device (GPU) with the host array named with an extra H on the end diff --git a/components/omega/doc/devGuide/HorzOperators.md b/components/omega/doc/devGuide/HorzOperators.md index 82efa155117c..a234f4d7efb6 100644 --- a/components/omega/doc/devGuide/HorzOperators.md +++ b/components/omega/doc/devGuide/HorzOperators.md @@ -39,6 +39,17 @@ Currently, the following operators are implemented: - `TangentialReconOnEdge` - `VectorReconOnCell` +`VectorReconOnCell` differs from the others in that it depends on +least-squares stencil and weight arrays (`NEdgesReconOnCell`, +`ReconStencilCell` and `ReconWeightsCell`) that are precomputed as a mesh +preprocessing step rather than by Omega. Constructing the operator on a +mesh whose file did not supply them (`HorzMesh::HasVectorRecon` is false) +is an error. It works on both spherical and planar meshes: on a sphere it +returns the local geographic (zonal and meridional) components, and on a +plane the Cartesian x and y components. It provides a single-layer form +and a form that takes a vertical index, for reconstructing one layer of +a multi-layer field. + Some tendency terms in the Omega PDE solver could in principle be constructed using these operators as building blocks. However, very often tendency terms require evaluation of slightly modified operators. Moreover, there is a diff --git a/components/omega/doc/userGuide/Decomp.md b/components/omega/doc/userGuide/Decomp.md index 938c368094ee..715618ca438a 100644 --- a/components/omega/doc/userGuide/Decomp.md +++ b/components/omega/doc/userGuide/Decomp.md @@ -39,8 +39,12 @@ An input mesh file must be provided that contains at a minimum - the total number of cells, edges and vertices (NCells, NEdges, NVertices) - the mesh connectivity contained in the arrays CellsOnCell, EdgesOnCell VerticesOnCell, CellsOnEdge, EdgesOnEdge, CellsOnVertex, EdgesOnVertex. -For spherical meshes, the vector reconstruction stencil arrays -NEdgesReconOnCell and ReconStencilCell are also required. +The vector reconstruction stencil arrays NEdgesReconOnCell and +ReconStencilCell are optional. They are precomputed as a mesh +preprocessing step, so only mesh files that have been through that step +contain them. A mesh without them is read normally, but reconstructing +vectors at cell centers (for example the zonal and meridional velocity +components) is then unavailable and requesting it is an error. Again, a full description of the mesh is given in the [Developer's Guide](#omega-dev-decomp). The file name for this input file is extracted from the HorzMeshIn input diff --git a/components/omega/src/base/Decomp.cpp b/components/omega/src/base/Decomp.cpp index c347a77d0feb..ac310035aafd 100644 --- a/components/omega/src/base/Decomp.cpp +++ b/components/omega/src/base/Decomp.cpp @@ -115,7 +115,7 @@ void readMesh( I4 &MaxCellsOnEdge, // max number of cells sharing edge I4 &VertexDegree, // number of cells/edges sharing vrtx I4 &MaxEdges2, // twice max number of edges on a cell - bool &OnSphere, // true if mesh is spherical + bool &HasVectorRecon, // true if mesh has vector reconstruction arrays std::vector &CellsOnCellInit, // cell neighbors for each cell std::vector &EdgesOnCellInit, // edge IDs for each cell edge std::vector &VerticesOnCellInit, // vertices around each cell @@ -224,19 +224,6 @@ void readMesh( I4 MaxEdgesOnEdge = MaxEdges2; // 2*MaxEdges, used below for // EdgesOnEdge/WeightsOnEdge offsets - // Determine whether the mesh is spherical or planar. This duplicates - // (temporarily) the OnSphere/on_a_sphere attribute parsing HorzMesh - // does via the full IOStream mechanism - here we just need a quick - // answer to decide whether the reconstruction stencil arrays (currently - // only generated for spherical meshes) are present in the file. - std::string OnSphereStr; - Err = IO::readMeta("OnSphere", OnSphereStr, MeshFileID, IO::GlobalID); - if (Err.isFail()) - Err = IO::readMeta("on_a_sphere", OnSphereStr, MeshFileID, IO::GlobalID); - std::transform(OnSphereStr.begin(), OnSphereStr.end(), OnSphereStr.begin(), - [](unsigned char c) { return std::tolower(c); }); - OnSphere = (OnSphereStr == "yes"); - // Create the linear decompositions for parallel IO // Determine the size of each block, divided as evenly as possible I4 NCellsChunk = (NCellsGlobal - 1) / NumTasks + 1; @@ -324,24 +311,20 @@ void readMesh( } // Create the parallel IO decompositions - IO::Rearranger Rearr = IO::RearrBox; - I4 OnCellDecomp = IO::createDecomp(IO::IOTypeI4, NDims, OnCellDims, - OnCellSize, OnCellOffset, Rearr); - I4 OnEdgeDecomp = IO::createDecomp(IO::IOTypeI4, NDims, OnEdgeDims, - OnEdgeSize, OnEdgeOffset, Rearr); - I4 OnEdgeDecomp2 = IO::createDecomp(IO::IOTypeI4, NDims, OnEdgeDims2, - OnEdgeSize2, OnEdgeOffset2, Rearr); - I4 OnVertexDecomp = IO::createDecomp(IO::IOTypeI4, NDims, OnVertexDims, - OnVertexSize, OnVertexOffset, Rearr); - I4 OnCellDecompScalar = -1; - I4 OnCellDecomp2 = -1; - if (OnSphere) { - OnCellDecompScalar = - IO::createDecomp(IO::IOTypeI4, 1, OnCellDimsScalar, OnCellSizeScalar, - OnCellOffsetScalar, Rearr); - OnCellDecomp2 = IO::createDecomp(IO::IOTypeI4, NDims, OnCellDims2, + IO::Rearranger Rearr = IO::RearrBox; + I4 OnCellDecomp = IO::createDecomp(IO::IOTypeI4, NDims, OnCellDims, + OnCellSize, OnCellOffset, Rearr); + I4 OnEdgeDecomp = IO::createDecomp(IO::IOTypeI4, NDims, OnEdgeDims, + OnEdgeSize, OnEdgeOffset, Rearr); + I4 OnEdgeDecomp2 = IO::createDecomp(IO::IOTypeI4, NDims, OnEdgeDims2, + OnEdgeSize2, OnEdgeOffset2, Rearr); + I4 OnVertexDecomp = IO::createDecomp(IO::IOTypeI4, NDims, OnVertexDims, + OnVertexSize, OnVertexOffset, Rearr); + I4 OnCellDecompScalar = + IO::createDecomp(IO::IOTypeI4, 1, OnCellDimsScalar, OnCellSizeScalar, + OnCellOffsetScalar, Rearr); + I4 OnCellDecomp2 = IO::createDecomp(IO::IOTypeI4, NDims, OnCellDims2, OnCellSize2, OnCellOffset2, Rearr); - } // Now read the connectivity arrays. Try reading under the new Omega // name convention and the older MPAS mesh names. @@ -444,23 +427,34 @@ void readMesh( // Vector reconstruction stencil - Omega-native fields, no legacy MPAS // name to fall back on. These are mesh-dependent, precomputed as a - // preprocessing step (least-squares pseudo-inverse). Only spherical - // meshes currently have these fields, so require them only in that case. - if (OnSphere) { - NEdgesReconOnCellInit.resize(OnCellSizeScalar); - ReconStencilCellInit.resize(OnCellSize2); - - VarName = "NEdgesReconOnCell"; - int NEdgesReconOnCellID; - Err = IO::readArray(&NEdgesReconOnCellInit[0], OnCellSizeScalar, VarName, - MeshFileID, OnCellDecompScalar, NEdgesReconOnCellID); - CHECK_ERROR_ABORT(Err, "Decomp: error reading NEdgesReconOnCell"); - + // preprocessing step (least-squares pseudo-inverse), so they are only + // present in mesh files that have been through that step. A mesh + // without them is not an error here - it simply cannot reconstruct + // vectors at cell centers, and the code that needs that capability + // (see VectorReconOnCell) aborts if HasVectorRecon is false. The first + // read failing is the detection mechanism, so the error it logs is + // expected for such meshes. + NEdgesReconOnCellInit.resize(OnCellSizeScalar); + ReconStencilCellInit.resize(OnCellSize2); + + VarName = "NEdgesReconOnCell"; + int NEdgesReconOnCellID; + Err = IO::readArray(&NEdgesReconOnCellInit[0], OnCellSizeScalar, VarName, + MeshFileID, OnCellDecompScalar, NEdgesReconOnCellID); + HasVectorRecon = !Err.isFail(); + + if (HasVectorRecon) { VarName = "ReconStencilCell"; int ReconStencilCellID; Err = IO::readArray(&ReconStencilCellInit[0], OnCellSize2, VarName, MeshFileID, OnCellDecomp2, ReconStencilCellID); CHECK_ERROR_ABORT(Err, "Decomp: error reading ReconStencilCell"); + } else { + LOG_INFO("Decomp: mesh file has no vector reconstruction arrays " + "(NEdgesReconOnCell); reconstruction of vectors at cell " + "centers will not be available"); + NEdgesReconOnCellInit.clear(); + ReconStencilCellInit.clear(); } // Initial decompositions are no longer needed so remove them now @@ -468,10 +462,8 @@ void readMesh( IO::destroyDecomp(OnEdgeDecomp); IO::destroyDecomp(OnEdgeDecomp2); IO::destroyDecomp(OnVertexDecomp); - if (OnSphere) { - IO::destroyDecomp(OnCellDecompScalar); - IO::destroyDecomp(OnCellDecomp2); - } + IO::destroyDecomp(OnCellDecompScalar); + IO::destroyDecomp(OnCellDecomp2); } // end readMesh @@ -580,7 +572,7 @@ Decomp::Decomp( HaloWidth = InHaloWidth; readMesh(FileID, InEnv, NCellsGlobal, NEdgesGlobal, NVerticesGlobal, - MaxEdges, MaxCellsOnEdge, VertexDegree, MaxEdges2, OnSphere, + MaxEdges, MaxCellsOnEdge, VertexDegree, MaxEdges2, HasVectorRecon, CellsOnCellInit, EdgesOnCellInit, VerticesOnCellInit, CellsOnEdgeInit, EdgesOnEdgeInit, VerticesOnEdgeInit, CellsOnVertexInit, EdgesOnVertexInit, NEdgesReconOnCellInit, @@ -635,9 +627,9 @@ Decomp::Decomp( // Redistribute the vector reconstruction stencil arrays to the same // final cell decomposition. This can happen as soon as CellID/CellLoc // are finalized above - it does not participate in defining the - // decomposition itself, unlike CellsOnCellInit. Only spherical meshes - // currently have these arrays. - if (OnSphere) { + // decomposition itself, unlike CellsOnCellInit. Skipped for meshes + // without the reconstruction arrays. + if (HasVectorRecon) { TimerFlag = Pacer::start("Decomp rearrange recon stencil", 2) && TimerFlag; rearrangeReconArrays(InEnv, NEdgesReconOnCellInit, ReconStencilCellInit); @@ -732,8 +724,8 @@ Decomp::Decomp( // ReconStencilCell - translated the same way as EdgesOnCell // NEdgeReconOnCellH is a count, not an ID, and needs no translation. - // Only spherical meshes currently have this array. - if (OnSphere) { + // Skipped for meshes without the reconstruction arrays. + if (HasVectorRecon) { for (int Cell = 0; Cell < NCellsSize; ++Cell) { for (int Edge = 0; Edge < MaxEdges2; ++Edge) { I4 GlobID = ReconStencilCellH(Cell, Edge); @@ -878,8 +870,8 @@ Decomp::Decomp( CellsOnVertex = createDeviceMirrorCopy(CellsOnVertexH); EdgesOnVertex = createDeviceMirrorCopy(EdgesOnVertexH); - // Only spherical meshes currently have the reconstruction stencil arrays - if (OnSphere) { + // Only meshes with reconstruction data have the stencil arrays + if (HasVectorRecon) { NEdgesReconOnCell = createDeviceMirrorCopy(NEdgesReconOnCellH); ReconStencilCell = createDeviceMirrorCopy(ReconStencilCellH); } diff --git a/components/omega/src/base/Decomp.h b/components/omega/src/base/Decomp.h index 47c7530fc5c5..0ac73214d2a8 100644 --- a/components/omega/src/base/Decomp.h +++ b/components/omega/src/base/Decomp.h @@ -278,9 +278,9 @@ class Decomp { // Vector reconstruction stencil (mesh-dependent, precomputed and stored // in the mesh file - see HorzMesh for the paired ReconWeightsCell) - bool OnSphere; ///< true if mesh is spherical (temporary local read of - ///< the OnSphere attribute - only spherical meshes - ///< currently have the reconstruction stencil below) + bool HasVectorRecon; ///< true if the mesh file supplied the vector + ///< reconstruction arrays below (and the paired + ///< ReconWeightsCell read by HorzMesh) Array1DI4 NEdgesReconOnCell; ///< Num of edges in reconstruction stencil HostArray1DI4 NEdgesReconOnCellH; ///< Num of edges in reconstruction stencil diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 78f041584e1d..dea53dac13b5 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -114,13 +114,12 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh NEdgesReconOnCell = MeshDecomp->NEdgesReconOnCell; ReconStencilCell = MeshDecomp->ReconStencilCell; - // OnSphere is needed early (before the mesh field definitions below) - // to decide whether to define the ReconWeightsCell field, so we - // use the value Decomp already read from the mesh file. It is - // redetermined below (with the rest of the sphere/plane attributes) - // once the full mesh stream is read, which is redundant but harmless - // since both come from the same file. - OnSphere = MeshDecomp->OnSphere; + // Whether the mesh file supplied the vector reconstruction arrays is + // needed early (before the mesh field definitions below) to decide + // whether to define the ReconWeightsCell field, so we take the answer + // from Decomp, which already attempted to read its half of those + // arrays from the same file. + HasVectorRecon = MeshDecomp->HasVectorRecon; // Create Omega Dimensions associated with this mesh createDimensions(MeshDecomp); @@ -477,7 +476,7 @@ void HorzMesh::completeReadArrays() { // the mesh dimension in the middle (as in Tracers' [Tracer, Cell, // Vert]), so exchange each R3 component as a 2D (Cell, MaxEdges2) // slice instead of the full 3D array. - if (OnSphere) { + if (HasVectorRecon) { for (int IComp = 0; IComp < 3; ++IComp) { auto ReconWeightsCellSlice = Kokkos::subview(ReconWeightsCell, Kokkos::ALL, IComp, Kokkos::ALL); @@ -512,7 +511,7 @@ void HorzMesh::completeReadArrays() { FCellH = createHostMirrorCopy(FCell); FEdgeH = createHostMirrorCopy(FEdge); FVertexH = createHostMirrorCopy(FVertex); - if (OnSphere) { + if (HasVectorRecon) { ReconWeightsCellH = createHostMirrorCopy(ReconWeightsCell); } @@ -1037,7 +1036,7 @@ void HorzMesh::defineMeshFields() { MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, FCell); - if (OnSphere) { + if (HasVectorRecon) { // Vector Reconstruction // NEdgesReconOnCell/ReconStencilCell come from Decomp (see // constructor), not read here. ReconWeightsCell is pure data and diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index c9eb3f2f755b..69ee24a86ee4 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -79,12 +79,14 @@ class HorzMesh { std::string MeshFileName; ///< name (and full path) of input mesh file // Global attributes - bool OnSphere; ///< true if mesh is spherical - bool OnPlane; ///< true if mesh is on a Cartesian plane - bool IsPeriodic; ///< true if planar mesh is periodic - Real SphereRadius; ///< radius (m) of sphere that mesh covers - Real XPeriod; ///< length (m) of periodicity in x direction - Real YPeriod; ///< length (m) of periodicity in y direction + bool OnSphere; ///< true if mesh is spherical + bool HasVectorRecon; ///< true if the mesh file supplied the vector + ///< reconstruction stencil and weights below + bool OnPlane; ///< true if mesh is on a Cartesian plane + bool IsPeriodic; ///< true if planar mesh is periodic + Real SphereRadius; ///< radius (m) of sphere that mesh covers + Real XPeriod; ///< length (m) of periodicity in x direction + Real YPeriod; ///< length (m) of periodicity in y direction // Sizes and global IDs // Note that all sizes are actual counts (1-based) so that loop extents diff --git a/components/omega/src/ocn/HorzOperators.cpp b/components/omega/src/ocn/HorzOperators.cpp index a530e8254dd1..05ef352f5225 100644 --- a/components/omega/src/ocn/HorzOperators.cpp +++ b/components/omega/src/ocn/HorzOperators.cpp @@ -72,9 +72,11 @@ VectorReconOnCell::VectorReconOnCell(HorzMesh const *Mesh) ReconStencilCell(Mesh->ReconStencilCell), ReconWeightsCell(Mesh->ReconWeightsCell), LatCell(Mesh->LatCell), LonCell(Mesh->LonCell) { - if (!Mesh->OnSphere) - ABORT_ERROR("VectorReconOnCell: reconstruction stencil/weights " - "are only available for spherical meshes"); + if (!Mesh->HasVectorRecon) + ABORT_ERROR("VectorReconOnCell: mesh {} has no vector reconstruction " + "data; the mesh file must supply NEdgesReconOnCell, " + "ReconStencilCell and ReconWeightsCell", + Mesh->MeshName); } } // namespace OMEGA diff --git a/components/omega/test/base/DecompTest.cpp b/components/omega/test/base/DecompTest.cpp index 72c5f2de5eb7..4cd4a78e2841 100644 --- a/components/omega/test/base/DecompTest.cpp +++ b/components/omega/test/base/DecompTest.cpp @@ -127,12 +127,12 @@ int main(int argc, char *argv[]) { ABORT_ERROR("DecompTest: Sum vertex ID test FAIL {} {}", SumVertices, RefSumVertices); - // Test the vector-reconstruction stencil arrays (spherical meshes - // only): each owned cell's count is in range, active entries are - // resolvable local edges, and padding columns carry the NEdgesAll - // sentinel (confirms no compaction, so columns stay aligned with - // ReconWeightsCell). - if (DefDecomp->OnSphere) { + // Test the vector-reconstruction stencil arrays, for meshes that + // supply them: each owned cell's count is in range, active entries + // are resolvable local edges, and padding columns carry the + // NEdgesAll sentinel (confirms no compaction, so columns stay + // aligned with ReconWeightsCell). + if (DefDecomp->HasVectorRecon) { I4 MaxEdges2 = DefDecomp->MaxEdges2; HostArray1DI4 NEdgesReconOnCellH = DefDecomp->NEdgesReconOnCellH; HostArray2DI4 ReconStencilCellH = DefDecomp->ReconStencilCellH; From fa33fbd4b3cfc53bfb4bdfc65d46ad452e415929 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 25 Aug 2026 13:48:17 -0500 Subject: [PATCH 3/7] Add VelocityZonalCell and VelocityMeridionalCell fields Add a VelocityReconAuxVars auxiliary variable class holding the zonal and meridional velocity components reconstructed at cell centers from the edge-normal velocity, and register both as fields in the AuxiliaryState field group so they can be requested by any IO stream. On a planar mesh the reconstructed vector already lies in the plane of the mesh, so the two components are the Cartesian x and y components. The reconstruction is written out here rather than calling VectorReconOnCell so that layers outside MinLayerCell/MaxLayerCell keep the fill value that Field::attachData wrote, which the operator, as a reference implementation, knows nothing about. Add AuxiliaryState::computeVelocityRecon to drive it. These are diagnostics that no tendency reads, so the computation is deliberately kept out of computeAll, which runs once per time stepper stage, and is over owned cells only. It aborts if the mesh supplied no reconstruction data. Nothing calls it yet. Co-Authored-By: Claude Opus 5 --- components/omega/src/ocn/AuxiliaryState.cpp | 35 ++++++++ components/omega/src/ocn/AuxiliaryState.h | 8 ++ .../auxiliaryVars/VelocityReconAuxVars.cpp | 78 +++++++++++++++++ .../ocn/auxiliaryVars/VelocityReconAuxVars.h | 87 +++++++++++++++++++ 4 files changed, 208 insertions(+) create mode 100644 components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.cpp create mode 100644 components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.h diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 9a37a0739d98..e4ed53152f92 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -29,6 +29,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux(stripDefault(Name), Mesh, VCoord), VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), SurfTracerRestAux(stripDefault(Name), Mesh, NTracers), + VelocityReconAux(stripDefault(Name), Mesh, VCoord), TracerAux(stripDefault(Name), Mesh, VCoord, NTracers), TimeStep(TimeStep) { @@ -45,6 +46,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux.registerFields(GroupName, AuxMeshName); VelocityDel2Aux.registerFields(GroupName, AuxMeshName); SurfTracerRestAux.registerFields(GroupName, AuxMeshName); + VelocityReconAux.registerFields(GroupName, AuxMeshName); TracerAux.registerFields(GroupName, AuxMeshName); } @@ -56,6 +58,7 @@ AuxiliaryState::~AuxiliaryState() { VorticityAux.unregisterFields(); VelocityDel2Aux.unregisterFields(); SurfTracerRestAux.unregisterFields(); + VelocityReconAux.unregisterFields(); TracerAux.unregisterFields(); FieldGroup::destroy(GroupName); @@ -323,6 +326,38 @@ void AuxiliaryState::computeAll(const OceanState *State, computeAll(State, TracerArray, TimeLevel, TimeLevel, ProjDt); } +// Compute the diagnostic zonal and meridional velocity components at cell +// centers. These are not used by the Omega equations, so this is kept out +// of computeAll (which runs once per time stepper stage) and is instead +// called once per time step. +void AuxiliaryState::computeVelocityRecon(const OceanState *State, + int VelTimeLevel) const { + + if (!Mesh->HasVectorRecon) + ABORT_ERROR("AuxiliaryState: {} and {} were requested but mesh {} has " + "no vector reconstruction data; the mesh file must supply " + "NEdgesReconOnCell, ReconStencilCell and ReconWeightsCell", + VelocityReconAux.VelocityZonalCell.label(), + VelocityReconAux.VelocityMeridionalCell.label(), + Mesh->MeshName); + + Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); + + OMEGA_SCOPE(LocVelocityReconAux, VelocityReconAux); + + Pacer::start("AuxState:computeVelocityRecon", 1); + + // Only owned cells are needed: these are diagnostics written to output, + // not inputs to any tendency that would read them in the halo. + parallelFor( + "velocityReconAuxState", {Mesh->NCellsOwned, VCoord->NVertLayers}, + KOKKOS_LAMBDA(int ICell, int K) { + LocVelocityReconAux.computeVarsOnCell(ICell, K, NormalVelEdge); + }); + + Pacer::stop("AuxState:computeVelocityRecon", 1); +} + // Create a non-default auxiliary state AuxiliaryState *AuxiliaryState::create(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 91047d6e160e..949c43bc2dda 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -15,6 +15,7 @@ #include "auxiliaryVars/SurfTracerRestAuxVars.h" #include "auxiliaryVars/TracerAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" +#include "auxiliaryVars/VelocityReconAuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" #include @@ -42,6 +43,7 @@ class AuxiliaryState { VorticityAuxVars VorticityAux; VelocityDel2AuxVars VelocityDel2Aux; SurfTracerRestAuxVars SurfTracerRestAux; + VelocityReconAuxVars VelocityReconAux; ~AuxiliaryState(); @@ -84,6 +86,12 @@ class AuxiliaryState { int ThickTimeLevel, int VelTimeLevel, const TimeInterval ProjDt) const; + // Compute the diagnostic zonal and meridional velocity components at + // cell centers. Unlike the auxiliary variables above, nothing in the + // Omega equations uses these, so they are computed once per time step + // rather than once per time stepper stage. + void computeVelocityRecon(const OceanState *State, int VelTimeLevel) const; + /// Compute all auxiliary variables based on an ocean state at a given time /// level void computeAll(const OceanState *State, const Array3DReal &TracerArray, diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.cpp new file mode 100644 index 000000000000..370ccdc87867 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.cpp @@ -0,0 +1,78 @@ +#include "VelocityReconAuxVars.h" +#include "DataTypes.h" +#include "Field.h" + +#include + +namespace OMEGA { + +VelocityReconAuxVars::VelocityReconAuxVars(const std::string &AuxStateSuffix, + const HorzMesh *Mesh, + const VertCoord *VCoord) + : VelocityZonalCell("VelocityZonalCell" + AuxStateSuffix, Mesh->NCellsSize, + VCoord->NVertLayers), + VelocityMeridionalCell("VelocityMeridionalCell" + AuxStateSuffix, + Mesh->NCellsSize, VCoord->NVertLayers), + OnSphere(Mesh->OnSphere), NEdgesReconOnCell(Mesh->NEdgesReconOnCell), + ReconStencilCell(Mesh->ReconStencilCell), + ReconWeightsCell(Mesh->ReconWeightsCell), LatCell(Mesh->LatCell), + LonCell(Mesh->LonCell), MinLayerCell(VCoord->MinLayerCell), + MaxLayerCell(VCoord->MaxLayerCell) {} + +void VelocityReconAuxVars::registerFields( + const std::string &AuxGroupName, // name of Auxiliary field group + const std::string &MeshName // name of horizontal mesh +) const { + + // Create fields + int NDims = 2; + std::vector DimNames(NDims); + std::string DimSuffix; + if (MeshName == "Default") { + DimSuffix = ""; + } else { + DimSuffix = MeshName; + } + + DimNames[0] = "NCells" + DimSuffix; + DimNames[1] = "NVertLayers"; + + // Zonal velocity on cells + auto VelocityZonalCellField = Field::create( + VelocityZonalCell.label(), // field name + "zonal velocity reconstructed at cell centers", // long name/describe + "m s^-1", // units + "eastward_sea_water_velocity", // CF standard Name + std::numeric_limits::lowest(), // min valid value + std::numeric_limits::max(), // max valid value + NDims, // number of dimensions + DimNames // dimension names + ); + + // Meridional velocity on cells + auto VelocityMeridionalCellField = Field::create( + VelocityMeridionalCell.label(), // field name + "meridional velocity reconstructed at cell centers", // long name + "m s^-1", // units + "northward_sea_water_velocity", // CF standard Name + std::numeric_limits::lowest(), // min valid value + std::numeric_limits::max(), // max valid value + NDims, // number of dimensions + DimNames // dimension names + ); + + // Add fields to FieldGroup + FieldGroup::addFieldToGroup(VelocityZonalCell.label(), AuxGroupName); + FieldGroup::addFieldToGroup(VelocityMeridionalCell.label(), AuxGroupName); + + // Attach data + VelocityZonalCellField->attachData(VelocityZonalCell); + VelocityMeridionalCellField->attachData(VelocityMeridionalCell); +} + +void VelocityReconAuxVars::unregisterFields() const { + Field::destroy(VelocityZonalCell.label()); + Field::destroy(VelocityMeridionalCell.label()); +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.h b/components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.h new file mode 100644 index 000000000000..1f94504c0bf9 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/VelocityReconAuxVars.h @@ -0,0 +1,87 @@ +#ifndef OMEGA_AUX_VELOCITY_RECON_H +#define OMEGA_AUX_VELOCITY_RECON_H + +#include "DataTypes.h" +#include "HorzMesh.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" + +#include + +namespace OMEGA { + +/// Zonal and meridional velocity components reconstructed at cell centers +/// from the edge-normal velocity, using the least-squares stencil and +/// weights supplied by the mesh file. On a planar mesh the reconstructed +/// vector already lies in the plane of the mesh, so the two components are +/// the Cartesian x and y components instead. +/// +/// These are diagnostic: nothing in the Omega equations uses them, so they +/// are computed once per time step rather than once per time stepper stage +/// (see AuxiliaryState::computeVelocityRecon). +/// +/// This holds the same reconstruction as the VectorReconOnCell horizontal +/// operator, which is its reference implementation and is what the +/// operator unit test exercises, with the valid layer range of a column +/// taken into account. +class VelocityReconAuxVars { + public: + Array2DReal VelocityZonalCell; + Array2DReal VelocityMeridionalCell; + + VelocityReconAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, + const VertCoord *VCoord); + + KOKKOS_FUNCTION void + computeVarsOnCell(int ICell, int K, const Array2DReal &NormalVelEdge) const { + + // Leave layers outside the valid range of this column at the fill + // value that Field::attachData wrote + if (K < MinLayerCell(ICell) || K > MaxLayerCell(ICell)) + return; + + // Accumulate the Cartesian components of the reconstructed vector + Real Ux = 0._Real, Uy = 0._Real, Uz = 0._Real; + + for (int J = 0; J < NEdgesReconOnCell(ICell); ++J) { + const I4 JEdge = ReconStencilCell(ICell, J); + const Real Field = NormalVelEdge(JEdge, K); + + Ux += ReconWeightsCell(ICell, 0, J) * Field; + Uy += ReconWeightsCell(ICell, 1, J) * Field; + Uz += ReconWeightsCell(ICell, 2, J) * Field; + } + + if (OnSphere) { + // cartesian to local geographic + const Real CLat = Kokkos::cos(LatCell(ICell)); + const Real SLat = Kokkos::sin(LatCell(ICell)); + const Real CLon = Kokkos::cos(LonCell(ICell)); + const Real SLon = Kokkos::sin(LonCell(ICell)); + + VelocityZonalCell(ICell, K) = -SLon * Ux + CLon * Uy; + VelocityMeridionalCell(ICell, K) = + -(CLon * Ux + SLon * Uy) * SLat + Uz * CLat; + } else { + VelocityZonalCell(ICell, K) = Ux; + VelocityMeridionalCell(ICell, K) = Uy; + } + } + + void registerFields(const std::string &AuxGroupName, + const std::string &MeshName) const; + void unregisterFields() const; + + private: + bool OnSphere; + Array1DI4 NEdgesReconOnCell; + Array2DI4 ReconStencilCell; + Array3DReal ReconWeightsCell; + Array1DReal LatCell; + Array1DReal LonCell; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + +} // namespace OMEGA +#endif From f499ed264d3ee4a76cdb8be0ca31e083837a8f91 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 25 Aug 2026 13:50:30 -0500 Subject: [PATCH 4/7] Compute the velocity components once per step, on demand Add IOStream::isFieldRequested, which reports whether a Field is in the contents of any defined stream, so a diagnostic that nothing will read or write can be skipped. Since validate replaces a group name in a stream's contents by its member Fields, this is only meaningful after the streams are validated; an unvalidated stream conservatively answers yes so an optional computation is done rather than skipped. Use it to gate AuxiliaryState::computeVelocityRecon, resolving the answer on the first call and caching it, since the streams are validated only after all Fields are defined and so well after the auxiliary state is constructed. Call it once per time step from both ocnRun loops, on the state the step just produced (time level 0), which is the state that updateTimeLevels has attached to the NormalVelocity field and so the one that any stream writing this step will write. Co-Authored-By: Claude Opus 5 --- components/omega/src/infra/IOStream.cpp | 21 +++++++++++++++++++++ components/omega/src/infra/IOStream.h | 11 +++++++++++ components/omega/src/ocn/AuxiliaryState.cpp | 17 +++++++++++++++-- components/omega/src/ocn/AuxiliaryState.h | 13 +++++++++++-- components/omega/src/ocn/OceanRun.cpp | 13 +++++++++++++ 5 files changed, 71 insertions(+), 4 deletions(-) diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp index 2a1cca616766..3c78a5a98f98 100644 --- a/components/omega/src/infra/IOStream.cpp +++ b/components/omega/src/infra/IOStream.cpp @@ -274,6 +274,27 @@ bool IOStream::validateAll() { } // End validateAll +//------------------------------------------------------------------------------ +// Determines whether a Field is in the contents of any defined stream. +bool IOStream::isFieldRequested(const std::string &FieldName) { + + for (auto Iter = AllStreams.begin(); Iter != AllStreams.end(); Iter++) { + std::shared_ptr ThisStream = Iter->second; + + // An unvalidated stream may still hold group names rather than the + // names of the group's member Fields, so we cannot tell whether it + // contains this Field. Assume that it might. + if (!ThisStream->Validated) + return true; + + if (ThisStream->Contents.find(FieldName) != ThisStream->Contents.end()) + return true; + } + + return false; + +} // End isFieldRequested + //------------------------------------------------------------------------------ // Reads a single stream if it is time. Error IOStream::read( diff --git a/components/omega/src/infra/IOStream.h b/components/omega/src/infra/IOStream.h index 7f5fe7d82c71..6cb02ba1e2b7 100644 --- a/components/omega/src/infra/IOStream.h +++ b/components/omega/src/infra/IOStream.h @@ -346,6 +346,17 @@ class IOStream { /// Returns true if all streams are valid. static bool validateAll(); + //--------------------------------------------------------------------------- + /// Determines whether a Field is in the contents of any defined stream, + /// so that a caller can skip computing a diagnostic that nothing will + /// read or write. Because validate replaces a group name in a stream's + /// contents by the names of the group's member Fields, this is only + /// meaningful after the streams have been validated. If any stream has + /// not been validated, the answer is conservatively true so that an + /// optional computation is done rather than skipped. + static bool isFieldRequested(const std::string &FieldName ///< [in] Field + ); + //--------------------------------------------------------------------------- /// Reads a stream if it is time. static Error read(const std::string &StreamName, ///< [in] Name of stream diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index e4ed53152f92..0d788db54be2 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -2,6 +2,7 @@ #include "Config.h" #include "Error.h" #include "Field.h" +#include "IOStream.h" #include "Logging.h" #include "Pacer.h" #include "Tendencies.h" @@ -329,9 +330,21 @@ void AuxiliaryState::computeAll(const OceanState *State, // Compute the diagnostic zonal and meridional velocity components at cell // centers. These are not used by the Omega equations, so this is kept out // of computeAll (which runs once per time stepper stage) and is instead -// called once per time step. +// called once per time step, and only if some IO stream asks for them. void AuxiliaryState::computeVelocityRecon(const OceanState *State, - int VelTimeLevel) const { + int VelTimeLevel) { + + if (!VelocityReconResolved) { + VelocityReconRequested = + IOStream::isFieldRequested( + VelocityReconAux.VelocityZonalCell.label()) || + IOStream::isFieldRequested( + VelocityReconAux.VelocityMeridionalCell.label()); + VelocityReconResolved = true; + } + + if (!VelocityReconRequested) + return; if (!Mesh->HasVectorRecon) ABORT_ERROR("AuxiliaryState: {} and {} were requested but mesh {} has " diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 949c43bc2dda..d686a9ed9090 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -89,8 +89,10 @@ class AuxiliaryState { // Compute the diagnostic zonal and meridional velocity components at // cell centers. Unlike the auxiliary variables above, nothing in the // Omega equations uses these, so they are computed once per time step - // rather than once per time stepper stage. - void computeVelocityRecon(const OceanState *State, int VelTimeLevel) const; + // rather than once per time stepper stage, and only if some IO stream + // asks for them. Not const because the answer to that question is + // resolved on the first call and cached. + void computeVelocityRecon(const OceanState *State, int VelTimeLevel); /// Compute all auxiliary variables based on an ocean state at a given time /// level @@ -114,6 +116,13 @@ class AuxiliaryState { VertAdv *VAdv; TimeInterval TimeStep; + /// Whether any IO stream asks for the reconstructed velocity + /// components, resolved on the first call to computeVelocityRecon. + /// This cannot be answered when the auxiliary state is constructed, + /// since the streams are validated only after all Fields are defined. + bool VelocityReconRequested = false; + bool VelocityReconResolved = false; + static AuxiliaryState *DefaultAuxState; static std::map> AllAuxStates; }; diff --git a/components/omega/src/ocn/OceanRun.cpp b/components/omega/src/ocn/OceanRun.cpp index 65653bfbb3ec..9a8f154a9b68 100644 --- a/components/omega/src/ocn/OceanRun.cpp +++ b/components/omega/src/ocn/OceanRun.cpp @@ -6,6 +6,7 @@ //===----------------------------------------------------------------------===// #include "Analysis.h" +#include "AuxiliaryState.h" #include "Forcing.h" #include "IOStream.h" #include "OceanDriver.h" @@ -65,6 +66,12 @@ int ocnRun(TimeInstant &CurrTime ///< [inout] current sim time Pacer::stop("Stepper:doStep", 1); } + // Compute diagnostics that are not needed by the time stepper, using + // the state the step just produced (time level 0). This is a no-op + // unless an IO stream asks for them. + AuxiliaryState *DefAuxState = AuxiliaryState::getDefault(); + DefAuxState->computeVelocityRecon(DefOceanState, 0); + // Compute analysis fields whose alarms are ringing Analysis *DefAnalysis = Analysis::getDefault(); DefAnalysis->computeAll(); @@ -135,6 +142,12 @@ int ocnRun(TimeInstant &CurrTime, ///< [inout] current sim time Pacer::stop("Stepper:doStep", 1); } + // Compute diagnostics that are not needed by the time stepper, using + // the state the step just produced (time level 0). This is a no-op + // unless an IO stream asks for them. + AuxiliaryState *DefAuxState = AuxiliaryState::getDefault(); + DefAuxState->computeVelocityRecon(DefOceanState, 0); + // Write any IOStreams with their alarms ringing IOStream::writeAll(OmegaClock); From c4c8b007e140455622e551509891b32ef0297c78 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 25 Aug 2026 13:54:41 -0500 Subject: [PATCH 5/7] Export the reconstructed surface velocity to the coupler Replace the hardcoded 1e-4 placeholder in OcnToCplFields::updateFields with the reconstructed zonal and meridional velocity in the surface layer, and drop the TODO it carried. Add AuxiliaryState::requireVelocityRecon so surface coupling can force the reconstruction to be computed every step regardless of what the IO streams ask for, and call it from SfcCoupling::init. That init also now checks the mesh supplies the reconstruction data, so a coupled run on a mesh without it fails at initialization rather than at the first time step. Initialize VertAdv and AuxiliaryState in SfcCouplingTest, which the coupling now depends on, and check the exported surface velocity against the reconstruction rather than skipping it. At NAccumSteps == 0 the running average is exactly the single sampled value, so the comparison is exact, as it already is for temperature and salinity. Co-Authored-By: Claude Opus 5 --- components/omega/src/ocn/AuxiliaryState.cpp | 7 +++ components/omega/src/ocn/AuxiliaryState.h | 5 ++ components/omega/src/ocn/SfcCoupling.cpp | 38 +++++++++++--- components/omega/test/ocn/SfcCouplingTest.cpp | 51 +++++++++++++++++-- 4 files changed, 90 insertions(+), 11 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 0d788db54be2..1ffe5a1fb6f9 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -327,6 +327,13 @@ void AuxiliaryState::computeAll(const OceanState *State, computeAll(State, TracerArray, TimeLevel, TimeLevel, ProjDt); } +// Force the reconstructed velocity components to be computed every time +// step, whatever the IO streams ask for. +void AuxiliaryState::requireVelocityRecon() { + VelocityReconRequested = true; + VelocityReconResolved = true; +} + // Compute the diagnostic zonal and meridional velocity components at cell // centers. These are not used by the Omega equations, so this is kept out // of computeAll (which runs once per time stepper stage) and is instead diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index d686a9ed9090..6821ea4056a0 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -94,6 +94,11 @@ class AuxiliaryState { // resolved on the first call and cached. void computeVelocityRecon(const OceanState *State, int VelTimeLevel); + // Force the reconstructed velocity components to be computed every time + // step even if no IO stream asks for them. Used by surface coupling, + // which needs them regardless of what is being written. + void requireVelocityRecon(); + /// Compute all auxiliary variables based on an ocean state at a given time /// level void computeAll(const OceanState *State, const Array3DReal &TracerArray, diff --git a/components/omega/src/ocn/SfcCoupling.cpp b/components/omega/src/ocn/SfcCoupling.cpp index 3c60fbe17c19..1f9057a3fe15 100644 --- a/components/omega/src/ocn/SfcCoupling.cpp +++ b/components/omega/src/ocn/SfcCoupling.cpp @@ -1,4 +1,5 @@ #include "SfcCoupling.h" +#include "AuxiliaryState.h" #include "Eos.h" #include "Error.h" #include "GlobalConstants.h" @@ -28,6 +29,22 @@ int SfcCoupling::init(const CouplingInitParams &CouplingInitParams) { OMEGA_REQUIRE(DefTimeStepper, "Null default TimeStepper pointer in SfcCoupling::init"); + // The surface velocity exported to the coupler is reconstructed at cell + // centers from the edge-normal velocity, so a coupled run needs a mesh + // that supplies the reconstruction stencil and weights. Fail here + // rather than at the first time step. + if (!DefHorzMesh->HasVectorRecon) + ABORT_ERROR("SfcCoupling: mesh {} has no vector reconstruction data, " + "which is needed for the surface velocity export; the mesh " + "file must supply NEdgesReconOnCell, ReconStencilCell and " + "ReconWeightsCell", + DefHorzMesh->MeshName); + + AuxiliaryState *DefAuxState = AuxiliaryState::getDefault(); + OMEGA_REQUIRE(DefAuxState, + "Null default AuxiliaryState pointer in SfcCoupling::init"); + DefAuxState->requireVelocityRecon(); + TimeInterval OcnTimeStep = DefTimeStepper->getTimeStep(); TimeInterval CplTimeStep = CouplingInitParams.CouplingTimeStep; @@ -317,7 +334,8 @@ void OcnToCplFields::updateFields(const OceanState *State, auto Salinity = Kokkos::subview(TracerArray, SalinityIdx, Kokkos::ALL, Kokkos::ALL); - VertCoord *DefVertCoord = VertCoord::getDefault(); + VertCoord *DefVertCoord = VertCoord::getDefault(); + AuxiliaryState *DefAuxState = AuxiliaryState::getDefault(); OMEGA_SCOPE(LocMinLayerCell, DefVertCoord->MinLayerCell); OMEGA_SCOPE(LocAvgSfcSalinity, AvgSfcSalinity); @@ -325,8 +343,12 @@ void OcnToCplFields::updateFields(const OceanState *State, OMEGA_SCOPE(LocAvgSfcVelZonal, AvgSfcVelocityZonal); OMEGA_SCOPE(LocAvgSfcVelMerid, AvgSfcVelocityMerid); - // TODO: Implement vector reconsturction for velocity field. - constexpr Real ConstSfcVelocity = 1e-4; + // The reconstruction is computed once per time step in ocnRun, before + // this is called (see AuxiliaryState::computeVelocityRecon). + OMEGA_SCOPE(LocVelocityZonal, + DefAuxState->VelocityReconAux.VelocityZonalCell); + OMEGA_SCOPE(LocVelocityMerid, + DefAuxState->VelocityReconAux.VelocityMeridionalCell); parallelFor( {NCellsOwned}, KOKKOS_LAMBDA(int ICell) { @@ -339,11 +361,13 @@ void OcnToCplFields::updateFields(const OceanState *State, LocAvgSfcSalinity(ICell) = updateAverage( LocAvgSfcSalinity(ICell), Salinity(ICell, KSfc), NAccumSteps); - LocAvgSfcVelZonal(ICell) = updateAverage( - LocAvgSfcVelZonal(ICell), ConstSfcVelocity, NAccumSteps); + LocAvgSfcVelZonal(ICell) = + updateAverage(LocAvgSfcVelZonal(ICell), + LocVelocityZonal(ICell, KSfc), NAccumSteps); - LocAvgSfcVelMerid(ICell) = updateAverage( - LocAvgSfcVelMerid(ICell), ConstSfcVelocity, NAccumSteps); + LocAvgSfcVelMerid(ICell) = + updateAverage(LocAvgSfcVelMerid(ICell), + LocVelocityMerid(ICell, KSfc), NAccumSteps); }); } diff --git a/components/omega/test/ocn/SfcCouplingTest.cpp b/components/omega/test/ocn/SfcCouplingTest.cpp index 5b389a21b6bb..532973d40607 100644 --- a/components/omega/test/ocn/SfcCouplingTest.cpp +++ b/components/omega/test/ocn/SfcCouplingTest.cpp @@ -1,4 +1,5 @@ #include "SfcCoupling.h" +#include "AuxiliaryState.h" #include "Config.h" #include "DataTypes.h" #include "Decomp.h" @@ -16,6 +17,7 @@ #include "OmegaKokkos.h" #include "Pacer.h" #include "TimeStepper.h" +#include "VertAdv.h" #include "VertCoord.h" #include "mpi.h" @@ -134,6 +136,11 @@ int initSfcCouplingTest(const std::string &MeshFile) { Forcing::init(); Tracers::init(); + // Needed by SfcCoupling, which exports the surface velocity + // reconstructed at cell centers by the auxiliary state + VertAdv::init(); + AuxiliaryState::init(); + return Err; } @@ -285,6 +292,10 @@ int testUpdateExportFields(const I4 NSteps) { } Tracers::copyToDevice(0); + // Mirror the run sequence: the reconstruction is refreshed each step + // before the export fields sample it + AuxiliaryState::getDefault()->computeVelocityRecon(DefState, 0); + DefCoupling->updateExportFields(DefState, Tracers::getAll(0)); ModelClock->advance(); @@ -379,6 +390,8 @@ int testExportToCoupler(const CouplingLayout Layout) { int TempIdx = CouplingParams.ExportIdxMap.at("So_t"); int SalinIdx = CouplingParams.ExportIdxMap.at("So_s"); int SshIdx = CouplingParams.ExportIdxMap.at("So_ssh"); + int VelUIdx = CouplingParams.ExportIdxMap.at("So_u"); + int VelVIdx = CouplingParams.ExportIdxMap.at("So_v"); DefCoupling->attachData(CplToOcnData.data(), OcnToCplData.data()); @@ -405,6 +418,25 @@ int testExportToCoupler(const CouplingLayout Layout) { } Tracers::copyToDevice(0); + // Give the edge-normal velocity a nonzero, cell-to-cell varying value + // and reconstruct it at cell centers, so the exported surface velocity + // is a real reconstruction. Halo edges are set here too, since the + // reconstruction stencil for an owned cell reaches beyond that cell. + AuxiliaryState *DefAuxState = AuxiliaryState::getDefault(); + HostArray2DReal NormalVelH = DefState->getNormalVelocityH(0); + for (int Edge = 0; Edge < NormalVelH.extent_int(0); Edge++) { + for (int K = 0; K < NormalVelH.extent_int(1); K++) { + NormalVelH(Edge, K) = 0.01_Real * static_cast((Edge % 7) + 1); + } + } + DefState->copyToDevice(0); + DefAuxState->computeVelocityRecon(DefState, 0); + + auto VelocityZonalH = + createHostMirrorCopy(DefAuxState->VelocityReconAux.VelocityZonalCell); + auto VelocityMeridH = createHostMirrorCopy( + DefAuxState->VelocityReconAux.VelocityMeridionalCell); + DefCoupling->updateExportFields(DefState, Tracers::getAll(0)); auto SshCellOwned = @@ -413,13 +445,22 @@ int testExportToCoupler(const CouplingLayout Layout) { DefCoupling->exportToCoupler(); - // Check 1: exportToCoupler properly packs into OcnToCplView. Velocity - // is skipped here: its averaging is a hardcoded stub pending real vector - // reconstruction (see OcnToCplFields::updateAverages), not yet - // meaningful to check. + // Check 1: exportToCoupler properly packs into OcnToCplView. At + // NAccumSteps == 0 the running average is exactly the single sampled + // value, so the exported velocity should match the reconstruction in + // the surface layer exactly. // copyToHost() converts temp to Kelvin (identity CT->PT w/ ConstantEos) int PackErr = 0; for (int Cell = 0; Cell < NCells; Cell++) { + int KSfc = DefVertCoord->MinLayerCellH(Cell); + if (OcnToCplData[flatIdx(Layout, Cell, VelUIdx, NCells, NExports)] != + VelocityZonalH(Cell, KSfc)) { + PackErr++; + } + if (OcnToCplData[flatIdx(Layout, Cell, VelVIdx, NCells, NExports)] != + VelocityMeridH(Cell, KSfc)) { + PackErr++; + } if (OcnToCplData[flatIdx(Layout, Cell, TempIdx, NCells, NExports)] != ExpectedTemp(Cell) + TkFrz) { PackErr++; @@ -524,6 +565,8 @@ int testEraseAndGet() { void finalizeSfcCouplingTest() { + AuxiliaryState::clear(); + VertAdv::clear(); Tracers::clear(); Forcing::clear(); OceanState::clear(); From 201799f24b0779e5775a6334cc9c56dd1d1fcd49 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 25 Aug 2026 13:55:47 -0500 Subject: [PATCH 6/7] Document the reconstructed velocity components Add VelocityZonalCell and VelocityMeridionalCell to the auxiliary variable list in the user guide, with a section on the three ways they differ from the other auxiliary variables: they are diagnostic and so computed once per step, they are only computed when a stream asks for them or Omega is coupled, and they need a mesh file that supplies the reconstruction stencil and weights. Document computeVelocityRecon and requireVelocityRecon in the auxiliary state developer guide, and IOStream::isFieldRequested in the IOStreams developer guide. Co-Authored-By: Claude Opus 5 --- .../omega/doc/devGuide/AuxiliaryState.md | 16 ++++++++++ components/omega/doc/devGuide/IOStreams.md | 14 +++++++- .../omega/doc/userGuide/AuxiliaryVariables.md | 32 +++++++++++++++++++ 3 files changed, 61 insertions(+), 1 deletion(-) diff --git a/components/omega/doc/devGuide/AuxiliaryState.md b/components/omega/doc/devGuide/AuxiliaryState.md index 18c3396e31a0..a5793f73b32a 100644 --- a/components/omega/doc/devGuide/AuxiliaryState.md +++ b/components/omega/doc/devGuide/AuxiliaryState.md @@ -41,6 +41,22 @@ given ocean state `State`, an array of tracers `TracerArray`, and time level `Ti AuxState.computeAll(State, TracerArray, TimeLevel); ``` +The reconstructed zonal and meridional velocity components are deliberately +not part of `computeAll`, since no tendency reads them and `computeAll` runs +once per time stepper stage. They are computed once per time step instead: +```c++ +AuxState.computeVelocityRecon(State, TimeLevel); +``` +This call does nothing unless the components are needed, which it decides +by asking `IOStream::isFieldRequested` whether any stream contains them. +A caller that needs them regardless of the streams, as surface coupling +does, can say so once during initialization: +```c++ +AuxState.requireVelocityRecon(); +``` +Computing them requires a mesh whose file supplied the reconstruction +stencil and weights; `computeVelocityRecon` aborts otherwise. + ## Removal of auxiliary states To erase a specific named auxiliary state use `erase` ```c++ diff --git a/components/omega/doc/devGuide/IOStreams.md b/components/omega/doc/devGuide/IOStreams.md index 84a5089b16d8..aa279938b06d 100644 --- a/components/omega/doc/devGuide/IOStreams.md +++ b/components/omega/doc/devGuide/IOStreams.md @@ -37,7 +37,19 @@ and the validation status can be checked with ``` All streams must be validated before use to make sure the Fields have been defined and the relevant data arrays have been attached to Fields and -are available to access. At the end of a simulation, IOStreams must be +are available to access. + +Once the streams have been validated, code that computes an optional +diagnostic can ask whether anything will actually read or write it: +```c++ + bool Requested = IOStream::isFieldRequested(FieldName); +``` +This returns true if the Field is in the contents of any defined stream. +Because validation is what replaces a group name in a stream's contents by +the names of the group's member Fields, the answer is only meaningful after +validation; if any stream is still unvalidated the answer is +conservatively true, so an optional computation is done rather than +skipped. At the end of a simulation, IOStreams must be finalized using ```c++ IOStream::finalize(ModelClock); diff --git a/components/omega/doc/userGuide/AuxiliaryVariables.md b/components/omega/doc/userGuide/AuxiliaryVariables.md index 3a12fd370931..634137dca163 100644 --- a/components/omega/doc/userGuide/AuxiliaryVariables.md +++ b/components/omega/doc/userGuide/AuxiliaryVariables.md @@ -31,9 +31,41 @@ The following auxiliary variables are currently available: | Del2TracersCell | laplacian of tracers on cells | SurfTracerRestoringDiffsCell | surface tracer restoring differences on cells | TracersMonthlySurfClimoCell | monthly climatology values to restore to for surface tracer on cells +| VelocityZonalCell | zonal velocity reconstructed at cell centers +| VelocityMeridionalCell | meridional velocity reconstructed at cell centers ## Kinetic energy on cells In [Ringler et al. (2010)](https://www.sciencedirect.com/science/article/pii/S0021999109006780), the cell-centered kinetic energy ($K_i$) is a geometry-weighted combination of the squared edge-normal velocities surrounding cell ($i$), constructed so that its discrete gradient enters the vector-invariant momentum equation as part of the Bernoulli-gradient term, ($-\nabla(K_i+\Phi_i)$). Although only one velocity component is stored at each C-grid edge, the differently oriented edges collectively represent the two-dimensional velocity magnitude; for uniform flow on an isotropic cell, the construction recovers ($K_i=\tfrac12|\mathbf{u}|^2$). Importantly, this definition is chosen for algebraic compatibility with the edge-based kinetic-energy norm and the discrete continuity equation, enabling the nondissipative momentum terms to conserve total energy to within time-discretization error rather than providing an arbitrary pointwise reconstruction of the velocity magnitude. It should be noted that this is the irrotational kinetic energy and is not the total kinetic energy. In principle, this implies that that $K_i$ could be an underestimate of the total kinetic energy. Even for a regular hexagon in an irrotational constant flow, this formulation underestimates the total kinetic energy by 13-16\% depeneding on the flow orientation. However, it has been shown with MPAS-Ocean at standard resolution (Icos30) that $K_i$ *exceeds* $K = 0.5 (u_{cell}^2 + v_{cell}^2)$, likely due to velocity noise at the grid scale that is filtered in the course of reconstructing velocities at cell-centers. This offers a justification for employing $K_i$ in other terms of the momentum equation such as bottom drag. + +## Reconstructed velocity components on cells + +Omega carries only the edge-normal component of the velocity, so the zonal +and meridional components at cell centers are reconstructed from it using +least-squares weights. These two variables differ from the others above in +three ways. + +They are diagnostic: nothing in the Omega equations reads them, so they are +computed once per model time step rather than once per time stepper stage. +They are also only computed when they are needed, that is when some IO +stream asks for one of them or when Omega is running coupled, since the +surface velocity exported to the coupler is taken from them. They are not +in the contents of any stream by default, so add them to a stream to have +them written: +```yaml + Contents: + - VelocityZonalCell + - VelocityMeridionalCell +``` + +They depend on the mesh file. The reconstruction stencil and weights +(`NEdgesReconOnCell`, `ReconStencilCell` and `ReconWeightsCell`) are +precomputed as a mesh preprocessing step rather than by Omega, so a mesh +file that has not been through that step cannot supply them. Requesting +these variables on such a mesh, or running coupled with one, is an error. + +On a planar mesh the reconstructed vector already lies in the plane of the +mesh, so the two variables hold the Cartesian x and y components instead of +zonal and meridional ones. From 67bd4e93fc3ac5b58232bb5cdc663d92be151c41 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 25 Aug 2026 16:08:04 -0500 Subject: [PATCH 7/7] Point the developer guide at the new planar test mesh PlanarPeriodic48x48.omega_vars.260825.nc is the first planar mesh to carry the vector reconstruction stencil and weights, so the CTests that use a planar mesh need it rather than the earlier datestamps. Co-Authored-By: Claude Opus 5 --- components/omega/doc/devGuide/QuickStart.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/omega/doc/devGuide/QuickStart.md b/components/omega/doc/devGuide/QuickStart.md index 38dd8e2475af..49531f3677d9 100644 --- a/components/omega/doc/devGuide/QuickStart.md +++ b/components/omega/doc/devGuide/QuickStart.md @@ -126,12 +126,12 @@ named files under the `test` directory. Appropriate mesh files can be downloaded from: - [Ocean Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/ocean.QU.240km.omega_vars.260807.nc) - [Global Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/cosine_bell_icos480.omega_vars.260807.nc) -- [Planar Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260720.nc) +- [Planar Mesh](https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260825.nc) ```sh cd test wget -O OmegaMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/ocean.QU.240km.omega_vars.260807.nc wget -O OmegaSphereMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/cosine_bell_icos480.omega_vars.260807.nc -wget -O OmegaPlanarMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260720.nc +wget -O OmegaPlanarMesh.nc https://web.lcrc.anl.gov/public/e3sm/polaris/ocean/omega_ctest/PlanarPeriodic48x48.omega_vars.260825.nc cd .. ```