From a10434ae9aa555ee502668835a1ba4f9cce0097c Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Wed, 4 Feb 2026 20:19:43 -0800 Subject: [PATCH 1/3] Fixes ssh calculation Moves SSH calculation to the VertCoord class and computes it as a single layer instead of a multi level field. --- components/omega/src/ocn/AuxiliaryState.cpp | 24 -------------- components/omega/src/ocn/Tendencies.cpp | 3 +- components/omega/src/ocn/TendencyTerms.h | 4 +-- components/omega/src/ocn/VertCoord.cpp | 31 +++++++++++++++++++ components/omega/src/ocn/VertCoord.h | 3 ++ .../auxiliaryVars/LayerThicknessAuxVars.cpp | 20 ------------ .../ocn/auxiliaryVars/LayerThicknessAuxVars.h | 7 ----- .../omega/test/ocn/TendencyTermsTest.cpp | 3 +- 8 files changed, 39 insertions(+), 56 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 4de2a47884f0..ff6037ed26b2 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -189,30 +189,6 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, }); Pacer::stop("AuxState:cellAuxState2", 2); - Pacer::start("AuxState:cellAuxState3", 2); - parallelForOuter( - "cellAuxState3", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(int 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) { - LocLayerThicknessAux.computeVarsOnCells( - ICell, KChunk, LayerThickCell, NormalVelEdge, - TimeStepSeconds); - }); - }); - Pacer::stop("AuxState:cellAuxState3", 2); - - Pacer::start("AuxState:computeVerticalVelocity", 2); - - const auto &FluxLayerThickEdge = LayerThicknessAux.FluxLayerThickEdge; - VAdv->computeVerticalVelocity(NormalVelEdge, FluxLayerThickEdge); - - Pacer::stop("AuxState:computeVerticalVelocity", 2); - Pacer::stop("AuxState:computeMomAux", 1); } diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index e07640507ec3..9d128a9d5f56 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -375,6 +375,7 @@ void Tendencies::computeVelocityTendenciesOnly( OMEGA_SCOPE(LocBottomDrag, BottomDrag); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + OMEGA_SCOPE(LocSshCell, VCoord->SshCell); Pacer::start("Tend:computeVelocityTendenciesOnly", 1); @@ -433,7 +434,7 @@ void Tendencies::computeVelocityTendenciesOnly( } // Compute sea surface height gradient - const Array2DReal &SSHCell = AuxState->LayerThicknessAux.SshCell; + const Array1DReal &SSHCell = LocSshCell; if (LocSSHGrad.Enabled) { Pacer::start("Tend:SSHGrad", 2); parallelForOuter( diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 5ccb1ec79ed3..3f4793e7af5c 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -173,7 +173,7 @@ class SSHGradOnEdge { /// The functor takes edge index, vertical chunk index, and array of /// layer thickness/SSH, outputs tendency array KOKKOS_FUNCTION void operator()(const Array2DReal &Tend, I4 IEdge, I4 KChunk, - const Array2DReal &SshCell) const { + const Array1DReal &SshCell) const { const I4 KStart = chunkStart(KChunk, MinLayerEdgeBot(IEdge)); const I4 KLen = chunkLength(KChunk, KStart, MaxLayerEdgeTop(IEdge)); @@ -184,7 +184,7 @@ class SSHGradOnEdge { for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; Tend(IEdge, K) -= EdgeMask(IEdge, K) * Gravity * - (SshCell(ICell1, K) - SshCell(ICell0, K)) * + (SshCell(ICell1) - SshCell(ICell0)) * InvDcEdge; } } diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 027006d1164c..3df0b8003b8d 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -163,6 +163,8 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord PressureMid = Array2DReal("PressureMid", NCellsSize, NVertLayers); ZInterface = Array2DReal("ZInterface", NCellsSize, NVertLayersP1); ZMid = Array2DReal("ZMid", NCellsSize, NVertLayers); + SshCell = Array1DReal("SshCell", NCellsSize); + GeopotentialMid = Array2DReal("GeopotentialMid", NCellsSize, NVertLayers); LayerThicknessTarget = Array2DReal("LayerThicknessTarget", NCellsSize, NVertLayers); @@ -174,6 +176,8 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord PressureMidH = createHostMirrorCopy(PressureMid); ZInterfaceH = createHostMirrorCopy(ZInterface); ZMidH = createHostMirrorCopy(ZMid); + SshCellH = createHostMirrorCopy(SshCellH); + GeopotentialMidH = createHostMirrorCopy(GeopotentialMid); LayerThicknessTargetH = createHostMirrorCopy(LayerThicknessTarget); RefLayerThicknessH = createHostMirrorCopy(RefLayerThickness); @@ -236,6 +240,8 @@ void VertCoord::defineFields() { PressMidFldName = "PressureMid"; ZInterfFldName = "ZInterface"; ZMidFldName = "ZMid"; + SshFldName = "SshCell"; + GeopotFldName = "GeopotentialMid"; LyrThickTargetFldName = "LayerThicknessTarget"; @@ -249,6 +255,7 @@ void VertCoord::defineFields() { ZMidFldName.append(Name); GeopotFldName.append(Name); LyrThickTargetFldName.append(Name); + SshFldName.append(Name); } // Create fields for VertCoord variables @@ -295,6 +302,19 @@ void VertCoord::defineFields() { DimNames // dimension names ); + auto SshField = Field::create( + SshFldName, // field name + "sea surface height at cell center", // long Name or description + "m", // units + "sea_surface_height", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValueReal, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dimension names + ); + + NDims = 2; DimNames.resize(NDims); DimNames[1] = "NVertLayersP1"; @@ -402,6 +422,7 @@ void VertCoord::defineFields() { VCoordGroup->addField(ZMidFldName); VCoordGroup->addField(GeopotFldName); VCoordGroup->addField(LyrThickTargetFldName); + VCoordGroup->addField(SshFldName); // Associate Field with data PressureInterfaceField->attachData(PressureInterface); @@ -410,6 +431,7 @@ void VertCoord::defineFields() { ZMidField->attachData(ZMid); GeopotentialMidField->attachData(GeopotentialMid); LayerThicknessTargetField->attachData(LayerThicknessTarget); + SshField->attachData(SshCell); } // end defineFields @@ -431,6 +453,7 @@ VertCoord::~VertCoord() { Field::destroy(ZMidFldName); Field::destroy(GeopotFldName); Field::destroy(LyrThickTargetFldName); + Field::destroy(SshFldName); FieldGroup::destroy(GroupName); } @@ -872,6 +895,7 @@ void VertCoord::computeZHeight( OMEGA_SCOPE(LocZInterf, ZInterface); OMEGA_SCOPE(LocZMid, ZMid); OMEGA_SCOPE(LocBotDepth, BottomDepth); + OMEGA_SCOPE(LocSshCell, SshCell); parallelForOuter( "computeZHeight", {NCellsAll}, @@ -891,6 +915,9 @@ void VertCoord::computeZHeight( LocZInterf(ICell, KLyr) = -LocBotDepth(ICell) + Accum; LocZMid(ICell, KLyr) = -LocBotDepth(ICell) + Accum - 0.5 * DZ; + if (KLyr == 0) { + LocSshCell(ICell) = LocZInterf(ICell,KLyr); + } } }); }); @@ -995,6 +1022,8 @@ void VertCoord::copyToHost() { deepCopy(PressureMidH, PressureMid); deepCopy(ZInterfaceH, ZInterface); deepCopy(ZMidH, ZMid); + deepCopy(SshCellH, SshCell); + deepCopy(GeopotentialMidH, GeopotentialMid); deepCopy(LayerThicknessTargetH, LayerThicknessTarget); deepCopy(RefLayerThicknessH, RefLayerThickness); @@ -1008,6 +1037,8 @@ void VertCoord::copyToDevice() { deepCopy(PressureMid, PressureMidH); deepCopy(ZInterface, ZInterfaceH); deepCopy(ZMid, ZMidH); + deepCopy(SshCell, SshCellH); + deepCopy(GeopotentialMid, GeopotentialMidH); deepCopy(LayerThicknessTarget, LayerThicknessTargetH); deepCopy(RefLayerThickness, RefLayerThicknessH); diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h index d39d3676b225..7de96b0da559 100644 --- a/components/omega/src/ocn/VertCoord.h +++ b/components/omega/src/ocn/VertCoord.h @@ -107,6 +107,7 @@ class VertCoord { Array2DReal ZMid; Array2DReal GeopotentialMid; Array2DReal LayerThicknessTarget; + Array1DReal SshCell; HostArray2DReal PressureInterfaceH; HostArray2DReal PressureMidH; @@ -114,6 +115,7 @@ class VertCoord { HostArray2DReal ZMidH; HostArray2DReal GeopotentialMidH; HostArray2DReal LayerThicknessTargetH; + HostArray1DReal SshCellH; // Vertical loop bounds Array1DI4 MinLayerCell; @@ -181,6 +183,7 @@ class VertCoord { std::string ZMidFldName; ///< Field name for midpoint Z height std::string GeopotFldName; ///< Field name for geopotential std::string LyrThickTargetFldName; ///< Field name for target thickness + std::string SshFldName; ///< Field name for sea surface height // methods diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp index 8a69aa7c8269..5d24bbf84329 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp @@ -12,14 +12,11 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix, Mesh->NEdgesSize, VCoord->NVertLayers), MeanLayerThickEdge("MeanLayerThickEdge" + AuxStateSuffix, Mesh->NEdgesSize, VCoord->NVertLayers), - SshCell("SshCell" + AuxStateSuffix, Mesh->NCellsSize, - VCoord->NVertLayers), ProvThickness("ProvThickness" + AuxStateSuffix, Mesh->NCellsSize, VCoord->NVertLayers), AreaCell(Mesh->AreaCell), DvEdge(Mesh->DvEdge), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), CellsOnEdge(Mesh->CellsOnEdge), - BottomDepth(VCoord->BottomDepth), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} @@ -69,20 +66,6 @@ void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName, DimNames // dimension names ); - // Sea surface height - DimNames[0] = "NCells" + DimSuffix; - auto SshCellField = Field::create( - SshCell.label(), // field name - "sea surface height at cell center", // long Name or description - "m", // units - "sea_surface_height", // CF standard Name - 0, // min valid value - std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries - NDims, // number of dimensions - DimNames // dimension names - ); - // Provisional Thickness auto ProvThicknessField = Field::create( ProvThickness.label(), // field name @@ -99,20 +82,17 @@ void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName, // Add fields to Aux field group FieldGroup::addFieldToGroup(FluxLayerThickEdge.label(), AuxGroupName); FieldGroup::addFieldToGroup(MeanLayerThickEdge.label(), AuxGroupName); - FieldGroup::addFieldToGroup(SshCell.label(), AuxGroupName); FieldGroup::addFieldToGroup(ProvThickness.label(), AuxGroupName); // Attach field data FluxLayerThickEdgeField->attachData(FluxLayerThickEdge); MeanLayerThickEdgeField->attachData(MeanLayerThickEdge); - SshCellField->attachData(SshCell); ProvThicknessField->attachData(ProvThickness); } void LayerThicknessAuxVars::unregisterFields() const { Field::destroy(FluxLayerThickEdge.label()); Field::destroy(MeanLayerThickEdge.label()); - Field::destroy(SshCell.label()); Field::destroy(ProvThickness.label()); } diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h index 68bb7da2296c..5c640eee40ea 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h @@ -16,7 +16,6 @@ class LayerThicknessAuxVars { public: Array2DReal FluxLayerThickEdge; Array2DReal MeanLayerThickEdge; - Array2DReal SshCell; Array2DReal ProvThickness; FluxThickEdgeOption FluxThickEdgeChoice; @@ -73,11 +72,6 @@ class LayerThicknessAuxVars { const int KStart = chunkStart(KChunk, MinLayerCell(ICell)); const int KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); - for (int KVec = 0; KVec < KLen; ++KVec) { - const int K = KStart + KVec; - SshCell(ICell, K) = LayerThickCell(ICell, K) - BottomDepth(ICell); - } - Real TmpProv[VecLength] = {0.}; Real DtInvAreaCell = Dt / AreaCell(ICell); @@ -109,7 +103,6 @@ class LayerThicknessAuxVars { Array2DI4 EdgesOnCell; Array2DReal EdgeSignOnCell; Array2DI4 CellsOnEdge; - Array1DReal BottomDepth; Array1DI4 MinLayerEdgeBot; Array1DI4 MaxLayerEdgeTop; Array1DI4 MinLayerCell; diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index c34fe574d10c..debf0fe5042d 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -533,8 +533,7 @@ int testSSHGrad(int NVertLayers, Real RTol) { }, ExactSSHGrad, EdgeComponent::Normal, Geom, Mesh, ExchangeHalos::No); - // Set input array - Array2DReal SSHCell("SSHCell", Mesh->NCellsSize, NVertLayers); + Array1DReal SSHCell("SSHCell", Mesh->NCellsSize); Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalar(X, Y); }, SSHCell, From 0141956ce62f66d5cbc7b2a397a8d4cbf2b2905e Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 5 Feb 2026 07:46:38 -0800 Subject: [PATCH 2/3] Fixes kmin and linting --- components/omega/src/ocn/VertCoord.cpp | 45 +++++++++++++------------- 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp index 3df0b8003b8d..28a57cb1cd67 100644 --- a/components/omega/src/ocn/VertCoord.cpp +++ b/components/omega/src/ocn/VertCoord.cpp @@ -160,10 +160,10 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord BottomDepth = Array1DReal("BottomDepth", NCellsSize); PressureInterface = Array2DReal("PressureInterface", NCellsSize, NVertLayersP1); - PressureMid = Array2DReal("PressureMid", NCellsSize, NVertLayers); - ZInterface = Array2DReal("ZInterface", NCellsSize, NVertLayersP1); - ZMid = Array2DReal("ZMid", NCellsSize, NVertLayers); - SshCell = Array1DReal("SshCell", NCellsSize); + PressureMid = Array2DReal("PressureMid", NCellsSize, NVertLayers); + ZInterface = Array2DReal("ZInterface", NCellsSize, NVertLayersP1); + ZMid = Array2DReal("ZMid", NCellsSize, NVertLayers); + SshCell = Array1DReal("SshCell", NCellsSize); GeopotentialMid = Array2DReal("GeopotentialMid", NCellsSize, NVertLayers); LayerThicknessTarget = @@ -172,11 +172,11 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord Array2DReal("RefLayerThickness", NCellsSize, NVertLayers); // Make host copies for device arrays not being read from file - PressureInterfaceH = createHostMirrorCopy(PressureInterface); - PressureMidH = createHostMirrorCopy(PressureMid); - ZInterfaceH = createHostMirrorCopy(ZInterface); - ZMidH = createHostMirrorCopy(ZMid); - SshCellH = createHostMirrorCopy(SshCellH); + PressureInterfaceH = createHostMirrorCopy(PressureInterface); + PressureMidH = createHostMirrorCopy(PressureMid); + ZInterfaceH = createHostMirrorCopy(ZInterface); + ZMidH = createHostMirrorCopy(ZMid); + SshCellH = createHostMirrorCopy(SshCellH); GeopotentialMidH = createHostMirrorCopy(GeopotentialMid); LayerThicknessTargetH = createHostMirrorCopy(LayerThicknessTarget); @@ -233,14 +233,14 @@ VertCoord *VertCoord::create( void VertCoord::defineFields() { // Set field names (append Name if not default) - MinLayerCellFldName = "MinLayerCell"; - MaxLayerCellFldName = "MaxLayerCell"; - BottomDepthFldName = "BottomDepth"; - PressInterfFldName = "PressureInterface"; - PressMidFldName = "PressureMid"; - ZInterfFldName = "ZInterface"; - ZMidFldName = "ZMid"; - SshFldName = "SshCell"; + MinLayerCellFldName = "MinLayerCell"; + MaxLayerCellFldName = "MaxLayerCell"; + BottomDepthFldName = "BottomDepth"; + PressInterfFldName = "PressureInterface"; + PressMidFldName = "PressureMid"; + ZInterfFldName = "ZInterface"; + ZMidFldName = "ZMid"; + SshFldName = "SshCell"; GeopotFldName = "GeopotentialMid"; LyrThickTargetFldName = "LayerThicknessTarget"; @@ -303,18 +303,17 @@ void VertCoord::defineFields() { ); auto SshField = Field::create( - SshFldName, // field name + SshFldName, // field name "sea surface height at cell center", // long Name or description "m", // units "sea_surface_height", // CF standard Name std::numeric_limits::min(), // min valid value std::numeric_limits::max(), // max valid value - FillValueReal, // scalar for undefined entries + FillValueReal, // scalar for undefined entries NDims, // number of dimensions - DimNames // dimension names + DimNames // dimension names ); - NDims = 2; DimNames.resize(NDims); DimNames[1] = "NVertLayersP1"; @@ -915,8 +914,8 @@ void VertCoord::computeZHeight( LocZInterf(ICell, KLyr) = -LocBotDepth(ICell) + Accum; LocZMid(ICell, KLyr) = -LocBotDepth(ICell) + Accum - 0.5 * DZ; - if (KLyr == 0) { - LocSshCell(ICell) = LocZInterf(ICell,KLyr); + if (KLyr == KMin) { + LocSshCell(ICell) = LocZInterf(ICell, KLyr); } } }); From 52001a1460312e78feaedc5d612c762fd4c96984 Mon Sep 17 00:00:00 2001 From: Luke Van Roekel Date: Thu, 5 Feb 2026 08:00:02 -0800 Subject: [PATCH 3/3] Linting fixes --- components/omega/src/ocn/Tendencies.cpp | 2 +- components/omega/src/ocn/TendencyTerms.h | 3 +-- components/omega/src/ocn/VertCoord.h | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 9d128a9d5f56..a438c1982560 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -434,7 +434,7 @@ void Tendencies::computeVelocityTendenciesOnly( } // Compute sea surface height gradient - const Array1DReal &SSHCell = LocSshCell; + const Array1DReal &SSHCell = LocSshCell; if (LocSSHGrad.Enabled) { Pacer::start("Tend:SSHGrad", 2); parallelForOuter( diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 3f4793e7af5c..988ff11d7c7c 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -184,8 +184,7 @@ class SSHGradOnEdge { for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; Tend(IEdge, K) -= EdgeMask(IEdge, K) * Gravity * - (SshCell(ICell1) - SshCell(ICell0)) * - InvDcEdge; + (SshCell(ICell1) - SshCell(ICell0)) * InvDcEdge; } } diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h index 7de96b0da559..d75c87e35eb1 100644 --- a/components/omega/src/ocn/VertCoord.h +++ b/components/omega/src/ocn/VertCoord.h @@ -183,7 +183,7 @@ class VertCoord { std::string ZMidFldName; ///< Field name for midpoint Z height std::string GeopotFldName; ///< Field name for geopotential std::string LyrThickTargetFldName; ///< Field name for target thickness - std::string SshFldName; ///< Field name for sea surface height + std::string SshFldName; ///< Field name for sea surface height // methods