From d969e535a4eda5cf0ceb1d603fbbe2cffd015ffd Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 12 Mar 2026 12:37:07 -0700 Subject: [PATCH 01/11] 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 eac5eeb2876b78cab46c2fb4750936dab382fe9e Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 12 Mar 2026 12:46:17 -0700 Subject: [PATCH 02/11] 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 e72719a3b2f4..d9a5fc223b27 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -62,6 +62,9 @@ void ForwardBackwardStepper::doStep( // u^{n+1} = u^{n} + R_u^{n+1} updateVelocityByTend(State, NextLevel, State, CurLevel, 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 04b7c89edbcb73419da4693ea4f6c473cdbea5f6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 12 Mar 2026 13:45:49 -0700 Subject: [PATCH 03/11] 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 aa977bcb67a6c8a86eac0828b1bf1232352d6365 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 15:59:55 -0700 Subject: [PATCH 04/11] 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 d9a5fc223b27..ae374fdd8c28 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -63,7 +63,8 @@ void ForwardBackwardStepper::doStep( updateVelocityByTend(State, NextLevel, State, CurLevel, 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 b5bbab2a156105d81e2bbd6bbd8883d3634d5573 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 16:01:30 -0700 Subject: [PATCH 05/11] 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 638acccf7e8d786c8db50c1869cc4bda7a8fec3b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 13 Mar 2026 16:44:43 -0700 Subject: [PATCH 06/11] 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 4ae3775007458f2c802ae2ef6f34a58c224bd2e4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 19 Mar 2026 09:01:02 -0700 Subject: [PATCH 07/11] 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 f038470fd2410852b32e20c4c941834cd7e22b21 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 8 Apr 2026 09:27:44 -0700 Subject: [PATCH 08/11] 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 93dc58e424d6b868e081f9d4078d3375a0237305 Mon Sep 17 00:00:00 2001 From: James Overfelt Date: Sat, 25 Apr 2026 10:20:10 -0700 Subject: [PATCH 09/11] 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 b7c670d43abe4026eaf113481f7a2e9163707408 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 29 Apr 2026 13:13:52 -0400 Subject: [PATCH 10/11] 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 ae374fdd8c28..5401679f23e3 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -40,6 +40,8 @@ void ForwardBackwardStepper::doStep( if (AuxState == nullptr) LOG_CRITICAL("Invalid AuxState"); + prescribeVelocity(State, CurLevel, State, CurLevel,SimTime); + // R_h^{n} = RHS_h(u^{n}, h^{n}, t^{n}) Tend->computeThicknessTendencies(State, AuxState, CurLevel, CurLevel, SimTime); @@ -47,6 +49,8 @@ void ForwardBackwardStepper::doStep( // h^{n+1} = h^{n} + R_h^{n} updateThicknessByTend(State, NextLevel, State, CurLevel, TimeStep); + prescribeThickness(State, CurLevel, State, CurLevel); + // R_phi^{n} = RHS_phi(u^{n}, h^{n}, phi^{n}, t^{n}) Tend->computeTracerTendencies(State, AuxState, CurTracerArray, CurLevel, CurLevel, SimTime); @@ -62,10 +66,6 @@ void ForwardBackwardStepper::doStep( // u^{n+1} = u^{n} + R_u^{n+1} updateVelocityByTend(State, NextLevel, State, CurLevel, 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 0575643dc1036da541d73b23f2aa985c27a48403 Mon Sep 17 00:00:00 2001 From: Hyun Kang Date: Wed, 29 Apr 2026 13:22:24 -0400 Subject: [PATCH 11/11] 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 5401679f23e3..09d237876254 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -40,7 +40,7 @@ void ForwardBackwardStepper::doStep( if (AuxState == nullptr) LOG_CRITICAL("Invalid AuxState"); - prescribeVelocity(State, CurLevel, State, CurLevel,SimTime); + prescribeVelocity(State, CurLevel, State, CurLevel, SimTime); // R_h^{n} = RHS_h(u^{n}, h^{n}, t^{n}) Tend->computeThicknessTendencies(State, AuxState, CurLevel, CurLevel, 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