Vertical Advection
1 Overview
The Vertical Advection module is responsible for the vertical transport of mass, momentum, and tracers within Omega. This implementation is consistent with Omega’s non-Boussinesq governing equations, where the prognostic vertical coordinate is a pseudo-height, effectively a normalized pressure coordinate. The methods compute the tendencies associated with the vertical advection terms in the governing equations, and are designed to allow for flexibility in numerical scheme configuration.
2 Requirements
2.1 Requirement: Compute transport pseudo-velocity
Compute the transport pseudo-velocity through vertical layer interfaces based on the divergence of the horizontal velocity field. Allow for possible movement of layer interfaces and tilt of layer interfaces
2.2 Requirement: Advection scheme flexibility
Support for multiple possible reconstruction schemes, allowing users to select among different orders of accuracy and upwinding at runtime through configuration options. Standard and monotonic flux-corrected transport (FCT) advection algorithms are available for the tendency calculations.
2.3 Requirement: Compute tendencies for mass, momentum, and tracers
Separate methods are neededto compute the vertical-advection tendencies
for mass, momentum, and tracers. Mass (LayerThickness) and horizontal
momentum (NormalVelocity) are stored as 2-D arrays. LayerThickness is
cell-centered, while NormalVelocity is defined on edges. These tendencies
use a first-order flux formulation that ensures strict conservation and energy
consistency with the continuity equation. Tracers is a 3-D, cell-centered
array. The tracer tendencies are computed using a configurable flux-form scheme
that supports multiple interface reconstruction options.
2.4 Requirement: Modular design
The vertical advection module must consist of classes and methods that are easily incorporated into different time-integration schemes and maintain compatibility with the Omega driver framework.
2.5 Desired: Monotonicity checks
The vertical advection module should include optional monotonicity checks to ensure that reconstructed tracer fields remain physically bounded and free of spurious oscillations when using the FCT scheme.
3 Algorithmic Formulation
The Omega V1 Governing Equations describe the evolution of mass, momentum, and tracers, including the contributions from vertical advection.
3.1 Pseudo-velocity
The vertical pseudo-velocity is derived from the divergence of the horizontal velocity \(d\tilde{w} / d\tilde{z} = - \nabla \cdot {\bf u}.\) In discrete form, the pseudo-velocity at the top interface of layer \(k\) in cell \(i\) is computed from the thickness-weighted divergence of horizontal velocity:
using the TRiSK formulation of the discrete divergence operator (see V0 design document). The boundary conditions at the top and bottom interfaces set the pseudo-velocity to zero.
3.2 Mass (Thickness)
The time evolution of the pseudo-thickness (represented in the code as
LayerThickness) due to vertical transport is
where \(\tilde{W}_{tr}\) represents the transport pseudo-velocity and accounts for the motion of the interface, contributions from the horizontal velocity through tilted interfaces, and surface source terms.
This first-order form ensures exact conservation of total column thickness and numerical stability.
3.3 Momentum (Velocity)
The momentum equation governs the evolution of the horizontal velocity field
(NormalVelocity),
where \(U = u - u_k\) removes the component of the horizontal velocity normal to a tilted interface. Square brackets \([\,]\) denote interpolation of quantities from cell centers to edges or from the middle of a layer to an interface.
3.4 Tracers
The contribution of vertical advection to tracer transport is expressed as
A variety of options are available for the scheme used to reconstruct the tracer values at interfaces \([\phi]_{k\pm1/2}\) (order of accuracy, centered vs. upwinded) which can be selected at runtime via configuration.
4 Design
4.1 Data types and parameters
4.1.1 Parameters
An enum class is used for the selection of the reconstruction scheme for the
fluxes at layer interfaces:
enum class VertFluxOption {
Second, /// 2nd-order centered
Third, /// 3rd-order upwind
Fourth /// 4th-order centered
};
Another enum class specifies the advection scheme:
enum class VertAdvOption {
Standard, /// standard non-monotonic scheme
FCT, /// flux-corrected transport scheme
};
4.1.2 Class/structs/data types
Separate functors are defined for each of the flux reconstruction options described in section 4.1.1, as well as for the low-order (first-order upwind) flux used in the flux-corrected transport scheme. These functors are invoked in parallel dispatches within the tracer tendency compute methods.
class VertFlux1stOrderUpwind {
public:
VertFlux1stOrderUpwind(const VertCoord *VCoord);
KOKKOS_FUNCTION void operator()(const Array3DReal &Flux, I4 ICell,
I4 KChunk, Array3DReal &Tracers,
Array2DReal &WTransp)
};
The VertAdv class provides the core functionality for vertical advection.
It defines the data structures, configuration parameters, and compute methods
needed for performing advective vertical transport. Instances of the flux
functors are created and managed within this class.
class VertAdv{
public:
// Logicals to enable/disable advection of specific fields
bool ThickVertAdvEnabled;
bool VelVertAdvEnabled;
bool TracerVertAdvEnabled;
// Enum options
VertFluxOption VertFluxChoice;
VertAdvOption VertAdvChoice;
// Core data arrays
Array2DReal VerticalTransport /// pseudo-velocity through top of layer
Array3DReal VerticalFlux /// fluxes at vertical interfaces
Array3DReal LowOrderVerticalFlux /// low order fluxes for FCT
// Compute methods
computeVerticalTransport(...);
computeThicknessVAdvTend(...);
computeVelocityVadvTend(...);
computeTracerVadvTend(...);
computeStdVAdvTend(...);
computeFCTVAdvTend(...);
// Flux functors
VertFlux1stOrderUpwind VFlux1stUpw
VertFlux2ndOrderCentered VFlux2ndCent
VertFlux3rdOrderUpwind VFlux3rdUpw
VertFlux4thOrderCentered VFlux4thCent
private:
// Mesh dimensions
I4 NVertLayers;
I4 NVertLayersP1;
I4 NCellsOwned;
I4 NCellsAll;
I4 NCellsSize;
I4 NEdgesOwned;
I4 NEdgesAll;
I4 NEdgesSize;
// Vertical loop bounds from VertCoord
Array1DI4 MinLayerCell;
Array1DI4 MaxLayerCell;
Array1DI4 MinLayerEdgeBot;
Array1DI4 MaxLayerEdgeTop;
// Arrays from HorzMesh
Array2DI4 CellsOnEdge;
Array2DI4 EdgesOnCell;
Array1DI4 NEdgesOnCell;
Array1DReal AreaCell;
Array1DReal DvEdge;
Array2DReal EdgeSignOnCell;
};
4.2 Methods
4.2.1 Creation, Destruction, Retrieval
Like other Omega components, instances of the VertAdv class will be stored in
a static map and a DefVAdv pointer will point to the default instance.
VertAdv *VertAdv::DefVAdv;
std::map<std::string, std::unique_ptr<VertAdv>> VertAdv::AllVAdv
Methods are required for construction, destruction, and retrieval. After
construction, object pointers are inserted in the map; upon destruction, they
are removed. Utility functions provide access to the default VertAdv instance
and any non-default instances.
An initialization method init is responsible for reading configuration
parameters and creating the default VertAdv object.
void VertAdv::init();
4.2.2 Pseudo-velocity computation
The quantities required by this method are either stored locally in the
VertAdv object, or are passed in through the State and AuxState objects.
This method computes the vertical pseudo-velocity based on the divergence of
the horizontal velocity field. A parallelForOuter construct iterates over the
cells in the horizontal mesh. Within each cell, a parallelForInner iterates
through the active layers to compute the divergence of horizontal velocity,
storing the result in scratch space. After the divergence is computed, a
parallelScanInner performs a prefix sum to obtain the vertical pseudo-velocity
at each layer interface. These pseudo-velocities are stored in the
VerticalTransport member array for later use in mass, momentum, and tracer
tendency calculations.
void VertAdv::computeVerticalTransport(const OceanState *State,
const AuxiliaryState *AuxState, int ThickTimeLevel,
int VelTimeLevel, TimeInterval Dt)
4.2.3 Tendency computations
The computeThicknessVAdvTend, computeVelocityVadvTend, and
computeTracerVadvTend methods are designed to be called by the
Tendencies class methods.
The computeThicknessVAdvTend method computes the vertical advection tendency
for pseudo-thickness (LayerThickness). Because the vertical pseudo-velocity
field (VerticalTransport) is already stored within the VertAdv object,
this method requires only the thickness tendency array as an argument. The
parallel execution of the computation is straightforward.
void VertAdv::computeThicknessVAdvTend(const Array2DReal &Tend) {
if (!ThickVertAdvEnabled) return;
OmegaScope(LocVertTransp, VerticalTransport);
parallelForOuter(...,
parallelForInner(...,
Tend(ICell, K) += LocVertTransp(ICell, K + 1) -
LocVertTransp(ICell, K);
);
);
}
The computeVelocityVadvTend method computes the vertical advection tendency
for horizontal velocity (NormalVelocity). The method requires the
NormalVelocity field from the State object and the FluxLayerThickEdge
field from the AuxState object as input. A parallelForOuter loop iterates
over the edges in the horizontal mesh, and a parallelForInner loop iterates
over the active layer interfaces. At each interface, NormalVelocity and
FluxLayerThickEdge are used to compute du/dz, and the VerticalTransport
is interpolated from cell centers to edges. The product of these, w du/dz is
averaged from interfaces to the middle of each layer and accumulated into the
velocity tendency.
void VertAdv::computeVelocityVAdvTend(const Array2DReal &Tend,
const Array2DReal &NormalVelocity, const Array2DReal &FluxLayerThickEdge)
The computeTracerVAdvTend method serves as an interface that dispatches to
either the standard advection algorithm, or the FCT scheme depending on
runtime configuration.
void VertAdv::computeTracerVadvTend(const Array3DReal &Tend,
const Array3DReal &Tracers)
The computeStdVAdvTend method performs the standard advection update by
looping over the cells and interfaces using nested parallelForOuter and
parallelForInner constructs. The flux functor chosen from the runtime
configuration is used to compute interface fluxes, which are then applied to
update the tracer tendencies.
The computeFCTVAdvTend method implements the flux-corrected transport scheme
following Zalesak 1979. At each interface, both a high-order flux (chosen via
configuration) and a low-order flux (1st-order upwind) are computed. These
fluxes are then blended and limited to ensure monotone transport.
5 Verification and Testing
Unit tests are required for each compute method to verify correct operation and numerical consistency. In addition, convergence tests should confirm the expected order of accuracy of the flux functors.