From 2790c93e30391d3424d5f7fd7f0ca6b35badfeca Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 23 Oct 2025 08:53:01 -0700 Subject: [PATCH 1/2] Add vertical advection design document --- components/omega/doc/design/VertAdv.md | 343 +++++++++++++++++++++++++ components/omega/doc/index.md | 1 + 2 files changed, 344 insertions(+) create mode 100644 components/omega/doc/design/VertAdv.md diff --git a/components/omega/doc/design/VertAdv.md b/components/omega/doc/design/VertAdv.md new file mode 100644 index 000000000000..8e500d73b6bc --- /dev/null +++ b/components/omega/doc/design/VertAdv.md @@ -0,0 +1,343 @@ +(omega-design-vert-adv)= +# 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](omega-design-governing-eqns-omega1) 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: + +$$ +\tilde{w}_{i,k-1/2} = \tilde{w}_{i,k+1/2} + +\frac{1}{A_i} \sum_{e\in EC(i)} n_{e,i} [\tilde{h}_{i,k}]_{e} u_{e,k} l_e, +$$ + +using the TRiSK formulation of the discrete divergence operator (see +[V0 design document](omega-design-shallow-water-omega0)). 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 + +$$ +\tilde{h}^{n+1}_{i,k} = \tilde{h}^{n}_{i,k} + +\Delta t \left( \left[\tilde{W}_{tr}\right]_{k+1/2} - +\left[\tilde{W}_{tr}\right]_{k-1/2} \right), +$$ + +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`), + +$$ +u^{n+1}_{e,k} = u^{n}_{e,k} + +\frac{\Delta t}{\left[\tilde{h}_{i,k} \right]_{e}} \left( \left[U \tilde{W}_{tr} \right]_{e,k+1/2} - +\left[U \tilde{W}_{tr} \right]_{e,k-1/2} \right), +$$ + +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 + +$$ +\tilde{h}^{n+1}_{i,k} \tilde{\phi}^{n+1}_{i,k} = \tilde{h}^{n}_{i,k} \tilde{\phi}^{n}_{i,k} + +\Delta t \left( \left[ \phi \tilde{W}_{tr} \right]_{k+1/2} - +\left[ \phi \tilde{W}_{tr} \right]_{k-1/2} \right). +$$ + +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: + +```c++ +enum class VertFluxOption { + Second, /// 2nd-order centered + Third, /// 3rd-order upwind + Fourth /// 4th-order centered +}; +``` + +Another `enum class` specifies the advection scheme: + +```c++ +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. + +```c++ +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. + +```c++ +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. + +```c++ +VertAdv *VertAdv::DefVAdv; +std::map> 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. + +```c++ +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. + +```c++ +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. + +```c++ +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. + +```c++ +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. + +```c++ +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. diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index c7469d2ff19b..4fba7d6603df 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -123,6 +123,7 @@ design/Timers design/TimeStepping design/Tracers design/TridiagonalSolver +design/VertAdv design/VertCoord design/VerticalMixingCoeff From 62859de45befb3ef63c9a017065012245f708fa3 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 23 Oct 2025 12:15:18 -0400 Subject: [PATCH 2/2] Linting fix --- components/omega/doc/design/VertAdv.md | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/components/omega/doc/design/VertAdv.md b/components/omega/doc/design/VertAdv.md index 8e500d73b6bc..dc5d4aa6f1ca 100644 --- a/components/omega/doc/design/VertAdv.md +++ b/components/omega/doc/design/VertAdv.md @@ -17,7 +17,7 @@ scheme configuration. 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 +layer interfaces and tilt of layer interfaces ### 2.2 Requirement: Advection scheme flexibility @@ -95,7 +95,7 @@ The momentum equation governs the evolution of the horizontal velocity field (`NormalVelocity`), $$ -u^{n+1}_{e,k} = u^{n}_{e,k} + +u^{n+1}_{e,k} = u^{n}_{e,k} + \frac{\Delta t}{\left[\tilde{h}_{i,k} \right]_{e}} \left( \left[U \tilde{W}_{tr} \right]_{e,k+1/2} - \left[U \tilde{W}_{tr} \right]_{e,k-1/2} \right), $$ @@ -129,7 +129,7 @@ fluxes at layer interfaces: ```c++ enum class VertFluxOption { - Second, /// 2nd-order centered + Second, /// 2nd-order centered Third, /// 3rd-order upwind Fourth /// 4th-order centered }; @@ -220,7 +220,7 @@ class VertAdv{ // Arrays from HorzMesh Array2DI4 CellsOnEdge; Array2DI4 EdgesOnCell; - Array1DI4 NEdgesOnCell; + Array1DI4 NEdgesOnCell; Array1DReal AreaCell; Array1DReal DvEdge; Array2DReal EdgeSignOnCell; @@ -254,7 +254,7 @@ 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. +`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 @@ -274,11 +274,11 @@ void VertAdv::computeVerticalTransport(const OceanState *State, #### 4.2.3 Tendency computations -The `computeThicknessVAdvTend`, `computeVelocityVadvTend`, and +The `computeThicknessVAdvTend`, `computeVelocityVadvTend`, and `computeTracerVadvTend` methods are designed to be called by the `Tendencies` class methods. -The `computeThicknessVAdvTend` method computes the vertical advection tendency +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 @@ -327,8 +327,8 @@ void VertAdv::computeTracerVadvTend(const Array3DReal &Tend, 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 +`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