Skip to content

Commit

Permalink
Add alternative element-local mapping (#34)
Browse files Browse the repository at this point in the history 
Co-authored-by: Andrés Rueda-Ramírez <[email protected]>
  • Loading branch information
@tristanmontoya @amrueda
tristanmontoya and amrueda authored Nov 6, 2024
1 parent ba74a0b commit 3b6f7c1
Show file tree
Hide file tree
Showing 5 changed files with 138 additions and 22 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,6 @@ run
run/*
out
out/*
docs/build
docs/build/*
.DS_Store
5 changes: 4 additions & 1 deletion examples/elixir_euler_spherical_advection_cartesian.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
111 changes: 97 additions & 14 deletions src/meshes/p4est_cubed_sphere.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
11 changes: 6 additions & 5 deletions src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
31 changes: 29 additions & 2 deletions test/test_spherical_advection.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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=[
Expand All @@ -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
Expand Down

0 comments on commit 3b6f7c1

Please sign in to comment.