From 53afa84178eda28f7787a09a2d598d758b525483 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 12 Mar 2026 12:37:07 -0700 Subject: [PATCH 01/12] Add prescribeState with one option to set the state at timeLevel2 to the state at timeLevel1 --- .../omega/src/timeStepping/TimeStepper.cpp | 58 +++++++++++++++++++ .../omega/src/timeStepping/TimeStepper.h | 25 ++++++++ 2 files changed, 83 insertions(+) diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 7ee026c384b7..b1e857bd1e24 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -460,6 +460,64 @@ void TimeStepper::updateStateByTend(OceanState *State1, int TimeLevel1, updateVelocityByTend(State1, TimeLevel1, State2, TimeLevel2, Coeff); } +//------------------------------------------------------------------------------ +// Reset state variables to their initial values +void TimeStepper::prescribeThickness(OceanState *State1, int TimeLevel1, + OceanState *State2, int TimeLevel2) const { + + Array2DReal LayerThick1 = State1->getLayerThickness(TimeLevel1); + Array2DReal LayerThick2 = State2->getLayerThickness(TimeLevel2); + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + + parallelForOuter( + "prescribeThickness", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + LayerThick1(ICell, K) = LayerThick2(ICell, K); + }); + }); +} + +//------------------------------------------------------------------------------ +void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, + OceanState *State2, int TimeLevel2) const { + + Array2DReal NormalVel1 = State1->getNormalVelocity(TimeLevel1); + Array2DReal NormalVel2 = State2->getNormalVelocity(TimeLevel2); + + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + + parallelForOuter( + "prescribeVelocity", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRange(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + NormalVel1(IEdge, K) = NormalVel2(IEdge, K); + }); + }); +} + +//------------------------------------------------------------------------------ +void TimeStepper::prescribeState(OceanState *State1, int TimeLevel1, + OceanState *State2, int TimeLevel2) const { + prescribeThickness(State1, TimeLevel1, State2, TimeLevel2); + prescribeVelocity(State1, TimeLevel1, State2, TimeLevel2); +} + //------------------------------------------------------------------------------ // Updates tracers // NextTracers = (CurTracers * LayerThickness2(TimeLevel2)) + diff --git a/components/omega/src/timeStepping/TimeStepper.h b/components/omega/src/timeStepping/TimeStepper.h index 6e4d581e2e7f..3696678fa0cb 100644 --- a/components/omega/src/timeStepping/TimeStepper.h +++ b/components/omega/src/timeStepping/TimeStepper.h @@ -204,6 +204,31 @@ class TimeStepper { TimeInterval Coeff ///< [in] time-related coeff for tendency ) const; + /// Resets the layer thickness at the working time level to the initial + /// condition stored at the reference time level. + void prescribeThickness( + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2 ///< [in] time level index of the destination data + ) const; + + /// Resets the velocity at the working time level to the initial condition. + void prescribeVelocity( + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2 ///< [in] time level index of the destination data + ) const; + + /// Reset thickness and velocity to their initial values + void prescribeState( + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2 ///< [in] time level index of the destination data + ) const; + /// Updates tracers /// NextTracers = (CurTracers * LayerThickness2(TimeLevel2)) + /// Coeff * TracersTend) / LayerThickness1(TimeLevel1) From c0d364a0a1c42e74a4d3976c79c50e20a3f6c1af Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 12 Mar 2026 12:46:17 -0700 Subject: [PATCH 02/12] Use prescribeState in time steppers --- components/omega/src/timeStepping/ForwardBackwardStepper.cpp | 3 +++ components/omega/src/timeStepping/RungeKutta2Stepper.cpp | 2 ++ components/omega/src/timeStepping/RungeKutta4Stepper.cpp | 3 +++ 3 files changed, 8 insertions(+) diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index 8717b46e3757..00da6d0a23e2 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -68,6 +68,9 @@ void ForwardBackwardStepper::doStep( updateTracersByTend(NextTracerArray, CurTracerArray, State, ThickNextLevel, State, ThickCurLevel, TimeStep); + prescribeThickness(State, NextLevel, State, CurLevel); + prescribeVelocity(State, NextLevel, State, CurLevel); + // Update time levels (New -> Old) of prognostic variables with halo // exchanges const MPI_Comm Comm = MeshHalo->getComm(); diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index 4bfefe1633eb..d34a3de804eb 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -56,6 +56,8 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state updateTracersByTend(NextTracerArray, CurTracerArray, State, NextLevel, State, CurLevel, TimeStep); + prescribeState(State, NextLevel, State, CurLevel); + // Update time levels (New -> Old) of prognostic variables with halo // exchanges const MPI_Comm Comm = MeshHalo->getComm(); diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index bfa8255c821d..6742a11263f7 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -88,6 +88,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, CurLevel, RKB[Stage] * TimeStep); + prescribeState(State, NextLevel, State, CurLevel); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); } else { // every other stage does: @@ -96,6 +97,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state // q^{n+1} += RKB[stage] * dt * R^{(s)} updateStateByTend(ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); + prescribeState(State, NextLevel, State, CurLevel); updateTracersByTend(ProvisTracers, CurTracerArray, ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); @@ -109,6 +111,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state CurLevel, CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, NextLevel, RKB[Stage] * TimeStep); + prescribeState(State, NextLevel, State, CurLevel); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); } } From 42db5da1a211fc7f0e66f81dab8038a525875a22 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 12 Mar 2026 13:45:49 -0700 Subject: [PATCH 03/12] Add a new config option to turn prescribeThickness, prescribeVelocity off by default --- components/omega/configs/Default.yml | 2 + .../omega/src/timeStepping/TimeStepper.cpp | 144 +++++++++++++----- .../omega/src/timeStepping/TimeStepper.h | 20 +++ 3 files changed, 129 insertions(+), 37 deletions(-) diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index 963ad7b01f16..c3ad529fa635 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -22,6 +22,8 @@ Omega: IODefaultFormat: pnetcdf State: NTimeLevels: 2 + PrescribeThicknessType: None + PrescribeVelocityType: None Advection: Coef3rdOrder: 0.25 FluxThicknessType: Center diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index b1e857bd1e24..2c7d44f1656b 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -12,7 +12,6 @@ #include "RungeKutta4Stepper.h" namespace OMEGA { - //------------------------------------------------------------------------------ // create the static class members // Default model time stepper @@ -20,6 +19,10 @@ TimeStepper *TimeStepper::DefaultTimeStepper = nullptr; // All defined time steppers std::map> TimeStepper::AllTimeSteppers; +PrescribeStateType TimeStepper::DefaultPrescribeThicknessMode = + PrescribeStateType::None; +PrescribeStateType TimeStepper::DefaultPrescribeVelocityMode = + PrescribeStateType::None; //------------------------------------------------------------------------------ // utility functions @@ -44,6 +47,33 @@ TimeStepperType getTimeStepperFromStr(const std::string &InString) { return TimeStepperChoice; } +PrescribeStateType getPrescribeThicknessTypeFromStr(const std::string &InString) { + + if (InString == "None") { + return PrescribeStateType::None; + } + if (InString == "Init") { + return PrescribeStateType::Init; + } + + ABORT_ERROR("PrescribeStateType should be 'None' or 'Init' for thickness but got {}", + InString); + return PrescribeStateType::Invalid; +} +PrescribeStateType getPrescribeVelocityTypeFromStr(const std::string &InString) { + + if (InString == "None") { + return PrescribeStateType::None; + } + if (InString == "Init") { + return PrescribeStateType::Init; + } + + ABORT_ERROR("PrescribeStateType should be 'None' or 'Init' for velocity but got {}", + InString); + return PrescribeStateType::Invalid; +} + //------------------------------------------------------------------------------ // Constructors and creation methods. @@ -119,6 +149,11 @@ TimeStepper *TimeStepper::create( ABORT_ERROR("Unknown time stepping method"); } + NewTimeStepper->PrescribeThicknessMode = + DefaultPrescribeThicknessMode; + NewTimeStepper->PrescribeVelocityMode = + DefaultPrescribeVelocityMode; + // Attach data pointers NewTimeStepper->attachData(InTend, InAuxState, InMesh, InVCoord, InMeshHalo); @@ -170,6 +205,11 @@ TimeStepper *TimeStepper::create( ABORT_ERROR("Unknown time stepping method"); } + NewTimeStepper->PrescribeThicknessMode = + DefaultPrescribeThicknessMode; + NewTimeStepper->PrescribeVelocityMode = + DefaultPrescribeVelocityMode; + // Store instance AllTimeSteppers.emplace(InName, NewTimeStepper); @@ -298,6 +338,22 @@ void TimeStepper::init1() { StopTime = StopTime2; } + Config StateConfig("State"); + Error StateErr = OmegaConfig->get(StateConfig); + if (StateErr.isSuccess()) { + std::string ThicknessMode; + if (StateConfig.get("PrescribeThicknessType", ThicknessMode).isSuccess()) { + TimeStepper::DefaultPrescribeThicknessMode = + getPrescribeThicknessTypeFromStr(ThicknessMode); + } + + std::string VelocityMode; + if (StateConfig.get("PrescribeVelocityType", VelocityMode).isSuccess()) { + TimeStepper::DefaultPrescribeVelocityMode = + getPrescribeVelocityTypeFromStr(VelocityMode); + } + } + // Now that all the inputs are defined, create the default time stepper // Use the partial creation function for only the time info. Data // pointers will be attached in phase 2 initialization @@ -465,50 +521,64 @@ void TimeStepper::updateStateByTend(OceanState *State1, int TimeLevel1, void TimeStepper::prescribeThickness(OceanState *State1, int TimeLevel1, OceanState *State2, int TimeLevel2) const { - Array2DReal LayerThick1 = State1->getLayerThickness(TimeLevel1); - Array2DReal LayerThick2 = State2->getLayerThickness(TimeLevel2); - - OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); - OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - - parallelForOuter( - "prescribeThickness", {Mesh->NCellsAll}, - KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { - const int KMin = MinLayerCell(ICell); - const int KMax = MaxLayerCell(ICell); - const int KRange = vertRange(KMin, KMax); + if (PrescribeThicknessMode == PrescribeStateType::None) { + return; + } - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - const int K = KMin + KChunk; - LayerThick1(ICell, K) = LayerThick2(ICell, K); - }); - }); + if (PrescribeThicknessMode == PrescribeStateType::Init) { + Array2DReal LayerThick1 = State1->getLayerThickness(TimeLevel1); + Array2DReal LayerThick2 = State2->getLayerThickness(TimeLevel2); + + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + + parallelForOuter( + "prescribeThickness", {Mesh->NCellsAll}, + KOKKOS_LAMBDA(int ICell, const TeamMember &Team) { + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + const int KRange = vertRange(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + LayerThick1(ICell, K) = LayerThick2(ICell, K); + }); + }); + return; + } } //------------------------------------------------------------------------------ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, OceanState *State2, int TimeLevel2) const { - Array2DReal NormalVel1 = State1->getNormalVelocity(TimeLevel1); - Array2DReal NormalVel2 = State2->getNormalVelocity(TimeLevel2); - - OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); - OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); - - parallelForOuter( - "prescribeVelocity", {Mesh->NEdgesAll}, - KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { - const int KMin = MinLayerEdgeBot(IEdge); - const int KMax = MaxLayerEdgeTop(IEdge); - const int KRange = vertRange(KMin, KMax); + if (PrescribeVelocityMode == PrescribeStateType::None) { + return; + } - parallelForInner( - Team, KRange, INNER_LAMBDA(int KChunk) { - const int K = KMin + KChunk; - NormalVel1(IEdge, K) = NormalVel2(IEdge, K); - }); - }); + if (PrescribeVelocityMode == PrescribeStateType::Init) { + Array2DReal NormalVel1 = State1->getNormalVelocity(TimeLevel1); + Array2DReal NormalVel2 = State2->getNormalVelocity(TimeLevel2); + + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + + parallelForOuter( + "prescribeVelocity", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRange(KMin, KMax); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + NormalVel1(IEdge, K) = NormalVel2(IEdge, K); + }); + }); + return; + } } //------------------------------------------------------------------------------ diff --git a/components/omega/src/timeStepping/TimeStepper.h b/components/omega/src/timeStepping/TimeStepper.h index 3696678fa0cb..92f631ad70af 100644 --- a/components/omega/src/timeStepping/TimeStepper.h +++ b/components/omega/src/timeStepping/TimeStepper.h @@ -62,6 +62,14 @@ enum class TimeStepperType { Invalid }; +/// An enum describing how a state variable should be prescribed from the +/// reference time level +enum class PrescribeStateType { + None, + Init, + Invalid +}; + //------------------------------------------------------------------------------ // Utility routine /// Translate string for time stepper type into enum @@ -69,6 +77,11 @@ TimeStepperType getTimeStepperFromStr( const std::string &InString ///< [in] choice of time stepping method ); +/// Translate string for prescribe state type into enum +PrescribeStateType getPrescribeStateTypeFromStr( + const std::string &InString ///< [in] choice of prescribe method +); + //------------------------------------------------------------------------------ /// A base class for Omega time steppers /// @@ -277,6 +290,10 @@ class TimeStepper { /// Time step TimeInterval TimeStep; + /// Prescribe state configuration + PrescribeStateType PrescribeThicknessMode; + PrescribeStateType PrescribeVelocityMode; + /// Start time TimeInstant StartTime; @@ -320,6 +337,9 @@ class TimeStepper { static TimeStepper *DefaultTimeStepper; /// All defined time steppers static std::map> AllTimeSteppers; + /// Prescribe modes parsed from config + static PrescribeStateType DefaultPrescribeThicknessMode; + static PrescribeStateType DefaultPrescribeVelocityMode; }; } // namespace OMEGA From d3fb28a39f60ea46ace2599acd27539306781bec Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 15:59:55 -0700 Subject: [PATCH 04/12] Add simTime as an argument to PrescribeState, PrescribeVelocity --- .../timeStepping/ForwardBackwardStepper.cpp | 3 ++- .../src/timeStepping/RungeKutta2Stepper.cpp | 2 +- .../src/timeStepping/RungeKutta4Stepper.cpp | 6 +++--- .../omega/src/timeStepping/TimeStepper.cpp | 8 +++++--- .../omega/src/timeStepping/TimeStepper.h | 18 ++++++++++-------- 5 files changed, 21 insertions(+), 16 deletions(-) diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index 00da6d0a23e2..b59e6c8b0580 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -69,7 +69,8 @@ void ForwardBackwardStepper::doStep( State, ThickCurLevel, TimeStep); prescribeThickness(State, NextLevel, State, CurLevel); - prescribeVelocity(State, NextLevel, State, CurLevel); + prescribeVelocity(State, NextLevel, State, CurLevel, + SimTime + TimeStep); // Update time levels (New -> Old) of prognostic variables with halo // exchanges diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index d34a3de804eb..bde140205cd7 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -56,7 +56,7 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state updateTracersByTend(NextTracerArray, CurTracerArray, State, NextLevel, State, CurLevel, TimeStep); - prescribeState(State, NextLevel, State, CurLevel); + prescribeState(State, NextLevel, State, CurLevel, SimTime + TimeStep); // Update time levels (New -> Old) of prognostic variables with halo // exchanges diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index 6742a11263f7..b6b3e94c6c41 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -88,7 +88,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, CurLevel, RKB[Stage] * TimeStep); - prescribeState(State, NextLevel, State, CurLevel); + prescribeState(State, NextLevel, State, CurLevel, StageTime); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); } else { // every other stage does: @@ -97,7 +97,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state // q^{n+1} += RKB[stage] * dt * R^{(s)} updateStateByTend(ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); - prescribeState(State, NextLevel, State, CurLevel); + prescribeState(State, NextLevel, State, CurLevel, StageTime); updateTracersByTend(ProvisTracers, CurTracerArray, ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); @@ -111,7 +111,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state CurLevel, CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, NextLevel, RKB[Stage] * TimeStep); - prescribeState(State, NextLevel, State, CurLevel); + prescribeState(State, NextLevel, State, CurLevel, StageTime); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); } } diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 2c7d44f1656b..28b911db43d5 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -551,7 +551,8 @@ void TimeStepper::prescribeThickness(OceanState *State1, int TimeLevel1, //------------------------------------------------------------------------------ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, - OceanState *State2, int TimeLevel2) const { + OceanState *State2, int TimeLevel2, + const TimeInstant &SimTime) const { if (PrescribeVelocityMode == PrescribeStateType::None) { return; @@ -583,9 +584,10 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, //------------------------------------------------------------------------------ void TimeStepper::prescribeState(OceanState *State1, int TimeLevel1, - OceanState *State2, int TimeLevel2) const { + OceanState *State2, int TimeLevel2, + const TimeInstant &SimTime) const { prescribeThickness(State1, TimeLevel1, State2, TimeLevel2); - prescribeVelocity(State1, TimeLevel1, State2, TimeLevel2); + prescribeVelocity(State1, TimeLevel1, State2, TimeLevel2, SimTime); } //------------------------------------------------------------------------------ diff --git a/components/omega/src/timeStepping/TimeStepper.h b/components/omega/src/timeStepping/TimeStepper.h index 92f631ad70af..f3e4e8cd98b2 100644 --- a/components/omega/src/timeStepping/TimeStepper.h +++ b/components/omega/src/timeStepping/TimeStepper.h @@ -228,18 +228,20 @@ class TimeStepper { /// Resets the velocity at the working time level to the initial condition. void prescribeVelocity( - OceanState *State1, ///< [out] destination state - int TimeLevel1, ///< [in] time level index of the reference data - OceanState *State2, ///< [in] source state (initial condition) - int TimeLevel2 ///< [in] time level index of the destination data + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2, ///< [in] time level index of the destination data + const TimeInstant &SimTime ///< [in] current simulation time ) const; /// Reset thickness and velocity to their initial values void prescribeState( - OceanState *State1, ///< [out] destination state - int TimeLevel1, ///< [in] time level index of the reference data - OceanState *State2, ///< [in] source state (initial condition) - int TimeLevel2 ///< [in] time level index of the destination data + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2, ///< [in] time level index of the destination data + const TimeInstant &SimTime ///< [in] current simulation time ) const; /// Updates tracers From 18b8a719a0251adb0fc188921c5ff1e4709ee11e Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 16:01:30 -0700 Subject: [PATCH 05/12] Fixup Init; Curr, Prev superimposed --- components/omega/src/timeStepping/TimeStepper.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 28b911db43d5..c8743082debb 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -575,7 +575,7 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { const int K = KMin + KChunk; - NormalVel1(IEdge, K) = NormalVel2(IEdge, K); + NormalVel2(IEdge, K) = NormalVel1(IEdge, K); }); }); return; From 9b388bdafb06fe29fbb4f9a4c7024e93871a579a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 16:44:43 -0700 Subject: [PATCH 06/12] Add non-divergent,divergent flow types --- .../omega/src/timeStepping/TimeStepper.cpp | 98 ++++++++++++++++++- .../omega/src/timeStepping/TimeStepper.h | 2 + 2 files changed, 97 insertions(+), 3 deletions(-) diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index c8743082debb..e11e1108e7ea 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -10,6 +10,7 @@ #include "ForwardBackwardStepper.h" #include "RungeKutta2Stepper.h" #include "RungeKutta4Stepper.h" +#include "Logging.h" namespace OMEGA { //------------------------------------------------------------------------------ @@ -64,12 +65,15 @@ PrescribeStateType getPrescribeVelocityTypeFromStr(const std::string &InString) if (InString == "None") { return PrescribeStateType::None; - } - if (InString == "Init") { + }else if (InString == "Init") { return PrescribeStateType::Init; + }else if (InString == "NonDivergent") { + return PrescribeStateType::NonDivergent; + }else if (InString == "Divergent") { + return PrescribeStateType::Divergent; } - ABORT_ERROR("PrescribeStateType should be 'None' or 'Init' for velocity but got {}", + ABORT_ERROR("PrescribeStateType should be 'None', 'Init', 'NonDivergent' or 'Divergent' for velocity but got {}", InString); return PrescribeStateType::Invalid; } @@ -579,6 +583,94 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, }); }); return; + } else if (PrescribeVelocityMode == PrescribeStateType::NonDivergent) { + Array2DReal NormalVel2 = State2->getNormalVelocity(TimeLevel2); + + OMEGA_SCOPE(LatEdge, Mesh->LatEdgeH); + OMEGA_SCOPE(LonEdge, Mesh->LonEdgeH); + OMEGA_SCOPE(AngleEdge, Mesh->AngleEdgeH); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBotH); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTopH); + + const Clock *ModelClock = StepClock.get(); + R8 ElapsedTimeSec; + TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getCurrentTime(); + ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds); + + const R8 Tau = 12. * Day2Sec; // 12 days in seconds + const R8 TSim = ElapsedTimeSec; + + parallelForOuter( + "prescribeVelocityNonDivergent", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRange(KMin, KMax); + + const R8 lon_p = + LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; + const R8 u = (10.0 / Tau) * + (Kokkos::pow(sin(lon_p), 2) * + sin(2.0 * LatEdge(IEdge)) * + cos(Pi * TSim / Tau) + + 2.0 * Pi * cos(LatEdge(IEdge))); + const R8 v = (10.0 / Tau) * sin(2.0 * lon_p) * + cos(LatEdge(IEdge)) * cos(Pi * TSim / Tau); + const R8 normalVel = REarth * ( + u * cos(AngleEdge(IEdge)) + v * sin(AngleEdge(IEdge))); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + NormalVel2(IEdge, K) = normalVel; + }); + }); + return; + } else if (PrescribeVelocityMode == PrescribeStateType::Divergent) { + Array2DReal NormalVel2 = State2->getNormalVelocity(TimeLevel2); + + OMEGA_SCOPE(LatEdge, Mesh->LatEdgeH); + OMEGA_SCOPE(LonEdge, Mesh->LonEdgeH); + OMEGA_SCOPE(AngleEdge, Mesh->AngleEdgeH); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBotH); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTopH); + + const Clock *ModelClock = StepClock.get(); + R8 ElapsedTimeSec; + TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getCurrentTime(); + ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds); + + const R8 Tau = 12. * Day2Sec; // 14 days in seconds + const R8 TSim = ElapsedTimeSec; + + parallelForOuter( + "prescribeVelocityDivergent", {Mesh->NEdgesAll}, + KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRange(KMin, KMax); + + const R8 lon_p = + LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; + const R8 u = (1.0 / Tau) * + (-10.0 / 2.0) * Kokkos::pow(sin(lon_p / 2), 2) * + sin(2.0 * LatEdge(IEdge)) * + Kokkos::pow(cos(LatEdge(IEdge)), 2) * + cos(Pi * TSim / Tau) + + 2.0 * Pi * cos(LatEdge(IEdge)); + const R8 v = (10.0 / (4 * Tau)) * + sin(lon_p) * Kokkos::pow(cos(LatEdge(IEdge)), 3) * + cos(Pi * TSim / Tau); + const R8 normalVel = REarth * ( + u * cos(AngleEdge(IEdge)) + v * sin(AngleEdge(IEdge))); + + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + const int K = KMin + KChunk; + NormalVel2(IEdge, K) = normalVel; + }); + }); + return; } } diff --git a/components/omega/src/timeStepping/TimeStepper.h b/components/omega/src/timeStepping/TimeStepper.h index f3e4e8cd98b2..da79672d4002 100644 --- a/components/omega/src/timeStepping/TimeStepper.h +++ b/components/omega/src/timeStepping/TimeStepper.h @@ -67,6 +67,8 @@ enum class TimeStepperType { enum class PrescribeStateType { None, Init, + NonDivergent, + Divergent, Invalid }; From f8d86b884c6b03b3003aed6a628c5f8f37e9d9b7 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 19 Mar 2026 09:01:02 -0700 Subject: [PATCH 07/12] Fixup nondiv,div velocity --- components/omega/src/timeStepping/TimeStepper.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index e11e1108e7ea..353e6398a439 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -609,8 +609,8 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const R8 lon_p = LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; - const R8 u = (10.0 / Tau) * - (Kokkos::pow(sin(lon_p), 2) * + const R8 u = (1 / Tau) * + (10.0 * Kokkos::pow(sin(lon_p), 2) * sin(2.0 * LatEdge(IEdge)) * cos(Pi * TSim / Tau) + 2.0 * Pi * cos(LatEdge(IEdge))); @@ -653,14 +653,14 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const R8 lon_p = LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; const R8 u = (1.0 / Tau) * - (-10.0 / 2.0) * Kokkos::pow(sin(lon_p / 2), 2) * + (-10.0 * Kokkos::pow(sin(lon_p / 2), 2) * sin(2.0 * LatEdge(IEdge)) * Kokkos::pow(cos(LatEdge(IEdge)), 2) * - cos(Pi * TSim / Tau) + + cos(Pi * TSim / Tau)) + 2.0 * Pi * cos(LatEdge(IEdge)); - const R8 v = (10.0 / (4 * Tau)) * - sin(lon_p) * Kokkos::pow(cos(LatEdge(IEdge)), 3) * - cos(Pi * TSim / Tau); + const R8 v = (10.0 / (2 * Tau)) * + (sin(lon_p) * Kokkos::pow(cos(LatEdge(IEdge)), 3) * + cos(Pi * TSim / Tau)); const R8 normalVel = REarth * ( u * cos(AngleEdge(IEdge)) + v * sin(AngleEdge(IEdge))); From f1755eaba567d4cbb0bfa69d1b080d80b74666bc Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 8 Apr 2026 09:27:44 -0700 Subject: [PATCH 08/12] Fix divergent parentheses --- components/omega/src/timeStepping/TimeStepper.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 353e6398a439..9dbce2f07504 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -653,14 +653,15 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const R8 lon_p = LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; const R8 u = (1.0 / Tau) * - (-10.0 * Kokkos::pow(sin(lon_p / 2), 2) * + (-5.0 * Kokkos::pow(sin(lon_p / 2), 2) * sin(2.0 * LatEdge(IEdge)) * Kokkos::pow(cos(LatEdge(IEdge)), 2) * - cos(Pi * TSim / Tau)) + - 2.0 * Pi * cos(LatEdge(IEdge)); - const R8 v = (10.0 / (2 * Tau)) * - (sin(lon_p) * Kokkos::pow(cos(LatEdge(IEdge)), 3) * - cos(Pi * TSim / Tau)); + cos(Pi * TSim / Tau) + + 2.0 * Pi * cos(LatEdge(IEdge))); + const R8 v = ((2.5 / Tau) * + sin(lon_p) * + Kokkos::pow(cos(LatEdge(IEdge)), 3) * + cos(Pi * TSim / Tau)); const R8 normalVel = REarth * ( u * cos(AngleEdge(IEdge)) + v * sin(AngleEdge(IEdge))); From 718bfdabf192f094e47fd3a2ccee7cdd44759bed Mon Sep 17 00:00:00 2001 From: James Overfelt Date: Sat, 25 Apr 2026 10:20:10 -0700 Subject: [PATCH 09/12] Fix ElapsedTimeInterval for time dependent prescribed velocities. ElapsedTimeInterval is the difference between the start time and the current time. But the wrong start time was used so the interval was being set to a sub-step time interval. --- components/omega/src/timeStepping/TimeStepper.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 9dbce2f07504..35f2ffd6f059 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -594,7 +594,7 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const Clock *ModelClock = StepClock.get(); R8 ElapsedTimeSec; - TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getCurrentTime(); + TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getStartTime(); ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds); const R8 Tau = 12. * Day2Sec; // 12 days in seconds @@ -637,7 +637,7 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const Clock *ModelClock = StepClock.get(); R8 ElapsedTimeSec; - TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getCurrentTime(); + TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getStartTime(); ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds); const R8 Tau = 12. * Day2Sec; // 14 days in seconds From e5289d42adca7ffe87fabe6ebfd3a1e4abbb2b74 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 29 Apr 2026 13:13:52 -0400 Subject: [PATCH 10/12] Fix prescribeState in time steppers - Updated the location and time information for prescribeState in the time steppers. --- .../omega/src/timeStepping/ForwardBackwardStepper.cpp | 8 ++++---- components/omega/src/timeStepping/RungeKutta2Stepper.cpp | 6 ++++-- components/omega/src/timeStepping/RungeKutta4Stepper.cpp | 7 ++++--- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index b59e6c8b0580..2458017d6542 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -45,6 +45,8 @@ void ForwardBackwardStepper::doStep( if (AuxState == nullptr) LOG_CRITICAL("Invalid AuxState"); + prescribeVelocity(State, CurLevel, State, CurLevel,SimTime); + // R_u^{n} = RHS_u(u^{n}, h^{n}, t^{n}) Tend->computeVelocityTendencies(State, AuxState, CurTracerArray, ThickCurLevel, VelCurLevel, TracerCurLevel, @@ -60,6 +62,8 @@ void ForwardBackwardStepper::doStep( // h^{n+1} = h^{n} + R_h^{n} updateThicknessByTend(State, ThickNextLevel, State, ThickCurLevel, TimeStep); + prescribeThickness(State, CurLevel, State, CurLevel); + // R_phi^{n} = RHS_phi(u^{n+1}, h^{n+1}, phi^{n}, t^{n}) Tend->computeTracerTendencies(State, AuxState, CurTracerArray, ThickNextLevel, VelNextLevel, SimTime); @@ -68,10 +72,6 @@ void ForwardBackwardStepper::doStep( updateTracersByTend(NextTracerArray, CurTracerArray, State, ThickNextLevel, State, ThickCurLevel, TimeStep); - prescribeThickness(State, NextLevel, State, CurLevel); - prescribeVelocity(State, NextLevel, State, CurLevel, - SimTime + TimeStep); - // Update time levels (New -> Old) of prognostic variables with halo // exchanges const MPI_Comm Comm = MeshHalo->getComm(); diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index bde140205cd7..5a325b45b502 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -34,6 +34,8 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); + prescribeState(State, CurLevel, State, CurLevel, SimTime); + // q = (h,u,phi) // R_q^{n} = RHS_q(u^{n}, h^{n}, phi^{n}, t^{n}) Tend->computeAllTendencies(State, AuxState, CurTracerArray, CurLevel, @@ -47,6 +49,8 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state State->exchangeHalo(NextLevel); MeshHalo->exchangeFullArrayHalo(NextTracerArray, OnCell); + prescribeState(State, NextLevel, State, CurLevel, SimTime + 0.5 * TimeStep); + // R_q^{n+0.5} = RHS_q(u^{n+0.5}, h^{n+0.5}, phi^{n+0.5}, t^{n+0.5}) Tend->computeAllTendencies(State, AuxState, NextTracerArray, NextLevel, NextLevel, NextLevel, SimTime + 0.5 * TimeStep); @@ -56,8 +60,6 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state updateTracersByTend(NextTracerArray, CurTracerArray, State, NextLevel, State, CurLevel, TimeStep); - prescribeState(State, NextLevel, State, CurLevel, SimTime + TimeStep); - // Update time levels (New -> Old) of prognostic variables with halo // exchanges const MPI_Comm Comm = MeshHalo->getComm(); diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index b6b3e94c6c41..945982a8a80c 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -76,6 +76,7 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); + TimeInstant ForcingStageTime = SimTime; for (int Stage = 0; Stage < NStages; ++Stage) { const TimeInstant StageTime = SimTime + RKC[Stage] * TimeStep; @@ -84,11 +85,11 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state // q^{n+1} = q^{n} + dt * RKB[0] * R^{(0)} if (Stage == 0) { weightTracers(NextTracerArray, CurTracerArray, State, CurLevel); + prescribeState(State, CurLevel, State, CurLevel, ForcingStageTime); Tend->computeAllTendencies(State, AuxState, CurTracerArray, CurLevel, CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, CurLevel, RKB[Stage] * TimeStep); - prescribeState(State, NextLevel, State, CurLevel, StageTime); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); } else { // every other stage does: @@ -97,7 +98,8 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state // q^{n+1} += RKB[stage] * dt * R^{(s)} updateStateByTend(ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); - prescribeState(State, NextLevel, State, CurLevel, StageTime); + ForcingStageTime += RKA[Stage] * TimeStep; + prescribeState(ProvisState, CurLevel, ProvisState, CurLevel, ForcingStageTime); updateTracersByTend(ProvisTracers, CurTracerArray, ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); @@ -111,7 +113,6 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state CurLevel, CurLevel, CurLevel, StageTime); updateStateByTend(State, NextLevel, State, NextLevel, RKB[Stage] * TimeStep); - prescribeState(State, NextLevel, State, CurLevel, StageTime); accumulateTracersUpdate(NextTracerArray, RKB[Stage] * TimeStep); } } From d3461a6892969a78501a718646bcf31370a5c557 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 29 Apr 2026 13:22:24 -0400 Subject: [PATCH 11/12] Fix linting issues --- .../timeStepping/ForwardBackwardStepper.cpp | 2 +- .../src/timeStepping/RungeKutta4Stepper.cpp | 7 +- .../omega/src/timeStepping/TimeStepper.cpp | 105 +++++++++--------- .../omega/src/timeStepping/TimeStepper.h | 72 ++++++------ 4 files changed, 89 insertions(+), 97 deletions(-) diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index 2458017d6542..47550e70debf 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -45,7 +45,7 @@ void ForwardBackwardStepper::doStep( if (AuxState == nullptr) LOG_CRITICAL("Invalid AuxState"); - prescribeVelocity(State, CurLevel, State, CurLevel,SimTime); + prescribeVelocity(State, CurLevel, State, CurLevel, SimTime); // R_u^{n} = RHS_u(u^{n}, h^{n}, t^{n}) Tend->computeVelocityTendencies(State, AuxState, CurTracerArray, diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index 945982a8a80c..b1cf5353d249 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -74,8 +74,8 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state const int CurLevel = 0; const int NextLevel = 1; - Array3DReal CurTracerArray = Tracers::getAll(CurLevel); - Array3DReal NextTracerArray = Tracers::getAll(NextLevel); + Array3DReal CurTracerArray = Tracers::getAll(CurLevel); + Array3DReal NextTracerArray = Tracers::getAll(NextLevel); TimeInstant ForcingStageTime = SimTime; for (int Stage = 0; Stage < NStages; ++Stage) { @@ -99,7 +99,8 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state updateStateByTend(ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); ForcingStageTime += RKA[Stage] * TimeStep; - prescribeState(ProvisState, CurLevel, ProvisState, CurLevel, ForcingStageTime); + prescribeState(ProvisState, CurLevel, ProvisState, CurLevel, + ForcingStageTime); updateTracersByTend(ProvisTracers, CurTracerArray, ProvisState, CurLevel, State, CurLevel, RKA[Stage] * TimeStep); diff --git a/components/omega/src/timeStepping/TimeStepper.cpp b/components/omega/src/timeStepping/TimeStepper.cpp index 35f2ffd6f059..b78319360618 100644 --- a/components/omega/src/timeStepping/TimeStepper.cpp +++ b/components/omega/src/timeStepping/TimeStepper.cpp @@ -8,9 +8,9 @@ #include "Config.h" #include "Error.h" #include "ForwardBackwardStepper.h" +#include "Logging.h" #include "RungeKutta2Stepper.h" #include "RungeKutta4Stepper.h" -#include "Logging.h" namespace OMEGA { //------------------------------------------------------------------------------ @@ -21,9 +21,9 @@ TimeStepper *TimeStepper::DefaultTimeStepper = nullptr; std::map> TimeStepper::AllTimeSteppers; PrescribeStateType TimeStepper::DefaultPrescribeThicknessMode = - PrescribeStateType::None; + PrescribeStateType::None; PrescribeStateType TimeStepper::DefaultPrescribeVelocityMode = - PrescribeStateType::None; + PrescribeStateType::None; //------------------------------------------------------------------------------ // utility functions @@ -48,7 +48,8 @@ TimeStepperType getTimeStepperFromStr(const std::string &InString) { return TimeStepperChoice; } -PrescribeStateType getPrescribeThicknessTypeFromStr(const std::string &InString) { +PrescribeStateType +getPrescribeThicknessTypeFromStr(const std::string &InString) { if (InString == "None") { return PrescribeStateType::None; @@ -57,23 +58,26 @@ PrescribeStateType getPrescribeThicknessTypeFromStr(const std::string &InString) return PrescribeStateType::Init; } - ABORT_ERROR("PrescribeStateType should be 'None' or 'Init' for thickness but got {}", - InString); + ABORT_ERROR( + "PrescribeStateType should be 'None' or 'Init' for thickness but got {}", + InString); return PrescribeStateType::Invalid; } -PrescribeStateType getPrescribeVelocityTypeFromStr(const std::string &InString) { +PrescribeStateType +getPrescribeVelocityTypeFromStr(const std::string &InString) { if (InString == "None") { return PrescribeStateType::None; - }else if (InString == "Init") { + } else if (InString == "Init") { return PrescribeStateType::Init; - }else if (InString == "NonDivergent") { + } else if (InString == "NonDivergent") { return PrescribeStateType::NonDivergent; - }else if (InString == "Divergent") { + } else if (InString == "Divergent") { return PrescribeStateType::Divergent; } - ABORT_ERROR("PrescribeStateType should be 'None', 'Init', 'NonDivergent' or 'Divergent' for velocity but got {}", + ABORT_ERROR("PrescribeStateType should be 'None', 'Init', 'NonDivergent' or " + "'Divergent' for velocity but got {}", InString); return PrescribeStateType::Invalid; } @@ -153,10 +157,8 @@ TimeStepper *TimeStepper::create( ABORT_ERROR("Unknown time stepping method"); } - NewTimeStepper->PrescribeThicknessMode = - DefaultPrescribeThicknessMode; - NewTimeStepper->PrescribeVelocityMode = - DefaultPrescribeVelocityMode; + NewTimeStepper->PrescribeThicknessMode = DefaultPrescribeThicknessMode; + NewTimeStepper->PrescribeVelocityMode = DefaultPrescribeVelocityMode; // Attach data pointers NewTimeStepper->attachData(InTend, InAuxState, InMesh, InVCoord, InMeshHalo); @@ -209,10 +211,8 @@ TimeStepper *TimeStepper::create( ABORT_ERROR("Unknown time stepping method"); } - NewTimeStepper->PrescribeThicknessMode = - DefaultPrescribeThicknessMode; - NewTimeStepper->PrescribeVelocityMode = - DefaultPrescribeVelocityMode; + NewTimeStepper->PrescribeThicknessMode = DefaultPrescribeThicknessMode; + NewTimeStepper->PrescribeVelocityMode = DefaultPrescribeVelocityMode; // Store instance AllTimeSteppers.emplace(InName, NewTimeStepper); @@ -346,7 +346,8 @@ void TimeStepper::init1() { Error StateErr = OmegaConfig->get(StateConfig); if (StateErr.isSuccess()) { std::string ThicknessMode; - if (StateConfig.get("PrescribeThicknessType", ThicknessMode).isSuccess()) { + if (StateConfig.get("PrescribeThicknessType", ThicknessMode) + .isSuccess()) { TimeStepper::DefaultPrescribeThicknessMode = getPrescribeThicknessTypeFromStr(ThicknessMode); } @@ -545,7 +546,7 @@ void TimeStepper::prescribeThickness(OceanState *State1, int TimeLevel1, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - const int K = KMin + KChunk; + const int K = KMin + KChunk; LayerThick1(ICell, K) = LayerThick2(ICell, K); }); }); @@ -555,8 +556,8 @@ void TimeStepper::prescribeThickness(OceanState *State1, int TimeLevel1, //------------------------------------------------------------------------------ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, - OceanState *State2, int TimeLevel2, - const TimeInstant &SimTime) const { + OceanState *State2, int TimeLevel2, + const TimeInstant &SimTime) const { if (PrescribeVelocityMode == PrescribeStateType::None) { return; @@ -578,7 +579,7 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - const int K = KMin + KChunk; + const int K = KMin + KChunk; NormalVel2(IEdge, K) = NormalVel1(IEdge, K); }); }); @@ -597,8 +598,8 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getStartTime(); ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds); - const R8 Tau = 12. * Day2Sec; // 12 days in seconds - const R8 TSim = ElapsedTimeSec; + const R8 Tau = 12. * Day2Sec; // 12 days in seconds + const R8 TSim = ElapsedTimeSec; parallelForOuter( "prescribeVelocityNonDivergent", {Mesh->NEdgesAll}, @@ -607,21 +608,19 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const int KMax = MaxLayerEdgeTop(IEdge); const int KRange = vertRange(KMin, KMax); - const R8 lon_p = - LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; - const R8 u = (1 / Tau) * - (10.0 * Kokkos::pow(sin(lon_p), 2) * - sin(2.0 * LatEdge(IEdge)) * - cos(Pi * TSim / Tau) + - 2.0 * Pi * cos(LatEdge(IEdge))); - const R8 v = (10.0 / Tau) * sin(2.0 * lon_p) * + const R8 lon_p = LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; + const R8 u = (1 / Tau) * (10.0 * Kokkos::pow(sin(lon_p), 2) * + sin(2.0 * LatEdge(IEdge)) * + cos(Pi * TSim / Tau) + + 2.0 * Pi * cos(LatEdge(IEdge))); + const R8 v = (10.0 / Tau) * sin(2.0 * lon_p) * cos(LatEdge(IEdge)) * cos(Pi * TSim / Tau); - const R8 normalVel = REarth * ( - u * cos(AngleEdge(IEdge)) + v * sin(AngleEdge(IEdge))); + const R8 normalVel = REarth * (u * cos(AngleEdge(IEdge)) + + v * sin(AngleEdge(IEdge))); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - const int K = KMin + KChunk; + const int K = KMin + KChunk; NormalVel2(IEdge, K) = normalVel; }); }); @@ -640,8 +639,8 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, TimeInterval ElapsedTimeInterval = SimTime - ModelClock->getStartTime(); ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds); - const R8 Tau = 12. * Day2Sec; // 14 days in seconds - const R8 TSim = ElapsedTimeSec; + const R8 Tau = 12. * Day2Sec; // 14 days in seconds + const R8 TSim = ElapsedTimeSec; parallelForOuter( "prescribeVelocityDivergent", {Mesh->NEdgesAll}, @@ -650,24 +649,22 @@ void TimeStepper::prescribeVelocity(OceanState *State1, int TimeLevel1, const int KMax = MaxLayerEdgeTop(IEdge); const int KRange = vertRange(KMin, KMax); - const R8 lon_p = - LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; - const R8 u = (1.0 / Tau) * - (-5.0 * Kokkos::pow(sin(lon_p / 2), 2) * - sin(2.0 * LatEdge(IEdge)) * - Kokkos::pow(cos(LatEdge(IEdge)), 2) * - cos(Pi * TSim / Tau) + - 2.0 * Pi * cos(LatEdge(IEdge))); - const R8 v = ((2.5 / Tau) * - sin(lon_p) * - Kokkos::pow(cos(LatEdge(IEdge)), 3) * - cos(Pi * TSim / Tau)); - const R8 normalVel = REarth * ( - u * cos(AngleEdge(IEdge)) + v * sin(AngleEdge(IEdge))); + const R8 lon_p = LonEdge(IEdge) - 2.0 * Pi * TSim / Tau; + const R8 u = + (1.0 / Tau) * (-5.0 * Kokkos::pow(sin(lon_p / 2), 2) * + sin(2.0 * LatEdge(IEdge)) * + Kokkos::pow(cos(LatEdge(IEdge)), 2) * + cos(Pi * TSim / Tau) + + 2.0 * Pi * cos(LatEdge(IEdge))); + const R8 v = + ((2.5 / Tau) * sin(lon_p) * + Kokkos::pow(cos(LatEdge(IEdge)), 3) * cos(Pi * TSim / Tau)); + const R8 normalVel = REarth * (u * cos(AngleEdge(IEdge)) + + v * sin(AngleEdge(IEdge))); parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - const int K = KMin + KChunk; + const int K = KMin + KChunk; NormalVel2(IEdge, K) = normalVel; }); }); diff --git a/components/omega/src/timeStepping/TimeStepper.h b/components/omega/src/timeStepping/TimeStepper.h index da79672d4002..4d591792b6b9 100644 --- a/components/omega/src/timeStepping/TimeStepper.h +++ b/components/omega/src/timeStepping/TimeStepper.h @@ -64,13 +64,7 @@ enum class TimeStepperType { /// An enum describing how a state variable should be prescribed from the /// reference time level -enum class PrescribeStateType { - None, - Init, - NonDivergent, - Divergent, - Invalid -}; +enum class PrescribeStateType { None, Init, NonDivergent, Divergent, Invalid }; //------------------------------------------------------------------------------ // Utility routine @@ -219,32 +213,32 @@ class TimeStepper { TimeInterval Coeff ///< [in] time-related coeff for tendency ) const; - /// Resets the layer thickness at the working time level to the initial - /// condition stored at the reference time level. - void prescribeThickness( - OceanState *State1, ///< [out] destination state - int TimeLevel1, ///< [in] time level index of the reference data - OceanState *State2, ///< [in] source state (initial condition) - int TimeLevel2 ///< [in] time level index of the destination data - ) const; - - /// Resets the velocity at the working time level to the initial condition. - void prescribeVelocity( - OceanState *State1, ///< [out] destination state - int TimeLevel1, ///< [in] time level index of the reference data - OceanState *State2, ///< [in] source state (initial condition) - int TimeLevel2, ///< [in] time level index of the destination data - const TimeInstant &SimTime ///< [in] current simulation time - ) const; - - /// Reset thickness and velocity to their initial values - void prescribeState( - OceanState *State1, ///< [out] destination state - int TimeLevel1, ///< [in] time level index of the reference data - OceanState *State2, ///< [in] source state (initial condition) - int TimeLevel2, ///< [in] time level index of the destination data - const TimeInstant &SimTime ///< [in] current simulation time - ) const; + /// Resets the layer thickness at the working time level to the initial + /// condition stored at the reference time level. + void prescribeThickness( + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2 ///< [in] time level index of the destination data + ) const; + + /// Resets the velocity at the working time level to the initial condition. + void prescribeVelocity( + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2, ///< [in] time level index of the destination data + const TimeInstant &SimTime ///< [in] current simulation time + ) const; + + /// Reset thickness and velocity to their initial values + void prescribeState( + OceanState *State1, ///< [out] destination state + int TimeLevel1, ///< [in] time level index of the reference data + OceanState *State2, ///< [in] source state (initial condition) + int TimeLevel2, ///< [in] time level index of the destination data + const TimeInstant &SimTime ///< [in] current simulation time + ) const; /// Updates tracers /// NextTracers = (CurTracers * LayerThickness2(TimeLevel2)) + @@ -294,9 +288,9 @@ class TimeStepper { /// Time step TimeInterval TimeStep; - /// Prescribe state configuration - PrescribeStateType PrescribeThicknessMode; - PrescribeStateType PrescribeVelocityMode; + /// Prescribe state configuration + PrescribeStateType PrescribeThicknessMode; + PrescribeStateType PrescribeVelocityMode; /// Start time TimeInstant StartTime; @@ -341,9 +335,9 @@ class TimeStepper { static TimeStepper *DefaultTimeStepper; /// All defined time steppers static std::map> AllTimeSteppers; - /// Prescribe modes parsed from config - static PrescribeStateType DefaultPrescribeThicknessMode; - static PrescribeStateType DefaultPrescribeVelocityMode; + /// Prescribe modes parsed from config + static PrescribeStateType DefaultPrescribeThicknessMode; + static PrescribeStateType DefaultPrescribeVelocityMode; }; } // namespace OMEGA From 8d2786081db5da72fc726fc5689400af719cceda Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 8 May 2026 08:32:16 -0600 Subject: [PATCH 12/12] Fixup prescribeThick,Vel calls for FB scheme Co-authored-by: Hyun (Hyun-Gyu) Kang <47987430+hyungyukang@users.noreply.github.com> --- .../omega/src/timeStepping/ForwardBackwardStepper.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index 47550e70debf..3f48b12c809e 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -45,7 +45,7 @@ void ForwardBackwardStepper::doStep( if (AuxState == nullptr) LOG_CRITICAL("Invalid AuxState"); - prescribeVelocity(State, CurLevel, State, CurLevel, SimTime); + prescribeVelocity(State, VelCurLevel, State, VelCurLevel, SimTime); // R_u^{n} = RHS_u(u^{n}, h^{n}, t^{n}) Tend->computeVelocityTendencies(State, AuxState, CurTracerArray, @@ -55,6 +55,11 @@ void ForwardBackwardStepper::doStep( // u^{n+1} = u^{n} + R_u^{n} updateVelocityByTend(State, VelNextLevel, State, VelCurLevel, TimeStep); + updateVelocityByTend(State, VelNextLevel, State, VelCurLevel, TimeStep); + + prescribeVelocity(State, VelNextLevel, State, VelCurLevel, + SimTime + TimeStep); + // R_h^{n} = RHS_h(u^{n+1}, h^{n}, t^{n}) Tend->computeThicknessTendencies(State, AuxState, ThickCurLevel, VelNextLevel, SimTime); @@ -62,7 +67,7 @@ void ForwardBackwardStepper::doStep( // h^{n+1} = h^{n} + R_h^{n} updateThicknessByTend(State, ThickNextLevel, State, ThickCurLevel, TimeStep); - prescribeThickness(State, CurLevel, State, CurLevel); + prescribeThickness(State, ThickNextLevel, State, ThickCurLevel); // R_phi^{n} = RHS_phi(u^{n+1}, h^{n+1}, phi^{n}, t^{n}) Tend->computeTracerTendencies(State, AuxState, CurTracerArray,