examples.md

Examples

1D Column examples

1D Wave Equation

The 1D linear wave example in examples/column/wave.jl demonstrates solving the following system of PDEs:

Equations and Discretizations

The first-order system for the 1D linear wave equation is:

$$\begin{aligned} \frac{\partial u}{\partial t} &= -\frac{\partial p}{\partial z} \\\ \frac{\partial p}{\partial t} &= -\frac{\partial u}{\partial z} \end{aligned}$$

This is discretized using finite difference operators with a staggered grid:

• u is defined at cell centers (CenterFiniteDifferenceSpace)
• p is defined at cell faces (FaceFiniteDifferenceSpace)

The discrete form is:

$$\begin{aligned} \frac{d u}{dt} &\approx - \nabla \cdot p \\\ \frac{d p}{dt} &\approx - \nabla u \end{aligned}$$

These are implemented in code using:

∂f = GradientC2F(left = SetValue(0.0), right = SetValue(0.0))
∂c = DivergenceF2C()

@. dp = -WVector(∂f(u))
@. du = -∂c(p)

with Dirichlet boundary conditions (SetValue(0.0)) applied at both ends.

Prognostic Variables

Variable	Description	Space
u	Scalar velocity	CenterFiniteDifferenceSpace
p	Face-centered pressure flux	FaceFiniteDifferenceSpace

Time Integration

The system is solved using the SSPRK33 3rd order Runge Kutta time integrator from the OrdinaryDiffEqSSPRK package:

prob = ODEProblem(tendency!, Y, (0.0, 4π))
sol = solve(prob, SSPRK33(), dt = 0.01, saveat = 0.1)

ClimaCore Operators Used

Operator	Purpose	Code Example
GradientC2F(...)	Computes the gradient from center (C) to face (F) locations	∂f = GradientC2F(...)
DivergenceF2C()	Computes the divergence from face (F) to center (C) locations	∂c = DivergenceF2C()
WVector(...)	Constructs a vertical vector field at faces	@. dp = -WVector(∂f(u))
SetValue(val)	Enforces a boundary condition by setting values at the boundary	GradientC2F(left = SetValue(0.0))

• GradientC2F: Computes ∂u/∂z, placing the result on face points.
• DivergenceF2C: Computes ∂p/∂z, placing the result on center points.
• WVector: Wraps scalars into 1D vertical vectors for compatibility with DivergenceF2C.
• SetValue(0.0): Applies fixed-value Dirichlet boundary conditions.

Problem Flow and Setup

1. Mesh and Domain: Create a 1D interval domain [0, 4π] with 30 elements.
2. Function Spaces: Construct staggered finite-difference spaces for u (cell centers) and p (cell faces).
3. Initial Conditions: Set u(z) = sin(z) and p(z) = 0.
4. Differential Operators: Apply GradientC2F and DivergenceF2C with Dirichlet boundaries.
5. Tendency Function: Compute time derivatives of u and p.
6. ODE Integration: Solve the PDE system using SSPRK33 with OrdinaryDiffEqSSPRK.
7. Output: Save the final state and animated wave evolution using Plots.jl.

2D Cartesian examples

Flux Limiters advection

The 2D Cartesian advection/transport example in examples/plane/limiters_advection.jl demonstrates the application of flux limiters in the horizontal direction, namely QuasiMonotoneLimiter, in a 2D Cartesian domain.

Equations and discretizations

Mass

Follows the continuity equation

$$\begin{equation} \frac{\partial}{\partial t} \rho = - \nabla \cdot(\rho \boldsymbol{u}) . \label{eq:2d-plane-advection-lim-continuity} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho \approx - wD[ \rho \boldsymbol{u}] . \label{eq:2d-plane-advection-lim-discrete-continuity} \end{equation}$$

Tracers

For the tracer concentration per unit mass q, the tracer density (scalar) \rho q follows the advection/transport equation

$$\begin{equation} \frac{\partial}{\partial t} \rho q = - \nabla \cdot(\rho q \boldsymbol{u}) + g(\rho, q). \label{eq:2d-plane-advection-lim-tracers} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho q \approx - wD[ \rho q \boldsymbol{u}] + g(\rho, q), \label{eq:2d-plane-advection-lim-discrete-tracers} \end{equation}$$

where g(\rho, q) = - \nu_4 [\nabla^4_h (\rho q)] represents the horizontal hyperdiffusion operator, with \nu_4 (measured in m^4/s) the hyperviscosity constant coefficient (set equal to zero by default in the example).

Currently tracers are only treated explicitly in the time discretization.

Prognostic variables

• \rho: density measured in kg/m³.
• \boldsymbol{u} velocity, a vector measured in m/s. Since this is a 2D problem, \boldsymbol{u} \equiv \boldsymbol{u}_h.
• \rho q: the tracer density scalar, where q is the tracer concentration per unit mass.

Differentiation operators

Because this is a purely 2D problem, there is no staggered vertical discretization, hence, there is no need of specifying variables at cell centers, faces or to reconstruct from faces to centers and vice versa.

• wD is the discrete horizontal weak spectral divergence, called wdiv in the example code.
• G is the discrete horizontal strong spectral gradient, called grad in the example code.

To discretize the hyperdiffusion operator, g(\rho, q) = - \nu_4 [\nabla^4 (\rho q)], in the horizontal direction, we compose the horizontal weak divergence, wD, and the horizontal gradient operator, G_h, twice, with an intermediate call to weighted_dss! between the two compositions, as in [g_2(\rho, g) \circ DSS(\rho, q) \circ g_1(\rho, q)], with:

• g_1(\rho, q) = wD(G_h(q))
• DSS(\rho, q) = DSS(g_1(\rho q))
• g_2(\rho, q) = -\nu_4 wD(\rho G_h(\rho q))
  • with \nu_4 the hyperviscosity coefficient.

Problem flow and set-up

This test case is set up in a Cartesian planar domain [-2 \pi, 2 \pi]^2, doubly periodic.

The flow was chosen to be a horizontal uniform rotation. Moreover, the flow is reversed halfway through the time period so that the tracer blobs go back to its initial configuration (using the same speed scaling constant which was derived to account for the distance travelled in all directions in a half period).

$$\begin{align} u &= -u_0 (y - c_y) \cos(\pi t / T_f) \nonumber \\\ v &= u_0 (x - c_x) \cos(\pi t / T_f) \label{eq:2d-plane-advection-lim-flow} \end{align}$$

where u_0 = \pi / 2 is the speed scaling factor to have the flow reversed halfway through the time period, \boldsymbol{c} = (c_x, c_y) is the center of the rotational flow, which coincides with the center of the domain, and T_f = 2 \pi is the final simulation time, which coincides with the temporal period to have a full rotation.

This example is set up to run with three possible initial conditions:

• cosine_bells
• gaussian_bells
• cylinders: two 2D slotted cylinders (test case available in the literature, cfr: GubaOpt2014).

Application of Flux Limiters

Because this is a fully 2D problem, the application of limiters does not affect the order of operations, which is implemented as follows:

1. Horizontal trasport with hyperdiffusion (with weak divergence wD)
2. Horizontal flux limiters
3. DSS

3D Cartesian examples

Flux Limiters advection

The 3D Cartesian advection/transport example in examples/hybrid/box/limiters_advection.jl demonstrates the application of flux limiters in the horizontal direction, namely QuasiMonotoneLimiter, in a hybrid Cartesian domain. It also demonstrates the usage of the high-order upwinding scheme in the vertical direction, called Upwind3rdOrderBiasedProductC2F.

Equations and discretizations

Mass

Follows the continuity equation

$$\begin{equation} \frac{\partial}{\partial t} \rho = - \nabla \cdot(\rho \boldsymbol{u}) . \label{eq:3d-box-advection-lim-continuity} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho \approx - D_h[ \rho (\boldsymbol{u}_h + I^c(\boldsymbol{u}_v))] - D^c_v[I^f(\rho \boldsymbol{u}_h)) + I^f(\rho) \boldsymbol{u}_v)] . \label{eq:3d-box-advection-lim-discrete-continuity} \end{equation}$$

Tracers

For the tracer concentration per unit mass q, the tracer density (scalar) \rho q follows the advection/transport equation

$$\begin{equation} \frac{\partial}{\partial t} \rho q = - \nabla \cdot(\rho q \boldsymbol{u}) + g(\rho, q). \label{eq:3d-box-advection-lim-tracers} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho q \approx - D_h[ \rho q (\boldsymbol{u}_h + I^c(\boldsymbol{u}_v))] - D^c_v\left[I^f(\rho q) U^f\left(I^f(\boldsymbol{u}_h) + \boldsymbol{u}_v, \frac{\rho q}{\rho} \right) \right] + g(\rho, q), \label{eq:3d-box-advection-lim-discrete-tracers} \end{equation}$$

where g(\rho, q) = - \nu_4 [\nabla^4_h (\rho q)] represents the horizontal hyperdiffusion operator, with \nu_4 (measured in m^4/s) the hyperviscosity constant coefficient.

Currently tracers are only treated explicitly in the time discretization.

Prognostic variables

• \rho: density measured in kg/m³. This is discretized at cell centers.
• \boldsymbol{u} velocity, a vector measured in m/s. This is discretized via \boldsymbol{u} = \boldsymbol{u}_h + \boldsymbol{u}_v where
  • \boldsymbol{u}_h = u_1 \boldsymbol{e}^1 + u_2 \boldsymbol{e}^2 is the projection onto horizontal covariant components (covariance here means with respect to the reference element), stored at cell centers.
  • \boldsymbol{u}_v = u_3 \boldsymbol{e}^3 is the projection onto the vertical covariant components, stored at cell faces.
• \rho q: the tracer density scalar, where q is the tracer concentration per unit mass, is stored at cell centers.

Operators

We make use of the following operators

Reconstructions

• I^c is the face-to-center reconstruction operator, called first_order_If2c in the example code.
• I^f is the center-to-face reconstruction operator, called first_order_Ic2f in the example code.
  • Currently this is just the arithmetic mean, but we will need to use a weighted version with stretched vertical grids.
• U^f is the center-to-face upwind product operator, called third_order_upwind_c2f in the example code
  • This operator is of third-order of accuracy (when used with a constant vertical velocity and some reduced, but still high-order for non constant vertical velocity).

Differentiation operators

• D_h is the discrete horizontal strong spectral divergence, called hdiv in the example code.
• wD_h is the discrete horizontal weak spectral divergence, called hwdiv in the example code.
• D^c_v is the face-to-center vertical divergence, called vdivf2c in the example code.
  • This example uses advective fluxes equal to zero at the top and bottom boundaries.
• G_h is the discrete horizontal spectral gradient, called hgrad in the example code.

To discretize the hyperdiffusion operator, g(\rho, q) = - \nu_4 [\nabla^4 (\rho q)], in the horizontal direction, we compose the horizontal weak divergence, wD_h, and the horizontal gradient operator, G_h, twice, with an intermediate call to weighted_dss! between the two compositions, as in [g_2(\rho, g) \circ DSS(\rho, q) \circ g_1(\rho, q)], with:

• g_1(\rho, q) = wD_h(G_h(q))
• DSS(\rho, q) = DSS(g_1(\rho q))
• g_2(\rho, q) = -\nu_4 wD_h(\rho G_h(\rho q))
  • with \nu_4 the hyperviscosity coefficient.

Application of Flux Limiters

!!! note

Since we use flux limiters that limit only operators defined in the spectral space (i.e., they are applied level-wise in the horizontal direction), the application of limiters has to follow a precise order in the sequence of operations that specifies the total tendency.

The order of operations should be the following:

1. Horizontal transport (with strong divergence D_h)
2. Horizontal Flux Limiters
3. Horizontal hyperdiffusion (with weak divergence wD_h)
4. Vertical transport
5. DSS

Problem flow and set-up

This test case is set up in a Cartesian (box) domain [-2 \pi, 2 \pi]^2 \times [0, 4 \pi] ~\textrm{m}^3, doubly periodic in the horizontal direction, but not in the vertical direction.

The flow was chosen to be a spiral, i.e., so to have a horizontal uniform rotation, and a vertical velocity \boldsymbol{u}_v \equiv w = 0 at the top and bottom boundaries, and \boldsymbol{u}_v \equiv w = 1 in the center of the domain. Moreover, the flow is reversed in all directions halfway through the time period so that the tracer blobs go back to its initial configuration (using the same speed scaling constant which was derived to account for the distance travelled in all directions in a half period).

$$\begin{align} u &= -u_0 (y - c_y) \cos(\pi t / T_f) \nonumber \\\ v &= u_0 (x - c_x) \cos(\pi t / T_f) \nonumber \\\ w &= u_0 \sin(\pi z / z_m) \cos(\pi t / T_f) \nonumber \label{eq:3d-box-advection-lim-flow} \end{align}$$

where u_0 = \pi / 2 is the speed scaling factor to have the flow reversed halfway through the time period, \boldsymbol{c} = (c_x, c_y) is the center of the rotational flow, which coincides with the center of the domain, z_m = 4 \pi is the maximum height of the domain, and T_f = 2 \pi is the final simulation time, which coincides with the temporal period to have a full rotation in the horizontal direction.

This example is set up to run with three possible initial conditions:

• cosine_bells
• gaussian_bells
• slotted_spheres: a slight modification of the 2D slotted cylinder test case available in the literature (cfr: GubaOpt2014).

Application of Flux Limiters

Because this is a Cartesian 3D problem, the application of limiters does not affect the order of operations, which is implemented as follows:

1. Horizontal transport + hyperdiffusion (with weak divergence wD_h)
2. Horizontal flux limiters
3. Vertical transport
4. DSS

2D Sphere examples

Flux Limiters advection

The 2D sphere advection/transport example in examples/sphere/limiters_advection.jl demonstrates the application of flux limiters in the horizontal direction, namely QuasiMonotoneLimiter, in a 2D spherical domain.

Equations and discretizations

Mass

Follows the continuity equation

$$\begin{equation} \frac{\partial}{\partial t} \rho = - \nabla \cdot(\rho \boldsymbol{u}) . \label{eq:2d-sphere-advection-lim-continuity} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho \approx - wD[ \rho \boldsymbol{u}] . \label{eq:2d-sphere-advection-lim-discrete-continuity} \end{equation}$$

Tracers

For the tracer concentration per unit mass q, the tracer density (scalar) \rho q follows the advection/transport equation

$$\begin{equation} \frac{\partial}{\partial t} \rho q = - \nabla \cdot(\rho q \boldsymbol{u}) + g(\rho, q). \label{eq:2d-sphere-advection-lim-tracers} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho q \approx - wD[ \rho q \boldsymbol{u}] + g(\rho, q), \label{eq:2d-sphere-advection-lim-discrete-tracers} \end{equation}$$

where g(\rho, q) = - \nu_4 [\nabla^4_h (\rho q)] represents the horizontal hyperdiffusion operator, with \nu_4 (measured in m^4/s) the hyperviscosity constant coefficient.

Currently tracers are only treated explicitly in the time discretization.

Prognostic variables

• \rho: density measured in kg/m³.
• \boldsymbol{u} velocity, a vector measured in m/s. Since this is a 2D problem, \boldsymbol{u} \equiv \boldsymbol{u}_h.
• \rho q: the tracer density scalar, where q is the tracer concentration per unit mass.

Differentiation operators

Because this is a purely 2D problem, there is no staggered vertical discretization, hence, there is no need of specifying variables at cell centers, faces or to reconstruct from faces to centers and vice versa.

• wD is the discrete horizontal weak spectral divergence, called wdiv in the example code.
• G is the discrete horizontal strong spectral gradient, called grad in the example code.

To discretize the hyperdiffusion operator, g(\rho, q) = - \nu_4 [\nabla^4 (\rho q)], in the horizontal direction, we compose the horizontal weak divergence, wD, and the horizontal gradient operator, G_h, twice, with an intermediate call to weighted_dss! between the two compositions, as in [g_2(\rho, g) \circ DSS(\rho, q) \circ g_1(\rho, q)], with:

• g_1(\rho, q) = wD(G_h(q))
• DSS(\rho, q) = DSS(g_1(\rho q))
• g_2(\rho, q) = -\nu_4 wD(\rho G_h(\rho q))
  • with \nu_4 the hyperviscosity coefficient.

Problem flow and set-up

This test case is set up in a Cartesian planar domain [-2 \pi, 2 \pi]^2, doubly periodic.

The flow was chosen to be a horizontal uniform rotation. Moreover, the flow is reversed halfway through the time period so that the tracer blobs go back to its initial configuration (using the same speed scaling constant which was derived to account for the distance travelled in all directions in a half period).

$$\begin{align} u &= k sin (\lambda)^2 sin (2 \phi) \cos(\pi t / T_f) + \frac{2 \pi}{T_f} \cos (\phi) \nonumber \\\ v &= k \sin (2 \lambda) \cos (\phi) \cos(\pi t / T_f) \label{eq:2d-sphere-lim-flow} \end{align}$$

where u_0 = 2 \pi R / T_f is the speed scaling factor to have the flow reversed halfway through the time period, T_f = 86400 * 12 (i.e., 12 days in seconds) is the final simulation time, which coincides with the temporal period to have a full rotation around the sphere of radius R.

This example is set up to run with three possible initial conditions:

• cosine_bells
• gaussian_bells
• cylinders: two 2D slotted cylinders (test case available in the literature, cfr: GubaOpt2014).

Application of Flux Limiters

Because this is a fully 2D problem, the application of limiters does not affect the order of operations, which is implemented as follows:

1. Horizontal trasport with hyperdiffusion (with weak divergence wD)
2. Horizontal flux limiters
3. DSS

Shallow-water equations

The shallow water equations in the so-called vector invariant form from Bao2014 are:

$$\begin{align} \frac{\partial h}{\partial t} + \nabla \cdot (h u) &= 0\\\ \frac{\partial u_i}{\partial t} + \nabla (\Phi + \tfrac{1}{2}\|u\|^2)_i &= (\boldsymbol{u} \times (f + \nabla \times \boldsymbol{u}))_i \label{eq:shallow-water} \end{align}$$

where f is the Coriolis term and \Phi = g(h+h_s), with g the gravitational accelration constant, h the (free) height of the fluid and h_s a non-uniform reference surface.

To the above set of equations, we allow the uset to add a hyperdiffusion operator, g(h, \boldsymbol{u}) = - \nu_4 [\nabla^4 (h, \boldsymbol{u})], with \nu_4 (measured in m^4/s) the hyperviscosity constant coefficient. In the hyperdiffusion expression, \nabla^4 represents a biharmonic operator, and it assumes a different formulation on curvilinear reference systems, depending on it being applied to a scalar field, such as h, or a vector field, such as \boldsymbol{u}.

The governing equations then become:

$$\begin{align} \frac{\partial h}{\partial t} + \nabla \cdot (h u) &= g(h, \boldsymbol{u})\\\ \frac{\partial u_i}{\partial t} + \nabla (\Phi + \tfrac{1}{2}\|u\|^2)_i &= (\boldsymbol{u} \times (f + \nabla \times \boldsymbol{u}))_i + g(h, \boldsymbol{u}) \label{eq:shallow-water-with-hyperdiff} \end{align}$$

Since this is a 2D problem (with related 2D vector field), the curl is defined to be

$$\begin{align} \omega^i = (\nabla \times u)^i &= \begin{cases} 0 &\text{ if \$i =1,2\$},\\\ \frac{1}{J} \left[ \frac{\partial u_2}{\partial \xi^1} - \frac{\partial u_1}{\partial \xi^2} \right] &\text{ if \$i=3\$}, \end{cases} \label{eq:2Dvorticity} \end{align}$$

where we have used the coordinate system in each 2D reference element, i.e., (\xi^1, \xi^2) \in [-1,1]\times[-1,1]. Similarly, if additionally v^1 = v^2 = 0, then

$$\begin{align} (\boldsymbol{u} \times \boldsymbol{v})_i = \begin{cases} J u^2 v^3 &\text{ if \$i=1\$},\\\ - J u^1 v^3 &\text{ if \$i=2\$},\\\ 0 &\text{ if \$i=3\$}. \end{cases} \end{align}$$

Hence, we can rewrite equations \eqref{eq:shallow-water} using the velocity representation in covariant coordinates, in this case u = u_1 \boldsymbol{b}^1 + u_2 \boldsymbol{b}^2 + 0\boldsymbol{b}^3, and g(h, \boldsymbol{u}) = 0 for simplicity, as:

$$\begin{align} \frac{\partial h}{\partial t} + \frac{1}{J}\frac{\partial}{\partial \xi^j}\Big(h J u^j\Big) &= 0\\\ \frac{\partial u_i}{\partial t} + \frac{\partial}{\partial \xi^i} (\Phi + \tfrac{1}{2}\|u\|^2) &= E_{ijk}u^j (f^k + \omega^k) . \label{eq:covariant-shallow-water} \end{align}$$

Prognostic variables

• h: scalar height field of the fluid, measured in m.
• \boldsymbol{u} velocity, a 2D vector measured in m/s.

Differentiation operators

Because this is a purely 2D problem, there is no staggered vertical discretization, hence, there is no need of specifying variables at cell centers, faces or to reconstruct from faces to centers and vice versa.

• D is the discrete horizontal strong spectral divergence, called div in the example code.
• wD is the discrete horizontal weak spectral divergence, called wdiv in the example code.
• G is the discrete horizontal strong spectral gradient, called grad in the example code.
• wG is the discrete horizontal weak spectral gradient, called wgrad in the example code.
• Curl is the discrete curl, called curl in the example code.
• wCurl is the discrete weak curl, called wcurl in the example code.

To discretize the hyperdiffusion operator, g(h, \boldsymbol{u}) = - \nu_4 [\nabla^4 (h, \boldsymbol{u})], in the horizontal direction, we compose the weak divergence, wD, and the gradient operator, G, twice, with an intermediate call to weighted_dss! between the two compositions, as in [g_2(h, \boldsymbol{u}) \circ DSS(h, \boldsymbol{u}) \circ g_1(h, \boldsymbol{u})]. Moreover, when g(h, \boldsymbol{u}) = - \nu_4 [\nabla^4 (h)], i.e., the operator is applied to a scalar field only, it is discretized composing the following operations:

• g_1(h) = wD(G(h))
• DSS(g_1(h))
• g_2(h) = -\nu_4 wD(G(h)) whereas, when the operator is applied to a vector field, i.e., g(h, \boldsymbol{u}) = - \nu_4 [\nabla^4 (\boldsymbol{u})], it is discretized as:
• g_1(h, \boldsymbol{u}) = wG(D(\boldsymbol{u})) - wCurl(Curl(\boldsymbol{u}))
• DSS(h, \boldsymbol{u}) = DSS(g_1(h, \boldsymbol{u}))
• g_2(h, \boldsymbol{u}) = -\nu_4 \left[ wG(D(\boldsymbol{u})) - wCurl(Curl(\boldsymbol{u})) \right]

In both cases, \nu_4 is the hyperviscosity coefficient.

Problem flow and set-up

This test case is set up on a 2D (surface) spherical domain represented by a cubed-sphere manifold.

This suite of examples contains five different test cases:

• One, invoked via the command-line argument steady_state, which reproduces Test Case 2 in Williamson1992. This test case gives the steady-state solution to the non-linear shallow water equations. It consists of a solid body rotation or zonal flow with the corresponding geostrophic height field. The Coriolis parameter is a function of latitude and longitude so the flow can be specified with the spherical coordinate poles not necessarily coincident with Earth's rotation axis. Hence, this test case can be run with a specified command-line argument for the angle \alpha that represents the angle between the north pole and the center of the top cube panel of the cubed-sphere geometry.
• A second one, invoked via the command-line argument steady_state_compact, reproduces Test Case 3 in Williamson1992. This test case gives the steady-state solution to the non-linear shallow water equations with nonlinear zonal geostrophic flow with compact support.
• A third one, invoked via the command-line argument mountain, reproduces Test Case 5 in Williamson1992. It represents a zonal flow over an isolated mountain, where the governing equations describe a global steady-state nonlinear zonal geostrophic flow, with a corresponding geostrophic height field over a non-uniform reference surface h_s.
• A fourth one, invoked via the command-line argument rossby_haurwitz, reproduces Test Case 6 in Williamson1992. It represents the solution of the nonlinear barotropic vorticity equation on the sphere.
• A fifth one, invoked via the command-line argument barotropic_instability, reproduces the test case in Galewsky2004 (also in Sec. 7.6 in Ullrich2010). This test case consists of a zonal jet with compact support at a latitude of 45°. A small height disturbance is then added, which causes the jet to become unstable and collapse into a highly vortical structure.

3D Sphere examples

Deformation Flow with Flux Limiters

The 3D sphere advection/transport example in examples/hybrid/sphere/deformation_flow.jl demonstrates the application of flux limiters in the horziontal direction, namely QuasiMonotoneLimiter, in a hybrid 3D spherical domain. It also demonstrates the usage of the flux-corrected transport in the vertical direction; by default, it uses FCTZalesak.

Equations and discretizations

The original test case (without limiters or flux-corrected transport) is specified in Section 1.1 of Ullrich2012DynamicalCM.

Mass

Follows the continuity equation

$$\begin{equation} \frac{\partial}{\partial t} \rho = - \nabla \cdot(\rho \boldsymbol{u}) . \label{eq:3d-sphere-lim-continuity} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho \approx - D_h[\rho \boldsymbol{u}^c] - D^c_v[I^f(\rho) \boldsymbol{u}^f]. \label{eq:3d-sphere-lim-discrete-continuity} \end{equation}$$

Tracers

This test case has five different tracer concentrations per unit mass q_i, hence five different tracer densities (scalar) \rho q_i. They all follow the same advection/transport equation

$$\begin{equation} \frac{\partial}{\partial t} \rho q = - \nabla \cdot(\rho q \boldsymbol{u}) + g(\rho, q). \label{eq:3d-sphere-lim-tracers} \end{equation}$$

This is discretized using the following

$$\begin{equation} \frac{\partial}{\partial t} \rho q \approx - D_h[ \rho q \boldsymbol{u}^c] - D^c_v\left[I^f(\rho q) * \boldsymbol{u}^f_h + FCT^f\left( \boldsymbol{u}^f_v, \frac{\rho q}{\rho} \right) \right] + g(\rho, q), \label{eq:3d-sphere-lim-discrete-tracers} \end{equation}$$

where g(\rho, q) = - \nu_4 [\nabla^4_h (\rho q)] represents the horizontal hyperdiffusion operator, with \nu_4 (measured in m^4/s) the hyperviscosity constant coefficient.

Currently tracers are only treated explicitly in the time discretization.

Prognostic variables

• \rho: density measured in kg/m³. This is discretized at cell centers.
• \boldsymbol{u} velocity, a vector measured in m/s. This is discretized via \boldsymbol{u} = \boldsymbol{u}_h + \boldsymbol{u}_v where
  • \boldsymbol{u}_h = u_1 \boldsymbol{e}^1 + u_2 \boldsymbol{e}^2 is the projection onto horizontal covariant components (covariance here means with respect to the reference element), stored at cell centers.
  • \boldsymbol{u}_v = u_3 \boldsymbol{e}^3 is the projection onto the vertical covariant components, stored at cell faces.
• \rho q_i: tracer density scalars, where q_i is a tracer concentration per unit mass, are stored at cell centers.

Operators

We make use of the following operators

Reconstructions

• I^c is the face-to-center reconstruction operator, called If2c in the example code.
• I^f is the center-to-face reconstruction operator, called Ic2f in the example code.
  • Currently this is just the arithmetic mean, but we will need to use a weighted version with stretched vertical grids.
• FCT^f denotes either the center-to-face upwind product operator (which represents no flux-corrected transport), the center-to-face Boris & Book FCT operator, or the center-to-face Zalesak FCT operator.

Differentiation operators

• D_h is the discrete horizontal strong spectral divergence, called hdiv in the example code.
• wD_h is the discrete horizontal weak spectral divergence, called hwdiv in the example code.
• D^c_v is the face-to-center vertical divergence, called vdivf2c in the example code.
  • This example uses advective fluxes equal to zero at the top and bottom boundaries.
• G_h is the discrete horizontal spectral gradient, called hgrad in the example code.

To discretize the hyperdiffusion operator for each tracer concentration, g(\rho, q_i) = - \nu_4 [\nabla^4 (\rho q_i)], in the horizontal direction, we compose the horizontal weak divergence, wD_h, and the horizontal gradient operator, G_h, twice, with an intermediate call to weighted_dss! between the two compositions, as in [g_2(\rho, g) \circ DSS(\rho, q) \circ g_1(\rho, q_i)], with:

• g_1(\rho, q_i) = wD_h(G_h(q_i))
• DSS(\rho, q_i) = DSS(g_1(\rho q_i))
• g_2(\rho, q_i) = -\nu_4 wD_h(\rho G_h(\rho q_i))
  • with \nu_4 the hyperviscosity coefficient.

Application of Flux Limiters

!!! note

Since we use flux limiters that limit only operators defined in the spectral space (i.e., they are applied level-wise in the horizontal direction), the application of limiters has to follow a precise order in the sequence of operations that specifies the total tendency.

The order of operations should be the following:

1. Horizontal transport (with strong divergence D_h)
2. Horizontal flux limiters
3. Horizontal hyperdiffusion (with weak divergence wD_h)
4. Vertical transport
5. DSS

Problem flow and set-up

This test case is set up in a 3D (shell) spherical domain where the elevation goes from z=0~\textrm{m} (i.e., from the radius of the sphere R = 6.37122 10^6~\textrm{m}) to z_{\textrm{top}} = 12000~\textrm{m}.

The flow (reversed halfway through the time period) is specified as \boldsymbol{u} = \boldsymbol{u}_a + \boldsymbol{u}_d, where the components are defined as follows:

$$\begin{align} u_a &= k \sin (\lambda')^2 \sin (2 \phi) \cos(\pi t / \tau) + \frac{2 \pi R}{\tau} \cos (\phi) \nonumber \\\ v_a &= k \sin (2 \lambda') \cos (\phi) \cos(\pi t / \tau) \nonumber \\\ u_d &= \frac{\omega_0 R}{ b / p_{\textrm{top}}} \cos (\lambda') \cos(\phi)^2 \cos(2 \pi t / \tau) \left[-exp \left( \frac{(p - p_0)}{ b p_{\textrm{top}}} \right) + exp \left( \frac{(p_{\textrm{top}} - p(zc))}{b p_{\textrm{top}}} \right) \right] \nonumber \label{eq:3d-sphere-lim-flow} \end{align}$$

where all values of the parameters can be found in Table 1.1 in the reference Ullrich2012DynamicalCM.