Skip to content

Commit

Permalink
Merge branch 'tm/element_local_mapping' into tm/covariant_form
Browse files Browse the repository at this point in the history
  • Loading branch information
@tristanmontoya
tristanmontoya committed Sep 10, 2024
2 parents 2d5e6c5 + bcb9fda commit f4f94fb
Show file tree
Hide file tree
Showing 3 changed files with 106 additions and 10 deletions.
5 changes: 4 additions & 1 deletion examples/elixir_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
82 changes: 74 additions & 8 deletions src/meshes/p4est_cubed_sphere.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
29 changes: 28 additions & 1 deletion test/test_spherical_advection.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 f4f94fb

Please sign in to comment.