From dc5a8da02e16879704a4c43005061427508ce16f Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 26 Jan 2026 10:23:14 -0800 Subject: [PATCH 01/15] starts adding surface tracer restoring --- cime | 2 +- components/omega/src/ocn/AuxiliaryState.cpp | 45 +++++++++- components/omega/src/ocn/AuxiliaryState.h | 2 + components/omega/src/ocn/Tendencies.cpp | 87 ++++++++++++++++++- components/omega/src/ocn/Tendencies.h | 1 + components/omega/src/ocn/TendencyTerms.cpp | 4 + components/omega/src/ocn/TendencyTerms.h | 26 ++++++ .../auxiliaryVars/SurfTracerRestAuxVars.cpp | 60 +++++++++++++ .../ocn/auxiliaryVars/SurfTracerRestAuxVars.h | 55 ++++++++++++ externals/ekat | 2 +- externals/scorpio | 2 +- 11 files changed, 278 insertions(+), 8 deletions(-) create mode 100644 components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp create mode 100644 components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h diff --git a/cime b/cime index bd710a378bbb..c64260ed94af 160000 --- a/cime +++ b/cime @@ -1 +1 @@ -Subproject commit bd710a378bbb5ad9efef5d958b63dab221e7b1fe +Subproject commit c64260ed94aff167a59635ed4a9ae5af16824e91 diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 4de2a47884f0..247ccf7d17ce 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -4,6 +4,7 @@ #include "Logging.h" #include "Pacer.h" #include "TimeStepper.h" +#include "Tendencies.h" namespace OMEGA { @@ -28,7 +29,8 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux(stripDefault(Name), Mesh, VCoord), VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), WindForcingAux(stripDefault(Name), Mesh), - TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) { + TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) { //, + // SurfTracerRestAux(stripDefault(Name), Mesh, VCoord, NTracers) { GroupName = "AuxiliaryState"; if (Name != "Default") { @@ -44,6 +46,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VelocityDel2Aux.registerFields(GroupName, AuxMeshName); WindForcingAux.registerFields(GroupName, AuxMeshName); TracerAux.registerFields(GroupName, AuxMeshName); + // SurfTracerRestAux.registerFields(GroupName, AuxMeshName); } // Destructor. Unregisters the fields with IOStreams and destroys this auxiliary @@ -55,6 +58,7 @@ AuxiliaryState::~AuxiliaryState() { VelocityDel2Aux.unregisterFields(); WindForcingAux.unregisterFields(); TracerAux.unregisterFields(); + // SurfTracerRestAux.unregisterFields(); FieldGroup::destroy(GroupName); } @@ -220,13 +224,17 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, void AuxiliaryState::computeAll(const OceanState *State, const Array3DReal &TracerArray, int ThickTimeLevel, int VelTimeLevel) const { - Array2DReal LayerThickCell = State->getLayerThickness(ThickTimeLevel); - Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); - + Array2DReal LayerThickCell; + Array2DReal NormalVelEdge; + State->getLayerThickness(LayerThickCell, ThickTimeLevel); + State->getNormalVelocity(NormalVelEdge, VelTimeLevel); + const int NTracers = TracerArray.extent_int(0); OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); OMEGA_SCOPE(LocTracerAux, TracerAux); + // OMEGA_SCOPE(LocSurfTracerRestAux, SurfTracerRestAux); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); @@ -274,6 +282,23 @@ void AuxiliaryState::computeAll(const OceanState *State, }); Pacer::stop("AuxState:cellAuxState4", 2); + // Pacer::start("AuxState:cellAuxState5", 2); + // parallelForOuter( + // "cellAuxState5", {NTracers, Mesh->NCellsAll}, + // KOKKOS_LAMBDA(int LTracer, 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) { + // LocSurfTracerRestAux.computeVarsOnCells( + // LTracer, ICell, KChunk, TracerArray, + // TracersMonthlySurfClimoArray); + // }); + // }); + // Pacer::stop("AuxState:cellAuxState5", 2); + Pacer::stop("AuxState:computeAll", 1); } @@ -393,6 +418,18 @@ void AuxiliaryState::readConfigOptions(Config *OmegaConfig) { } else { ABORT_ERROR("AuxiliaryState: Unknown InterpType requested"); } + + // Config SurfRestConfig("SurfaceRestoring"); + // Err += OmegaConfig->get(SurfRestConfig); + + // Err += SurfRestConfig.get("MaxDiff", SurfTracerRestAux.MaxDiff); + // CHECK_ERROR_ABORT( + // Err, "AuxiliaryState: MaxDiff not found in SurfaceRestoringConfig"); + // Err += SurfRestConfig.get("PistonVelocity", + // Tendencies::getDefault()->SurfaceTracerRestoring.PistonVelocity); + // CHECK_ERROR_ABORT( + // Err, "AuxiliaryState: PistonVelocity not found in + // SurfaceRestoringConfig"); } //------------------------------------------------------------------------------ diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index f3c14be20ac3..94ded5b70eda 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -16,6 +16,7 @@ #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" #include "auxiliaryVars/WindForcingAuxVars.h" +// #include "auxiliaryVars/SurfTracerRestAuxVars.h" #include #include @@ -42,6 +43,7 @@ class AuxiliaryState { VorticityAuxVars VorticityAux; VelocityDel2AuxVars VelocityDel2Aux; WindForcingAuxVars WindForcingAux; + // SurfTracerRestAuxVars SurfTracerRestAux; ~AuxiliaryState(); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 3c3aacb8dde5..fadf3162040b 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -248,9 +248,70 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config CHECK_ERROR_ABORT(Err, "Tendencies: EddyDiff4 not found in TendConfig"); } + Err += TendConfig.get("PressureGradTendencyEnable", this->PGrad->Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("SurfaceTracerRestoringEnable", + this->SurfaceTracerRestoring.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: SurfaceTracerRestoringEnable not found in TendConfig"); + if (this->SurfaceTracerRestoring.Enabled) { + Config SurfRestConfig("SurfaceRestoring"); + Err += OmegaConfig->get(SurfRestConfig); + Err += SurfRestConfig.get("PistonVelocity", + this->SurfaceTracerRestoring.PistonVelocity); + CHECK_ERROR_ABORT( + Err, + "Tendencies: PistonVelocity not found in SurfaceRestoringConfig"); + + std::vector TracersToRestore; + SurfRestConfig.get("TracersToRestore", TracersToRestore); + + // Enable restoring for specified individual tracers + I4 NumInvalidTracers = 0; + std::vector TracerIdsToRestoreVec; + for (const auto &TracerName : TracersToRestore) { + I4 TracerIndex = -1; + Tracers::getIndex(TracerIndex, TracerName); + if (TracerIndex == -1) { + LOG_ERROR("Tendencies: Tracer {} in TracersToRestore is not " + "defined", + TracerName); + NumInvalidTracers++; + } else { + bool IsDuplicate = false; + for (const auto ExistingTracerIndex : TracerIdsToRestoreVec) { + if (ExistingTracerIndex == TracerIndex) { + IsDuplicate = true; + break; + } + } + if (!IsDuplicate) { + TracerIdsToRestoreVec.push_back(TracerIndex); + } + } + } + if (NumInvalidTracers > 0) { + ABORT_ERROR("Tendencies: {} invalid tracer(s) in TracersToRestore " + "configuration", + NumInvalidTracers); + } + + this->SurfaceTracerRestoring.NTracersToRestore = + static_cast(TracerIdsToRestoreVec.size()); + this->SurfaceTracerRestoring.TracerIdsToRestore = Array1DI4( + "TracerIdsToRestore", this->SurfaceTracerRestoring.NTracersToRestore); + HostArray1DI4 TracerIdsToRestoreH( + "TracerIdsToRestoreH", + this->SurfaceTracerRestoring.NTracersToRestore); + for (I4 R = 0; R < this->SurfaceTracerRestoring.NTracersToRestore; ++R) { + TracerIdsToRestoreH(R) = TracerIdsToRestoreVec[R]; + } + deepCopy(this->SurfaceTracerRestoring.TracerIdsToRestore, + TracerIdsToRestoreH); + } } //------------------------------------------------------------------------------ @@ -325,6 +386,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VelocityHyperDiff(Mesh, VCoord), WindForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), + TracerDiffusion(Mesh, VCoord), SurfaceTracerRestoring(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad) { @@ -653,6 +715,7 @@ void Tendencies::computeTracerTendenciesOnly( OMEGA_SCOPE(LocTracerHorzAdv, TracerHorzAdv); OMEGA_SCOPE(LocTracerDiffusion, TracerDiffusion); OMEGA_SCOPE(LocTracerHyperDiff, TracerHyperDiff); + OMEGA_SCOPE(LocSurfaceTracerRestoring, SurfaceTracerRestoring); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); @@ -746,8 +809,8 @@ void Tendencies::computeTracerTendenciesOnly( Pacer::stop("Tend:tracerHyperDiff", 2); } + // compute tracer tendency from vertical advection Pacer::start("Tend:computeTracerVAdvTend", 2); - // compute tracer tendencies from vertical advection Array2DReal ThicknessForVAdv; if (VAdv->VertAdvChoice == VertAdvOption::Standard) { ThicknessForVAdv = State->getLayerThickness(ThickTimeLevel); @@ -758,6 +821,28 @@ void Tendencies::computeTracerTendenciesOnly( TimeStep); Pacer::stop("Tend:computeTracerVAdvTend", 2); + // compute tracer surface restoring + const Array2DReal &TracersMonthlySurfClimo = + AuxState->SurfTracerRestAux.TracersMonthlySurfClimoCell; + const I4 NTracersToRestore = LocSurfaceTracerRestoring.NTracersToRestore; + const auto &TracerIdsToRestore = + LocSurfaceTracerRestoring.TracerIdsToRestore; + if (LocSurfaceTracerRestoring.Enabled && NTracersToRestore > 0) { + Pacer::start("Tend:surfaceTracerRestoring", 2); + parallelFor( + {NTracersToRestore, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int R, int ICell) { + const int KMin = MinLayerCell(ICell); + const int L = TracerIdsToRestore(R); + LocSurfaceTracerRestoring(LocTracerTend, L, ICell, KMin, + TracersMonthlySurfClimo, TracerArray); + }); + Pacer::stop("Tend:surfaceTracerRestoring", 2); + } else if (LocSurfaceTracerRestoring.Enabled && NTracersToRestore == 0) { + ABORT_ERROR("Tendencies: SurfaceTracerRestoring is enabled but " + "TracerIdsToRestore is empty"); + } + 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 d5687f9274a6..96651eeacf7a 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -71,6 +71,7 @@ class Tendencies { TracerHorzAdvOnCell TracerHorzAdv; TracerDiffOnCell TracerDiffusion; TracerHyperDiffOnCell TracerHyperDiff; + SurfaceTracerRestoringOnCell SurfaceTracerRestoring; std::string Name; diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 74dbcd033649..35efb216a34d 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -111,6 +111,9 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, MinLayerEdgeBot(VCoord->MinLayerEdgeBot), MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} +SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell( + const HorzMesh *Mesh) {} + void TracerHorzAdvOnCell::init() { const HorzMesh *Mesh = this->HorzontalMesh; const auto MaxEdges2 = Mesh->MaxEdges2; @@ -133,6 +136,7 @@ void TracerHorzAdvOnCell::init() { {NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { masksAndCoefficients(IEdge); }); Kokkos::fence(); } + } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 55a108e8d3f5..8e92b842646e 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -570,5 +570,31 @@ class TracerHyperDiffOnCell { Array1DI4 MaxLayerEdgeTop; }; +/// Surface tracer restoring term +class SurfaceTracerRestoringOnCell { + public: + bool Enabled; + Real PistonVelocity = 1.585e-5; + + /// constructor declaration + SurfaceTracerRestoringOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + /// The functor takes the cell index and the array for the tracer surface + /// restoring values, outputs tendency array + KOKKOS_FUNCTION void + operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, + const Array2DReal &SurfTracerRestValuesCell) const { + if (KChunk == 0) { + const I4 K = MinLayerCell(ICell); + + Tend(L, ICell, K) += + PistonVelocity * SurfTracerRestValuesCell(L, ICell); + } + } + + private: + Array1DI4 MinLayerCell; +}; + } // namespace OMEGA #endif diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp new file mode 100644 index 000000000000..852efd5da553 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp @@ -0,0 +1,60 @@ +#include "SurfTracerRestAuxVars.h" +#include "DataTypes.h" +#include "Field.h" + +#include + +namespace OMEGA { + +SurfTracerRestAuxVars::SurfTracerRestAuxVars(const std::string &AuxStateSuffix, + const HorzMesh *Mesh, + const VertCoord *VCoord, + const I4 NTracers) + : SurfTracerRestValuesCell("SurfTracerRestValuesCell" + AuxStateSuffix, + NTracers, Mesh->NCellsSize), + MinLayerCell(VCoord->MinLayerCell) {} + +void SurfTracerRestAuxVars::registerFields( + const std::string &AuxGroupName, // name of Auxiliary field group + const std::string &MeshName // name of horizontal mesh +) const { + + // Create fields + const Real FillValue = -9.99e30; + int NDims = 2; + std::vector DimNames(NDims); + std::string DimSuffix; + if (MeshName == "Default") { + DimSuffix = ""; + } else { + DimSuffix = MeshName; + } + + // Zonal wind stress + DimNames[0] = "NTracers"; + DimNames[1] = "NCells" + DimSuffix; + auto SurfTracerRestValuesCellField = + Field::create(SurfTracerRestValuesCell.label(), // field name + "tracer surface restoring value", // long name/describe + "tracer units", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValue, // scalar for undefined entries + NDims, // number of dimensions + DimNames // dim names + ); + + // Add fields to FieldGroup + FieldGroup::addFieldToGroup(SurfTracerRestValuesCell.label(), AuxGroupName); + + // Attach data + SurfTracerRestValuesCellField->attachData( + SurfTracerRestValuesCell); +} + +void SurfTracerRestAuxVars::unregisterFields() const { + Field::destroy(SurfTracerRestValuesCell.label()); +} + +} // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h new file mode 100644 index 000000000000..a43472dd11f4 --- /dev/null +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h @@ -0,0 +1,55 @@ +#ifndef OMEGA_AUX_SURF_REST_H +#define OMEGA_AUX_SURF_REST_H + +#include "DataTypes.h" +#include "HorzMesh.h" +#include "HorzOperators.h" +#include "OmegaKokkos.h" +#include "VertCoord.h" + +#include + +namespace OMEGA { + +class SurfTracerRestAuxVars { + public: + Array2DReal SurfTracerRestValuesCell; + Real MaxDiff = 100.0; // maximum allowed difference for restoring + /// Need to add under sea ice restoring option when that is available + + SurfTracerRestAuxVars(const std::string &AuxStateSuffix, + const HorzMesh *Mesh, const VertCoord *VCoord, + const I4 NTracers); + + KOKKOS_FUNCTION void + computeVarsOnCells(int L, int ICell, int KChunk, + const Array3DReal &TracerCell, + const Array2DReal &TracersMonthlySurfClimoCell) const { + + if (KChunk == 0) { + const I4 K = MinLayerCell(ICell); + + Real DeltaTracer; + if ((TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K)) > + MaxDiff) { + SurfTracerRestValuesCell(L, ICell) = MaxDiff; + } else if ((TracersMonthlySurfClimoCell(L, ICell) - + TracerCell(L, ICell, K)) < -MaxDiff) { + SurfTracerRestValuesCell(L, ICell) = -MaxDiff; + } else { + SurfTracerRestValuesCell(L, ICell) = + TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K); + } + } + } + + void registerFields(const std::string &AuxGroupName, + const std::string &MeshName) const; + void unregisterFields() const; + + private: + Array1DI4 MinLayerCell; +}; + +} // namespace OMEGA +#endif diff --git a/externals/ekat b/externals/ekat index d13cf66f3543..7156c944701f 160000 --- a/externals/ekat +++ b/externals/ekat @@ -1 +1 @@ -Subproject commit d13cf66f3543df10445cd1fd35edac3165cc2833 +Subproject commit 7156c944701fd19fa5e27ffcb86819f12851ae67 diff --git a/externals/scorpio b/externals/scorpio index a486b499e811..9df851035a89 160000 --- a/externals/scorpio +++ b/externals/scorpio @@ -1 +1 @@ -Subproject commit a486b499e8113f92c136f020cc0e7a68c664f01a +Subproject commit 9df851035a89ba1ca0581ccdbb88fcfc0ec38992 From 1717bfa1bae76ffbf1d19b7904ac0fad448de487 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Wed, 4 Mar 2026 08:33:23 -0800 Subject: [PATCH 02/15] completes working surface forcing --- components/omega/configs/Default.yml | 4 + components/omega/src/ocn/AuxiliaryState.cpp | 50 ++++------ components/omega/src/ocn/AuxiliaryState.h | 4 +- components/omega/src/ocn/Tendencies.cpp | 12 ++- .../auxiliaryVars/SurfTracerRestAuxVars.cpp | 5 +- .../ocn/auxiliaryVars/SurfTracerRestAuxVars.h | 32 +++---- .../omega/test/ocn/AuxiliaryVarsTest.cpp | 91 +++++++++++++++++++ .../omega/test/ocn/TendencyTermsTest.cpp | 78 ++++++++++++++++ 8 files changed, 218 insertions(+), 58 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index eb393b4d0982..7d166110429a 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -31,6 +31,9 @@ Omega: VerticalTracerFluxOrder: 3 WindStress: InterpType: Isotropic + SurfaceRestoring: + MaxDiff: 100.0 + PistonVelocity: 1.585e-5 VertCoord: MovementWeightType: Uniform PressureGrad: @@ -46,6 +49,7 @@ Omega: ViscDel4: 1.2e11 DivFactor: 1.0 WindForcingTendencyEnable: false + SurfaceTracerRestoringEnable: true BottomDragTendencyEnable: false BottomDragCoeff: 0.0 TracerHorzAdvTendencyEnable: true diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 247ccf7d17ce..6031f7898658 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -29,8 +29,8 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux(stripDefault(Name), Mesh, VCoord), VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), WindForcingAux(stripDefault(Name), Mesh), - TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) { //, - // SurfTracerRestAux(stripDefault(Name), Mesh, VCoord, NTracers) { + TracerAux(stripDefault(Name), Mesh, VCoord, NTracers), + SurfTracerRestAux(stripDefault(Name), Mesh, VCoord, NTracers) { GroupName = "AuxiliaryState"; if (Name != "Default") { @@ -46,7 +46,7 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VelocityDel2Aux.registerFields(GroupName, AuxMeshName); WindForcingAux.registerFields(GroupName, AuxMeshName); TracerAux.registerFields(GroupName, AuxMeshName); - // SurfTracerRestAux.registerFields(GroupName, AuxMeshName); + SurfTracerRestAux.registerFields(GroupName, AuxMeshName); } // Destructor. Unregisters the fields with IOStreams and destroys this auxiliary @@ -58,7 +58,7 @@ AuxiliaryState::~AuxiliaryState() { VelocityDel2Aux.unregisterFields(); WindForcingAux.unregisterFields(); TracerAux.unregisterFields(); - // SurfTracerRestAux.unregisterFields(); + SurfTracerRestAux.unregisterFields(); FieldGroup::destroy(GroupName); } @@ -233,7 +233,7 @@ void AuxiliaryState::computeAll(const OceanState *State, OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); OMEGA_SCOPE(LocTracerAux, TracerAux); - // OMEGA_SCOPE(LocSurfTracerRestAux, SurfTracerRestAux); + OMEGA_SCOPE(LocSurfTracerRestAux, SurfTracerRestAux); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); @@ -282,22 +282,13 @@ void AuxiliaryState::computeAll(const OceanState *State, }); Pacer::stop("AuxState:cellAuxState4", 2); - // Pacer::start("AuxState:cellAuxState5", 2); - // parallelForOuter( - // "cellAuxState5", {NTracers, Mesh->NCellsAll}, - // KOKKOS_LAMBDA(int LTracer, 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) { - // LocSurfTracerRestAux.computeVarsOnCells( - // LTracer, ICell, KChunk, TracerArray, - // TracersMonthlySurfClimoArray); - // }); - // }); - // Pacer::stop("AuxState:cellAuxState5", 2); + Pacer::start("AuxState:cellAuxState5", 2); + parallelFor( + "cellAuxState5", {NTracers, Mesh->NCellsAll}, + KOKKOS_LAMBDA(int LTracer, int ICell) { + LocSurfTracerRestAux.computeVarsOnCells(LTracer, ICell, TracerArray); + }); + Pacer::stop("AuxState:cellAuxState5", 2); Pacer::stop("AuxState:computeAll", 1); } @@ -419,17 +410,12 @@ void AuxiliaryState::readConfigOptions(Config *OmegaConfig) { ABORT_ERROR("AuxiliaryState: Unknown InterpType requested"); } - // Config SurfRestConfig("SurfaceRestoring"); - // Err += OmegaConfig->get(SurfRestConfig); - - // Err += SurfRestConfig.get("MaxDiff", SurfTracerRestAux.MaxDiff); - // CHECK_ERROR_ABORT( - // Err, "AuxiliaryState: MaxDiff not found in SurfaceRestoringConfig"); - // Err += SurfRestConfig.get("PistonVelocity", - // Tendencies::getDefault()->SurfaceTracerRestoring.PistonVelocity); - // CHECK_ERROR_ABORT( - // Err, "AuxiliaryState: PistonVelocity not found in - // SurfaceRestoringConfig"); + Config SurfRestConfig("SurfaceRestoring"); + Err += OmegaConfig->get(SurfRestConfig); + + Err += SurfRestConfig.get("MaxDiff", this->SurfTracerRestAux.MaxDiff); + CHECK_ERROR_ABORT( + Err, "AuxiliaryState: MaxDiff not found in SurfaceRestoringConfig"); } //------------------------------------------------------------------------------ diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h index 94ded5b70eda..1228f06962d4 100644 --- a/components/omega/src/ocn/AuxiliaryState.h +++ b/components/omega/src/ocn/AuxiliaryState.h @@ -12,11 +12,11 @@ #include "VertCoord.h" #include "auxiliaryVars/KineticAuxVars.h" #include "auxiliaryVars/LayerThicknessAuxVars.h" +#include "auxiliaryVars/SurfTracerRestAuxVars.h" #include "auxiliaryVars/TracerAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" #include "auxiliaryVars/WindForcingAuxVars.h" -// #include "auxiliaryVars/SurfTracerRestAuxVars.h" #include #include @@ -43,7 +43,7 @@ class AuxiliaryState { VorticityAuxVars VorticityAux; VelocityDel2AuxVars VelocityDel2Aux; WindForcingAuxVars WindForcingAux; - // SurfTracerRestAuxVars SurfTracerRestAux; + SurfTracerRestAuxVars SurfTracerRestAux; ~AuxiliaryState(); diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index fadf3162040b..b3e6de666024 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -142,12 +142,17 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config ) { Error Err; // error code + // we need access to the top-level Omega config when pulling values + // from other sections such as SurfaceRestoring + Config *OmegaConfig = Config::getOmegaConfig(); + Config TendConfig("Tendencies"); Err += OmegaConfig->get(TendConfig); CHECK_ERROR_ABORT(Err, "Tendencies: Tendencies group not found in Config"); Err += TendConfig.get("ThicknessFluxTendencyEnable", this->ThicknessFluxDiv.Enabled); + CHECK_ERROR_ABORT( Err, "Tendencies: ThicknessFluxTendencyEnable not found in TendConfig"); @@ -252,7 +257,7 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config Err += TendConfig.get("PressureGradTendencyEnable", this->PGrad->Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); - + Err += TendConfig.get("SurfaceTracerRestoringEnable", this->SurfaceTracerRestoring.Enabled); CHECK_ERROR_ABORT( @@ -838,11 +843,8 @@ void Tendencies::computeTracerTendenciesOnly( TracersMonthlySurfClimo, TracerArray); }); Pacer::stop("Tend:surfaceTracerRestoring", 2); - } else if (LocSurfaceTracerRestoring.Enabled && NTracersToRestore == 0) { - ABORT_ERROR("Tendencies: SurfaceTracerRestoring is enabled but " - "TracerIdsToRestore is empty"); } - + Pacer::stop("Tend:computeTracerTendenciesOnly", 1); } // end tracer tendency compute diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp index 852efd5da553..3837c131f9e0 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp @@ -12,6 +12,9 @@ SurfTracerRestAuxVars::SurfTracerRestAuxVars(const std::string &AuxStateSuffix, const I4 NTracers) : SurfTracerRestValuesCell("SurfTracerRestValuesCell" + AuxStateSuffix, NTracers, Mesh->NCellsSize), + TracersMonthlySurfClimoCell("TracersMonthlySurfClimoCell" + + AuxStateSuffix, + NTracers, Mesh->NCellsSize), MinLayerCell(VCoord->MinLayerCell) {} void SurfTracerRestAuxVars::registerFields( @@ -30,7 +33,7 @@ void SurfTracerRestAuxVars::registerFields( DimSuffix = MeshName; } - // Zonal wind stress + // Tracers surface restoring values DimNames[0] = "NTracers"; DimNames[1] = "NCells" + DimSuffix; auto SurfTracerRestValuesCellField = diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h index a43472dd11f4..2411eecaf75e 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h @@ -14,6 +14,7 @@ namespace OMEGA { class SurfTracerRestAuxVars { public: Array2DReal SurfTracerRestValuesCell; + Array2DReal TracersMonthlySurfClimoCell; Real MaxDiff = 100.0; // maximum allowed difference for restoring /// Need to add under sea ice restoring option when that is available @@ -22,24 +23,19 @@ class SurfTracerRestAuxVars { const I4 NTracers); KOKKOS_FUNCTION void - computeVarsOnCells(int L, int ICell, int KChunk, - const Array3DReal &TracerCell, - const Array2DReal &TracersMonthlySurfClimoCell) const { - - if (KChunk == 0) { - const I4 K = MinLayerCell(ICell); - - Real DeltaTracer; - if ((TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K)) > - MaxDiff) { - SurfTracerRestValuesCell(L, ICell) = MaxDiff; - } else if ((TracersMonthlySurfClimoCell(L, ICell) - - TracerCell(L, ICell, K)) < -MaxDiff) { - SurfTracerRestValuesCell(L, ICell) = -MaxDiff; - } else { - SurfTracerRestValuesCell(L, ICell) = - TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K); - } + computeVarsOnCells(int L, int ICell, const Array3DReal &TracerCell) const { + + const I4 K = MinLayerCell(ICell); + + if ((TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K)) > + MaxDiff) { + SurfTracerRestValuesCell(L, ICell) = MaxDiff; + } else if ((TracersMonthlySurfClimoCell(L, ICell) - + TracerCell(L, ICell, K)) < -MaxDiff) { + SurfTracerRestValuesCell(L, ICell) = -MaxDiff; + } else { + SurfTracerRestValuesCell(L, ICell) = + TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K); } } diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index 44957762f953..b1d57f19fc6b 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -15,6 +15,7 @@ #include "VertCoord.h" #include "auxiliaryVars/KineticAuxVars.h" #include "auxiliaryVars/LayerThicknessAuxVars.h" +#include "auxiliaryVars/SurfTracerRestAuxVars.h" #include "auxiliaryVars/TracerAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" @@ -68,6 +69,8 @@ struct TestSetupPlane { ErrorMeasures ExpectedNormalStressErrors = {0.0033910709836867704, 0.0039954090464502795}; + ErrorMeasures ExpectedSurfTracerRestErrors = {0, 0}; + KOKKOS_FUNCTION Real layerThickness(Real X, Real Y) const { return 2 + std::cos(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } @@ -196,6 +199,8 @@ struct TestSetupSphere { ErrorMeasures ExpectedNormalStressErrors = {0.0038588958862868362, 0.003813760171030077}; + ErrorMeasures ExpectedSurfTracerRestErrors = {0, 0}; + KOKKOS_FUNCTION Real layerThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); } @@ -777,6 +782,90 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, return Err; } + +int testSurfTracerRestAuxVars(int NVertLayers, int NTracers, Real RTol) { + int Err = 0; + TestSetup Setup; + + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + SurfTracerRestAuxVars SurfTracerRestAux("", Mesh, VCoord, NTracers); + SurfTracerRestAux.MaxDiff = 0.5; // choose a small max‑difference so all + // three branches (>MaxDiff, <-MaxDiff, + // inside) are exercised. + + Array2DReal InputField("InputField", Mesh->NCellsSize, 1); + deepCopy(InputField, 0.0); + + // fill monthly climo with a simple analytic function that produces + // positive/negative variations, e.g. tracer+velocityx: + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { + return Setup.tracer(X, Y) + Setup.velocityX(X, Y); + }, + InputField, Geom, Mesh, OnCell); + + parallelFor( + {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { + SurfTracerRestAux.TracersMonthlySurfClimoCell(L, ICell) = + InputField(ICell, 0); + }); + + Array3DReal TracersOnCell("TracersOnCell", NTracers, Mesh->NCellsSize, + NVertLayers); + deepCopy(TracersOnCell, 0.0); + // fill the tracer array and add a trivial vertical dependence so that + // the surface value differs from deeper layers + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, + TracersOnCell, Geom, Mesh, OnCell); + parallelFor( + {NTracers, Mesh->NCellsSize, NVertLayers}, + KOKKOS_LAMBDA(int L, int ICell, int K) { + TracersOnCell(L, ICell, K) += 0.04_Real * K; + }); + + // compute the numerical result + parallelFor( + {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { + SurfTracerRestAux.computeVarsOnCells(L, ICell, TracersOnCell); + }); + // const auto &NumSurfRest = SurfTracerRestAux.SurfTracerRestValuesCell; + + // build the exact result by re‑implementing the class logic, including + // lookup of MinLayerCell + Array1DReal ExactSurfRest("ExactSurfRest", Mesh->NCellsOwned); + Array1DReal NumSurfRest("NumSurfRest", Mesh->NCellsOwned); + const auto &MinLayerCell = VCoord->MinLayerCell; + parallelFor( + {Mesh->NCellsOwned}, KOKKOS_LAMBDA(int ICell) { + const int K = MinLayerCell(ICell); + const Real diff = + SurfTracerRestAux.TracersMonthlySurfClimoCell(0, ICell) - + TracersOnCell(0, ICell, K); + if (diff > SurfTracerRestAux.MaxDiff) { + ExactSurfRest(ICell) = SurfTracerRestAux.MaxDiff; + } else if (diff < -SurfTracerRestAux.MaxDiff) { + ExactSurfRest(ICell) = -SurfTracerRestAux.MaxDiff; + } else { + ExactSurfRest(ICell) = diff; + } + NumSurfRest(ICell) = + SurfTracerRestAux.SurfTracerRestValuesCell(0, ICell); + }); + + ErrorMeasures SurfRestErrors; + Err += + computeErrors(SurfRestErrors, NumSurfRest, ExactSurfRest, Mesh, OnCell); + Err += checkErrors("AuxVarsTest", "SurfTracerRest", SurfRestErrors, + Setup.ExpectedSurfTracerRestErrors, RTol); + + if (Err == 0) { + LOG_INFO("AuxVarsTest: SurfTracerRestAuxVars PASS"); + } + return Err; +} //------------------------------------------------------------------------------ // The initialization routine for aux vars testing int initAuxVarsTest(const std::string &mesh) { @@ -847,6 +936,8 @@ int auxVarsTest(const std::string &mesh = DefaultMeshFile) { Err += testWindForcingAuxVars(RTol); + Err += testSurfTracerRestAuxVars(NVertLayers, NTracers, RTol); + if (Err == 0) { LOG_INFO("AuxVarsTest: Successful completion"); } diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 980a7d759dcb..93d8003e50cd 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -58,6 +58,7 @@ struct TestSetupPlane { 0.00290978146207349032}; ErrorMeasures ExpectedTrDel4Errors = {0.00508833446725232875, 0.00523080740758275625}; + ErrorMeasures ExpectedSurfTrRestErrors = {0, 0}; ErrorMeasures ExpectedWindForcingErrors = {0, 0}; ErrorMeasures ExpectedBottomDragErrors = {0.033848740052302935, 0.01000133508329411}; @@ -197,6 +198,7 @@ struct TestSetupSphere { 0.005105510870642706}; ErrorMeasures ExpectedTrDel4Errors = {0.0008646345116716073, 0.0007118574326665881}; + ErrorMeasures ExpectedSurfTrRestErrors = {0, 0}; ErrorMeasures ExpectedWindForcingErrors = {0, 0}; ErrorMeasures ExpectedBottomDragErrors = {0.0015333449035655053, 0.0014897009917655022}; @@ -1011,6 +1013,80 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) { return Err; } // end testTracerHyperDiffOnCell +int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { + + I4 Err = 0; + TestSetup Setup; + + const auto Mesh = HorzMesh::getDefault(); + const auto VCoord = VertCoord::getDefault(); + + // Compute exact result: zero everywhere except at surface layer where it + // equals PistonVelocity * input field + Array3DReal ExactSurfRest("ExactSurfRest", NTracers, Mesh->NCellsOwned, + NVertLayers); + // Initialize exact array to zero + deepCopy(ExactSurfRest, 0); + + // Create temporary array for input field (scalarA) + Array2DReal InputField("InputField", Mesh->NCellsSize, 1); + // Set input array (surface tracer restoring values at each cell) + Array2DReal SurfTracerRestValuesCell("SurfTracerRestValuesCell", NTracers, + Mesh->NCellsSize, NVertLayers); + deepCopy(InputField, 0); + deepCopy(SurfTracerRestValuesCell, 0); + + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, + InputField, Geom, Mesh, OnCell); + + // Set exact array at surface layer only + // The functor modifies Tend(L, ICell, MinLayerCell(ICell)) when KChunk==0 + // We need access to MinLayerCell to set exact values correctly + // For now, use a straightforward approach: set all MinLayer indices to 0 + // (topmost layer) in the numerical computation loop + SurfaceTracerRestoringOnCell SurfRestOnC(Mesh, VCoord); + + // Set the exact array at the minimum layer for each cell + // For this test we assume the minimum layer index is 0 (surface). The + // exact result will therefore only be non-zero at KLayer==0 and equal to + // PistonVelocity * input field. Perform the computation directly on the + // device rather than staging data on the host. + parallelFor( + {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { + SurfTracerRestValuesCell(L, ICell) = InputField(ICell, 0); + ExactSurfRest(L, ICell, 0) = + SurfRestOnC.PistonVelocity * InputField(ICell, 0); + }); + + // Compute numerical result: initialize to zero and apply functor + Array3DReal NumSurfRest("NumSurfRest", NTracers, Mesh->NCellsOwned, + NVertLayers); + deepCopy(NumSurfRest, 0); + + // Loop over all KChunk values to verify functor only modifies at KChunk==0 + parallelFor( + {NTracers, Mesh->NCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(int L, int ICell, int KLayer) { + SurfRestOnC(NumSurfRest, L, ICell, KLayer, SurfTracerRestValuesCell); + }); + + // Compute errors + ErrorMeasures SurfRestErrors; + Err += + computeErrors(SurfRestErrors, NumSurfRest, ExactSurfRest, Mesh, OnCell); + + // Check error values (expect zero error) + Err += checkErrors("TendencyTermsTest", "SurfaceTracerRestoring", + SurfRestErrors, Setup.ExpectedSurfTrRestErrors, RTol); + + if (Err == 0) { + LOG_INFO("TendencyTermsTest: SurfaceTracerRestoring PASS"); + } + + return Err; +} // end testSurfaceTracerRestoringOnCell + void initTendTest(const std::string &MeshFile, int NVertLayers) { Error Err; @@ -1091,6 +1167,8 @@ int tendencyTermsTest(const std::string &MeshFile = DefaultMeshFile) { Err += testTracerHyperDiffOnCell(NVertLayers, NTracers, RTol); + Err += testSurfaceTracerRestoringOnCell(NVertLayers, NTracers, RTol); + if (Err == 0) { LOG_INFO("TendencyTermsTest: Successful completion"); } From f142ae5f0204d21e4d012905ad4fe2893bf0ae56 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Wed, 4 Mar 2026 08:51:53 -0800 Subject: [PATCH 03/15] rebase --- components/omega/src/ocn/AuxiliaryState.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 6031f7898658..93e8c6979bb6 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -224,11 +224,9 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, void AuxiliaryState::computeAll(const OceanState *State, const Array3DReal &TracerArray, int ThickTimeLevel, int VelTimeLevel) const { - Array2DReal LayerThickCell; - Array2DReal NormalVelEdge; - State->getLayerThickness(LayerThickCell, ThickTimeLevel); - State->getNormalVelocity(NormalVelEdge, VelTimeLevel); - + Array2DReal LayerThickCell = State->getLayerThickness(ThickTimeLevel); + Array2DReal NormalVelEdge = State->getNormalVelocity(VelTimeLevel); + const int NTracers = TracerArray.extent_int(0); OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); From 864cf69092228549b14069bdcf87d8fcc025b2bc Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Wed, 4 Mar 2026 09:37:41 -0800 Subject: [PATCH 04/15] updates docs --- components/omega/doc/devGuide/AuxiliaryVariables.md | 2 ++ components/omega/doc/devGuide/TendencyTerms.md | 1 + components/omega/doc/userGuide/AuxiliaryVariables.md | 2 ++ components/omega/doc/userGuide/TendencyTerms.md | 2 ++ 4 files changed, 7 insertions(+) diff --git a/components/omega/doc/devGuide/AuxiliaryVariables.md b/components/omega/doc/devGuide/AuxiliaryVariables.md index d75a8f2b3f84..e899e20e2812 100644 --- a/components/omega/doc/devGuide/AuxiliaryVariables.md +++ b/components/omega/doc/devGuide/AuxiliaryVariables.md @@ -104,3 +104,5 @@ The following auxiliary variable groups are currently implemented: | WindForcingAuxVars | ZonalStressCell || || MeridStressCell || || NormalStressEdge || +| SurfTracerRestAuxVars | SurfTracerRestValuesCell || +|| TracersMonthlySurfClimoCell || diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index e665fb03ed76..5d8770ab74d6 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -41,3 +41,4 @@ implemented: - `TracerHighOrderHorzAdvOnCell` - `TracerDiffOnCell` - `TracerHyperDiffOnCell` +- `SurfaceTracerRestoringOnCell` \ No newline at end of file diff --git a/components/omega/doc/userGuide/AuxiliaryVariables.md b/components/omega/doc/userGuide/AuxiliaryVariables.md index 4f6fbc18c616..ca441e0b7e6c 100644 --- a/components/omega/doc/userGuide/AuxiliaryVariables.md +++ b/components/omega/doc/userGuide/AuxiliaryVariables.md @@ -32,3 +32,5 @@ The following auxiliary variables are currently available: | ZonalStressCell | zonal component of wind stress on cells | MeridStressCell | meridional component of wind stress on cells | NormalStressEdge | normal component of wind stress on edge +| SurfTracerRestValuesCell | surface tracer restoring values on cells +| TracersMonthlySurfClimoCell | monthly climatology value for surface tracer on cells diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 5757e98ebe0f..66035b3ed135 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -20,6 +20,7 @@ tendency terms are currently implemented: | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers | WindForcingOnEdge | forcing by wind stress, defined on edges | BottomDragOnEdge | bottom drag, defined on edges +| SurfaceTracerRestoringOnCell | surface tracer restoring, defined on cells Among the internal data stored by each functor is a `bool` which can enable or disable the contribution of that particular term to the tendency. These flags @@ -54,6 +55,7 @@ the currently available tendency terms: | WindForcingOnEdge | WindForcingTendencyEnable | enable/disable term | BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term | | BottomDragCoeff | bottom drag coefficient +| SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term ## Second Order Horizontal Advection Algorithm From 8777f79b383939e82add36097a7c0537e9a7cdce Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Fri, 20 Mar 2026 14:40:54 -0700 Subject: [PATCH 05/15] adds ability to apply surface restoring to specified tracers or tracer groups --- components/omega/configs/Default.yml | 2 + components/omega/src/ocn/TendencyTerms.h | 4 +- .../omega/test/ocn/TendencyTermsTest.cpp | 116 ++++++++++-------- 3 files changed, 68 insertions(+), 54 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 7d166110429a..f8a701ad140b 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -34,6 +34,8 @@ Omega: SurfaceRestoring: MaxDiff: 100.0 PistonVelocity: 1.585e-5 + TracerGroupsToRestore: [Base] + TracersToRestore: [] VertCoord: MovementWeightType: Uniform PressureGrad: diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 8e92b842646e..69f7e39ec69d 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -575,6 +575,7 @@ class SurfaceTracerRestoringOnCell { public: bool Enabled; Real PistonVelocity = 1.585e-5; + Array1DI4 TracerRestoreMask; /// constructor declaration SurfaceTracerRestoringOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); @@ -584,7 +585,8 @@ class SurfaceTracerRestoringOnCell { KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, const Array2DReal &SurfTracerRestValuesCell) const { - if (KChunk == 0) { + if (KChunk == 0 && + (TracerRestoreMask.extent(0) == 0 || TracerRestoreMask(L) == 1)) { const I4 K = MinLayerCell(ICell); Tend(L, ICell, K) += diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 93d8003e50cd..69f585d0e90d 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -1021,64 +1021,74 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { const auto Mesh = HorzMesh::getDefault(); const auto VCoord = VertCoord::getDefault(); - // Compute exact result: zero everywhere except at surface layer where it - // equals PistonVelocity * input field - Array3DReal ExactSurfRest("ExactSurfRest", NTracers, Mesh->NCellsOwned, - NVertLayers); - // Initialize exact array to zero - deepCopy(ExactSurfRest, 0); - - // Create temporary array for input field (scalarA) - Array2DReal InputField("InputField", Mesh->NCellsSize, 1); - // Set input array (surface tracer restoring values at each cell) - Array2DReal SurfTracerRestValuesCell("SurfTracerRestValuesCell", NTracers, - Mesh->NCellsSize, NVertLayers); - deepCopy(InputField, 0); - deepCopy(SurfTracerRestValuesCell, 0); + if (NTracers < 3) { + LOG_ERROR("TendencyTermsTest: SurfaceTracerRestoring requires at least 3 " + "tracers, found {}", + NTracers); + return 1; + } - Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, - InputField, Geom, Mesh, OnCell); - - // Set exact array at surface layer only - // The functor modifies Tend(L, ICell, MinLayerCell(ICell)) when KChunk==0 - // We need access to MinLayerCell to set exact values correctly - // For now, use a straightforward approach: set all MinLayer indices to 0 - // (topmost layer) in the numerical computation loop - SurfaceTracerRestoringOnCell SurfRestOnC(Mesh, VCoord); - - // Set the exact array at the minimum layer for each cell - // For this test we assume the minimum layer index is 0 (surface). The - // exact result will therefore only be non-zero at KLayer==0 and equal to - // PistonVelocity * input field. Perform the computation directly on the - // device rather than staging data on the host. - parallelFor( - {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { - SurfTracerRestValuesCell(L, ICell) = InputField(ICell, 0); - ExactSurfRest(L, ICell, 0) = - SurfRestOnC.PistonVelocity * InputField(ICell, 0); - }); + const char *CaseLabels[3] = {"SurfaceTracerRestoringSalinityOnly", + "SurfaceTracerRestoringTemperatureOnly", + "SurfaceTracerRestoringTempSaltDebug"}; - // Compute numerical result: initialize to zero and apply functor - Array3DReal NumSurfRest("NumSurfRest", NTracers, Mesh->NCellsOwned, - NVertLayers); - deepCopy(NumSurfRest, 0); + const I4 CaseMasks[3][3] = { + {0, 1, 0}, + {1, 0, 0}, + {1, 1, 1}, + }; - // Loop over all KChunk values to verify functor only modifies at KChunk==0 - parallelFor( - {NTracers, Mesh->NCellsOwned, NVertLayers}, - KOKKOS_LAMBDA(int L, int ICell, int KLayer) { - SurfRestOnC(NumSurfRest, L, ICell, KLayer, SurfTracerRestValuesCell); - }); + for (int Case = 0; Case < 3; ++Case) { - // Compute errors - ErrorMeasures SurfRestErrors; - Err += - computeErrors(SurfRestErrors, NumSurfRest, ExactSurfRest, Mesh, OnCell); + Array3DReal ExactSurfRest("ExactSurfRest", NTracers, Mesh->NCellsOwned, + NVertLayers); + Array2DReal InputField("InputField", Mesh->NCellsSize, 1); + Array2DReal SurfTracerRestValuesCell("SurfTracerRestValuesCell", NTracers, + Mesh->NCellsSize, NVertLayers); + Array3DReal NumSurfRest("NumSurfRest", NTracers, Mesh->NCellsOwned, + NVertLayers); - // Check error values (expect zero error) - Err += checkErrors("TendencyTermsTest", "SurfaceTracerRestoring", - SurfRestErrors, Setup.ExpectedSurfTrRestErrors, RTol); + deepCopy(ExactSurfRest, 0); + deepCopy(InputField, 0); + deepCopy(SurfTracerRestValuesCell, 0); + deepCopy(NumSurfRest, 0); + + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, + InputField, Geom, Mesh, OnCell); + + SurfaceTracerRestoringOnCell SurfRestOnC(Mesh, VCoord); + SurfRestOnC.TracerRestoreMask = Array1DI4("TracerRestoreMask", NTracers); + + HostArray1DI4 TracerRestoreMaskH("TracerRestoreMaskH", NTracers); + deepCopy(TracerRestoreMaskH, 0); + for (int L = 0; L < 3; ++L) { + TracerRestoreMaskH(L) = CaseMasks[Case][L]; + } + deepCopy(SurfRestOnC.TracerRestoreMask, TracerRestoreMaskH); + + parallelFor( + {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { + SurfTracerRestValuesCell(L, ICell) = InputField(ICell, 0); + ExactSurfRest(L, ICell, 0) = SurfRestOnC.TracerRestoreMask(L) * + SurfRestOnC.PistonVelocity * + InputField(ICell, 0); + }); + + parallelFor( + {NTracers, Mesh->NCellsOwned, NVertLayers}, + KOKKOS_LAMBDA(int L, int ICell, int KLayer) { + SurfRestOnC(NumSurfRest, L, ICell, KLayer, + SurfTracerRestValuesCell); + }); + + ErrorMeasures SurfRestErrors; + Err += computeErrors(SurfRestErrors, NumSurfRest, ExactSurfRest, Mesh, + OnCell); + + Err += checkErrors("TendencyTermsTest", CaseLabels[Case], SurfRestErrors, + Setup.ExpectedSurfTrRestErrors, RTol); + } if (Err == 0) { LOG_INFO("TendencyTermsTest: SurfaceTracerRestoring PASS"); From 4367a5e8ffac49e0fdbf1d614d25a042208c7b85 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 23 Mar 2026 09:17:28 -0700 Subject: [PATCH 06/15] removed group tracer optioin in config and checks if maxdiff is pos --- components/omega/configs/Default.yml | 1 - components/omega/src/ocn/AuxiliaryState.cpp | 4 ++++ components/omega/src/ocn/Tendencies.cpp | 2 +- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index f8a701ad140b..748f289c04a8 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -34,7 +34,6 @@ Omega: SurfaceRestoring: MaxDiff: 100.0 PistonVelocity: 1.585e-5 - TracerGroupsToRestore: [Base] TracersToRestore: [] VertCoord: MovementWeightType: Uniform diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 93e8c6979bb6..ecd45fb0bcf7 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -414,6 +414,10 @@ void AuxiliaryState::readConfigOptions(Config *OmegaConfig) { Err += SurfRestConfig.get("MaxDiff", this->SurfTracerRestAux.MaxDiff); CHECK_ERROR_ABORT( Err, "AuxiliaryState: MaxDiff not found in SurfaceRestoringConfig"); + if (this->SurfTracerRestAux.MaxDiff <= 0) { + ABORT_ERROR("AuxiliaryState: MaxDiff must be positive, got {}", + this->SurfTracerRestAux.MaxDiff); + } } //------------------------------------------------------------------------------ diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index b3e6de666024..a56216be42a5 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -257,7 +257,7 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config Err += TendConfig.get("PressureGradTendencyEnable", this->PGrad->Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); - + Err += TendConfig.get("SurfaceTracerRestoringEnable", this->SurfaceTracerRestoring.Enabled); CHECK_ERROR_ABORT( From 4189450194175019513fbee01c39080089835585 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Tue, 24 Mar 2026 11:08:01 -0700 Subject: [PATCH 07/15] adds more documentation and renames vars --- components/omega/configs/Default.yml | 4 +- .../omega/doc/devGuide/AuxiliaryVariables.md | 7 +- components/omega/doc/devGuide/Forcing.md | 79 +++++++++++++++++ .../omega/doc/devGuide/TendencyTerms.md | 8 +- components/omega/doc/index.md | 5 ++ .../omega/doc/userGuide/AuxiliaryVariables.md | 9 +- components/omega/doc/userGuide/Forcing.md | 84 +++++++++++++++++++ .../omega/doc/userGuide/TendencyTerms.md | 6 +- components/omega/src/ocn/AuxiliaryState.cpp | 9 ++ components/omega/src/ocn/TendencyTerms.h | 4 +- .../auxiliaryVars/SurfTracerRestAuxVars.cpp | 51 +++++++---- .../ocn/auxiliaryVars/SurfTracerRestAuxVars.h | 8 +- .../omega/test/ocn/AuxiliaryVarsTest.cpp | 3 +- .../omega/test/ocn/TendencyTermsTest.cpp | 11 +-- 14 files changed, 251 insertions(+), 37 deletions(-) create mode 100644 components/omega/doc/devGuide/Forcing.md create mode 100644 components/omega/doc/userGuide/Forcing.md diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 748f289c04a8..d1222da0dc39 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -32,9 +32,9 @@ Omega: WindStress: InterpType: Isotropic SurfaceRestoring: - MaxDiff: 100.0 + TracersToRestore: [Temperature, Salinity] PistonVelocity: 1.585e-5 - TracersToRestore: [] + MaxDiff: 100.0 VertCoord: MovementWeightType: Uniform PressureGrad: diff --git a/components/omega/doc/devGuide/AuxiliaryVariables.md b/components/omega/doc/devGuide/AuxiliaryVariables.md index e899e20e2812..de96e3186921 100644 --- a/components/omega/doc/devGuide/AuxiliaryVariables.md +++ b/components/omega/doc/devGuide/AuxiliaryVariables.md @@ -104,5 +104,10 @@ The following auxiliary variable groups are currently implemented: | WindForcingAuxVars | ZonalStressCell || || MeridStressCell || || NormalStressEdge || -| SurfTracerRestAuxVars | SurfTracerRestValuesCell || +| SurfTracerRestAuxVars | SurfTracerRestoringDiffsCell || || TracersMonthlySurfClimoCell || + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-dev-forcing). diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md new file mode 100644 index 000000000000..d1cdbf284ebd --- /dev/null +++ b/components/omega/doc/devGuide/Forcing.md @@ -0,0 +1,79 @@ +(omega-dev-forcing)= + +# Forcing + +This page describes design and implementation details for forcing-related +pathways in Omega, currently this includes: + +- Wind forcing +- Surface tracer restoring + +## Wind forcing design + +### Wind forcing data flow + +1. External fields provide: + - `WindStressZonal` + - `WindStressMeridional` +2. Auxiliary-state compute builds `NormalStressEdge` from those fields. +3. Tendency term applies wind-stress forcing to edge-normal velocity tendency. + +### Wind forcing key classes/components + +- `WindForcingAuxVars` + - Stores wind-stress cell fields and computed `NormalStressEdge` + - Applies configured interpolation choice (`InterpType`) +- `AuxiliaryState::computeMomAux` + - Calls `WindForcingAuxVars::computeVarsOnEdge` +- `WindForcingOnEdge` tendency term + - Adds contribution proportional to normal stress and inverse layer + thickness in the surface layer + +### Wind forcing config coupling + +- `Omega.WindStress.InterpType` + - mapped to `InterpCellToEdgeOption` +- `Omega.Tendencies.WindForcingTendencyEnable` + - gates execution of wind forcing tendency kernel + +## Surface tracer restoring design + +### Surface tracer restoring data flow + +1. External fields provide target values: `TracersMonthlySurfClimoCell` +2. Auxiliary-state compute forms bounded restoring differences: + - `SurfTracerRestoringDiffsCell = clamp(target - tracer_surface, +/- MaxDiff)` +3. Tendency term applies restoring only at surface layer and only for tracers + enabled by restore mask. + +### Surface tracer restoring key classes/components + +- `SurfTracerRestAuxVars` + - Inputs: `TracersMonthlySurfClimoCell`, tracer state array + - Output: `SurfTracerRestoringDiffsCell` + - Uses `MinLayerCell` to select surface layer index +- `SurfaceTracerRestoringOnCell` tendency term + - Applies `PistonVelocity * SurfTracerRestoringDiffsCell` at surface +- `Tendencies` + - Builds restore mask from `SurfaceRestoring.TracersToRestore` + - Calls surface restoring term when enabled + +### Surface tracer restoring config coupling + +- `Omega.SurfaceRestoring.MaxDiff` + - bound for restoring differences +- `Omega.SurfaceRestoring.PistonVelocity` + - tendency scaling +- `Omega.SurfaceRestoring.TracersToRestore` + - tracer-level enable list +- `Omega.Tendencies.SurfaceTracerRestoringEnable` + - gates restoring tendency execution + + ## Notes + +- If a tracer is not listed in `TracersToRestore`, no restoring tendency is + applied to that tracer. +- `MaxDiff` must be positive. A runtime check will error out if not. +- It is assumed that the incoming `TracersMonthlySurfClimoCell` values have the correct units appropriate for Omega. If this is not true, conversion should be implemented. +- Surface tracer restoring is active everywhere if enabled. A flag to turn it off under sea ice will need to be added in later development if this feature is desired. +- A global sum has not been implemented for the surface tracer restoring, but should be added in later development. diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 5d8770ab74d6..01641a60685a 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -41,4 +41,10 @@ implemented: - `TracerHighOrderHorzAdvOnCell` - `TracerDiffOnCell` - `TracerHyperDiffOnCell` -- `SurfaceTracerRestoringOnCell` \ No newline at end of file +- `SurfaceTracerRestoringOnCell` + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-dev-forcing). +>>>>>>> 41ae670f82 (adds more documentation and renames vars) diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 220965ff1fdb..aac646b407ad 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -53,6 +53,7 @@ userGuide/PGrad userGuide/Timing userGuide/VerticalMixingCoeff userGuide/VertAdv +userGuide/Forcing ``` ```{toctree} @@ -97,7 +98,11 @@ devGuide/VertCoord devGuide/PGrad devGuide/Timing devGuide/VerticalMixingCoeff +<<<<<<< HEAD devGuide/VertAdv +======= +devGuide/Forcing +>>>>>>> 41ae670f82 (adds more documentation and renames vars) ``` ```{toctree} diff --git a/components/omega/doc/userGuide/AuxiliaryVariables.md b/components/omega/doc/userGuide/AuxiliaryVariables.md index ca441e0b7e6c..f6e377c388f9 100644 --- a/components/omega/doc/userGuide/AuxiliaryVariables.md +++ b/components/omega/doc/userGuide/AuxiliaryVariables.md @@ -32,5 +32,10 @@ The following auxiliary variables are currently available: | ZonalStressCell | zonal component of wind stress on cells | MeridStressCell | meridional component of wind stress on cells | NormalStressEdge | normal component of wind stress on edge -| SurfTracerRestValuesCell | surface tracer restoring values on cells -| TracersMonthlySurfClimoCell | monthly climatology value for surface tracer on cells +| SurfTracerRestoringDiffsCell | surface tracer restoring differences on cells +| TracersMonthlySurfClimoCell | monthly climatology values to restore to for surface tracer on cells + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-user-forcing). diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md new file mode 100644 index 000000000000..5c973801a2a2 --- /dev/null +++ b/components/omega/doc/userGuide/Forcing.md @@ -0,0 +1,84 @@ +(omega-user-forcing)= + +# Forcing + +This page documents the user-facing configuration and behavior for current forcing in Omega: + +- Wind forcing +- Surface tracer restoring + +## Wind forcing + +Wind forcing adds momentum tendency from surface wind stress. + +### Wind forcing configuration + +Wind forcing behavior is controlled by two configuration blocks: + +```yaml +Omega: + WindStress: + InterpType: Isotropic + + Tendencies: + WindForcingTendencyEnable: true +``` + +- `WindStress.InterpType` + - `Isotropic`: isotropic cell-to-edge interpolation for wind stress + - `Anisotropic`: anisotropic interpolation option +- `Tendencies.WindForcingTendencyEnable`: switch to enable wind forcing tendency + +### Required input fields + +Wind forcing uses auxiliary wind-stress fields: + +- `WindStressZonal` +- `WindStressMeridional` + +These are used to form edge-normal stress (`NormalStressEdge`) that enters +momentum tendencies. + +## Surface tracer restoring + +Surface tracer restoring applies a piston-velocity tendency, or damping, at the ocean +surface for selected tracers. This is implemented to mitigate drifts in chosen tracers +(most often salinity) by nudging the model's simulated tracer values towards observed climatological values. +This process prevents oceanic regimes from shifting away from reality due to errors in surface freshwater +forcing (in the case of salinity restoring). Currently, it is applied everywhere when enabled. + +### Surface tracer restoring configuration + +Surface tracer restoring is controlled by two configuration blocks: + +```yaml +Omega: + SurfaceRestoring: + TracersToRestore: [Temperature, Salinity] + PistonVelocity: 1.585e-5 + MaxDiff: 100.0 + + Tendencies: + SurfaceTracerRestoringEnable: true +``` + +- `TracersToRestore`: list of tracer names that restoring is applied to +- `PistonVelocity`: restoring rate coefficient +- `MaxDiff`: cap on restoring difference magnitude +- `SurfaceTracerRestoringEnable`: switch to enable surface tracer restoring + +### Restoring target fields + +Surface restoring uses auxiliary fields: + +- `TracersMonthlySurfClimoCell`: restoring target climatological values +- `SurfTracerRestoringDiffsCell`: computed target-minus-state differences + +The restoring tendency is computed at the surface layer only and is limited by +`MaxDiff`. + +## Notes + +- If a tracer is not listed in `TracersToRestore`, no restoring tendency is + applied to that tracer. +- `MaxDiff` must be positive. A runtime check will error out if not. diff --git a/components/omega/doc/userGuide/TendencyTerms.md b/components/omega/doc/userGuide/TendencyTerms.md index 66035b3ed135..4dbcfcb13519 100644 --- a/components/omega/doc/userGuide/TendencyTerms.md +++ b/components/omega/doc/userGuide/TendencyTerms.md @@ -57,7 +57,6 @@ the currently available tendency terms: | | BottomDragCoeff | bottom drag coefficient | SurfaceTracerRestoringOnCell | SurfaceTracerRestoringEnable | enable/disable term - ## Second Order Horizontal Advection Algorithm The horizontal advection is done independently within each ocean layer @@ -138,3 +137,8 @@ is about order 1.7 as shown in {numref}`tracer-higher-order-convergence`: :width: 600 px Tracer higer order convergence example of a cosine bell advected on a sphere showing an order 1.71 convergence rate ``` + +## See Also + +Additional information on forcing (currently wind forcing and surface tracer +restoring) is detailed in [](omega-user-forcing). diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index ecd45fb0bcf7..a306d12151cf 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -431,6 +431,15 @@ I4 AuxiliaryState::exchangeHalo() { Err += MeshHalo->exchangeFullArrayHalo(WindForcingAux.MeridStressCell, OnCell); + const I4 NTracers = + SurfTracerRestAux.TracersMonthlySurfClimoCell.extent_int(0); + for (I4 LTracer = 0; LTracer < NTracers; ++LTracer) { + auto TracerSurfClimoCell = + Kokkos::subview(SurfTracerRestAux.TracersMonthlySurfClimoCell, + LTracer, Kokkos::ALL()); + Err += MeshHalo->exchangeFullArrayHalo(TracerSurfClimoCell, OnCell); + } + return Err; } // end exchangeHalo diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 69f7e39ec69d..63d570a9b39d 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -584,13 +584,13 @@ class SurfaceTracerRestoringOnCell { /// restoring values, outputs tendency array KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, - const Array2DReal &SurfTracerRestValuesCell) const { + const Array2DReal &SurfTracerRestoringDiffsCell) const { if (KChunk == 0 && (TracerRestoreMask.extent(0) == 0 || TracerRestoreMask(L) == 1)) { const I4 K = MinLayerCell(ICell); Tend(L, ICell, K) += - PistonVelocity * SurfTracerRestValuesCell(L, ICell); + PistonVelocity * SurfTracerRestoringDiffsCell(L, ICell); } } diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp index 3837c131f9e0..5d61813f38b7 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp @@ -10,8 +10,9 @@ SurfTracerRestAuxVars::SurfTracerRestAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, const VertCoord *VCoord, const I4 NTracers) - : SurfTracerRestValuesCell("SurfTracerRestValuesCell" + AuxStateSuffix, - NTracers, Mesh->NCellsSize), + : SurfTracerRestoringDiffsCell("SurfTracerRestoringDiffsCell" + + AuxStateSuffix, + NTracers, Mesh->NCellsSize), TracersMonthlySurfClimoCell("TracersMonthlySurfClimoCell" + AuxStateSuffix, NTracers, Mesh->NCellsSize), @@ -34,30 +35,46 @@ void SurfTracerRestAuxVars::registerFields( } // Tracers surface restoring values - DimNames[0] = "NTracers"; - DimNames[1] = "NCells" + DimSuffix; - auto SurfTracerRestValuesCellField = - Field::create(SurfTracerRestValuesCell.label(), // field name - "tracer surface restoring value", // long name/describe - "tracer units", // units - "", // CF standard Name - std::numeric_limits::min(), // min valid value - std::numeric_limits::max(), // max valid value + DimNames[0] = "NTracers"; + DimNames[1] = "NCells" + DimSuffix; + auto SurfTracerRestoringDiffsCellField = Field::create( + SurfTracerRestoringDiffsCell.label(), // field name + "surface tracer restoring differences", // long name/describe + "tracer units", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value + FillValue, // scalar for undefined entries + NDims, // number of dimensions + DimNames); // dim names + + auto TracersMonthlySurfClimoCellField = + Field::create(TracersMonthlySurfClimoCell.label(), + "monthly surface tracer climatology", // long name/describe + "tracer units", // units + "", // CF standard Name + std::numeric_limits::min(), // min valid value + std::numeric_limits::max(), // max valid value FillValue, // scalar for undefined entries NDims, // number of dimensions - DimNames // dim names - ); + DimNames); // dim names // Add fields to FieldGroup - FieldGroup::addFieldToGroup(SurfTracerRestValuesCell.label(), AuxGroupName); + FieldGroup::addFieldToGroup(SurfTracerRestoringDiffsCell.label(), + AuxGroupName); + FieldGroup::addFieldToGroup(TracersMonthlySurfClimoCell.label(), + AuxGroupName); // Attach data - SurfTracerRestValuesCellField->attachData( - SurfTracerRestValuesCell); + SurfTracerRestoringDiffsCellField->attachData( + SurfTracerRestoringDiffsCell); + TracersMonthlySurfClimoCellField->attachData( + TracersMonthlySurfClimoCell); } void SurfTracerRestAuxVars::unregisterFields() const { - Field::destroy(SurfTracerRestValuesCell.label()); + Field::destroy(SurfTracerRestoringDiffsCell.label()); + Field::destroy(TracersMonthlySurfClimoCell.label()); } } // namespace OMEGA diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h index 2411eecaf75e..6b77c2237699 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h @@ -13,7 +13,7 @@ namespace OMEGA { class SurfTracerRestAuxVars { public: - Array2DReal SurfTracerRestValuesCell; + Array2DReal SurfTracerRestoringDiffsCell; Array2DReal TracersMonthlySurfClimoCell; Real MaxDiff = 100.0; // maximum allowed difference for restoring /// Need to add under sea ice restoring option when that is available @@ -29,12 +29,12 @@ class SurfTracerRestAuxVars { if ((TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K)) > MaxDiff) { - SurfTracerRestValuesCell(L, ICell) = MaxDiff; + SurfTracerRestoringDiffsCell(L, ICell) = MaxDiff; } else if ((TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K)) < -MaxDiff) { - SurfTracerRestValuesCell(L, ICell) = -MaxDiff; + SurfTracerRestoringDiffsCell(L, ICell) = -MaxDiff; } else { - SurfTracerRestValuesCell(L, ICell) = + SurfTracerRestoringDiffsCell(L, ICell) = TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K); } } diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index b1d57f19fc6b..f62856065bf3 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -831,7 +831,6 @@ int testSurfTracerRestAuxVars(int NVertLayers, int NTracers, Real RTol) { {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { SurfTracerRestAux.computeVarsOnCells(L, ICell, TracersOnCell); }); - // const auto &NumSurfRest = SurfTracerRestAux.SurfTracerRestValuesCell; // build the exact result by re‑implementing the class logic, including // lookup of MinLayerCell @@ -852,7 +851,7 @@ int testSurfTracerRestAuxVars(int NVertLayers, int NTracers, Real RTol) { ExactSurfRest(ICell) = diff; } NumSurfRest(ICell) = - SurfTracerRestAux.SurfTracerRestValuesCell(0, ICell); + SurfTracerRestAux.SurfTracerRestoringDiffsCell(0, ICell); }); ErrorMeasures SurfRestErrors; diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 69f585d0e90d..835b05becd77 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -1043,14 +1043,15 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { Array3DReal ExactSurfRest("ExactSurfRest", NTracers, Mesh->NCellsOwned, NVertLayers); Array2DReal InputField("InputField", Mesh->NCellsSize, 1); - Array2DReal SurfTracerRestValuesCell("SurfTracerRestValuesCell", NTracers, - Mesh->NCellsSize, NVertLayers); + Array2DReal SurfTracerRestoringDiffsCell("SurfTracerRestoringDiffsCell", + NTracers, Mesh->NCellsSize, + NVertLayers); Array3DReal NumSurfRest("NumSurfRest", NTracers, Mesh->NCellsOwned, NVertLayers); deepCopy(ExactSurfRest, 0); deepCopy(InputField, 0); - deepCopy(SurfTracerRestValuesCell, 0); + deepCopy(SurfTracerRestoringDiffsCell, 0); deepCopy(NumSurfRest, 0); Err += setScalar( @@ -1069,7 +1070,7 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { parallelFor( {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { - SurfTracerRestValuesCell(L, ICell) = InputField(ICell, 0); + SurfTracerRestoringDiffsCell(L, ICell) = InputField(ICell, 0); ExactSurfRest(L, ICell, 0) = SurfRestOnC.TracerRestoreMask(L) * SurfRestOnC.PistonVelocity * InputField(ICell, 0); @@ -1079,7 +1080,7 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { {NTracers, Mesh->NCellsOwned, NVertLayers}, KOKKOS_LAMBDA(int L, int ICell, int KLayer) { SurfRestOnC(NumSurfRest, L, ICell, KLayer, - SurfTracerRestValuesCell); + SurfTracerRestoringDiffsCell); }); ErrorMeasures SurfRestErrors; From c785491a02f56728ef741f32a87060d9873ff1dd Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Tue, 24 Mar 2026 14:13:42 -0700 Subject: [PATCH 08/15] removed layer loop from surface restoring calc --- components/omega/doc/devGuide/Forcing.md | 10 +++++-- components/omega/doc/userGuide/Forcing.md | 5 ++++ components/omega/src/ocn/Tendencies.cpp | 3 +- components/omega/src/ocn/TendencyTerms.h | 16 +++------- .../omega/test/ocn/TendencyTermsTest.cpp | 21 +++++--------- .../test/timeStepping/TimeStepperTest.cpp | 29 ++++++++++--------- 6 files changed, 39 insertions(+), 45 deletions(-) diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index d1cdbf284ebd..649fff6fff65 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -55,8 +55,10 @@ pathways in Omega, currently this includes: - `SurfaceTracerRestoringOnCell` tendency term - Applies `PistonVelocity * SurfTracerRestoringDiffsCell` at surface - `Tendencies` - - Builds restore mask from `SurfaceRestoring.TracersToRestore` - - Calls surface restoring term when enabled + - Builds `TracerRestoreMask` from `SurfaceRestoring.TracersToRestore` + - Applies tracer-selection logic at call site in + `computeTracerTendenciesOnly` + - Aborts if restoring is enabled but no restore mask is available ### Surface tracer restoring config coupling @@ -65,7 +67,7 @@ pathways in Omega, currently this includes: - `Omega.SurfaceRestoring.PistonVelocity` - tendency scaling - `Omega.SurfaceRestoring.TracersToRestore` - - tracer-level enable list + - tracer-level enable list used to build `TracerRestoreMask` - `Omega.Tendencies.SurfaceTracerRestoringEnable` - gates restoring tendency execution @@ -73,6 +75,8 @@ pathways in Omega, currently this includes: - If a tracer is not listed in `TracersToRestore`, no restoring tendency is applied to that tracer. +- If restoring is enabled but no restore mask is available at tendency + compute-time, Omega aborts with an error. - `MaxDiff` must be positive. A runtime check will error out if not. - It is assumed that the incoming `TracersMonthlySurfClimoCell` values have the correct units appropriate for Omega. If this is not true, conversion should be implemented. - Surface tracer restoring is active everywhere if enabled. A flag to turn it off under sea ice will need to be added in later development if this feature is desired. diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index 5c973801a2a2..a7179e392c7c 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -67,6 +67,9 @@ Omega: - `MaxDiff`: cap on restoring difference magnitude - `SurfaceTracerRestoringEnable`: switch to enable surface tracer restoring +When restoring is enabled, Omega builds an internal `TracerRestoreMask` from +`TracersToRestore` and applies restoring only to tracers whose mask value is 1. + ### Restoring target fields Surface restoring uses auxiliary fields: @@ -81,4 +84,6 @@ The restoring tendency is computed at the surface layer only and is limited by - If a tracer is not listed in `TracersToRestore`, no restoring tendency is applied to that tracer. +- If surface restoring is enabled but no restore mask is available at tendency + compute-time, Omega aborts with an error. - `MaxDiff` must be positive. A runtime check will error out if not. diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index a56216be42a5..0fc59d173c44 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -391,8 +391,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VelocityHyperDiff(Mesh, VCoord), WindForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord), TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), - TracerDiffusion(Mesh, VCoord), SurfaceTracerRestoring(Mesh, VCoord), - CustomThicknessTend(InCustomThicknessTend), + SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad) { // Tendency arrays diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 63d570a9b39d..5a06f172eb76 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -578,24 +578,16 @@ class SurfaceTracerRestoringOnCell { Array1DI4 TracerRestoreMask; /// constructor declaration - SurfaceTracerRestoringOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + SurfaceTracerRestoringOnCell(const HorzMesh *Mesh); /// The functor takes the cell index and the array for the tracer surface /// restoring values, outputs tendency array KOKKOS_FUNCTION void - operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk, + operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KMin, const Array2DReal &SurfTracerRestoringDiffsCell) const { - if (KChunk == 0 && - (TracerRestoreMask.extent(0) == 0 || TracerRestoreMask(L) == 1)) { - const I4 K = MinLayerCell(ICell); - - Tend(L, ICell, K) += - PistonVelocity * SurfTracerRestoringDiffsCell(L, ICell); - } + Tend(L, ICell, KMin) += + PistonVelocity * SurfTracerRestoringDiffsCell(L, ICell); } - - private: - Array1DI4 MinLayerCell; }; } // namespace OMEGA diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 835b05becd77..01bc834630d6 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -1058,29 +1058,22 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, InputField, Geom, Mesh, OnCell); - SurfaceTracerRestoringOnCell SurfRestOnC(Mesh, VCoord); - SurfRestOnC.TracerRestoreMask = Array1DI4("TracerRestoreMask", NTracers); - - HostArray1DI4 TracerRestoreMaskH("TracerRestoreMaskH", NTracers); - deepCopy(TracerRestoreMaskH, 0); - for (int L = 0; L < 3; ++L) { - TracerRestoreMaskH(L) = CaseMasks[Case][L]; - } - deepCopy(SurfRestOnC.TracerRestoreMask, TracerRestoreMaskH); + SurfaceTracerRestoringOnCell SurfRestOnC(Mesh); parallelFor( {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { SurfTracerRestoringDiffsCell(L, ICell) = InputField(ICell, 0); - ExactSurfRest(L, ICell, 0) = SurfRestOnC.TracerRestoreMask(L) * + ExactSurfRest(L, ICell, 0) = CaseMasks[Case][L] * SurfRestOnC.PistonVelocity * InputField(ICell, 0); }); parallelFor( - {NTracers, Mesh->NCellsOwned, NVertLayers}, - KOKKOS_LAMBDA(int L, int ICell, int KLayer) { - SurfRestOnC(NumSurfRest, L, ICell, KLayer, - SurfTracerRestoringDiffsCell); + {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { + if (CaseMasks[Case][L] == 1) { + SurfRestOnC(NumSurfRest, L, ICell, 0, + SurfTracerRestoringDiffsCell); + } }); ErrorMeasures SurfRestErrors; diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index aa2e2e89bbe2..27d4cdb75c67 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -231,20 +231,21 @@ int initTimeStepperTest(const std::string &mesh) { } // Disable all other tendencies - TestTendencies->ThicknessFluxDiv.Enabled = false; - TestTendencies->PotentialVortHAdv.Enabled = false; - TestTendencies->KEGrad.Enabled = false; - TestTendencies->SSHGrad.Enabled = false; - TestTendencies->VelocityDiffusion.Enabled = false; - TestTendencies->VelocityHyperDiff.Enabled = false; - TestTendencies->TracerHorzAdv.Enabled = false; - TestTendencies->TracerDiffusion.Enabled = false; - TestTendencies->TracerHyperDiff.Enabled = false; - TestTendencies->WindForcing.Enabled = false; - TestTendencies->BottomDrag.Enabled = false; - DefVAdv->ThickVertAdvEnabled = false; - DefVAdv->VelVertAdvEnabled = false; - DefVAdv->TracerVertAdvEnabled = false; + TestTendencies->ThicknessFluxDiv.Enabled = false; + TestTendencies->PotentialVortHAdv.Enabled = false; + TestTendencies->KEGrad.Enabled = false; + TestTendencies->SSHGrad.Enabled = false; + TestTendencies->VelocityDiffusion.Enabled = false; + TestTendencies->VelocityHyperDiff.Enabled = false; + TestTendencies->TracerHorzAdv.Enabled = false; + TestTendencies->TracerDiffusion.Enabled = false; + TestTendencies->TracerHyperDiff.Enabled = false; + TestTendencies->WindForcing.Enabled = false; + TestTendencies->SurfaceTracerRestoring.Enabled = false; + TestTendencies->BottomDrag.Enabled = false; + DefVAdv->ThickVertAdvEnabled = false; + DefVAdv->VelVertAdvEnabled = false; + DefVAdv->TracerVertAdvEnabled = false; return Err; } From acba12c77e6285f183adfab5a8aec3196dfbc723 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 30 Mar 2026 14:14:53 -0700 Subject: [PATCH 09/15] adds suggestions from pr review --- components/omega/src/ocn/AuxiliaryState.cpp | 36 ++------ components/omega/src/ocn/TendencyTerms.h | 15 +++- .../auxiliaryVars/SurfTracerRestAuxVars.cpp | 28 +----- .../ocn/auxiliaryVars/SurfTracerRestAuxVars.h | 26 +----- .../omega/test/ocn/AuxiliaryVarsTest.cpp | 90 ------------------- .../omega/test/ocn/TendencyTermsTest.cpp | 89 +++++++++++++----- 6 files changed, 89 insertions(+), 195 deletions(-) diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index a306d12151cf..0db44aa79239 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -29,8 +29,8 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux(stripDefault(Name), Mesh, VCoord), VelocityDel2Aux(stripDefault(Name), Mesh, VCoord), WindForcingAux(stripDefault(Name), Mesh), - TracerAux(stripDefault(Name), Mesh, VCoord, NTracers), - SurfTracerRestAux(stripDefault(Name), Mesh, VCoord, NTracers) { + SurfTracerRestAux(stripDefault(Name), Mesh, NTracers), + TracerAux(stripDefault(Name), Mesh, VCoord, NTracers) { GroupName = "AuxiliaryState"; if (Name != "Default") { @@ -45,8 +45,8 @@ AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, VorticityAux.registerFields(GroupName, AuxMeshName); VelocityDel2Aux.registerFields(GroupName, AuxMeshName); WindForcingAux.registerFields(GroupName, AuxMeshName); - TracerAux.registerFields(GroupName, AuxMeshName); SurfTracerRestAux.registerFields(GroupName, AuxMeshName); + TracerAux.registerFields(GroupName, AuxMeshName); } // Destructor. Unregisters the fields with IOStreams and destroys this auxiliary @@ -57,8 +57,8 @@ AuxiliaryState::~AuxiliaryState() { VorticityAux.unregisterFields(); VelocityDel2Aux.unregisterFields(); WindForcingAux.unregisterFields(); - TracerAux.unregisterFields(); SurfTracerRestAux.unregisterFields(); + TracerAux.unregisterFields(); FieldGroup::destroy(GroupName); } @@ -74,7 +74,6 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, OMEGA_SCOPE(LocVorticityAux, VorticityAux); OMEGA_SCOPE(LocVelocityDel2Aux, VelocityDel2Aux); OMEGA_SCOPE(LocWindForcingAux, WindForcingAux); - OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerVertexBot, VCoord->MinLayerVertexBot); @@ -231,7 +230,6 @@ void AuxiliaryState::computeAll(const OceanState *State, OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); OMEGA_SCOPE(LocTracerAux, TracerAux); - OMEGA_SCOPE(LocSurfTracerRestAux, SurfTracerRestAux); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); @@ -280,14 +278,6 @@ void AuxiliaryState::computeAll(const OceanState *State, }); Pacer::stop("AuxState:cellAuxState4", 2); - Pacer::start("AuxState:cellAuxState5", 2); - parallelFor( - "cellAuxState5", {NTracers, Mesh->NCellsAll}, - KOKKOS_LAMBDA(int LTracer, int ICell) { - LocSurfTracerRestAux.computeVarsOnCells(LTracer, ICell, TracerArray); - }); - Pacer::stop("AuxState:cellAuxState5", 2); - Pacer::stop("AuxState:computeAll", 1); } @@ -407,17 +397,6 @@ void AuxiliaryState::readConfigOptions(Config *OmegaConfig) { } else { ABORT_ERROR("AuxiliaryState: Unknown InterpType requested"); } - - Config SurfRestConfig("SurfaceRestoring"); - Err += OmegaConfig->get(SurfRestConfig); - - Err += SurfRestConfig.get("MaxDiff", this->SurfTracerRestAux.MaxDiff); - CHECK_ERROR_ABORT( - Err, "AuxiliaryState: MaxDiff not found in SurfaceRestoringConfig"); - if (this->SurfTracerRestAux.MaxDiff <= 0) { - ABORT_ERROR("AuxiliaryState: MaxDiff must be positive, got {}", - this->SurfTracerRestAux.MaxDiff); - } } //------------------------------------------------------------------------------ @@ -431,12 +410,13 @@ I4 AuxiliaryState::exchangeHalo() { Err += MeshHalo->exchangeFullArrayHalo(WindForcingAux.MeridStressCell, OnCell); + // Performing halo exchange on individual tracers because full halo exchange + // on a 2D array assumes the first dimension is the vertical const I4 NTracers = SurfTracerRestAux.TracersMonthlySurfClimoCell.extent_int(0); for (I4 LTracer = 0; LTracer < NTracers; ++LTracer) { - auto TracerSurfClimoCell = - Kokkos::subview(SurfTracerRestAux.TracersMonthlySurfClimoCell, - LTracer, Kokkos::ALL()); + auto TracerSurfClimoCell = Kokkos::subview( + SurfTracerRestAux.TracersMonthlySurfClimoCell, LTracer, Kokkos::ALL); Err += MeshHalo->exchangeFullArrayHalo(TracerSurfClimoCell, OnCell); } diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 5a06f172eb76..19387931efbf 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -574,8 +574,11 @@ class TracerHyperDiffOnCell { class SurfaceTracerRestoringOnCell { public: bool Enabled; - Real PistonVelocity = 1.585e-5; - Array1DI4 TracerRestoreMask; + Real PistonVelocity = 1.585e-5; ///< piston velocity + Real MaxDiff = 100.0; ///< maximum allowed restoring difference + I4 NTracersToRestore = 0; ///< number of tracers to restore + Array1DI4 TracerIdsToRestore; ///< tracer IDs to restore + /// Need to add under sea ice restoring option when that is available /// constructor declaration SurfaceTracerRestoringOnCell(const HorzMesh *Mesh); @@ -584,9 +587,13 @@ class SurfaceTracerRestoringOnCell { /// restoring values, outputs tendency array KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KMin, - const Array2DReal &SurfTracerRestoringDiffsCell) const { + const Array2DReal &TracersMonthlySurfClimoCell, + const Array3DReal &TracerCell) const { + Tend(L, ICell, KMin) += - PistonVelocity * SurfTracerRestoringDiffsCell(L, ICell); + PistonVelocity * Kokkos::clamp(TracersMonthlySurfClimoCell(L, ICell) - + TracerCell(L, ICell, KMin), + -MaxDiff, MaxDiff); } }; diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp index 5d61813f38b7..9a0219ff8204 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.cpp @@ -8,15 +8,10 @@ namespace OMEGA { SurfTracerRestAuxVars::SurfTracerRestAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh, - const VertCoord *VCoord, const I4 NTracers) - : SurfTracerRestoringDiffsCell("SurfTracerRestoringDiffsCell" + - AuxStateSuffix, - NTracers, Mesh->NCellsSize), - TracersMonthlySurfClimoCell("TracersMonthlySurfClimoCell" + + : TracersMonthlySurfClimoCell("TracersMonthlySurfClimoCell" + AuxStateSuffix, - NTracers, Mesh->NCellsSize), - MinLayerCell(VCoord->MinLayerCell) {} + NTracers, Mesh->NCellsSize) {} void SurfTracerRestAuxVars::registerFields( const std::string &AuxGroupName, // name of Auxiliary field group @@ -34,20 +29,6 @@ void SurfTracerRestAuxVars::registerFields( DimSuffix = MeshName; } - // Tracers surface restoring values - DimNames[0] = "NTracers"; - DimNames[1] = "NCells" + DimSuffix; - auto SurfTracerRestoringDiffsCellField = Field::create( - SurfTracerRestoringDiffsCell.label(), // field name - "surface tracer restoring differences", // long name/describe - "tracer units", // units - "", // CF standard Name - std::numeric_limits::min(), // min valid value - std::numeric_limits::max(), // max valid value - FillValue, // scalar for undefined entries - NDims, // number of dimensions - DimNames); // dim names - auto TracersMonthlySurfClimoCellField = Field::create(TracersMonthlySurfClimoCell.label(), "monthly surface tracer climatology", // long name/describe @@ -60,20 +41,15 @@ void SurfTracerRestAuxVars::registerFields( DimNames); // dim names // Add fields to FieldGroup - FieldGroup::addFieldToGroup(SurfTracerRestoringDiffsCell.label(), - AuxGroupName); FieldGroup::addFieldToGroup(TracersMonthlySurfClimoCell.label(), AuxGroupName); // Attach data - SurfTracerRestoringDiffsCellField->attachData( - SurfTracerRestoringDiffsCell); TracersMonthlySurfClimoCellField->attachData( TracersMonthlySurfClimoCell); } void SurfTracerRestAuxVars::unregisterFields() const { - Field::destroy(SurfTracerRestoringDiffsCell.label()); Field::destroy(TracersMonthlySurfClimoCell.label()); } diff --git a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h index 6b77c2237699..c5bfc32bb43d 100644 --- a/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h +++ b/components/omega/src/ocn/auxiliaryVars/SurfTracerRestAuxVars.h @@ -13,38 +13,14 @@ namespace OMEGA { class SurfTracerRestAuxVars { public: - Array2DReal SurfTracerRestoringDiffsCell; Array2DReal TracersMonthlySurfClimoCell; - Real MaxDiff = 100.0; // maximum allowed difference for restoring - /// Need to add under sea ice restoring option when that is available SurfTracerRestAuxVars(const std::string &AuxStateSuffix, - const HorzMesh *Mesh, const VertCoord *VCoord, - const I4 NTracers); - - KOKKOS_FUNCTION void - computeVarsOnCells(int L, int ICell, const Array3DReal &TracerCell) const { - - const I4 K = MinLayerCell(ICell); - - if ((TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K)) > - MaxDiff) { - SurfTracerRestoringDiffsCell(L, ICell) = MaxDiff; - } else if ((TracersMonthlySurfClimoCell(L, ICell) - - TracerCell(L, ICell, K)) < -MaxDiff) { - SurfTracerRestoringDiffsCell(L, ICell) = -MaxDiff; - } else { - SurfTracerRestoringDiffsCell(L, ICell) = - TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, K); - } - } + const HorzMesh *Mesh, const I4 NTracers); void registerFields(const std::string &AuxGroupName, const std::string &MeshName) const; void unregisterFields() const; - - private: - Array1DI4 MinLayerCell; }; } // namespace OMEGA diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp index f62856065bf3..44957762f953 100644 --- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp +++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp @@ -15,7 +15,6 @@ #include "VertCoord.h" #include "auxiliaryVars/KineticAuxVars.h" #include "auxiliaryVars/LayerThicknessAuxVars.h" -#include "auxiliaryVars/SurfTracerRestAuxVars.h" #include "auxiliaryVars/TracerAuxVars.h" #include "auxiliaryVars/VelocityDel2AuxVars.h" #include "auxiliaryVars/VorticityAuxVars.h" @@ -69,8 +68,6 @@ struct TestSetupPlane { ErrorMeasures ExpectedNormalStressErrors = {0.0033910709836867704, 0.0039954090464502795}; - ErrorMeasures ExpectedSurfTracerRestErrors = {0, 0}; - KOKKOS_FUNCTION Real layerThickness(Real X, Real Y) const { return 2 + std::cos(TwoPi * X / Lx) * std::cos(TwoPi * Y / Ly); } @@ -199,8 +196,6 @@ struct TestSetupSphere { ErrorMeasures ExpectedNormalStressErrors = {0.0038588958862868362, 0.003813760171030077}; - ErrorMeasures ExpectedSurfTracerRestErrors = {0, 0}; - KOKKOS_FUNCTION Real layerThickness(Real Lon, Real Lat) const { return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4)); } @@ -782,89 +777,6 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell, return Err; } - -int testSurfTracerRestAuxVars(int NVertLayers, int NTracers, Real RTol) { - int Err = 0; - TestSetup Setup; - - const auto Mesh = HorzMesh::getDefault(); - const auto VCoord = VertCoord::getDefault(); - - SurfTracerRestAuxVars SurfTracerRestAux("", Mesh, VCoord, NTracers); - SurfTracerRestAux.MaxDiff = 0.5; // choose a small max‑difference so all - // three branches (>MaxDiff, <-MaxDiff, - // inside) are exercised. - - Array2DReal InputField("InputField", Mesh->NCellsSize, 1); - deepCopy(InputField, 0.0); - - // fill monthly climo with a simple analytic function that produces - // positive/negative variations, e.g. tracer+velocityx: - Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { - return Setup.tracer(X, Y) + Setup.velocityX(X, Y); - }, - InputField, Geom, Mesh, OnCell); - - parallelFor( - {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { - SurfTracerRestAux.TracersMonthlySurfClimoCell(L, ICell) = - InputField(ICell, 0); - }); - - Array3DReal TracersOnCell("TracersOnCell", NTracers, Mesh->NCellsSize, - NVertLayers); - deepCopy(TracersOnCell, 0.0); - // fill the tracer array and add a trivial vertical dependence so that - // the surface value differs from deeper layers - Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.tracer(X, Y); }, - TracersOnCell, Geom, Mesh, OnCell); - parallelFor( - {NTracers, Mesh->NCellsSize, NVertLayers}, - KOKKOS_LAMBDA(int L, int ICell, int K) { - TracersOnCell(L, ICell, K) += 0.04_Real * K; - }); - - // compute the numerical result - parallelFor( - {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { - SurfTracerRestAux.computeVarsOnCells(L, ICell, TracersOnCell); - }); - - // build the exact result by re‑implementing the class logic, including - // lookup of MinLayerCell - Array1DReal ExactSurfRest("ExactSurfRest", Mesh->NCellsOwned); - Array1DReal NumSurfRest("NumSurfRest", Mesh->NCellsOwned); - const auto &MinLayerCell = VCoord->MinLayerCell; - parallelFor( - {Mesh->NCellsOwned}, KOKKOS_LAMBDA(int ICell) { - const int K = MinLayerCell(ICell); - const Real diff = - SurfTracerRestAux.TracersMonthlySurfClimoCell(0, ICell) - - TracersOnCell(0, ICell, K); - if (diff > SurfTracerRestAux.MaxDiff) { - ExactSurfRest(ICell) = SurfTracerRestAux.MaxDiff; - } else if (diff < -SurfTracerRestAux.MaxDiff) { - ExactSurfRest(ICell) = -SurfTracerRestAux.MaxDiff; - } else { - ExactSurfRest(ICell) = diff; - } - NumSurfRest(ICell) = - SurfTracerRestAux.SurfTracerRestoringDiffsCell(0, ICell); - }); - - ErrorMeasures SurfRestErrors; - Err += - computeErrors(SurfRestErrors, NumSurfRest, ExactSurfRest, Mesh, OnCell); - Err += checkErrors("AuxVarsTest", "SurfTracerRest", SurfRestErrors, - Setup.ExpectedSurfTracerRestErrors, RTol); - - if (Err == 0) { - LOG_INFO("AuxVarsTest: SurfTracerRestAuxVars PASS"); - } - return Err; -} //------------------------------------------------------------------------------ // The initialization routine for aux vars testing int initAuxVarsTest(const std::string &mesh) { @@ -935,8 +847,6 @@ int auxVarsTest(const std::string &mesh = DefaultMeshFile) { Err += testWindForcingAuxVars(RTol); - Err += testSurfTracerRestAuxVars(NVertLayers, NTracers, RTol); - if (Err == 0) { LOG_INFO("AuxVarsTest: Successful completion"); } diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 01bc834630d6..5d117fd171a0 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -36,6 +36,7 @@ #include #include +#include using namespace OMEGA; @@ -1018,8 +1019,7 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { I4 Err = 0; TestSetup Setup; - const auto Mesh = HorzMesh::getDefault(); - const auto VCoord = VertCoord::getDefault(); + const auto Mesh = HorzMesh::getDefault(); if (NTracers < 3) { LOG_ERROR("TendencyTermsTest: SurfaceTracerRestoring requires at least 3 " @@ -1028,52 +1028,97 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { return 1; } + // Test multiple cases with different combinations of tracers being restored const char *CaseLabels[3] = {"SurfaceTracerRestoringSalinityOnly", "SurfaceTracerRestoringTemperatureOnly", "SurfaceTracerRestoringTempSaltDebug"}; - const I4 CaseMasks[3][3] = { - {0, 1, 0}, - {1, 0, 0}, - {1, 1, 1}, + const std::vector> CaseTracerIds = { + {1}, // Salinity Only + {0}, // Temperature Only + {0, 1, 2}, // Temperature, Salinity, and a Debug Tracer }; + // Loop over cases, computing exact result, numerical result, and + // errors for each. for (int Case = 0; Case < 3; ++Case) { Array3DReal ExactSurfRest("ExactSurfRest", NTracers, Mesh->NCellsOwned, NVertLayers); Array2DReal InputField("InputField", Mesh->NCellsSize, 1); - Array2DReal SurfTracerRestoringDiffsCell("SurfTracerRestoringDiffsCell", - NTracers, Mesh->NCellsSize, - NVertLayers); + Array2DReal TracersMonthlySurfClimoCell("TracersMonthlySurfClimoCell", + NTracers, Mesh->NCellsSize); + Array3DReal TracersOnCell("TracersOnCell", NTracers, Mesh->NCellsSize, + NVertLayers); Array3DReal NumSurfRest("NumSurfRest", NTracers, Mesh->NCellsOwned, NVertLayers); deepCopy(ExactSurfRest, 0); deepCopy(InputField, 0); - deepCopy(SurfTracerRestoringDiffsCell, 0); + deepCopy(TracersMonthlySurfClimoCell, 0); + deepCopy(TracersOnCell, 0); deepCopy(NumSurfRest, 0); + // Set Input Field values. Use a combination of scalarB and vectorX to + // ensure the full surface tracer restoring logic is exercised (i.e., + // cases where the difference is >, <, and within the max-difference + // threshold). Err += setScalar( - KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarA(X, Y); }, + KOKKOS_LAMBDA(Real X, Real Y) { + return Setup.scalarB(X, Y) + Setup.vectorX(X, Y); + }, InputField, Geom, Mesh, OnCell); - SurfaceTracerRestoringOnCell SurfRestOnC(Mesh); + // Set TracersOnCell values (use scalarB for simplicity, but could be + // any field). + Err += setScalar( + KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y); }, + TracersOnCell, Geom, Mesh, OnCell); + parallelFor( + {NTracers, Mesh->NCellsSize, NVertLayers}, + KOKKOS_LAMBDA(int L, int ICell, int K) { + TracersOnCell(L, ICell, K) += 0.04_Real * K; + }); + SurfaceTracerRestoringOnCell SurfRestOnC(Mesh); + SurfRestOnC.MaxDiff = 0.5; // choose a small max‑difference so all + // three branches (>MaxDiff, <-MaxDiff, + // inside) are exercised. + SurfRestOnC.PistonVelocity = 1.585e-5; + + // Build host-selected tracer IDs for restoring and copy to device. + const I4 NTracersToRestore = static_cast(CaseTracerIds[Case].size()); + Array1DI4 TracerIdsToRestore("TracerIdsToRestore", NTracersToRestore); + HostArray1DI4 TracerIdsToRestoreH("TracerIdsToRestoreH", + NTracersToRestore); + for (I4 R = 0; R < NTracersToRestore; ++R) { + TracerIdsToRestoreH(R) = CaseTracerIds[Case][R]; + } + deepCopy(TracerIdsToRestore, TracerIdsToRestoreH); + + // Compute exact result using the same logic as + // SurfaceTracerRestoringOnCell, but iterating + // selected tracer IDs to match restoring implementation. parallelFor( - {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { - SurfTracerRestoringDiffsCell(L, ICell) = InputField(ICell, 0); - ExactSurfRest(L, ICell, 0) = CaseMasks[Case][L] * - SurfRestOnC.PistonVelocity * - InputField(ICell, 0); + {NTracersToRestore, Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int R, int ICell) { + const int L = TracerIdsToRestore(R); + TracersMonthlySurfClimoCell(L, ICell) = InputField(ICell, 0); + const Real Diff = TracersMonthlySurfClimoCell(L, ICell) - + TracersOnCell(L, ICell, 0); + + ExactSurfRest(L, ICell, 0) = + SurfRestOnC.PistonVelocity * + Kokkos::clamp(Diff, -SurfRestOnC.MaxDiff, SurfRestOnC.MaxDiff); }); + // Compute numerical result parallelFor( - {NTracers, Mesh->NCellsOwned}, KOKKOS_LAMBDA(int L, int ICell) { - if (CaseMasks[Case][L] == 1) { - SurfRestOnC(NumSurfRest, L, ICell, 0, - SurfTracerRestoringDiffsCell); - } + {NTracersToRestore, Mesh->NCellsOwned}, + KOKKOS_LAMBDA(int R, int ICell) { + const int L = TracerIdsToRestore(R); + SurfRestOnC(NumSurfRest, L, ICell, 0, TracersMonthlySurfClimoCell, + TracersOnCell); }); ErrorMeasures SurfRestErrors; From 56be44b7f6fdb082c02fd202708e6b39779ff289 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 30 Mar 2026 15:02:56 -0700 Subject: [PATCH 10/15] rebases and updates docs --- components/omega/doc/devGuide/AuxiliaryVariables.md | 3 +-- components/omega/doc/devGuide/Forcing.md | 13 +++++++------ components/omega/doc/devGuide/TendencyTerms.md | 1 - components/omega/doc/userGuide/Forcing.md | 6 +++--- 4 files changed, 11 insertions(+), 12 deletions(-) diff --git a/components/omega/doc/devGuide/AuxiliaryVariables.md b/components/omega/doc/devGuide/AuxiliaryVariables.md index de96e3186921..bd43757a42e3 100644 --- a/components/omega/doc/devGuide/AuxiliaryVariables.md +++ b/components/omega/doc/devGuide/AuxiliaryVariables.md @@ -104,8 +104,7 @@ The following auxiliary variable groups are currently implemented: | WindForcingAuxVars | ZonalStressCell || || MeridStressCell || || NormalStressEdge || -| SurfTracerRestAuxVars | SurfTracerRestoringDiffsCell || -|| TracersMonthlySurfClimoCell || +| SurfTracerRestAuxVars | TracersMonthlySurfClimoCell || ## See Also diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 649fff6fff65..00534b0a448d 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -44,7 +44,7 @@ pathways in Omega, currently this includes: 2. Auxiliary-state compute forms bounded restoring differences: - `SurfTracerRestoringDiffsCell = clamp(target - tracer_surface, +/- MaxDiff)` 3. Tendency term applies restoring only at surface layer and only for tracers - enabled by restore mask. + selected from `SurfaceRestoring.TracersToRestore`. ### Surface tracer restoring key classes/components @@ -55,10 +55,11 @@ pathways in Omega, currently this includes: - `SurfaceTracerRestoringOnCell` tendency term - Applies `PistonVelocity * SurfTracerRestoringDiffsCell` at surface - `Tendencies` - - Builds `TracerRestoreMask` from `SurfaceRestoring.TracersToRestore` + - Parses `SurfaceRestoring.TracersToRestore` and resolves tracer indices + - Builds `TracerIdsToRestore` and `NTracersToRestore` - Applies tracer-selection logic at call site in `computeTracerTendenciesOnly` - - Aborts if restoring is enabled but no restore mask is available + - Aborts if restoring is enabled but no tracer IDs are available ### Surface tracer restoring config coupling @@ -67,15 +68,15 @@ pathways in Omega, currently this includes: - `Omega.SurfaceRestoring.PistonVelocity` - tendency scaling - `Omega.SurfaceRestoring.TracersToRestore` - - tracer-level enable list used to build `TracerRestoreMask` + - tracer-level enable list used to build `TracerIdsToRestore` - `Omega.Tendencies.SurfaceTracerRestoringEnable` - gates restoring tendency execution - ## Notes +## Notes - If a tracer is not listed in `TracersToRestore`, no restoring tendency is applied to that tracer. -- If restoring is enabled but no restore mask is available at tendency +- If restoring is enabled but no tracer IDs are available at tendency compute-time, Omega aborts with an error. - `MaxDiff` must be positive. A runtime check will error out if not. - It is assumed that the incoming `TracersMonthlySurfClimoCell` values have the correct units appropriate for Omega. If this is not true, conversion should be implemented. diff --git a/components/omega/doc/devGuide/TendencyTerms.md b/components/omega/doc/devGuide/TendencyTerms.md index 01641a60685a..08005b1c00be 100644 --- a/components/omega/doc/devGuide/TendencyTerms.md +++ b/components/omega/doc/devGuide/TendencyTerms.md @@ -47,4 +47,3 @@ implemented: Additional information on forcing (currently wind forcing and surface tracer restoring) is detailed in [](omega-dev-forcing). ->>>>>>> 41ae670f82 (adds more documentation and renames vars) diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index a7179e392c7c..e1df87fcdf3b 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -67,8 +67,8 @@ Omega: - `MaxDiff`: cap on restoring difference magnitude - `SurfaceTracerRestoringEnable`: switch to enable surface tracer restoring -When restoring is enabled, Omega builds an internal `TracerRestoreMask` from -`TracersToRestore` and applies restoring only to tracers whose mask value is 1. +When restoring is enabled, Omega resolves `TracersToRestore` into an internal +list of tracer IDs and applies restoring only to tracers in that list. ### Restoring target fields @@ -84,6 +84,6 @@ The restoring tendency is computed at the surface layer only and is limited by - If a tracer is not listed in `TracersToRestore`, no restoring tendency is applied to that tracer. -- If surface restoring is enabled but no restore mask is available at tendency +- If surface restoring is enabled but no tracer IDs are available at tendency compute-time, Omega aborts with an error. - `MaxDiff` must be positive. A runtime check will error out if not. From c6ddfac11d764ac3f2872a59c0105dbc1401c2b6 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 30 Mar 2026 15:42:15 -0700 Subject: [PATCH 11/15] removes redundent call of OmegaConfig --- components/omega/src/ocn/Tendencies.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 0fc59d173c44..030277596901 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -142,17 +142,12 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config ) { Error Err; // error code - // we need access to the top-level Omega config when pulling values - // from other sections such as SurfaceRestoring - Config *OmegaConfig = Config::getOmegaConfig(); - Config TendConfig("Tendencies"); Err += OmegaConfig->get(TendConfig); CHECK_ERROR_ABORT(Err, "Tendencies: Tendencies group not found in Config"); Err += TendConfig.get("ThicknessFluxTendencyEnable", this->ThicknessFluxDiv.Enabled); - CHECK_ERROR_ABORT( Err, "Tendencies: ThicknessFluxTendencyEnable not found in TendConfig"); @@ -253,7 +248,6 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config CHECK_ERROR_ABORT(Err, "Tendencies: EddyDiff4 not found in TendConfig"); } - Err += TendConfig.get("PressureGradTendencyEnable", this->PGrad->Enabled); CHECK_ERROR_ABORT( Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig"); From 4f79be503f1fc07681499aedff5cac67317b0585 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 30 Mar 2026 16:02:15 -0700 Subject: [PATCH 12/15] reverts submodule changes --- cime | 2 +- components/omega/configs/Default.yml | 2 +- components/omega/doc/index.md | 3 --- components/omega/src/ocn/AuxiliaryState.cpp | 4 ++-- components/omega/src/ocn/TendencyTerms.cpp | 3 +-- externals/ekat | 2 +- externals/scorpio | 2 +- 7 files changed, 7 insertions(+), 11 deletions(-) diff --git a/cime b/cime index c64260ed94af..bd710a378bbb 160000 --- a/cime +++ b/cime @@ -1 +1 @@ -Subproject commit c64260ed94aff167a59635ed4a9ae5af16824e91 +Subproject commit bd710a378bbb5ad9efef5d958b63dab221e7b1fe diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index d1222da0dc39..40a5fa8ded40 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -50,7 +50,7 @@ Omega: ViscDel4: 1.2e11 DivFactor: 1.0 WindForcingTendencyEnable: false - SurfaceTracerRestoringEnable: true + SurfaceTracerRestoringEnable: false BottomDragTendencyEnable: false BottomDragCoeff: 0.0 TracerHorzAdvTendencyEnable: true diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index aac646b407ad..6712bb3e7f8d 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -98,11 +98,8 @@ devGuide/VertCoord devGuide/PGrad devGuide/Timing devGuide/VerticalMixingCoeff -<<<<<<< HEAD devGuide/VertAdv -======= devGuide/Forcing ->>>>>>> 41ae670f82 (adds more documentation and renames vars) ``` ```{toctree} diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp index 0db44aa79239..30efe9b5ffc5 100644 --- a/components/omega/src/ocn/AuxiliaryState.cpp +++ b/components/omega/src/ocn/AuxiliaryState.cpp @@ -3,8 +3,8 @@ #include "Field.h" #include "Logging.h" #include "Pacer.h" -#include "TimeStepper.h" #include "Tendencies.h" +#include "TimeStepper.h" namespace OMEGA { @@ -74,6 +74,7 @@ void AuxiliaryState::computeMomAux(const OceanState *State, int ThickTimeLevel, OMEGA_SCOPE(LocVorticityAux, VorticityAux); OMEGA_SCOPE(LocVelocityDel2Aux, VelocityDel2Aux); OMEGA_SCOPE(LocWindForcingAux, WindForcingAux); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerVertexBot, VCoord->MinLayerVertexBot); @@ -230,7 +231,6 @@ void AuxiliaryState::computeAll(const OceanState *State, OMEGA_SCOPE(LocLayerThicknessAux, LayerThicknessAux); OMEGA_SCOPE(LocTracerAux, TracerAux); - OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp index 35efb216a34d..27b19e9f7a68 100644 --- a/components/omega/src/ocn/TendencyTerms.cpp +++ b/components/omega/src/ocn/TendencyTerms.cpp @@ -113,7 +113,7 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh, SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell( const HorzMesh *Mesh) {} - + void TracerHorzAdvOnCell::init() { const HorzMesh *Mesh = this->HorzontalMesh; const auto MaxEdges2 = Mesh->MaxEdges2; @@ -136,7 +136,6 @@ void TracerHorzAdvOnCell::init() { {NEdgesAll}, KOKKOS_LAMBDA(int IEdge) { masksAndCoefficients(IEdge); }); Kokkos::fence(); } - } // end namespace OMEGA //===----------------------------------------------------------------------===// diff --git a/externals/ekat b/externals/ekat index 7156c944701f..d13cf66f3543 160000 --- a/externals/ekat +++ b/externals/ekat @@ -1 +1 @@ -Subproject commit 7156c944701fd19fa5e27ffcb86819f12851ae67 +Subproject commit d13cf66f3543df10445cd1fd35edac3165cc2833 diff --git a/externals/scorpio b/externals/scorpio index 9df851035a89..a486b499e811 160000 --- a/externals/scorpio +++ b/externals/scorpio @@ -1 +1 @@ -Subproject commit 9df851035a89ba1ca0581ccdbb88fcfc0ec38992 +Subproject commit a486b499e8113f92c136f020cc0e7a68c664f01a From 8fa3a02603444660f36ec16651d0f8534d2b3341 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 6 Apr 2026 10:15:46 -0700 Subject: [PATCH 13/15] removes maxdiff and adds suggestions from pr review --- components/omega/configs/Default.yml | 1 - components/omega/doc/devGuide/Forcing.md | 23 ++++++++----------- components/omega/doc/userGuide/Forcing.md | 5 +--- components/omega/src/ocn/OceanInit.cpp | 9 ++++++++ components/omega/src/ocn/TendencyTerms.h | 6 ++--- .../omega/test/ocn/TendencyTermsTest.cpp | 11 ++------- 6 files changed, 23 insertions(+), 32 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 40a5fa8ded40..549a6fd3bec3 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -34,7 +34,6 @@ Omega: SurfaceRestoring: TracersToRestore: [Temperature, Salinity] PistonVelocity: 1.585e-5 - MaxDiff: 100.0 VertCoord: MovementWeightType: Uniform PressureGrad: diff --git a/components/omega/doc/devGuide/Forcing.md b/components/omega/doc/devGuide/Forcing.md index 00534b0a448d..5d2d52f20d6b 100644 --- a/components/omega/doc/devGuide/Forcing.md +++ b/components/omega/doc/devGuide/Forcing.md @@ -40,11 +40,9 @@ pathways in Omega, currently this includes: ### Surface tracer restoring data flow -1. External fields provide target values: `TracersMonthlySurfClimoCell` -2. Auxiliary-state compute forms bounded restoring differences: - - `SurfTracerRestoringDiffsCell = clamp(target - tracer_surface, +/- MaxDiff)` -3. Tendency term applies restoring only at surface layer and only for tracers - selected from `SurfaceRestoring.TracersToRestore`. +1. External fields provide target values: `TracersMonthlySurfClimoCell` (values and units should match the state variables) +2. Auxiliary-state compute forms restoring differences: `SurfTracerRestoringDiffsCell = target - tracer_surface` +3. Tendency term applies restoring only at surface layer and only for tracers selected from `SurfaceRestoring.TracersToRestore`. ### Surface tracer restoring key classes/components @@ -63,8 +61,6 @@ pathways in Omega, currently this includes: ### Surface tracer restoring config coupling -- `Omega.SurfaceRestoring.MaxDiff` - - bound for restoring differences - `Omega.SurfaceRestoring.PistonVelocity` - tendency scaling - `Omega.SurfaceRestoring.TracersToRestore` @@ -74,11 +70,10 @@ pathways in Omega, currently this includes: ## Notes -- If a tracer is not listed in `TracersToRestore`, no restoring tendency is - applied to that tracer. -- If restoring is enabled but no tracer IDs are available at tendency - compute-time, Omega aborts with an error. -- `MaxDiff` must be positive. A runtime check will error out if not. -- It is assumed that the incoming `TracersMonthlySurfClimoCell` values have the correct units appropriate for Omega. If this is not true, conversion should be implemented. +- If a tracer is not listed in `TracersToRestore`, no restoring tendency is applied to that tracer. +- If restoring is enabled but no tracer IDs are available at tendency compute-time, Omega aborts with an error. +- It is assumed that the incoming `TracersMonthlySurfClimoCell` fields (values and units) match the Omega state variables (i.e. conservative temperature and absolute salinity for Teos-10). If not, a pre-processing conversion should be implemented. - Surface tracer restoring is active everywhere if enabled. A flag to turn it off under sea ice will need to be added in later development if this feature is desired. -- A global sum has not been implemented for the surface tracer restoring, but should be added in later development. +- Unlike MPAS-Ocean, a `MaxDiff` clamping is not applied here. This check should instead be implemented in Ocean Validate when that is available. +- A global scaling to ensure zero-sum has not been implemented for the surface tracer restoring, but should be added in later development. +- At this stage, there is no temporal interpolation applied to the restoring targets, the raw `TracersMonthlySurfClimoCell` snapshot is used. diff --git a/components/omega/doc/userGuide/Forcing.md b/components/omega/doc/userGuide/Forcing.md index e1df87fcdf3b..469693958b41 100644 --- a/components/omega/doc/userGuide/Forcing.md +++ b/components/omega/doc/userGuide/Forcing.md @@ -56,7 +56,6 @@ Omega: SurfaceRestoring: TracersToRestore: [Temperature, Salinity] PistonVelocity: 1.585e-5 - MaxDiff: 100.0 Tendencies: SurfaceTracerRestoringEnable: true @@ -64,7 +63,6 @@ Omega: - `TracersToRestore`: list of tracer names that restoring is applied to - `PistonVelocity`: restoring rate coefficient -- `MaxDiff`: cap on restoring difference magnitude - `SurfaceTracerRestoringEnable`: switch to enable surface tracer restoring When restoring is enabled, Omega resolves `TracersToRestore` into an internal @@ -78,7 +76,7 @@ Surface restoring uses auxiliary fields: - `SurfTracerRestoringDiffsCell`: computed target-minus-state differences The restoring tendency is computed at the surface layer only and is limited by -`MaxDiff`. +the configured `PistonVelocity` and target-minus-state difference. ## Notes @@ -86,4 +84,3 @@ The restoring tendency is computed at the surface layer only and is limited by applied to that tracer. - If surface restoring is enabled but no tracer IDs are available at tendency compute-time, Omega aborts with an error. -- `MaxDiff` must be positive. A runtime check will error out if not. diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index a22e3b3b0667..02b9f9010d78 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -139,6 +139,15 @@ int initOmegaModules(MPI_Comm Comm) { PressureGrad::init(); Eos::init(); Tendencies::init(); + + // Validate SurfaceTracerRestoring configuration + Tendencies *DefTend = Tendencies::getDefault(); + if (DefTend->SurfaceTracerRestoring.Enabled && + DefTend->SurfaceTracerRestoring.NTracersToRestore == 0) { + ABORT_ERROR("OceanInit: SurfaceTracerRestoring is enabled but " + "TracersToRestore is empty"); + } + TimeStepper::init2(); Err = OceanState::init(); diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h index 19387931efbf..4a1d7cbe3968 100644 --- a/components/omega/src/ocn/TendencyTerms.h +++ b/components/omega/src/ocn/TendencyTerms.h @@ -575,7 +575,6 @@ class SurfaceTracerRestoringOnCell { public: bool Enabled; Real PistonVelocity = 1.585e-5; ///< piston velocity - Real MaxDiff = 100.0; ///< maximum allowed restoring difference I4 NTracersToRestore = 0; ///< number of tracers to restore Array1DI4 TracerIdsToRestore; ///< tracer IDs to restore /// Need to add under sea ice restoring option when that is available @@ -591,9 +590,8 @@ class SurfaceTracerRestoringOnCell { const Array3DReal &TracerCell) const { Tend(L, ICell, KMin) += - PistonVelocity * Kokkos::clamp(TracersMonthlySurfClimoCell(L, ICell) - - TracerCell(L, ICell, KMin), - -MaxDiff, MaxDiff); + PistonVelocity * + (TracersMonthlySurfClimoCell(L, ICell) - TracerCell(L, ICell, KMin)); } }; diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp index 5d117fd171a0..a63154561c8b 100644 --- a/components/omega/test/ocn/TendencyTermsTest.cpp +++ b/components/omega/test/ocn/TendencyTermsTest.cpp @@ -1060,9 +1060,7 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { deepCopy(NumSurfRest, 0); // Set Input Field values. Use a combination of scalarB and vectorX to - // ensure the full surface tracer restoring logic is exercised (i.e., - // cases where the difference is >, <, and within the max-difference - // threshold). + // ensure the full surface tracer restoring logic is exercised. Err += setScalar( KOKKOS_LAMBDA(Real X, Real Y) { return Setup.scalarB(X, Y) + Setup.vectorX(X, Y); @@ -1081,9 +1079,6 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { }); SurfaceTracerRestoringOnCell SurfRestOnC(Mesh); - SurfRestOnC.MaxDiff = 0.5; // choose a small max‑difference so all - // three branches (>MaxDiff, <-MaxDiff, - // inside) are exercised. SurfRestOnC.PistonVelocity = 1.585e-5; // Build host-selected tracer IDs for restoring and copy to device. @@ -1107,9 +1102,7 @@ int testSurfaceTracerRestoringOnCell(int NVertLayers, int NTracers, Real RTol) { const Real Diff = TracersMonthlySurfClimoCell(L, ICell) - TracersOnCell(L, ICell, 0); - ExactSurfRest(L, ICell, 0) = - SurfRestOnC.PistonVelocity * - Kokkos::clamp(Diff, -SurfRestOnC.MaxDiff, SurfRestOnC.MaxDiff); + ExactSurfRest(L, ICell, 0) = SurfRestOnC.PistonVelocity * Diff; }); // Compute numerical result From a9e2fa446132aff09884039f0c71d90385947804 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 6 Apr 2026 12:11:19 -0700 Subject: [PATCH 14/15] linting --- 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 030277596901..76554da242ef 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -837,7 +837,7 @@ void Tendencies::computeTracerTendenciesOnly( }); Pacer::stop("Tend:surfaceTracerRestoring", 2); } - + Pacer::stop("Tend:computeTracerTendenciesOnly", 1); } // end tracer tendency compute From 761b24dadb773c7af5bfb5dc1f77f29dbdacd800 Mon Sep 17 00:00:00 2001 From: Kat Smith Date: Mon, 6 Apr 2026 14:29:31 -0700 Subject: [PATCH 15/15] adds suggestion from pr review --- components/omega/src/ocn/Tendencies.cpp | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index 76554da242ef..7788bda4363d 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -302,14 +302,9 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config static_cast(TracerIdsToRestoreVec.size()); this->SurfaceTracerRestoring.TracerIdsToRestore = Array1DI4( "TracerIdsToRestore", this->SurfaceTracerRestoring.NTracersToRestore); - HostArray1DI4 TracerIdsToRestoreH( - "TracerIdsToRestoreH", - this->SurfaceTracerRestoring.NTracersToRestore); - for (I4 R = 0; R < this->SurfaceTracerRestoring.NTracersToRestore; ++R) { - TracerIdsToRestoreH(R) = TracerIdsToRestoreVec[R]; - } deepCopy(this->SurfaceTracerRestoring.TracerIdsToRestore, - TracerIdsToRestoreH); + HostArray1DI4(TracerIdsToRestoreVec.data(), + TracerIdsToRestoreVec.size())); } }