add EC formulation

examples/elixir_shallowwater_covariant_cubed_sphere_EC.jl

@@ -0,0 +1,75 @@
+###############################################################################
+# DGSEM for the shallow water equations on the cubed sphere
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Parameters
+
+initial_condition = initial_condition_geostrophic_balance
+polydeg = 3
+cells_per_dimension = 5
+n_saves = 10
+
+###############################################################################
+# Spatial discretization
+
+tspan = (0.0, 1.0 * SECONDS_PER_DAY)
+
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = polydeg,
+                              initial_refinement_level = 0,
+                              element_local_mapping = true)
+
+equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION,
+                                             EARTH_ROTATION_RATE,
+                                             global_coordinate_system = GlobalSphericalCoordinates())
+
+# Create DG solver with polynomial degree = p
+volume_flux = (flux_split_covariant, flux_nonconservative_split_covariant)
+surface_flux = (flux_split_covariant, flux_nonconservative_split_covariant)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+initial_condition_transformed = transform_initial_condition(initial_condition, equations)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
+                                    source_terms = source_terms_split_covariant)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0 to T
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 50,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(dt = (tspan[2] - tspan[1]) / n_saves,
+                                     solution_variables = cons2cons)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 0.4)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 100.0, save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

src/TrixiAtmo.jl

@@ -37,7 +37,8 @@ export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D,
 export GlobalCartesianCoordinates, GlobalSphericalCoordinates
 
 export flux_chandrashekar, flux_LMARS, flux_split_covariant, flux_nonconservative_weak_form,
-       flux_nonconservative_split_covariant
+       flux_nonconservative_split_covariant, flux_split_covariant_lax_friedrichs,
+       source_terms_split_covariant
 
 export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal,
        lake_at_rest_error, source_terms_lagrange_multiplier, source_terms_weak_form,

src/equations/covariant_shallow_water.jl

@@ -1,10 +1,14 @@
 @muladd begin
 #! format: noindent
 
+"""
+    CovariantShallowWaterEquations2D{GlobalCoordinateSystem} <:  
+        AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3}
+"""
 struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} <:
        AbstractCovariantEquations{2, 3, GlobalCoordinateSystem, 3}
-    gravity::RealT
-    rotation_rate::RealT
+    gravity::RealT  # acceleration due to gravity
+    rotation_rate::RealT  # rotation rate for Coriolis term 
     global_coordinate_system::GlobalCoordinateSystem
     function CovariantShallowWaterEquations2D(gravity::RealT,
                                               rotation_rate::RealT;
@@ -14,35 +18,37 @@ struct CovariantShallowWaterEquations2D{GlobalCoordinateSystem, RealT <: Real} <
     end
 end
 
+# Our implementation of flux-differencing formulation uses nonconservative terms, but the 
+# standard weak form does not. To handle both options, we have defined a dummy kernel for 
+# the nonconservative terms that does nothing when VolumeIntegralWeakForm is used with a 
+# nonconservative system.
+Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True()
+
 # The conservative variables are the height and contravariant momentum components
 function Trixi.varnames(::typeof(cons2cons), ::CovariantShallowWaterEquations2D)
     return ("h", "h_vcon1", "h_vcon2")
 end
 
-# Convenience function to extract the velocity
-function velocity(u, ::CovariantShallowWaterEquations2D)
-    return SVector(u[2] / u[1], u[3] / u[1])
+# The primitive variables are the height and contravariant velocity components
+function Trixi.varnames(::typeof(cons2prim), ::CovariantShallowWaterEquations2D)
+    return ("h", "vcon1", "vcon2")
 end
 
-# Convenience function to extract the momentum
-function momentum(u, ::CovariantShallowWaterEquations2D)
-    return SVector(u[2], u[3])
+# The change of variables contravariant2global converts the two local contravariant vector 
+# components u[2] and u[3] to the three global vector components specified by 
+# equations.global_coordinate_system (e.g. spherical or Cartesian). This transformation 
+# works for both primitive and conservative variables, although varnames refers 
+# specifically to transformations from conservative variables.
+function Trixi.varnames(::typeof(contravariant2global),
+                        ::CovariantShallowWaterEquations2D)
+    return ("h", "h_vglo1", "h_vglo2", "h_vglo3")
 end
 
-# Our implementation of flux-differencing formulation uses nonconservative terms, but the 
-# standard weak form does not. To handle both options, we have defined a dummy kernel for 
-# the nonconservative terms that does nothing when VolumeIntegralWeakForm is used with a 
-# nonconservative system.
-Trixi.have_nonconservative_terms(::CovariantShallowWaterEquations2D) = True()
-
-# Entropy function (total energy per unit volume)
-@inline function Trixi.entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D)
-    h, h_vcon1, h_vcon2 = u
-    Gcov = metric_covariant(aux_vars, equations)
-    vcon = SVector(h_vcon1 / h, h_vcon2 / h)
-    vcov = Gcov * vcon
-    return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2)
-end
+# Convenience functions to extract variables from state vector
+@inline waterheight(u, ::CovariantShallowWaterEquations2D) = u[1]
+@inline velocity(u, ::CovariantShallowWaterEquations2D) = SVector(u[2] / u[1],
+                                                                  u[3] / u[1])
+@inline momentum(u, ::CovariantShallowWaterEquations2D) = SVector(u[2], u[3])
 
 @inline function Trixi.cons2prim(u, aux_vars, ::CovariantShallowWaterEquations2D)
     h, h_vcon1, h_vcon2 = u
@@ -54,21 +60,12 @@ end
     return SVector(h, h * vcon1, h * vcon2)
 end
 
-# Entropy variables (partial derivatives of entropy with respect to conservative variables)
 @inline function Trixi.cons2entropy(u, aux_vars,
                                     equations::CovariantShallowWaterEquations2D)
-    h, h_vcon1, h_vcon2 = u
-    Gcov = metric_covariant(aux_vars, equations)
-    vcon = SVector(h_vcon1 / h, h_vcon2 / h)
-    vcov = Gcov * vcon
-    w1 = equations.gravity * h - 0.5f0 * dot(vcov, vcon)
-    return SVector{3}(w1, vcov[1], vcov[2])
-end
-
-# Height and three global momentum components
-function Trixi.varnames(::typeof(contravariant2global),
-                        ::CovariantShallowWaterEquations2D)
-    return ("h", "h_vglo1", "h_vglo2", "h_vglo3")
+    h = waterheight(u, equations)
+    vcon = velocity(u, equations)
+    vcov = metric_covariant(aux_vars, equations) * vcon
+    return SVector{3}(equations.gravity * h - 0.5f0 * dot(vcov, vcon), vcov[1], vcov[2])
 end
 
 # Convert contravariant momentum components to the global coordinate system
@@ -84,6 +81,15 @@ end
     vcon1, vcon2 = basis_contravariant(aux_vars, equations) * SVector(u[2], u[3], u[4])
     return SVector(u[1], vcon1, vcon2)
 end
+
+# Entropy function (total energy per unit volume)
+@inline function Trixi.entropy(u, aux_vars, equations::CovariantShallowWaterEquations2D)
+    h = waterheight(u, equations)
+    vcon = velocity(u, equations)
+    vcov = metric_covariant(aux_vars, equations) * vcon
+    return 0.5f0 * (dot(vcov, vcon) + equations.gravity * h^2)
+end
+
 # The flux for the covariant form takes in the element container and node/element indices
 # in order to give the flux access to the geometric information
 @inline function Trixi.flux(u, aux_vars, orientation::Integer,
@@ -171,16 +177,15 @@ end
 # Geometric and Coriolis sources for a rotating sphere with VolumeIntegralWeakForm
 @inline function source_terms_weak_form(u, x, t, aux_vars,
                                         equations::CovariantShallowWaterEquations2D)
-    h, h_vcon1, h_vcon2 = u
-
     # Geometric variables
     Gcon = metric_contravariant(aux_vars, equations)
     (Gamma1, Gamma2) = christoffel_symbols(aux_vars, equations)
     J = area_element(aux_vars, equations)
 
     # Physical variables
-    h_vcon = SVector{2}(h_vcon1, h_vcon2)
-    v_con = SVector{2}(h_vcon1 / h, h_vcon2 / h)
+    h = waterheight(u, equations)
+    h_vcon = momentum(u, equations)
+    v_con = velocity(u, equations)
 
     # Doubly-contravariant flux tensor
     T = h_vcon * v_con' + 0.5f0 * equations.gravity * h^2 * Gcon
@@ -192,8 +197,8 @@ end
     s_geo = SVector(sum(Gamma1 .* T), sum(Gamma2 .* T))
 
     # Combined source terms
-    source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon1 - Gcon[1, 1] * h_vcon2)
-    source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon1 - Gcon[2, 1] * h_vcon2)
+    source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon[1] - Gcon[1, 1] * h_vcon[2])
+    source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon[1] - Gcon[2, 1] * h_vcon[2])
 
     # Do not scale by Jacobian since apply_jacobian! is called before this
     return SVector(zero(eltype(u)), -source_1, -source_2)
@@ -202,7 +207,6 @@ end
 # Geometric and Coriolis sources for a rotating sphere with VolumeIntegralFluxDifferencing
 @inline function source_terms_split_covariant(u, x, t, aux_vars,
                                               equations::CovariantShallowWaterEquations2D)
-    h, h_vcon1, h_vcon2 = u
 
     # Geometric variables
     Gcov = metric_covariant(aux_vars, equations)
@@ -211,8 +215,8 @@ end
     J = area_element(aux_vars, equations)
 
     # Physical variables
-    h_vcon = SVector{2}(h_vcon1, h_vcon2)
-    v_con = SVector{2}(h_vcon1 / h, h_vcon2 / h)
+    h_vcon = momentum(u, equations)
+    v_con = velocity(u, equations)
 
     # Doubly-contravariant and mixed inertial flux tensors
     h_vcon_vcon = h_vcon * v_con'
@@ -226,8 +230,8 @@ end
              Gcon * (Gamma1 * h_vcov_vcon[1, :] + Gamma2 * h_vcov_vcon[2, :]))
 
     # Combined source terms
-    source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon1 - Gcon[1, 1] * h_vcon2)
-    source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon1 - Gcon[2, 1] * h_vcon2)
+    source_1 = s_geo[1] + f * J * (Gcon[1, 2] * h_vcon[1] - Gcon[1, 1] * h_vcon[2])
+    source_2 = s_geo[2] + f * J * (Gcon[2, 2] * h_vcon[1] - Gcon[2, 1] * h_vcon[2])
 
     # Do not scale by Jacobian since apply_jacobian! is called before this
     return SVector(zero(eltype(u)), -source_1, -source_2)
@@ -246,23 +250,19 @@ end
     Gcon_ll = metric_contravariant(aux_vars_ll, equations)
     Gcon_rr = metric_contravariant(aux_vars_rr, equations)
 
-    phi_ll = max(h_ll * equations.gravity, 0)
-    phi_rr = max(h_rr * equations.gravity, 0)
-
     return max(abs(h_vcon_ll[orientation] / h_ll) +
-               sqrt(Gcon_ll[orientation, orientation] * phi_ll),
+               sqrt(Gcon_ll[orientation, orientation] * h_ll * equations.gravity),
                abs(h_vcon_rr[orientation] / h_rr) +
-               sqrt(Gcon_rr[orientation, orientation] * phi_rr))
+               sqrt(Gcon_rr[orientation, orientation] * h_rr * equations.gravity))
 end
 
 # Maximum wave speeds with respect to the covariant basis
 @inline function Trixi.max_abs_speeds(u, aux_vars,
                                       equations::CovariantShallowWaterEquations2D)
     h, h_vcon1, h_vcon2 = u
     Gcon = metric_contravariant(aux_vars, equations)
-    phi = max(h * equations.gravity, 0)
-    return abs(h_vcon1 / h) + sqrt(Gcon[1, 1] * phi),
-           abs(h_vcon2 / h) + sqrt(Gcon[2, 2] * phi)
+    return abs(h_vcon1 / h) + sqrt(Gcon[1, 1] * h * equations.gravity),
+           abs(h_vcon2 / h) + sqrt(Gcon[2, 2] * h * equations.gravity)
 end
 
 # Rossby-Haurwitz wave

src/equations/equations.jl

@@ -184,7 +184,7 @@ end
             "basis_contravariant[1,2]",     # e₁ ⋅ a²
             "basis_contravariant[2,2]",     # e₂ ⋅ a²
             "basis_contravariant[3,2]",     # e₃ ⋅ a²
-            "area_element",                 # √G = √(G₁₁G₂₂ - G₁₂G₂₁) = ||a₁ × a₂||
+            "area_element",                 # J = √(G₁₁G₂₂ - G₁₂G₂₁) = ||a₁ × a₂||
             "metric_covariant[1,1]",        # G₁₁
             "metric_covariant[1,2]",        # G₁₂ = G₂₁
             "metric_covariant[2,2]",        # G₂₂
@@ -216,7 +216,7 @@ end
                          aux_vars[11], aux_vars[12])
 end
 
-# Extract the area element √G = (det(AᵀA))^(1/2) from the auxiliary variables
+# Extract the area element J = (det(AᵀA))^(1/2) from the auxiliary variables
 @inline function area_element(aux_vars, ::AbstractCovariantEquations{2})
     return aux_vars[13]
 end
@@ -274,10 +274,9 @@ end
                                                               aux_vars_rr,
                                                               orientation_or_normal_direction,
                                                               equations::AbstractCovariantEquations)
-    sqrtG = area_element(aux_vars_ll, equations)
     λ = dissipation.max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr,
                                   orientation_or_normal_direction, equations)
-    return -0.5f0 * sqrtG * λ * (u_rr - u_ll)
+    return -0.5f0 * area_element(aux_vars_ll, equations) * λ * (u_rr - u_ll)
 end
 
 # Define the two-point nonconservative flux for weak form to be a no-op

src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl

@@ -123,6 +123,56 @@ end
     Trixi.weak_form_kernel!(du, u, element, mesh, False(), equations, dg, cache, alpha)
 end
 
+# Non-conservative flux differencing kernel which uses contravariant flux components, 
+# passing the geometric information contained in the auxiliary variables to the flux 
+# function
+@inline function Trixi.flux_differencing_kernel!(du, u, element, mesh::P4estMesh{2},
+                                                 nonconservative_terms::True,
+                                                 equations::AbstractCovariantEquations{2},
+                                                 volume_flux, dg::DGSEM, cache,
+                                                 alpha = true)
+    (; derivative_split) = dg.basis
+    (; aux_node_vars) = cache.auxiliary_variables
+    symmetric_flux, nonconservative_flux = volume_flux
+
+    # Apply the symmetric flux as usual
+    Trixi.flux_differencing_kernel!(du, u, element, mesh, False(), equations,
+                                    symmetric_flux, dg, cache, alpha)
+
+    for j in eachnode(dg), i in eachnode(dg)
+        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)
+
+        # ξ¹ direction
+        integral_contribution = zero(u_node)
+        for ii in eachnode(dg)
+            u_node_ii = Trixi.get_node_vars(u, equations, dg, ii, j, element)
+            aux_node_ii = get_node_aux_vars(aux_node_vars, equations, dg, ii, j,
+                                            element)
+            flux1 = nonconservative_flux(u_node, u_node_ii, aux_node, aux_node_ii, 1,
+                                         equations)
+            integral_contribution = integral_contribution +
+                                    derivative_split[i, ii] * flux1
+        end
+
+        # ξ² direction
+        for jj in eachnode(dg)
+            u_node_jj = Trixi.get_node_vars(u, equations, dg, i, jj, element)
+            aux_node_jj = get_node_aux_vars(aux_node_vars, equations, dg, i, jj,
+                                            element)
+            flux2 = nonconservative_flux(u_node, u_node_jj, aux_node, aux_node_jj, 2,
+                                         equations)
+            integral_contribution = integral_contribution +
+                                    derivative_split[j, jj] * flux2
+        end
+
+        Trixi.multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution,
+                                         equations, dg, i, j, element)
+    end
+
+    return nothing
+end
+
 # Pointwise interface flux in local coordinates for problems without nonconservative terms
 @inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2},
                                             nonconservative_terms::False,