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/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/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/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 .. ``` 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. 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/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 9a37a0739d98..1ffe5a1fb6f9 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" @@ -29,6 +30,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 +47,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 +59,7 @@ AuxiliaryState::~AuxiliaryState() { VorticityAux.unregisterFields(); VelocityDel2Aux.unregisterFields(); SurfTracerRestAux.unregisterFields(); + VelocityReconAux.unregisterFields(); TracerAux.unregisterFields(); FieldGroup::destroy(GroupName); @@ -323,6 +327,57 @@ 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 +// called once per time step, and only if some IO stream asks for them. +void AuxiliaryState::computeVelocityRecon(const OceanState *State, + 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 " + "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..6821ea4056a0 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,19 @@ 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, 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); + + // 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, @@ -106,6 +121,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/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/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/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); 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/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 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; 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"); } 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();