Add centered pressure gradient

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]

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");
+```

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

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

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();
 

components/omega/src/ocn/OceanFinal.cpp

@@ -7,13 +7,15 @@
 
 #include "AuxiliaryState.h"
 #include "Decomp.h"
+#include "Eos.h"
 #include "Field.h"
 #include "Halo.h"
 #include "HorzMesh.h"
 #include "IO.h"
 #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();

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();