streamline tests and improve readability

trixi-framework · Sep 10, 2024 · 4324e62 · 4324e62
1 parent d2cd693
commit 4324e62
Show file tree

Hide file tree

Showing 5 changed files with 48 additions and 129 deletions.
diff --git a/examples/elixir_spherical_advection_covariant.jl b/examples/elixir_spherical_advection_covariant.jl
@@ -38,9 +38,11 @@ end
 polydeg = 3
 cells_per_dimension = 2
 
+element_local_mapping = true
+
 mesh = P4estMeshCubedSphere2D(5, EARTH_RADIUS, polydeg = polydeg,
                               initial_refinement_level = 0,
-                              element_local_mapping = true)
+                              element_local_mapping = element_local_mapping)
 
 equations = CovariantLinearAdvectionEquation2D()
 
@@ -50,7 +52,7 @@ initial_condition = initial_condition_advection_earth
 volume_integral = VolumeIntegralWeakForm()
 
 # Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux
-solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs, 
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs,
                volume_integral = volume_integral)
 
 # A semidiscretization collects data structures and functions for the spatial discretization

diff --git a/examples/elixir_spherical_advection_covariant_flux_differencing.jl b/examples/elixir_spherical_advection_covariant_flux_differencing.jl
diff --git a/src/meshes/p4est_cubed_sphere.jl b/src/meshes/p4est_cubed_sphere.jl
@@ -2,31 +2,32 @@
 #! format: noindent
 
 """
-        P4estMeshCubedSphere2D(trees_per_face_dimension, radius;
-                                polydeg, RealT=Float64,
-                                initial_refinement_level=0, unsaved_changes=true,
-                                p4est_partition_allow_for_coarsening=true)
-
-    Build a "Cubed Sphere" mesh as a 2D `P4estMesh` with
-    `6 * trees_per_face_dimension^2` trees.
-
-    The mesh will have no boundaries.
-
-    # Arguments
-
-      - `trees_per_face_dimension::Integer`: the number of trees in the two local dimensions of
-        each face.
-      - `radius::Integer`: the radius of the sphere.
-      - `polydeg::Integer`: polynomial degree used to store the geometry of the mesh.
-        The mapping will be approximated by an interpolation polynomial
-        of the specified degree for each tree.
-      - `RealT::Type`: the type that should be used for coordinates.
-      - `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts.
-      - `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file.
-      - `p4est_partition_allow_for_coarsening::Bool`: Must be `true` when using AMR to make mesh adaptivity
-        independent of domain partitioning. Should be `false` for static meshes
-        to permit more fine-grained partitioning.
-    """
+    P4estMeshCubedSphere2D(trees_per_face_dimension, radius;
+                            polydeg, RealT=Float64,
+                            initial_refinement_level=0, unsaved_changes=true,
+                            p4est_partition_allow_for_coarsening=true,
+                            element_local_mapping=false)
+
+Build a "Cubed Sphere" mesh as a 2D `P4estMesh` with
+`6 * trees_per_face_dimension^2` trees.
+
+The mesh will have no boundaries.
+
+# Arguments
+- `trees_per_face_dimension::Integer`: the number of trees in the two local dimensions of
+                                       each face.
+- `radius::Integer`: the radius of the sphere.
+- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh.
+                      The mapping will be approximated by an interpolation polynomial
+                      of the specified degree for each tree.
+- `RealT::Type`: the type that should be used for coordinates.
+- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts.
+- `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file.
+- `p4est_partition_allow_for_coarsening::Bool`: Must be `true` when using AMR to make mesh adaptivity
+                                                independent of domain partitioning. Should be `false` for static meshes
+                                                to permit more fine-grained partitioning.
+- `element_local_mapping::Bool`: option to use the element-local mapping from Guba et al. (see https://doi.org/10.5194/gmd-7-2803-2014, Appendix A).
+"""
 function P4estMeshCubedSphere2D(trees_per_face_dimension, radius;
                                 polydeg, RealT = Float64,
                                 initial_refinement_level = 0,

diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d_cartesian.jl
@@ -125,13 +125,14 @@ function init_elements_2d_manifold_in_3d!(elements,
                                           mesh::Union{P4estMesh{2}, T8codeMesh{2}},
                                           equations::AbstractEquations{3},
                                           basis::LobattoLegendreBasis)
-    (; node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian) = elements
+    (; node_coordinates, jacobian_matrix,
+    contravariant_vectors, inverse_jacobian) = elements
 
     calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis)
 
     for element in 1:Trixi.ncells(mesh)
 
-        # Compute Jacobian matrix as usual using Giraldo's cross-product form
+        # Compute Jacobian matrix as usual
         Trixi.calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis)
 
         # Compute contravariant vectors
@@ -231,9 +232,9 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                               jacobian_matrix, node_coordinates,
                                               basis::LobattoLegendreBasis)
     for j in eachnode(basis), i in eachnode(basis)
-        r = sqrt(node_coordinates[1, i, j, element]^2 +
-                 node_coordinates[2, i, j, element]^2 +
-                 node_coordinates[3, i, j, element]^2)
+        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
@@ -248,7 +249,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[l, 2, i, j,
                                                                           element] *
                                                           node_coordinates[m, i, j,
-                                                                           element]) / r
+                                                                           element]) /
+                                                         radius
 
             contravariant_vectors[n, 2, i, j, element] = (jacobian_matrix[l, 1, i, j,
                                                                           element] *
@@ -258,7 +260,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[m, 1, i, j,
                                                                           element] *
                                                           node_coordinates[l, i, j,
-                                                                           element]) / r
+                                                                           element]) /
+                                                         radius
 
             contravariant_vectors[n, 3, i, j, element] = (jacobian_matrix[m, 1, i, j,
                                                                           element] *
@@ -268,7 +271,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr
                                                           jacobian_matrix[m, 2, i, j,
                                                                           element] *
                                                           jacobian_matrix[l, 1, i, j,
-                                                                          element]) / r
+                                                                          element]) /
+                                                         radius
         end
     end
 

diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl
@@ -23,7 +23,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples")
                             1.37453400e-02,
                             0.00000000e+00,
                             4.09322999e-01,
-                        ], element_local_mapping=false)
+                        ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -34,7 +34,6 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples")
     end
 end
 
-
 @trixiatmo_testset "Spherical advection, covariant weak form" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_spherical_advection_covariant.jl"),
@@ -58,9 +57,11 @@ end
     end
 end
 
-@trixiatmo_testset "Spherical advection, covariant flux-differencing form" begin
+# The covariant flux-differencing form should be equivalent to the weak form when the 
+# arithmetic mean is used as the two-point flux
+@trixiatmo_testset "Spherical advection, covariant flux-differencing, arithmetic mean" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
-                                 "elixir_spherical_advection_covariant_flux_differencing.jl"),
+                                 "elixir_spherical_advection_covariant.jl"),
                         l2=[
                             1.00016205e+00,
                             0.00000000e+00,
@@ -70,7 +71,7 @@ end
                             1.42451839e+01,
                             0.00000000e+00,
                             0.00000000e+00,
-                        ])
+                        ], volume_integral=VolumeIntegralFluxDifferencing(flux_central))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let