From d55b4e4d5570643f716afb3f2d334d82b7a9f836 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 2 Feb 2026 16:37:19 -0800 Subject: [PATCH 01/19] changes shear mixing to LMD94 formulation and adds ri smoothing --- components/omega/configs/Default.yml | 2 +- .../omega/doc/devGuide/VerticalMixingCoeff.md | 14 +- .../doc/userGuide/VerticalMixingCoeff.md | 15 +- components/omega/src/ocn/VertMix.cpp | 98 ++++++++++-- components/omega/test/ocn/VertMixTest.cpp | 142 ++++++++++++++++++ 5 files changed, 245 insertions(+), 26 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index ce0de75afebf..bc526ba3c738 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -92,7 +92,7 @@ Omega: BaseShearValue: 0.005 RiCrit: 0.7 Exponent: 3.0 - RiSmoothLoops: 2 + RiSmoothLoops: 3 IOStreams: HorzMeshIn: UsePointerFile: false diff --git a/components/omega/doc/devGuide/VerticalMixingCoeff.md b/components/omega/doc/devGuide/VerticalMixingCoeff.md index a59151ef27cc..5ff40bc0505d 100644 --- a/components/omega/doc/devGuide/VerticalMixingCoeff.md +++ b/components/omega/doc/devGuide/VerticalMixingCoeff.md @@ -40,7 +40,7 @@ VertMix: BaseShearValue: 0.005 RiCrit: 0.7 Exponent: 3.0 - RiSmoothLoops: 2 + RiSmoothLoops: 3 ``` ## Class Structure @@ -63,12 +63,12 @@ VertMix: - `ConvDiff`: Convective mixing coefficient (m²/s; Default: 1.0) - `ConvTriggerBVF`: Trigger threshold for convective mixing (Default: 0.0) -3. Shear 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) +3. Shear Mixing: + - `EnableShearMix`: Flag to enable/disable shear mixing (Default: True) + - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) + - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) + - `ShearExponent`: Exponent parameter number for the LMB94 formulation (Default: 3.0) + - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 3) ## Core Functionality (Vertical Mixing Coefficient Calculation) diff --git a/components/omega/doc/userGuide/VerticalMixingCoeff.md b/components/omega/doc/userGuide/VerticalMixingCoeff.md index 2f97a169d8fe..a530b6781acc 100644 --- a/components/omega/doc/userGuide/VerticalMixingCoeff.md +++ b/components/omega/doc/userGuide/VerticalMixingCoeff.md @@ -6,8 +6,9 @@ The vertical mixing module in Omega handles the parameterization of unresolved v processes in the ocean. It calculates vertical diffusivity and viscosity coefficients that determine how properties (like momentum, heat, salt, and biogeochemical tracers) mix vertically in the ocean model. Currently, Omega offers three different mixing processes within the water column: (1) a constant -background mixing value, (2) a convective instability mixing value, and (3) a Richardson number-dependent shear-instability-driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 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. +background mixing value, (2) a convective instability mixing value, and (3) a Richardson number +dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These are linearly additive and are describe a bit +more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear mixing calculation. The user-configurable options are the following parameters in the yaml configuration file: @@ -25,14 +26,14 @@ VertMix: 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 + RiSmoothLoops: 3 # Number of Richardson number smoothing loops ``` ## Vertical Mixing Processes/Types ### 1. Background Mixing -A constant background mixing value that represents small-scale mixing processes not explicitly resolved 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). +A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior. ### 2. Convective Mixing @@ -50,14 +51,14 @@ This is different than some current implementations (i.e. in MPAS-Ocean and the ### 3. Shear-Instability-Driven Mixing -Mixing induced by vertical pseudo-velocity shear, implemented using the LMD94 scheme, through the gradient Richardson number (ratio of buoyancy to shear). +Mixing induced by vertical velocity shear, implemented using the LMD94 scheme, through the gradient Richardson number (ratio of buoyancy to shear). $$ \nu_{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 +\nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 < Ri_g < Ri_{crit}\\ +0.0 \quad \text{ if } Ri_{crit} < Ri_g \end{cases} $$ diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 815dab5e6814..72690a1eec98 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -201,6 +201,7 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, OMEGA_SCOPE(LocBackVisc, BackVisc); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(NVertLayers, VCoord->NVertLayers); /// First, initialize VertDiff and VertVisc to background values parallelForOuter( @@ -237,22 +238,97 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, + KChunk, BruntVaisalaFreqSq); LocComputeGradRichardsonNum( LocGradRichNum, ICell, KChunk, NormalVelocity, TangentialVelocity, BruntVaisalaFreqSq); }); + LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); + LocGradRichNum(ICell, NVertLayers) = + LocGradRichNum(ICell, NVertLayers - 1); + }); + deepCopy(GradRichNumSmoothed, LocGradRichNum); + for (int SmoothLoop = 0; + SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + if (SmoothLoop == 0) + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + LocGradRichNum); + else + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + GradRichNumSmoothed); + }); + }); + } + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); - 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); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, + KChunk, GradRichNumSmoothed); + }); + }); + } else if (LocComputeVertMixShear.Enabled) { + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeGradRichardsonNum( + LocGradRichNum, ICell, KChunk, NormalVelocity, + TangentialVelocity, BruntVaisalaFreqSq); + }); + LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); + LocGradRichNum(ICell, NVertLayers) = + LocGradRichNum(ICell, NVertLayers - 1); + }); + deepCopy(GradRichNumSmoothed, LocGradRichNum); + for (int SmoothLoop = 0; + SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + if (SmoothLoop == 0) + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + LocGradRichNum); + else + LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, + GradRichNumSmoothed); + }); + }); + } + parallelForOuter( + "VertMix-ConvPlusShear", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, + KChunk, GradRichNumSmoothed); }); }); /// Smooth Richardson number with 1-2-1 filter the number of times diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index d1bf97015798..c25338cbabc8 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -205,6 +205,9 @@ void testGradRichNum() { TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, BruntVaisalaFreqSqCell); + const auto &MinLayerCell = VCoord->MinLayerCell; + const auto &MaxLayerCell = VCoord->MaxLayerCell; + /// Check all array values against expected value int NumMismatches = 0; OMEGA_SCOPE(GradRichNum, TestVertMix->GradRichNum); @@ -259,6 +262,15 @@ void testGradRichNum() { // If test fails, print bad values and abort if (NumMismatches != 0) { + auto GradRichNumH = createHostMirrorCopy(GradRichNum); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (!isApprox(GradRichNumH(I, K), RiExpValue, RTol)) + LOG_ERROR("TestVertMix: GradRichNum Bad Value: " + "GradRichNum({},{}) = {}; Expected {}", + I, K, GradRichNumH(I, K), RiExpValue); + } + } ABORT_ERROR("TestVertMix: GradRichNum FAIL with {} bad values", NumMismatches); } else { @@ -352,6 +364,28 @@ void testOneTwoOneFilter() { // If test fails, print bad values and abort if (NumMismatches != 0) { + auto GradRichNumH = createHostMirrorCopy(GradRichNum); + auto GradRichNumSmoothedH = createHostMirrorCopy(GradRichNumSmoothed); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K >= 1 && K < NVertLayers - 1) { + // Interior layers should be smoothed to 0.0 + if (!isApprox(GradRichNumSmoothedH(I, K), 0.0, RTol)) + LOG_ERROR("TestVertMix: GradRichNumSmoothed Bad Value: " + "GradRichNumSmoothed({},{}) = {}; Expected {}", + I, K, GradRichNumSmoothedH(I, K), 0.0); + } else { + // Boundary layers (K==0 or K==NVertLayers-1) should be copied + // from input + if (!isApprox(GradRichNumSmoothedH(I, K), GradRichNumH(I, K), + RTol)) + LOG_ERROR("TestVertMix: GradRichNumSmoothed Bad Value: " + "GradRichNumSmoothed({},{}) = {}; Expected {}", + I, K, GradRichNumSmoothedH(I, K), + GradRichNumH(I, K)); + } + } + } ABORT_ERROR("TestVertMix: GradRichNumSmoothed FAIL with {} bad values", NumMismatches); } else { @@ -575,6 +609,30 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { + auto VertViscOutH = createHostMirrorCopy(VertViscOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } else if (K < 30) { + // Interior layers + if (!isApprox(VertViscOutH(I, K), VertConvExp, RTol)) + LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), VertConvExp); + } else { + // Interior layers + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } + } + } ABORT_ERROR("TestVertMixConv: VertVisc FAIL with {} bad values", NumMismatches); } else { @@ -614,6 +672,30 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { + auto VertDiffOutH = createHostMirrorCopy(VertDiffOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } else if (K < 30) { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), VertConvExp, RTol)) + LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), VertConvExp); + } else { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } + } + } ABORT_ERROR("TestVertMixConv: VertDiff FAIL with {} bad values", NumMismatches); } else { @@ -713,6 +795,36 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { + auto VertViscOutH = createHostMirrorCopy(VertViscOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } else if (K < 20) { + // Interior layers + if (!isApprox(VertViscOutH(I, K), VertShearBaseExp, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), VertShearBaseExp); + } else if (K >= 20 && K < 40) { + // Interior layers + if (!isApprox(VertViscOutH(I, K), VertShearExp, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), VertShearExp); + } else { + // Interior layers + if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " + "VertVisc({},{}) = {}; Expected {}", + I, K, VertViscOutH(I, K), 0.0_Real); + } + } + } ABORT_ERROR("TestVertMixShear: VertVisc FAIL with {} bad values", NumMismatches); } else { @@ -756,6 +868,36 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { + auto VertDiffOutH = createHostMirrorCopy(VertDiffOut); + for (int I = 0; I < NCellsAll; ++I) { + for (int K = 0; K < NVertLayers; ++K) { + if (K == 0) { + // Surface should be 0.0 + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } else if (K < 20) { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), VertShearBaseExp, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), VertShearBaseExp); + } else if (K >= 20 && K < 40) { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), VertShearExp, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), VertShearExp); + } else { + // Interior layers + if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) + LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " + "VertDiff({},{}) = {}; Expected {}", + I, K, VertDiffOutH(I, K), 0.0_Real); + } + } + } ABORT_ERROR("TestVertMixShear: VertDiff FAIL with {} bad values", NumMismatches); } else { From cc101033c44423d97da8dafa78f9133273ec7322 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Thu, 21 May 2026 09:52:13 -0700 Subject: [PATCH 02/19] changes gradrichnum, vertvisc, and vertdiff to nvertlayersp1 dim --- components/omega/src/ocn/VertMix.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 72690a1eec98..4c08d607bcb6 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -329,6 +329,18 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, Team, KRange, INNER_LAMBDA(int KChunk) { LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, KChunk, GradRichNumSmoothed); + + 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); }); }); /// Smooth Richardson number with 1-2-1 filter the number of times From aa2ce7cd9e08fe9af8a9dd906cb75fd2a2891020 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 4 May 2026 20:51:31 -0700 Subject: [PATCH 03/19] Add ALE projection time intervals to time steppers --- components/omega/src/timeStepping/RungeKutta4Stepper.h | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.h b/components/omega/src/timeStepping/RungeKutta4Stepper.h index 16efa1acd9ba..d4dbe8f8072b 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.h +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.h @@ -46,10 +46,7 @@ class RungeKutta4Stepper : public TimeStepper { Real RKB[NStages]; Real RKC[NStages]; - // Projection factors for each RK stage. These are the RKC coefficients - // shifted by one RK stage and represent the fraction of the timestep over - // which fields are projected from the old state to the state at the end of a - // RK stage. + // Projection factors Real RKProj[NStages]; }; From fdd426f693faa4840b44d026803e81948ea3a21d Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 11 May 2026 00:46:30 -0400 Subject: [PATCH 04/19] Implement vertical mixing tendency - Implemented implicit vertical mixing using TriDiagDiffSolver. - Added applyVertMixImplicit to VertMix. - Applied vertical mixing once at the end of each time step. - Note: Tangential velocity is temporarily handled in VertMix. --- components/omega/configs/Default.yml | 2 + components/omega/src/ocn/OceanInit.cpp | 88 ++++++ components/omega/src/ocn/OceanState.h | 10 +- components/omega/src/ocn/Tendencies.cpp | 36 ++- components/omega/src/ocn/Tendencies.h | 15 + components/omega/src/ocn/VertMix.cpp | 281 +++++++++++++++++- components/omega/src/ocn/VertMix.h | 140 +++++++++ .../timeStepping/ForwardBackwardStepper.cpp | 10 + .../src/timeStepping/RungeKutta2Stepper.cpp | 10 + .../src/timeStepping/RungeKutta4Stepper.cpp | 11 + components/omega/test/infra/IOStreamTest.cpp | 5 + components/omega/test/ocn/StateTest.cpp | 5 + components/omega/test/ocn/TendenciesTest.cpp | 6 +- .../test/timeStepping/TimeStepperTest.cpp | 10 +- 14 files changed, 617 insertions(+), 12 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index bc526ba3c738..92299ae6f259 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -66,6 +66,8 @@ Omega: VelocityVertAdvTendencyEnable: true TracerVertAdvTendencyEnable: true PressureGradTendencyEnable: true + VelVertMixTendencyEnable: true + TracerVertMixTendencyEnable: true ManufacturedSolution: WavelengthX: 5.0e6 WavelengthY: 4.33013e6 diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index ceabc0ede033..1348f197c747 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -31,6 +31,7 @@ #include "Tracers.h" #include "VertAdv.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" @@ -106,6 +107,44 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator TimeStepper *DefStepper = TimeStepper::getDefault(); Clock *ModelClock = DefStepper->getClock(); + // Initialize IOStreams - this does not yet validate the contents + // of each file, only creates streams from Config + IOStream::init(ModelClock); + + IO::init(Comm); + Field::init(ModelClock); + Decomp::init(); + + Err = Halo::init(); + if (Err != 0) { + ABORT_ERROR("ocnInit: Error initializing default halo"); + } + + HorzMesh::init(); + VertCoord::init(); + Tracers::init(); + VertAdv::init(); + AuxiliaryState::init(); + Eos::init(); + PressureGrad::init(); + VertMix::init(); + Tendencies::init(); + + // Validate SurfaceTracerRestoring configuration + Tendencies *DefTend = Tendencies::getDefault(); + if (DefTend->SurfaceTracerRestoring.Enabled && + DefTend->SurfaceTracerRestoring.NTracersToRestore == 0) { + ABORT_ERROR("OceanInit: SurfaceTracerRestoring is enabled but " + "TracersToRestore is empty"); + } + + TimeStepper::init2(); + + Err = OceanState::init(); + if (Err != 0) { + ABORT_ERROR("ocnInit: Error initializing default state"); + } + // Now that all fields have been defined, validate all the streams // contents bool StreamsValid = IOStream::validateAll(); @@ -164,6 +203,55 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator } Tracers::copyToHost(CurTimeLevel); +//////////////////////////////////////////////// +//////////////////////////////////////////////// + + // For initial GeometricThick -> PseudoThick + HorzMesh *Mesh = HorzMesh::getDefault(); + VertCoord *VCoord = VertCoord::getDefault(); + Eos *EosInstance = Eos::getInstance(); + + // get temperature and salinity + I4 ConservTempIdx; + I4 AbsSalinityIdx; + Array3DReal TracerArray = Tracers::getAll(0); + + Tracers::getIndex(ConservTempIdx, "Temperature"); + Tracers::getIndex(AbsSalinityIdx, "Salinity"); + + const auto ConservTemp = + Kokkos::subview(TracerArray, ConservTempIdx, Kokkos::ALL, Kokkos::ALL); + const auto AbsSalinity = + Kokkos::subview(TracerArray, AbsSalinityIdx, Kokkos::ALL, Kokkos::ALL); + + Array2DReal PressureMidDbar("PressureMidDbar", VCoord->PressureMid.layout()); + deepCopy(PressureMidDbar, 0.0); + EosInstance->computeSpecVol(ConservTemp, AbsSalinity, PressureMidDbar); + const auto &SpecVol = EosInstance->SpecVol; + + // get layer thickness + Array2DReal LayerThickCell = DefState->getLayerThickness(0); + OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocNCellsAll, Mesh->NCellsAll); + + parallelForOuter( + "GtoP", {LocNCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = LocMinLayerCell(ICell); + const int KMax = LocMaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + LayerThickCell(ICell,K) = LayerThickCell(ICell,K) / (SpecVol(ICell,K)*1026.0); + + }); + }); + +//////////////////////////////////////////////// +////////////////////////////////////////////////// + return Err; } // end ocnInit diff --git a/components/omega/src/ocn/OceanState.h b/components/omega/src/ocn/OceanState.h index 3632449ccc32..59f6845c3305 100644 --- a/components/omega/src/ocn/OceanState.h +++ b/components/omega/src/ocn/OceanState.h @@ -48,11 +48,6 @@ class OceanState { OceanState(const OceanState &) = delete; OceanState(OceanState &&) = delete; - // Current time index - // this index is circular so that it returns to index 0 - // if it is over max index - I4 CurTimeIndex; ///< Time dimension array index for current level - /// Get the current time level index associated with a time level I4 getTimeIndex(const I4 TimeLevel) const; @@ -63,6 +58,11 @@ class OceanState { std::string Name; + // Current time index + // this index is circular so that it returns to index 0 + // if it is over max index + I4 CurTimeIndex; ///< Time dimension array index for current level + // Sizes and global IDs // Note that all sizes are actual counts (1-based) so that loop extents // should always use the 0:NCellsXX-1 form. diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index ba0df5dd692a..69c4fa22b27d 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -20,6 +20,7 @@ #include "TimeStepper.h" #include "Tracers.h" #include "VertAdv.h" +#include "VertMix.h" #include namespace OMEGA { @@ -39,6 +40,7 @@ void Tendencies::init() { TimeStepper *DefTimeStepper = TimeStepper::getDefault(); Eos *DefEos = Eos::getInstance(); PressureGrad *DefPGrad = PressureGrad::getDefault(); + VertMix *DefVertMix = VertMix::getInstance(); I4 NTracers = Tracers::getNumTracers(); @@ -80,7 +82,8 @@ void Tendencies::init() { // Ceate default tendencies Tendencies::DefaultTendencies = create( "Default", DefHorzMesh, DefVertCoord, DefVertAdv, DefPGrad, DefEos, - NTracers, TimeStep, &TendConfig, CustomThickTend, CustomVelTend); + DefVertMix, NTracers, TimeStep, &TendConfig, CustomThickTend, + CustomVelTend); DefaultTendencies->readConfig(OmegaConfig); @@ -317,6 +320,17 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config HostArray1DI4(TracerIdsToRestoreVec.data(), TracerIdsToRestoreVec.size())); } + + Err += TendConfig.get("VelVertMixTendencyEnable", + this->VMix->VelVertMixSetup.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: VelVertMixTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("TracerVertMixTendencyEnable", + this->VMix->TracerVertMixSetup.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: TracerVertMixTendencyEnable not found in TendConfig"); + } //------------------------------------------------------------------------------ @@ -380,6 +394,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VertAdv *VAdv, ///< [in] Vertical advection PressureGrad *PGrad, ///< [in] Pressure gradient Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStepIn, ///< [in] Time step Config *Options, ///< [in] Configuration options @@ -393,7 +408,8 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), - CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad) { + CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad), + VMix(VMix) { // Tendency arrays PseudoThicknessTend = Array2DReal("PseudoThicknessTend", Mesh->NCellsSize, @@ -418,13 +434,27 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VertAdv *VAdv, ///< [in] Vertical advection PressureGrad *PGrad, ///< [in] Pressure gradient Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStepIn, ///< [in] Time step Config *Options) ///< [in] Configuration options - : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, NTracersIn, + : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, VMix, NTracersIn, TimeStepIn, Options, CustomTendencyType{}, CustomTendencyType{}) {} +Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies + const HorzMesh *Mesh, ///< [in] Horizontal mesh + VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection + PressureGrad *PGrad, ///< [in] Pressure gradient + Eos *EqState, ///< [in] Equation of state + int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStepIn, ///< [in] Time step + Config *Options) ///< [in] Configuration options + : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, VertMix::getInstance(), + NTracersIn, TimeStepIn, Options, CustomTendencyType{}, + CustomTendencyType{}) {} + //------------------------------------------------------------------------------ // Compute tendencies for the pseudo-thickness equation void Tendencies::computePseudoThicknessTendenciesOnly( diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 5d26ca041458..d2aa85ffd2ef 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -41,6 +41,7 @@ #include "TimeMgr.h" #include "VertAdv.h" #include "VertCoord.h" +#include "VertMix.h" #include #include @@ -169,6 +170,7 @@ class Tendencies { VertAdv *VAdv, ///< [in] Vertical advection PressureGrad *PGrad, ///< [in] Pressure gradient Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStep, ///< [in] Time step Config *Options, ///< [in] Configuration options @@ -186,6 +188,18 @@ class Tendencies { Config *Options ///< [in] Configuration options ); + Tendencies(const std::string &Name, ///< [in] Name for tendencies + const HorzMesh *Mesh, ///< [in] Horizontal mesh + VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection + PressureGrad *PGrad, ///< [in] Pressure gradient + Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing + int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStep, ///< [in] Time step + Config *Options ///< [in] Configuration options + ); + void defineFields(); // forbid copy and move construction @@ -199,6 +213,7 @@ class Tendencies { CustomTendencyType CustomVelocityTend; Eos *EqState; ///< Pointer to equation of state PressureGrad *PGrad; ///< Pointer to pressure gradient + VertMix *VMix; ///< Pointer to vertical mixing I4 NTracers; ///< Number of tracers TimeInterval TimeStep; ///< Time step diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 4c08d607bcb6..c6e6749d3866 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -12,7 +12,12 @@ #include "VertMix.h" #include "DataTypes.h" +#include "Eos.h" +#include "GlobalConstants.h" #include "HorzMesh.h" +#include "HorzOperators.h" +#include "TriDiagSolvers.h" +#include "TimeStepper.h" namespace OMEGA { @@ -36,6 +41,19 @@ GradRichardsonNum::GradRichardsonNum(const HorzMesh *Mesh, OneTwoOneFilter::OneTwoOneFilter(const VertCoord *VCoord) : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} +VelVertMixSetupOnEdge::VelVertMixSetupOnEdge(const HorzMesh *Mesh, + const VertCoord *VCoord) + : Enabled(false), LocRhoSw(RhoSw), NVertLayers(VCoord->NVertLayers), + CellsOnEdge(Mesh->CellsOnEdge), EdgeMask(VCoord->EdgeMask), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} + +TracerVertMixSetupOnCell::TracerVertMixSetupOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) + : Enabled(false), LocRhoSw(RhoSw), NVertLayers(VCoord->NVertLayers), + 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 @@ -43,7 +61,8 @@ VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object ) : ComputeVertMixConv(VCoord), ComputeVertMixShear(VCoord), ComputeGradRichardsonNum(Mesh, VCoord), ComputeOneTwoOneFilter(VCoord), - Name(Name), Mesh(Mesh), VCoord(VCoord) { + Name(Name), Mesh(Mesh), VCoord(VCoord), VelVertMixSetup(Mesh, VCoord), + TracerVertMixSetup(Mesh, VCoord) { VertDiff = Array2DReal("VertDiff", Mesh->NCellsSize, VCoord->NVertLayersP1); VertVisc = Array2DReal("VertVisc", Mesh->NCellsSize, VCoord->NVertLayersP1); GradRichNum = @@ -51,6 +70,10 @@ VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object GradRichNumSmoothed = Array2DReal("GradRichNumSmoothed", Mesh->NCellsSize, VCoord->NVertLayersP1); + //TODO: Temporary handling of TangentialVelocity + TangentialVelocity = Array2DReal("TangentialVelocity", Mesh->NEdgesSize, + VCoord->NVertLayers); + defineFields(); } @@ -498,4 +521,260 @@ void VertMix::defineFields() { } // end defineIOFields + +// Dpply implicit velocity vertical mixing +void VertMix::applyVelVertMixImplicit( + OceanState *State, ///< [in] State variables + const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + int ThickTimeLevel, ///< [in] Time level + int VelTimeLevel ///< [in] Time level +) { + + OMEGA_SCOPE(LocNEdgesAll, Mesh->NEdgesAll); + OMEGA_SCOPE(LocVelVertMixSetup, VelVertMixSetup); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + + const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; + + // Compute velocity vertical mixing + if (LocVelVertMixSetup.Enabled) { + Pacer::start("Tend:velocityVertMix", 1); + + Eos *EosInstance = Eos::getInstance(); + VertMix *VertMixInstance = VertMix::getInstance(); + + if (!EosInstance) { + LOG_WARN("Eos has not been initialized. Skipping calculation of " + "VelVertMix tendency"); + } else if (!VertMixInstance) { + LOG_WARN("VertMix has not been initialized. Skipping calculation of " + "VelVertMix tendency"); + } else { + + // Obtain TimeStep + const auto *DefTimeStepper = TimeStepper::getDefault(); + const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); + R8 DT; + TimeStep.get(DT, TimeUnits::Seconds); + + const auto &SpecVol = EosInstance->SpecVol; + const auto &VertVisc = VertMixInstance->VertVisc; + const auto &LayerThickEdge = + AuxState->LayerThicknessAux.MeanLayerThickEdge; + + const I4 NVertLayers = VCoord->NVertLayers; + TeamPolicy Policy = + TriDiagDiffSolver::makeTeamPolicy(Mesh->NEdgesAll, NVertLayers); + + Kokkos::parallel_for( + Policy, KOKKOS_LAMBDA(const TeamMember &Team) { + const int IStart = Team.league_rank() * VecLength; + + TriDiagDiffScratch Scratch(Team, NVertLayers); + + Kokkos::parallel_for( + TeamThreadRange(Team, NVertLayers), [=](int K) { + for (int IVec = 0; IVec < VecLength; ++IVec) { + const int IEdge = IStart + IVec; + + if (IEdge >= LocNEdgesAll) { + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + Real G, H, X; + LocVelVertMixSetup(IEdge, K, DT, SpecVol, + LayerThickEdge, VertVisc, + NormalVelEdge, G, H, X); + Scratch.G(K, IVec) = G; + Scratch.H(K, IVec) = H; + Scratch.X(K, IVec) = X; + } + }); + + Team.team_barrier(); + TriDiagDiffSolver::solve(Team, Scratch); + Team.team_barrier(); + + Kokkos::parallel_for( + TeamThreadRange(Team, NVertLayers), [=](int K) { + for (int IVec = 0; IVec < VecLength; ++IVec) { + const int IEdge = IStart + IVec; + if (IEdge < LocNEdgesAll && + K >= MinLayerEdgeBot(IEdge) && + K <= MaxLayerEdgeTop(IEdge)) { + NormalVelEdge(IEdge, K) = Scratch.X(K, IVec); + } + } + }); + }); + } + Pacer::stop("Tend:velocityVertMix", 1); + } +} + +// Apply implicit tracer vertical mixing +void VertMix::applyTracerVertMixImplicit( + OceanState *State, ///< [in] State variables + const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + Array3DReal &TracerArray, ///< [in] Tracer array + int NTracers, ///< [in] Number of tracers + int ThickTimeLevel, ///< [in] Time level + int VelTimeLevel ///< [in] Time level +) { + + OMEGA_SCOPE(LocNCellsAll, Mesh->NCellsAll); + OMEGA_SCOPE(LocTracerVertMixSetup, TracerVertMixSetup); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + + const Array2DReal &LayerThickCell = State->LayerThickness[ThickTimeLevel]; + + if (LocTracerVertMixSetup.Enabled) { + Pacer::start("Tend:tracerVertMix", 1); + + Eos *EosInstance = Eos::getInstance(); + VertMix *VertMixInstance = VertMix::getInstance(); + + if (!EosInstance) { + LOG_WARN("Eos has not been initialized. Skipping calculation of " + "PresGradZ tendency"); + } else if (!VertMixInstance) { + LOG_WARN("VertMix has not been initialized. Skipping calculation of " + "VelVertMix tendency"); + } else { + + // Obtain TimeStep + const auto *DefTimeStepper = TimeStepper::getDefault(); + const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); + R8 DT; + TimeStep.get(DT, TimeUnits::Seconds); + + const auto &SpecVol = EosInstance->SpecVol; + const auto &VertDiff = VertMixInstance->VertDiff; + + const I4 NVertLayers = VCoord->NVertLayers; + TeamPolicy Policy = + TriDiagDiffSolver::makeTeamPolicy(Mesh->NCellsAll, NVertLayers); + + for (int LT = 0; LT < NTracers; ++LT) { + const I4 L = LT; + + Kokkos::parallel_for( + Policy, KOKKOS_LAMBDA(const TeamMember &Team) { + const int IStart = Team.league_rank() * VecLength; + + TriDiagDiffScratch Scratch(Team, NVertLayers); + + Kokkos::parallel_for( + TeamThreadRange(Team, NVertLayers), [=](int K) { + for (int IVec = 0; IVec < VecLength; ++IVec) { + const int ICell = IStart + IVec; + + if (ICell >= LocNCellsAll) { + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + Real G, H, X; + LocTracerVertMixSetup(L, ICell, K, DT, SpecVol, + LayerThickCell, VertDiff, + TracerArray, G, H, X); + Scratch.G(K, IVec) = G; + Scratch.H(K, IVec) = H; + Scratch.X(K, IVec) = X; + } + }); + + Team.team_barrier(); + TriDiagDiffSolver::solve(Team, Scratch); + Team.team_barrier(); + + Kokkos::parallel_for( + TeamThreadRange(Team, NVertLayers), [=](int K) { + for (int IVec = 0; IVec < VecLength; ++IVec) { + const int ICell = IStart + IVec; + if (ICell < LocNCellsAll && + K >= MinLayerCell(ICell) && + K <= MaxLayerCell(ICell)) { + TracerArray(L, ICell, K) = Scratch.X(K, IVec); + } + } + }); + }); + + } // for LT + } + Pacer::stop("Tend:tracerVertMix", 1); + } + +} // end all tendency compute + +/// Apply implicit vertical mixing to velocities and tracers +void VertMix::applyVertMixImplicit(OceanState *State, AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, + int TimeLevel) { + + // get NormalVelocity + Array2DReal NormalVelEdge = State->getNormalVelocity(TimeLevel); + + // get temperature and salinity + I4 ConservTempIdx; + I4 AbsSalinityIdx; + Tracers::getIndex(ConservTempIdx, "Temperature"); + Tracers::getIndex(AbsSalinityIdx, "Salinity"); + + const auto ConservTemp = + Kokkos::subview(TracerArray, ConservTempIdx, Kokkos::ALL, Kokkos::ALL); + const auto AbsSalinity = + Kokkos::subview(TracerArray, AbsSalinityIdx, Kokkos::ALL, Kokkos::ALL); + + // TODO: Temporary handling of computation of tangential velocity + // Compute tangential velocity + OMEGA_SCOPE(MinLayerEdgeTop, VCoord->MinLayerEdgeTop); + OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot); + OMEGA_SCOPE(LocTangentialVelocity, TangentialVelocity); + + TangentialReconOnEdge TanReconEdge(Mesh); + + parallelForOuter({Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeTop(IEdge); + const int KMax = MaxLayerEdgeBot(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + TanReconEdge(LocTangentialVelocity, IEdge, KChunk, + NormalVelEdge); + }); + }); + + + // Update Pressure, SpecVol + AuxState->computeMomVertAux(State, TracerArray, TimeLevel, + TimeLevel); + + // Compute Brunt-Vaisala frequency squared + Eos *EqState = Eos::getInstance(); + EqState->computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, + VCoord->PressureInterface, + EqState->SpecVol); + + // Compute vertical mixing coefficients + computeVertMix(NormalVelEdge, TangentialVelocity, + EqState->BruntVaisalaFreqSq); + + // Apply implicit mixing to velocities + applyVelVertMixImplicit(State, AuxState, TimeLevel, TimeLevel); + + // Apply implicit mixing to tracers + applyTracerVertMixImplicit(State, AuxState, TracerArray, NTracers, + TimeLevel, TimeLevel); +} + } // namespace OMEGA diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index 9466eb67d7ec..eb7879d7f015 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -228,6 +228,127 @@ class OneTwoOneFilter { Array1DI4 MaxLayerCell; }; +/// Velocity vertical mixing +class VelVertMixSetupOnEdge { + public: + bool Enabled; + Real LocRhoSw; + + /// constructor declaration + VelVertMixSetupOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); + + /// The functor takes edge index, vertical chunk index, and arrays for + /// layer specific volume, layer thickness on edge, + /// interface pressure, and outputs tendency array + KOKKOS_FUNCTION void operator()(I4 IEdge, I4 K, Real DT, + const Array2DReal &SpecVol, + const Array2DReal &LayerThickEdge, + const Array2DReal &VertVisc, + const Array2DReal &NormalVelEdge, Real &G, + Real &H, Real &X) const { + + const I4 JCell0 = CellsOnEdge(IEdge, 0); + const I4 JCell1 = CellsOnEdge(IEdge, 1); + + const I4 KMin = MinLayerEdgeBot(IEdge); + const I4 KMax = MaxLayerEdgeTop(IEdge); + + G = 0.0_Real; + H = 1.0_Real; + + if (K < KMin || K > KMax) { + X = 0.0_Real; + } else { + if (K < KMax) { + const Real LayerThickEdgeK = LayerThickEdge(IEdge, K); + const Real LayerThickEdgeKp1 = LayerThickEdge(IEdge, K + 1); + + const Real LayerThickEdgeTop = + 0.5_Real * (LayerThickEdgeK + LayerThickEdgeKp1); + + // Interpolation from cell center to top using + // the two-point linear interpolation + const Real SpecVolEdgeTop = + (0.5_Real * (SpecVol(JCell0, K) + SpecVol(JCell1, K)) * + LayerThickEdgeKp1 + + 0.5_Real * (SpecVol(JCell0, K + 1) + SpecVol(JCell1, K + 1)) * + LayerThickEdgeK) / + (LayerThickEdgeK + LayerThickEdgeKp1); + + const Real ViscAlphaEdgeTop = + 0.5_Real * (VertVisc(JCell0, K + 1) + VertVisc(JCell1, K + 1)) / + (LocRhoSw * SpecVolEdgeTop); + + G = DT * ViscAlphaEdgeTop / + (LayerThickEdgeTop * LayerThickEdge(IEdge, K)); + } + X = NormalVelEdge(IEdge, K); + } + } + + private: + I4 NVertLayers; + Array2DI4 CellsOnEdge; + Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; +}; + +// Tracer vertical mixing term +class TracerVertMixSetupOnCell { + public: + bool Enabled; + Real LocRhoSw; + + TracerVertMixSetupOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + KOKKOS_FUNCTION void operator()(I4 L, I4 ICell, I4 K, Real DT, + const Array2DReal &SpecVol, + const Array2DReal &LayerThickCell, + const Array2DReal &VertDiff, + const Array3DReal &TracersOnCell, Real &G, + Real &H, Real &X) const { + + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + + G = 0.0_Real; + H = 1.0_Real; + + if (K < KMin || K > KMax) { + X = 0.0_Real; + } else { + if (K < KMax) { + const Real LayerThickCellK = LayerThickCell(ICell, K); + const Real LayerThickCellKp1 = LayerThickCell(ICell, K + 1); + + const Real LayerThickCellTop = + 0.5_Real * (LayerThickCellK + LayerThickCellKp1); + + // Interpolation from cell center to top using + // the two-point linear interpolation + const Real SpecVolCellTop = + (SpecVol(ICell, K) * LayerThickCellKp1 + + SpecVol(ICell, K + 1) * LayerThickCellK) / + (LayerThickCellK + LayerThickCellKp1); + + const Real DiffAlphaCellTop = + VertDiff(ICell, K + 1) / (LocRhoSw * SpecVolCellTop); + + G = DT * DiffAlphaCellTop / + (LayerThickCellTop * LayerThickCell(ICell, K)); + } + X = TracersOnCell(L, ICell, K); + } + } + + private: + I4 NVertLayers; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + + /// Class for Vertical Mixing Coefficient (VertMix) calculations class VertMix { public: @@ -243,6 +364,9 @@ class VertMix { Array2DReal GradRichNumSmoothed; ///< Smoothed Gradient Richardson number field + //TODO: Temporary handling of TangentialVelocity + Array2DReal TangentialVelocity; ///< Tangential velocity + std::string VertDiffFldName; ///< Field name for vertical diffusivity std::string VertViscFldName; ///< Field name for vertical viscosity std::string @@ -264,6 +388,9 @@ class VertMix { ///< calculation OneTwoOneFilter ComputeOneTwoOneFilter; ///< Functor for 1-2-1 filtering + VelVertMixSetupOnEdge VelVertMixSetup; + TracerVertMixSetupOnCell TracerVertMixSetup; + /// Compute vertical diffusivity and viscosity for all cells/layers void computeVertMix(const Array2DReal &NormalVelocity, const Array2DReal &TangentialVelocity, @@ -272,6 +399,19 @@ class VertMix { /// Initialize VertMix from config and mesh static void init(); + void applyVelVertMixImplicit(OceanState *State, + const AuxiliaryState *AuxState, + int ThickTimeLevel, int VelTimeLevel); + void applyTracerVertMixImplicit(OceanState *State, + const AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, + int ThickTimeLevel, int VelTimeLevel); + + /// Apply implicit vertical mixing to velocities and tracers + void applyVertMixImplicit(OceanState *State, AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, + int TimeLevel); + private: /// Private constructor VertMix(const std::string &Name, const HorzMesh *Mesh, diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index 11ddc664d1eb..e20d534152ff 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -6,6 +6,7 @@ #include "ForwardBackwardStepper.h" #include "Pacer.h" +#include "VertMix.h" namespace OMEGA { @@ -37,6 +38,8 @@ void ForwardBackwardStepper::doStep( const int ThickNextLevel = 1; const int TracerNextLevel = 1; + int NTracers = Tracers::getNumTracers(); + Array3DReal CurTracerArray = Tracers::getAll(TracerCurLevel); Array3DReal NextTracerArray = Tracers::getAll(TracerNextLevel); @@ -86,6 +89,13 @@ void ForwardBackwardStepper::doStep( Tracers::updateTimeLevels(); Pacer::stop("ForwardBackward:haloExch", 3); + // Apply implicit vertical mixing + if (VMix->VelVertMixSetup.Enabled or + VMix->TracerVertMixSetup.Enabled) { + VMix->applyVertMixImplicit(State, AuxState, NextTracerArray, + NTracers, State->CurTimeIndex); + } + validateOceanState(State, AuxState, VertCoord::getDefault(), 0); // Advance the clock and update the simulation time diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index 372415859102..3d43c744a8fa 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -6,6 +6,7 @@ #include "RungeKutta2Stepper.h" #include "Pacer.h" +#include "VertMix.h" namespace OMEGA { @@ -31,6 +32,8 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state const int CurLevel = 0; const int NextLevel = 1; + int NTracers = Tracers::getNumTracers(); + Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); @@ -70,6 +73,13 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state Tracers::updateTimeLevels(); Pacer::stop("RK2:haloExch", 3); + // Apply implicit vertical mixing + if (VMix->VelVertMixSetup.Enabled or + VMix->TracerVertMixSetup.Enabled) { + VMix->applyVertMixImplicit(State, AuxState, NextTracerArray, + NTracers, State->CurTimeIndex); + } + validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); // Advance the clock and update the simulation time diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index 5ad5386e28f9..925164b5c3ce 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -6,6 +6,7 @@ #include "RungeKutta4Stepper.h" #include "Pacer.h" +#include "VertMix.h" namespace OMEGA { @@ -78,11 +79,14 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state const int CurLevel = 0; const int NextLevel = 1; + int NTracers = Tracers::getNumTracers(); Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); TimeInstant ForcingStageTime = SimTime; + VertMix *VMix = VertMix::getInstance(); + for (int Stage = 0; Stage < NStages; ++Stage) { const TimeInstant StageTime = SimTime + RKC[Stage] * TimeStep; // first stage does: @@ -135,6 +139,13 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state Tracers::updateTimeLevels(); Pacer::stop("RK4:haloExch", 3); + // Apply implicit vertical mixing + if (VMix->VelVertMixSetup.Enabled or + VMix->TracerVertMixSetup.Enabled) { + VMix->applyVertMixImplicit(State, AuxState, NextTracerArray, + NTracers, State->CurTimeIndex); + } + validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); // Advance the clock and update the simulation time diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index 1678a08470f6..9fc87452298d 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -31,6 +31,7 @@ #include "TimeStepper.h" #include "Tracers.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include #include @@ -108,6 +109,9 @@ void initIOStreamTest(Clock *&ModelClock // Model clock PressureGrad::init(); + // Initialize VertMix + VertMix::init(); + // Intialize Tendencies Tendencies::init(); @@ -229,6 +233,7 @@ int main(int argc, char **argv) { // Clean up environments TimeStepper::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Tendencies::clear(); Tracers::clear(); OceanState::clear(); diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp index 252056ef2074..fa3fa8f5fb1c 100644 --- a/components/omega/test/ocn/StateTest.cpp +++ b/components/omega/test/ocn/StateTest.cpp @@ -28,6 +28,7 @@ #include "Tendencies.h" #include "TimeStepper.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include @@ -95,6 +96,9 @@ void initStateTest() { // Initialize pressure gradient PressureGrad::init(); + // Create tendencies + VertMix::init(); + // Create tendencies Tendencies::init(); @@ -408,6 +412,7 @@ int main(int argc, char *argv[]) { Tracers::clear(); AuxiliaryState::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Tendencies::clear(); TimeStepper::clear(); HorzMesh::clear(); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 5c8bb8d3d3f0..5268d0a436ea 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -22,6 +22,7 @@ #include "Pacer.h" #include "TimeStepper.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include @@ -142,6 +143,7 @@ int initTendenciesTest(const std::string &mesh) { PressureGrad::init(); Eos::init(); Forcing::init(); + VertMix::init(); int StateErr = OceanState::init(); if (StateErr != 0) { @@ -177,6 +179,7 @@ int testTendencies() { const auto VAdv = VertAdv::getDefault(); const auto PGrad = PressureGrad::getDefault(); const auto EqState = Eos::getInstance(); + const auto VMix = VertMix::getInstance(); VCoord->NVertLayers = 12; // test creation of another tendencies @@ -188,7 +191,7 @@ int testTendencies() { int NTracersTest = 3; Tendencies::create("TestTendencies", Mesh, VCoord, VAdv, PGrad, EqState, - NTracersTest, ZeroTimeStep, &TendConfig); + VMix, NTracersTest, ZeroTimeStep, &TendConfig); // test retrievel of another tendencies if (Tendencies::get("TestTendencies")) { @@ -340,6 +343,7 @@ void finalizeTendenciesTest() { Forcing::clear(); Tracers::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Eos::destroyInstance(); AuxiliaryState::clear(); OceanState::clear(); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 81f6bcef5938..ec458273c351 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -37,6 +37,7 @@ #include "TimeMgr.h" #include "Tracers.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include @@ -187,6 +188,7 @@ int initTimeStepperTest(const std::string &mesh) { AuxiliaryState::init(); Eos::init(); PressureGrad::init(); + VertMix::init(); Tendencies::init(); // finish initializing default time stepper @@ -205,6 +207,7 @@ int initTimeStepperTest(const std::string &mesh) { auto *DefHalo = Halo::getDefault(); auto *DefEos = Eos::getInstance(); auto *DefPGrad = PressureGrad::getDefault(); + auto *DefVMix = VertMix::getInstance(); int NTracers = Tracers::getNumTracers(); const int NTimeLevels = 2; @@ -233,8 +236,8 @@ int initTimeStepperTest(const std::string &mesh) { // Creating non-default tendencies with custom velocity tendencies auto *TestTendencies = Tendencies::create( "TestTendencies", DefMesh, DefVertCoord, DefVAdv, DefPGrad, DefEos, - NTracers, ZeroTimeStep, &Options, Tendencies::CustomTendencyType{}, - DecayVelocityTendency{}); + DefVMix, NTracers, ZeroTimeStep, &Options, + Tendencies::CustomTendencyType{}, DecayVelocityTendency{}); if (!TestTendencies) { Err++; LOG_ERROR("TimeStepperTest: error creating test tendencies"); @@ -256,6 +259,8 @@ int initTimeStepperTest(const std::string &mesh) { DefVAdv->ThickVertAdvEnabled = false; DefVAdv->VelVertAdvEnabled = false; DefVAdv->TracerVertAdvEnabled = false; + DefVMix->VelVertMixSetup.Enabled = false; + DefVMix->TracerVertMixSetup.Enabled = false; return Err; } @@ -295,6 +300,7 @@ void finalizeTimeStepperTest() { Tracers::clear(); TimeStepper::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Eos::destroyInstance(); Tendencies::clear(); AuxiliaryState::clear(); From 2589d47806a6e294e71f16c769b831022580e92f Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Sun, 17 May 2026 22:54:20 -0400 Subject: [PATCH 05/19] Update VertMixImplicit code and linting --- components/omega/src/ocn/OceanInit.cpp | 49 ------------------- components/omega/src/ocn/Tendencies.cpp | 17 +++---- components/omega/src/ocn/VertMix.cpp | 48 +++++++++--------- components/omega/src/ocn/VertMix.h | 20 ++++---- .../timeStepping/ForwardBackwardStepper.cpp | 9 ++-- .../src/timeStepping/RungeKutta2Stepper.cpp | 9 ++-- .../src/timeStepping/RungeKutta4Stepper.cpp | 9 ++-- 7 files changed, 54 insertions(+), 107 deletions(-) diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index 1348f197c747..e7fe9938401e 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -203,55 +203,6 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator } Tracers::copyToHost(CurTimeLevel); -//////////////////////////////////////////////// -//////////////////////////////////////////////// - - // For initial GeometricThick -> PseudoThick - HorzMesh *Mesh = HorzMesh::getDefault(); - VertCoord *VCoord = VertCoord::getDefault(); - Eos *EosInstance = Eos::getInstance(); - - // get temperature and salinity - I4 ConservTempIdx; - I4 AbsSalinityIdx; - Array3DReal TracerArray = Tracers::getAll(0); - - Tracers::getIndex(ConservTempIdx, "Temperature"); - Tracers::getIndex(AbsSalinityIdx, "Salinity"); - - const auto ConservTemp = - Kokkos::subview(TracerArray, ConservTempIdx, Kokkos::ALL, Kokkos::ALL); - const auto AbsSalinity = - Kokkos::subview(TracerArray, AbsSalinityIdx, Kokkos::ALL, Kokkos::ALL); - - Array2DReal PressureMidDbar("PressureMidDbar", VCoord->PressureMid.layout()); - deepCopy(PressureMidDbar, 0.0); - EosInstance->computeSpecVol(ConservTemp, AbsSalinity, PressureMidDbar); - const auto &SpecVol = EosInstance->SpecVol; - - // get layer thickness - Array2DReal LayerThickCell = DefState->getLayerThickness(0); - OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell); - OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell); - OMEGA_SCOPE(LocNCellsAll, Mesh->NCellsAll); - - parallelForOuter( - "GtoP", {LocNCellsAll}, - KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { - const int KMin = LocMinLayerCell(ICell); - const int KMax = LocMaxLayerCell(ICell); - const int KRange = vertRange(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - const int K = KMin + KChunk; - LayerThickCell(ICell,K) = LayerThickCell(ICell,K) / (SpecVol(ICell,K)*1026.0); - - }); - }); - -//////////////////////////////////////////////// -////////////////////////////////////////////////// - return Err; } // end ocnInit diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 69c4fa22b27d..c0feea0d8652 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -40,7 +40,7 @@ void Tendencies::init() { TimeStepper *DefTimeStepper = TimeStepper::getDefault(); Eos *DefEos = Eos::getInstance(); PressureGrad *DefPGrad = PressureGrad::getDefault(); - VertMix *DefVertMix = VertMix::getInstance(); + VertMix *DefVertMix = VertMix::getInstance(); I4 NTracers = Tracers::getNumTracers(); @@ -80,10 +80,10 @@ void Tendencies::init() { TimeInterval TimeStep = DefTimeStepper->getTimeStep(); // Ceate default tendencies - Tendencies::DefaultTendencies = create( - "Default", DefHorzMesh, DefVertCoord, DefVertAdv, DefPGrad, DefEos, - DefVertMix, NTracers, TimeStep, &TendConfig, CustomThickTend, - CustomVelTend); + Tendencies::DefaultTendencies = + create("Default", DefHorzMesh, DefVertCoord, DefVertAdv, DefPGrad, + DefEos, DefVertMix, NTracers, TimeStep, &TendConfig, + CustomThickTend, CustomVelTend); DefaultTendencies->readConfig(OmegaConfig); @@ -330,7 +330,6 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config this->VMix->TracerVertMixSetup.Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: TracerVertMixTendencyEnable not found in TendConfig"); - } //------------------------------------------------------------------------------ @@ -451,9 +450,9 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStepIn, ///< [in] Time step Config *Options) ///< [in] Configuration options - : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, VertMix::getInstance(), - NTracersIn, TimeStepIn, Options, CustomTendencyType{}, - CustomTendencyType{}) {} + : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, + VertMix::getInstance(), NTracersIn, TimeStepIn, Options, + CustomTendencyType{}, CustomTendencyType{}) {} //------------------------------------------------------------------------------ // Compute tendencies for the pseudo-thickness equation diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index c6e6749d3866..8cb9fe776c82 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -16,8 +16,8 @@ #include "GlobalConstants.h" #include "HorzMesh.h" #include "HorzOperators.h" -#include "TriDiagSolvers.h" #include "TimeStepper.h" +#include "TriDiagSolvers.h" namespace OMEGA { @@ -51,8 +51,7 @@ VelVertMixSetupOnEdge::VelVertMixSetupOnEdge(const HorzMesh *Mesh, TracerVertMixSetupOnCell::TracerVertMixSetupOnCell(const HorzMesh *Mesh, const VertCoord *VCoord) : Enabled(false), LocRhoSw(RhoSw), NVertLayers(VCoord->NVertLayers), - MinLayerCell(VCoord->MinLayerCell), - MaxLayerCell(VCoord->MaxLayerCell) {} + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} /// Constructor for VertMix VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object @@ -70,9 +69,9 @@ VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object GradRichNumSmoothed = Array2DReal("GradRichNumSmoothed", Mesh->NCellsSize, VCoord->NVertLayersP1); - //TODO: Temporary handling of TangentialVelocity - TangentialVelocity = Array2DReal("TangentialVelocity", Mesh->NEdgesSize, - VCoord->NVertLayers); + // TODO: Temporary handling of TangentialVelocity + TangentialVelocity = + Array2DReal("TangentialVelocity", Mesh->NEdgesSize, VCoord->NVertLayers); defineFields(); } @@ -521,7 +520,6 @@ void VertMix::defineFields() { } // end defineIOFields - // Dpply implicit velocity vertical mixing void VertMix::applyVelVertMixImplicit( OceanState *State, ///< [in] State variables @@ -535,7 +533,8 @@ void VertMix::applyVelVertMixImplicit( OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); - const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; + const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; + const Array2DReal &LayerThickCell = State->LayerThickness[ThickTimeLevel]; // Compute velocity vertical mixing if (LocVelVertMixSetup.Enabled) { @@ -560,8 +559,6 @@ void VertMix::applyVelVertMixImplicit( const auto &SpecVol = EosInstance->SpecVol; const auto &VertVisc = VertMixInstance->VertVisc; - const auto &LayerThickEdge = - AuxState->LayerThicknessAux.MeanLayerThickEdge; const I4 NVertLayers = VCoord->NVertLayers; TeamPolicy Policy = @@ -587,7 +584,7 @@ void VertMix::applyVelVertMixImplicit( Real G, H, X; LocVelVertMixSetup(IEdge, K, DT, SpecVol, - LayerThickEdge, VertVisc, + LayerThickCell, VertVisc, NormalVelEdge, G, H, X); Scratch.G(K, IVec) = G; Scratch.H(K, IVec) = H; @@ -716,9 +713,9 @@ void VertMix::applyTracerVertMixImplicit( } // end all tendency compute /// Apply implicit vertical mixing to velocities and tracers -void VertMix::applyVertMixImplicit(OceanState *State, AuxiliaryState *AuxState, - Array3DReal &TracerArray, int NTracers, - int TimeLevel) { +void VertMix::VertMixImplicit(OceanState *State, AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, + int TimeLevel) { // get NormalVelocity Array2DReal NormalVelEdge = State->getNormalVelocity(TimeLevel); @@ -734,6 +731,9 @@ void VertMix::applyVertMixImplicit(OceanState *State, AuxiliaryState *AuxState, const auto AbsSalinity = Kokkos::subview(TracerArray, AbsSalinityIdx, Kokkos::ALL, Kokkos::ALL); + // get an instance of equation of state + Eos *EqState = Eos::getInstance(); + // TODO: Temporary handling of computation of tangential velocity // Compute tangential velocity OMEGA_SCOPE(MinLayerEdgeTop, VCoord->MinLayerEdgeTop); @@ -742,8 +742,8 @@ void VertMix::applyVertMixImplicit(OceanState *State, AuxiliaryState *AuxState, TangentialReconOnEdge TanReconEdge(Mesh); - parallelForOuter({Mesh->NEdgesAll}, - KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { const int KMin = MinLayerEdgeTop(IEdge); const int KMax = MaxLayerEdgeBot(IEdge); const int KRange = vertRangeChunked(KMin, KMax); @@ -754,27 +754,23 @@ void VertMix::applyVertMixImplicit(OceanState *State, AuxiliaryState *AuxState, }); }); - // Update Pressure, SpecVol - AuxState->computeMomVertAux(State, TracerArray, TimeLevel, - TimeLevel); + AuxState->computeMomVertAux(State, TracerArray, TimeLevel, TimeLevel); // Compute Brunt-Vaisala frequency squared - Eos *EqState = Eos::getInstance(); - EqState->computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, - VCoord->PressureInterface, - EqState->SpecVol); + EqState->computeBruntVaisalaFreqSq( + ConservTemp, AbsSalinity, VCoord->PressureInterface, EqState->SpecVol); // Compute vertical mixing coefficients - computeVertMix(NormalVelEdge, TangentialVelocity, + computeVertMix(NormalVelEdge, LocTangentialVelocity, EqState->BruntVaisalaFreqSq); // Apply implicit mixing to velocities applyVelVertMixImplicit(State, AuxState, TimeLevel, TimeLevel); // Apply implicit mixing to tracers - applyTracerVertMixImplicit(State, AuxState, TracerArray, NTracers, - TimeLevel, TimeLevel); + applyTracerVertMixImplicit(State, AuxState, TracerArray, NTracers, TimeLevel, + TimeLevel); } } // namespace OMEGA diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index eb7879d7f015..903262b95567 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -242,7 +242,7 @@ class VelVertMixSetupOnEdge { /// interface pressure, and outputs tendency array KOKKOS_FUNCTION void operator()(I4 IEdge, I4 K, Real DT, const Array2DReal &SpecVol, - const Array2DReal &LayerThickEdge, + const Array2DReal &LayerThickCell, const Array2DReal &VertVisc, const Array2DReal &NormalVelEdge, Real &G, Real &H, Real &X) const { @@ -260,8 +260,11 @@ class VelVertMixSetupOnEdge { X = 0.0_Real; } else { if (K < KMax) { - const Real LayerThickEdgeK = LayerThickEdge(IEdge, K); - const Real LayerThickEdgeKp1 = LayerThickEdge(IEdge, K + 1); + const Real LayerThickEdgeK = 0.5_Real * (LayerThickCell(JCell0, K) + + LayerThickCell(JCell1, K)); + const Real LayerThickEdgeKp1 = + 0.5_Real * + (LayerThickCell(JCell0, K + 1) + LayerThickCell(JCell1, K + 1)); const Real LayerThickEdgeTop = 0.5_Real * (LayerThickEdgeK + LayerThickEdgeKp1); @@ -279,8 +282,7 @@ class VelVertMixSetupOnEdge { 0.5_Real * (VertVisc(JCell0, K + 1) + VertVisc(JCell1, K + 1)) / (LocRhoSw * SpecVolEdgeTop); - G = DT * ViscAlphaEdgeTop / - (LayerThickEdgeTop * LayerThickEdge(IEdge, K)); + G = DT * ViscAlphaEdgeTop / (LayerThickEdgeTop * LayerThickEdgeK); } X = NormalVelEdge(IEdge, K); } @@ -348,7 +350,6 @@ class TracerVertMixSetupOnCell { Array1DI4 MaxLayerCell; }; - /// Class for Vertical Mixing Coefficient (VertMix) calculations class VertMix { public: @@ -364,7 +365,7 @@ class VertMix { Array2DReal GradRichNumSmoothed; ///< Smoothed Gradient Richardson number field - //TODO: Temporary handling of TangentialVelocity + // TODO: Temporary handling of TangentialVelocity Array2DReal TangentialVelocity; ///< Tangential velocity std::string VertDiffFldName; ///< Field name for vertical diffusivity @@ -408,9 +409,8 @@ class VertMix { int ThickTimeLevel, int VelTimeLevel); /// Apply implicit vertical mixing to velocities and tracers - void applyVertMixImplicit(OceanState *State, AuxiliaryState *AuxState, - Array3DReal &TracerArray, int NTracers, - int TimeLevel); + void VertMixImplicit(OceanState *State, AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, int TimeLevel); private: /// Private constructor diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index e20d534152ff..e1b46d79f2fc 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -43,6 +43,8 @@ void ForwardBackwardStepper::doStep( Array3DReal CurTracerArray = Tracers::getAll(TracerCurLevel); Array3DReal NextTracerArray = Tracers::getAll(TracerNextLevel); + VertMix *VMix = VertMix::getInstance(); + if (State == nullptr) LOG_CRITICAL("Invalid State"); if (AuxState == nullptr) @@ -90,10 +92,9 @@ void ForwardBackwardStepper::doStep( Pacer::stop("ForwardBackward:haloExch", 3); // Apply implicit vertical mixing - if (VMix->VelVertMixSetup.Enabled or - VMix->TracerVertMixSetup.Enabled) { - VMix->applyVertMixImplicit(State, AuxState, NextTracerArray, - NTracers, State->CurTimeIndex); + if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { + VMix->VertMixImplicit(State, AuxState, NextTracerArray, NTracers, + VelNextLevel); } validateOceanState(State, AuxState, VertCoord::getDefault(), 0); diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index 3d43c744a8fa..387e24f79185 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -37,6 +37,8 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); + VertMix *VMix = VertMix::getInstance(); + prescribeState(State, CurLevel, State, CurLevel, SimTime); // q = (h,u,phi) @@ -74,10 +76,9 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state Pacer::stop("RK2:haloExch", 3); // Apply implicit vertical mixing - if (VMix->VelVertMixSetup.Enabled or - VMix->TracerVertMixSetup.Enabled) { - VMix->applyVertMixImplicit(State, AuxState, NextTracerArray, - NTracers, State->CurTimeIndex); + if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { + VMix->VertMixImplicit(State, AuxState, NextTracerArray, NTracers, + NextLevel); } validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index 925164b5c3ce..2b6ddc96c1e4 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -79,7 +79,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state const int CurLevel = 0; const int NextLevel = 1; - int NTracers = Tracers::getNumTracers(); + int NTracers = Tracers::getNumTracers(); Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); @@ -140,10 +140,9 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state Pacer::stop("RK4:haloExch", 3); // Apply implicit vertical mixing - if (VMix->VelVertMixSetup.Enabled or - VMix->TracerVertMixSetup.Enabled) { - VMix->applyVertMixImplicit(State, AuxState, NextTracerArray, - NTracers, State->CurTimeIndex); + if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { + VMix->VertMixImplicit(State, AuxState, NextTracerArray, NTracers, + NextLevel); } validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); From d730ab465f02a6b49687e3a92cbc3e96754e6f96 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 18 May 2026 04:56:32 -0400 Subject: [PATCH 06/19] Update codes and add comments --- components/omega/src/ocn/VertMix.cpp | 37 ++++++++++++------- .../timeStepping/ForwardBackwardStepper.cpp | 5 ++- .../src/timeStepping/RungeKutta2Stepper.cpp | 5 ++- .../src/timeStepping/RungeKutta4Stepper.cpp | 5 ++- 4 files changed, 32 insertions(+), 20 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 8cb9fe776c82..2b8b36bf745d 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -520,7 +520,7 @@ void VertMix::defineFields() { } // end defineIOFields -// Dpply implicit velocity vertical mixing +// Apply implicit velocity vertical mixing void VertMix::applyVelVertMixImplicit( OceanState *State, ///< [in] State variables const AuxiliaryState *AuxState, ///< [in] Auxilary state variables @@ -567,9 +567,12 @@ void VertMix::applyVelVertMixImplicit( Kokkos::parallel_for( Policy, KOKKOS_LAMBDA(const TeamMember &Team) { const int IStart = Team.league_rank() * VecLength; + const int ILen = + Kokkos::max(0, Kokkos::min(VecLength, LocNEdgesAll - IStart)); TriDiagDiffScratch Scratch(Team, NVertLayers); + // Construct a tri-diag diffusion matrix and RHS Kokkos::parallel_for( TeamThreadRange(Team, NVertLayers), [=](int K) { for (int IVec = 0; IVec < VecLength; ++IVec) { @@ -583,31 +586,34 @@ void VertMix::applyVelVertMixImplicit( } Real G, H, X; - LocVelVertMixSetup(IEdge, K, DT, SpecVol, - LayerThickCell, VertVisc, - NormalVelEdge, G, H, X); + LocVelVertMixSetup(IEdge, K, DT, SpecVol, LayerThickCell, + VertVisc, NormalVelEdge, G, H, X); + Scratch.G(K, IVec) = G; Scratch.H(K, IVec) = H; Scratch.X(K, IVec) = X; } }); + // Solve the tri-diag diffusion system Team.team_barrier(); TriDiagDiffSolver::solve(Team, Scratch); Team.team_barrier(); + // Store the solution vector X Kokkos::parallel_for( TeamThreadRange(Team, NVertLayers), [=](int K) { - for (int IVec = 0; IVec < VecLength; ++IVec) { + for (int IVec = 0; IVec < ILen; ++IVec) { const int IEdge = IStart + IVec; - if (IEdge < LocNEdgesAll && - K >= MinLayerEdgeBot(IEdge) && + + if (K >= MinLayerEdgeBot(IEdge) && K <= MaxLayerEdgeTop(IEdge)) { NormalVelEdge(IEdge, K) = Scratch.X(K, IVec); } } }); }); + } Pacer::stop("Tend:velocityVertMix", 1); } @@ -657,15 +663,16 @@ void VertMix::applyTracerVertMixImplicit( TeamPolicy Policy = TriDiagDiffSolver::makeTeamPolicy(Mesh->NCellsAll, NVertLayers); - for (int LT = 0; LT < NTracers; ++LT) { - const I4 L = LT; - + for (int L = 0; L < NTracers; ++L) { Kokkos::parallel_for( Policy, KOKKOS_LAMBDA(const TeamMember &Team) { const int IStart = Team.league_rank() * VecLength; + const int ILen = + Kokkos::max(0, Kokkos::min(VecLength, LocNCellsAll - IStart)); TriDiagDiffScratch Scratch(Team, NVertLayers); + // Construct a tri-diag diffusion matrix and RHS Kokkos::parallel_for( TeamThreadRange(Team, NVertLayers), [=](int K) { for (int IVec = 0; IVec < VecLength; ++IVec) { @@ -688,16 +695,18 @@ void VertMix::applyTracerVertMixImplicit( } }); + // Solve the tri-diag diffusion system Team.team_barrier(); TriDiagDiffSolver::solve(Team, Scratch); Team.team_barrier(); + // Store the solution vector X Kokkos::parallel_for( TeamThreadRange(Team, NVertLayers), [=](int K) { - for (int IVec = 0; IVec < VecLength; ++IVec) { + for (int IVec = 0; IVec < ILen; ++IVec) { const int ICell = IStart + IVec; - if (ICell < LocNCellsAll && - K >= MinLayerCell(ICell) && + + if (K >= MinLayerCell(ICell) && K <= MaxLayerCell(ICell)) { TracerArray(L, ICell, K) = Scratch.X(K, IVec); } @@ -705,7 +714,7 @@ void VertMix::applyTracerVertMixImplicit( }); }); - } // for LT + } // for L } Pacer::stop("Tend:tracerVertMix", 1); } diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index e1b46d79f2fc..bcb0d9463835 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -92,9 +92,10 @@ void ForwardBackwardStepper::doStep( Pacer::stop("ForwardBackward:haloExch", 3); // Apply implicit vertical mixing + CurTracerArray = Tracers::getAll(VelCurLevel); if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { - VMix->VertMixImplicit(State, AuxState, NextTracerArray, NTracers, - VelNextLevel); + VMix->VertMixImplicit(State, AuxState, CurTracerArray, NTracers, + State->CurTimeIndex); } validateOceanState(State, AuxState, VertCoord::getDefault(), 0); diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index 387e24f79185..16fc8fa36c0d 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -76,9 +76,10 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state Pacer::stop("RK2:haloExch", 3); // Apply implicit vertical mixing + CurTracerArray = Tracers::getAll(CurLevel); if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { - VMix->VertMixImplicit(State, AuxState, NextTracerArray, NTracers, - NextLevel); + VMix->VertMixImplicit(State, AuxState, CurTracerArray, NTracers, + State->CurTimeIndex); } validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index 2b6ddc96c1e4..3b21f972854e 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -140,9 +140,10 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state Pacer::stop("RK4:haloExch", 3); // Apply implicit vertical mixing + CurTracerArray = Tracers::getAll(CurLevel); if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { - VMix->VertMixImplicit(State, AuxState, NextTracerArray, NTracers, - NextLevel); + VMix->VertMixImplicit(State, AuxState, CurTracerArray, NTracers, + State->CurTimeIndex); } validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); From 5e932a668913d25d13e03867be2c147e0879c995 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 18 May 2026 09:37:25 -0400 Subject: [PATCH 07/19] Fix linting issues --- components/omega/src/ocn/VertMix.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 2b8b36bf745d..128a4aa91696 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -567,8 +567,8 @@ void VertMix::applyVelVertMixImplicit( Kokkos::parallel_for( Policy, KOKKOS_LAMBDA(const TeamMember &Team) { const int IStart = Team.league_rank() * VecLength; - const int ILen = - Kokkos::max(0, Kokkos::min(VecLength, LocNEdgesAll - IStart)); + const int ILen = Kokkos::max( + 0, Kokkos::min(VecLength, LocNEdgesAll - IStart)); TriDiagDiffScratch Scratch(Team, NVertLayers); @@ -586,8 +586,9 @@ void VertMix::applyVelVertMixImplicit( } Real G, H, X; - LocVelVertMixSetup(IEdge, K, DT, SpecVol, LayerThickCell, - VertVisc, NormalVelEdge, G, H, X); + LocVelVertMixSetup(IEdge, K, DT, SpecVol, + LayerThickCell, VertVisc, + NormalVelEdge, G, H, X); Scratch.G(K, IVec) = G; Scratch.H(K, IVec) = H; @@ -613,7 +614,6 @@ void VertMix::applyVelVertMixImplicit( } }); }); - } Pacer::stop("Tend:velocityVertMix", 1); } @@ -667,8 +667,8 @@ void VertMix::applyTracerVertMixImplicit( Kokkos::parallel_for( Policy, KOKKOS_LAMBDA(const TeamMember &Team) { const int IStart = Team.league_rank() * VecLength; - const int ILen = - Kokkos::max(0, Kokkos::min(VecLength, LocNCellsAll - IStart)); + const int ILen = Kokkos::max( + 0, Kokkos::min(VecLength, LocNCellsAll - IStart)); TriDiagDiffScratch Scratch(Team, NVertLayers); From cd32dc96b76589ee8ce54915da588e9b3a365bf6 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 27 May 2026 00:54:52 -0400 Subject: [PATCH 08/19] Change LayerThick to PseudoThick in VertMix tend --- components/omega/src/ocn/VertMix.cpp | 10 +++---- components/omega/src/ocn/VertMix.h | 41 ++++++++++++++-------------- 2 files changed, 26 insertions(+), 25 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 128a4aa91696..b186f01ae0b6 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -533,8 +533,8 @@ void VertMix::applyVelVertMixImplicit( OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); - const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; - const Array2DReal &LayerThickCell = State->LayerThickness[ThickTimeLevel]; + const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; + const Array2DReal &PseudoThickCell = State->PseudoThickness[ThickTimeLevel]; // Compute velocity vertical mixing if (LocVelVertMixSetup.Enabled) { @@ -587,7 +587,7 @@ void VertMix::applyVelVertMixImplicit( Real G, H, X; LocVelVertMixSetup(IEdge, K, DT, SpecVol, - LayerThickCell, VertVisc, + PseudoThickCell, VertVisc, NormalVelEdge, G, H, X); Scratch.G(K, IVec) = G; @@ -634,7 +634,7 @@ void VertMix::applyTracerVertMixImplicit( OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - const Array2DReal &LayerThickCell = State->LayerThickness[ThickTimeLevel]; + const Array2DReal &PseudoThickCell = State->PseudoThickness[ThickTimeLevel]; if (LocTracerVertMixSetup.Enabled) { Pacer::start("Tend:tracerVertMix", 1); @@ -687,7 +687,7 @@ void VertMix::applyTracerVertMixImplicit( Real G, H, X; LocTracerVertMixSetup(L, ICell, K, DT, SpecVol, - LayerThickCell, VertDiff, + PseudoThickCell, VertDiff, TracerArray, G, H, X); Scratch.G(K, IVec) = G; Scratch.H(K, IVec) = H; diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index 903262b95567..d9c0dc187e67 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -242,7 +242,7 @@ class VelVertMixSetupOnEdge { /// interface pressure, and outputs tendency array KOKKOS_FUNCTION void operator()(I4 IEdge, I4 K, Real DT, const Array2DReal &SpecVol, - const Array2DReal &LayerThickCell, + const Array2DReal &PseudoThickCell, const Array2DReal &VertVisc, const Array2DReal &NormalVelEdge, Real &G, Real &H, Real &X) const { @@ -260,29 +260,30 @@ class VelVertMixSetupOnEdge { X = 0.0_Real; } else { if (K < KMax) { - const Real LayerThickEdgeK = 0.5_Real * (LayerThickCell(JCell0, K) + - LayerThickCell(JCell1, K)); - const Real LayerThickEdgeKp1 = + const Real PseudoThickEdgeK = 0.5_Real * - (LayerThickCell(JCell0, K + 1) + LayerThickCell(JCell1, K + 1)); + (PseudoThickCell(JCell0, K) + PseudoThickCell(JCell1, K)); + const Real PseudoThickEdgeKp1 = + 0.5_Real * (PseudoThickCell(JCell0, K + 1) + + PseudoThickCell(JCell1, K + 1)); - const Real LayerThickEdgeTop = - 0.5_Real * (LayerThickEdgeK + LayerThickEdgeKp1); + const Real PseudoThickEdgeTop = + 0.5_Real * (PseudoThickEdgeK + PseudoThickEdgeKp1); // Interpolation from cell center to top using // the two-point linear interpolation const Real SpecVolEdgeTop = (0.5_Real * (SpecVol(JCell0, K) + SpecVol(JCell1, K)) * - LayerThickEdgeKp1 + + PseudoThickEdgeKp1 + 0.5_Real * (SpecVol(JCell0, K + 1) + SpecVol(JCell1, K + 1)) * - LayerThickEdgeK) / - (LayerThickEdgeK + LayerThickEdgeKp1); + PseudoThickEdgeK) / + (PseudoThickEdgeK + PseudoThickEdgeKp1); const Real ViscAlphaEdgeTop = 0.5_Real * (VertVisc(JCell0, K + 1) + VertVisc(JCell1, K + 1)) / (LocRhoSw * SpecVolEdgeTop); - G = DT * ViscAlphaEdgeTop / (LayerThickEdgeTop * LayerThickEdgeK); + G = DT * ViscAlphaEdgeTop / (PseudoThickEdgeTop * PseudoThickEdgeK); } X = NormalVelEdge(IEdge, K); } @@ -306,7 +307,7 @@ class TracerVertMixSetupOnCell { KOKKOS_FUNCTION void operator()(I4 L, I4 ICell, I4 K, Real DT, const Array2DReal &SpecVol, - const Array2DReal &LayerThickCell, + const Array2DReal &PseudoThickCell, const Array2DReal &VertDiff, const Array3DReal &TracersOnCell, Real &G, Real &H, Real &X) const { @@ -321,24 +322,24 @@ class TracerVertMixSetupOnCell { X = 0.0_Real; } else { if (K < KMax) { - const Real LayerThickCellK = LayerThickCell(ICell, K); - const Real LayerThickCellKp1 = LayerThickCell(ICell, K + 1); + const Real PseudoThickCellK = PseudoThickCell(ICell, K); + const Real PseudoThickCellKp1 = PseudoThickCell(ICell, K + 1); - const Real LayerThickCellTop = - 0.5_Real * (LayerThickCellK + LayerThickCellKp1); + const Real PseudoThickCellTop = + 0.5_Real * (PseudoThickCellK + PseudoThickCellKp1); // Interpolation from cell center to top using // the two-point linear interpolation const Real SpecVolCellTop = - (SpecVol(ICell, K) * LayerThickCellKp1 + - SpecVol(ICell, K + 1) * LayerThickCellK) / - (LayerThickCellK + LayerThickCellKp1); + (SpecVol(ICell, K) * PseudoThickCellKp1 + + SpecVol(ICell, K + 1) * PseudoThickCellK) / + (PseudoThickCellK + PseudoThickCellKp1); const Real DiffAlphaCellTop = VertDiff(ICell, K + 1) / (LocRhoSw * SpecVolCellTop); G = DT * DiffAlphaCellTop / - (LayerThickCellTop * LayerThickCell(ICell, K)); + (PseudoThickCellTop * PseudoThickCell(ICell, K)); } X = TracersOnCell(L, ICell, K); } From a3a6f2a3d757c2069a609d205f8faf8004add394 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 27 May 2026 09:53:51 -0400 Subject: [PATCH 09/19] Fix bugs after rebase --- .../omega/doc/devGuide/VerticalMixingCoeff.md | 14 +- .../doc/userGuide/VerticalMixingCoeff.md | 15 +- components/omega/src/ocn/VertMix.cpp | 96 +----------- components/omega/test/ocn/VertMixTest.cpp | 143 +----------------- 4 files changed, 20 insertions(+), 248 deletions(-) diff --git a/components/omega/doc/devGuide/VerticalMixingCoeff.md b/components/omega/doc/devGuide/VerticalMixingCoeff.md index 5ff40bc0505d..a59151ef27cc 100644 --- a/components/omega/doc/devGuide/VerticalMixingCoeff.md +++ b/components/omega/doc/devGuide/VerticalMixingCoeff.md @@ -40,7 +40,7 @@ VertMix: BaseShearValue: 0.005 RiCrit: 0.7 Exponent: 3.0 - RiSmoothLoops: 3 + RiSmoothLoops: 2 ``` ## Class Structure @@ -63,12 +63,12 @@ VertMix: - `ConvDiff`: Convective mixing coefficient (m²/s; Default: 1.0) - `ConvTriggerBVF`: Trigger threshold for convective mixing (Default: 0.0) -3. Shear Mixing: - - `EnableShearMix`: Flag to enable/disable shear mixing (Default: True) - - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) - - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) - - `ShearExponent`: Exponent parameter number for the LMB94 formulation (Default: 3.0) - - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 3) +3. Shear Instability Driven Mixing: + - `EnableShearMix`: Flag to enable/disable shear instability driven mixing (Default: True) + - `BaseShearValue`: Base values of maximum viscosity and diffusivity for the LMD94 interior shear instability driven mixing formulation (Default: 0.005) + - `RiCrit`: Critical Richardson number for the LMD94 formulation (Default: 0.7) + - `ShearExponent`: Exponent parameter number for the LMD94 formulation (Default: 3.0) + - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 2) ## Core Functionality (Vertical Mixing Coefficient Calculation) diff --git a/components/omega/doc/userGuide/VerticalMixingCoeff.md b/components/omega/doc/userGuide/VerticalMixingCoeff.md index a530b6781acc..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 [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These are linearly additive and are describe a bit -more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear mixing calculation. +background mixing value, (2) a convective instability mixing value, and (3) a Richardson number-dependent shear-instability-driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 shear instability driven mixing parameterization. These are linearly additive and are describe a bit +more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear instability driven mixing calculation. The user-configurable options are the following parameters in the yaml configuration file: @@ -26,14 +25,14 @@ VertMix: BaseShearValue: 0.005 # Base viscosity/diffusivity value RiCrit: 0.7 # Critical Richardson number Exponent: 3.0 # Richardson number exponent - RiSmoothLoops: 3 # Number of Richardson number smoothing loops + RiSmoothLoops: 2 # Number of Richardson number smoothing loops ``` ## Vertical Mixing Processes/Types ### 1. Background Mixing -A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior. +A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior and is assumed roughly equivalent to the globally averaged interior mixing from all sources (e.g., from internal wave breaking). ### 2. Convective Mixing @@ -51,14 +50,14 @@ This is different than some current implementations (i.e. in MPAS-Ocean and the ### 3. Shear-Instability-Driven Mixing -Mixing induced by vertical velocity shear, implemented using the LMD94 scheme, through the gradient Richardson number (ratio of buoyancy to shear). +Mixing induced by vertical pseudo-velocity shear, implemented using the LMD94 scheme, through the gradient Richardson number (ratio of buoyancy to shear). $$ \nu_{shear} = \kappa_{shear} = \begin{cases} \nu_o \quad \text{ if } Ri_g < 0\\ -\nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 < Ri_g < Ri_{crit}\\ -0.0 \quad \text{ if } Ri_{crit} < Ri_g +\nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 \leq Ri_g < Ri_{crit}\\ +0.0 \quad \text{ if } Ri_{crit} \leq Ri_g \end{cases} $$ diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index b186f01ae0b6..7066210c5746 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -223,7 +223,7 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, OMEGA_SCOPE(LocBackVisc, BackVisc); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - OMEGA_SCOPE(NVertLayers, VCoord->NVertLayers); + // OMEGA_SCOPE(NVertLayers, VCoord->NVertLayers); /// First, initialize VertDiff and VertVisc to background values parallelForOuter( @@ -260,97 +260,10 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, - KChunk, BruntVaisalaFreqSq); LocComputeGradRichardsonNum( LocGradRichNum, ICell, KChunk, NormalVelocity, TangentialVelocity, BruntVaisalaFreqSq); }); - LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); - LocGradRichNum(ICell, NVertLayers) = - LocGradRichNum(ICell, NVertLayers - 1); - }); - deepCopy(GradRichNumSmoothed, LocGradRichNum); - for (int SmoothLoop = 0; - SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - if (SmoothLoop == 0) - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - LocGradRichNum); - else - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - GradRichNumSmoothed); - }); - }); - } - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, - KChunk, GradRichNumSmoothed); - }); - }); - } else if (LocComputeVertMixShear.Enabled) { - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeGradRichardsonNum( - LocGradRichNum, ICell, KChunk, NormalVelocity, - TangentialVelocity, BruntVaisalaFreqSq); - }); - LocGradRichNum(ICell, 0) = LocGradRichNum(ICell, 1); - LocGradRichNum(ICell, NVertLayers) = - LocGradRichNum(ICell, NVertLayers - 1); - }); - deepCopy(GradRichNumSmoothed, LocGradRichNum); - for (int SmoothLoop = 0; - SmoothLoop < LocComputeVertMixShear.RiSmoothLoops; ++SmoothLoop) { - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - if (SmoothLoop == 0) - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - LocGradRichNum); - else - LocOneTwoOneFilter(GradRichNumSmoothed, ICell, KChunk, - GradRichNumSmoothed); - }); - }); - } - parallelForOuter( - "VertMix-ConvPlusShear", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(I4 ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRangeChunked(KMin, KMax); - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - LocComputeVertMixShear(LocVertDiff, LocVertVisc, ICell, - KChunk, GradRichNumSmoothed); teamBarrier(Team); @@ -617,7 +530,7 @@ void VertMix::applyVelVertMixImplicit( } Pacer::stop("Tend:velocityVertMix", 1); } -} +} // applyVelVertMixImplicit // Apply implicit tracer vertical mixing void VertMix::applyTracerVertMixImplicit( @@ -719,7 +632,7 @@ void VertMix::applyTracerVertMixImplicit( Pacer::stop("Tend:tracerVertMix", 1); } -} // end all tendency compute +} // applyTracerVertMixImplicit /// Apply implicit vertical mixing to velocities and tracers void VertMix::VertMixImplicit(OceanState *State, AuxiliaryState *AuxState, @@ -780,6 +693,7 @@ void VertMix::VertMixImplicit(OceanState *State, AuxiliaryState *AuxState, // Apply implicit mixing to tracers applyTracerVertMixImplicit(State, AuxState, TracerArray, NTracers, TimeLevel, TimeLevel); -} + +} // VertMixImplicit } // namespace OMEGA diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index c25338cbabc8..75bd1d2ed67b 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -128,6 +128,7 @@ void testGradRichNum() { OMEGA_SCOPE(DcEdge, Mesh->DcEdge); OMEGA_SCOPE(DvEdge, Mesh->DvEdge); OMEGA_SCOPE(CellsOnCell, Mesh->CellsOnCell); + // OMEGA_SCOPE(CellsOnEdge, Mesh->CellsOnEdge); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); @@ -205,9 +206,6 @@ void testGradRichNum() { TestVertMix->computeVertMix(NormalVelEdge, TangVelEdge, BruntVaisalaFreqSqCell); - const auto &MinLayerCell = VCoord->MinLayerCell; - const auto &MaxLayerCell = VCoord->MaxLayerCell; - /// Check all array values against expected value int NumMismatches = 0; OMEGA_SCOPE(GradRichNum, TestVertMix->GradRichNum); @@ -262,15 +260,6 @@ void testGradRichNum() { // If test fails, print bad values and abort if (NumMismatches != 0) { - auto GradRichNumH = createHostMirrorCopy(GradRichNum); - for (int I = 0; I < NCellsAll; ++I) { - for (int K = 0; K < NVertLayers; ++K) { - if (!isApprox(GradRichNumH(I, K), RiExpValue, RTol)) - LOG_ERROR("TestVertMix: GradRichNum Bad Value: " - "GradRichNum({},{}) = {}; Expected {}", - I, K, GradRichNumH(I, K), RiExpValue); - } - } ABORT_ERROR("TestVertMix: GradRichNum FAIL with {} bad values", NumMismatches); } else { @@ -364,28 +353,6 @@ void testOneTwoOneFilter() { // If test fails, print bad values and abort if (NumMismatches != 0) { - auto GradRichNumH = createHostMirrorCopy(GradRichNum); - auto GradRichNumSmoothedH = createHostMirrorCopy(GradRichNumSmoothed); - for (int I = 0; I < NCellsAll; ++I) { - for (int K = 0; K < NVertLayers; ++K) { - if (K >= 1 && K < NVertLayers - 1) { - // Interior layers should be smoothed to 0.0 - if (!isApprox(GradRichNumSmoothedH(I, K), 0.0, RTol)) - LOG_ERROR("TestVertMix: GradRichNumSmoothed Bad Value: " - "GradRichNumSmoothed({},{}) = {}; Expected {}", - I, K, GradRichNumSmoothedH(I, K), 0.0); - } else { - // Boundary layers (K==0 or K==NVertLayers-1) should be copied - // from input - if (!isApprox(GradRichNumSmoothedH(I, K), GradRichNumH(I, K), - RTol)) - LOG_ERROR("TestVertMix: GradRichNumSmoothed Bad Value: " - "GradRichNumSmoothed({},{}) = {}; Expected {}", - I, K, GradRichNumSmoothedH(I, K), - GradRichNumH(I, K)); - } - } - } ABORT_ERROR("TestVertMix: GradRichNumSmoothed FAIL with {} bad values", NumMismatches); } else { @@ -609,30 +576,6 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { - auto VertViscOutH = createHostMirrorCopy(VertViscOut); - for (int I = 0; I < NCellsAll; ++I) { - for (int K = 0; K < NVertLayers; ++K) { - if (K == 0) { - // Surface should be 0.0 - if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " - "VertVisc({},{}) = {}; Expected {}", - I, K, VertViscOutH(I, K), 0.0_Real); - } else if (K < 30) { - // Interior layers - if (!isApprox(VertViscOutH(I, K), VertConvExp, RTol)) - LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " - "VertVisc({},{}) = {}; Expected {}", - I, K, VertViscOutH(I, K), VertConvExp); - } else { - // Interior layers - if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixConv: VertVisc Bad Value: " - "VertVisc({},{}) = {}; Expected {}", - I, K, VertViscOutH(I, K), 0.0_Real); - } - } - } ABORT_ERROR("TestVertMixConv: VertVisc FAIL with {} bad values", NumMismatches); } else { @@ -672,30 +615,6 @@ void testConvVertMix() { NumMismatches); if (NumMismatches != 0) { - auto VertDiffOutH = createHostMirrorCopy(VertDiffOut); - for (int I = 0; I < NCellsAll; ++I) { - for (int K = 0; K < NVertLayers; ++K) { - if (K == 0) { - // Surface should be 0.0 - if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " - "VertDiff({},{}) = {}; Expected {}", - I, K, VertDiffOutH(I, K), 0.0_Real); - } else if (K < 30) { - // Interior layers - if (!isApprox(VertDiffOutH(I, K), VertConvExp, RTol)) - LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " - "VertDiff({},{}) = {}; Expected {}", - I, K, VertDiffOutH(I, K), VertConvExp); - } else { - // Interior layers - if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixConv: VertDiff Bad Value: " - "VertDiff({},{}) = {}; Expected {}", - I, K, VertDiffOutH(I, K), 0.0_Real); - } - } - } ABORT_ERROR("TestVertMixConv: VertDiff FAIL with {} bad values", NumMismatches); } else { @@ -795,36 +714,6 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { - auto VertViscOutH = createHostMirrorCopy(VertViscOut); - for (int I = 0; I < NCellsAll; ++I) { - for (int K = 0; K < NVertLayers; ++K) { - if (K == 0) { - // Surface should be 0.0 - if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " - "VertVisc({},{}) = {}; Expected {}", - I, K, VertViscOutH(I, K), 0.0_Real); - } else if (K < 20) { - // Interior layers - if (!isApprox(VertViscOutH(I, K), VertShearBaseExp, RTol)) - LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " - "VertVisc({},{}) = {}; Expected {}", - I, K, VertViscOutH(I, K), VertShearBaseExp); - } else if (K >= 20 && K < 40) { - // Interior layers - if (!isApprox(VertViscOutH(I, K), VertShearExp, RTol)) - LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " - "VertVisc({},{}) = {}; Expected {}", - I, K, VertViscOutH(I, K), VertShearExp); - } else { - // Interior layers - if (!isApprox(VertViscOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixShear: VertVisc Bad Value: " - "VertVisc({},{}) = {}; Expected {}", - I, K, VertViscOutH(I, K), 0.0_Real); - } - } - } ABORT_ERROR("TestVertMixShear: VertVisc FAIL with {} bad values", NumMismatches); } else { @@ -868,36 +757,6 @@ void testShearVertMix() { NumMismatches); if (NumMismatches != 0) { - auto VertDiffOutH = createHostMirrorCopy(VertDiffOut); - for (int I = 0; I < NCellsAll; ++I) { - for (int K = 0; K < NVertLayers; ++K) { - if (K == 0) { - // Surface should be 0.0 - if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " - "VertDiff({},{}) = {}; Expected {}", - I, K, VertDiffOutH(I, K), 0.0_Real); - } else if (K < 20) { - // Interior layers - if (!isApprox(VertDiffOutH(I, K), VertShearBaseExp, RTol)) - LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " - "VertDiff({},{}) = {}; Expected {}", - I, K, VertDiffOutH(I, K), VertShearBaseExp); - } else if (K >= 20 && K < 40) { - // Interior layers - if (!isApprox(VertDiffOutH(I, K), VertShearExp, RTol)) - LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " - "VertDiff({},{}) = {}; Expected {}", - I, K, VertDiffOutH(I, K), VertShearExp); - } else { - // Interior layers - if (!isApprox(VertDiffOutH(I, K), 0.0_Real, RTol)) - LOG_ERROR("TestVertMixShear: VertDiff Bad Value: " - "VertDiff({},{}) = {}; Expected {}", - I, K, VertDiffOutH(I, K), 0.0_Real); - } - } - } ABORT_ERROR("TestVertMixShear: VertDiff FAIL with {} bad values", NumMismatches); } else { From 45fa17544b26de9b8eff88d1afd138699d584ced Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 27 May 2026 10:28:34 -0400 Subject: [PATCH 10/19] Fix tridiag policy change in applyVertMix --- components/omega/src/ocn/VertMix.cpp | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 7066210c5746..fcf84b7445bd 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -474,11 +474,10 @@ void VertMix::applyVelVertMixImplicit( const auto &VertVisc = VertMixInstance->VertVisc; const I4 NVertLayers = VCoord->NVertLayers; - TeamPolicy Policy = - TriDiagDiffSolver::makeTeamPolicy(Mesh->NEdgesAll, NVertLayers); + auto LConfig = TriDiagSolver::makeLaunchConfig(Mesh->NEdgesAll, NVertLayers); - Kokkos::parallel_for( - Policy, KOKKOS_LAMBDA(const TeamMember &Team) { + parallelForOuter( + LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { const int IStart = Team.league_rank() * VecLength; const int ILen = Kokkos::max( 0, Kokkos::min(VecLength, LocNEdgesAll - IStart)); @@ -486,8 +485,7 @@ void VertMix::applyVelVertMixImplicit( TriDiagDiffScratch Scratch(Team, NVertLayers); // Construct a tri-diag diffusion matrix and RHS - Kokkos::parallel_for( - TeamThreadRange(Team, NVertLayers), [=](int K) { + parallelForInner(Team, NVertLayers, [=](int K) { for (int IVec = 0; IVec < VecLength; ++IVec) { const int IEdge = IStart + IVec; @@ -573,12 +571,11 @@ void VertMix::applyTracerVertMixImplicit( const auto &VertDiff = VertMixInstance->VertDiff; const I4 NVertLayers = VCoord->NVertLayers; - TeamPolicy Policy = - TriDiagDiffSolver::makeTeamPolicy(Mesh->NCellsAll, NVertLayers); + auto LConfig = TriDiagSolver::makeLaunchConfig(Mesh->NCellsAll, NVertLayers); for (int L = 0; L < NTracers; ++L) { - Kokkos::parallel_for( - Policy, KOKKOS_LAMBDA(const TeamMember &Team) { + parallelForOuter( + LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { const int IStart = Team.league_rank() * VecLength; const int ILen = Kokkos::max( 0, Kokkos::min(VecLength, LocNCellsAll - IStart)); @@ -586,8 +583,7 @@ void VertMix::applyTracerVertMixImplicit( TriDiagDiffScratch Scratch(Team, NVertLayers); // Construct a tri-diag diffusion matrix and RHS - Kokkos::parallel_for( - TeamThreadRange(Team, NVertLayers), [=](int K) { + parallelForInner(Team, NVertLayers, [=](int K) { for (int IVec = 0; IVec < VecLength; ++IVec) { const int ICell = IStart + IVec; From fb349bcea3bca71ba1b19a53d022b91b38bc5d40 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 27 May 2026 10:41:14 -0400 Subject: [PATCH 11/19] Fix linting issue --- components/omega/src/ocn/VertMix.cpp | 83 ++++++++++++++-------------- 1 file changed, 42 insertions(+), 41 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index fcf84b7445bd..6f8f1c763355 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -474,7 +474,8 @@ void VertMix::applyVelVertMixImplicit( const auto &VertVisc = VertMixInstance->VertVisc; const I4 NVertLayers = VCoord->NVertLayers; - auto LConfig = TriDiagSolver::makeLaunchConfig(Mesh->NEdgesAll, NVertLayers); + auto LConfig = + TriDiagSolver::makeLaunchConfig(Mesh->NEdgesAll, NVertLayers); parallelForOuter( LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { @@ -486,26 +487,25 @@ void VertMix::applyVelVertMixImplicit( // Construct a tri-diag diffusion matrix and RHS parallelForInner(Team, NVertLayers, [=](int K) { - for (int IVec = 0; IVec < VecLength; ++IVec) { - const int IEdge = IStart + IVec; - - if (IEdge >= LocNEdgesAll) { - Scratch.G(K, IVec) = 0._Real; - Scratch.H(K, IVec) = 1._Real; - Scratch.X(K, IVec) = 0._Real; - continue; - } - - Real G, H, X; - LocVelVertMixSetup(IEdge, K, DT, SpecVol, - PseudoThickCell, VertVisc, - NormalVelEdge, G, H, X); - - Scratch.G(K, IVec) = G; - Scratch.H(K, IVec) = H; - Scratch.X(K, IVec) = X; - } - }); + for (int IVec = 0; IVec < VecLength; ++IVec) { + const int IEdge = IStart + IVec; + + if (IEdge >= LocNEdgesAll) { + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + Real G, H, X; + LocVelVertMixSetup(IEdge, K, DT, SpecVol, PseudoThickCell, + VertVisc, NormalVelEdge, G, H, X); + + Scratch.G(K, IVec) = G; + Scratch.H(K, IVec) = H; + Scratch.X(K, IVec) = X; + } + }); // Solve the tri-diag diffusion system Team.team_barrier(); @@ -571,7 +571,8 @@ void VertMix::applyTracerVertMixImplicit( const auto &VertDiff = VertMixInstance->VertDiff; const I4 NVertLayers = VCoord->NVertLayers; - auto LConfig = TriDiagSolver::makeLaunchConfig(Mesh->NCellsAll, NVertLayers); + auto LConfig = + TriDiagSolver::makeLaunchConfig(Mesh->NCellsAll, NVertLayers); for (int L = 0; L < NTracers; ++L) { parallelForOuter( @@ -584,25 +585,25 @@ void VertMix::applyTracerVertMixImplicit( // Construct a tri-diag diffusion matrix and RHS parallelForInner(Team, NVertLayers, [=](int K) { - for (int IVec = 0; IVec < VecLength; ++IVec) { - const int ICell = IStart + IVec; - - if (ICell >= LocNCellsAll) { - Scratch.G(K, IVec) = 0._Real; - Scratch.H(K, IVec) = 1._Real; - Scratch.X(K, IVec) = 0._Real; - continue; - } - - Real G, H, X; - LocTracerVertMixSetup(L, ICell, K, DT, SpecVol, - PseudoThickCell, VertDiff, - TracerArray, G, H, X); - Scratch.G(K, IVec) = G; - Scratch.H(K, IVec) = H; - Scratch.X(K, IVec) = X; - } - }); + for (int IVec = 0; IVec < VecLength; ++IVec) { + const int ICell = IStart + IVec; + + if (ICell >= LocNCellsAll) { + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + Real G, H, X; + LocTracerVertMixSetup(L, ICell, K, DT, SpecVol, + PseudoThickCell, VertDiff, + TracerArray, G, H, X); + Scratch.G(K, IVec) = G; + Scratch.H(K, IVec) = H; + Scratch.X(K, IVec) = X; + } + }); // Solve the tri-diag diffusion system Team.team_barrier(); From ce54fde870e144727efb3a96a6db8009840d5afc Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 1 Jun 2026 21:07:13 -0400 Subject: [PATCH 12/19] Fix bugs after rebase --- components/omega/configs/Default.yml | 2 +- components/omega/src/ocn/OceanInit.cpp | 39 +------------------ components/omega/src/ocn/Tendencies.cpp | 13 ------- components/omega/src/ocn/Tendencies.h | 11 ------ components/omega/src/ocn/VertMix.cpp | 1 - .../src/timeStepping/RungeKutta4Stepper.h | 5 ++- components/omega/test/ocn/VertMixTest.cpp | 1 - 7 files changed, 6 insertions(+), 66 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 92299ae6f259..08164bdaa373 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -94,7 +94,7 @@ Omega: BaseShearValue: 0.005 RiCrit: 0.7 Exponent: 3.0 - RiSmoothLoops: 3 + RiSmoothLoops: 2 IOStreams: HorzMeshIn: UsePointerFile: false diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index e7fe9938401e..c0393e692add 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -107,44 +107,6 @@ int ocnInit(MPI_Comm Comm ///< [in] ocean MPI communicator TimeStepper *DefStepper = TimeStepper::getDefault(); Clock *ModelClock = DefStepper->getClock(); - // Initialize IOStreams - this does not yet validate the contents - // of each file, only creates streams from Config - IOStream::init(ModelClock); - - IO::init(Comm); - Field::init(ModelClock); - Decomp::init(); - - Err = Halo::init(); - if (Err != 0) { - ABORT_ERROR("ocnInit: Error initializing default halo"); - } - - HorzMesh::init(); - VertCoord::init(); - Tracers::init(); - VertAdv::init(); - AuxiliaryState::init(); - Eos::init(); - PressureGrad::init(); - VertMix::init(); - Tendencies::init(); - - // Validate SurfaceTracerRestoring configuration - Tendencies *DefTend = Tendencies::getDefault(); - if (DefTend->SurfaceTracerRestoring.Enabled && - DefTend->SurfaceTracerRestoring.NTracersToRestore == 0) { - ABORT_ERROR("OceanInit: SurfaceTracerRestoring is enabled but " - "TracersToRestore is empty"); - } - - TimeStepper::init2(); - - Err = OceanState::init(); - if (Err != 0) { - ABORT_ERROR("ocnInit: Error initializing default state"); - } - // Now that all fields have been defined, validate all the streams // contents bool StreamsValid = IOStream::validateAll(); @@ -239,6 +201,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) { AuxiliaryState::init(); Eos::init(); PressureGrad::init(); + VertMix::init(); Tendencies::init(); // Validate SurfaceTracerRestoring configuration diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index c0feea0d8652..089daa8d0c3b 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -441,19 +441,6 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies TimeStepIn, Options, CustomTendencyType{}, CustomTendencyType{}) {} -Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies - const HorzMesh *Mesh, ///< [in] Horizontal mesh - VertCoord *VCoord, ///< [in] Vertical coordinate - VertAdv *VAdv, ///< [in] Vertical advection - PressureGrad *PGrad, ///< [in] Pressure gradient - Eos *EqState, ///< [in] Equation of state - int NTracersIn, ///< [in] Number of tracers - TimeInterval TimeStepIn, ///< [in] Time step - Config *Options) ///< [in] Configuration options - : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, - VertMix::getInstance(), NTracersIn, TimeStepIn, Options, - CustomTendencyType{}, CustomTendencyType{}) {} - //------------------------------------------------------------------------------ // Compute tendencies for the pseudo-thickness equation void Tendencies::computePseudoThicknessTendenciesOnly( diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index d2aa85ffd2ef..ac826f915a85 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -177,17 +177,6 @@ class Tendencies { CustomTendencyType InCustomThicknessTend, CustomTendencyType InCustomVelocityTend); - Tendencies(const std::string &Name, ///< [in] Name for tendencies - const HorzMesh *Mesh, ///< [in] Horizontal mesh - VertCoord *VCoord, ///< [in] Vertical coordinate - VertAdv *VAdv, ///< [in] Vertical advection - PressureGrad *PGrad, ///< [in] Pressure gradient - Eos *EqState, ///< [in] Equation of state - int NTracersIn, ///< [in] Number of tracers - TimeInterval TimeStep, ///< [in] Time step - Config *Options ///< [in] Configuration options - ); - Tendencies(const std::string &Name, ///< [in] Name for tendencies const HorzMesh *Mesh, ///< [in] Horizontal mesh VertCoord *VCoord, ///< [in] Vertical coordinate diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 6f8f1c763355..4262a5932e34 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -223,7 +223,6 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, OMEGA_SCOPE(LocBackVisc, BackVisc); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - // OMEGA_SCOPE(NVertLayers, VCoord->NVertLayers); /// First, initialize VertDiff and VertVisc to background values parallelForOuter( diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.h b/components/omega/src/timeStepping/RungeKutta4Stepper.h index d4dbe8f8072b..16efa1acd9ba 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.h +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.h @@ -46,7 +46,10 @@ class RungeKutta4Stepper : public TimeStepper { Real RKB[NStages]; Real RKC[NStages]; - // Projection factors + // Projection factors for each RK stage. These are the RKC coefficients + // shifted by one RK stage and represent the fraction of the timestep over + // which fields are projected from the old state to the state at the end of a + // RK stage. Real RKProj[NStages]; }; diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index 75bd1d2ed67b..d1bf97015798 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -128,7 +128,6 @@ void testGradRichNum() { OMEGA_SCOPE(DcEdge, Mesh->DcEdge); OMEGA_SCOPE(DvEdge, Mesh->DvEdge); OMEGA_SCOPE(CellsOnCell, Mesh->CellsOnCell); - // OMEGA_SCOPE(CellsOnEdge, Mesh->CellsOnEdge); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); From 9b82c47227e92295c6782a6faab940418c0029c0 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 1 Jun 2026 21:18:37 -0400 Subject: [PATCH 13/19] Fix GPU compile error --- components/omega/src/ocn/VertMix.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 4262a5932e34..fcf2abfa996b 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -512,8 +512,7 @@ void VertMix::applyVelVertMixImplicit( Team.team_barrier(); // Store the solution vector X - Kokkos::parallel_for( - TeamThreadRange(Team, NVertLayers), [=](int K) { + parallelForInner(Team, NVertLayers, [=](int K) { for (int IVec = 0; IVec < ILen; ++IVec) { const int IEdge = IStart + IVec; @@ -610,8 +609,7 @@ void VertMix::applyTracerVertMixImplicit( Team.team_barrier(); // Store the solution vector X - Kokkos::parallel_for( - TeamThreadRange(Team, NVertLayers), [=](int K) { + parallelForInner(Team, NVertLayers, [=](int K) { for (int IVec = 0; IVec < ILen; ++IVec) { const int ICell = IStart + IVec; From 79c0a77bdcf1fe93cb928141a7d324fadeff28e1 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 1 Jun 2026 21:23:05 -0400 Subject: [PATCH 14/19] Fix linting issue --- components/omega/src/ocn/VertMix.cpp | 36 ++++++++++++++-------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index fcf2abfa996b..79323e9eff50 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -513,15 +513,15 @@ void VertMix::applyVelVertMixImplicit( // Store the solution vector X parallelForInner(Team, NVertLayers, [=](int K) { - for (int IVec = 0; IVec < ILen; ++IVec) { - const int IEdge = IStart + IVec; - - if (K >= MinLayerEdgeBot(IEdge) && - K <= MaxLayerEdgeTop(IEdge)) { - NormalVelEdge(IEdge, K) = Scratch.X(K, IVec); - } - } - }); + for (int IVec = 0; IVec < ILen; ++IVec) { + const int IEdge = IStart + IVec; + + if (K >= MinLayerEdgeBot(IEdge) && + K <= MaxLayerEdgeTop(IEdge)) { + NormalVelEdge(IEdge, K) = Scratch.X(K, IVec); + } + } + }); }); } Pacer::stop("Tend:velocityVertMix", 1); @@ -610,15 +610,15 @@ void VertMix::applyTracerVertMixImplicit( // Store the solution vector X parallelForInner(Team, NVertLayers, [=](int K) { - for (int IVec = 0; IVec < ILen; ++IVec) { - const int ICell = IStart + IVec; - - if (K >= MinLayerCell(ICell) && - K <= MaxLayerCell(ICell)) { - TracerArray(L, ICell, K) = Scratch.X(K, IVec); - } - } - }); + for (int IVec = 0; IVec < ILen; ++IVec) { + const int ICell = IStart + IVec; + + if (K >= MinLayerCell(ICell) && + K <= MaxLayerCell(ICell)) { + TracerArray(L, ICell, K) = Scratch.X(K, IVec); + } + } + }); }); } // for L From 5ea6e16f4c4334b1e651f57cd9216b145306feb5 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Thu, 4 Jun 2026 12:57:37 -0400 Subject: [PATCH 15/19] Fix bugs for non-flat bottom topography and minor updates --- components/omega/src/ocn/VertMix.cpp | 47 ++++++++--- components/omega/src/ocn/VertMix.h | 119 +++++++++++++-------------- 2 files changed, 95 insertions(+), 71 deletions(-) diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 79323e9eff50..acd486a0c448 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -472,7 +472,7 @@ void VertMix::applyVelVertMixImplicit( const auto &SpecVol = EosInstance->SpecVol; const auto &VertVisc = VertMixInstance->VertVisc; - const I4 NVertLayers = VCoord->NVertLayers; + const int NVertLayers = VCoord->NVertLayers; auto LConfig = TriDiagSolver::makeLaunchConfig(Mesh->NEdgesAll, NVertLayers); @@ -490,6 +490,18 @@ void VertMix::applyVelVertMixImplicit( const int IEdge = IStart + IVec; if (IEdge >= LocNEdgesAll) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + + if (K < KMin || K > KMax) { + // Fill values Scratch.G(K, IVec) = 0._Real; Scratch.H(K, IVec) = 1._Real; Scratch.X(K, IVec) = 0._Real; @@ -497,8 +509,9 @@ void VertMix::applyVelVertMixImplicit( } Real G, H, X; - LocVelVertMixSetup(IEdge, K, DT, SpecVol, PseudoThickCell, - VertVisc, NormalVelEdge, G, H, X); + LocVelVertMixSetup(IEdge, K, KMin, KMax, DT, SpecVol, + PseudoThickCell, VertVisc, + NormalVelEdge, G, H, X); Scratch.G(K, IVec) = G; Scratch.H(K, IVec) = H; @@ -568,7 +581,7 @@ void VertMix::applyTracerVertMixImplicit( const auto &SpecVol = EosInstance->SpecVol; const auto &VertDiff = VertMixInstance->VertDiff; - const I4 NVertLayers = VCoord->NVertLayers; + const int NVertLayers = VCoord->NVertLayers; auto LConfig = TriDiagSolver::makeLaunchConfig(Mesh->NCellsAll, NVertLayers); @@ -587,6 +600,18 @@ void VertMix::applyTracerVertMixImplicit( const int ICell = IStart + IVec; if (ICell >= LocNCellsAll) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + + if (K < KMin || K > KMax) { + // Fill values Scratch.G(K, IVec) = 0._Real; Scratch.H(K, IVec) = 1._Real; Scratch.X(K, IVec) = 0._Real; @@ -594,9 +619,9 @@ void VertMix::applyTracerVertMixImplicit( } Real G, H, X; - LocTracerVertMixSetup(L, ICell, K, DT, SpecVol, - PseudoThickCell, VertDiff, - TracerArray, G, H, X); + LocTracerVertMixSetup(L, ICell, K, KMin, KMax, DT, + SpecVol, PseudoThickCell, + VertDiff, TracerArray, G, H, X); Scratch.G(K, IVec) = G; Scratch.H(K, IVec) = H; Scratch.X(K, IVec) = X; @@ -652,16 +677,16 @@ void VertMix::VertMixImplicit(OceanState *State, AuxiliaryState *AuxState, // TODO: Temporary handling of computation of tangential velocity // Compute tangential velocity - OMEGA_SCOPE(MinLayerEdgeTop, VCoord->MinLayerEdgeTop); - OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); OMEGA_SCOPE(LocTangentialVelocity, TangentialVelocity); TangentialReconOnEdge TanReconEdge(Mesh); parallelForOuter( {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { - const int KMin = MinLayerEdgeTop(IEdge); - const int KMax = MaxLayerEdgeBot(IEdge); + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); const int KRange = vertRangeChunked(KMin, KMax); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index d9c0dc187e67..8ea2d311e6f2 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -240,7 +240,7 @@ class VelVertMixSetupOnEdge { /// The functor takes edge index, vertical chunk index, and arrays for /// layer specific volume, layer thickness on edge, /// interface pressure, and outputs tendency array - KOKKOS_FUNCTION void operator()(I4 IEdge, I4 K, Real DT, + KOKKOS_FUNCTION void operator()(I4 IEdge, I4 K, I4 KMin, I4 KMax, Real DT, const Array2DReal &SpecVol, const Array2DReal &PseudoThickCell, const Array2DReal &VertVisc, @@ -250,42 +250,41 @@ class VelVertMixSetupOnEdge { const I4 JCell0 = CellsOnEdge(IEdge, 0); const I4 JCell1 = CellsOnEdge(IEdge, 1); - const I4 KMin = MinLayerEdgeBot(IEdge); - const I4 KMax = MaxLayerEdgeTop(IEdge); - + // Fill values G = 0.0_Real; H = 1.0_Real; - - if (K < KMin || K > KMax) { - X = 0.0_Real; - } else { - if (K < KMax) { - const Real PseudoThickEdgeK = - 0.5_Real * - (PseudoThickCell(JCell0, K) + PseudoThickCell(JCell1, K)); - const Real PseudoThickEdgeKp1 = - 0.5_Real * (PseudoThickCell(JCell0, K + 1) + - PseudoThickCell(JCell1, K + 1)); - - const Real PseudoThickEdgeTop = - 0.5_Real * (PseudoThickEdgeK + PseudoThickEdgeKp1); - - // Interpolation from cell center to top using - // the two-point linear interpolation - const Real SpecVolEdgeTop = - (0.5_Real * (SpecVol(JCell0, K) + SpecVol(JCell1, K)) * - PseudoThickEdgeKp1 + - 0.5_Real * (SpecVol(JCell0, K + 1) + SpecVol(JCell1, K + 1)) * - PseudoThickEdgeK) / - (PseudoThickEdgeK + PseudoThickEdgeKp1); - - const Real ViscAlphaEdgeTop = - 0.5_Real * (VertVisc(JCell0, K + 1) + VertVisc(JCell1, K + 1)) / - (LocRhoSw * SpecVolEdgeTop); - - G = DT * ViscAlphaEdgeTop / (PseudoThickEdgeTop * PseudoThickEdgeK); - } - X = NormalVelEdge(IEdge, K); + X = 0.0_Real; + + const Real PseudoThickEdgeK = + 0.5_Real * (PseudoThickCell(JCell0, K) + PseudoThickCell(JCell1, K)); + + // Row-scaled conservative form. + // Unknown is u^{n+1}. + H = PseudoThickEdgeK; + X = PseudoThickEdgeK * NormalVelEdge(IEdge, K); + + if (K < KMax) { + const Real PseudoThickEdgeKp1 = + 0.5_Real * + (PseudoThickCell(JCell0, K + 1) + PseudoThickCell(JCell1, K + 1)); + + const Real PseudoThickEdgeTop = + 0.5_Real * (PseudoThickEdgeK + PseudoThickEdgeKp1); + + // Interpolation from cell center to top using + // the two-point linear interpolation + const Real SpecVolEdgeTop = + (0.5_Real * (SpecVol(JCell0, K) + SpecVol(JCell1, K)) * + PseudoThickEdgeKp1 + + 0.5_Real * (SpecVol(JCell0, K + 1) + SpecVol(JCell1, K + 1)) * + PseudoThickEdgeK) / + (PseudoThickEdgeK + PseudoThickEdgeKp1); + + const Real ViscAlphaEdgeTop = + 0.5_Real * (VertVisc(JCell0, K + 1) + VertVisc(JCell1, K + 1)) / + (LocRhoSw * SpecVolEdgeTop); + + G = DT * ViscAlphaEdgeTop / PseudoThickEdgeTop; } } @@ -305,43 +304,43 @@ class TracerVertMixSetupOnCell { TracerVertMixSetupOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); - KOKKOS_FUNCTION void operator()(I4 L, I4 ICell, I4 K, Real DT, - const Array2DReal &SpecVol, + KOKKOS_FUNCTION void operator()(I4 L, I4 ICell, I4 K, I4 KMin, I4 KMax, + Real DT, const Array2DReal &SpecVol, const Array2DReal &PseudoThickCell, const Array2DReal &VertDiff, const Array3DReal &TracersOnCell, Real &G, Real &H, Real &X) const { - const I4 KMin = MinLayerCell(ICell); - const I4 KMax = MaxLayerCell(ICell); - + // Fill values G = 0.0_Real; H = 1.0_Real; + X = 0.0_Real; - if (K < KMin || K > KMax) { - X = 0.0_Real; - } else { - if (K < KMax) { - const Real PseudoThickCellK = PseudoThickCell(ICell, K); - const Real PseudoThickCellKp1 = PseudoThickCell(ICell, K + 1); + const Real PseudoThickCellK = PseudoThickCell(ICell, K); - const Real PseudoThickCellTop = - 0.5_Real * (PseudoThickCellK + PseudoThickCellKp1); + // Row-scaled conservative form. + // Unknown is phi^{n+1}. + H = PseudoThickCellK; + X = PseudoThickCellK * TracersOnCell(L, ICell, K); - // Interpolation from cell center to top using - // the two-point linear interpolation - const Real SpecVolCellTop = - (SpecVol(ICell, K) * PseudoThickCellKp1 + - SpecVol(ICell, K + 1) * PseudoThickCellK) / - (PseudoThickCellK + PseudoThickCellKp1); + if (K < KMax) { - const Real DiffAlphaCellTop = - VertDiff(ICell, K + 1) / (LocRhoSw * SpecVolCellTop); + const Real PseudoThickCellKp1 = PseudoThickCell(ICell, K + 1); - G = DT * DiffAlphaCellTop / - (PseudoThickCellTop * PseudoThickCell(ICell, K)); - } - X = TracersOnCell(L, ICell, K); + const Real PseudoThickCellTop = + 0.5_Real * (PseudoThickCellK + PseudoThickCellKp1); + + // Interpolation from cell center to top using + // the two-point linear interpolation + const Real SpecVolCellTop = + (SpecVol(ICell, K) * PseudoThickCellKp1 + + SpecVol(ICell, K + 1) * PseudoThickCellK) / + (PseudoThickCellK + PseudoThickCellKp1); + + const Real DiffAlphaCellTop = + VertDiff(ICell, K + 1) / (LocRhoSw * SpecVolCellTop); + + G = DT * DiffAlphaCellTop / PseudoThickCellTop; } } From cdfbcbf7fcaa6d0d4ffb01112a7acd2b7f28fcec Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Mon, 29 Jun 2026 23:26:15 -0400 Subject: [PATCH 16/19] =?UTF-8?q?Address=20reviewers=E2=80=99=20comments?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Renaming cell-interface variables from 'Top' to 'Bot'. --- components/omega/src/ocn/VertMix.h | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index 8ea2d311e6f2..008476125c99 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -268,23 +268,23 @@ class VelVertMixSetupOnEdge { 0.5_Real * (PseudoThickCell(JCell0, K + 1) + PseudoThickCell(JCell1, K + 1)); - const Real PseudoThickEdgeTop = + const Real PseudoThickEdgeBot = 0.5_Real * (PseudoThickEdgeK + PseudoThickEdgeKp1); // Interpolation from cell center to top using // the two-point linear interpolation - const Real SpecVolEdgeTop = + const Real SpecVolEdgeBot = (0.5_Real * (SpecVol(JCell0, K) + SpecVol(JCell1, K)) * PseudoThickEdgeKp1 + 0.5_Real * (SpecVol(JCell0, K + 1) + SpecVol(JCell1, K + 1)) * PseudoThickEdgeK) / (PseudoThickEdgeK + PseudoThickEdgeKp1); - const Real ViscAlphaEdgeTop = + const Real ViscAlphaEdgeBot = 0.5_Real * (VertVisc(JCell0, K + 1) + VertVisc(JCell1, K + 1)) / - (LocRhoSw * SpecVolEdgeTop); + (LocRhoSw * SpecVolEdgeBot); - G = DT * ViscAlphaEdgeTop / PseudoThickEdgeTop; + G = DT * ViscAlphaEdgeBot / PseudoThickEdgeBot; } } @@ -327,20 +327,20 @@ class TracerVertMixSetupOnCell { const Real PseudoThickCellKp1 = PseudoThickCell(ICell, K + 1); - const Real PseudoThickCellTop = + const Real PseudoThickCellBot = 0.5_Real * (PseudoThickCellK + PseudoThickCellKp1); // Interpolation from cell center to top using // the two-point linear interpolation - const Real SpecVolCellTop = + const Real SpecVolCellBot = (SpecVol(ICell, K) * PseudoThickCellKp1 + SpecVol(ICell, K + 1) * PseudoThickCellK) / (PseudoThickCellK + PseudoThickCellKp1); - const Real DiffAlphaCellTop = - VertDiff(ICell, K + 1) / (LocRhoSw * SpecVolCellTop); + const Real DiffAlphaCellBot = + VertDiff(ICell, K + 1) / (LocRhoSw * SpecVolCellBot); - G = DT * DiffAlphaCellTop / PseudoThickCellTop; + G = DT * DiffAlphaCellBot / PseudoThickCellBot; } } From a14740a4b6d71bec039290eec5114e9777826167 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Tue, 30 Jun 2026 10:25:18 -0400 Subject: [PATCH 17/19] =?UTF-8?q?Address=20reviewers=E2=80=99=20comments?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Validate the VertMix tendency once at initialization instead of during each time step. - Abort the model if vertical mixing tendencies are enabled but Eos or VertMix has not been initialized. --- components/omega/src/ocn/Tendencies.cpp | 14 ++ components/omega/src/ocn/VertMix.cpp | 270 +++++++++++------------- 2 files changed, 141 insertions(+), 143 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 089daa8d0c3b..52a4c4c33eca 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -321,6 +321,7 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config TracerIdsToRestoreVec.size())); } + // Validate VertMix tendency Err += TendConfig.get("VelVertMixTendencyEnable", this->VMix->VelVertMixSetup.Enabled); CHECK_ERROR_ABORT( @@ -330,6 +331,19 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config this->VMix->TracerVertMixSetup.Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: TracerVertMixTendencyEnable not found in TendConfig"); + + if (this->VMix->VelVertMixSetup.Enabled || + this->VMix->TracerVertMixSetup.Enabled) { + + if (!this->EqState) { + ABORT_ERROR("Tendencies: Eos must be initialized when" + "vertical mixing tendencies are enabled"); + } + if (!this->VMix) { + ABORT_ERROR("Tendencies: VertMix must be initialized when" + "vertical mixing tendencies are enabled"); + } + } } //------------------------------------------------------------------------------ diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index acd486a0c448..8421efd0fe8f 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -455,41 +455,135 @@ void VertMix::applyVelVertMixImplicit( Eos *EosInstance = Eos::getInstance(); VertMix *VertMixInstance = VertMix::getInstance(); - if (!EosInstance) { - LOG_WARN("Eos has not been initialized. Skipping calculation of " - "VelVertMix tendency"); - } else if (!VertMixInstance) { - LOG_WARN("VertMix has not been initialized. Skipping calculation of " - "VelVertMix tendency"); - } else { - - // Obtain TimeStep - const auto *DefTimeStepper = TimeStepper::getDefault(); - const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); - R8 DT; - TimeStep.get(DT, TimeUnits::Seconds); - - const auto &SpecVol = EosInstance->SpecVol; - const auto &VertVisc = VertMixInstance->VertVisc; - - const int NVertLayers = VCoord->NVertLayers; - auto LConfig = - TriDiagSolver::makeLaunchConfig(Mesh->NEdgesAll, NVertLayers); + // Obtain TimeStep + const auto *DefTimeStepper = TimeStepper::getDefault(); + const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); + R8 DT; + TimeStep.get(DT, TimeUnits::Seconds); + const auto &SpecVol = EosInstance->SpecVol; + const auto &VertVisc = VertMixInstance->VertVisc; + + const int NVertLayers = VCoord->NVertLayers; + auto LConfig = + TriDiagSolver::makeLaunchConfig(Mesh->NEdgesAll, NVertLayers); + + parallelForOuter( + LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { + const int IStart = Team.league_rank() * VecLength; + const int ILen = + Kokkos::max(0, Kokkos::min(VecLength, LocNEdgesAll - IStart)); + + TriDiagDiffScratch Scratch(Team, NVertLayers); + + // Construct a tri-diag diffusion matrix and RHS + parallelForInner(Team, NVertLayers, [=](int K) { + for (int IVec = 0; IVec < VecLength; ++IVec) { + const int IEdge = IStart + IVec; + + if (IEdge >= LocNEdgesAll) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + + if (K < KMin || K > KMax) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + Real G, H, X; + LocVelVertMixSetup(IEdge, K, KMin, KMax, DT, SpecVol, + PseudoThickCell, VertVisc, NormalVelEdge, + G, H, X); + + Scratch.G(K, IVec) = G; + Scratch.H(K, IVec) = H; + Scratch.X(K, IVec) = X; + } + }); + + // Solve the tri-diag diffusion system + Team.team_barrier(); + TriDiagDiffSolver::solve(Team, Scratch); + Team.team_barrier(); + + // Store the solution vector X + parallelForInner(Team, NVertLayers, [=](int K) { + for (int IVec = 0; IVec < ILen; ++IVec) { + const int IEdge = IStart + IVec; + + if (K >= MinLayerEdgeBot(IEdge) && + K <= MaxLayerEdgeTop(IEdge)) { + NormalVelEdge(IEdge, K) = Scratch.X(K, IVec); + } + } + }); + }); + + Pacer::stop("Tend:velocityVertMix", 1); + } +} // applyVelVertMixImplicit + +// Apply implicit tracer vertical mixing +void VertMix::applyTracerVertMixImplicit( + OceanState *State, ///< [in] State variables + const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + Array3DReal &TracerArray, ///< [in] Tracer array + int NTracers, ///< [in] Number of tracers + int ThickTimeLevel, ///< [in] Time level + int VelTimeLevel ///< [in] Time level +) { + + OMEGA_SCOPE(LocNCellsAll, Mesh->NCellsAll); + OMEGA_SCOPE(LocTracerVertMixSetup, TracerVertMixSetup); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + + const Array2DReal &PseudoThickCell = State->PseudoThickness[ThickTimeLevel]; + + if (LocTracerVertMixSetup.Enabled) { + Pacer::start("Tend:tracerVertMix", 1); + + Eos *EosInstance = Eos::getInstance(); + VertMix *VertMixInstance = VertMix::getInstance(); + + // Obtain TimeStep + const auto *DefTimeStepper = TimeStepper::getDefault(); + const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); + R8 DT; + TimeStep.get(DT, TimeUnits::Seconds); + + const auto &SpecVol = EosInstance->SpecVol; + const auto &VertDiff = VertMixInstance->VertDiff; + + const int NVertLayers = VCoord->NVertLayers; + auto LConfig = + TriDiagSolver::makeLaunchConfig(Mesh->NCellsAll, NVertLayers); + + for (int L = 0; L < NTracers; ++L) { parallelForOuter( LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { const int IStart = Team.league_rank() * VecLength; const int ILen = Kokkos::max( - 0, Kokkos::min(VecLength, LocNEdgesAll - IStart)); + 0, Kokkos::min(VecLength, LocNCellsAll - IStart)); TriDiagDiffScratch Scratch(Team, NVertLayers); // Construct a tri-diag diffusion matrix and RHS parallelForInner(Team, NVertLayers, [=](int K) { for (int IVec = 0; IVec < VecLength; ++IVec) { - const int IEdge = IStart + IVec; + const int ICell = IStart + IVec; - if (IEdge >= LocNEdgesAll) { + if (ICell >= LocNCellsAll) { // Fill values Scratch.G(K, IVec) = 0._Real; Scratch.H(K, IVec) = 1._Real; @@ -497,8 +591,8 @@ void VertMix::applyVelVertMixImplicit( continue; } - const int KMin = MinLayerEdgeBot(IEdge); - const int KMax = MaxLayerEdgeTop(IEdge); + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); if (K < KMin || K > KMax) { // Fill values @@ -509,10 +603,9 @@ void VertMix::applyVelVertMixImplicit( } Real G, H, X; - LocVelVertMixSetup(IEdge, K, KMin, KMax, DT, SpecVol, - PseudoThickCell, VertVisc, - NormalVelEdge, G, H, X); - + LocTracerVertMixSetup(L, ICell, K, KMin, KMax, DT, + SpecVol, PseudoThickCell, VertDiff, + TracerArray, G, H, X); Scratch.G(K, IVec) = G; Scratch.H(K, IVec) = H; Scratch.X(K, IVec) = X; @@ -527,127 +620,18 @@ void VertMix::applyVelVertMixImplicit( // Store the solution vector X parallelForInner(Team, NVertLayers, [=](int K) { for (int IVec = 0; IVec < ILen; ++IVec) { - const int IEdge = IStart + IVec; + const int ICell = IStart + IVec; - if (K >= MinLayerEdgeBot(IEdge) && - K <= MaxLayerEdgeTop(IEdge)) { - NormalVelEdge(IEdge, K) = Scratch.X(K, IVec); + if (K >= MinLayerCell(ICell) && + K <= MaxLayerCell(ICell)) { + TracerArray(L, ICell, K) = Scratch.X(K, IVec); } } }); }); - } - Pacer::stop("Tend:velocityVertMix", 1); - } -} // applyVelVertMixImplicit -// Apply implicit tracer vertical mixing -void VertMix::applyTracerVertMixImplicit( - OceanState *State, ///< [in] State variables - const AuxiliaryState *AuxState, ///< [in] Auxilary state variables - Array3DReal &TracerArray, ///< [in] Tracer array - int NTracers, ///< [in] Number of tracers - int ThickTimeLevel, ///< [in] Time level - int VelTimeLevel ///< [in] Time level -) { - - OMEGA_SCOPE(LocNCellsAll, Mesh->NCellsAll); - OMEGA_SCOPE(LocTracerVertMixSetup, TracerVertMixSetup); - OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); - OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - - const Array2DReal &PseudoThickCell = State->PseudoThickness[ThickTimeLevel]; + } // for L - if (LocTracerVertMixSetup.Enabled) { - Pacer::start("Tend:tracerVertMix", 1); - - Eos *EosInstance = Eos::getInstance(); - VertMix *VertMixInstance = VertMix::getInstance(); - - if (!EosInstance) { - LOG_WARN("Eos has not been initialized. Skipping calculation of " - "PresGradZ tendency"); - } else if (!VertMixInstance) { - LOG_WARN("VertMix has not been initialized. Skipping calculation of " - "VelVertMix tendency"); - } else { - - // Obtain TimeStep - const auto *DefTimeStepper = TimeStepper::getDefault(); - const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); - R8 DT; - TimeStep.get(DT, TimeUnits::Seconds); - - const auto &SpecVol = EosInstance->SpecVol; - const auto &VertDiff = VertMixInstance->VertDiff; - - const int NVertLayers = VCoord->NVertLayers; - auto LConfig = - TriDiagSolver::makeLaunchConfig(Mesh->NCellsAll, NVertLayers); - - for (int L = 0; L < NTracers; ++L) { - parallelForOuter( - LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { - const int IStart = Team.league_rank() * VecLength; - const int ILen = Kokkos::max( - 0, Kokkos::min(VecLength, LocNCellsAll - IStart)); - - TriDiagDiffScratch Scratch(Team, NVertLayers); - - // Construct a tri-diag diffusion matrix and RHS - parallelForInner(Team, NVertLayers, [=](int K) { - for (int IVec = 0; IVec < VecLength; ++IVec) { - const int ICell = IStart + IVec; - - if (ICell >= LocNCellsAll) { - // Fill values - Scratch.G(K, IVec) = 0._Real; - Scratch.H(K, IVec) = 1._Real; - Scratch.X(K, IVec) = 0._Real; - continue; - } - - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - - if (K < KMin || K > KMax) { - // Fill values - Scratch.G(K, IVec) = 0._Real; - Scratch.H(K, IVec) = 1._Real; - Scratch.X(K, IVec) = 0._Real; - continue; - } - - Real G, H, X; - LocTracerVertMixSetup(L, ICell, K, KMin, KMax, DT, - SpecVol, PseudoThickCell, - VertDiff, TracerArray, G, H, X); - Scratch.G(K, IVec) = G; - Scratch.H(K, IVec) = H; - Scratch.X(K, IVec) = X; - } - }); - - // Solve the tri-diag diffusion system - Team.team_barrier(); - TriDiagDiffSolver::solve(Team, Scratch); - Team.team_barrier(); - - // Store the solution vector X - parallelForInner(Team, NVertLayers, [=](int K) { - for (int IVec = 0; IVec < ILen; ++IVec) { - const int ICell = IStart + IVec; - - if (K >= MinLayerCell(ICell) && - K <= MaxLayerCell(ICell)) { - TracerArray(L, ICell, K) = Scratch.X(K, IVec); - } - } - }); - }); - - } // for L - } Pacer::stop("Tend:tracerVertMix", 1); } From 31bac0ccf1f02babd8083780dc8edd468145e7da Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 1 Jul 2026 09:33:21 -0400 Subject: [PATCH 18/19] =?UTF-8?q?Address=20reviewers=E2=80=99=20comments?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Improve the inline documentation in VertMixSetup --- components/omega/src/ocn/VertMix.h | 78 ++++++++++++++++++++++++++---- 1 file changed, 68 insertions(+), 10 deletions(-) diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index 008476125c99..7e8312fa06f8 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -234,12 +234,8 @@ class VelVertMixSetupOnEdge { bool Enabled; Real LocRhoSw; - /// constructor declaration VelVertMixSetupOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); - /// The functor takes edge index, vertical chunk index, and arrays for - /// layer specific volume, layer thickness on edge, - /// interface pressure, and outputs tendency array KOKKOS_FUNCTION void operator()(I4 IEdge, I4 K, I4 KMin, I4 KMax, Real DT, const Array2DReal &SpecVol, const Array2DReal &PseudoThickCell, @@ -247,20 +243,51 @@ class VelVertMixSetupOnEdge { const Array2DReal &NormalVelEdge, Real &G, Real &H, Real &X) const { - const I4 JCell0 = CellsOnEdge(IEdge, 0); - const I4 JCell1 = CellsOnEdge(IEdge, 1); + // Implicit vertical mixing is solved using the diffusion-form + // tridiagonal solver, `TriDiagDiffSolver`, provided by `TriDiagSolvers`. + // + // * Governing tridiagonal system: + // + // -G_{k-1} Y_{k-1} + // + (G_{k-1} + G_k + H_k) Y_k + // - G_k Y_{k+1} + // = X_k + // + // G: coefficient contributing to both diagonal and off-diagonal entries + // H: coefficient contributing only to the diagonal entry + // X: RHS vector + // + // * Matrix structure: + // + // [ D G 0 ... 0 ][ ] [ ] : k = 0 + // [ G D G ... 0 ][ ] [ ] + // [ 0 G D G ... 0 ][ ] [ ] + // [ ... ][Y] = [X] + // [ 0 ... G D G 0 ][ ] [ ] + // [ 0 ... 0 G D G ][ ] [ ] + // [ 0 ... 0 0 G D ][ ] [ ] : k = NVertLayers - 1 + // + // where D = G_{k-1} + G_k + H_k + // + // G_k = [DT * VertVisc_k^top / (Rho_0 * SpecVol_k^top)] + // / PseudoThick_k^top + // H_k = PseudoThick_k + // X_k = PseudoThick_k * NormVel_k // Fill values G = 0.0_Real; H = 1.0_Real; X = 0.0_Real; + const I4 JCell0 = CellsOnEdge(IEdge, 0); + const I4 JCell1 = CellsOnEdge(IEdge, 1); + const Real PseudoThickEdgeK = 0.5_Real * (PseudoThickCell(JCell0, K) + PseudoThickCell(JCell1, K)); - // Row-scaled conservative form. - // Unknown is u^{n+1}. H = PseudoThickEdgeK; + + // Unknown is NormVel^{n+1}. X = PseudoThickEdgeK * NormalVelEdge(IEdge, K); if (K < KMax) { @@ -311,6 +338,37 @@ class TracerVertMixSetupOnCell { const Array3DReal &TracersOnCell, Real &G, Real &H, Real &X) const { + // Implicit vertical mixing is solved using the diffusion-form + // tridiagonal solver, `TriDiagDiffSolver`, provided by `TriDiagSolvers`. + // + // * Governing tridiagonal system: + // + // -G_{k-1} Y_{k-1} + // + (G_{k-1} + G_k + H_k) Y_k + // - G_k Y_{k+1} + // = X_k + // + // G: coefficient contributing to both diagonal and off-diagonal entries + // H: coefficient contributing only to the diagonal entry + // X: RHS vector + // + // * Matrix structure: + // + // [ D G 0 ... 0 ][ ] [ ] : k = 0 + // [ G D G ... 0 ][ ] [ ] + // [ 0 G D G ... 0 ][ ] [ ] + // [ ... ][Y] = [X] + // [ 0 ... G D G 0 ][ ] [ ] + // [ 0 ... 0 G D G ][ ] [ ] + // [ 0 ... 0 0 G D ][ ] [ ] : k = NVertLayers - 1 + // + // where D = G_{k-1} + G_k + H_k + // + // G_k = [DT * VertDiff_k^top / (Rho_0 * SpecVol_k^top)] + // / PseudoThick_k^top + // H_k = PseudoThick_k + // X_k = PseudoThick_k * Phi_k + // Fill values G = 0.0_Real; H = 1.0_Real; @@ -318,9 +376,9 @@ class TracerVertMixSetupOnCell { const Real PseudoThickCellK = PseudoThickCell(ICell, K); - // Row-scaled conservative form. - // Unknown is phi^{n+1}. H = PseudoThickCellK; + + // Unknown is Phi^{n+1}. X = PseudoThickCellK * TracersOnCell(L, ICell, K); if (K < KMax) { From 102dadaee4d94946aae07f14997606a3747d5f55 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 1 Jul 2026 13:27:23 -0400 Subject: [PATCH 19/19] Resolve CTest issues after rebasing --- components/omega/test/infra/IOStreamTest.cpp | 4 ++++ components/omega/test/ocn/StateTest.cpp | 6 +++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index 9fc87452298d..b3126ad21a64 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -15,6 +15,7 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Eos.h" #include "Error.h" #include "Field.h" #include "Forcing.h" @@ -109,6 +110,9 @@ void initIOStreamTest(Clock *&ModelClock // Model clock PressureGrad::init(); + // Initialize Eos + Eos::init(); + // Initialize VertMix VertMix::init(); diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp index fa3fa8f5fb1c..09caff2e6d73 100644 --- a/components/omega/test/ocn/StateTest.cpp +++ b/components/omega/test/ocn/StateTest.cpp @@ -13,6 +13,7 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Eos.h" #include "Error.h" #include "Field.h" #include "Halo.h" @@ -96,7 +97,10 @@ void initStateTest() { // Initialize pressure gradient PressureGrad::init(); - // Create tendencies + // Initialize equation of state + Eos::init(); + + // Initialize vertical mixing VertMix::init(); // Create tendencies