diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 111e4c59ee88..61fc232141fc 100644 --- a/components/omega/configs/Default.yml +++ b/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 diff --git a/components/omega/doc/design/RayleighDamping.md b/components/omega/doc/design/RayleighDamping.md new file mode 100644 index 000000000000..069b1ad8636e --- /dev/null +++ b/components/omega/doc/design/RayleighDamping.md @@ -0,0 +1,527 @@ +(omega-design-rayleigh-damping)= +# Rayleigh Damping + +**Table of Contents** +1. [Overview](#1-overview) +2. [Requirements](#2-requirements) +3. [Algorithmic Formulation](#3-algorithmic-formulation) +4. [Design](#4-design) +5. [Verification and Testing](#5-verification-and-testing) + +## 1 Overview + +An ocean initial condition assembled from observational climatologies is not in +dynamical balance with the model's discrete equations. When such a state is used +to start a simulation, the adjustment radiates large-amplitude barotropic and +inertia-gravity waves. These waves carry velocities far larger than the +climatological circulation, and they force the model to run with a very small +time step (or to fail outright on the CFL condition) until they have propagated +away or dissipated. + +The standard remedy, used routinely to produce E3SM ocean initial conditions +with MPAS-Ocean, is a *dynamic adjustment* spin-up: the model is run for a short +period with a strong linear drag on momentum, the so-called Rayleigh damping +term, which removes the transient wave energy on a prescribed time scale. The +spin-up proceeds as a sequence of adjustment segments, each run as a separate +job. Across the sequence the damping coefficient is reduced and the time step is +increased, until the damping is small enough to be turned off and the time step +has reached its production value. At that point the state is dynamically +balanced and can be integrated forward. + +The staging of the time step is as important as the staging of the coefficient. +Running the entire sequence at the small initial time step would take far longer +than necessary, and it would leave the stability of the production time step +untested until the adjustment was over. Each segment therefore serves double +duty: it damps transients, and it demonstrates that the state tolerates the +next, longer time step. + +Omega currently has no such capability, which is a blocker for generating E3SM +initial conditions (see +[Omega issue #495](https://github.com/E3SM-Project/Omega/issues/495)). This +document describes the addition of Rayleigh damping to Omega's momentum +equation, following the formulation used by MPAS-Ocean but expressed in Omega's +non-Boussinesq, pseudo-thickness-weighted form and configured with Omega naming +conventions. + +The scope here is deliberately narrow: a single, spatially uniform damping +coefficient applied to the normal velocity at all active layers, treated +implicitly in time. Two variants present in MPAS-Ocean are explicitly *not* +carried over: + +- `config_Rayleigh_damping_depth_variable`, which divides the coefficient by the + total column thickness. This is used by the subgrid wetting-and-drying work + rather than by dynamic adjustment, and its coefficient has different units + (m s$^{-1}$ rather than s$^{-1}$), which would make a single config option + ambiguous. +- `config_Rayleigh_bottom_damping_coeff`, which applies the damping only in the + bottom layer. That term is physically a linear bottom friction, not Rayleigh + damping, and belongs with Omega's bottom drag rather than here. If it is + needed later, it should be added as a new `Type` under the existing + `BottomDragTendency` config group. + +## 2 Requirements + +### 2.1 Requirement: Linear damping of momentum with a configurable time scale + +Omega shall provide a term that damps the normal velocity toward zero at a rate +set by a single, run-time configurable inverse time scale $c_R$ with units of +s$^{-1}$. + +Dynamic adjustment requires control over how fast transient wave energy is +removed. A damping time scale of a few hours ($c_R \sim 10^{-4}$ s$^{-1}$) is +short enough to suppress barotropic and near-inertial transients within a few +days of simulation, while a subsequent segment with a weaker coefficient allows +the slower baroclinic adjustment to proceed without being over-damped. + +### 2.2 Requirement: Damping is applied to all active layers of the momentum equation only + +The damping shall act on the normal velocity in every active layer between the +top and bottom of each edge column, and shall not act on pseudo-thickness or on +tracers. + +The transients being removed are momentum modes. Damping thickness or tracers +would relax the model toward a state that is not the intended initial condition +and would violate conservation of mass and tracer content. Applying the damping +over the whole column rather than only near the bottom removes barotropic and +baroclinic transients at comparable rates, which is the behavior needed for +adjustment. + +### 2.3 Requirement: Unconditional stability for large damping coefficients + +The discrete damping term shall be stable and monotone for arbitrarily large +values of $c_R \Delta t$, and shall not produce sign reversal or growth of the +velocity. + +The whole point of the term is to allow the time steps to be as *large* as +possible during adjustment. An explicit treatment would impose +$c_R \Delta t < 1$ for monotone decay, which couples the very quantity being +tuned ($c_R$) to the quantity being maximized ($\Delta t$) and defeats the +purpose. An implicit treatment removes that coupling entirely. + +### 2.4 Requirement: Disabled by default with no effect on existing results + +Rayleigh damping shall be off by default, and when off, model results shall be +bit-for-bit identical to results without the feature. + +Rayleigh damping is a spin-up tool, not part of the physics of a production +simulation. Leaving it enabled in a science run would artificially dissipate the +circulation. Bit-for-bit neutrality when disabled keeps existing test baselines +valid. + +### 2.5 Requirement: Configuration errors are detected at initialization + +If Rayleigh damping is enabled in a configuration where it cannot be applied, +Omega shall abort during initialization with an informative message rather than +silently ignoring the setting. + +A silently dropped damping term would produce an adjustment run that appears to +succeed but yields an unbalanced initial condition, and the failure would only +be discovered much later in the E3SM initialization workflow. + +### 2.6 Desired: Time-varying damping coefficient + +It is desirable for the damping coefficient to decrease smoothly over the course +of an adjustment segment, for example linearly or exponentially between a +starting and an ending value. + +The value of such a ramp is limited, and it is worth being precise about what it +would and would not buy. It would *not* replace the sequence of adjustment jobs: +as described in Section 1, each segment also raises the time step, and Omega has +no mechanism for changing the time step mid-run and no intention of adding one. +The sequence is therefore intrinsic to the workflow regardless of how $c_R$ is +specified. What a ramp would offer is the removal of the abrupt change in $c_R$ +at each segment boundary, replacing a step with a continuous transition. That is +a nice-to-have rather than a motivating capability, which is why it is listed as +desired rather than required. + +The design accommodates it at no cost: the coefficient is a plain scalar member +of the setup functor, so a ramp can later be implemented by updating that scalar +once per time step, with no change to the tridiagonal setup or to the config +group structure. + +## 3 Algorithmic Formulation + +### 3.1 Continuous form + +Rayleigh damping adds a linear drag to the momentum equation. In terms of the +normal velocity at an edge, $u_e$, the added tendency is + +$$ +\left. \frac{\partial u_e}{\partial t} \right|_{\mathrm{Rayleigh}} = -c_R u_e , +$$ + +where $c_R \ge 0$ has units of s$^{-1}$ and $1 / c_R$ is the e-folding time +scale of the damping. In isolation this gives pure exponential decay, +$u_e(t) = u_e(0) e^{-c_R t}$, at every wavenumber and every vertical mode. The +term is scale-selective in neither space nor time: it damps the mean circulation +along with the transients. That is acceptable for adjustment precisely because +the adjustment run is short and the transients are much larger in amplitude than +the circulation being preserved. + +The term is not energy conserving by construction. It removes kinetic energy at +a rate + +$$ +\left. \frac{\partial}{\partial t} \frac{u_e^2}{2} \right|_{\mathrm{Rayleigh}} += -c_R u_e^2 , +$$ + +and that energy is not returned to any other reservoir in the model. + +### 3.2 Layered, pseudo-thickness-weighted form + +Omega's discrete momentum equation is non-Boussinesq and is advanced in +pseudo-thickness-weighted form, with $\left[\tilde{h}_k\right]_e$ the +pseudo-thickness of layer $k$ interpolated to edge $e$. Grouping the Rayleigh +term with the vertical momentum stress divergence, the portion of the momentum +equation that is treated implicitly is + +$$ +\left[\tilde{h}_k\right]_e \frac{\partial u_{e,k}}{\partial t} += \left[ \tilde{\nu}_v \frac{\partial u}{\partial \tilde{z}} \right]_{e,k} +- \left[ \tilde{\nu}_v \frac{\partial u}{\partial \tilde{z}} \right]_{e,k+1} +- c_R \left[\tilde{h}_k\right]_e u_{e,k} . +$$ + +Because the damping rate is independent of thickness, the Rayleigh term picks up +the same $\left[\tilde{h}_k\right]_e$ factor as the time derivative. This is the +non-Boussinesq counterpart of the MPAS-Ocean form, in which the equivalent +system is divided through by the layer thickness and the term appears simply as +$c_R$ on the diagonal. + +### 3.3 Discrete implicit form + +Omega already solves the vertical momentum stress divergence implicitly with a +backward-Euler step, using the diffusion-form tridiagonal solver +`TriDiagDiffSolver` (see the [Tridiagonal Solver](TridiagonalSolver) design +document). The system solved at each edge for the updated normal velocity +$Y_k = u^{n+1}_{e,k}$ 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 , +$$ + +where $G_k$ holds the off-diagonal viscous coupling between layers $k$ and +$k+1$, $H_k$ holds contributions that appear only on the diagonal, and $X_k$ is +the right-hand side. Without damping, + +$$ +H_k = \left[\tilde{h}_k\right]_e ++ \delta_{k,k_{\mathrm{bot}}} \, \Delta t \, C_D \, \frac{\rho}{\rho_0} \left| \mathbf{u}_e \right| , +\qquad +X_k = \left[\tilde{h}_k\right]_e u^n_{e,k} , +$$ + +with the second term in $H_k$ the optional implicit bottom drag, present only in +the bottom layer. + +Backward-Euler treatment of the Rayleigh term contributes +$\Delta t \, c_R \left[\tilde{h}_k\right]_e u^{n+1}_{e,k}$ to the left-hand side, +so the *only* change to the system is an additional diagonal contribution in +every active layer: + +$$ +H_k = \left[\tilde{h}_k\right]_e \left( 1 + \Delta t \, c_R \right) ++ \delta_{k,k_{\mathrm{bot}}} \, \Delta t \, C_D \, \frac{\rho}{\rho_0} \left| \mathbf{u}_e \right| . +$$ + +Neither $G_k$ nor $X_k$ is modified, and the matrix keeps its symmetric, +diagonally dominant structure, so the existing solver applies unchanged. + +In the limit of no vertical viscosity and no bottom drag, the system decouples +and reduces to + +$$ +u^{n+1}_{e,k} = \frac{u^n_{e,k}}{1 + \Delta t \, c_R} , +$$ + +which is the analytic solution of the damping equation accurate to first order +in $\Delta t$, is strictly decreasing in magnitude, preserves sign, and is +bounded for any $\Delta t \, c_R \ge 0$. This satisfies the stability +requirement in Section 2.3, and is the relation the unit test in Section 5.1 +checks. It is also identical to the single-active-layer branch of MPAS-Ocean's +`ocn_vel_vmix_tend_implicit_rayleigh`, after accounting for the fact that +MPAS-Ocean divides its system by layer thickness. + +### 3.4 Typical coefficient values + +For reference, the coefficients used in MPAS-Ocean dynamic adjustment sequences +are of order $10^{-4}$ s$^{-1}$ (a damping time scale of about 2.8 hours) for +the first, most aggressive segment, decreasing by roughly an order of magnitude +over several segments before damping is switched off. Omega's default value +follows this convention, but the feature is disabled by default so the default +value is inert. + +## 4 Design + +The design adds no new class and no new field. Omega already assembles the +implicit velocity solve in one place, the `VelVertMixSetupOnEdge` functor in +`VertMix.h`, and that functor already carries an optional diagonal contribution +for implicit bottom drag. Rayleigh damping is a second such contribution, so it +is added the same way: two new public members on the functor, a few lines in +`operator()`, config reading in `VertMix::init()`, and a consistency check in +`Tendencies`. + +### 4.1 Data types and parameters + +#### 4.1.1 Parameters + +A new top-level `RayleighDamping` group is added to the Omega config, placed +after the `VertMix` group in `components/omega/configs/Default.yml`: + +```yaml + RayleighDamping: + Enable: false + DampingCoeff: 1.0e-4 +``` + +- `Enable` (logical, default `false`). Turns Rayleigh damping on. When `false`, + the damping term is not evaluated at all and results are bit-for-bit + unchanged, satisfying requirement 2.4. +- `DampingCoeff` (real, default `1.0e-4`, units s$^{-1}$). The coefficient + $c_R$. Any non-negative value is allowed; a negative value is a configuration + error and aborts at initialization. The default corresponds to a damping time + scale of about 2.8 hours, the typical value for the first segment of a dynamic + adjustment, but has no effect unless `Enable` is `true`. + +A top-level group is used rather than a subgroup of `Tendencies` because +Rayleigh damping is a self-contained spin-up capability rather than one of the +physical tendency terms, and because a top-level group is the natural place to +add the ramp parameters of requirement 2.6 later. + +Note that `DampingCoeff` is read only at initialization. Changing it therefore +requires a new job segment. This is not a limitation in practice: the time step +is also fixed for the duration of a run, so an adjustment sequence must be +staged as separate jobs in any case. + +Because the term is applied inside the implicit velocity vertical mixing solve, +enabling it requires the existing option + +```yaml + Tendencies: + VelVertMixTendencyEnable: true +``` + +which is the default. This is the same constraint that already applies to +`BottomDragTendency` with `Mode: Implicit`. + +#### 4.1.2 Class/structs/data types + +No new class is introduced. Two public members are added to +`VelVertMixSetupOnEdge` in `components/omega/src/ocn/VertMix.h`, alongside the +existing bottom drag members: + +```c++ +class VelVertMixSetupOnEdge { + public: + bool Enabled; + bool ImplicitBottomDragEnabled; + Real BottomDragCoeff; + bool RayleighDampingEnabled; ///< Enable Rayleigh damping flag + Real RayleighDampingCoeff; ///< Rayleigh damping coefficient (s^-1) + Real LocRhoSw; + ... +}; +``` + +Both are initialized in the constructor initializer list in `VertMix.cpp` to +`false` and `0.0_Real` respectively, so that a `VertMix` instance created +without reading a config performs no damping. + +### 4.2 Methods + +No new public method is added. Four existing pieces of code change. + +#### 4.2.1 `VelVertMixSetupOnEdge::operator()` + +The functor already computes the diagonal contribution `H` for layer `K` at edge +`IEdge`. The damping contribution is added immediately after the +pseudo-thickness is assigned to `H` and before the bottom drag contribution: + +```c++ + 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; + } +``` + +The functor is only invoked for layers between `KMin` and `KMax` at each edge; +layers outside that range are filled with the inactive-column values +(`G = 0`, `H = 1`, `X = 0`) by the caller in `VertMix::applyVelVertMixImplicit`. +The damping therefore automatically applies to exactly the active layers, +satisfying requirement 2.2, and needs no additional masking. The off-diagonal +term `G` and the right-hand side `X` are untouched. + +Guarding the addition with `RayleighDampingEnabled` rather than relying on a +zero coefficient mirrors the existing bottom drag pattern, avoids the extra +work in the common case, and makes the disabled path trivially bit-for-bit. + +#### 4.2.2 `VertMix::init()` + +`VertMix::init()` already reads the `VertMix` config group and its `Background`, +`Convective`, and `Shear` subgroups. It gains a block that reads the new +top-level group, following the same error-handling and logging idiom: + +```c++ + /// 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); + } +``` + +`VertMix::init()` is called from `OceanInit.cpp` before `Tendencies::init()`, +so the flag is set before the consistency check below runs. + +#### 4.2.3 `Tendencies` configuration check + +`Tendencies` reads `VelVertMixTendencyEnable` into +`VMix->VelVertMixSetup.Enabled` and already aborts if implicit bottom drag is +requested without it. An analogous check for Rayleigh damping is added +immediately after, satisfying requirement 2.5: + +```c++ + if (this->VMix->VelVertMixSetup.RayleighDampingEnabled && + !this->VMix->VelVertMixSetup.Enabled) { + ABORT_ERROR("Tendencies: RayleighDamping Enable requires " + "VelVertMixTendencyEnable to be true"); + } +``` + +#### 4.2.4 Documentation + +The `\ConfigInput` block in the header comment of `VertMix.h` is extended with +the new group, and the user guide gains a short section describing the option +and its use in dynamic adjustment. A new entry `design/RayleighDamping` is added +to the design document list in `components/omega/doc/index.md`. + +### 4.3 Items explicitly unaffected + +- **Halo exchanges.** The damping is a pointwise operation on each edge column + and introduces no new horizontal stencil, so no additional halo exchange is + needed. The implicit solve is already performed over all edges including + halos. +- **Restart and I/O.** No new prognostic or diagnostic field is created, so the + restart contents and all I/O streams are unchanged. A restart written during + an adjustment segment is read normally by a subsequent segment with a + different coefficient. +- **Explicit tendency path.** `Tendencies::computeVelocityTendencies` and the + functors in `TendencyTerms.h` are untouched. +- **Time steppers.** The implicit vertical mixing step is applied by the time + steppers through `VertMix::VertMixImplicit`, so all time steppers that call it + pick up the damping with no change. +- **Performance.** The added cost is one multiply-add per active edge layer + inside a kernel that is already memory-bound, and only when the feature is + enabled. + +## 5 Verification and Testing + +### 5.1 Test: single-layer analytic decay + +In `components/omega/test/ocn/VertMixTest.cpp`, add a test that configures a +column with vertical viscosity and bottom drag set to zero, Rayleigh damping +enabled, and a uniform nonzero initial normal velocity. After a single call to +`VertMix::applyVelVertMixImplicit`, the normal velocity in every active layer at +every edge must equal $u^n / (1 + \Delta t \, c_R)$ to within round-off of the +working precision. + +This is the exact analytic solution of the discrete system in the decoupled +limit, so the test verifies both the presence of the term and the correctness of +its pseudo-thickness weighting: an unweighted or doubly weighted term would give +a thickness-dependent answer and fail. + +List which requirements it tests: + - tests requirement 2.1 + - tests requirement 2.2 + +### 5.2 Test: disabled damping is bit-for-bit neutral + +Run the existing `VertMixTest` cases with `Enable: false` and confirm the +velocity and tracer results are bit-for-bit identical to the values expected +before this feature was added. Separately, confirm that with `Enable: true` and +`DampingCoeff: 0.0` the results are also bit-for-bit identical, which verifies +that the added diagonal contribution vanishes exactly rather than approximately. + +List which requirements it tests: + - tests requirement 2.4 + +### 5.3 Test: stability at large coefficient + +Repeat the configuration of Section 5.1 with $c_R \Delta t = 10^3$, far outside +the range an explicit scheme could tolerate. The resulting velocity must be +reduced in magnitude by the expected factor, must retain the sign of the initial +velocity, and must be finite. The test is repeated with both signs of the +initial velocity. + +List which requirements it tests: + - tests requirement 2.3 + +### 5.4 Test: combined with vertical viscosity and bottom drag + +Configure a multi-layer column with nonuniform pseudo-thickness, nonzero +vertical viscosity, implicit bottom drag enabled, and Rayleigh damping enabled. +Assemble the same tridiagonal system independently on the host with a +straightforward Thomas-algorithm implementation in the test, and require the +solver result to match to round-off. + +This verifies that the damping contribution is added to the diagonal only, that +it does not perturb the off-diagonal coupling or the right-hand side, and that +it composes correctly with the bottom drag term that shares the same diagonal. + +List which requirements it tests: + - tests requirement 2.1 + - tests requirement 2.2 + +### 5.5 Test: configuration error handling + +Confirm that a configuration with `RayleighDamping: Enable: true` and +`Tendencies: VelVertMixTendencyEnable: false` aborts during initialization with +the expected message, and that a negative `DampingCoeff` aborts during +`VertMix::init()`. + +List which requirements it tests: + - tests requirement 2.5 + +### 5.6 Test: dynamic adjustment in a realistic configuration + +As a system-level check, run a global Omega configuration started from a +climatological initial condition, first without damping and then with +`DampingCoeff: 1.0e-4`. The damped run must show the maximum normal velocity and +the globally averaged kinetic energy decaying over the first few simulated days, +and must remain stable at a time step at which the undamped run either fails or +requires a substantially smaller step. This is the end-to-end demonstration that +the feature serves its purpose, and is the natural place for a Polaris dynamic +adjustment task. + +List which requirements it tests: + - tests requirement 2.1 + - tests requirement 2.3 diff --git a/components/omega/doc/devGuide/RayleighDamping.md b/components/omega/doc/devGuide/RayleighDamping.md new file mode 100644 index 000000000000..9f738d71f7e7 --- /dev/null +++ b/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. diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 5e29df3996f4..0f94a46f4628 100644 --- a/components/omega/doc/index.md +++ b/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 diff --git a/components/omega/doc/userGuide/RayleighDamping.md b/components/omega/doc/userGuide/RayleighDamping.md new file mode 100644 index 000000000000..43a31d9d5cf2 --- /dev/null +++ b/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. diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 2ac8c6db4769..468051e5c9f6 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/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( diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index cb0a563e1fc8..25643b39ab7f 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/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) diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index 03c97e852132..26bb56b8c05d 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/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); diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index 7ed3601c068e..0aad5345662f 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -4,15 +4,19 @@ /// \brief Test driver for OMEGA Vertical Mixing Coefficients /// /// This driver tests that VertMix can be called and returns expected values -/// of diffusivity, viscosity and Brunt-Vaisala frequency +/// of diffusivity, viscosity and Brunt-Vaisala frequency, and that the +/// implicit velocity vertical mixing solve, including Rayleigh damping, +/// produces the expected velocities /// //===-----------------------------------------------------------------------===/ #include "VertMix.h" +#include "AuxiliaryState.h" #include "Config.h" #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Eos.h" #include "Field.h" #include "FillValues.h" #include "HorzMesh.h" @@ -21,12 +25,20 @@ #include "KPPMix.h" #include "Logging.h" #include "MachEnv.h" +#include "OceanState.h" #include "OceanTestCommon.h" #include "OmegaKokkos.h" #include "Pacer.h" +#include "TimeStepper.h" +#include "Tracers.h" +#include "VertAdv.h" #include "VertCoord.h" #include "mpi.h" +#include +#include +#include + using namespace OMEGA; /// Test constants and expected values @@ -79,6 +91,12 @@ void initVertMixTest() { Config("Omega"); Config::readAll("omega.yml"); + /// Initialize the default time stepper, which provides the calendar, the + /// model clock for stream IO and the time step used by the implicit + /// vertical mixing solve + TimeStepper::init1(); + Clock *ModelClock = TimeStepper::getDefault()->getClock(); + // Initialize the IO system IO::init(DefComm); @@ -88,13 +106,6 @@ void initVertMixTest() { /// Initialize Halo Halo::init(); - /// Create dummy model clock for stream IO - Calendar::init("No Leap"); - TimeInstant StartTime(0, 1, 1, 0, 0, 0.0); - TimeInterval TimeStep(1, TimeUnits::Hours); - Clock ModelClockTmp(StartTime, TimeStep); - Clock *ModelClock = &ModelClockTmp; - /// Initialize IO streams for mesh IO Field::init(ModelClock); IOStream::init(ModelClock); @@ -105,6 +116,13 @@ void initVertMixTest() { /// Initialize vertical coordinate VertCoord::init(false); + /// Initialize the tracers, auxiliary state and equation of state needed + /// by the implicit vertical mixing solve + Tracers::init(); + VertAdv::init(); + AuxiliaryState::init(); + Eos::init(); + /// Initialize VertMix VertMix::init(); @@ -1144,10 +1162,369 @@ void testInteriorSnapshotForMatchBoth() { } } +//------------------------------------------------------------------------------ +// Tests of Rayleigh damping in the implicit velocity vertical mixing solve + +/// Name of the ocean state used by the implicit vertical mixing tests +const std::string ImplicitStateName = "VertMixImplicitTest"; + +/// Test input values for the implicit vertical mixing tests +constexpr Real ImplicitNV = 0.3; // Normal velocity scale in m/s +constexpr Real ImplicitDrag = 1.0e-3; // Bottom drag coefficient +constexpr Real ImplicitViscSc = 1.0e-3; // Vertical viscosity scale in m^2/s + +/// Relative tolerance for comparing implicit solve results with expected +/// values, which differ only by round-off +const Real SolveRTol = 1000 * std::numeric_limits::epsilon(); + +/// Initial normal velocity, which alternates in sign between edges so that +/// both signs are tested +KOKKOS_INLINE_FUNCTION Real initNormalVel(I4 IEdge, I4 K) { + const Real Sign = (IEdge % 2 == 0) ? 1.0_Real : -1.0_Real; + return Sign * ImplicitNV * (1.0_Real + 0.01_Real * K); +} + +/// Return the time step (s) used by VertMix::applyVelVertMixImplicit +R8 getImplicitTimeStep() { + R8 DT; + TimeStepper::getDefault()->getTimeStep().get(DT, TimeUnits::Seconds); + return DT; +} + +/// Set the state, auxiliary state, specific volume and viscosity used by the +/// implicit velocity vertical mixing solve. Pseudo-thickness is nonuniform in +/// both the horizontal and the vertical. If UseVisc is false, the vertical +/// viscosity is zero so that the layers decouple. +OceanState *initImplicitTestState(bool UseVisc) { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + auto *MeshHalo = Halo::getDefault(); + auto *AuxState = AuxiliaryState::getDefault(); + VertMix *VMix = VertMix::getInstance(); + Eos *EqState = Eos::getInstance(); + const I4 NLayers = VCoord->NVertLayers; + + OceanState *State = OceanState::get(ImplicitStateName); + if (!State) { + State = OceanState::create(ImplicitStateName, Mesh, MeshHalo, NLayers, 2); + if (!State) + ABORT_ERROR("TestVertMixImplicit: error creating test state"); + } + + Array2DReal PseudoThick = State->getPseudoThickness(0); + Array2DReal NormalVel = State->getNormalVelocity(0); + Array2DReal SpecVol = EqState->SpecVol; + Array2DReal VertVisc = VMix->VertVisc; + Array2DReal KECell = AuxState->KineticAux.KineticEnergyCell; + + parallelFor( + "initImplicitCells", {Mesh->NCellsAll, NLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + PseudoThick(ICell, K) = + 1.0_Real + 0.5_Real * K + 0.1_Real * (ICell % 7); + SpecVol(ICell, K) = (1.0_Real + 1.0e-3_Real * (K % 3)) / RhoSw; + KECell(ICell, K) = 0.01_Real * (1 + (ICell + K) % 3); + }); + + parallelFor( + "initImplicitEdges", {Mesh->NEdgesAll, NLayers}, + KOKKOS_LAMBDA(I4 IEdge, I4 K) { + NormalVel(IEdge, K) = initNormalVel(IEdge, K); + }); + + if (UseVisc) { + parallelFor( + "initImplicitVisc", {Mesh->NCellsAll, NLayers + 1}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + VertVisc(ICell, K) = ImplicitViscSc * (1 + (ICell + K) % 5); + }); + } else { + deepCopy(VertVisc, 0.0_Real); + } + + return State; +} + +/// Return a host copy of an array. Unlike createHostMirrorCopy, the result +/// never aliases the array, even when host and device memory are the same. +HostArray2DReal copyToHost(const Array2DReal &Arr) { + const auto ArrMirror = createHostMirrorCopy(Arr); + HostArray2DReal ArrH("ArrH", Arr.extent(0), Arr.extent(1)); + deepCopy(ArrH, ArrMirror); + return ArrH; +} + +/// Run the implicit velocity vertical mixing solve with the given options +/// and return a host copy of the resulting normal velocity +HostArray2DReal runImplicitVelVertMix(bool UseVisc, bool BottomDragEnabled, + bool RayleighEnabled, + Real RayleighCoeff) { + VertMix *VMix = VertMix::getInstance(); + auto *State = initImplicitTestState(UseVisc); + + VMix->VelVertMixSetup.Enabled = true; + VMix->VelVertMixSetup.ImplicitBottomDragEnabled = BottomDragEnabled; + VMix->VelVertMixSetup.BottomDragCoeff = ImplicitDrag; + VMix->VelVertMixSetup.RayleighDampingEnabled = RayleighEnabled; + VMix->VelVertMixSetup.RayleighDampingCoeff = RayleighCoeff; + + VMix->applyVelVertMixImplicit(State, AuxiliaryState::getDefault(), 0, 0); + + // Restore the default (disabled) state of the optional terms + VMix->VelVertMixSetup.ImplicitBottomDragEnabled = false; + VMix->VelVertMixSetup.BottomDragCoeff = 0.0_Real; + VMix->VelVertMixSetup.RayleighDampingEnabled = false; + VMix->VelVertMixSetup.RayleighDampingCoeff = 0.0_Real; + + return copyToHost(State->getNormalVelocity(0)); +} + +/// Independently assemble the tridiagonal system for the implicit velocity +/// solve on the host and solve it with the Thomas algorithm, starting from +/// the initial normal velocity. Assumes that initImplicitTestState has been +/// called with the same UseVisc. +HostArray2DReal computeReferenceVel(bool BottomDragEnabled, R8 RayleighCoeff) { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + const I4 NLayers = VCoord->NVertLayers; + const R8 DT = getImplicitTimeStep(); + auto *AuxState = AuxiliaryState::getDefault(); + const auto *State = OceanState::get(ImplicitStateName); + const auto HH = createHostMirrorCopy(State->getPseudoThickness(0)); + const auto SpecVH = createHostMirrorCopy(Eos::getInstance()->SpecVol); + const auto ViscH = createHostMirrorCopy(VertMix::getInstance()->VertVisc); + const auto KEH = + createHostMirrorCopy(AuxState->KineticAux.KineticEnergyCell); + const auto CellsOnEdgeH = createHostMirrorCopy(Mesh->CellsOnEdge); + const auto KMinH = createHostMirrorCopy(VCoord->MinLayerEdgeBot); + const auto KMaxH = createHostMirrorCopy(VCoord->MaxLayerEdgeTop); + + // The state's normal velocity has already been updated by the implicit + // solve, so the initial velocity is recomputed here instead + HostArray2DReal RefVel("RefVel", Mesh->NEdgesSize, NLayers); + for (I4 IEdge = 0; IEdge < Mesh->NEdgesOwned; ++IEdge) { + for (I4 K = 0; K < NLayers; ++K) { + RefVel(IEdge, K) = initNormalVel(IEdge, K); + } + } + + std::vector Lower(NLayers), Diag(NLayers), Upper(NLayers), Rhs(NLayers); + + for (I4 IEdge = 0; IEdge < Mesh->NEdgesOwned; ++IEdge) { + const I4 C0 = CellsOnEdgeH(IEdge, 0); + const I4 C1 = CellsOnEdgeH(IEdge, 1); + const I4 KMin = KMinH(IEdge); + const I4 KMax = KMaxH(IEdge); + + // Skip columns with no active layers + if (KMax < KMin) + continue; + + auto edgeAvg = [&](const auto &Arr, I4 K) { + return 0.5 * (R8(Arr(C0, K)) + R8(Arr(C1, K))); + }; + + // Coupling between layers K and K+1 through the interface between them + std::vector Coupling(NLayers, 0.0); + for (I4 K = KMin; K < KMax; ++K) { + const R8 HK = edgeAvg(HH, K); + const R8 HKp1 = edgeAvg(HH, K + 1); + const R8 HBot = 0.5 * (HK + HKp1); + const R8 AlpB = + (edgeAvg(SpecVH, K) * HKp1 + edgeAvg(SpecVH, K + 1) * HK) / + (HK + HKp1); + const R8 Visc = edgeAvg(ViscH, K + 1); + Coupling[K] = DT * Visc / (R8(RhoSw) * AlpB) / HBot; + } + + for (I4 K = KMin; K <= KMax; ++K) { + const R8 HK = edgeAvg(HH, K); + R8 DiagOnly = HK * (1.0 + DT * RayleighCoeff); + if (BottomDragEnabled && K == KMax) { + const R8 VelMag = std::sqrt(R8(KEH(C0, K)) + R8(KEH(C1, K))); + DiagOnly += DT * R8(ImplicitDrag) * VelMag / + (R8(RhoSw) * edgeAvg(SpecVH, K)); + } + const R8 GAbove = (K > KMin) ? Coupling[K - 1] : 0.0; + const R8 GBelow = (K < KMax) ? Coupling[K] : 0.0; + Lower[K] = -GAbove; + Upper[K] = -GBelow; + Diag[K] = GAbove + GBelow + DiagOnly; + Rhs[K] = HK * R8(initNormalVel(IEdge, K)); + } + + // Thomas algorithm: forward elimination and back substitution + for (I4 K = KMin + 1; K <= KMax; ++K) { + const R8 W = Lower[K] / Diag[K - 1]; + Diag[K] -= W * Upper[K - 1]; + Rhs[K] -= W * Rhs[K - 1]; + } + R8 Below = Rhs[KMax] / Diag[KMax]; + RefVel(IEdge, KMax) = Below; + for (I4 K = KMax - 1; K >= KMin; --K) { + Below = (Rhs[K] - Upper[K] * Below) / Diag[K]; + RefVel(IEdge, K) = Below; + } + } + + return RefVel; +} + +/// Count mismatches between the implicit solve result and the expected +/// velocity on owned edges. Active layers must match to within SolveRTol, +/// and inactive layers must be unchanged from the initial velocity. +int countImplicitMismatches(const HostArray2DReal &Vel, + const HostArray2DReal &ExpectedVel) { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + const auto KMinH = createHostMirrorCopy(VCoord->MinLayerEdgeBot); + const auto KMaxH = createHostMirrorCopy(VCoord->MaxLayerEdgeTop); + + int NumMismatches = 0; + for (I4 IEdge = 0; IEdge < Mesh->NEdgesOwned; ++IEdge) { + for (I4 K = 0; K < VCoord->NVertLayers; ++K) { + if (K < KMinH(IEdge) || K > KMaxH(IEdge)) { + if (Vel(IEdge, K) != initNormalVel(IEdge, K)) + NumMismatches++; + } else if (!std::isfinite(Vel(IEdge, K)) || + !isApprox(Vel(IEdge, K), ExpectedVel(IEdge, K), + SolveRTol)) { + NumMismatches++; + } + } + } + return NumMismatches; +} + +/// With no vertical viscosity and no bottom drag, the layers decouple and +/// the implicit damping gives u^{n+1} = u^n / (1 + DT * Coeff) in every +/// active layer. The analytic factor is independent of pseudo-thickness, so +/// this also checks the pseudo-thickness weighting of the damping term. For +/// large DT * Coeff, this also checks that the decay is stable, monotone, +/// and preserves the sign of the velocity. +void testRayleighDampingDecay(const std::string &TestName, Real Coeff) { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + const R8 DT = getImplicitTimeStep(); + const R8 Factor = 1.0 / (1.0 + DT * Coeff); + + HostArray2DReal Vel = runImplicitVelVertMix(false, false, true, Coeff); + + HostArray2DReal ExpectedVel("ExpectedVel", Vel.extent(0), Vel.extent(1)); + int NumSignErrors = 0; + for (I4 IEdge = 0; IEdge < Mesh->NEdgesOwned; ++IEdge) { + for (I4 K = 0; K < VCoord->NVertLayers; ++K) { + const Real U0 = initNormalVel(IEdge, K); + ExpectedVel(IEdge, K) = Factor * U0; + if (Vel(IEdge, K) * U0 <= 0.0_Real || + std::abs(Vel(IEdge, K)) > std::abs(U0)) + NumSignErrors++; + } + } + + const int NumMismatches = countImplicitMismatches(Vel, ExpectedVel); + + if (NumMismatches != 0 || NumSignErrors != 0) { + ABORT_ERROR("TestVertMixRayleigh: {} FAIL with DT * Coeff = {}, " + "{} mismatches and {} sign or growth errors", + TestName, DT * Coeff, NumMismatches, NumSignErrors); + } else { + LOG_INFO("TestVertMixRayleigh: {} PASS", TestName); + } +} + +/// Disabled Rayleigh damping, and enabled damping with a zero coefficient, +/// must give bit-for-bit identical results, and both must match the system +/// without damping +void testRayleighDampingNeutral() { + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + HostArray2DReal DisabledVel = + runImplicitVelVertMix(true, true, false, 1.0e-4_Real); + HostArray2DReal RefVel = computeReferenceVel(true, 0.0); + HostArray2DReal ZeroCoeffVel = + runImplicitVelVertMix(true, true, true, 0.0_Real); + + int NumDiffs = 0; + for (I4 IEdge = 0; IEdge < Mesh->NEdgesOwned; ++IEdge) { + for (I4 K = 0; K < VCoord->NVertLayers; ++K) { + if (DisabledVel(IEdge, K) != ZeroCoeffVel(IEdge, K)) + NumDiffs++; + } + } + + const int NumMismatches = countImplicitMismatches(DisabledVel, RefVel); + + if (NumDiffs != 0 || NumMismatches != 0) { + ABORT_ERROR("TestVertMixRayleigh: Neutral FAIL with {} differences " + "from zero coefficient and {} mismatches from reference", + NumDiffs, NumMismatches); + } else { + LOG_INFO("TestVertMixRayleigh: Neutral PASS"); + } +} + +/// Rayleigh damping combined with nonuniform vertical viscosity and implicit +/// bottom drag must match an independent host solve of the same system +void testRayleighDampingCombined() { + const Real Coeff = 1.0e-3_Real; + + HostArray2DReal Vel = runImplicitVelVertMix(true, true, true, Coeff); + HostArray2DReal RefVel = computeReferenceVel(true, Coeff); + + // Also make sure the damping changed the result + HostArray2DReal UndampedVel = + runImplicitVelVertMix(true, true, false, Coeff); + + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + const auto KMinH = createHostMirrorCopy(VCoord->MinLayerEdgeBot); + const auto KMaxH = createHostMirrorCopy(VCoord->MaxLayerEdgeTop); + int NumUnchanged = 0; + for (I4 IEdge = 0; IEdge < Mesh->NEdgesOwned; ++IEdge) { + for (I4 K = KMinH(IEdge); K <= KMaxH(IEdge); ++K) { + if (!(std::abs(Vel(IEdge, K)) < std::abs(UndampedVel(IEdge, K)))) + NumUnchanged++; + } + } + + const int NumMismatches = countImplicitMismatches(Vel, RefVel); + + if (NumMismatches != 0 || NumUnchanged != 0) { + ABORT_ERROR("TestVertMixRayleigh: Combined FAIL with {} mismatches " + "from reference and {} layers not damped", + NumMismatches, NumUnchanged); + } else { + LOG_INFO("TestVertMixRayleigh: Combined PASS"); + } +} + +/// Run all of the Rayleigh damping tests +void testRayleighDamping() { + const R8 DT = getImplicitTimeStep(); + + // Typical coefficient for the first segment of a dynamic adjustment + testRayleighDampingDecay("Decay", 1.0e-4_Real); + + // DT * Coeff = 1000, far outside the range an explicit scheme could + // tolerate + testRayleighDampingDecay("LargeCoeff", Real(1.0e3 / DT)); + + testRayleighDampingNeutral(); + testRayleighDampingCombined(); +} + /// Finalize and clean up all test infrastructure void finalizeVertMixTest() { KPPMix::destroyInstance(); VertMix::destroyInstance(); + Eos::destroyInstance(); + OceanState::clear(); + AuxiliaryState::clear(); + VertAdv::clear(); + Tracers::clear(); + TimeStepper::clear(); HorzMesh::clear(); Halo::clear(); VertCoord::clear(); @@ -1168,6 +1545,8 @@ void finalizeVertMixTest() { // for only shear // --> next tests the linear superposition of the // background, convective, and shear contributions +// --> finally tests Rayleigh damping in the implicit +// velocity vertical mixing solve void vertMixTest() { // initialize vertical mix and other infrastructure @@ -1182,6 +1561,9 @@ void vertMixTest() { testTotalVertMix(); testInteriorSnapshotForMatchBoth(); + // test Rayleigh damping in the implicit velocity vertical mixing solve + testRayleighDamping(); + // clean up finalizeVertMixTest(); }