Add Vertical Advection design document

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<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.
+
+```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.

components/omega/doc/index.md

@@ -123,6 +123,7 @@ design/Timers
 design/TimeStepping
 design/Tracers
 design/TridiagonalSolver
+design/VertAdv
 design/VertCoord
 design/VerticalMixingCoeff