diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml index ce0de75afebf..08164bdaa373 100644 --- a/components/omega/configs/Default.yml +++ b/components/omega/configs/Default.yml @@ -66,6 +66,8 @@ Omega: VelocityVertAdvTendencyEnable: true TracerVertAdvTendencyEnable: true PressureGradTendencyEnable: true + VelVertMixTendencyEnable: true + TracerVertMixTendencyEnable: true ManufacturedSolution: WavelengthX: 5.0e6 WavelengthY: 4.33013e6 diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp index ceabc0ede033..c0393e692add 100644 --- a/components/omega/src/ocn/OceanInit.cpp +++ b/components/omega/src/ocn/OceanInit.cpp @@ -31,6 +31,7 @@ #include "Tracers.h" #include "VertAdv.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" @@ -200,6 +201,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) { AuxiliaryState::init(); Eos::init(); PressureGrad::init(); + VertMix::init(); Tendencies::init(); // Validate SurfaceTracerRestoring configuration diff --git a/components/omega/src/ocn/OceanState.h b/components/omega/src/ocn/OceanState.h index 3632449ccc32..59f6845c3305 100644 --- a/components/omega/src/ocn/OceanState.h +++ b/components/omega/src/ocn/OceanState.h @@ -48,11 +48,6 @@ class OceanState { OceanState(const OceanState &) = delete; OceanState(OceanState &&) = delete; - // Current time index - // this index is circular so that it returns to index 0 - // if it is over max index - I4 CurTimeIndex; ///< Time dimension array index for current level - /// Get the current time level index associated with a time level I4 getTimeIndex(const I4 TimeLevel) const; @@ -63,6 +58,11 @@ class OceanState { std::string Name; + // Current time index + // this index is circular so that it returns to index 0 + // if it is over max index + I4 CurTimeIndex; ///< Time dimension array index for current level + // Sizes and global IDs // Note that all sizes are actual counts (1-based) so that loop extents // should always use the 0:NCellsXX-1 form. diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp index ba0df5dd692a..52a4c4c33eca 100644 --- a/components/omega/src/ocn/Tendencies.cpp +++ b/components/omega/src/ocn/Tendencies.cpp @@ -20,6 +20,7 @@ #include "TimeStepper.h" #include "Tracers.h" #include "VertAdv.h" +#include "VertMix.h" #include namespace OMEGA { @@ -39,6 +40,7 @@ void Tendencies::init() { TimeStepper *DefTimeStepper = TimeStepper::getDefault(); Eos *DefEos = Eos::getInstance(); PressureGrad *DefPGrad = PressureGrad::getDefault(); + VertMix *DefVertMix = VertMix::getInstance(); I4 NTracers = Tracers::getNumTracers(); @@ -78,9 +80,10 @@ void Tendencies::init() { TimeInterval TimeStep = DefTimeStepper->getTimeStep(); // Ceate default tendencies - Tendencies::DefaultTendencies = create( - "Default", DefHorzMesh, DefVertCoord, DefVertAdv, DefPGrad, DefEos, - NTracers, TimeStep, &TendConfig, CustomThickTend, CustomVelTend); + Tendencies::DefaultTendencies = + create("Default", DefHorzMesh, DefVertCoord, DefVertAdv, DefPGrad, + DefEos, DefVertMix, NTracers, TimeStep, &TendConfig, + CustomThickTend, CustomVelTend); DefaultTendencies->readConfig(OmegaConfig); @@ -317,6 +320,30 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config HostArray1DI4(TracerIdsToRestoreVec.data(), TracerIdsToRestoreVec.size())); } + + // Validate VertMix tendency + Err += TendConfig.get("VelVertMixTendencyEnable", + this->VMix->VelVertMixSetup.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: VelVertMixTendencyEnable not found in TendConfig"); + + Err += TendConfig.get("TracerVertMixTendencyEnable", + this->VMix->TracerVertMixSetup.Enabled); + CHECK_ERROR_ABORT( + Err, "Tendencies: TracerVertMixTendencyEnable not found in TendConfig"); + + if (this->VMix->VelVertMixSetup.Enabled || + this->VMix->TracerVertMixSetup.Enabled) { + + if (!this->EqState) { + ABORT_ERROR("Tendencies: Eos must be initialized when" + "vertical mixing tendencies are enabled"); + } + if (!this->VMix) { + ABORT_ERROR("Tendencies: VertMix must be initialized when" + "vertical mixing tendencies are enabled"); + } + } } //------------------------------------------------------------------------------ @@ -380,6 +407,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VertAdv *VAdv, ///< [in] Vertical advection PressureGrad *PGrad, ///< [in] Pressure gradient Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStepIn, ///< [in] Time step Config *Options, ///< [in] Configuration options @@ -393,7 +421,8 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord), TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh), CustomThicknessTend(InCustomThicknessTend), - CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad) { + CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad), + VMix(VMix) { // Tendency arrays PseudoThicknessTend = Array2DReal("PseudoThicknessTend", Mesh->NCellsSize, @@ -418,10 +447,11 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies VertAdv *VAdv, ///< [in] Vertical advection PressureGrad *PGrad, ///< [in] Pressure gradient Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStepIn, ///< [in] Time step Config *Options) ///< [in] Configuration options - : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, NTracersIn, + : Tendencies(Name_, Mesh, VCoord, VAdv, PGrad, EqState, VMix, NTracersIn, TimeStepIn, Options, CustomTendencyType{}, CustomTendencyType{}) {} diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h index 5d26ca041458..ac826f915a85 100644 --- a/components/omega/src/ocn/Tendencies.h +++ b/components/omega/src/ocn/Tendencies.h @@ -41,6 +41,7 @@ #include "TimeMgr.h" #include "VertAdv.h" #include "VertCoord.h" +#include "VertMix.h" #include #include @@ -169,6 +170,7 @@ class Tendencies { VertAdv *VAdv, ///< [in] Vertical advection PressureGrad *PGrad, ///< [in] Pressure gradient Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStep, ///< [in] Time step Config *Options, ///< [in] Configuration options @@ -181,6 +183,7 @@ class Tendencies { VertAdv *VAdv, ///< [in] Vertical advection PressureGrad *PGrad, ///< [in] Pressure gradient Eos *EqState, ///< [in] Equation of state + VertMix *VMix, ///< [in] Vertical mixing int NTracersIn, ///< [in] Number of tracers TimeInterval TimeStep, ///< [in] Time step Config *Options ///< [in] Configuration options @@ -199,6 +202,7 @@ class Tendencies { CustomTendencyType CustomVelocityTend; Eos *EqState; ///< Pointer to equation of state PressureGrad *PGrad; ///< Pointer to pressure gradient + VertMix *VMix; ///< Pointer to vertical mixing I4 NTracers; ///< Number of tracers TimeInterval TimeStep; ///< Time step diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index 815dab5e6814..1e2d002c2a9e 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -12,7 +12,12 @@ #include "VertMix.h" #include "DataTypes.h" +#include "Eos.h" +#include "GlobalConstants.h" #include "HorzMesh.h" +#include "HorzOperators.h" +#include "TimeStepper.h" +#include "TriDiagSolvers.h" namespace OMEGA { @@ -36,6 +41,18 @@ GradRichardsonNum::GradRichardsonNum(const HorzMesh *Mesh, OneTwoOneFilter::OneTwoOneFilter(const VertCoord *VCoord) : MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} +VelVertMixSetupOnEdge::VelVertMixSetupOnEdge(const HorzMesh *Mesh, + const VertCoord *VCoord) + : Enabled(false), LocRhoSw(RhoSw), NVertLayers(VCoord->NVertLayers), + CellsOnEdge(Mesh->CellsOnEdge), EdgeMask(VCoord->EdgeMask), + MinLayerEdgeBot(VCoord->MinLayerEdgeBot), + MaxLayerEdgeTop(VCoord->MaxLayerEdgeTop) {} + +TracerVertMixSetupOnCell::TracerVertMixSetupOnCell(const HorzMesh *Mesh, + const VertCoord *VCoord) + : Enabled(false), LocRhoSw(RhoSw), NVertLayers(VCoord->NVertLayers), + MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {} + /// Constructor for VertMix VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object const HorzMesh *Mesh, ///< [in] Horizontal mesh @@ -43,7 +60,8 @@ VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object ) : ComputeVertMixConv(VCoord), ComputeVertMixShear(VCoord), ComputeGradRichardsonNum(Mesh, VCoord), ComputeOneTwoOneFilter(VCoord), - Name(Name), Mesh(Mesh), VCoord(VCoord) { + Name(Name), Mesh(Mesh), VCoord(VCoord), VelVertMixSetup(Mesh, VCoord), + TracerVertMixSetup(Mesh, VCoord) { VertDiff = Array2DReal("VertDiff", Mesh->NCellsSize, VCoord->NVertLayersP1); VertVisc = Array2DReal("VertVisc", Mesh->NCellsSize, VCoord->NVertLayersP1); GradRichNum = @@ -51,6 +69,10 @@ VertMix::VertMix(const std::string &Name, ///< [in] Name for VertMix object GradRichNumSmoothed = Array2DReal("GradRichNumSmoothed", Mesh->NCellsSize, VCoord->NVertLayersP1); + // TODO: Temporary handling of TangentialVelocity + TangentialVelocity = + Array2DReal("TangentialVelocity", Mesh->NEdgesSize, VCoord->NVertLayers); + defineFields(); } @@ -410,4 +432,273 @@ void VertMix::defineFields() { } // end defineIOFields +// Apply implicit velocity vertical mixing +void VertMix::applyVelVertMixImplicit( + OceanState *State, ///< [in] State variables + const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + int ThickTimeLevel, ///< [in] Time level + int VelTimeLevel ///< [in] Time level +) { + + OMEGA_SCOPE(LocNEdgesAll, Mesh->NEdgesAll); + OMEGA_SCOPE(LocVelVertMixSetup, VelVertMixSetup); + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + + const Array2DReal &NormalVelEdge = State->NormalVelocity[VelTimeLevel]; + const Array2DReal &PseudoThickCell = State->PseudoThickness[ThickTimeLevel]; + + // Compute velocity vertical mixing + if (LocVelVertMixSetup.Enabled) { + Pacer::start("Tend:velocityVertMix", 1); + + Eos *EosInstance = Eos::getInstance(); + VertMix *VertMixInstance = VertMix::getInstance(); + + // Obtain TimeStep + const auto *DefTimeStepper = TimeStepper::getDefault(); + const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); + R8 DT; + TimeStep.get(DT, TimeUnits::Seconds); + + const auto &SpecVol = EosInstance->SpecVol; + const auto &VertVisc = VertMixInstance->VertVisc; + + const int NVertLayers = VCoord->NVertLayers; + const int LocVecLength = VecLength; + auto LConfig = + TriDiagSolver::makeLaunchConfig(Mesh->NEdgesAll, NVertLayers); + + parallelForOuter( + LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { + const int IStart = Team.league_rank() * LocVecLength; + const int ILen = Kokkos::max( + 0, Kokkos::min(LocVecLength, LocNEdgesAll - IStart)); + + TriDiagDiffScratch Scratch(Team, NVertLayers); + + // Construct a tri-diag diffusion matrix and RHS + parallelForInner(Team, NVertLayers, [=](int K) { + for (int IVec = 0; IVec < LocVecLength; ++IVec) { + const int IEdge = IStart + IVec; + + if (IEdge >= LocNEdgesAll) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + + if (K < KMin || K > KMax) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + Real G, H, X; + LocVelVertMixSetup(IEdge, K, KMin, KMax, DT, SpecVol, + PseudoThickCell, VertVisc, NormalVelEdge, + G, H, X); + + Scratch.G(K, IVec) = G; + Scratch.H(K, IVec) = H; + Scratch.X(K, IVec) = X; + } + }); + + // Solve the tri-diag diffusion system + Team.team_barrier(); + TriDiagDiffSolver::solve(Team, Scratch); + Team.team_barrier(); + + // Store the solution vector X + parallelForInner(Team, NVertLayers, [=](int K) { + for (int IVec = 0; IVec < ILen; ++IVec) { + const int IEdge = IStart + IVec; + + if (K >= MinLayerEdgeBot(IEdge) && + K <= MaxLayerEdgeTop(IEdge)) { + NormalVelEdge(IEdge, K) = Scratch.X(K, IVec); + } + } + }); + }); + + Pacer::stop("Tend:velocityVertMix", 1); + } +} // applyVelVertMixImplicit + +// Apply implicit tracer vertical mixing +void VertMix::applyTracerVertMixImplicit( + OceanState *State, ///< [in] State variables + const AuxiliaryState *AuxState, ///< [in] Auxilary state variables + Array3DReal &TracerArray, ///< [in] Tracer array + int NTracers, ///< [in] Number of tracers + int ThickTimeLevel, ///< [in] Time level + int VelTimeLevel ///< [in] Time level +) { + + OMEGA_SCOPE(LocNCellsAll, Mesh->NCellsAll); + OMEGA_SCOPE(LocTracerVertMixSetup, TracerVertMixSetup); + OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); + OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); + + const Array2DReal &PseudoThickCell = State->PseudoThickness[ThickTimeLevel]; + + if (LocTracerVertMixSetup.Enabled) { + Pacer::start("Tend:tracerVertMix", 1); + + Eos *EosInstance = Eos::getInstance(); + VertMix *VertMixInstance = VertMix::getInstance(); + + // Obtain TimeStep + const auto *DefTimeStepper = TimeStepper::getDefault(); + const TimeInterval TimeStep = DefTimeStepper->getTimeStep(); + R8 DT; + TimeStep.get(DT, TimeUnits::Seconds); + + const auto &SpecVol = EosInstance->SpecVol; + const auto &VertDiff = VertMixInstance->VertDiff; + + const int NVertLayers = VCoord->NVertLayers; + auto LConfig = + TriDiagSolver::makeLaunchConfig(Mesh->NCellsAll, NVertLayers); + const int LocVecLength = VecLength; + + for (int L = 0; L < NTracers; ++L) { + parallelForOuter( + LConfig, KOKKOS_LAMBDA(int, const TeamMember &Team) { + const int IStart = Team.league_rank() * LocVecLength; + const int ILen = Kokkos::max( + 0, Kokkos::min(LocVecLength, LocNCellsAll - IStart)); + + TriDiagDiffScratch Scratch(Team, NVertLayers); + + // Construct a tri-diag diffusion matrix and RHS + parallelForInner(Team, NVertLayers, [=](int K) { + for (int IVec = 0; IVec < LocVecLength; ++IVec) { + const int ICell = IStart + IVec; + + if (ICell >= LocNCellsAll) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + const int KMin = MinLayerCell(ICell); + const int KMax = MaxLayerCell(ICell); + + if (K < KMin || K > KMax) { + // Fill values + Scratch.G(K, IVec) = 0._Real; + Scratch.H(K, IVec) = 1._Real; + Scratch.X(K, IVec) = 0._Real; + continue; + } + + Real G, H, X; + LocTracerVertMixSetup(L, ICell, K, KMin, KMax, DT, + SpecVol, PseudoThickCell, VertDiff, + TracerArray, G, H, X); + Scratch.G(K, IVec) = G; + Scratch.H(K, IVec) = H; + Scratch.X(K, IVec) = X; + } + }); + + // Solve the tri-diag diffusion system + Team.team_barrier(); + TriDiagDiffSolver::solve(Team, Scratch); + Team.team_barrier(); + + // Store the solution vector X + parallelForInner(Team, NVertLayers, [=](int K) { + for (int IVec = 0; IVec < ILen; ++IVec) { + const int ICell = IStart + IVec; + + if (K >= MinLayerCell(ICell) && + K <= MaxLayerCell(ICell)) { + TracerArray(L, ICell, K) = Scratch.X(K, IVec); + } + } + }); + }); + + } // for L + + Pacer::stop("Tend:tracerVertMix", 1); + } + +} // applyTracerVertMixImplicit + +/// Apply implicit vertical mixing to velocities and tracers +void VertMix::VertMixImplicit(OceanState *State, AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, + int TimeLevel) { + + // get NormalVelocity + Array2DReal NormalVelEdge = State->getNormalVelocity(TimeLevel); + + // get temperature and salinity + I4 ConservTempIdx; + I4 AbsSalinityIdx; + Tracers::getIndex(ConservTempIdx, "Temperature"); + Tracers::getIndex(AbsSalinityIdx, "Salinity"); + + const auto ConservTemp = + Kokkos::subview(TracerArray, ConservTempIdx, Kokkos::ALL, Kokkos::ALL); + const auto AbsSalinity = + Kokkos::subview(TracerArray, AbsSalinityIdx, Kokkos::ALL, Kokkos::ALL); + + // get an instance of equation of state + Eos *EqState = Eos::getInstance(); + + // TODO: Temporary handling of computation of tangential velocity + // Compute tangential velocity + OMEGA_SCOPE(MinLayerEdgeBot, VCoord->MinLayerEdgeBot); + OMEGA_SCOPE(MaxLayerEdgeTop, VCoord->MaxLayerEdgeTop); + OMEGA_SCOPE(LocTangentialVelocity, TangentialVelocity); + + TangentialReconOnEdge TanReconEdge(Mesh); + + parallelForOuter( + {Mesh->NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int KMin = MinLayerEdgeBot(IEdge); + const int KMax = MaxLayerEdgeTop(IEdge); + const int KRange = vertRangeChunked(KMin, KMax); + parallelForInner( + Team, KRange, INNER_LAMBDA(int KChunk) { + TanReconEdge(LocTangentialVelocity, IEdge, KChunk, + NormalVelEdge); + }); + }); + + // Update Pressure, SpecVol + AuxState->computeMomVertAux(State, TracerArray, TimeLevel, TimeLevel); + + // Compute Brunt-Vaisala frequency squared + EqState->computeBruntVaisalaFreqSq( + ConservTemp, AbsSalinity, VCoord->PressureInterface, EqState->SpecVol); + + // Compute vertical mixing coefficients + computeVertMix(NormalVelEdge, LocTangentialVelocity, + EqState->BruntVaisalaFreqSq); + + // Apply implicit mixing to velocities + applyVelVertMixImplicit(State, AuxState, TimeLevel, TimeLevel); + + // Apply implicit mixing to tracers + applyTracerVertMixImplicit(State, AuxState, TracerArray, NTracers, TimeLevel, + TimeLevel); + +} // VertMixImplicit + } // namespace OMEGA diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index 9466eb67d7ec..7e8312fa06f8 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -228,6 +228,186 @@ class OneTwoOneFilter { Array1DI4 MaxLayerCell; }; +/// Velocity vertical mixing +class VelVertMixSetupOnEdge { + public: + bool Enabled; + Real LocRhoSw; + + VelVertMixSetupOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord); + + KOKKOS_FUNCTION void operator()(I4 IEdge, I4 K, I4 KMin, I4 KMax, Real DT, + const Array2DReal &SpecVol, + const Array2DReal &PseudoThickCell, + const Array2DReal &VertVisc, + const Array2DReal &NormalVelEdge, Real &G, + Real &H, Real &X) const { + + // Implicit vertical mixing is solved using the diffusion-form + // tridiagonal solver, `TriDiagDiffSolver`, provided by `TriDiagSolvers`. + // + // * Governing tridiagonal system: + // + // -G_{k-1} Y_{k-1} + // + (G_{k-1} + G_k + H_k) Y_k + // - G_k Y_{k+1} + // = X_k + // + // G: coefficient contributing to both diagonal and off-diagonal entries + // H: coefficient contributing only to the diagonal entry + // X: RHS vector + // + // * Matrix structure: + // + // [ D G 0 ... 0 ][ ] [ ] : k = 0 + // [ G D G ... 0 ][ ] [ ] + // [ 0 G D G ... 0 ][ ] [ ] + // [ ... ][Y] = [X] + // [ 0 ... G D G 0 ][ ] [ ] + // [ 0 ... 0 G D G ][ ] [ ] + // [ 0 ... 0 0 G D ][ ] [ ] : k = NVertLayers - 1 + // + // where D = G_{k-1} + G_k + H_k + // + // G_k = [DT * VertVisc_k^top / (Rho_0 * SpecVol_k^top)] + // / PseudoThick_k^top + // H_k = PseudoThick_k + // X_k = PseudoThick_k * NormVel_k + + // Fill values + G = 0.0_Real; + H = 1.0_Real; + X = 0.0_Real; + + const I4 JCell0 = CellsOnEdge(IEdge, 0); + const I4 JCell1 = CellsOnEdge(IEdge, 1); + + const Real PseudoThickEdgeK = + 0.5_Real * (PseudoThickCell(JCell0, K) + PseudoThickCell(JCell1, K)); + + H = PseudoThickEdgeK; + + // Unknown is NormVel^{n+1}. + X = PseudoThickEdgeK * NormalVelEdge(IEdge, K); + + if (K < KMax) { + const Real PseudoThickEdgeKp1 = + 0.5_Real * + (PseudoThickCell(JCell0, K + 1) + PseudoThickCell(JCell1, K + 1)); + + const Real PseudoThickEdgeBot = + 0.5_Real * (PseudoThickEdgeK + PseudoThickEdgeKp1); + + // Interpolation from cell center to top using + // the two-point linear interpolation + const Real SpecVolEdgeBot = + (0.5_Real * (SpecVol(JCell0, K) + SpecVol(JCell1, K)) * + PseudoThickEdgeKp1 + + 0.5_Real * (SpecVol(JCell0, K + 1) + SpecVol(JCell1, K + 1)) * + PseudoThickEdgeK) / + (PseudoThickEdgeK + PseudoThickEdgeKp1); + + const Real ViscAlphaEdgeBot = + 0.5_Real * (VertVisc(JCell0, K + 1) + VertVisc(JCell1, K + 1)) / + (LocRhoSw * SpecVolEdgeBot); + + G = DT * ViscAlphaEdgeBot / PseudoThickEdgeBot; + } + } + + private: + I4 NVertLayers; + Array2DI4 CellsOnEdge; + Array2DReal EdgeMask; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLayerEdgeTop; +}; + +// Tracer vertical mixing term +class TracerVertMixSetupOnCell { + public: + bool Enabled; + Real LocRhoSw; + + TracerVertMixSetupOnCell(const HorzMesh *Mesh, const VertCoord *VCoord); + + KOKKOS_FUNCTION void operator()(I4 L, I4 ICell, I4 K, I4 KMin, I4 KMax, + Real DT, const Array2DReal &SpecVol, + const Array2DReal &PseudoThickCell, + const Array2DReal &VertDiff, + const Array3DReal &TracersOnCell, Real &G, + Real &H, Real &X) const { + + // Implicit vertical mixing is solved using the diffusion-form + // tridiagonal solver, `TriDiagDiffSolver`, provided by `TriDiagSolvers`. + // + // * Governing tridiagonal system: + // + // -G_{k-1} Y_{k-1} + // + (G_{k-1} + G_k + H_k) Y_k + // - G_k Y_{k+1} + // = X_k + // + // G: coefficient contributing to both diagonal and off-diagonal entries + // H: coefficient contributing only to the diagonal entry + // X: RHS vector + // + // * Matrix structure: + // + // [ D G 0 ... 0 ][ ] [ ] : k = 0 + // [ G D G ... 0 ][ ] [ ] + // [ 0 G D G ... 0 ][ ] [ ] + // [ ... ][Y] = [X] + // [ 0 ... G D G 0 ][ ] [ ] + // [ 0 ... 0 G D G ][ ] [ ] + // [ 0 ... 0 0 G D ][ ] [ ] : k = NVertLayers - 1 + // + // where D = G_{k-1} + G_k + H_k + // + // G_k = [DT * VertDiff_k^top / (Rho_0 * SpecVol_k^top)] + // / PseudoThick_k^top + // H_k = PseudoThick_k + // X_k = PseudoThick_k * Phi_k + + // Fill values + G = 0.0_Real; + H = 1.0_Real; + X = 0.0_Real; + + const Real PseudoThickCellK = PseudoThickCell(ICell, K); + + H = PseudoThickCellK; + + // Unknown is Phi^{n+1}. + X = PseudoThickCellK * TracersOnCell(L, ICell, K); + + if (K < KMax) { + + const Real PseudoThickCellKp1 = PseudoThickCell(ICell, K + 1); + + const Real PseudoThickCellBot = + 0.5_Real * (PseudoThickCellK + PseudoThickCellKp1); + + // Interpolation from cell center to top using + // the two-point linear interpolation + const Real SpecVolCellBot = + (SpecVol(ICell, K) * PseudoThickCellKp1 + + SpecVol(ICell, K + 1) * PseudoThickCellK) / + (PseudoThickCellK + PseudoThickCellKp1); + + const Real DiffAlphaCellBot = + VertDiff(ICell, K + 1) / (LocRhoSw * SpecVolCellBot); + + G = DT * DiffAlphaCellBot / PseudoThickCellBot; + } + } + + private: + I4 NVertLayers; + Array1DI4 MinLayerCell; + Array1DI4 MaxLayerCell; +}; + /// Class for Vertical Mixing Coefficient (VertMix) calculations class VertMix { public: @@ -243,6 +423,9 @@ class VertMix { Array2DReal GradRichNumSmoothed; ///< Smoothed Gradient Richardson number field + // TODO: Temporary handling of TangentialVelocity + Array2DReal TangentialVelocity; ///< Tangential velocity + std::string VertDiffFldName; ///< Field name for vertical diffusivity std::string VertViscFldName; ///< Field name for vertical viscosity std::string @@ -264,6 +447,9 @@ class VertMix { ///< calculation OneTwoOneFilter ComputeOneTwoOneFilter; ///< Functor for 1-2-1 filtering + VelVertMixSetupOnEdge VelVertMixSetup; + TracerVertMixSetupOnCell TracerVertMixSetup; + /// Compute vertical diffusivity and viscosity for all cells/layers void computeVertMix(const Array2DReal &NormalVelocity, const Array2DReal &TangentialVelocity, @@ -272,6 +458,18 @@ class VertMix { /// Initialize VertMix from config and mesh static void init(); + void applyVelVertMixImplicit(OceanState *State, + const AuxiliaryState *AuxState, + int ThickTimeLevel, int VelTimeLevel); + void applyTracerVertMixImplicit(OceanState *State, + const AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, + int ThickTimeLevel, int VelTimeLevel); + + /// Apply implicit vertical mixing to velocities and tracers + void VertMixImplicit(OceanState *State, AuxiliaryState *AuxState, + Array3DReal &TracerArray, int NTracers, int TimeLevel); + private: /// Private constructor VertMix(const std::string &Name, const HorzMesh *Mesh, diff --git a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp index 11ddc664d1eb..bcb0d9463835 100644 --- a/components/omega/src/timeStepping/ForwardBackwardStepper.cpp +++ b/components/omega/src/timeStepping/ForwardBackwardStepper.cpp @@ -6,6 +6,7 @@ #include "ForwardBackwardStepper.h" #include "Pacer.h" +#include "VertMix.h" namespace OMEGA { @@ -37,9 +38,13 @@ void ForwardBackwardStepper::doStep( const int ThickNextLevel = 1; const int TracerNextLevel = 1; + int NTracers = Tracers::getNumTracers(); + Array3DReal CurTracerArray = Tracers::getAll(TracerCurLevel); Array3DReal NextTracerArray = Tracers::getAll(TracerNextLevel); + VertMix *VMix = VertMix::getInstance(); + if (State == nullptr) LOG_CRITICAL("Invalid State"); if (AuxState == nullptr) @@ -86,6 +91,13 @@ void ForwardBackwardStepper::doStep( Tracers::updateTimeLevels(); Pacer::stop("ForwardBackward:haloExch", 3); + // Apply implicit vertical mixing + CurTracerArray = Tracers::getAll(VelCurLevel); + if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { + VMix->VertMixImplicit(State, AuxState, CurTracerArray, NTracers, + State->CurTimeIndex); + } + validateOceanState(State, AuxState, VertCoord::getDefault(), 0); // Advance the clock and update the simulation time diff --git a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp index 372415859102..16fc8fa36c0d 100644 --- a/components/omega/src/timeStepping/RungeKutta2Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta2Stepper.cpp @@ -6,6 +6,7 @@ #include "RungeKutta2Stepper.h" #include "Pacer.h" +#include "VertMix.h" namespace OMEGA { @@ -31,9 +32,13 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state const int CurLevel = 0; const int NextLevel = 1; + int NTracers = Tracers::getNumTracers(); + Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); + VertMix *VMix = VertMix::getInstance(); + prescribeState(State, CurLevel, State, CurLevel, SimTime); // q = (h,u,phi) @@ -70,6 +75,13 @@ void RungeKutta2Stepper::doStep(OceanState *State, // model state Tracers::updateTimeLevels(); Pacer::stop("RK2:haloExch", 3); + // Apply implicit vertical mixing + CurTracerArray = Tracers::getAll(CurLevel); + if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { + VMix->VertMixImplicit(State, AuxState, CurTracerArray, NTracers, + State->CurTimeIndex); + } + validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); // Advance the clock and update the simulation time diff --git a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp index 5ad5386e28f9..3b21f972854e 100644 --- a/components/omega/src/timeStepping/RungeKutta4Stepper.cpp +++ b/components/omega/src/timeStepping/RungeKutta4Stepper.cpp @@ -6,6 +6,7 @@ #include "RungeKutta4Stepper.h" #include "Pacer.h" +#include "VertMix.h" namespace OMEGA { @@ -78,11 +79,14 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state const int CurLevel = 0; const int NextLevel = 1; + int NTracers = Tracers::getNumTracers(); Array3DReal CurTracerArray = Tracers::getAll(CurLevel); Array3DReal NextTracerArray = Tracers::getAll(NextLevel); TimeInstant ForcingStageTime = SimTime; + VertMix *VMix = VertMix::getInstance(); + for (int Stage = 0; Stage < NStages; ++Stage) { const TimeInstant StageTime = SimTime + RKC[Stage] * TimeStep; // first stage does: @@ -135,6 +139,13 @@ void RungeKutta4Stepper::doStep(OceanState *State, // model state Tracers::updateTimeLevels(); Pacer::stop("RK4:haloExch", 3); + // Apply implicit vertical mixing + CurTracerArray = Tracers::getAll(CurLevel); + if (VMix->VelVertMixSetup.Enabled or VMix->TracerVertMixSetup.Enabled) { + VMix->VertMixImplicit(State, AuxState, CurTracerArray, NTracers, + State->CurTimeIndex); + } + validateOceanState(State, AuxState, VertCoord::getDefault(), CurLevel); // Advance the clock and update the simulation time diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp index 1678a08470f6..b3126ad21a64 100644 --- a/components/omega/test/infra/IOStreamTest.cpp +++ b/components/omega/test/infra/IOStreamTest.cpp @@ -15,6 +15,7 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Eos.h" #include "Error.h" #include "Field.h" #include "Forcing.h" @@ -31,6 +32,7 @@ #include "TimeStepper.h" #include "Tracers.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include #include @@ -108,6 +110,12 @@ void initIOStreamTest(Clock *&ModelClock // Model clock PressureGrad::init(); + // Initialize Eos + Eos::init(); + + // Initialize VertMix + VertMix::init(); + // Intialize Tendencies Tendencies::init(); @@ -229,6 +237,7 @@ int main(int argc, char **argv) { // Clean up environments TimeStepper::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Tendencies::clear(); Tracers::clear(); OceanState::clear(); diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp index 252056ef2074..09caff2e6d73 100644 --- a/components/omega/test/ocn/StateTest.cpp +++ b/components/omega/test/ocn/StateTest.cpp @@ -13,6 +13,7 @@ #include "DataTypes.h" #include "Decomp.h" #include "Dimension.h" +#include "Eos.h" #include "Error.h" #include "Field.h" #include "Halo.h" @@ -28,6 +29,7 @@ #include "Tendencies.h" #include "TimeStepper.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include @@ -95,6 +97,12 @@ void initStateTest() { // Initialize pressure gradient PressureGrad::init(); + // Initialize equation of state + Eos::init(); + + // Initialize vertical mixing + VertMix::init(); + // Create tendencies Tendencies::init(); @@ -408,6 +416,7 @@ int main(int argc, char *argv[]) { Tracers::clear(); AuxiliaryState::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Tendencies::clear(); TimeStepper::clear(); HorzMesh::clear(); diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp index 5c8bb8d3d3f0..5268d0a436ea 100644 --- a/components/omega/test/ocn/TendenciesTest.cpp +++ b/components/omega/test/ocn/TendenciesTest.cpp @@ -22,6 +22,7 @@ #include "Pacer.h" #include "TimeStepper.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include @@ -142,6 +143,7 @@ int initTendenciesTest(const std::string &mesh) { PressureGrad::init(); Eos::init(); Forcing::init(); + VertMix::init(); int StateErr = OceanState::init(); if (StateErr != 0) { @@ -177,6 +179,7 @@ int testTendencies() { const auto VAdv = VertAdv::getDefault(); const auto PGrad = PressureGrad::getDefault(); const auto EqState = Eos::getInstance(); + const auto VMix = VertMix::getInstance(); VCoord->NVertLayers = 12; // test creation of another tendencies @@ -188,7 +191,7 @@ int testTendencies() { int NTracersTest = 3; Tendencies::create("TestTendencies", Mesh, VCoord, VAdv, PGrad, EqState, - NTracersTest, ZeroTimeStep, &TendConfig); + VMix, NTracersTest, ZeroTimeStep, &TendConfig); // test retrievel of another tendencies if (Tendencies::get("TestTendencies")) { @@ -340,6 +343,7 @@ void finalizeTendenciesTest() { Forcing::clear(); Tracers::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Eos::destroyInstance(); AuxiliaryState::clear(); OceanState::clear(); diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp index 81f6bcef5938..ec458273c351 100644 --- a/components/omega/test/timeStepping/TimeStepperTest.cpp +++ b/components/omega/test/timeStepping/TimeStepperTest.cpp @@ -37,6 +37,7 @@ #include "TimeMgr.h" #include "Tracers.h" #include "VertCoord.h" +#include "VertMix.h" #include "mpi.h" #include @@ -187,6 +188,7 @@ int initTimeStepperTest(const std::string &mesh) { AuxiliaryState::init(); Eos::init(); PressureGrad::init(); + VertMix::init(); Tendencies::init(); // finish initializing default time stepper @@ -205,6 +207,7 @@ int initTimeStepperTest(const std::string &mesh) { auto *DefHalo = Halo::getDefault(); auto *DefEos = Eos::getInstance(); auto *DefPGrad = PressureGrad::getDefault(); + auto *DefVMix = VertMix::getInstance(); int NTracers = Tracers::getNumTracers(); const int NTimeLevels = 2; @@ -233,8 +236,8 @@ int initTimeStepperTest(const std::string &mesh) { // Creating non-default tendencies with custom velocity tendencies auto *TestTendencies = Tendencies::create( "TestTendencies", DefMesh, DefVertCoord, DefVAdv, DefPGrad, DefEos, - NTracers, ZeroTimeStep, &Options, Tendencies::CustomTendencyType{}, - DecayVelocityTendency{}); + DefVMix, NTracers, ZeroTimeStep, &Options, + Tendencies::CustomTendencyType{}, DecayVelocityTendency{}); if (!TestTendencies) { Err++; LOG_ERROR("TimeStepperTest: error creating test tendencies"); @@ -256,6 +259,8 @@ int initTimeStepperTest(const std::string &mesh) { DefVAdv->ThickVertAdvEnabled = false; DefVAdv->VelVertAdvEnabled = false; DefVAdv->TracerVertAdvEnabled = false; + DefVMix->VelVertMixSetup.Enabled = false; + DefVMix->TracerVertMixSetup.Enabled = false; return Err; } @@ -295,6 +300,7 @@ void finalizeTimeStepperTest() { Tracers::clear(); TimeStepper::clear(); PressureGrad::clear(); + VertMix::destroyInstance(); Eos::destroyInstance(); Tendencies::clear(); AuxiliaryState::clear();