From 23dd90b18df217f6d8d4c04cd098f1f4a1ae9089 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 2 Jul 2026 07:22:59 -0500 Subject: [PATCH 01/25] Fix biharmonic velocity mixing at boundary vertices Compute the Laplacian of relative vorticity (Del2RelVortVertex) over the full vertex valid range [MinLayerVertexTop, MaxLayerVertexBot] and clamp each edge contribution to that edge's valid range [MinLayerEdgeTop, MaxLayerEdgeBot], matching VorticityAuxVars::computeVarsOnVertex. Previously Del2RelVortVertex was computed only over the narrower [MinLayerVertexBot, MaxLayerVertexTop] range (layers where every surrounding cell is active) and summed Del2Edge with no per-edge clamp. The biharmonic velocity tendency (VelocityHyperDiffOnEdge) reads Del2RelVortVertex over the edge range up to MaxLayerEdgeTop, which for a deep edge sharing a vertex with a shallower cell exceeds MaxLayerVertexTop. Those boundary-vertex layers were never computed, so on partial-bottom meshes the biharmonic term was under-computed there. This is a latent correctness bug on its own; on the fill-values branch the uncomputed layers hold FillValueReal, so the term instead blew up to ~1e45, corrupting NormalVelocity and, through the flux divergence and vertical advection, PseudoThickness and the tracers (surfaced by the new state validation in DRIVER_TEST). The wider range gives boundary vertices a valid, generally non-zero value from their active edges; inactive edges contribute zero via the per-edge clamp and Del2Edge's zeroed boundary band, so no fill value is read. Co-Authored-By: Claude Opus 4.8 --- components/omega/src/ocn/AuxiliaryState.cpp | 10 ++++--- .../ocn/auxiliaryVars/VelocityDel2AuxVars.cpp | 6 ++-- .../ocn/auxiliaryVars/VelocityDel2AuxVars.h | 30 ++++++++++++++----- 3 files changed, 32 insertions(+), 14 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index c2705a3b9cc6..7ac27c57adf8 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -117,10 +117,8 @@ void AuxiliaryState::computeMomAux(const OceanState *State, OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - OMEGA_SCOPE(MinLayerVertexBot, VCoord->MinLayerVertexBot); OMEGA_SCOPE(MinLayerVertexTop, VCoord->MinLayerVertexTop); OMEGA_SCOPE(MaxLayerVertexBot, VCoord->MaxLayerVertexBot); - OMEGA_SCOPE(MaxLayerVertexTop, VCoord->MaxLayerVertexTop); OMEGA_SCOPE(MinLayerEdgeTop, VCoord->MinLayerEdgeTop); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot); @@ -204,8 +202,12 @@ void AuxiliaryState::computeMomAux(const OceanState *State, parallelForOuter( "vertexAuxState2", {Mesh->NVerticesAll}, KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) { - const int KMin = MinLayerVertexBot(IVertex); - const int KMax = MaxLayerVertexTop(IVertex); + // Del2RelVortVertex is computed over the full vertex valid range + // [MinLayerVertexTop, MaxLayerVertexBot] so that boundary-vertex + // layers read by the biharmonic velocity tendency are valid rather + // than fill values (see VelocityDel2AuxVars::computeVarsOnVertex). + const int KMin = MinLayerVertexTop(IVertex); + const int KMax = MaxLayerVertexBot(IVertex); const int KRange = vertRangeChunked(KMin, KMax); parallelForInner( diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp index ffb8046ad5af..500473d18a7f 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp @@ -24,8 +24,10 @@ VelocityDel2AuxVars::VelocityDel2AuxVars(const std::string &AuxStateSuffix, AreaTriangle(Mesh->AreaTriangle), VertexDegree(Mesh->VertexDegree), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), - MinLayerVertexBot(VCoord->MinLayerVertexBot), - MaxLayerVertexTop(VCoord->MaxLayerVertexTop), + MinLayerEdgeTop(VCoord->MinLayerEdgeTop), + MaxLayerEdgeBot(VCoord->MaxLayerEdgeBot), + MinLayerVertexTop(VCoord->MinLayerVertexTop), + MaxLayerVertexBot(VCoord->MaxLayerVertexBot), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h index b6f05f9c0d8b..3e9e57782536 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h @@ -75,8 +75,17 @@ class VelocityDel2AuxVars { } KOKKOS_FUNCTION void computeVarsOnVertex(int IVertex, int KChunk) const { - const int KStart = chunkStart(KChunk, MinLayerVertexBot(IVertex)); - const int KLen = chunkLength(KChunk, KStart, MaxLayerVertexTop(IVertex)); + // Compute over the full vertex valid range [MinLayerVertexTop, + // MaxLayerVertexBot] so that boundary-vertex layers (where only some + // surrounding cells are active) receive a valid value. Each edge's + // contribution is clamped to that edge's valid range [MinLayerEdgeTop, + // MaxLayerEdgeBot], where Del2Edge has been computed or zeroed; this + // matches VorticityAuxVars::computeVarsOnVertex and avoids reading + // uninitialized (fill-value) layers of Del2Edge for deeper edges. + const int KStartVertex = chunkStart(KChunk, MinLayerVertexTop(IVertex)); + const int KLenVertex = + chunkLength(KChunk, KStartVertex, MaxLayerVertexBot(IVertex)); + const int KEndVertex = KStartVertex + KLenVertex - 1; const Real InvAreaTriangle = 1._Real / AreaTriangle(IVertex); @@ -84,16 +93,19 @@ class VelocityDel2AuxVars { for (int J = 0; J < VertexDegree; ++J) { const int JEdge = EdgesOnVertex(IVertex, J); - for (int KVec = 0; KVec < KLen; ++KVec) { - const int K = KStart + KVec; + const int KStartEdge = + Kokkos::max(KStartVertex, MinLayerEdgeTop(JEdge)); + const int KEndEdge = Kokkos::min(KEndVertex, MaxLayerEdgeBot(JEdge)); + for (int K = KStartEdge; K <= KEndEdge; ++K) { + const int KVec = K - KStartVertex; Del2RelVortVertexTmp[KVec] += InvAreaTriangle * DcEdge(JEdge) * EdgeSignOnVertex(IVertex, J) * Del2Edge(JEdge, K); } } - for (int KVec = 0; KVec < KLen; ++KVec) { - const int K = KStart + KVec; + for (int KVec = 0; KVec < KLenVertex; ++KVec) { + const int K = KStartVertex + KVec; Del2RelVortVertex(IVertex, K) = Del2RelVortVertexTmp[KVec]; } } @@ -118,8 +130,10 @@ class VelocityDel2AuxVars { I4 VertexDegree; Array1DI4 MinLayerEdgeBot; Array1DI4 MaxLayerEdgeTop; - Array1DI4 MinLayerVertexBot; - Array1DI4 MaxLayerVertexTop; + Array1DI4 MinLayerEdgeTop; + Array1DI4 MaxLayerEdgeBot; + Array1DI4 MinLayerVertexTop; + Array1DI4 MaxLayerVertexBot; Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; }; From fbaf74184539df049307b583469b1619c9d1174e Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 05:04:21 -0500 Subject: [PATCH 02/25] Add a design document for fixing fill values --- components/omega/doc/design/FillValues.md | 442 ++++++++++++++++++++++ components/omega/doc/index.md | 1 + 2 files changed, 443 insertions(+) create mode 100644 components/omega/doc/design/FillValues.md diff --git a/components/omega/doc/design/FillValues.md b/components/omega/doc/design/FillValues.md new file mode 100644 index 000000000000..f5395438d4d6 --- /dev/null +++ b/components/omega/doc/design/FillValues.md @@ -0,0 +1,442 @@ +(omega-design-fill-values)= +# Fill Values + +**Table of Contents** +1. [Overview](#1-overview) +2. [Requirements](#2-requirements) +3. [Algorithmic Formulation](#3-algorithmic-formulation) +4. [Design](#4-design) +5. [Verification and Testing](#5-verification-and-testing) + +## 1 Overview + +Omega supports output of ocean fields to NetCDF files via I/O streams. Each field +carries a declared fill value stored as metadata under the `_FillValue` attribute, +which marks entries that are undefined or outside the active domain. Currently, fill +value constants are declared locally in each module with inconsistent values +(typically `-9.99e30` for real-valued fields, `-999` for integers). These values do +not match the NetCDF-C standard fill values that analysis tools and visualization +software commonly expect. Additionally, the Kokkos arrays backing these fields are +never explicitly initialized with fill values: inactive ocean layers (model layers +below `MaxLayerCell`) hold uninitialized memory rather than a well-defined sentinel, +and the distinction between "no data" and "zero" is lost. + +This design standardizes fill values across all Omega fields to match the NetCDF-C +standard, ensures that field arrays are automatically initialized with the correct +fill value when data is attached to a `Field`, and verifies through CTests that the +resulting behavior is correct for inactive layers and for edge/vertex fields at +domain boundaries. + +## 2 Requirements + +### 2.1 Requirement: Centralized fill value constants + +Fill value constants must be defined in a single location in Omega that is accessible +to all modules. Scatter of locally-declared constants has produced inconsistencies +(e.g., `-9.99e30` vs `-9.99E+30`) and makes future maintenance error-prone. + +### 2.2 Requirement: Standard fill values matching NetCDF-C + +Fill value constants must exactly match the NetCDF-C standard values (`NC_FILL_*`) +as exposed by the SCORPIO library (`PIO_FILL_*` in `pio.h`). Using standard values +ensures compatibility with NetCDF-aware tools (ncview, Xarray, NCO, etc.) that +recognize and handle fill values automatically. The required values are: + +| Type | Omega name | Value | +|------|-----------------|-------------------------------| +| I4 | `FillValueI4` | -2147483647 (`PIO_FILL_INT`) | +| I8 | `FillValueI8` | `PIO_FILL_INT64` | +| R4 | `FillValueR4` | 9.9692099683868690e+36f (`PIO_FILL_FLOAT`) | +| R8 | `FillValueR8` | 9.9692099683868690e+36 (`PIO_FILL_DOUBLE`) | +| Real | `FillValueReal` | Alias for `FillValueR4` or `FillValueR8` depending on build | + +### 2.3 Requirement: Automatic array initialization at attach time + +When `Field::attachData()` is called, the attached Kokkos array must be +automatically initialized with the standard fill value for the array's element +type before any compute routine runs. The fill value is deduced from +`typename T::value_type` at `attachData()` time using the constants in +`FillValues.h` (`FillValueI4`, `FillValueI8`, `FillValueR4`, or `FillValueR8`), +so the caller does not need to specify it. This guarantees that inactive ocean +layers (below `MaxLayerCell`) and any other undefined entries contain a +well-defined sentinel in output, without requiring each module to manually +perform the fill. + +### 2.4 Requirement: No fill value argument in Field::create() + +The `FillValue` argument is removed from `Field::create()`. Because the fill +value is always one of the standard constants determined by the array element +type, specifying it at field creation time is redundant. The fill value and the +`_FillValue` metadata entry are set automatically when `attachData()` is called. + +### 2.5 Requirement: Inactive layers contain fill values in output + +After compute routines execute, model layers with index `k > MaxLayerCell(ICell)` +must contain the declared fill value in output. With automatic initialization at +attach time (Requirement 2.3) and compute routines that only write to active layers, +this is satisfied without changes to individual compute functions. + +### 2.6 Desired: Valid values at domain boundary edges and vertices + +Edge and vertex fields at valid boundaries adjacent to land or to bathymetry must +contain computed valid values, not fill values. Switching from zero-initialization +to fill-value initialization must not silently break any boundary computation that +previously assumed zero defaults. This mirrors a problem encountered in MPAS-Ocean +when initializing `NormalVelocity` with a fill value rather than zero: code that +implicitly assumed zero at land boundaries produced incorrect results. During +implementation, all compute routines for edge/vertex fields must be audited to +confirm that every valid boundary entry receives an explicit computed value. As an +example of what such an audit can reveal: `TracerHorzAdvOnCell` accumulated its +internal working array `HighOrderFlxHorz` over a cell's full active layer range +without clamping to each contributing edge's active range, reading fill values from +edges shallower than the cell — see Section 4.5 for the resolution. + +### 2.7 Requirement: CTests + +CTests must verify Requirements 2.1–2.3, 2.5, and Desired Requirement 2.6 for +`NormalVelocity`. + +## 3 Algorithmic Formulation + +No new numerical algorithm is introduced. The key behavioral change is in +`Field::attachData()`, which gains a fill step equivalent to: + +```c++ +Kokkos::deep_copy(InDataArray, fill_scalar); +``` + +where `fill_scalar` is the typed fill value extracted from the field's `FieldMeta` +map (stored as `std::any` under key `"_FillValue"`), and the element type +`typename T::value_type` of the Kokkos view determines the `std::any_cast` type. +`Kokkos::deep_copy` with a scalar source broadcasts that value to every element of +the view and dispatches correctly for both host and device memory spaces. + +## 4 Design + +### 4.1 Data types and parameters + +#### 4.1.1 Constants: FillValue definitions in FillValues.h + +The fill value constants are defined in a new lightweight header +`components/omega/src/base/FillValues.h`. The values match the NetCDF-C `NC_FILL_*` +constants (which are identical to the SCORPIO `PIO_FILL_*` constants), but are +written as numeric literals so that the header can be included from any context +without pulling in `pio.h` or `mpi.h`, which carry strict include-order requirements. + +A type-indexed variable template `FillValue` is the primary definition. Its +primary template is intentionally left undefined so that instantiation with an +unsupported type produces a link error. Explicit specializations are provided for +each supported scalar type: + +```c++ +template +constexpr T FillValue; // primary template intentionally undefined + +template <> inline constexpr I4 FillValue = -2147483647; // NC_FILL_INT +template <> inline constexpr I8 FillValue = -9223372036854775806LL; // NC_FILL_INT64 +template <> inline constexpr R4 FillValue = 9.9692099683868690e+36f; // NC_FILL_FLOAT +template <> inline constexpr R8 FillValue = 9.9692099683868690e+36; // NC_FILL_DOUBLE +``` + +Named aliases are provided for readability in comparison and test code: + +```c++ +constexpr I4 FillValueI4 = FillValue; +constexpr I8 FillValueI8 = FillValue; +constexpr R4 FillValueR4 = FillValue; +constexpr R8 FillValueR8 = FillValue; +#if defined(SINGLE_PRECISION) +constexpr Real FillValueReal = FillValue; +#else +constexpr Real FillValueReal = FillValue; +#endif +``` + +`FillValues.h` is included from `Field.h`. Because every module that creates or +uses fields already includes `Field.h` (directly or transitively), no additional +include directives are needed in module source files. + +#### 4.1.2 Class/struct changes + +No new classes or structs are introduced. The changes are confined to: +- new header `FillValues.h` with constants (Section 4.1.1) +- `FillValues.h` included from `Field.h` +- the `attachData` template methods in `Field.h` / `Field.cpp` (Section 4.2) +- removal of local fill value declarations in module files (Section 4.3) + +### 4.2 Methods + +#### 4.2.1 Field::attachData() auto-fill + +The `attachData()` instance method and the `attachFieldData()` static helper +in `components/omega/src/infra/Field.h` gain a private fill step controlled by an +optional `FillOnAttach` parameter (default `true`): + +```c++ +template +void attachData(const T &InDataArray, bool FillOnAttach = true) { + OMEGA_ASSERT(isKokkosArray, + "Field::attachData requires Kokkos array as input"); + DataArray = std::make_shared(InDataArray); + DataType = checkArrayType(); + MemLoc = findArrayMemLoc(); + if (FillOnAttach) + fillWithValue(InDataArray); // auto-fill with type-deduced FillValue +} + +template +static void attachFieldData(const std::string &FieldName, + const T &InDataArray, bool FillOnAttach = true); +``` + +The default `FillOnAttach = true` preserves the fill behavior for the normal +initialization workflow (fresh arrays attached before any data is written). Pass +`FillOnAttach = false` when re-attaching an existing data-filled array where the +fill must not overwrite computed values — see Section 4.2.3 for the canonical use case. + +The private helper `fillWithValue()` performs the following steps: + +1. Derives the element type as `typename T::non_const_value_type` (using the + non-const form so that const view types resolve to the same underlying type). +2. Looks up the fill value via `FillValue`, a constexpr variable template + specialization from `FillValues.h`. Supporting a new type requires only adding + one specialization there; no changes to `fillWithValue` are needed. +3. Stores the typed fill value in `FieldMeta["_FillValue"]` (the NetCDF/CF + standard attribute name), making it available to IO and metadata queries. +4. Calls `Kokkos::deep_copy(InDataArray, scalar)` to broadcast the scalar to every + element of the view. This works for host views (`HostArray*`), device views + (`Array*`), and all dimensionalities (1D–5D) without additional dispatch. +5. Only the view being attached is filled; no assumption is made about whether the + other memory space (host or device) has been attached. + +Ordering invariants preserved by this design: +- `Field::create()` is always called before `attachData()`, so the field exists + in `AllFields` before fill time. +- The `"_FillValue"` metadata entry is set at `attachData()` time, which precedes + any IO write operation in the initialization workflow. +- Restart reads happen after `attachData()` and overwrite fill entries with data + from the restart file. +- Compute routines run after initialization and overwrite active entries, leaving + inactive entries (below `MaxLayerCell`) at the fill value. + +#### 4.2.2 No per-module manual initialization required + +Because `attachData()` handles initialization automatically, no module needs to +call an explicit fill routine. Modules may still call `Kokkos::deep_copy` or +similar to set specific values, but this is independent of fill-value management. + +#### 4.2.3 Ordering invariant and anti-patterns + +**Invariant**: `attachData()` with `FillOnAttach = true` (the default) must always +be called *before* any real values are written to the array. Any data written before +such a call is silently overwritten by the auto-fill step and is dead code. + +Three anti-patterns discovered and corrected during implementation: + +**Class A — Wasted pre-fill**: A `deepCopy` (or equivalent write) is applied to an +array before `defineFields()` calls `attachData()` on it. Because `attachData()` +overwrites the entire array, the earlier write has no effect. + +```c++ +// WRONG: initial value is overwritten by attachData() +deepCopy(SpecVol, 1.0_Real / RhoSw); +defineFields(); // calls attachData(), fills SpecVol with FillValueReal +``` + +**Class B — Compute-time re-initialization**: A `deepCopy(arr, 0)` (or any scalar) +at the top of a recurring compute method resets the entire array on every timestep, +overwriting the correct fill values in inactive layers with zero. + +```c++ +// WRONG: destroys fill values in inactive layers every timestep +void Foo::compute(...) { + deepCopy(LocArr, 0); + parallelFor({NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + // writes only active layers KMin..KMax + }); +} +``` + +The correct pattern for Classes A and B is to do nothing: `attachData()` sets the +fill value once at initialization, and compute kernels that iterate only over active +layers naturally leave inactive layers at the fill value. + +**Class C — Time-level pointer update**: `attachData()` is called each timestep to +re-point the IO field to the newly-current time-level array (e.g., in +`OceanState::updateTimeLevels()` and `Tracers::updateTimeLevels()`). The array +already contains the just-computed state. Calling `attachData()` with the default +`FillOnAttach = true` here would silently destroy the computed values before IO +writes them. This is the most dangerous anti-pattern because the model runs to +completion without aborting, but all output fields contain fill values. + +```c++ +// WRONG: destroys computed state before IO can write it +void OceanState::updateTimeLevels() { + CurTimeIndex = (CurTimeIndex + 1) % NTimeLevels; + Field::attachFieldData(NormalVelocityFldName, + NormalVelocity[CurTimeIndex]); // fills! +} + +// CORRECT: pointer update only, computed values preserved +void OceanState::updateTimeLevels() { + CurTimeIndex = (CurTimeIndex + 1) % NTimeLevels; + Field::attachFieldData(NormalVelocityFldName, + NormalVelocity[CurTimeIndex], false); +} +``` + +A related instance occurs in the IO write path when the time-coordinate array is +populated before `attachFieldData("time", OutTime)` is called. The correct pattern +is to attach first (which fills the single-element array with the fill value), then +assign the elapsed time to `OutTime(0)`. Because Kokkos host views share backing +memory, the assignment is visible through the stored pointer before the field is +written to the file. + +### 4.3 Module-level changes: removing FillValue from Field::create() calls + +The `FillValue` parameter is removed from `Field::create()`. All ~50 call sites +across Omega drop the fill value argument. The pattern is: + +```c++ +// BEFORE — explicit fill value argument +Field::create(FieldName, Description, Units, StdName, + ValidMin, ValidMax, FillValueReal, NDims, DimNames); + +// AFTER — fill value is deduced automatically from the array type at attachData() +Field::create(FieldName, Description, Units, StdName, + ValidMin, ValidMax, NDims, DimNames); +``` + +The `Tracers::define()` helper in `src/ocn/Tracers.h` / `Tracers.cpp` and +`TracerDefs.inc` also have their `FillValue` parameter and arguments removed, +since `Tracers::define()` calls `Field::create()` internally. + +### 4.4 Boundary edge and vertex handling + +During implementation, all compute routines that populate fields defined on mesh +edges or vertices must be audited to confirm that every valid boundary entry — +including edges and vertices adjacent to land cells and edges at the base of the +water column — receives an explicit computed value rather than relying on +zero-initialization. Where such reliance is found, an explicit assignment must be +added. The CTest in Section 5.5 detects regressions after the change. + +### 4.5 Edge flux field masking + +Edge fields that represent fluxes normal to the edge face — such as +`NormalVelocityTend`, `VelocityDel2Aux.Del2Edge`, and `NormalVelocity` itself — +must carry zeros (not fill values) in the boundary layer sub-ranges +`[MinLayerEdgeTop, MinLayerEdgeBot)` and `(MaxLayerEdgeTop, MaxLayerEdgeBot]`, +where one neighboring cell is active but the other is not. Two helper methods +are added to `VertCoord` (`components/omega/src/ocn/VertCoord.h/.cpp`): + +**`zeroEdgeField(Array2DReal &Arr, I4 NEdgesAll)`** — zeros all layers in +`[MinLayerEdgeTop, MaxLayerEdgeBot]`. Used before recomputing a flux-type edge +field each time step so that boundary edges show 0, not `FillValueReal`. Layers +outside the valid range retain their fill value from `attachData()`. Called for: +- `NormalVelocityTend` at the start of `Tendencies::computeVelocityTendenciesOnly()` +- `VelocityDel2Aux.Del2Edge` before the `edgeAuxState2` kernel in + `AuxiliaryState::computeMomAux()`, because the Del4 hyperdiffusion kernel + accumulates `Del2Edge` over all edges sharing a cell or vertex; fill values in + boundary layers of a neighbouring edge would corrupt `Del2DivCell` and + `Del2RelVortVertex`. + +**`applyEdgeLayerMask(Array2DReal &Arr, I4 NEdgesAll)`** — enforces the full +three-zone mask on an edge field after IC or restart read: layers outside +`[MinLayerEdgeTop, MaxLayerEdgeBot]` are set to `FillValueReal`; layers inside +`[MinLayerEdgeTop, MaxLayerEdgeBot]` but outside the active range +`[MinLayerEdgeBot, MaxLayerEdgeTop]` are set to 0; active layers are left +unchanged. Called in `ocnInit` (`OceanInit.cpp`) for +`NormalVelocity[CurTimeLevel]` after `exchangeHalo` and before `copyToHost`. +`NormalVelocity` is the only edge-based state field read from IC/restart; all +other IC fields (pseudo-thickness, tracers, VertCoord fields) are cell-based +and need no edge masking. + +Fields that are *not* fluxes — `MeanPseudoThickEdge`, `FluxPseudoThickEdge` — +are genuinely undefined at boundary edges (interpolating thickness when one +neighboring cell is land has no physical meaning) and therefore correctly retain +`FillValueReal` in those layers; no zeroing is applied to them. +`VorticityAux` on edges is already computed over the full valid range +`[MinLayerEdgeTop, MaxLayerEdgeBot]`, so it requires no additional masking. + +**`HighOrderFlxHorz` in `TracerHorzAdvOnCell`** — the 3-D working array +`(NTracers, NEdgesSize, NVertLayers)` that accumulates horizontal tracer flux +contributions per edge requires separate treatment. The edge-pass overload of +`TracerHorzAdvOnCell::operator()` writes `HighOrderFlxHorz(L, IEdge, K)` only +for layers `K ∈ [0, MaxLayerEdgeTop(IEdge)]`; layers beyond that index retain +their fill values from `attachData()`. The cell-pass overload then reads +`HighOrderFlxHorz` over the cell's own active range `[MinLayerCell, MaxLayerCell]`, +which for cells adjacent to shallower edges extends past `MaxLayerEdgeTop` of +those edges — reading fill values and corrupting active-layer tracer tendencies. +The fix is to add `MinLayerEdgeBot` and `MaxLayerEdgeTop` as members of +`TracerHorzAdvOnCell` and clamp the cell-pass inner loop to +`[max(KStart, MinLayerEdgeBot(IEdge)), min(KEnd, MaxLayerEdgeTop(IEdge)+1))`. +This follows the identical pattern already used by `PseudoThicknessFluxDivOnCell`, +`TracerDiffOnCell`, `TracerAuxVars::computeVarsOnCells`, and every other +cell-level flux-accumulation kernel, and avoids the need to zero +`HighOrderFlxHorz` before the edge pass each timestep. + +## 5 Verification and Testing + +A new test `test/infra/FillValueTest.cpp` is added and registered as +`FILL_VALUE_TEST` with 8 MPI tasks using `add_omega_test()` in the test +`CMakeLists.txt`. The test follows the standard Omega test pattern: initialize +`MachEnv`, logging, IO, `Decomp`, and `Dimension` before running assertions. + +### 5.1 Test: fill constant values + +Verify that each Omega fill value constant exactly equals its NetCDF-C counterpart: + +```c++ +TstEval("FillValueI4 == NC_FILL_INT", FillValueI4, (I4)NC_FILL_INT, Err); +TstEval("FillValueI8 == NC_FILL_INT64", FillValueI8, (I8)NC_FILL_INT64, Err); +TstEval("FillValueR4 == NC_FILL_FLOAT", FillValueR4, (R4)NC_FILL_FLOAT, Err); +TstEval("FillValueR8 == NC_FILL_DOUBLE",FillValueR8, (R8)NC_FILL_DOUBLE, Err); +``` + +Tests requirements: 2.1, 2.2, 2.7 + +### 5.2 Test: attachData auto-fill + +Create a `Field`, allocate a Kokkos host or device array whose elements are all set +to a distinct sentinel (e.g., 0), then call `attachData()`. Verify that after the +call every element of the array equals `FillValueR8` (or the appropriate type). +Repeat for I4, R4, and R8 types. + +Tests requirements: 2.3, 2.7 + +### 5.3 Test: inactive layers contain fill values after compute + +Using the smallest test mesh (the `planar` mesh used by `VertCoordTest`), +run `VertCoord::computePressure()` or `computeGeomZHeight()`. For each cell, verify +that all entries at depth index `k > MaxLayerCell(ICell)` equal `FillValueR8`. + +Tests requirements: 2.5, 2.7 + +### 5.4 Test: NetCDF output contains correct fill value attribute + +Write a field to a test NetCDF file via `IOStream`. Read back the file's `_FillValue` +attribute for the variable and verify it equals `NC_FILL_DOUBLE` (or `NC_FILL_FLOAT` +for reduced-precision fields). Also read back the array entries for inactive layers +and verify they equal the fill value. + +Tests requirements: 2.5, 2.7 + +### 5.5 Test: NormalVelocity three-zone fill/zero/real pattern + +After full state initialization, IC read, halo exchange, and +`VertCoord::applyEdgeLayerMask()`, verify that `NormalVelocity` satisfies the +three-zone invariant for every owned edge: + +- **Fully inactive layers** (outside `[MinLayerEdgeTop, MaxLayerEdgeBot]`): + equal `FillValueReal`. +- **Boundary layers** (inside `[MinLayerEdgeTop, MaxLayerEdgeBot]` but outside + `[MinLayerEdgeBot, MaxLayerEdgeTop]`): equal `0`. +- **Active layers** (`[MinLayerEdgeBot, MaxLayerEdgeTop]`): not equal to + `FillValueReal` (for a quiescent IC the value is `0`, which is valid). + +This guards against two regressions: (a) the original zero-initialization +assumption (boundary layers silently filled with zero before `attachData()` +was changed) and (b) IC files that store zeros in inactive layers, which +`applyEdgeLayerMask` must replace with `FillValueReal`. + +Tests requirements: 2.5, 2.6, 2.7 diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 8d1b38501df8..eb8504db8047 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -115,6 +115,7 @@ design/Decomp design/Driver design/EOS design/Error +design/FillValues design/Halo design/HorzMeshClass design/Logging From 985b8e6b7e23592f46922835c5a3314445853640 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 05:39:36 -0500 Subject: [PATCH 03/25] Add fill value constant This merge adds 5 constexpr fill value constants (FillValueI4, FillValueI8, FillValueR4, FillValueR8, FillValueReal) in namespace OMEGA, wrapping PIO_FILL_* from SCORPIO. --- components/omega/src/base/FillValues.h | 32 ++++++++++++++++++++++++++ components/omega/src/base/IO.h | 1 + 2 files changed, 33 insertions(+) create mode 100644 components/omega/src/base/FillValues.h diff --git a/components/omega/src/base/FillValues.h b/components/omega/src/base/FillValues.h new file mode 100644 index 000000000000..1db933d96d3e --- /dev/null +++ b/components/omega/src/base/FillValues.h @@ -0,0 +1,32 @@ +#ifndef OMEGA_FILL_VALUES_H +#define OMEGA_FILL_VALUES_H +//===-- base/FillValues.h - standard fill value constants -------*- C++ -*-===// +/// +/// \file +/// \brief Standard fill value constants matching NetCDF-C NC_FILL_* values. +/// +/// These values exactly match the SCORPIO PIO_FILL_* constants (which wrap +/// NC_FILL_* from netcdf.h). Defined here as literals so this header can be +/// included without pulling in pio.h or mpi.h, which have strict include-order +/// requirements that make them unsafe to include from generic headers. +// +//===----------------------------------------------------------------------===// + +#include "DataTypes.h" + +namespace OMEGA { + +constexpr I4 FillValueI4 = -2147483647; ///< NC_FILL_INT +constexpr I8 FillValueI8 = -9223372036854775806LL; ///< NC_FILL_INT64 +constexpr R4 FillValueR4 = 9.9692099683868690e+36f; ///< NC_FILL_FLOAT +constexpr R8 FillValueR8 = 9.9692099683868690e+36; ///< NC_FILL_DOUBLE +#if defined(SINGLE_PRECISION) +constexpr Real FillValueReal = FillValueR4; +#else +constexpr Real FillValueReal = FillValueR8; +#endif + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// +#endif // OMEGA_FILL_VALUES_H diff --git a/components/omega/src/base/IO.h b/components/omega/src/base/IO.h index a8ec08e73f4a..4aa606558500 100644 --- a/components/omega/src/base/IO.h +++ b/components/omega/src/base/IO.h @@ -364,6 +364,7 @@ void writeNDVar(void *Variable, ///< [in] variable to be written ); } // end namespace IO + } // end namespace OMEGA //===----------------------------------------------------------------------===// From e96a067ac3b6c93c26c9e68848d9a1b65f652f83 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 05:41:51 -0500 Subject: [PATCH 04/25] Fill fields with fill values during attachData() This merge adds private fillWithValue() template helper (using if constexpr on T::value_type) and a call to it at the end of attachData(). Arrays are now auto-filled with the declared fill value at attach time. --- components/omega/src/infra/Field.h | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/components/omega/src/infra/Field.h b/components/omega/src/infra/Field.h index 313883184edc..3ab67b54d579 100644 --- a/components/omega/src/infra/Field.h +++ b/components/omega/src/infra/Field.h @@ -19,6 +19,7 @@ #include "DataTypes.h" #include "Dimension.h" #include "Error.h" +#include "FillValues.h" #include "Logging.h" #include "OmegaKokkos.h" #include "TimeMgr.h" @@ -26,6 +27,7 @@ #include #include #include +#include namespace OMEGA { @@ -82,6 +84,25 @@ class Field { /// various types and cast to the appropriate type when needed. std::shared_ptr DataArray; + /// Fills every element of InDataArray with the field's declared fill value. + /// Called automatically by attachData() so inactive entries are initialized + /// to a well-defined sentinel before any compute routine runs. + template void fillWithValue(const T &InDataArray) { + using ValType = typename T::value_type; + auto AnyFill = FieldMeta.at("FillValue"); + ValType FillVal; + if constexpr (std::is_same_v) + FillVal = std::any_cast(AnyFill); + else if constexpr (std::is_same_v) + FillVal = std::any_cast(AnyFill); + else if constexpr (std::is_same_v) + // Fill values are stored as R8 (Real) in legacy code; cast safely + FillVal = static_cast(std::any_cast(AnyFill)); + else + FillVal = std::any_cast(AnyFill); + Kokkos::deep_copy(InDataArray, FillVal); + } + public: //--------------------------------------------------------------------------- // Initialization @@ -314,6 +335,10 @@ class Field { // Determine type and location DataType = checkArrayType(); MemLoc = findArrayMemLoc(); + + // Initialize every element to the declared fill value so inactive + // entries (e.g. layers below MaxLayerCell) are well-defined in output + fillWithValue(InDataArray); }; //--------------------------------------------------------------------------- From 0a7aec591689716e10537b890b70a615bcca782b Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 05:44:48 -0500 Subject: [PATCH 05/25] Update fill values in all fields This merge removes all local fill value declarations and replaces references with the centralized constants. --- components/omega/src/infra/Field.cpp | 2 +- components/omega/src/ocn/Eos.cpp | 32 +-- components/omega/src/ocn/HorzMesh.cpp | 215 +++++++++--------- components/omega/src/ocn/OceanState.cpp | 8 +- components/omega/src/ocn/Tendencies.cpp | 6 +- components/omega/src/ocn/TendencyTerms.cpp | 4 +- components/omega/src/ocn/TracerDefs.inc | 10 +- components/omega/src/ocn/VertAdv.cpp | 3 +- components/omega/src/ocn/VertCoord.cpp | 29 +-- components/omega/src/ocn/VertMix.cpp | 15 +- .../src/ocn/auxiliaryVars/KineticAuxVars.cpp | 11 +- .../auxiliaryVars/PseudoThicknessAuxVars.cpp | 9 +- .../auxiliaryVars/SurfTracerRestAuxVars.cpp | 9 +- .../src/ocn/auxiliaryVars/TracerAuxVars.cpp | 5 +- .../ocn/auxiliaryVars/VelocityDel2AuxVars.cpp | 17 +- .../ocn/auxiliaryVars/VorticityAuxVars.cpp | 21 +- 16 files changed, 185 insertions(+), 211 deletions(-) diff --git a/components/omega/src/infra/Field.cpp b/components/omega/src/infra/Field.cpp index 5c24a9efd01d..921ea7885aa0 100644 --- a/components/omega/src/infra/Field.cpp +++ b/components/omega/src/infra/Field.cpp @@ -54,7 +54,7 @@ void Field::init(const Clock *ModelClock // [in] default model clock std::string CalName = CalendarCFName[CalKind]; std::vector DimNamesTmp; // empty dim names vector std::shared_ptr TimeField = - create("time", "time", UnitString, "time", 0.0, 1.e20, -9.99e30, 0, + create("time", "time", UnitString, "time", 0.0, 1.e20, FillValueR8, 0, DimNamesTmp, true, true); TimeField->addMetadata("calendar", CalName); } diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 8101e6bab42f..d2ecf014abb0 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -49,9 +49,6 @@ Eos::Eos(const std::string &Name, ///< [in] Name for eos object BruntVaisalaFreqSq = Array2DReal("BruntVaisalaFreqSq", Mesh->NCellsSize, VCoord->NVertLayersP1); - deepCopy(SpecVol, 1.0_Real / RhoSw); - deepCopy(SpecVolDisplaced, 1.0_Real / RhoSw); - defineFields(); } @@ -142,8 +139,6 @@ void Eos::computeSpecVol(const Array2DReal &ConservTemp, OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - deepCopy(LocSpecVol, 0); /// Initialize local specific volume to zero - I4 KDisp = 0; /// No displacement in this case /// Dispatch to the correct EOS calculation @@ -208,9 +203,6 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp, OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - deepCopy(LocSpecVolDisplaced, - 0); /// Initialize local specific volume to zero - /// Dispatch to the correct EOS calculation /// If EosChoice is Linear, the displaced specific /// volume is the same as the specific volume @@ -273,9 +265,6 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - deepCopy(LocBruntVaisalaFreqSq, - 0); /// Initialize local Brunt-Vaisala frequency to zero - /// Dispatch to the correct squared Brunt-Vaisala frequency calculation if (EosChoice == EosType::LinearEos) { /// If Linear EOS, use linear squared Brunt-Vaisala frequency calculation @@ -368,8 +357,7 @@ void Eos::defineFields() { } /// Create fields for state variables - const Real FillValue = -9.99e30; - int NDims = 2; + int NDims = 2; std::vector DimNames(NDims); DimNames[0] = "NCells"; DimNames[1] = "NVertLayers"; @@ -382,9 +370,9 @@ void Eos::defineFields() { "sea_water_specific_volume", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValue, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + FillValueReal, // Scalar used for undefined entries + NDims, // Number of dimensions + DimNames // Dimension names ); /// Create and register the displaced specific volume field auto SpecVolDisplacedField = @@ -395,9 +383,9 @@ void Eos::defineFields() { "sea_water_specific_volume_displaced", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValue, // Scalar used for undefined entried - NDims, // Number of dimensions - DimNames // Dimension names + FillValueReal, // Scalar used for undefined entries + NDims, // Number of dimensions + DimNames // Dimension names ); // Brunt-Vaisala frequency is located at interfaces @@ -411,9 +399,9 @@ void Eos::defineFields() { "sea_water_brunt_vaisala_frequency_squared", // CF-ish Name std::numeric_limits::min(), // Min valid value std::numeric_limits::max(), // Max valid value - FillValue, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + FillValueReal, // Scalar used for undefined entries + NDims, // Number of dimensions + DimNames // Dimension names ); // Create a field group for the eos-specific state fields diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 9e47b81e350a..6fde75403f2c 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -559,10 +559,10 @@ void HorzMesh::defineMeshFields() { "x", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, XCell); @@ -576,10 +576,10 @@ void HorzMesh::defineMeshFields() { "y", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, YCell); @@ -593,10 +593,10 @@ void HorzMesh::defineMeshFields() { "z", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, ZCell); @@ -610,10 +610,10 @@ void HorzMesh::defineMeshFields() { "latitude", // CF standard Name -3.1415927, // min valid value 3.1415927, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LatCell); @@ -624,13 +624,13 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Longitude coordinates of cell centers", // long Name "radians", // units - "longitude", // CF standard Name - 0.0, // min valid value - 6.28319, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "longitude", // CF standard Name + 0.0, // min valid value + 6.28319, // max valid value + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LonCell); @@ -646,10 +646,10 @@ void HorzMesh::defineMeshFields() { "x", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, XEdge); @@ -663,10 +663,10 @@ void HorzMesh::defineMeshFields() { "y", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, YEdge); @@ -680,10 +680,10 @@ void HorzMesh::defineMeshFields() { "z", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, ZEdge); @@ -697,10 +697,10 @@ void HorzMesh::defineMeshFields() { "latitude", // CF standard Name -3.1415927, // min valid value 3.1415927, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LatEdge); @@ -714,10 +714,10 @@ void HorzMesh::defineMeshFields() { "longitude", // CF standard Name 0.0, // min valid value 6.28319, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LonEdge); @@ -733,10 +733,10 @@ void HorzMesh::defineMeshFields() { "x", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, XVertex); @@ -750,10 +750,10 @@ void HorzMesh::defineMeshFields() { "y", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, YVertex); @@ -767,10 +767,10 @@ void HorzMesh::defineMeshFields() { "z", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, ZVertex); @@ -781,13 +781,13 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Latitude coordinates of cell vertices", // long Name "radians", // units - "latitude", // CF standard Name - -3.1415927, // min valid value - 3.1415927, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "latitude", // CF standard Name + -3.1415927, // min valid value + 3.1415927, // max valid value + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LatVertex); @@ -798,33 +798,34 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Longitude coordinates of cell vertices", // long Name "radians", // units - "longitude", // CF standard Name - 0.0, // min valid value - 6.28319, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "longitude", // CF standard Name + 0.0, // min valid value + 6.28319, // max valid value + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LonVertex); // Other mesh properties // Mesh areas, lengths, and angles - DimNames[0] = "NCells"; - FieldName = "AreaCell"; - AreaCell = Array1DReal("AreaCell", NCellsSize); // allocate space - auto AreaCellField = Field::create(FieldName, // field name - "Area of each cell (m^2)", // long name - "m2", // units - "cell_area", // CF standard name - 0.0, // min valid value - 9.99E+30, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent - ); + DimNames[0] = "NCells"; + FieldName = "AreaCell"; + AreaCell = Array1DReal("AreaCell", NCellsSize); // allocate space + auto AreaCellField = + Field::create(FieldName, // field name + "Area of each cell (m^2)", // long name + "m2", // units + "cell_area", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, AreaCell); @@ -838,7 +839,7 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -853,13 +854,13 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Area of each triangle in the dual grid (m^2)", // lng name "m2", // units - "cell_area", // CF standard name - 0.0, // min valid value - 9.99E+30, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "cell_area", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, AreaTriangle); @@ -875,7 +876,7 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -893,7 +894,7 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -906,15 +907,15 @@ void HorzMesh::defineMeshFields() { auto AngleEdgeField = Field::create(FieldName, // field name "Angle the edge normal makes with local eastward direction" - " (radians)", // long name - "radians", // units - "", // CF standard name - -3.1415927, // min valid value - 3.1415927, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + " (radians)", // long name + "radians", // units + "", // CF standard name + -3.1415927, // min valid value + 3.1415927, // max valid value + FillValueReal, // scalar for undefined entries + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, AngleEdge); @@ -934,7 +935,7 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -955,7 +956,7 @@ void HorzMesh::defineMeshFields() { "", // CF standard name -1.0, // min valid value 1.0, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -976,7 +977,7 @@ void HorzMesh::defineMeshFields() { "coriolis_parameter", // CF standard name -1.0E-3, // min valid value 1.0E-3, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -994,7 +995,7 @@ void HorzMesh::defineMeshFields() { "coriolis_parameter", // CF standard name -1.0E-3, // min valid value 1.0E-3, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -1012,7 +1013,7 @@ void HorzMesh::defineMeshFields() { "coriolis_parameter", // CF standard name -1.0E-3, // min valid value 1.0E-3, // max valid value - -9.99E+30, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index e696d06a5a7e..6ee850207f95 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/components/omega/src/ocn/OceanState.cpp @@ -202,9 +202,9 @@ void OceanState::defineFields() { "sea_water_velocity", // CF standard Name -9.99E+10, // min valid value 9.99E+10, // max valid value - -9.99E+30, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names ); DimNames[0] = "NCells"; @@ -215,7 +215,7 @@ void OceanState::defineFields() { "cell_thickness", // CF standard Name 0.0, // min valid value 9.99E+30, // max valid value - -9.99E+30, // scalar used for undefined entries + FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index ba0df5dd692a..77851168871b 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -338,7 +338,7 @@ void Tendencies::defineFields() { auto PseudoThicknessTendField = Field::create(PseudoThicknessTendFieldName, "Pseudo-thickness tendency", "m/s", "cell_thickness_tendency", -9.99E+10, 9.99E+10, - -9.99E+30, NDims, DimNamesThickness); + FillValueReal, NDims, DimNamesThickness); NDims = 3; std::vector DimNamesTracer(NDims); DimNamesTracer[0] = "NTracers"; @@ -346,7 +346,7 @@ void Tendencies::defineFields() { DimNamesTracer[2] = "NVertLayers"; auto TracerTendField = Field::create( TracerTendFieldName, "Tracer tendency", "kg/m^3/s", "tracer_tendency", - -9.99E+10, 9.99E+10, -9.99E+30, NDims, DimNamesTracer); + -9.99E+10, 9.99E+10, FillValueReal, NDims, DimNamesTracer); NDims = 2; std::vector DimNamesVelocity(NDims); DimNamesVelocity[0] = "NEdges"; @@ -354,7 +354,7 @@ void Tendencies::defineFields() { auto NormalVelocityTendField = Field::create(NormalVelocityTendFieldName, "Normal velocity tendency", "m/s^2", "sea_water_velocity_tendency", -9.99E+10, - 9.99E+10, -9.99E+30, NDims, DimNamesVelocity); + 9.99E+10, FillValueReal, NDims, DimNamesVelocity); std::string TendGroupName = "Tendencies"; if (Name != "Default") { diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 4d7ce499941f..5b142888142c 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -88,9 +88,7 @@ TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh, Mesh->NEdgesAll, VCoord->NVertLayers), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell), - DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) { - deepCopy(HighOrderFlxHorz, 0); -} + DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell) {} TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) diff --git a/components/omega/src/ocn/TracerDefs.inc b/components/omega/src/ocn/TracerDefs.inc index fb2326e9c203..03134e75da5d 100644 --- a/components/omega/src/ocn/TracerDefs.inc +++ b/components/omega/src/ocn/TracerDefs.inc @@ -27,16 +27,16 @@ static void defineAllTracers() { "sea_water_conservative_temperature", ///< [in] CF standard Name -273.15, ///< [in] min valid field value 100.0, ///< [in] max valid field value - 1.e33, ///< [in] value of undef entries + FillValueReal, ///< [in] value of undef entries IndxTemp); ///< [out] (optional) static index define("Salinity", "Absolute Salinity", "g kg-1", - "sea_water_absolute_salinity", 0.0, 50.0, 1.e33, IndxSalt); + "sea_water_absolute_salinity", 0.0, 50.0, FillValueReal, IndxSalt); - define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, 1.e33, + define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, FillValueReal, IndxDebug1); - define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, 1.e33, + define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, FillValueReal, IndxDebug2); - define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, 1.e33, + define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, FillValueReal, IndxDebug3); } diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 04e363a0ad3a..dee2b4e2d90b 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -128,8 +128,7 @@ void VertAdv::defineFields() { LowOrderVertFluxFldName.append(Name); } - const Real FillValueReal = -9.99e30; - I4 NDims = 2; + I4 NDims = 2; std::vector DimNames(NDims); DimNames[0] = "NCells"; DimNames[1] = "NVertLayersP1"; diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 146ceb414efc..536a9bfa9a3d 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -247,9 +247,7 @@ void VertCoord::defineFields() { } // Create fields for VertCoord variables - const I4 FillValueI4 = -999; - const Real FillValueReal = -9.99e30; - int NDims = 1; + int NDims = 1; std::vector DimNames(NDims); DimNames[0] = "NCells"; @@ -524,15 +522,6 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) { // uninitialized and will need to be initialized explicitly if needed. if (ReadStream) { - I4 FillValueI4 = -1; - Real FillValueReal = -999._Real; - - deepCopy(MinLayerCell, FillValueI4); - deepCopy(MaxLayerCell, FillValueI4); - deepCopy(BottomGeomDepth, FillValueReal); - deepCopy(RefPseudoThickness, FillValueReal); - deepCopy(VertCoordMovementWeights, FillValueReal); - // Fetch input stream and validate std::string StreamName = "InitialVertCoord"; if (Name != "Default") { @@ -591,7 +580,8 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) { Accum += LocBottomGeomDepth(I); }, Sum3); - if (Sum3 < 0.) { + if (Sum3 / BottomGeomDepth.extent_int(0) >= + 0.5_Real * FillValueReal) { ABORT_ERROR("VertCoord: Error reading BottomGeomDepth from {}", StreamName); } @@ -603,7 +593,9 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) { Accum += LocPseudoThick(I, J); }, Sum4); - if (Sum4 < 0.) { + I4 N4 = RefPseudoThickness.extent_int(0) * + RefPseudoThickness.extent_int(1); + if (Sum4 / N4 >= 0.5_Real * FillValueReal) { ABORT_ERROR("VertCoord: Error reading RefPseudoThickness " "from {}", StreamName); @@ -620,12 +612,15 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) { } }, Sum5, NumNonZero); - if (Sum5 < 0.) { + I4 N5 = VertCoordMovementWeights.extent_int(0); + if (Sum5 / N5 >= 0.5_Real * FillValueReal) { + // Stream did not populate weights — use default + deepCopy(VertCoordMovementWeights, 1._Real); + } else if (Sum5 < 0.) { ABORT_ERROR("VertCoord: Error reading VertCoordMovementWeights " "from {}", StreamName); - } - if (NumNonZero == 0) { + } else if (NumNonZero == 0) { // TODO: ABORT_ERROR when all weights equal 0 deepCopy(VertCoordMovementWeights, 1._Real); } diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 815dab5e6814..01c923e81b23 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -332,8 +332,7 @@ void VertMix::defineFields() { } /// Create fields for state variables - const Real FillValue = -9.99e30; - int NDims = 2; + int NDims = 2; std::vector DimNames(NDims); DimNames[0] = "NCells"; DimNames[1] = "NVertLayersP1"; @@ -347,9 +346,9 @@ void VertMix::defineFields() { "vertical_diffusivity", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValue, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + FillValueReal, // Scalar used for undefined entries + NDims, // Number of dimensions + DimNames // Dimension names ); /// Create and register the VertVisc field auto VertViscField = @@ -360,9 +359,9 @@ void VertMix::defineFields() { "vertical_viscosity", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValue, // Scalar used for undefined entried - NDims, // Number of dimensions - DimNames // Dimension names + FillValueReal, // Scalar used for undefined entries + NDims, // Number of dimensions + DimNames // Dimension names ); /// Create and register the GradRichNum field auto GradRichNumField = diff --git a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp index 5e5bb08d4714..8c76aabaa922 100644 --- a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp @@ -23,8 +23,7 @@ void KineticAuxVars::registerFields( ) const { // Create fields - const Real FillValue = -9.99e30; - int NDims = 2; + int NDims = 2; std::vector DimNames(NDims); std::string DimSuffix; if (MeshName == "Default") { @@ -43,9 +42,9 @@ void KineticAuxVars::registerFields( "specific_kinetic_energy_of_sea_water", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries - 2, // number of dimensions - DimNames // dim names + FillValueReal, // scalar for undefined entries + 2, // number of dimensions + DimNames // dim names ); // Velocity divergence on cells @@ -56,7 +55,7 @@ void KineticAuxVars::registerFields( "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries + FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp index f22e80968447..3e96b4940e78 100644 --- a/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp @@ -25,8 +25,7 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { // Create/define fields - const Real FillValue = -9.99e30; - int NDims = 2; + int NDims = 2; std::vector DimNames(NDims); std::string DimSuffix; if (MeshName == "Default") { @@ -47,7 +46,7 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -61,7 +60,7 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries + FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -74,7 +73,7 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp index 191c22e05043..37e81f770ef0 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp @@ -19,8 +19,7 @@ void SurfTracerRestAuxVars::registerFields( ) const { // Create fields - const Real FillValue = -9.99e30; - int NDims = 2; + int NDims = 2; std::vector DimNames(NDims); std::string DimSuffix; if (MeshName == "Default") { @@ -39,9 +38,9 @@ void SurfTracerRestAuxVars::registerFields( "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries - NDims, // number of dimensions - DimNames); // dim names + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames); // dim names // Add fields to FieldGroup FieldGroup::addFieldToGroup(TracersMonthlySurfClimoCell.label(), diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp index 9c2645441877..252bc64e3605 100644 --- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp @@ -21,8 +21,7 @@ void TracerAuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { // Create fields - const Real FillValue = -9.99e30; - int NDims = 3; + int NDims = 3; std::vector DimNames(NDims); std::string DimSuffix; if (MeshName == "Default") { @@ -43,7 +42,7 @@ void TracerAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp index 500473d18a7f..6ff83e743f78 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp @@ -34,8 +34,7 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { // Create/define fields - const Real FillValue = -9.99e30; - int NDims = 2; + int NDims = 2; std::vector DimNames(NDims); DimNames[1] = "NVertLayers"; std::string DimSuffix; @@ -54,7 +53,7 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undef entries + FillValueReal, // scalar for undef entries NDims, // number of dimensions DimNames // dimension names ); @@ -69,9 +68,9 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries - NDims, // number of dimensions - DimNames // dimension names + FillValueReal, // scalar used for undefined entries + NDims, // number of dimensions + DimNames // dimension names ); // Del2 of relative vorticity on vertices @@ -83,9 +82,9 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries - NDims, // number of dimensions - DimNames // dimension names + FillValueReal, // scalar used for undefined entries + NDims, // number of dimensions + DimNames // dimension names ); // Add fields to Aux Field group diff --git a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp index fabb2d53e0a1..d84666303c6e 100644 --- a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp @@ -35,8 +35,7 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, const std::string &MeshName) const { // Create fields with metadata - const Real FillValue = -9.99e30; - int NDims = 2; + int NDims = 2; std::vector DimNames(NDims); std::string DimSuffix; if (MeshName == "Default") { @@ -56,9 +55,9 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "ocean_relative_vorticity", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names ); // Normalized relative vorticity on vertices @@ -69,7 +68,7 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries + FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -83,7 +82,7 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries + FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -98,7 +97,7 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries + FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -112,9 +111,9 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used for undefined entries - NDims, // number of dimensions - DimNames // dimension names + FillValueReal, // scalar used for undefined entries + NDims, // number of dimensions + DimNames // dimension names ); // Add fields to Aux field group From f9f83d71c5727ef50e39f7337126009a65582e62 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 11:23:56 -0500 Subject: [PATCH 06/25] Determine fill values automatically by array type --- components/omega/src/base/FillValues.h | 25 +++++++++++++++++++------ components/omega/src/infra/Field.cpp | 9 +++------ components/omega/src/infra/Field.h | 24 +++++++++--------------- components/omega/src/infra/IOStream.cpp | 16 ++++++++-------- 4 files changed, 39 insertions(+), 35 deletions(-) diff --git a/components/omega/src/base/FillValues.h b/components/omega/src/base/FillValues.h index 1db933d96d3e..20ee9593e64b 100644 --- a/components/omega/src/base/FillValues.h +++ b/components/omega/src/base/FillValues.h @@ -16,14 +16,27 @@ namespace OMEGA { -constexpr I4 FillValueI4 = -2147483647; ///< NC_FILL_INT -constexpr I8 FillValueI8 = -9223372036854775806LL; ///< NC_FILL_INT64 -constexpr R4 FillValueR4 = 9.9692099683868690e+36f; ///< NC_FILL_FLOAT -constexpr R8 FillValueR8 = 9.9692099683868690e+36; ///< NC_FILL_DOUBLE +/// Type-indexed fill value. Instantiating FillValue for an unsupported +/// type T produces a link error (primary template is intentionally undefined). +template constexpr T FillValue; + +template <> inline constexpr I4 FillValue = -2147483647; ///< NC_FILL_INT +template <> +inline constexpr I8 FillValue = -9223372036854775806LL; ///< NC_FILL_INT64 +template <> +inline constexpr R4 FillValue = 9.9692099683868690e+36f; ///< NC_FILL_FLOAT +template <> +inline constexpr R8 FillValue = 9.9692099683868690e+36; ///< NC_FILL_DOUBLE + +/// Named aliases for readability in comparison and test code. +constexpr I4 FillValueI4 = FillValue; +constexpr I8 FillValueI8 = FillValue; +constexpr R4 FillValueR4 = FillValue; +constexpr R8 FillValueR8 = FillValue; #if defined(SINGLE_PRECISION) -constexpr Real FillValueReal = FillValueR4; +constexpr Real FillValueReal = FillValue; #else -constexpr Real FillValueReal = FillValueR8; +constexpr Real FillValueReal = FillValue; #endif } // end namespace OMEGA diff --git a/components/omega/src/infra/Field.cpp b/components/omega/src/infra/Field.cpp index 921ea7885aa0..d2e7095b5f79 100644 --- a/components/omega/src/infra/Field.cpp +++ b/components/omega/src/infra/Field.cpp @@ -54,8 +54,8 @@ void Field::init(const Clock *ModelClock // [in] default model clock std::string CalName = CalendarCFName[CalKind]; std::vector DimNamesTmp; // empty dim names vector std::shared_ptr TimeField = - create("time", "time", UnitString, "time", 0.0, 1.e20, FillValueR8, 0, - DimNamesTmp, true, true); + create("time", "time", UnitString, "time", 0.0, 1.e20, 0, DimNamesTmp, + true, true); TimeField->addMetadata("calendar", CalName); } @@ -85,8 +85,7 @@ Field::create(const std::string &FieldName, // [in] Name of variable/field const std::string &StdName, // [in] CF standard Name const std::any ValidMin, // [in] min valid field value const std::any ValidMax, // [in] max valid field value - const std::any FillValue, // [in] scalar for undefined entries - const int NumDims, // [in] number of dimensions + const int NumDims, // [in] number of dimensions const std::vector &Dimensions, // [in] dim names bool TimeDependent, // [in] flag for time dependent field bool RetainPrecision // [in] flag to retain full prec in IO @@ -118,8 +117,6 @@ Field::create(const std::string &FieldName, // [in] Name of variable/field ThisField->FieldMeta["valid_min"] = ValidMin; ThisField->FieldMeta["ValidMax"] = ValidMax; ThisField->FieldMeta["valid_max"] = ValidMax; - ThisField->FieldMeta["FillValue"] = FillValue; - ThisField->FieldMeta["_FillValue"] = FillValue; // Set the time-dependent flag ThisField->TimeDependent = TimeDependent; diff --git a/components/omega/src/infra/Field.h b/components/omega/src/infra/Field.h index 3ab67b54d579..69463a98d187 100644 --- a/components/omega/src/infra/Field.h +++ b/components/omega/src/infra/Field.h @@ -84,23 +84,18 @@ class Field { /// various types and cast to the appropriate type when needed. std::shared_ptr DataArray; - /// Fills every element of InDataArray with the field's declared fill value. + /// Fills every element of InDataArray with the standard fill value for the + /// array's element type, then records that value in the field metadata. /// Called automatically by attachData() so inactive entries are initialized /// to a well-defined sentinel before any compute routine runs. template void fillWithValue(const T &InDataArray) { - using ValType = typename T::value_type; - auto AnyFill = FieldMeta.at("FillValue"); - ValType FillVal; - if constexpr (std::is_same_v) - FillVal = std::any_cast(AnyFill); - else if constexpr (std::is_same_v) - FillVal = std::any_cast(AnyFill); - else if constexpr (std::is_same_v) - // Fill values are stored as R8 (Real) in legacy code; cast safely - FillVal = static_cast(std::any_cast(AnyFill)); - else - FillVal = std::any_cast(AnyFill); - Kokkos::deep_copy(InDataArray, FillVal); + using ValType = typename T::non_const_value_type; + constexpr ValType FillVal = FillValue; + // _FillValue is the NetCDF/CF standard attribute recognized by analysis + // tools (ncview, Xarray, NCO, ...); store only that, no non-standard + // alias + FieldMeta["_FillValue"] = FillVal; + Kokkos::deep_copy((InDataArray, FillVal); } public: @@ -133,7 +128,6 @@ class Field { const std::string &StdName, ///< [in] CF standard Name const std::any ValidMin, ///< [in] min valid field value const std::any ValidMax, ///< [in] max valid field value - const std::any FillValue, ///< [in] scalar for undefined entries const int NumDims, ///< [in] number of dimensions const std::vector &Dimensions, ///< [in] dim names const bool TimeDependent = true, ///< [in] opt flag for unlim time diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp index 95e082480ae9..9cab382c805e 100644 --- a/components/omega/src/infra/IOStream.cpp +++ b/components/omega/src/infra/IOStream.cpp @@ -1019,8 +1019,8 @@ void IOStream::writeFieldData( DataI4.resize(LocSize); DataPtr = DataI4.data(); // get fill value - Err += FieldPtr->getMetadata("FillValue", FillValI4); - CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}", + Err += FieldPtr->getMetadata("_FillValue", FillValI4); + CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}", FieldName); FillValPtr = &FillValI4; @@ -1161,8 +1161,8 @@ void IOStream::writeFieldData( DataI8.resize(LocSize); DataPtr = DataI8.data(); // Get fill value - Err += FieldPtr->getMetadata("FillValue", FillValI8); - CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}", + Err += FieldPtr->getMetadata("_FillValue", FillValI8); + CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}", FieldName); FillValPtr = &FillValI8; @@ -1302,8 +1302,8 @@ void IOStream::writeFieldData( DataR4.resize(LocSize); DataPtr = DataR4.data(); // Get fill value - Err += FieldPtr->getMetadata("FillValue", FillValR4); - CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}", + Err += FieldPtr->getMetadata("_FillValue", FillValR4); + CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}", FieldName); FillValPtr = &FillValR4; @@ -1441,8 +1441,8 @@ void IOStream::writeFieldData( case ArrayDataType::R8: // Get fill value - Err += FieldPtr->getMetadata("FillValue", FillValR8); - CHECK_ERROR_ABORT(Err, "Error retrieving FillValue for Field {}", + Err += FieldPtr->getMetadata("_FillValue", FillValR8); + CHECK_ERROR_ABORT(Err, "Error retrieving _FillValue for Field {}", FieldName); DataR8.resize(LocSize); if (ReducePrecision and !RetainPrecision) { From 515b930ec3837a8e3997ede82606595acdbbd2fe Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 11:25:20 -0500 Subject: [PATCH 07/25] Drop manually assigned fill values --- components/omega/src/ocn/Eos.cpp | 15 +- components/omega/src/ocn/HorzMesh.cpp | 189 ++++++++---------- components/omega/src/ocn/OceanState.cpp | 6 +- components/omega/src/ocn/Tendencies.cpp | 6 +- components/omega/src/ocn/TracerDefs.inc | 12 +- components/omega/src/ocn/Tracers.cpp | 9 +- components/omega/src/ocn/Tracers.h | 1 - components/omega/src/ocn/VertAdv.cpp | 8 +- components/omega/src/ocn/VertCoord.cpp | 20 +- components/omega/src/ocn/VertMix.cpp | 10 +- .../src/ocn/auxiliaryVars/KineticAuxVars.cpp | 6 +- .../auxiliaryVars/PseudoThicknessAuxVars.cpp | 3 - .../auxiliaryVars/SurfTracerRestAuxVars.cpp | 5 +- .../src/ocn/auxiliaryVars/TracerAuxVars.cpp | 1 - .../ocn/auxiliaryVars/VelocityDel2AuxVars.cpp | 11 +- .../ocn/auxiliaryVars/VorticityAuxVars.cpp | 13 +- 16 files changed, 122 insertions(+), 193 deletions(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index d2ecf014abb0..16ea4da9f5a2 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -370,9 +370,8 @@ void Eos::defineFields() { "sea_water_specific_volume", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValueReal, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + NDims, // Number of dimensions + DimNames // Dimension names ); /// Create and register the displaced specific volume field auto SpecVolDisplacedField = @@ -383,9 +382,8 @@ void Eos::defineFields() { "sea_water_specific_volume_displaced", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValueReal, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + NDims, // Number of dimensions + DimNames // Dimension names ); // Brunt-Vaisala frequency is located at interfaces @@ -399,9 +397,8 @@ void Eos::defineFields() { "sea_water_brunt_vaisala_frequency_squared", // CF-ish Name std::numeric_limits::min(), // Min valid value std::numeric_limits::max(), // Max valid value - FillValueReal, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + NDims, // Number of dimensions + DimNames // Dimension names ); // Create a field group for the eos-specific state fields diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 6fde75403f2c..bb8b909aad42 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -559,10 +559,9 @@ void HorzMesh::defineMeshFields() { "x", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, XCell); @@ -576,10 +575,9 @@ void HorzMesh::defineMeshFields() { "y", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, YCell); @@ -593,10 +591,9 @@ void HorzMesh::defineMeshFields() { "z", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, ZCell); @@ -610,10 +607,9 @@ void HorzMesh::defineMeshFields() { "latitude", // CF standard Name -3.1415927, // min valid value 3.1415927, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LatCell); @@ -624,13 +620,12 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Longitude coordinates of cell centers", // long Name "radians", // units - "longitude", // CF standard Name - 0.0, // min valid value - 6.28319, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "longitude", // CF standard Name + 0.0, // min valid value + 6.28319, // max valid value + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LonCell); @@ -646,10 +641,9 @@ void HorzMesh::defineMeshFields() { "x", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, XEdge); @@ -663,10 +657,9 @@ void HorzMesh::defineMeshFields() { "y", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, YEdge); @@ -680,10 +673,9 @@ void HorzMesh::defineMeshFields() { "z", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, ZEdge); @@ -697,10 +689,9 @@ void HorzMesh::defineMeshFields() { "latitude", // CF standard Name -3.1415927, // min valid value 3.1415927, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LatEdge); @@ -714,10 +705,9 @@ void HorzMesh::defineMeshFields() { "longitude", // CF standard Name 0.0, // min valid value 6.28319, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LonEdge); @@ -733,10 +723,9 @@ void HorzMesh::defineMeshFields() { "x", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, XVertex); @@ -750,10 +739,9 @@ void HorzMesh::defineMeshFields() { "y", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, YVertex); @@ -767,10 +755,9 @@ void HorzMesh::defineMeshFields() { "z", // CF standard Name 0.0, // min valid value 7.0E+6, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, ZVertex); @@ -781,13 +768,12 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Latitude coordinates of cell vertices", // long Name "radians", // units - "latitude", // CF standard Name - -3.1415927, // min valid value - 3.1415927, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "latitude", // CF standard Name + -3.1415927, // min valid value + 3.1415927, // max valid value + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LatVertex); @@ -798,34 +784,31 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Longitude coordinates of cell vertices", // long Name "radians", // units - "longitude", // CF standard Name - 0.0, // min valid value - 6.28319, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "longitude", // CF standard Name + 0.0, // min valid value + 6.28319, // max valid value + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, LonVertex); // Other mesh properties // Mesh areas, lengths, and angles - DimNames[0] = "NCells"; - FieldName = "AreaCell"; - AreaCell = Array1DReal("AreaCell", NCellsSize); // allocate space - auto AreaCellField = - Field::create(FieldName, // field name - "Area of each cell (m^2)", // long name - "m2", // units - "cell_area", // CF standard name - 0.0, // min valid value - 9.99E+30, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent - ); + DimNames[0] = "NCells"; + FieldName = "AreaCell"; + AreaCell = Array1DReal("AreaCell", NCellsSize); // allocate space + auto AreaCellField = Field::create(FieldName, // field name + "Area of each cell (m^2)", // long name + "m2", // units + "cell_area", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent + ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, AreaCell); @@ -839,7 +822,6 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -854,13 +836,12 @@ void HorzMesh::defineMeshFields() { Field::create(FieldName, // field name "Area of each triangle in the dual grid (m^2)", // lng name "m2", // units - "cell_area", // CF standard name - 0.0, // min valid value - 9.99E+30, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + "cell_area", // CF standard name + 0.0, // min valid value + 9.99E+30, // max valid value + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, AreaTriangle); @@ -876,7 +857,6 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -894,7 +874,6 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -907,15 +886,14 @@ void HorzMesh::defineMeshFields() { auto AngleEdgeField = Field::create(FieldName, // field name "Angle the edge normal makes with local eastward direction" - " (radians)", // long name - "radians", // units - "", // CF standard name - -3.1415927, // min valid value - 3.1415927, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // num of dimensions - DimNames, // dimension names - false // not time dependent + " (radians)", // long name + "radians", // units + "", // CF standard name + -3.1415927, // min valid value + 3.1415927, // max valid value + NDims, // num of dimensions + DimNames, // dimension names + false // not time dependent ); MeshGroupIn->addField(FieldName); Field::attachFieldData(FieldName, AngleEdge); @@ -935,7 +913,6 @@ void HorzMesh::defineMeshFields() { "", // CF standard name 0.0, // min valid value 9.99E+30, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -956,7 +933,6 @@ void HorzMesh::defineMeshFields() { "", // CF standard name -1.0, // min valid value 1.0, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -977,7 +953,6 @@ void HorzMesh::defineMeshFields() { "coriolis_parameter", // CF standard name -1.0E-3, // min valid value 1.0E-3, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -995,7 +970,6 @@ void HorzMesh::defineMeshFields() { "coriolis_parameter", // CF standard name -1.0E-3, // min valid value 1.0E-3, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent @@ -1013,7 +987,6 @@ void HorzMesh::defineMeshFields() { "coriolis_parameter", // CF standard name -1.0E-3, // min valid value 1.0E-3, // max valid value - FillValueReal, // scalar for undefined entries NDims, // num of dimensions DimNames, // dimension names false // not time dependent diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index 6ee850207f95..35b607fa8149 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/components/omega/src/ocn/OceanState.cpp @@ -202,9 +202,8 @@ void OceanState::defineFields() { "sea_water_velocity", // CF standard Name -9.99E+10, // min valid value 9.99E+10, // max valid value - FillValueReal, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); DimNames[0] = "NCells"; @@ -215,7 +214,6 @@ void OceanState::defineFields() { "cell_thickness", // CF standard Name 0.0, // min valid value 9.99E+30, // max valid value - FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 77851168871b..6d477f3f6274 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -338,7 +338,7 @@ void Tendencies::defineFields() { auto PseudoThicknessTendField = Field::create(PseudoThicknessTendFieldName, "Pseudo-thickness tendency", "m/s", "cell_thickness_tendency", -9.99E+10, 9.99E+10, - FillValueReal, NDims, DimNamesThickness); + NDims, DimNamesThickness); NDims = 3; std::vector DimNamesTracer(NDims); DimNamesTracer[0] = "NTracers"; @@ -346,7 +346,7 @@ void Tendencies::defineFields() { DimNamesTracer[2] = "NVertLayers"; auto TracerTendField = Field::create( TracerTendFieldName, "Tracer tendency", "kg/m^3/s", "tracer_tendency", - -9.99E+10, 9.99E+10, FillValueReal, NDims, DimNamesTracer); + -9.99E+10, 9.99E+10, NDims, DimNamesTracer); NDims = 2; std::vector DimNamesVelocity(NDims); DimNamesVelocity[0] = "NEdges"; @@ -354,7 +354,7 @@ void Tendencies::defineFields() { auto NormalVelocityTendField = Field::create(NormalVelocityTendFieldName, "Normal velocity tendency", "m/s^2", "sea_water_velocity_tendency", -9.99E+10, - 9.99E+10, FillValueReal, NDims, DimNamesVelocity); + 9.99E+10, NDims, DimNamesVelocity); std::string TendGroupName = "Tendencies"; if (Name != "Default") { diff --git a/components/omega/src/ocn/TracerDefs.inc b/components/omega/src/ocn/TracerDefs.inc index 03134e75da5d..adafaa945f11 100644 --- a/components/omega/src/ocn/TracerDefs.inc +++ b/components/omega/src/ocn/TracerDefs.inc @@ -27,16 +27,12 @@ static void defineAllTracers() { "sea_water_conservative_temperature", ///< [in] CF standard Name -273.15, ///< [in] min valid field value 100.0, ///< [in] max valid field value - FillValueReal, ///< [in] value of undef entries IndxTemp); ///< [out] (optional) static index define("Salinity", "Absolute Salinity", "g kg-1", - "sea_water_absolute_salinity", 0.0, 50.0, FillValueReal, IndxSalt); + "sea_water_absolute_salinity", 0.0, 50.0, IndxSalt); - define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, FillValueReal, - IndxDebug1); - define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, FillValueReal, - IndxDebug2); - define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, FillValueReal, - IndxDebug3); + define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, IndxDebug1); + define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, IndxDebug2); + define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, IndxDebug3); } diff --git a/components/omega/src/ocn/Tracers.cpp b/components/omega/src/ocn/Tracers.cpp index efa7b00b635f..bb7735b3b7cc 100644 --- a/components/omega/src/ocn/Tracers.cpp +++ b/components/omega/src/ocn/Tracers.cpp @@ -177,8 +177,7 @@ void Tracers::init() { //--------------------------------------------------------------------------- I4 Tracers::define(const std::string &Name, const std::string &Description, const std::string &Units, const std::string &StdName, - const Real ValidMin, const Real ValidMax, - const Real FillValue, I4 &Index) { + const Real ValidMin, const Real ValidMax, I4 &Index) { // Do nothing if this tracer is not selected if (TracerIndexes.find(Name) == TracerIndexes.end()) { @@ -201,9 +200,9 @@ I4 Tracers::define(const std::string &Name, const std::string &Description, // create a tracer field std::string TracerFieldName = Name; - auto TracerField = Field::create(TracerFieldName, Description, Units, - StdName, ValidMin, ValidMax, FillValue, - TracerDimNames.size(), TracerDimNames); + auto TracerField = + Field::create(TracerFieldName, Description, Units, StdName, ValidMin, + ValidMax, TracerDimNames.size(), TracerDimNames); if (!TracerField) { LOG_ERROR("Tracers: Tracer field '{}' is not created", TracerFieldName); return -2; diff --git a/components/omega/src/ocn/Tracers.h b/components/omega/src/ocn/Tracers.h index 6480880a6235..a1ee899056b8 100644 --- a/components/omega/src/ocn/Tracers.h +++ b/components/omega/src/ocn/Tracers.h @@ -81,7 +81,6 @@ class Tracers { const std::string &StdName, ///< [in] CF standard Name const Real ValidMin, ///< [in] min valid field value const Real ValidMax, ///< [in] max valid field value - const Real FillValue, ///< [in] value for undef entries I4 &Index = IndxInvalid ///< [out] (optional) index value ); diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index dee2b4e2d90b..fc3fd4d99a11 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -142,7 +142,6 @@ void VertAdv::defineFields() { "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -155,9 +154,8 @@ void VertAdv::defineFields() { "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); NDims = 3; @@ -173,7 +171,6 @@ void VertAdv::defineFields() { "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -186,7 +183,6 @@ void VertAdv::defineFields() { "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 536a9bfa9a3d..451e5218dd9a 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -258,9 +258,8 @@ void VertCoord::defineFields() { "", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValueI4, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); auto MaxLayerCellField = Field::create( @@ -270,9 +269,8 @@ void VertCoord::defineFields() { "", // CF standard Name -1, // min valid value std::numeric_limits::max(), // max valid value - FillValueI4, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); auto BottomGeomDepthField = Field::create( @@ -283,7 +281,6 @@ void VertCoord::defineFields() { "sea_floor_depth_below_geoid", // CF standard Name 0.0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -295,7 +292,6 @@ void VertCoord::defineFields() { "sea_surface_height", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -312,7 +308,6 @@ void VertCoord::defineFields() { "", // CF standard Name 0.0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -329,7 +324,6 @@ void VertCoord::defineFields() { "", // CF standard Name 0.0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -346,7 +340,6 @@ void VertCoord::defineFields() { "sea_water_pressure", // CF standard Name 0.0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -358,7 +351,6 @@ void VertCoord::defineFields() { "height", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -372,7 +364,6 @@ void VertCoord::defineFields() { "sea_water_pressure", // CF standard Name 0.0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -384,7 +375,6 @@ void VertCoord::defineFields() { "height", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -396,7 +386,6 @@ void VertCoord::defineFields() { "geopotential", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -409,7 +398,6 @@ void VertCoord::defineFields() { "", // CF standard Name 0.0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 01c923e81b23..b6ea017406fc 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -346,9 +346,8 @@ void VertMix::defineFields() { "vertical_diffusivity", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValueReal, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + NDims, // Number of dimensions + DimNames // Dimension names ); /// Create and register the VertVisc field auto VertViscField = @@ -359,9 +358,8 @@ void VertMix::defineFields() { "vertical_viscosity", // CF-ish Name 0.0, // Min valid value std::numeric_limits::max(), // Max valid value - FillValueReal, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + NDims, // Number of dimensions + DimNames // Dimension names ); /// Create and register the GradRichNum field auto GradRichNumField = diff --git a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp index 8c76aabaa922..5347699acdcc 100644 --- a/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/KineticAuxVars.cpp @@ -42,9 +42,8 @@ void KineticAuxVars::registerFields( "specific_kinetic_energy_of_sea_water", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries - 2, // number of dimensions - DimNames // dim names + 2, // number of dimensions + DimNames // dim names ); // Velocity divergence on cells @@ -55,7 +54,6 @@ void KineticAuxVars::registerFields( "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp index 3e96b4940e78..2a608c176dc1 100644 --- a/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/PseudoThicknessAuxVars.cpp @@ -46,7 +46,6 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -60,7 +59,6 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -73,7 +71,6 @@ void PseudoThicknessAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name 0, // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp index 37e81f770ef0..cc075b8ab634 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp @@ -38,9 +38,8 @@ void SurfTracerRestAuxVars::registerFields( "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries - NDims, // number of dimensions - DimNames); // dim names + NDims, // number of dimensions + DimNames); // dim names // Add fields to FieldGroup FieldGroup::addFieldToGroup(TracersMonthlySurfClimoCell.label(), diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp index 252bc64e3605..6cbaffdf4dce 100644 --- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp @@ -42,7 +42,6 @@ void TracerAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries NDims, // number of dimensions DimNames // dimension names ); diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp index 6ff83e743f78..e67da008c01b 100644 --- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp @@ -53,7 +53,6 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undef entries NDims, // number of dimensions DimNames // dimension names ); @@ -68,9 +67,8 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); // Del2 of relative vorticity on vertices @@ -82,9 +80,8 @@ void VelocityDel2AuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); // Add fields to Aux Field group diff --git a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp index d84666303c6e..062aeac80cf6 100644 --- a/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/VorticityAuxVars.cpp @@ -55,9 +55,8 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "ocean_relative_vorticity", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); // Normalized relative vorticity on vertices @@ -68,7 +67,6 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -82,7 +80,6 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -97,7 +94,6 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries NDims, // number of dimensions DimNames // dimension names ); @@ -111,9 +107,8 @@ void VorticityAuxVars::registerFields(const std::string &AuxGroupName, "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar used for undefined entries - NDims, // number of dimensions - DimNames // dimension names + NDims, // number of dimensions + DimNames // dimension names ); // Add fields to Aux field group From f4013d09bd757c4ffd1ac5a6524182b2f6e60a90 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 11:33:18 -0500 Subject: [PATCH 08/25] Make detection of fill values in VertCoord reads less brittle --- components/omega/src/ocn/VertCoord.cpp | 95 ++++++++++++++------------ 1 file changed, 53 insertions(+), 42 deletions(-) diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 451e5218dd9a..d4e160437ad3 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -535,82 +535,93 @@ void VertCoord::setStreamArrays(const bool ReadStream, Halo *MeshHalo) { Err = IOStream::read(StreamName); if (Err.isFail()) { LOG_INFO("VertCoord: Error while reading {} stream", StreamName); - I4 Sum1 = 0; + I4 NFill1 = 0; parallelReduce( {MinLayerCell.extent_int(0)}, - KOKKOS_LAMBDA(int I, int &Accum) { - Accum += LocMinLayerCell(I); + KOKKOS_LAMBDA(int I, I4 &Count) { + if (LocMinLayerCell(I) == FillValueI4) + ++Count; }, - Sum1); - if (Sum1 < 0) { + NFill1); + if (NFill1 > 0) { LOG_INFO("VertCoord: Error reading MinLayerCell from {}, " "using MinLayerCell = 0", StreamName); deepCopy(MinLayerCell, 1); } - I4 Sum2 = 0; + I4 NFill2 = 0; parallelReduce( {MaxLayerCell.extent_int(0)}, - KOKKOS_LAMBDA(int I, int &Accum) { - Accum += LocMaxLayerCell(I); + KOKKOS_LAMBDA(int I, I4 &Count) { + if (LocMaxLayerCell(I) == FillValueI4) + ++Count; }, - Sum2); - if (Sum2 < 0) { + NFill2); + if (NFill2 > 0) { LOG_INFO("VertCoord: Error reading MaxLayerCell from {}, " "using MaxLayerCell = NVertLayers - 1", StreamName); deepCopy(MaxLayerCell, NVertLayers); } - Real Sum3 = 0.; + I4 NFill3 = 0; parallelReduce( {BottomGeomDepth.extent_int(0)}, - KOKKOS_LAMBDA(int I, Real &Accum) { - Accum += LocBottomGeomDepth(I); + KOKKOS_LAMBDA(int I, I4 &Count) { + if (LocBottomGeomDepth(I) == FillValueReal) + ++Count; }, - Sum3); - if (Sum3 / BottomGeomDepth.extent_int(0) >= - 0.5_Real * FillValueReal) { + NFill3); + if (NFill3 > 0) { ABORT_ERROR("VertCoord: Error reading BottomGeomDepth from {}", StreamName); } - Real Sum4 = 0.; + I4 NFill4 = 0; parallelReduce( {RefPseudoThickness.extent_int(0), RefPseudoThickness.extent_int(1)}, - KOKKOS_LAMBDA(int I, int J, Real &Accum) { - Accum += LocPseudoThick(I, J); + KOKKOS_LAMBDA(int I, int J, I4 &Count) { + if (LocPseudoThick(I, J) == FillValueReal) + ++Count; }, - Sum4); - I4 N4 = RefPseudoThickness.extent_int(0) * - RefPseudoThickness.extent_int(1); - if (Sum4 / N4 >= 0.5_Real * FillValueReal) { + NFill4); + if (NFill4 > 0) { ABORT_ERROR("VertCoord: Error reading RefPseudoThickness " "from {}", StreamName); } - Real Sum5 = 0.; - I4 NumNonZero = 0; + I4 NFill5 = 0; parallelReduce( {VertCoordMovementWeights.extent_int(0)}, - KOKKOS_LAMBDA(int I, Real &Accum, I4 &NonZeroCount) { - Real W = LocVCoordMvmtWgts(I); - Accum += W; - if (W != 0) { - NonZeroCount += 1; - } + KOKKOS_LAMBDA(int I, I4 &Count) { + if (LocVCoordMvmtWgts(I) == FillValueReal) + ++Count; }, - Sum5, NumNonZero); - I4 N5 = VertCoordMovementWeights.extent_int(0); - if (Sum5 / N5 >= 0.5_Real * FillValueReal) { + NFill5); + if (NFill5 > 0) { // Stream did not populate weights — use default deepCopy(VertCoordMovementWeights, 1._Real); - } else if (Sum5 < 0.) { - ABORT_ERROR("VertCoord: Error reading VertCoordMovementWeights " - "from {}", - StreamName); - } else if (NumNonZero == 0) { - // TODO: ABORT_ERROR when all weights equal 0 - deepCopy(VertCoordMovementWeights, 1._Real); + } else { + Real Sum5 = 0.; + I4 NumNonZero = 0; + parallelReduce( + {VertCoordMovementWeights.extent_int(0)}, + KOKKOS_LAMBDA(int I, Real &Accum, I4 &NonZeroCount) { + Real W = LocVCoordMvmtWgts(I); + Accum += W; + if (W != 0) { + NonZeroCount += 1; + } + }, + Sum5, NumNonZero); + if (Sum5 < 0.) { + ABORT_ERROR( + "VertCoord: Error reading VertCoordMovementWeights " + "from {}", + StreamName); + } else if (NumNonZero == 0) { + // TODO: ABORT_ERROR when all weights equal 0 + deepCopy(VertCoordMovementWeights, 1._Real); + } } } } else { @@ -1014,7 +1025,7 @@ void VertCoord::computeGeomZHeight( Team, KRange, INNER_LAMBDA(int K, Real &Accum, bool IsFinal) { const I4 KLyr = KMax - K; Real DZ = RhoSw * SpecVol(ICell, KLyr) * - PseudoThickness(ICell, KLyr); + PseudoThickness(ICell, KLyr); Accum += DZ; if (IsFinal) { LocZInterf(ICell, KLyr) = -LocBotGeomDepth(ICell) + Accum; From 775af6a59048e14bd0a53cb5843553c7035ad7b3 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 12:42:19 -0500 Subject: [PATCH 09/25] Add FillOnAttach argument to attachData() In some cases we definitely don't want to overwrite fields we attach with fill values. Fix 2 such cases in Tracers and OceanState --- components/omega/src/infra/Field.h | 15 ++++++++++----- components/omega/src/ocn/OceanState.cpp | 7 ++++--- components/omega/src/ocn/Tracers.cpp | 3 ++- 3 files changed, 16 insertions(+), 9 deletions(-) diff --git a/components/omega/src/infra/Field.h b/components/omega/src/infra/Field.h index 69463a98d187..a94d0eb01d47 100644 --- a/components/omega/src/infra/Field.h +++ b/components/omega/src/infra/Field.h @@ -318,7 +318,8 @@ class Field { /// templated based on the array data type so a template argument with /// the proper array data type (eg ) must be supplied. template - void attachData(const T &InDataArray ///< [in] Array with data to attach + void attachData(const T &InDataArray, ///< [in] Array with data to attach + bool FillOnAttach = true ///< [in] fill array with fill value ) { OMEGA_ASSERT(isKokkosArray, "Field::attachData requires Kokkos array as input"); @@ -331,8 +332,11 @@ class Field { MemLoc = findArrayMemLoc(); // Initialize every element to the declared fill value so inactive - // entries (e.g. layers below MaxLayerCell) are well-defined in output - fillWithValue(InDataArray); + // entries (e.g. layers below MaxLayerCell) are well-defined in output. + // Pass FillOnAttach=false when re-attaching an existing data-filled array + // (e.g. time-level pointer updates) to preserve the computed values. + if (FillOnAttach) + fillWithValue(InDataArray); }; //--------------------------------------------------------------------------- @@ -344,7 +348,8 @@ class Field { template static void attachFieldData(const std::string &FieldName, ///< [in] Name of Field - const T &InDataArray ///< [in] Array with data to attach + const T &InDataArray, ///< [in] Array with data to attach + bool FillOnAttach = true ///< [in] fill array with fill value ) { OMEGA_ASSERT(isKokkosArray, "Field::attachFieldData requires Kokkos array as input"); @@ -354,7 +359,7 @@ class Field { // Retrieve the field auto ThisField = AllFields[FieldName]; // Attach the data array - ThisField->attachData(InDataArray); + ThisField->attachData(InDataArray, FillOnAttach); } else { // field has not yet been defined ABORT_ERROR("Field: error attaching data to {}. Field not defined", diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index 35b607fa8149..b74e4df1d0c5 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/components/omega/src/ocn/OceanState.cpp @@ -318,11 +318,12 @@ void OceanState::updateTimeLevels() { // Update current time index for pseudo-thickness and normal velocity CurTimeIndex = (CurTimeIndex + 1) % NTimeLevels; - // Update IOField data associations + // Update IOField data associations. Pass false to avoid overwriting the + // just-computed state with fill values; this is a pointer update only. Field::attachFieldData(NormalVelocityFldName, - NormalVelocity[CurTimeIndex]); + NormalVelocity[CurTimeIndex], false); Field::attachFieldData(PseudoThicknessFldName, - PseudoThickness[CurTimeIndex]); + PseudoThickness[CurTimeIndex], false); } // end updateTimeLevels diff --git a/components/omega/src/ocn/Tracers.cpp b/components/omega/src/ocn/Tracers.cpp index bb7735b3b7cc..c84fa5b33454 100644 --- a/components/omega/src/ocn/Tracers.cpp +++ b/components/omega/src/ocn/Tracers.cpp @@ -446,7 +446,8 @@ void Tracers::updateTimeLevels() { Array2DReal TracerSubview = Kokkos::subview( TracerArrays[CurTimeIndex], TracerIndex, Kokkos::ALL, Kokkos::ALL); - TracerField->attachData(TracerSubview); + // Pass false to avoid overwriting just-computed tracer values with fill. + TracerField->attachData(TracerSubview, false); } return; From 3d2c737673405b61abe37ac52bd0a7c65b54e17d Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 12:44:00 -0500 Subject: [PATCH 10/25] Reorder 2 cases so we attach first, then fill in the array --- components/omega/src/infra/IOStream.cpp | 6 ++- components/omega/test/infra/FieldTest.cpp | 54 +++++++++++------------ 2 files changed, 31 insertions(+), 29 deletions(-) diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp index 9cab382c805e..873322fca62f 100644 --- a/components/omega/src/infra/IOStream.cpp +++ b/components/omega/src/infra/IOStream.cpp @@ -2482,13 +2482,15 @@ void IOStream::writeStream( if (Multiframe and !FileAlarm.isRinging()) FileTime = *(FileAlarm.getRingTimePrev()); - // Update the time field with elapsed time since simulation start + // Update the time field with elapsed time since simulation start. + // Attach before setting the value so the fill-on-attach does not + // overwrite the elapsed time (Kokkos host views share backing memory). TimeInterval ElapsedTime = SimTime - StartTime; R8 ElapsedTimeR8; ElapsedTime.get(ElapsedTimeR8, TimeUnits::Seconds); HostArray1DR8 OutTime("OutTime", 1); - OutTime(0) = ElapsedTimeR8; Field::attachFieldData("time", OutTime); + OutTime(0) = ElapsedTimeR8; // Reset alarms and flags if (OnStartup) diff --git a/components/omega/test/infra/FieldTest.cpp b/components/omega/test/infra/FieldTest.cpp index f3032cea8a77..530a5436f900 100644 --- a/components/omega/test/infra/FieldTest.cpp +++ b/components/omega/test/infra/FieldTest.cpp @@ -322,6 +322,33 @@ void initFieldTest() { Array5DR4 Data5DR4("Test5DR4", NCellsSize, NVertLayers, NTracers, NTime, NStuff); + // Attach arrays before populating so the fill-on-attach does not overwrite + // the reference values written below. Kokkos views share backing memory, so + // assignments through the local handle are visible via the stored pointer. + Test1DI4H->attachData(Data1DI4H); + Test1DI8H->attachData(Data1DI8H); + Test1DR4H->attachData(Data1DR4H); + Test1DR8H->attachData(Data1DR8H); + + Field::attachFieldData("Test2DI4H", Data2DI4H); + Field::attachFieldData("Test2DI8H", Data2DI8H); + Field::attachFieldData("Test2DR4H", Data2DR4H); + Field::attachFieldData("Test2DR8H", Data2DR8H); + + Field::attachFieldData("Test1DI4", Data1DI4); + Field::attachFieldData("Test1DI8", Data1DI8); + Field::attachFieldData("Test1DR4", Data1DR4); + Field::attachFieldData("Test1DR8", Data1DR8); + + Test2DI4->attachData(Data2DI4); + Test2DI8->attachData(Data2DI8); + Test2DR4->attachData(Data2DR4); + Test2DR8->attachData(Data2DR8); + + Test3DI4->attachData(Data3DI4); + Test4DI8->attachData(Data4DI8); + Test5DR4->attachData(Data5DR4); + // Host arrays vertical vector for (int K = 0; K < NVertLayers; ++K) { Data1DR8H(K) = RefR8 + K; @@ -388,33 +415,6 @@ void initFieldTest() { } }); - // Attach data arrays - // Use member function for some and name interface for others - - Test1DI4H->attachData(Data1DI4H); - Test1DI8H->attachData(Data1DI8H); - Test1DR4H->attachData(Data1DR4H); - Test1DR8H->attachData(Data1DR8H); - - Field::attachFieldData("Test2DI4H", Data2DI4H); - Field::attachFieldData("Test2DI8H", Data2DI8H); - Field::attachFieldData("Test2DR4H", Data2DR4H); - Field::attachFieldData("Test2DR8H", Data2DR8H); - - Field::attachFieldData("Test1DI4", Data1DI4); - Field::attachFieldData("Test1DI8", Data1DI8); - Field::attachFieldData("Test1DR4", Data1DR4); - Field::attachFieldData("Test1DR8", Data1DR8); - - Test2DI4->attachData(Data2DI4); - Test2DI8->attachData(Data2DI8); - Test2DR4->attachData(Data2DR4); - Test2DR8->attachData(Data2DR8); - - Test3DI4->attachData(Data3DI4); - Test4DI8->attachData(Data4DI8); - Test5DR4->attachData(Data5DR4); - // End of init } // End initialization Fields From 4505f492222d866dfaa5bbaff3e2d537fd925f75 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 17:14:53 -0500 Subject: [PATCH 11/25] Fill fluxes on edges first with zeros over full physical range We don't want fill values at every boundary edge becausee these are also valid. Instead, we want zeros except for edges between two inactive cells (e.g. both adjacent cells are below bathymetry) --- components/omega/src/ocn/AuxiliaryState.cpp | 6 +++ components/omega/src/ocn/Tendencies.cpp | 15 +++----- components/omega/src/ocn/VertCoord.cpp | 42 ++++++++++++++++++++- components/omega/src/ocn/VertCoord.h | 13 +++++++ 4 files changed, 66 insertions(+), 10 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 7ac27c57adf8..7c0bc3d71ca5 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -167,6 +167,12 @@ void AuxiliaryState::computeMomAux(const OceanState *State, const auto &VelocityDivCell = KineticAux.VelocityDivCell; const auto &RelVortVertex = VorticityAux.RelVortVertex; + // Zero Del2Edge at boundary layers before accumulating over neighbouring + // edges in computeVarsOnCell/Vertex. Del4 hyperdiffusion reads Del2Edge at + // all edges sharing a cell or vertex; fill values in boundary layers would + // corrupt Del2DivCell and Del2RelVortVertex. + VCoord->zeroEdgeField(VelocityDel2Aux.Del2Edge, Mesh->NEdgesAll); + Pacer::start("AuxState:edgeAuxState2", 2); parallelForOuter( "edgeAuxState2", {Mesh->NEdgesAll}, diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 6d477f3f6274..eca1252826b7 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -517,15 +517,12 @@ void Tendencies::computeVelocityTendenciesOnly( Pacer::start("Tend:computeVelocityTendenciesOnly", 1); - parallelForOuter( - {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { - const int KMin = MinLayerEdgeBot(IEdge); - const int KMax = MaxLayerEdgeTop(IEdge); - - parallelForInner( - Team, Range{KMin, KMax}, - INNER_LAMBDA(int K) { LocNormalVelocityTend(IEdge, K) = 0; }); - }); + // Zero NVTend over all layers where at least one neighboring cell is active: + // [MinLayerEdgeTop, MaxLayerEdgeBot]. This includes boundary edges adjacent + // to land or bathymetry steps, which have no computed tendency but should + // show 0, not FillValueReal. Fully inactive layers (outside this range) + // keep their FillValueReal from attachData(). + VCoord->zeroEdgeField(NormalVelocityTend, Mesh->NEdgesAll); // Compute potential vorticity horizontal advection const Array2DReal &FluxPseudoThickEdge = diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index d4e160437ad3..f22077026e3b 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -1025,7 +1025,7 @@ void VertCoord::computeGeomZHeight( Team, KRange, INNER_LAMBDA(int K, Real &Accum, bool IsFinal) { const I4 KLyr = KMax - K; Real DZ = RhoSw * SpecVol(ICell, KLyr) * - PseudoThickness(ICell, KLyr); + PseudoThickness(ICell, KLyr); Accum += DZ; if (IsFinal) { LocZInterf(ICell, KLyr) = -LocBotGeomDepth(ICell) + Accum; @@ -1180,6 +1180,46 @@ VertCoord *VertCoord::get(const std::string Name ///< [in] Name of VertCoord } // end get VertCoord +//------------------------------------------------------------------------------ +// Zero all layers in [MinLayerEdgeTop, MaxLayerEdgeBot] of an edge field. +void VertCoord::zeroEdgeField(Array2DReal Arr, I4 NEdgesAll) const { + OMEGA_SCOPE(LocMinLayerEdgeTop, MinLayerEdgeTop); + OMEGA_SCOPE(LocMaxLayerEdgeBot, MaxLayerEdgeBot); + parallelForOuter( + {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KTop = LocMinLayerEdgeTop(IEdge); + const int KBot = Kokkos::max(0, LocMaxLayerEdgeBot(IEdge)); + parallelForInner( + Team, Range{KTop, KBot}, + INNER_LAMBDA(int K) { Arr(IEdge, K) = 0._Real; }); + }); +} + +//------------------------------------------------------------------------------ +// Enforce 3-zone masking on an edge field after IC or restart read. +void VertCoord::applyEdgeLayerMask(Array2DReal Arr, I4 NEdgesAll) const { + OMEGA_SCOPE(LocMinLayerEdgeTop, MinLayerEdgeTop); + OMEGA_SCOPE(LocMaxLayerEdgeBot, MaxLayerEdgeBot); + OMEGA_SCOPE(LocMinLayerEdgeBot, MinLayerEdgeBot); + OMEGA_SCOPE(LocMaxLayerEdgeTop, MaxLayerEdgeTop); + I4 LocNVertLayers = NVertLayers; + parallelForOuter( + {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KTop = LocMinLayerEdgeTop(IEdge); + const int KBot = Kokkos::max(0, LocMaxLayerEdgeBot(IEdge)); + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + parallelForInner( + Team, Range{0, LocNVertLayers - 1}, INNER_LAMBDA(int K) { + if (K < KTop || K > KBot) + Arr(IEdge, K) = FillValueReal; + else if (K < KMin || K > KMax) + Arr(IEdge, K) = 0._Real; + // else: active layer - keep IC/restart value + }); + }); +} + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h index 1e478a04c373..848833e9c9fc 100644 --- a/components/omega/src/ocn/VertCoord.h +++ b/components/omega/src/ocn/VertCoord.h @@ -242,6 +242,19 @@ class VertCoord { /// Initialize computational masks void setMasks(); + /// Zero all layers in [MinLayerEdgeTop, MaxLayerEdgeBot] of an edge field. + /// Call before computing a flux-type edge field so boundary edges show 0, + /// not FillValueReal. Layers outside the valid range retain their fill + /// value. + void zeroEdgeField(Array2DReal Arr, I4 NEdgesAll) const; + + /// Enforce 3-zone masking on an edge field after IC or restart read: + /// K < MinLayerEdgeTop or K > MaxLayerEdgeBot -> FillValueReal + /// [MinLayerEdgeTop, MaxLayerEdgeBot] but outside + /// [MinLayerEdgeBot, MaxLayerEdgeTop] -> 0 + /// [MinLayerEdgeBot, MaxLayerEdgeTop] -> unchanged (IC/restart values) + void applyEdgeLayerMask(Array2DReal Arr, I4 NEdgesAll) const; + /// Sums the mass thickness times g from the top layer down, starting with /// the surface pressure void computePressure( From cc58ac209d31b1e6438f0e44d43b7c154fb978f2 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 05:46:55 -0500 Subject: [PATCH 12/25] Add fill-value CTests --- components/omega/test/CMakeLists.txt | 11 + components/omega/test/infra/FieldTest.cpp | 38 +-- components/omega/test/ocn/FillValueTest.cpp | 344 ++++++++++++++++++++ 3 files changed, 374 insertions(+), 19 deletions(-) create mode 100644 components/omega/test/ocn/FillValueTest.cpp diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index b30760fea72b..adb8d8fda1d3 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -523,6 +523,17 @@ add_omega_test( "-n;8" ) +################## +# Fill Value test +################## + +add_omega_test( + FILL_VALUE_TEST + testFillValue.exe + ocn/FillValueTest.cpp + "-n;8" +) + ########################## # Vert Mix test ########################## diff --git a/components/omega/test/infra/FieldTest.cpp b/components/omega/test/infra/FieldTest.cpp index 530a5436f900..9be90e8038ec 100644 --- a/components/omega/test/infra/FieldTest.cpp +++ b/components/omega/test/infra/FieldTest.cpp @@ -158,44 +158,44 @@ void initFieldTest() { DimNames[0] = "NCells"; auto Test1DI4H = Field::create("Test1DI4H", "Test 1DI4 field on host", "Units1DI4H", - "var_name_1DI4", 0, 100000, 999, 1, DimNames); + "var_name_1DI4", 0, 100000, 1, DimNames); DimNames[0] = "NEdges"; auto Test1DI8H = Field::create("Test1DI8H", "Test 1DI8 field on host", "Units1DI8H", - "var_name_1DI8", 0, 100000, 999, 1, DimNames); + "var_name_1DI8", 0, 100000, 1, DimNames); DimNames[0] = "NVertices"; auto Test1DR4H = Field::create("Test1DR4H", "Test 1DR4 field on host", "Units1DR4H", - "var_name_1DR4", 0.0, 100000.0, 999, 1, DimNames); + "var_name_1DR4", 0.0, 100000.0, 1, DimNames); DimNames[0] = "NVertLayers"; auto Test1DR8H = Field::create("Test1DR8H", "Test 1DR8 field on host", "Units1DR8H", - "var_name_1DR8", 0.0, 100000.0, 999, 1, DimNames); + "var_name_1DR8", 0.0, 100000.0, 1, DimNames); // 1D Fields on device DimNames[0] = "NCells"; auto Test1DI4 = Field::create("Test1DI4", "Test 1DI4 field on device", "Units1DI4", - "var_name_1DI4", 0, 100000, 999, 1, DimNames); + "var_name_1DI4", 0, 100000, 1, DimNames); DimNames[0] = "NEdges"; auto Test1DI8 = Field::create("Test1DI8", "Test 1DI8 field on device", "Units1DI8", - "var_name_1DI8", 0, 100000, 999, 1, DimNames); + "var_name_1DI8", 0, 100000, 1, DimNames); DimNames[0] = "NVertices"; auto Test1DR4 = Field::create("Test1DR4", "Test 1DR4 field on device", "Units1DR4", - "var_name_1DR4", 0.0, 100000.0, 999, 1, DimNames); + "var_name_1DR4", 0.0, 100000.0, 1, DimNames); DimNames[0] = "NVertLayers"; auto Test1DR8 = Field::create("Test1DR8", "Test 1DR8 field on device", "Units1DR8", - "var_name_1DR8", 0.0, 100000.0, 999, 1, DimNames); + "var_name_1DR8", 0.0, 100000.0, 1, DimNames); // 2D Fields on host DimNames.resize(2); @@ -204,25 +204,25 @@ void initFieldTest() { DimNames[1] = "NVertLayers"; auto Test2DI4H = Field::create("Test2DI4H", "Test 2DI4 field on host", "Units2DI4H", - "var_name_2DI4", 0, 100000, 999, 2, DimNames); + "var_name_2DI4", 0, 100000, 2, DimNames); DimNames[0] = "NEdges"; DimNames[1] = "NVertLayers"; auto Test2DI8H = Field::create("Test2DI8H", "Test 2DI8 field on host", "Units2DI8H", - "var_name_2DI8", 0, 100000, 999, 2, DimNames); + "var_name_2DI8", 0, 100000, 2, DimNames); DimNames[0] = "NVertices"; DimNames[1] = "NVertLayers"; auto Test2DR4H = Field::create("Test2DR4H", "Test 2DR4 field on host", "Units2DR4H", - "var_name_2DR4", 0.0, 100000.0, 999, 2, DimNames); + "var_name_2DR4", 0.0, 100000.0, 2, DimNames); DimNames[0] = "NCells"; DimNames[1] = "NVertLayers"; auto Test2DR8H = Field::create("Test2DR8H", "Test 2DR8 field on host", "Units2DR8H", - "var_name_2DR8", 0.0, 100000.0, 999, 2, DimNames); + "var_name_2DR8", 0.0, 100000.0, 2, DimNames); // 2D Fields on device @@ -230,25 +230,25 @@ void initFieldTest() { DimNames[1] = "NVertLayers"; auto Test2DI4 = Field::create("Test2DI4", "Test 2DI4 field on device", "Units2DI4", - "var_name_2DI4", 0, 100000, 999, 2, DimNames); + "var_name_2DI4", 0, 100000, 2, DimNames); DimNames[0] = "NEdges"; DimNames[1] = "NVertLayers"; auto Test2DI8 = Field::create("Test2DI8", "Test 2DI8 field on device", "Units2DI8", - "var_name_2DI8", 0, 100000, 999, 2, DimNames); + "var_name_2DI8", 0, 100000, 2, DimNames); DimNames[0] = "NVertices"; DimNames[1] = "NVertLayers"; auto Test2DR4 = Field::create("Test2DR4", "Test 2DR4 field on device", "Units2DR4", - "var_name_2DR4", 0.0, 100000.0, 999, 2, DimNames); + "var_name_2DR4", 0.0, 100000.0, 2, DimNames); DimNames[0] = "NEdges"; DimNames[1] = "NVertLayers"; auto Test2DR8 = Field::create("Test2DR8", "Test 2DR8 field on device", "Units2DR8", - "var_name_2DR8", 0.0, 100000.0, 999, 2, DimNames); + "var_name_2DR8", 0.0, 100000.0, 2, DimNames); // Higher dimension fields on device @@ -258,7 +258,7 @@ void initFieldTest() { DimNames[2] = "NTracers"; auto Test3DI4 = Field::create("Test3DI4", "Test 3DI4 field on device", "Units3DI4", - "var_name_3DI4", 0, 100000, 999, 3, DimNames); + "var_name_3DI4", 0, 100000, 3, DimNames); DimNames.resize(4); DimNames[0] = "NCells"; @@ -267,7 +267,7 @@ void initFieldTest() { DimNames[3] = "NTime"; auto Test4DI8 = Field::create("Test4DI8", "Test 4DI8 field on device", "Units4DI8", - "var_name_4DI8", 0, 100000, 999, 4, DimNames); + "var_name_4DI8", 0, 100000, 4, DimNames); DimNames.resize(5); DimNames[0] = "NCells"; @@ -277,7 +277,7 @@ void initFieldTest() { DimNames[4] = "NStuff"; auto Test5DR4 = Field::create("Test5DR4", "Test 5DR4 field on device", "Units5DR4", - "var_name_5DR4", 0, 100000, 999, 5, DimNames); + "var_name_5DR4", 0, 100000, 5, DimNames); /// Create two field groups for 1D and 2D fields auto FieldGroup1D = FieldGroup::create("FieldGroup1D"); diff --git a/components/omega/test/ocn/FillValueTest.cpp b/components/omega/test/ocn/FillValueTest.cpp new file mode 100644 index 000000000000..def12bcb4987 --- /dev/null +++ b/components/omega/test/ocn/FillValueTest.cpp @@ -0,0 +1,344 @@ +//===-- Test driver for Omega fill value standardization --------*- C++ -*-===// +// +/// \file +/// \brief Test driver for Omega fill value constants and auto-fill behavior +/// +/// Tests: +/// 1. FillValueI4/I8/R4/R8 constants match the NetCDF-C NC_FILL_* values +/// 2. Field::attachData() auto-fills the array with the declared fill value +/// 3. Inactive layers (k >= MaxLayerCell) contain FillValueReal after +/// VertCoord initialization and compute +/// 4. NormalVelocity has the correct 3-zone pattern after IC read and +/// masking: +/// fully-inactive layers = FillValueReal, boundary layers = 0, active +/// layers +/// != FillValueReal (quiescent IC gives 0, which is valid) +// +//===----------------------------------------------------------------------===// + +#include "AuxiliaryState.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "Error.h" +#include "Field.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "IO.h" +#include "IOStream.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "PGrad.h" +#include "Pacer.h" +#include "Tendencies.h" +#include "TimeMgr.h" +#include "TimeStepper.h" +#include "Tracers.h" +#include "VertCoord.h" +#include "mpi.h" + +#include + +using namespace OMEGA; + +//------------------------------------------------------------------------------ +void logNormalVelocityFillDiagnostics(const char *Stage) { + + auto *DefEnv = MachEnv::getDefault(); + auto *DefMesh = HorzMesh::getDefault(); + auto *DefVertCoord = VertCoord::getDefault(); + auto *DefState = OceanState::getDefault(); + + I4 NEdgesOwned = DefMesh->NEdgesOwned; + I4 NVertLayers = DefVertCoord->NVertLayers; + + auto NVH = createHostMirrorCopy(DefState->NormalVelocity[0]); + const auto &MinLayerEdgeTopH = DefVertCoord->MinLayerEdgeTopH; + const auto &MaxLayerEdgeBotH = DefVertCoord->MaxLayerEdgeBotH; + const auto &MinLayerEdgeBotH = DefVertCoord->MinLayerEdgeBotH; + const auto &MaxLayerEdgeTopH = DefVertCoord->MaxLayerEdgeTopH; + + int ErrActiveLocal = 0; + int NLogged = 0; + for (int IEdge = 0; IEdge < NEdgesOwned; ++IEdge) { + int KTop = MinLayerEdgeTopH(IEdge); + int KBot = std::max(0, (int)MaxLayerEdgeBotH(IEdge)); + int KMin = MinLayerEdgeBotH(IEdge); + int KMax = MaxLayerEdgeTopH(IEdge); + + for (int K = 0; K < NVertLayers; ++K) { + Real Val = NVH(IEdge, K); + bool Active = (KMax >= 0 && K >= KMin && K <= KMax); + + if (Active && Val == FillValueReal) { + ++ErrActiveLocal; + if (NLogged < 5) { + LOG_INFO("FillValueTest: {} active fill sample: IEdge={} K={} " + "KTop={} KBot={} KMin={} KMax={} Val={}", + Stage, IEdge, K, KTop, KBot, KMin, KMax, Val); + ++NLogged; + } + } + } + } + + int ErrActiveGlobal = 0; + MPI_Allreduce(&ErrActiveLocal, &ErrActiveGlobal, 1, MPI_INT, MPI_SUM, + DefEnv->getComm()); + LOG_INFO("FillValueTest: {} active NormalVelocity fill count: local={} " + "global={}", + Stage, ErrActiveLocal, ErrActiveGlobal); +} + +//------------------------------------------------------------------------------ +// Full initialization needed for OceanState (matches StateTest.cpp pattern) +void initFillValueTest() { + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + initLogging(DefEnv); + + Config("Omega"); + Config::readAll("omega.yml"); + + TimeStepper::init1(); + TimeStepper *DefStepper = TimeStepper::getDefault(); + Clock *ModelClock = DefStepper->getClock(); + + IO::init(DefComm); + + Field::init(ModelClock); + IOStream::init(ModelClock); + + Decomp::init(); + + int Err = Halo::init(); + if (Err != 0) + ABORT_ERROR("FillValueTest: error initializing default halo"); + + HorzMesh::init(ModelClock); + VertCoord::init(); + Tracers::init(); + AuxiliaryState::init(); + PressureGrad::init(); + Tendencies::init(); + TimeStepper::init2(); + + Err = OceanState::init(); + if (Err != 0) + ABORT_ERROR("FillValueTest: error initializing OceanState"); + + bool StreamsValid = IOStream::validateAll(); + if (!StreamsValid) + ABORT_ERROR("FillValueTest: error validating IO streams"); +} + +//------------------------------------------------------------------------------ +int main(int argc, char *argv[]) { + + Error ErrAll; + + MPI_Init(&argc, &argv); + Kokkos::initialize(); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + { + initFillValueTest(); + + // Read initial state so NormalVelocity has real values (quiescent IC: + // NV=0 for all active layers). Mirrors the pattern in ocnInit / + // StateTest. + { + Clock *ModelClock = TimeStepper::getDefault()->getClock(); + Metadata ReqMeta; + Error Err1 = IOStream::read("InitialState", ModelClock, ReqMeta); + CHECK_ERROR_ABORT(Err1, "FillValueTest: error reading initial state"); + logNormalVelocityFillDiagnostics("after InitialState read"); + + OceanState::getDefault()->exchangeHalo(0); + logNormalVelocityFillDiagnostics("after NormalVelocity halo exchange"); + + // Apply 3-zone mask after IC read: boundary layers -> 0, + // inactive layers -> FillValueReal, active layers unchanged. + VertCoord::getDefault()->applyEdgeLayerMask( + OceanState::getDefault()->NormalVelocity[0], + HorzMesh::getDefault()->NEdgesAll); + logNormalVelocityFillDiagnostics("after edge layer mask"); + + OceanState::getDefault()->copyToHost(0); + } + + auto *DefMesh = HorzMesh::getDefault(); + auto *DefVertCoord = VertCoord::getDefault(); + auto *DefState = OceanState::getDefault(); + + // ------------------------------------------------------------------ + // Test 1: FillValue constants match NetCDF-C NC_FILL_* values + // ------------------------------------------------------------------ + { + int Err = 0; + if (FillValueI4 != static_cast(NC_FILL_INT)) + ++Err; + if (FillValueI8 != static_cast(NC_FILL_INT64)) + ++Err; + if (static_cast(FillValueR4) != + static_cast(NC_FILL_FLOAT)) + ++Err; + if (FillValueR8 != static_cast(NC_FILL_DOUBLE)) + ++Err; + + if (Err == 0) { + LOG_INFO("FillValueTest: fill constant values PASS"); + } else { + ErrAll += Error(ErrorCode::Fail, + "FillValueTest: fill constant values FAIL"); + } + } + + // ------------------------------------------------------------------ + // Test 2: attachData() auto-fills the array with the declared fill value + // ------------------------------------------------------------------ + { + // Create a small non-distributed test dimension and field + I4 NTest = 10; + Dimension::create("FVTestDim", NTest); + auto TestField = + Field::create("FVTestField", "fill value test", "1", "", -1.0e40, + 1.0e40, 1, {"FVTestDim"}, false); + + // Allocate a host array set to 0.0 (not the fill value) + HostArray1DR8 TestArr("FVTestArr", NTest); + Kokkos::deep_copy(TestArr, 0.0); + + // attachData() should automatically fill with FillValueR8 + TestField->attachData(TestArr); + + int Err = 0; + for (int I = 0; I < NTest; ++I) { + if (TestArr(I) != FillValueR8) + ++Err; + } + + if (Err == 0) { + LOG_INFO("FillValueTest: attachData auto-fill PASS"); + } else { + ErrAll += Error(ErrorCode::Fail, + "FillValueTest: attachData auto-fill FAIL"); + } + } + + // ------------------------------------------------------------------ + // Test 3: Inactive layers (k >= MaxLayerCell) contain FillValueReal + // in a VertCoord field (GeomZMid) after initialization + // ------------------------------------------------------------------ + { + I4 NCellsOwned = DefMesh->NCellsOwned; + I4 NVertLayers = DefVertCoord->NVertLayers; + + auto GeomZMidH = createHostMirrorCopy(DefVertCoord->GeomZMid); + + int Err = 0; + for (int ICell = 0; ICell < NCellsOwned; ++ICell) { + I4 KMax = DefVertCoord->MaxLayerCellH(ICell); + for (int K = KMax; K < NVertLayers; ++K) { + if (GeomZMidH(ICell, K) != FillValueReal) + ++Err; + } + } + + if (Err == 0) { + LOG_INFO("FillValueTest: inactive layers contain fill value PASS"); + } else { + ErrAll += + Error(ErrorCode::Fail, + "FillValueTest: inactive layers contain fill value FAIL"); + } + } + + // ------------------------------------------------------------------ + // Test 4: NormalVelocity has the correct 3-zone pattern after IC read + // and applyEdgeLayerMask(): + // - Fully inactive layers (outside [MinLayerEdgeTop, MaxLayerEdgeBot]): + // FillValueReal + // - Boundary layers ([MinLayerEdgeTop, MaxLayerEdgeBot] but outside + // [MinLayerEdgeBot, MaxLayerEdgeTop]): 0 + // - Active layers ([MinLayerEdgeBot, MaxLayerEdgeTop]): + // != FillValueReal (quiescent IC gives 0, which is valid) + // ------------------------------------------------------------------ + { + I4 NEdgesOwned = DefMesh->NEdgesOwned; + I4 NVertLayers = DefVertCoord->NVertLayers; + + auto NVH = createHostMirrorCopy(DefState->NormalVelocity[0]); + const auto &MinLayerEdgeTopH = DefVertCoord->MinLayerEdgeTopH; + const auto &MaxLayerEdgeBotH = DefVertCoord->MaxLayerEdgeBotH; + const auto &MinLayerEdgeBotH = DefVertCoord->MinLayerEdgeBotH; + const auto &MaxLayerEdgeTopH = DefVertCoord->MaxLayerEdgeTopH; + + int ErrInactive = 0; // non-fill in fully-inactive zone (bad) + int ErrBoundary = 0; // non-zero in boundary zone (bad) + int ErrActive = 0; // fill value in active zone (bad) + + for (int IEdge = 0; IEdge < NEdgesOwned; ++IEdge) { + int KTop = MinLayerEdgeTopH(IEdge); + int KBot = std::max(0, (int)MaxLayerEdgeBotH(IEdge)); + int KMin = MinLayerEdgeBotH(IEdge); + int KMax = MaxLayerEdgeTopH(IEdge); // -1 for land-adjacent edges + + for (int K = 0; K < NVertLayers; ++K) { + Real Val = NVH(IEdge, K); + bool valid = (K >= KTop && K <= KBot); + bool active = (KMax >= 0 && K >= KMin && K <= KMax); + bool boundary = valid && !active; + + if (!valid && Val != FillValueReal) + ++ErrInactive; + if (boundary && Val != 0._Real) + ++ErrBoundary; + if (active && Val == FillValueReal) + ++ErrActive; + } + } + + if (ErrInactive == 0 && ErrBoundary == 0 && ErrActive == 0) { + LOG_INFO( + "FillValueTest: NormalVelocity fill/zero/real pattern PASS"); + } else { + ErrAll += Error( + ErrorCode::Fail, + "FillValueTest: NormalVelocity fill/zero/real pattern FAIL" + " (ErrInactive={}, ErrBoundary={}, ErrActive={})", + ErrInactive, ErrBoundary, ErrActive); + } + } + + // ------------------------------------------------------------------ + // Finalize + // ------------------------------------------------------------------ + OceanState::clear(); + Tendencies::clear(); + PressureGrad::clear(); + AuxiliaryState::clear(); + Tracers::clear(); + IOStream::finalize(); + TimeStepper::clear(); + HorzMesh::clear(); + VertCoord::clear(); + Field::clear(); + Dimension::clear(); + Halo::clear(); + Decomp::clear(); + MachEnv::removeAll(); + } + Pacer::finalize(); + Kokkos::finalize(); + int RetVal = ErrAll.isFail() ? 1 : 0; + MPI_Finalize(); + + return RetVal; +} +//===----------------------------------------------------------------------===// From 48e65e60f06865ac174fdbff110143653d42963c Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jun 2026 06:43:05 -0500 Subject: [PATCH 13/25] Update the docs --- components/omega/doc/devGuide/Field.md | 18 ++++++++++++++---- components/omega/doc/devGuide/Tracers.md | 18 ++++++++++++------ components/omega/doc/userGuide/Tracers.md | 18 ++++++++++++------ 3 files changed, 38 insertions(+), 16 deletions(-) diff --git a/components/omega/doc/devGuide/Field.md b/components/omega/doc/devGuide/Field.md index 85244430b133..62f2ad295d73 100644 --- a/components/omega/doc/devGuide/Field.md +++ b/components/omega/doc/devGuide/Field.md @@ -37,7 +37,6 @@ Fields are created with standard metadata using StdName, ///< [in] CF standard Name (string) ValidMin, ///< [in] min valid field value (same type as ValidMax, ///< [in] max valid field value field data) - FillValue, ///< [in] scalar used for undefined entries NumDims, ///< [in] number of dimensions (int) Dimensions, ///< [in] dim names (vector of strings) TimeDependent ///< [in] (opt, default true) if time varying @@ -49,9 +48,13 @@ not exist, an empty string can be provided. This is uncommon for most fields since the CF conventions maintain a fairly complete list, but can be the case for some intermediate calculations or unique analyses. If there is no restriction on valid range, an appropriately large range should be provided for -the data type. Similarly, if a FillValue is not being used, a very unique -number should be supplied to prevent accidentally treating valid data as a -FillValue. The optional TimeDependent argument can be omitted and is assumed +the data type. The fill value is not specified when creating a field; it is +automatically deduced from the element type of the array passed to +`attachData()` and set to the corresponding standard constant from `FillValues.h` +(`FillValueI4`, `FillValueI8`, `FillValueR4`, or `FillValueR8`). These match +the NetCDF-C `NC_FILL_*` standard values recognized by analysis tools such as +ncview, Xarray, and NCO. The optional TimeDependent argument can be omitted +and is assumed to be true by default. Fields with this attribute will be output with the unlimited time dimension added. Time should not be added explicitly in the dimension list since it will be added during I/O. Fields that do not change @@ -113,6 +116,13 @@ is provided using the field name: InDataArray // [in] Array with data to attach ); ``` +When `attachData` is called, every element of the attached array is automatically +initialized to the field's declared fill value. This ensures that inactive entries +(e.g., ocean layers below `MaxLayerCell`) hold a well-defined NetCDF-standard +sentinel in output without any explicit initialization in module code. Subsequent +compute calls overwrite active entries; the fill value persists only where no valid +data is written. + Note that the data is assumed to reside in only one location so if a mirror array exists (eg if replicated on host and device), a separate Field may be needed. However, it is is better to define only one location and allow the diff --git a/components/omega/doc/devGuide/Tracers.md b/components/omega/doc/devGuide/Tracers.md index 5bb55c1e158a..f455df08e880 100644 --- a/components/omega/doc/devGuide/Tracers.md +++ b/components/omega/doc/devGuide/Tracers.md @@ -24,17 +24,23 @@ static void defineAllTracers() { "sea_water_potential_temperature", ///< [in] CF standard Name -273.15, ///< [in] min valid field value 100.0, ///< [in] max valid field value - 1.e33, ///< [in] value for undef entries + FillValueReal, ///< [in] value for undef entries IndxTemp); ///< [out] (optional) static index - define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0, 1.e33, - IndxSalt); - define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, 1.e33); - define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, 1.e33); - define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, 1.e33); + define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0, + FillValueReal, IndxSalt); + define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, FillValueReal); + define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, FillValueReal); + define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, FillValueReal); } ``` +The `FillValueReal` constant is defined in `FillValues.h` and matches the +NetCDF-C standard fill value (`NC_FILL_DOUBLE` or `NC_FILL_FLOAT` depending on +the build precision). It is automatically in scope in `TracerDefs.inc` via +`Field.h`. New tracers should always use `FillValueReal` (or `FillValueI4` / +`FillValueI8` for integer fields) rather than hardcoded literals. + To add a new tracer, simply call the `define` function with the appropriate arguments. Index argument is optional one that allows to access the tracer data using the given tracer index variable. diff --git a/components/omega/doc/userGuide/Tracers.md b/components/omega/doc/userGuide/Tracers.md index d3fefb14af1e..18ba1c93a430 100644 --- a/components/omega/doc/userGuide/Tracers.md +++ b/components/omega/doc/userGuide/Tracers.md @@ -26,17 +26,23 @@ static void defineAllTracers() { "sea_water_potential_temperature", ///< [in] CF standard Name -273.15, ///< [in] min valid field value 100.0, ///< [in] max valid field value - 1.e33, ///< [in] value for undef entries + FillValueReal, ///< [in] value for undef entries IndxTemp); ///< [out] (optional) static index - define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0, 1.e33, - IndxSalt); - define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, 1.e33); - define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, 1.e33); - define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, 1.e33); + define("Salt", "Salinity", "psu", "sea_water_salinity", 0.0, 50.0, + FillValueReal, IndxSalt); + define("Debug1", "Debug Tracer 1", "none", "none", 0.0, 100.0, FillValueReal); + define("Debug2", "Debug Tracer 2", "none", "none", 0.0, 100.0, FillValueReal); + define("Debug3", "Debug Tracer 3", "none", "none", 0.0, 100.0, FillValueReal); } ``` +The `FillValueReal` constant is defined in `FillValues.h` and matches the +NetCDF-C standard fill value (`NC_FILL_DOUBLE` or `NC_FILL_FLOAT` depending on +the build precision). It is automatically in scope in `TracerDefs.inc` via +`Field.h`. New tracers should always use `FillValueReal` (or `FillValueI4` / +`FillValueI8` for integer fields) rather than hardcoded literals. + To add a new tracer, simply call the `define` function with the appropriate arguments. Index argument is optional one that allows to access the tracer data using the given tracer index variable. From 02cf466ab63744eb2925337fbabc8f97927027d6 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Fri, 5 Jun 2026 05:50:18 -0500 Subject: [PATCH 14/25] Enforce fill and zero values in NormalVelocity from IC Rather than assuming that the NormalVelocity read from an initial condition has correct masking, we ensure that it has fill values for inactive layers, zeros at boundary edges (both adjacent to land and to bathymetry), and leave its values unchanged for active, non-boundary layers. --- components/omega/src/ocn/OceanInit.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index ceabc0ede033..acd5093a6dbc 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -149,6 +149,15 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator OceanState *DefState = OceanState::getDefault(); I4 CurTimeLevel = 0; DefState->exchangeHalo(CurTimeLevel); + + // Enforce 3-zone edge mask on NormalVelocity: fully-inactive layers get + // FillValueReal, boundary layers (one active neighbor) get 0, active layers + // keep their IC/restart value. + VertCoord *DefVCoord = VertCoord::getDefault(); + HorzMesh *DefMesh = HorzMesh::getDefault(); + DefVCoord->applyEdgeLayerMask(DefState->getNormalVelocity(CurTimeLevel), + DefMesh->NEdgesAll); + DefState->copyToHost(CurTimeLevel); Forcing *DefForcing = Forcing::getDefault(); From 74375167f2e3c052f9ad749c5f6071c1a4e67106 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sat, 6 Jun 2026 06:22:04 -0500 Subject: [PATCH 15/25] Fix bounds on DivHU computation to avoid inactive layers --- components/omega/src/ocn/VertAdv.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index fc3fd4d99a11..0a29ec7f6b17 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -379,6 +379,7 @@ void VertAdv::computeVerticalPseudoVelocity( OMEGA_SCOPE(LocEOnC, Mesh->EdgesOnCell); OMEGA_SCOPE(LocDvE, Mesh->DvEdge); OMEGA_SCOPE(LocESOnC, Mesh->EdgeSignOnCell); + OMEGA_SCOPE(LocMaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); // Loop over all cells owned by the task parallelForOuter( @@ -402,12 +403,16 @@ void VertAdv::computeVerticalPseudoVelocity( const I4 KLen = chunkLength(KChunk, KStart, KMax); for (int J = 0; J < LocNEOnC(ICell); ++J) { - const I4 JEdge = LocEOnC(ICell, J); + const I4 JEdge = LocEOnC(ICell, J); + const I4 MaxKEdge = LocMaxLayerEdgeTop(JEdge); for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; - DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) * - FluxPseudoThickEdge(JEdge, K) * - NormalVelocity(JEdge, K) * InvAreaCell; + if (K <= MaxKEdge) { + DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) * + FluxPseudoThickEdge(JEdge, K) * + NormalVelocity(JEdge, K) * + InvAreaCell; + } } } for (int KVec = 0; KVec < KLen; ++KVec) { From 56192c71dcc8a0e88cb10d23d4f0f8f01363edb5 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sat, 6 Jun 2026 11:00:58 -0500 Subject: [PATCH 16/25] Mask NormalVelocity in time stepper We need to ensure that it is zero at active boundary edges. --- components/omega/src/timeStepping/TimeStepper.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 22a6c71a843b..334e3bb3a0fe 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -501,6 +501,10 @@ void TimeStepper::updateVelocityByTend(OceanState *State1, int TimeLevel1, CoeffSeconds * NormalVelTend(IEdge, K); }); }); + // Zero boundary layers; updateVelocityByTend only writes [KMin,KMax], so + // layers outside that range retain fill values that would corrupt + // tendencies. + VCoord->applyEdgeLayerMask(NormalVel1, Mesh->NEdgesAll); } //------------------------------------------------------------------------------ From 25c9a3274967fdf883a5ddabe1c26c05250746ab Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 22 Jun 2026 21:42:13 +0200 Subject: [PATCH 17/25] Apply suggestions from code review Co-authored-by: Maciej Waruszewski --- components/omega/src/base/FillValues.h | 14 +++++--------- components/omega/src/infra/Field.h | 2 +- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/components/omega/src/base/FillValues.h b/components/omega/src/base/FillValues.h index 20ee9593e64b..21b6caf9ed1b 100644 --- a/components/omega/src/base/FillValues.h +++ b/components/omega/src/base/FillValues.h @@ -29,15 +29,11 @@ template <> inline constexpr R8 FillValue = 9.9692099683868690e+36; ///< NC_FILL_DOUBLE /// Named aliases for readability in comparison and test code. -constexpr I4 FillValueI4 = FillValue; -constexpr I8 FillValueI8 = FillValue; -constexpr R4 FillValueR4 = FillValue; -constexpr R8 FillValueR8 = FillValue; -#if defined(SINGLE_PRECISION) -constexpr Real FillValueReal = FillValue; -#else -constexpr Real FillValueReal = FillValue; -#endif +constexpr I4 FillValueI4 = FillValue; +constexpr I8 FillValueI8 = FillValue; +constexpr R4 FillValueR4 = FillValue; +constexpr R8 FillValueR8 = FillValue; +constexpr Real FillValueReal = FillValue; } // end namespace OMEGA diff --git a/components/omega/src/infra/Field.h b/components/omega/src/infra/Field.h index a94d0eb01d47..6683ec9d1bc8 100644 --- a/components/omega/src/infra/Field.h +++ b/components/omega/src/infra/Field.h @@ -95,7 +95,7 @@ class Field { // tools (ncview, Xarray, NCO, ...); store only that, no non-standard // alias FieldMeta["_FillValue"] = FillVal; - Kokkos::deep_copy((InDataArray, FillVal); + deepCopy(InDataArray, FillVal); } public: From 3961f3d3cba9e554a5cf44dca9417d79e02a0684 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 22 Jun 2026 15:35:29 -0500 Subject: [PATCH 18/25] Apply edge layer mask to NormalVelocity in OceanState constructor --- components/omega/src/ocn/OceanState.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index b74e4df1d0c5..0e8ffe682cfe 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/components/omega/src/ocn/OceanState.cpp @@ -115,6 +115,14 @@ OceanState::OceanState( NEdgesSize, NVertLayers); } + // Apply edge layer mask to NormalVelocity for all time levels so that + // fully-inactive layers are initialized with FillValueReal and boundary + // layers with 0. + VertCoord *DefVCoord = VertCoord::getDefault(); + for (int I = 0; I < NTimeLevels; I++) { + DefVCoord->applyEdgeLayerMask(NormalVelocity[I], NEdgesAll); + } + // Register fields and metadata for IO defineFields(); From 2052fa28cd5a0cdf0789dfdc280662f673f1ceda Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 22 Jun 2026 15:47:27 -0500 Subject: [PATCH 19/25] Create applyCellLayerMask, applyVertexLayerMask --- components/omega/src/ocn/VertCoord.cpp | 41 ++++++++++++++++++++++++++ components/omega/src/ocn/VertCoord.h | 13 ++++++++ 2 files changed, 54 insertions(+) diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index f22077026e3b..07cab988b895 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -1195,6 +1195,47 @@ void VertCoord::zeroEdgeField(Array2DReal Arr, I4 NEdgesAll) const { }); } +//------------------------------------------------------------------------------ +// Enforce masking on a cell field after IC or restart read. +void VertCoord::applyCellLayerMask(Array2DReal Arr, I4 NCellsAll) const { + OMEGA_SCOPE(LocMinLayerCell, MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell); + I4 LocNVertLayers = NVertLayers; + parallelForOuter( + {NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + parallelForInner( + Team, Range{0, LocNVertLayers - 1}, INNER_LAMBDA(int K) { + if (K < KMin || K > KMax) + Arr(ICell, K) = FillValueReal; + // else: active layer - keep IC/restart value + }); + }); +} + +//------------------------------------------------------------------------------ +// Enforce masking on a vertex field after IC or restart read. +void VertCoord::applyVertexLayerMask(Array2DReal Arr, I4 NVerticesAll) const { + OMEGA_SCOPE(LocMinLayerVertexTop, MinLayerVertexTop); + OMEGA_SCOPE(LocMaxLayerVertexBot, MaxLayerVertexBot); + I4 LocNVertLayers = NVertLayers; + parallelForOuter( + {NVerticesAll}, KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) { + const int KMin = LocMinLayerVertexTop(IVertex); + const int KMax = LocMaxLayerVertexBot(IVertex); + parallelForInner( + Team, Range{0, LocNVertLayers - 1}, INNER_LAMBDA(int K) { + if (K < KMin || K > KMax) + Arr(IVertex, K) = FillValueReal; + // else: active layer - keep IC/restart value + // Unlike edges, boundary vertices (those with one or more + // active surrounding cells) hold valid, generally non-zero + // data and are not zeroed. + }); + }); +} + //------------------------------------------------------------------------------ // Enforce 3-zone masking on an edge field after IC or restart read. void VertCoord::applyEdgeLayerMask(Array2DReal Arr, I4 NEdgesAll) const { diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h index 848833e9c9fc..977c931813ca 100644 --- a/components/omega/src/ocn/VertCoord.h +++ b/components/omega/src/ocn/VertCoord.h @@ -255,6 +255,19 @@ class VertCoord { /// [MinLayerEdgeBot, MaxLayerEdgeTop] -> unchanged (IC/restart values) void applyEdgeLayerMask(Array2DReal Arr, I4 NEdgesAll) const; + /// Enforce masking on a cell field after IC or restart read: + /// K < MinLayerCell(ICell) or K > MaxLayerCell(ICell) -> FillValueReal + /// [MinLayerCell(ICell), MaxLayerCell(ICell)] -> unchanged (IC/restart + /// values) + void applyCellLayerMask(Array2DReal Arr, I4 NCellsAll) const; + + /// Enforce masking on a vertex field after IC or restart read. Unlike + /// edges, a boundary vertex with one or more active surrounding cells holds + /// valid (generally non-zero) data, so there is no zeroed boundary zone: + /// K < MinLayerVertexTop or K > MaxLayerVertexBot -> FillValueReal + /// [MinLayerVertexTop, MaxLayerVertexBot] -> unchanged (IC/restart values) + void applyVertexLayerMask(Array2DReal Arr, I4 NVerticesAll) const; + /// Sums the mass thickness times g from the top layer down, starting with /// the surface pressure void computePressure( From e145479a8d17096b3fd05e758423886f94a9e83d Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 22 Jun 2026 16:29:19 -0500 Subject: [PATCH 20/25] Use single method in OceanRun for applying fill values to OceanState variables --- components/omega/src/ocn/OceanInit.cpp | 11 ++++------- components/omega/src/ocn/OceanState.cpp | 16 ++++++++++++++++ components/omega/src/ocn/OceanState.h | 4 ++++ 3 files changed, 24 insertions(+), 7 deletions(-) diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index acd5093a6dbc..494cdb9b3776 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -150,13 +150,10 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator I4 CurTimeLevel = 0; DefState->exchangeHalo(CurTimeLevel); - // Enforce 3-zone edge mask on NormalVelocity: fully-inactive layers get - // FillValueReal, boundary layers (one active neighbor) get 0, active layers - // keep their IC/restart value. - VertCoord *DefVCoord = VertCoord::getDefault(); - HorzMesh *DefMesh = HorzMesh::getDefault(); - DefVCoord->applyEdgeLayerMask(DefState->getNormalVelocity(CurTimeLevel), - DefMesh->NEdgesAll); + // Enforce layer masks on state and tracer variables: fully-inactive layers + // get FillValueReal, boundary layers get 0, active layers keep their + // IC/restart value. + DefState->applyLayerMasks(CurTimeLevel); DefState->copyToHost(CurTimeLevel); diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index 0e8ffe682cfe..42764e97af36 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/components/omega/src/ocn/OceanState.cpp @@ -17,6 +17,7 @@ #include "MachEnv.h" #include "OmegaKokkos.h" #include "TimeStepper.h" +#include "Tracers.h" namespace OMEGA { @@ -301,6 +302,21 @@ void OceanState::copyToHost(const I4 TimeLevel) { deepCopy(NormalVelocityH[TimeIndex], NormalVelocity[TimeIndex]); } // end copyToHost +//------------------------------------------------------------------------------ +// Apply layer masks to NormalVelocity, PseudoThickness, Temperature, and +// Salinity at the given time level after reading an IC or restart file. +void OceanState::applyLayerMasks(const I4 TimeLevel) { + + VertCoord *DefVCoord = VertCoord::getDefault(); + DefVCoord->applyEdgeLayerMask(getNormalVelocity(TimeLevel), NEdgesAll); + DefVCoord->applyCellLayerMask(getPseudoThickness(TimeLevel), NCellsAll); + DefVCoord->applyCellLayerMask(Tracers::getByName(TimeLevel, "Temperature"), + NCellsAll); + DefVCoord->applyCellLayerMask(Tracers::getByName(TimeLevel, "Salinity"), + NCellsAll); + +} // end applyLayerMasks + //------------------------------------------------------------------------------ // Perform state halo exchange // TimeLevel == [1:new, 0:current, -1:previous, -2:two times ago, ...] diff --git a/components/omega/src/ocn/OceanState.h b/components/omega/src/ocn/OceanState.h index 3632449ccc32..fe3891e05b7c 100644 --- a/components/omega/src/ocn/OceanState.h +++ b/components/omega/src/ocn/OceanState.h @@ -122,6 +122,10 @@ class OceanState { /// Get normal velocity host array at given time level HostArray2DReal getNormalVelocityH(const I4 TimeLevel) const; + /// Apply layer masks to NormalVelocity, PseudoThickness, Temperature, and + /// Salinity after reading an IC or restart file + void applyLayerMasks(const I4 TimeLevel); + /// Exchange halo void exchangeHalo(const I4 TimeLevel); From 5e06ba37c7234d6916d138294573575f87c130f9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 25 Jun 2026 17:59:43 -0500 Subject: [PATCH 21/25] Fixup rebase on vmix shear changes --- components/omega/src/ocn/VertMix.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index b6ea017406fc..2b988e8ddb19 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -369,9 +369,8 @@ void VertMix::defineFields() { "sea_water_gradient_richardson_number", // CF-ish Name std::numeric_limits::min(), // Min valid value std::numeric_limits::max(), // Max valid value - FillValue, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + NDims, // Number of dimensions + DimNames // Dimension names ); /// Create and register the GradRichNumSmoothed field auto GradRichNumSmoothedField = Field::create( @@ -381,9 +380,8 @@ void VertMix::defineFields() { "sea_water_gradient_richardson_number_smoothed", // CF-ish Name std::numeric_limits::min(), // Min valid value std::numeric_limits::max(), // Max valid value - FillValue, // Scalar used for undefined entries - NDims, // Number of dimensions - DimNames // Dimension names + NDims, // Number of dimensions + DimNames // Dimension names ); // Create a field group for the vertmix-specific state fields From 35bb0b0fc1f3c49ec34221e84abd17bacfbd1bd4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 26 Jun 2026 10:06:09 -0500 Subject: [PATCH 22/25] Fixup rebase on surface stress changes --- components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp b/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp index 4375f79f63e3..f61cee3ccc26 100644 --- a/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp +++ b/components/omega/src/ocn/forcingVars/SfcStressForcingVars.cpp @@ -36,7 +36,6 @@ void SfcStressForcingVars::registerFields( "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined NDims, // number of dimensions DimNames // dim names ); @@ -48,7 +47,6 @@ void SfcStressForcingVars::registerFields( "", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValue, // scalar used undefined NDims, // number of dimensions DimNames // dimension names ); From 0e9e2342dd0881b369ee6f26aea0f1441bf26fc7 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 30 Jun 2026 06:48:33 -0500 Subject: [PATCH 23/25] Add unit tests for applyCellLayerMask and applyVertexLayerMask Add Test 5 and Test 6 to FillValueTest.cpp covering the cell and vertex layer-mask methods added in 27c02f57f7. Each test fills a synthetic field with a sentinel value, applies the mask, and verifies the resulting zones: cells expect FillValueReal in inactive layers and the sentinel in active layers; vertices expect the 3-zone pattern (fill / 0 / sentinel) matching the existing edge-mask test. Co-Authored-By: Claude Opus 4.8 --- components/omega/test/ocn/FillValueTest.cpp | 113 ++++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/components/omega/test/ocn/FillValueTest.cpp b/components/omega/test/ocn/FillValueTest.cpp index def12bcb4987..d0b268a9e721 100644 --- a/components/omega/test/ocn/FillValueTest.cpp +++ b/components/omega/test/ocn/FillValueTest.cpp @@ -13,6 +13,11 @@ /// fully-inactive layers = FillValueReal, boundary layers = 0, active /// layers /// != FillValueReal (quiescent IC gives 0, which is valid) +/// 5. applyCellLayerMask sets inactive cell layers to FillValueReal and +/// leaves active layers ([MinLayerCell, MaxLayerCell]) unchanged +/// 6. applyVertexLayerMask sets inactive vertex layers to FillValueReal and +/// leaves active layers ([MinLayerVertexTop, MaxLayerVertexBot]) +/// unchanged (no zeroed boundary zone, unlike edges) // //===----------------------------------------------------------------------===// @@ -316,6 +321,114 @@ int main(int argc, char *argv[]) { } } + // ------------------------------------------------------------------ + // Test 5: applyCellLayerMask sets inactive cell layers to FillValueReal + // and leaves active layers ([MinLayerCell, MaxLayerCell]) + // unchanged. + // A synthetic cell field is filled with a sentinel value (!= 0 and + // != FillValueReal); after masking: + // - Inactive layers (K < MinLayerCell or K > MaxLayerCell): + // FillValueReal + // - Active layers ([MinLayerCell, MaxLayerCell]): sentinel (unchanged) + // ------------------------------------------------------------------ + { + I4 NCellsOwned = DefMesh->NCellsOwned; + I4 NVertLayers = DefVertCoord->NVertLayers; + + const Real Sentinel = 42.0_Real; + + Array2DReal SynthCell("SynthCell", DefMesh->NCellsSize, NVertLayers); + deepCopy(SynthCell, Sentinel); + DefVertCoord->applyCellLayerMask(SynthCell, DefMesh->NCellsAll); + + auto SynthCellH = createHostMirrorCopy(SynthCell); + const auto &MinLayerCellH = DefVertCoord->MinLayerCellH; + const auto &MaxLayerCellH = DefVertCoord->MaxLayerCellH; + + int ErrInactive = 0; // non-fill in inactive zone (bad) + int ErrActive = 0; // not sentinel in active zone (bad) + + for (int ICell = 0; ICell < NCellsOwned; ++ICell) { + int KMin = MinLayerCellH(ICell); + int KMax = MaxLayerCellH(ICell); // -1 for fully-land cells + + for (int K = 0; K < NVertLayers; ++K) { + Real Val = SynthCellH(ICell, K); + bool active = (K >= KMin && K <= KMax); + + if (active && Val != Sentinel) + ++ErrActive; + if (!active && Val != FillValueReal) + ++ErrInactive; + } + } + + if (ErrInactive == 0 && ErrActive == 0) { + LOG_INFO("FillValueTest: applyCellLayerMask pattern PASS"); + } else { + ErrAll += Error(ErrorCode::Fail, + "FillValueTest: applyCellLayerMask pattern FAIL" + " (ErrInactive={}, ErrActive={})", + ErrInactive, ErrActive); + } + } + + // ------------------------------------------------------------------ + // Test 6: applyVertexLayerMask sets inactive vertex layers to + // FillValueReal and leaves active layers + // ([MinLayerVertexTop, MaxLayerVertexBot]) unchanged. + // Unlike edges, a boundary vertex with one or more active surrounding + // cells holds valid (generally non-zero) data, so there is no zeroed + // boundary zone. + // A synthetic vertex field is filled with a sentinel value (!= 0 and + // != FillValueReal); after masking: + // - Inactive layers (outside [MinLayerVertexTop, MaxLayerVertexBot]): + // FillValueReal + // - Active layers ([MinLayerVertexTop, MaxLayerVertexBot]): sentinel + // (unchanged) + // ------------------------------------------------------------------ + { + I4 NVerticesOwned = DefMesh->NVerticesOwned; + I4 NVertLayers = DefVertCoord->NVertLayers; + + const Real Sentinel = 42.0_Real; + + Array2DReal SynthVtx("SynthVtx", DefMesh->NVerticesSize, NVertLayers); + deepCopy(SynthVtx, Sentinel); + DefVertCoord->applyVertexLayerMask(SynthVtx, DefMesh->NVerticesAll); + + auto SynthVtxH = createHostMirrorCopy(SynthVtx); + const auto &MinLayerVertexTopH = DefVertCoord->MinLayerVertexTopH; + const auto &MaxLayerVertexBotH = DefVertCoord->MaxLayerVertexBotH; + + int ErrInactive = 0; // non-fill in inactive zone (bad) + int ErrActive = 0; // not sentinel in active zone (bad) + + for (int IVertex = 0; IVertex < NVerticesOwned; ++IVertex) { + int KMin = MinLayerVertexTopH(IVertex); + int KMax = MaxLayerVertexBotH(IVertex); // -1 for invalid vertices + + for (int K = 0; K < NVertLayers; ++K) { + Real Val = SynthVtxH(IVertex, K); + bool active = (K >= KMin && K <= KMax); + + if (active && Val != Sentinel) + ++ErrActive; + if (!active && Val != FillValueReal) + ++ErrInactive; + } + } + + if (ErrInactive == 0 && ErrActive == 0) { + LOG_INFO("FillValueTest: applyVertexLayerMask pattern PASS"); + } else { + ErrAll += Error(ErrorCode::Fail, + "FillValueTest: applyVertexLayerMask pattern FAIL" + " (ErrInactive={}, ErrActive={})", + ErrInactive, ErrActive); + } + } + // ------------------------------------------------------------------ // Finalize // ------------------------------------------------------------------ From 621e04f53fb7e4e44b35716023526be9a36c5e8f Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 30 Jun 2026 06:50:43 -0500 Subject: [PATCH 24/25] Update FillValues design doc to match implementation Document applyCellLayerMask and applyVertexLayerMask alongside the edge helpers in Section 4.5, and describe the OceanState::applyLayerMasks driver called from ocnInit. Correct stale facts in Section 5: the test lives at test/ocn/FillValueTest.cpp, uses direct comparison with Error accumulation, and now includes cell and vertex mask cases (5.5, 5.6); remove the unimplemented NetCDF attribute test. Update Requirement 2.7 for the expanded cell/vertex coverage. Co-Authored-By: Claude Opus 4.8 --- components/omega/doc/design/FillValues.md | 126 +++++++++++++++------- 1 file changed, 89 insertions(+), 37 deletions(-) diff --git a/components/omega/doc/design/FillValues.md b/components/omega/doc/design/FillValues.md index f5395438d4d6..20a6958a1ed8 100644 --- a/components/omega/doc/design/FillValues.md +++ b/components/omega/doc/design/FillValues.md @@ -93,8 +93,10 @@ edges shallower than the cell — see Section 4.5 for the resolution. ### 2.7 Requirement: CTests -CTests must verify Requirements 2.1–2.3, 2.5, and Desired Requirement 2.6 for -`NormalVelocity`. +CTests must verify Requirements 2.1–2.3, 2.5, and Desired Requirement 2.6. The +three-zone boundary behavior (2.6) is exercised for `NormalVelocity` via the edge +mask; the two-zone cell and vertex masks are verified with synthetic cell and +vertex fields. ## 3 Algorithmic Formulation @@ -345,11 +347,36 @@ three-zone mask on an edge field after IC or restart read: layers outside `[MinLayerEdgeTop, MaxLayerEdgeBot]` are set to `FillValueReal`; layers inside `[MinLayerEdgeTop, MaxLayerEdgeBot]` but outside the active range `[MinLayerEdgeBot, MaxLayerEdgeTop]` are set to 0; active layers are left -unchanged. Called in `ocnInit` (`OceanInit.cpp`) for -`NormalVelocity[CurTimeLevel]` after `exchangeHalo` and before `copyToHost`. -`NormalVelocity` is the only edge-based state field read from IC/restart; all -other IC fields (pseudo-thickness, tracers, VertCoord fields) are cell-based -and need no edge masking. +unchanged. + +**`applyCellLayerMask(Array2DReal &Arr, I4 NCellsAll)`** — enforces the +two-zone mask on a cell field after IC or restart read: layers outside the +active range `[MinLayerCell, MaxLayerCell]` are set to `FillValueReal`; active +layers are left unchanged. Cell fields have no boundary (zero) zone because a +cell column is either active or inactive at a given layer — there is no +neighbor-dependent partial-activity range as there is for edges and vertices. + +**`applyVertexLayerMask(Array2DReal &Arr, I4 NVerticesAll)`** — enforces the +two-zone mask on a vertex field: layers outside `[MinLayerVertexTop, +MaxLayerVertexBot]` are set to `FillValueReal`; active layers are left +unchanged. Unlike an edge, a vertex does **not** have a zeroed boundary zone: a +boundary vertex with one or more active surrounding cells holds valid, generally +non-zero data (for example, relative vorticity at such a vertex is computed from +its active surrounding edges — see `VorticityAuxVars::computeVarsOnVertex`, which +loops over the full `[MinLayerVertexTop, MaxLayerVertexBot]` range). Zeroing that +boundary band would discard real signal, so the whole valid range is kept. There +is currently no vertex-based state field read from IC/restart, so this method has +no production caller yet; it is provided for completeness and is exercised by the +CTest (Section 5.6). + +All three `apply*LayerMask` methods use the inclusive `[Min, Max]` layer +convention consistent with `computeGeomZHeight`/`computePressure`. They are +driven through `OceanState::applyLayerMasks(TimeLevel)` +(`src/ocn/OceanState.cpp`), which is called from `ocnInit` (`OceanInit.cpp`) +after `exchangeHalo` and before `copyToHost`. `applyLayerMasks` applies the +edge mask to `NormalVelocity` and the cell mask to `PseudoThickness`, +`Temperature`, and `Salinity` — the state and tracer fields read from +IC/restart. Fields that are *not* fluxes — `MeanPseudoThickEdge`, `FluxPseudoThickEdge` — are genuinely undefined at boundary edges (interpolating thickness when one @@ -377,51 +404,44 @@ cell-level flux-accumulation kernel, and avoids the need to zero ## 5 Verification and Testing -A new test `test/infra/FillValueTest.cpp` is added and registered as -`FILL_VALUE_TEST` with 8 MPI tasks using `add_omega_test()` in the test -`CMakeLists.txt`. The test follows the standard Omega test pattern: initialize -`MachEnv`, logging, IO, `Decomp`, and `Dimension` before running assertions. +A test `test/ocn/FillValueTest.cpp` is added and registered as `FILL_VALUE_TEST` +with 8 MPI tasks (building `testFillValue.exe`) using `add_omega_test()` in the +test `CMakeLists.txt`. The test performs a full ocean initialization (matching +the `StateTest` pattern: `MachEnv`, logging, Config, `TimeStepper`, IO, `Field`, +`IOStream`, `Decomp`, `Halo`, `HorzMesh`, `VertCoord`, `Tracers`, +`AuxiliaryState`, `PressureGrad`, `Tendencies`, `OceanState`), reads the initial +state, exchanges halos, and applies `applyEdgeLayerMask` to `NormalVelocity` +before running its assertions. Each test counts mismatches and, on any failure, +accumulates an `Error(ErrorCode::Fail, ...)` into a top-level `Error` whose +final state sets the process return code. ### 5.1 Test: fill constant values -Verify that each Omega fill value constant exactly equals its NetCDF-C counterpart: - -```c++ -TstEval("FillValueI4 == NC_FILL_INT", FillValueI4, (I4)NC_FILL_INT, Err); -TstEval("FillValueI8 == NC_FILL_INT64", FillValueI8, (I8)NC_FILL_INT64, Err); -TstEval("FillValueR4 == NC_FILL_FLOAT", FillValueR4, (R4)NC_FILL_FLOAT, Err); -TstEval("FillValueR8 == NC_FILL_DOUBLE",FillValueR8, (R8)NC_FILL_DOUBLE, Err); -``` +Verify that each Omega fill value constant exactly equals its NetCDF-C +counterpart (`FillValueI4 == NC_FILL_INT`, `FillValueI8 == NC_FILL_INT64`, +`FillValueR4 == NC_FILL_FLOAT`, `FillValueR8 == NC_FILL_DOUBLE`), comparing +directly against the `` `NC_FILL_*` macros. Tests requirements: 2.1, 2.2, 2.7 ### 5.2 Test: attachData auto-fill -Create a `Field`, allocate a Kokkos host or device array whose elements are all set -to a distinct sentinel (e.g., 0), then call `attachData()`. Verify that after the -call every element of the array equals `FillValueR8` (or the appropriate type). -Repeat for I4, R4, and R8 types. +Create a `Field`, allocate a Kokkos host array whose elements are all set to a +distinct sentinel (`0`), then call `attachData()`. Verify that after the call +every element of the array equals `FillValueR8`. Tests requirements: 2.3, 2.7 -### 5.3 Test: inactive layers contain fill values after compute - -Using the smallest test mesh (the `planar` mesh used by `VertCoordTest`), -run `VertCoord::computePressure()` or `computeGeomZHeight()`. For each cell, verify -that all entries at depth index `k > MaxLayerCell(ICell)` equal `FillValueR8`. - -Tests requirements: 2.5, 2.7 - -### 5.4 Test: NetCDF output contains correct fill value attribute +### 5.3 Test: inactive layers contain fill values -Write a field to a test NetCDF file via `IOStream`. Read back the file's `_FillValue` -attribute for the variable and verify it equals `NC_FILL_DOUBLE` (or `NC_FILL_FLOAT` -for reduced-precision fields). Also read back the array entries for inactive layers -and verify they equal the fill value. +After `VertCoord` initialization, verify that the cell field `GeomZMid` carries +`FillValueReal` in every inactive layer (`k >= MaxLayerCell(ICell)`) for all +owned cells. This confirms that the auto-fill applied at `attachData()` time +survives initialization for inactive layers. Tests requirements: 2.5, 2.7 -### 5.5 Test: NormalVelocity three-zone fill/zero/real pattern +### 5.4 Test: NormalVelocity three-zone fill/zero/real pattern After full state initialization, IC read, halo exchange, and `VertCoord::applyEdgeLayerMask()`, verify that `NormalVelocity` satisfies the @@ -440,3 +460,35 @@ was changed) and (b) IC files that store zeros in inactive layers, which `applyEdgeLayerMask` must replace with `FillValueReal`. Tests requirements: 2.5, 2.6, 2.7 + +### 5.5 Test: applyCellLayerMask two-zone pattern + +Allocate a synthetic cell field, fill it with a sentinel value (≠ 0 and +≠ `FillValueReal`), apply `VertCoord::applyCellLayerMask()`, and verify for every +owned cell that: + +- **Inactive layers** (outside `[MinLayerCell, MaxLayerCell]`): equal + `FillValueReal`. +- **Active layers** (`[MinLayerCell, MaxLayerCell]`): retain the sentinel value. + +Using a sentinel rather than IC data gives an exact check that the mask +overwrites only the inactive layers. + +Tests requirements: 2.5, 2.7 + +### 5.6 Test: applyVertexLayerMask two-zone pattern + +Allocate a synthetic vertex field, fill it with the same sentinel value, apply +`VertCoord::applyVertexLayerMask()`, and verify for every owned vertex that: + +- **Inactive layers** (outside `[MinLayerVertexTop, MaxLayerVertexBot]`): equal + `FillValueReal`. +- **Active layers** (`[MinLayerVertexTop, MaxLayerVertexBot]`): retain the + sentinel value. + +Vertices have no zeroed boundary zone (a boundary vertex holds valid data), so +the whole valid range is checked for the sentinel. Because no vertex-based state +field is read from IC/restart, this synthetic test is the primary coverage for +`applyVertexLayerMask`. + +Tests requirements: 2.5, 2.7 From b5695022a35475315c9db6b7d6e629873a835d17 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 30 Jun 2026 06:51:55 -0500 Subject: [PATCH 25/25] Document VertCoord layer masking in dev and user guides Add a Layer masking subsection to the VertCoord devGuide describing zeroEdgeField, applyEdgeLayerMask, applyCellLayerMask, and applyVertexLayerMask, their zones, the inclusive active-layer convention, and the applyLayerMasks driver. Add a userGuide note explaining that inactive layers carry the fill value and boundary layers carry 0 in output. Co-Authored-By: Claude Opus 4.8 --- components/omega/doc/devGuide/VertCoord.md | 22 +++++++++++++++++++++ components/omega/doc/userGuide/VertCoord.md | 9 +++++++++ 2 files changed, 31 insertions(+) diff --git a/components/omega/doc/devGuide/VertCoord.md b/components/omega/doc/devGuide/VertCoord.md index 90bb3e017c84..f3da8dbbd0b3 100644 --- a/components/omega/doc/devGuide/VertCoord.md +++ b/components/omega/doc/devGuide/VertCoord.md @@ -128,3 +128,25 @@ This `parallel_for` iterates over vertical chunks to facilitate vectorization on The vector length on GPUs is set to 1 to maximize parallelism. The `computeGeopotential` method uses hierarchical parallelism in a very similar way to `computeTargetThickness`, except that it doesn't require a column sum. It has an outer `parallel_for` that splits horizontal cells into teams and an inner `parallel_for` that does vertical computations in chunks. + +### Layer masking + +Field arrays are auto-filled with `FillValueReal` when they are attached to a `Field` (see the [Fill Values design](../design/FillValues.md)). +Compute kernels only write to active layers, so inactive layers naturally retain the fill value. +After an initial-condition or restart read, however, state and tracer arrays come from the input file and may hold zeros (or arbitrary values) in inactive and boundary layers. +`VertCoord` provides a family of helper methods that enforce the correct layer pattern on such arrays. +All of them use the inclusive `[Min, Max]` active-layer convention consistent with `computeGeomZHeight`/`computePressure`, and use hierarchical parallelism (an outer `parallelForOuter` over the horizontal dimension and an inner `parallelForInner` over the vertical layers). + +| Method | Element | Zones applied | +|--------|---------|---------------| +| `zeroEdgeField(Arr, NEdgesAll)` | edge | sets `[MinLayerEdgeTop, MaxLayerEdgeBot]` to 0; layers outside retain their fill value | +| `applyEdgeLayerMask(Arr, NEdgesAll)` | edge | three-zone: outside `[MinLayerEdgeTop, MaxLayerEdgeBot]` → `FillValueReal`; inside that but outside `[MinLayerEdgeBot, MaxLayerEdgeTop]` → 0; active layers unchanged | +| `applyCellLayerMask(Arr, NCellsAll)` | cell | two-zone: outside `[MinLayerCell, MaxLayerCell]` → `FillValueReal`; active layers unchanged | +| `applyVertexLayerMask(Arr, NVerticesAll)` | vertex | two-zone: outside `[MinLayerVertexTop, MaxLayerVertexBot]` → `FillValueReal`; active layers unchanged | + +Cell fields have no boundary (zero) zone because a cell column is either active or inactive at a given layer, with no neighbor-dependent partial-activity range. +Vertex fields also have no zeroed boundary zone, but for a different reason: a boundary vertex with one or more active surrounding cells holds valid, generally non-zero data (for example, relative vorticity computed from its active surrounding edges over the full `[MinLayerVertexTop, MaxLayerVertexBot]` range), so the entire valid range is kept rather than zeroed. +This is unlike a flux-type edge field, where the normal flux through a face bordering land genuinely vanishes. +The flux-type edge fields handled by `zeroEdgeField` (for example `NormalVelocityTend` and `VelocityDel2Aux.Del2Edge`) are zeroed before being recomputed each time step so that boundary edges carry 0 rather than `FillValueReal`. +The `apply*LayerMask` methods are driven through `OceanState::applyLayerMasks`, which is called once from `ocnInit` after the IC/restart read: it applies the edge mask to `NormalVelocity` and the cell mask to `PseudoThickness`, `Temperature`, and `Salinity`. +`applyVertexLayerMask` has no production caller yet because no vertex-based state field is read from the IC/restart; it is provided for completeness and verified by the fill-value CTest. diff --git a/components/omega/doc/userGuide/VertCoord.md b/components/omega/doc/userGuide/VertCoord.md index 21d197efd6c6..7f9eb1b0b59d 100644 --- a/components/omega/doc/userGuide/VertCoord.md +++ b/components/omega/doc/userGuide/VertCoord.md @@ -46,3 +46,12 @@ Multiple instances of the vertical coordinate class can be created and accessed | VertCoordMovementWeights | weights to specify how total column thickness changes are distributed across layers | - | | RefPseudoThickness | reference pseudo-thickness used to distribute total column thickness changes | m | | BottomGeomDepth | positive down distance from the reference geoid to the bottom | m | + +### Inactive and boundary layers in output + +Each output field declares a fill value (matching the NetCDF-C standard `_FillValue`) that marks entries outside the active domain. +In a cell column only the active layers (`[MinLayerCell, MaxLayerCell]`) carry valid data; deeper, inactive layers contain the fill value. +Edge and vertex fields have an intermediate boundary range where some, but not all, of the surrounding cells are active. +For flux-type edge fields (such as `NormalVelocity`), the flux through a face bordering land vanishes, so those boundary layers carry `0`; layers fully outside the valid range carry the fill value. +Vertex fields, by contrast, keep their computed (generally non-zero) values throughout the valid range, since a boundary vertex still has active surrounding cells contributing to quantities like vorticity; only layers outside the valid range carry the fill value. +After an initial-condition or restart read, `VertCoord` enforces this pattern on the state and tracer fields so that analysis tools (ncview, Xarray, NCO, etc.) recognize and mask the inactive entries automatically.