diff --git a/components/omega/doc/design/StateValidation.md b/components/omega/doc/design/StateValidation.md new file mode 100644 index 000000000000..01bcf131a905 --- /dev/null +++ b/components/omega/doc/design/StateValidation.md @@ -0,0 +1,273 @@ + + +(omega-design-state-validation)= + +# Ocean State Validation + +## 1 Overview + +The Ocean State Validation module provides a mechanism for checking the +physical plausibility of the ocean prognostic state and selected auxiliary +variables at runtime. After each timestep (or at user-defined intervals) the +model can call `validateOceanState` to scan every owned cell and vertical layer +for NaN (Not-a-Number) values and values that lie outside a pre-defined +physically meaningful range. Any detected anomaly is reported through the +OMEGA logging infrastructure as a critical error together with a full stack +backtrace; the run is then terminated via `MPI_Abort` so that corrupted output +is not silently written to disk. + +## 2 Requirements + +### 2.1 Requirement: Check for NaN values + +The validation function must detect NaN values in the ocean state fields. +NaNs indicate a numerical instability or a programming error, and their +presence in the prognostic state makes continued time-stepping meaningless. +Every element of each validated array must be tested. + +### 2.2 Requirement: Check for out-of-bounds values + +In addition to NaN detection, the validation function must check that every +element of each validated field lies within a physically plausible range. +Values outside these ranges indicate catastrophic model failure and should +halt the simulation before corrupted output is written to disk. + +### 2.3 Requirement: Validate LayerThickness + +`LayerThickness` must be validated for each owned cell over all vertical +layers. The valid range is $[10^{-10},\, 1000]$ m. +Negative or near-zero layer thicknesses indicate numerical collapse of the +column and must be caught immediately. + +### 2.4 Requirement: Validate KineticEnergyCell + +`KineticEnergyCell` from the kinematic auxiliary state must be validated for +each owned cell over all vertical layers. The valid range is $[0,\, 10]$ +m$^2$ s$^{-2}$. Negative kinetic energies are unphysical, and values +exceeding 10 m$^2$ s$^{-2}$ correspond to current speeds above +$\sim 4.5$ m s$^{-1}$, which are unrealistic for the open ocean. + +### 2.5 Requirement: Validate Temperature tracer + +Ocean Conservative Temperature must be validated for each owned cell over +all vertical layers. The valid range is $[-10,\, 50]$ °C. +This broad range accommodates all realistic oceanographic regimes including +polar and hydrothermal vent environments. + +### 2.6 Requirement: Validate Salinity tracer + +Ocean Absolute Salinity must be validated for each owned cell over all +vertical layers. The valid range is $[-2,\, 60]$ g kg$^{-1}$. +Values below $-2$ g kg$^{-1}$ are unphysical and values above 60 g kg$^{-1}$ +are outside the valid domain of the TEOS-10 equation of state. + +### 2.7 Requirement: GPU/CPU portability + +All validation kernels must execute on both CPU and GPU hardware using the +Kokkos parallel programming model and therefore must be expressed as Kokkos +parallel reductions. + +### 2.8 Requirement: Informative error reporting + +On detection of any failure the module must log a critical-level message that +identifies the field name, the nature of the problem (NaN or out-of-bounds), +and the number of offending elements. After all fields are checked the +module must additionally print a stack backtrace to assist with debugging, +then abort the run via `MPI_Abort`. + +### 2.9 Requirement: Graceful handling of absent tracers + +If the Temperature or Salinity tracer is not present in the tracer registry +(e.g. in configurations that do not activate active tracers) the +corresponding check must be skipped silently rather than causing an error. + +### 2.10 Desired: Configurable valid ranges + +In the future it may be desirable to allow the user to override the default +valid ranges through the OMEGA configuration system (e.g. for idealised +process studies that intentionally use non-oceanic parameter values). + +### 2.11 Desired: Configurable validation frequency + +In the future it may be desirable to allow the user to control whether +validation is performed every timestep, every N timesteps, or only at +specific points in the run (e.g. after restart reads). + +## 3 Algorithmic Formulation + +No complex numerical algorithms are required. Each field is checked with two +independent `parallelReduce` passes over the domain, restricted to active +cells where `CellMask(i, k) > 0` (inactive cells, such as land cells, are +skipped): + +1. **NaN pass** – counts elements for which `Kokkos::isnan(val)` is `true`. +2. **Bounds pass** – counts elements that are finite yet lie outside + $[\text{MinVal},\, \text{MaxVal}]$. + +Separating the two passes avoids potentially undefined behaviour when +comparing NaN values with `<` or `>`. + +For a 2-D field $f_{i,k}$ with $i \in [0,\, N_\text{cells})$ and +$k \in [0,\, N_\text{vert})$ the two counts are: + +$$ +N_\text{NaN} = \sum_{i,k} M_{i,k}\,\mathbf{1}[\,\text{isnan}(f_{i,k})\,] +$$ + +$$ +N_\text{OOB} = \sum_{i,k} M_{i,k}\,\mathbf{1}[\,\lnot\,\text{isnan}(f_{i,k}) + \land (f_{i,k} < f_\text{min} \lor f_{i,k} > f_\text{max})\,] +$$ + +where $M_{i,k}$ is `CellMask(i, k)` (1 for active cell-layer, 0 for +inactive). + +For a 3-D tracer array $T_{n,i,k}$ the same expressions are applied at a +fixed tracer index $n$. + +The validation traverses only the `NCellsOwned` cells (excluding halo cells) +to avoid double-counting in the parallel decomposition. + +## 4 Design + +The module is implemented as a free function (`validateOceanState`) plus two +file-local helper functions. It does not introduce a class or persistent state. + +### 4.1 Data types and parameters + +#### 4.1.1 Parameters + +The valid ranges for each field are compile-time constants embedded in the +implementation: + +| Field | MinVal | MaxVal | +|----------------------|-----------|--------| +| `LayerThickness` | 1×10⁻¹⁰ | 1000 | +| `KineticEnergyCell` | 0 | 10 | +| `Temperature` | −10 | 50 | +| `Salinity` | −2 | 60 | + +#### 4.1.2 Class/structs/data types + +No new classes or data types are introduced. The module uses the existing +`OceanState`, `AuxiliaryState`, `VertCoord`, and `Tracers` types from the +OMEGA ocean component. + +### 4.2 Methods + +#### 4.2.1 `checkOceanState` (public) + +Performs all field checks and returns the total count of errors found: + +```c++ +I4 checkOceanState(const OceanState *State, + const AuxiliaryState *AuxState, + const VertCoord *VCoord, + I4 TimeLevel); +``` + +Checks all fields described in Section 2, skipping inactive cells +(`CellMask == 0`). Logs critical messages for each type of error. Returns +the total number of errors as an `I4`; returns 0 if all checks pass. Does +**not** abort. Suitable for calling from tests. + +#### 4.2.2 `validateOceanState` (public) + +Production entry-point that aborts on failure: + +```c++ +void validateOceanState(const OceanState *State, + const AuxiliaryState *AuxState, + const VertCoord *VCoord, + I4 TimeLevel); +``` + +Calls `checkOceanState` and aborts via `MPI_Abort` if the return value is +greater than zero. + +#### 4.2.3 `checkArray2D` (file-local helper) + +```c++ +static std::pair checkArray2D(const Array2DReal &Arr, + I4 NRows, I4 NCols, + Real MinVal, Real MaxVal, + bool CheckMin, + const Array2DReal &CellMask); +``` + +Performs the NaN and bounds counts for a 2-D device array over the first +`NRows` rows and `NCols` columns, restricted to active cells +(`CellMask(Row, Col) > 0`). When `CheckMin` is `false` only the upper +bound is enforced (not needed for any current field but useful for future +extension). Returns `{NaNCount, OutOfRangeCount}`. + +#### 4.2.4 `checkTracerArray` (file-local helper) + +```c++ +static std::pair checkTracerArray(const Array3DReal &Tracers3D, + I4 TracerIdx, + I4 NCells, I4 NVert, + Real MinVal, Real MaxVal, + const Array2DReal &CellMask); +``` + +Performs the NaN and bounds counts for a single tracer slice (identified by +`TracerIdx`) of the 3-D tracer array, restricted to active cells +(`CellMask(Cell, K) > 0`). Returns `{NaNCount, OutOfRangeCount}`. + +#### 4.2.5 `abortWithMessage` (file-local helper) + +```c++ +static void abortWithMessage(const std::string &Msg); +``` + +Logs `Msg` at critical severity, prints a stack backtrace using `cpptrace`, +then calls `MPI_Abort` on the default OMEGA communicator with error code +`ErrorCode::Critical`. + +## 5 Verification and Testing + +### 5.1 Test: Valid state passes without abort + +A unit test constructs a minimal OMEGA environment (MachEnv, Decomp, +HorzMesh, VertCoord, Tracers, OceanState, AuxiliaryState) using the standard +test mesh `OmegaMesh.nc`. All state arrays are filled with physically +plausible values: + +- `LayerThickness` = 100 m (valid range [1×10⁻¹⁰, 1000]) +- `NormalVelocity` = 0 m s⁻¹ (not directly validated but required for + `KineticEnergyCell` to be zero) +- `Temperature` = 10 °C (valid range [−10, 50]) +- `Salinity` = 35 g kg⁻¹ (valid range [−2, 60]) + +`AuxiliaryState::computeAll` is called to populate `KineticEnergyCell` +before `validateOceanState` is invoked. + +The test passes if `validateOceanState` returns without calling `MPI_Abort`. + +Tests requirements: 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.9. + +### 5.2 Negative tests: invalid values are detected + +The public `checkOceanState` function is used for negative tests so that +errors can be detected without triggering `MPI_Abort`. Each sub-test: +1. Resets the state to valid values via `restoreValidState`. +2. Injects a single type of invalid value (NaN or OOB) into one field using + a `parallelFor` kernel that overwrites all owned cell-layer entries. +3. Calls `checkOceanState` and verifies a non-zero error count is returned. + +The following sub-tests are implemented: + +| Sub-test | Injected value | Field | +|-------------------------------|------------------------|--------------------| +| `testNaNLayerThickness` | NaN | LayerThickness | +| `testOOBHighLayerThickness` | 2000 m (> max 1000 m) | LayerThickness | +| `testOOBLowLayerThickness` | −1 m (< min 1×10⁻¹⁰) | LayerThickness | +| `testNaNKineticEnergy` | NaN | KineticEnergyCell | +| `testOOBKineticEnergy` | 9999 J kg⁻¹ (> max 10) | KineticEnergyCell | +| `testNaNTemperature` | NaN | Temperature | +| `testOOBTemperature` | 9999 °C (> max 50) | Temperature | +| `testNaNSalinity` | NaN | Salinity | +| `testOOBSalinity` | 9999 g kg⁻¹ (> max 60) | Salinity | + +Tests requirements: 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8. diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 6712bb3e7f8d..8d1b38501df8 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -125,6 +125,7 @@ design/IO design/IOStreams design/Reductions design/State +design/StateValidation design/SubmesoscaleEddies design/Tendency design/Tendencies diff --git a/components/omega/src/ocn/StateValidation.cpp b/components/omega/src/ocn/StateValidation.cpp new file mode 100644 index 000000000000..84e6f382d9f1 --- /dev/null +++ b/components/omega/src/ocn/StateValidation.cpp @@ -0,0 +1,246 @@ +//===-- ocn/StateValidation.cpp - ocean state validation --------*- C++ -*-===// +// +// Validates ocean state fields by checking for NaN values and +// out-of-bounds conditions. Any failure triggers a critical error log with +// backtrace and MPI_Abort on the local communicator. +// +//===----------------------------------------------------------------------===// + +#include "StateValidation.h" + +#include "AuxiliaryState.h" +#include "DataTypes.h" +#include "Error.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "Tracers.h" +#include "VertCoord.h" +#include "mpi.h" + +#include +#include +#include +#include + +namespace OMEGA { + +//------------------------------------------------------------------------------ +// Helper: abort on the local Omega communicator with a message and backtrace +static void abortWithMessage(const std::string &Msg) { + LOG_CRITICAL("{}", Msg); + cpptrace::generate_trace().print(); + MPI_Comm Comm = MachEnv::getDefault()->getComm(); + MPI_Abort(Comm, static_cast(ErrorCode::Critical)); +} + +//------------------------------------------------------------------------------ +// Helper: count NaN entries and out-of-range entries in a 2-D Real device +// array over the first NCells/NEdges rows and NVert columns, restricted to +// active cells (CellMask > 0). +// Returns {NaNCount, OutOfRangeCount}. +static std::pair checkArray2D(const Array2DReal &Arr, I4 NRows, + I4 NCols, Real MinVal, Real MaxVal, + bool CheckMin, + const Array2DReal &CellMask) { + I4 NaNCount = 0; + I4 OutOfRangeCount = 0; + + parallelReduce( + "CheckNaN", {NRows, NCols}, + KOKKOS_LAMBDA(int Row, int Col, int &Accum) { + if (CellMask(Row, Col) == 0) { + return; + } + Real Val = Arr(Row, Col); + if (Kokkos::isnan(Val)) { + ++Accum; + } + }, + NaNCount); + + parallelReduce( + "CheckBounds", {NRows, NCols}, + KOKKOS_LAMBDA(int Row, int Col, int &Accum) { + if (CellMask(Row, Col) == 0) { + return; + } + Real Val = Arr(Row, Col); + if (!Kokkos::isnan(Val)) { + if (Val > MaxVal) { + ++Accum; + } else if (CheckMin && Val < MinVal) { + ++Accum; + } + } + }, + OutOfRangeCount); + + return {NaNCount, OutOfRangeCount}; +} + +//------------------------------------------------------------------------------ +// Helper: count NaN and out-of-range entries for a single tracer (row = cell, +// col = vert) extracted from the 3-D tracer array at the given tracer index, +// restricted to active cells (CellMask > 0). +static std::pair checkTracerArray(const Array3DReal &Tracers3D, + I4 TracerIdx, I4 NCells, I4 NVert, + Real MinVal, Real MaxVal, + const Array2DReal &CellMask) { + I4 NaNCount = 0; + I4 OutOfRangeCount = 0; + + parallelReduce( + "CheckTracerNaN", {NCells, NVert}, + KOKKOS_LAMBDA(int Cell, int K, int &Accum) { + if (CellMask(Cell, K) == 0) { + return; + } + Real Val = Tracers3D(TracerIdx, Cell, K); + if (Kokkos::isnan(Val)) { + ++Accum; + } + }, + NaNCount); + + parallelReduce( + "CheckTracerBounds", {NCells, NVert}, + KOKKOS_LAMBDA(int Cell, int K, int &Accum) { + if (CellMask(Cell, K) == 0) { + return; + } + Real Val = Tracers3D(TracerIdx, Cell, K); + if (!Kokkos::isnan(Val)) { + if (Val < MinVal || Val > MaxVal) { + ++Accum; + } + } + }, + OutOfRangeCount); + + return {NaNCount, OutOfRangeCount}; +} + +//------------------------------------------------------------------------------ +/// Check ocean state fields for NaN and out-of-bounds conditions. +/// Only active cells (where CellMask > 0) are checked. +/// Returns the total count of errors found; does not abort. +I4 checkOceanState(const OceanState *State, const AuxiliaryState *AuxState, + const VertCoord *VCoord, I4 TimeLevel) { + + I4 TotalErrors = 0; + + const Array2DReal &CellMask = VCoord->CellMask; + + // ------------------------------------------------------------------------- + // LayerThickness: valid range [1e-10, 1000] + // ------------------------------------------------------------------------- + { + Array2DReal LayerThick = State->getLayerThickness(TimeLevel); + auto [NaNs, OOB] = + checkArray2D(LayerThick, State->NCellsOwned, State->NVertLayers, + static_cast(1e-10), static_cast(1000.0), + /*CheckMin=*/true, CellMask); + + if (NaNs > 0) { + LOG_CRITICAL( + "StateValidation: LayerThickness contains {} NaN value(s)", NaNs); + TotalErrors += NaNs; + } + if (OOB > 0) { + LOG_CRITICAL("StateValidation: LayerThickness has {} value(s) outside " + "valid range [1e-10, 1000]", + OOB); + TotalErrors += OOB; + } + } + + // ------------------------------------------------------------------------- + // KineticEnergyCell: valid range [0, 10] + // ------------------------------------------------------------------------- + { + const Array2DReal &KE = AuxState->KineticAux.KineticEnergyCell; + auto [NaNs, OOB] = + checkArray2D(KE, State->NCellsOwned, State->NVertLayers, + static_cast(0.0), static_cast(10.0), + /*CheckMin=*/true, CellMask); + + if (NaNs > 0) { + LOG_CRITICAL( + "StateValidation: KineticEnergyCell contains {} NaN value(s)", + NaNs); + TotalErrors += NaNs; + } + if (OOB > 0) { + LOG_CRITICAL( + "StateValidation: KineticEnergyCell has {} value(s) outside " + "valid range [0, 10]", + OOB); + TotalErrors += OOB; + } + } + + // ------------------------------------------------------------------------- + // Temperature tracer: valid range [-10, 50] + // ------------------------------------------------------------------------- + if (Tracers::IndxTemp != -1) { + Array3DReal AllTracers = Tracers::getAll(TimeLevel); + auto [NaNs, OOB] = checkTracerArray( + AllTracers, Tracers::IndxTemp, State->NCellsOwned, State->NVertLayers, + static_cast(-10.0), static_cast(50.0), CellMask); + + if (NaNs > 0) { + LOG_CRITICAL("StateValidation: Temperature contains {} NaN value(s)", + NaNs); + TotalErrors += NaNs; + } + if (OOB > 0) { + LOG_CRITICAL("StateValidation: Temperature has {} value(s) outside " + "valid range [-10, 50]", + OOB); + TotalErrors += OOB; + } + } + + // ------------------------------------------------------------------------- + // Salinity tracer: valid range [-2, 60] + // ------------------------------------------------------------------------- + if (Tracers::IndxSalt != -1) { + Array3DReal AllTracers = Tracers::getAll(TimeLevel); + auto [NaNs, OOB] = checkTracerArray( + AllTracers, Tracers::IndxSalt, State->NCellsOwned, State->NVertLayers, + static_cast(-2.0), static_cast(60.0), CellMask); + + if (NaNs > 0) { + LOG_CRITICAL("StateValidation: Salinity contains {} NaN value(s)", + NaNs); + TotalErrors += NaNs; + } + if (OOB > 0) { + LOG_CRITICAL("StateValidation: Salinity has {} value(s) outside " + "valid range [-2, 60]", + OOB); + TotalErrors += OOB; + } + } + + return TotalErrors; +} + +//------------------------------------------------------------------------------ +/// Validate ocean state fields for NaN and out-of-bounds conditions. +/// Only active cells (where CellMask > 0) are checked. +/// Aborts via MPI_Abort on failure. +void validateOceanState(const OceanState *State, const AuxiliaryState *AuxState, + const VertCoord *VCoord, I4 TimeLevel) { + + if (checkOceanState(State, AuxState, VCoord, TimeLevel) > 0) { + abortWithMessage("StateValidation: Ocean state validation failed. " + "See critical messages above for details."); + } +} + +} // namespace OMEGA + +//===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/StateValidation.h b/components/omega/src/ocn/StateValidation.h new file mode 100644 index 000000000000..bf49ba74f63d --- /dev/null +++ b/components/omega/src/ocn/StateValidation.h @@ -0,0 +1,70 @@ +#ifndef OMEGA_STATEVALIDATION_H +#define OMEGA_STATEVALIDATION_H +//===-- ocn/StateValidation.h - ocean state validation ----------*- C++ -*-===// +// +/// \file +/// \brief Declares state validation functions for ocean state validation +/// +/// Provides two functions: +/// - checkOceanState: checks for NaN and out-of-bounds conditions and +/// returns the total error count without aborting. Suitable for testing. +/// - validateOceanState: calls checkOceanState and aborts via MPI_Abort if +/// any errors are found. +// +//===----------------------------------------------------------------------===// + +#include "AuxiliaryState.h" +#include "OceanState.h" +#include "Tracers.h" +#include "VertCoord.h" + +namespace OMEGA { + +/// Check ocean state fields for NaN values and out-of-bounds conditions +/// without aborting, returning the total count of errors found. +/// +/// Only active ocean cells (where CellMask > 0) are checked. Critical log +/// messages are emitted for each type of error found. +/// +/// The following fields are validated: +/// - LayerThickness : [1e-10, 1000] (from OceanState) +/// - KineticEnergyCell : [0, 10] +/// (from AuxiliaryState::KineticAux) +/// - Temperature tracer : [-10, 50] (from Tracers) +/// - Salinity tracer : [-2, 60] (from Tracers) +/// +/// \param[in] State Ocean state to validate +/// \param[in] AuxState Auxiliary state containing KineticEnergyCell +/// \param[in] VCoord Vertical coordinate containing the CellMask +/// \param[in] TimeLevel Time level index to validate (typically 0 = current) +/// \return I4 total count of errors found across all checked fields (0 = valid) +I4 checkOceanState(const OceanState *State, const AuxiliaryState *AuxState, + const VertCoord *VCoord, I4 TimeLevel); + +/// Check ocean state fields for NaN values and out-of-bounds conditions. +/// +/// Calls checkOceanState and aborts via MPI_Abort if any errors are found. +/// Only active ocean cells (where CellMask > 0) are checked. +/// +/// The following fields are validated: +/// - LayerThickness : [1e-10, 1000] (from OceanState) +/// - KineticEnergyCell : [0, 10] +/// (from AuxiliaryState::KineticAux) +/// - Temperature tracer : [-10, 50] (from Tracers) +/// - Salinity tracer : [-2, 60] (from Tracers) +/// +/// If any check fails a critical error is logged with an informative message +/// and a stack backtrace, and the run is aborted via MPI_Abort on the +/// communicator obtained from the default MachEnv. +/// +/// \param[in] State Ocean state to validate +/// \param[in] AuxState Auxiliary state containing KineticEnergyCell +/// \param[in] VCoord Vertical coordinate containing the CellMask +/// \param[in] TimeLevel Time level index to validate (typically 0 = current) +void validateOceanState(const OceanState *State, const AuxiliaryState *AuxState, + const VertCoord *VCoord, I4 TimeLevel); + +} // namespace OMEGA + +//===----------------------------------------------------------------------===// +#endif // defined OMEGA_STATEVALIDATION_H diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 301a04ede49c..4d3830103691 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -516,3 +516,14 @@ add_omega_test( ocn/VertAdvTest.cpp "-n;1" ) + +########################## +# State Validation test +########################## + +add_omega_test( + STATE_VALIDATION_TEST + testStateValidation.exe + ocn/StateValidationTest.cpp + "-n;8" +) diff --git a/components/omega/test/ocn/StateValidationTest.cpp b/components/omega/test/ocn/StateValidationTest.cpp new file mode 100644 index 000000000000..cbff1cc956c0 --- /dev/null +++ b/components/omega/test/ocn/StateValidationTest.cpp @@ -0,0 +1,413 @@ +//===-- Test driver for OMEGA StateValidation -----------------------*- C++ +//-*-===/ +// +/// \file +/// \brief Test driver for ocean state validation +/// +/// Tests both the positive path (valid state passes) and the negative paths +/// (NaN and out-of-bounds values in each checked field are detected) of +/// validateOceanState / checkOceanState. Checked fields are: +/// - LayerThickness, KineticEnergyCell, Temperature, Salinity. +// +//===-----------------------------------------------------------------------===/ + +#include "StateValidation.h" +#include "AuxiliaryState.h" +#include "Config.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 "Pacer.h" +#include "TimeStepper.h" +#include "Tracers.h" +#include "VertAdv.h" +#include "VertCoord.h" +#include "mpi.h" + +#include +#include + +using namespace OMEGA; + +//------------------------------------------------------------------------------ +// Initialize the Omega subsystems required for state validation testing + +int initStateValidationTest(const std::string &MeshFile) { + int Err = 0; + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + initLogging(DefEnv); + LOG_INFO("------ StateValidation unit tests ------"); + + Config("Omega"); + Config::readAll("omega.yml"); + + TimeStepper::init1(); + + IO::init(DefComm); + Decomp::init(MeshFile); + + IOStream::init(); + + int HaloErr = Halo::init(); + if (HaloErr != 0) { + Err++; + LOG_ERROR("StateValidationTest: error initializing default halo"); + } + + HorzMesh::init(); + VertCoord::init(); + Tracers::init(); + + int StateErr = OceanState::init(); + if (StateErr != 0) { + Err++; + LOG_ERROR("StateValidationTest: error initializing default state"); + } + + VertAdv::init(); + + return Err; +} + +//------------------------------------------------------------------------------ +// Fill state and auxiliary/tracer arrays with known-valid values + +static int fillValidState() { + int Err = 0; + + auto *State = OceanState::getDefault(); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + const int NEdges = State->NEdgesAll; + + // LayerThickness: fill with 100 m (valid range [1e-10, 1000]) + Array2DReal LayerThick = State->getLayerThickness(0); + parallelFor( + "FillLayerThick", {NCells, NVert}, + KOKKOS_LAMBDA(int ICell, int K) { LayerThick(ICell, K) = 100.0; }); + + // NormalVelocity: fill with 0 (not checked, but needed for AuxState) + Array2DReal NormalVel = State->getNormalVelocity(0); + parallelFor( + "FillNormalVel", {NEdges, NVert}, + KOKKOS_LAMBDA(int IEdge, int K) { NormalVel(IEdge, K) = 0.0; }); + + // Exchange halos so auxiliary state computations are consistent + State->exchangeHalo(0); + + // Tracers: fill Temperature with 10 C and Salinity with 35 g/kg + // Use deepCopy with individual tracer subviews + if (Tracers::getNumTracers() > 0) { + // Temperature = 10.0 (valid: -10 to 50) + if (Tracers::IndxTemp != -1) { + Array2DReal TempArr = Tracers::getByIndex(0, Tracers::IndxTemp); + deepCopy(TempArr, static_cast(10.0)); + } + + // Salinity = 35.0 (valid: -2 to 60) + if (Tracers::IndxSalt != -1) { + Array2DReal SaltArr = Tracers::getByIndex(0, Tracers::IndxSalt); + deepCopy(SaltArr, static_cast(35.0)); + } + } + + return Err; +} + +//------------------------------------------------------------------------------ +// Restore the default state to valid values and recompute the auxiliary state + +static void restoreValidState(OceanState *State, AuxiliaryState *AuxState) { + fillValidState(); + Array3DReal AllTracers = Tracers::getAll(0); + AuxState->computeAll(State, AllTracers, 0); +} + +//------------------------------------------------------------------------------ +// Positive test: validate a clean, valid state — expects 0 errors + +static int testValidState(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + LOG_INFO("StateValidationTest: Testing validation on valid state"); + validateOceanState(State, AuxState, VCoord, 0); + LOG_INFO("StateValidationTest: Valid state validation PASS"); + return 0; +} + +//------------------------------------------------------------------------------ +// Negative tests: inject an invalid value, verify checkOceanState returns > 0, +// then restore valid state. +// +// Returns 0 on success (i.e. the error was caught), 1 otherwise. + +static int expectErrors(const char *TestName, OceanState *State, + AuxiliaryState *AuxState, VertCoord *VCoord) { + I4 Errs = checkOceanState(State, AuxState, VCoord, 0); + if (Errs == 0) { + LOG_ERROR("StateValidationTest: {} - expected errors but got none", + TestName); + return 1; + } + LOG_INFO("StateValidationTest: {} PASS (caught {} error(s))", TestName, + Errs); + return 0; +} + +// --- LayerThickness --- + +static int testNaNLayerThickness(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + restoreValidState(State, AuxState); + const Real NaN = std::numeric_limits::quiet_NaN(); + Array2DReal LT = State->getLayerThickness(0); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + parallelFor( + "InjectNaNLayerThick", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { LT(I, K) = NaN; }); + return expectErrors("NaN in LayerThickness", State, AuxState, VCoord); +} + +static int testOOBHighLayerThickness(OceanState *State, + AuxiliaryState *AuxState, + VertCoord *VCoord) { + restoreValidState(State, AuxState); + Array2DReal LT = State->getLayerThickness(0); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + // 2000 m is above the valid max of 1000 m + parallelFor( + "InjectOOBHighLayerThick", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { LT(I, K) = 2000.0; }); + return expectErrors("OOB-high in LayerThickness", State, AuxState, VCoord); +} + +static int testOOBLowLayerThickness(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + restoreValidState(State, AuxState); + Array2DReal LT = State->getLayerThickness(0); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + // -1.0 m is below the valid min of 1e-10 m + parallelFor( + "InjectOOBLowLayerThick", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { LT(I, K) = -1.0; }); + return expectErrors("OOB-low in LayerThickness", State, AuxState, VCoord); +} + +// --- KineticEnergyCell --- + +static int testNaNKineticEnergy(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + restoreValidState(State, AuxState); + const Real NaN = std::numeric_limits::quiet_NaN(); + Array2DReal KE = AuxState->KineticAux.KineticEnergyCell; + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + parallelFor( + "InjectNaNKE", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { KE(I, K) = NaN; }); + return expectErrors("NaN in KineticEnergyCell", State, AuxState, VCoord); +} + +static int testOOBKineticEnergy(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + restoreValidState(State, AuxState); + Array2DReal KE = AuxState->KineticAux.KineticEnergyCell; + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + // 9999 J/kg is above the valid max of 10 J/kg + parallelFor( + "InjectOOBKE", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { KE(I, K) = 9999.0; }); + return expectErrors("OOB in KineticEnergyCell", State, AuxState, VCoord); +} + +// --- Temperature tracer --- + +static int testNaNTemperature(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + if (Tracers::IndxTemp == -1) + return 0; // tracer not active; skip + restoreValidState(State, AuxState); + const Real NaN = std::numeric_limits::quiet_NaN(); + Array2DReal TempArr = Tracers::getByIndex(0, Tracers::IndxTemp); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + parallelFor( + "InjectNaNTemp", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { TempArr(I, K) = NaN; }); + return expectErrors("NaN in Temperature", State, AuxState, VCoord); +} + +static int testOOBTemperature(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + if (Tracers::IndxTemp == -1) + return 0; // tracer not active; skip + restoreValidState(State, AuxState); + Array2DReal TempArr = Tracers::getByIndex(0, Tracers::IndxTemp); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + // 9999 C is above the valid max of 50 C + parallelFor( + "InjectOOBTemp", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { TempArr(I, K) = 9999.0; }); + return expectErrors("OOB in Temperature", State, AuxState, VCoord); +} + +// --- Salinity tracer --- + +static int testNaNSalinity(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + if (Tracers::IndxSalt == -1) + return 0; // tracer not active; skip + restoreValidState(State, AuxState); + const Real NaN = std::numeric_limits::quiet_NaN(); + Array2DReal SaltArr = Tracers::getByIndex(0, Tracers::IndxSalt); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + parallelFor( + "InjectNaNSalt", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { SaltArr(I, K) = NaN; }); + return expectErrors("NaN in Salinity", State, AuxState, VCoord); +} + +static int testOOBSalinity(OceanState *State, AuxiliaryState *AuxState, + VertCoord *VCoord) { + if (Tracers::IndxSalt == -1) + return 0; // tracer not active; skip + restoreValidState(State, AuxState); + Array2DReal SaltArr = Tracers::getByIndex(0, Tracers::IndxSalt); + const int NCells = State->NCellsAll; + const int NVert = State->NVertLayers; + // 9999 g/kg is above the valid max of 60 g/kg + parallelFor( + "InjectOOBSalt", {NCells, NVert}, + KOKKOS_LAMBDA(int I, int K) { SaltArr(I, K) = 9999.0; }); + return expectErrors("OOB in Salinity", State, AuxState, VCoord); +} + +//------------------------------------------------------------------------------ +// Run all state validation tests + +int testStateValidation() { + int Err = 0; + + // Initialize the auxiliary state (needed for KineticEnergyCell) + AuxiliaryState::init(); + auto *DefAuxState = AuxiliaryState::getDefault(); + + if (!DefAuxState) { + Err++; + LOG_ERROR("StateValidationTest: Default AuxiliaryState not found"); + return Err; + } + + auto *DefState = OceanState::getDefault(); + if (!DefState) { + Err++; + LOG_ERROR("StateValidationTest: Default OceanState not found"); + return Err; + } + + VertCoord *VCoord = VertCoord::getDefault(); + + // Fill state arrays with valid values and compute auxiliary state + fillValidState(); + { + Array3DReal AllTracers = Tracers::getAll(0); + DefAuxState->computeAll(DefState, AllTracers, 0); + } + + // ---- Positive test ---- + Err += testValidState(DefState, DefAuxState, VCoord); + + // ---- Negative tests: each injects a bad value and verifies detection ---- + + // LayerThickness + Err += testNaNLayerThickness(DefState, DefAuxState, VCoord); + Err += testOOBHighLayerThickness(DefState, DefAuxState, VCoord); + Err += testOOBLowLayerThickness(DefState, DefAuxState, VCoord); + + // KineticEnergyCell + Err += testNaNKineticEnergy(DefState, DefAuxState, VCoord); + Err += testOOBKineticEnergy(DefState, DefAuxState, VCoord); + + // Temperature tracer + Err += testNaNTemperature(DefState, DefAuxState, VCoord); + Err += testOOBTemperature(DefState, DefAuxState, VCoord); + + // Salinity tracer + Err += testNaNSalinity(DefState, DefAuxState, VCoord); + Err += testOOBSalinity(DefState, DefAuxState, VCoord); + + AuxiliaryState::clear(); + return Err; +} + +//------------------------------------------------------------------------------ +// Finalize Omega objects + +void finalizeStateValidationTest() { + Tracers::clear(); + OceanState::clear(); + VertAdv::clear(); + VertCoord::clear(); + HorzMesh::clear(); + Field::clear(); + Dimension::clear(); + TimeStepper::clear(); + Halo::clear(); + Decomp::clear(); + MachEnv::removeAll(); +} + +//------------------------------------------------------------------------------ +// Main entry point + +int main(int argc, char *argv[]) { + int RetVal = 0; + + MPI_Init(&argc, &argv); + Kokkos::initialize(argc, argv); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + + { + int Err = initStateValidationTest("OmegaMesh.nc"); + if (Err != 0) { + LOG_CRITICAL("StateValidationTest: Error during initialization"); + } else { + RetVal += testStateValidation(); + } + finalizeStateValidationTest(); + } + + if (RetVal == 0) + LOG_INFO("------ StateValidation unit tests successful ------"); + + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + if (RetVal >= 256) + RetVal = 255; + + return RetVal; + +} // end of main +//===-----------------------------------------------------------------------===/