diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 745278aece02..47167f1d4272 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -86,9 +86,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: 2 IOStreams: HorzMeshIn: UsePointerFile: false diff --git a/components/omega/doc/devGuide/VerticalMixingCoeff.md b/components/omega/doc/devGuide/VerticalMixingCoeff.md index 2b3fa0a86fa2..a59151ef27cc 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 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 (vertical) filter to produce `GradRichNumSmoothed` which is used by the LMD94 shear instability driven 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: 2 ``` ## 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 @@ -60,11 +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) - - `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) +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) @@ -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 cb9b4524f49e..2f97a169d8fe 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 [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. +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: @@ -23,48 +22,49 @@ 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: 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 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 and is assumed roughly equivalent to the globally averaged interior mixing from all sources (e.g., from internal wave breaking). ### 2. Convective Mixing -Enhanced convective adjustment mixing that occurs in statically unstable regions of the water column to parameterize convection and homogenize properties. In Omega this is mixing is defaulted to occur when the squared Brunt Vaisala Frequency is less than 0.0 (unstable), and is off when equal to (neutral) or greater than (stable) 0.0. +Enhanced convective adjustment mixing that occurs in statically unstable regions of the water column to parameterize convection and homogenize properties. In Omega this mixing is defaulted to occur when the squared Brunt Vaisala Frequency is less than 0.0 (unstable), and is off when equal to (neutral) or greater than (stable) 0.0. The convective diffusivity is thus added to the background diffusivity ($\kappa_b$) to form the total diffusivity coefficient $\kappa$. $$ \kappa = \begin{cases} -+\kappa_b + \kappa_{conv} \quad \text{ if } N^2 < N^2_{crit}\\ -+\kappa_b \quad \text{ if } N^2 \geq N^2_{crit} +\kappa_b + \kappa_{conv} \quad \text{ if } N^2 < N^2_{crit}\\ +\kappa_b \quad \text{ if } N^2 \geq N^2_{crit} \end{cases} $$ 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 pseudo-velocity shear, implemented using the Pacanowski-Philander scheme, through the gradient Richardson number (ratio of buoyancy to shear). +Mixing induced by vertical pseudo-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 \leq Ri_g < Ri_{crit}\\ +0.0 \quad \text{ if } Ri_{crit} \leq 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_g$ is defined as: $$ -Ri = \frac{N^2}{\left|\frac{\partial \mathbf{U}}{\partial z}\right|^2}\,, +Ri_g = \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_g$ 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$, $\left| \frac{\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri_g$ of `k` are all defined at the cell center, top interface of layer `k`. $N^2$, $\nu$ and $\kappa$ are set to zero at the surface. diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index a0035dc9b65a..a530549d493d 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -680,9 +680,6 @@ class Teos10BruntVaisalaFreqSq { /// Linear squared Brunt-Vaisala frequency calculator class LinearBruntVaisalaFreqSq { public: - /// Coefficients from LinearEos (overwritten by config file if set there) - Real RhoT0S0 = 1000.0; ///< Reference density (kg m^-3) at (T,S)=(0,0) - /// constructor declaration LinearBruntVaisalaFreqSq(const VertCoord *VCoord); @@ -702,7 +699,7 @@ class LinearBruntVaisalaFreqSq { /// K-1 and K Do not need to use displaced specific volume here /// since only the linear EOS is used with this BVF formulation. BruntVaisalaFreqSq(ICell, K) = - -(Gravity / RhoT0S0) * + -(Gravity / RhoSw) * ((1.0_Real / SpecVol(ICell, K - 1)) - (1.0_Real / SpecVol(ICell, K))) / (GeomZMid(ICell, K - 1) - GeomZMid(ICell, K)); diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index c64fae2bafea..815dab5e6814 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -16,24 +16,40 @@ namespace OMEGA { -ShearMix::ShearMix(const HorzMesh *Mesh, const VertCoord *VCoord) - : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell), - GeomZMid(VCoord->GeomZMid), 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), GeomZMid(VCoord->GeomZMid), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), + CellsOnCell(Mesh->CellsOnCell), MinLayerCell(VCoord->MinLayerCell), + MaxLayerCell(VCoord->MaxLayerCell), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), + MaxLayerEdgeBot(VCoord->MaxLayerEdgeBot), 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) { - VertDiff = Array2DReal("VertDiff", Mesh->NCellsSize, VCoord->NVertLayers); - VertVisc = Array2DReal("VertVisc", Mesh->NCellsSize, VCoord->NVertLayers); + : ComputeVertMixConv(VCoord), ComputeVertMixShear(VCoord), + ComputeGradRichardsonNum(Mesh, VCoord), ComputeOneTwoOneFilter(VCoord), + Name(Name), Mesh(Mesh), VCoord(VCoord) { + VertDiff = Array2DReal("VertDiff", Mesh->NCellsSize, VCoord->NVertLayersP1); + VertVisc = Array2DReal("VertVisc", Mesh->NCellsSize, VCoord->NVertLayersP1); + GradRichNum = + Array2DReal("GradRichNum", Mesh->NCellsSize, VCoord->NVertLayersP1); + GradRichNumSmoothed = Array2DReal("GradRichNumSmoothed", Mesh->NCellsSize, + VCoord->NVertLayersP1); defineFields(); } @@ -54,7 +70,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 +152,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,57 +182,118 @@ 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(LocGradRichNumSmoothed, + GradRichNumSmoothed); /// 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(LocBackDiff, BackDiff); + OMEGA_SCOPE(LocBackVisc, BackVisc); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - /// Initialize VertDiff and VertVisc to background values - deepCopy(LocVertDiff, BackDiff); - deepCopy(LocVertVisc, BackVisc); - - /// Dispatch to the correct VertMix calculation - if (LocComputeVertMixShear.Enabled && LocComputeVertMixConv.Enabled) { + /// First, initialize VertDiff and VertVisc to background values + parallelForOuter( + "VertMix-BackAndRich", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int KStart = chunkStart(KChunk, KMin); + const int KLen = chunkLength(KChunk, KStart, KMax); + for (int KVec = 0; KVec < KLen; ++KVec) { + const int K = KStart + KVec; + LocVertDiff(ICell, K) = LocBackDiff; + LocVertVisc(ICell, K) = LocBackVisc; + LocGradRichNum(ICell, K) = + LocComputeGradRichardsonNum.RiInitValue; + LocGradRichNumSmoothed(ICell, K) = + LocComputeGradRichardsonNum.RiInitValue; + } + }); + }); + /// Second, compute shear mixing if enabled + if (LocComputeVertMixShear.Enabled) { + /// Compute Richardson number parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + "VertMix-ComputeRi", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); + const int KMin = MinLayerCell(ICell) + 1; const int KMax = MaxLayerCell(ICell); const int KRange = vertRangeChunked(KMin, KMax); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixShear( - LocVertDiff, LocVertVisc, ICell, KChunk, NormalVelocity, + LocComputeGradRichardsonNum( + LocGradRichNum, ICell, KChunk, NormalVelocity, TangentialVelocity, BruntVaisalaFreqSq); - LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, - KChunk, BruntVaisalaFreqSq); + }); + + teamBarrier(Team); + + // Fill Richardson number at vertical boundaries using the + // closest valid value. This is equivalent to doing one-sided + // differencing at the boundary. + Kokkos::single( + PerTeam(Team), INNER_LAMBDA() { + LocGradRichNum(ICell, MinLayerCell(ICell)) = + LocGradRichNum(ICell, KMin); + LocGradRichNum(ICell, MaxLayerCell(ICell) + 1) = + LocGradRichNum(ICell, KMax); }); }); - } else if (LocComputeVertMixShear.Enabled) { + /// Smooth Richardson number with 1-2-1 filter the number of times + /// specified by RiSmoothLoops + deepCopy(LocGradRichNumSmoothed, LocGradRichNum); + for (int SmoothLoop = 0; + SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { + parallelForOuter( + "VertMix-RiSmooth", {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) { + LocOneTwoOneFilter(LocGradRichNumSmoothed, ICell, KChunk, + LocGradRichNumSmoothed); + }); + }); + } + /// Compute shear mixing using smoothed Richardson number parallelForOuter( - "VertMix-ShearOnly", {Mesh->NCellsAll}, + "VertMix-Shear", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); + const int KMin = MinLayerCell(ICell) + 1; const int KMax = MaxLayerCell(ICell); const int KRange = vertRangeChunked(KMin, KMax); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixShear( - LocVertDiff, LocVertVisc, ICell, KChunk, NormalVelocity, - TangentialVelocity, BruntVaisalaFreqSq); + LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, + KChunk, LocGradRichNumSmoothed); }); }); - } else if (LocComputeVertMixConv.Enabled) { + } + /// Third, compute convective mixing if 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 KMin = MinLayerCell(ICell) + 1; const int KMax = MaxLayerCell(ICell); const int KRange = vertRangeChunked(KMin, KMax); @@ -221,24 +303,32 @@ 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; - }); } + /// Finally, zero viscosity/diffusivity at surface and bottom boundaries + parallelFor( + "VertMix-Boundaries", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; + LocVertDiff(ICell, KMin) = 0.0_Real; + LocVertVisc(ICell, KMin) = 0.0_Real; + LocVertDiff(ICell, KMax) = 0.0_Real; + LocVertVisc(ICell, KMax) = 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 @@ -246,7 +336,7 @@ void VertMix::defineFields() { int NDims = 2; std::vector DimNames(NDims); DimNames[0] = "NCells"; - DimNames[1] = "NVertLayers"; + DimNames[1] = "NVertLayersP1"; /// Create and register the Diffusivity field auto VertDiffField = @@ -274,6 +364,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 +399,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 261ecd5ca8b6..9466eb67d7ec 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -38,19 +38,15 @@ class ConvectiveMix { operator()(Array2DReal VertDiff, Array2DReal VertVisc, I4 ICell, I4 KChunk, const Array2DReal &BruntVaisalaFreqSq) const { - const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell)); + const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1); const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); 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 (BruntVaisalaFreqSq(ICell, K) < ConvTriggerBVF) { + VertDiff(ICell, K) += ConvDiff; + VertVisc(ICell, K) += ConvDiff; } } } @@ -65,68 +61,169 @@ 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 LMD94 + Real ShearRiCrit = 0.7; ///< Critical Richardson number of LMD94 Real ShearExponent = - 2.0; /// Exponent value used in Pacanowski and Philander (1981) Eqs (1). + 3.0; /// Exponent value used interior shear mixing calculation of LMD94 + I4 RiSmoothLoops = 2; ///< 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 &GradRichNumSmoothed) const { + + const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1); + const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + + 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; + } + } + } + + private: + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + +/// Class for Gradient Richardson Number calculation +class GradRichardsonNum { + public: + Real RiInitValue = 100.0_Real; ///< Initial Richardson number value + + /// 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 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1); 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; - if (K == 0) { - 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); - } - Real DelZMid = GeomZMid(ICell, K - 1) - GeomZMid(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); + GradRichNumNorm[KVec] = 1.0e-12_Real; + GradRichNumTmp[KVec] = RiInitValue; + } + + 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; + I4 K1 = K - 1; + I4 K2 = K; + + // Skip this edge contribution if it would access + // invalid edge velocity levels. + if (K1 > MaxLayerEdgeBot(JEdge) || K2 > MaxLayerEdgeBot(JEdge)) + continue; + + Real DNormVel = + NormalVelocity(JEdge, K1) - NormalVelocity(JEdge, K2); + Real DTanVel = + TangentialVelocity(JEdge, K1) - TangentialVelocity(JEdge, K2); + Real DzEdge = + 0.5_Real * (GeomZMid(ICell, K1) + GeomZMid(JCell, K1) - + (GeomZMid(ICell, K2) + GeomZMid(JCell, K2))); + Real ShearSquared = + (DNormVel * DNormVel + DTanVel * DTanVel) / (DzEdge * DzEdge); + Real RiEdge = + Kokkos::max(0.0_Real, + 0.5_Real * (BruntVaisalaFreqSq(ICell, K2) + + BruntVaisalaFreqSq(JCell, K2))) / + (ShearSquared + 1.0e-12_Real); + + Real Weight = 0.25_Real * DcEdge(JEdge) * DvEdge(JEdge); + GradRichNumNorm[KVec] += Weight; + GradRichNumTmp[KVec] += Weight * RiEdge; } } + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + GradRichNum(ICell, K) = GradRichNumTmp[KVec] / GradRichNumNorm[KVec]; + } } private: - Array1DReal DcEdge; - Array1DReal DvEdge; - Array1DReal AreaCell; Array2DReal GeomZMid; - Array1DI4 NEdgesOnCell; Array2DI4 EdgesOnCell; + Array2DI4 CellsOnCell; + Array2DI4 CellsOnEdge; + Array1DI4 NEdgesOnCell; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; + Array1DI4 MaxLayerEdgeBot; + 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) + 1); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + if (K > MinLayerCell(ICell) && K < MaxLayerCell(ICell)) { + // 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); + } + } + } + + private: Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; }; @@ -140,13 +237,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 +259,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/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 2e54a23306b7..0e7eb2c5be76 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -49,7 +49,7 @@ const Real TeosBVFExpValue = 0.020913834194283325; // Expected value for TEOS-10 squared Brunt-Vaisala // frequency const Real LinearBVFExpValue = - 0.017834796542017275; // Expected value for Linear squared Brunt-Vaisala + 0.017382842633545097; // Expected value for Linear squared Brunt-Vaisala // frequency const Real GswBVFExpValue = 0.02081197958166906; // Expected value from GSW-C library diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index 98ad51d04051..d1bf97015798 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -28,33 +28,35 @@ using namespace OMEGA; /// Test constants and expected values -constexpr int NVertLayers = 60; +constexpr int NVertLayers = 60; +constexpr int NVertLayersP1 = NVertLayers + 1; /// 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. @@ -110,14 +112,264 @@ 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(); + auto *MeshHalo = Halo::getDefault(); + VCoord->NVertLayers = NVertLayers; + VCoord->NVertLayersP1 = NVertLayersP1; + I4 NCellsSize = Mesh->NCellsSize; + I4 NEdgesAll = Mesh->NEdgesAll; + OMEGA_SCOPE(GeomZMid, VCoord->GeomZMid); + OMEGA_SCOPE(NEdgesOnCell, Mesh->NEdgesOnCell); + OMEGA_SCOPE(EdgesOnCell, Mesh->EdgesOnCell); + OMEGA_SCOPE(AreaCell, Mesh->AreaCell); + OMEGA_SCOPE(DcEdge, Mesh->DcEdge); + OMEGA_SCOPE(DvEdge, Mesh->DvEdge); + OMEGA_SCOPE(CellsOnCell, Mesh->CellsOnCell); + 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 NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesAll, NVertLayers); + auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesAll, NVertLayers); + auto BruntVaisalaFreqSqCell = + Array2DReal("BruntVaisalaFreqSqCell", NCellsSize, NVertLayersP1); + /// Use deep copy to initialize results + deepCopy(NormalVelEdge, NV); + deepCopy(TangVelEdge, TV); + deepCopy(BruntVaisalaFreqSqCell, BVFP); + deepCopy(TestVertMix->GradRichNum, 0.0); + + parallelFor( + "populateArrays", {Mesh->NCellsAll, NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + GeomZMid(ICell, K) = -K; + NEdgesOnCell(ICell) = 5; + AreaCell(ICell) = 3.6e10_Real; + }); + + // Also test guard is working for skipping layer edge contribution + // if it would access invalid edge velocity levels + parallelFor( + "setMinMax", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + MinLayerCell(ICell) = 0; + if (ICell % 2 == 0) { + MaxLayerCell(ICell) = 10; + } else { + MaxLayerCell(ICell) = NVertLayers - 2; + } + }); + + // Refresh edge layer ranges after overriding MaxLayerCell. + VCoord->minMaxLayerEdge(MeshHalo); + + // Rebind GradRichardsonNum so the functor captures refreshed edge ranges. + TestVertMix->ComputeGradRichardsonNum = GradRichardsonNum(Mesh, VCoord); + + // Recapture after minMaxLayerEdge since it reallocates edge layer views. + OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot); + + // filling CellsOnCell with simple mapping for this test + parallelFor( + "populateArrays", {Mesh->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) { + if (K > MaxLayerEdgeBot(IEdge)) { + NormalVelEdge(IEdge, K) = + -9.99e30; // Fill value to cause failure if used + TangVelEdge(IEdge, K) = + -9.99e30; // Fill value to cause failure if used + } else { + 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); + + /// 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) + 1; + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + + // Match the production K1/K2 logic to determine whether + // this layer has any valid edge contributions. + I4 K1 = K - 1; + I4 K2 = K; + + if (K == MaxLayerCell(ICell) + 1) { + K1 = K - 2; + K2 = K - 1; + } + + bool HasValidEdge = false; + for (int J = 0; J < NEdgesOnCell(ICell); ++J) { + const I4 JEdge = EdgesOnCell(ICell, J); + if (K1 <= MaxLayerEdgeBot(JEdge) && + K2 <= MaxLayerEdgeBot(JEdge)) { + HasValidEdge = true; + break; + } + } + + if (HasValidEdge) { + if (!isApprox(GradRichNum(ICell, K), RiExpValue, RTol)) + InnerCount++; + } else { + // With all edges skipped, GradRichNum remains at sentinel + // scale (~1e14) from initialization in the functor. + if (!(GradRichNum(ICell, K) > 1.0e10_Real)) + InnerCount++; + } + }, + NumMismatchesCol); + + Kokkos::single(PerTeam(Team), + [&]() { OuterCount += NumMismatchesCol; }); + }, + NumMismatches); + + // If test fails, print bad values and abort + if (NumMismatches != 0) { + 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; + VCoord->NVertLayersP1 = NVertLayersP1; + I4 NCellsSize = Mesh->NCellsSize; + I4 NChunks = VCoord->NVertLayers / VecLength; + 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", NCellsSize, NVertLayersP1); + auto GradRichNum = Array2DReal("GradRichNum", NCellsSize, NVertLayersP1); + /// 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", {Mesh->NCellsAll, NVertLayersP1}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + if (K % 2 == 0) { + GradRichNum(ICell, K) = 1.0; + } else { + GradRichNum(ICell, K) = -1.0; + } + }); + + parallelFor( + "setMinMax", {Mesh->NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + MinLayerCell(ICell) = 0; + MaxLayerCell(ICell) = NVertLayers - 1; + }); + + // Apply the 1-2-1 filter to each cell + OMEGA_SCOPE(ComputeOneTwoOneFilter, TestVertMix->ComputeOneTwoOneFilter); + parallelFor( + "ApplyOneTwoOneFilter", {Mesh->NCellsAll, NChunks}, + KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { + 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) + 1; + const int KRange = vertRange(KMin, KMax); + parallelReduceInner( + Team, KRange, + INNER_LAMBDA(int KOff, int &InnerCount) { + const int K = KMin + KOff; + if (K > MinLayerCell(ICell) && K < MaxLayerCell(ICell)) { + // 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) 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) { + 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(); - const auto VCoord = VertCoord::getDefault(); - VCoord->NVertLayers = NVertLayers; - I4 NCellsSize = Mesh->NCellsSize; - I4 NEdgesSize = Mesh->NEdgesSize; - I4 NEdgesAll = Mesh->NEdgesAll; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + VCoord->NVertLayersP1 = NVertLayersP1; + I4 NCellsSize = Mesh->NCellsSize; + I4 NEdgesSize = Mesh->NEdgesSize; + I4 NEdgesAll = Mesh->NEdgesAll; OMEGA_SCOPE(GeomZMid, VCoord->GeomZMid); /// Get VertMix instance to test @@ -127,7 +379,7 @@ void testBackVertMix() { auto NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesSize, NVertLayers); auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesSize, NVertLayers); auto BruntVaisalaFreqSqCell = - Array2DReal("BruntVaisalaFreqSqCell", NCellsSize, NVertLayers); + Array2DReal("BruntVaisalaFreqSqCell", NCellsSize, NVertLayersP1); /// Use deep copy initialize with reference or zero values deepCopy(NormalVelEdge, NV); @@ -161,21 +413,21 @@ 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}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { - // Surface layer should be zero + if (K == KMin || K == KMax) { + // Surface and bottom layers should be zero if (BackVertVisc(ICell, K) != 0.0_Real) InnerCount++; } else { @@ -191,26 +443,29 @@ 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}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { - // Surface layer should be zero + if (K == KMin || K == KMax) { + // Surface and bottom layers should be zero if (BackVertDiff(ICell, K) != 0.0_Real) InnerCount++; } else { @@ -231,82 +486,84 @@ 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(); - VCoord->NVertLayers = NVertLayers; - I4 NCellsSize = Mesh->NCellsSize; - I4 NEdgesAll = Mesh->NEdgesAll; - OMEGA_SCOPE(GeomZMid, VCoord->GeomZMid); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + VCoord->NVertLayersP1 = NVertLayersP1; + I4 NCellsSize = Mesh->NCellsSize; + 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", NCellsSize, NVertLayers); + auto BruntVaisalaFreqSqIn = + Array2DReal("BruntVaisalaFreqSqIn", NCellsSize, NVertLayersP1); + auto VertDiffOut = + Array2DReal("VertDiffOut", Mesh->NCellsAll, NVertLayersP1); + auto VertViscOut = + Array2DReal("VertViscOut", Mesh->NCellsAll, NVertLayersP1); /// 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", {Mesh->NCellsAll, NVertLayers}, - KOKKOS_LAMBDA(I4 ICell, I4 K) { GeomZMid(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; + "populateArrays", {Mesh->NCellsAll, NVertLayersP1}, + 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); + OMEGA_SCOPE(ComputeVertMixConv, TestVertMix->ComputeVertMixConv); + parallelFor( + "ApplyVertMixConv", {Mesh->NCellsAll, NChunks}, + KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { + 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}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { - // Surface layer should be zero - if (ConvVertVisc(ICell, K) != 0.0_Real) + if (K == KMin || K == KMax) { + // Surface and bottom layers should be zero + 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++; } }, @@ -318,32 +575,34 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ConvVertViscH = createHostMirrorCopy(ConvVertVisc); - ABORT_ERROR("TestVertMixConv: VertVisc FAIL, " - "expected {}, got {} with {} mismatches", - VertDiffConvExp, ConvVertViscH(1, 1), NumMismatches); + 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}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { - // Surface layer should be zero - if (ConvVertDiff(ICell, K) != 0.0_Real) + if (K == KMin || K == KMax) { + // Surface and bottom layers should be zero + 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++; } }, @@ -355,10 +614,10 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ConvVertDiffH = createHostMirrorCopy(ConvVertDiff); - ABORT_ERROR("TestVertMixConv: VertDiff FAIL, " - "expected {}, got {} with {} mismatches", - VertDiffConvExp, ConvVertDiffH(1, 1), NumMismatches); + ABORT_ERROR("TestVertMixConv: VertDiff FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMixConv: VertDiff PASS"); } return; @@ -366,90 +625,83 @@ void testConvVertMix() { void testShearVertMix() { /// Get mesh and coordinate info - const auto Mesh = HorzMesh::getDefault(); - const auto VCoord = VertCoord::getDefault(); - VCoord->NVertLayers = NVertLayers; - I4 NCellsSize = Mesh->NCellsSize; - I4 NEdgesAll = Mesh->NEdgesAll; - OMEGA_SCOPE(GeomZMid, VCoord->GeomZMid); - OMEGA_SCOPE(NEdgesOnCell, Mesh->NEdgesOnCell); - OMEGA_SCOPE(AreaCell, Mesh->AreaCell); - OMEGA_SCOPE(DcEdge, Mesh->DcEdge); - OMEGA_SCOPE(DvEdge, Mesh->DvEdge); + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + VCoord->NVertLayersP1 = NVertLayersP1; + I4 NCellsSize = Mesh->NCellsSize; + 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", NCellsSize, NVertLayers); + auto GradRichNumSmoothedIn = + Array2DReal("GradRichNumSmoothedIn", NCellsSize, NVertLayersP1); + auto VertDiffOut = + Array2DReal("VertDiffOut", Mesh->NCellsAll, NVertLayersP1); + auto VertViscOut = + Array2DReal("VertViscOut", Mesh->NCellsAll, NVertLayersP1); /// 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", {Mesh->NCellsAll, NVertLayers}, + "populateArrays", {Mesh->NCellsAll, NVertLayersP1}, KOKKOS_LAMBDA(I4 ICell, I4 K) { - GeomZMid(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 + OMEGA_SCOPE(ComputeVertMixShear, TestVertMix->ComputeVertMixShear); + 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", {Mesh->NCellsAll, NChunks}, + KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { + 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}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { - if (ShearVertVisc(ICell, K) != 0.0_Real) + if (K == KMin || K == KMax) { + // Surface and bottom layers should be zero + 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++; } }, @@ -461,38 +713,38 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ShearVertViscH = createHostMirrorCopy(ShearVertVisc); - ABORT_ERROR("TestVertMixShear: VertVisc FAIL, " - "expected {}, got {} with {} mismatches", - VertViscShearExp, ShearVertViscH(1, 1), NumMismatches); + 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}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { - if (ShearVertDiff(ICell, K) != 0.0_Real) + if (K == KMin || K == KMax) { + // Surface and bottom layers should be zero + 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++; } }, @@ -504,10 +756,10 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { - auto ShearVertDiffH = createHostMirrorCopy(ShearVertDiff); - ABORT_ERROR("TestVertMixShear: VertDiff FAIL, " - "expected {}, got {} with {} mismatches", - VertDiffShearExp, ShearVertDiffH(1, 1), NumMismatches); + ABORT_ERROR("TestVertMixShear: VertDiff FAIL with {} bad values", + NumMismatches); + } else { + LOG_INFO("TestVertMixShear: VertDiff PASS"); } return; @@ -516,16 +768,18 @@ void testShearVertMix() { /// Test vertical mixing coefficients calculation for all cells/layers void testTotalVertMix() { /// Get mesh and coordinate info - const auto Mesh = HorzMesh::getDefault(); - const auto VCoord = VertCoord::getDefault(); - VCoord->NVertLayers = NVertLayers; - I4 NCellsSize = Mesh->NCellsSize; - I4 NEdgesAll = Mesh->NEdgesAll; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + VCoord->NVertLayersP1 = NVertLayersP1; + I4 NCellsSize = Mesh->NCellsSize; + I4 NEdgesAll = Mesh->NEdgesAll; OMEGA_SCOPE(GeomZMid, VCoord->GeomZMid); 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(); @@ -534,7 +788,7 @@ void testTotalVertMix() { auto NormalVelEdge = Array2DReal("NormalVelEdge", NEdgesAll, NVertLayers); auto TangVelEdge = Array2DReal("TangVelEdge", NEdgesAll, NVertLayers); auto BruntVaisalaFreqSqCell = - Array2DReal("BruntVaisalaFreqSqCell", NCellsSize, NVertLayers); + Array2DReal("BruntVaisalaFreqSqCell", NCellsSize, NVertLayersP1); /// Use deep copy to initialize with the ref value deepCopy(NormalVelEdge, NV); @@ -544,6 +798,7 @@ void testTotalVertMix() { deepCopy(BruntVaisalaFreqSqCell, BVFP); deepCopy(TestVertMix->VertDiff, 0.0); deepCopy(TestVertMix->VertVisc, 0.0); + deepCopy(TestVertMix->GradRichNumSmoothed, 0.0); parallelFor( "populateArrays", {Mesh->NCellsAll, NVertLayers}, @@ -553,6 +808,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", {Mesh->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) { @@ -583,17 +849,18 @@ void testTotalVertMix() { KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { + if (K == KMin || K == KMax) { + // Surface and bottom layers should be zero if (VertDiffP(ICell, K) != 0.0_Real) InnerCount++; // K = 1 should have ref value - } else if (K == 1) { + } else if (K == KMin + 1) { if (!isApprox(VertDiffP(ICell, K), VertDiffExpValueP, RTol)) InnerCount++; // otherwise check for invalid values @@ -616,6 +883,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 @@ -625,17 +894,17 @@ void testTotalVertMix() { KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { + if (K == KMin || K == KMax) { if (VertViscP(ICell, K) != 0.0_Real) InnerCount++; // K = 1 should have ref value - } else if (K == 1) { + } else if (K == KMin + 1) { if (!isApprox(VertViscP(ICell, K), VertViscExpValueP, RTol)) InnerCount++; // otherwise check for invalid values @@ -658,6 +927,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 @@ -678,17 +949,17 @@ void testTotalVertMix() { KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { + if (K == KMin || K == KMax) { if (VertDiffN(ICell, K) != 0.0_Real) InnerCount++; // K = 1 should have ref value - } else if (K == 1) { + } else if (K == KMin + 1) { if (!isApprox(VertDiffN(ICell, K), VertDiffExpValueN, RTol)) InnerCount++; // otherwise check for invalid values @@ -711,6 +982,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 @@ -720,17 +993,17 @@ void testTotalVertMix() { KOKKOS_LAMBDA(int ICell, const TeamMember &Team, int &OuterCount) { int NumMismatchesCol; const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); + const int KMax = MaxLayerCell(ICell) + 1; const int KRange = vertRange(KMin, KMax); parallelReduceInner( Team, KRange, INNER_LAMBDA(int KOff, int &InnerCount) { const int K = KMin + KOff; - if (K == 0) { + if (K == KMin || K == KMax) { if (VertViscN(ICell, K) != 0.0_Real) InnerCount++; // K = 1 should have ref value - } else if (K == 1) { + } else if (K == KMin + 1) { if (!isApprox(VertViscN(ICell, K), VertViscExpValueN, RTol)) InnerCount++; // otherwise check for invalid values @@ -753,6 +1026,8 @@ void testTotalVertMix() { ABORT_ERROR("TestVertMix: VertViscNegative FAIL, " "expected {}, got {} with {} mismatches", VertViscExpValueN, VertViscNH(1, 1), NumMismatches); + } else { + LOG_INFO("TestVertMixTotal: VertViscNeg PASS"); } return; @@ -771,12 +1046,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() { @@ -785,6 +1062,8 @@ void vertMixTest() { initVertMixTest(); // test each vertical mix option + testGradRichNum(); + testOneTwoOneFilter(); testBackVertMix(); testConvVertMix(); testShearVertMix();