Add test for spherical advection in Cartesian coordinates (#30)

* add tests for spherical advection in Cartesian coords * revert change to moist bubble case * revert change to moist bubble case * hard-code dimension of node coordinates to avoid type instability * metric terms for shell of radius != 1
trixi-framework · Sep 3, 2024 · 15062dc · 15062dc
1 parent d8b0d31
commit 15062dc
Show file tree

Hide file tree

Showing 4 changed files with 56 additions and 6 deletions.
diff --git a/examples/elixir_euler_spherical_advection_cartesian.jl b/examples/elixir_euler_spherical_advection_cartesian.jl
@@ -6,6 +6,8 @@ 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

diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl
@@ -150,7 +150,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
     # Hanging nodes will cause holes in the mesh if its polydeg is higher
     # than the polydeg of the solver.
     @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh"
-
+    @assert size(mesh.tree_node_coordinates, 1)==3 "Shell must use 3D node coordinates"
     calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis.nodes)
 end
 
@@ -162,7 +162,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
     # places and the additional information passed to the compiler makes them faster
     # than native `Array`s.
     tmp1 = Trixi.StrideArray(undef, real(mesh),
-                             Trixi.StaticInt(size(mesh.tree_node_coordinates, 1)),
+                             Trixi.StaticInt(3),
                              Trixi.static_length(nodes),
                              Trixi.static_length(mesh.nodes))
     matrix1 = Trixi.StrideArray(undef, real(mesh),
@@ -173,7 +173,6 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
 
     # Macros from `p4est`
     p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL
-    p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL
     p4est_quadrant_len(l) = 1 << (Trixi.P4EST_MAXLEVEL - l)
 
     trees = Trixi.unsafe_wrap_sc(Trixi.p4est_tree_t, mesh.p4est.trees)
@@ -208,6 +207,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates,
     return node_coordinates
 end
 
+# This only works for a sphere
 function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArray{<:Any,
                                                                                    5},
                                               element,
@@ -216,6 +216,10 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
     @unpack derivative_matrix = basis
 
     for j in eachnode(basis), i in eachnode(basis)
+        radius = sqrt(node_coordinates[1, i, j, element]^2 +
+                      node_coordinates[2, i, j, element]^2 +
+                      node_coordinates[3, i, j, element]^2)
+
         for n in 1:3
             # (n, m, l) cyclic
             m = (n % 3) + 1
@@ -229,7 +233,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[l, 2, i, j,
                                                                           element] *
                                                           node_coordinates[m, i, j,
-                                                                           element])
+                                                                           element]) /
+                                                         radius
 
             contravariant_vectors[n, 2, i, j, element] = (jacobian_matrix[l, 1, i, j,
                                                                           element] *
@@ -239,7 +244,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[m, 1, i, j,
                                                                           element] *
                                                           node_coordinates[l, i, j,
-                                                                           element])
+                                                                           element]) /
+                                                         radius
 
             contravariant_vectors[n, 3, i, j, element] = (jacobian_matrix[m, 1, i, j,
                                                                           element] *
@@ -249,7 +255,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[m, 2, i, j,
                                                                           element] *
                                                           jacobian_matrix[l, 1, i, j,
-                                                                          element])
+                                                                          element]) /
+                                                         radius
         end
     end
 

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -16,4 +16,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3)
     @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "moist_euler"
         include("test_2d_moist_euler.jl")
     end
+
+    @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection"
+        include("test_spherical_advection.jl")
+    end
 end
diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl
@@ -0,0 +1,37 @@
+module TestSphericalAdvection
+
+using Test
+using TrixiAtmo
+
+include("test_trixiatmo.jl")
+
+EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples")
+
+@trixiatmo_testset "elixir_euler_spherical_advection_cartesian" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_spherical_advection_cartesian.jl"),
+                        l2=[
+                            8.44505073e-03,
+                            8.23414117e-03,
+                            1.84210648e-03,
+                            0.00000000e+00,
+                            6.44302430e-02,
+                        ],
+                        linf=[
+                            9.48950488e-02,
+                            9.64811952e-02,
+                            1.37453400e-02,
+                            0.00000000e+00,
+                            4.09322999e-01,
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100
+    end
+end
+
+end # module