Skip to content

Commit 8fc6203

Browse files
xylarxylar
committed
Add a design document for fixing fill values
1 parent 1a7c657 commit 8fc6203

2 files changed

Lines changed: 271 additions & 0 deletions

File tree

‎components/omega/doc/design/FillValues.md‎

Lines changed: 270 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,270 @@
1+
(omega-design-fill-values)=
2+
# Fill Values
3+
4+
**Table of Contents**
5+
1. [Overview](#1-overview)
6+
2. [Requirements](#2-requirements)
7+
3. [Algorithmic Formulation](#3-algorithmic-formulation)
8+
4. [Design](#4-design)
9+
5. [Verification and Testing](#5-verification-and-testing)
10+
11+
## 1 Overview
12+
13+
Omega supports output of ocean fields to NetCDF files via I/O streams. Each field
14+
carries a declared fill value stored as metadata under the `_FillValue` attribute,
15+
which marks entries that are undefined or outside the active domain. Currently, fill
16+
value constants are declared locally in each module with inconsistent values
17+
(typically `-9.99e30` for real-valued fields, `-999` for integers). These values do
18+
not match the NetCDF-C standard fill values that analysis tools and visualization
19+
software commonly expect. Additionally, the Kokkos arrays backing these fields are
20+
never explicitly initialized with fill values: inactive ocean layers (model layers
21+
below `MaxLayerCell`) hold uninitialized memory rather than a well-defined sentinel,
22+
and the distinction between "no data" and "zero" is lost.
23+
24+
This design standardizes fill values across all Omega fields to match the NetCDF-C
25+
standard, ensures that field arrays are automatically initialized with the correct
26+
fill value when data is attached to a `Field`, and verifies through CTests that the
27+
resulting behavior is correct for inactive layers and for edge/vertex fields at
28+
domain boundaries.
29+
30+
## 2 Requirements
31+
32+
### 2.1 Requirement: Centralized fill value constants
33+
34+
Fill value constants must be defined in a single location in Omega that is accessible
35+
to all modules. Scatter of locally-declared constants has produced inconsistencies
36+
(e.g., `-9.99e30` vs `-9.99E+30`) and makes future maintenance error-prone.
37+
38+
### 2.2 Requirement: Standard fill values matching NetCDF-C
39+
40+
Fill value constants must exactly match the NetCDF-C standard values (`NC_FILL_*`)
41+
as exposed by the SCORPIO library (`PIO_FILL_*` in `pio.h`). Using standard values
42+
ensures compatibility with NetCDF-aware tools (ncview, Xarray, NCO, etc.) that
43+
recognize and handle fill values automatically. The required values are:
44+
45+
| Type | Omega name | Value |
46+
|------|-----------------|-------------------------------|
47+
| I4 | `FillValueI4` | -2147483647 (`PIO_FILL_INT`) |
48+
| I8 | `FillValueI8` | `PIO_FILL_INT64` |
49+
| R4 | `FillValueR4` | 9.9692099683868690e+36f (`PIO_FILL_FLOAT`) |
50+
| R8 | `FillValueR8` | 9.9692099683868690e+36 (`PIO_FILL_DOUBLE`) |
51+
| Real | `FillValueReal` | Alias for `FillValueR4` or `FillValueR8` depending on build |
52+
53+
### 2.3 Requirement: Automatic array initialization at attach time
54+
55+
When `Field::attachData()` is called, the attached Kokkos array must be
56+
automatically initialized with the field's declared fill value before any compute
57+
routine runs. This guarantees that inactive ocean layers (below `MaxLayerCell`) and
58+
any other undefined entries contain a well-defined sentinel in output, without
59+
requiring each module to manually perform the fill.
60+
61+
### 2.4 Requirement: Standardized Field::create() calls
62+
63+
All calls to `Field::create()` across Omega must use the centralized fill value
64+
constants rather than locally declared values. This eliminates the scattered local
65+
declarations and ensures consistency.
66+
67+
### 2.5 Requirement: Inactive layers contain fill values in output
68+
69+
After compute routines execute, model layers with index `k > MaxLayerCell(ICell)`
70+
must contain the declared fill value in output. With automatic initialization at
71+
attach time (Requirement 2.3) and compute routines that only write to active layers,
72+
this is satisfied without changes to individual compute functions.
73+
74+
### 2.6 Desired: Valid values at domain boundary edges and vertices
75+
76+
Edge and vertex fields at valid boundaries adjacent to land or to bathymetry must
77+
contain computed valid values, not fill values. Switching from zero-initialization
78+
to fill-value initialization must not silently break any boundary computation that
79+
previously assumed zero defaults. This mirrors a problem encountered in MPAS-Ocean
80+
when initializing `NormalVelocity` with a fill value rather than zero: code that
81+
implicitly assumed zero at land boundaries produced incorrect results. During
82+
implementation, all compute routines for edge/vertex fields must be audited to
83+
confirm that every valid boundary entry receives an explicit computed value.
84+
85+
### 2.7 Requirement: CTests
86+
87+
CTests must verify Requirements 2.1–2.3, 2.5, and Desired Requirement 2.6 for
88+
`NormalVelocity`.
89+
90+
## 3 Algorithmic Formulation
91+
92+
No new numerical algorithm is introduced. The key behavioral change is in
93+
`Field::attachData<T>()`, which gains a fill step equivalent to:
94+
95+
```c++
96+
Kokkos::deep_copy(InDataArray, fill_scalar);
97+
```
98+
99+
where `fill_scalar` is the typed fill value extracted from the field's `FieldMeta`
100+
map (stored as `std::any` under key `"FillValue"`), and the element type
101+
`typename T::value_type` of the Kokkos view determines the `std::any_cast` type.
102+
`Kokkos::deep_copy` with a scalar source broadcasts that value to every element of
103+
the view and dispatches correctly for both host and device memory spaces.
104+
105+
## 4 Design
106+
107+
### 4.1 Data types and parameters
108+
109+
#### 4.1.1 Constants: FillValue definitions in IO.h
110+
111+
The five fill value constants are added to
112+
`components/omega/src/base/IO.h`, which already includes `pio.h` from SCORPIO.
113+
`IO.h` is the appropriate home because it already owns all PIO type mappings
114+
(`IOTypeI4`, `IOTypeR8`, etc.) and `Field.h` includes it transitively, making the
115+
constants available throughout Omega without new include directives in most modules.
116+
117+
```c++
118+
// Standard fill values — wrap SCORPIO PIO_FILL_* constants (which wrap NC_FILL_*)
119+
constexpr I4 FillValueI4 = PIO_FILL_INT; // -2147483647
120+
constexpr I8 FillValueI8 = PIO_FILL_INT64; // per NC_FILL_INT64
121+
constexpr R4 FillValueR4 = PIO_FILL_FLOAT; // 9.9692099683868690e+36
122+
constexpr R8 FillValueR8 = PIO_FILL_DOUBLE; // 9.9692099683868690e+36
123+
#if defined(SINGLE_PRECISION)
124+
constexpr Real FillValueReal = FillValueR4;
125+
#else
126+
constexpr Real FillValueReal = FillValueR8;
127+
#endif
128+
```
129+
130+
#### 4.1.2 Class/struct changes
131+
132+
No new classes or structs are introduced. The changes are confined to:
133+
- constants added to `IO.h` (Section 4.1.1)
134+
- the `attachData` template methods in `Field.h` / `Field.cpp` (Section 4.2)
135+
- removal of local fill value declarations in module files (Section 4.3)
136+
137+
### 4.2 Methods
138+
139+
#### 4.2.1 Field::attachData() auto-fill
140+
141+
The `attachData<T>()` instance method and the `attachFieldData<T>()` static helper
142+
in `components/omega/src/infra/Field.h` / `Field.cpp` gain a private fill step
143+
after the array is stored:
144+
145+
```c++
146+
template <typename T>
147+
void attachData(const T &InDataArray) {
148+
OMEGA_ASSERT(isKokkosArray<T>,
149+
"Field::attachData requires Kokkos array as input");
150+
DataArray = std::make_shared<T>(InDataArray);
151+
DataType = checkArrayType<T>();
152+
MemLoc = findArrayMemLoc<T>();
153+
fillWithValue<T>(InDataArray); // auto-fill with declared FillValue
154+
}
155+
```
156+
157+
The private helper `fillWithValue<T>()` performs the following steps:
158+
159+
1. Reads the `"FillValue"` entry from `FieldMeta` as `std::any`.
160+
2. Uses `typename T::value_type` to `std::any_cast` to the correct scalar type
161+
(I4, I8, R4, or R8).
162+
3. Calls `Kokkos::deep_copy(InDataArray, scalar)` to broadcast the scalar to every
163+
element of the view. This works for host views (`HostArray*`), device views
164+
(`Array*`), and all dimensionalities (1D–5D) without additional dispatch.
165+
4. Only the view being attached is filled; no assumption is made about whether the
166+
other memory space (host or device) has been attached.
167+
168+
Ordering invariants preserved by this design:
169+
- `Field::create()` is always called before `attachData()`, so the fill value is
170+
always present in `FieldMeta` at fill time.
171+
- Restart reads happen after `attachData()` in the existing initialization workflow
172+
and overwrite fill entries with data from the restart file.
173+
- Compute routines run after initialization and overwrite active entries, leaving
174+
inactive entries (below `MaxLayerCell`) at the fill value.
175+
176+
#### 4.2.2 No per-module manual initialization required
177+
178+
Because `attachData()` handles initialization automatically, no module needs to
179+
call an explicit fill routine. Modules may still call `Kokkos::deep_copy` or
180+
similar to set specific values, but this is independent of fill-value management.
181+
182+
### 4.3 Module-level changes: standardized Field::create() calls
183+
184+
All modules that declare local fill value constants are updated to remove those
185+
declarations and reference the centralized constants from `IO.h`. The pattern is:
186+
187+
```c++
188+
// BEFORE — local declaration (approximately 20 files)
189+
const Real FillValue = -9.99e30;
190+
// or
191+
const I4 FillValueI4 = -999;
192+
193+
// AFTER — use centralized constants from IO.h (no local declaration needed)
194+
// FillValueReal, FillValueI4, etc. are already in scope via IO.h
195+
```
196+
197+
Files that require updates (representative list; full list determined during
198+
implementation by searching for local fill value patterns):
199+
200+
- `src/ocn/VertMix.cpp`, `VertAdv.cpp`, `VertCoord.cpp`
201+
- `src/ocn/Eos.cpp`, `Tendencies.cpp`, `OceanState.cpp`, `HorzMesh.cpp`
202+
- `src/ocn/WindForcingAuxVars.cpp`, `VorticityAuxVars.cpp`, `KineticAuxVars.cpp`
203+
- `src/ocn/TracerAuxVars.cpp`, `PseudoThicknessAuxVars.cpp`,
204+
`SurfTracerRestAuxVars.cpp`, `VelocityDel2AuxVars.cpp`
205+
206+
### 4.4 Boundary edge and vertex handling
207+
208+
During implementation, all compute routines that populate fields defined on mesh
209+
edges or vertices must be audited to confirm that every valid boundary entry —
210+
including edges and vertices adjacent to land cells and edges at the base of the
211+
water column — receives an explicit computed value rather than relying on
212+
zero-initialization. Where such reliance is found, an explicit assignment must be
213+
added. The CTest in Section 5.5 detects regressions after the change.
214+
215+
## 5 Verification and Testing
216+
217+
A new test `test/infra/FillValueTest.cpp` is added and registered as
218+
`FILL_VALUE_TEST` with 8 MPI tasks using `add_omega_test()` in the test
219+
`CMakeLists.txt`. The test follows the standard Omega test pattern: initialize
220+
`MachEnv`, logging, IO, `Decomp`, and `Dimension` before running assertions.
221+
222+
### 5.1 Test: fill constant values
223+
224+
Verify that each Omega fill value constant exactly equals its NetCDF-C counterpart:
225+
226+
```c++
227+
TstEval<I4>("FillValueI4 == NC_FILL_INT", FillValueI4, (I4)NC_FILL_INT, Err);
228+
TstEval<I8>("FillValueI8 == NC_FILL_INT64", FillValueI8, (I8)NC_FILL_INT64, Err);
229+
TstEval<R4>("FillValueR4 == NC_FILL_FLOAT", FillValueR4, (R4)NC_FILL_FLOAT, Err);
230+
TstEval<R8>("FillValueR8 == NC_FILL_DOUBLE",FillValueR8, (R8)NC_FILL_DOUBLE, Err);
231+
```
232+
233+
Tests requirements: 2.1, 2.2, 2.7
234+
235+
### 5.2 Test: attachData auto-fill
236+
237+
Create a `Field`, allocate a Kokkos host or device array whose elements are all set
238+
to a distinct sentinel (e.g., 0), then call `attachData()`. Verify that after the
239+
call every element of the array equals `FillValueR8` (or the appropriate type).
240+
Repeat for I4, R4, and R8 types.
241+
242+
Tests requirements: 2.3, 2.7
243+
244+
### 5.3 Test: inactive layers contain fill values after compute
245+
246+
Using the smallest test mesh (the `planar` mesh used by `VertCoordTest`),
247+
run `VertCoord::computePressure()` or `computeGeomZHeight()`. For each cell, verify
248+
that all entries at depth index `k > MaxLayerCell(ICell)` equal `FillValueR8`.
249+
250+
Tests requirements: 2.5, 2.7
251+
252+
### 5.4 Test: NetCDF output contains correct fill value attribute
253+
254+
Write a field to a test NetCDF file via `IOStream`. Read back the file's `_FillValue`
255+
attribute for the variable and verify it equals `NC_FILL_DOUBLE` (or `NC_FILL_FLOAT`
256+
for reduced-precision fields). Also read back the array entries for inactive layers
257+
and verify they equal the fill value.
258+
259+
Tests requirements: 2.5, 2.7
260+
261+
### 5.5 Test: NormalVelocity at land boundaries contains valid values
262+
263+
After full state initialization on a test mesh that contains at least one land cell,
264+
verify that `NormalVelocity` at all valid boundary edges adjacent to land contains
265+
a computed value (expected to be zero for a quiescent initial state) and not the
266+
fill value. This guards against the class of regression where switching from
267+
zero-initialization to fill-value initialization silently corrupts boundary
268+
computations.
269+
270+
Tests requirements: 2.6, 2.7

‎components/omega/doc/index.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,7 @@ design/Decomp
115115
design/Driver
116116
design/EOS
117117
design/Error
118+
design/FillValues
118119
design/Halo
119120
design/HorzMeshClass
120121
design/Logging

0 commit comments

Comments
 (0)