Merge branch 'main' into arr/shallow_water_cartesian_topography

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.28.1
+        uses: crate-ci/typos@v1.29.0

Project.toml

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

examples/elixir_shallowwater_cubed_sphere_shell_advection.jl

@@ -51,9 +51,7 @@ mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = 3,
                               #initial_refinement_level = 0, # Comment to use cells_per_dimension in the convergence test
                               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 = transform_to_cartesian(initial_condition, equations)
+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,

examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl

@@ -46,9 +46,7 @@ mesh = P4estMeshQuadIcosahedron2D(cells_per_dimension[1], EARTH_RADIUS,
                                   #initial_refinement_level = 0,
                                   polydeg = polydeg)
 
-# Convert initial condition given in terms of zonal and meridional velocity components to 
-# one given in terms of Cartesian momentum components
-initial_condition_transformed = transform_to_cartesian(initial_condition, equations)
+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,

examples/elixir_spherical_advection_covariant.jl → examples/elixir_spherical_advection_covariant_cubed_sphere.jl

@@ -1,32 +1,31 @@
 ###############################################################################
 # DGSEM for the linear advection equation on the cubed sphere
 ###############################################################################
+# To run a convergence test, use
+# convergence_test("../examples/elixir_spherical_advection_covariant_cubed_sphere.jl", 4, cells_per_dimension = (3,3))
 
 using OrdinaryDiffEq, Trixi, TrixiAtmo
 
 ###############################################################################
 # Spatial discretization
 
-cells_per_dimension = 5
+cells_per_dimension = (5, 5)
 initial_condition = initial_condition_gaussian
 
-equations = CovariantLinearAdvectionEquation2D()
+equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = GlobalCartesianCoordinates())
 
 # 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. For the covariant form to work 
 # properly, we currently need polydeg to equal that of the solver, 
-# initial_refinement_level = 0, and element_local_mapping = true.
-mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS,
+# initial_refinement_level = 0 (default), and element_local_mapping = true.
+mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS,
                               polydeg = Trixi.polydeg(solver),
-                              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 = transform_to_contravariant(initial_condition, equations)
+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)
@@ -48,7 +47,7 @@ analysis_callback = AnalysisCallback(semi, interval = 10,
 
 # The SaveSolutionCallback allows to save the solution to a file in regular intervals
 save_solution = SaveSolutionCallback(interval = 10,
-                                     solution_variables = contravariant2spherical)
+                                     solution_variables = contravariant2global)
 
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
 stepsize_callback = StepsizeCallback(cfl = 0.7)

examples/elixir_spherical_advection_covariant_quad_icosahedron.jl

@@ -0,0 +1,66 @@
+###############################################################################
+# DGSEM for the linear advection equation on the cubed sphere
+###############################################################################
+# To run a convergence test, use
+# convergence_test("../examples/elixir_spherical_advection_covariant_quad_icosahedron.jl", 4, cells_per_dimension = (1,1))
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Spatial discretization
+
+cells_per_dimension = (2, 2)
+initial_condition = initial_condition_gaussian
+
+equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = GlobalCartesianCoordinates())
+
+# 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. For the covariant form to work 
+# properly, we currently need polydeg to equal that of the solver, and
+# initial_refinement_level = 0 (default)
+mesh = P4estMeshQuadIcosahedron2D(cells_per_dimension[1], EARTH_RADIUS,
+                                  polydeg = Trixi.polydeg(solver))
+
+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)
+
+###############################################################################
+# 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 = contravariant2global)
+
+# 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

@@ -15,7 +15,7 @@ using Printf: @sprintf
 using Static: True, False
 using StrideArrays: PtrArray
 using StaticArrayInterface: static_size
-using LinearAlgebra: norm, dot
+using LinearAlgebra: cross, norm, dot, det
 using Reexport: @reexport
 using LoopVectorization: @turbo
 using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
@@ -32,6 +32,7 @@ include("callbacks_step/callbacks_step.jl")
 
 export CompressibleMoistEulerEquations2D, ShallowWaterEquations3D,
        CovariantLinearAdvectionEquation2D
+export GlobalCartesianCoordinates, GlobalSphericalCoordinates
 
 export flux_chandrashekar, flux_LMARS
 
@@ -43,9 +44,9 @@ export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProdu
        MetricTermsInvariantCurl
 export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
        EARTH_ROTATION_RATE, SECONDS_PER_DAY
-export spherical2contravariant, contravariant2spherical, spherical2cartesian,
-       transform_to_cartesian, transform_to_contravariant
-export initial_condition_gaussian
+export global2contravariant, contravariant2global, spherical2cartesian,
+       transform_initial_condition
+export initial_condition_gaussian, initial_condition_gaussian_cartesian
 
 export examples_dir
 end # module TrixiAtmo

src/callbacks_step/analysis_covariant.jl

@@ -1,7 +1,44 @@
 @muladd begin
 #! format: noindent
 
-# Entropy time derivative which uses auxiliary variables
+# For the covariant form, we want to integrate using the exact area element 
+# √G = (det(AᵀA))^(1/2), which is stored in cache.auxiliary_variables, not the approximate 
+# area element used in the Cartesian formulation, which stored in cache.elements.
+function Trixi.integrate_via_indices(func::Func, u,
+                                     mesh::Union{StructuredMesh{2},
+                                                 StructuredMeshView{2},
+                                                 UnstructuredMesh2D, P4estMesh{2},
+                                                 T8codeMesh{2}},
+                                     equations::AbstractCovariantEquations{2},
+                                     dg::DGSEM, cache, args...;
+                                     normalize = true) where {Func}
+    (; weights) = dg.basis
+    (; aux_node_vars) = cache.auxiliary_variables
+
+    # Initialize integral with zeros of the right shape
+    integral = zero(func(u, 1, 1, 1, equations, dg, args...))
+    total_volume = zero(real(mesh))
+
+    # Use quadrature to numerically integrate over entire domain
+    for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
+            sqrtG = area_element(aux_node, equations)
+            integral += weights[i] * weights[j] * sqrtG *
+                        func(u, i, j, element, equations, dg, args...)
+            total_volume += weights[i] * weights[j] * sqrtG
+        end
+    end
+
+    # Normalize with total volume
+    if normalize
+        integral = integral / total_volume
+    end
+
+    return integral
+end
+
+# Entropy time derivative for cons2entropy function which depends on auxiliary variables
 function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t,
                        mesh::P4estMesh{2},
                        equations::AbstractCovariantEquations{2}, dg::DG, cache)
@@ -28,7 +65,6 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
     (; weights) = dg.basis
     (; node_coordinates) = cache.elements
     (; aux_node_vars) = cache.auxiliary_variables
-
     # Set up data structures
     l2_error = zero(func(Trixi.get_node_vars(u, equations, dg, 1, 1, 1), equations))
     linf_error = copy(l2_error)
@@ -51,7 +87,7 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
             u_numerical = Trixi.get_node_vars(u, equations, dg, i, j, element)
             diff = func(u_exact, equations) - func(u_numerical, equations)
 
-            # For the L2 error, integrate with respect to volume element stored in aux vars 
+            # For the L2 error, integrate with respect to area element stored in aux vars 
             sqrtG = area_element(aux_node, equations)
             l2_error += diff .^ 2 * (weights[i] * weights[j] * sqrtG)
 

src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl

@@ -42,7 +42,9 @@ end
 
 function Trixi.save_solution_file(u, time, dt, timestep,
                                   mesh::P4estMesh{2},
-                                  equations::AbstractEquations{3}, dg::DG, cache,
+                                  equations::Union{AbstractEquations{3},
+                                                   AbstractCovariantEquations{2}},
+                                  dg::DG, cache,
                                   solution_callback,
                                   element_variables = Dict{Symbol, Any}(),
                                   node_variables = Dict{Symbol, Any}();

src/callbacks_step/save_solution_covariant.jl

@@ -1,13 +1,17 @@
+@muladd begin
+#! format: noindent
+
 # Convert to another set of variables using a solution_variables function that depends on 
 # auxiliary variables
 function convert_variables(u, solution_variables, mesh::P4estMesh{2},
-                           equations, dg, cache)
+                           equations::AbstractCovariantEquations{2}, dg, cache)
     (; aux_node_vars) = cache.auxiliary_variables
 
     # Extract the number of solution variables to be output 
     # (may be different than the number of conservative variables) 
     n_vars = length(solution_variables(Trixi.get_node_vars(u, equations, dg, 1, 1, 1),
-                                       get_node_aux_vars(aux_node_vars, equations, dg, 1, 1,
+                                       get_node_aux_vars(aux_node_vars, equations, dg,
+                                                         1, 1,
                                                          1), equations))
     # Allocate storage for output
     data = Array{eltype(u)}(undef, n_vars, nnodes(dg), nnodes(dg), nelements(dg, cache))
@@ -25,29 +29,4 @@ function convert_variables(u, solution_variables, mesh::P4estMesh{2},
     end
     return data
 end
-
-# Specialized save_solution_file method that supports a solution_variables function which 
-# depends on auxiliary variables. The conversion must be defined as solution_variables(u, 
-# aux_vars, equations), and an additional method must be defined as solution_variables(u,
-# equations) = u, such that no conversion is done when auxiliary variables are not provided.
-function Trixi.save_solution_file(u_ode, t, dt, iter,
-                                  semi::SemidiscretizationHyperbolic{<:Trixi.AbstractMesh,
-                                                                     <:AbstractCovariantEquations},
-                                  solution_callback,
-                                  element_variables = Dict{Symbol, Any}(),
-                                  node_variables = Dict{Symbol, Any}();
-                                  system = "")
-    mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(semi)
-    u = Trixi.wrap_array_native(u_ode, semi)
-
-    # Perform the variable conversion at each node
-    data = convert_variables(u, solution_callback.solution_variables, mesh, equations,
-                             solver, cache)
-
-    # Call the existing Trixi.save_solution_file, which will use solution_variables(u, 
-    # equations). Since the variables are already converted above, we therefore require 
-    # solution_variables(u, equations) = u.
-    Trixi.save_solution_file(data, t, dt, iter, mesh, equations, solver, cache,
-                             solution_callback, element_variables,
-                             node_variables, system = system)
-end
+end # muladd