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.27.0
+        uses: crate-ci/typos@v1.28.1

.github/workflows/ci.yml

@@ -66,7 +66,7 @@ jobs:
 
       - name: Upload coverage data (Codecov)
         if: ${{ matrix.os == 'ubuntu-latest' }}
-        uses: codecov/codecov-action@v4
+        uses: codecov/codecov-action@v5
         with:
           files: lcov.info
           token: ${{ secrets.CODECOV_TOKEN }}

Project.toml

@@ -8,6 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
@@ -20,6 +21,7 @@ LinearAlgebra = "1"
 LoopVectorization = "0.12.118"
 MuladdMacro = "0.2.2"
 NLsolve = "4.5.1"
+Printf = "1"
 Reexport = "1.0"
 Static = "0.8, 1"
 StaticArrayInterface = "1.5.1"

examples/elixir_shallowwater_cubed_sphere_shell_advection.jl

@@ -3,12 +3,22 @@ using OrdinaryDiffEq
 using Trixi
 using TrixiAtmo
 
+# To run a convergence test, we have two options:
+# 1. Use the p4est variable initial_refinement_level to refine the grid:
+#    - To do this, line 46 ("initial_refinement_level = 0") must NOT be a comment
+#    - Call convergence_test("../examples/elixir_shallowwater_cubed_sphere_shell_advection.jl", 4, initial_refinement_level = 0)
+#    - NOT OPTIMAL: Good convergence the first iterations, but then it stagnates. Reason: The geometry does not improve with refinement.
+# 2. Use the variable trees_per_face_dimension of P4estMeshCubedSphere2D
+#    - To do this, line 46 ("initial_refinement_level = 0") MUST BE commented/removed.
+#    - Call convergence_test("../examples/elixir_shallowwater_cubed_sphere_shell_advection.jl", 4, cells_per_dimension = (3,3))
+#    - OPTIMAL convergence of polydeg + 1. Reason: The geometry improves with refinement.
+
 ###############################################################################
 # semidiscretization of the linear advection equation
 
-cells_per_dimension = 5
 initial_condition = initial_condition_gaussian
 polydeg = 3
+cells_per_dimension = (5, 5)
 
 # We use the ShallowWaterEquations3D equations structure but modify the rhs! function to
 # convert it to a variable-coefficient advection equation
@@ -37,8 +47,8 @@ function source_terms_convert_to_linear_advection(u, du, x, t,
 end
 
 # Create a 2D cubed sphere mesh the size of the Earth
-mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = 3,
-                              initial_refinement_level = 0,
+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 
@@ -59,12 +69,12 @@ ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY))
 summary_callback = SummaryCallback()
 
 # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
-analysis_callback = AnalysisCallback(semi, interval = 10,
+analysis_callback = AnalysisCallback(semi, interval = 100,
                                      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,
+save_solution = SaveSolutionCallback(interval = 100,
                                      solution_variables = cons2prim)
 
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step

examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl

@@ -0,0 +1,91 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+
+# To run a convergence test, we have two options:
+# 1. Use the p4est variable initial_refinement_level to refine the grid:
+#    - To do this, line 46 ("initial_refinement_level = 0") must NOT be a comment
+#    - Call convergence_test("../examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl", 4, initial_refinement_level = 0)
+#    - NOT OPTIMAL: Good convergence the first iterations, but then it stagnates. Reason: The geometry does not improve with refinement.
+# 2. Use the variable trees_per_face_dimension of P4estMeshQuadIcosahedron2D
+#    - To do this, line 46 ("initial_refinement_level = 0") MUST BE commented/removed
+#    - Call convergence_test("../examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl", 4, cells_per_dimension = (1,1))
+#    - OPTIMAL convergence of polydeg + 1. Reason: The geometry improves with refinement.
+
+###############################################################################
+# semidiscretization of the linear advection equation
+initial_condition = initial_condition_gaussian
+polydeg = 3
+cells_per_dimension = (2, 2)
+
+# 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
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
+
+# 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::ShallowWaterEquations3D,
+                                                  normal_direction)
+    v1 = u[2] / u[1]
+    v2 = u[3] / u[1]
+    v3 = u[4] / u[1]
+
+    s2 = du[1] * v1 - du[2]
+    s3 = du[1] * v2 - du[3]
+    s4 = du[1] * v3 - du[4]
+
+    return SVector(0.0, s2, s3, s4, 0.0)
+end
+
+# Create a 2D quad-based icosahedral mesh the size of the Earth
+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)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver,
+                                    source_terms = source_terms_convert_to_linear_advection)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+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 = 100,
+                                     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 = 100,
+                                     solution_variables = cons2prim)
+
+# 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, # 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

@@ -48,7 +48,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 = cons2cons)
+                                     solution_variables = contravariant2spherical)
 
 # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
 stepsize_callback = StepsizeCallback(cfl = 0.7)

src/TrixiAtmo.jl

@@ -11,6 +11,7 @@ module TrixiAtmo
 using Reexport: @reexport
 using Trixi
 using MuladdMacro: @muladd
+using Printf: @sprintf
 using Static: True, False
 using StrideArrays: PtrArray
 using StaticArrayInterface: static_size
@@ -36,7 +37,8 @@ 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 P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct,
+       MetricTermsInvariantCurl
 export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
        EARTH_ROTATION_RATE, SECONDS_PER_DAY
 export spherical2contravariant, contravariant2spherical, spherical2cartesian,

src/callbacks_step/analysis_covariant.jl

@@ -11,13 +11,13 @@ function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t,
     Trixi.integrate_via_indices(u, mesh, equations, dg, cache,
                                 du) do u, i, j, element, equations, dg, du_node
         # Get auxiliary variables, solution variables, and time derivative at given node
-        a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
+        aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
         u_node = Trixi.get_node_vars(u, equations, dg, i, j, element)
         du_node = Trixi.get_node_vars(du, equations, dg, i, j, element)
 
         # compute ∂S/∂u ⋅ ∂u/∂t, where the entropy variables ∂S/∂u depend on the solution 
         # and auxiliary variables
-        dot(cons2entropy(u_node, a_node, equations), du_node)
+        dot(cons2entropy(u_node, aux_node, equations), du_node)
     end
 end
 
@@ -44,15 +44,15 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
 
             # Convert exact solution into contravariant components using geometric
             # information stored in aux vars
-            a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
-            u_exact = initial_condition(x_node, t, a_node, equations)
+            aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
+            u_exact = initial_condition(x_node, t, aux_node, equations)
 
             # Compute the difference as usual
             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 
-            sqrtG = area_element(a_node, equations)
+            sqrtG = area_element(aux_node, equations)
             l2_error += diff .^ 2 * (weights[i] * weights[j] * sqrtG)
 
             # Compute Linf error as usual

src/callbacks_step/callbacks_step.jl

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

src/callbacks_step/save_solution_covariant.jl

@@ -0,0 +1,53 @@
+# 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)
+    (; 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,
+                                                         1), equations))
+    # Allocate storage for output
+    data = Array{eltype(u)}(undef, n_vars, nnodes(dg), nnodes(dg), nelements(dg, cache))
+
+    # Loop over all nodes and convert variables, passing in auxiliary variables
+    Trixi.@threaded for element in eachelement(dg, cache)
+        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)
+            u_node_converted = solution_variables(u_node, aux_node, equations)
+            for v in 1:n_vars
+                data[v, i, j, element] = u_node_converted[v]
+            end
+        end
+    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

src/callbacks_step/stepsize_dg2d.jl

@@ -59,8 +59,8 @@ function Trixi.max_dt(u, t, mesh::P4estMesh{2}, constant_speed::False,
         max_lambda1 = max_lambda2 = zero(max_scaled_speed)
         for j in eachnode(dg), i in eachnode(dg)
             u_node = Trixi.get_node_vars(u, equations, dg, i, j, element)
-            a_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
-            lambda1, lambda2 = Trixi.max_abs_speeds(u_node, a_node, equations)
+            aux_node = get_node_aux_vars(aux_node_vars, equations, dg, i, j, element)
+            lambda1, lambda2 = Trixi.max_abs_speeds(u_node, aux_node, equations)
             max_lambda1 = max(max_lambda1, lambda1)
             max_lambda2 = max(max_lambda2, lambda2)
         end

src/equations/covariant_advection.jl

@@ -47,6 +47,26 @@ function velocity(u, ::CovariantLinearAdvectionEquation2D)
     return SVector(u[2], u[3])
 end
 
+# Convert contravariant velocity/momentum components to zonal and meridional components
+@inline function contravariant2spherical(u, aux_vars,
+                                         equations::CovariantLinearAdvectionEquation2D)
+    vlon, vlat = basis_covariant(aux_vars, equations) * velocity(u, equations)
+    return SVector(u[1], vlon, vlat)
+end
+
+# Convert zonal and meridional velocity/momentum components to contravariant components
+@inline function spherical2contravariant(u_spherical, aux_vars,
+                                         equations::CovariantLinearAdvectionEquation2D)
+    vcon1, vcon2 = basis_covariant(aux_vars, equations) \
+                   velocity(u_spherical, equations)
+    return SVector(u_spherical[1], vcon1, vcon2)
+end
+
+function Trixi.varnames(::typeof(contravariant2spherical),
+                        ::CovariantLinearAdvectionEquation2D)
+    return ("scalar", "vlon", "vlat")
+end
+
 # We will define the "entropy variables" here to just be the scalar variable in the first 
 # slot, with zeros in the second and third positions
 @inline function Trixi.cons2entropy(u, aux_vars,
@@ -91,19 +111,4 @@ end
                                       equations::CovariantLinearAdvectionEquation2D)
     return abs.(velocity(u, equations))
 end
-
-# Convert contravariant velocity/momentum components to zonal and meridional components
-@inline function contravariant2spherical(u, aux_vars,
-                                         equations::CovariantLinearAdvectionEquation2D)
-    vlon, vlat = basis_covariant(aux_vars, equations) * velocity(u, equations)
-    return SVector(u[1], vlon, vlat)
-end
-
-# Convert zonal and meridional velocity/momentum components to contravariant components
-@inline function spherical2contravariant(u_spherical, aux_vars,
-                                         equations::CovariantLinearAdvectionEquation2D)
-    vcon1, vcon2 = basis_covariant(aux_vars, equations) \
-                   velocity(u_spherical, equations)
-    return SVector(u_spherical[1], vcon1, vcon2)
-end
 end # @muladd

src/equations/equations.jl

@@ -48,6 +48,13 @@ dispatching on the return type.
 """
 @inline have_aux_node_vars(::AbstractEquations) = False()
 
+# cons2cons method which takes in unused aux_vars variable
+@inline Trixi.cons2cons(u, aux_vars, equations) = u
+
+# If no auxiliary variables are passed into the conversion to spherical coordinates, do not 
+# do any conversion.
+@inline contravariant2spherical(u, equations) = u
+
 # For the covariant form, the auxiliary variables are the the NDIMS^2 entries of the 
 # covariant basis matrix
 @inline have_aux_node_vars(::AbstractCovariantEquations) = True()

src/equations/reference_data.jl

@@ -52,7 +52,7 @@ for Cartesian and covariant formulations.
 @inline function initial_condition_gaussian(x, t)
     RealT = eltype(x)
 
-    a = EARTH_RADIUS  # radius of the sphere in metres
+    a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) #radius of the sphere
     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
@@ -68,7 +68,14 @@ for Cartesian and covariant formulations.
     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)
+    r_xy = x[1]^2 + x[2]^2
+    if r_xy == 0.0
+        lon = pi / 2
+    else
+        lon = atan(x[2], x[1])
+    end
+
+    lat = asin(x[3] / a)
     vlon = V * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha))
     vlat = -V * sin(lon) * sin(alpha)
 

src/meshes/meshes.jl

@@ -1,2 +1,2 @@
-export P4estMeshCubedSphere2D
 include("p4est_cubed_sphere.jl")
+include("p4est_icosahedron_quads.jl")