diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl index 0749c83..0e21dbd 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl +++ b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl @@ -3,11 +3,21 @@ using OrdinaryDiffEq using Trixi using TrixiAtmo +# To run a convergence test, we have two options: +# 1. Use the p4est variable initial_refinement_level to refine the grid: +# - To do this, line 46 ("initial_refinement_level = 0") must NOT be a comment +# - Call convergence_test("../examples/elixir_shallowwater_cubed_sphere_shell_advection.jl", 4, initial_refinement_level = 0) +# - NOT OPTIMAL: Good convergence the first iterations, but then it stagnates. Reason: The geometry does not improve with refinement. +# 2. Use the variable trees_per_face_dimension of P4estMeshCubedSphere2D +# - To do this, line 46 ("initial_refinement_level = 0") MUST BE commented/removed. +# - Call convergence_test("../examples/elixir_shallowwater_cubed_sphere_shell_advection.jl", 4, cells_per_dimension = (3,3)) +# - OPTIMAL convergence of polydeg + 1. Reason: The geometry improves with refinement. + ############################################################################### # semidiscretization of the linear advection equation -cells_per_dimension = 5 initial_condition = initial_condition_gaussian +cells_per_dimension = (5, 5) # We use the ShallowWaterEquations3D equations structure but modify the rhs! function to # convert it to a variable-coefficient advection equation @@ -32,8 +42,8 @@ function source_terms_convert_to_linear_advection(u, du, x, t, end # Create a 2D cubed sphere mesh the size of the Earth -mesh = P4estMeshCubedSphere2D(cells_per_dimension, EARTH_RADIUS, polydeg = 3, - initial_refinement_level = 0, +mesh = P4estMeshCubedSphere2D(cells_per_dimension[1], EARTH_RADIUS, polydeg = 3, + #initial_refinement_level = 0, # Comment to use cells_per_dimension in the convergence test element_local_mapping = false) # Convert initial condition given in terms of zonal and meridional velocity components to @@ -54,12 +64,12 @@ ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY)) summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback = AnalysisCallback(semi, interval = 10, +analysis_callback = AnalysisCallback(semi, interval = 100, save_analysis = true, extra_analysis_errors = (:conservation_error,)) # The SaveSolutionCallback allows to save the solution to a file in regular intervals -save_solution = SaveSolutionCallback(interval = 10, +save_solution = SaveSolutionCallback(interval = 100, solution_variables = cons2prim) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step diff --git a/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl b/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl new file mode 100644 index 0000000..a9ff9dc --- /dev/null +++ b/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl @@ -0,0 +1,91 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiAtmo + +# To run a convergence test, we have two options: +# 1. Use the p4est variable initial_refinement_level to refine the grid: +# - To do this, line 46 ("initial_refinement_level = 0") must NOT be a comment +# - Call convergence_test("../examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl", 4, initial_refinement_level = 0) +# - NOT OPTIMAL: Good convergence the first iterations, but then it stagnates. Reason: The geometry does not improve with refinement. +# 2. Use the variable trees_per_face_dimension of P4estMeshQuadIcosahedron2D +# - To do this, line 46 ("initial_refinement_level = 0") MUST BE commented/removed +# - Call convergence_test("../examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl", 4, cells_per_dimension = (1,1)) +# - OPTIMAL convergence of polydeg + 1. Reason: The geometry improves with refinement. + +############################################################################### +# semidiscretization of the linear advection equation +initial_condition = initial_condition_gaussian +polydeg = 3 +cells_per_dimension = (2, 2) + +# We use the ShallowWaterEquations3D equations structure but modify the rhs! function to +# convert it to a variable-coefficient advection equation +equations = ShallowWaterEquations3D(gravity_constant = 0.0) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) + +# Source term function to transform the Euler equations into a linear advection equation with variable advection velocity +function source_terms_convert_to_linear_advection(u, du, x, t, + equations::ShallowWaterEquations3D, + normal_direction) + v1 = u[2] / u[1] + v2 = u[3] / u[1] + v3 = u[4] / u[1] + + s2 = du[1] * v1 - du[2] + s3 = du[1] * v2 - du[3] + s4 = du[1] * v3 - du[4] + + return SVector(0.0, s2, s3, s4, 0.0) +end + +# Create a 2D quad-based icosahedral mesh the size of the Earth +mesh = P4estMeshQuadIcosahedron2D(cells_per_dimension[1], EARTH_RADIUS, + #initial_refinement_level = 0, + polydeg = polydeg) + +# Convert initial condition given in terms of zonal and meridional velocity components to +# one given in terms of Cartesian momentum components +initial_condition_transformed = transform_to_cartesian(initial_condition, equations) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, solver, + source_terms = source_terms_convert_to_linear_advection) + +############################################################################### +# ODE solvers, callbacks etc. + +ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY)) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 100, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.7) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 4fbe100..6803370 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -37,7 +37,8 @@ export flux_chandrashekar, flux_LMARS export velocity, waterheight, pressure, energy_total, energy_kinetic, energy_internal, lake_at_rest_error, source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier! -export P4estCubedSphere2D, MetricTermsCrossProduct, MetricTermsInvariantCurl +export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct, + MetricTermsInvariantCurl export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION, EARTH_ROTATION_RATE, SECONDS_PER_DAY export spherical2contravariant, contravariant2spherical, spherical2cartesian, diff --git a/src/equations/reference_data.jl b/src/equations/reference_data.jl index e7a59fb..a3b6883 100644 --- a/src/equations/reference_data.jl +++ b/src/equations/reference_data.jl @@ -52,7 +52,7 @@ for Cartesian and covariant formulations. @inline function initial_condition_gaussian(x, t) RealT = eltype(x) - a = EARTH_RADIUS # radius of the sphere in metres + a = sqrt(x[1]^2 + x[2]^2 + x[3]^2) #radius of the sphere V = convert(RealT, 2π) * a / (12 * SECONDS_PER_DAY) # speed of rotation alpha = convert(RealT, π / 4) # angle of rotation h_0 = 1000.0f0 # bump height in metres @@ -68,7 +68,14 @@ for Cartesian and covariant formulations. h = h_0 * exp(-b_0 * ((x[1] - x_0[1])^2 + (x[2] - x_0[2])^2 + (x[3] - x_0[3])^2)) # get zonal and meridional components of the velocity - lon, lat = atan(x[2], x[1]), asin(x[3] / a) + r_xy = x[1]^2 + x[2]^2 + if r_xy == 0.0 + lon = pi / 2 + else + lon = atan(x[2], x[1]) + end + + lat = asin(x[3] / a) vlon = V * (cos(lat) * cos(alpha) + sin(lat) * cos(lon) * sin(alpha)) vlat = -V * sin(lon) * sin(alpha) diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index 7012004..8bd9c4b 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -1,2 +1,2 @@ -export P4estMeshCubedSphere2D include("p4est_cubed_sphere.jl") +include("p4est_icosahedron_quads.jl") diff --git a/src/meshes/p4est_icosahedron_quads.jl b/src/meshes/p4est_icosahedron_quads.jl new file mode 100644 index 0000000..f15677c --- /dev/null +++ b/src/meshes/p4est_icosahedron_quads.jl @@ -0,0 +1,752 @@ +@muladd begin +#! format: noindent + +# Using element_local mapping +""" + P4estMeshQuadIcosahedron2D(trees_per_face_dimension, radius; + polydeg, RealT=Float64, + initial_refinement_level=0, unsaved_changes=true, + p4est_partition_allow_for_coarsening=true) + +Build a quad-based icosahedral mesh as a 2D `P4estMesh` with +`60 * trees_per_face_dimension^2` trees (20 triangular faces of the icosahedron, +each subdivided into 3 parent quads, each of which subdivided into `trees_per_face_dimension^2` trees). + +The node coordinates of the trees will be obtained using the element-local mapping from +Appendix A of [Guba et al. (2014)](https://doi.org/10.5194/gmd-7-2803-2014). +See [`P4estMeshCubedSphere2D`](@ref) for more information about the element-local mapping. + +The mesh will have no boundaries. + +# Arguments +- `trees_per_face_dimension::Integer`: the number of trees in the two local dimensions of + each parent quad. +- `radius::Integer`: the radius of the sphere. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. + 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. +- `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. + +!!! warning + Adaptivity and MPI parallelization are not yet supported for equations in covariant + form, and we require `initial_refinement_level = 0` for such cases. Furthermore, the + calculation of the metric terms for the covariant form currently requires `polydeg` to + be equal to the polynomial degree of the solver. +!!! +""" +function P4estMeshQuadIcosahedron2D(trees_per_face_dimension, radius; + polydeg, RealT = Float64, + initial_refinement_level = 0, + unsaved_changes = true, + p4est_partition_allow_for_coarsening = true) + connectivity = connectivity_icosahedron_2D(trees_per_face_dimension) + + n_trees = 60 * trees_per_face_dimension^2 # 20 triangles subdivided into 3 quads each + + basis = LobattoLegendreBasis(RealT, polydeg) + nodes = basis.nodes + + tree_node_coordinates = Array{RealT, 4}(undef, 3, + ntuple(_ -> length(nodes), 2)..., + n_trees) + calc_tree_node_coordinates_quad_icosahedron_local!(tree_node_coordinates, nodes, + trees_per_face_dimension, radius) + + p4est = Trixi.new_p4est(connectivity, initial_refinement_level) + + boundary_names = fill(Symbol("---"), 2 * 2, n_trees) + + return P4estMesh{2}(p4est, tree_node_coordinates, nodes, + boundary_names, "", unsaved_changes, + p4est_partition_allow_for_coarsening) +end + +# Fig 1: Illustration of the unfolded icosahedron with the numbering of the triangular faces +# ,'|. +# ,' | `. +# ,' 4 | 3 `. +# ,'_____ |_____`. +# \ /\ / +# \ 5 / \ 2 / +# \ / 1 \ / +# \/_______/_______________________________ +# \ /\ /\ /\ /\ /\ +# \ 6 / \ 7 / \ 8 / \ 9 / \ 10 / \ +# \ / 11 \ / 12 \ / 13 \ / 14 \ / 15 \ +# \/_______/_______/_______/_______/______\ +# /\ /\ +# / \ 20 / \ +# / 19 \ / 16 \ +# /______\/______\ +# . | ,' +# `.18 | 17 ,' +# `. | ,' +# `|,' +# +# Each triangle is subdivided into 3 quadrilaterals with a local (ξ,η)-coordinate system. +# +# Fig 2: +# /\ +# / \ +# / \ +# / \ +# / \ +# / \ +# / 3 \ +# / \ +# / η ξ \ +# / ⎻⎼⎼⎽⎽ ⎽⎽⎼⎼⎻ \ +# /‾⎺⎺⎻⎻⎼⎼⎽⎽ ⎽⎽⎼⎼⎻⎻⎺⎺‾ \ +# / | \ +# / | \ +# / | \ +# / η | ↑η \ +# / / 1 | | 2 \ +# / / | | \ +# / / | | \ +# / -------->ξ | └------->ξ \ +# /__________________|___________________\ +# +# Each of those quadrilaterlas is subdivided into trees_per_face_dimension^2 trees. +# +# We use the following numbering for the 12 vertices of the icosahedron +# Fig 3: +# 5 +# ,'|. +# ,' | `. +# ,' | `. +# 6,'_____1|_____`.4 +# \ /\ / +# \ / \ / +# \ / \ / +# 2\/______3/______4_______5_______6_______2 +# \ /\ /\ /\ /\ /\ +# \ / \ / \ / \ / \ / \ +# \ / \ / \ / \ / \ / \ +# \/_______/_______/_______/_______/______\7 +# 7 8 9 10 11\ /\ +# / \ / \ +# / \ / \ +# /______\/______\ +# 10 . |12 ,'8 +# `. | ,' +# `. | ,' +# `|,' +# 9 + +# Function to compute the vertices' coordinates of an icosahedron inscribed in a sphere of radius `radius` +function calc_node_coordinates_icosahedron_vertices(radius; RealT = Float64) + vertices = Array{RealT, 2}(undef, 3, 12) + + vertices[:, 1] = [0, 0, 1] + vertices[:, 2] = [ + -sqrt(1 / 10 * (5 - sqrt(5))), + -1 / 2 - 1 / (2 * sqrt(5)), + 1 / sqrt(5) + ] + vertices[:, 3] = [ + sqrt(1 / 10 * (5 - sqrt(5))), + -1 / 2 - 1 / (2 * sqrt(5)), + 1 / sqrt(5) + ] + vertices[:, 4] = [ + sqrt(1 / 10 * (5 + sqrt(5))), + 1 / 2 - 1 / (2 * sqrt(5)), + 1 / sqrt(5) + ] + vertices[:, 5] = [0, 2 / sqrt(5), 1 / sqrt(5)] + vertices[:, 6] = [ + -sqrt(1 / 10 * (5 + sqrt(5))), + 1 / 2 - 1 / (2 * sqrt(5)), + 1 / sqrt(5) + ] + vertices[:, 7] = [0, -2 / sqrt(5), -1 / sqrt(5)] + vertices[:, 8] = [ + sqrt(1 / 10 * (5 + sqrt(5))), + 1 / (2 * sqrt(5)) - 1 / 2, + -1 / sqrt(5) + ] + vertices[:, 9] = [ + sqrt(1 / 10 * (5 - sqrt(5))), + 1 / 2 + 1 / (2 * sqrt(5)), + -1 / sqrt(5) + ] + vertices[:, 10] = [ + -sqrt(1 / 10 * (5 - sqrt(5))), + 1 / 2 + 1 / (2 * sqrt(5)), + -1 / sqrt(5) + ] + vertices[:, 11] = [ + -sqrt(1 / 10 * (5 + sqrt(5))), + 1 / (2 * sqrt(5)) - 1 / 2, + -1 / sqrt(5) + ] + vertices[:, 12] = [0, 0, -1] + + return vertices * radius +end + +# Index map for the corner vertices of the triangular faces on the icosahedron (see Fig 1 and Fig 3) +# We use a counter-clockwise numbering +const icosahedron_triangle_vertices_idx_map = ([2, 3, 1], # Triangle 1 + [3, 4, 1], # Triangle 2 + [4, 5, 1], # Triangle 3 + [5, 6, 1], # Triangle 4 + [6, 2, 1], # Triangle 5 + [3, 2, 7], # Triangle 6 + [4, 3, 8], # Triangle 7 + [5, 4, 9], # Triangle 8 + [6, 5, 10], # Triangle 9 + [2, 6, 11], # Triangle 10 + [7, 8, 3], # Triangle 11 + [8, 9, 4], # Triangle 12 + [9, 10, 5], # Triangle 13 + [10, 11, 6], # Triangle 14 + [11, 7, 2], # Triangle 15 + [8, 7, 12], # Triangle 16 + [9, 8, 12], # Triangle 17 + [10, 9, 12], # Triangle 18 + [11, 10, 12], # Triangle 19 + [7, 11, 12]) + +# We use a local numbering to obtain the quad vertices of each triangular face +# +# Fig 4: Local quad vertex numbering for a triangular face (corner vertices of the triangular face in parenthesis) +# 5 (3) +# /\ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# 6/ 4\ +# /⎺⎻⎼⎽ ⎽⎼⎻⎺ \ +# / ⎺⎻⎼⎽7 ⎽⎼⎻⎺ \ +# / ⎺| \ +# / | \ +# / | \ +# / | \ +# / | \ +# / | \ +# / | \ +# 1/_________________2|__________________3\ +# (1) (2) + +# Index map for the vertices of each quad on the triangular faces of the icosahedron (see Fig 4) +const icosahedron_quad_vertices_idx_map = ([1, 2, 7, 6], # Quad 1 + [2, 3, 4, 7], # Quad 2 + [7, 4, 5, 6]) # Quad 3 +end + +# Function to initialize the p4est connectivity for the icosahedral grid. +# For reference, see Fig 1 and Fig 2 above. +function connectivity_icosahedron_2D(trees_per_face_dimension) + num_triangles = 20 + trees_per_triangle = 3 + n_cells = trees_per_face_dimension + + linear_indices = LinearIndices((n_cells, n_cells, trees_per_triangle, + num_triangles)) + + # Vertices represent the coordinates of the forest. This is used by `p4est` + # to write VTK files. + # Trixi.jl doesn't use the coordinates from `p4est`, so the vertices can be empty. + n_vertices = 0 + n_trees = n_cells * n_cells * trees_per_triangle * num_triangles + + # No corner connectivity is needed + n_corners = 0 + vertices = Trixi.C_NULL + tree_to_vertex = Trixi.C_NULL + + tree_to_tree = Array{Trixi.p4est_topidx_t, 2}(undef, 4, n_trees) + tree_to_face = Array{Int8, 2}(undef, 4, n_trees) + + # Connectivities for the first 5 (1:5) and the last 5 (16:20) triangular faces + # Notes: + # - We subtract 1 because `p4est` uses zero-based indexing + # - We use circshift to do a circular shift in the triangle list of the 5 triangles that + # are connected to each pole (useful due to periodicity) + # - We use triangle_offset_base to store the offset for the connections along the base of the + # triangles. For example, triangle 1 is connected to 6 (offset = +5), and triangle 16 + # is connected to 11 (offset = -5) + # - We use triangle_offset_13 to store the offset for the connections along the edge of the + # triangles where quads 1 and 3 are, and triangle_offset_23 to store the offset for the connections + # along the edge of the triangles where quads 1 and 2 are. For example, triangle 1 is connected to + # triangle 5 along the side where quads 1 and 3 are (triangle_offset_13 = +4), and triangle 1 is also + # connected to triangle 2 along the side where quads 2 and 3 are (triangle_offset_23 = 1) + triangle_list = Vector(101:105) + triangle_list_1 = Vector(1:5) + triangle_list_2 = Vector(16:20) + # initialize triangle_offset_base with a random value (scopes!) + triangle_offset_base = 1000 + + # Loop over the triangles + for triangle in [1:5; 16:20] + if triangle <= 5 + triangle_offset_base = 5 + triangle_list = triangle_list_1 + triangle_offset_13 = 4 + triangle_offset_23 = 1 + offset = 0 + else # triangle >=16 + triangle_offset_base = -5 + triangle_list = triangle_list_2 + triangle_offset_13 = 1 + triangle_offset_23 = 4 + offset = 15 + end + + # Quad 1 + ######## + for cell_y in 1:n_cells, cell_x in 1:n_cells + tree = linear_indices[cell_x, cell_y, 1, triangle] + + # Negative xi-direction + if cell_x > 1 # Connect to tree at the same quad + tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, 1, + triangle] - 1 + tree_to_face[1, tree] = 1 + else + tree_to_tree[1, tree] = linear_indices[end, cell_y, 2, + circshift(triangle_list, + -triangle_offset_13)[triangle - offset]] - + 1 + tree_to_face[1, tree] = 1 + end + + # Positive xi-direction + if cell_x < n_cells # Connect to tree at the same quad + tree_to_tree[2, tree] = linear_indices[cell_x + 1, cell_y, 1, + triangle] - 1 + tree_to_face[2, tree] = 0 + else + tree_to_tree[2, tree] = linear_indices[1, cell_y, 2, triangle] - 1 + tree_to_face[2, tree] = 0 + end + + # Negative eta-direction + if cell_y > 1 # Connect to tree at the same quad + tree_to_tree[3, tree] = linear_indices[cell_x, cell_y - 1, 1, + triangle] - 1 + tree_to_face[3, tree] = 3 + else + tree_to_tree[3, tree] = linear_indices[n_cells - cell_x + 1, 1, 2, + triangle + triangle_offset_base] - + 1 + tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + end + + # Positive eta-direction + if cell_y < n_cells # Connect to tree at the same quad + tree_to_tree[4, tree] = linear_indices[cell_x, cell_y + 1, 1, + triangle] - 1 + tree_to_face[4, tree] = 2 + else + tree_to_tree[4, tree] = linear_indices[1, n_cells - cell_x + 1, 3, + triangle] - 1 + tree_to_face[4, tree] = 4 # first face dimensions are oppositely oriented, add 4 + end + end + + # Quad 2 + ######## + for cell_y in 1:n_cells, cell_x in 1:n_cells + tree = linear_indices[cell_x, cell_y, 2, triangle] + + # Negative xi-direction + if cell_x > 1 # Connect to tree at the same quad + tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, 2, + triangle] - 1 + tree_to_face[1, tree] = 1 + else + tree_to_tree[1, tree] = linear_indices[end, cell_y, 1, triangle] - 1 + tree_to_face[1, tree] = 1 + end + + # Positive xi-direction + if cell_x < n_cells # Connect to tree at the same quad + tree_to_tree[2, tree] = linear_indices[cell_x + 1, cell_y, 2, + triangle] - 1 + tree_to_face[2, tree] = 0 + else + tree_to_tree[2, tree] = linear_indices[1, cell_y, 1, + circshift(triangle_list, + -triangle_offset_23)[triangle - offset]] - + 1 + tree_to_face[2, tree] = 0 + end + + # Negative eta-direction + if cell_y > 1 # Connect to tree at the same quad + tree_to_tree[3, tree] = linear_indices[cell_x, cell_y - 1, 2, + triangle] - 1 + tree_to_face[3, tree] = 3 + else + tree_to_tree[3, tree] = linear_indices[n_cells - cell_x + 1, 1, 1, + triangle + triangle_offset_base] - + 1 + tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + end + + # Positive eta-direction + if cell_y < n_cells # Connect to tree at the same quad + tree_to_tree[4, tree] = linear_indices[cell_x, cell_y + 1, 2, + triangle] - 1 + tree_to_face[4, tree] = 2 + else + tree_to_tree[4, tree] = linear_indices[cell_x, 1, 3, triangle] - 1 + tree_to_face[4, tree] = 2 + end + end + + # Quad 3 + ######## + for cell_y in 1:n_cells, cell_x in 1:n_cells + tree = linear_indices[cell_x, cell_y, 3, triangle] + + # Negative xi-direction + if cell_x > 1 # Connect to tree at the same quad + tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, 3, + triangle] - 1 + tree_to_face[1, tree] = 1 + else + tree_to_tree[1, tree] = linear_indices[n_cells - cell_y + 1, end, 1, + triangle] - 1 + tree_to_face[1, tree] = 7 # first face dimensions are oppositely oriented, add 4 + end + + # Positive xi-direction + if cell_x < n_cells # Connect to tree at the same quad + tree_to_tree[2, tree] = linear_indices[cell_x + 1, cell_y, 3, + triangle] - 1 + tree_to_face[2, tree] = 0 + else + tree_to_tree[2, tree] = linear_indices[cell_y, end, 3, + circshift(triangle_list, + -triangle_offset_23)[triangle - offset]] - + 1 + tree_to_face[2, tree] = 3 + end + + # Negative eta-direction + if cell_y > 1 # Connect to tree at the same quad + tree_to_tree[3, tree] = linear_indices[cell_x, cell_y - 1, 3, + triangle] - 1 + tree_to_face[3, tree] = 3 + else + tree_to_tree[3, tree] = linear_indices[cell_x, end, 2, triangle] - 1 + tree_to_face[3, tree] = 3 + end + + # Positive eta-direction + if cell_y < n_cells # Connect to tree at the same quad + tree_to_tree[4, tree] = linear_indices[cell_x, cell_y + 1, 3, + triangle] - 1 + tree_to_face[4, tree] = 2 + else + tree_to_tree[4, tree] = linear_indices[end, cell_x, 3, + circshift(triangle_list, + -triangle_offset_13)[triangle - offset]] - + 1 + tree_to_face[4, tree] = 1 + end + end + end + + # Connectivities for the triangular faces 6:15 + # Notes + # - We subtract 1 because `p4est` uses zero-based indexing + # - We use triangle_connectivity_n* to store the triangle connectivity along edges + # - We use circshift to do a circular shift in the triangle list of triangles 6:15 + # (useful due to periodicity) + # - We use triangle_offset_base to store the offset for the connections along the base of the + # triangles. For example, triangle 6 is connected to 1 (offset = -5), and triangle 11 + # is connected to 16 (offset = +5) + triangle_list = Vector(6:15) + # Triangle connectivity along edges with orientation north-east (/) + triangle_connectivity_ne = circshift(triangle_list, -5) + # Triangle connectivity along edges with orientation north-west (\) + triangle_connectivity_nw = [15; 11:14; 7:10; 6] + # initialize triangle_offset_base (scopes!) + triangle_offset_base = 1000 + + # Loop over the triangles + for triangle in 6:15 + if triangle <= 10 + triangle_offset_base = -5 + else # triangle >=11 + triangle_offset_base = 5 + end + + # Quad 1 + ######## + for cell_y in 1:n_cells, cell_x in 1:n_cells + tree = linear_indices[cell_x, cell_y, 1, triangle] + + # Negative xi-direction + if cell_x > 1 # Connect to tree at the same quad + tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, 1, + triangle] - 1 + tree_to_face[1, tree] = 1 + else + tree_to_tree[1, tree] = linear_indices[n_cells - cell_y + 1, end, 3, + triangle_connectivity_ne[triangle - 5]] - + 1 + tree_to_face[1, tree] = 7 # first face dimensions are oppositely oriented, add 4 + end + + # Positive xi-direction + if cell_x < n_cells # Connect to tree at the same quad + tree_to_tree[2, tree] = linear_indices[cell_x + 1, cell_y, 1, + triangle] - 1 + tree_to_face[2, tree] = 0 + else + tree_to_tree[2, tree] = linear_indices[1, cell_y, 2, triangle] - 1 + tree_to_face[2, tree] = 0 + end + + # Negative eta-direction + if cell_y > 1 # Connect to tree at the same quad + tree_to_tree[3, tree] = linear_indices[cell_x, cell_y - 1, 1, + triangle] - 1 + tree_to_face[3, tree] = 3 + else + tree_to_tree[3, tree] = linear_indices[n_cells - cell_x + 1, 1, 2, + triangle + triangle_offset_base] - + 1 + tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + end + + # Positive eta-direction + if cell_y < n_cells # Connect to tree at the same quad + tree_to_tree[4, tree] = linear_indices[cell_x, cell_y + 1, 1, + triangle] - 1 + tree_to_face[4, tree] = 2 + else + tree_to_tree[4, tree] = linear_indices[1, n_cells - cell_x + 1, 3, + triangle] - 1 + tree_to_face[4, tree] = 4 # first face dimensions are oppositely oriented, add 4 + end + end + + # Quad 2 + ######## + for cell_y in 1:n_cells, cell_x in 1:n_cells + tree = linear_indices[cell_x, cell_y, 2, triangle] + + # Negative xi-direction + if cell_x > 1 # Connect to tree at the same quad + tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, 2, + triangle] - 1 + tree_to_face[1, tree] = 1 + else + tree_to_tree[1, tree] = linear_indices[end, cell_y, 1, triangle] - 1 + tree_to_face[1, tree] = 1 + end + + # Positive xi-direction + if cell_x < n_cells # Connect to tree at the same quad + tree_to_tree[2, tree] = linear_indices[cell_x + 1, cell_y, 2, + triangle] - 1 + tree_to_face[2, tree] = 0 + else + tree_to_tree[2, tree] = linear_indices[end, n_cells - cell_y + 1, 3, + triangle_connectivity_nw[triangle - 5]] - + 1 + tree_to_face[2, tree] = 5 # first face dimensions are oppositely oriented, add 4 + end + + # Negative eta-direction + if cell_y > 1 # Connect to tree at the same quad + tree_to_tree[3, tree] = linear_indices[cell_x, cell_y - 1, 2, + triangle] - 1 + tree_to_face[3, tree] = 3 + else + tree_to_tree[3, tree] = linear_indices[n_cells - cell_x + 1, 1, 1, + triangle + triangle_offset_base] - + 1 + tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + end + + # Positive eta-direction + if cell_y < n_cells # Connect to tree at the same quad + tree_to_tree[4, tree] = linear_indices[cell_x, cell_y + 1, 2, + triangle] - 1 + tree_to_face[4, tree] = 2 + else + tree_to_tree[4, tree] = linear_indices[cell_x, 1, 3, triangle] - 1 + tree_to_face[4, tree] = 2 + end + end + + # Quad 3 + ######## + for cell_y in 1:n_cells, cell_x in 1:n_cells + tree = linear_indices[cell_x, cell_y, 3, triangle] + + # Negative xi-direction + if cell_x > 1 # Connect to tree at the same quad + tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, 3, + triangle] - 1 + tree_to_face[1, tree] = 1 + else + tree_to_tree[1, tree] = linear_indices[n_cells - cell_y + 1, end, 1, + triangle] - 1 + tree_to_face[1, tree] = 7 # first face dimensions are oppositely oriented, add 4 + end + + # Positive xi-direction + if cell_x < n_cells # Connect to tree at the same quad + tree_to_tree[2, tree] = linear_indices[cell_x + 1, cell_y, 3, + triangle] - 1 + tree_to_face[2, tree] = 0 + else + tree_to_tree[2, tree] = linear_indices[end, n_cells - cell_y + 1, 2, + triangle_connectivity_nw[triangle - 5]] - + 1 + tree_to_face[2, tree] = 5 # first face dimensions are oppositely oriented, add 4 + end + + # Negative eta-direction + if cell_y > 1 # Connect to tree at the same quad + tree_to_tree[3, tree] = linear_indices[cell_x, cell_y - 1, 3, + triangle] - 1 + tree_to_face[3, tree] = 3 + else + tree_to_tree[3, tree] = linear_indices[cell_x, end, 2, triangle] - 1 + tree_to_face[3, tree] = 3 + end + + # Positive eta-direction + if cell_y < n_cells # Connect to tree at the same quad + tree_to_tree[4, tree] = linear_indices[cell_x, cell_y + 1, 3, + triangle] - 1 + tree_to_face[4, tree] = 2 + else + tree_to_tree[4, tree] = linear_indices[1, n_cells - cell_x + 1, 1, + triangle_connectivity_ne[triangle - 5]] - + 1 + tree_to_face[4, tree] = 4 # first face dimensions are oppositely oriented, add 4 + end + end + end + + tree_to_corner = Trixi.C_NULL + # `p4est` docs: "in trivial cases it is just a pointer to a p4est_topix value of 0." + # We don't need corner connectivity, so this is a trivial case. + ctt_offset = zeros(Trixi.p4est_topidx_t, 1) + + corner_to_tree = Trixi.C_NULL + corner_to_corner = Trixi.C_NULL + + connectivity = Trixi.p4est_connectivity_new_copy(n_vertices, n_trees, n_corners, + vertices, tree_to_vertex, + tree_to_tree, tree_to_face, + tree_to_corner, ctt_offset, + corner_to_tree, corner_to_corner) + + @assert Trixi.p4est_connectivity_is_valid(connectivity) == 1 + + return connectivity +end + +# Calculate physical coordinates of each node of a 2D quad-icosahedral mesh using the +# element-local mapping from Guba et al. (see https://doi.org/10.5194/gmd-7-2803-2014, +# Appendix A). +function calc_tree_node_coordinates_quad_icosahedron_local!(node_coordinates::AbstractArray{<:Any, + 4}, + nodes, + trees_per_face_dimension, + radius) + num_triangles = 20 + trees_per_triangle = 3 + n_cells = trees_per_face_dimension + + linear_indices = LinearIndices((n_cells, n_cells, trees_per_triangle, + num_triangles)) + + triangle_vertices = Array{eltype(node_coordinates), 2}(undef, 3, 7) + icosahedron_vertices = calc_node_coordinates_icosahedron_vertices(radius) + + # Get cell length in reference mesh + dx = 2 / n_cells + dy = 2 / n_cells + + # Loop over all the triangles of the mesh + for triangle in 1:num_triangles + calc_node_coordinates_triangle_vertices!(triangle_vertices, + icosahedron_vertices, + radius, triangle) + + # Loop over each parent quad in each triangle + for local_tree in 1:3 + idx = icosahedron_quad_vertices_idx_map[local_tree] + + # Vertices of the parent quad + v1_quad = triangle_vertices[:, idx[1]] + v2_quad = triangle_vertices[:, idx[2]] + v3_quad = triangle_vertices[:, idx[3]] + v4_quad = triangle_vertices[:, idx[4]] + + # Loop over the cells/trees in each parent quad + for cell_y in 1:n_cells, cell_x in 1:n_cells + tree = linear_indices[cell_x, cell_y, local_tree, triangle] + + x_0 = -1 + (cell_x - 1) * dx + y_0 = -1 + (cell_y - 1) * dy + x_1 = -1 + cell_x * dx + y_1 = -1 + cell_y * dy + + # Obtain the coordinates of the corner nodes for the tree + v1 = local_mapping(x_0, y_0, v1_quad, v2_quad, v3_quad, v4_quad, radius) + v2 = local_mapping(x_1, y_0, v1_quad, v2_quad, v3_quad, v4_quad, radius) + v3 = local_mapping(x_1, y_1, v1_quad, v2_quad, v3_quad, v4_quad, radius) + v4 = local_mapping(x_0, y_1, v1_quad, v2_quad, v3_quad, v4_quad, radius) + + 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 +end + +function calc_node_coordinates_triangle_vertices!(triangle_vertices, + icosahedron_vertices, + radius, triangle) + # Retrieve triangle vertices + corners_triangle = icosahedron_vertices[:, + icosahedron_triangle_vertices_idx_map[triangle]] + + triangle_vertices[:, 1] = corners_triangle[:, 1] + + v2_bilinear = 0.5 * (corners_triangle[:, 1] + corners_triangle[:, 2]) + triangle_vertices[:, 2] = radius * v2_bilinear / norm(v2_bilinear) + + triangle_vertices[:, 3] = corners_triangle[:, 2] + + v4_bilinear = 0.5 * (corners_triangle[:, 2] + corners_triangle[:, 3]) + triangle_vertices[:, 4] = radius * v4_bilinear / norm(v4_bilinear) + + triangle_vertices[:, 5] = corners_triangle[:, 3] + + v6_bilinear = 0.5 * (corners_triangle[:, 1] + corners_triangle[:, 3]) + triangle_vertices[:, 6] = radius * v6_bilinear / norm(v6_bilinear) + + v7_bilinear = (corners_triangle[:, 1] + corners_triangle[:, 2] + + corners_triangle[:, 3]) / 3 + triangle_vertices[:, 7] = radius * v7_bilinear / norm(v7_bilinear) +end diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 0d6b8e3..72dd5e9 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -7,21 +7,49 @@ include("test_trixiatmo.jl") EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") -@trixiatmo_testset "Spherical advection, Cartesian weak form, LLF surface flux" begin +@trixiatmo_testset "Spherical advection (cubed sphere), Cartesian weak form, LLF surface flux" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_cubed_sphere_shell_advection.jl"), l2=[ - 0.7963221028128544, - 20.317489255709344, - 8.811382597221147, - 20.317512894705885, + 0.796321633853675, + 20.317829852384286, + 8.810001095524816, + 20.317829852393054, 0.0 ], linf=[ - 10.872101730749478, - 289.65159632622454, - 95.16499199537975, - 289.65159632621726, + 10.872101731709677, + 289.6515963524798, + 95.1288712006542, + 289.65159635247255, + 0.0 + ]) + # and small reference values + # 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 (quad icosahedron), Cartesian weak form, LLF surface flux" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_quad_icosahedron_shell_advection.jl"), + l2=[ + 0.4570227714893776, + 11.807355540221533, + 4.311881740776439, + 11.807355540221321, + 0.0 + ], + linf=[ + 13.591965583197862, + 364.7641889538827, + 93.69731833991227, + 364.76418895387906, 0.0 ]) # and small reference values @@ -39,17 +67,17 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_cubed_sphere_shell_advection.jl"), l2=[ - 0.8933384470346987, - 22.84334163690791, - 9.776600357597692, - 22.843351937152637, + 0.893342967294085, + 22.84887991899238, + 9.75885058673732, + 22.84887991899168, 0.0 ], linf=[ - 14.289456304666146, - 380.6958334078372, - 120.59259301636666, - 380.6958334078299, + 14.289456304679561, + 380.6958334082083, + 120.59259301648672, + 380.695833408201, 0.0 ], element_local_mapping=true) # and small reference values