diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index d46de8c24e95..d71f8538d4b8 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -34,6 +34,8 @@ Omega: VertCoord: Density0: 1026.0 MovementWeightType: Uniform + PressureGrad: + PressureGradType: Centered Tendencies: ThicknessFluxTendencyEnable: true PVTendencyEnable: true @@ -57,6 +59,7 @@ Omega: ThicknessVertAdvTendencyEnable: true VelocityVertAdvTendencyEnable: true TracerVertAdvTendencyEnable: true + PressureGradTendencyEnable: true Tracers: Base: [Temperature, Salinity] Debug: [Debug1, Debug2, Debug3] diff --git a/components/omega/doc/devGuide/PGrad.md b/components/omega/doc/devGuide/PGrad.md new file mode 100644 index 000000000000..cfbe45fcabd8 --- /dev/null +++ b/components/omega/doc/devGuide/PGrad.md @@ -0,0 +1,185 @@ +(omega-dev-pgrad)= + +# Pressure Gradient (PGrad) + +Omega includes a `PressureGrad` class that computes horizontal pressure gradient +tendencies for the non-Boussinesq momentum equation. The implementation supports a +centered difference scheme as the default, with a placeholder for future high-order +methods. The class follows the same factory pattern used by other Omega modules. + +## PressureGradType enum + +An enumeration of the available pressure gradient schemes is defined in `PGrad.h`: + +```c++ +enum class PressureGradType { Centered, HighOrder1, HighOrder2 }; +``` + +This is used to select which pressure gradient method is applied at runtime. + +## Initialization + +An instance of `PressureGrad` requires both a [`HorzMesh`](#omega-dev-horz-mesh) and +a [`VertCoord`](#omega-dev-vert-coord), so these classes and all of their dependencies +must be initialized before `PressureGrad` can be initialized. The static method: + +```c++ +OMEGA::PressureGrad::init(); +``` + +initializes the default `PressureGrad` instance using the default `HorzMesh` and +`VertCoord` instances and the global Omega configuration. A pointer to the default +instance can be retrieved at any time using: + +```c++ +OMEGA::PressureGrad* DefPGrad = OMEGA::PressureGrad::getDefault(); +``` + +## Creating additional instances + +Additional named instances can be created using: + +```c++ +OMEGA::PressureGrad* MyPGrad = + OMEGA::PressureGrad::create("MyPGrad", Mesh, VCoord, Options); +``` + +where `Mesh` is a pointer to a `HorzMesh`, `VCoord` is a pointer to a `VertCoord`, +and `Options` is a pointer to the `Config`. A named instance can be retrieved later +using: + +```c++ +OMEGA::PressureGrad* MyPGrad = OMEGA::PressureGrad::get("MyPGrad"); +``` + +## Constructor behavior + +The constructor reads the `PressureGrad` section from the configuration, stores +references to mesh and vertical coordinate data needed for computation, and enables +the appropriate functor based on the configured `PressureGradType`. It also allocates +placeholder arrays for tidal potential and self-attraction/loading, which are intended +to be populated by a future tidal forcing module. Currently these arrays are initialized +to zero. + +## Computing the pressure gradient + +To compute pressure gradient tendencies and accumulate them into a tendency array: + +```c++ +PGrad->computePressureGrad(Tend, State, VCoord, EqState, TimeLevel); +``` + +where: +- `Tend` is a 2D array `(NEdgesAll × NVertLayers)` that the pressure gradient + tendency is accumulated into +- `State` is the current `OceanState`, from which layer thickness is extracted + at the given `TimeLevel` +- `VCoord` provides pressure, interface height, and geopotential fields +- `EqState` provides the specific volume field +- `TimeLevel` selects which time level of the state to use + +The method uses hierarchical Kokkos parallelism: an outer `parallelForOuter` loop +iterates over edges, and an inner `parallelForInner` loop iterates over vertical +chunks. The appropriate functor is dispatched based on `PressureGradChoice`. + +## Functors + +### PressureGradCentered + +This functor implements a centered difference approximation of the pressure gradient +tendency. For each edge, it first computes the layer-invariant tidal and +self-attraction/loading contribution: + +``` +GradGeoPot = grad(TidalPotential) + grad(SelfAttractionLoading) +``` + +Then, for each vertical layer `K`, it computes three terms: + +1. **Montgomery potential gradient**: The average of the horizontal gradients of the + Montgomery potential ($\alpha p + g z$) at the top (interface `K`) and bottom + (interface `K+1`) of the layer. This compactly represents the combined effect + of the pressure gradient and the geopotential contribution from tilted coordinate + surfaces. + +2. **Specific volume correction**: A correction term equal to the edge-averaged + pressure at mid-layer multiplied by the horizontal gradient of specific volume. + This accounts for horizontal density variations that are not captured by the + Montgomery potential form. + +3. **Tidal and geopotential forcing** (`GradGeoPot`): The external geopotential + contribution from tidal forcing and self-attraction/loading, applied uniformly + across all layers at an edge. + +The tendency update for each layer is: + +``` +Tend(IEdge, K) += EdgeMask(IEdge, K) * (-GradMontPot + PGradAlpha - GradGeoPot) +``` + +where `EdgeMask` is applied to enforce land boundary conditions. The functor operator +signature is: + +```c++ +KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, + const Array2DReal &PressureMid, + const Array2DReal &PressureInterface, + const Array2DReal &ZInterface, + const Array1DReal &TidalPotential, + const Array1DReal &SelfAttractionLoading, + const Array2DReal &SpecVol) const; +``` + +### PressureGradHighOrder + +This functor is a placeholder for a future high-order pressure gradient implementation +suitable for ice shelf cavities and complex bathymetry. Currently it performs no +computation (a no-op). + +## Configuration + +The pressure gradient type is selected in the input YAML file: + +```yaml +PressureGrad: + PressureGradType: 'centered' +``` + +Valid options for `PressureGradType` are: +- `'centered'` or `'Centered'`: centered difference approximation (default) +- `'HighOrder1'`: first high-order method (placeholder, future implementation) + +If an unrecognized value is provided, the implementation falls back to the centered +scheme and logs an informational message. + +## Data members + +The `PressureGrad` class stores the following key data: + +| Member | Type | Description | +| ------ | ---- | ----------- | +| `NEdgesAll` | `I4` | Total number of edges including halo | +| `NChunks` | `I4` | Number of vertical chunks for vectorization | +| `NVertLayers` | `I4` | Number of vertical layers | +| `NVertLayersP1` | `I4` | Number of vertical layers plus one | +| `MinLayerEdgeBot` | `Array1DI4` | Minimum active layer index for each edge | +| `MaxLayerEdgeTop` | `Array1DI4` | Maximum active layer index for each edge | +| `TidalPotential` | `Array1DReal` | Tidal potential (placeholder, currently zero) | +| `SelfAttractionLoading` | `Array1DReal` | Self-attraction and loading term (placeholder, currently zero) | +| `CenteredPGrad` | `PressureGradCentered` | Centered pressure gradient functor | +| `HighOrderPGrad` | `PressureGradHighOrder` | High-order pressure gradient functor | +| `PressureGradChoice` | `PressureGradType` | Selected pressure gradient method | + +## Removal + +To remove all `PressureGrad` instances: + +```c++ +OMEGA::PressureGrad::clear(); +``` + +To remove a specific named instance: + +```c++ +OMEGA::PressureGrad::erase("MyPGrad"); +``` diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 1dd11e400a43..41bd611140b3 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -49,6 +49,7 @@ userGuide/Reductions userGuide/Tracers userGuide/TridiagonalSolvers userGuide/VertCoord +userGuide/PGrad userGuide/Timing userGuide/VerticalMixingCoeff userGuide/VertAdv @@ -93,6 +94,7 @@ devGuide/Reductions devGuide/Tracers devGuide/TridiagonalSolvers devGuide/VertCoord +devGuide/PGrad devGuide/Timing devGuide/VerticalMixingCoeff devGuide/VertAdv diff --git a/components/omega/doc/userGuide/PGrad.md b/components/omega/doc/userGuide/PGrad.md new file mode 100644 index 000000000000..0fc7bb5968a2 --- /dev/null +++ b/components/omega/doc/userGuide/PGrad.md @@ -0,0 +1,68 @@ +(omega-user-pgrad)= + +# Pressure Gradient + +The pressure gradient term in the momentum equation represents the force per unit +mass due to horizontal variations in pressure and geopotential. This term is +essential for capturing the dynamics of +ocean circulation, including both barotropic and baroclinic motions. + +## Physical Background + +In the layered non-Boussinesq momentum equation solved in Omega, the pressure +gradient tendency for each edge and layer includes three contributions: + +1. **Montgomery potential gradient**: The horizontal gradient of the Montgomery + potential ($\alpha p + g z$), averaged across the top and bottom interfaces of + each layer. The Montgomery potential combines the pressure gradient and the + geopotential, and its gradient along coordinate surfaces accounts for both the + direct pressure force and the effect of tilted layer interfaces that arise when + using a general vertical coordinate. + +2. **Specific volume correction**: A correction term proportional to the gradient + of specific volume (inverse density) at each edge. This term ensures that + horizontal density variations between the two cells sharing an edge are properly + represented in the pressure gradient force. + +3. **External geopotential forcing**: Contributions from the tidal potential and + the self-attraction and loading (SAL) terms. These represent gravitational + forcing from astronomical tides and the deformation of the solid Earth and ocean + surface in response to the ocean mass distribution. These terms are currently set + to zero and will be provided by a future tidal forcing module. + +## Configuration Options + +The pressure gradient method is configured in the input YAML file under the +`PressureGrad` section: + +```yaml +PressureGrad: + PressureGradType: 'centered' +``` + +### Available Methods + +**Centered Difference** (`'centered'` or `'Centered'`) +- Computes the pressure gradient using a centered finite-difference approximation + of the Montgomery potential gradient and specific volume correction +- Suitable for global ocean simulations without ice shelf cavities +- Default and currently the only fully implemented option + +**High-Order** (`'HighOrder1'`) +- Placeholder for a future high-order pressure gradient method based on volume + integral formulations +- Intended for simulations with ice shelf cavities and steep bathymetry where the + centered scheme may be inaccurate +- Not yet implemented; selecting this option produces zero pressure gradient tendency + +## Dependencies + +The pressure gradient calculation requires the following Omega components to be +initialized first: + +- [**Horizontal Mesh**](omega-user-horz-mesh): provides mesh geometry including + distances between cell centers and edge connectivity +- [**Vertical Coordinate**](omega-user-vert-coord): provides pressure at layer + mid-points and interfaces, interface heights ($z$), and geopotential +- [**Equation of State**](omega-user-eos): provides the specific volume field +- [**Ocean State**](omega-user-state): provides the current layer thicknesses diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 707297f01113..06c9695f7267 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -44,7 +44,6 @@ class Teos10Eos { I4 KDisp) const { Real SpecVolPCoeffs[6 * VecLength]; - const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); @@ -63,15 +62,15 @@ class Teos10Eos { if (KDisp == 0) { // No displacement SpecVol(ICell, K) = - calcRefProfile(Pressure(ICell, K)) + - calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, K)); + calcRefProfile(Pressure(ICell, K) * Pa2Db) + + calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, K) * Pa2Db); } else { // Displacement, use the displaced pressure I4 KTmp = Kokkos::min(K + KDisp, MaxLayerCell(ICell)); KTmp = Kokkos::max(MinLayerCell(ICell), KTmp); SpecVol(ICell, K) = - calcRefProfile(Pressure(ICell, KTmp)) + - calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, KTmp)); + calcRefProfile(Pressure(ICell, KTmp) * Pa2Db) + + calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, KTmp) * Pa2Db); } } } @@ -312,15 +311,15 @@ class Teos10BruntVaisalaFreqSq { Real PInt = 0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1)); Real SpInt = 0.5_Real * (SpecVol(ICell, K) + SpecVol(ICell, K - 1)); - Real AlphaInt = calcAlpha(SaInt, CtInt, PInt, SpInt); - Real BetaInt = calcBeta(SaInt, CtInt, PInt, SpInt); + Real AlphaInt = calcAlpha(SaInt, CtInt, PInt * Pa2Db, SpInt); + Real BetaInt = calcBeta(SaInt, CtInt, PInt * Pa2Db, SpInt); Real DSa = AbsSalinity(ICell, K) - AbsSalinity(ICell, K - 1); Real DCt = ConservTemp(ICell, K) - ConservTemp(ICell, K - 1); Real DP = Pressure(ICell, K) - Pressure(ICell, K - 1); BruntVaisalaFreqSq(ICell, K) = Gravity * Gravity * (BetaInt * DSa - AlphaInt * DCt) / - (SpInt * Db2Pa * DP); + (SpInt * DP); } } } @@ -568,7 +567,6 @@ class Eos { const Array2DReal &AbsSalinity, const Array2DReal &Pressure, const Array2DReal &SpecVol); - /// Initialize EOS from config and mesh static void init(); diff --git a/components/omega/src/ocn/OceanFinal.cpp b/components/omega/src/ocn/OceanFinal.cpp index 3715e95baef6..a518db52a175 100644 --- a/components/omega/src/ocn/OceanFinal.cpp +++ b/components/omega/src/ocn/OceanFinal.cpp @@ -7,6 +7,7 @@ #include "AuxiliaryState.h" #include "Decomp.h" +#include "Eos.h" #include "Field.h" #include "Halo.h" #include "HorzMesh.h" @@ -14,6 +15,7 @@ #include "MachEnv.h" #include "OceanDriver.h" #include "OceanState.h" +#include "PGrad.h" #include "Tendencies.h" #include "TimeMgr.h" #include "TimeStepper.h" @@ -33,6 +35,8 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time // clean up all objects Tracers::clear(); TimeStepper::clear(); + PressureGrad::clear(); + Eos::destroyInstance(); Tendencies::clear(); AuxiliaryState::clear(); OceanState::clear(); diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index b5c21aa1508f..a22e3b3b0667 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -11,6 +11,7 @@ #include "Config.h" #include "DataTypes.h" #include "Decomp.h" +#include "Eos.h" #include "Error.h" #include "Field.h" #include "Halo.h" @@ -21,6 +22,7 @@ #include "MachEnv.h" #include "OceanDriver.h" #include "OceanState.h" +#include "PGrad.h" #include "Pacer.h" #include "Tendencies.h" #include "TimeMgr.h" @@ -134,6 +136,8 @@ int initOmegaModules(MPI_Comm Comm) { Tracers::init(); VertAdv::init(); AuxiliaryState::init(); + PressureGrad::init(); + Eos::init(); Tendencies::init(); TimeStepper::init2(); diff --git a/components/omega/src/ocn/OceanState.cpp b/components/omega/src/ocn/OceanState.cpp index aef68a66ef51..4d3e7ce8919a 100644 --- a/components/omega/src/ocn/OceanState.cpp +++ b/components/omega/src/ocn/OceanState.cpp @@ -44,6 +44,9 @@ int OceanState::init() { LOG_ERROR("TimeStepper needs to be initialized before OceanState"); } int NTimeLevels = DefTimeStepper->getNTimeLevels(); + LOG_INFO("OceanState: Initializing default state with {} vertical layers " + "and {} time levels", + NVertLayers, NTimeLevels); if (NTimeLevels < 2) { LOG_ERROR("OceanState: the number of time level is lower than 2"); diff --git a/components/omega/src/ocn/PGrad.cpp b/components/omega/src/ocn/PGrad.cpp new file mode 100644 index 000000000000..33f3320c7606 --- /dev/null +++ b/components/omega/src/ocn/PGrad.cpp @@ -0,0 +1,237 @@ +//===-- ocn/PGrad.cpp - Pressure Gradient Term -----------------*- C++ -*-===// +// +// Implements the PGrad manager and two discretizations: Centered and +// HighOrder. +// +//===----------------------------------------------------------------------===// + +#include "PGrad.h" +#include "Eos.h" +#include "Error.h" +#include "Field.h" +#include "HorzMesh.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" + +namespace OMEGA { + +PressureGrad *PressureGrad::DefaultPGrad = nullptr; +std::map> PressureGrad::AllPGrads; + +//------------------------------------------------------------------------------ +// Initialize the PressureGrad. Assumes that HorzMesh and VertCoord have already +// been initialized. +void PressureGrad::init() { + + // Retrieve default mesh and vertical coordinate + HorzMesh *DefMesh = HorzMesh::getDefault(); + VertCoord *DefVCoord = VertCoord::getDefault(); + + // Retrieve omega config + Config *OmegaConfig = Config::getOmegaConfig(); + + // Create the default PressureGrad and set pointer to it + PressureGrad::DefaultPGrad = + PressureGrad::create("Default", DefMesh, DefVCoord, OmegaConfig); + +} // end init + +//------------------------------------------------------------------------------ +// Create a new PressureGrad object and add it to the map +PressureGrad * +PressureGrad::create(const std::string &Name, /// [in] Name for PressureGrad + const HorzMesh *Mesh, ///< [in] Horizontal mesh + const VertCoord *VCoord, ///< [in] Vertical coordinate + Config *Options) { ///< [in] Configuration options + + // Check to see if a PressureGrad of the same name already exists and + // if so, exit with an error + if (AllPGrads.find(Name) != AllPGrads.end()) { + LOG_ERROR("Attempted to create a PressureGrad with name {} but a " + "PressureGrad of that name already exists", + Name); + return nullptr; + } + + // create a new PressureGrad on the heap and put it in a map of + // unique_ptrs, which will manage its lifetime + auto *NewPGrad = new PressureGrad(Mesh, VCoord, Options); + AllPGrads.emplace(Name, NewPGrad); + + return NewPGrad; + +} // end create + +//------------------------------------------------------------------------------ +// Get the default pressure gradient instance +PressureGrad *PressureGrad::getDefault() { + + return DefaultPGrad; + +} // end get default + +//------------------------------------------------------------------------------ +// Constructor for PressureGrad +PressureGrad::PressureGrad( + const HorzMesh *Mesh, ///< [in] Horizontal mesh + const VertCoord *VCoord, ///< [in] Vertical coordinate + Config *Options) ///< [in] Configuration options + : MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), CenteredPGrad(Mesh, VCoord), + HighOrderPGrad(Mesh, VCoord) { + + // store mesh sizes + NEdgesAll = Mesh->NEdgesAll; + NEdgesOwned = Mesh->NEdgesOwned; + NVertLayers = VCoord->NVertLayers; + NVertLayersP1 = NVertLayers + 1; + + // Read config options for PressureGrad type + // and enable the appropriate functor + Config PGradConfig("PressureGrad"); + Error Err; + Err += Options->get(PGradConfig); + CHECK_ERROR_ABORT(Err, + "PressureGrad: PressureGrad group not found in Config"); + std::string PGradTypeStr; + Err += PGradConfig.get("PressureGradType", PGradTypeStr); + + if (PGradTypeStr == "centered" || PGradTypeStr == "Centered") { + PressureGradChoice = PressureGradType::Centered; + this->CenteredPGrad.Enabled = true; + } else if (PGradTypeStr == "HighOrder1") { + PressureGradChoice = PressureGradType::HighOrder1; + this->HighOrderPGrad.Enabled = true; + } else { + LOG_INFO( + "PGrad: Unknown PressureGradType in config, defaulting to centered"); + } + + // Temporary: initialization of tidal potential and SAL + TidalPotential = Array1DReal("TidalPotential", Mesh->NCellsSize); + SelfAttractionLoading = + Array1DReal("SelfAttractionLoading", Mesh->NCellsSize); + +} // end constructor + +//------------------------------------------------------------------------------ +// Destructor for PressureGrad +PressureGrad::~PressureGrad() { + + // No operations needed, Kokkos arrays removed when no longer in scope + +} // end destructor + +//------------------------------------------------------------------------------ +// Remove PressureGrad instances before exit +void PressureGrad::clear() { + + AllPGrads.clear(); + + DefaultPGrad = nullptr; // prevent dangling pointe'r + +} // end clear + +//------------------------------------------------------------------------------ +// Remove PressureGrad from list by name +void PressureGrad::erase(const std::string &Name) { + + AllPGrads.erase(Name); + +} // end erase + +//------------------------------------------------------------------------------ +// Get pressure gradient instance by name +PressureGrad *PressureGrad::get(const std::string &Name ///< [in] Name of +) { + + auto it = AllPGrads.find(Name); + + if (it != AllPGrads.end()) { + return it->second.get(); + } else { + LOG_ERROR("PressureGrad::get: Attempt to retrieve non-existent " + "PressureGrad:"); + LOG_ERROR("{} has not been defined or has been removed", Name); + return nullptr; + } + +} // end get pressure gradient + +//------------------------------------------------------------------------------ +// Compute pressure gradient tendencies and add into Tend array +void PressureGrad::computePressureGrad(Array2DReal &Tend, + const Array2DReal &PressureMid, + const Array2DReal &PressureInterface, + const Array2DReal &SpecVol, + const Array2DReal &ZInterface, + const Array2DReal &LayerThick) const { + + OMEGA_SCOPE(LocCenteredPGrad, CenteredPGrad); + OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad); + OMEGA_SCOPE(LocMinLayerEdgeBot, MinLayerEdgeBot); + OMEGA_SCOPE(LocMaxLayerEdgeTop, MaxLayerEdgeTop); + OMEGA_SCOPE(LocTidalPotential, TidalPotential); + OMEGA_SCOPE(LocSelfAttractionLoading, SelfAttractionLoading); + + if (PressureGradChoice == PressureGradType::Centered) { + + // computes centered geopotential and pressure gradient tendency + parallelForOuter( + "pgrad-centered", {NEdgesAll}, + KOKKOS_LAMBDA(I4 IEdge, const TeamMember &Team) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocCenteredPGrad(Tend, IEdge, KChunk, PressureMid, + PressureInterface, ZInterface, + LocTidalPotential, + LocSelfAttractionLoading, SpecVol); + }); + }); + + } else { + + // computes high-order geopotential and pressure gradient tendency + parallelForOuter( + "pgrad-highorder", {NEdgesAll}, + KOKKOS_LAMBDA(I4 IEdge, const TeamMember &Team) { + const int KMin = LocMinLayerEdgeBot(IEdge); + const int KMax = LocMaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocHighOrderPGrad(Tend, IEdge, KChunk, PressureMid, + PressureInterface, ZInterface, + LocTidalPotential, + LocSelfAttractionLoading, SpecVol); + }); + }); + } +} // end compute pressure gradient + +//------------------------------------------------------------------------------ +// Constructor for centered pressure gradient functor +PressureGradCentered::PressureGradCentered( + const HorzMesh *Mesh, ///< [in] Horizontal mesh + const VertCoord *VCoord ///< [in] Vertical coordinate + ) + : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge), + EdgeMask(VCoord->EdgeMask), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} + +//------------------------------------------------------------------------------ +// Constructor for high order pressure gradient functor +PressureGradHighOrder::PressureGradHighOrder( + const HorzMesh *Mesh, ///< [in] Horizontal mesh + const VertCoord *VCoord ///< [in] Vertical coordinate + ) + : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge), + EdgeMask(VCoord->EdgeMask), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/PGrad.h b/components/omega/src/ocn/PGrad.h new file mode 100644 index 000000000000..9c64b4b86e82 --- /dev/null +++ b/components/omega/src/ocn/PGrad.h @@ -0,0 +1,204 @@ +#ifndef OMEGA_PGRAD_H +#define OMEGA_PGRAD_H +//===-- ocn/PGrad.h - Pressure Gradient -----------------*- C++ -*-===// +/// +/// Implements the PressureGrad class which provides a centered and +/// high-order pressure gradient option and dispatches computations to +/// functor objects. This follows the patterns used in Eos.h/Eos.cpp. +// +//===----------------------------------------------------------------------===// + +#include "Config.h" +#include "Eos.h" +#include "GlobalConstants.h" +#include "HorzMesh.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" +#include + +namespace OMEGA { + +enum class PressureGradType { Centered, HighOrder1, HighOrder2 }; + +// Centered pressure gradient functor +class PressureGradCentered { + public: + bool Enabled; + + // constructor declaration + PressureGradCentered(const HorzMesh *Mesh, ///< [in] Horizontal mesh + const VertCoord *VCoord ///< [in] Vertical coordinate + ); + + // Compute centered pressure gradient contribution for given edge and + // vertical chunk. This appends results into the Tend array (in-place). + KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, + const Array2DReal &PressureMid, + const Array2DReal &PressureInterface, + const Array2DReal &ZInterface, + const Array1DReal &TidalPotential, + const Array1DReal &SelfAttractionLoading, + const Array2DReal &SpecVol) const { + + const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); + + const I4 ICell0 = CellsOnEdge(IEdge, 0); + const I4 ICell1 = CellsOnEdge(IEdge, 1); + const Real InvDcEdge = 1.0_Real / DcEdge(IEdge); + + Real GradGeoPot = + (TidalPotential(ICell1) - TidalPotential(ICell0)) * InvDcEdge + + (SelfAttractionLoading(ICell1) - SelfAttractionLoading(ICell0)) * + InvDcEdge; + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + Real MontPotCell0K = + PressureInterface(ICell0, K) * SpecVol(ICell0, K) + + Gravity * ZInterface(ICell0, K); + Real MontPotCell1K = + PressureInterface(ICell1, K) * SpecVol(ICell1, K) + + Gravity * ZInterface(ICell1, K); + Real GradMontPotK = (MontPotCell1K - MontPotCell0K) * InvDcEdge; + + Real MontPotCell0Kp1 = + PressureInterface(ICell0, K + 1) * SpecVol(ICell0, K) + + Gravity * ZInterface(ICell0, K + 1); + Real MontPotCell1Kp1 = + PressureInterface(ICell1, K + 1) * SpecVol(ICell1, K) + + Gravity * ZInterface(ICell1, K + 1); + Real GradMontPotKp1 = (MontPotCell1Kp1 - MontPotCell0Kp1) * InvDcEdge; + Real GradMontPot = 0.5_Real * (GradMontPotK + GradMontPotKp1); + + Real PGradAlpha = + 0.5_Real * (PressureMid(ICell1, K) + PressureMid(ICell0, K)) * + (SpecVol(ICell1, K) - SpecVol(ICell0, K)) * InvDcEdge; + Tend(IEdge, K) += + EdgeMask(IEdge, K) * (-GradMontPot + PGradAlpha - GradGeoPot); + } + } + + private: + Array2DI4 CellsOnEdge; + Array1DReal DcEdge; + Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; +}; + +// High-order pressure gradient functor (placeholder) +class PressureGradHighOrder { + public: + bool Enabled; + + // constructor declaration + PressureGradHighOrder(const HorzMesh *Mesh, ///< [in] Horizontal mesh + const VertCoord *VCoord ///< [in] Vertical coordinate + ); + + KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, + const Array2DReal &PressureMid, + const Array2DReal &PressureInterface, + const Array2DReal &ZInterface, + const Array1DReal &TidalPotential, + const Array1DReal &SelfAttractionLoading, + const Array2DReal &SpecVol) const { + + // Placeholder: for now, no-op (future high-order implementation) + const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + Tend(IEdge, K) += 0.0_Real; + } + } + + private: + Array2DI4 CellsOnEdge; + Array1DReal DcEdge; + Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; +}; + +// Pressure gradient manager class +class PressureGrad { + public: + // Flag to indicate if pressure gradient term is enabled + bool Enabled; + + // Initialize the default instance + static void init(); + + // Create a new pressure gradient object and add to map + static PressureGrad *create(const std::string &Name, const HorzMesh *Mesh, + const VertCoord *VCoord, Config *Options); + + // Get the default instance + static PressureGrad *getDefault(); + + // Get instance by name + static PressureGrad * + get(const std::string &Name ///< [in] Name of PressureGrad + ); + + // Deallocates arrays and deletes instance + static void clear(); + + // Remove pressure gradient object by name + static void erase(const std::string &Name ///< [in] Name of PressureGrad + ); + + // Destructor + ~PressureGrad(); + + // Compute pressure gradient tendencies and add into Tend array + void computePressureGrad(Array2DReal &Tend, const Array2DReal &PressureMid, + const Array2DReal &PressureInterface, + const Array2DReal &SpecVol, + const Array2DReal &ZInterface, + const Array2DReal &LayerThick) const; + + private: + // Construct a new pressure gradient object + PressureGrad(const HorzMesh *Mesh, const VertCoord *VCoord, Config *Options); + + // forbid copy and move construction + PressureGrad(const PressureGrad &) = delete; + PressureGrad(PressureGrad &&) = delete; + + // Pointer to default pressure gradient object + static PressureGrad *DefaultPGrad; + + // Mesh-related sizes + I4 NEdgesAll = 0; + I4 NEdgesOwned = 0; + I4 NVertLayers = 0; + I4 NVertLayersP1 = 0; + + // Data required for computation (stored copies of VCoord arrays) + Array1DI4 MinLayerEdgeBot; ///< min vertical layer on each edge + Array1DI4 MaxLayerEdgeTop; ///< max vertical layer on each edge + + // Temporary: to be moveed to tidal forcing module in future + Array1DReal TidalPotential; ///< Tidal potential for tidal forcing + Array1DReal + SelfAttractionLoading; ///< Self attraction and loading for tidal forcing + + // Instances of functors + PressureGradCentered CenteredPGrad; + PressureGradHighOrder HighOrderPGrad; + + // Choice from config + PressureGradType PressureGradChoice = PressureGradType::Centered; + + // Map of all pressure gradient objects by name + static std::map> AllPGrads; + +}; // end class PressureGrad + +} // namespace OMEGA +#endif diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index b228c829082b..23259f36c379 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -10,7 +10,11 @@ #include "Tendencies.h" #include "CustomTendencyTerms.h" +#include "Eos.h" #include "Error.h" +#include "Field.h" +#include "OceanState.h" +#include "PGrad.h" #include "Pacer.h" #include "TimeStepper.h" #include "Tracers.h" @@ -32,6 +36,8 @@ void Tendencies::init() { VertCoord *DefVertCoord = VertCoord::getDefault(); VertAdv *DefVertAdv = VertAdv::getDefault(); TimeStepper *DefTimeStepper = TimeStepper::getDefault(); + Eos *DefEos = Eos::getInstance(); + PressureGrad *DefPGrad = PressureGrad::getDefault(); I4 NTracers = Tracers::getNumTracers(); @@ -71,9 +77,9 @@ void Tendencies::init() { TimeInterval TimeStep = DefTimeStepper->getTimeStep(); // Ceate default tendencies - Tendencies::DefaultTendencies = - create("Default", DefHorzMesh, DefVertCoord, DefVertAdv, NTracers, - TimeStep, &TendConfig, CustomThickTend, CustomVelTend); + Tendencies::DefaultTendencies = create( + "Default", DefHorzMesh, DefVertCoord, DefVertAdv, DefPGrad, DefEos, + NTracers, TimeStep, &TendConfig, CustomThickTend, CustomVelTend); DefaultTendencies->readConfig(OmegaConfig); @@ -238,17 +244,76 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config Err += TendConfig.get("EddyDiff4", this->TracerHyperDiff.EddyDiff4); CHECK_ERROR_ABORT(Err, "Tendencies: EddyDiff4 not found in TendConfig"); } + + Err += TendConfig.get("PressureGradTendencyEnable", this->PGrad->Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); } +//------------------------------------------------------------------------------ +// Define fields associated with tendencies +void Tendencies::defineFields() { + std::string LayerThicknessTendFieldName = "LayerThicknessTend"; + std::string NormalVelocityTendFieldName = "NormalVelocityTend"; + std::string TracerTendFieldName = "TracerTend"; + if (Name != "Default") { + LayerThicknessTendFieldName.append(Name); + NormalVelocityTendFieldName.append(Name); + TracerTendFieldName.append(Name); + } + + int NDims = 2; + std::vector DimNamesThickness(NDims); + DimNamesThickness[0] = "NCells"; + DimNamesThickness[1] = "NVertLayers"; + auto LayerThicknessTendField = + Field::create(LayerThicknessTendFieldName, "Layer thickness tendency", + "m/s", "cell_thickness_tendency", -9.99E+10, 9.99E+10, + -9.99E+30, NDims, DimNamesThickness); + NDims = 3; + std::vector DimNamesTracer(NDims); + DimNamesTracer[0] = "NTracers"; + DimNamesTracer[1] = "NCells"; + 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); + NDims = 2; + std::vector DimNamesVelocity(NDims); + DimNamesVelocity[0] = "NEdges"; + DimNamesVelocity[1] = "NVertLayers"; + 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); + + std::string TendGroupName = "Tendencies"; + if (Name != "Default") { + TendGroupName.append(Name); + } + auto TendGroup = FieldGroup::create(TendGroupName); + + TendGroup->addField(LayerThicknessTendFieldName); + TendGroup->addField(NormalVelocityTendFieldName); + TendGroup->addField(TracerTendFieldName); + + LayerThicknessTendField->attachData(LayerThicknessTend); + NormalVelocityTendField->attachData(NormalVelocityTend); + TracerTendField->attachData(TracerTend); + +} // end defineFields + //------------------------------------------------------------------------------ // Construct a new group of tendencies -Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies - const HorzMesh *Mesh, ///< [in] Horizontal mesh - const VertCoord *VCoord, ///< [in] Vertical coordinate - VertAdv *VAdv, ///< [in] Vertical advection - int NTracersIn, ///< [in] Number of tracers - TimeInterval TimeStepIn, ///< [in] Time step - Config *Options, ///< [in] Configuration options +Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies + const HorzMesh *Mesh, ///< [in] Horizontal mesh + VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection + PressureGrad *PGrad, ///< [in] Pressure gradient + Eos *EqState, ///< [in] Equation of state + int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStepIn, ///< [in] Time step + Config *Options, ///< [in] Configuration options CustomTendencyType InCustomThicknessTend, CustomTendencyType InCustomVelocityTend) : Mesh(Mesh), VCoord(VCoord), VAdv(VAdv), ThicknessFluxDiv(Mesh, VCoord), @@ -258,7 +323,7 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies BottomDrag(Mesh, VCoord), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend), - CustomVelocityTend(InCustomVelocityTend) { + CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad) { // Tendency arrays LayerThicknessTend = @@ -268,20 +333,27 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies TracerTend = Array3DReal("TracerTend", NTracersIn, Mesh->NCellsSize, VCoord->NVertLayers); + Name = Name_; + NTracers = NTracersIn; TimeStep = TimeStepIn; + defineFields(); + } // end constructor -Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies - const HorzMesh *Mesh, ///< [in] Horizontal mesh - const VertCoord *VCoord, ///< [in] Vertical coordinate - VertAdv *VAdv, ///< [in] Vertical advection - int NTracersIn, ///< [in] Number of tracers - TimeInterval TimeStepIn, ///< [in] Time step - Config *Options) ///< [in] Configuration options - : Tendencies(Name, Mesh, VCoord, VAdv, NTracersIn, TimeStepIn, Options, - CustomTendencyType{}, CustomTendencyType{}) {} +Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies + const HorzMesh *Mesh, ///< [in] Horizontal mesh + VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection + PressureGrad *PGrad, ///< [in] Pressure gradient + Eos *EqState, ///< [in] Equation of state + int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStepIn, ///< [in] Time step + Config *Options) ///< [in] Configuration options + : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, NTracersIn, + TimeStepIn, Options, CustomTendencyType{}, + CustomTendencyType{}) {} //------------------------------------------------------------------------------ // Compute tendencies for layer thickness equation @@ -357,8 +429,10 @@ void Tendencies::computeThicknessTendenciesOnly( void Tendencies::computeVelocityTendenciesOnly( const OceanState *State, ///< [in] State variables const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + const Array3DReal &TracerArray, ///< [in] Tracer array int ThickTimeLevel, ///< [in] Time level int VelTimeLevel, ///< [in] Time level + int TracerTimeLevel, ///< [in] Time level TimeInstant Time ///< [in] Time ) { @@ -529,6 +603,37 @@ void Tendencies::computeVelocityTendenciesOnly( Pacer::stop("Tend:customVelocityTend", 2); } + // Compute pressure gradient + if (PGrad->Enabled) { + + // Temporary handling of surface pressure + Array1DReal SurfacePressure("SurfacePressure", Mesh->NCellsSize); + deepCopy(SurfacePressure, 0.0_Real); + + Pacer::start("Tend:pressureGradTerm", 2); + Array2DReal LayerThick = State->getLayerThickness(ThickTimeLevel); + VCoord->computePressure(LayerThick, SurfacePressure); + + const auto &PressureMid = VCoord->PressureMid; + const auto &PressureInterface = VCoord->PressureInterface; + Array2DReal Temp = Kokkos::subview(TracerArray, Tracers::IndxTemp, + Kokkos::ALL, Kokkos::ALL); + Array2DReal Salinity = Kokkos::subview(TracerArray, Tracers::IndxSalt, + Kokkos::ALL, Kokkos::ALL); + EqState->computeSpecVol(Temp, Salinity, PressureMid); + + // Temporary: ensure vertical geometric/geopotential fields are updated + // for pressure-gradient tendency calculations. + const auto &SpecVol = EqState->SpecVol; + VCoord->computeZHeight(LayerThick, SpecVol); + + const auto &ZInterface = VCoord->ZInterface; + PGrad->computePressureGrad(LocNormalVelocityTend, PressureMid, + PressureInterface, SpecVol, ZInterface, + LayerThick); + Pacer::stop("Tend:pressureGradTerm", 2); + } + Pacer::stop("Tend:computeVelocityTendenciesOnly", 1); } // end velocity tendency compute @@ -696,15 +801,17 @@ void Tendencies::computeThicknessTendencies( void Tendencies::computeVelocityTendencies( const OceanState *State, ///< [in] State variables const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + const Array3DReal &TracerArray, ///< [in] Tracer array int ThickTimeLevel, ///< [in] Time level int VelTimeLevel, ///< [in] Time level + int TracerTimeLevel, ///< [in] Time level TimeInstant Time ///< [in] Time ) { Pacer::start("Tend:computeVelocityTendencies", 1); AuxState->computeMomAux(State, ThickTimeLevel, VelTimeLevel); - computeVelocityTendenciesOnly(State, AuxState, ThickTimeLevel, VelTimeLevel, - Time); + computeVelocityTendenciesOnly(State, AuxState, TracerArray, ThickTimeLevel, + VelTimeLevel, TracerTimeLevel, Time); Pacer::stop("Tend:computeVelocityTendencies", 1); } @@ -759,13 +866,15 @@ void Tendencies::computeAllTendencies( const Array3DReal &TracerArray, ///< [in] Tracer array int ThickTimeLevel, ///< [in] Time level int VelTimeLevel, ///< [in] Time level + int TracerTimeLevel, ///< [in] Time level TimeInstant Time ///< [in] Time ) { AuxState->computeAll(State, TracerArray, ThickTimeLevel, VelTimeLevel); + computeThicknessTendenciesOnly(State, AuxState, ThickTimeLevel, VelTimeLevel, Time); - computeVelocityTendenciesOnly(State, AuxState, ThickTimeLevel, VelTimeLevel, - Time); + computeVelocityTendenciesOnly(State, AuxState, TracerArray, ThickTimeLevel, + VelTimeLevel, TracerTimeLevel, Time); computeTracerTendenciesOnly(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel, Time); } // end all tendency compute diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 9bc4f0584b5a..4ba53a313e82 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -33,11 +33,14 @@ #include "AuxiliaryState.h" #include "Config.h" +#include "Eos.h" #include "HorzMesh.h" #include "OceanState.h" +#include "PGrad.h" #include "TendencyTerms.h" #include "TimeMgr.h" #include "VertAdv.h" +#include "VertCoord.h" #include #include @@ -69,6 +72,8 @@ class Tendencies { TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; + std::string Name; + // Methods to compute tendency groups void computeThicknessTendencies(const OceanState *State, const AuxiliaryState *AuxState, @@ -76,8 +81,9 @@ class Tendencies { TimeInstant Time); void computeVelocityTendencies(const OceanState *State, const AuxiliaryState *AuxState, + const Array3DReal &TracerArray, int ThickTimeLevel, int VelTimeLevel, - TimeInstant Time); + int TracerTimeLevel, TimeInstant Time); void computeTracerTendencies(const OceanState *State, const AuxiliaryState *AuxState, const Array3DReal &TracerArray, @@ -86,15 +92,17 @@ class Tendencies { void computeAllTendencies(const OceanState *State, const AuxiliaryState *AuxState, const Array3DReal &TracerArray, int ThickTimeLevel, - int VelTimeLevel, TimeInstant Time); + int VelTimeLevel, int TracerTimeLevel, + TimeInstant Time); void computeThicknessTendenciesOnly(const OceanState *State, const AuxiliaryState *AuxState, int ThickTimeLevel, int VelTimeLevel, TimeInstant Time); void computeVelocityTendenciesOnly(const OceanState *State, const AuxiliaryState *AuxState, + const Array3DReal &TracerArray, int ThickTimeLevel, int VelTimeLevel, - TimeInstant Time); + int TracerTimeLevel, TimeInstant Time); void computeTracerTendenciesOnly(const OceanState *State, const AuxiliaryState *AuxState, const Array3DReal &TracerArray, @@ -150,8 +158,10 @@ class Tendencies { // Construct a new tendency object Tendencies(const std::string &Name, ///< [in] Name for tendencies const HorzMesh *Mesh, ///< [in] Horizontal mesh - const VertCoord *VCoord, ///< [in] Vertical coordinate + VertCoord *VCoord, ///< [in] Vertical coordinate VertAdv *VAdv, ///< [in] Vertical advection + PressureGrad *PGrad, ///< [in] Pressure gradient + Eos *EqState, ///< [in] Equation of state int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStep, ///< [in] Time step Config *Options, ///< [in] Configuration options @@ -160,22 +170,30 @@ class Tendencies { Tendencies(const std::string &Name, ///< [in] Name for tendencies const HorzMesh *Mesh, ///< [in] Horizontal mesh - const VertCoord *VCoord, ///< [in] Vertical coordinate + VertCoord *VCoord, ///< [in] Vertical coordinate VertAdv *VAdv, ///< [in] Vertical advection + PressureGrad *PGrad, ///< [in] Pressure gradient + Eos *EqState, ///< [in] Equation of state int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStep, ///< [in] Time step Config *Options ///< [in] Configuration options ); + void defineFields(); + // forbid copy and move construction Tendencies(const Tendencies &) = delete; Tendencies(Tendencies &&) = delete; - const HorzMesh *Mesh; ///< Pointer to horizontal mesh - const VertCoord *VCoord; ///< Pointer to vertical coordinate - VertAdv *VAdv; ///< Pointer to vertical advection - I4 NTracers; ///< Number of tracers - TimeInterval TimeStep; ///< Time step + const HorzMesh *Mesh; ///< Pointer to horizontal mesh + VertCoord *VCoord; ///< Pointer to vertical coordinate + VertAdv *VAdv; ///< Pointer to vertical advection + CustomTendencyType CustomThicknessTend; + CustomTendencyType CustomVelocityTend; + Eos *EqState; ///< Pointer to equation of state + PressureGrad *PGrad; ///< Pointer to pressure gradient + I4 NTracers; ///< Number of tracers + TimeInterval TimeStep; ///< Time step // Pointer to default tendencies static Tendencies *DefaultTendencies; @@ -183,9 +201,6 @@ class Tendencies { // Map of all tendency objects static std::map> AllTendencies; - CustomTendencyType CustomThicknessTend; - CustomTendencyType CustomVelocityTend; - }; // end class Tendencies } // namespace OMEGA diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index 12c53a59a436..e72719a3b2f4 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -56,8 +56,8 @@ void ForwardBackwardStepper::doStep( CurLevel, TimeStep); // R_u^{n+1} = RHS_u(u^{n}, h^{n+1}, t^{n+1}) - Tend->computeVelocityTendencies(State, AuxState, NextLevel, CurLevel, - SimTime + TimeStep); + Tend->computeVelocityTendencies(State, AuxState, NextTracerArray, NextLevel, + CurLevel, NextLevel, SimTime + TimeStep); // u^{n+1} = u^{n} + R_u^{n+1} updateVelocityByTend(State, NextLevel, State, CurLevel, TimeStep); diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index f4b697529db0..47ffc54b71fe 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -37,7 +37,7 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state // q = (h,u,phi) // R_q^{n} = RHS_q(u^{n}, h^{n}, phi^{n}, t^{n}) Tend->computeAllTendencies(State, AuxState, CurTracerArray, CurLevel, - CurLevel, SimTime); + CurLevel, CurLevel, SimTime); // q^{n+0.5} = q^{n} + 0.5*dt*R_q^{n} updateStateByTend(State, NextLevel, State, CurLevel, 0.5 * TimeStep); @@ -49,7 +49,7 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state // R_q^{n+0.5} = RHS_q(u^{n+0.5}, h^{n+0.5}, phi^{n+0.5}, t^{n+0.5}) Tend->computeAllTendencies(State, AuxState, NextTracerArray, NextLevel, - NextLevel, SimTime + 0.5 * TimeStep); + NextLevel, NextLevel, SimTime + 0.5 * TimeStep); // q^{n+1} = q^{n} + dt*R_q^{n+0.5} updateStateByTend(State, NextLevel, State, CurLevel, TimeStep); diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index 1b882e7631b9..1160a4ab7213 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -85,7 +85,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state if (Stage == 0) { weightTracers(NextTracerArray, CurTracerArray, State, CurLevel); Tend->computeAllTendencies(State, AuxState, CurTracerArray, CurLevel, - CurLevel, StageTime); + CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, CurLevel, RKB[Stage] * TimeStep); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); @@ -111,7 +111,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state Pacer::stop("RK4:haloExchProvis", 3); Tend->computeAllTendencies(ProvisState, AuxState, ProvisTracers, - CurLevel, CurLevel, StageTime); + CurLevel, CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, NextLevel, RKB[Stage] * TimeStep); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 271af1c40d46..301a04ede49c 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -346,6 +346,17 @@ add_omega_test( "-n;8" ) +################## +# Pressure gradient test +################## + +add_omega_test( + PGRAD_TEST + testPGrad.exe + ocn/PGradTest.cpp + "-n;1" +) + ################## # State test ################## diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index ff89d97f71b5..451ca4293edf 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -23,7 +23,9 @@ #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" @@ -103,6 +105,11 @@ void initIOStreamTest(Clock *&ModelClock // Model clock if (TmpErr != 0) ABORT_ERROR("IOStreamTest: Error initializing OceanState"); + PressureGrad::init(); + + // Intialize Tendencies + Tendencies::init(); + // Initialize Tracers Tracers::init(); @@ -220,6 +227,8 @@ int main(int argc, char **argv) { // Clean up environments TimeStepper::clear(); + PressureGrad::clear(); + Tendencies::clear(); Tracers::clear(); OceanState::clear(); AuxiliaryState::clear(); diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 5202594abedb..aaf7b1d96e80 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -51,9 +51,9 @@ const Real GswBVFExpValue = 0.02081197958166906; // Expected value from GSW-C library /// Test input values -const Real Sa = 30.0; // Absolute Salinity in g/kg -const Real Ct = 10.0; // Conservative Temperature in degC -const Real P = 1000.0; // Pressure in dbar +const Real Sa = 30.0; // Absolute Salinity in g/kg +const Real Ct = 10.0; // Conservative Temperature in degC +const Real P = 1000.0 * Db2Pa; // Pressure in Pa const I4 KDisp = 1; // Displate parcel to K=1 for TEOS-10 eos const Real RTol = 1e-10; // Relative tolerance for isApprox checks @@ -540,17 +540,17 @@ void testBruntVaisalaFreqSqTeos10() { ZMid(ICell, 1) = -993.1071379053125_Real; SArray(ICell, 1) = Sa; TArray(ICell, 1) = Ct + 10.0_Real; - PArray(ICell, 1) = P + 1.0_Real; + PArray(ICell, 1) = (P * Pa2Db + 1.0_Real) * Db2Pa; } else if (K == 2) { ZMid(ICell, 2) = -994.0968821072275_Real; SArray(ICell, 2) = Sa + 1.0_Real; TArray(ICell, 2) = Ct + 5.0_Real; - PArray(ICell, K) = P + 2.0_Real; + PArray(ICell, K) = (P * Pa2Db + 2.0_Real) * Db2Pa; } else { // fill rest with valid junk to avoid Nans and Inf ZMid(ICell, K) = -994.0968821072275_Real - 0.1_Real * K; SArray(ICell, K) = Sa + 1.0_Real + 0.1_Real * K; TArray(ICell, K) = Ct + 5.0_Real - 0.01_Real * K; - PArray(ICell, K) = P + 2.0_Real + 0.1_Real * K; + PArray(ICell, K) = (P * Pa2Db + 2.0_Real + 0.1_Real * K) * Db2Pa; } }); @@ -647,7 +647,7 @@ void checkValueGswcSpecVol() { const Real RTol = 1e-10; /// Get specific volume from GSW-C library - double SpecVol = gsw_specvol(Sa, Ct, P); + double SpecVol = gsw_specvol(Sa, Ct, P * Pa2Db); /// Check the value against the expected TEOS-10 value bool Check = isApprox(SpecVol, TeosSVExpValue, RTol); if (!Check) { @@ -667,8 +667,9 @@ void checkValueGswcN2() { // Input arrays: length nz+1 double Salt[4] = {Sa - 1.0, Sa, Sa + 1.0}; // Absolute Salinity (g/kg) double Temp[4] = {Ct + 15.0, Ct + 10.0, - Ct + 5.0}; // Conservative Temperature (deg C) - double Press[4] = {P, P + 1.0, P + 2.0}; // Pressure (dbar) + Ct + 5.0}; // Conservative Temperature (deg C) + double Press[4] = {P * Pa2Db, P * Pa2Db + 1.0, + P * Pa2Db + 2.0}; // Pressure (dbar) // Latitude (degrees north) double Latitude[4] = {0.0, 0.0, 0.0}; diff --git a/components/omega/test/ocn/PGradTest.cpp b/components/omega/test/ocn/PGradTest.cpp new file mode 100644 index 000000000000..09887a6d9ce6 --- /dev/null +++ b/components/omega/test/ocn/PGradTest.cpp @@ -0,0 +1,334 @@ +//===-- Test driver for OMEGA Pressure Gradient (PGrad) --------------*- +// C++-*-===/ +// +/// \file +/// \brief Test driver for PressureGrad module +// +//===----------------------------------------------------------------------===/ + +#include "PGrad.h" + +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "Eos.h" +#include "Error.h" +#include "Field.h" +#include "GlobalConstants.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 "TimeStepper.h" +#include "Tracers.h" +#include "VertCoord.h" +#include "mpi.h" + +using namespace OMEGA; + +void initPGradTest() { + + Error Err; + int Err1; + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + // Initialize the Logging system + initLogging(DefEnv); + + // Read default config if present + Config("Omega"); + Config::readAll("omega.yml"); + + // First step of time stepper initialization needed for IOstream + TimeStepper::init1(); + + // Initialize the IO system + IO::init(DefComm); + + // Create the default decomposition (initializes the decomposition) + Decomp::init(); + + // Initialize streams + IOStream::init(); + + // Initialize the default halo + Err1 = Halo::init(); + if (Err1 != 0) { + LOG_ERROR("PGrad: error initializing default halo"); + Err += Error(ErrorCode::Fail, "PGrad: error initializing default halo"); + } + + // Initialize the default mesh + HorzMesh::init(); + + // Initialize the default vertical coordinate + VertCoord::init(); + + // Initialize the equation of state + Eos::init(); + + // Initialize ocean state + OceanState::init(); + + // Initialize tracers + Tracers::init(); + + CHECK_ERROR_ABORT(Err, "PGrad: error during initialization"); +} + +int main(int argc, char *argv[]) { + int RetVal = 0; + + MPI_Init(&argc, &argv); + Kokkos::initialize(); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + + { + initPGradTest(); + // Initialize default PressureGrad + PressureGrad::init(); + + HorzMesh *DefMesh = HorzMesh::getDefault(); + VertCoord *VCoord = VertCoord::getDefault(); + OceanState *DefState = OceanState::getDefault(); + Eos *DefEos = Eos::getInstance(); + + // create arrays: Tend on edges, Pressure/Geopotential/SpecVol on cells + Array2DReal Tend("Tend", DefMesh->NEdgesSize, VCoord->NVertLayers); + Array2DReal SpecVolOld("SpecVolOld", DefMesh->NCellsSize, + VCoord->NVertLayers); + Array2DReal PressureMidOld("PressureMidOld", DefMesh->NCellsSize, + VCoord->NVertLayers); + Array1DReal SurfacePressure("SurfacePressure", DefMesh->NCellsSize); + + I4 NEdgesAll = DefMesh->NEdgesAll; + I4 NCellsAll = DefMesh->NCellsAll; + + I4 NVertLayers = 60; + Real DC = 30000.0_Real; + I4 NRefinements = 4; + HostArray1DReal Rmse("Rmse", NRefinements); + for (int Refinement = 0; Refinement < NRefinements; ++Refinement) { + + LOG_INFO("PGradTest: Starting refinement level {}", Refinement); + VCoord->NVertLayers = NVertLayers; + VCoord->NVertLayersP1 = NVertLayers + 1; + + auto &MinLayerCell = VCoord->MinLayerCell; + auto &MaxLayerCell = VCoord->MaxLayerCell; + parallelFor( + {NCellsAll}, KOKKOS_LAMBDA(int i) { + MinLayerCell(i) = 0; + MaxLayerCell(i) = NVertLayers - 1; + }); + + auto &MinLayerEdgeBot = VCoord->MinLayerEdgeBot; + auto &MaxLayerEdgeTop = VCoord->MaxLayerEdgeTop; + parallelFor( + {NEdgesAll}, KOKKOS_LAMBDA(int i) { + MinLayerEdgeBot(i) = 0; + MaxLayerEdgeTop(i) = NVertLayers - 1; + }); + + auto &CellsOnEdge = DefMesh->CellsOnEdge; + auto &DcEdge = DefMesh->DcEdge; + parallelFor( + {NEdgesAll}, KOKKOS_LAMBDA(int i) { + CellsOnEdge(i, 0) = 0; + CellsOnEdge(i, 1) = 1; + DcEdge(i) = DC; + }); + + // Fetch reference desnity from Config + Real Density0; + Density0 = RhoSw; + + I4 TimeLevel = 0; + + // get state and tracer arrays + Array2DReal LayerThick = DefState->getLayerThickness(TimeLevel); + Array2DReal Temp = Tracers::getByName(TimeLevel, "Temperature"); + Array2DReal Salinity = Tracers::getByName(TimeLevel, "Salinity"); + + // set Z interface and mid-point locations + Real ZBottom = -1000.0_Real; + Real DZ = 2.0_Real * (-ZBottom / NVertLayers); + auto &BottomDepth = VCoord->BottomDepth; + auto &ZInterface = VCoord->ZInterface; + auto &ZMid = VCoord->ZMid; + Real TiltFactor = 0.495_Real; + parallelFor( + {NCellsAll}, KOKKOS_LAMBDA(int i) { + ZInterface(i, NVertLayers) = ZBottom; + SurfacePressure(i) = 0.0_Real; + BottomDepth(i) = 0.0_Real; + for (int k = NVertLayers - 1; k >= 0; --k) { + Real X = (k + i) % 2; + Real Dz = (2.0_Real * TiltFactor - 1.0_Real) * X * DZ + + (1.0_Real - TiltFactor) * + DZ; // staggered layer thickness + ZInterface(i, k) = ZInterface(i, k + 1) + Dz; + LayerThick(i, k) = ZInterface(i, k) - ZInterface(i, k + 1); + ZMid(i, k) = + 0.5_Real * (ZInterface(i, k) + ZInterface(i, k + 1)); + BottomDepth(i) += Dz; + } + }); + + LOG_INFO("NVertLayers = {}", NVertLayers); + LOG_INFO("dC = {}", DC); + DefState->copyToHost(0); + HostArray2DReal LayerThickH = DefState->getLayerThicknessH(TimeLevel); + for (int i = 0; i < 2; ++i) { + for (int k = 0; k < 2; ++k) { + LOG_INFO("LayerThick({}, {}) = {}", i, k, LayerThickH(i, k)); + } + } + + // set simple temperature and salinity profiles + auto &SpecVol = DefEos->SpecVol; + parallelFor( + {NCellsAll, NVertLayers}, KOKKOS_LAMBDA(int i, int k) { + Real T0 = 30.0; + Real TB = 5.0; + Real S0 = 30.0; + Real SB = 40.0; + + Real Phi0 = (ZMid(i, k) - ZBottom) / (-ZBottom); + Real PhiB = 1.0_Real - Phi0; + + Temp(i, k) = T0 * Phi0 + TB * PhiB; + Salinity(i, k) = S0 * Phi0 + SB * PhiB; + SpecVol(i, k) = 1.0_Real / Density0; + SpecVolOld(i, k) = SpecVol(i, k); + }); + + // Iterate to converge LayerThick, SpecVol, PressureMid + auto &PressureMid = VCoord->PressureMid; + VCoord->computePressure(LayerThick, SurfacePressure); + deepCopy(PressureMidOld, PressureMid); + for (int Iteration = 0; Iteration < 15; ++Iteration) { + + // compute specific volume from EOS + VCoord->computePressure(LayerThick, SurfacePressure); + DefEos->computeSpecVol(Temp, Salinity, PressureMid); + + // compute psuedo thickness from specific volume + parallelFor( + {NCellsAll, NVertLayers}, KOKKOS_LAMBDA(int i, int k) { + LayerThick(i, k) = 1.0_Real / (SpecVol(i, k) * Density0) * + (ZInterface(i, k) - ZInterface(i, k + 1)); + }); + + // compute difference from previous iteration + Real MaxValue = 0.0_Real; + parallelReduce( + {NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(int i, int k, Real &max) { + Real Diff = Kokkos::abs(SpecVol(i, k) - SpecVolOld(i, k)); + if (Diff > max) + max = Diff; + }, + Kokkos::Max(MaxValue)); + + // check convergence + if (MaxValue < 1e-12_Real) { + LOG_INFO("converged: max diff = {}", MaxValue); + break; + } else { + parallelFor( + {NCellsAll, NVertLayers}, KOKKOS_LAMBDA(int i, int k) { + SpecVolOld(i, k) = SpecVol(i, k); + }); + } + } + + // compute pressure once more with converged LayerThick + VCoord->computePressure(LayerThick, SurfacePressure); + + // compute z levels + VCoord->computeZHeight(LayerThick, SpecVol); + + // get PressureGrad instance + PressureGrad *DefPGrad = PressureGrad::getDefault(); + if (!DefPGrad) { + LOG_INFO("PGrad: default instance not created by init"); + } + + // compute pressure gradient + deepCopy(Tend, 0.0_Real); + + const auto &PressureInterface = VCoord->PressureInterface; + DefPGrad->computePressureGrad(Tend, PressureMid, PressureInterface, + SpecVol, ZInterface, LayerThick); + + // compute errors + Real MaxValue = 0.0_Real; + parallelReduce( + {NEdgesAll, NVertLayers - 2}, + KOKKOS_LAMBDA(int i, int k, Real &max) { + Real Val = Kokkos::abs(Tend(i, k + 1)); + if (Val > max) + max = Val; + }, + Kokkos::Max(MaxValue)); + Real SumValue = 0.0_Real; + parallelReduce( + {NEdgesAll, NVertLayers - 2}, + KOKKOS_LAMBDA(int i, int k, Real &LSum) { + LSum += Tend(i, k + 1) * Tend(i, k + 1); + }, + Kokkos::Sum(SumValue)); + Real RmseVal = std::sqrt(SumValue / (NEdgesAll * (NVertLayers - 2))); + Rmse(Refinement) = RmseVal; + + LOG_INFO("refinement level {}: max |Tend| = {}, average Tend = {}", + Refinement, MaxValue, RmseVal); + + // coarsen for next iteration + DC = DC * 2.0_Real; + NVertLayers = NVertLayers / 2; + + } // refinement loop + + // Test for second order convergence + // resolution (dC) increases in refimenent loop + if (Rmse(0) < Rmse(NRefinements - 1) / pow(4.0_Real, NRefinements - 1)) { + RetVal = 0; + } else { + RetVal = 1; + } + + // cleanup + PressureGrad::clear(); + IOStream::finalize(); + TimeStepper::clear(); + Tracers::clear(); + Eos::destroyInstance(); + OceanState::clear(); + VertCoord::clear(); + Field::clear(); + Dimension::clear(); + HorzMesh::clear(); + Halo::clear(); + Decomp::clear(); + MachEnv::removeAll(); + } + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + if (RetVal >= 256) + RetVal = 255; + return RetVal; +} diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp index 0c2f631db277..e93a32ad2030 100644 --- a/components/omega/test/ocn/StateTest.cpp +++ b/components/omega/test/ocn/StateTest.cpp @@ -23,7 +23,9 @@ #include "MachEnv.h" #include "OceanState.h" #include "OmegaKokkos.h" +#include "PGrad.h" #include "Pacer.h" +#include "Tendencies.h" #include "TimeStepper.h" #include "VertCoord.h" #include "mpi.h" @@ -90,6 +92,9 @@ void initStateTest() { // Initialize Aux State variables AuxiliaryState::init(); + // Initialize pressure gradient + PressureGrad::init(); + // Create tendencies Tendencies::init(); @@ -399,6 +404,7 @@ int main(int argc, char *argv[]) { OceanState::clear(); Tracers::clear(); AuxiliaryState::clear(); + PressureGrad::clear(); Tendencies::clear(); TimeStepper::clear(); HorzMesh::clear(); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 40afa13e531f..c87828f6bccf 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -1,9 +1,11 @@ #include "Tendencies.h" #include "AuxiliaryState.h" #include "Config.h" +#include "CustomTendencyTerms.h" #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Eos.h" #include "Error.h" #include "Field.h" #include "GlobalConstants.h" @@ -15,6 +17,7 @@ #include "MachEnv.h" #include "OceanTestCommon.h" #include "OmegaKokkos.h" +#include "PGrad.h" #include "Pacer.h" #include "TimeStepper.h" #include "VertCoord.h" @@ -57,6 +60,11 @@ int initState() { auto *VCoord = VertCoord::getDefault(); auto *State = OceanState::getDefault(); + // Define tendency fields + int NDims = 2; + std::vector DimNamesThickness(NDims); + DimNamesThickness[0] = "NCells"; + Array2DReal LayerThickCell = State->getLayerThickness(0); Array2DReal NormalVelEdge = State->getNormalVelocity(0); @@ -124,6 +132,8 @@ int initTendenciesTest(const std::string &mesh) { VertCoord::init(); Tracers::init(); VertAdv::init(); + PressureGrad::init(); + Eos::init(); int StateErr = OceanState::init(); if (StateErr != 0) { @@ -138,6 +148,7 @@ int initTendenciesTest(const std::string &mesh) { int testTendencies() { int Err = 0; + Error Err1; // test initialization Tendencies::init(); @@ -156,13 +167,20 @@ int testTendencies() { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); const auto VAdv = VertAdv::getDefault(); + const auto PGrad = PressureGrad::getDefault(); + const auto EqState = Eos::getInstance(); VCoord->NVertLayers = 12; + // test creation of another tendencies TimeInterval ZeroTimeStep; // Zero-length time step placeholder Config *Options = Config::getOmegaConfig(); - Tendencies::create("TestTendencies", Mesh, VCoord, VAdv, 3, ZeroTimeStep, - Options); + Config TendConfig("Tendencies"); + Err1 = Options->get(TendConfig); + int NTracersTest = 3; + + Tendencies::create("TestTendencies", Mesh, VCoord, VAdv, PGrad, EqState, + NTracersTest, ZeroTimeStep, &TendConfig); // test retrievel of another tendencies if (Tendencies::get("TestTendencies")) { @@ -195,9 +213,11 @@ int testTendencies() { Array3DReal TracerArray = Tracers::getAll(0); int ThickTimeLevel = 0; int VelTimeLevel = 0; + int TracerTimeLevel = 0; TimeInstant Time; DefTendencies->computeAllTendencies(State, AuxState, TracerArray, - ThickTimeLevel, VelTimeLevel, Time); + ThickTimeLevel, VelTimeLevel, + TracerTimeLevel, Time); // check that everything got computed correctly int NCellsOwned = Mesh->NCellsOwned; @@ -235,6 +255,8 @@ int testTendencies() { void finalizeTendenciesTest() { Tracers::clear(); + PressureGrad::clear(); + Eos::destroyInstance(); AuxiliaryState::clear(); OceanState::clear(); VertAdv::clear(); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 95aca280fa42..0de465c17fe4 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -19,6 +19,7 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Eos.h" #include "Error.h" #include "Field.h" #include "Halo.h" @@ -28,6 +29,7 @@ #include "MachEnv.h" #include "OceanState.h" #include "OmegaKokkos.h" +#include "PGrad.h" #include "Pacer.h" #include "TendencyTerms.h" #include "TimeMgr.h" @@ -173,6 +175,8 @@ int initTimeStepperTest(const std::string &mesh) { auto *DefVAdv = VertAdv::getDefault(); AuxiliaryState::init(); + Eos::init(); + PressureGrad::init(); Tendencies::init(); // finish initializing default time stepper @@ -187,8 +191,10 @@ int initTimeStepperTest(const std::string &mesh) { // Creating non-default state and auxiliary state to use only one vertical // layer - auto *DefMesh = HorzMesh::getDefault(); - auto *DefHalo = Halo::getDefault(); + auto *DefMesh = HorzMesh::getDefault(); + auto *DefHalo = Halo::getDefault(); + auto *DefEos = Eos::getInstance(); + auto *DefPGrad = PressureGrad::getDefault(); int NTracers = Tracers::getNumTracers(); const int NTimeLevels = 2; @@ -216,8 +222,9 @@ int initTimeStepperTest(const std::string &mesh) { // Creating non-default tendencies with custom velocity tendencies auto *TestTendencies = Tendencies::create( - "TestTendencies", DefMesh, DefVertCoord, DefVAdv, NTracers, ZeroTimeStep, - &Options, Tendencies::CustomTendencyType{}, DecayVelocityTendency{}); + "TestTendencies", DefMesh, DefVertCoord, DefVAdv, DefPGrad, DefEos, + NTracers, ZeroTimeStep, &Options, Tendencies::CustomTendencyType{}, + DecayVelocityTendency{}); if (!TestTendencies) { Err++; LOG_ERROR("TimeStepperTest: error creating test tendencies"); @@ -275,6 +282,8 @@ void timeLoop(TimeInstant TimeStart, Real TimeEnd) { void finalizeTimeStepperTest() { Tracers::clear(); TimeStepper::clear(); + PressureGrad::clear(); + Eos::destroyInstance(); Tendencies::clear(); AuxiliaryState::clear(); OceanState::clear();