Add 2D covariant form using P4estMesh

Project.toml

@@ -25,5 +25,5 @@ Static = "0.8, 1"
 StaticArrayInterface = "1.5.1"
 StaticArrays = "1"
 StrideArrays = "0.1.28"
-Trixi = "0.7, 0.8, 0.9"
+Trixi = "0.9"
 julia = "1.9"

docs/Project.toml

@@ -1,7 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
-TrixiAtmo = "c9ed1054-d170-44a9-8ee2-d5566f5d1389"
 
 [compat]
 Documenter = "1.3"

docs/make.jl

@@ -1,7 +1,14 @@
-using TrixiAtmo
 using Documenter
 using DocumenterInterLinks
 
+# Fix for https://github.com/trixi-framework/Trixi.jl/issues/668
+if (get(ENV, "CI", nothing) != "true") &&
+   (get(ENV, "TRIXI_DOC_DEFAULT_ENVIRONMENT", nothing) != "true")
+    push!(LOAD_PATH, dirname(@__DIR__))
+end
+
+using TrixiAtmo
+
 # Provide external links to the Trixi.jl docs (project root and inventory file)
 links = InterLinks("Trixi" => ("https://trixi-framework.github.io/Trixi.jl/stable/",
                                "https://trixi-framework.github.io/Trixi.jl/stable/objects.inv"))

examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl

@@ -45,12 +45,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     v0 = 1.0 # Velocity at the "equator"
     alpha = 0.0 #pi / 4
     v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha))
-    v_lat = -v0 * sin(phi) * sin(alpha)
+    vlat = -v0 * sin(phi) * sin(alpha)
 
     # Transform to Cartesian coordinate system
-    v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long
-    v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long
-    v3 = sin(colat) * v_lat
+    v1 = -cos(colat) * cos(phi) * vlat - sin(phi) * v_long
+    v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long
+    v3 = sin(colat) * vlat
 
     return prim2cons(SVector(rho, v1, v2, v3, 0), equations)
 end

examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl

@@ -47,12 +47,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     v0 = 1.0 # Velocity at the "equator"
     alpha = 0.0 #pi / 4
     v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha))
-    v_lat = -v0 * sin(phi) * sin(alpha)
+    vlat = -v0 * sin(phi) * sin(alpha)
 
     # Transform to Cartesian coordinate system
-    v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long
-    v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long
-    v3 = sin(colat) * v_lat
+    v1 = -cos(colat) * cos(phi) * vlat - sin(phi) * v_long
+    v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long
+    v3 = sin(colat) * vlat
 
     return prim2cons(SVector(rho, v1, v2, v3, 0), equations)
 end

examples/elixir_euler_spherical_advection_cartesian.jl → examples/elixir_shallowwater_cubed_sphere_shell_advection.jl

@@ -6,58 +6,19 @@ using TrixiAtmo
 ###############################################################################
 # semidiscretization of the linear advection equation
 
-# We use the Euler equations structure but modify the rhs! function to convert it to a
-# variable-coefficient advection equation
-equations = CompressibleEulerEquations3D(1.4)
+cells_per_dimension = 5
+initial_condition = initial_condition_gaussian
+
+# We use the ShallowWaterEquations3D equations structure but modify the rhs! function to
+# convert it to a variable-coefficient advection equation
+equations = ShallowWaterEquations3D(gravity_constant = 0.0)
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
-polydeg = 3
-solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
-
-# Initial condition for a Gaussian density profile with constant pressure
-# and the velocity of a rotating solid body
-function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEquations3D)
-    # Gaussian density
-    rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2))
-    # Constant pressure
-    p = 1.0
-
-    # Spherical coordinates for the point x
-    if sign(x[2]) == 0.0
-        signy = 1.0
-    else
-        signy = sign(x[2])
-    end
-    # Co-latitude
-    colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2))
-    # Latitude (auxiliary variable)
-    lat = -colat + 0.5 * pi
-    # Longitude
-    r_xy = sqrt(x[1]^2 + x[2]^2)
-    if r_xy == 0.0
-        phi = pi / 2
-    else
-        phi = signy * acos(x[1] / r_xy)
-    end
-
-    # Compute the velocity of a rotating solid body
-    # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system)
-    v0 = 1.0 # Velocity at the "equator"
-    alpha = 0.0 #pi / 4
-    v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha))
-    v_lat = -v0 * sin(phi) * sin(alpha)
-
-    # Transform to Cartesian coordinate system
-    v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long
-    v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long
-    v3 = sin(colat) * v_lat
-
-    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
-end
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
-# Source term function to transform the Euler equations into the linear advection equations with variable advection velocity
+# Source term function to transform the Euler equations into a linear advection equation with variable advection velocity
 function source_terms_convert_to_linear_advection(u, du, x, t,
-                                                  equations::CompressibleEulerEquations3D,
+                                                  equations::ShallowWaterEquations3D,
                                                   normal_direction)
     v1 = u[2] / u[1]
     v2 = u[3] / u[1]
@@ -66,29 +27,27 @@ function source_terms_convert_to_linear_advection(u, du, x, t,
     s2 = du[1] * v1 - du[2]
     s3 = du[1] * v2 - du[3]
     s4 = du[1] * v3 - du[4]
-    s5 = 0.5 * (s2 * v1 + s3 * v2 + s4 * v3) - du[5]
 
-    return SVector(0.0, s2, s3, s4, s5)
+    return SVector(0.0, s2, s3, s4, 0.0)
 end
 
-initial_condition = initial_condition_advection_sphere
-
-element_local_mapping = false
-
-mesh = P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg,
+# Create a 2D cubed sphere mesh the size of the Earth
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = 3,
                               initial_refinement_level = 0,
-                              element_local_mapping = element_local_mapping)
+                              element_local_mapping = false)
+
+# Convert initial condition given in terms of zonal and meridional velocity components to 
+# one given in terms of Cartesian momentum components
+initial_condition_transformed = spherical2cartesian(initial_condition, equations)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
                                     source_terms = source_terms_convert_to_linear_advection)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-# Create ODE problem with time span from 0.0 to π
-tspan = (0.0, pi)
-ode = semidiscretize(semi, tspan)
+ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY))
 
 # At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
 # and resets the timers
@@ -97,8 +56,7 @@ summary_callback = SummaryCallback()
 # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
 analysis_callback = AnalysisCallback(semi, interval = 10,
                                      save_analysis = true,
-                                     extra_analysis_errors = (:conservation_error,),
-                                     extra_analysis_integrals = (Trixi.density,))
+                                     extra_analysis_errors = (:conservation_error,))
 
 # The SaveSolutionCallback allows to save the solution to a file in regular intervals
 save_solution = SaveSolutionCallback(interval = 10,

examples/elixir_shallowwater_cubed_sphere_shell_standard.jl

@@ -41,12 +41,12 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     v0 = 1.0 # Velocity at the "equator"
     alpha = 0.0 #pi / 4
     v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha))
-    v_lat = -v0 * sin(phi) * sin(alpha)
+    vlat = -v0 * sin(phi) * sin(alpha)
 
     # Transform to Cartesian coordinate system
-    v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long
-    v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long
-    v3 = sin(colat) * v_lat
+    v1 = -cos(colat) * cos(phi) * vlat - sin(phi) * v_long
+    v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long
+    v3 = sin(colat) * vlat
 
     return prim2cons(SVector(rho, v1, v2, v3, 0), equations)
 end

examples/elixir_spherical_advection_covariant.jl

@@ -0,0 +1,65 @@
+###############################################################################
+# DGSEM for the linear advection equation on the cubed sphere
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Spatial discretization
+
+cells_per_dimension = 5
+initial_condition = initial_condition_gaussian
+
+equations = CovariantLinearAdvectionEquation2D()
+
+# Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs,
+               volume_integral = VolumeIntegralWeakForm())
+
+# Create a 2D cubed sphere mesh the size of the Earth
+mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = 3,
+                              initial_refinement_level = 0,
+                              element_local_mapping = true)
+
+# Convert initial condition given in terms of zonal and meridional velocity components to 
+# one given in terms of contravariant velocity components
+initial_condition_transformed = spherical2contravariant(initial_condition, equations)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0 to T
+ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY))
+
+# 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 = 10,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = cons2cons)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+# 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 = 1.0, save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

src/TrixiAtmo.jl

@@ -8,32 +8,38 @@ See also: [trixi-framework/TrixiAtmo.jl](https://github.com/trixi-framework/Trix
 """
 module TrixiAtmo
 
+using Reexport: @reexport
 using Trixi
 using MuladdMacro: @muladd
 using Static: True, False
 using StrideArrays: PtrArray
 using StaticArrayInterface: static_size
-using LinearAlgebra: norm
+using LinearAlgebra: norm, dot
 using Reexport: @reexport
 using LoopVectorization: @turbo
-@reexport using StaticArrays: SVector
+@reexport using StaticArrays: SVector, SMatrix
 
 include("auxiliary/auxiliary.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")
 include("semidiscretization/semidiscretization.jl")
 include("solvers/solvers.jl")
 include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
-include("callbacks_step/stepsize_dg2d.jl")
+include("callbacks_step/callbacks_step.jl")
 
-export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D
+export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D,
+       CovariantLinearAdvectionEquation2D
 
 export flux_chandrashekar, flux_LMARS
 
 export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal,
        lake_at_rest_error, source_terms_lagrange_multiplier,
        clean_solution_lagrange_multiplier!
 export P4estCubedSphere2D, MetricTermsCrossProduct, MetricTermsInvariantCurl
-export examples_dir
+export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
+       EARTH_ROTATION_RATE, SECONDS_PER_DAY
+export spherical2contravariant, contravariant2spherical, spherical2cartesian
+export initial_condition_gaussian
 
+export examples_dir
 end # module TrixiAtmo

src/auxiliary/auxiliary.jl

@@ -7,6 +7,7 @@ read-only and should *not* be modified.
 To find out which files are available, use, e.g., `readdir`:
 
 # Examples
+
 ```@example
 readdir(examples_dir())
 ```

src/callbacks_step/analysis_covariant.jl

@@ -0,0 +1,53 @@
+@muladd begin
+#! format: noindent
+
+function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t,
+                       mesh::P4estMesh{2},
+                       equations::AbstractCovariantEquations{2}, dg::DG, cache)
+    # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ
+    Trixi.integrate_via_indices(u, mesh, equations, dg, cache,
+                                du) do u, i, j, element, equations, dg, du
+        u_node = Trixi.get_node_vars(u, equations, dg, i, j, element)
+        du_node = Trixi.get_node_vars(du, equations, dg, i, j, element)
+        dot(cons2entropy(u_node, equations, cache, (i, j), element), du_node)
+    end
+end
+
+function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
+                                equations::AbstractCovariantEquations{2},
+                                initial_condition, dg::DGSEM, cache, cache_analysis)
+    (; weights) = dg.basis
+    (; node_coordinates) = cache.elements
+
+    # Set up data structures
+    l2_error = zero(func(Trixi.get_node_vars(u, equations, dg, 1, 1, 1), equations))
+    linf_error = copy(l2_error)
+    total_volume = zero(real(mesh))
+
+    # Iterate over all elements for error calculations
+    for element in eachelement(dg, cache)
+
+        # Calculate errors at each volume quadrature node
+        for j in eachnode(dg), i in eachnode(dg)
+            x = Trixi.get_node_coords(node_coordinates, equations, dg, i, j, element)
+
+            u_exact = initial_condition(x, t, equations, cache, (i, j), element)
+
+            u_numerical = Trixi.get_node_vars(u, equations, dg, i, j, element)
+
+            diff = func(u_exact, equations) - func(u_numerical, equations)
+
+            J = volume_element(cache.elements, (i, j), element)
+
+            l2_error += diff .^ 2 * (weights[i] * weights[j] * J)
+            linf_error = @. max(linf_error, abs(diff))
+            total_volume += weights[i] * weights[j] * J
+        end
+    end
+
+    # For L2 error, divide by total volume
+    l2_error = @. sqrt(l2_error / total_volume)
+
+    return l2_error, linf_error
+end
+end # muladd

src/callbacks_step/callbacks_step.jl

@@ -0,0 +1,2 @@
+include("analysis_covariant.jl")
+include("stepsize_dg2d.jl")