Add centered pressure gradient implimentation

components/omega/configs/Default.yml

@@ -23,13 +23,19 @@ Omega:
   State:
     NTimeLevels: 2
   Advection:
+    Coef3rdOrder: 0.25
     FluxThicknessType: Center
+    HorzTracerFluxOrder: 2
     FluxTracerType: Center
+    VerticalTracerFluxLimiterEnable: true
+    VerticalTracerFluxOrder: 3
   WindStress:
     InterpType: Isotropic
   VertCoord:
     Density0: 1026.0
     MovementWeightType: Uniform
+  PressureGrad:
+    PressureGradType: Centered
   Tendencies:
     ThicknessFluxTendencyEnable: true
     PVTendencyEnable: true
@@ -50,6 +56,10 @@ Omega:
     EddyDiff4: 0.0
     UseCustomTendency: false
     ManufacturedSolutionTendency: false
+    ThicknessVertAdvTendencyEnable: true
+    VelocityVertAdvTendencyEnable: true
+    TracerVertAdvTendencyEnable: true
+    PressureGradTendencyEnable: true
   Tracers:
     Base: [Temperature, Salinity]
     Debug: [Debug1, Debug2, Debug3]

components/omega/doc/design/Tendencies.md

@@ -39,14 +39,23 @@ retrieval.
 ```c++
 class Tendencies{
  public:
+   // Arrays for accumulating tendencies
    Array2DReal LayerThicknessTend;
    Array2DReal NormalVelocityTend;
+   Array3DReal TracerTend;
+   // Instances of tendency terms
    ThicknessFluxDivOnCell ThicknessFluxDiv;
    PotentialVortHAdvOnEdge PotientialVortHAdv;
    KEGradOnEdge KEGrad;
    SSHGradOnEdge SSHGrad;
    VelocityDiffusionOnEdge VelocityDiffusion;
    VelocityHyperDiffOnEdge VelocityHyperDiff;
+   WindForcingOnEdge WindForcing;
+   BottomDragOnEdge BottomDrag;
+   TracerDiffOnCell TracerDiffusion;
+   TracerHyperDiffOnCell TracerHyperDiff;
+   TracerHorzAdvOnCell TracerHorzAdv;
+   TracerHighOrderHorzAdvOnCell TracerHighOrderHorzAdv;
  private:
    static Tendencies *DefaultTendencies;
    static std::map<std::string, std::unique_ptr<Tendencies>> AllTendencies;
@@ -79,29 +88,29 @@ Tendencies* Tendencies::create(const std::string &Name, const HorzMesh *Mesh, in
 #### 4.2.2 Initialization
 The init method will create the default tendencies and return an error code:
 ```c++
-int Tendencies::init();
+static void Tendencies::init();
 ```
 
 #### 4.2.3 Retrieval
 There will be methods for getting the default and non-default tendencies instances:
 ```c++
-Tendencies *Tendencies::getDefault();
-Tendencies *Tendencies::get(const std::string &Name);
+static Tendencies *Tendencies::getDefault();
+static Tendencies *Tendencies::get(const std::string &Name);
 ```
 
 #### 4.2.4 Computation
 The 'computeAll' method will compute and accumulate all layer thickness and normal velocity tendencies using the ocean
 state and auxiliary state from a given time level:
 ```c++
-void Tendencies::computeAllTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel);
+void Tendencies::computeAllTendencies(const OceanState *State, const AuxilaryState *AuxState, const Array3DReal &TracerArray, int TimeLevel, int VelTimeLevel, TimeInstant Time);
 ```
 The layer thickness tendencies will be computed with a method:
 ```c++
-void Tendencies::computeThicknessTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel);
+void Tendencies::computeThicknessTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel, int VelTimeLevel, TimeInstant Time);
 ```
 The normal velocity tendencies will be computed with a method:
 ```c++
-void Tendencies::computeVelocityTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel);
+void Tendencies::computeVelocityTendencies(const OceanState *State, const AuxilaryState *AuxState, int TimeLevel, int VelTimeLevel, TimeInstant Time);
 ```
 
 
@@ -111,8 +120,8 @@ when they are no longer in scope. The erase method
 will remove a named tendencies instance, whereas the clear method will remove all of
 them. Both will call the destructor in the process.
 ```c++
-void Tendencies::erase(const std::string &Name);
-void Tendencies::clear();
+static void Tendencies::erase(const std::string &Name);
+static void Tendencies::clear();
 ```
 
 ## 5 Verification and Testing

components/omega/doc/design/Tendency.md

@@ -39,24 +39,25 @@ The class contains private member variables for any constant data that are defin
 Each tendency term will have a `bool` variable which can be set to enable/disable the computation of the tendency.
 
 ```c++
-class ThicknessFluxDivergenceOnCell {
+class ThicknessFluxDivOnCell {
   public:
 
     bool enabled = false;
 
-    ThicknessFluxDivergenceOnCell(const HorzMesh *Mesh, Config *Options);
+    ThicknessFluxDivOnCell(const HorzMesh *Mesh, Config *Options);
 
     KOKKOS_FUNCTION void operator()(Array2DReal &Tend
                                     int ICell,
                                     int KChunk,
-                                    const Array2DReal &ThicknessFlux);
+                                    const Array2DReal &ThicknessFlux,
+                                    const Array2DReal &NormalVelEdge);
 
   private:
     Array1DI4 NEdgesOnCell;
     Array2DI4 EdgesOnCell;
-    Array1DR8 DvEdge;
-    Array1DR8 AreaCell;
-    Array2DR8 EdgeSignOnCell;
+    Array1DReal DvEdge;
+    Array1DReal AreaCell;
+    Array2DReal EdgeSignOnCell;
 };
 ```
 
@@ -66,31 +67,38 @@ class ThicknessFluxDivergenceOnCell {
 Tendency functor constructors are responsible for initializing the private member variables:
 
 ```c++
-ThicknessFluxDivergenceOnCell::ThicknessFluxDivergenceOnCell(HorzMesh const *Mesh, Config *Options)
+ThicknessFluxDivOnCell::ThicknessFluxDivOnCell(HorzMesh const *Mesh)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell),
-      EdgeSignOnCell(Mesh->EdgeSignOnCell) {
-
-    enabled = Options->get("ThicknessFluxTendencyEnable")
-}
+      EdgeSignOnCell(Mesh->EdgeSignOnCell) { }
 ```
 
 #### 4.2.2 operator
 The operator method implements the tendency computation for a chunk of vertical layers at a given horizontal mesh location.
 The inner loop over a chunk of vertical layers enables CPU vectorization.
 
 ```c++
-KOKKOS_FUNCTION void ThicknessFluxDivergenceOnCell::operator()(Array2DReal &Tend
-                                                               int ICell,
-                                                               int KChunk,
-                                                               const Array2DReal &ThicknessFlux)  const {
-   const Real InvAreaCell = 1. / AreaCell(iCell);
-   for (int J = 0; J < NEdgesOnCell(ICell); ++J) {
-      const int JEdge = EdgesOnCell(ICell, J);
-      for (int K = KChunk * VecLength; K < (KChunk + 1) * VecLength; ++K) {
-         Tend(ICell, K) -= DvEdge(JEdge) * EdgeSignOnCell(ICell, J) * ThicknessFlux(JEdge, K) * InvAreaCell;
-      }
-   }
+KOKKOS_FUNCTION void ThicknessFluxDivOnCell::operator()(Array2DReal &Tend
+                                                        int ICell,
+                                                        int KChunk,
+                                                        const Array2DReal &ThicknessFlux,
+                                                        const Array2DReal &NormalVelEdge)  const {
+    const I4 KStart        = KChunk * VecLength;
+    const Real InvAreaCell = 1._Real / AreaCell(ICell);
+    Real DivTmp[VecLength] = {0};
+    for (int J = 0; J < NEdgesOnCell(ICell); ++J) {
+       const I4 JEdge = EdgesOnCell(ICell, J);
+       for (int KVec = 0; KVec < VecLength; ++KVec) {
+          const I4 K = KStart + KVec;
+          DivTmp[KVec] -= DvEdge(JEdge) * EdgeSignOnCell(ICell, J) *
+                          ThicknessFlux(JEdge, K) * NormalVelEdge(JEdge, K) *
+                          InvAreaCell;
+       }
+    }
+    for (int KVec = 0; KVec < VecLength; ++KVec) {
+       const I4 K = KStart + KVec;
+       Tend(ICell, K) -= DivTmp[KVec];
+    }
 }
 ```
 

components/omega/doc/devGuide/OceanState.md

@@ -40,38 +40,38 @@ For now, member variables that are host arrays have variable names are appended
 `H`.  Array variable names not ending in `H` are device arrays.  For a given time level,
 host to device array is performed via:
 ```c++
-Err = State->copyToDevice(TimeLevel);
+State->copyToDevice(TimeLevel);
 ```
 and a copy from device to host is performed by:
 ```c++
-Err = State->copyToHost(TimeLevel);
+State->copyToHost(TimeLevel);
 ```
+These functions validate the `TimeLevel` via `OMEGA_REQUIRE` and will abort if invalid.
 Eventually, the host arrays will be eliminated when `IO` and `Halo` are extended to
 handle host <-> device transfers.
 
 A time level update to advance the solution at the end of a model timestep is done by:
 ```c++
-Err = State->updateTimeLevels();
+State->updateTimeLevels();
 ```
 This shifts the time level indices within the `State` instance. A halo exchange is
 also performed on these arrays and the IOFields data pointer is attached
 to the current time level.
 
 The arrays associated with a given time level can be accessed with the functions:
 ```c++
-Array2DReal LayerThick;
-Err = State->getLayerThickness(LayerThick, TimeLevel);
-Array2DReal NormVel;
-Err = State->getNormalVelocity(NormVel, TimeLevel);
+Array2DReal LayerThick = State->getLayerThickness(TimeLevel);
+Array2DReal NormVel = State->getNormalVelocity(TimeLevel);
 ```
 for the device arrays and
 ```c++
-HostArray2DReal LayerThickH;
-Err = State->getLayerThicknessH(LayerThickH, TimeLevel);
-HostArray2DReal NormVelH;
-Err = State->getNormalVelocityH(NormVelH, TimeLevel);
+HostArray2DReal LayerThickH = State->getLayerThicknessH(TimeLevel);
+HostArray2DReal NormVelH = State->getNormalVelocityH(TimeLevel);
 ```
-for the host arrays. The time level convention is:
+for the host arrays. These functions return the arrays directly and will abort
+via `OMEGA_REQUIRE` if an invalid `TimeLevel` is provided.
+
+The time level convention is:
 | time level | `TimeLevel` |
 |------------|-------------|
 | New | 1 |