Add velocity-component reconstruction for I/O

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++

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

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

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);

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 ..
 ```
 

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.

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

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<I4> &CellsOnCellInit,       // cell neighbors for each cell
     std::vector<I4> &EdgesOnCellInit,       // edge IDs for each cell edge
     std::vector<I4> &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,34 +427,43 @@ 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
    IO::destroyDecomp(OnCellDecomp);
    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);
    }

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