From 0752c95b31684ca249fd0c02dea565d555623fd9 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 2 Feb 2026 16:37:19 -0800 Subject: [PATCH 1/4] changes shear mixing to LMD94 formulation and adds ri smoothing --- components/omega/configs/Default.yml | 7 +- components/omega/doc/devGuide/EOS.md | 2 +- .../omega/doc/devGuide/VerticalMixingCoeff.md | 29 +- .../doc/userGuide/VerticalMixingCoeff.md | 31 +- components/omega/src/ocn/VertMix.cpp | 192 +++++- components/omega/src/ocn/VertMix.h | 200 ++++-- components/omega/test/ocn/VertMixTest.cpp | 584 ++++++++++++++---- 7 files changed, 817 insertions(+), 228 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 93b4fe7f7bf3..be6b649820df 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -69,9 +69,10 @@ Omega: TriggerBVF: 0.0 Shear: Enable: true - NuZero: 0.005 - Alpha: 5.0 - Exponent: 2.0 + BaseShearValue: 0.005 + RiCrit: 0.7 + Exponent: 3.0 + RiSmoothLoops: 3 IOStreams: InitialVertCoord: UsePointerFile: false diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 0762727a25b5..9aeaec4fff23 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -7,7 +7,7 @@ and `BruntVaisalaFreqSq`. Current EOS options are a linear EOS or an EOS compute 75 term expansion from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566). If `SpecVolDisplaced` is calculated with the linear EOS option, it will be equal to `SpecVol` as there is no pressure/depth dependence for the linear EOS. `SpecVolDisplaced` computes specific volume -adiabatically displaced to `K + KDisp` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreqSq`, as +adiabatically displaced to `K + KDisp` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreq`, as `SpecVol` is an input for the `BruntVaisalaFreqSq` calculation. If the linear EOS option is used, then the `BruntVaisalaFreqSq` is calculated using linear coefficients. If the TEOS-10 option is used, the `BruntVaisalaFreqSq` is calculated with non-linear coefficients according to the [TEOS-10 toolbox](https://www.teos-10.org/software.htm). Note: two assumption for ease of computation and efficiency have been made diff --git a/components/omega/doc/devGuide/VerticalMixingCoeff.md b/components/omega/doc/devGuide/VerticalMixingCoeff.md index 2b3fa0a86fa2..79736261bc72 100644 --- a/components/omega/doc/devGuide/VerticalMixingCoeff.md +++ b/components/omega/doc/devGuide/VerticalMixingCoeff.md @@ -2,12 +2,12 @@ # Vertical Mixing Coefficients -Omega includes a `VertMix` class that provides functions that compute `VertDiff` and `VertVisc`, the -vertical diffusivity and viscosity, where both are defined at the center of the cell and top of the layer. +Omega includes a `VertMix` class that provides functions that compute `VertDiff`, `VertVisc`, `GradRichNum`, and `GradRichNumSmoothed`, the +vertical diffusivity and viscosity, the gradient Richardson number, a smoothed gradient Richardson number, where all are defined at the center of the cell and top of the layer. Currently the values of `VertDiff` and `VertVisc` are calculated using the linear combination of three options: (1) a constant background mixing value, (2) a convective instability mixing value, and (3) a Richardson -number dependent shear mixing value from the [Pacanowski and Philander (1981)](https://journals.ametsoc.org/view/journals/phoc/11/11/1520-0485_1981_011_1443_povmin_2_0_co_2.xml) parameterization. These options are linearly additive. In the future, additional additive options will be implemented, such as the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). For both the convective and shear mixing values `BruntVaisalaFreqSq` is needed, which -is calculated by the `EOS` class. +number dependent shear mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These options are linearly additive. In the future, additional additive options will be implemented, such as the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). For both the convective and shear mixing values `BruntVaisalaFreq` is needed, which +is calculated by the `EOS` class. `GradRichNum` is smoothed using a 1-2-1 filter to produce `GradRichNumSmoothed` which is used by the LMD94 shear mixing formulation. ## Initialization and Usage @@ -37,9 +37,10 @@ VertMix: TriggerBVF: 0.0 Shear: Enable: true - NuZero: 0.005 - Alpha: 5.0 - Exponent: 2.0 + BaseShearValue: 0.005 + RiCrit: 0.7 + Exponent: 3.0 + RiSmoothLoops: 3 ``` ## Class Structure @@ -48,6 +49,8 @@ VertMix: - `VertDiff`: 2D array storing vertical diffusivity coefficients (m²/s) - `VertVisc`: 2D array storing vertical viscosity coefficients (m²/s) +- `GradRichNum`: 2D array storing the gradient Richardson number +- `GradRichNumSmoothed`: 2D array storing the smoothed gradient Richardson number ### Mixing Parameters @@ -62,9 +65,10 @@ VertMix: 3. Shear Mixing: - `EnableShearMix`: Flag to enable/disable shear mixing (Default: True) - - `ShearNuZero`: Base coefficient for Pacanowski-Philander scheme (Default: 0.005) - - `ShearAlpha`: Alpha parameter for P-P scheme (Default: 5.0) - - `ShearExponent`: Exponent parameter for P-P scheme (Default: 2.0) + - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) + - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) + - `ShearExponent`: Exponent parameter number for the LMB94 formulation (Default: 3.0) + - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 3) ## Core Functionality (Vertical Mixing Coefficient Calculation) @@ -75,8 +79,3 @@ void computeVertMix(const Array2DReal &NormalVelocity, const Array2DReal &TangentialVelocity, const Array2DReal &BruntVaisalaFreqSq); ``` - -This method combines the effects of: -- Background mixing (constant coefficients) -- Convective mixing (triggered by static instability) -- Shear instability driven mixing (Pacanowski-Philander scheme; to be changed to [Large et al., 1994](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) shear mixing scheme in a later development) diff --git a/components/omega/doc/userGuide/VerticalMixingCoeff.md b/components/omega/doc/userGuide/VerticalMixingCoeff.md index a7f92931f1dd..129033aa07f0 100644 --- a/components/omega/doc/userGuide/VerticalMixingCoeff.md +++ b/components/omega/doc/userGuide/VerticalMixingCoeff.md @@ -7,8 +7,8 @@ processes in the ocean. It calculates vertical diffusivity and viscosity coeffic determine how properties (like momentum, heat, salt, and biogeochemical tracers) mix vertically in the ocean model. Currently, Omega offers three different mixing processes within the water column: (1) a constant background mixing value, (2) a convective instability mixing value, and (3) a Richardson number -dependent shear instability driven mixing value from the [Pacanowski and Philander (1981)](https://journals.ametsoc.org/view/journals/phoc/11/11/1520-0485_1981_011_1443_povmin_2_0_co_2.xml) parameterization. These are linearly additive and are describe a bit -more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. +dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These are linearly additive and are describe a bit +more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear mixing calculation. The user-configurable options are the following parameters in the yaml configuration file: @@ -23,17 +23,17 @@ VertMix: TriggerBVF: 0.0 # Squared Brunt-Vaisala frequency threshold Shear: Enable: true # Enables the shear-instability driven mixing option - NuZero: 1.0e-2 # Base viscosity coefficient - Alpha: 5.0 # Stability parameter - Exponent: 2 # Richardson number exponent + BaseShearValue: 0.005 # Base viscosity/diffusivity value + RiCrit: 0.7 # Critical Richardson number + Exponent: 3.0 # Richardson number exponent + RiSmoothLoops: 3 # Number of Richardson number smoothing loops ``` ## Vertical Mixing Processes/Types ### 1. Background Mixing -A constant background mixing value that represents small-scale mixing processes not explicitly resolved by the model. Typically, this is chosen to represent low values of vertical mixing -happening in the ocean's stratified interior. +A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior. ### 2. Convective Mixing @@ -51,20 +51,21 @@ This is different than some current implementations (i.e. in MPAS-Ocean and the ### 3. Shear Mixing -Mixing induced by vertical velocity shear, implemented using the Pacanowski-Philander scheme, through the gradient Richardson number (ratio of buoyancy to shear). +Mixing induced by vertical velocity shear, implemented using the LMD94 scheme, through the gradient Richardson number (ratio of buoyancy to shear). $$ -\nu = \frac{\nu_o}{(1+\alpha Ri)^n} + \nu_b\,, -$$ - -$$ -\kappa = \frac{\nu}{(1+\alpha Ri)} + \kappa_b\,. +\nu_{shear}\,, \kappa_{shear} = = +\begin{cases} +\nu_o \quad \text{ if } Ri_g < 0\\ +\nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 < Ri_g < Ri_{crit}\\ +0.0 \quad \text{ if } Ri_{crit} < Ri_g +\end{cases} $$ -where $Ri$ is defined as: +where $\nu_o$, $Ri_{crit}$, and $p$ are constant parameters set in the `VertMix` section of the yaml file (`BaseShearValue`, `RiCrit`, and `Exponent` under the `Shear` header). $Ri$ is defined as: $$ Ri = \frac{N^2}{\left|\frac{\partial \mathbf{U}}{\partial z}\right|^2}\,, $$ -where $\nu_o$, $\alpha$, $n$, $\nu_b$, $\kappa_b$ are constant parameters set in the `VertMix` section of the yaml file (`NuZero`, `Alpha`, `Exponent` under the `Shear` header and `Viscosity`, `Diffusivity` under the `Background` header). $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $N^2$, $\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri$ of `K` are all defined at the cell center, top interface of layer `K`. $N^2$, $\nu_{shear}$ and $\kappa_{shear}$ are set to zero for the surface layer. In a later development, the shear mixing option will be changed to the interior shear mixing formulation in [Large et al., 1994](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). +where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 filter before being used to calculate the shear. $N^2$, $\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri$ of `K` are all defined at the cell center, top interface of layer `K`. $N^2$, $\nu_{shear}$ and $\kappa_{shear}$ are set to zero for the surface layer. diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index f0e0db9a7055..04e785d0fd8d 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -16,24 +16,39 @@ namespace OMEGA { -ShearMix::ShearMix(const HorzMesh *Mesh, const VertCoord *VCoord) - : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), - ZMid(VCoord->ZMid), NEdgesOnCell(Mesh->NEdgesOnCell), - AreaCell(Mesh->AreaCell), EdgesOnCell(Mesh->EdgesOnCell), - DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge) {} +ShearMix::ShearMix(const VertCoord *VCoord) + : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} ConvectiveMix::ConvectiveMix(const VertCoord *VCoord) : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} +GradRichardsonNum::GradRichardsonNum(const HorzMesh *Mesh, + const VertCoord *VCoord) + : NVertLayers(VCoord->NVertLayers), ZMid(VCoord->ZMid), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), + CellsOnCell(Mesh->CellsOnCell), MinLayerCell(VCoord->MinLayerCell), + MaxLayerCell(VCoord->MaxLayerCell), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), DcEdge(Mesh->DcEdge), + DvEdge(Mesh->DvEdge) {} + +OneTwoOneFilter::OneTwoOneFilter(const VertCoord *VCoord) + : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} + /// Constructor for VertMix VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object const HorzMesh *Mesh, ///< [in] Horizontal mesh const VertCoord *VCoord ///< [in] Vertical coordinate ) - : ComputeVertMixConv(VCoord), ComputeVertMixShear(Mesh, VCoord), Name(Name), - Mesh(Mesh), VCoord(VCoord) { + : ComputeVertMixConv(VCoord), ComputeVertMixShear(VCoord), + ComputeGradRichardsonNum(Mesh, VCoord), ComputeOneTwoOneFilter(VCoord), + Name(Name), Mesh(Mesh), VCoord(VCoord) { VertDiff = Array2DReal("VertDiff", Mesh->NCellsAll, VCoord->NVertLayers); VertVisc = Array2DReal("VertVisc", Mesh->NCellsAll, VCoord->NVertLayers); + GradRichNum = + Array2DReal("GradRichNum", Mesh->NCellsAll, VCoord->NVertLayers); + GradRichNumSmoothed = + Array2DReal("GradRichNumSmoothed", Mesh->NCellsAll, VCoord->NVertLayers); defineFields(); } @@ -54,7 +69,7 @@ void VertMix::destroyInstance() { } /// Initializes the VertMix (Vertical Mixing Coefficients) class and its -/// options. it ASSUMES that HorzMesh was initialized and initializes the +/// options. It ASSUMES that HorzMesh was initialized and initializes the /// VertMix class by using the default mesh, reading the config file, and /// setting parameters for the background, convective, and/or shear mixing /// routines. Returns 0 on success, or an error code if any required option is @@ -136,22 +151,27 @@ void VertMix::init() { LOG_INFO("VertMix::init: Shear mixing is disabled."); } else { LOG_INFO("VertMix::init: Shear mixing is enabled."); - Err += ShearConfig.get("NuZero", - DefVertMix->ComputeVertMixShear.ShearNuZero); - CHECK_ERROR_ABORT( - Err, - "VertMix::init: Parameter Shear:NuZero not found in ShearConfig"); + Err += ShearConfig.get("BaseShearValue", + DefVertMix->ComputeVertMixShear.BaseShearValue); + CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter Shear:BaseShearValue " + "not found in ShearConfig"); - Err += - ShearConfig.get("Alpha", DefVertMix->ComputeVertMixShear.ShearAlpha); + Err += ShearConfig.get("RiCrit", + DefVertMix->ComputeVertMixShear.ShearRiCrit); CHECK_ERROR_ABORT( - Err, "VertMix::init: Parameter Shear:Alpha not found in ShearConfig"); + Err, + "VertMix::init: Parameter Shear:RiCrit not found in ShearConfig"); Err += ShearConfig.get("Exponent", DefVertMix->ComputeVertMixShear.ShearExponent); CHECK_ERROR_ABORT( Err, "VertMix::init: Parameter Shear:Exponent not found in ShearConfig"); + + Err += ShearConfig.get("RiSmoothLoops", + DefVertMix->ComputeVertMixShear.RiSmoothLoops); + CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter Shear:RiSmoothLoops not " + "found in ShearConfig"); } } // end init @@ -161,18 +181,28 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, const Array2DReal &BruntVaisalaFreqSq) { OMEGA_SCOPE(LocVertDiff, VertDiff); /// Create a local view for computation OMEGA_SCOPE(LocVertVisc, VertVisc); /// Create a local view for computation + OMEGA_SCOPE(LocGradRichNum, + GradRichNum); /// Local view for computation OMEGA_SCOPE( LocComputeVertMixConv, ComputeVertMixConv); /// Local view for Convective VertMix computation OMEGA_SCOPE( LocComputeVertMixShear, ComputeVertMixShear); /// Local view for Shear VertMix computation + OMEGA_SCOPE( + LocComputeGradRichardsonNum, + ComputeGradRichardsonNum); /// Local view for GradRichNum computation + OMEGA_SCOPE(LocOneTwoOneFilter, + ComputeOneTwoOneFilter); /// Local view for 1-2-1 filter OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(NVertLayers, VCoord->NVertLayers); /// Initialize VertDiff and VertVisc to background values deepCopy(LocVertDiff, BackDiff); deepCopy(LocVertVisc, BackVisc); + deepCopy(LocGradRichNum, 100.0_Real); + deepCopy(GradRichNumSmoothed, 100.0_Real); /// Dispatch to the correct VertMix calculation if (LocComputeVertMixShear.Enabled && LocComputeVertMixConv.Enabled) { @@ -185,16 +215,52 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixShear( - LocVertDiff, LocVertVisc, ICell, KChunk, NormalVelocity, - TangentialVelocity, BruntVaisalaFreqSq); LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, KChunk, BruntVaisalaFreqSq); + LocComputeGradRichardsonNum( + LocGradRichNum, ICell, KChunk, NormalVelocity, + TangentialVelocity, BruntVaisalaFreqSq); + }); + LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); + LocGradRichNum(ICell, NVertLayers) = + LocGradRichNum(ICell, NVertLayers - 1); + }); + deepCopy(GradRichNumSmoothed, LocGradRichNum); + for (int SmoothLoop = 0; + SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + if (SmoothLoop == 0) + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + LocGradRichNum); + else + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + GradRichNumSmoothed); + }); + }); + } + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, + KChunk, GradRichNumSmoothed); }); }); } else if (LocComputeVertMixShear.Enabled) { parallelForOuter( - "VertMix-ShearOnly", {Mesh->NCellsAll}, + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { const int KMin = MinLayerCell(ICell); const int KMax = MaxLayerCell(ICell); @@ -202,10 +268,45 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixShear( - LocVertDiff, LocVertVisc, ICell, KChunk, NormalVelocity, + LocComputeGradRichardsonNum( + LocGradRichNum, ICell, KChunk, NormalVelocity, TangentialVelocity, BruntVaisalaFreqSq); }); + LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); + LocGradRichNum(ICell, NVertLayers) = + LocGradRichNum(ICell, NVertLayers - 1); + }); + deepCopy(GradRichNumSmoothed, LocGradRichNum); + for (int SmoothLoop = 0; + SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + if (SmoothLoop == 0) + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + LocGradRichNum); + else + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + GradRichNumSmoothed); + }); + }); + } + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, + KChunk, GradRichNumSmoothed); + }); }); } else if (LocComputeVertMixConv.Enabled) { parallelForOuter( @@ -221,24 +322,27 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, KChunk, BruntVaisalaFreqSq); }); }); - } else { - parallelFor( - "VertMix-Background", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { - LocVertDiff(ICell, 0) = 0.0_Real; - LocVertVisc(ICell, 0) = 0.0_Real; - }); } + parallelFor( + "VertMix-Surface", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + LocVertDiff(ICell, 0) = 0.0_Real; + LocVertVisc(ICell, 0) = 0.0_Real; + }); } /// Define IO fields and metadata for output void VertMix::defineFields() { /// Set field names (append Name if not default) - VertDiffFldName = "VertDiff"; - VertViscFldName = "VertVisc"; + VertDiffFldName = "VertDiff"; + VertViscFldName = "VertVisc"; + GradRichNumFldName = "GradRichNum"; + GradRichNumSmoothedFldName = "GradRichNumSmoothed"; if (Name != "Default") { VertDiffFldName.append(Name); VertViscFldName.append(Name); + GradRichNumFldName.append(Name); + GradRichNumSmoothedFldName.append(Name); } /// Create fields for state variables @@ -274,6 +378,30 @@ void VertMix::defineFields() { NDims, // Number of dimensions DimNames // Dimension names ); + /// Create and register the GradRichNum field + auto GradRichNumField = + Field::create(GradRichNumFldName, // Field name + "Gradient Richardson number", // Long Name + "dimensionless", // Units + "sea_water_gradient_richardson_number", // CF-ish Name + std::numeric_limits::min(), // Min valid value + std::numeric_limits::max(), // Max valid value + FillValue, // Scalar used for undefined entries + NDims, // Number of dimensions + DimNames // Dimension names + ); + /// Create and register the GradRichNumSmoothed field + auto GradRichNumSmoothedField = Field::create( + GradRichNumSmoothedFldName, // Field name + "Smoothed Gradient Richardson number", // Long Name + "dimensionless", // Units + "sea_water_gradient_richardson_number_smoothed", // CF-ish Name + std::numeric_limits::min(), // Min valid value + std::numeric_limits::max(), // Max valid value + FillValue, // Scalar used for undefined entries + NDims, // Number of dimensions + DimNames // Dimension names + ); // Create a field group for the vertmix-specific state fields VertMixGroupName = "VertMix"; @@ -285,10 +413,14 @@ void VertMix::defineFields() { // Add fields to the VertMix group VertMixGroup->addField(VertDiffFldName); VertMixGroup->addField(VertViscFldName); + VertMixGroup->addField(GradRichNumFldName); + VertMixGroup->addField(GradRichNumSmoothedFldName); // Attach Kokkos views to the fields VertDiffField->attachData(VertDiff); VertViscField->attachData(VertVisc); + GradRichNumField->attachData(GradRichNum); + GradRichNumSmoothedField->attachData(GradRichNumSmoothed); } // end defineIOFields diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index bd3abcd1b185..ac22f8e72ab7 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -65,21 +65,19 @@ class ShearMix { bool Enabled = true; ///< Enable shear mixing flag // Shear mixing parameters - Real ShearNuZero = - 0.005; ///< Numerator of Pacanowski and Philander (1981) Eq (1). - Real ShearAlpha = 5.0; ///< Alpha value used in Pacanowski and Philander - ///< (1981) Eqs (1) and (2). + Real BaseShearValue = 0.005; ///< Base shear vertical viscosity and + ///< diffusivity (m^2 s^-1) of LMG94 + Real ShearRiCrit = 0.7; ///< Critical Richardson number of LMG94 Real ShearExponent = - 2.0; /// Exponent value used in Pacanowski and Philander (1981) Eqs (1). + 3.0; /// Exponent value used interior shear mixing calculation of LMG94 + I4 RiSmoothLoops = 3; ///< Number of smoothing loops for Richardson number /// Constructor for ShearMix - ShearMix(const HorzMesh *Mesh, const VertCoord *VCoord); + ShearMix(const VertCoord *VCoord); KOKKOS_FUNCTION void operator()(Array2DReal VertDiff, Array2DReal VertVisc, I4 ICell, I4 KChunk, - const Array2DReal &NormalVelocity, - const Array2DReal &TangentialVelocity, - const Array2DReal &BruntVaisalaFreqSq) const { + const Array2DReal &GradRichNumSmoothed) const { const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); @@ -90,43 +88,144 @@ class ShearMix { VertVisc(ICell, K) = 0.0_Real; VertDiff(ICell, K) = 0.0_Real; } else { - Real ShearViscVal = 0.0; - Real InvAreaCell = 1.0_Real / AreaCell(ICell); - Real ShearSquared = 0.0; - for (int J = 0; J < NEdgesOnCell(ICell); ++J) { - I4 JEdge = EdgesOnCell(ICell, J); - Real Factor = - 0.5_Real * DcEdge(JEdge) * DvEdge(JEdge) * InvAreaCell; - Real DelNormVel = - NormalVelocity(JEdge, K - 1) - NormalVelocity(JEdge, K); - Real DelTangVel = TangentialVelocity(JEdge, K - 1) - - TangentialVelocity(JEdge, K); - ShearSquared = ShearSquared + Factor * (DelNormVel * DelNormVel + - DelTangVel * DelTangVel); + if (GradRichNumSmoothed(ICell, K) < 0.0_Real) { + VertDiff(ICell, K) += BaseShearValue; + VertVisc(ICell, K) += BaseShearValue; + } else if (GradRichNumSmoothed(ICell, K) >= 0.0_Real && + GradRichNumSmoothed(ICell, K) < ShearRiCrit) { + VertDiff(ICell, K) += + Kokkos::pow( + 1.0_Real - + (GradRichNumSmoothed(ICell, K) / ShearRiCrit) * + (GradRichNumSmoothed(ICell, K) / ShearRiCrit), + ShearExponent) * + BaseShearValue; + VertVisc(ICell, K) += + Kokkos::pow( + 1.0_Real - + (GradRichNumSmoothed(ICell, K) / ShearRiCrit) * + (GradRichNumSmoothed(ICell, K) / ShearRiCrit), + ShearExponent) * + BaseShearValue; + } else { + VertDiff(ICell, K) += 0.0_Real; + VertVisc(ICell, K) += 0.0_Real; } - Real DelZMid = ZMid(ICell, K - 1) - ZMid(ICell, K); - ShearSquared = ShearSquared / (DelZMid * DelZMid); - - Real RichardsonNum = BruntVaisalaFreqSq(ICell, K) / - Kokkos::max(1.0e-12_Real, ShearSquared); - - ShearViscVal = - ShearNuZero / Kokkos::pow(1.0_Real + ShearAlpha * RichardsonNum, - ShearExponent); - VertVisc(ICell, K) += ShearViscVal; - VertDiff(ICell, K) += - VertVisc(ICell, K) / (1.0_Real + ShearAlpha * RichardsonNum); } } } private: - Array1DReal DcEdge; - Array1DReal DvEdge; - Array1DReal AreaCell; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + +/// Class for Gradient Richardson Number calculation +class GradRichardsonNum { + public: + /// constructor declaration + GradRichardsonNum(const HorzMesh *Mesh, const VertCoord *VCoord); + + // The functor takes the full arrays of Richardson number (inout), + // the index ICell, and normal and tangential velocities as inputs, + // and outputs the Richardson number. + KOKKOS_FUNCTION void + operator()(Array2DReal GradRichNum, I4 ICell, I4 KChunk, + const Array2DReal &NormalVelocity, + const Array2DReal &TangentialVelocity, + const Array2DReal &BruntVaisalaFreqSq) const { + + const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + + Real GradRichNumNorm[VecLength]; + Real GradRichNumTmp[VecLength]; + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + GradRichNumNorm[K] = 1.0e-12_Real; + GradRichNumTmp[K] = 100.0_Real; + } + + for (int J = 0; J < NEdgesOnCell(ICell); ++J) { + I4 JEdge = EdgesOnCell(ICell, J); + I4 JCell = CellsOnCell(ICell, J); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + if (K < 1 || K >= NVertLayers) + continue; + + Real DNormVel = + NormalVelocity(JEdge, K - 1) - NormalVelocity(JEdge, K); + Real DTanVel = + TangentialVelocity(JEdge, K - 1) - TangentialVelocity(JEdge, K); + Real DzEdge = 0.5_Real * (ZMid(ICell, K - 1) + ZMid(JCell, K - 1) - + ZMid(ICell, K) - ZMid(JCell, K)); + Real ShearSquared = + (DNormVel * DNormVel + DTanVel * DTanVel) / (DzEdge * DzEdge); + Real RiEdge = + Kokkos::max(0.0_Real, + 0.5_Real * (BruntVaisalaFreqSq(ICell, K) + + BruntVaisalaFreqSq(JCell, K))) / + (ShearSquared + 1.0e-12_Real); + + Real Weight = 0.25_Real * DcEdge(JEdge) * DvEdge(JEdge); + GradRichNumNorm[K] = GradRichNumNorm[K] + Weight; + GradRichNumTmp[K] = GradRichNumTmp[K] + Weight * RiEdge; + } + } + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + GradRichNum(ICell, K) = GradRichNumTmp[K] / GradRichNumNorm[K]; + } + } + + private: Array2DReal ZMid; - Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; + Array2DI4 CellsOnCell; + Array2DI4 CellsOnEdge; + Array1DI4 NEdgesOnCell; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; + Array1DReal DcEdge; + Array1DReal DvEdge; + I4 NVertLayers; + I4 NCellsAll; +}; + +/// Class for Gradient Richardson Number calculation +class OneTwoOneFilter { + public: + /// constructor declaration + OneTwoOneFilter(const VertCoord *VCoord); + // The functor takes the full arrays of Richardson number (inout), + // the index ICell, and normal and tangential velocities as inputs, + // and outputs the Richardson number. + KOKKOS_FUNCTION void operator()(Array2DReal VarOut, I4 ICell, I4 KChunk, + const Array2DReal &VarIn) const { + + const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + if (K < MinLayerCell(ICell) || K > MaxLayerCell(ICell)) { + VarOut(ICell, K) = VarIn(ICell, K); + } else { + // apply 1-2-1 filter + VarOut(ICell, K) = + (VarIn(ICell, K - 1) + 2.0_Real * VarIn(ICell, K) + + VarIn(ICell, K + 1)) / + 4.0_Real; + } + } + } + + private: Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; }; @@ -140,13 +239,20 @@ class VertMix { /// Destroy instance (frees Kokkos views) static void destroyInstance(); - Array2DReal VertDiff; ///< Vertical diffusivity field (m^2 s^-1) - Array2DReal VertVisc; ///< Vertical viscosity field (m^2 s^-1) - - std::string VertDiffFldName; ///< Field name for vertical diffusivity - std::string VertViscFldName; ///< Field name for vertical viscosity - std::string VertMixGroupName; ///< VertMix group name (for config) - std::string Name; ///< Name of this VertMix instance + Array2DReal VertDiff; ///< Vertical diffusivity field (m^2 s^-1) + Array2DReal VertVisc; ///< Vertical viscosity field (m^2 s^-1) + Array2DReal GradRichNum; ///< Gradient Richardson number field + Array2DReal + GradRichNumSmoothed; ///< Smoothed Gradient Richardson number field + + std::string VertDiffFldName; ///< Field name for vertical diffusivity + std::string VertViscFldName; ///< Field name for vertical viscosity + std::string + GradRichNumFldName; ///< Field name for gradient Richardson number + std::string GradRichNumSmoothedFldName; ///< Field name for smoothed gradient + ///< Richardson number + std::string VertMixGroupName; ///< VertMix group name (for config) + std::string Name; ///< Name of this VertMix instance // Background mixing parameters Real BackDiff = 1.0e-5; ///< Background vertical diffusivity (m^2 s^-1) @@ -155,6 +261,10 @@ class VertMix { ConvectiveMix ComputeVertMixConv; ///< Functor for Convective VertMix calculation ShearMix ComputeVertMixShear; ///< Functor for Shear VertMix calculation + GradRichardsonNum + ComputeGradRichardsonNum; ///< Functor for Gradient Richardson Number + ///< calculation + OneTwoOneFilter ComputeOneTwoOneFilter; ///< Functor for 1-2-1 filtering /// Compute vertical diffusivity and viscosity for all cells/layers void computeVertMix(const Array2DReal &NormalVelocity, diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index 9c0cf2101867..905ec1801674 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -31,29 +31,30 @@ constexpr int NVertLayers = 60; /// Values to test against const Real VertDiffExpValueN = - 1.0393923290498872; // Expected value for diffusivity for positive BVF + 1.00501; // Expected value for diffusivity for positive BVF const Real VertViscExpValueN = - 1.0198269656595984; // Expected value for viscosity for positive BVF + 1.0051; // Expected value for viscosity for positive BVF const Real VertDiffExpValueP = - 0.0015685660274841844; // Expected value for diffusivity for negative BVF + 0.003882748571163051; // Expected value for diffusivity for negative BVF const Real VertViscExpValueP = - 0.002332474675614262; // Expected value for viscosity for negative BVF + 0.003972748571163051; // Expected value for viscosity for negative BVF const Real VertDiffBackExp = 1.0e-5; // Expected value for background diffusivity const Real VertViscBackExp = 1.0e-4; // Expected value for background viscosity -const Real VertDiffConvExp = +const Real VertConvExp = 1.0; // Expected value for convective diffusivity/viscosity -const Real VertDiffShearExp = - 0.039183698912901; // Expected value for shear diffusivity -const Real VertViscShearExp = - 0.01972696565959843; // Expected value for shear viscosity +const Real VertShearExp = + 0.00387274859; // Expected value for shear diffusivity/viscosity +const Real VertShearBaseExp = + 0.005; // Expected value for shear diffusivity/viscosity +const Real RiExpValue = 0.2; // Expected value for gradient Richardson number /// Test input values const Real BVFP = 0.1; // Positive Brunt-Vaisala frequency in s^-2 const Real BVFN = -0.1; // Negative Brunt-Vaisala frequency in s^-2 const Real NV = 1.0; // Normal velocity in m/s const Real TV = 1.0; // Tangential velocity in m/s -const Real RTol = 1e-8; // Relative tolerance for isApprox checks +const Real RTol = 1e-7; // Relative tolerance for isApprox checks /// The initialization routine for VertMix testing. It calls various /// init routines, including the creation of the default decomposition. @@ -98,6 +99,225 @@ void initVertMixTest() { ABORT_ERROR("VertMixTest: VertMix retrieval FAIL"); } +void testGradRichNum() { + /// Get mesh and coordinate info + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + I4 NCellsAll = Mesh->NCellsAll; + I4 NEdgesAll = Mesh->NEdgesAll; + OMEGA_SCOPE(ZMid, VCoord->ZMid); + OMEGA_SCOPE(NEdgesOnCell, Mesh->NEdgesOnCell); + OMEGA_SCOPE(AreaCell, Mesh->AreaCell); + OMEGA_SCOPE(DcEdge, Mesh->DcEdge); + OMEGA_SCOPE(DvEdge, Mesh->DvEdge); + OMEGA_SCOPE(CellsOnCell, Mesh->CellsOnCell); + + /// Get VertMix instance to test + VertMix *TestVertMix = VertMix::getInstance(); + + /// Create and fill ocean state arrays + auto NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesAll, NVertLayers); + auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesAll, NVertLayers); + auto BruntVaisalaFreqSqCell = + Array2DReal("BruntVaisalaFreqSqCell", NCellsAll, NVertLayers); + /// Use deep copy to initialize results + deepCopy(NormalVelEdge, NV); + deepCopy(TangVelEdge, TV); + deepCopy(BruntVaisalaFreqSqCell, BVFP); + deepCopy(TestVertMix->GradRichNum, 0.0); + + parallelFor( + "populateArrays", {NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + ZMid(ICell, K) = -K; + NEdgesOnCell(ICell) = 5; + AreaCell(ICell) = 3.6e10_Real; + }); + + // filling CellsOnCell with simple mapping for this test + parallelFor( + "populateArrays", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + CellsOnCell(ICell, 0) = ICell; + CellsOnCell(ICell, 1) = ICell; + CellsOnCell(ICell, 2) = ICell; + CellsOnCell(ICell, 3) = ICell; + CellsOnCell(ICell, 4) = ICell; + }); + + parallelFor( + "populateArrays", {NEdgesAll, NVertLayers}, + KOKKOS_LAMBDA(I4 IEdge, I4 K) { + NormalVelEdge(IEdge, K) = NormalVelEdge(IEdge, K) + 0.5 * K; + TangVelEdge(IEdge, K) = TangVelEdge(IEdge, K) + 0.5 * K; + DcEdge(IEdge) = 2.0e5_Real; + DvEdge(IEdge) = 1.45e5_Real; + }); + + /// Compute gradient Richardson number + TestVertMix->ComputeVertMixShear.Enabled = true; + TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, + BruntVaisalaFreqSqCell); + + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + + /// Check all array values against expected value + int NumMismatches = 0; + OMEGA_SCOPE(GradRichNum, TestVertMix->GradRichNum); + parallelReduceOuter( + "CheckGradRichNum", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (!isApprox(GradRichNum(ICell, K), RiExpValue, RTol)) + InnerCount++; + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); + }, + NumMismatches); + + // If test fails, print bad values and abort + if (NumMismatches != 0) { + auto GradRichNumH = createHostMirrorCopy(GradRichNum); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (!isApprox(GradRichNumH(I, K), RiExpValue, RTol)) + LOG_ERROR("TestVertMix: GradRichNum Bad Value: " + "GradRichNum({},{}) = {}; Expected {}", + I, K, GradRichNumH(I, K), RiExpValue); + } + } + ABORT_ERROR("TestVertMix: GradRichNum FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMix: GradRichNum PASS"); + } + + return; +} + +void testOneTwoOneFilter() { + /// Get mesh and coordinate info + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + I4 NCellsAll = Mesh->NCellsAll; + I4 NChunks = VCoord->NVertLayers / VecLength; + OMEGA_SCOPE(ZMid, VCoord->ZMid); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + + /// Get VertMix instance to test + VertMix *TestVertMix = VertMix::getInstance(); + + /// Create and fill ocean state arrays + auto GradRichNumSmoothed = + Array2DReal("GradRichNumSmoothed", NCellsAll, NVertLayers); + auto GradRichNum = Array2DReal("GradRichNum", NCellsAll, NVertLayers); + /// Use deep copy to initialize results + deepCopy(GradRichNumSmoothed, 1.0); + deepCopy(GradRichNum, 1.0); + + // Populate GradRichNum with alternating +1.0 and -1.0 values in vertical + // GradRichNumSmoothed should smooth these to 0.0 + parallelFor( + "populateArrays", {NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + if (K > 0) { + GradRichNum(ICell, K) = -1.0 * GradRichNum(ICell, K - 1); + } + }); + + parallelFor( + "setMinMax", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + MinLayerCell(ICell) = 1; // avoid K=0 boundary + MaxLayerCell(ICell) = VCoord->NVertLayers - 2; // avoid top boundary + }); + + // Apply the 1-2-1 filter to each cell + parallelFor( + "ApplyOneTwoOneFilter", {NCellsAll, NChunks}, + KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { + TestVertMix->ComputeOneTwoOneFilter(GradRichNumSmoothed, ICell, + KChunk, GradRichNum); + }); + + /// Check all array values against expected value + int NumMismatches = 0; + parallelReduceOuter( + "CheckGradRichNum", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { + int NumMismatchesCol; + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (K >= 1 && K < NVertLayers - 1) { + // Interior layers should be smoothed to 0.0 + if (!isApprox(GradRichNumSmoothed(ICell, K), 0.0_Real, + RTol)) + InnerCount++; + } else { + // Boundary layers (K==0 or K==NVertLayers-1) should be the + // same as input + if (!isApprox(GradRichNumSmoothed(ICell, K), + GradRichNum(ICell, K), RTol)) + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); + }, + NumMismatches); + + // If test fails, print bad values and abort + if (NumMismatches != 0) { + auto GradRichNumH = createHostMirrorCopy(GradRichNum); + auto GradRichNumSmoothedH = createHostMirrorCopy(GradRichNumSmoothed); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K >= 1 && K < NVertLayers - 1) { + // Interior layers should be smoothed to 0.0 + if (!isApprox(GradRichNumSmoothedH(I, K), 0.0, RTol)) + LOG_ERROR("TestVertMix: GradRichNumSmoothed Bad Value: " + "GradRichNumSmoothed({},{}) = {}; Expected {}", + I, K, GradRichNumSmoothedH(I, K), 0.0); + } else { + // Boundary layers (K==0 or K==NVertLayers-1) should be copied + // from input + if (!isApprox(GradRichNumSmoothedH(I, K), GradRichNumH(I, K), + RTol)) + LOG_ERROR("TestVertMix: GradRichNumSmoothed Bad Value: " + "GradRichNumSmoothed({},{}) = {}; Expected {}", + I, K, GradRichNumSmoothedH(I, K), + GradRichNumH(I, K)); + } + } + } + ABORT_ERROR("TestVertMix: GradRichNumSmoothed FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMix: GradRichNumSmoothed PASS"); + } + + return; +} + void testBackVertMix() { // Get mesh and coordinate info const auto Mesh = HorzMesh::getDefault(); @@ -150,7 +370,7 @@ void testBackVertMix() { Array2DReal BackVertVisc = TestVertMix->VertVisc; Array2DReal BackVertDiff = TestVertMix->VertDiff; - /// Check total Visc against linear addition of components + /// Check Visc against expected value int NumMismatches = 0; parallelReduceOuter( "CheckVertMixMatrix-BackgroundVisc", {Mesh->NCellsAll}, @@ -180,12 +400,15 @@ void testBackVertMix() { }, NumMismatches); - if (NumMismatches != 0) + if (NumMismatches != 0) { ABORT_ERROR("TestVertMixBack: VertVisc FAIL, " "expected {}, got {} with {} mismatches", VertViscBackExp, BackVertVisc(1, 1), NumMismatches); + } else { + LOG_INFO("TestVertMixBack: VertVisc PASS"); + } - /// Check total Diff against linear addition of components + /// Check Diff against expected value NumMismatches = 0; parallelReduceOuter( "CheckVertMixMatrix-BackgroundDiff", {Mesh->NCellsAll}, @@ -220,13 +443,14 @@ void testBackVertMix() { ABORT_ERROR("TestVertMixBack: VertDiff FAIL, " "expected {}, got {} with {} mismatches", VertDiffBackExp, BackVertDiffH(1, 1), NumMismatches); + } else { + LOG_INFO("TestVertMixBack: VertDiff PASS"); } return; } void testConvVertMix() { - // Get mesh and coordinate info const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); @@ -234,49 +458,46 @@ void testConvVertMix() { I4 NCellsAll = Mesh->NCellsAll; I4 NEdgesAll = Mesh->NEdgesAll; OMEGA_SCOPE(ZMid, VCoord->ZMid); + I4 NChunks = VCoord->NVertLayers / VecLength; /// Get VertMix instance to test VertMix *TestVertMix = VertMix::getInstance(); /// Create and fill ocean state arrays - auto NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesAll, NVertLayers); - auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesAll, NVertLayers); - auto BruntVaisalaFreqSqCell = - Array2DReal("BruntVaisalaFreqSqCell", NCellsAll, NVertLayers); + auto BruntVaisalaFreqSqIn = + Array2DReal("BruntVaisalaFreqSqIn", NCellsAll, NVertLayers); + auto VertDiffOut = Array2DReal("VertDiffOut", NCellsAll, NVertLayers); + auto VertViscOut = Array2DReal("VertViscOut", NCellsAll, NVertLayers); /// Use deep copy to initialize with the ref value - deepCopy(NormalVelEdge, NV); - deepCopy(TangVelEdge, TV); - deepCopy(BruntVaisalaFreqSqCell, BVFN); - deepCopy(TestVertMix->VertDiff, 0.0); - deepCopy(TestVertMix->VertVisc, 0.0); + deepCopy(BruntVaisalaFreqSqIn, 0.0); + deepCopy(VertDiffOut, 0.0); + deepCopy(VertViscOut, 0.0); + // Populate arrays: positive BVF in lower half (conv off), + // negative in upper half (conv on) parallelFor( "populateArrays", {NCellsAll, NVertLayers}, - KOKKOS_LAMBDA(I4 ICell, I4 K) { ZMid(ICell, K) = -K; }); - - parallelFor( - "populateArrays", {NEdgesAll, NVertLayers}, - KOKKOS_LAMBDA(I4 IEdge, I4 K) { - NormalVelEdge(IEdge, K) = NormalVelEdge(IEdge, K) + 0.5 * K; - TangVelEdge(IEdge, K) = TangVelEdge(IEdge, K) + 0.5 * K; + KOKKOS_LAMBDA(I4 ICell, I4 K) { + if (K < 30) { + BruntVaisalaFreqSqIn(ICell, K) = -0.2; + } else { + BruntVaisalaFreqSqIn(ICell, K) = 0.2; + } }); /// Compute only convective vertical viscosity and diffusivity - TestVertMix->BackDiff = 0.0; - TestVertMix->BackVisc = 0.0; - TestVertMix->ComputeVertMixConv.Enabled = true; - TestVertMix->ComputeVertMixShear.Enabled = false; - TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, - BruntVaisalaFreqSqCell); + parallelFor( + "ApplyVertMixConv", {NCellsAll, NChunks}, + KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { + TestVertMix->ComputeVertMixConv(VertDiffOut, VertViscOut, ICell, + KChunk, BruntVaisalaFreqSqIn); + }); const auto &MinLayerCell = VCoord->MinLayerCell; const auto &MaxLayerCell = VCoord->MaxLayerCell; - Array2DReal ConvVertVisc = TestVertMix->VertVisc; - Array2DReal ConvVertDiff = TestVertMix->VertDiff; - - /// Check total Visc against linear addition of components + /// Check Visc against expected value int NumMismatches = 0; parallelReduceOuter( "CheckVertMixMatrix-ConvectiveVisc", {Mesh->NCellsAll}, @@ -291,11 +512,13 @@ void testConvVertMix() { const int K = KMin + KOff; if (K == 0) { // Surface layer should be zero - if (ConvVertVisc(ICell, K) != 0.0_Real) + if (VertViscOut(ICell, K) != 0.0_Real) + InnerCount++; + } else if (K < 30) { + if (!isApprox(VertViscOut(ICell, K), VertConvExp, RTol)) InnerCount++; } else { - if (!isApprox(ConvVertVisc(ICell, K), VertDiffConvExp, - RTol)) + if (!isApprox(VertViscOut(ICell, K), 0.0_Real, RTol)) InnerCount++; } }, @@ -307,13 +530,37 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ConvVertViscH = createHostMirrorCopy(ConvVertVisc); - ABORT_ERROR("TestVertMixConv: VertVisc FAIL, " - "expected {}, got {} with {} mismatches", - VertDiffConvExp, ConvVertViscH(1, 1), NumMismatches); + auto VertViscOutH = createHostMirrorCopy(VertViscOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } else if (K < 30) { + // Interior layers + if (!isApprox(VertViscOutH(I, K), VertConvExp, RTol)) + LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), VertConvExp); + } else { + // Interior layers + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } + } + } + ABORT_ERROR("TestVertMixConv: VertVisc FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMixConv: VertVisc PASS"); } - /// Check total Diff against linear addition of components + /// Check Diff against expected value NumMismatches = 0; parallelReduceOuter( "CheckVertMixMatrix-ConvectiveDiff", {Mesh->NCellsAll}, @@ -328,11 +575,13 @@ void testConvVertMix() { const int K = KMin + KOff; if (K == 0) { // Surface layer should be zero - if (ConvVertDiff(ICell, K) != 0.0_Real) + if (VertDiffOut(ICell, K) != 0.0_Real) + InnerCount++; + } else if (K < 30) { + if (!isApprox(VertDiffOut(ICell, K), VertConvExp, RTol)) InnerCount++; } else { - if (!isApprox(ConvVertDiff(ICell, K), VertDiffConvExp, - RTol)) + if (!isApprox(VertDiffOut(ICell, K), 0.0_Real, RTol)) InnerCount++; } }, @@ -344,10 +593,34 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ConvVertDiffH = createHostMirrorCopy(ConvVertDiff); - ABORT_ERROR("TestVertMixConv: VertDiff FAIL, " - "expected {}, got {} with {} mismatches", - VertDiffConvExp, ConvVertDiffH(1, 1), NumMismatches); + auto VertDiffOutH = createHostMirrorCopy(VertDiffOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } else if (K < 30) { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), VertConvExp, RTol)) + LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), VertConvExp); + } else { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } + } + } + ABORT_ERROR("TestVertMixConv: VertDiff FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMixConv: VertDiff PASS"); } return; @@ -359,61 +632,50 @@ void testShearVertMix() { const auto VCoord = VertCoord::getDefault(); VCoord->NVertLayers = NVertLayers; I4 NCellsAll = Mesh->NCellsAll; - I4 NEdgesAll = Mesh->NEdgesAll; - OMEGA_SCOPE(ZMid, VCoord->ZMid); - OMEGA_SCOPE(NEdgesOnCell, Mesh->NEdgesOnCell); - OMEGA_SCOPE(AreaCell, Mesh->AreaCell); - OMEGA_SCOPE(DcEdge, Mesh->DcEdge); - OMEGA_SCOPE(DvEdge, Mesh->DvEdge); + I4 NChunks = VCoord->NVertLayers / VecLength; /// Get VertMix instance to test VertMix *TestVertMix = VertMix::getInstance(); /// Create and fill ocean state arrays - auto NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesAll, NVertLayers); - auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesAll, NVertLayers); - auto BruntVaisalaFreqSqCell = - Array2DReal("BruntVaisalaFreqSqCell", NCellsAll, NVertLayers); + auto GradRichNumSmoothedIn = + Array2DReal("GradRichNumSmoothedIn", NCellsAll, NVertLayers); + auto VertDiffOut = Array2DReal("VertDiffOut", NCellsAll, NVertLayers); + auto VertViscOut = Array2DReal("VertViscOut", NCellsAll, NVertLayers); /// Use Kokkos::deep_copy to fill the entire view with the ref value - deepCopy(NormalVelEdge, NV); - deepCopy(TangVelEdge, TV); - deepCopy(BruntVaisalaFreqSqCell, BVFN); - deepCopy(TestVertMix->VertDiff, 0.0); - deepCopy(TestVertMix->VertVisc, 0.0); + deepCopy(GradRichNumSmoothedIn, 0.0); + deepCopy(VertDiffOut, 0.0); + deepCopy(VertViscOut, 0.0); + // Populate arrays: negative Ri in upper third (base shear value), + // positive in middle third (altered shear value), large positive + // in lower third (no shear) parallelFor( "populateArrays", {NCellsAll, NVertLayers}, KOKKOS_LAMBDA(I4 ICell, I4 K) { - ZMid(ICell, K) = -K; - NEdgesOnCell(ICell) = 5; - AreaCell(ICell) = 3.6e10_Real; + if (K < 20) { + GradRichNumSmoothedIn(ICell, K) = -0.2; + } else if (K >= 20 && K < 40) { + GradRichNumSmoothedIn(ICell, K) = 0.2; + } else { + GradRichNumSmoothedIn(ICell, K) = 10.0; + } }); + /// Compute only shear vertical viscosity and diffusivity + TestVertMix->ComputeVertMixShear.ShearExponent = 3.0; parallelFor( - "populateArrays", {NEdgesAll, NVertLayers}, - KOKKOS_LAMBDA(I4 IEdge, I4 K) { - NormalVelEdge(IEdge, K) = NormalVelEdge(IEdge, K) + 0.5 * K; - TangVelEdge(IEdge, K) = TangVelEdge(IEdge, K) + 0.5 * K; - DcEdge(IEdge) = 2.0e5_Real; - DvEdge(IEdge) = 1.45e5_Real; + "ApplyVertMixShear", {NCellsAll, NChunks}, + KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { + TestVertMix->ComputeVertMixShear(VertDiffOut, VertViscOut, ICell, + KChunk, GradRichNumSmoothedIn); }); - /// Compute only shear vertical viscosity and diffusivity - TestVertMix->BackDiff = 0.0; - TestVertMix->BackVisc = 0.0; - TestVertMix->ComputeVertMixConv.Enabled = false; - TestVertMix->ComputeVertMixShear.Enabled = true; - TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, - BruntVaisalaFreqSqCell); - const auto &MinLayerCell = VCoord->MinLayerCell; const auto &MaxLayerCell = VCoord->MaxLayerCell; - Array2DReal ShearVertVisc = TestVertMix->VertVisc; - Array2DReal ShearVertDiff = TestVertMix->VertDiff; - - /// Check total Visc against linear addition of components + /// Check Visc against expected value int NumMismatches = 0; parallelReduceOuter( "CheckVertMixMatrix-ShearVisc", {Mesh->NCellsAll}, @@ -427,18 +689,17 @@ void testShearVertMix() { INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; if (K == 0) { - if (ShearVertVisc(ICell, K) != 0.0_Real) + if (VertViscOut(ICell, K) != 0.0_Real) InnerCount++; - // K = 1 should have ref value - } else if (K == 1) { - if (!isApprox(ShearVertVisc(ICell, K), VertViscShearExp, + } else if (K < 20) { + if (!isApprox(VertViscOut(ICell, K), VertShearBaseExp, RTol)) InnerCount++; - // otherwise check for invalid values + } else if (K >= 20 && K < 40) { + if (!isApprox(VertViscOut(ICell, K), VertShearExp, RTol)) + InnerCount++; } else { - if (ShearVertVisc(ICell, K) == 0.0 or - Kokkos::isnan(ShearVertVisc(ICell, K)) or - Kokkos::isinf(ShearVertVisc(ICell, K))) + if (!isApprox(VertViscOut(ICell, K), 0.0_Real, RTol)) InnerCount++; } }, @@ -450,13 +711,43 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ShearVertViscH = createHostMirrorCopy(ShearVertVisc); - ABORT_ERROR("TestVertMixShear: VertVisc FAIL, " - "expected {}, got {} with {} mismatches", - VertViscShearExp, ShearVertViscH(1, 1), NumMismatches); + auto VertViscOutH = createHostMirrorCopy(VertViscOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } else if (K < 20) { + // Interior layers + if (!isApprox(VertViscOutH(I, K), VertShearBaseExp, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), VertShearBaseExp); + } else if (K >= 20 && K < 40) { + // Interior layers + if (!isApprox(VertViscOutH(I, K), VertShearExp, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), VertShearExp); + } else { + // Interior layers + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } + } + } + ABORT_ERROR("TestVertMixShear: VertVisc FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMixShear: VertVisc PASS"); } - /// Check total Diff against linear addition of components + /// Check Diff against expected value NumMismatches = 0; parallelReduceOuter( "CheckVertMixMatrix-ShearVisc", {Mesh->NCellsAll}, @@ -470,18 +761,17 @@ void testShearVertMix() { INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; if (K == 0) { - if (ShearVertDiff(ICell, K) != 0.0_Real) + if (VertDiffOut(ICell, K) != 0.0_Real) InnerCount++; - // K = 1 should have ref value - } else if (K == 1) { - if (!isApprox(ShearVertDiff(ICell, K), VertDiffShearExp, + } else if (K < 20) { + if (!isApprox(VertDiffOut(ICell, K), VertShearBaseExp, RTol)) InnerCount++; - // otherwise check for invalid values + } else if (K >= 20 && K < 40) { + if (!isApprox(VertDiffOut(ICell, K), VertShearExp, RTol)) + InnerCount++; } else { - if (ShearVertDiff(ICell, K) == 0.0 or - Kokkos::isnan(ShearVertDiff(ICell, K)) or - Kokkos::isinf(ShearVertDiff(ICell, K))) + if (!isApprox(VertDiffOut(ICell, K), 0.0_Real, RTol)) InnerCount++; } }, @@ -493,10 +783,40 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ShearVertDiffH = createHostMirrorCopy(ShearVertDiff); - ABORT_ERROR("TestVertMixShear: VertDiff FAIL, " - "expected {}, got {} with {} mismatches", - VertDiffShearExp, ShearVertDiffH(1, 1), NumMismatches); + auto VertDiffOutH = createHostMirrorCopy(VertDiffOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } else if (K < 20) { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), VertShearBaseExp, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), VertShearBaseExp); + } else if (K >= 20 && K < 40) { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), VertShearExp, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), VertShearExp); + } else { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } + } + } + ABORT_ERROR("TestVertMixShear: VertDiff FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMixShear: VertDiff PASS"); } return; @@ -515,6 +835,7 @@ void testTotalVertMix() { OMEGA_SCOPE(AreaCell, Mesh->AreaCell); OMEGA_SCOPE(DcEdge, Mesh->DcEdge); OMEGA_SCOPE(DvEdge, Mesh->DvEdge); + OMEGA_SCOPE(CellsOnCell, Mesh->CellsOnCell); /// Get VertMix instance to test VertMix *TestVertMix = VertMix::getInstance(); @@ -533,6 +854,7 @@ void testTotalVertMix() { deepCopy(BruntVaisalaFreqSqCell, BVFP); deepCopy(TestVertMix->VertDiff, 0.0); deepCopy(TestVertMix->VertVisc, 0.0); + deepCopy(TestVertMix->GradRichNumSmoothed, 0.0); parallelFor( "populateArrays", {NCellsAll, NVertLayers}, @@ -542,6 +864,17 @@ void testTotalVertMix() { AreaCell(ICell) = 3.6e10_Real; }); + // current mesh has some CellsOnCell value > NCellsAll, so + // filling CellsOnCell with simple mapping for this test + parallelFor( + "populateArrays", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + CellsOnCell(ICell, 0) = ICell; + CellsOnCell(ICell, 1) = ICell; + CellsOnCell(ICell, 2) = ICell; + CellsOnCell(ICell, 3) = ICell; + CellsOnCell(ICell, 4) = ICell; + }); + parallelFor( "populateArrays", {NEdgesAll, NVertLayers}, KOKKOS_LAMBDA(I4 IEdge, I4 K) { @@ -564,6 +897,7 @@ void testTotalVertMix() { OMEGA_SCOPE(VertDiffP, TestVertMix->VertDiff); OMEGA_SCOPE(VertViscP, TestVertMix->VertVisc); + OMEGA_SCOPE(GradRichNumSmoothed, TestVertMix->GradRichNumSmoothed); /// Check all VertDiff array values against expected value int NumMismatches = 0; @@ -605,6 +939,8 @@ void testTotalVertMix() { ABORT_ERROR("TestVertMixTotal: VertDiffPositive FAIL, " "expected {}, got {} with {} mismatches", VertDiffExpValueP, VertDiffPH(1, 1), NumMismatches); + } else { + LOG_INFO("TestVertMixTotal: VertDiffPos PASS"); } /// Check all VertVisc array values against expected value @@ -647,6 +983,8 @@ void testTotalVertMix() { ABORT_ERROR("TestVertMixTotal: VertViscPositive FAIL, " "expected {}, got {} with {} mismatches", VertViscExpValueP, VertViscPH(1, 1), NumMismatches); + } else { + LOG_INFO("TestVertMixTotal: VertViscPos PASS"); } // Now test with negative BVF @@ -700,6 +1038,8 @@ void testTotalVertMix() { ABORT_ERROR("TestVertMix: VertDiffNegative FAIL, " "expected {}, got {} with {} mismatches", VertDiffExpValueN, VertDiffNH(1, 1), NumMismatches); + } else { + LOG_INFO("TestVertMixTotal: VertDiffNeg PASS"); } /// Check all VertVisc array values against expected value @@ -742,6 +1082,8 @@ void testTotalVertMix() { ABORT_ERROR("TestVertMix: VertViscNegative FAIL, " "expected {}, got {} with {} mismatches", VertViscExpValueN, VertViscNH(1, 1), NumMismatches); + } else { + LOG_INFO("TestVertMixTotal: VertViscNeg PASS"); } return; @@ -760,12 +1102,14 @@ void finalizeVertMixTest() { } // the main tests (all in one to have the same log): -// --> one tests the vertical diffusivity and viscosity +// --> one tests the gradient richardson number calculation +// --> next tests the 1-2-1 filter (smoothing) +// --> next tests the vertical diffusivity and viscosity // with only background on // --> next tests the vertical diffusivity and viscosity -// with only convective on +// for only convective // --> next tests the vertical diffusivity and viscosity -// with only shear on +// for only shear // --> next tests the linear superposition of the // background, convective, and shear contributions void vertMixTest() { @@ -774,6 +1118,8 @@ void vertMixTest() { initVertMixTest(); // test each vertical mix option + testGradRichNum(); + testOneTwoOneFilter(); testBackVertMix(); testConvVertMix(); testShearVertMix(); From d21deb7f4cf90a9241dcc46e93c0927ccdcee8e2 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Tue, 17 Feb 2026 14:53:31 -0800 Subject: [PATCH 2/4] changes for gpu --- components/omega/src/ocn/VertMix.cpp | 84 ++++---------------- components/omega/src/ocn/VertMix.h | 95 ++++++++++++----------- components/omega/test/ocn/VertMixTest.cpp | 31 +++++--- 3 files changed, 87 insertions(+), 123 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 04e785d0fd8d..ea6b221a44d1 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -183,6 +183,8 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, OMEGA_SCOPE(LocVertVisc, VertVisc); /// Create a local view for computation OMEGA_SCOPE(LocGradRichNum, GradRichNum); /// Local view for computation + OMEGA_SCOPE(LocGradRichNumSmoothed, + GradRichNumSmoothed); /// Local view for computation OMEGA_SCOPE( LocComputeVertMixConv, ComputeVertMixConv); /// Local view for Convective VertMix computation @@ -202,12 +204,12 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, deepCopy(LocVertDiff, BackDiff); deepCopy(LocVertVisc, BackVisc); deepCopy(LocGradRichNum, 100.0_Real); - deepCopy(GradRichNumSmoothed, 100.0_Real); + deepCopy(LocGradRichNumSmoothed, 100.0_Real); /// Dispatch to the correct VertMix calculation - if (LocComputeVertMixShear.Enabled && LocComputeVertMixConv.Enabled) { + if (LocComputeVertMixShear.Enabled) { parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + "VertMix-ComputeRi", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { const int KMin = MinLayerCell(ICell); const int KMax = MaxLayerCell(ICell); @@ -215,21 +217,16 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, - KChunk, BruntVaisalaFreqSq); LocComputeGradRichardsonNum( LocGradRichNum, ICell, KChunk, NormalVelocity, TangentialVelocity, BruntVaisalaFreqSq); }); - LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); - LocGradRichNum(ICell, NVertLayers) = - LocGradRichNum(ICell, NVertLayers - 1); }); - deepCopy(GradRichNumSmoothed, LocGradRichNum); + deepCopy(LocGradRichNumSmoothed, LocGradRichNum); for (int SmoothLoop = 0; SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + "VertMix-RiSmooth", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { const int KMin = MinLayerCell(ICell); const int KMax = MaxLayerCell(ICell); @@ -237,80 +234,31 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { if (SmoothLoop == 0) - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - LocGradRichNum); + LocOneTwoOneFilter(LocGradRichNumSmoothed, ICell, + KChunk, LocGradRichNum); else - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - GradRichNumSmoothed); + LocOneTwoOneFilter(LocGradRichNumSmoothed, ICell, + KChunk, LocGradRichNumSmoothed); }); }); } parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, - KChunk, GradRichNumSmoothed); - }); - }); - } else if (LocComputeVertMixShear.Enabled) { - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + "VertMix-Shear", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { const int KMin = MinLayerCell(ICell); const int KMax = MaxLayerCell(ICell); const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeGradRichardsonNum( - LocGradRichNum, ICell, KChunk, NormalVelocity, - TangentialVelocity, BruntVaisalaFreqSq); - }); - LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); - LocGradRichNum(ICell, NVertLayers) = - LocGradRichNum(ICell, NVertLayers - 1); - }); - deepCopy(GradRichNumSmoothed, LocGradRichNum); - for (int SmoothLoop = 0; - SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - if (SmoothLoop == 0) - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - LocGradRichNum); - else - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - GradRichNumSmoothed); - }); - }); - } - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, - KChunk, GradRichNumSmoothed); + KChunk, LocGradRichNumSmoothed); }); }); - } else if (LocComputeVertMixConv.Enabled) { + } + if (LocComputeVertMixConv.Enabled) { parallelForOuter( - "VertMix-ConvOnly", {Mesh->NCellsAll}, + "VertMix-Conv", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { const int KMin = MinLayerCell(ICell); const int KMax = MaxLayerCell(ICell); diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index ac22f8e72ab7..b8c79472bfdc 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -43,14 +43,12 @@ class ConvectiveMix { for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; - if (K == 0) { - VertVisc(ICell, K) = 0.0_Real; - VertDiff(ICell, K) = 0.0_Real; - } else { - if (BruntVaisalaFreqSq(ICell, K) < ConvTriggerBVF) { - VertDiff(ICell, K) += ConvDiff; - VertVisc(ICell, K) += ConvDiff; - } + if (K <= MinLayerCell(ICell) || K >= MaxLayerCell(ICell)) + continue; + + if (BruntVaisalaFreqSq(ICell, K) < ConvTriggerBVF) { + VertDiff(ICell, K) += ConvDiff; + VertVisc(ICell, K) += ConvDiff; } } } @@ -84,33 +82,31 @@ class ShearMix { for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; - if (K == 0) { - VertVisc(ICell, K) = 0.0_Real; - VertDiff(ICell, K) = 0.0_Real; + if (K <= MinLayerCell(ICell) || K >= MaxLayerCell(ICell)) + continue; + + if (GradRichNumSmoothed(ICell, K) < 0.0_Real) { + VertDiff(ICell, K) += BaseShearValue; + VertVisc(ICell, K) += BaseShearValue; + } else if (GradRichNumSmoothed(ICell, K) >= 0.0_Real && + GradRichNumSmoothed(ICell, K) < ShearRiCrit) { + VertDiff(ICell, K) += + Kokkos::pow( + 1.0_Real - + (GradRichNumSmoothed(ICell, K) / ShearRiCrit) * + (GradRichNumSmoothed(ICell, K) / ShearRiCrit), + ShearExponent) * + BaseShearValue; + VertVisc(ICell, K) += + Kokkos::pow( + 1.0_Real - + (GradRichNumSmoothed(ICell, K) / ShearRiCrit) * + (GradRichNumSmoothed(ICell, K) / ShearRiCrit), + ShearExponent) * + BaseShearValue; } else { - if (GradRichNumSmoothed(ICell, K) < 0.0_Real) { - VertDiff(ICell, K) += BaseShearValue; - VertVisc(ICell, K) += BaseShearValue; - } else if (GradRichNumSmoothed(ICell, K) >= 0.0_Real && - GradRichNumSmoothed(ICell, K) < ShearRiCrit) { - VertDiff(ICell, K) += - Kokkos::pow( - 1.0_Real - - (GradRichNumSmoothed(ICell, K) / ShearRiCrit) * - (GradRichNumSmoothed(ICell, K) / ShearRiCrit), - ShearExponent) * - BaseShearValue; - VertVisc(ICell, K) += - Kokkos::pow( - 1.0_Real - - (GradRichNumSmoothed(ICell, K) / ShearRiCrit) * - (GradRichNumSmoothed(ICell, K) / ShearRiCrit), - ShearExponent) * - BaseShearValue; - } else { - VertDiff(ICell, K) += 0.0_Real; - VertVisc(ICell, K) += 0.0_Real; - } + VertDiff(ICell, K) += 0.0_Real; + VertVisc(ICell, K) += 0.0_Real; } } } @@ -152,21 +148,30 @@ class GradRichardsonNum { I4 JCell = CellsOnCell(ICell, J); for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; - if (K < 1 || K >= NVertLayers) - continue; + I4 K1 = K - 1; + I4 K2 = K; + + if (K <= MinLayerCell(ICell)) { + K1 = K; + K2 = K + 1; + } + if (K >= MaxLayerCell(ICell)) { + K1 = K - 2; + K2 = K - 1; + } Real DNormVel = - NormalVelocity(JEdge, K - 1) - NormalVelocity(JEdge, K); + NormalVelocity(JEdge, K1) - NormalVelocity(JEdge, K2); Real DTanVel = - TangentialVelocity(JEdge, K - 1) - TangentialVelocity(JEdge, K); - Real DzEdge = 0.5_Real * (ZMid(ICell, K - 1) + ZMid(JCell, K - 1) - - ZMid(ICell, K) - ZMid(JCell, K)); + TangentialVelocity(JEdge, K1) - TangentialVelocity(JEdge, K2); + Real DzEdge = 0.5_Real * (ZMid(ICell, K1) + ZMid(JCell, K1) - + (ZMid(ICell, K2) + ZMid(JCell, K2))); Real ShearSquared = (DNormVel * DNormVel + DTanVel * DTanVel) / (DzEdge * DzEdge); Real RiEdge = Kokkos::max(0.0_Real, - 0.5_Real * (BruntVaisalaFreqSq(ICell, K) + - BruntVaisalaFreqSq(JCell, K))) / + 0.5_Real * (BruntVaisalaFreqSq(ICell, K2) + + BruntVaisalaFreqSq(JCell, K2))) / (ShearSquared + 1.0e-12_Real); Real Weight = 0.25_Real * DcEdge(JEdge) * DvEdge(JEdge); @@ -213,14 +218,14 @@ class OneTwoOneFilter { for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; - if (K < MinLayerCell(ICell) || K > MaxLayerCell(ICell)) { - VarOut(ICell, K) = VarIn(ICell, K); - } else { + if (K > MinLayerCell(ICell) && K < MaxLayerCell(ICell) - 1) { // apply 1-2-1 filter VarOut(ICell, K) = (VarIn(ICell, K - 1) + 2.0_Real * VarIn(ICell, K) + VarIn(ICell, K + 1)) / 4.0_Real; + } else { + VarOut(ICell, K) = VarIn(ICell, K); } } } diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index 905ec1801674..abd1f9a9ba8e 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -112,6 +112,9 @@ void testGradRichNum() { OMEGA_SCOPE(DcEdge, Mesh->DcEdge); OMEGA_SCOPE(DvEdge, Mesh->DvEdge); OMEGA_SCOPE(CellsOnCell, Mesh->CellsOnCell); + OMEGA_SCOPE(CellsOnEdge, Mesh->CellsOnEdge); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); /// Get VertMix instance to test VertMix *TestVertMix = VertMix::getInstance(); @@ -135,6 +138,12 @@ void testGradRichNum() { AreaCell(ICell) = 3.6e10_Real; }); + parallelFor( + "setMinMax", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + MinLayerCell(ICell) = 0; + MaxLayerCell(ICell) = VCoord->NVertLayers; + }); + // filling CellsOnCell with simple mapping for this test parallelFor( "populateArrays", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { @@ -159,8 +168,8 @@ void testGradRichNum() { TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, BruntVaisalaFreqSqCell); - const auto &MinLayerCell = VCoord->MinLayerCell; - const auto &MaxLayerCell = VCoord->MaxLayerCell; + // const auto &MinLayerCell = VCoord->MinLayerCell; + // const auto &MaxLayerCell = VCoord->MaxLayerCell; /// Check all array values against expected value int NumMismatches = 0; @@ -233,15 +242,17 @@ void testOneTwoOneFilter() { parallelFor( "populateArrays", {NCellsAll, NVertLayers}, KOKKOS_LAMBDA(I4 ICell, I4 K) { - if (K > 0) { - GradRichNum(ICell, K) = -1.0 * GradRichNum(ICell, K - 1); + if (K % 2 == 0) { + GradRichNum(ICell, K) = 1.0; + } else { + GradRichNum(ICell, K) = -1.0; } }); parallelFor( "setMinMax", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { - MinLayerCell(ICell) = 1; // avoid K=0 boundary - MaxLayerCell(ICell) = VCoord->NVertLayers - 2; // avoid top boundary + MinLayerCell(ICell) = 0; + MaxLayerCell(ICell) = VCoord->NVertLayers - 1; }); // Apply the 1-2-1 filter to each cell @@ -265,14 +276,14 @@ void testOneTwoOneFilter() { Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K >= 1 && K < NVertLayers - 1) { + if (K > MinLayerCell(ICell) && K < MaxLayerCell(ICell) - 1) { // Interior layers should be smoothed to 0.0 if (!isApprox(GradRichNumSmoothed(ICell, K), 0.0_Real, RTol)) InnerCount++; } else { - // Boundary layers (K==0 or K==NVertLayers-1) should be the - // same as input + // Boundary layers (K==0 or K==NVertLayers - 1) should be + // the same as input if (!isApprox(GradRichNumSmoothed(ICell, K), GradRichNum(ICell, K), RTol)) InnerCount++; @@ -291,7 +302,7 @@ void testOneTwoOneFilter() { auto GradRichNumSmoothedH = createHostMirrorCopy(GradRichNumSmoothed); for (int I = 0; I < NCellsAll; ++I) { for (int K = 0; K < NVertLayers; ++K) { - if (K >= 1 && K < NVertLayers - 1) { + if (K > MinLayerCell(I) && K < MaxLayerCell(I) - 1) { // Interior layers should be smoothed to 0.0 if (!isApprox(GradRichNumSmoothedH(I, K), 0.0, RTol)) LOG_ERROR("TestVertMix: GradRichNumSmoothed Bad Value: " From 3e301a3dbb212ffb306d2ca0167024e9cfcd03f5 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 16 Mar 2026 14:32:04 -0700 Subject: [PATCH 3/4] adds suggestions from pr review --- components/omega/configs/Default.yml | 2 +- .../omega/doc/devGuide/VerticalMixingCoeff.md | 18 +++++++-------- .../doc/userGuide/VerticalMixingCoeff.md | 19 ++++++++-------- components/omega/src/ocn/VertMix.h | 22 ++++++++++--------- components/omega/test/ocn/VertMixTest.cpp | 2 +- 5 files changed, 32 insertions(+), 31 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index be6b649820df..0e750a020bf4 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -72,7 +72,7 @@ Omega: BaseShearValue: 0.005 RiCrit: 0.7 Exponent: 3.0 - RiSmoothLoops: 3 + RiSmoothLoops: 2 IOStreams: InitialVertCoord: UsePointerFile: false diff --git a/components/omega/doc/devGuide/VerticalMixingCoeff.md b/components/omega/doc/devGuide/VerticalMixingCoeff.md index 79736261bc72..22b3b6c5eed3 100644 --- a/components/omega/doc/devGuide/VerticalMixingCoeff.md +++ b/components/omega/doc/devGuide/VerticalMixingCoeff.md @@ -6,8 +6,8 @@ Omega includes a `VertMix` class that provides functions that compute `VertDiff` vertical diffusivity and viscosity, the gradient Richardson number, a smoothed gradient Richardson number, where all are defined at the center of the cell and top of the layer. Currently the values of `VertDiff` and `VertVisc` are calculated using the linear combination of three options: (1) a constant background mixing value, (2) a convective instability mixing value, and (3) a Richardson -number dependent shear mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These options are linearly additive. In the future, additional additive options will be implemented, such as the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). For both the convective and shear mixing values `BruntVaisalaFreq` is needed, which -is calculated by the `EOS` class. `GradRichNum` is smoothed using a 1-2-1 filter to produce `GradRichNumSmoothed` which is used by the LMD94 shear mixing formulation. +number dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear instability driven mixing parameterization. These options are linearly additive. In the future, additional additive options will be implemented, such as the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). For both the convective and shear instability driven mixing values `BruntVaisalaFreqSq` is needed, which +is calculated by the `EOS` class. `GradRichNum` is smoothed using a 1-2-1 filter to produce `GradRichNumSmoothed` which is used by the LMD94 shear instability driven mixing formulation. ## Initialization and Usage @@ -40,7 +40,7 @@ VertMix: BaseShearValue: 0.005 RiCrit: 0.7 Exponent: 3.0 - RiSmoothLoops: 3 + RiSmoothLoops: 2 ``` ## Class Structure @@ -63,12 +63,12 @@ VertMix: - `ConvDiff`: Convective mixing coefficient (m²/s; Default: 1.0) - `ConvTriggerBVF`: Trigger threshold for convective mixing (Default: 0.0) -3. Shear Mixing: - - `EnableShearMix`: Flag to enable/disable shear mixing (Default: True) - - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) - - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) - - `ShearExponent`: Exponent parameter number for the LMB94 formulation (Default: 3.0) - - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 3) +3. Shear Instability Driven Mixing: + - `EnableShearMix`: Flag to enable/disable shear instability driven mixing (Default: True) + - `BaseShearValue`: Base values of maximum viscosity and diffusivity for the LMD94 interior shear instability driven mixing formulation (Default: 0.005) + - `RiCrit`: Critical Richardson number for the LMD94 formulation (Default: 0.7) + - `ShearExponent`: Exponent parameter number for the LMD94 formulation (Default: 3.0) + - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 2) ## Core Functionality (Vertical Mixing Coefficient Calculation) diff --git a/components/omega/doc/userGuide/VerticalMixingCoeff.md b/components/omega/doc/userGuide/VerticalMixingCoeff.md index 129033aa07f0..a66d5d28f017 100644 --- a/components/omega/doc/userGuide/VerticalMixingCoeff.md +++ b/components/omega/doc/userGuide/VerticalMixingCoeff.md @@ -6,9 +6,8 @@ The vertical mixing module in Omega handles the parameterization of unresolved v processes in the ocean. It calculates vertical diffusivity and viscosity coefficients that determine how properties (like momentum, heat, salt, and biogeochemical tracers) mix vertically in the ocean model. Currently, Omega offers three different mixing processes within the water column: (1) a constant -background mixing value, (2) a convective instability mixing value, and (3) a Richardson number -dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These are linearly additive and are describe a bit -more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear mixing calculation. +background mixing value, (2) a convective instability mixing value, and (3) a Richardson number-dependent shear-instability-driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 shear instability driven mixing parameterization. These are linearly additive and are describe a bit +more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear instability driven mixing calculation. The user-configurable options are the following parameters in the yaml configuration file: @@ -26,14 +25,14 @@ VertMix: BaseShearValue: 0.005 # Base viscosity/diffusivity value RiCrit: 0.7 # Critical Richardson number Exponent: 3.0 # Richardson number exponent - RiSmoothLoops: 3 # Number of Richardson number smoothing loops + RiSmoothLoops: 2 # Number of Richardson number smoothing loops ``` ## Vertical Mixing Processes/Types ### 1. Background Mixing -A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior. +A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior and is assumed roughly equivalent to the globally averaged interior mixing from all sources (e.g., from internal wave breaking). ### 2. Convective Mixing @@ -49,16 +48,16 @@ $$ This is different than some current implementations (i.e. in MPAS-Ocean and the CVMix package), where convective adjustment occurs both with unstable and neutral conditions ($N^2 \leq N^2_{crit}$). $\kappa_{conv}$ and $N^2_{crit}$ are constant parameters set in the `VertMix` section of the yaml file (`Diffusivity` and `TriggerBVF` under the `Convective` header). -### 3. Shear Mixing +### 3. Shear-Instability-Driven Mixing Mixing induced by vertical velocity shear, implemented using the LMD94 scheme, through the gradient Richardson number (ratio of buoyancy to shear). $$ -\nu_{shear}\,, \kappa_{shear} = = +\nu_{shear} = = \kappa_{shear} = = \begin{cases} \nu_o \quad \text{ if } Ri_g < 0\\ -\nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 < Ri_g < Ri_{crit}\\ -0.0 \quad \text{ if } Ri_{crit} < Ri_g +\nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 \leq Ri_g < Ri_{crit}\\ +0.0 \quad \text{ if } Ri_{crit} \leq Ri_g \end{cases} $$ @@ -68,4 +67,4 @@ $$ Ri = \frac{N^2}{\left|\frac{\partial \mathbf{U}}{\partial z}\right|^2}\,, $$ -where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 filter before being used to calculate the shear. $N^2$, $\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri$ of `K` are all defined at the cell center, top interface of layer `K`. $N^2$, $\nu_{shear}$ and $\kappa_{shear}$ are set to zero for the surface layer. +where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 (vertical) filter before being used to calculate the shear-instability-driven mixing. $N^2$, $\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri$ of `K` are all defined at the cell center, top interface of layer `K`. $N^2$, $\nu_{shear}$ and $\kappa_{shear}$ are set to zero at the surface. diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index b8c79472bfdc..ae8de1c76868 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -68,7 +68,7 @@ class ShearMix { Real ShearRiCrit = 0.7; ///< Critical Richardson number of LMG94 Real ShearExponent = 3.0; /// Exponent value used interior shear mixing calculation of LMG94 - I4 RiSmoothLoops = 3; ///< Number of smoothing loops for Richardson number + I4 RiSmoothLoops = 2; ///< Number of smoothing loops for Richardson number /// Constructor for ShearMix ShearMix(const VertCoord *VCoord); @@ -138,9 +138,8 @@ class GradRichardsonNum { Real GradRichNumTmp[VecLength]; for (int KVec = 0; KVec < KLen; ++KVec) { - const I4 K = KStart + KVec; - GradRichNumNorm[K] = 1.0e-12_Real; - GradRichNumTmp[K] = 100.0_Real; + GradRichNumNorm[KVec] = 1.0e-12_Real; + GradRichNumTmp[KVec] = 100.0_Real; } for (int J = 0; J < NEdgesOnCell(ICell); ++J) { @@ -151,11 +150,14 @@ class GradRichardsonNum { I4 K1 = K - 1; I4 K2 = K; - if (K <= MinLayerCell(ICell)) { + if (K < MinLayerCell(ICell) || K > MaxLayerCell(ICell)) + continue; + + if (K == MinLayerCell(ICell)) { K1 = K; K2 = K + 1; } - if (K >= MaxLayerCell(ICell)) { + if (K == MaxLayerCell(ICell)) { K1 = K - 2; K2 = K - 1; } @@ -174,15 +176,15 @@ class GradRichardsonNum { BruntVaisalaFreqSq(JCell, K2))) / (ShearSquared + 1.0e-12_Real); - Real Weight = 0.25_Real * DcEdge(JEdge) * DvEdge(JEdge); - GradRichNumNorm[K] = GradRichNumNorm[K] + Weight; - GradRichNumTmp[K] = GradRichNumTmp[K] + Weight * RiEdge; + Real Weight = 0.25_Real * DcEdge(JEdge) * DvEdge(JEdge); + GradRichNumNorm[KVec] = GradRichNumNorm[KVec] + Weight; + GradRichNumTmp[KVec] = GradRichNumTmp[KVec] + Weight * RiEdge; } } for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; - GradRichNum(ICell, K) = GradRichNumTmp[K] / GradRichNumNorm[K]; + GradRichNum(ICell, K) = GradRichNumTmp[KVec] / GradRichNumNorm[KVec]; } } diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index abd1f9a9ba8e..daacda1cea56 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -141,7 +141,7 @@ void testGradRichNum() { parallelFor( "setMinMax", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { MinLayerCell(ICell) = 0; - MaxLayerCell(ICell) = VCoord->NVertLayers; + MaxLayerCell(ICell) = VCoord->NVertLayers - 1; }); // filling CellsOnCell with simple mapping for this test From bcecf5b86a52d19ade337cc9e3d74cc77de4a3d8 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Tue, 17 Mar 2026 13:43:51 -0700 Subject: [PATCH 4/4] fixes gpu test hang --- components/omega/doc/devGuide/EOS.md | 2 +- components/omega/src/ocn/VertMix.h | 3 --- components/omega/test/ocn/VertMixTest.cpp | 12 +++++------- 3 files changed, 6 insertions(+), 11 deletions(-) diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 9aeaec4fff23..0762727a25b5 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -7,7 +7,7 @@ and `BruntVaisalaFreqSq`. Current EOS options are a linear EOS or an EOS compute 75 term expansion from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566). If `SpecVolDisplaced` is calculated with the linear EOS option, it will be equal to `SpecVol` as there is no pressure/depth dependence for the linear EOS. `SpecVolDisplaced` computes specific volume -adiabatically displaced to `K + KDisp` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreq`, as +adiabatically displaced to `K + KDisp` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreqSq`, as `SpecVol` is an input for the `BruntVaisalaFreqSq` calculation. If the linear EOS option is used, then the `BruntVaisalaFreqSq` is calculated using linear coefficients. If the TEOS-10 option is used, the `BruntVaisalaFreqSq` is calculated with non-linear coefficients according to the [TEOS-10 toolbox](https://www.teos-10.org/software.htm). Note: two assumption for ease of computation and efficiency have been made diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index ae8de1c76868..23219fcf177d 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -104,9 +104,6 @@ class ShearMix { (GradRichNumSmoothed(ICell, K) / ShearRiCrit), ShearExponent) * BaseShearValue; - } else { - VertDiff(ICell, K) += 0.0_Real; - VertVisc(ICell, K) += 0.0_Real; } } } diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index daacda1cea56..d0502605d982 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -141,7 +141,7 @@ void testGradRichNum() { parallelFor( "setMinMax", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { MinLayerCell(ICell) = 0; - MaxLayerCell(ICell) = VCoord->NVertLayers - 1; + MaxLayerCell(ICell) = NVertLayers - 1; }); // filling CellsOnCell with simple mapping for this test @@ -252,7 +252,7 @@ void testOneTwoOneFilter() { parallelFor( "setMinMax", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { MinLayerCell(ICell) = 0; - MaxLayerCell(ICell) = VCoord->NVertLayers - 1; + MaxLayerCell(ICell) = NVertLayers - 1; }); // Apply the 1-2-1 filter to each cell @@ -334,8 +334,6 @@ void testBackVertMix() { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); VCoord->NVertLayers = NVertLayers; - I4 NCellsSize = Mesh->NCellsSize; - I4 NEdgesSize = Mesh->NEdgesSize; I4 NCellsAll = Mesh->NCellsAll; I4 NEdgesAll = Mesh->NEdgesAll; OMEGA_SCOPE(ZMid, VCoord->ZMid); @@ -344,10 +342,10 @@ void testBackVertMix() { VertMix *TestVertMix = VertMix::getInstance(); /// Create and fill ocean state arrays - auto NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesSize, NVertLayers); - auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesSize, NVertLayers); + auto NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesAll, NVertLayers); + auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesAll, NVertLayers); auto BruntVaisalaFreqSqCell = - Array2DReal("BruntVaisalaFreqSqCell", NCellsSize, NVertLayers); + Array2DReal("BruntVaisalaFreqSqCell", NCellsAll, NVertLayers); /// Use deep copy initialize with reference or zero values deepCopy(NormalVelEdge, NV);