covariant weak form with tests

examples/elixir_spherical_advection_cartesian.jl

@@ -1,13 +1,13 @@
-###############################################################################
-# DGSEM for the linear advection equation on the cubed using Lagrange multiplier approach
-###############################################################################
 
 using OrdinaryDiffEq
 using Trixi
+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)
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
@@ -101,7 +101,8 @@ end
 
 initial_condition = initial_condition_advection_sphere
 
-mesh = P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, initial_refinement_level = 0)
+mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg,
+                                        initial_refinement_level = 0)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
@@ -126,24 +127,27 @@ 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,))
 
 # The SaveSolutionCallback allows to save the solution to a file in regular intervals
 save_solution = SaveSolutionCallback(interval = 10,
                                      solution_variables = cons2prim)
 
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
-# stepsize_callback = StepsizeCallback(cfl = 1.4)
+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)
+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_semi_custom, CarpenterKennedy2N54(williamson_condition = false),
-            dt = pi * 1e-3, adaptive = false, save_everystep = false, callback = callbacks);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
 
 # Print the timer summary
 summary_callback()

examples/elixir_spherical_advection_covariant.jl

@@ -8,7 +8,7 @@ using OrdinaryDiffEq, Trixi, TrixiAtmo
 # Problem definition
 
 function initial_condition_advection_sphere(x, t, ::CovariantLinearAdvectionEquation2D)
-    
+
     # Gaussian density
     rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2))
 

src/TrixiAtmo.jl

@@ -17,6 +17,8 @@ using StrideArrays: StrideArray, StaticInt, PtrArray
 using LinearAlgebra: norm, dot
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 
+@reexport using StaticArrays: SVector
+
 include("auxiliary/auxiliary.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")

src/equations/covariant_advection.jl

@@ -30,7 +30,8 @@ end
 """
     Compute a given contravariant flux component
 """
-@inline function Trixi.flux(u, orientation::Integer, ::CovariantLinearAdvectionEquation2D)
+@inline function Trixi.flux(u, orientation::Integer,
+                            ::CovariantLinearAdvectionEquation2D)
     z = zero(eltype(u))
     return SVector(u[orientation + 1] * u[1], z, z)
 end
@@ -62,13 +63,12 @@ end
     return abs(u[2]), abs(u[3])
 end
 
-
 """
     Convert from 4-vector of scalar plus 3 Cartesian components to 3-vector of scalar plus 2 contravariant velocity/momentum components
 """
-@inline function cartesian2contravariant(u_cartesian, ::CovariantLinearAdvectionEquation2D,
+@inline function cartesian2contravariant(u_cartesian,
+                                         ::CovariantLinearAdvectionEquation2D,
                                          i, j, element, cache)
-
     (; contravariant_vectors, inverse_jacobian) = cache.elements
 
     Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, j,
@@ -89,13 +89,11 @@ end
 """
 @inline function contravariant2cartesian(u_node, ::CovariantLinearAdvectionEquation2D,
                                          i, j, element, cache)
-
     A11, A21, A31, A12, A22, A32 = view(cache.elements.jacobian_matrix, :, :, i, j,
                                         element)
     return SVector(u_node[1],
                    A11 * u_node[2] + A12 * u_node[3],
                    A21 * u_node[2] + A22 * u_node[3],
                    A31 * u_node[2] + A32 * u_node[3])
 end
-
 end # @muladd