diff --git a/.gitignore b/.gitignore index 8a53047..c5ec89b 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,6 @@ run run/* out out/* +docs/build +docs/build/* .DS_Store diff --git a/examples/elixir_euler_spherical_advection_cartesian.jl b/examples/elixir_euler_spherical_advection_cartesian.jl index 9a0e7b5..b887955 100644 --- a/examples/elixir_euler_spherical_advection_cartesian.jl +++ b/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) diff --git a/src/meshes/p4est_cubed_sphere.jl b/src/meshes/p4est_cubed_sphere.jl index 07f0d9b..ca4cd1e 100644 --- a/src/meshes/p4est_cubed_sphere.jl +++ b/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 diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl index 248523c..8aa6583 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl @@ -32,11 +32,12 @@ end 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 diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 31f7a29..c9e9300 100644 --- a/test/test_spherical_advection.jl +++ b/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