From 5af4160b7695310078e85b70aefb4fef90e9fbac Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Sat, 27 Dec 2025 12:43:53 -0800 Subject: [PATCH 01/42] VertAdv module definition; creation, destruction, retrieval methods --- components/omega/src/ocn/VertAdv.cpp | 207 +++++++++++++++++++++++++++ components/omega/src/ocn/VertAdv.h | 138 ++++++++++++++++++ 2 files changed, 345 insertions(+) create mode 100644 components/omega/src/ocn/VertAdv.cpp create mode 100644 components/omega/src/ocn/VertAdv.h diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp new file mode 100644 index 000000000000..4ae65d1d94e3 --- /dev/null +++ b/components/omega/src/ocn/VertAdv.cpp @@ -0,0 +1,207 @@ +//===-- ocn/VertAdv.cpp - vertical advection --------------------*- C++ -*-===// +// +//===----------------------------------------------------------------------===// + +#include "Tracers.h" +#include "VertAdv.h" + +namespace OMEGA { + +// create static class members +VertAdv *VertAdv::DefaultVertAdv = nullptr; +std::map> VertAdv::AllVertAdvs; + +//------------------------------------------------------------------------------ +// create the default VertAdv, requires prior initialization of default +// HorzMesh and VertCoord +void VertAdv::init() { + + auto Mesh = HorzMesh::getDefault(); + auto VCoord = VertCoord::getDefault(); + + Config *OmegaConfig = Config::getOmegaConfig(); + + VertAdv::DefaultVertAdv = create("Default", Mesh, VCoord, OmegaConfig); + +} // end init + +//------------------------------------------------------------------------------ +// constructor +VertAdv::VertAdv(const std::string &Name_, ///< [in] name for new VertAdv + const HorzMesh *InMesh, ///< [in] associated HorzMesh + const VertCoord *InVCoord, ///< [in] associated VertCoord + Config *Options ///< [in] configuration options +) { + + // Read options from Config object + readConfigOptions(Options); + + // Store name suffix + Name = Name_; + + // Store pointers to Mesh and VertCoord objects + Mesh = InMesh; + VCoord = InVCoord; + + // Store dimension sizes + NVertLayers = VCoord->NVertLayers; + NVertLayersP1 = VCoord->NVertLayersP1; + NCellsOwned = Mesh->NCellsOwned; + NCellsAll = Mesh->NCellsAll; + NCellsSize = Mesh->NCellsSize; + NEdgesOwned = Mesh->NEdgesOwned; + NEdgesAll = Mesh->NEdgesAll; + NEdgesSize = Mesh->NEdgesSize; + NTracers = Tracers::getNumTracers(); + + // Allocate member arrays + VerticalVelocity = + Array2DReal("VerticalVelocity", NCellsSize, NVertLayersP1); + TotalVerticalVelocity = + Array2DReal("TotalVerticalVelocity", NCellsSize, NVertLayersP1); + VertFlux = Array3DReal("VertFlux", NTracers, NCellsSize, NVertLayersP1); + + // Low-order flux array only needed for flux-corrected transport + if (VertAdvChoice == VertAdvOption::FCT) { + LowOrderVertFlux = + Array3DReal("LowOrderVertFlux", NTracers, NCellsSize, NVertLayersP1); + } + +} // end constructor + +//------------------------------------------------------------------------------ +// create a new VertAdv instance +VertAdv *VertAdv::create(const std::string &Name, ///< [in] name for new VertAdv + const HorzMesh *Mesh, ///< [in] associated HorzMesh + const VertCoord *VCoord, ///< [in] associated VertCoord + Config *Options ///< [in] configuration options +) { + // Check to see if a VertAdv of the same name already exists and, if so, + // exit with an error + if (AllVertAdvs.find(Name) != AllVertAdvs.end()) { + LOG_ERROR("Attempted to create a VertAdv with name {} but a VertAdv " + "of that name already exists", + Name); + return nullptr; + } + + // create a new VertAdv on the heap and put it in a map of unique_ptrs, + // which will manage its lifetime + auto *NewVertAdv = new VertAdv(Name, Mesh, VCoord, Options); + AllVertAdvs.emplace(Name, NewVertAdv); + + return NewVertAdv; + +} // end create + +//------------------------------------------------------------------------------ +// destructor +VertAdv::~VertAdv() {} // end destructor + +//------------------------------------------------------------------------------ +// Removes all VertAdvs to clean up before exit +void VertAdv::clear() { + + AllVertAdvs.clear(); // removes all VertAdvs from the list and in the + // process, calls the destructors for each + +} // end clear + +//------------------------------------------------------------------------------ +// Removes a VertAdv from map by name +void VertAdv::erase(std::string Name) { + AllVertAdvs.erase(Name); // removes the VertAdv from the list and in the + // process, calls the destructor +} // end erase + +//------------------------------------------------------------------------------ +// Get default VertAdv +VertAdv *VertAdv::getDefault() { return VertAdv::DefaultVertAdv; } + +//------------------------------------------------------------------------------ +// Get VertAdv by name +VertAdv *VertAdv::get(const std::string Name ///< [in] Name of VertAdv +) { + + // look for an instance of this name + auto it = AllVertAdvs.find(Name); + + // if found, return the VertAdv pointer + if (it != AllVertAdvs.end()) { + return it->second.get(); + + // otherwise print error and return null pointer + } else { + LOG_ERROR("VertAdv::get: Attempt to retrieve non-existant VertAdv:"); + LOG_ERROR("{} has not been defined or has been removed", Name); + return nullptr; + } + +} // end get VertAdv + +// Read and set config options +void VertAdv::readConfigOptions(Config *OmegaConfig) { + + Error Err; // Error code + + Config TendConfig("Tendencies"); + Err += OmegaConfig->get(TendConfig); + CHECK_ERROR_ABORT(Err, "VertAdv: Tendencies group not in Config"); + + Err += TendConfig.get("ThicknessVertAdvTendencyEnable", + this->ThickVertAdvEnabled); + CHECK_ERROR_ABORT( + Err, "VertAdv: ThicknessVertAdvTendencyEnable not found in TendConfig"); + + Err += + TendConfig.get("VelocityVertAdvTendencyEnable", this->VelVertAdvEnabled); + CHECK_ERROR_ABORT( + Err, "VertAdv: ThicknessVertAdvTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("TracerVertAdvTendencyEnable", + this->TracerVertAdvEnabled); + CHECK_ERROR_ABORT( + Err, "VertAdv: ThicknessVertAdvTendencyEnable not found in TendConfig"); + + Config AdvectConfig("Advection"); + Err += OmegaConfig->get(AdvectConfig); + CHECK_ERROR_ABORT(Err, "VertAdv: Advection group not in Config"); + + bool FluxLimiterOn; + Err += AdvectConfig.get("VerticalTracerFluxLimiterEnabled", FluxLimiterOn); + CHECK_ERROR_ABORT( + Err, + "VertAdv: VerticalTracerFluxLimiterEnabled not found in AdvectConfig"); + if (FluxLimiterOn) { + VertAdvChoice = VertAdvOption::FCT; + } else { + VertAdvChoice = VertAdvOption::Standard; + } + + I4 VertFluxOrder; + Err += AdvectConfig.get("VerticalTracerFluxOrder", VertFluxOrder); + CHECK_ERROR_ABORT( + Err, "VertAdv: VerticalTracerFluxOrder not found in AdvectConfig"); + + switch (VertFluxOrder) { + case (2): + VertFluxChoice = VertFluxOption::Second; + break; + case (3): + VertFluxChoice = VertFluxOption::Third; + break; + case (4): + VertFluxChoice = VertFluxOption::Fourth; + break; + default: + ABORT_ERROR("VertAdv: Invalid option for VerticalTracerFluxOrder found " + "in AdvectConfig. Must be 2, 3, or 4"); + } + + Err += AdvectConfig.get("Coef3rdOrder", Coef3rdOrder); + CHECK_ERROR_ABORT(Err, "VertAdv: Coef3rdOrder not found in AdvectConfig"); + +} // end readConfigOptions + +} // end namespace OMEGA +//===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h new file mode 100644 index 000000000000..04f780fca123 --- /dev/null +++ b/components/omega/src/ocn/VertAdv.h @@ -0,0 +1,138 @@ +#ifndef OMEGA_VERTADV_H +#define OMEGA_VERTADV_H +//===-- ocn/VertAdv.h - vertical advection ----------------------*- C++ -*-===// +// +//===----------------------------------------------------------------------===// + +#include "AuxiliaryState.h" +#include "Config.h" +#include "DataTypes.h" +#include "Error.h" +#include "HorzMesh.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "TimeMgr.h" +#include "VertCoord.h" + +#include +#include +#include + +namespace OMEGA { + +// Enum for flux types +enum class VertFluxOption { Second, Third, Fourth }; +// Enum for advection scheme type +enum class VertAdvOption { Standard, FCT }; + +class VertAdv { + + public: + // Logicals to enable/disable advection of specific fields + bool ThickVertAdvEnabled = false; + bool VelVertAdvEnabled = false; + bool TracerVertAdvEnabled = false; + + // enums storing options chosen by user + VertFluxOption VertFluxChoice; + VertAdvOption VertAdvChoice; + + // Mesh dimensions + I4 NVertLayers; + I4 NVertLayersP1; + I4 NCellsOwned; + I4 NCellsAll; + I4 NCellsSize; + I4 NEdgesOwned; + I4 NEdgesAll; + I4 NEdgesSize; + + // Number of tracers + I4 NTracers; + + // Coefficient for blending 3rd-order and 4th-order reconstruction of tracers + Real Coef3rdOrder; + + // VertAdv instance name + std::string Name; + + Array2DReal VerticalVelocity; // pseudovelocity through top of cell + Array2DReal TotalVerticalVelocity; // transport velocity through top of Cell + Array3DReal VertFlux; // fluxes at vertical interfaces + Array3DReal LowOrderVertFlux; // low-order fluxes for FCT + + // public methods + + // Initialize the default VertAdv instance + static void init(); + + /// Creates a new vertical advection object by calling the constructor and + /// puts it in the AllVertAdvs map. + static VertAdv *create(const std::string &Name, // [in] name for new VertAdv + const HorzMesh *Mesh, // [in] associated HorzMesh + const VertCoord *VCoord, // [in] associated VertCoord + Config *Options // [in] configuration options + ); + + /// Destructor - deallocates all memory and deletes a VertAdv + ~VertAdv(); + + /// Deallocates arrays + static void clear(); + + /// Remove a VertAdv by name + static void erase(std::string InName); + + /// Retrieve the default VertAdv + static VertAdv *getDefault(); + + /// Retreive a VertAdv by name + static VertAdv *get(std::string name); + + /// Read and set config options + void readConfigOptions(Config *Options); + + private: + // Vertical loop bounds from VertCoord + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; + + // Arrays from HorzMesh + Array2DI4 CellsOnEdge; + Array2DI4 EdgesOnCell; + Array1DI4 NEdgesOnCell; + Array1DReal AreaCell; + Array1DReal DvEdge; + Array2DReal EdgeSignOnCell; + + // small number to avoid numerical difficulties in computing + // scale factors for FCT. + Real Eps = 1e-10; + + const HorzMesh *Mesh; + const VertCoord *VCoord; + + static VertAdv *DefaultVertAdv; + static std::map> AllVertAdvs; + + // private methods + + /// construct a new vertical coordinate object + VertAdv(const std::string &Name, // [in] Name for new VertAdv + const HorzMesh *Mesh, // [in] associated HorzMesh + const VertCoord *VCoord, // [in] associated VertCoord + Config *Options // [in] configuration options + ); + + // Forbid copy and move construction + VertAdv(const VertAdv &) = delete; + VertAdv(VertAdv &&) = delete; + +}; // end class VertAdv + +} // end namespace OMEGA + +//===----------------------------------------------------------------------===// +#endif // defined OMEGA_VERTICALADV_H From f09b3ae1dd62e18af531b1900fa955b4b77f918c Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Sat, 27 Dec 2025 13:27:46 -0800 Subject: [PATCH 02/42] Add optional per-team scratch space to parallelForOuter --- components/omega/src/infra/OmegaKokkos.h | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/components/omega/src/infra/OmegaKokkos.h b/components/omega/src/infra/OmegaKokkos.h index 5a3eb0b82499..a7ab8bca89da 100644 --- a/components/omega/src/infra/OmegaKokkos.h +++ b/components/omega/src/infra/OmegaKokkos.h @@ -78,6 +78,8 @@ using ScratchMemSpace = ExecSpace::scratch_memory_space; using Kokkos::MemoryUnmanaged; using Kokkos::PerTeam; using Kokkos::TeamThreadRange; +using RealScratchArray = + Kokkos::View; /// team_size for hierarchical parallelism #ifdef OMEGA_TARGET_DEVICE @@ -321,7 +323,8 @@ KOKKOS_INLINE_FUNCTION void teamBarrier(const TeamMember &Team) { // parallelForOuter: with label template inline void parallelForOuter(const std::string &Label, - const int (&UpperBounds)[N], F &&Functor) { + const int (&UpperBounds)[N], F &&Functor, + int ScratchValsPerTeam = 0) { auto LinFunctor = LinearIdxWrapper{std::forward(Functor), UpperBounds}; int LinBound = 1; @@ -330,6 +333,12 @@ inline void parallelForOuter(const std::string &Label, } auto Policy = TeamPolicy(LinBound, OMEGA_TEAMSIZE); + + if (ScratchValsPerTeam > 0) { + Policy.set_scratch_size( + 0, Kokkos::PerTeam(ScratchValsPerTeam * sizeof(Real))); + } + Kokkos::parallel_for( Label, Policy, KOKKOS_LAMBDA(const TeamMember &Team) { const int TeamId = Team.league_rank(); @@ -339,8 +348,10 @@ inline void parallelForOuter(const std::string &Label, // parallelForOuter: without label template -inline void parallelForOuter(const int (&UpperBounds)[N], F &&Functor) { - parallelForOuter("", UpperBounds, std::forward(Functor)); +inline void parallelForOuter(const int (&UpperBounds)[N], F &&Functor, + int ScratchValsPerTeam = 0) { + parallelForOuter("", UpperBounds, std::forward(Functor), + ScratchValsPerTeam); } // parallelForInner From ff5fa57c2933708aac7430e2ab3cd5614ac65ce8 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Sat, 27 Dec 2025 14:46:39 -0800 Subject: [PATCH 03/42] Add computeVerticalVelocity --- components/omega/src/ocn/VertAdv.cpp | 91 +++++++++++++++++++++++++++- components/omega/src/ocn/VertAdv.h | 7 +++ 2 files changed, 95 insertions(+), 3 deletions(-) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 4ae65d1d94e3..426fc15ed288 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -2,8 +2,8 @@ // //===----------------------------------------------------------------------===// -#include "Tracers.h" #include "VertAdv.h" +#include "Tracers.h" namespace OMEGA { @@ -16,8 +16,8 @@ std::map> VertAdv::AllVertAdvs; // HorzMesh and VertCoord void VertAdv::init() { - auto Mesh = HorzMesh::getDefault(); - auto VCoord = VertCoord::getDefault(); + auto Mesh = HorzMesh::getDefault(); + auto VCoord = VertCoord::getDefault(); Config *OmegaConfig = Config::getOmegaConfig(); @@ -139,6 +139,7 @@ VertAdv *VertAdv::get(const std::string Name ///< [in] Name of VertAdv } // end get VertAdv +//------------------------------------------------------------------------------ // Read and set config options void VertAdv::readConfigOptions(Config *OmegaConfig) { @@ -203,5 +204,89 @@ void VertAdv::readConfigOptions(Config *OmegaConfig) { } // end readConfigOptions +//------------------------------------------------------------------------------ +// Compute VerticalVelocity and TotalVerticalVelocity from the horizontal +// velocity (NormalVelocity) and the layer thickness used for fluxes through +// edges (FluxLayerThickEdge) +void VertAdv::computeVerticalVelocity(const Array2DReal &NormalVelocity, + const Array2DReal &FluxLayerThickEdge) { + + OMEGA_SCOPE(LocNVertLayers, NVertLayers); + OMEGA_SCOPE(LocAreaCell, Mesh->AreaCell); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocNEOnC, Mesh->NEdgesOnCell); + OMEGA_SCOPE(LocEOnC, Mesh->EdgesOnCell); + OMEGA_SCOPE(LocDvE, Mesh->DvEdge); + OMEGA_SCOPE(LocESOnC, Mesh->EdgeSignOnCell); + + // Loop over all cells owned by the task + parallelForOuter( + "computeVerticalVelocity", {NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + RealScratchArray DivHU(Team.team_scratch(0), LocNVertLayers); + + const Real InvAreaCell = 1._Real / LocAreaCell(ICell); + + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + I4 KRange = vertRangeChunked(KMin, KMax); + + // Compute thickness-weighted divergence of horizontal velocity + // in each layer + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + Real DivHUTmp[VecLength] = {0}; + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int J = 0; J < LocNEOnC(ICell); ++J) { + const I4 JEdge = LocEOnC(ICell, J); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + DivHUTmp[KVec] -= LocDvE(JEdge) * LocESOnC(ICell, J) * + FluxLayerThickEdge(JEdge, K) * + NormalVelocity(JEdge, K) * InvAreaCell; + } + } + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + DivHU(K) = DivHUTmp[KVec]; + } + }); + + Team.team_barrier(); + + // Set velocity through top and bottom interfaces to zero + Kokkos::single( + PerTeam(Team), INNER_LAMBDA() { + VerticalVelocity(ICell, KMin) = 0.; + VerticalVelocity(ICell, KMax + 1) = 0.; + }); + KRange = vertRange(KMin + 1, KMax); + + // Prefix sum of divergence to determine velocity through each + // interface + parallelScanInner( + Team, KRange, INNER_LAMBDA(int K, Real &Accum, bool IsFinal) { + const I4 KRev = KMax - K; + Accum -= DivHU(KRev); + + if (IsFinal) { + VerticalVelocity(ICell, KRev) = Accum; + } + }); + }, + NVertLayers); + + // TODO: currently assuming TotalVerticalVelocity = VerticalVelocity, i.e. + // purely from divergence of horizontal velocity. Need to add optional + // corrections to transport velocity from other contributions, e.g. + // movement of vertical interfaces, contribution of horizontal velocity + // through tilted interface. + deepCopy(TotalVerticalVelocity, VerticalVelocity); + +} // end computeVerticalVelocity + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index 04f780fca123..74516aa6350e 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -92,6 +92,13 @@ class VertAdv { /// Read and set config options void readConfigOptions(Config *Options); + /// Determine transport due to vertical advection from divergence of + /// horizontal advection. + void computeVerticalVelocity( + const Array2DReal &NormalVelocity, + const Array2DReal &FluxLayerThickEdge + ); + private: // Vertical loop bounds from VertCoord Array1DI4 MinLayerCell; From d9deb075aca302c11ca04ed701b47057b1a79211 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Dec 2025 09:14:03 -0800 Subject: [PATCH 04/42] Add computeThicknessVAdvTend --- components/omega/src/ocn/VertAdv.cpp | 44 ++++++++++++++++++++++++++-- components/omega/src/ocn/VertAdv.h | 9 ++++-- 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 426fc15ed288..12d1a9e6d79c 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -208,8 +208,10 @@ void VertAdv::readConfigOptions(Config *OmegaConfig) { // Compute VerticalVelocity and TotalVerticalVelocity from the horizontal // velocity (NormalVelocity) and the layer thickness used for fluxes through // edges (FluxLayerThickEdge) -void VertAdv::computeVerticalVelocity(const Array2DReal &NormalVelocity, - const Array2DReal &FluxLayerThickEdge) { +void VertAdv::computeVerticalVelocity( + const Array2DReal &NormalVelocity, /// [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge /// [in] layer thickness at edges +) { OMEGA_SCOPE(LocNVertLayers, NVertLayers); OMEGA_SCOPE(LocAreaCell, Mesh->AreaCell); @@ -288,5 +290,43 @@ void VertAdv::computeVerticalVelocity(const Array2DReal &NormalVelocity, } // end computeVerticalVelocity +//------------------------------------------------------------------------------ +// Compute thickness tendency due to vertical advection +void VertAdv::computeThicknessVAdvTend( + const Array2DReal &ThickTend // [inout] thickness tendency +) { + + // Return if vertical advection thickness tendency not enabled + if (!ThickVertAdvEnabled) + return; + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocTotVertVelocity, TotalVerticalVelocity); + + // Loop over every owned cell, pseudo thickness tendency is simply + // difference in pseudo velocity between bottom and top interface for + // each layer + parallelForOuter( + "computeThicknessVAdvTend", {NCellsOwned}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + ThickTend(ICell, K) += LocTotVertVelocity(ICell, K + 1) - + LocTotVertVelocity(ICell, K); + } + }); + }); + +} // end computeThicknessVAdvTend + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index 74516aa6350e..f1be78443b5a 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -95,8 +95,13 @@ class VertAdv { /// Determine transport due to vertical advection from divergence of /// horizontal advection. void computeVerticalVelocity( - const Array2DReal &NormalVelocity, - const Array2DReal &FluxLayerThickEdge + const Array2DReal &NormalVelocity, // [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge // [in] layer thickness at edges + ); + + /// Compute pseudo thickness tendency due to vertical advection + void computeThicknessVAdvTend( + const Array2DReal &ThickTend // [inout] thickness tendency ); private: From 93876a34eba7d655dce799b28675be3bed89f8e4 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Dec 2025 09:56:15 -0800 Subject: [PATCH 05/42] Add computeVelocityVAdvTend --- components/omega/src/ocn/VertAdv.cpp | 87 +++++++++++++++++++++++++++- components/omega/src/ocn/VertAdv.h | 7 +++ 2 files changed, 91 insertions(+), 3 deletions(-) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 12d1a9e6d79c..b2345d7ae4c3 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -209,8 +209,8 @@ void VertAdv::readConfigOptions(Config *OmegaConfig) { // velocity (NormalVelocity) and the layer thickness used for fluxes through // edges (FluxLayerThickEdge) void VertAdv::computeVerticalVelocity( - const Array2DReal &NormalVelocity, /// [in] horizontal velocity - const Array2DReal &FluxLayerThickEdge /// [in] layer thickness at edges + const Array2DReal &NormalVelocity, ///< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges ) { OMEGA_SCOPE(LocNVertLayers, NVertLayers); @@ -293,7 +293,7 @@ void VertAdv::computeVerticalVelocity( //------------------------------------------------------------------------------ // Compute thickness tendency due to vertical advection void VertAdv::computeThicknessVAdvTend( - const Array2DReal &ThickTend // [inout] thickness tendency + const Array2DReal &ThickTend ///< [inout] thickness tendency ) { // Return if vertical advection thickness tendency not enabled @@ -328,5 +328,86 @@ void VertAdv::computeThicknessVAdvTend( } // end computeThicknessVAdvTend +//------------------------------------------------------------------------------ +// Compute velocity tendency due to vertical advection +void VertAdv::computeVelocityVAdvTend( + const Array2DReal &VelTend, ///< [inout] horizontal velocity tendency + const Array2DReal &NormalVelocity, ///< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges +) { + + // Return if vertical advection velocity tendency not enabled + if (!VelVertAdvEnabled) + return; + + OMEGA_SCOPE(LocCOnE, Mesh->CellsOnEdge); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + OMEGA_SCOPE(LocTotVertVelocity, TotalVerticalVelocity); + OMEGA_SCOPE(EdgeMask, VCoord->EdgeMask); + OMEGA_SCOPE(LocNVertLayersP1, NVertLayersP1); + + // Loop over every owned edge + parallelForOuter( + "computeVelocityVAdvTend", {NEdgesOwned}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const I4 Cell1 = LocCOnE(IEdge, 0); + const I4 Cell2 = LocCOnE(IEdge, 1); + const I4 KMin = MinLayerEdgeBot(IEdge); + const I4 KMax = MaxLayerEdgeTop(IEdge); + I4 KRange = vertRangeChunked(KMin + 1, KMax); + + // Allocate scratch space for W times Du/Dz at vertical interfaces + // between edges + RealScratchArray WDuDzEdge(Team.team_scratch(0), LocNVertLayersP1); + + // Flux is zero at top and bottom + Kokkos::single( + PerTeam(Team), INNER_LAMBDA() { + WDuDzEdge(Kmin) = 0._Real; + WDuDzEdge(KMax + 1) = 0._Real; + }); + + // Average vertical velocities from cell centers to edges and multiply + // by derivative of horizontal velocity to obtain flux of horizontal + // velocity at the vertical interfaces between edges + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 1); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + const Real WAvg = 0.5_Real * (LocTotVertVelocity(Cell1, K) + + LocTotVertVelocity(Cell2, K)); + WDuDzEdge(K) = + WAvg * + (NormalVelocity(IEdge, K - 1) - + NormalVelocity(IEdge, K)) / + (0.5_Real * (FluxLayerThickEdge(IEdge, K - 1) + + FluxLayerThickEdge(IEdge, K))); + } + }); + + Team.team_barrier(); + + KRange = vertRangeChunked(KMin, KMax); + // Average W*Du/Dz from interfaces to layer midpoints + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + VelTend(IEdge, K) -= EdgeMask(IEdge, K) * 0.5_Real * + (WDuDzEdge(K) + WDuDzEdge(K + 1)); + } + }); + }, + NVertLayersP1); + +} // end computeVelocityVAdvTend + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index f1be78443b5a..9a56b1d398a6 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -104,6 +104,13 @@ class VertAdv { const Array2DReal &ThickTend // [inout] thickness tendency ); + /// Compute velocity tendency due to vertical advection + void computeVelocityVAdvTend( + const Array2DReal &VelTend, // [inout] horizontal velocity tendency + const Array2DReal &NormalVelocity, // [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge // [in] layer thickness at edges + ); + private: // Vertical loop bounds from VertCoord Array1DI4 MinLayerCell; From 43fa0c82f718bb4d8e81ad31b168b458ebf7079c Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Mon, 29 Dec 2025 14:06:08 -0800 Subject: [PATCH 06/42] Add computeTracerVAdvTend dispatch function --- components/omega/src/ocn/VertAdv.cpp | 26 ++++++++++++++++++++++++++ components/omega/src/ocn/VertAdv.h | 9 +++++++++ 2 files changed, 35 insertions(+) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index b2345d7ae4c3..130ef4824277 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -409,5 +409,31 @@ void VertAdv::computeVelocityVAdvTend( } // end computeVelocityVAdvTend +//------------------------------------------------------------------------------ +// Compute tracer tendency due to vertical advection, TimeStep is only needed +// as an arugement for flux-corrected transport +void VertAdv::computeTracerVAdvTend( + const Array3DReal &TracerTend, ///< [inout] tracer tendencies + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &LayerThickness, ///< [in] layer thickness + const TimeInterval TimeStep ///< [in] (optional) time step +) { + + computeVerticalFluxes(Tracers, LayerThickness); + + // dispatch to appropriate algorithm based on configuration settings + switch (VertAdvChoice) { + case VertAdvOption::Standard: + computeStdVAdvTend(TracerTend); + break; + case VertAdvOption::FCT: + R8 Dt; + TimeStep.get(Dt, TimeUnits::Seconds); + computeFCTVAdvTend(TracerTend, Tracers, LayerThickness, Dt); + break; + } + +} // end computeTracerVAdvTend + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index 9a56b1d398a6..a310b13ad044 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -111,6 +111,15 @@ class VertAdv { const Array2DReal &FluxLayerThickEdge // [in] layer thickness at edges ); + /// Compute tracer tendency due to vertical advection, TimeStep is only + /// needed as an argument for flux-corrected transport + void computeTracerVAdvTend( + const Array3DReal &TracerTend, // [inout] tracer tendencies + const Array3DReal &Tracers, // [in] tracer array + const Array2DReal &LayerThickness, // [in] layer thickness + const TimeInterval TimeStep = TimeInterval() // [in] (optional) time step + ); + private: // Vertical loop bounds from VertCoord Array1DI4 MinLayerCell; From 56ca68608728d7bdff5721cf281f3f1a91378d55 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 1 Jan 2026 10:48:58 -0800 Subject: [PATCH 07/42] Add computeVerticalFluxes --- components/omega/src/ocn/VertAdv.cpp | 167 ++++++++++++++++++++++++++- components/omega/src/ocn/VertAdv.h | 6 + 2 files changed, 172 insertions(+), 1 deletion(-) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 130ef4824277..fefe32c31c1d 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -364,7 +364,7 @@ void VertAdv::computeVelocityVAdvTend( // Flux is zero at top and bottom Kokkos::single( PerTeam(Team), INNER_LAMBDA() { - WDuDzEdge(Kmin) = 0._Real; + WDuDzEdge(KMin) = 0._Real; WDuDzEdge(KMax + 1) = 0._Real; }); @@ -435,5 +435,170 @@ void VertAdv::computeTracerVAdvTend( } // end computeTracerVAdvTend +//------------------------------------------------------------------------------ +// Compute tracer fluxes due to vertical advection, the particular scheme used +// is chosen via configuration settings +void VertAdv::computeVerticalFluxes( + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &LayerThickness ///< [in] layer thickness +) { + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocTotVertVel, TotalVerticalVelocity); + OMEGA_SCOPE(LocVertFlux, VertFlux); + OMEGA_SCOPE(LocLowOrderVertFlux, LowOrderVertFlux); + OMEGA_SCOPE(LocCoef3rdOrder, Coef3rdOrder); + + // Compute the fluxes used for the standard VAdv scheme, or the high-order + // fluxes used for the FCT VAdv scheme, store in VertFlux member array + switch (VertFluxChoice) { + // 2nd-order centered fluxes + case VertFluxOption::Second: + parallelForOuter( + "computeVerticalFluxes-Second", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 2, KMax - 1); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 2); + const I4 KLen = chunkLength(KChunk, KStart, KMax - 1); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + const Real InvLayerThickSum = + 1._Real / (LayerThickness(ICell, K - 1) + + LayerThickness(ICell, K)); + const Real VerticalWeightK = + LayerThickness(ICell, K - 1) * InvLayerThickSum; + const Real VerticalWeightKm1 = + LayerThickness(ICell, K) * InvLayerThickSum; + LocVertFlux(L, ICell, K) = + LocTotVertVel(ICell, K) * + (VerticalWeightK * Tracers(L, ICell, K) + + VerticalWeightKm1 * Tracers(L, ICell, K - 1)); + } + }); + }); + break; + // 3rd-order upwind fluxes + case VertFluxOption::Third: + parallelForOuter( + "computeVerticalFluxes-Third", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 2, KMax - 1); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 2); + const I4 KLen = chunkLength(KChunk, KStart, KMax - 1); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + LocVertFlux(L, ICell, K) = + (LocTotVertVel(ICell, K) * + (7._Real * (Tracers(L, ICell, K) + + Tracers(L, ICell, K - 1)) - + (Tracers(L, ICell, K + 1) + + Tracers(L, ICell, K - 2))) - + LocCoef3rdOrder * + std::abs(LocTotVertVel(ICell, K)) * + ((Tracers(L, ICell, K + 1) - + Tracers(L, ICell, K - 2)) - + 3._Real * (Tracers(L, ICell, K) - + Tracers(L, ICell, K - 1)))) / + 12._Real; + } + }); + }); + break; + // 4th-order centered fluxes + case VertFluxOption::Fourth: + parallelForOuter( + "computeVerticalFluxes-Fourth", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 2, KMax - 1); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 2); + const I4 KLen = chunkLength(KChunk, KStart, KMax - 1); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + LocVertFlux(L, ICell, K) = + LocTotVertVel(ICell, K) * + (7._Real * (Tracers(L, ICell, K) + + Tracers(L, ICell, K - 1)) - + (Tracers(L, ICell, K + 1) + + Tracers(L, ICell, K - 2))) / + 12._Real; + } + }); + }); + break; + } + + // Handle fluxes on interfaces near top and bottom layers. Fluxes are 0 on + // top-most (KMin) and bottom-most (KMax + 1) interfaces, use second order + // for fluxes on next-to-top (KMin + 1) and next-to-bottom (KMax) interfaces. + parallelFor( + "computeVerticalFluxes-TopBot", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + for (int K : {KMin, KMax + 1}) { + LocVertFlux(L, ICell, K) = 0._Real; + } + for (int K : {KMin + 1, KMax}) { + const Real InvLayerThickSum = + 1._Real / + (LayerThickness(ICell, K - 1) + LayerThickness(ICell, K)); + const Real VerticalWeightK = + LayerThickness(ICell, K - 1) * InvLayerThickSum; + const Real VerticalWeightKm1 = + LayerThickness(ICell, K) * InvLayerThickSum; + LocVertFlux(L, ICell, K) = + LocTotVertVel(ICell, K) * + (VerticalWeightK * Tracers(L, ICell, K) + + VerticalWeightKm1 * Tracers(L, ICell, K - 1)); + } + }); + + // If using FCT scheme, compute 1st-order upwind fluxes for low-order and + // remove low-order flux from high-order flux + if (VertAdvChoice == VertAdvOption::FCT) { + parallelForOuter( + "computeVerticalFluxes-LowOrder", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin + 1, KMax); + + for (int K : {KMin, KMax + 1}) { + LocLowOrderVertFlux(L, ICell, K) = 0._Real; + } + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 1); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + LocLowOrderVertFlux(L, ICell, K) = + Kokkos::min(0._Real, LocTotVertVel(ICell, K)) * + Tracers(L, ICell, K - 1) + + Kokkos::max(0._Real, LocTotVertVel(ICell, K)) * + Tracers(L, ICell, K); + + LocVertFlux(L, ICell, K) -= + LocLowOrderVertFlux(L, ICell, K); + } + }); + }); + } + +} // end computeVerticalFluxes + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index a310b13ad044..3683a76bd05e 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -120,6 +120,12 @@ class VertAdv { const TimeInterval TimeStep = TimeInterval() // [in] (optional) time step ); + /// Compute tracer fluxes due to vertical advection + void computeVerticalFluxes( + const Array3DReal &Tracers, // [in] tracer array + const Array2DReal &LayerThickness // [in] layer thickness + ); + private: // Vertical loop bounds from VertCoord Array1DI4 MinLayerCell; From 42609bf1e8f3232997c16c696cfdbf7f16558505 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 1 Jan 2026 11:41:44 -0800 Subject: [PATCH 08/42] Add computeStdVAdvTend --- components/omega/src/ocn/VertAdv.cpp | 33 ++++++++++++++++++++++++++++ components/omega/src/ocn/VertAdv.h | 6 +++++ 2 files changed, 39 insertions(+) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index fefe32c31c1d..0f19d365add7 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -600,5 +600,38 @@ void VertAdv::computeVerticalFluxes( } // end computeVerticalFluxes +//------------------------------------------------------------------------------ +// Compute tracer tendencies due to vertical advection using standard advection +// scheme +void VertAdv::computeStdVAdvTend( + const Array3DReal &TracerTend ///< [inout] tracer tendencies +) { + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocVertFlux, VertFlux); + + // Loop over owned cells, tracer tendency in each layer is computed from + // difference between fluxes through bottom and top interfaces + parallelForOuter( + "computeStdVAdvTend", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + const I4 KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + TracerTend(L, ICell, K) += + LocVertFlux(L, ICell, K + 1) - LocVertFlux(L, ICell, K); + } + }); + }); + +} // end computeStdVAdvTend + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index 3683a76bd05e..ee7b198ab672 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -126,6 +126,12 @@ class VertAdv { const Array2DReal &LayerThickness // [in] layer thickness ); + /// Compute tracer tendencies due to vertical advection using standard + /// advection scheme + void + computeStdVAdvTend(const Array3DReal &TracerTend // [inout] tracer tendencies + ); + private: // Vertical loop bounds from VertCoord Array1DI4 MinLayerCell; From bd522d7bf5618ce625a9341ba15540eb6632e98a Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 1 Jan 2026 16:19:09 -0800 Subject: [PATCH 09/42] Add computeFCTVAdvTend --- components/omega/src/ocn/VertAdv.cpp | 153 +++++++++++++++++++++++++++ components/omega/src/ocn/VertAdv.h | 9 ++ 2 files changed, 162 insertions(+) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 0f19d365add7..e7735ac90181 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -633,5 +633,158 @@ void VertAdv::computeStdVAdvTend( } // end computeStdVAdvTend +//------------------------------------------------------------------------------ +// Compute tracer tendencies due to vertical advection using flux-corrected +// transport scheme. ProvThickness input is provisional layer thickness after +// horizontal thickness flux +void VertAdv::computeFCTVAdvTend( + const Array3DReal &TracerTend, ///< [inout] tracer tendencies + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &ProvThickness, ///< [in] provisional layer thickness + const Real Dt ///< [in] time step +) { + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + OMEGA_SCOPE(LocTotVertVel, TotalVerticalVelocity); + OMEGA_SCOPE(LocVertFlux, VertFlux); + OMEGA_SCOPE(LocLowOrderVertFlux, LowOrderVertFlux); + OMEGA_SCOPE(LocNVertLayers, NVertLayers); + OMEGA_SCOPE(LocEps, Eps); + + parallelForOuter( + "computeFCTVAdvTend", {NTracers, NCellsOwned}, + KOKKOS_LAMBDA(int L, int ICell, const TeamMember &Team) { + const I4 KMin = MinLayerCell(ICell); + const I4 KMax = MaxLayerCell(ICell); + I4 KRange = vertRangeChunked(KMin, KMax); + + RealScratchArray InvNewProvThick(Team.team_scratch(0), + LocNVertLayers); + RealScratchArray WorkTend(Team.team_scratch(0), LocNVertLayers); + RealScratchArray FlxIn(Team.team_scratch(0), LocNVertLayers); + RealScratchArray FlxOut(Team.team_scratch(0), LocNVertLayers); + RealScratchArray RescaledFlux(Team.team_scratch(0), + LocNVertLayers + 1); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + InvNewProvThick(K) = + 1._Real / (ProvThickness(ICell, K) + + Dt * (LocTotVertVel(ICell, K + 1) - + LocTotVertVel(ICell, K))); + Real TracerMax; + Real TracerMin; + // Determine bounds on tracer from neighbor values for + // limiting + if (K == KMin) { + TracerMax = Kokkos::max(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1)); + TracerMin = Kokkos::min(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1)); + } else if (K == KMax) { + TracerMax = Kokkos::max(Tracers(L, ICell, K - 1), + Tracers(L, ICell, K)); + TracerMin = Kokkos::min(Tracers(L, ICell, K - 1), + Tracers(L, ICell, K)); + } else { + TracerMax = + Kokkos::max(Tracers(L, ICell, K - 1), + Kokkos::max(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1))); + TracerMin = + Kokkos::min(Tracers(L, ICell, K - 1), + Kokkos::min(Tracers(L, ICell, K), + Tracers(L, ICell, K + 1))); + } + + // Accumulate upwind flux in WorkTend + WorkTend(K) = LocLowOrderVertFlux(L, ICell, K + 1) - + LocLowOrderVertFlux(L, ICell, K); + // Accumulate remaining high-order flux into layer + FlxIn(K) = + Kokkos::max(0._Real, LocVertFlux(L, ICell, K + 1)) - + Kokkos::min(0._Real, LocVertFlux(L, ICell, K)); + // Accumulate remaining high-order flux out of layer + FlxOut(K) = + Kokkos::min(0._Real, LocVertFlux(L, ICell, K + 1)) - + Kokkos::max(0._Real, LocVertFlux(L, ICell, K)); + // Build scale factors to limit flux for FCT using the + // bounds determined above and bounds on newly updated + // values. Factors are stored in FlxIn and FlxOut + // scratch space. + Real TracerMinNew = + (Tracers(L, ICell, K) * ProvThickness(ICell, K) + + Dt * (WorkTend(K) + FlxOut(K))) * + InvNewProvThick(K); + Real TracerMaxNew = + (Tracers(L, ICell, K) * ProvThickness(ICell, K) + + Dt * (WorkTend(K) + FlxIn(K))) * + InvNewProvThick(K); + Real TracerUpwindNew = + (Tracers(L, ICell, K) * ProvThickness(ICell, K) + + Dt * WorkTend(K)) * + InvNewProvThick(K); + Real ScaleFactor = + (TracerMax - TracerUpwindNew) / + (TracerMaxNew - TracerUpwindNew + LocEps); + FlxIn(K) = + Kokkos::min(1._Real, Kokkos::max(0._Real, ScaleFactor)); + ScaleFactor = (TracerUpwindNew - TracerMin) / + (TracerUpwindNew - TracerMinNew + LocEps); + FlxOut(K) = + Kokkos::min(1._Real, Kokkos::max(0._Real, ScaleFactor)); + } + }); + + Team.team_barrier(); + + KRange = vertRangeChunked(KMin + 1, KMax); + + // Rescale the high-order vertical flux + RescaledFlux(KMin) = 0._Real; + RescaledFlux(KMax + 1) = 0._Real; + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin + 1); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + + RescaledFlux(K) = + Kokkos::max(0._Real, LocVertFlux(L, ICell, K)) * + Kokkos::min(FlxOut(K), FlxIn(K - 1)) + + Kokkos::min(0._Real, LocVertFlux(L, ICell, K)) * + Kokkos::min(FlxOut(K - 1), FlxIn(K)); + } + }); + + Team.team_barrier(); + + // Accumulate total FCT vertical advection tendency + KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const I4 KStart = chunkStart(KChunk, KMin); + const I4 KLen = chunkLength(KChunk, KStart, KMax); + + for (int KVec = 0; KVec < KLen; ++KVec) { + const I4 K = KStart + KVec; + WorkTend(K) += RescaledFlux(K + 1) - RescaledFlux(K); + TracerTend(L, ICell, K) += WorkTend(K); + } + }); + // TODO: Monotonicity and diagnostic checks + }, + 5 * NVertLayers + 1); + +} // end computeFTCVAdvTend + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index ee7b198ab672..cbe808720e8c 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -132,6 +132,15 @@ class VertAdv { computeStdVAdvTend(const Array3DReal &TracerTend // [inout] tracer tendencies ); + /// Compute tracer tendencies due to vertical advection using flux-corrected + /// transport scheme + void computeFCTVAdvTend( + const Array3DReal &TracerTend, // [inout] tracer tendencies + const Array3DReal &Tracers, // [in] tracer array + const Array2DReal &ProvThickness, // [in] provisional layer thickness + const Real Dt // [in] time step + ); + private: // Vertical loop bounds from VertCoord Array1DI4 MinLayerCell; From 421ab5044ff33a4483f1e33adeba650614f383c5 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Tue, 6 Jan 2026 11:24:16 -0800 Subject: [PATCH 10/42] Add ProvThickness to auxvars (follow-up restructuring needed) --- components/omega/src/ocn/AuxiliaryState.cpp | 5 +-- .../auxiliaryVars/LayerThicknessAuxVars.cpp | 7 ++-- .../ocn/auxiliaryVars/LayerThicknessAuxVars.h | 32 +++++++++++++++++-- 3 files changed, 37 insertions(+), 7 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 9643fe264cd5..ef3cea652d64 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -196,8 +196,9 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - LocLayerThicknessAux.computeVarsOnCells(ICell, KChunk, - LayerThickCell); + LocLayerThicknessAux.computeVarsOnCells( + ICell, KChunk, LayerThickCell, NormalVelEdge, 0._Real); + // TODO: make timestep available to this call }); }); Pacer::stop("AuxState:cellAuxState3", 2); diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp index cb937863528c..872b68a82b9d 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp @@ -14,8 +14,11 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix, Mesh->NEdgesSize, VCoord->NVertLayers), SshCell("SshCell" + AuxStateSuffix, Mesh->NCellsSize, VCoord->NVertLayers), - CellsOnEdge(Mesh->CellsOnEdge), BottomDepth(Mesh->BottomDepth), - MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + ProvThickness("ProvThickness" + AuxStateSuffix, Mesh->NCellsSize, + VCoord->NVertLayers), + AreaCell(Mesh->AreaCell), DvEdge(Mesh->AreaCell), + EdgeSignOnCell(Mesh->EdgeSignOnCell), CellsOnEdge(Mesh->CellsOnEdge), + BottomDepth(Mesh->BottomDepth), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h index 0620d546fd8f..a415bfa02ba3 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h @@ -17,6 +17,7 @@ class LayerThicknessAuxVars { Array2DReal FluxLayerThickEdge; Array2DReal MeanLayerThickEdge; Array2DReal SshCell; + Array2DReal ProvThickness; FluxThickEdgeOption FluxThickEdgeChoice; @@ -63,9 +64,10 @@ class LayerThicknessAuxVars { } } - KOKKOS_FUNCTION void - computeVarsOnCells(int ICell, int KChunk, - const Array2DReal &LayerThickCell) const { + KOKKOS_FUNCTION void computeVarsOnCells(int ICell, int KChunk, + const Array2DReal &LayerThickCell, + const Array2DReal &NormalVelEdge, + const Real Dt) const { // Temporary for stacked shallow water const int KStart = chunkStart(KChunk, MinLayerCell(ICell)); @@ -76,6 +78,25 @@ class LayerThicknessAuxVars { SshCell(ICell, K) = LayerThickCell(ICell, K) - BottomDepth(ICell); } + Real TmpProv[VecLength] = {0.}; + + Real DtInvAreaCell = Dt / AreaCell(ICell); + for (int J = 0; J < NEdgesOnCell(ICell); ++J) { + const I4 JEdge = EdgesOnCell(ICell, J); + const Real Factor = + DtInvAreaCell * DvEdge(JEdge) * EdgeSignOnCell(ICell, J); + for (int KVec = 0; KVec < KLen; ++KVec) { + const int K = KStart + KVec; + TmpProv[KVec] += + Factor * FluxLayerThickEdge(JEdge, K) * NormalVelEdge(JEdge, K); + } + } + + for (int KVec = 0; KVec < KLen; ++KVec) { + const int K = KStart + KVec; + ProvThickness(ICell, K) = LayerThickCell(ICell, K) + TmpProv[KVec]; + } + /* Real TotalThickness = 0.0; for (int K = 0; K < NVertLayers; K++) { @@ -91,6 +112,11 @@ class LayerThicknessAuxVars { void unregisterFields() const; private: + Array1DReal AreaCell; + Array1DReal DvEdge; + Array1DI4 NEdgesOnCell; + Array2DI4 EdgesOnCell; + Array2DReal EdgeSignOnCell; Array2DI4 CellsOnEdge; Array1DReal BottomDepth; Array1DI4 MinLayerEdgeBot; From 725210c1370f1a9c22d065baea959470f5120772 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 8 Jan 2026 15:54:07 -0800 Subject: [PATCH 11/42] Add new config variables to default config --- components/omega/configs/Default.yml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index d592acd1359a..7718ac3675e3 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -23,8 +23,11 @@ Omega: State: NTimeLevels: 2 Advection: + Coef3rdOrder: 0.25 FluxThicknessType: Center FluxTracerType: Center + VerticalTracerFluxLimiterEnabled: true + VerticalTracerFluxOrder: 3 WindStress: InterpType: Isotropic VertCoord: @@ -50,6 +53,9 @@ Omega: EddyDiff4: 0.0 UseCustomTendency: false ManufacturedSolutionTendency: false + ThicknessVertAdvTendencyEnable: true + VelocityVertAdvTendencyEnable: true + TracerVertAdvTendencyEnable: true Tracers: Base: [Temperature, Salinity] Debug: [Debug1, Debug2, Debug3] From 980fef8bc7a67a9c07e8fc1388ff6673d20911bd Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 8 Jan 2026 15:58:15 -0800 Subject: [PATCH 12/42] Add VertAdv unit tests --- components/omega/test/CMakeLists.txt | 11 + components/omega/test/ocn/VertAdvTest.cpp | 437 ++++++++++++++++++++++ 2 files changed, 448 insertions(+) create mode 100644 components/omega/test/ocn/VertAdvTest.cpp diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt index 2f4b182f3ed0..3aeb12df79e9 100644 --- a/components/omega/test/CMakeLists.txt +++ b/components/omega/test/CMakeLists.txt @@ -476,3 +476,14 @@ add_omega_test( ocn/VertMixTest.cpp "-n;4" ) + +################## +# VAdv test +################## + +add_omega_test( + VERTADV_TEST + testVertAdv.exe + ocn/VertAdvTest.cpp + "-n;1" +) diff --git a/components/omega/test/ocn/VertAdvTest.cpp b/components/omega/test/ocn/VertAdvTest.cpp new file mode 100644 index 000000000000..71f7e448f2e7 --- /dev/null +++ b/components/omega/test/ocn/VertAdvTest.cpp @@ -0,0 +1,437 @@ +//===-- Test driver for OMEGA Vertical Advection ----------------*- C++ -*-===/ +// +/// \file +/// \brief Test driver for OMEGA VertAdv class +/// +/// +// +//===-----------------------------------------------------------------------===/ + +#include "VertAdv.h" +#include "Config.h" +#include "DataTypes.h" +#include "Decomp.h" +#include "GlobalConstants.h" +#include "HorzMesh.h" +#include "IO.h" +#include "IOStream.h" +#include "Logging.h" +#include "OceanState.h" +#include "OmegaKokkos.h" +#include "TimeStepper.h" +#include "VertCoord.h" + +#include "Pacer.h" +#include +#include +#include + +using namespace OMEGA; + +constexpr Real Alpha = 0.197_Real; // random phase angle +constexpr Real Length = 1024._Real; + +// Test function used for tendTest +KOKKOS_FUNCTION Real testFunc(const Real X) { + Real Arg = Pi * (X + Alpha) / Length; + return std::sin(Arg) * std::sin(Arg) * std::cos(Arg) * std::cos(Arg); +} + +// Function used for analytical derivative in tendTest +KOKKOS_FUNCTION Real testDeriv(const Real X) { + Real C = Pi / Length; + return (C / 2._Real) * std::sin(4._Real * C * (X + Alpha)); +} + +// For a given resolution, use computeStdVAdvTend to compute tracer +// tendencies, then calculate L2 error +Real tendTest(const int NLayers, VertAdv *VAdv) { + + // Setup test for this resolution + Real Delta = Length / NLayers; + + VAdv->NVertLayers = NLayers; + VAdv->TotalVerticalVelocity = + Array2DReal("TotalVerticaVelocity", 1, NLayers + 1); + Array1DReal XLayer("XLayer", NLayers); + OMEGA_SCOPE(LocTotVertVel, VAdv->TotalVerticalVelocity); + Array2DReal LayerThick("LayerThickness", 1, NLayers); + Array3DReal TracerArray("TracerArray", 1, 1, NLayers); + Array3DReal Tend("TracerArray", 1, 1, NLayers); + + // Set uniform velocity throughout domain + deepCopy(VAdv->TotalVerticalVelocity, 1._Real); + // Set velocities at top and bottom interfaces to 0. + parallelFor( + {1}, KOKKOS_LAMBDA(const int &) { + LocTotVertVel(0, 0) = 0._Real; + LocTotVertVel(0, NLayers) = 0._Real; + }); + // Set X values, LayerThickness, and Tracer values throughout domain + parallelFor( + {NLayers}, KOKKOS_LAMBDA(int K) { + XLayer(K) = Delta * (K + 0.5_Real); + LayerThick(0, K) = Delta; + TracerArray(0, 0, K) = testFunc(XLayer(K)); + }); + + // Compute tracer tendencies + VAdv->computeTracerVAdvTend(Tend, TracerArray, LayerThick); + + HostArray3DReal TendH = createHostMirrorCopy(Tend); + HostArray1DReal XLayerH = createHostMirrorCopy(XLayer); + Real SumSq = 0._Real; + int Count = 0; + + // Compute errors at interior of range, ignoring layers near top and bottom + // where lower order is used + for (int K = 2; K < NLayers - 2; ++K) { + Real Expected = testDeriv(XLayerH(K)); + Real Diff = TendH(0, 0, K) / Delta - Expected; + SumSq += Diff * Diff; + Count++; + } + Real L2 = std::sqrt(SumSq / Count); + + return L2; +} + +// Initialize needed modules +void initVertAdvTest() { + + I4 Err; + + MachEnv::init(MPI_COMM_WORLD); + MachEnv *DefEnv = MachEnv::getDefault(); + MPI_Comm DefComm = DefEnv->getComm(); + + // Initialize the Logging system + initLogging(DefEnv); + + // Open config file + Config("Omega"); + Config::readAll("omega.yml"); + + // First step of time stepper initialization needed for IOstream + TimeStepper::init1(); + + // Initialize the IO system + IO::init(DefComm); + + // Create the default decomposition (initializes the decomposition) + Decomp::init(); + + // Initialize streams + IOStream::init(); + + // Initialize the default halo + Err = Halo::init(); + if (Err != 0) + ABORT_ERROR("VertAdvTest: error initializing default halo"); + + // Initialize the default mesh + HorzMesh::init(); + + // Initialize the default vertical coordinate + VertCoord::init(false); + + // Initialize tracers + Tracers::init(); + + // Initialize the default vertical advection + VertAdv::init(); +} + +// Clean-up modules +void finalizeVertAdvTest() { + + IOStream::finalize(); + Tracers::clear(); + VertAdv::clear(); + VertCoord::clear(); + TimeStepper::clear(); + HorzMesh::clear(); + Field::clear(); + Dimension::clear(); + Halo::clear(); + Decomp::clear(); + MachEnv::removeAll(); +} + +int main(int argc, char *argv[]) { + + // Initialize error code + Error ErrAll; + + // Initialize the global MPI environment + MPI_Init(&argc, &argv); + Kokkos::initialize(); + Pacer::initialize(MPI_COMM_WORLD); + Pacer::setPrefix("Omega:"); + { + initVertAdvTest(); + + auto DefMesh = HorzMesh::getDefault(); + auto DefDecomp = Decomp::getDefault(); + auto DefVertAdv = VertAdv::getDefault(); + auto DefVCoord = VertCoord::getDefault(); + + // Only need to test within one column or along one edge and with one + // tracer + DefVertAdv->NCellsOwned = 1; + DefVertAdv->NCellsAll = 1; + DefVertAdv->NEdgesOwned = 1; + DefVertAdv->NEdgesAll = 1; + DefVertAdv->NTracers = 1; + + I4 NVertLayers = DefVertAdv->NVertLayers; + I4 NVertLayersP1 = DefVertAdv->NVertLayersP1; + I4 NTracers = DefVertAdv->NTracers; + + Array2DReal FluxLayerThickEdge("FluxLayerThickEdge", + DefDecomp->NEdgesSize, NVertLayers); + Array2DReal NormalVelEdge("NormalVelEdge", DefDecomp->NEdgesSize, + NVertLayers); + + const I4 ICell0 = 0; + const I4 NEdgesOnCell0 = DefDecomp->NEdgesOnCell(ICell0); + + // Test for computeVerticalVelocity + I4 Err = 0; + + OMEGA_SCOPE(LocEOnC, DefDecomp->EdgesOnCell); + OMEGA_SCOPE(LocESOnC, DefMesh->EdgeSignOnCell); + + // Set unifrom values for layer thickness and velocity on edges of column + parallelFor( + {NVertLayers}, KOKKOS_LAMBDA(int K) { + for (int J = 0; J < NEdgesOnCell0; ++J) { + const I4 JEdge = LocEOnC(ICell0, J); + FluxLayerThickEdge(JEdge, K) = 1._Real; + NormalVelEdge(JEdge, K) = 1._Real * LocESOnC(ICell0, J); + } + }); + + // Compute vertical velocity + DefVertAdv->computeVerticalVelocity(NormalVelEdge, FluxLayerThickEdge); + + HostArray2DReal VertVelH = + createHostMirrorCopy(DefVertAdv->VerticalVelocity); + + Real Tol = 1e-10; + + // Divergence is equal in each layer. Expected vertical velocity through + // each interface is sum over negative divergence starting from bottom + // layer + Real Perim0 = 0.; + for (int J = 0; J < NEdgesOnCell0; ++J) { + Perim0 += DefMesh->DvEdge(DefMesh->EdgesOnCell(ICell0, J)); + } + Real InvAreaCell0 = 1._Real / DefMesh->AreaCell(ICell0); + + for (int K = 1; K < NVertLayersP1; ++K) { + Real Expected = (NVertLayers - K) * Perim0 * InvAreaCell0; + Real Diff = std::abs(VertVelH(ICell0, K) - Expected); + if (Diff > Tol) { + ++Err; + } + } + + if (Err != 0) { + ErrAll += Error(ErrorCode::Fail, + "VertAdvTest: computeVerticalVelocity FAIL"); + } + + // Test for computeThicknessVAdvTend + Err = 0; + + Array2DReal TendCell2D("TendCell2D", DefMesh->NCellsSize, NVertLayers); + // Compute thickness tendencies for each layer using velocities from + // previous tests + DefVertAdv->computeThicknessVAdvTend(TendCell2D); + HostArray2DReal TendCell2DH = createHostMirrorCopy(TendCell2D); + + // Expected thickness tendency for each layer is difference between + // velocity at bottom and top interface, i.e. the divergence + for (int K = 1; K < NVertLayers; ++K) { + Real Expected = -Perim0 * InvAreaCell0; + Real Diff = std::abs(TendCell2DH(ICell0, K) - Expected); + if (Diff > Tol) { + ++Err; + } + } + + if (Err != 0) { + ErrAll += Error(ErrorCode::Fail, + "VertAdvTest: computeThicknessVAdvTend FAIL"); + } + + // Test for computeVelocityVAdvTend + Err = 0; + + // Set uniform vertical velocity in each column along the edge and + // uniform layer thickness at the edge. Horizontal velocity is set + // to a periodic function + Array2DReal TendEdge2D("TendEdge2D", DefMesh->NEdgesSize, NVertLayers); + const I4 IEdge0 = 0; + OMEGA_SCOPE(LocCOnE, DefDecomp->CellsOnEdge); + OMEGA_SCOPE(LocVertVel, DefVertAdv->TotalVerticalVelocity); + parallelFor( + {NVertLayers}, KOKKOS_LAMBDA(int K) { + NormalVelEdge(IEdge0, K) = sin(Pi * K / NVertLayers); + FluxLayerThickEdge(IEdge0, K) = 1._Real; + const I4 Cell1 = LocCOnE(IEdge0, 0); + const I4 Cell2 = LocCOnE(IEdge0, 1); + LocVertVel(Cell1, K) = 1._Real; + LocVertVel(Cell2, K) = 1._Real; + }); + + // Compute velocity tendencies for each layer + DefVertAdv->computeVelocityVAdvTend(TendEdge2D, NormalVelEdge, + FluxLayerThickEdge); + + Tol = 1._Real / (NVertLayers * NVertLayers); + + HostArray2DReal TendEdge2DH = createHostMirrorCopy(TendEdge2D); + + // Expected velocity tendency is derivative of initial distribution + for (int K = 0; K < NVertLayers; ++K) { + Real Expected = (Pi / NVertLayers) * cos(Pi * K / NVertLayers); + if (K == 0 or K == NVertLayers - 1) + Expected /= 2._Real; + Real Diff = std::abs(TendEdge2DH(IEdge0, K) - Expected); + if (Diff > Tol) { + ++Err; + } + } + + if (Err != 0) { + ErrAll += Error(ErrorCode::Fail, + "VertAdvTest: computeVelocityVAdvTend FAIL"); + } + + // Tests for computeTracerVAdvTend + + DefVertAdv->VertAdvChoice = VertAdvOption::Standard; + + // Setup for convergence tests with computeStdVAdvTend + constexpr I4 NLayersArray[] = {128, 256, 512, 1024}; + constexpr VertFluxOption AllOrders[] = {VertFluxOption::Second, + VertFluxOption::Third, + VertFluxOption::Fourth}; + static const std::map OrderOfAcc = { + {VertFluxOption::Second, 2}, + {VertFluxOption::Third, 3}, + {VertFluxOption::Fourth, 4}}; + static const std::map OrderStr = { + {VertFluxOption::Second, "2nd"}, + {VertFluxOption::Third, "3rd"}, + {VertFluxOption::Fourth, "4th"}}; + + // Set Coef3rdOrder to 1 so VertFluxOption::Third is calculated with + // purely 3rd-order and not a blend of 3rd- and 4th-order fluxes + DefVertAdv->Coef3rdOrder = 1.; + + Tol = 0.1_Real; + + // Compute L2 error for computeStdVAdvTend with each VertFluxOption over a + // set of resolutions. Compare successive errors for each Order to confirm + // expected order of accuracy. + for (VertFluxOption Order : AllOrders) { + Err = 0; + DefVertAdv->VertFluxChoice = Order; + std::vector L2Errors; + for (I4 NLayers : NLayersArray) { + DefVCoord->MaxLayerCell(0) = NLayers - 1; + Real L2Err = tendTest(NLayers, DefVertAdv); + L2Errors.push_back(L2Err); + } + Real ExpectedRat = std::pow(2._Real, OrderOfAcc.at(Order)); + for (I4 I = 0; I < L2Errors.size() - 1; ++I) { + Real Rat = L2Errors[I] / L2Errors[I + 1]; + Real RelErr = std::abs(Rat - ExpectedRat) / ExpectedRat; + if (RelErr > Tol) { + ++Err; + } + } + if (Err != 0) { + ErrAll += + Error(ErrorCode::Fail, + "VertAdvTest: computeStdVAdvTend with {} order FAIL", + OrderStr.at(Order)); + } + } + + // Monotonicity check for computeFCTVAdvTend + + NVertLayers = 256; + DefVertAdv->VertAdvChoice = VertAdvOption::FCT; + DefVertAdv->VertFluxChoice = VertFluxOption::Third; + DefVertAdv->Coef3rdOrder = 0.25; + std::vector TestVel = {1._Real, -1._Real}; + std::vector StartIdx = {NVertLayers / 2, NVertLayers / 4}; + + // Test a top-hat function for a uniform velocity in both directions + for (int IMono = 0; IMono < 2; ++IMono) { + Err = 0; + Array3DReal TracerArray("TracerArray", 1, 1, NVertLayers); + Array2DReal LayerThick("LayerThick", 1, NVertLayers); + Array3DReal TendCell3D("TendCell3D", 1, 1, NVertLayers); + deepCopy(LayerThick, 1._Real); + DefVCoord->MaxLayerCell(0) = NVertLayers - 1; + DefVertAdv->NVertLayers = NVertLayers; + LocVertVel = Array2DReal("TotalVerticaVelocity", 1, NVertLayers + 1); + deepCopy(LocVertVel, TestVel[IMono]); + parallelFor( + {1}, KOKKOS_LAMBDA(const int &) { + LocVertVel(0, 0) = 0._Real; + LocVertVel(0, NVertLayers) = 1._Real; + }); + deepCopy(TracerArray, 0._Real); + const I4 KStart = StartIdx[IMono]; + const I4 KRange = NVertLayers / 4; + parallelFor( + {KRange}, KOKKOS_LAMBDA(const int K) { + TracerArray(0, 0, KStart + K) = 1._Real; + }); + Real Dt = 0.5_Real; + TimeInterval TimeStep(Dt, TimeUnits::Seconds); + + const I4 NTSteps = 100; + + // Integrate 100 time steps and confirm monotonicity + for (int N = 0; N < NTSteps; ++N) { + deepCopy(TendCell3D, 0._Real); + DefVertAdv->computeTracerVAdvTend(TendCell3D, TracerArray, + LayerThick, TimeStep); + parallelFor( + {NVertLayers}, KOKKOS_LAMBDA(const int K) { + TracerArray(0, 0, K) += Dt * TendCell3D(0, 0, K); + }); + } + HostArray3DReal TracerH = createHostMirrorCopy(TracerArray); + HostArray3DReal TendH = createHostMirrorCopy(TendCell3D); + for (int K = 0; K < NVertLayers; ++K) { + if (TracerH(0, 0, K) < 0._Real or TracerH(0, 0, K) > 1._Real) { + ++Err; + } + } + if (Err != 0) { + ErrAll += + Error(ErrorCode::Fail, + "VertAdvTest: computeFCTVAdvTend with velocity = {} " + "monotonicity FAIL", + TestVel[IMono]); + } + } + + finalizeVertAdvTest(); + } + Pacer::finalize(); + Kokkos::finalize(); + MPI_Finalize(); + + CHECK_ERROR_ABORT(ErrAll, "VertAdv unit tests FAIL"); + + return 0; +} From 9e9d60696f5c45cfd53560bc3ae0a0c5273f0be7 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Fri, 9 Jan 2026 14:11:45 -0800 Subject: [PATCH 13/42] Init some previously missed pointers and add Field for ProvThickness --- .../ocn/auxiliaryVars/LayerThicknessAuxVars.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp index 872b68a82b9d..506a328749e0 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp @@ -17,6 +17,7 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix, ProvThickness("ProvThickness" + AuxStateSuffix, Mesh->NCellsSize, VCoord->NVertLayers), AreaCell(Mesh->AreaCell), DvEdge(Mesh->AreaCell), + NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), CellsOnEdge(Mesh->CellsOnEdge), BottomDepth(Mesh->BottomDepth), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), @@ -81,21 +82,37 @@ void LayerThicknessAuxVars::registerFields(const std::string &AuxGroupName, DimNames // dimension names ); + // Provisional Thickness + auto ProvThicknessField = Field::create( + ProvThickness.label(), // field name + "layer thickness after horizontal flux", // long Name or description + "m", // units + "", // 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 + ); + // 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()); } } // namespace OMEGA From d859c845df05868bd4413b0168c8f3c1810b55c0 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Fri, 9 Jan 2026 14:14:50 -0800 Subject: [PATCH 14/42] Comments, bug fix, and add Fields for VertAdv arrays --- components/omega/src/ocn/VertAdv.cpp | 184 ++++++++++++++++++++++----- components/omega/src/ocn/VertAdv.h | 95 +++++++++----- 2 files changed, 212 insertions(+), 67 deletions(-) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index e7735ac90181..97b8ab813e67 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -1,10 +1,18 @@ //===-- ocn/VertAdv.cpp - vertical advection --------------------*- C++ -*-===// // +// The VertAdv class contains methods for calculating the vertical velocity and +// the tendencies of thickness, horizontal velocity, and tracers due to vertical +// advective transport, and member variables for storing the vertical velocity +// and the fluxes of tracers at vertical interfaces in each column. +// //===----------------------------------------------------------------------===// #include "VertAdv.h" +#include "Field.h" #include "Tracers.h" +#include + namespace OMEGA { // create static class members @@ -13,7 +21,7 @@ std::map> VertAdv::AllVertAdvs; //------------------------------------------------------------------------------ // create the default VertAdv, requires prior initialization of default -// HorzMesh and VertCoord +// HorzMesh, VertCoord, and Tracers void VertAdv::init() { auto Mesh = HorzMesh::getDefault(); @@ -27,10 +35,10 @@ void VertAdv::init() { //------------------------------------------------------------------------------ // constructor -VertAdv::VertAdv(const std::string &Name_, ///< [in] name for new VertAdv - const HorzMesh *InMesh, ///< [in] associated HorzMesh - const VertCoord *InVCoord, ///< [in] associated VertCoord - Config *Options ///< [in] configuration options +VertAdv::VertAdv(const std::string &Name_, //< [in] name for new VertAdv + const HorzMesh *InMesh, //< [in] associated HorzMesh + const VertCoord *InVCoord, //< [in] associated VertCoord + Config *Options //< [in] configuration options ) { // Read options from Config object @@ -67,14 +75,16 @@ VertAdv::VertAdv(const std::string &Name_, ///< [in] name for new VertAdv Array3DReal("LowOrderVertFlux", NTracers, NCellsSize, NVertLayersP1); } + defineFields(); + } // end constructor //------------------------------------------------------------------------------ // create a new VertAdv instance -VertAdv *VertAdv::create(const std::string &Name, ///< [in] name for new VertAdv - const HorzMesh *Mesh, ///< [in] associated HorzMesh - const VertCoord *VCoord, ///< [in] associated VertCoord - Config *Options ///< [in] configuration options +VertAdv *VertAdv::create(const std::string &Name, //< [in] name for new VertAdv + const HorzMesh *Mesh, //< [in] associated HorzMesh + const VertCoord *VCoord, //< [in] associated VertCoord + Config *Options //< [in] configuration options ) { // Check to see if a VertAdv of the same name already exists and, if so, // exit with an error @@ -94,9 +104,117 @@ VertAdv *VertAdv::create(const std::string &Name, ///< [in] name for new VertAdv } // end create +void VertAdv::defineFields() { + + // set field names (append Name if not default) + VerticalVelocityFldName = "VerticalVelocity"; + TotalVertVelocityFldName = "TotalVerticalVelocity"; + VertFluxFldName = "VertFlux"; + LowOrderVertFluxFldName = "LowOrderVertFlux"; + + if (Name != "Default") { + VerticalVelocityFldName.append(Name); + TotalVertVelocityFldName.append(Name); + VertFluxFldName.append(Name); + LowOrderVertFluxFldName.append(Name); + } + + const Real FillValueReal = -9.99e30; + I4 NDims = 2; + std::vector DimNames(NDims); + DimNames[0] = "NCells"; + DimNames[1] = "NVertLayersP1"; + + auto VerticalVelocityField = Field::create( + VerticalVelocityFldName, // field name + "Vertical pseudovelocity across a pseudoheight surface", // long name or + // description + "m s^-1", // units + "", // 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 + ); + + auto TotalVertVelocityField = + Field::create(TotalVertVelocityFldName, // field name + "Total vertical pseudovelocity across a moving, tilted " + "pseudoheight surface", // long name or description + "m s^-1", // units + "", // 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 = 3; + DimNames.resize(NDims); + DimNames[2] = "NTracers"; + + auto VertFluxField = Field::create( + VertFluxFldName, // field name + "Vertical flux of tracers across a pseudoheight surface", // long name or + // description + "", // units + "", // 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 + ); + + auto LowOrderVertFluxField = + Field::create(LowOrderVertFluxFldName, // field name + "Low-order vertical flux of tracers across a pseudoheight " + "surface", // long name or description + "", // units + "", // 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 + ); + + GroupName = "VertAdv"; + if (Name != "Default") { + GroupName.append(Name); + } + + // Create a field group for VertAdv fields + auto VertAdvGroup = FieldGroup::create(GroupName); + + VertAdvGroup->addField(VerticalVelocityFldName); + VertAdvGroup->addField(TotalVertVelocityFldName); + VertAdvGroup->addField(VertFluxFldName); + VertAdvGroup->addField(LowOrderVertFluxFldName); + + // Associate Fields with data + VerticalVelocityField->attachData(VerticalVelocity); + TotalVertVelocityField->attachData(TotalVerticalVelocity); + VertFluxField->attachData(VertFlux); + LowOrderVertFluxField->attachData(LowOrderVertFlux); + +} // end defineFields + //------------------------------------------------------------------------------ // destructor -VertAdv::~VertAdv() {} // end destructor +VertAdv::~VertAdv() { + + if (FieldGroup::exists(GroupName)) { + Field::destroy(VerticalVelocityFldName); + Field::destroy(TotalVertVelocityFldName); + Field::destroy(VertFluxFldName); + Field::destroy(LowOrderVertFluxFldName); + FieldGroup::destroy(GroupName); + } + +} // end destructor //------------------------------------------------------------------------------ // Removes all VertAdvs to clean up before exit @@ -120,7 +238,7 @@ VertAdv *VertAdv::getDefault() { return VertAdv::DefaultVertAdv; } //------------------------------------------------------------------------------ // Get VertAdv by name -VertAdv *VertAdv::get(const std::string Name ///< [in] Name of VertAdv +VertAdv *VertAdv::get(const std::string Name //< [in] Name of VertAdv ) { // look for an instance of this name @@ -209,10 +327,11 @@ void VertAdv::readConfigOptions(Config *OmegaConfig) { // velocity (NormalVelocity) and the layer thickness used for fluxes through // edges (FluxLayerThickEdge) void VertAdv::computeVerticalVelocity( - const Array2DReal &NormalVelocity, ///< [in] horizontal velocity - const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges + const Array2DReal &NormalVelocity, //< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge //< [in] layer thickness at edges ) { + OMEGA_SCOPE(LocVertVel, VerticalVelocity); OMEGA_SCOPE(LocNVertLayers, NVertLayers); OMEGA_SCOPE(LocAreaCell, Mesh->AreaCell); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); @@ -262,8 +381,8 @@ void VertAdv::computeVerticalVelocity( // Set velocity through top and bottom interfaces to zero Kokkos::single( PerTeam(Team), INNER_LAMBDA() { - VerticalVelocity(ICell, KMin) = 0.; - VerticalVelocity(ICell, KMax + 1) = 0.; + LocVertVel(ICell, KMin) = 0.; + LocVertVel(ICell, KMax + 1) = 0.; }); KRange = vertRange(KMin + 1, KMax); @@ -275,7 +394,7 @@ void VertAdv::computeVerticalVelocity( Accum -= DivHU(KRev); if (IsFinal) { - VerticalVelocity(ICell, KRev) = Accum; + LocVertVel(ICell, KRev) = Accum; } }); }, @@ -293,7 +412,7 @@ void VertAdv::computeVerticalVelocity( //------------------------------------------------------------------------------ // Compute thickness tendency due to vertical advection void VertAdv::computeThicknessVAdvTend( - const Array2DReal &ThickTend ///< [inout] thickness tendency + const Array2DReal &ThickTend //< [inout] thickness tendency ) { // Return if vertical advection thickness tendency not enabled @@ -331,9 +450,9 @@ void VertAdv::computeThicknessVAdvTend( //------------------------------------------------------------------------------ // Compute velocity tendency due to vertical advection void VertAdv::computeVelocityVAdvTend( - const Array2DReal &VelTend, ///< [inout] horizontal velocity tendency - const Array2DReal &NormalVelocity, ///< [in] horizontal velocity - const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges + const Array2DReal &VelTend, //< [inout] horizontal velocity tendency + const Array2DReal &NormalVelocity, //< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge //< [in] layer thickness at edges ) { // Return if vertical advection velocity tendency not enabled @@ -413,15 +532,16 @@ void VertAdv::computeVelocityVAdvTend( // Compute tracer tendency due to vertical advection, TimeStep is only needed // as an arugement for flux-corrected transport void VertAdv::computeTracerVAdvTend( - const Array3DReal &TracerTend, ///< [inout] tracer tendencies - const Array3DReal &Tracers, ///< [in] tracer array - const Array2DReal &LayerThickness, ///< [in] layer thickness - const TimeInterval TimeStep ///< [in] (optional) time step + const Array3DReal &TracerTend, //< [inout] tracer tendencies + const Array3DReal &Tracers, //< [in] tracer array + const Array2DReal &LayerThickness, //< [in] layer thickness + const TimeInterval TimeStep //< [in] (optional) time step ) { + // Compute tracer fluxes at the interfaces computeVerticalFluxes(Tracers, LayerThickness); - // dispatch to appropriate algorithm based on configuration settings + // Dispatch to appropriate algorithm based on configuration settings switch (VertAdvChoice) { case VertAdvOption::Standard: computeStdVAdvTend(TracerTend); @@ -439,8 +559,8 @@ void VertAdv::computeTracerVAdvTend( // Compute tracer fluxes due to vertical advection, the particular scheme used // is chosen via configuration settings void VertAdv::computeVerticalFluxes( - const Array3DReal &Tracers, ///< [in] tracer array - const Array2DReal &LayerThickness ///< [in] layer thickness + const Array3DReal &Tracers, //< [in] tracer array + const Array2DReal &LayerThickness //< [in] layer thickness ) { OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); @@ -604,7 +724,7 @@ void VertAdv::computeVerticalFluxes( // Compute tracer tendencies due to vertical advection using standard advection // scheme void VertAdv::computeStdVAdvTend( - const Array3DReal &TracerTend ///< [inout] tracer tendencies + const Array3DReal &TracerTend //< [inout] tracer tendencies ) { OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); @@ -638,10 +758,10 @@ void VertAdv::computeStdVAdvTend( // transport scheme. ProvThickness input is provisional layer thickness after // horizontal thickness flux void VertAdv::computeFCTVAdvTend( - const Array3DReal &TracerTend, ///< [inout] tracer tendencies - const Array3DReal &Tracers, ///< [in] tracer array - const Array2DReal &ProvThickness, ///< [in] provisional layer thickness - const Real Dt ///< [in] time step + const Array3DReal &TracerTend, //< [inout] tracer tendencies + const Array3DReal &Tracers, //< [in] tracer array + const Array2DReal &ProvThickness, //< [in] provisional layer thickness + const Real Dt //< [in] time step ) { OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index cbe808720e8c..5a173fdd7fb7 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -1,6 +1,13 @@ #ifndef OMEGA_VERTADV_H #define OMEGA_VERTADV_H -//===-- ocn/VertAdv.h - vertical advection ----------------------*- C++ -*-===// +//===-- ocn/VertAdv.h - vertical advection definitions ----------*- C++ -*-===// +// +/// \file +/// \brief Contains methods and variables for vertical advection +/// +/// This header defines the VertAdv class which contains methods and variables +/// for calculating the vertical velocity and for vertical transport of +/// thickness, velocity, and tracers. // //===----------------------------------------------------------------------===// @@ -50,16 +57,26 @@ class VertAdv { // Number of tracers I4 NTracers; - // Coefficient for blending 3rd-order and 4th-order reconstruction of tracers + // Coefficient for blending 3rd-order and 4th-order reconstruction of + // tracers with VertFluxOption::Third. Coef3rdOrder = 1 gives purely + // 3rd-order, Coef3rdOrder = 0 gives purely 4th-order. Real Coef3rdOrder; - // VertAdv instance name + Array2DReal VerticalVelocity; ///< pseudovelocity through top of cell + Array2DReal + TotalVerticalVelocity; ///< transport velocity through top of Cell + Array3DReal VertFlux; ///< fluxes at vertical interfaces + Array3DReal LowOrderVertFlux; ///< low-order fluxes for FCT + + // VertAdv instance name and group name std::string Name; + std::string GroupName; - Array2DReal VerticalVelocity; // pseudovelocity through top of cell - Array2DReal TotalVerticalVelocity; // transport velocity through top of Cell - Array3DReal VertFlux; // fluxes at vertical interfaces - Array3DReal LowOrderVertFlux; // low-order fluxes for FCT + // Field names + std::string VerticalVelocityFldName; + std::string TotalVertVelocityFldName; + std::string VertFluxFldName; + std::string LowOrderVertFluxFldName; // public methods @@ -68,10 +85,11 @@ class VertAdv { /// Creates a new vertical advection object by calling the constructor and /// puts it in the AllVertAdvs map. - static VertAdv *create(const std::string &Name, // [in] name for new VertAdv - const HorzMesh *Mesh, // [in] associated HorzMesh - const VertCoord *VCoord, // [in] associated VertCoord - Config *Options // [in] configuration options + static VertAdv * + create(const std::string &Name, ///< [in] name for new VertAdv + const HorzMesh *Mesh, ///< [in] associated HorzMesh + const VertCoord *VCoord, ///< [in] associated VertCoord + Config *Options ///< [in] configuration options ); /// Destructor - deallocates all memory and deletes a VertAdv @@ -95,50 +113,54 @@ class VertAdv { /// Determine transport due to vertical advection from divergence of /// horizontal advection. void computeVerticalVelocity( - const Array2DReal &NormalVelocity, // [in] horizontal velocity - const Array2DReal &FluxLayerThickEdge // [in] layer thickness at edges + const Array2DReal &NormalVelocity, ///< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges ); /// Compute pseudo thickness tendency due to vertical advection void computeThicknessVAdvTend( - const Array2DReal &ThickTend // [inout] thickness tendency + const Array2DReal &ThickTend ///< [inout] thickness tendency ); /// Compute velocity tendency due to vertical advection void computeVelocityVAdvTend( - const Array2DReal &VelTend, // [inout] horizontal velocity tendency - const Array2DReal &NormalVelocity, // [in] horizontal velocity - const Array2DReal &FluxLayerThickEdge // [in] layer thickness at edges + const Array2DReal &VelTend, ///< [inout] horizontal velocity tendency + const Array2DReal &NormalVelocity, ///< [in] horizontal velocity + const Array2DReal &FluxLayerThickEdge ///< [in] layer thickness at edges ); - /// Compute tracer tendency due to vertical advection, TimeStep is only - /// needed as an argument for flux-corrected transport + /// Compute tracer tendency due to vertical advection, The LayerThickness at + /// the beginning of the time step is passed for standard advection, while + /// the provisional thickness including horizontal thickness flux is passed + /// for flux-corrected transport. The TimeStep is only needed as an argument + /// for flux-corrected transport void computeTracerVAdvTend( - const Array3DReal &TracerTend, // [inout] tracer tendencies - const Array3DReal &Tracers, // [in] tracer array - const Array2DReal &LayerThickness, // [in] layer thickness - const TimeInterval TimeStep = TimeInterval() // [in] (optional) time step + const Array3DReal &TracerTend, ///< [inout] tracer tendencies + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &LayerThickness, ///< [in] layer thickness + const TimeInterval TimeStep = + TimeInterval() ///< [in] (optional) time step ); /// Compute tracer fluxes due to vertical advection void computeVerticalFluxes( - const Array3DReal &Tracers, // [in] tracer array - const Array2DReal &LayerThickness // [in] layer thickness + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &LayerThickness ///< [in] layer thickness ); /// Compute tracer tendencies due to vertical advection using standard /// advection scheme - void - computeStdVAdvTend(const Array3DReal &TracerTend // [inout] tracer tendencies + void computeStdVAdvTend( + const Array3DReal &TracerTend ///< [inout] tracer tendencies ); /// Compute tracer tendencies due to vertical advection using flux-corrected /// transport scheme void computeFCTVAdvTend( - const Array3DReal &TracerTend, // [inout] tracer tendencies - const Array3DReal &Tracers, // [in] tracer array - const Array2DReal &ProvThickness, // [in] provisional layer thickness - const Real Dt // [in] time step + const Array3DReal &TracerTend, ///< [inout] tracer tendencies + const Array3DReal &Tracers, ///< [in] tracer array + const Array2DReal &ProvThickness, ///< [in] provisional layer thickness + const Real Dt ///< [in] time step ); private: @@ -169,12 +191,15 @@ class VertAdv { // private methods /// construct a new vertical coordinate object - VertAdv(const std::string &Name, // [in] Name for new VertAdv - const HorzMesh *Mesh, // [in] associated HorzMesh - const VertCoord *VCoord, // [in] associated VertCoord - Config *Options // [in] configuration options + VertAdv(const std::string &Name, ///< [in] Name for new VertAdv + const HorzMesh *Mesh, ///< [in] associated HorzMesh + const VertCoord *VCoord, ///< [in] associated VertCoord + Config *Options ///< [in] configuration options ); + /// define field metadata + void defineFields(); + // Forbid copy and move construction VertAdv(const VertAdv &) = delete; VertAdv(VertAdv &&) = delete; From f811db199f90ddda72b95a6fc91dc130e0f70e94 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Fri, 9 Jan 2026 14:25:52 -0800 Subject: [PATCH 15/42] Some bug fixes for running on GPU --- components/omega/test/ocn/VertAdvTest.cpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/components/omega/test/ocn/VertAdvTest.cpp b/components/omega/test/ocn/VertAdvTest.cpp index 71f7e448f2e7..18b1bea12f9f 100644 --- a/components/omega/test/ocn/VertAdvTest.cpp +++ b/components/omega/test/ocn/VertAdvTest.cpp @@ -186,7 +186,6 @@ int main(int argc, char *argv[]) { I4 NVertLayers = DefVertAdv->NVertLayers; I4 NVertLayersP1 = DefVertAdv->NVertLayersP1; - I4 NTracers = DefVertAdv->NTracers; Array2DReal FluxLayerThickEdge("FluxLayerThickEdge", DefDecomp->NEdgesSize, NVertLayers); @@ -194,7 +193,7 @@ int main(int argc, char *argv[]) { NVertLayers); const I4 ICell0 = 0; - const I4 NEdgesOnCell0 = DefDecomp->NEdgesOnCell(ICell0); + const I4 NEdgesOnCell0 = DefDecomp->NEdgesOnCellH(ICell0); // Test for computeVerticalVelocity I4 Err = 0; @@ -225,9 +224,9 @@ int main(int argc, char *argv[]) { // layer Real Perim0 = 0.; for (int J = 0; J < NEdgesOnCell0; ++J) { - Perim0 += DefMesh->DvEdge(DefMesh->EdgesOnCell(ICell0, J)); + Perim0 += DefMesh->DvEdgeH(DefMesh->EdgesOnCellH(ICell0, J)); } - Real InvAreaCell0 = 1._Real / DefMesh->AreaCell(ICell0); + Real InvAreaCell0 = 1._Real / DefMesh->AreaCellH(ICell0); for (int K = 1; K < NVertLayersP1; ++K) { Real Expected = (NVertLayers - K) * Perim0 * InvAreaCell0; @@ -333,6 +332,7 @@ int main(int argc, char *argv[]) { DefVertAdv->Coef3rdOrder = 1.; Tol = 0.1_Real; + OMEGA_SCOPE(LocMaxLyrCell, DefVCoord->MaxLayerCell); // Compute L2 error for computeStdVAdvTend with each VertFluxOption over a // set of resolutions. Compare successive errors for each Order to confirm @@ -342,8 +342,10 @@ int main(int argc, char *argv[]) { DefVertAdv->VertFluxChoice = Order; std::vector L2Errors; for (I4 NLayers : NLayersArray) { - DefVCoord->MaxLayerCell(0) = NLayers - 1; - Real L2Err = tendTest(NLayers, DefVertAdv); + parallelFor( + {1}, + KOKKOS_LAMBDA(const int &) { LocMaxLyrCell(0) = NLayers - 1; }); + Real L2Err = tendTest(NLayers, DefVertAdv); L2Errors.push_back(L2Err); } Real ExpectedRat = std::pow(2._Real, OrderOfAcc.at(Order)); @@ -370,6 +372,9 @@ int main(int argc, char *argv[]) { DefVertAdv->Coef3rdOrder = 0.25; std::vector TestVel = {1._Real, -1._Real}; std::vector StartIdx = {NVertLayers / 2, NVertLayers / 4}; + parallelFor( + {1}, + KOKKOS_LAMBDA(const int &) { LocMaxLyrCell(0) = NVertLayers - 1; }); // Test a top-hat function for a uniform velocity in both directions for (int IMono = 0; IMono < 2; ++IMono) { @@ -378,8 +383,7 @@ int main(int argc, char *argv[]) { Array2DReal LayerThick("LayerThick", 1, NVertLayers); Array3DReal TendCell3D("TendCell3D", 1, 1, NVertLayers); deepCopy(LayerThick, 1._Real); - DefVCoord->MaxLayerCell(0) = NVertLayers - 1; - DefVertAdv->NVertLayers = NVertLayers; + DefVertAdv->NVertLayers = NVertLayers; LocVertVel = Array2DReal("TotalVerticaVelocity", 1, NVertLayers + 1); deepCopy(LocVertVel, TestVel[IMono]); parallelFor( From 5b40aa4b086f60eaf9d71755045fb2122756309b Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Sat, 10 Jan 2026 15:24:00 -0500 Subject: [PATCH 16/42] Bug fixes: BottomDepth no longer in HorzMesh and typo fixes --- .../omega/src/ocn/CustomTendencyTerms.cpp | 5 +++-- components/omega/src/ocn/HorzMesh.cpp | 10 --------- components/omega/src/ocn/HorzMesh.h | 9 -------- .../auxiliaryVars/LayerThicknessAuxVars.cpp | 4 ++-- components/omega/test/ocn/HorzMeshTest.cpp | 22 ------------------- 5 files changed, 5 insertions(+), 45 deletions(-) diff --git a/components/omega/src/ocn/CustomTendencyTerms.cpp b/components/omega/src/ocn/CustomTendencyTerms.cpp index ff9c17da6435..79ab03630912 100644 --- a/components/omega/src/ocn/CustomTendencyTerms.cpp +++ b/components/omega/src/ocn/CustomTendencyTerms.cpp @@ -10,6 +10,7 @@ #include "Config.h" #include "GlobalConstants.h" #include "TimeStepper.h" +#include "VertCoord.h" namespace OMEGA { @@ -76,8 +77,8 @@ void ManufacturedSolution::init() { /// This test case assumes that the restingThickness is horizontally uniform /// and that only one vertical level is used so only one set of indices is /// used here. - HorzMesh *DefHorzMesh = HorzMesh::getDefault(); - R8 H0 = DefHorzMesh->BottomDepthH(0); + VertCoord *DefVCoord = VertCoord::getDefault(); + R8 H0 = DefVCoord->BottomDepthH(0); // Define and compute common constants R8 Kx = TwoPi / WavelengthX; // Wave in X-dir diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp index 0c626b49905e..0494273df9cf 100644 --- a/components/omega/src/ocn/HorzMesh.cpp +++ b/components/omega/src/ocn/HorzMesh.cpp @@ -111,9 +111,6 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh // Read x/y/z and lon/lat coordinates for cells, edges, and vertices readCoordinates(); - // Read the cell-centered bottom depth - readBottomDepth(); - // Read the mesh areas, lengths, and angles readMeasurements(); @@ -439,12 +436,6 @@ void HorzMesh::readCoordinates() { } // end readCoordinates -//------------------------------------------------------------------------------ -// Read the cell-centered bottom depth -void HorzMesh::readBottomDepth() { - readCellArray(BottomDepthH, "bottomDepth"); -} // end readDepth - //------------------------------------------------------------------------------ // Read the mesh areas (cell, triangle, and kite), // lengths (between centers and vertices), and edge angles @@ -600,7 +591,6 @@ void HorzMesh::copyToDevice() { AngleEdge = createDeviceMirrorCopy(AngleEdgeH); WeightsOnEdge = createDeviceMirrorCopy(WeightsOnEdgeH); FVertex = createDeviceMirrorCopy(FVertexH); - BottomDepth = createDeviceMirrorCopy(BottomDepthH); FEdge = createDeviceMirrorCopy(FEdgeH); XCell = createDeviceMirrorCopy(XCellH); YCell = createDeviceMirrorCopy(YCellH); diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h index 0b244c0e17a9..424c34e3a5b9 100644 --- a/components/omega/src/ocn/HorzMesh.h +++ b/components/omega/src/ocn/HorzMesh.h @@ -42,16 +42,12 @@ class HorzMesh { void readCoordinates(); - void readBottomDepth(); - void readMeasurements(); void readWeights(); void readCoriolis(); - // void computeEdgeSign(); - void copyToDevice(); // int computeMesh(); @@ -232,11 +228,6 @@ class HorzMesh { Array1DReal FVertex; ///< Coriolis parameter at vertices (radians s^-1) HostArray1DReal FVertexH; ///< Coriolis parameter at vertices (radians s^-1) - // Depth - - Array1DReal BottomDepth; ///< Depth of the bottom of the ocean (m) - HostArray1DReal BottomDepthH; ///< Depth of the bottom of the ocean (m) - // Edge sign Array2DReal EdgeSignOnCell; ///< Sign of vector connecting cells diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp index 506a328749e0..ba1a0cdde4a8 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp @@ -16,10 +16,10 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix, VCoord->NVertLayers), ProvThickness("ProvThickness" + AuxStateSuffix, Mesh->NCellsSize, VCoord->NVertLayers), - AreaCell(Mesh->AreaCell), DvEdge(Mesh->AreaCell), + AreaCell(Mesh->AreaCell), DvEdge(Mesh->DvEdge), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), CellsOnEdge(Mesh->CellsOnEdge), - BottomDepth(Mesh->BottomDepth), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + BottomDepth(VCoord->BottomDepth), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp index 0c218325b100..f178b3e8f171 100644 --- a/components/omega/test/ocn/HorzMeshTest.cpp +++ b/components/omega/test/ocn/HorzMeshTest.cpp @@ -368,28 +368,6 @@ int main(int argc, char *argv[]) { LOG_INFO("HorzMeshTest: Vertex lon/lat test PASS"); } - // Test bounds of bathymetry - // Find minimum and maximum values of the bottom depth - // and compares to reasonable values - // Tests that the bottom depth has been read in correctly - R8 MaxBathy = -1e10; - R8 MinBathy = 1e10; - for (int Cell = 0; Cell < LocCells; Cell++) { - if (Mesh->BottomDepthH(Cell) < MinBathy) { - MinBathy = Mesh->BottomDepthH(Cell); - } - if (Mesh->BottomDepthH(Cell) > MaxBathy) { - MaxBathy = Mesh->BottomDepthH(Cell); - } - } - - if ((MinBathy > 0) && (MaxBathy < 11000.0)) { - LOG_INFO("HorzMeshTest: Bathy min/max test PASS"); - } else { - RetVal += 1; - LOG_INFO("HorzMeshTest: Bathy min/max test FAIL"); - } - // Test cell areas // Find the global sum of all the local cell areas // and compares to reasonable value for Earth's ocean area From 9bdf05d844bc3a0efbb28749c00177cd4cc26992 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Sun, 11 Jan 2026 08:37:50 -0800 Subject: [PATCH 17/42] Make host copies of VertAdv arrays --- components/omega/src/ocn/VertAdv.cpp | 30 ++++++++++++++++++++++++++++ components/omega/src/ocn/VertAdv.h | 12 +++++++++++ 2 files changed, 42 insertions(+) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 97b8ab813e67..886b839025e1 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -69,10 +69,16 @@ VertAdv::VertAdv(const std::string &Name_, //< [in] name for new VertAdv Array2DReal("TotalVerticalVelocity", NCellsSize, NVertLayersP1); VertFlux = Array3DReal("VertFlux", NTracers, NCellsSize, NVertLayersP1); + // Allocate host copies + VerticalVelocityH = createHostMirrorCopy(VerticalVelocity); + TotalVerticalVelocityH = createHostMirrorCopy(TotalVerticalVelocity); + VertFluxH = createHostMirrorCopy(VertFlux); + // Low-order flux array only needed for flux-corrected transport if (VertAdvChoice == VertAdvOption::FCT) { LowOrderVertFlux = Array3DReal("LowOrderVertFlux", NTracers, NCellsSize, NVertLayersP1); + LowOrderVertFluxH = createHostMirrorCopy(LowOrderVertFlux); } defineFields(); @@ -232,6 +238,30 @@ void VertAdv::erase(std::string Name) { // process, calls the destructor } // end erase +//------------------------------------------------------------------------------ +// Perform deepCopy for each variable array from device to host +void VertAdv::copyToHost() { + + deepCopy(VerticalVelocityH, VerticalVelocity); + deepCopy(TotalVerticalVelocityH, TotalVerticalVelocity); + deepCopy(VertFluxH, VertFlux); + if (VertAdvChoice == VertAdvOption::FCT) { + deepCopy(LowOrderVertFluxH, LowOrderVertFlux); + } +} + +//------------------------------------------------------------------------------ +// Perform deepCopy for each variable array from host to device +void VertAdv::copyToDevice() { + + deepCopy(VerticalVelocity, VerticalVelocityH); + deepCopy(TotalVerticalVelocity, TotalVerticalVelocityH); + deepCopy(VertFlux, VertFluxH); + if (VertAdvChoice == VertAdvOption::FCT) { + deepCopy(LowOrderVertFlux, LowOrderVertFluxH); + } +} + //------------------------------------------------------------------------------ // Get default VertAdv VertAdv *VertAdv::getDefault() { return VertAdv::DefaultVertAdv; } diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index 5a173fdd7fb7..95bd30968913 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -68,6 +68,12 @@ class VertAdv { Array3DReal VertFlux; ///< fluxes at vertical interfaces Array3DReal LowOrderVertFlux; ///< low-order fluxes for FCT + // Arrays on host + HostArray2DReal VerticalVelocityH; + HostArray2DReal TotalVerticalVelocityH; + HostArray3DReal VertFluxH; + HostArray3DReal LowOrderVertFluxH; + // VertAdv instance name and group name std::string Name; std::string GroupName; @@ -92,6 +98,12 @@ class VertAdv { Config *Options ///< [in] configuration options ); + /// Copy member arrays from device to host + void copyToHost(); + + /// Copy member arrays from host to device + void copyToDevice(); + /// Destructor - deallocates all memory and deletes a VertAdv ~VertAdv(); From 7f8f0efff69a0f492b50c99b5a83146fd7100a93 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Sun, 11 Jan 2026 09:46:39 -0800 Subject: [PATCH 18/42] Linting fix --- .../omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp index ba1a0cdde4a8..8a69aa7c8269 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.cpp @@ -19,7 +19,8 @@ LayerThicknessAuxVars::LayerThicknessAuxVars(const std::string &AuxStateSuffix, AreaCell(Mesh->AreaCell), DvEdge(Mesh->DvEdge), NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell), EdgeSignOnCell(Mesh->EdgeSignOnCell), CellsOnEdge(Mesh->CellsOnEdge), - BottomDepth(VCoord->BottomDepth), MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + BottomDepth(VCoord->BottomDepth), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop), MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} From d6257968366c3707a8c97e5773c2786f039e2583 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Tue, 3 Feb 2026 12:01:18 -0800 Subject: [PATCH 19/42] Resize flux arrays in VertAdvTest --- components/omega/test/ocn/VertAdvTest.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/components/omega/test/ocn/VertAdvTest.cpp b/components/omega/test/ocn/VertAdvTest.cpp index 18b1bea12f9f..816ea9183192 100644 --- a/components/omega/test/ocn/VertAdvTest.cpp +++ b/components/omega/test/ocn/VertAdvTest.cpp @@ -53,6 +53,7 @@ Real tendTest(const int NLayers, VertAdv *VAdv) { VAdv->NVertLayers = NLayers; VAdv->TotalVerticalVelocity = Array2DReal("TotalVerticaVelocity", 1, NLayers + 1); + VAdv->VertFlux = Array3DReal("VertFlux", 1, 1, NLayers + 1); Array1DReal XLayer("XLayer", NLayers); OMEGA_SCOPE(LocTotVertVel, VAdv->TotalVerticalVelocity); Array2DReal LayerThick("LayerThickness", 1, NLayers); @@ -366,7 +367,9 @@ int main(int argc, char *argv[]) { // Monotonicity check for computeFCTVAdvTend - NVertLayers = 256; + NVertLayers = 256; + DefVertAdv->LowOrderVertFlux = + Array3DReal("LowOrderVertFlux", 1, 1, NVertLayers + 1); DefVertAdv->VertAdvChoice = VertAdvOption::FCT; DefVertAdv->VertFluxChoice = VertFluxOption::Third; DefVertAdv->Coef3rdOrder = 0.25; From ec916bfdc6ab3982e92597adef2d78c04dd6dde8 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 5 Feb 2026 08:32:23 -0800 Subject: [PATCH 20/42] Rename config option from Enabled to Enable --- components/omega/configs/Default.yml | 2 +- components/omega/src/ocn/VertAdv.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 7718ac3675e3..2f19e41c06dd 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -26,7 +26,7 @@ Omega: Coef3rdOrder: 0.25 FluxThicknessType: Center FluxTracerType: Center - VerticalTracerFluxLimiterEnabled: true + VerticalTracerFluxLimiterEnable: true VerticalTracerFluxOrder: 3 WindStress: InterpType: Isotropic diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 886b839025e1..72bad626b94f 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -317,7 +317,7 @@ void VertAdv::readConfigOptions(Config *OmegaConfig) { CHECK_ERROR_ABORT(Err, "VertAdv: Advection group not in Config"); bool FluxLimiterOn; - Err += AdvectConfig.get("VerticalTracerFluxLimiterEnabled", FluxLimiterOn); + Err += AdvectConfig.get("VerticalTracerFluxLimiterEnable", FluxLimiterOn); CHECK_ERROR_ABORT( Err, "VertAdv: VerticalTracerFluxLimiterEnabled not found in AdvectConfig"); From 2ea99f182a363ab8bbc2aca73ebda2de5c90a4bb Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 5 Feb 2026 08:52:38 -0800 Subject: [PATCH 21/42] Remove unnecessary header --- components/omega/src/ocn/VertAdv.h | 1 - 1 file changed, 1 deletion(-) diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index 95bd30968913..cae279ffe1fe 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -11,7 +11,6 @@ // //===----------------------------------------------------------------------===// -#include "AuxiliaryState.h" #include "Config.h" #include "DataTypes.h" #include "Error.h" From 60bc8bf6b9d810fabfb793abdcc94b6e762a156a Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 5 Feb 2026 08:56:24 -0800 Subject: [PATCH 22/42] Skip vertical tracer advection when not enabled --- components/omega/src/ocn/VertAdv.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 72bad626b94f..d97f9da6bdcb 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -568,6 +568,10 @@ void VertAdv::computeTracerVAdvTend( const TimeInterval TimeStep //< [in] (optional) time step ) { + // Return if vertical advection tracer tendency not enabled + if (!TracerVertAdvEnabled) + return; + // Compute tracer fluxes at the interfaces computeVerticalFluxes(Tracers, LayerThickness); From f777c610d82c18919ae57d8714e39d4470a1e9bf Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Tue, 10 Feb 2026 15:33:20 -0800 Subject: [PATCH 23/42] Add TimeStep to AuxState for ProvThickness --- components/omega/src/ocn/AuxiliaryState.cpp | 20 +++++++++++++------ components/omega/src/ocn/AuxiliaryState.h | 7 +++++-- components/omega/src/ocn/VertAdv.cpp | 2 +- .../omega/test/ocn/AuxiliaryStateTest.cpp | 3 ++- .../test/timeStepping/TimeStepperTest.cpp | 4 +++- 5 files changed, 25 insertions(+), 11 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index ef3cea652d64..2b33505507d5 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -3,6 +3,7 @@ #include "Field.h" #include "Logging.h" #include "Pacer.h" +#include "TimeStepper.h" namespace OMEGA { @@ -19,7 +20,7 @@ static std::string stripDefault(const std::string &Name) { // fields with IOStreams AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, const VertCoord *VCoord, - int NTracers) + int NTracers, TimeInterval TimeStep) : Mesh(Mesh), MeshHalo(MeshHalo), VCoord(VCoord), Name(stripDefault(Name)), KineticAux(stripDefault(Name), Mesh, VCoord), LayerThicknessAux(stripDefault(Name), Mesh, VCoord), @@ -82,6 +83,9 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, OMEGA_SCOPE(MaxLayerEdgeBot, VCoord->MaxLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + R8 TimeStepSeconds; + TimeStep.get(TimeStepSeconds, TimeUnits::Seconds); + Pacer::start("AuxState:computeMomAux", 1); Pacer::start("AuxState:vertexAuxState1", 2); @@ -197,8 +201,8 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { LocLayerThicknessAux.computeVarsOnCells( - ICell, KChunk, LayerThickCell, NormalVelEdge, 0._Real); - // TODO: make timestep available to this call + ICell, KChunk, LayerThickCell, NormalVelEdge, + TimeStepSeconds); }); }); Pacer::stop("AuxState:cellAuxState3", 2); @@ -274,7 +278,8 @@ void AuxiliaryState::computeAll(const OceanState *State, AuxiliaryState *AuxiliaryState::create(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, const VertCoord *VCoord, - const int NTracers) { + const int NTracers, + TimeInterval TimeStep) { if (AllAuxStates.find(Name) != AllAuxStates.end()) { LOG_ERROR("Attempted to create a new AuxiliaryState with name {} but it " "already exists", @@ -283,7 +288,7 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name, } auto *NewAuxState = - new AuxiliaryState(Name, Mesh, MeshHalo, VCoord, NTracers); + new AuxiliaryState(Name, Mesh, MeshHalo, VCoord, NTracers, TimeStep); AllAuxStates.emplace(Name, NewAuxState); return NewAuxState; @@ -295,11 +300,14 @@ void AuxiliaryState::init() { const HorzMesh *DefMesh = HorzMesh::getDefault(); Halo *DefHalo = Halo::getDefault(); const VertCoord *DefVCoord = VertCoord::getDefault(); + const TimeStepper *DefTimeStepper = TimeStepper::getDefault(); int NTracers = Tracers::getNumTracers(); + TimeInterval TimeStep = DefTimeStepper->getTimeStep(); AuxiliaryState::DefaultAuxState = - AuxiliaryState::create("Default", DefMesh, DefHalo, DefVCoord, NTracers); + AuxiliaryState::create("Default", DefMesh, DefHalo, DefVCoord, NTracers, + TimeStep); Config *OmegaConfig = Config::getOmegaConfig(); DefaultAuxState->readConfigOptions(OmegaConfig); diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 313cc157f0f3..5b997d93e3c7 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -6,6 +6,7 @@ #include "Halo.h" #include "HorzMesh.h" #include "OceanState.h" +#include "TimeMgr.h" #include "Tracers.h" #include "VertCoord.h" #include "auxiliaryVars/KineticAuxVars.h" @@ -51,7 +52,7 @@ class AuxiliaryState { // Create a non-default auxiliary state static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, const VertCoord *VCoord, - int NTracers); + int NTracers, TimeInterval TimeStep); /// Get the default auxiliary state static AuxiliaryState *getDefault(); @@ -84,7 +85,7 @@ class AuxiliaryState { private: AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, - const VertCoord *VCoord, int NTracers); + const VertCoord *VCoord, int NTracers, TimeInterval TimeStep); AuxiliaryState(const AuxiliaryState &) = delete; AuxiliaryState(AuxiliaryState &&) = delete; @@ -92,6 +93,8 @@ class AuxiliaryState { const HorzMesh *Mesh; Halo *MeshHalo; const VertCoord *VCoord; + TimeInterval TimeStep; + static AuxiliaryState *DefaultAuxState; static std::map> AllAuxStates; }; diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index d97f9da6bdcb..c59b77505ab6 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -320,7 +320,7 @@ void VertAdv::readConfigOptions(Config *OmegaConfig) { Err += AdvectConfig.get("VerticalTracerFluxLimiterEnable", FluxLimiterOn); CHECK_ERROR_ABORT( Err, - "VertAdv: VerticalTracerFluxLimiterEnabled not found in AdvectConfig"); + "VertAdv: VerticalTracerFluxLimiterEnable not found in AdvectConfig"); if (FluxLimiterOn) { VertAdvChoice = VertAdvOption::FCT; } else { diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index 543f074d5f24..ea77c74c6e6b 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -150,8 +150,9 @@ int testAuxState() { const auto *Mesh = HorzMesh::getDefault(); auto *MeshHalo = Halo::getDefault(); const auto *VCoord = VertCoord::getDefault(); + TimeInterval TimeStep; // test creation of another auxiliary state - AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, VCoord, 3); + AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, VCoord, 3, TimeStep); // test retrievel of another if (AuxiliaryState::get("AnotherAuxState")) { diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 8188bac63f60..f286fe83ba90 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -204,8 +204,10 @@ int initTimeStepperTest(const std::string &mesh) { LOG_ERROR("TimeStepperTest: error creating test state"); } + TimeInterval ZeroTimeStep; auto *TestAuxState = AuxiliaryState::create("TestAuxState", DefMesh, DefHalo, - DefVertCoord, NTracers); + DefVertCoord, NTracers, + ZeroTimeStep); Config *OmegaConfig = Config::getOmegaConfig(); TestAuxState->readConfigOptions(OmegaConfig); From f7fe5964c1975f315df508ba9efff8b034b7e6d8 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Tue, 10 Feb 2026 17:33:42 -0800 Subject: [PATCH 24/42] Remove commented out ssh calculation --- .../omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h | 9 --------- 1 file changed, 9 deletions(-) diff --git a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h index a415bfa02ba3..68bb7da2296c 100644 --- a/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/LayerThicknessAuxVars.h @@ -96,15 +96,6 @@ class LayerThicknessAuxVars { const int K = KStart + KVec; ProvThickness(ICell, K) = LayerThickCell(ICell, K) + TmpProv[KVec]; } - - /* - Real TotalThickness = 0.0; - for (int K = 0; K < NVertLayers; K++) { - TotalThickness += LayerThickCell(ICell, K); - } - - SshCell(ICell) = TotalThickness - BottomDepth(ICell); - */ } void registerFields(const std::string &AuxGroupName, From 91d6523e9ad652e5cef45bbeb5d886a2b0ba1e5b Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Tue, 10 Feb 2026 20:19:06 -0800 Subject: [PATCH 25/42] Add LayerThicknessAux.computeVarsOnCells to AuxState::computeAll This PR adds ProvThickness to LayerThicknessAux, which is needed by the tracer equations. PR #335 removes SSH from LayerThicknessAux and removes computeVarsOnCells from computeMomAux. --- components/omega/src/ocn/AuxiliaryState.cpp | 34 ++++++++++++++++----- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 2b33505507d5..4807232af675 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -221,16 +221,37 @@ void AuxiliaryState::computeAll(const OceanState *State, const int NTracers = TracerArray.extent_int(0); + OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); OMEGA_SCOPE(LocTracerAux, TracerAux); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + R8 TimeStepSeconds; + TimeStep.get(TimeStepSeconds, TimeUnits::Seconds); + Pacer::start("AuxState:computeAll", 1); computeMomAux(State, ThickTimeLevel, VelTimeLevel); + 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:edgeAuxState4", 2); parallelForOuter( "edgeAuxState4", {NTracers, Mesh->NEdgesAll}, @@ -297,17 +318,16 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name, // Create the default auxiliary state. Assumes that HorzMesh, VertCoord and // Halo have been initialized. void AuxiliaryState::init() { - const HorzMesh *DefMesh = HorzMesh::getDefault(); - Halo *DefHalo = Halo::getDefault(); - const VertCoord *DefVCoord = VertCoord::getDefault(); + const HorzMesh *DefMesh = HorzMesh::getDefault(); + Halo *DefHalo = Halo::getDefault(); + const VertCoord *DefVCoord = VertCoord::getDefault(); const TimeStepper *DefTimeStepper = TimeStepper::getDefault(); - int NTracers = Tracers::getNumTracers(); + int NTracers = Tracers::getNumTracers(); TimeInterval TimeStep = DefTimeStepper->getTimeStep(); - AuxiliaryState::DefaultAuxState = - AuxiliaryState::create("Default", DefMesh, DefHalo, DefVCoord, NTracers, - TimeStep); + AuxiliaryState::DefaultAuxState = AuxiliaryState::create( + "Default", DefMesh, DefHalo, DefVCoord, NTracers, TimeStep); Config *OmegaConfig = Config::getOmegaConfig(); DefaultAuxState->readConfigOptions(OmegaConfig); From 3d3f1a37c71301b98651989d093e4078417ce078 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 07:23:23 -0800 Subject: [PATCH 26/42] Add VertAdv to Tendencies --- components/omega/src/ocn/Tendencies.cpp | 14 +++++++++----- components/omega/src/ocn/Tendencies.h | 5 ++++- components/omega/test/ocn/TendenciesTest.cpp | 4 +++- .../omega/test/timeStepping/TimeStepperTest.cpp | 11 +++++++---- 4 files changed, 23 insertions(+), 11 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 7cf8a7a401a8..418e9b0e5ab0 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -13,6 +13,7 @@ #include "Error.h" #include "Pacer.h" #include "Tracers.h" +#include "VertAdv.h" namespace OMEGA { @@ -27,6 +28,7 @@ void Tendencies::init() { HorzMesh *DefHorzMesh = HorzMesh::getDefault(); VertCoord *DefVertCoord = VertCoord::getDefault(); + VertAdv *DefVertAdv = VertAdv::getDefault(); I4 NTracers = Tracers::getNumTracers(); @@ -65,8 +67,8 @@ void Tendencies::init() { // Ceate default tendencies Tendencies::DefaultTendencies = - create("Default", DefHorzMesh, DefVertCoord, NTracers, &TendConfig, - CustomThickTend, CustomVelTend); + create("Default", DefHorzMesh, DefVertCoord, DefVertAdv, NTracers, + &TendConfig, CustomThickTend, CustomVelTend); DefaultTendencies->readTendConfig(&TendConfig); @@ -214,11 +216,12 @@ void Tendencies::readTendConfig( Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies const HorzMesh *Mesh, ///< [in] Horizontal mesh const VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers Config *Options, ///< [in] Configuration options CustomTendencyType InCustomThicknessTend, CustomTendencyType InCustomVelocityTend) - : Mesh(Mesh), VCoord(VCoord), ThicknessFluxDiv(Mesh, VCoord), + : Mesh(Mesh), VCoord(VCoord), VAdv(VAdv), ThicknessFluxDiv(Mesh, VCoord), PotientialVortHAdv(Mesh, VCoord), KEGrad(Mesh, VCoord), SSHGrad(Mesh, VCoord), VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord), WindForcing(Mesh, VCoord), @@ -242,10 +245,11 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies const HorzMesh *Mesh, ///< [in] Horizontal mesh const VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers Config *Options) ///< [in] Configuration options - : Tendencies(Name, Mesh, VCoord, NTracersIn, Options, CustomTendencyType{}, - CustomTendencyType{}) {} + : Tendencies(Name, Mesh, VCoord, VAdv, NTracersIn, Options, + CustomTendencyType{}, CustomTendencyType{}) {} //------------------------------------------------------------------------------ // Compute tendencies for layer thickness equation diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 851ce1fc646c..c8c0d84b5348 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -37,7 +37,7 @@ #include "OceanState.h" #include "TendencyTerms.h" #include "TimeMgr.h" -#include "VertCoord.h" +#include "VertAdv.h" #include #include @@ -151,6 +151,7 @@ class Tendencies { Tendencies(const std::string &Name, ///< [in] Name for tendencies const HorzMesh *Mesh, ///< [in] Horizontal mesh const VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers Config *Options, ///< [in] Configuration options CustomTendencyType InCustomThicknessTend, @@ -159,6 +160,7 @@ class Tendencies { Tendencies(const std::string &Name, ///< [in] Name for tendencies const HorzMesh *Mesh, ///< [in] Horizontal mesh const VertCoord *VCoord, ///< [in] Vertical coordinate + VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers Config *Options ///< [in] Configuration options ); @@ -169,6 +171,7 @@ class Tendencies { const HorzMesh *Mesh; ///< Pointer to horizontal mesh const VertCoord *VCoord; ///< Pointer to vertical coordinate + VertAdv *VAdv; ///< Pointer to vertical advection I4 NTracers; ///< Number of tracers // Pointer to default tendencies diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 4595bfea6d88..38533c112579 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -126,6 +126,7 @@ int initTendenciesTest(const std::string &mesh) { HorzMesh::init(); VertCoord::init(); Tracers::init(); + VertAdv::init(); // VertCoord::init2(); @@ -159,10 +160,11 @@ int testTendencies() { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); + const auto VAdv = VertAdv::getDefault(); VCoord->NVertLayers = 12; // test creation of another tendencies Config *Options = Config::getOmegaConfig(); - Tendencies::create("TestTendencies", Mesh, VCoord, 3, Options); + Tendencies::create("TestTendencies", Mesh, VCoord, VAdv, 3, Options); // test retrievel of another tendencies if (Tendencies::get("TestTendencies")) { diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index f286fe83ba90..b30419faf180 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -177,6 +177,10 @@ int initTimeStepperTest(const std::string &mesh) { auto *DefVertCoord = VertCoord::getDefault(); Tracers::init(); + + VertAdv::init(); + auto *DefVertAdv = VertAdv::getDefault(); + AuxiliaryState::init(); Tendencies::init(); @@ -205,9 +209,8 @@ int initTimeStepperTest(const std::string &mesh) { } TimeInterval ZeroTimeStep; - auto *TestAuxState = AuxiliaryState::create("TestAuxState", DefMesh, DefHalo, - DefVertCoord, NTracers, - ZeroTimeStep); + auto *TestAuxState = AuxiliaryState::create( + "TestAuxState", DefMesh, DefHalo, DefVertCoord, NTracers, ZeroTimeStep); Config *OmegaConfig = Config::getOmegaConfig(); TestAuxState->readConfigOptions(OmegaConfig); @@ -221,7 +224,7 @@ int initTimeStepperTest(const std::string &mesh) { // Creating non-default tendencies with custom velocity tendencies auto *TestTendencies = Tendencies::create( - "TestTendencies", DefMesh, DefVertCoord, NTracers, &Options, + "TestTendencies", DefMesh, DefVertCoord, DefVertAdv, NTracers, &Options, Tendencies::CustomTendencyType{}, DecayVelocityTendency{}); if (!TestTendencies) { Err++; From 9616d8f7f720904eb70f00cbeda364c848f3d8a6 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 07:36:27 -0800 Subject: [PATCH 27/42] Disable VertAdv tends in TimeStepperTest --- components/omega/test/timeStepping/TimeStepperTest.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index b30419faf180..6b530eb491f8 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -243,6 +243,9 @@ int initTimeStepperTest(const std::string &mesh) { TestTendencies->TracerHyperDiff.Enabled = false; TestTendencies->WindForcing.Enabled = false; TestTendencies->BottomDrag.Enabled = false; + DefVertAdv->ThickVertAdvEnabled = false; + DefVertAdv->VelVertAdvEnabled = false; + DefVertAdv->TracerVertAdvEnabled = false; return Err; } @@ -283,6 +286,7 @@ void finalizeTimeStepperTest() { Tendencies::clear(); AuxiliaryState::clear(); OceanState::clear(); + VertAdv::clear(); VertCoord::clear(); HorzMesh::clear(); Dimension::clear(); From 9a5ef3abe577a0043b767e7b367eeaf387c6d69b Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 08:39:03 -0800 Subject: [PATCH 28/42] Add VertAdv to AuxState --- components/omega/src/ocn/AuxiliaryState.cpp | 20 ++++++++++--------- components/omega/src/ocn/AuxiliaryState.h | 8 ++++++-- .../omega/test/ocn/AuxiliaryStateTest.cpp | 7 ++++++- components/omega/test/ocn/TendenciesTest.cpp | 1 + .../test/timeStepping/TimeStepperTest.cpp | 17 ++++++++-------- 5 files changed, 33 insertions(+), 20 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 4807232af675..7330315248fb 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -20,9 +20,10 @@ static std::string stripDefault(const std::string &Name) { // fields with IOStreams AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, const VertCoord *VCoord, - int NTracers, TimeInterval TimeStep) - : Mesh(Mesh), MeshHalo(MeshHalo), VCoord(VCoord), Name(stripDefault(Name)), - KineticAux(stripDefault(Name), Mesh, VCoord), + VertAdv *VAdv, int NTracers, + TimeInterval TimeStep) + : Mesh(Mesh), MeshHalo(MeshHalo), VCoord(VCoord), VAdv(VAdv), + Name(stripDefault(Name)), KineticAux(stripDefault(Name), Mesh, VCoord), LayerThicknessAux(stripDefault(Name), Mesh, VCoord), VorticityAux(stripDefault(Name), Mesh, VCoord), VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), @@ -298,7 +299,7 @@ void AuxiliaryState::computeAll(const OceanState *State, // Create a non-default auxiliary state AuxiliaryState *AuxiliaryState::create(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, - const VertCoord *VCoord, + const VertCoord *VCoord, VertAdv *VAdv, const int NTracers, TimeInterval TimeStep) { if (AllAuxStates.find(Name) != AllAuxStates.end()) { @@ -308,26 +309,27 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name, return nullptr; } - auto *NewAuxState = - new AuxiliaryState(Name, Mesh, MeshHalo, VCoord, NTracers, TimeStep); + auto *NewAuxState = new AuxiliaryState(Name, Mesh, MeshHalo, VCoord, VAdv, + NTracers, TimeStep); AllAuxStates.emplace(Name, NewAuxState); return NewAuxState; } -// Create the default auxiliary state. Assumes that HorzMesh, VertCoord and -// Halo have been initialized. +// Create the default auxiliary state. Assumes that HorzMesh, VertCoord, +// VertAdv, and Halo have been initialized. void AuxiliaryState::init() { const HorzMesh *DefMesh = HorzMesh::getDefault(); Halo *DefHalo = Halo::getDefault(); const VertCoord *DefVCoord = VertCoord::getDefault(); + VertAdv *DefVAdv = VertAdv::getDefault(); const TimeStepper *DefTimeStepper = TimeStepper::getDefault(); int NTracers = Tracers::getNumTracers(); TimeInterval TimeStep = DefTimeStepper->getTimeStep(); AuxiliaryState::DefaultAuxState = AuxiliaryState::create( - "Default", DefMesh, DefHalo, DefVCoord, NTracers, TimeStep); + "Default", DefMesh, DefHalo, DefVCoord, DefVAdv, NTracers, TimeStep); Config *OmegaConfig = Config::getOmegaConfig(); DefaultAuxState->readConfigOptions(OmegaConfig); diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 5b997d93e3c7..60a414d39d35 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -8,6 +8,7 @@ #include "OceanState.h" #include "TimeMgr.h" #include "Tracers.h" +#include "VertAdv.h" #include "VertCoord.h" #include "auxiliaryVars/KineticAuxVars.h" #include "auxiliaryVars/LayerThicknessAuxVars.h" @@ -52,7 +53,8 @@ class AuxiliaryState { // Create a non-default auxiliary state static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, const VertCoord *VCoord, - int NTracers, TimeInterval TimeStep); + VertAdv *Vadv, int NTracers, + TimeInterval TimeStep); /// Get the default auxiliary state static AuxiliaryState *getDefault(); @@ -85,7 +87,8 @@ class AuxiliaryState { private: AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, - const VertCoord *VCoord, int NTracers, TimeInterval TimeStep); + const VertCoord *VCoord, VertAdv *Vadv, int NTracers, + TimeInterval TimeStep); AuxiliaryState(const AuxiliaryState &) = delete; AuxiliaryState(AuxiliaryState &&) = delete; @@ -93,6 +96,7 @@ class AuxiliaryState { const HorzMesh *Mesh; Halo *MeshHalo; const VertCoord *VCoord; + VertAdv *VAdv; TimeInterval TimeStep; static AuxiliaryState *DefaultAuxState; diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp index ea77c74c6e6b..9d42cd7bc589 100644 --- a/components/omega/test/ocn/AuxiliaryStateTest.cpp +++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp @@ -127,6 +127,8 @@ int initAuxStateTest(const std::string &mesh) { LOG_ERROR("AuxStateTest: error initializing default state"); } + VertAdv::init(); + return Err; } @@ -150,9 +152,11 @@ int testAuxState() { const auto *Mesh = HorzMesh::getDefault(); auto *MeshHalo = Halo::getDefault(); const auto *VCoord = VertCoord::getDefault(); + auto *VAdv = VertAdv::getDefault(); TimeInterval TimeStep; // test creation of another auxiliary state - AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, VCoord, 3, TimeStep); + AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, VCoord, VAdv, 3, + TimeStep); // test retrievel of another if (AuxiliaryState::get("AnotherAuxState")) { @@ -324,6 +328,7 @@ int testAuxState() { void finalizeAuxStateTest() { Tracers::clear(); OceanState::clear(); + VertAdv::clear(); VertCoord::clear(); HorzMesh::clear(); Field::clear(); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 38533c112579..63177b8721b3 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -240,6 +240,7 @@ void finalizeTendenciesTest() { Tracers::clear(); AuxiliaryState::clear(); OceanState::clear(); + VertAdv::clear(); VertCoord::clear(); HorzMesh::clear(); Field::clear(); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 6b530eb491f8..a0c71c91ee1b 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -179,7 +179,7 @@ int initTimeStepperTest(const std::string &mesh) { Tracers::init(); VertAdv::init(); - auto *DefVertAdv = VertAdv::getDefault(); + auto *DefVAdv = VertAdv::getDefault(); AuxiliaryState::init(); Tendencies::init(); @@ -208,9 +208,10 @@ int initTimeStepperTest(const std::string &mesh) { LOG_ERROR("TimeStepperTest: error creating test state"); } - TimeInterval ZeroTimeStep; - auto *TestAuxState = AuxiliaryState::create( - "TestAuxState", DefMesh, DefHalo, DefVertCoord, NTracers, ZeroTimeStep); + TimeInterval ZeroTimeStep; // Zero-length time step placeholder + auto *TestAuxState = + AuxiliaryState::create("TestAuxState", DefMesh, DefHalo, DefVertCoord, + DefVAdv, NTracers, ZeroTimeStep); Config *OmegaConfig = Config::getOmegaConfig(); TestAuxState->readConfigOptions(OmegaConfig); @@ -224,7 +225,7 @@ int initTimeStepperTest(const std::string &mesh) { // Creating non-default tendencies with custom velocity tendencies auto *TestTendencies = Tendencies::create( - "TestTendencies", DefMesh, DefVertCoord, DefVertAdv, NTracers, &Options, + "TestTendencies", DefMesh, DefVertCoord, DefVAdv, NTracers, &Options, Tendencies::CustomTendencyType{}, DecayVelocityTendency{}); if (!TestTendencies) { Err++; @@ -243,9 +244,9 @@ int initTimeStepperTest(const std::string &mesh) { TestTendencies->TracerHyperDiff.Enabled = false; TestTendencies->WindForcing.Enabled = false; TestTendencies->BottomDrag.Enabled = false; - DefVertAdv->ThickVertAdvEnabled = false; - DefVertAdv->VelVertAdvEnabled = false; - DefVertAdv->TracerVertAdvEnabled = false; + DefVAdv->ThickVertAdvEnabled = false; + DefVAdv->VelVertAdvEnabled = false; + DefVAdv->TracerVertAdvEnabled = false; return Err; } From 3c06217e6807dfa9006bd786d98494b77c0cd3e9 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 11:32:13 -0800 Subject: [PATCH 29/42] Add computeVerticalVelocity to computeMomAux --- components/omega/src/ocn/AuxiliaryState.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 7330315248fb..1d55e1c96eac 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -208,6 +208,13 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, }); 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); } From d2359cefda65523c5ff38798c33227f8864d8548 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 11:34:09 -0800 Subject: [PATCH 30/42] Add VertAdv init and clear to OceanDriver methods --- components/omega/src/ocn/OceanFinal.cpp | 1 + components/omega/src/ocn/OceanInit.cpp | 2 ++ 2 files changed, 3 insertions(+) diff --git a/components/omega/src/ocn/OceanFinal.cpp b/components/omega/src/ocn/OceanFinal.cpp index bcadc327a324..3715e95baef6 100644 --- a/components/omega/src/ocn/OceanFinal.cpp +++ b/components/omega/src/ocn/OceanFinal.cpp @@ -36,6 +36,7 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time Tendencies::clear(); AuxiliaryState::clear(); OceanState::clear(); + VertAdv::clear(); VertCoord::clear(); Dimension::clear(); Field::clear(); diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index 884a9cc51d7b..150d7def405b 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -26,6 +26,7 @@ #include "TimeMgr.h" #include "TimeStepper.h" #include "Tracers.h" +#include "VertAdv.h" #include "VertCoord.h" #include "mpi.h" @@ -131,6 +132,7 @@ int initOmegaModules(MPI_Comm Comm) { HorzMesh::init(); VertCoord::init(); Tracers::init(); + VertAdv::init(); AuxiliaryState::init(); Tendencies::init(); TimeStepper::init2(); From 1a2cb4093f040e6031b4c764d7bbb01c3f0acf5f Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 11:59:28 -0800 Subject: [PATCH 31/42] Add computeThicknessVAdvTend to computeThicknessTendenciesOnly --- components/omega/src/ocn/Tendencies.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 418e9b0e5ab0..711a55cc91df 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -305,6 +305,11 @@ void Tendencies::computeThicknessTendenciesOnly( Pacer::stop("Tend:thicknessFluxDiv", 2); } + Pacer::start("Tend:computeThicknessVAdvTend", 2); + // Compute thickness tendency from vertical advection + VAdv->computeThicknessVAdvTend(LayerThicknessTend); + Pacer::stop("Tend:computeThicknessVAdvTend", 2); + if (CustomThicknessTend) { Pacer::start("Tend:customThicknessTend", 2); CustomThicknessTend(LocLayerThicknessTend, State, AuxState, From f4beb4efc9f154b73749031c9ff9c1f353857f5d Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 12:58:13 -0800 Subject: [PATCH 32/42] Add computeVelocityVAdvTend to computeVelocityTendenciesOnly --- components/omega/src/ocn/Tendencies.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 711a55cc91df..da858d42cabf 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -457,6 +457,12 @@ void Tendencies::computeVelocityTendenciesOnly( Pacer::stop("Tend:velocityHyperDiff", 2); } + Pacer::start("Tend:computeVelocityVAdvTend", 2); + // Compute velocity tendency from vertical advection + VAdv->computeVelocityVAdvTend(NormalVelocityTend, NormalVelEdge, + FluxLayerThickEdge); + Pacer::stop("Tend:computeVelocityVAdvTend", 2); + // Compute wind forcing const auto &NormalStressEdge = AuxState->WindForcingAux.NormalStressEdge; const auto &MeanLayerThickEdge = From a1a84d79bf3d160d56ccc92bb29917a96ae77a6e Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 15:55:39 -0800 Subject: [PATCH 33/42] Add computeTracerVAdvTend to computeTracerTendenciesOnly Also add TimeStep to Tendencies --- components/omega/src/ocn/Tendencies.cpp | 33 +++++++++++++++---- components/omega/src/ocn/Tendencies.h | 3 ++ components/omega/test/ocn/TendenciesTest.cpp | 5 ++- .../test/timeStepping/TimeStepperTest.cpp | 4 +-- 4 files changed, 35 insertions(+), 10 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index da858d42cabf..0ff4242484b4 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -12,6 +12,7 @@ #include "CustomTendencyTerms.h" #include "Error.h" #include "Pacer.h" +#include "TimeStepper.h" #include "Tracers.h" #include "VertAdv.h" @@ -21,14 +22,15 @@ Tendencies *Tendencies::DefaultTendencies = nullptr; std::map> Tendencies::AllTendencies; //------------------------------------------------------------------------------ -// Initialize the tendencies. Assumes that HorzMesh and VertCoord has alread -// been initialized. +// Initialize the tendencies. Assumes that HorzMesh, VertCoord, VertAdv, and +// TimeStepper has already been initialized. void Tendencies::init() { Error Err; // error code - HorzMesh *DefHorzMesh = HorzMesh::getDefault(); - VertCoord *DefVertCoord = VertCoord::getDefault(); - VertAdv *DefVertAdv = VertAdv::getDefault(); + HorzMesh *DefHorzMesh = HorzMesh::getDefault(); + VertCoord *DefVertCoord = VertCoord::getDefault(); + VertAdv *DefVertAdv = VertAdv::getDefault(); + TimeStepper *DefTimeStepper = TimeStepper::getDefault(); I4 NTracers = Tracers::getNumTracers(); @@ -65,10 +67,12 @@ void Tendencies::init() { } // end if UseCustomTendency + TimeInterval TimeStep = DefTimeStepper->getTimeStep(); + // Ceate default tendencies Tendencies::DefaultTendencies = create("Default", DefHorzMesh, DefVertCoord, DefVertAdv, NTracers, - &TendConfig, CustomThickTend, CustomVelTend); + TimeStep, &TendConfig, CustomThickTend, CustomVelTend); DefaultTendencies->readTendConfig(&TendConfig); @@ -218,6 +222,7 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies const VertCoord *VCoord, ///< [in] Vertical coordinate VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStepIn, ///< [in] Time step Config *Options, ///< [in] Configuration options CustomTendencyType InCustomThicknessTend, CustomTendencyType InCustomVelocityTend) @@ -239,6 +244,7 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies VCoord->NVertLayers); NTracers = NTracersIn; + TimeStep = TimeStepIn; } // end constructor @@ -247,8 +253,9 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies const VertCoord *VCoord, ///< [in] Vertical coordinate VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStepIn, ///< [in] Time step Config *Options) ///< [in] Configuration options - : Tendencies(Name, Mesh, VCoord, VAdv, NTracersIn, Options, + : Tendencies(Name, Mesh, VCoord, VAdv, NTracersIn, TimeStepIn, Options, CustomTendencyType{}, CustomTendencyType{}) {} //------------------------------------------------------------------------------ @@ -597,6 +604,18 @@ void Tendencies::computeTracerTendenciesOnly( Pacer::stop("Tend:tracerHyperDiff", 2); } + Pacer::start("Tend:computeTracerVAdvTend", 2); + // compute tracer tendencies from vertical advection + Array2DReal ThicknessForVAdv; + if (VAdv->VertAdvChoice == VertAdvOption::Standard) { + State->getLayerThickness(ThicknessForVAdv, ThickTimeLevel); + } else if (VAdv->VertAdvChoice == VertAdvOption::FCT) { + ThicknessForVAdv = AuxState->LayerThicknessAux.ProvThickness; + } + VAdv->computeTracerVAdvTend(TracerTend, TracerArray, ThicknessForVAdv, + TimeStep); + Pacer::stop("Tend:computeTracerVAdvTend", 2); + Pacer::stop("Tend:computeTracerTendenciesOnly", 1); } // end tracer tendency compute diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index c8c0d84b5348..12a103efd98c 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -153,6 +153,7 @@ class Tendencies { const VertCoord *VCoord, ///< [in] Vertical coordinate VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStep, ///< [in] Time step Config *Options, ///< [in] Configuration options CustomTendencyType InCustomThicknessTend, CustomTendencyType InCustomVelocityTend); @@ -162,6 +163,7 @@ class Tendencies { const VertCoord *VCoord, ///< [in] Vertical coordinate VertAdv *VAdv, ///< [in] Vertical advection int NTracersIn, ///< [in] Number of tracers + TimeInterval TimeStep, ///< [in] Time step Config *Options ///< [in] Configuration options ); @@ -173,6 +175,7 @@ class Tendencies { const VertCoord *VCoord; ///< Pointer to vertical coordinate VertAdv *VAdv; ///< Pointer to vertical advection I4 NTracers; ///< Number of tracers + TimeInterval TimeStep; ///< Time step // Pointer to default tendencies static Tendencies *DefaultTendencies; diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 63177b8721b3..f511fef80bf1 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -163,8 +163,11 @@ int testTendencies() { const auto VAdv = VertAdv::getDefault(); VCoord->NVertLayers = 12; // test creation of another tendencies + + TimeInterval ZeroTimeStep; // Zero-length time step placeholder Config *Options = Config::getOmegaConfig(); - Tendencies::create("TestTendencies", Mesh, VCoord, VAdv, 3, Options); + Tendencies::create("TestTendencies", Mesh, VCoord, VAdv, 3, ZeroTimeStep, + Options); // test retrievel of another tendencies if (Tendencies::get("TestTendencies")) { diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index a0c71c91ee1b..c6a168a1b1af 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -225,8 +225,8 @@ int initTimeStepperTest(const std::string &mesh) { // Creating non-default tendencies with custom velocity tendencies auto *TestTendencies = Tendencies::create( - "TestTendencies", DefMesh, DefVertCoord, DefVAdv, NTracers, &Options, - Tendencies::CustomTendencyType{}, DecayVelocityTendency{}); + "TestTendencies", DefMesh, DefVertCoord, DefVAdv, NTracers, ZeroTimeStep, + &Options, Tendencies::CustomTendencyType{}, DecayVelocityTendency{}); if (!TestTendencies) { Err++; LOG_ERROR("TimeStepperTest: error creating test tendencies"); From 594604bd1899167e5d0e85e3c825537942db79a3 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Wed, 11 Feb 2026 17:16:15 -0800 Subject: [PATCH 34/42] Make TimeStep nonoptional for computeTracerVAdvTend --- components/omega/src/ocn/VertAdv.cpp | 2 +- components/omega/src/ocn/VertAdv.h | 3 +-- components/omega/test/ocn/VertAdvTest.cpp | 3 ++- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index c59b77505ab6..30ea70cf2c8e 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -565,7 +565,7 @@ void VertAdv::computeTracerVAdvTend( const Array3DReal &TracerTend, //< [inout] tracer tendencies const Array3DReal &Tracers, //< [in] tracer array const Array2DReal &LayerThickness, //< [in] layer thickness - const TimeInterval TimeStep //< [in] (optional) time step + const TimeInterval TimeStep //< [in] time step ) { // Return if vertical advection tracer tendency not enabled diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index cae279ffe1fe..f8daaa26d38f 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -149,8 +149,7 @@ class VertAdv { const Array3DReal &TracerTend, ///< [inout] tracer tendencies const Array3DReal &Tracers, ///< [in] tracer array const Array2DReal &LayerThickness, ///< [in] layer thickness - const TimeInterval TimeStep = - TimeInterval() ///< [in] (optional) time step + const TimeInterval TimeStep ///< [in] time step ); /// Compute tracer fluxes due to vertical advection diff --git a/components/omega/test/ocn/VertAdvTest.cpp b/components/omega/test/ocn/VertAdvTest.cpp index 816ea9183192..b8bcc0e22863 100644 --- a/components/omega/test/ocn/VertAdvTest.cpp +++ b/components/omega/test/ocn/VertAdvTest.cpp @@ -59,6 +59,7 @@ Real tendTest(const int NLayers, VertAdv *VAdv) { Array2DReal LayerThick("LayerThickness", 1, NLayers); Array3DReal TracerArray("TracerArray", 1, 1, NLayers); Array3DReal Tend("TracerArray", 1, 1, NLayers); + TimeInterval ZeroTimeStep; // Set uniform velocity throughout domain deepCopy(VAdv->TotalVerticalVelocity, 1._Real); @@ -77,7 +78,7 @@ Real tendTest(const int NLayers, VertAdv *VAdv) { }); // Compute tracer tendencies - VAdv->computeTracerVAdvTend(Tend, TracerArray, LayerThick); + VAdv->computeTracerVAdvTend(Tend, TracerArray, LayerThick, ZeroTimeStep); HostArray3DReal TendH = createHostMirrorCopy(Tend); HostArray1DReal XLayerH = createHostMirrorCopy(XLayer); From 913d9e38a2dcf2fabf501fd5a11b18400ad3b1a5 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 12 Feb 2026 08:47:13 -0800 Subject: [PATCH 35/42] Add VertAdv devGuide and userGuide --- components/omega/doc/devGuide/VertAdv.md | 103 ++++++++++++++++++++++ components/omega/doc/index.md | 2 + components/omega/doc/userGuide/VertAdv.md | 47 ++++++++++ 3 files changed, 152 insertions(+) create mode 100644 components/omega/doc/devGuide/VertAdv.md create mode 100644 components/omega/doc/userGuide/VertAdv.md diff --git a/components/omega/doc/devGuide/VertAdv.md b/components/omega/doc/devGuide/VertAdv.md new file mode 100644 index 000000000000..822c40ee13b6 --- /dev/null +++ b/components/omega/doc/devGuide/VertAdv.md @@ -0,0 +1,103 @@ +(omega-dev-vert-adv)= + +## Vertical Advection + +The `VertAdv` class contains methods and arrays for calculating vertical +velocities and the tendencies of thickness, horizontal velocity, and tracers due +to vertical advection. + +### Initialization + +The default `VertAdv` instance is created by calling `VertAdv::init()` method. +The default `HorzMesh`, `VertCoord`, and `Tracers` objects must be initialized +first. The default instance is retrieved by: +``` +VertAdv *DefVertAdv = VertAdv::getDefault(); +``` + +Additional instances can be created with the `create` method, which call the +constructor and places the new instance in a map identified with the label +`Name`: +``` +VertAdv *VertAdv::create(const std::string &Name, ///< [in] name for new VertAdv + const HorzMesh *Mesh, ///< [in] associated HorzMesh + const VertCoord *VCoord, ///< [in] associated VertCoord + Config *Options ///< [in] configuration options +) +``` +This instance is retrieved with the get method: +``` +VertAdv *NewVertAdv = VertAdv::get("Name"); +``` +The constructor reads configuration options, stores some info from the +`HorzMesh`, `VertCoord`, and `Tracers` objects, allocates needed arrays, +and defines `Field` metadata. + +### Variables + +The following arrays are stored as members of the `VertAdv` class. + +| Variable Name | Type | Dimensions | +| ------------- | ---- | ---------- | +| VerticalVelocity | Real | NCellsSize, NVertLayersP1 | +| TotalVerticalVelocity | Real | NCellsSize, NVertLayersP1 | +| VertFlux | Real | NTracers, NCellsSize, NVertLayersP1 | +| LowOrderVertFlux | Real | NTracers, NCellsSize, NVertLayersP1 | + +The `VerticalVelocity` represents the raw vertical pseudovelocity through the +top interfaces of cell layers, computed from the divergence of the horizontal +velocity. The `TotalVerticalVelocity` is the transport velocity and includes +corrections applied to the `VerticalVelocity`. The `VertFlux` and +`LowOrderVertFlux` store tracer fluxes at layer interfaces. The specific +numerical algorithms to compute these fluxes are chosen by the user through +configuration options. + +### Methods + +The `VerticalVelocity` and `TotalVerticalVelocity` arrays are computed by +calling the `computeVerticalVelocity` method: +``` +VertAdv::computeVerticalVelocity(NormalVelocity, FluxLayerThickEdge); +``` +This method takes as input the `NormalVelocity` field from the `OceanState` +object, and the `FluxLayerThickEdge` field from the `AuxiliaryState`. At +present, `TotalVerticalVelocity` is equivalent to `VerticalVelocity`; +additional corrections will be added in subsequent updates. The +`TotalVerticalVelocity` array is used by each of the tendency methods, +therefore `computeVerticalVelocity` must be called before any tendency +computations. + +There are three methods for computing vertical advection tendencies of +thickness, horizontal velocity and tracers that are called from the time +stepper: `computeThicknessVAdvTend`, `computeVelocityVAdvTend`, and +`computeTracerVAdvTend`. Each method updates the corresponding tendency array +in place (inout). + +Only the thickness tendency array is passed to the `computeThicknessVAdvTend` +method: +``` +VertAdv::computeThicknessVAdvTend(ThickTend); +``` + +The `computeVelocityVAdvTend` method requires the `NormalVelocity` and +`FluxLayerThickEdge` fields, in addition to the velocity tendency array: +``` +VertAdv::computeVelocityVAdvTend(VelTend, NormalVelocity, FluxLayerThickEdge); +``` + +The tracer vertical advection tendency depends on the configured settings. +The `computeTracerVAdvTend` method takes as arguments the full 3D array of +tracers, a thickness array (selected based on the configured advection +algorithm) and a `TimeInterval` representing the time step, along with the +tracer tendency array: +``` +VertAdv::computeTracerVAdvTend(TracerTend, TracerArray, Thickness, TimeStep); +``` +For the standard advection algorithm, `LayerThickness` from the `OceanState` is +passed as the thickness argument; for the FCT algorithm, the `ProvThickness` +from the `AuxiliaryState` is used instead. + +This method first calls `computeVerticalFluxes` to compute the tracer fluxes at +layer interfaces, using the order of accuracy specified in the configuration +file. It then dispatches to the selected advection algorithm (standard or FCT) +to compute the tracer tendencies. diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 121d65a0cbfd..2cbd0facd2e4 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -51,6 +51,7 @@ userGuide/TridiagonalSolvers userGuide/VertCoord userGuide/Timing userGuide/VerticalMixingCoeff +userGuide/VertAdv ``` ```{toctree} @@ -94,6 +95,7 @@ devGuide/TridiagonalSolvers devGuide/VertCoord devGuide/Timing devGuide/VerticalMixingCoeff +devGuide/VertAdv ``` ```{toctree} diff --git a/components/omega/doc/userGuide/VertAdv.md b/components/omega/doc/userGuide/VertAdv.md new file mode 100644 index 000000000000..38511a201956 --- /dev/null +++ b/components/omega/doc/userGuide/VertAdv.md @@ -0,0 +1,47 @@ +(omega-user-vert-adv)= + +## Vertical Advection + +### Overview + +In Omega, vertical advection represents the transport of mass, momentum, and +tracer quantities along the vertical coordinate due to vertical motion within a +water column. In the context of Omega's +[governing equations](omega-design-governing-eqns-omega1), vertical advection +contributes to the tendencies of pseudo-thickness, horizontal velocity, and +tracer fields. Vertical motion is inferred from the divergence of horizontal +flow through the continuity equation and is necessary for accurately +representing vertical transport in three-dimensional ocean simulations. + +The `VertAdv` class contains variables and functions for computing the vertical +velocity and the tendencies of `LayerThickness`, `NormalVelocity`, and +`Tracers`. The algorithms used are determined by options specified in the +configuration file. + +### Configuration Options + +Within the `Tendencies` group of the configuration file, flags are used to +enable or disable the contribution of vertical advection to the tendencies of +thickness, velocity, and tracers: +| Configuration name | Options | +| ------------------ | ------- | +| ThicknessVertAdvTendencyEnable | true / false | +| VelocityVertAdvTendencyEnable | true / false | +| TracerVertAdvTendencyEnable | true / false | + +Within the `Advection` group, options are provided to configure the algorithms +used to compute the tracer tendencies: +| Configuration name | Description | Options | +| ------------------ | ----------- | ------- | +| VerticalTracerFluxOrder | order of accuracy used for the tracer flux calculation | 2, 3, or 4 | +| VerticalTracerFluxLimiterEnable | selects the tracer tendency algorithm | true / false | +| Coef3rdOrder | coefficient used for blending 3rd- and 4th-order fluxes when `VerticalTracerFluxOrder == 3` | 0.0 - 1.0 | + +For `VerticalTracerFluxLimiterEnable`, `true` enables flux-corrected transport, +and `false` selects the standard algorithm. For `Coef3rdOrder`, `1` corresponds +to purely 3rd-order, `0` purely 4th-order. + +The `VertAdv` class provides four fields that can be written to output by +adding them to the contents of an output stream in the configuration file: +`VerticalVelocity`, `TotalVerticalVelocity`, `VertFlux`, and `LowOrderVertFlux`. +These fields collectively make up the `VertAdv` output group. From bf7f3f7ed3ff24450fcf451dde33a2ed9c2e68a8 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Fri, 13 Feb 2026 18:47:30 -0600 Subject: [PATCH 36/42] Add early return to VertAdv methods when NVertLayers == 1 --- components/omega/src/ocn/VertAdv.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 30ea70cf2c8e..21e1c94456bc 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -361,6 +361,10 @@ void VertAdv::computeVerticalVelocity( const Array2DReal &FluxLayerThickEdge //< [in] layer thickness at edges ) { + // Return if mesh only has a single vertical layer + if (NVertLayers == 1) + return; + OMEGA_SCOPE(LocVertVel, VerticalVelocity); OMEGA_SCOPE(LocNVertLayers, NVertLayers); OMEGA_SCOPE(LocAreaCell, Mesh->AreaCell); @@ -449,6 +453,10 @@ void VertAdv::computeThicknessVAdvTend( if (!ThickVertAdvEnabled) return; + // Return if mesh only has a single vertical layer + if (NVertLayers == 1) + return; + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(LocTotVertVelocity, TotalVerticalVelocity); @@ -489,6 +497,10 @@ void VertAdv::computeVelocityVAdvTend( if (!VelVertAdvEnabled) return; + // Return if mesh only has a single vertical layer + if (NVertLayers == 1) + return; + OMEGA_SCOPE(LocCOnE, Mesh->CellsOnEdge); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); @@ -572,6 +584,10 @@ void VertAdv::computeTracerVAdvTend( if (!TracerVertAdvEnabled) return; + // Return if mesh only has a single vertical layer + if (NVertLayers == 1) + return; + // Compute tracer fluxes at the interfaces computeVerticalFluxes(Tracers, LayerThickness); From e813b29ce9f3c3042f9675db9d6db793970bce6d Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 19 Feb 2026 17:20:58 -0600 Subject: [PATCH 37/42] Fix typos in comments and docs --- components/omega/doc/devGuide/VertAdv.md | 2 +- components/omega/src/ocn/AuxiliaryState.h | 4 ++-- components/omega/src/ocn/VertAdv.cpp | 2 +- components/omega/src/ocn/VertAdv.h | 2 +- components/omega/test/ocn/VertAdvTest.cpp | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/components/omega/doc/devGuide/VertAdv.md b/components/omega/doc/devGuide/VertAdv.md index 822c40ee13b6..672a25ef550e 100644 --- a/components/omega/doc/devGuide/VertAdv.md +++ b/components/omega/doc/devGuide/VertAdv.md @@ -15,7 +15,7 @@ first. The default instance is retrieved by: VertAdv *DefVertAdv = VertAdv::getDefault(); ``` -Additional instances can be created with the `create` method, which call the +Additional instances can be created with the `create` method, which calls the constructor and places the new instance in a map identified with the label `Name`: ``` diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 60a414d39d35..f3c14be20ac3 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -53,7 +53,7 @@ class AuxiliaryState { // Create a non-default auxiliary state static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, const VertCoord *VCoord, - VertAdv *Vadv, int NTracers, + VertAdv *VAdv, int NTracers, TimeInterval TimeStep); /// Get the default auxiliary state @@ -87,7 +87,7 @@ class AuxiliaryState { private: AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo, - const VertCoord *VCoord, VertAdv *Vadv, int NTracers, + const VertCoord *VCoord, VertAdv *VAdv, int NTracers, TimeInterval TimeStep); AuxiliaryState(const AuxiliaryState &) = delete; diff --git a/components/omega/src/ocn/VertAdv.cpp b/components/omega/src/ocn/VertAdv.cpp index 21e1c94456bc..343136b10d34 100644 --- a/components/omega/src/ocn/VertAdv.cpp +++ b/components/omega/src/ocn/VertAdv.cpp @@ -572,7 +572,7 @@ void VertAdv::computeVelocityVAdvTend( //------------------------------------------------------------------------------ // Compute tracer tendency due to vertical advection, TimeStep is only needed -// as an arugement for flux-corrected transport +// as an argument for flux-corrected transport void VertAdv::computeTracerVAdvTend( const Array3DReal &TracerTend, //< [inout] tracer tendencies const Array3DReal &Tracers, //< [in] tracer array diff --git a/components/omega/src/ocn/VertAdv.h b/components/omega/src/ocn/VertAdv.h index f8daaa26d38f..6d74b91835ce 100644 --- a/components/omega/src/ocn/VertAdv.h +++ b/components/omega/src/ocn/VertAdv.h @@ -115,7 +115,7 @@ class VertAdv { /// Retrieve the default VertAdv static VertAdv *getDefault(); - /// Retreive a VertAdv by name + /// Retrieve a VertAdv by name static VertAdv *get(std::string name); /// Read and set config options diff --git a/components/omega/test/ocn/VertAdvTest.cpp b/components/omega/test/ocn/VertAdvTest.cpp index b8bcc0e22863..e60a6586aed0 100644 --- a/components/omega/test/ocn/VertAdvTest.cpp +++ b/components/omega/test/ocn/VertAdvTest.cpp @@ -203,7 +203,7 @@ int main(int argc, char *argv[]) { OMEGA_SCOPE(LocEOnC, DefDecomp->EdgesOnCell); OMEGA_SCOPE(LocESOnC, DefMesh->EdgeSignOnCell); - // Set unifrom values for layer thickness and velocity on edges of column + // Set uniform values for layer thickness and velocity on edges of column parallelFor( {NVertLayers}, KOKKOS_LAMBDA(int K) { for (int J = 0; J < NEdgesOnCell0; ++J) { From beeec710d302a66a8689c9886b625ea2deb2054c Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 19 Feb 2026 18:40:25 -0600 Subject: [PATCH 38/42] Update HorzMesh user doc --- components/omega/doc/userGuide/HorzMesh.md | 1 - 1 file changed, 1 deletion(-) diff --git a/components/omega/doc/userGuide/HorzMesh.md b/components/omega/doc/userGuide/HorzMesh.md index 7d50db5c529e..2160b51c82ae 100644 --- a/components/omega/doc/userGuide/HorzMesh.md +++ b/components/omega/doc/userGuide/HorzMesh.md @@ -23,7 +23,6 @@ This includes the following variables: | XCell, YCell, ZCell | Cartesian coordinates of cell centers | m | | XEdge, YEdge, ZEdge | Cartesian coordinates of edge centers | m | | XVertex, YVertex, ZVertex | Cartesian coordinates of vertices | m | -| BottomDepth | Depth of the bottom of the ocean at cell centers | m | | FCell, FEdge, FVertex | Coriolis parameter at cell centers/edges/vertices | radians/s | | LonCell, LatCell | Longitude/latitude coordinates of cell centers | radians | | LonEdge, LatEdge | Longitude/latitude coordinates of edge centers | radians | From 44bd58a3c54f7a6739ea955903df92ed80d802a5 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 19 Feb 2026 18:41:13 -0600 Subject: [PATCH 39/42] Add sanity check for bottom depth --- components/omega/test/ocn/VertCoordTest.cpp | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/components/omega/test/ocn/VertCoordTest.cpp b/components/omega/test/ocn/VertCoordTest.cpp index 6c4648e7ac36..8fe474f41018 100644 --- a/components/omega/test/ocn/VertCoordTest.cpp +++ b/components/omega/test/ocn/VertCoordTest.cpp @@ -104,6 +104,26 @@ int main(int argc, char *argv[]) { I4 VertexDegree = DefMesh->VertexDegree; I4 NVertLayers = DefVertCoord->NVertLayers; + // Rest bottom depth successful read + R8 MaxBathy = -1e10; + R8 MinBathy = 1e10; + + for (int ICell = 0; ICell < NCellsOwned; ++ICell) { + if (DefVertCoord->BottomDepthH(ICell) < MinBathy) { + MinBathy = DefVertCoord->BottomDepthH(ICell); + } + if (DefVertCoord->BottomDepthH(ICell) > MaxBathy) { + MaxBathy = DefVertCoord->BottomDepthH(ICell); + } + } + + if ((MinBathy > 0) && (MaxBathy < 11000.)) { + LOG_INFO("VertCoordTest: Bathy min/max test PASS"); + } else { + ErrAll += Error( + ErrorCode::Fail, "VertCoordTest: Bathy min/max test FAIL"); + } + // Tests for computePressure Array2DReal LayerThickness("LayerThickness", NCellsSize, NVertLayers); From 6a34375b91f9aef907a97bad5b112ba82c16a96c Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Thu, 19 Feb 2026 18:50:16 -0600 Subject: [PATCH 40/42] Fix linting fail --- components/omega/test/ocn/VertCoordTest.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/omega/test/ocn/VertCoordTest.cpp b/components/omega/test/ocn/VertCoordTest.cpp index 8fe474f41018..a1289cb5fef1 100644 --- a/components/omega/test/ocn/VertCoordTest.cpp +++ b/components/omega/test/ocn/VertCoordTest.cpp @@ -120,8 +120,8 @@ int main(int argc, char *argv[]) { if ((MinBathy > 0) && (MaxBathy < 11000.)) { LOG_INFO("VertCoordTest: Bathy min/max test PASS"); } else { - ErrAll += Error( - ErrorCode::Fail, "VertCoordTest: Bathy min/max test FAIL"); + ErrAll += + Error(ErrorCode::Fail, "VertCoordTest: Bathy min/max test FAIL"); } // Tests for computePressure From bda6af4e956ee8798d5730ce5435d632a91848c4 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Fri, 20 Feb 2026 09:30:52 -0600 Subject: [PATCH 41/42] Rename TracerTend to LocTracerTend in function call --- components/omega/src/ocn/Tendencies.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 0ff4242484b4..6f41c50a3126 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -612,7 +612,7 @@ void Tendencies::computeTracerTendenciesOnly( } else if (VAdv->VertAdvChoice == VertAdvOption::FCT) { ThicknessForVAdv = AuxState->LayerThicknessAux.ProvThickness; } - VAdv->computeTracerVAdvTend(TracerTend, TracerArray, ThicknessForVAdv, + VAdv->computeTracerVAdvTend(LocTracerTend, TracerArray, ThicknessForVAdv, TimeStep); Pacer::stop("Tend:computeTracerVAdvTend", 2); From 10eec4853c5d5cff072890e5aa88c2af006dedd2 Mon Sep 17 00:00:00 2001 From: Brian O'Neill Date: Tue, 24 Feb 2026 08:21:43 -0600 Subject: [PATCH 42/42] Remove commented-out line --- components/omega/test/ocn/TendenciesTest.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index f511fef80bf1..3a5418830e43 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -128,8 +128,6 @@ int initTendenciesTest(const std::string &mesh) { Tracers::init(); VertAdv::init(); - // VertCoord::init2(); - int StateErr = OceanState::init(); if (StateErr != 0) { Err++;