Add Rayleigh damping

components/omega/configs/Default.yml

@@ -145,6 +145,9 @@ Omega:
       IceFractionThresholdForMinimumOSBL: 0.15
       MinimumOSBLUnderSeaIce: 5.0
       DebugDiagnostics: false
+  RayleighDamping:
+    Enable: false
+    DampingCoeff: 1.0e-4
   IOStreams:
     HorzMeshIn:
       UsePointerFile: false

components/omega/doc/devGuide/RayleighDamping.md

@@ -0,0 +1,43 @@
+(omega-dev-rayleigh-damping)=
+
+# Rayleigh Damping
+
+Rayleigh damping (see the [user guide](#omega-user-rayleigh-damping)) is not a separate class or tendency term. It is a diagonal contribution to the tridiagonal system solved by the implicit velocity vertical mixing in `VertMix::applyVelVertMixImplicit`.
+
+## Implementation
+
+The `VelVertMixSetupOnEdge` functor in `VertMix.h`, which sets up the tridiagonal system for each edge column, has two public members for Rayleigh damping, alongside the members for implicit bottom drag:
+
+```c++
+bool RayleighDampingEnabled; ///< Enable Rayleigh damping flag
+Real RayleighDampingCoeff;   ///< Rayleigh damping coefficient (s^-1)
+```
+
+The constructor sets them to `false` and `0`. They are read from the `RayleighDamping` config group in `VertMix::init()`, which aborts if `DampingCoeff` is negative.
+
+The tridiagonal system solved for the new normal velocity $Y_k$ in each active layer is
+
+$$
+-G_{k-1} Y_{k-1} + \left( G_{k-1} + G_k + H_k \right) Y_k - G_k Y_{k+1} = X_k .
+$$
+
+The backward-Euler Rayleigh term adds $\Delta t \, c_R \, \tilde{h}_k$ to the diagonal-only coefficient $H_k$, where $\tilde{h}_k$ is the pseudo-thickness of layer $k$ interpolated to the edge. With bottom drag, $H_k$ becomes
+
+$$
+H_k = \tilde{h}_k \left( 1 + \Delta t \, c_R \right)
++ \delta_{k,k_{\mathrm{bot}}} \, \Delta t \, C_D \, \frac{\rho}{\rho_0} \left| \mathbf{u} \right| .
+$$
+
+The off-diagonal coefficients $G_k$ and the right-hand side $X_k = \tilde{h}_k u^n_k$ are unchanged. The functor is only called for active layers, so damping is applied to exactly those layers. The contribution is guarded by `RayleighDampingEnabled`, so the disabled path is unchanged.
+
+Because the damping is part of the implicit solve, `Tendencies` aborts at initialization if Rayleigh damping is enabled but `VelVertMixTendencyEnable` is `false`, as it already does for implicit bottom drag. `VertMix::init()` must therefore be called before `Tendencies::init()`.
+
+All time steppers apply the damping through `VertMix::VertMixImplicit`. No fields, halo exchanges or I/O are added.
+
+## Testing
+
+`VertMixTest` calls `VertMix::applyVelVertMixImplicit` with a test `OceanState` and checks:
+
+- With no vertical viscosity or bottom drag, every active layer decays by exactly $1/(1 + \Delta t \, c_R)$, both for a typical coefficient and for $c_R \Delta t = 10^3$. Velocity must keep its sign and not grow. Inactive layers must be unchanged.
+- With viscosity and bottom drag, disabled damping and enabled damping with $c_R = 0$ give bit-for-bit identical results.
+- With nonuniform viscosity, implicit bottom drag and damping all enabled, the result matches a host Thomas-algorithm solve of the system above.

components/omega/doc/index.md

@@ -55,6 +55,7 @@ userGuide/PGrad
 userGuide/Timing
 userGuide/VerticalMixingCoeff
 userGuide/KPPMix
+userGuide/RayleighDamping
 userGuide/SubmesoEddies
 userGuide/VertAdv
 userGuide/Forcing
@@ -109,6 +110,7 @@ devGuide/PGrad
 devGuide/Timing
 devGuide/VerticalMixingCoeff
 devGuide/KPPMix
+devGuide/RayleighDamping
 devGuide/SubmesoEddies
 devGuide/VertAdv
 devGuide/Forcing
@@ -142,6 +144,7 @@ design/Metadata
 design/PGrad
 design/IO
 design/IOStreams
+design/RayleighDamping
 design/Reductions
 design/State
 design/StateValidation

components/omega/doc/userGuide/RayleighDamping.md

@@ -0,0 +1,34 @@
+(omega-user-rayleigh-damping)=
+
+# Rayleigh Damping
+
+Rayleigh damping adds a linear drag on the normal velocity to the momentum equation:
+
+$$
+\left. \frac{\partial u}{\partial t} \right|_{\mathrm{Rayleigh}} = -c_R u ,
+$$
+
+where $c_R$ (s$^{-1}$) is the damping coefficient and $1/c_R$ is the e-folding time scale of the damping. The damping acts on every active layer of every edge, but not on pseudo-thickness or tracers.
+
+Rayleigh damping is a spin-up tool, not part of the physics of a production simulation. It is used for *dynamic adjustment*: a short run, starting from an initial condition that is not in dynamical balance (e.g. one built from observational climatologies), in which strong damping removes the large barotropic and inertia-gravity waves generated as the model adjusts. A dynamic adjustment is typically run as a sequence of jobs (restarting from the previous one) in which the damping coefficient is reduced and the time step increased from one job to the next until damping can be turned off and the time step has reached its production value. Coefficients of order $10^{-4}$ s$^{-1}$ (a time scale of about 2.8 hours) are typical for the first, most aggressive job.
+
+Rayleigh damping is configured with the `RayleighDamping` group in the yaml configuration file:
+
+```yaml
+RayleighDamping:
+  Enable: false        # Enables Rayleigh damping
+  DampingCoeff: 1.0e-4 # Damping coefficient (s^-1)
+```
+
+Rayleigh damping is disabled by default, in which case `DampingCoeff` has no effect and results are bit-for-bit identical to results without the feature. `DampingCoeff` must be non-negative and is read only at initialization, so changing it requires a new job.
+
+The damping term is treated implicitly in time as part of the implicit vertical mixing of velocity, so it is stable for any value of $c_R \Delta t$. Rayleigh damping therefore requires implicit velocity vertical mixing to be enabled:
+
+```yaml
+Tendencies:
+  VelVertMixTendencyEnable: true
+```
+
+which is the default. Omega aborts during initialization if Rayleigh damping is enabled and `VelVertMixTendencyEnable` is `false`.
+
+See the [Rayleigh damping design document](#omega-design-rayleigh-damping) for the details of the formulation.

components/omega/src/ocn/Tendencies.cpp

@@ -474,6 +474,12 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config
                   "VelVertMixTendencyEnable to be true");
    }
 
+   if (this->VMix->VelVertMixSetup.RayleighDampingEnabled &&
+       !this->VMix->VelVertMixSetup.Enabled) {
+      ABORT_ERROR("Tendencies: RayleighDamping Enable requires "
+                  "VelVertMixTendencyEnable to be true");
+   }
+
    Err += TendConfig.get("TracerVertMixTendencyEnable",
                          this->TracerVertMixTendencyEnable);
    CHECK_ERROR_ABORT(

components/omega/src/ocn/VertMix.cpp

@@ -46,7 +46,8 @@ OneTwoOneFilter::OneTwoOneFilter(const VertCoord *VCoord)
 VelVertMixSetupOnEdge::VelVertMixSetupOnEdge(const HorzMesh *Mesh,
                                              const VertCoord *VCoord)
     : Enabled(false), ImplicitBottomDragEnabled(false),
-      BottomDragCoeff(0.0_Real), LocRhoSw(RhoSw),
+      BottomDragCoeff(0.0_Real), RayleighDampingEnabled(false),
+      RayleighDampingCoeff(0.0_Real), LocRhoSw(RhoSw),
       NVertLayers(VCoord->NVertLayers), CellsOnEdge(Mesh->CellsOnEdge),
       EdgeMask(VCoord->EdgeMask), MinLayerEdgeBot(VCoord->MinLayerEdgeBot),
       MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {}
@@ -208,6 +209,37 @@ void VertMix::init() {
       CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter Shear:RiSmoothLoops not "
                              "found in ShearConfig");
    }
+
+   /// Get RayleighDamping group from Omega config
+   Config RayleighConfig("RayleighDamping");
+   Err += OmegaConfig->get(RayleighConfig);
+   CHECK_ERROR_ABORT(
+       Err, "VertMix::init: RayleighDamping group not found in Config");
+
+   Err += RayleighConfig.get(
+       "Enable", DefVertMix->VelVertMixSetup.RayleighDampingEnabled);
+   CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter RayleighDamping:Enable not "
+                          "found in RayleighConfig");
+
+   if (!DefVertMix->VelVertMixSetup.RayleighDampingEnabled) {
+      LOG_INFO("VertMix::init: Rayleigh damping is disabled.");
+   } else {
+      Err += RayleighConfig.get(
+          "DampingCoeff", DefVertMix->VelVertMixSetup.RayleighDampingCoeff);
+      CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter "
+                             "RayleighDamping:DampingCoeff not found in "
+                             "RayleighConfig");
+
+      if (DefVertMix->VelVertMixSetup.RayleighDampingCoeff < 0.0_Real) {
+         ABORT_ERROR("VertMix::init: RayleighDamping:DampingCoeff must be "
+                     "non-negative but got {}",
+                     DefVertMix->VelVertMixSetup.RayleighDampingCoeff);
+      }
+
+      LOG_INFO("VertMix::init: Rayleigh damping is enabled with coefficient "
+               "{} s^-1.",
+               DefVertMix->VelVertMixSetup.RayleighDampingCoeff);
+   }
 } // end init
 
 /// Compute diffusivity and viscosity for all cells/layers (no displacement)

components/omega/src/ocn/VertMix.h

@@ -234,6 +234,8 @@ class VelVertMixSetupOnEdge {
    bool Enabled;
    bool ImplicitBottomDragEnabled;
    Real BottomDragCoeff;
+   bool RayleighDampingEnabled; ///< Enable Rayleigh damping flag
+   Real RayleighDampingCoeff;   ///< Rayleigh damping coefficient (s^-1)
    Real LocRhoSw;
 
    VelVertMixSetupOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord);
@@ -274,7 +276,8 @@ class VelVertMixSetupOnEdge {
       //
       //   G_k = [DT * VertVisc_k^top / (Rho_0 * SpecVol_k^top)]
       //         / PseudoThick_k^top
-      //   H_k = PseudoThick_k, plus implicit bottom-drag contribution
+      //   H_k = PseudoThick_k * (1 + DT * RayleighDampingCoeff), plus
+      //         implicit bottom-drag contribution
       //   X_k = PseudoThick_k * NormVel_k
 
       // Fill values
@@ -290,6 +293,13 @@ class VelVertMixSetupOnEdge {
 
       H = PseudoThickEdgeK;
 
+      // Rayleigh damping: a linear drag on momentum applied in every active
+      // layer. Treated implicitly, so it is unconditionally stable for any
+      // DT * RayleighDampingCoeff.
+      if (RayleighDampingEnabled) {
+         H += DT * RayleighDampingCoeff * PseudoThickEdgeK;
+      }
+
       // Unknown is NormVel^{n+1}.
       X = PseudoThickEdgeK * NormalVelEdge(IEdge, K);