Add 2D covariant form using P4estMesh

.gitignore

@@ -3,4 +3,6 @@ run
 run/*
 out
 out/*
+docs/build
+docs/build/*
 .DS_Store

Project.toml

@@ -23,5 +23,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"

examples/elixir_euler_spherical_advection_cartesian.jl → examples/elixir_spherical_advection_cartesian.jl

@@ -101,8 +101,11 @@ end
 
 initial_condition = initial_condition_advection_sphere
 
+element_local_mapping = false
+
 mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg,
-                                        initial_refinement_level = 0)
+                                        initial_refinement_level = 0,
+                                        element_local_mapping = element_local_mapping)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/elixir_spherical_advection_covariant.jl

@@ -0,0 +1,98 @@
+###############################################################################
+# DGSEM for the linear advection equation on the cubed sphere
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Problem definition
+
+function initial_condition_advection_earth(x, t, ::CovariantLinearAdvectionEquation2D)
+    RealT = eltype(x)
+
+    # set parameters
+    a = EARTH_RADIUS  # radius of the sphere in metres
+    V = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY)  # speed of rotation
+    alpha = convert(RealT, π / 4)  # angle of rotation
+    h_0 = 1000.0f0  # bump height in metres
+    b_0 = 5.0f0 / (a^2)  # bump width
+    lon_0, lat_0 = convert(RealT, 3π / 2), 0.0f0  # initial bump location
+
+    # convert initial position to Cartesian coordinates
+    x_0 = SVector(a * cos(lat_0) * cos(lon_0), a * cos(lat_0) * sin(lon_0), a * sin(lat_0))
+
+    # compute Gaussian bump profile
+    h = h_0 * exp(-b_0 * ((x[1] - x_0[1])^2 + (x[2] - x_0[2])^2 + (x[3] - x_0[3])^2))
+
+    # get zonal and meridional components of the velocity
+    lon, lat = atan(x[2], x[1]), asin(x[3] / a)
+    v_lon = V * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha))
+    v_lat = -V * sin(lon) * sin(alpha)
+
+    return SVector(h, v_lon, v_lat)
+end
+
+###############################################################################
+# Spatial discretization
+
+polydeg = 3
+cells_per_dimension = 2
+
+element_local_mapping = true
+
+mesh = P4estMeshCubedSphere2D(5, EARTH_RADIUS, polydeg = polydeg,
+                              initial_refinement_level = 0,
+                              element_local_mapping = element_local_mapping)
+
+equations = CovariantLinearAdvectionEquation2D()
+
+initial_condition = initial_condition_advection_earth
+
+# Local Lax-Friedrichs surface flux
+surface_flux = flux_lax_friedrichs
+
+# Standard weak-form volume integral
+volume_integral = VolumeIntegralWeakForm()
+
+# Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, 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.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 = 1.0, save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

src/TrixiAtmo.jl

@@ -8,25 +8,29 @@ 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 Reexport: @reexport
-@reexport using StaticArrays: SVector
+using StrideArrays: StrideArray, StaticInt, PtrArray
+using LinearAlgebra: norm, dot
+
+@reexport using StaticArrays: SVector, SMatrix
 
 include("auxiliary/auxiliary.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")
 include("solvers/solvers.jl")
 include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
+include("callbacks_step/callbacks_step.jl")
 
-export CompressibleMoistEulerEquations2D
+export CompressibleMoistEulerEquations2D, CovariantLinearAdvectionEquation2D
 
 export flux_chandrashekar, flux_LMARS
 
-export examples_dir
+export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
+       EARTH_ROTATION_RATE, SECONDS_PER_DAY
 
+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,51 @@
+
+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.elements, 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 = spherical2contravariant(initial_condition(x, t, equations),
+                                              equations, cache.elements, 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

src/callbacks_step/callbacks_step.jl

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

src/equations/covariant_advection.jl

@@ -0,0 +1,99 @@
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3, 3}
+
+A variable-coefficient linear advection equation can be defined on a two-dimensional
+manifold $S \subset \mathbb{R}^3$ as
+```math
+\partial_t h + \nabla_S \cdot (h \vec{v}) = 0.
+```
+We treat this problem as a system of equations in which the first variable is the scalar
+conserved quantity $h$, and the second two are the contravariant components $v^1$ and $v^2$ 
+used in the expansion 
+```math
+\vec{v} = v^1 \vec{a}_1 + v^2 \vec{a}_2,
+```
+which are spatially varying but remain constant in time (i.e. no flux or dissipation is 
+applied to such variables). The resulting system is then given on the reference element as 
+```math
+J \frac{\partial}{\partial t}\left[\begin{array}{c} h \\ v^1 \\ v^2 \end{array}\right] +
+\frac{\partial}{\partial \xi^1} \left[\begin{array}{c} J h v^1 \\ 0 \\ 0 \end{array}\right]
++ 
+\frac{\partial}{\partial \xi^2} \left[\begin{array}{c} J h v^2 \\ 0 \\ 0 \end{array}\right] 
+= \left[\begin{array}{c} 0 \\ 0 \\ 0 \end{array}\right],
+```
+where $J$ is the area element (see the documentation for [`P4estElementContainerCovariant`](@ref)).
+!!! note
+    The initial condition should be prescribed as $[h, u, v]^{\mathrm{T}}$ in terms of the 
+    global velocity components $u$ and $v$ (i.e. zonal and meridional components in the 
+    case of a spherical shell). The transformation to local contravariant components $v^1$ 
+    and $v^2$ is handled internally within `Trixi.compute_coefficients!`.
+"""
+struct CovariantLinearAdvectionEquation2D <: AbstractCovariantEquations{2, 3, 3} end
+
+function Trixi.varnames(::typeof(cons2cons), ::CovariantLinearAdvectionEquation2D)
+    return ("scalar", "v_con_1", "v_con_2")
+end
+
+# Compute the entropy variables (requires element container and indices)
+@inline function Trixi.cons2entropy(u, equations::CovariantLinearAdvectionEquation2D,
+                                    elements, i, j, element)
+    z = zero(eltype(u))
+    return SVector{3}(u[1], z, z)
+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, orientation::Integer,
+                            ::CovariantLinearAdvectionEquation2D,
+                            elements, i, j, element)
+    z = zero(eltype(u))
+    J = volume_element(elements, i, j, element)
+    return SVector(J * u[1] * u[orientation + 1], z, z)
+end
+
+# Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity 
+# components, as they should remain unchanged in time
+@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
+                                                              orientation_or_normal_direction,
+                                                              equations::CovariantLinearAdvectionEquation2D,
+                                                              elements, i_ll, j_ll,
+                                                              i_rr, j_rr, element)
+    z = zero(eltype(u_ll))
+    J = volume_element(elements, i_ll, j_ll, element)
+    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
+                                  equations, elements, i_ll, j_ll, i_rr, j_rr, element)
+    return -0.5f0 * J * λ * SVector(u_rr[1] - u_ll[1], z, z)
+end
+
+# Maximum wave speed with respect to the a specific orientation
+@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
+                                           ::CovariantLinearAdvectionEquation2D,
+                                           elements, i_ll, j_ll, i_rr, j_rr, element)
+    return max(abs(u_ll[orientation + 1]), abs(u_rr[orientation + 1]))
+end
+
+# Maximum wave speeds in each direction for CFL calculation
+@inline function Trixi.max_abs_speeds(u, ::CovariantLinearAdvectionEquation2D,
+                                      elements, i, j, element)
+    return abs(u[2]), abs(u[3])
+end
+
+# Convert contravariant velocity/momentum components to zonal and meridional components
+@inline function contravariant2spherical(u::SVector{3},
+                                         ::CovariantLinearAdvectionEquation2D,
+                                         elements, i, j, element)
+    v_lon, v_lat = contravariant2spherical(u[2], u[3], elements, i, j, element)
+    return SVector(u[1], v_lon, v_lat)
+end
+
+# Convert zonal and meridional velocity/momentum components to contravariant components
+@inline function spherical2contravariant(u::SVector{3},
+                                         ::CovariantLinearAdvectionEquation2D,
+                                         elements, i, j, element)
+    v_con_1, v_con_2 = spherical2contravariant(u[2], u[3], elements, i, j, element)
+    return SVector(u[1], v_con_1, v_con_2)
+end
+end # @muladd

src/equations/equations.jl

@@ -1,5 +1,70 @@
+@muladd begin
+#! format: noindent
+
 using Trixi: AbstractEquations
 
+# Physical constants
+const EARTH_RADIUS = 6.37122  # m
+const EARTH_GRAVITATIONAL_ACCELERATION = 9.80616  # m/s²
+const EARTH_ROTATION_RATE = 7.292e-5  # rad/s
+const SECONDS_PER_DAY = 8.64e4
+
+@doc raw"""
+    AbstractCovariantEquations{NDIMS, 
+                               NDIMS_AMBIENT, 
+                               NVARS} <: AbstractEquations{NDIMS, NVARS} 
+
+Abstract type used to dispatch on systems of equations in covariant form, in which fluxes
+and prognostic variables are stored and computed in terms of their contravariant components 
+defining their expansions in terms of the local covariant tangent basis. The type parameter
+`NDIMS` denotes the dimension of the manifold on which the equations are solved, while
+`NDIMS_AMBIENT` is the dimension of the ambient space in which such a manifold is embedded. 
+
+!!! note 
+    Components of vector-valued fields should be prescibed within the global coordinate 
+    system (i.e. zonal and meridional components in the case of a spherical shell). 
+    When dispatched on this type, the function `Trixi.compute_coefficients!` will 
+    internally use the `covariant_basis` field of the container type 
+    [`P4estElementContainerCovariant`](@ref) to obtain the local contravariant components
+    used in the solver. 
+"""
+abstract type AbstractCovariantEquations{NDIMS,
+                                         NDIMS_AMBIENT,
+                                         NVARS} <: AbstractEquations{NDIMS, NVARS} end
+
+# Numerical flux plus dissipation which passes node/element indices and cache. 
+# We assume that u_ll and u_rr have been transformed into the same local coordinate system.
+@inline function (numflux::Trixi.FluxPlusDissipation)(u_ll, u_rr,
+                                                      orientation_or_normal_direction,
+                                                      equations::AbstractCovariantEquations{2},
+                                                      elements, i_ll, j_ll, i_rr, j_rr,
+                                                      element)
+
+    # The flux and dissipation need to be defined for the specific equation set
+    flux = numflux.numerical_flux(u_ll, u_rr, orientation_or_normal_direction,
+                                  equations,
+                                  elements, i_ll, j_ll, i_rr, j_rr, element)
+    diss = numflux.dissipation(u_ll, u_rr, orientation_or_normal_direction, equations,
+                               elements, i_ll, j_ll, i_rr, j_rr, element)
+    return flux + diss
+end
+
+# Central flux which passes node/element indices and cache. 
+# We assume that u_ll and u_rr have been transformed into the same local coordinate system.
+@inline function Trixi.flux_central(u_ll, u_rr,
+                                    orientation_or_normal_direction,
+                                    equations::AbstractCovariantEquations{2},
+                                    elements, i_ll, j_ll, i_rr, j_rr, element)
+    flux_ll = Trixi.flux(u_ll, orientation_or_normal_direction, equations,
+                         elements, i_ll, j_ll, element)
+    flux_rr = Trixi.flux(u_rr, orientation_or_normal_direction, equations,
+                         elements, i_rr, j_rr, element)
+
+    return 0.5f0 * (flux_ll + flux_rr)
+end
+
+include("covariant_advection.jl")
 abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <:
               AbstractEquations{NDIMS, NVARS} end
 include("compressible_moist_euler_2d_lucas.jl")
+end # @muladd