Merge pull request #223 from mark-petersen/omega/design-doc-v1-eqns

Omega V1 Governing Equations Design Document

Author: mark-petersen

components/omega/doc/conf.py

@@ -41,6 +41,7 @@
     "sphinx.ext.mathjax",
 ]
 
+numfig = True
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
 
@@ -73,6 +74,7 @@
 myst_dmath_double_inline = True
 myst_enable_checkboxes = True
 suppress_warnings = ['myst.header']
+myst_heading_start_level = 1
 
 # -- HTML output -------------------------------------------------
 

components/omega/doc/design/OmegaV0ShallowWater.md

@@ -109,7 +109,7 @@ Next we replace the non-linear advection term with the right-hand side of the ve
 $$
 \boldsymbol{u} \cdot \nabla \boldsymbol{u} &= (\nabla \times \boldsymbol{u}) \times \boldsymbol{u} + \nabla \frac{|\boldsymbol{u}|^2}{2} \\
 &= \{\boldsymbol{k} \cdot (\nabla \times \boldsymbol{u})\} \boldsymbol{k} \times \boldsymbol{u} + \nabla \frac{|\boldsymbol{u}|^2}{2} \\
-&= \omega \boldsymbol{u}^{\perp} + \nabla K.
+&= \zeta \boldsymbol{u}^{\perp} + \nabla K.
 $$
 
 The governing equations for Omega-0 in continuous form are then
@@ -141,10 +141,10 @@ The drag consists of simple Rayleigh drag for spin-up as well as quadratic botto
 
 Here $q$ is the potential vorticity, so that the term
 $$
-q\left(h\boldsymbol{u}^{\perp}\right) = \frac{\omega + f}{h}\left(h\boldsymbol{u}^{\perp}\right)
- = \omega \boldsymbol{u}^{\perp} + f \boldsymbol{u}^{\perp}
+q\left(h\boldsymbol{u}^{\perp}\right) = \frac{\zeta + f}{h}\left(h\boldsymbol{u}^{\perp}\right)
+ = \zeta \boldsymbol{u}^{\perp} + f \boldsymbol{u}^{\perp}
 $$
-is composed of the rotational part of the advection $\omega \boldsymbol{u}^{\perp}$ and the Coriolis term $f \boldsymbol{u}^{\perp}$. This term is discussed in [Ringler et al. 2010](https://www.sciencedirect.com/science/article/pii/S0021999109006780) sections 2.1 and 2.2.
+is composed of the rotational part of the advection $\zeta \boldsymbol{u}^{\perp}$ and the Coriolis term $f \boldsymbol{u}^{\perp}$. This term is discussed in [Ringler et al. 2010](https://www.sciencedirect.com/science/article/pii/S0021999109006780) sections 2.1 and 2.2.
 
 The thickness equation (2) is derived from conservation of mass for a fluid with constant density, which reduces to conservation of volume. The model domain uses fixed horizontal cells with horizontal areas that are constant in time, so the area drops out and only the layer thickness $h$ remains as the prognostic variable.
 
@@ -191,6 +191,7 @@ $$
 \mathcal{F}_e =  C_W \frac{(u_W - u_e)\left|u_W - u_e\right|}{[h_i]_e}
 $$
 
+(32-variable-definitions)=
 ### 3.2 Variable Definitions
 
 Table 1. Definition of variables
@@ -215,7 +216,7 @@ Table 1. Definition of variables
 | $K$                 | kinetic energy              | m$^2$/s$^2$  | cell     | KineticEnergyCell  | $K = \left\| {\boldsymbol u} \right\|^2 / 2$ |
 | $l_e $      | edge length (vertex span) | m    | edge     | DvEdge  |                                                              |
 | $n_{e,i}$ | edge normal sign | unitless | edge | EdgeSignOnCell| also by index ordering of CellsOnEdge |
-| $q$                 | potential vorticity         | 1/m/s    | vertex   |   | $q = \eta/h = \left(\omega+f\right)/h$                                                 |
+| $q$                 | potential vorticity         | 1/m/s    | vertex   |   | $q = \left(\zeta+f\right)/h$                                                 |
 | $Ra$                | Rayleigh drag coefficient   | 1/s      | constant |   |                                                              |
 | $t$                 | time                        | s        | none     |   |                                                              |
 | $t_{e,v}$ | edge tangential sign | unitless | edge | EdgeSignOnVertex | also by index ordering of VerticesOnEdge |
@@ -225,13 +226,12 @@ Table 1. Definition of variables
 | ${\boldsymbol u}_W$ | wind velocity               | m/s      | edge     |   |                                                              |
 | $w_{e,e'}$ | tangential edge weights| unitless | edge     |   | weights defined in [Ringler et al. 2010](https://www.sciencedirect.com/science/article/pii/S0021999109006780) eqn 24 |
 |${\tilde w}_{e,e'}$ | normalized tangential edge weights| unitless | edge     | WeightsOnEdge   | normalized weights, ${\tilde w}_{e,e'} = \tfrac{l_{e'}}{d_e} w_{e,e'}$ |
-| $\eta$              | absolute vorticity          | 1/s      | vertex   |   | $\eta=\omega + f$ |
 | $\kappa_2$          | tracer diffusion            | m$^2$/s    | cell     |   |                                                              |
 | $\kappa_4$          | biharmonic tracer diffusion | m$^4$/s    | cell     |   |                                                              |
 | $\nu_2$             | viscosity                   | m$^2$/s    | edge     |   |                                                              |
 | $\nu_4$             | biharmonic viscosity        | m$^4$/s    | edge     |   |                                                              |
 | $\phi$              | tracer                      | varies | cell     |   | units may be kg/m$^3$ or similar |
-| $\omega$             | relative vorticity          | 1/s      | vertex   |  RelativeVorticity | $\omega={\boldsymbol k} \cdot \left( \nabla \times {\boldsymbol u}\right)$ |
+| $\zeta$             | relative vorticity          | 1/s      | vertex   |  RelativeVorticity | $\zeta={\boldsymbol k} \cdot \left( \nabla \times {\boldsymbol u}\right)$ |
 
 <!--- Note: Table created with [markdown table generator](https://www.tablesgenerator.com/markdown_tables) and original [google sheet](https://docs.google.com/spreadsheets/d/1rz-QXDiwfemq5NpSR1XsvomI7aSKQ1myTNweCY4afcE/edit#gid=0). --->
 
@@ -274,11 +274,12 @@ The definitions of geometric variables may be found in
 [Ringler et al. 2010](https://www.sciencedirect.com/science/article/pii/S0021999109006780), and
 [Thuburn and Cotter 2012](https://epubs.siam.org/doi/10.1137/110850293).
 
-### 3.2 Operator Formulation
+(33-operator-formulation)=
+### 3.3 Operator Formulation
 
 The TRiSK formulation of the discrete operators are as follows. See [Bishnu et al. 2023](https://gmd.copernicus.org/articles/16/5539/2023) section 4.1 and Figure 1 for a description and documentation of convergence rates, as well as [Bishnu et al. 2021](https://doi.org/10.5281/zenodo.7439539). All TRiSK spatial operators show second-order convergence on a uniform hexagon grid, except for the curl on vertices, which is first order. The curl interpolated from vertices to cell centers regains second order convergence. The rates of convergence are typically less than second order on nonuniform meshes, including spherical meshes.
 
-#### 3.2.1. Divergence
+#### 3.3.1 Divergence
 The divergence operator maps a vector field's edge normal component $F_e$ to a cell center ([Ringler et al. 2010](https://www.sciencedirect.com/science/article/pii/S0021999109006780) eqn 21),
 
 $$
@@ -302,7 +303,7 @@ D_i
 $$
 
 
-#### 3.2.2. Gradient
+#### 3.3.2 Gradient
 
 The gradient operator maps a cell-centered scalar to an edge-normal vector component
 ([Ringler et al. 2010](https://www.sciencedirect.com/science/article/pii/S0021999109006780) eqn 22),
@@ -322,7 +323,7 @@ $$
 $$
 where $\left\{i_1, i_2\right\}$ are the cells neighboring edge $e$.  The indices $\left\{i_1, i_2\right\}$ are ordered such that the normal vector ${\bf n}$ points from cell $i_1$ to cell $i_2$. In the code ${\bf n}$ points from `CellsOnEdge(IEdge, 0)` to `CellsOnEdge(IEdge, 1)`.
 
-#### 3.2.3. Curl
+#### 3.3.3 Curl
 The curl operator maps a vector's edge normal component $u_e$ to a scalar field at vertex
 ([Ringler et al. 2010](https://www.sciencedirect.com/science/article/pii/S0021999109006780) eqn 23),
 
@@ -337,7 +338,7 @@ where ${\hat A}_v$ is the area of the dual mesh cell $v$, i.e. the triangle surr
 ${\bf x}_v$, then $t_{e,v}=1$. The summation is over $e\in EV(v)$, which are the edges that terminate at vertex $v$. There are *always* three edges that terminate at each vertex for Voronoi Tessellations, and four edges on each vertex for quadrilateral meshes, unless neighboring cells are missing for land boundaries. Again, the subscript $v$ is dropped and a vertex is assumed for the location of the curl.
 
 
-#### 3.2.4. Perpendicular vector component
+#### 3.3.4 Perpendicular vector component
 
 The perpendicular component a vector field is defined as
 $$
@@ -370,34 +371,34 @@ u_e^\perp = \sum_{e'\in ECP(e)} {\tilde w}_{e,e'} \, u_{e'}.
 $$
 This simple operation can be seen in MPAS-Ocean in the subroutine `ocn_diagnostic_solve_vortVel` in the file `mpas_ocn_diagnostics.F`.
 
-#### 3.2.5. Perpendicular Gradient
+#### 3.3.5 Perpendicular Gradient
 The gradient of a scalar at the middle of an edge, pointing tangentially along the edge (from one vertex to the other), is sometimes used. For example, the del2 formulation requires the perpendicular gradient of vorticity. This is called the perpendicular gradient because the standard gradient is normal to the edge.
 
 The perpendicular gradient maps a scalar at vertices to an edge-tangential vector component,
 $$
-\left(\nabla \omega_v\right)_e^\perp = \frac{1}{l_e} \sum_{v\in VE(e)} -t_{e,v}\omega_v.
+\left(\nabla \zeta_v\right)_e^\perp = \frac{1}{l_e} \sum_{v\in VE(e)} -t_{e,v}\zeta_v.
 $$
 Like the standard gradient, in practice we drop the sign indicator $t_{e,v}$ and rewrite it using the index ordering,
 $$
-\left(\nabla \omega_v\right)_e^\perp = \frac{\omega_{v2} - \omega_{v1}}{l_e},
+\left(\nabla \zeta_v\right)_e^\perp = \frac{\zeta_{v2} - \zeta_{v1}}{l_e},
 $$
 where the positive vector ${\bf n}^\perp$ is 90$^o$ to the left of ${\bf n}$. The indices are ordered such that ${\bf n}^\perp$ points from $v_1$ to $v_2$, which corresponds to `VerticesOnEdge(IEdge, 0)` and `VerticesOnEdge(IEdge, 1)` in the code.
 
-#### 3.2.6. Cell to Edge Interpolation
+#### 3.3.6 Cell to Edge Interpolation
 The mid-point average of a scalar from cell centers to the adjoining edge is
 $$
 [h_i]_e = \frac{1}{2} \sum_{i\in CE(e)} h_i.
 $$
 In a Voronoi tessellation the edge is defined to be at the mid-point between the two cell centers, so this formula is identical to an area-weighted average using the triangles between edge $e$ and the neighboring cell centers.
 
-#### 3.2.7. Vertex to Edge Interpolation
+#### 3.3.7 Vertex to Edge Interpolation
 The mid-point average of a scalar from vertices to the middle of the connecting edge is
 $$
-[\omega_v]_e = \frac{1}{2} \sum_{v\in VE(e)} \omega_v.
+[\zeta_v]_e = \frac{1}{2} \sum_{v\in VE(e)} \zeta_v.
 $$
 The is a distance-weighted average, since the edge quantity lives at the mid-point between the vertices. One could alternatively compute an area-weighted average using the dual-mesh cell area ${\hat A}_v$ surrounding each vertex, but that is not done and would result in very small differences.
 
-#### 3.2.8. Cell to Vertex Interpolation
+#### 3.3.8 Cell to Vertex Interpolation
 The area-weighted average of a scalar at a vertex from the three surrounding cells is
 $$
 [h_i]_v = \frac{1}{{\hat A}_v} \sum_{i\in CV(v)} h_i {\tilde A}_{v,i}.
@@ -408,21 +409,22 @@ $$
 {\hat A}_v = \sum_{i\in CV(v)} {\tilde A}_{v,i}.
 $$
 
-#### 3.2.9. Vertex to Cell Interpolation
+#### 3.3.9 Vertex to Cell Interpolation
 The area-weighted average of a scalar at a cell from the surrounding vertices is
 $$
 [h_v]_i = \frac{1}{A_i} \sum_{v\in VC(i)} h_v {\tilde A}_{v,i}
 $$
 
-#### 3.2.10. Vector from Edge to Cell
+#### 3.3.10 Vector from Edge to Cell
 
 The prognostic velocity variable on the edge is the edge normal velocity, $u_e$. The tangential velocity $u_e^\perp$ is computed diagnostically from $u_e$. In addition, the full vector may be computed at the cell center from the edge normal velocities $u_e$ of the edges on that cell. That is done in MPAS with radial basis functions, and is explained in this [previous design document](https://github.com/MPAS-Dev/MPAS-Documents/blob/master/shared/rbf_design/rbf.pdf).
 
-### 3.3 Momentum Terms
+(34-momentum-terms)=
+### 3.4 Momentum Terms
 
 The computation of each term in (4-6) is now described in detail, along with alternative formulations.
 
-#### 3.3.1. Kinetic energy gradient
+#### 3.4.1 Kinetic energy gradient
 
 The kinetic energy gradient term is the non-rotational part of the nonlinear advection. Fundamentally, it maps the prognostic edge-normal velocity from edges back to an edge scalar quantity. We use the standard gradient formulation from cell center to edge,
 
@@ -465,81 +467,84 @@ $$
 for the final kinetic energy at the cell center. Note that addition of $[K_v]_i$ enlarges the stencil. One could also use $u_e^{\perp}$ and compute the kinetic energy at the edge itself in order to enlarge the stencil, but that method is not used here.
 See [Calandrini et al. 2021](https://www.sciencedirect.com/science/article/pii/S146350032100161X) section 2.3 for more information.
 
-#### 3.3.2. Potential vorticity term
+#### 3.4.2 Potential vorticity term
 
 The potential vorticity term, $q\left(h\boldsymbol{u}^{\perp}\right)$, includes the rotational part of the advection. It may be computed in two ways,
 
 $$
-\text{Option 1:  } \left[ \frac{\omega_v +f_v}{[h_i]_v}\right]_e\left([h_i]_e u_e^{\perp}\right)
+\text{Option 1:  } \left[ \frac{\zeta_v +f_v}{[h_i]_v}\right]_e\left([h_i]_e u_e^{\perp}\right)
 $$
 
 $$
-\text{Option 2:  }([\omega_v]_e +f_e) u_e^\perp
+\text{Option 2:  }([\zeta_v]_e +f_e) u_e^\perp
 $$
 
 The first computes the potential vorticity $q_v$ at the vertex and interpolates that quantity to the edge, which is what is done in MPAS-Ocean. One may also cancel the thickness $h$ (ignoring the interpolated locations) and use option 2. Additional interpolation options and results are presented in [Calandrini et al. 2021](https://www.sciencedirect.com/science/article/pii/S146350032100161X)
 
-#### 3.3.3. Sea surface height gradient
+#### 3.4.3 Sea surface height gradient
 The sea surface height (SSH) gradient uses the standard gradient formulation from cell center to edge,
 
 $$
 -g\nabla(h_i-b_i) &= -g\frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i}(h_i-b_i)\\
 &= -g\frac{(h_{i2}-b_{i2}) - (h_{i1}-b_{i1}) }{d_e}.
 $$
 
-#### 3.3.4. Del2 momentum dissipation
+(344-del2-momentum-dissipation)=
+#### 3.4.4 Del2 momentum dissipation
 The Del2, or Laplacian operator, viscous momentum dissipation maps edge-normal velocity back to the edge-normal component of the Laplacian. In TRiSK this is done with the vorticity-divergence formulation, which may be written as
 
 $$
 \nabla^2 {\bf u} &=  \nabla \left( \nabla\cdot{\bf u}\right) - \nabla \times \left( \nabla \times {\bf u} \right) \\
- &=  \nabla {\bf D} - \nabla \times \omega \\
- &=  \nabla {\bf D} - {\bf k} \times \nabla \omega \\
- &=  \nabla {\bf D} - \nabla^\perp \omega.
+ &=  \nabla {\bf D} - \nabla \times \zeta \\
+ &=  \nabla {\bf D} - {\bf k} \times \nabla \zeta \\
+ &=  \nabla {\bf D} - \nabla^\perp \zeta.
 $$
 This formulation is also mentioned in [Gassmann 2011](https://linkinghub.elsevier.com/retrieve/pii/S0021999111000325) (equation 17), [Gassman 2018](https://onlinelibrary.wiley.com/doi/10.1002/qj.3294) (equation 44) and [Lapolli et al. 2024](https://www.sciencedirect.com/science/article/pii/S1463500324000222) (section 4.5.2).
 
 For our discretization, the full Del2 term is written as
 $$
-\nu_2 \nabla^2 u_e &= \nu_2 \left( \nabla D_i - \nabla^\perp \omega_v\right) \\
-&= \nu_2 \left( \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i}D_i - \frac{1}{l_e} \sum_{v\in VE(e)} -t_{e,v}\omega_v\right) \\
-&= \nu_2 \left( \frac{D_{i2} - D_{i1}}{d_e} - \frac{\omega_{v2} - \omega_{v1}}{l_e} \right),
+\nu_2 \nabla^2 u_e &= \nu_2 \left( \nabla D_i - \nabla^\perp \zeta_v\right) \\
+&= \nu_2 \left( \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i}D_i - \frac{1}{l_e} \sum_{v\in VE(e)} -t_{e,v}\zeta_v\right) \\
+&= \nu_2 \left( \frac{D_{i2} - D_{i1}}{d_e} - \frac{\zeta_{v2} - \zeta_{v1}}{l_e} \right),
 $$
 where the ordering of indices $\{i_i, i_2\}$ and $\{v_1, v_2\}$ are explained in the gradient operator sections above.
 
 An alternative formulation for the Del2 dissipation on an unstructured mesh is presented in section 4.2 and Appendix B of
 [Gassman 2018](https://onlinelibrary.wiley.com/doi/10.1002/qj.3294).
 
-#### 3.3.5. Del4 momentum dissipation
+(345-del4-momentum-dissipation)=
+#### 3.4.5 Del4 momentum dissipation
 The Del4, or biharmonic, momentum dissipation also maps edge-normal velocity back to the edge-normal component of the Laplacian. This is done with two applications of the Del2 operator above.
 
 $$
 - \nu_4 \nabla^4 u_e = - \nu_4 \nabla^2 \left( \nabla^2 u_e \right)
 $$
 
-#### 3.3.6. Rayleigh Drag
+#### 3.4.6 Rayleigh Drag
 Rayleigh drag is simply linear drag, applied to all levels. It is typically only used during the spin-up process to damp large velocities during the initial adjustment. The Rayleigh coefficient $Ra$ is simply a scalar constant,
 
 $$
 - Ra \: u_e.
 $$
 
-#### 3.3.7. Bottom drag
+#### 3.4.7 Bottom drag
 Bottom drag is more relevant to layered models than to shallow water systems, but is included here for completeness. It is a quadratic drag applied to the edge velocity,
 
 $$
 - C_D \frac{u_e\left|u_e\right|}{[h_i]_e}.
 $$
 
-#### 3.3.8. Wind forcing
+#### 3.4.8 Wind forcing
 Wind forcing has the same form as the bottom drag, but the forcing is the difference between the current velocity and the wind $u_W$, interpolated and projected to the edge normal direction,
 
 $$
 - C_W \frac{(u_W - u_e)\left|u_W - u_e\right|}{[h_i]_e}.
 $$
 
-### 3.4 Thickness and Tracer Terms
+(35-thickness-and-tracer-terms)=
+### 3.5 Thickness and Tracer Terms
 
-#### 3.4.1. Tracer advection
+#### 3.5.1 Tracer advection
 
 There are many schemes available for tracer advection. Simple schemes include centered advection and upstream. MPAS-Ocean uses Flux Corrected Transport, which is fourth order under normal conditions, and reduces to third order to preserve monotonicity.
 
@@ -558,7 +563,8 @@ Note that the thickness advection is identical to the tracer advection when $\ph
 
 More details of the tracer advection scheme will be given in a future design document.
 
-#### 3.4.2. Del2 tracer diffusion
+(352-del2-tracer-diffusion)=
+#### 3.5.2 Del2 tracer diffusion
 Tracer diffusion is applied with a Laplacian operator on the cell-centered tracer phi, and the product of the operator is also at the cell center. The Laplacian may be written as the divergence of the gradient,
 
 $$
@@ -567,23 +573,24 @@ $$
 \kappa_2 h_i \nabla \cdot \left( \nabla \phi_i \right).
 $$
 
-and the stencils in Section 3.2 are used. Here $\kappa_2$ is the del2 diffusion coefficient, and the operator is thickness-weighted by $h_i$. The position of a gradient is assumed to be at an edge, and a divergence at the cell center. This could be written explicitly as
+and the stencils in Section 3.3 are used. Here $\kappa_2$ is the del2 diffusion coefficient, and the operator is thickness-weighted by $h_i$. The position of a gradient is assumed to be at an edge, and a divergence at the cell center. This could be written explicitly as
 
 $$
 \kappa_2 h_i \left( \nabla^2 \phi_i \right)_i
 =
 \kappa_2 h_i \left( \nabla \cdot \left( \nabla \phi_i \right)_e \right)_i.
 $$
 
-#### 3.4.3. Del4 tracer diffusion
+(353-del4-tracer-diffusion)=
+#### 3.5.3 Del4 tracer diffusion
 The del4 tracer diffusion is simply the Laplacian operator applied twice,
 
 $$
 - \kappa_4 h_i \nabla^4 \phi_i
 =
 - \kappa_4 h_i \nabla^2 \left( \nabla^2 \phi_i \right).
 $$
-The del2 operator using the divergence of the gradient in the last section is used, with the simple stencils from Section 3.2.
+The del2 operator using the divergence of the gradient in the last section is used, with the simple stencils from Section 3.3.
 
 ## 4 Design