Add alternative element-local mapping

.gitignore

@@ -3,4 +3,6 @@ run
 run/*
 out
 out/*
+docs/build
+docs/build/*
 .DS_Store

examples/elixir_euler_spherical_advection_cartesian.jl

@@ -101,8 +101,11 @@ end
 
 initial_condition = initial_condition_advection_sphere
 
+element_local_mapping = false
+
 mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg,
-                                        initial_refinement_level = 0)
+                                        initial_refinement_level = 0,
+                                        element_local_mapping = element_local_mapping)
 
 # A semidiscretization collects data structures and functions for the spatial discretization
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

src/meshes/p4est_cubed_sphere.jl

@@ -5,7 +5,8 @@
     P4estMeshCubedSphere2D(trees_per_face_dimension, radius;
                             polydeg, RealT=Float64,
                             initial_refinement_level=0, unsaved_changes=true,
-                            p4est_partition_allow_for_coarsening=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.
@@ -20,16 +21,28 @@ The mesh will have no boundaries.
                       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.
+- `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.
+- `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 alternative element-local mapping from
+  Appendix A of [Guba et al. (2014)](https://doi.org/10.5194/gmd-7-2803-2014). If set to 
+  `true`, the four corner vertex positions for each element will be obtained through an 
+  equiangular gnomonic projection [(Ronchi et al. 1996)](https://doi.org/10.1006/jcph.1996.
+  0047), and the tree node coordinates within the element will be obtained by first using 
+  a bilinear mapping based on the four corner vertices, and then projecting the bilinearly 
+  mapped nodes onto the spherical surface by normalizing the resulting Cartesian 
+  coordinates and scaling by  `radius`. If set to `false`, the equiangular gnomonic 
+  projection will be used for all node coordinates.
 """
 function P4estMeshCubedSphere2D(trees_per_face_dimension, radius;
                                 polydeg, RealT = Float64,
-                                initial_refinement_level = 0, unsaved_changes = true,
-                                p4est_partition_allow_for_coarsening = true)
+                                initial_refinement_level = 0,
+                                unsaved_changes = true,
+                                p4est_partition_allow_for_coarsening = true,
+                                element_local_mapping = false)
     connectivity = connectivity_cubed_sphere_2D(trees_per_face_dimension)
 
     n_trees = 6 * trees_per_face_dimension^2
@@ -40,8 +53,14 @@ function P4estMeshCubedSphere2D(trees_per_face_dimension, radius;
     tree_node_coordinates = Array{RealT, 4}(undef, 3,
                                             ntuple(_ -> length(nodes), 2)...,
                                             n_trees)
-    calc_tree_node_coordinates_cubed_sphere_2D!(tree_node_coordinates, nodes,
-                                                trees_per_face_dimension, radius)
+    if element_local_mapping
+        calc_tree_node_coordinates_cubed_sphere_local!(tree_node_coordinates, nodes,
+                                                       trees_per_face_dimension, radius)
+    else
+        calc_tree_node_coordinates_cubed_sphere_standard!(tree_node_coordinates, nodes,
+                                                          trees_per_face_dimension,
+                                                          radius)
+    end
 
     p4est = Trixi.new_p4est(connectivity, initial_refinement_level)
 
@@ -272,11 +291,14 @@ function connectivity_cubed_sphere_2D(trees_per_face_dimension)
     return connectivity
 end
 
-# Calculate physical coordinates of each node of a 2D cubed sphere mesh.
-function calc_tree_node_coordinates_cubed_sphere_2D!(node_coordinates::AbstractArray{<:Any,
-                                                                                     4},
-                                                     nodes, trees_per_face_dimension,
-                                                     radius)
+# Calculate physical coordinates of each node of a 2D cubed sphere mesh by mapping 
+# the LGL quadrature nodes onto a Cartesian element grid, then using the equiangular 
+# gnomonic mapping to place the nodes on the spherical surface
+function calc_tree_node_coordinates_cubed_sphere_standard!(node_coordinates::AbstractArray{<:Any,
+                                                                                           4},
+                                                           nodes,
+                                                           trees_per_face_dimension,
+                                                           radius)
     n_cells_x = n_cells_y = trees_per_face_dimension
 
     linear_indices = LinearIndices((n_cells_x, n_cells_y, 6))
@@ -309,4 +331,65 @@ function calc_tree_node_coordinates_cubed_sphere_2D!(node_coordinates::AbstractA
         end
     end
 end
+
+# Calculate physical coordinates of each node of a 2D cubed sphere mesh using the
+# element-local mapping from Guba et al. (see https://doi.org/10.5194/gmd-7-2803-2014,
+# Appendix A). We thank Oswald Knoth for bringing this paper to our attention.
+function calc_tree_node_coordinates_cubed_sphere_local!(node_coordinates::AbstractArray{<:Any,
+                                                                                        4},
+                                                        nodes, trees_per_face_dimension,
+                                                        radius)
+    n_cells_x = n_cells_y = trees_per_face_dimension
+
+    linear_indices = LinearIndices((n_cells_x, n_cells_y, 6))
+
+    # Get cell length in reference mesh
+    dx = 2 / n_cells_x
+    dy = 2 / n_cells_y
+
+    for direction in 1:6
+        for cell_y in 1:n_cells_y, cell_x in 1:n_cells_x
+            tree = linear_indices[cell_x, cell_y, direction]
+
+            x_offset = -1 + (cell_x - 1) * dx + dx / 2
+            y_offset = -1 + (cell_y - 1) * dy + dy / 2
+            z_offset = 0
+
+            # Vertices for bilinear mapping
+            v1 = Trixi.cubed_sphere_mapping(x_offset - dx / 2,
+                                            y_offset - dx / 2,
+                                            z_offset, radius, 0, direction)
+
+            v2 = Trixi.cubed_sphere_mapping(x_offset + dx / 2,
+                                            y_offset - dx / 2,
+                                            z_offset, radius, 0, direction)
+
+            v3 = Trixi.cubed_sphere_mapping(x_offset + dx / 2,
+                                            y_offset + dx / 2,
+                                            z_offset, radius, 0, direction)
+
+            v4 = Trixi.cubed_sphere_mapping(x_offset - dx / 2,
+                                            y_offset + dx / 2,
+                                            z_offset, radius, 0, direction)
+
+            for j in eachindex(nodes), i in eachindex(nodes)
+                # node_coordinates are the mapped reference node coordinates
+                node_coordinates[:, i, j, tree] .= local_mapping(nodes[i], nodes[j],
+                                                                 v1, v2, v3, v4, radius)
+            end
+        end
+    end
+end
+
+# Local mapping from the reference quadrilateral to a spherical patch based on four vertex
+# positions on the sphere, provided in Cartesian coordinates
+@inline function local_mapping(xi1, xi2, v1, v2, v3, v4, radius)
+
+    # Construct a bilinear mapping based on the four corner vertices
+    x_bilinear = 0.25f0 * ((1 - xi1) * (1 - xi2) * v1 + (1 + xi1) * (1 - xi2) * v2 +
+                  (1 + xi1) * (1 + xi2) * v3 + (1 - xi1) * (1 + xi2) * v4)
+
+    # Project the mapped local coordinates onto the sphere
+    return radius * x_bilinear / norm(x_bilinear)
+end
 end # @muladd

src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl

@@ -32,11 +32,12 @@
     size(elements.node_coordinates, ndims(elements) + 2)
 end
 @inline Base.ndims(::P4estElementContainerPtrArray{NDIMS}) where {NDIMS} = NDIMS
-@inline function Base.eltype(::P4estElementContainerPtrArray{NDIMS, RealT, uEltype}) where {
-                                                                                            NDIMS,
-                                                                                            RealT,
-                                                                                            uEltype
-                                                                                            }
+@inline function Base.eltype(::P4estElementContainerPtrArray{NDIMS,
+                                                             RealT, uEltype}) where {
+                                                                                     NDIMS,
+                                                                                     RealT,
+                                                                                     uEltype
+                                                                                     }
     uEltype
 end
 

test/test_spherical_advection.jl

@@ -7,7 +7,7 @@ include("test_trixiatmo.jl")
 
 EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples")
 
-@trixiatmo_testset "elixir_euler_spherical_advection_cartesian" begin
+@trixiatmo_testset "Spherical advection, Cartesian form with standard mapping" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_spherical_advection_cartesian.jl"),
                         l2=[
@@ -23,7 +23,34 @@ 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
+        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
+
+@trixiatmo_testset "Spherical advection, Cartesian form with element-local mapping" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_spherical_advection_cartesian.jl"),
+                        l2=[
+                            8.97501204e-03,
+                            8.74316738e-03,
+                            1.99380928e-03,
+                            0.00000000e+00,
+                            6.45266612e-02,
+                        ],
+                        linf=[
+                            1.01363241e-01,
+                            1.03082705e-01,
+                            1.41424723e-02,
+                            0.00000000e+00,
+                            4.10471741e-01,
+                        ], element_local_mapping=true)
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let