diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index b2b1a7c37b16..4bc3703f0b8a 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -59,6 +59,19 @@ Omega: DRhoDT: -0.2 DRhoDS: 0.8 RhoT0S0: 1000.0 + VertMix: + Background: + Diffusivity: 1.0e-5 + Viscosity: 1.0e-4 + Convective: + Enable: true + Diffusivity: 1.0 + TriggerBVF: 0.0 + Shear: + Enable: true + NuZero: 0.005 + Alpha: 5.0 + Exponent: 2.0 IOStreams: InitialVertCoord: UsePointerFile: false diff --git a/components/omega/doc/design/images/ocean.jpg b/components/omega/doc/design/images/ocean.jpg new file mode 100644 index 000000000000..70607c7b7a41 Binary files /dev/null and b/components/omega/doc/design/images/ocean.jpg differ diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 2cef706e9382..015308b04bac 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -2,12 +2,13 @@ # Equation of State (EOS) -Omega includes an `Eos` class that provides functions that compute `SpecVol` and `SpecVolDisplaced`. -Current EOS options are a linear EOS or an EOS computed using the TEOS-10 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 +Omega includes an `Eos` class that provides functions that compute `SpecVol`, `SpecVolDisplaced`, +and `BruntVaisalaFreq`. Current EOS options are a linear EOS or an EOS computed using the TEOS-10 +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`. +adiabatically displaced to `K + KDisp`. Note: `SpecVol` must be calculated before `BruntVaisalaFreq`, as +`SpecVol` is an input for the `BruntVaisalaFreq` calculation. ## Eos type @@ -49,9 +50,16 @@ and pressure arrays and displaced vertical index level, do Eos.computeSpecVolDisp(ConsrvTemp, AbsSalinity, Pressure, KDisp); ``` -where `KDisp` is the number of `k` levels you want to displace each specific volume level to. +where `KDisp` is the number of `k` layers you want to displace each specific volume layer to. For example, to displace each level to one below, set `KDisp = 1`. +To compute `BruntVaisalaFreq` for a particular set of temperature, salinity, pressure, and specific +volume arrays, do + +```c++ +Eos.computeBruntVaisalaFreq(ConservTemp, AbsSalinity, Pressure, SpecVol); +``` + ## Removal of Eos To clear the Eos instance do: diff --git a/components/omega/doc/devGuide/VerticalMixingCoeff.md b/components/omega/doc/devGuide/VerticalMixingCoeff.md new file mode 100644 index 000000000000..40639c9cb882 --- /dev/null +++ b/components/omega/doc/devGuide/VerticalMixingCoeff.md @@ -0,0 +1,86 @@ +(omega-dev-vertmix) = + +# 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 top of cell centers. 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 `BruntVaisalaFreq` is needed, which +is calculated by the `EOS` class. + +## Initialization and Usage + +The primary class `VertMix` is implemented using the Singleton pattern to ensure a single instance manages all vertical mixing calculations. + +```c++ +// Initialize VertMix +VertMix* vertMix = VertMix::getInstance("Default", mesh, nLevels); + +// Compute mixing coefficients +vertMix->computeVertMix(VertDiff, VertVisc, NormalVelocity, + TangentialVelocity, BruntVaisalaFreq); +``` + +## Configuration + +The initialization process reads parameters from the yaml configuration file with the following structure and +default values: + +```yaml +VertMix: + Background: + Viscosity: 1e-4 + Diffusivity: 1e-5 + Convective: + Enable: true + Diffusivity: 1.0 + TriggerBVF: 0.0 + Shear: + Enable: true + NuZero: 0.005 + Alpha: 5.0 + Exponent: 2.0 +``` + +## Class Structure + +### Core Data Members + +- `VertDiff`: 2D array storing vertical diffusivity coefficients (m²/s) +- `VertVisc`: 2D array storing vertical viscosity coefficients (m²/s) + +### Mixing Parameters + +1. Background Mixing: + - `BackDiff`: Background vertical diffusivity (m²/s; Default: 1e-5) + - `BackVisc`: Background vertical viscosity (m²/s: Default: 1e-4) + +2. Convective Mixing: + - `EnableConvMix`: Flag to enable/disable convective mixing (Default: True) + - `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) + +## Core Functionality (Vertical Mixing Coefficient Calculation) + +The main computation is handled by: + +```cpp +void computeVertMix(Array2DReal VertDiff, + Array2DReal VertVisc, + const Array2DReal &NormalVelocity, + const Array2DReal &TangentialVelocity, + const Array2DReal &BruntVaisalaFreq); +``` + +This method combines the effects of: +- Background mixing (constant coefficients) +- Convective mixing (triggered by static instability) +- Shear instability driven mixing (Pacanowski-Philander scheme) diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index c7469d2ff19b..e58f9739dbce 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -50,6 +50,7 @@ userGuide/Tracers userGuide/TridiagonalSolvers userGuide/VertCoord userGuide/Timing +userGuide/VerticalMixingCoeff ``` ```{toctree} @@ -90,6 +91,7 @@ devGuide/Tracers devGuide/TridiagonalSolvers devGuide/VertCoord devGuide/Timing +devGuide/VerticalMixingCoeff ``` ```{toctree} diff --git a/components/omega/doc/userGuide/EOS.md b/components/omega/doc/userGuide/EOS.md index 3a7ca6b40097..1816d26b7b2b 100644 --- a/components/omega/doc/userGuide/EOS.md +++ b/components/omega/doc/userGuide/EOS.md @@ -4,7 +4,7 @@ The equation of state (EOS) for the ocean describes the relationship between specific volume of seawater (in $\textrm{m}^3/\textrm{kg}$; the reciprocal of density) and temperature (in $^{\circ}\textrm{C}$), salinity (in $\textrm{g/kg}$), and pressure (in $\textrm{dbar}$). Through the hydrostatic balance (which relates density/specific volume gradients to pressure gradients), the equation of state provides a connection between active tracers (temperature and salinity) and the fluid dynamics. -Two choices of EOS are provided by Omega: a linear EOS and a TEOS-10 EOS. The linear EOS simplifies the relationship by excluding the influence of pressure and using constant expansion/contraction coefficients, making the specific volume a simple linear function of temperature and salinity. However, this option is only recommended for simpler idealized test cases as its accuracy is not sufficient for real ocean simulations. The TEOS-10 EOS is a 75-term polynomial expression from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566) that approximates the [Thermodynamic Equation of Seawater 2010](https://www.teos-10.org/pubs/TEOS-10_Manual.pdf) , but in a less complex and more computationally efficient manner, and is the preferred EOS for real ocean simulations in Omega. +Two choices of EOS are provided by Omega: a linear EOS and a TEOS-10 EOS. The linear EOS simplifies the relationship by excluding the influence of pressure and using constant expansion/contraction coefficients, making the specific volume a simple linear function of temperature and salinity. However, this option is only recommended for simpler idealized test cases as its accuracy is not sufficient for real ocean simulations. The TEOS-10 EOS is a 75-term polynomial expression from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566) that approximates the [Thermodynamic Equation of Seawater 2010](https://www.teos-10.org/pubs/TEOS-10_Manual.pdf), but in a less complex and more computationally efficient manner, and is the preferred EOS for real ocean simulations in Omega. The user-configurable options are: `EosType` (choose either `Linear` or `Teos-10`), as well as the parameters needed for the linear EOS. @@ -19,4 +19,12 @@ Eos: where `DRhoDT` is the thermal expansion coefficient ($\textrm{kg}/(\textrm{m}^3 \cdot ^{\circ}\textrm{C})$), `DRhoDS` is the saline contraction coefficient ($\textrm{kg}/\textrm{m}^3$), and `RhoT0S0` is the reference density at (T,S)=(0,0) (in $\textrm{kg}/\textrm{m}^3$). -In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` is also calculated by the EOS. This calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels, capturing the effects of pressure/depth changes. This is primarily used to calculate quantities for determining the water column stability (i.e. the stratification) and the vertical mixing coefficients (viscosity and diffusivity). Note: when using the linear EOS, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth. +In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` and `BruntVaisalaFreq` are also calculated by the EOS. + +## Displaced Specific Volume + +The `Eos` class calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels, capturing the effects of pressure/depth changes. This is primarily used to calculate quantities for determining the water column stability (i.e. the stratification) and the vertical mixing coefficients (viscosity and diffusivity). Note: when using the `Linear` EOS option, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth. + +## Brunt Vaisala Frequency + +The `Eos` class also calculates the Brunt Vaisala Frequency, which is used by the `VertMix` class to determine water column stability for both the convective adjustment and shear-instability driven mixing. When using the `Linear` EOS option, the `BruntVaisalaFreq` is calculated using with linear coefficients. When using the `Teos-10` EOS option, the `BruntVaisalaFreq` is calculated with non-linear coefficients according to the TOES-10. For both options, `SpecVol` must be calculated prior to calculating `BruntVaisalaFreq`, as it is an input to the `BruntVaisalaFreq` calculation. diff --git a/components/omega/doc/userGuide/VerticalMixingCoeff.md b/components/omega/doc/userGuide/VerticalMixingCoeff.md new file mode 100644 index 000000000000..cb2410386dbd --- /dev/null +++ b/components/omega/doc/userGuide/VerticalMixingCoeff.md @@ -0,0 +1,68 @@ +(omega-user-vertmix)= + +# Vertical Mixing Coefficients + +The vertical mixing module in Omega handles the parameterization of unresolved vertical mixing +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. + +The user-configurable options are the following parameters in the yaml configuration file: + +```yaml +VertMix: + Background: + Viscosity: 1.0e-4 # Background vertical viscosity (m²/s) + Diffusivity: 1.0e-5 # Background vertical diffusivity (m²/s) + Convective: + Enable: true # Enables the convective-induced mixing option + Diffusivity: 1.0 # Convective mixing coefficient (m²/s) + TriggerBVF: -1e-4 # 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 +``` + +## 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. + +### 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 Brunt Vaisala Frequency is less than 0.0 (unstable), and is off when equal to (neutral) or greater than (stable) 0.0. + +$$ +\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} +\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}$). + +### 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). + +$$ +\nu = \frac{\nu_o}{(1+\alpha Ri)^n} + \nu_b\,, +$$ + +$$ +\kappa = \frac{\nu}{(1+\alpha Ri)} + \kappa_b\,. +$$ + +where $Ri$ is defined as: + +$$ +Ri = \frac{N^2}{\left|\frac{\partial \mathbf{U}}{\partial z}\right|^2}\,, +$$ diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index b8d1751c24c3..d6669a038b1d 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -14,25 +14,36 @@ namespace OMEGA { /// Constructor for Teos10Eos -Teos10Eos::Teos10Eos(int NVertLayers) : NVertLayers(NVertLayers) { - SpecVolPCoeffs = Array2DReal("SpecVolPCoeffs", 6, VecLength); -} +Teos10Eos::Teos10Eos(const VertCoord *VCoord) + : NVertLayers(VCoord->NVertLayers) {} /// Constructor for LinearEos -LinearEos::LinearEos() {} +LinearEos::LinearEos(const VertCoord *VCoord) + : NVertLayers(VCoord->NVertLayers) {} + +/// Constructor for Teos10 Brunt-Vaisala frequency +Teos10BruntVaisalaFreq::Teos10BruntVaisalaFreq(const VertCoord *VCoord) + : NVertLayers(VCoord->NVertLayers) {} + +/// Constructor for Linear Brunt-Vaisala frequency +LinearBruntVaisalaFreq::LinearBruntVaisalaFreq(const VertCoord *VCoord) + : NVertLayers(VCoord->NVertLayers), ZMid(VCoord->ZMid) {} /// Constructor for Eos Eos::Eos(const std::string &Name_, ///< [in] Name for eos object const HorzMesh *Mesh, ///< [in] Horizontal mesh - int NVertLayers ///< [in] Number of vertical layers + const VertCoord *VCoord ///< [in] Vertical coordinate ) - : ComputeSpecVolTeos10(NVertLayers) { - SpecVol = Array2DReal("SpecVol", Mesh->NCellsAll, NVertLayers); + : ComputeSpecVolLinear(VCoord), ComputeBruntVaisalaFreqLinear(VCoord), + ComputeBruntVaisalaFreqTeos10(VCoord), ComputeSpecVolTeos10(VCoord) { + SpecVol = Array2DReal("SpecVol", Mesh->NCellsAll, VCoord->NVertLayers); SpecVolDisplaced = - Array2DReal("SpecVolDisplaced", Mesh->NCellsAll, NVertLayers); + Array2DReal("SpecVolDisplaced", Mesh->NCellsAll, VCoord->NVertLayers); + BruntVaisalaFreq = + Array2DReal("BruntVaisalaFreq", Mesh->NCellsAll, VCoord->NVertLayers); // Array dimension lengths NCellsAll = Mesh->NCellsAll; - NChunks = NVertLayers / VecLength; + NChunks = VCoord->NVertLayers / VecLength; Name = Name_; defineFields(); @@ -60,8 +71,8 @@ void Eos::destroyInstance() { void Eos::init() { if (!Instance) { - Instance = new Eos("Default", HorzMesh::getDefault(), - VertCoord::getDefault()->NVertLayers); + Instance = + new Eos("Default", HorzMesh::getDefault(), VertCoord::getDefault()); } Error Err; // error code @@ -176,19 +187,58 @@ void Eos::computeSpecVolDisp(const Array2DReal &ConservTemp, } } +/// Compute Brunt-Vaisala frequency (NSquared) for all cells/layers +void Eos::computeBruntVaisalaFreq(const Array2DReal &ConservTemp, + const Array2DReal &AbsSalinity, + const Array2DReal &Pressure, + const Array2DReal &SpecVol) { + OMEGA_SCOPE(LocBruntVaisalaFreq, + BruntVaisalaFreq); /// Local view for computation + OMEGA_SCOPE( + LocComputeBruntVaisalaFreqLinear, + ComputeBruntVaisalaFreqLinear); /// Local view for linear computation + OMEGA_SCOPE( + LocComputeBruntVaisalaFreqTeos10, + ComputeBruntVaisalaFreqTeos10); /// Local view for TEOS-10 computation + deepCopy(LocBruntVaisalaFreq, + 0); /// Initialize local Brunt-Vaisala frequency to zero + + /// Dispatch to the correct Brunt-Vaisala frequency calculation + if (EosChoice == EosType::LinearEos) { + /// If Linear EOS, use linear Brunt-Vaisala frequency calculation + parallelFor( + "bvf-linear", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + LocComputeBruntVaisalaFreqLinear(LocBruntVaisalaFreq, ICell, + SpecVol); + }); + } else if (EosChoice == EosType::Teos10Eos) { + /// If TEOS-10 EOS, use TEOS-10 Brunt-Vaisala frequency calculation + parallelFor( + "bvf-teos10", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + /// Compute Brunt-Vaisala frequency + LocComputeBruntVaisalaFreqTeos10(LocBruntVaisalaFreq, ICell, + ConservTemp, AbsSalinity, + Pressure, SpecVol); + }); + } +} + /// Define IO fields and metadata for output void Eos::defineFields() { /// Set field names (append Name if not default) SpecVolFldName = "SpecVol"; SpecVolDisplacedFldName = "SpecVolDisplaced"; + BruntVaisalaFreqFldName = "BruntVaisalaFreq"; if (Name != "Default") { SpecVolFldName.append(Name); SpecVolDisplacedFldName.append(Name); + BruntVaisalaFreqFldName.append(Name); } /// Create fields for state variables - int NDims = 2; + const Real FillValue = -9.99e30; + int NDims = 2; std::vector DimNames(NDims); DimNames[0] = "NCells"; DimNames[1] = "NVertLayers"; @@ -200,8 +250,8 @@ void Eos::defineFields() { "m3 kg-1", // Units "sea_water_specific_volume", // CF-ish Name 0.0, // Min valid value - 9.99E+30, // Max valid value - -9.99E+30, // Scalar used for undefined entries + std::numeric_limits::max(), // Max valid value + FillValue, // Scalar used for undefined entries NDims, // Number of dimensions DimNames // Dimension names ); @@ -213,8 +263,20 @@ void Eos::defineFields() { "m3 kg-1", // Units "sea_water_specific_volume_displaced", // CF-ish Name 0.0, // Min valid value - 9.99E+30, // Max valid value - -9.99E+30, // Scalar used for undefined entried + std::numeric_limits::max(), // Max valid value + FillValue, // Scalar used for undefined entried + NDims, // Number of dimensions + DimNames // Dimension names + ); + /// Create and register the BruntVaisalaFreq field + auto BruntVaisalaFreqField = + Field::create(BruntVaisalaFreqFldName, // Field name + "Brunt-Vaisala frequency squared", // Long Name + "s-2", // Units + "sea_water_brunt_vaisala_frequency_squared", // 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 ); @@ -229,10 +291,12 @@ void Eos::defineFields() { // Add fields to the EOS group EosGroup->addField(SpecVolDisplacedFldName); EosGroup->addField(SpecVolFldName); + EosGroup->addField(BruntVaisalaFreqFldName); // Attach Kokkos views to the fields SpecVolDisplacedField->attachData(SpecVolDisplaced); SpecVolField->attachData(SpecVol); + BruntVaisalaFreqField->attachData(BruntVaisalaFreq); } // end defineIOFields diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index b9eb0f949c31..8c019c7b6ecd 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -29,10 +29,8 @@ enum class EosType { /// TEOS10 75-term Polynomial Equation of State class Teos10Eos { public: - Array2DReal SpecVolPCoeffs; - /// constructor declaration - Teos10Eos(int NVertLayers); + Teos10Eos(const VertCoord *VCoord); // The functor takes the full arrays of specific volume (inout), // the indices ICell and KChunk, and the ocean tracers (conservative) @@ -44,41 +42,43 @@ class Teos10Eos { const Array2DReal &Pressure, I4 KDisp) const { - OMEGA_SCOPE(LocSpecVolPCoeffs, SpecVolPCoeffs); + Real SpecVolPCoeffs[6 * VecLength]; const I4 KStart = KChunk * VecLength; for (int KVec = 0; KVec < VecLength; ++KVec) { const I4 K = KStart + KVec; - + if (K >= NVertLayers) + continue; /// Calculate the local specific volume polynomial pressure - /// coefficients - calcPCoeffs(LocSpecVolPCoeffs, KVec, ConservTemp(ICell, K), + /// coefficients with cell center values + calcPCoeffs(SpecVolPCoeffs, KVec, ConservTemp(ICell, K), AbsSalinity(ICell, K)); /// Calculate the specific volume at the given pressure - /// If KDisp is 0, we use the current pressure, otherwise we use the - /// displaced pressure (K + KDisp) + /// If KDisp is 0, we use the current pressure, otherwise + /// we use the displaced pressure (K + KDisp) /// Note: KDisp is only used for TEOS-10, for Linear EOS it /// is always 0. if (KDisp == 0) { // No displacement SpecVol(ICell, K) = calcRefProfile(Pressure(ICell, K)) + - calcDelta(LocSpecVolPCoeffs, KVec, Pressure(ICell, K)); + calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, K)); } else { // Displacement, use the displaced pressure I4 KTmp = Kokkos::min(K + KDisp, NVertLayers - 1); KTmp = Kokkos::max(0, KTmp); SpecVol(ICell, K) = calcRefProfile(Pressure(ICell, KTmp)) + - calcDelta(LocSpecVolPCoeffs, KVec, Pressure(ICell, KTmp)); + calcDelta(SpecVolPCoeffs, KVec, Pressure(ICell, KTmp)); } } } /// TEOS-10 helpers /// Calculate pressure polynomial coefficients for TEOS-10 - KOKKOS_FUNCTION void calcPCoeffs(Array2DReal SpecVolPCoeffs, const I4 K, - const Real Ct, const Real Sa) const { + KOKKOS_FUNCTION void calcPCoeffs(Real (&SpecVolPCoeffs)[6 * VecLength], + const I4 K, const Real Ct, + const Real Sa) const { constexpr Real SAu = 40.0 * 35.16504 / 35.0; constexpr Real CTu = 40.0; constexpr Real DeltaS = 24.0; @@ -161,18 +161,18 @@ class Teos10Eos { constexpr Real V014 = 3.1454099902e-07; constexpr Real V005 = 4.2369007180e-09; - SpecVolPCoeffs(5, K) = V005; - SpecVolPCoeffs(4, K) = V014 * Tt + V104 * Ss + V004; - SpecVolPCoeffs(3, K) = + SpecVolPCoeffs[5 + 6 * K] = V005; + SpecVolPCoeffs[4 + 6 * K] = V014 * Tt + V104 * Ss + V004; + SpecVolPCoeffs[3 + 6 * K] = (V023 * Tt + V113 * Ss + V013) * Tt + (V203 * Ss + V103) * Ss + V003; - SpecVolPCoeffs(2, K) = + SpecVolPCoeffs[2 + 6 * K] = (((V042 * Tt + V132 * Ss + V032) * Tt + (V222 * Ss + V122) * Ss + V022) * Tt + ((V312 * Ss + V212) * Ss + V112) * Ss + V012) * Tt + (((V402 * Ss + V302) * Ss + V202) * Ss + V102) * Ss + V002; - SpecVolPCoeffs(1, K) = + SpecVolPCoeffs[1 + 6 * K] = ((((V051 * Tt + V141 * Ss + V041) * Tt + (V231 * Ss + V131) * Ss + V031) * Tt + @@ -182,7 +182,7 @@ class Teos10Eos { Tt + ((((V501 * Ss + V401) * Ss + V301) * Ss + V201) * Ss + V101) * Ss + V001; - SpecVolPCoeffs(0, K) = + SpecVolPCoeffs[0 + 6 * K] = (((((V060 * Tt + V150 * Ss + V050) * Tt + (V240 * Ss + V140) * Ss + V040) * Tt + @@ -197,25 +197,25 @@ class Teos10Eos { V100) * Ss + V000; - - // could insert a check here (abs(value)> 0 value or = NVertLayers) + continue; SpecVol(ICell, K) = 1.0_Real / (RhoT0S0 + (DRhodT * ConservTemp(ICell, K) + DRhodS * AbsSalinity(ICell, K))); } } + + private: + I4 NVertLayers; +}; + +/// Functor for calculating the Brunt-Vaisala frequency using TEOS-10 +class Teos10BruntVaisalaFreq { + public: + /// Constructor for BruntVaisalaFreq + Teos10BruntVaisalaFreq(const VertCoord *VCoord); + + // The functor takes the full arrays of Brunt-Vaisala frequency (inout), + // the indix ICell, and the ocean tracers (conservative) temperature, + // (absolute) salinity, pressure, specific volume as inputs, and outputs + // the Brunt-Vaisala frequency. + KOKKOS_FUNCTION void operator()(Array2DReal BruntVaisalaFreq, I4 ICell, + const Array2DReal &ConservTemp, + const Array2DReal &AbsSalinity, + const Array2DReal &Pressure, + const Array2DReal &SpecVol) const { + + Real Gravity = 9.80616; + Real Db2Pa = 1.0e4; + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // No Brunt-Vaisala frequency at surface + BruntVaisalaFreq(ICell, K) = 0.0_Real; + } else { + // Calculate Brunt-Vaisala frequency + Real CtInt = + 0.5_Real * (ConservTemp(ICell, K) + ConservTemp(ICell, K - 1)); + Real SaInt = + 0.5_Real * (AbsSalinity(ICell, K) + AbsSalinity(ICell, K - 1)); + Real PInt = + 0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1)); + Real SpInt = 0.5_Real * (SpecVol(ICell, K) + SpecVol(ICell, K - 1)); + Real AlphaInt = calcAlpha(SaInt, CtInt, PInt, SpInt); + Real BetaInt = calcBeta(SaInt, CtInt, PInt, SpInt); + Real DSa = AbsSalinity(ICell, K) - AbsSalinity(ICell, K - 1); + Real DCt = ConservTemp(ICell, K) - ConservTemp(ICell, K - 1); + Real DP = Pressure(ICell, K) - Pressure(ICell, K - 1); + + BruntVaisalaFreq(ICell, K) = Gravity * Gravity * + ((1.0_Real / SpInt) / (Db2Pa * DP)) * + (BetaInt * DSa - AlphaInt * DCt); + } + } + } + + /// Calculate alpha values for the Brunt-Vaisala frequency + KOKKOS_FUNCTION Real calcAlpha(Real Sa, Real Ct, Real P, Real Sp) const { + constexpr Real Factor = 0.0248826675584615; + constexpr Real Offset = 5.971840214030754e-1; + constexpr Real Pa2Db = 1.0e-4; + Real Ss = Kokkos::sqrt(Factor * Sa + Offset); + Real Tt = 0.025_Real * Ct; + Real Pp = P * Pa2Db; + + constexpr Real A000 = -1.56497346750e-5; + constexpr Real A001 = 1.85057654290e-5; + constexpr Real A002 = -1.17363867310e-6; + constexpr Real A003 = -3.65270065530e-7; + constexpr Real A004 = 3.14540999020e-7; + constexpr Real A010 = 5.55242129680e-5; + constexpr Real A011 = -2.34332137060e-5; + constexpr Real A012 = 4.26100574800e-6; + constexpr Real A013 = 5.73918103180e-7; + constexpr Real A020 = -4.95634777770e-5; + constexpr Real A021 = 2.37838968519e-5; + constexpr Real A022 = -1.38397620111e-6; + constexpr Real A030 = 2.76445290808e-5; + constexpr Real A031 = -1.36408749928e-5; + constexpr Real A032 = -2.53411666056e-7; + constexpr Real A040 = -4.02698077700e-6; + constexpr Real A041 = 2.53683834070e-6; + constexpr Real A050 = 1.23258565608e-6; + constexpr Real A100 = 3.50095997640e-5; + constexpr Real A101 = -9.56770881560e-6; + constexpr Real A102 = -5.56991545570e-6; + constexpr Real A103 = -2.72956962370e-7; + constexpr Real A110 = -7.48716846880e-5; + constexpr Real A111 = -4.73566167220e-7; + constexpr Real A112 = 7.82747741600e-7; + constexpr Real A120 = 7.24244384490e-5; + constexpr Real A121 = -1.03676320965e-5; + constexpr Real A122 = 2.32856664276e-8; + constexpr Real A130 = -3.50383492616e-5; + constexpr Real A131 = 5.18268711320e-6; + constexpr Real A140 = -1.65263794500e-6; + constexpr Real A200 = -4.35926785610e-5; + constexpr Real A201 = 1.11008347650e-5; + constexpr Real A202 = 5.46207488340e-6; + constexpr Real A210 = 7.18156455200e-5; + constexpr Real A211 = 5.85666925900e-6; + constexpr Real A212 = -1.31462208134e-6; + constexpr Real A220 = -4.30608991440e-5; + constexpr Real A221 = 9.49659182340e-7; + constexpr Real A230 = 1.74814722392e-5; + constexpr Real A300 = 3.45324618280e-5; + constexpr Real A301 = -9.84471178440e-6; + constexpr Real A302 = -1.35441856270e-6; + constexpr Real A310 = -3.73971683740e-5; + constexpr Real A311 = -9.76522784000e-7; + constexpr Real A320 = 6.85899736680e-6; + constexpr Real A400 = -1.19594097880e-5; + constexpr Real A401 = 2.59092252600e-6; + constexpr Real A410 = 7.71906784880e-6; + constexpr Real A500 = 1.38645945810e-6; + + Real Rval = + A000 + + Ss * (A100 + Ss * (A200 + Ss * (A300 + Ss * (A400 + A500 * Ss)))) + + Tt * (A010 + Ss * (A110 + Ss * (A210 + Ss * (A310 + A410 * Ss))) + + Tt * (A020 + Ss * (A120 + Ss * (A220 + A320 * Ss)) + + Tt * (A030 + Ss * (A130 + A230 * Ss) + + Tt * (A040 + A140 * Ss + A050 * Tt)))) + + Pp * (A001 + Ss * (A101 + Ss * (A201 + Ss * (A301 + A401 * Ss))) + + Tt * (A011 + Ss * (A111 + Ss * (A211 + A311 * Ss)) + + Tt * (A021 + Ss * (A121 + A221 * Ss) + + Tt * (A031 + A131 * Ss + A041 * Tt))) + + Pp * (A002 + Ss * (A102 + Ss * (A202 + A302 * Ss)) + + Tt * (A012 + Ss * (A112 + A212 * Ss) + + Tt * (A022 + A122 * Ss + A032 * Tt)) + + Pp * (A003 + A103 * Ss + A013 * Tt + A004 * Pp))); + + return 0.025_Real * Rval / Sp; + } + + /// Calculate alpha values for the Brunt-Vaisala frequency + KOKKOS_FUNCTION Real calcBeta(Real Sa, Real Ct, Real P, Real Sp) const { + constexpr Real Factor = 0.0248826675584615; + constexpr Real Offset = 5.971840214030754e-1; + constexpr Real Pa2Db = 1.0e-4; + Real Ss = Kokkos::sqrt(Factor * Sa + Offset); + Real Tt = 0.025_Real * Ct; + Real Pp = P * Pa2Db; + + constexpr Real B000 = -3.10389819760e-4; + constexpr Real B003 = 3.63101885150e-7; + constexpr Real B004 = -1.11471254230e-7; + constexpr Real B010 = 3.50095997640e-5; + constexpr Real B013 = -2.72956962370e-7; + constexpr Real B020 = -3.74358423440e-5; + constexpr Real B030 = 2.41414794830e-5; + constexpr Real B040 = -8.75958731540e-6; + constexpr Real B050 = -3.30527589000e-7; + constexpr Real B100 = 1.33856134076e-3; + constexpr Real B103 = 3.34926075600e-8; + constexpr Real B110 = -8.71853571220e-5; + constexpr Real B120 = 7.18156455200e-5; + constexpr Real B130 = -2.87072660960e-5; + constexpr Real B140 = 8.74073611960e-6; + constexpr Real B200 = -2.55143801811e-3; + constexpr Real B210 = 1.03597385484e-4; + constexpr Real B220 = -5.60957525610e-5; + constexpr Real B230 = 6.85899736680e-6; + constexpr Real B300 = 2.32344279772e-3; + constexpr Real B310 = -4.78376391520e-5; + constexpr Real B320 = 1.54381356976e-5; + constexpr Real B400 = -1.05461852535e-3; + constexpr Real B410 = 6.93229729050e-6; + constexpr Real B500 = 1.91594743830e-4; + constexpr Real B001 = 2.42624687470e-5; + constexpr Real B011 = -9.56770881560e-6; + constexpr Real B021 = -2.36783083610e-7; + constexpr Real B031 = -3.45587736550e-6; + constexpr Real B041 = 1.29567177830e-6; + constexpr Real B101 = -6.95849219480e-5; + constexpr Real B111 = 2.22016695300e-5; + constexpr Real B121 = 5.85666925900e-6; + constexpr Real B131 = 6.33106121560e-7; + constexpr Real B201 = 1.12412331915e-4; + constexpr Real B211 = -2.95341353532e-5; + constexpr Real B221 = -1.46478417600e-6; + constexpr Real B301 = -6.92888744480e-5; + constexpr Real B311 = 1.03636901040e-5; + constexpr Real B401 = 1.54637136265e-5; + constexpr Real B002 = -5.84844329840e-7; + constexpr Real B012 = -5.56991545570e-6; + constexpr Real B022 = 3.91373870800e-7; + constexpr Real B032 = 7.76188880920e-9; + constexpr Real B102 = -9.62445031940e-6; + constexpr Real B112 = 1.09241497668e-5; + constexpr Real B122 = -1.31462208134e-6; + constexpr Real B202 = 1.47789320994e-5; + constexpr Real B212 = -4.06325568810e-6; + constexpr Real B302 = -7.12478989080e-6; + + Real Rval = + B000 + + Ss * (B100 + Ss * (B200 + Ss * (B300 + Ss * (B400 + B500 * Ss)))) + + Tt * (B010 + Ss * (B110 + Ss * (B210 + Ss * (B310 + B410 * Ss))) + + Tt * (B020 + Ss * (B120 + Ss * (B220 + B320 * Ss)) + + Tt * (B030 + Ss * (B130 + B230 * Ss) + + Tt * (B040 + B140 * Ss + B050 * Tt)))) + + Pp * (B001 + Ss * (B101 + Ss * (B201 + Ss * (B301 + B401 * Ss))) + + Tt * (B011 + Ss * (B111 + Ss * (B211 + B311 * Ss)) + + Tt * (B021 + Ss * (B121 + B221 * Ss) + + Tt * (B031 + B131 * Ss + B041 * Tt))) + + Pp * (B002 + Ss * (B102 + Ss * (B202 + B302 * Ss)) + + Tt * (B012 + Ss * (B112 + B212 * Ss) + + Tt * (B022 + B122 * Ss + B032 * Tt)) + + Pp * (B003 + B103 * Ss + B013 * Tt + B004 * Pp))); + + return -0.5_Real * Rval * Factor / (Sp * Ss); + } + + private: + I4 NVertLayers; +}; + +/// Linear Brunt-Vaisala frequency calculator +class LinearBruntVaisalaFreq { + 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 + LinearBruntVaisalaFreq(const VertCoord *VCoord); + + // The functor takes the full arrays of Brunt-Vaisala frequency (inout), + // the index ICell, and the specific volume and layer thickness as inputs, + // and outputs the Brunt-Vaisala frequency. + KOKKOS_FUNCTION void operator()(Array2DReal BruntVaisalaFreq, I4 ICell, + const Array2DReal &SpecVol) const { + Real Gravity = 9.80616; // gravitational acceleration + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + /// No Brunt-Vaisala frequency at the top level + BruntVaisalaFreq(ICell, K) = 0.0_Real; + } else { + /// Calculate Brunt-Vaisala frequency at mid-point between K-1 and K + /// Do not need to use displaced specific volume here since only the + /// linear EOS is used with this BVF formulation. + BruntVaisalaFreq(ICell, K) = -(Gravity / RhoT0S0) * + ((1.0_Real / SpecVol(ICell, K - 1)) - + (1.0_Real / SpecVol(ICell, K))) / + (ZMid(ICell, K - 1) - ZMid(ICell, K)); + } + } + } + + private: + Array2DReal ZMid; + I4 NVertLayers; }; /// Class for Equation of State (EOS) calculations @@ -278,12 +525,15 @@ class Eos { static void destroyInstance(); EosType EosChoice; ///< Current EOS type in use - Array2DReal SpecVol; ///< Specific volume field + Array2DReal SpecVol; ///< Specific volume field at level centers Array2DReal SpecVolDisplaced; ///< Displaced specific volume field + Array2DReal BruntVaisalaFreq; ///< Brunt-Vaisala frequency field std::string SpecVolFldName; ///< Field name for specific volume std::string SpecVolDisplacedFldName; ///< Field name for displaced specific volume + std::string + BruntVaisalaFreqFldName; ///< Field name for Brunt-Vaisala frequency std::string EosGroupName; ///< EOS group name (for config) std::string Name; ///< Name of this EOS instance @@ -297,12 +547,18 @@ class Eos { const Array2DReal &AbsSalinity, const Array2DReal &Pressure, I4 KDisp); + /// Compute Brunt-Vaisala frequency for all cells/layers + void computeBruntVaisalaFreq(const Array2DReal &ConservTemp, + const Array2DReal &AbsSalinity, + const Array2DReal &Pressure, + const Array2DReal &SpecVol); + /// Initialize EOS from config and mesh static void init(); private: /// Private constructor - Eos(const std::string &Name, const HorzMesh *Mesh, int NVertLayers); + Eos(const std::string &Name, const HorzMesh *Mesh, const VertCoord *VCoord); /// Private destructor ~Eos(); @@ -316,11 +572,17 @@ class Eos { Eos(Eos &&) = delete; Eos &operator=(Eos &&) = delete; - I4 NCellsAll; ///< Number of horizontal cells - I4 NChunks; ///< Number of vertical chunks (for vectorization) + I4 NCellsAll; ///< Number of horizontal cells + I4 NChunks; ///< Number of vertical chunks (for vectorization) + I4 NVertLayers; ///< Number of vertical layers + const HorzMesh *Mesh; Teos10Eos ComputeSpecVolTeos10; ///< TEOS-10 specific volume calculator LinearEos ComputeSpecVolLinear; ///< Linear specific volume calculator + Teos10BruntVaisalaFreq + ComputeBruntVaisalaFreqTeos10; ///< TEOS-10 Brunt-Vaisala calculator + LinearBruntVaisalaFreq + ComputeBruntVaisalaFreqLinear; ///< Linear Brunt-Vaisala calculator // Define fields and metadata void defineFields(); diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index f08b18458e6d..a49c252c9624 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -643,6 +643,7 @@ void HorzMesh::copyToDevice() { YCell = createDeviceMirrorCopy(YCellH); XEdge = createDeviceMirrorCopy(XEdgeH); YEdge = createDeviceMirrorCopy(YEdgeH); + LatCell = createDeviceMirrorCopy(LatCellH); } // end copyToDevice diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index a1a73a0db83a..b0664647b295 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -163,6 +163,8 @@ class HorzMesh { HostArray1DReal ZCellH; ///< Z Coordinates of cell centers (m) HostArray1DReal LonCellH; ///< Longitude location of cell centers (radians) + + Array1DReal LatCell; ///< Latitude location of cell centers (radians) HostArray1DReal LatCellH; ///< Latitude location of cell centers (radians) Array1DReal XEdge; ///< X Coordinate of edge midpoints (m) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp new file mode 100644 index 000000000000..ebeddef19dbf --- /dev/null +++ b/components/omega/src/ocn/VertMix.cpp @@ -0,0 +1,270 @@ +//===-- ocn/VertMix.cpp - Vertical Mixing Coefficients -----------*- C++ +//-*-===// +// +// The VertMix class is responsible for managing the calculation of the +// vertical diffusivity and viscosity needed for the vertical mixing. +// It currently has background, convective, and shear mixing options, and +// they can be additively combined, depending on configuration options. It +// contains arrays that store the vertical top-of-cell diffusivity and +// viscosity values for each cell and vertical level. +// +//===----------------------------------------------------------------------===// + +#include "VertMix.h" +#include "DataTypes.h" +#include "HorzMesh.h" + +namespace OMEGA { + +ShearMix::ShearMix(const HorzMesh *Mesh, const VertCoord *VCoord) + : NVertLayers(VCoord->NVertLayers), ZMid(VCoord->ZMid), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), + DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AreaCell(Mesh->AreaCell) {} + +ConvectiveMix::ConvectiveMix(const VertCoord *VCoord) + : NVertLayers(VCoord->NVertLayers) {} + +/// 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) { + VertDiff = Array2DReal("VertDiff", Mesh->NCellsAll, VCoord->NVertLayers + 1); + VertVisc = Array2DReal("VertVisc", Mesh->NCellsAll, VCoord->NVertLayers + 1); + NCellsAll = Mesh->NCellsAll; + NChunks = VCoord->NVertLayers / VecLength; + Name = Name_; + + defineFields(); +} + +/// Destructor for VertMix +VertMix::~VertMix() {} + +/// Instance management +VertMix *VertMix::Instance = nullptr; + +/// Get instance of VertMix +VertMix *VertMix::getInstance() { return Instance; } + +/// Destroy instance of VertMix +void VertMix::destroyInstance() { + delete Instance; + Instance = nullptr; +} + +/// Initializes the VertMix (Vertical Mixing Coefficients) class and its +/// 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 +/// missing. +void VertMix::init() { + + if (!Instance) { + Instance = new VertMix("Default", HorzMesh::getDefault(), + VertCoord::getDefault()); + } + + Error Err; // error code + + /// Retrieve default VertMix + VertMix *vertmix = VertMix::getInstance(); + + /// Get VertMixConfig group from Omega config + Config *OmegaConfig = Config::getOmegaConfig(); + Config VertMixConfig("VertMix"); + Err += OmegaConfig->get(VertMixConfig); + CHECK_ERROR_ABORT(Err, "VertMix::init: VertMix group not found in Config"); + + /// Get Background from VertMixConfig + /// and set associated parameters + Config BackConfig("Background"); + Err += VertMixConfig.get(BackConfig); + CHECK_ERROR_ABORT( + Err, "VertMix::init: Background subgroup not found in VertMixConfig"); + + /// Get diffusivity and viscosity parameters + Err += BackConfig.get("Viscosity", vertmix->BackVisc); + CHECK_ERROR_ABORT( + Err, + "VertMix::init: Parameter Background:Viscosity not found in BackConfig"); + + Err += BackConfig.get("Diffusivity", vertmix->BackDiff); + CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter Background:Diffusivity not " + "found in BackConfig"); + + /// Get Convective from VertMixConfig + Config ConvConfig("Convective"); + Err += VertMixConfig.get(ConvConfig); + CHECK_ERROR_ABORT( + Err, "VertMix::init: Convective subgroup not found in VertMixConfig"); + + /// Get convective diffusivity and viscosity parameters + Err += ConvConfig.get("Enable", vertmix->ComputeVertMixConv.Enabled); + CHECK_ERROR_ABORT( + Err, + "VertMix::init: Parameter Convective:Enable not found in ConvConfig"); + + if (!vertmix->ComputeVertMixConv.Enabled) { + LOG_INFO("VertMix::init: Convective mixing is disabled."); + } else { + LOG_INFO("VertMix::init: Convective mixing is enabled."); + Err += + ConvConfig.get("Diffusivity", vertmix->ComputeVertMixConv.ConvDiff); + CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter Convective:Diffusivity " + "not found in ConvConfig"); + + Err += ConvConfig.get("TriggerBVF", + vertmix->ComputeVertMixConv.ConvTriggerBVF); + CHECK_ERROR_ABORT(Err, "VertMix::init: Parameter Convective:TriggerBVF " + "not found in ConvConfig"); + } + + /// Get Shear from VertMixConfig + Config ShearConfig("Shear"); + Err += VertMixConfig.get(ShearConfig); + CHECK_ERROR_ABORT( + Err, "VertMix::init: Shear subgroup not found in VertMixConfig"); + + /// Get shear diffusivity and viscosity parameters + Err += ShearConfig.get("Enable", vertmix->ComputeVertMixShear.Enabled); + CHECK_ERROR_ABORT( + Err, "VertMix::init: Parameter Shear:Enable not found in ShearConfig"); + + if (!vertmix->ComputeVertMixShear.Enabled) { + LOG_INFO("VertMix::init: Shear mixing is disabled."); + } else { + LOG_INFO("VertMix::init: Shear mixing is enabled."); + Err += + ShearConfig.get("NuZero", vertmix->ComputeVertMixShear.ShearNuZero); + CHECK_ERROR_ABORT( + Err, + "VertMix::init: Parameter Shear:NuZero not found in ShearConfig"); + + Err += ShearConfig.get("Alpha", vertmix->ComputeVertMixShear.ShearAlpha); + CHECK_ERROR_ABORT( + Err, "VertMix::init: Parameter Shear:Alpha not found in ShearConfig"); + + Err += ShearConfig.get("Exponent", + vertmix->ComputeVertMixShear.ShearExponent); + CHECK_ERROR_ABORT( + Err, + "VertMix::init: Parameter Shear:Exponent not found in ShearConfig"); + } +} // end init + +/// Compute diffusivity and viscosity for all cells/layers (no displacement) +void VertMix::computeVertMix(const Array2DReal &NormalVelocity, + const Array2DReal &TangentialVelocity, + const Array2DReal &BruntVaisalaFreq) { + OMEGA_SCOPE(LocVertDiff, VertDiff); /// Create a local view for computation + OMEGA_SCOPE(LocVertVisc, VertVisc); /// Create a local view for computation + OMEGA_SCOPE( + LocComputeVertMixConv, + ComputeVertMixConv); /// Local view for Convective VertMix computation + OMEGA_SCOPE( + LocComputeVertMixShear, + ComputeVertMixShear); /// Local view for Shear VertMix computation + + /// Initialize VertDiff and VertVisc to background values + deepCopy(LocVertDiff, BackDiff); + deepCopy(LocVertVisc, BackVisc); + + /// Dispatch to the correct VertMix calculation + if (LocComputeVertMixShear.Enabled && LocComputeVertMixConv.Enabled) { + parallelFor( + "VertMix-ConvPlusShear", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, + BruntVaisalaFreq); + LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, + NormalVelocity, TangentialVelocity, + BruntVaisalaFreq); + }); + } else if (LocComputeVertMixShear.Enabled) { + parallelFor( + "VertMix-ShearOnly", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, + NormalVelocity, TangentialVelocity, + BruntVaisalaFreq); + }); + } else if (LocComputeVertMixConv.Enabled) { + parallelFor( + "VertMix-ConvOnly", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, + BruntVaisalaFreq); + }); + } else { + parallelFor( + "VertMix-Background", {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) { + LocVertDiff(ICell, 0) = 0.0_Real; + LocVertVisc(ICell, 0) = 0.0_Real; + }); + } + Kokkos::fence(); +} + +/// Define IO fields and metadata for output +void VertMix::defineFields() { + + /// Set field names (append Name if not default) + VertDiffFldName = "VertDiff"; + VertViscFldName = "VertVisc"; + if (Name != "Default") { + VertDiffFldName.append(Name); + VertViscFldName.append(Name); + } + + /// Create fields for state variables + const Real FillValue = -9.99e30; + int NDims = 2; + std::vector DimNames(NDims); + DimNames[0] = "NCells"; + DimNames[1] = "NVertLayers"; + + /// Create and register the Diffusivity field + auto VertDiffField = + Field::create(VertDiffFldName, // Field name + "Vertical diffusivity at center and" + " top of cell", // Long Name + "m2 s-1", // Units + "vertical_diffusivity", // CF-ish Name + 0.0, // 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 VertVisc field + auto VertViscField = + Field::create(VertViscFldName, // Field name + "Vertical viscosity at center and" + " top of cell", // Long Name + "m2 s-1", // Units + "vertical_viscosity", // CF-ish Name + 0.0, // Min valid value + std::numeric_limits::max(), // Max valid value + FillValue, // Scalar used for undefined entried + NDims, // Number of dimensions + DimNames // Dimension names + ); + + // Create a field group for the vertmix-specific state fields + VertMixGroupName = "VertMix"; + if (Name != "Default") { + VertMixGroupName.append(Name); + } + auto VertMixGroup = FieldGroup::create(VertMixGroupName); + + // Add fields to the VertMix group + VertMixGroup->addField(VertDiffFldName); + VertMixGroup->addField(VertViscFldName); + + // Attach Kokkos views to the fields + VertDiffField->attachData(VertDiff); + VertViscField->attachData(VertVisc); + +} // end defineIOFields + +} // namespace OMEGA diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h new file mode 100644 index 000000000000..9b50e08461b1 --- /dev/null +++ b/components/omega/src/ocn/VertMix.h @@ -0,0 +1,183 @@ +#ifndef OMEGA_VERTMIX_H +#define OMEGA_VERTMIX_H +//===-- ocn/VertMix.cpp - Vertical Mixing Coefficients -----------*- C++ +//-*-===// +// +/// \file +/// \brief Contains functors for calculating vertical diffusivity and viscosity +/// +/// This header defines functors to be called by the time-stepping scheme +/// to calculate the vertical diffusivity and viscosity +// +//===----------------------------------------------------------------------===// + +#include "AuxiliaryState.h" +#include "Config.h" +#include "HorzMesh.h" +#include "MachEnv.h" +#include "OmegaKokkos.h" +#include "TimeMgr.h" +#include "VertCoord.h" +#include + +namespace OMEGA { + +class ConvectiveMix { + public: + bool Enabled = true; ///< Enable convective mixing flag + + // Convective mixing parameters + Real ConvDiff = + 1.0; ///< Convective vertical viscosity and diffusivity (m^2 s^-1) + Real ConvTriggerBVF = 0.0; ///< Reference density (kg m^-3) at (T,S)=(0,0) + + /// Constructor for ConvectiveMix + ConvectiveMix(const VertCoord *VCoord); + + KOKKOS_FUNCTION void operator()(Array2DReal VertDiff, Array2DReal VertVisc, + I4 ICell, + const Array2DReal &BruntVaisalaFreq) const { + for (int K = 1; K < NVertLayers; ++K) { + if (K == 0) { + VertVisc(ICell, K) = 0.0_Real; + VertDiff(ICell, K) = 0.0_Real; + } else { + if (BruntVaisalaFreq(ICell, K) < ConvTriggerBVF) { + VertDiff(ICell, K) = VertDiff(ICell, K) + ConvDiff; + VertVisc(ICell, K) = VertVisc(ICell, K) + ConvDiff; + } + } + } + } + + private: + I4 NVertLayers; +}; + +class ShearMix { + public: + 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 ShearExponent = + 2.0; /// Exponent value used in Pacanowski and Philander (1981) Eqs (1). + + /// Constructor for ShearMix + ShearMix(const HorzMesh *Mesh, const VertCoord *VCoord); + + KOKKOS_FUNCTION void operator()(Array2DReal VertDiff, Array2DReal VertVisc, + I4 ICell, const Array2DReal &NormalVelocity, + const Array2DReal &TangentialVelocity, + const Array2DReal &BruntVaisalaFreq) const { + + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + VertVisc(ICell, K) = 0.0_Real; + VertDiff(ICell, K) = 0.0_Real; + } else { + Real ShearViscVal = 0.0; + Real InvAreaCell = 1.0 / 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 = ZMid(ICell, K - 1) - ZMid(ICell, K); + ShearSquared = ShearSquared / (DelZMid * DelZMid); + + Real RichardsonNum = BruntVaisalaFreq(ICell, K) / + Kokkos::max(1.0e-12_Real, ShearSquared); + + ShearViscVal = + ShearNuZero / Kokkos::pow(1.0_Real + ShearAlpha * RichardsonNum, + ShearExponent); + VertVisc(ICell, K) = VertVisc(ICell, K) + ShearViscVal; + VertDiff(ICell, K) = + VertDiff(ICell, K) + + ShearViscVal / (1.0_Real + ShearAlpha * RichardsonNum); + } + } + } + + private: + Array1DReal DcEdge; + Array1DReal DvEdge; + Array1DReal AreaCell; + Array2DReal ZMid; + Array1DI4 NEdgesOnCell; + Array2DI4 EdgesOnCell; + I4 NVertLayers; +}; + +/// Class for Vertical Mixing Coefficient (VertMix) calculations +class VertMix { + public: + /// Get instance of VertMix + static VertMix *getInstance(); + + /// 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 + + // Background mixing parameters + Real BackDiff = 1.0e-5; ///< Background vertical diffusivity (m^2 s^-1) + Real BackVisc = 1.0e-4; ///< Background vertical viscosity (m^2 s^-1) + + ConvectiveMix + ComputeVertMixConv; ///< Functor for Convective VertMix calculation + ShearMix ComputeVertMixShear; ///< Functor for Shear VertMix calculation + + /// Compute vertical diffusivity and viscosity for all cells/layers + void computeVertMix(const Array2DReal &NormalVelocity, + const Array2DReal &TangentialVelocity, + const Array2DReal &BruntVaisalaFreq); + + /// Initialize VertMix from config and mesh + static void init(); + + private: + /// Private constructor + VertMix(const std::string &Name, const HorzMesh *Mesh, + const VertCoord *VCoord); + + /// Private destructor + ~VertMix(); + + static VertMix *Instance; ///< Instance pointer + + // Delete copy and move constructors and assignment operators + VertMix(const VertMix &) = delete; + VertMix &operator=(const VertMix &) = delete; + VertMix(VertMix &&) = delete; + VertMix &operator=(VertMix &&) = delete; + + I4 NCellsAll; ///< Number of horizontal cells + I4 NEdgesAll; ///< Number of horizontal edges + I4 NChunks; ///< Number of vertical chunks (for vectorization) + I4 NVertLayers; ///< Number of vertical layers + + // Define fields and metadata + void defineFields(); + +}; // End class VertMix + +} // namespace OMEGA +#endif diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 4a1ec2bf1a9a..12906e8dbd05 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -458,3 +458,14 @@ add_omega_test( ocn/VertCoordTest.cpp "-n;8" ) + +########################## +# Vert Mix test +########################## + +add_omega_test( + VERTMIX_TEST + testVertMix.exe + ocn/VertMixTest.cpp + "-n;4" +) diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 028d136e9c16..91feb0b96f5c 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -37,16 +37,24 @@ using namespace OMEGA; constexpr int NVertLayers = 60; /// Published values (TEOS-10 and linear) to test against -const Real TeosExpValue = 0.0009732819628; // Expected value for TEOS-10 eos +const Real TeosSVExpValue = + 0.0009732819628; // Expected value for TEOS-10 specific volume const Real LinearExpValue = - 0.0009784735812133072; // Expected value for Linear eos + 0.0009784735812133072; // Expected value for Linear specific volume +const Real TeosBVFExpValue = + 0.020911744281268286; // Expected value for TEOS-10 Brunt-Vaisala frequency +const Real LinearBVFExpValue = + 0.017833905406889013; // Expected value for Linear Brunt-Vaisala frequency +const Real GswBVFExpValue = + 0.02081197958166906; // Expected value from GSW-C library /// Test input values -double Sa = 30.0; // Absolute Salinity in g/kg -double Ct = 10.0; // Conservative Temperature in degC -double P = 1000.0; // Pressure in dbar -const I4 KDisp = 1; // Displace parcel to K=1 for TEOS-10 eos -const Real RTol = 1e-10; // Relative tolerance for isApprox checks +const Real Sa = 30.0; // Absolute Salinity in g/kg +const Real Ct = 10.0; // Conservative Temperature in degC +const Real P = 1000.0; // Pressure in dbar + +const I4 KDisp = 1; // Displace parcel to K=1 for TEOS-10 eos +const Real RTol = 1e-10; // Relative tolerance for isApprox checks /// The initialization routine for Eos testing. It calls various /// init routines, including the creation of the default decomposition. @@ -76,6 +84,8 @@ I4 initEosTest(const std::string &mesh) { /// Initialize vertical coordinate (phase 1) VertCoord::init1(); + auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; /// Initialize mesh HorzMesh::init(); @@ -98,16 +108,21 @@ I4 initEosTest(const std::string &mesh) { /// Test Linear EOS calculation for all cells/layers int testEosLinear() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::LinearEos; /// Create and fill ocean state arrays - Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLayers); - Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLayers); - Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLayers); + Array2DReal SArray = + Array2DReal("SArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal TArray = + Array2DReal("TArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal PArray = + Array2DReal("PArray", Mesh->NCellsAll, VCoord->NVertLayers); /// Use Kokkos::deep_copy to fill the entire view with the ref value deepCopy(SArray, Sa); deepCopy(TArray, Ct); @@ -121,7 +136,7 @@ int testEosLinear() { int numMismatches = 0; Array2DReal SpecVol = TestEos->SpecVol; parallelReduce( - "CheckSpecVolMatrix-linear", {Mesh->NCellsAll, NVertLayers}, + "CheckSpecVolMatrix-linear", {Mesh->NCellsAll, VCoord->NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { if (!isApprox(SpecVol(i, j), LinearExpValue, RTol)) { localCount++; @@ -145,16 +160,21 @@ int testEosLinear() { /// Test Linear EOS calculation with vertical displacement int testEosLinearDisplaced() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::LinearEos; /// Create and fill ocean state arrays - Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLayers); - Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLayers); - Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLayers); + Array2DReal SArray = + Array2DReal("SArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal TArray = + Array2DReal("TArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal PArray = + Array2DReal("PArray", Mesh->NCellsAll, VCoord->NVertLayers); /// Use Kokkos::deep_copy to fill the entire view with the ref value deepCopy(SArray, Sa); deepCopy(TArray, Ct); @@ -168,7 +188,7 @@ int testEosLinearDisplaced() { int numMismatches = 0; Array2DReal SpecVolDisplaced = TestEos->SpecVolDisplaced; parallelReduce( - "CheckSpecVolDispMatrix-linear", {Mesh->NCellsAll, NVertLayers}, + "CheckSpecVolDispMatrix-linear", {Mesh->NCellsAll, VCoord->NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { if (!isApprox(SpecVolDisplaced(i, j), LinearExpValue, RTol)) { localCount++; @@ -190,18 +210,112 @@ int testEosLinearDisplaced() { return Err; } +/// Test linear Brunt-Vaisala frequency calculation for all cells/layers +int testBruntVaisalaFreqLinear() { + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::LinearEos; + + /// Create and fill ocean state arrays + Array2DReal SArray = + Array2DReal("SArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal TArray = + Array2DReal("TArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal PArray = + Array2DReal("PArray", Mesh->NCellsAll, VCoord->NVertLayers); + /// Use Kokkos::deep_copy to fill the entire view with the ref value + // deepCopy(SArray, Sa); + // deepCopy(TArray, Ct); + // deepCopy(PArray, P); + deepCopy(TestEos->SpecVol, 0.0); + deepCopy(TestEos->BruntVaisalaFreq, 0.0); + + // parallelFor( + // "populateArrays", {Mesh->NCellsAll, NVertLayers}, + // KOKKOS_LAMBDA(I4 ICell, I4 K) { + // VCoord->ZMid(ICell,K) = -K; + // SArray(ICell,K) = Sa + 0.1*K; + // TArray(ICell,K) = Ct + 0.25*K; + // PArray(ICell,K) = P + 0.25*K; + // }); + + parallelFor( + "populateArrays", {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + VCoord->ZMid(ICell, 0) = -992.1173890198451; + VCoord->ZMid(ICell, 1) = -993.1071379053125; + VCoord->ZMid(ICell, 2) = -994.0968821072275; + SArray(ICell, 0) = Sa - 1.0; + SArray(ICell, 1) = Sa; + SArray(ICell, 2) = Sa + 1.0; + TArray(ICell, 0) = Ct + 15.0; + TArray(ICell, 1) = Ct + 10.0; + TArray(ICell, 2) = Ct + 5.0; + PArray(ICell, 0) = P; + PArray(ICell, 1) = P + 1.0; + PArray(ICell, 2) = P + 2.0; + }); + + /// Compute specific volume first + TestEos->computeSpecVol(TArray, SArray, PArray); + Array2DReal SpecVol = TestEos->SpecVol; + + /// Compute Brunt-Vaisala frequency + TestEos->computeBruntVaisalaFreq(TArray, SArray, PArray, SpecVol); + Array2DReal BruntVaisalaFreq = TestEos->BruntVaisalaFreq; + auto BruntVaisalaFreqH = createHostMirrorCopy(BruntVaisalaFreq); + /// Check all array values against expected value + int numMismatches = 0; + if (!isApprox(BruntVaisalaFreqH(1, 1), LinearBVFExpValue, RTol)) { + numMismatches = 1; + } else { + numMismatches = 0; + } + + // int numMismatches = 0; + // parallelReduce( + // "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, + // KOKKOS_LAMBDA(int i, int j, int &localCount) { + // if (!isApprox(BruntVaisalaFreq(i, j), LinearBVFExpValue, RTol)) { + // localCount++; + // } + // }, + // numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("EosTest: Linear BruntVaisalaFreq isApprox FAIL, " + "expected {}, got {} with {} mismatches", + LinearBVFExpValue, BruntVaisalaFreqH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("EosTest BruntVaisalaFreqCalc Linear: PASS"); + } + + return Err; +} + /// Test TEOS-10 EOS calculation for all cells/layers int testEosTeos10() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::Teos10Eos; /// Create and fill ocean state arrays - Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLayers); - Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLayers); - Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLayers); + Array2DReal SArray = + Array2DReal("SArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal TArray = + Array2DReal("TArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal PArray = + Array2DReal("PArray", Mesh->NCellsAll, VCoord->NVertLayers); /// Use Kokkos::deep_copy to fill the entire view with the ref value deepCopy(SArray, Sa); deepCopy(TArray, Ct); @@ -215,9 +329,9 @@ int testEosTeos10() { int numMismatches = 0; Array2DReal SpecVol = TestEos->SpecVol; parallelReduce( - "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, + "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, VCoord->NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { - if (!isApprox(SpecVol(i, j), TeosExpValue, RTol)) { + if (!isApprox(SpecVol(i, j), TeosSVExpValue, RTol)) { localCount++; } }, @@ -228,7 +342,7 @@ int testEosTeos10() { Err++; LOG_ERROR("EosTest: TEOS SpecVol isApprox FAIL, " "expected {}, got {} with {} mismatches", - TeosExpValue, SpecVolH(1, 1), numMismatches); + TeosSVExpValue, SpecVolH(1, 1), numMismatches); } if (Err == 0) { LOG_INFO("EosTest SpecVolCalc TEOS-10: PASS"); @@ -239,16 +353,21 @@ int testEosTeos10() { /// Test TEOS-10 EOS calculation with vertical displacement int testEosTeos10Displaced() { - int Err = 0; - const auto *Mesh = HorzMesh::getDefault(); + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::Teos10Eos; /// Create and fill ocean state arrays - Array2DReal SArray = Array2DReal("SArray", Mesh->NCellsAll, NVertLayers); - Array2DReal TArray = Array2DReal("TArray", Mesh->NCellsAll, NVertLayers); - Array2DReal PArray = Array2DReal("PArray", Mesh->NCellsAll, NVertLayers); + Array2DReal SArray = + Array2DReal("SArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal TArray = + Array2DReal("TArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal PArray = + Array2DReal("PArray", Mesh->NCellsAll, VCoord->NVertLayers); /// Use Kokkos::deep_copy to fill the entire view with the ref value deepCopy(SArray, Sa); deepCopy(TArray, Ct); @@ -262,9 +381,9 @@ int testEosTeos10Displaced() { int numMismatches = 0; Array2DReal SpecVolDisplaced = TestEos->SpecVolDisplaced; parallelReduce( - "CheckSpecVolDispMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, + "CheckSpecVolDispMatrix-Teos", {Mesh->NCellsAll, VCoord->NVertLayers}, KOKKOS_LAMBDA(int i, int j, int &localCount) { - if (!isApprox(SpecVolDisplaced(i, j), TeosExpValue, RTol)) { + if (!isApprox(SpecVolDisplaced(i, j), TeosSVExpValue, RTol)) { localCount++; } }, @@ -275,7 +394,7 @@ int testEosTeos10Displaced() { Err++; LOG_ERROR("EosTest: TEOS SpecVolDisp isApprox FAIL, " "expected {}, got {} with {} mismatches", - TeosExpValue, SpecVolDisplacedH(1, 1), numMismatches); + TeosSVExpValue, SpecVolDisplacedH(1, 1), numMismatches); } if (Err == 0) { LOG_INFO("EosTest SpecVolCalcDisp TEOS-10: PASS"); @@ -284,6 +403,93 @@ int testEosTeos10Displaced() { return Err; } +/// Test TEOS-10 Brunt-Vaisala frequency calculation for all cells/layer +int testBruntVaisalaFreqTeos10() { + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::Teos10Eos; + + /// Create and fill ocean state arrays + Array2DReal SArray = + Array2DReal("SArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal TArray = + Array2DReal("TArray", Mesh->NCellsAll, VCoord->NVertLayers); + Array2DReal PArray = + Array2DReal("PArray", Mesh->NCellsAll, VCoord->NVertLayers); + /// Use Kokkos::deep_copy to fill the entire view with the ref value + // deepCopy(SArray, Sa); + // deepCopy(TArray, Ct); + // deepCopy(PArray, P); + deepCopy(TestEos->BruntVaisalaFreq, 0.0); + deepCopy(TestEos->SpecVol, 0.0); + + // parallelFor( + // "populateArrays", {Mesh->NCellsAll, NVertLayers}, + // KOKKOS_LAMBDA(I4 ICell, I4 K) { + // VCoord->ZMid(ICell,K) = -K; + // SArray(ICell,K) = Sa + 0.1*K; + // TArray(ICell,K) = Ct + 0.25*K; + // PArray(ICell,K) = P + 0.25*K; + // }); + + parallelFor( + "populateArrays", {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + VCoord->ZMid(ICell, 0) = -992.1173890198451; + VCoord->ZMid(ICell, 1) = -993.1071379053125; + VCoord->ZMid(ICell, 2) = -994.0968821072275; + SArray(ICell, 0) = Sa - 1.0; + SArray(ICell, 1) = Sa; + SArray(ICell, 2) = Sa + 1.0; + TArray(ICell, 0) = Ct + 15.0; + TArray(ICell, 1) = Ct + 10.0; + TArray(ICell, 2) = Ct + 5.0; + PArray(ICell, 0) = P; + PArray(ICell, 1) = P + 1.0; + PArray(ICell, 2) = P + 2.0; + }); + + /// Compute specific volume first + TestEos->computeSpecVol(TArray, SArray, PArray); + Array2DReal SpecVol = TestEos->SpecVol; + + /// Compute Brunt-Vaisala frequency + TestEos->computeBruntVaisalaFreq(TArray, SArray, PArray, SpecVol); + Array2DReal BruntVaisalaFreq = TestEos->BruntVaisalaFreq; + auto BruntVaisalaFreqH = createHostMirrorCopy(BruntVaisalaFreq); + /// Check all array values against expected value + int numMismatches = 0; + if (!isApprox(BruntVaisalaFreqH(1, 1), TeosBVFExpValue, RTol)) { + numMismatches = 1; + } else { + numMismatches = 0; + } + // parallelReduce( + // "CheckSpecVolMatrix-Teos", {Mesh->NCellsAll, NVertLayers}, + // KOKKOS_LAMBDA(int i, int j, int &localCount) { + // if (!isApprox(BruntVaisalaFreq(i, j), TeosBVFExpValue, RTol)) { + // localCount++; + // } + // }, + // numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("EosTest: TEOS BruntVaisalaFreq isApprox FAIL, " + "expected {}, got {} with {} mismatches", + TeosBVFExpValue, BruntVaisalaFreqH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("EosTest BruntVaisalaFreqCalc TEOS-10: PASS"); + } + + return Err; +} + /// Finalize and clean up all test infrastructure void finalizeEosTest() { HorzMesh::clear(); @@ -302,12 +508,12 @@ int checkValueGswcSpecVol() { /// Get specific volume from GSW-C library double SpecVol = gsw_specvol(Sa, Ct, P); /// Check the value against the expected TEOS-10 value - bool Check = isApprox(SpecVol, TeosExpValue, RTol); + bool Check = isApprox(SpecVol, TeosSVExpValue, RTol); if (!Check) { Err++; LOG_ERROR( "checkValueGswcSpecVol: SpecVol isApprox FAIL, expected {}, got {}", - TeosExpValue, SpecVol); + TeosSVExpValue, SpecVol); } if (Err == 0) { LOG_INFO("checkValueGswcSpecVol: PASS"); @@ -315,15 +521,54 @@ int checkValueGswcSpecVol() { return Err; } -// the main test (all in one to have the same log) +/// Test that the external GSW-C library returns the expected N2 +int checkValueGswcN2() { + int Err = 0; + const Real RTol = 1e-10; + + // Number of intervals (nz) + int Nz = 2; + + // Input arrays: length nz+1 + double Salt[4] = {Sa - 1.0, Sa, Sa + 1.0}; // Absolute Salinity (g/kg) + double Temp[4] = {Ct + 15.0, Ct + 10.0, + Ct + 5.0}; // Conservative Temperature (deg C) + double Press[4] = {P, P + 1.0, P + 2.0}; // Pressure (dbar) + + // Latitude (degrees north) + double Latitude[4] = {0.0, 0.0, 0.0}; + + // Output arrays: length nz + double N2[Nz]; // Brunt–Väisälä frequency squared + double PMid[Nz]; // Midpoint pressure + + /// Get specific volume from GSW-C library + gsw_nsquared(Salt, Temp, Press, Latitude, Nz, N2, PMid); + + /// Check the value against the expected TEOS-10 value + bool Check = isApprox(N2[0], GswBVFExpValue, RTol); + if (!Check) { + Err++; + LOG_ERROR("checkValueGswcN2: N2 isApprox FAIL, expected {}, got {}", + GswBVFExpValue, N2[0]); + } + if (Err == 0) { + LOG_INFO("checkValueGswcN2: PASS"); + } + return Err; +} + +// the main tests (all in one to have the same log): // Single value test: // --> test calls the external GSW-C library // and compares the specific volume to the published value // Full array tests: // --> one tests the value on a Eos with linear option // --> next checks the value on a Eos with linear displaced option +// --> next checks the value of the linear Brunt Vaisala Freq. calculation // --> next checks the value on a Eos with TEOS-10 option // --> next checks the value on a Eos with TEOS-10 displaced option +// --> last checks the value of the TOES-10 Brunt Vaisala Freq. calculation int eosTest(const std::string &MeshFile = "OmegaMesh.nc") { int Err = initEosTest(MeshFile); if (Err != 0) { @@ -333,12 +578,15 @@ int eosTest(const std::string &MeshFile = "OmegaMesh.nc") { LOG_INFO("Single value checks:"); Err += checkValueGswcSpecVol(); + Err += checkValueGswcN2(); LOG_INFO("Full array checks:"); Err += testEosLinear(); Err += testEosLinearDisplaced(); + Err += testBruntVaisalaFreqLinear(); Err += testEosTeos10(); Err += testEosTeos10Displaced(); + Err += testBruntVaisalaFreqTeos10(); if (Err == 0) { LOG_INFO("EosTest: Successful completion"); diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp new file mode 100644 index 000000000000..2e49ba60b585 --- /dev/null +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -0,0 +1,658 @@ +//===-- Test driver for OMEGA Vertical Mixing Coefficients ---------------*- C++ +//-*-===/ +// +/// \file +/// \brief Test driver for OMEGA Vertical Mixing Coefficients +/// +/// This driver tests that VertMix can be called and returns expected values +/// of diffusivity, viscosity and Brunt-Vaisala frequency +/// +//===-----------------------------------------------------------------------===/ + +#include "VertMix.h" +#include "Config.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "Dimension.h" +#include "IO.h" +#include "Logging.h" +#include "MachEnv.h" +#include "OceanTestCommon.h" +#include "OmegaKokkos.h" +#include "Pacer.h" +#include "mpi.h" + +// added for debug +#include "AuxiliaryState.h" +#include "Field.h" +#include "Halo.h" +#include "HorzMesh.h" +#include "VertCoord.h" + +using namespace OMEGA; + +/// Test constants and expected values +constexpr int NVertLayers = 60; + +/// Values to test against +const Real VertDiffExpValueN = + 1.0391936989129011; // Expected value for diffusivity for positive BVF +const Real VertViscExpValueN = + 1.0198269656595984; // Expected value for viscosity for positive BVF +const Real VertDiffExpValueP = + 0.0015017457509864886; // Expected value for diffusivity for negative BVF +const Real VertViscExpValueP = + 0.002332474675614262; // 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 = + 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 + +/// 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 + +/// The initialization routine for VertMix testing. It calls various +/// init routines, including the creation of the default decomposition. +I4 initVertMixTest(const std::string &mesh) { + + I4 Err = 0; + + /// Initialize the Machine Environment class - this also creates + /// the default MachEnv. Then retrieve the default environment and + /// some needed data members. + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + /// Initialize logging + initLogging(DefEnv); + + /// Open and read config file + Config("Omega"); + Config::readAll("omega.yml"); + + // Initialize the IO system + IO::init(DefComm); + + /// Initialize decomposition + Decomp::init(); + + // Initialize the vertical coordinate (phase 1) + VertCoord::init1(); + + /// Initialize mesh + HorzMesh::init(); + + /// Initialize VertMix + VertMix::init(); + + /// Retrieve VertMix + VertMix *DefVertMix = VertMix::getInstance(); + if (DefVertMix) { + LOG_INFO("VertMixTest: VertMix retrieval PASS"); + } else { + Err++; + LOG_INFO("VertMixTest: VertMix retrieval FAIL"); + return -1; + } + + return Err; +} + +int testBackVertMix() { + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + /// Get VertMix instance to test + VertMix *TestVertMix = VertMix::getInstance(); + + /// Create and fill ocean state arrays + Array2DReal NormalVelEdge = + Array2DReal("NormalVelEdge", Mesh->NEdgesSize, VCoord->NVertLayers); + Array2DReal TangVelEdge = + Array2DReal("TangVelEdge", Mesh->NEdgesSize, VCoord->NVertLayers); + Array2DReal BruntVaisalaFreqCell = Array2DReal( + "BruntVaisalaFreqCell", Mesh->NCellsSize, VCoord->NVertLayers); + /// Use Kokkos::deep_copy to fill the entire view with the ref value + deepCopy(NormalVelEdge, NV); + deepCopy(TangVelEdge, TV); + deepCopy(BruntVaisalaFreqCell, BVFN); + deepCopy(TestVertMix->VertDiff, 0.0); + deepCopy(TestVertMix->VertVisc, 0.0); + + parallelFor( + "populateArrays", {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { VCoord->ZMid(ICell, K) = -K; }); + + parallelFor( + "populateArrays", {Mesh->NEdgesAll, VCoord->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; + }); + + /// Compute only background vertical viscosity and diffusivity + TestVertMix->BackDiff = 1.0e-5; + TestVertMix->BackVisc = 1.0e-4; + TestVertMix->ComputeVertMixConv.Enabled = false; + TestVertMix->ComputeVertMixShear.Enabled = false; + TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, + BruntVaisalaFreqCell); + Array2DReal BackVertVisc = TestVertMix->VertVisc; + Array2DReal BackVertDiff = TestVertMix->VertDiff; + + /// Check total Visc against linear addition of components + int numMismatches = 0; + parallelReduce( + "CheckVertMixMatrix-BackgroundVisc", + {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (j == 0) { + // Surface layer should be zero + if (!isApprox(BackVertVisc(i, j), 0.0_Real, RTol)) { + localCount++; + } + } else { + if (!isApprox(BackVertVisc(i, j), VertViscBackExp, RTol)) { + localCount++; + } + } + return; + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixBack: VertVisc isApprox FAIL, " + "expected {}, got {} mismatches", + VertViscBackExp, numMismatches); + } + if (Err == 0) { + LOG_INFO("VertMixTest TestVertMixBack-Visc: PASS"); + } + + /// Check total Diff against linear addition of components + numMismatches = 0; + parallelReduce( + "CheckVertMixMatrix-BackgroundDiff", + {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (j == 0) { + // Surface layer should be zero + if (!isApprox(BackVertDiff(i, j), 0.0_Real, RTol)) { + localCount++; + } + } else { + if (!isApprox(BackVertDiff(i, j), VertDiffBackExp, RTol)) { + localCount++; + } + } + return; + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixBack: VertDiff isApprox FAIL, " + "expected {}, got {} mismatches", + VertDiffBackExp, numMismatches); + } + if (Err == 0) { + LOG_INFO("VertMixTest TestVertMixBack-Diff: PASS"); + } + + return Err; +} + +int testConvVertMix() { + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + /// Get VertMix instance to test + VertMix *TestVertMix = VertMix::getInstance(); + + /// Create and fill ocean state arrays + Array2DReal NormalVelEdge = + Array2DReal("NormalVelEdge", Mesh->NEdgesAll, VCoord->NVertLayers); + Array2DReal TangVelEdge = + Array2DReal("TangVelEdge", Mesh->NEdgesAll, VCoord->NVertLayers); + Array2DReal BruntVaisalaFreqCell = Array2DReal( + "BruntVaisalaFreqCell", Mesh->NCellsAll, VCoord->NVertLayers); + /// Use Kokkos::deep_copy to fill the entire view with the ref value + deepCopy(NormalVelEdge, NV); + deepCopy(TangVelEdge, TV); + deepCopy(BruntVaisalaFreqCell, BVFN); + deepCopy(TestVertMix->VertDiff, 0.0); + deepCopy(TestVertMix->VertVisc, 0.0); + + parallelFor( + "populateArrays", {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { VCoord->ZMid(ICell, K) = -K; }); + + parallelFor( + "populateArrays", {Mesh->NEdgesAll, VCoord->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; + }); + + /// 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, + BruntVaisalaFreqCell); + Array2DReal ConvVertVisc = TestVertMix->VertVisc; + Array2DReal ConvVertDiff = TestVertMix->VertDiff; + + /// Check total Visc against linear addition of components + int numMismatches = 0; + parallelReduce( + "CheckVertMixMatrix-ConvectiveVisc", + {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (j == 0) { + // Surface layer should be zero + if (!isApprox(ConvVertVisc(i, j), 0.0_Real, RTol)) { + localCount++; + } + } else { + if (!isApprox(ConvVertVisc(i, j), VertDiffConvExp, RTol)) { + localCount++; + } + } + return; + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixConv: VertVisc isApprox FAIL, " + "expected {}, got {} mismatches", + VertDiffConvExp, numMismatches); + } + if (Err == 0) { + LOG_INFO("VertMixTest TestVertMixConv-Visc: PASS"); + } + + /// Check total Diff against linear addition of components + numMismatches = 0; + parallelReduce( + "CheckVertMixMatrix-ConvectiveDiff", + {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(int i, int j, int &localCount) { + if (j == 0) { + // Surface layer should be zero + if (!isApprox(ConvVertDiff(i, j), 0.0_Real, RTol)) { + localCount++; + } + } else { + if (!isApprox(ConvVertDiff(i, j), VertDiffConvExp, RTol)) { + localCount++; + } + } + return; + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixConv: VertDiff isApprox FAIL, " + "expected {}, got {} mismatches", + VertDiffConvExp, numMismatches); + } + if (Err == 0) { + LOG_INFO("VertMixTest TestVertMixConv-Diff: PASS"); + } + + return Err; +} + +int testShearVertMix() { + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + /// Get VertMix instance to test + VertMix *TestVertMix = VertMix::getInstance(); + + /// Create and fill ocean state arrays + Array2DReal NormalVelEdge = + Array2DReal("NormalVelEdge", Mesh->NEdgesAll, VCoord->NVertLayers); + Array2DReal TangVelEdge = + Array2DReal("TangVelEdge", Mesh->NEdgesAll, VCoord->NVertLayers); + Array2DReal BruntVaisalaFreqCell = Array2DReal( + "BruntVaisalaFreqCell", Mesh->NCellsAll, VCoord->NVertLayers); + /// Use Kokkos::deep_copy to fill the entire view with the ref value + deepCopy(NormalVelEdge, NV); + deepCopy(TangVelEdge, TV); + deepCopy(BruntVaisalaFreqCell, BVFN); + deepCopy(TestVertMix->VertDiff, 0.0); + deepCopy(TestVertMix->VertVisc, 0.0); + + parallelFor( + "populateArrays", {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + VCoord->ZMid(ICell, K) = -K; + Mesh->NEdgesOnCell(ICell) = 5; + Mesh->AreaCell(ICell) = 3.6e10_Real; + }); + + parallelFor( + "populateArrays", {Mesh->NEdgesAll, VCoord->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; + Mesh->DcEdge(IEdge) = 2.0e5_Real; + Mesh->DvEdge(IEdge) = 1.45e5_Real; + }); + + /// 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, + BruntVaisalaFreqCell); + Array2DReal ShearVertVisc = TestVertMix->VertVisc; + Array2DReal ShearVertDiff = TestVertMix->VertDiff; + + /// Check total Visc against linear addition of components + int numMismatches = 0; + Kokkos::fence(); + auto ShearVertViscH = createHostMirrorCopy(ShearVertVisc); + parallelReduce( + "CheckVertMixMatrix-ShearVisc", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int i, int &localCount) { + // Surface layer should be zero + if (!isApprox(ShearVertVisc(i, 0), 0.0_Real, RTol)) { + localCount++; + } + if (!isApprox(ShearVertVisc(i, 1), VertViscShearExp, RTol)) { + localCount++; + } + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixShear: VertVisc isApprox FAIL, " + "expected {}, got {} with {} mismatches", + VertViscShearExp, ShearVertViscH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("VertMixTest TestVertMixShear-Visc: PASS"); + } + + /// Check total Diff against linear addition of components + numMismatches = 0; + auto ShearVertDiffH = createHostMirrorCopy(ShearVertDiff); + parallelReduce( + "CheckVertMixMatrix-ShearDiff", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int i, int &localCount) { + // Surface layer should be zero + if (!isApprox(ShearVertDiff(i, 0), 0.0_Real, RTol)) { + localCount++; + } + if (!isApprox(ShearVertDiff(i, 1), VertDiffShearExp, RTol)) { + localCount++; + } + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixShear: VertDiff isApprox FAIL, " + "expected {}, got {} with {} mismatches", + VertDiffShearExp, ShearVertDiffH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("VertMixTest TestVertMixShear-Diff: PASS"); + } + + return Err; +} + +/// Test vertical mixing coefficients calculation for all cells/layers +int testTotalVertMix() { + int Err = 0; + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + VCoord->NVertLayers = NVertLayers; + /// Get VertMix instance to test + VertMix *TestVertMix = VertMix::getInstance(); + + /// Create and fill ocean state arrays + Array2DReal NormalVelEdge = + Array2DReal("NormalVelEdge", Mesh->NEdgesAll, VCoord->NVertLayers); + Array2DReal TangVelEdge = + Array2DReal("TangVelEdge", Mesh->NEdgesAll, VCoord->NVertLayers); + Array2DReal BruntVaisalaFreqCell = Array2DReal( + "BruntVaisalaFreqCell", Mesh->NCellsAll, VCoord->NVertLayers); + /// Use Kokkos::deep_copy to fill the entire view with the ref value + deepCopy(NormalVelEdge, NV); + deepCopy(TangVelEdge, TV); + + // Test with positive BVF first + deepCopy(BruntVaisalaFreqCell, BVFP); + deepCopy(TestVertMix->VertDiff, 0.0); + deepCopy(TestVertMix->VertVisc, 0.0); + + parallelFor( + "populateArrays", {Mesh->NCellsAll, VCoord->NVertLayers}, + KOKKOS_LAMBDA(I4 ICell, I4 K) { + VCoord->ZMid(ICell, K) = -K; + Mesh->NEdgesOnCell(ICell) = 5; + Mesh->AreaCell(ICell) = 3.6e10_Real; + }); + + parallelFor( + "populateArrays", {Mesh->NEdgesAll, VCoord->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; + Mesh->DcEdge(IEdge) = 2.0e5_Real; + Mesh->DvEdge(IEdge) = 1.45e5_Real; + }); + + /// Compute vertical viscosity and diffusivity + TestVertMix->BackDiff = 1.0e-5; + TestVertMix->BackVisc = 1.0e-4; + TestVertMix->ComputeVertMixConv.Enabled = true; + TestVertMix->ComputeVertMixShear.Enabled = true; + TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, + BruntVaisalaFreqCell); + Array2DReal VertDiffP = TestVertMix->VertDiff; + Array2DReal VertViscP = TestVertMix->VertVisc; + + /// Check all VertDiff array values against expected value + int numMismatches = 0; + auto VertDiffPH = createHostMirrorCopy(VertDiffP); + parallelReduce( + "CheckVertMixMatrix-TotalPosDiff", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int i, int &localCount) { + // Surface layer should be zero + if (!isApprox(VertDiffP(i, 0), 0.0_Real, RTol)) { + localCount++; + } + if (!isApprox(VertDiffP(i, 1), VertDiffExpValueP, RTol)) { + localCount++; + } + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixTotal: VertDiffPositive isApprox FAIL, " + "expected {}, got {} with {} mismatches", + VertDiffExpValueP, VertDiffPH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("TestVertMixTotal VertDiffPositive: PASS"); + } + + /// Check all VertVisc array values against expected value + numMismatches = 0; + auto VertViscPH = createHostMirrorCopy(VertViscP); + parallelReduce( + "CheckVertMixMatrix-TotalPosVisc", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int i, int &localCount) { + // Surface layer should be zero + if (!isApprox(VertViscP(i, 0), 0.0_Real, RTol)) { + localCount++; + } + if (!isApprox(VertViscP(i, 1), VertViscExpValueP, RTol)) { + localCount++; + } + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMixTotal: VertViscPositive isApprox FAIL, " + "expected {}, got {} with {} mismatches", + VertViscExpValueP, VertViscPH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("TestVertMixTotal VertViscPositive: PASS"); + } + + // Now test with negative BVF + deepCopy(BruntVaisalaFreqCell, BVFN); + deepCopy(TestVertMix->VertDiff, 0.0); + deepCopy(TestVertMix->VertVisc, 0.0); + + /// Compute vertical viscosity and diffusivity + TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, + BruntVaisalaFreqCell); + Array2DReal VertDiffN = TestVertMix->VertDiff; + Array2DReal VertViscN = TestVertMix->VertVisc; + + /// Check all VertDiff array values against expected value + numMismatches = 0; + auto VertDiffNH = createHostMirrorCopy(VertDiffN); + parallelReduce( + "CheckVertMixMatrix-TotalNegDiff", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int i, int &localCount) { + // Surface layer should be zero + if (!isApprox(VertDiffN(i, 0), 0.0_Real, RTol)) { + localCount++; + } + if (!isApprox(VertDiffN(i, 1), VertDiffExpValueN, RTol)) { + localCount++; + } + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMix: VertDiffNegative isApprox FAIL, " + "expected {}, got {} with {} mismatches", + VertDiffExpValueN, VertDiffNH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("TestVertMix VertDiffNegative: PASS"); + } + + /// Check all VertVisc array values against expected value + numMismatches = 0; + auto VertViscNH = createHostMirrorCopy(VertViscN); + parallelReduce( + "CheckVertMixMatrix-TotalNegVisc", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int i, int &localCount) { + // Surface layer should be zero + if (!isApprox(VertViscN(i, 0), 0.0_Real, RTol)) { + localCount++; + } + if (!isApprox(VertViscN(i, 1), VertViscExpValueN, RTol)) { + localCount++; + } + }, + numMismatches); + + if (numMismatches != 0) { + Err++; + LOG_ERROR("TestVertMix: VertViscNegative isApprox FAIL, " + "expected {}, got {} with {} mismatches", + VertViscExpValueN, VertViscNH(1, 1), numMismatches); + } + if (Err == 0) { + LOG_INFO("TestVertMix VertViscNegative : PASS"); + } + + return Err; +} + +/// Finalize and clean up all test infrastructure +void finalizeVertMixTest() { + HorzMesh::clear(); + VertCoord::clear(); + Decomp::clear(); + Field::clear(); + Dimension::clear(); + MachEnv::removeAll(); +} + +// the main tests (all in one to have the same log): +// --> one tests the vertical diffusivity and viscosity +// with only background on +// --> next tests the vertical diffusivity and viscosity +// with only convective on +// --> next tests the vertical diffusivity and viscosity +// with only shear on +// --> next tests the linear superposition of the +// background, convective, and shear contributions +int vertMixTest(const std::string &MeshFile = "OmegaMesh.nc") { + int Err = initVertMixTest(MeshFile); + if (Err != 0) { + LOG_CRITICAL("VertMix: Error initializing"); + } + const auto &Mesh = HorzMesh::getDefault(); + + Err += testBackVertMix(); + Err += testConvVertMix(); + Err += testShearVertMix(); + Err += testTotalVertMix(); + + if (Err == 0) { + LOG_INFO("VertMix: Successful completion"); + } + finalizeVertMixTest(); + return Err; +} + +// The test driver for VertMix testing +int main(int argc, char *argv[]) { + + int RetVal = 0; + + MPI_Init(&argc, &argv); + Kokkos::initialize(argc, argv); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + + RetVal += vertMixTest(); + + VertMix::destroyInstance(); + Kokkos::finalize(); + MPI_Finalize(); + + if (RetVal >= 256) + RetVal = 255; + + return RetVal; + +} // end of main +//===-----------------------------------------------------------------------===/