Merge pull request #410 from hyungyukang/omega/vmix-hookup

This PR adds the velocity and tracer vertical mixing tendency terms. The implicit vertical mixing is solved by using TriDiagDiffSolver and is applied once at the end of each time step.

    - Implemented implicit vertical mixing using TriDiagDiffSolver.
    - Added VertMixImplicit to VertMix.
    - Applied vertical mixing once at the end of each time step.
    - Note: Tangential velocity, which is required by ComputeGradRichardsonNum in VertMix, is temporarily handled in VertMix.

Author: sbrus89

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

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

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.

components/omega/src/ocn/Tendencies.cpp

@@ -20,6 +20,7 @@
 #include "TimeStepper.h"
 #include "Tracers.h"
 #include "VertAdv.h"
+#include "VertMix.h"
 #include <string>
 
 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{}) {}
 

components/omega/src/ocn/Tendencies.h

@@ -41,6 +41,7 @@
 #include "TimeMgr.h"
 #include "VertAdv.h"
 #include "VertCoord.h"
+#include "VertMix.h"
 
 #include <functional>
 #include <memory>
@@ -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