add more comments and docs, fix notation

Author: tristanmontoya

src/callbacks_step/save_solution_covariant.jl

@@ -43,13 +43,13 @@ end
     u_node = Trixi.get_node_vars(u, equations, dg, i, j, element)
     aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
     relative_vorticity = calc_vorticity_node(u, equations, dg, cache, i, j, element)
-    b = bottom_topography(aux_node, equations)
+    h_s = bottom_topography(aux_node, equations)
     primitive_global = contravariant2global(cons2prim(u_node, aux_node, equations),
                                             aux_node, equations)
-    return SVector(primitive_global..., b, relative_vorticity)
+    return SVector(primitive_global..., h_s, relative_vorticity)
 end
 
-# Calculate relative vorticity ζ = ∂₁v₂ - ∂₂v₁ for equations in covariant form
+# Calculate relative vorticity ζ = (∂₁v₂ - ∂₂v₁)/J for equations in covariant form
 @inline function calc_vorticity_node(u, equations::AbstractCovariantEquations{2},
                                      dg::DGSEM, cache, i, j, element)
     (; derivative_matrix) = dg.basis

src/equations/covariant_advection.jl

@@ -51,7 +51,7 @@ end
 # The conservative variables are the scalar conserved quantity and two contravariant 
 # velocity components.
 function Trixi.varnames(::typeof(cons2cons), ::CovariantLinearAdvectionEquation2D)
-    return ("scalar", "vcon1", "vcon2")
+    return ("h", "vcon1", "vcon2")
 end
 
 # Convenience function to extract the velocity
@@ -62,9 +62,9 @@ end
 # Convert contravariant velocity components to the global coordinate system
 @inline function contravariant2global(u, aux_vars,
                                       equations::CovariantLinearAdvectionEquation2D)
-    vglo1, vglo2, vglo3 = basis_covariant(aux_vars, equations) *
-                          velocity_contravariant(u, equations)
-    return SVector(u[1], vglo1, vglo2, vglo3)
+    v1, v2, v3 = basis_covariant(aux_vars, equations) *
+                 velocity_contravariant(u, equations)
+    return SVector(u[1], v1, v2, v3)
 end
 
 # Convert velocity components in the global coordinate system to contravariant components
@@ -77,7 +77,7 @@ end
 # Scalar conserved quantity and three global velocity components
 function Trixi.varnames(::typeof(contravariant2global),
                         ::CovariantLinearAdvectionEquation2D)
-    return ("scalar", "vglo1", "vglo2", "vglo3")
+    return ("h", "v1", "v2", "v3")
 end
 
 # We will define the "entropy variables" here to just be the scalar variable in the first 

src/equations/covariant_shallow_water.jl

@@ -57,10 +57,10 @@ components $G_{ab}$ as
   Journal of Computational Physics 271:224-243. 
   [DOI: 10.1016/j.jcp.2013.11.033](https://doi.org/10.1016/j.jcp.2013.11.033)
 
-``` note
-    For problems with variable bottom topography and entropy-stable schemes, 
-    [SplitCovariantShallowWaterEquations2D](@ref) should be used.
-```
+!!! note
+    When solving problems with variable bottom topography as well as when using
+    entropy-stable schemes, [SplitCovariantShallowWaterEquations2D](@ref) should be used
+    instead.
 """
 struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} <:
        AbstractCovariantShallowWaterEquations2D{GlobalCoordinateSystem}
@@ -95,7 +95,7 @@ end
 # specifically to transformations from conservative variables.
 function Trixi.varnames(::typeof(contravariant2global),
                         ::AbstractCovariantShallowWaterEquations2D)
-    return ("h", "h_vglo1", "h_vglo2", "h_vglo3")
+    return ("h", "h_v1", "h_v2", "h_v3")
 end
 
 # Convenience functions to extract physical variables from state vector
@@ -112,35 +112,34 @@ end
 @inline function Trixi.cons2prim(u, aux_vars,
                                  equations::AbstractCovariantShallowWaterEquations2D)
     h, h_vcon1, h_vcon2 = u
-    b = bottom_topography(aux_vars, equations)
-    H = h + b
-    return SVector(H, h_vcon1 / h, h_vcon2 / h)
+    h_s = bottom_topography(aux_vars, equations)
+    return SVector(h + h_s, h_vcon1 / h, h_vcon2 / h)
 end
 
 @inline function Trixi.prim2cons(u, aux_vars,
                                  equations::AbstractCovariantShallowWaterEquations2D)
     H, vcon1, vcon2 = u
-    b = bottom_topography(aux_vars, equations)
-    h = H - b
-    return SVector(h, h * vcon1, h * vcon2)
+    h_s = bottom_topography(aux_vars, equations)
+    return SVector(H - h_s, h * vcon1, h * vcon2)
 end
 
+# Entropy variables are w = (g(h+hₛ) - (v₁v¹ + v₂v²)/2, v₁, v₂)ᵀ
 @inline function Trixi.cons2entropy(u, aux_vars,
                                     equations::AbstractCovariantShallowWaterEquations2D)
     h = waterheight(u, equations)
-    b = bottom_topography(aux_vars, equations)
+    h_s = bottom_topography(aux_vars, equations)
     vcon = velocity_contravariant(u, equations)
     vcov = metric_covariant(aux_vars, equations) * vcon
-    return SVector{3}(equations.gravity * (h + b) - 0.5f0 * dot(vcov, vcon), vcov[1],
-                      vcov[2])
+    return SVector{3}(equations.gravity * (h + h_s) - 0.5f0 * dot(vcov, vcon),
+                      vcov[1], vcov[2])
 end
 
 # Convert contravariant momentum components to the global coordinate system
 @inline function contravariant2global(u, aux_vars,
                                       equations::AbstractCovariantShallowWaterEquations2D)
-    h_vglo1, h_vglo2, h_vglo3 = basis_covariant(aux_vars, equations) *
-                                momentum_contravariant(u, equations)
-    return SVector(waterheight(u, equations), h_vglo1, h_vglo2, h_vglo3)
+    h_v1, h_v2, h_v3 = basis_covariant(aux_vars, equations) *
+                       momentum_contravariant(u, equations)
+    return SVector(waterheight(u, equations), h_v1, h_v2, h_v3)
 end
 
 # Convert momentum components in the global coordinate system to contravariant components
@@ -151,15 +150,15 @@ end
     return SVector(u[1], h_vcon1, h_vcon2)
 end
 
-# Entropy function (total energy per unit volume)
+# Entropy function (total energy) given by S = (h(v₁v¹ + v₂v²) + gh² + ghhₛ)/2
 @inline function Trixi.entropy(u, aux_vars,
                                equations::AbstractCovariantShallowWaterEquations2D)
     h = waterheight(u, equations)
-    b = bottom_topography(aux_vars, equations)
+    h_s = bottom_topography(aux_vars, equations)
     vcon = velocity_contravariant(u, equations)
     vcov = metric_covariant(aux_vars, equations) * vcon
     return 0.5f0 * (h * dot(vcov, vcon) + equations.gravity * h^2) +
-           equations.gravity * h * b
+           equations.gravity * h * h_s
 end
 
 # Flux as a function of the state vector u, as well as the auxiliary variables aux_vars, 

src/equations/covariant_shallow_water_split.jl

@@ -111,8 +111,8 @@ end
     Gcon_ll = metric_contravariant(aux_vars_ll, equations)
     Gcov_rr = metric_covariant(aux_vars_rr, equations)
     J_ll = area_element(aux_vars_ll, equations)
-    b_jump = bottom_topography(aux_vars_rr, equations) -
-             bottom_topography(aux_vars_ll, equations)
+    h_s_jump = bottom_topography(aux_vars_rr, equations) -
+               bottom_topography(aux_vars_ll, equations)
 
     # Physical variables
     h_ll = waterheight(u_ll, equations)
@@ -122,7 +122,8 @@ end
     vcov_rr = Gcov_rr * vcon_rr
 
     geometric_term = 0.5f0 * h_vcon_ll[orientation] * (Gcon_ll * vcov_rr - vcon_rr)
-    pressure_term = equations.gravity * Gcon_ll[:, orientation] * h_ll * (h_rr + b_jump)
+    pressure_term = equations.gravity * Gcon_ll[:, orientation] * h_ll *
+                    (h_rr + h_s_jump)
 
     return SVector(zero(eltype(u_ll)), J_ll * (geometric_term[1] + pressure_term[1]),
                    J_ll * (geometric_term[2] + pressure_term[2]))

src/equations/equations.jl

@@ -172,7 +172,8 @@ end
            nvars_christoffel
 end
 
-# Return auxiliary variable names for 2D covariant form
+# Return auxiliary variable names for 2D covariant form. Note that e₁, e₂, e₃ are the
+# global basis vectors, which may be Cartesian or spherical
 @inline function auxvarnames(::AbstractCovariantEquations{2})
     return ("basis_covariant[1,1]",         # e₁ ⋅ a₁
             "basis_covariant[2,1]",         # e₂ ⋅ a₁
@@ -193,7 +194,7 @@ end
             "metric_contravariant[1,1]",    # G¹¹
             "metric_contravariant[1,2]",    # G¹² = G²¹
             "metric_contravariant[2,2]",    # G²²
-            "bottom_topography",            # b
+            "bottom_topography",            # hₛ
             "christoffel_symbols[1][1,1]",  # Γ¹₁₁
             "christoffel_symbols[1][1,2]",  # Γ¹₁₂ = Γ¹₂₁
             "christoffel_symbols[1][2,2]",  # Γ¹₂₂
@@ -236,7 +237,7 @@ end
                          aux_vars[18], aux_vars[19])
 end
 
-# Extract the bottom topography b from the auxiliary variables
+# Extract the bottom topography hₛ from the auxiliary variables
 @inline function bottom_topography(aux_vars, ::AbstractCovariantEquations{2})
     return aux_vars[20]
 end

src/equations/reference_data.jl

@@ -65,7 +65,7 @@ This problem is adapted from Case 1 of the test suite described in the following
 
     # Convert primitive variables from Cartesian coordinates to the chosen global 
     # coordinate system, which depends on the equation type
-    return cartesian2global(SVector(h, vx, vy, vz, 0.0f0), x, equations)
+    return cartesian2global(SVector(h, vx, vy, vz, zero(RealT)), x, equations)
 end
 
 @doc raw"""
@@ -104,7 +104,8 @@ test suite described in the following paper:
 
     # Convert primitive variables from spherical coordinates to the chosen global 
     # coordinate system, which depends on the equation type
-    return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations)
+    return spherical2global(SVector(h, vlon, vlat, zero(RealT), zero(RealT)), x,
+                            equations)
 end
 
 @doc raw"""
@@ -126,12 +127,12 @@ $g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$, $\Omega = 7.292 \times 10^{-5} \ \mathr
 and $h_0 = 8000 \ \mathrm{m}$ and defining the functions 
 ```math
 \begin{aligned}
-A(\theta) &:= \frac{\omega}{2}(2 \Omega+\omega) \cos^2 \theta + 
+A(\theta) &= \frac{\omega}{2}(2 \Omega+\omega) \cos^2 \theta + 
 \frac{1}{4} K^2 \cos^{2 R} \theta\Big((R+1) \cos^2\theta +\left(2 R^2-R-2\right) - 
 \big(2 R^2 / \cos^2 \theta\big) \Big), \\
-B(\theta) &:= \frac{2(\Omega+\omega) K}{(R+1)(R+2)} \cos ^R \theta\big((R^2+2 R+2) - 
+B(\theta) &= \frac{2(\Omega+\omega) K}{(R+1)(R+2)} \cos ^R \theta\big((R^2+2 R+2) - 
 (R+1)^2 \cos^2 \theta\big), \\
-C(\theta) &:=  \frac{1}{4} K^2 \cos^{2 R} \theta\big((R+1) \cos^2 \theta-(R+2)\big),
+C(\theta) &=  \frac{1}{4} K^2 \cos^{2 R} \theta\big((R+1) \cos^2 \theta-(R+2)\big),
 \end{aligned}
 ```
 the initial height field is given by
@@ -147,9 +148,8 @@ This problem corresponds to Case 6 of the test suite described in the following
 """
 @inline function initial_condition_rossby_haurwitz(x, t, equations)
     RealT = eltype(x)
-    a = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
-    lon = atan(x[2], x[1])
-    lat = asin(x[3] / a)
+    a = sqrt(x[1]^2 + x[2]^2 + x[3]^2)  # radius of the sphere
+    lon, lat = atan(x[2], x[1]), asin(x[3] / a)
 
     h_0 = 8.0f3
     omega = 7.848f-6
@@ -176,7 +176,8 @@ This problem corresponds to Case 6 of the test suite described in the following
 
     # Convert primitive variables from spherical coordinates to the chosen global 
     # coordinate system, which depends on the equation type
-    return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations)
+    return spherical2global(SVector(h, vlon, vlat, zero(RealT), zero(RealT)), x,
+                            equations)
 end
 
 @doc raw"""
@@ -222,7 +223,8 @@ test suite described in the following paper:
 
     # Convert primitive variables from spherical coordinates to the chosen global 
     # coordinate system, which depends on the equation type
-    return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations)
+    return spherical2global(SVector(h, vlon, vlat, zero(RealT), zero(RealT)), x,
+                            equations)
 end
 
 # Bottom topography function to pass as auxiliary_field keyword argument in constructor for 
@@ -279,12 +281,12 @@ H(\vec{x},t) =
 \vec{\varphi}(\vec{c},t) \cdot \vec{x}/\lVert \vec{x} \rVert \big)^2 +
 (\vec{\Omega} \cdot \vec{x})^2 + 2k_1\right),
 ```
-where we take $v_0 = 2\pi a / (12 \ \mathrm{days})$ and 
-$k_1 = 133681 \ \mathrm{m}^2/\mathrm{s}^2$, and we use $\lVert \cdot \rVert$ to denote the 
-Euclidean norm. To use this test case with [`SplitCovariantShallowWaterEquations2D`](@ref), 
-the keyword argument `auxiliary_field = bottom_topography_unsteady_solid_body_rotation` 
-should be passed into the `SemidiscretizationHyperbolic` constructor. This analytical 
-solution was derived in the following paper:
+where $v_0 = 2\pi a / (12 \ \mathrm{days})$, $k_1 = 133681 \ \mathrm{m}^2/\mathrm{s}^2$,
+and $\lVert \cdot \rVert$ denotes the Euclidean norm. To use this test case with
+[`SplitCovariantShallowWaterEquations2D`](@ref), the keyword argument
+`auxiliary_field = bottom_topography_unsteady_solid_body_rotation` should be passed into
+the `SemidiscretizationHyperbolic` constructor. This analytical solution was derived in the
+following paper:
 - M. Läuter, D. Handorf, and K. Dethloff (2005). Unsteady analytical solutions of the 
   spherical shallow water equations. Journal of Computational Physics 210:535–553.
   [DOI: 10.1016/j.jcp.2005.04.022](https://doi.org/10.1016/j.jcp.2005.04.022)
@@ -320,7 +322,7 @@ solution was derived in the following paper:
 
     # Convert primitive variables from Cartesian coordinates to the chosen global 
     # coordinate system, which depends on the equation type
-    return cartesian2global(SVector(H, v[1], v[2], v[3]), x, equations)
+    return cartesian2global(SVector(H, v[1], v[2], v[3], zero(RealT)), x, equations)
 end
 
 # Bottom topography function to pass as auxiliary_field keyword argument in constructor for 
@@ -332,8 +334,34 @@ end
 @doc raw"""
     initial_condition_barotropic_instability(x, t, equations)
 
-Barotrotropic instability initiated by a perturbation applied to a mid-latitude jet. A full
-description of this case is provided in the following paper:
+Barotrotropic instability initiated by a perturbation applied to a mid-latitude jet. The  velocity field is a purely zonal flow, given as a function of the latitude $\theta$ as
+```math
+v_\lambda(\theta) = \begin{cases}
+u_0 \exp(-4 / (\theta_1 - \theta_0)^{-2})\exp((\theta - \theta_0)^{-1}
+(\theta - \theta_1)^{-1}), & \quad \theta_0 < \theta < \theta_1, \\
+0 & \quad \text{otherwise},
+\end{cases}
+```
+where $u_0 = 80 \ \mathrm{m}/\mathrm{s}$, $\theta_0 = \pi/7$, and 
+$\theta_1 = \pi/2 - \theta_0$. The background geopotential height field is given by 
+```math
+h_0(\theta) = 10158 \ \mathrm{m} - 
+\frac{a}{g} \int_{-\pi/2}^\theta v_\lambda(\theta')\big(2\Omega\sin\theta' + 
+v_\lambda(\theta')\tan\theta' / a \big)\, \mathrm{d}\theta',
+``` 
+where $a = 6.37122 \times 10^3\ \mathrm{m}$ is the Earth's radius, 
+$g = 9.80616 \ \mathrm{m}/\mathrm{s}^2$ is the Earth's gravitational acceleration, and
+$\Omega = 7.292 \times 10^{-5}\ \mathrm{s}^{-1}$ is the Earth's rotation rate. The 
+perturbation is then added to obtain the following geopotential height field:
+```math
+h(\lambda, \theta) = \begin{cases}
+h_0(\theta) + \delta h \cos\theta \exp(-(\lambda/\alpha)^2) 
+\exp(-((\theta_2 -\theta)/\beta)^2), & \quad -\pi < \lambda < \pi,\\
+h_0(\theta), & \quad \text{otherwise},
+\end{cases}
+```
+where $\lambda$ is the longitude coordinate, and we take $\alpha = 1/3$, $\beta = 1/15$, 
+and $\delta h = 120 \ \mathrm{m}$. This problem was proposed in the following paper:
 - J. Galewsky, R. K. Scott, and L. M. Polvani (2004). An initial-value problem for
   testing numerical models of the global shallow-water equations. Tellus A 56.5:429–440.
   [DOI: 10.3402/tellusa.v56i5.14436](https://doi.org/10.3402/tellusa.v56i5.14436)
@@ -366,7 +394,8 @@ description of this case is provided in the following paper:
 
     # Convert primitive variables from spherical coordinates to the chosen global 
     # coordinate system, which depends on the equation type
-    return spherical2global(SVector(h, vlon, vlat, zero(RealT)), x, equations)
+    return spherical2global(SVector(h, vlon, vlat, zero(RealT), zero(RealT)), x,
+                            equations)
 end
 
 @inline function galewsky_velocity(lat, u_0, lat_0, lat_1)