Higher order horizontal advection copied from mpas-ocean to Omega.

components/omega/configs/Default.yml

@@ -44,6 +44,7 @@ Omega:
     BottomDragTendencyEnable: false
     BottomDragCoeff: 0.0
     TracerHorzAdvTendencyEnable: true
+    TracerHorzAdvTendencyOrder: 1
     TracerDiffTendencyEnable: true
     EddyDiff2: 10.0
     TracerHyperDiffTendencyEnable: true

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/TendencyTerms.md

@@ -35,8 +35,9 @@ implemented:
 - `SSHGradOnEdge`
 - `VelocityDiffusionOnEdge`
 - `VelocityHyperDiffOnEdge`
+- `WindForcingOnEdge`
+- `BottomDragOnEdge`
 - `TracerHorzAdvOnCell`
+- `TracerHighOrderHorzAdvOnCell`
 - `TracerDiffOnCell`
 - `TracerHyperDiffOnCell`
-- `WindForcingOnEdge`
-- `BottomDragOnEdge`

components/omega/doc/userGuide/TendencyTerms.md

@@ -15,6 +15,7 @@ tendency terms are currently implemented:
 | VelocityDiffusionOnEdge | Laplacian horizontal mixing, defined on edges
 | VelocityHyperDiffOnEdge | biharmonic horizontal mixing, defined on edges
 | TracerHorzAdvOnCell | horizontal advection of thickness-weighted tracers
+| TracerHighOrderHorzAdvOnCell | second order horizontal advection of thickness-weighted tracers
 | TracerDiffOnCell | horizontal diffusion of thickness-weighted tracers
 | TracerHyperDiffOnCell | biharmonic horizontal mixing of thickness-weighted tracers
 | WindForcingOnEdge | forcing by wind stress, defined on edges
@@ -43,10 +44,95 @@ the currently available tendency terms:
 | | ViscDel4 | horizontal biharmonic mixing coefficient for normal velocity
 | | DivFactor | scale factor for the divergence term
 | TracerHorzAdvOnCell | TracerHorzAdvTendencyEnable | enable/disable term
+| | TracerHorzAdvTendencyOrder | 1 for standard linear advection
+| TracerHighOrderHorzAdvOnCell | TracerHorzAdvTendencyEnable | enable/disable term
+| | TracerHorzAdvTendencyOrder | 2 for second order advection algorithm
 | TracerDiffOnCell | TracerDiffTendencyEnable | enable/disable term
 | | EddyDiff2 | horizontal diffusion coefficient
 | TracerHyperDiffOnCell | TracerHyperDiffTendencyEnable | enable/disable term
 | | EddyDiff4 | biharmonic horizontal mixing coeffienct for tracers
 | WindForcingOnEdge | WindForcingTendencyEnable | enable/disable term
 | BottomDragOnEdge | BottomDragTendencyEnable | enable/disable term
 | | BottomDragCoeff | bottom drag coefficient
+
+
+## Second Order Horizontal Advection Algorithm
+
+The horizontal advection is done independently within each ocean layer
+so a two dimensional scheme is required to advect mixing ratios.
+The second order horizontal advection scheme is described in the paper
+
+William C. Skamarock and Almut Gassmann,
+"Conservative Transport Schemes for Spherical Geodesic Grids: High-Order Flux Operators for ODE-Based Time Integration"
+Monthly Weather Review, Vol 139: Issue 9, pp 2962–2975, 2011. doi.org/10.1175/MWR-D-10-05056.1
+
+and only a summary of a few equations from that paper are reproduced here.
+
+The conservative form of the scalar transport equation
+within a fluid is given in equation (1) of Skamarock and Gassmann,
+$$
+ \frac{\partial(\rho\psi)}{\partial t} = -\nabla\cdot\mathbf{V}\rho\psi
+$$
+where $\rho$ is the fluid density, $\psi$ is the mixing ratio, and $\mathbf{V}$ is
+the fluid velocity. This is simplified by assuming that $\rho$ is constant and cancels out.
+The advection scheme for a finite volume method involves approximating this equation along an element
+edge at a certain time step and multiplying by the time step length to get an approximate scalar transport
+from element to element during the time step.
+The left hand side of this equation is a time derivative and for Omega this is handled by standard Runge—Kutta methods and
+not discussed here but can be found in the Skamarock and Gassmann paper.
+The right hand side of this equation is the spatial derivative that needs to be evaluated by the Omega horizontal advection scheme.
+Omega ocean uses unstructured Voronoi meshes and a
+finite volume scheme with $\psi$ defined on cell centers. A first order approximation
+to the spatial derivative of the left hand termacross a single edge of an element would be the difference of $\psi$ at
+element centroids divided by the distance from cell centroid to cell centroid.  This is expressed in equation (5) of the paper.
+For edge $i$ between elements $i+1/2$ and $i-1/2$,
+$$
+ \frac{\partial(u\psi_i)}{\partial x} = \frac{1}{\Delta x}[F_{i+1/2}(u\psi) - F_{i-1/2}(u\psi)] + O(\Delta x^2).
+$$
+Where $u$ is $\mathbf{V}\cdot\mathbf{n}$ where $\mathbf{n}$ is the unit normal along the edge being evaluated and $F$ is the evaluation
+of $\psi$ for the elements on each side of the edge being evaluated.
+This is used when the TracerHorzAdvOnCell user option is active and TracerHorzAdvTendencyOrder is 1.
+To make the comparison to higher order methods explicit, this first order method is equivalent to defining $\psi$ as a linear function,
+$\psi = c_0 + c_x x + c_y y$, between the two elements sharing the edge $i$ and taking the directed derivative along the edge.
+
+
+For higher order flux calculations, higher order derivatives of $\psi$ are needed and for the user option of
+TracerHorzAdvOnCell along with TracerHorzAdvTendencyOrder set to 2, a quadratic approximation of $\psi$ is used,
+$$
+\psi = c_0 + c_x x + c_y y + c_{xx} x^2 + c_{xy} xy + c_{yy} y^2,
+$$
+defined for each edge in the mesh. The six coefficients in this quadratic equation are calculated from
+a least squares fit of $\psi$ in a neighborhood of elements around an edge. The
+least squares fit is defined in equation (12) of the paper,
+$$
+ \mathbf{f = (P^T P)^{-1}P^T s = Bs}
+$$
+where the vector $\mathbf{f} = [c_0,c_x,c_y,c_{xx},c_{xy},c_{yy}]$ and $\mathbf{P}$ is an $m\times 6$ matrix with each row,
+$(1,x_i,y_i,x_i^2,x_iy_i,y_i^2)$, based on the cell center coordinates where there is one row for each of the $m$ elements in the neighborhood.
+This results in a $6\times m$ matrix $\mathbf{B}$. The values of $\mathbf{s}$ are the mixing ratios,
+$\mathbf{s}=[\psi_0, \psi_1,...,\psi_m]$, of the elements in the neighborhood at the time step being evaluated.
+Note that $\mathbf{B}$ depends on the mesh geometry only, so only needs to
+be computed once for each edge during initialization. Since only the second order terms are used to determine the
+second order derivatives needed for the higher order flux corrections to the lower order linear flux given above
+and the derivatives are taken along a line connecting
+cell centers, the amount of data that needs to be stored is minimized and $\mathbf{B}$ reduces to a vector.
+
+The second order derivatives of $\psi$ are then combined as shown in equation (11) of Skamarock and Gassmann to get
+a second order horizontal transport algorithm.
+
+The computation of these second order fluxes is complicated by having to do these computations on a sphere.
+Figure 2 from the Skamarock and Gassmann paper describs the slight modifications needed to compute the entries
+of the $\mathbf{P}$ matrix for a spherical geometry where the distances need to be measured as arc lengths along
+the surface.
+
+In practice, the convergence of this higher order algorithm on a sphere is slightly below
+second order. For the advection of a $\psi$ in the shape of a cosine bell
+being advected around a sphere and compared with the analytic solution, the spatial convergence for grid refinement
+is about order 1.7 as shown in  {numref}`tracer-higher-order-convergence`:
+
+```{figure} images/higher_order_tracer_convergence_on_sphere.jpeg
+:name:  tracer-higher-order-convergence
+:align: center
+:width: 600 px
+Tracer higer order convergence example of a cosine bell advected on a sphere showing an order 1.71 convergence rate
+```

components/omega/src/infra/Config.cpp

@@ -33,17 +33,7 @@ int Config::NumReadGroups       = 0;
 MPI_Comm Config::ConfigComm;
 Config Config::ConfigAll;
 
-//------------------------------------------------------------------------------
-// Constructor that creates configuration with a given name and an emtpy
-// YAML node that will be filled later with either a readAll or get call.
-Config::Config(const std::string &InName // [in] name of config, node
-) {
-
-   // If this is the first Config created, set variables for reading the
-   // input configuration file. In particular, to avoid too many MPI tasks
-   // reading the same input stream, we divide the tasks into groups and
-   // only one group at a time loads the full configuration from a stream
-
+void Config::Initialize() {
    if (NotInitialized) {
       // Determine MPI variables
       MachEnv *DefEnv = MachEnv::getDefault();
@@ -58,6 +48,18 @@ Config::Config(const std::string &InName // [in] name of config, node
 
       NotInitialized = false; // now initialized for future calls
    }
+}
+//------------------------------------------------------------------------------
+// Constructor that creates configuration with a given name and an emtpy
+// YAML node that will be filled later with either a readAll or get call.
+Config::Config(const std::string &InName // [in] name of config, node
+) {
+
+   // If this is the first Config created, set variables for reading the
+   // input configuration file. In particular, to avoid too many MPI tasks
+   // reading the same input stream, we divide the tasks into groups and
+   // only one group at a time loads the full configuration from a stream
+   Initialize();
 
    // Set the name - this is also used as the name of the root YAML node
    // (a map node) in this configuration.

components/omega/src/infra/Config.h

@@ -59,6 +59,8 @@ class Config {
  public:
    // Methods
 
+   static void Initialize();
+
    /// Constructor that creates a Configuration with a given name and an
    /// empty YAML::Node. The node will be filled later with either a read
    /// or get operation.

components/omega/src/ocn/CustomTendencyTerms.h

@@ -22,8 +22,8 @@
 //
 //===----------------------------------------------------------------------===//
 
-#include "HorzMesh.h"
-#include "TendencyTerms.h"
+#include "AuxiliaryState.h"
+#include "OceanState.h"
 #include "TimeMgr.h"
 
 namespace OMEGA {

components/omega/src/ocn/HorzMesh.cpp

@@ -17,7 +17,6 @@
 #include "Error.h"
 #include "IO.h"
 #include "Logging.h"
-#include "MachEnv.h"
 #include "OmegaKokkos.h"
 
 namespace OMEGA {
@@ -46,19 +45,21 @@ void HorzMesh::init() {
 
 HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
                    Decomp *MeshDecomp       //< [in] Decomp for the new mesh
-) {
-
+                   )
+    : CellID(MeshDecomp->CellID) ///< global cell ID for each local cell
+{
    MeshName = Name;
 
    // Retrieve mesh files name from Decomp
    MeshFileName = MeshDecomp->MeshFileName;
 
    // Retrieve mesh cell/edge/vertex totals from Decomp
-   NCellsHalo  = MeshDecomp->NCellsHalo;
-   NCellsHaloH = MeshDecomp->NCellsHaloH;
-   NCellsOwned = MeshDecomp->NCellsOwned;
-   NCellsAll   = MeshDecomp->NCellsAll;
-   NCellsSize  = MeshDecomp->NCellsSize;
+   NCellsHalo   = MeshDecomp->NCellsHalo;
+   NCellsHaloH  = MeshDecomp->NCellsHaloH;
+   NCellsOwned  = MeshDecomp->NCellsOwned;
+   NCellsAll    = MeshDecomp->NCellsAll;
+   NCellsSize   = MeshDecomp->NCellsSize;
+   NCellsGlobal = MeshDecomp->NCellsGlobal;
 
    NEdgesHalo     = MeshDecomp->NEdgesHalo;
    NEdgesHaloH    = MeshDecomp->NEdgesHaloH;
@@ -67,6 +68,7 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
    NEdgesSize     = MeshDecomp->NEdgesSize;
    MaxCellsOnEdge = MeshDecomp->MaxCellsOnEdge;
    MaxEdges       = MeshDecomp->MaxEdges;
+   NEdgesGlobal   = MeshDecomp->NEdgesGlobal;
 
    NVerticesHalo  = MeshDecomp->NVerticesHalo;
    NVerticesHaloH = MeshDecomp->NVerticesHaloH;