diff --git a/examples/elixir_spherical_advection_cartesian.jl b/examples/elixir_spherical_advection_cartesian.jl index 9a0e7b5..b887955 100644 --- a/examples/elixir_spherical_advection_cartesian.jl +++ b/examples/elixir_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 88d4352..855c1ba 100644 --- a/src/meshes/p4est_cubed_sphere.jl +++ b/src/meshes/p4est_cubed_sphere.jl @@ -29,8 +29,10 @@ """ 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 @@ -41,8 +43,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_cartesian!(tree_node_coordinates, nodes, + trees_per_face_dimension, + radius) + end p4est = Trixi.new_p4est(connectivity, initial_refinement_level) @@ -274,10 +282,11 @@ function connectivity_cubed_sphere_2D(trees_per_face_dimension) 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) +function calc_tree_node_coordinates_cubed_sphere_cartesian!(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)) @@ -310,4 +319,61 @@ 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. +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 planar 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 + +@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/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 025891a..c11ebdf 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -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