diff --git a/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl b/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl index 1da3358..b3ed994 100644 --- a/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl +++ b/examples/elixir_shallowwater_quad_icosahedron_shell_advection.jl @@ -3,17 +3,28 @@ 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 = 3, surface_flux = flux_lax_friedrichs) +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, @@ -31,8 +42,9 @@ function source_terms_convert_to_linear_advection(u, du, x, t, end # Create a 2D cubed sphere mesh the size of the Earth -mesh = P4estMeshQuadIcosahedron2D(EARTH_RADIUS, polydeg = 3, - initial_refinement_level = 1) +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 @@ -52,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/src/meshes/p4est_icosahedron_quads.jl b/src/meshes/p4est_icosahedron_quads.jl index 902d53f..b84478c 100644 --- a/src/meshes/p4est_icosahedron_quads.jl +++ b/src/meshes/p4est_icosahedron_quads.jl @@ -2,14 +2,14 @@ #! format: noindent # Using element_local mapping -function P4estMeshQuadIcosahedron2D(radius; +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() + connectivity = connectivity_icosahedron_2D(trees_per_face_dimension) - n_trees = 60 # 20 triangles subdivided into 3 quads each + n_trees = 60 * trees_per_face_dimension^2 # 20 triangles subdivided into 3 quads each basis = LobattoLegendreBasis(RealT, polydeg) nodes = basis.nodes @@ -18,7 +18,7 @@ function P4estMeshQuadIcosahedron2D(radius; ntuple(_ -> length(nodes), 2)..., n_trees) calc_tree_node_coordinates_quad_icosahedron_local!(tree_node_coordinates, nodes, - radius) + trees_per_face_dimension, radius) p4est = Trixi.new_p4est(connectivity, initial_refinement_level) @@ -29,7 +29,7 @@ function P4estMeshQuadIcosahedron2D(radius; p4est_partition_allow_for_coarsening) end -function connectivity_icosahedron_2D() +function connectivity_icosahedron_2D(trees_per_face_dimension) # Illustration of the unfolded icosahedron with the numbering of the triangular faces # # ,'|. @@ -75,17 +75,21 @@ function connectivity_icosahedron_2D() # / / | | \ # / -------->ξ | └------->ξ \ # /__________________|___________________\ + # + # Each of those quadrilaterlas is subdivided into trees_per_face_dimension^2 trees num_triangles = 20 trees_per_triangle = 3 + n_cells = trees_per_face_dimension - linear_indices = LinearIndices((trees_per_triangle, num_triangles)) + 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 = trees_per_triangle * num_triangles + n_trees = n_cells * n_cells * trees_per_triangle * num_triangles # No corner connectivity is needed n_corners = 0 @@ -132,75 +136,159 @@ function connectivity_icosahedron_2D() # Quad 1 ######## - tree = linear_indices[1, triangle] - - # Negative xi-direction - tree_to_tree[1, tree] = linear_indices[2, - circshift(triangle_list, - -triangle_offset_13)[triangle - offset]] - - 1 - tree_to_face[1, tree] = 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 - tree_to_tree[2, tree] = linear_indices[2, triangle] - 1 - tree_to_face[2, tree] = 0 + # 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 - tree_to_tree[3, tree] = linear_indices[2, triangle + triangle_offset_base] - 1 - tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + # 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 - tree_to_tree[4, tree] = linear_indices[3, triangle] - 1 - tree_to_face[4, tree] = 4 # first face dimensions are oppositely oriented, add 4 + # 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 ######## - tree = linear_indices[2, triangle] - - # Negative xi-direction - tree_to_tree[1, tree] = linear_indices[1, triangle] - 1 - tree_to_face[1, tree] = 1 + 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 - tree_to_tree[2, tree] = linear_indices[1, - circshift(triangle_list, - -triangle_offset_23)[triangle - offset]] - - 1 - tree_to_face[2, tree] = 0 + # 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 - tree_to_tree[3, tree] = linear_indices[1, triangle + triangle_offset_base] - 1 - tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + # 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 - tree_to_tree[4, tree] = linear_indices[3, triangle] - 1 - tree_to_face[4, tree] = 2 + # 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 ######## - tree = linear_indices[3, triangle] - - # Negative xi-direction - tree_to_tree[1, tree] = linear_indices[1, triangle] - 1 - tree_to_face[1, tree] = 7 # first face dimensions are oppositely oriented, add 4 - - # Positive xi-direction - tree_to_tree[2, tree] = linear_indices[3, - circshift(triangle_list, - -triangle_offset_23)[triangle - offset]] - - 1 - tree_to_face[2, tree] = 3 - - # Negative eta-direction - tree_to_tree[3, tree] = linear_indices[2, triangle] - 1 - tree_to_face[3, tree] = 3 - - # Positive eta-direction - tree_to_tree[4, tree] = linear_indices[3, - circshift(triangle_list, - -triangle_offset_13)[triangle - offset]] - - 1 - tree_to_face[4, tree] = 1 + 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 @@ -230,71 +318,155 @@ function connectivity_icosahedron_2D() # Quad 1 ######## - tree = linear_indices[1, triangle] - - # Negative xi-direction - tree_to_tree[1, tree] = linear_indices[3, - triangle_connectivity_ne[triangle - 5]] - - 1 - tree_to_face[1, tree] = 7 # first face dimensions are oppositely oriented, add 4 + 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 - tree_to_tree[2, tree] = linear_indices[2, triangle] - 1 - tree_to_face[2, tree] = 0 + # 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 - tree_to_tree[3, tree] = linear_indices[2, triangle + triangle_offset_base] - 1 - tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + # 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 - tree_to_tree[4, tree] = linear_indices[3, triangle] - 1 - tree_to_face[4, tree] = 4 # first face dimensions are oppositely oriented, add 4 + # 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 ######## - tree = linear_indices[2, triangle] - - # Negative xi-direction - tree_to_tree[1, tree] = linear_indices[1, triangle] - 1 - tree_to_face[1, tree] = 1 + 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 - tree_to_tree[2, tree] = linear_indices[3, - triangle_connectivity_nw[triangle - 5]] - - 1 - tree_to_face[2, tree] = 5 # first face dimensions are oppositely oriented, add 4 + # 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 - tree_to_tree[3, tree] = linear_indices[1, triangle + triangle_offset_base] - 1 - tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + # 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 - tree_to_tree[4, tree] = linear_indices[3, triangle] - 1 - tree_to_face[4, tree] = 2 + # 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 ######## - tree = linear_indices[3, triangle] - - # Negative xi-direction - tree_to_tree[1, tree] = linear_indices[1, triangle] - 1 - tree_to_face[1, tree] = 7 # first face dimensions are oppositely oriented, add 4 - - # Positive xi-direction - tree_to_tree[2, tree] = linear_indices[2, - triangle_connectivity_nw[triangle - 5]] - - 1 - tree_to_face[2, tree] = 5 # first face dimensions are oppositely oriented, add 4 - - # Negative eta-direction - tree_to_tree[3, tree] = linear_indices[2, triangle] - 1 - tree_to_face[3, tree] = 3 - - # Positive eta-direction - tree_to_tree[4, tree] = linear_indices[1, - triangle_connectivity_ne[triangle - 5]] - - 1 - tree_to_face[4, tree] = 4 # first face dimensions are oppositely oriented, add 4 + 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 @@ -321,34 +493,55 @@ end # Appendix A). function calc_tree_node_coordinates_quad_icosahedron_local!(node_coordinates::AbstractArray{<:Any, 4}, - nodes, radius) + nodes, + trees_per_face_dimension, + radius) num_triangles = 20 trees_per_triangle = 3 + n_cells = trees_per_face_dimension - linear_indices = LinearIndices((trees_per_triangle, num_triangles)) + 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_icosahedron_vertices(radius) + # Get cell length in reference mesh + dx = 2 / n_cells + dy = 2 / n_cells + for triangle in 1:num_triangles calc_icosahedron_triangle_vertices!(triangle_vertices, icosahedron_vertices, radius, triangle) for local_tree in 1:3 - tree = linear_indices[local_tree, triangle] - idx = tree_vertices_idx(local_tree) # Vertices for bilinear mapping - v1 = triangle_vertices[:, idx[1]] - v2 = triangle_vertices[:, idx[2]] - v3 = triangle_vertices[:, idx[3]] - v4 = triangle_vertices[:, idx[4]] - - 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) + v1_quad = triangle_vertices[:, idx[1]] + v2_quad = triangle_vertices[:, idx[2]] + v3_quad = triangle_vertices[:, idx[3]] + v4_quad = triangle_vertices[:, idx[4]] + + 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 + + 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 diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 9f82897..68eff25 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -39,17 +39,17 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_quad_icosahedron_shell_advection.jl"), l2=[ - 1.8407728405357193, - 33.896971556215654, - 14.365741070429909, - 33.895196216486305, + 0.4570227714893776, + 11.807355540221533, + 4.311881740776439, + 11.807355540221321, 0.0 ], linf=[ - 34.15213600054404, - 601.2676604837784, - 208.32483362278194, - 601.2676604838791, + 13.591965583197862, + 364.7641889538827, + 93.69731833991227, + 364.76418895387906, 0.0 ]) # and small reference values