From 24a37931781a09fccb7157554c4ac10374d6a9ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 26 Jul 2024 15:07:30 +0200 Subject: [PATCH] Added containers and 2D p4est mesh from my spherical shell implementation in Trixi.jl (cherry-picked from f45378e) Co-authored-by: Tristan Montoya --- src/TrixiAtmo.jl | 2 + src/meshes/meshes.jl | 2 + src/meshes/p4est_cubed_sphere.jl | 312 ++++++++++++++++++ .../dgsem_p4est_covariant/containers.jl | 212 ++++++++++++ src/solvers/solvers.jl | 1 + 5 files changed, 529 insertions(+) create mode 100644 src/meshes/meshes.jl create mode 100644 src/meshes/p4est_cubed_sphere.jl create mode 100644 src/solvers/dgsem_p4est_covariant/containers.jl create mode 100644 src/solvers/solvers.jl diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 61b3702..c3f068e 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -18,6 +18,8 @@ using Reexport: @reexport include("auxiliary/auxiliary.jl") include("equations/equations.jl") +include("meshes/meshes.jl") +include("solvers/solvers.jl") export CompressibleMoistEulerEquations2D diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl new file mode 100644 index 0000000..7012004 --- /dev/null +++ b/src/meshes/meshes.jl @@ -0,0 +1,2 @@ +export P4estMeshCubedSphere2D +include("p4est_cubed_sphere.jl") diff --git a/src/meshes/p4est_cubed_sphere.jl b/src/meshes/p4est_cubed_sphere.jl new file mode 100644 index 0000000..07f0d9b --- /dev/null +++ b/src/meshes/p4est_cubed_sphere.jl @@ -0,0 +1,312 @@ +@muladd begin +#! format: noindent + +""" + P4estMeshCubedSphere2D(trees_per_face_dimension, radius; + polydeg, RealT=Float64, + initial_refinement_level=0, unsaved_changes=true, + p4est_partition_allow_for_coarsening=true) + +Build a "Cubed Sphere" mesh as a 2D `P4estMesh` with +`6 * trees_per_face_dimension^2` trees. + +The mesh will have no boundaries. + +# Arguments +- `trees_per_face_dimension::Integer`: the number of trees in the two local dimensions of + each face. +- `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. +""" +function P4estMeshCubedSphere2D(trees_per_face_dimension, radius; + polydeg, RealT = Float64, + initial_refinement_level = 0, unsaved_changes = true, + p4est_partition_allow_for_coarsening = true) + connectivity = connectivity_cubed_sphere_2D(trees_per_face_dimension) + + n_trees = 6 * trees_per_face_dimension^2 + + 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_cubed_sphere_2D!(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 + +# Connectivity of a 2D cubed sphere +function connectivity_cubed_sphere_2D(trees_per_face_dimension) + n_cells_x = n_cells_y = trees_per_face_dimension + + linear_indices = LinearIndices((trees_per_face_dimension, trees_per_face_dimension, + 6)) + + # 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 = 6 * n_cells_x * n_cells_y + + # 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) + + # Illustration of the local coordinates of each face. ξ and η are the first + # local coordinates of each face. The third local coordinate ζ is always + # pointing outwards, which yields a right-handed coordinate system for each face. + # ┌────────────────────────────────────────────────────┐ + # ╱│ ╱│ + # ╱ │ ξ <───┐ ╱ │ + # ╱ │ ╱ ╱ │ + # ╱ │ 4 (+y) V ╱ │ + # ╱ │ η ╱ │ + # ╱ │ ╱ │ + # ╱ │ ╱ │ + # ╱ │ ╱ │ + # ╱ │ ╱ │ + # ╱ │ 5 (-z) η ╱ │ + # ╱ │ ↑ ╱ │ + # ╱ │ │ ╱ │ + # ╱ │ ξ <───┘ ╱ │ + # ┌────────────────────────────────────────────────────┐ 2 (+x) │ + # │ │ │ │ + # │ │ │ ξ │ + # │ │ │ ↑ │ + # │ 1 (-x) │ │ │ │ + # │ │ │ │ │ + # │ ╱│ │ │ ╱ │ + # │ V │ │ │ V │ + # │ η ↓ │ │ η │ + # │ ξ └──────────────────────────────────────│─────────────┘ + # │ ╱ η 6 (+z) │ ╱ + # │ ╱ ↑ │ ╱ + # │ ╱ │ │ ╱ + # │ ╱ └───> ξ │ ╱ + # │ ╱ │ ╱ + # │ ╱ │ ╱ Global coordinates: + # │ ╱ │ ╱ y + # │ ╱ ┌───> ξ │ ╱ ↑ + # │ ╱ ╱ │ ╱ │ + # │ ╱ V 3 (-y) │ ╱ │ + # │ ╱ η │ ╱ └─────> x + # │ ╱ │ ╱ ╱ + # │╱ │╱ V + # └────────────────────────────────────────────────────┘ z + 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] + + # Subtract 1 because `p4est` uses zero-based indexing + # Negative x-direction + if cell_x > 1 # Connect to tree at the same face + tree_to_tree[1, tree] = linear_indices[cell_x - 1, cell_y, + direction] - 1 + tree_to_face[1, tree] = 1 + elseif direction == 1 # This is the -x face + target = 4 + tree_to_tree[1, tree] = linear_indices[end, cell_y, target] - 1 + tree_to_face[1, tree] = 1 + elseif direction == 2 # This is the +x face + target = 3 + tree_to_tree[1, tree] = linear_indices[end, cell_y, target] - 1 + tree_to_face[1, tree] = 1 + elseif direction == 3 # This is the -y face + target = 1 + tree_to_tree[1, tree] = linear_indices[end, cell_y, target] - 1 + tree_to_face[1, tree] = 1 + elseif direction == 4 # This is the +y face + target = 2 + tree_to_tree[1, tree] = linear_indices[end, cell_y, target] - 1 + tree_to_face[1, tree] = 1 + elseif direction == 5 # This is the -z face + target = 2 + tree_to_tree[1, tree] = linear_indices[cell_y, 1, target] - 1 + tree_to_face[1, tree] = 2 + else # direction == 6, this is the +z face + target = 1 + tree_to_tree[1, tree] = linear_indices[end - cell_y + 1, end, + target] - 1 + tree_to_face[1, tree] = 7 # first face dimensions are oppositely oriented, add 4 + end + + # Positive x-direction + if cell_x < n_cells_x # Connect to tree at the same face + tree_to_tree[2, tree] = linear_indices[cell_x + 1, cell_y, + direction] - 1 + tree_to_face[2, tree] = 0 + elseif direction == 1 # This is the -x face + target = 3 + tree_to_tree[2, tree] = linear_indices[1, cell_y, target] - 1 + tree_to_face[2, tree] = 0 + elseif direction == 2 # This is the +x face + target = 4 + tree_to_tree[2, tree] = linear_indices[1, cell_y, target] - 1 + tree_to_face[2, tree] = 0 + elseif direction == 3 # This is the -y face + target = 2 + tree_to_tree[2, tree] = linear_indices[1, cell_y, target] - 1 + tree_to_face[2, tree] = 0 + elseif direction == 4 # This is the +y face + target = 1 + tree_to_tree[2, tree] = linear_indices[1, cell_y, target] - 1 + tree_to_face[2, tree] = 0 + elseif direction == 5 # This is the -z face + target = 1 + tree_to_tree[2, tree] = linear_indices[end - cell_y + 1, 1, + target] - 1 + tree_to_face[2, tree] = 6 # first face dimensions are oppositely oriented, add 4 + else # direction == 6, this is the +z face + target = 2 + tree_to_tree[2, tree] = linear_indices[cell_y, end, target] - 1 + tree_to_face[2, tree] = 3 + end + + # Negative y-direction + if cell_y > 1 # Connect to tree at the same face + tree_to_tree[3, tree] = linear_indices[cell_x, cell_y - 1, + direction] - 1 + tree_to_face[3, tree] = 3 + elseif direction == 1 + target = 5 + tree_to_tree[3, tree] = linear_indices[end, end - cell_x + 1, + target] - 1 + tree_to_face[3, tree] = 5 # first face dimensions are oppositely oriented, add 4 + elseif direction == 2 + target = 5 + tree_to_tree[3, tree] = linear_indices[1, cell_x, target] - 1 + tree_to_face[3, tree] = 0 + elseif direction == 3 + target = 5 + tree_to_tree[3, tree] = linear_indices[end - cell_x + 1, 1, + target] - 1 + tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + elseif direction == 4 + target = 5 + tree_to_tree[3, tree] = linear_indices[cell_x, end, target] - 1 + tree_to_face[3, tree] = 3 + elseif direction == 5 + target = 3 + tree_to_tree[3, tree] = linear_indices[end - cell_x + 1, 1, + target] - 1 + tree_to_face[3, tree] = 6 # first face dimensions are oppositely oriented, add 4 + else # direction == 6 + target = 3 + tree_to_tree[3, tree] = linear_indices[cell_x, end, target] - 1 + tree_to_face[3, tree] = 3 + end + + # Positive y-direction + if cell_y < n_cells_y # Connect to tree at the same face + tree_to_tree[4, tree] = linear_indices[cell_x, cell_y + 1, + direction] - 1 + tree_to_face[4, tree] = 2 + elseif direction == 1 + target = 6 + tree_to_tree[4, tree] = linear_indices[1, end - cell_x + 1, + target] - 1 + tree_to_face[4, tree] = 4 # first face dimensions are oppositely oriented, add 4 + elseif direction == 2 + target = 6 + tree_to_tree[4, tree] = linear_indices[end, cell_x, target] - 1 + tree_to_face[4, tree] = 1 + elseif direction == 3 + target = 6 + tree_to_tree[4, tree] = linear_indices[cell_x, 1, target] - 1 + tree_to_face[4, tree] = 2 + elseif direction == 4 + target = 6 + tree_to_tree[4, tree] = linear_indices[end - cell_x + 1, end, + target] - 1 + tree_to_face[4, tree] = 7 # first face dimensions are oppositely oriented, add 4 + elseif direction == 5 + target = 4 + tree_to_tree[4, tree] = linear_indices[cell_x, 1, target] - 1 + tree_to_face[4, tree] = 2 + else # direction == 6 + target = 4 + tree_to_tree[4, tree] = linear_indices[end - cell_x + 1, end, + target] - 1 + tree_to_face[4, tree] = 7 # 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 cubed sphere mesh. +function calc_tree_node_coordinates_cubed_sphere_2D!(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 + + for j in eachindex(nodes), i in eachindex(nodes) + # node_coordinates are the mapped reference node coordinates + node_coordinates[:, i, j, tree] .= Trixi.cubed_sphere_mapping(x_offset + + dx / 2 * + nodes[i], + y_offset + + dy / 2 * + nodes[j], + z_offset, + radius, + 0, + direction) + end + end + end +end +end # @muladd diff --git a/src/solvers/dgsem_p4est_covariant/containers.jl b/src/solvers/dgsem_p4est_covariant/containers.jl new file mode 100644 index 0000000..4a52f36 --- /dev/null +++ b/src/solvers/dgsem_p4est_covariant/containers.jl @@ -0,0 +1,212 @@ +@muladd begin +#! format: noindent + +# Create element container and initialize element data +function Trixi.init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, + T8codeMesh{NDIMS, RealT}}, + equations::AbstractCovariantEquations2D, + basis, + ::Type{uEltype}) where {NDIMS, RealT <: Real, + uEltype <: Real} + nelements = Trixi.ncells(mesh) + + ndims_spa = size(mesh.tree_node_coordinates, 1) + + _node_coordinates = Vector{RealT}(undef, + ndims_spa * nnodes(basis)^NDIMS * nelements) + node_coordinates = Trixi.unsafe_wrap(Array, pointer(_node_coordinates), + (ndims_spa, + ntuple(_ -> nnodes(basis), NDIMS)..., + nelements)) + + _jacobian_matrix = Vector{RealT}(undef, + ndims_spa * NDIMS * nnodes(basis)^NDIMS * + nelements) + jacobian_matrix = Trixi.unsafe_wrap(Array, pointer(_jacobian_matrix), + (ndims_spa, NDIMS, + ntuple(_ -> nnodes(basis), NDIMS)..., + nelements)) + + _contravariant_vectors = Vector{RealT}(undef, + ndims_spa^2 * nnodes(basis)^NDIMS * + nelements) + contravariant_vectors = Trixi.unsafe_wrap(Array, pointer(_contravariant_vectors), + (ndims_spa, ndims_spa, + ntuple(_ -> nnodes(basis), NDIMS)..., + nelements)) + + _inverse_jacobian = Vector{RealT}(undef, nnodes(basis)^NDIMS * nelements) + inverse_jacobian = Trixi.unsafe_wrap(Array, pointer(_inverse_jacobian), + (ntuple(_ -> nnodes(basis), NDIMS)..., + nelements)) + + _surface_flux_values = Vector{uEltype}(undef, + nvariables(equations) * + nnodes(basis)^(NDIMS - 1) * (NDIMS * 2) * + nelements) + surface_flux_values = Trixi.unsafe_wrap(Array, pointer(_surface_flux_values), + (nvariables(equations), + ntuple(_ -> nnodes(basis), NDIMS - 1)..., + NDIMS * 2, nelements)) + + elements = Trixi.P4estElementContainer{NDIMS, RealT, uEltype, NDIMS + 1, NDIMS + 2, + NDIMS + 3}(node_coordinates, jacobian_matrix, + contravariant_vectors, + inverse_jacobian, + surface_flux_values, + _node_coordinates, + _jacobian_matrix, + _contravariant_vectors, + _inverse_jacobian, + _surface_flux_values) + + init_elements_covariant!(elements, mesh, basis) + + return elements +end + +function init_elements_covariant!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + basis::LobattoLegendreBasis) + (; node_coordinates, jacobian_matrix, + contravariant_vectors, inverse_jacobian) = elements + + calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis) + + for element in 1:Trixi.ncells(mesh) + + # Compute Jacobian matrix as usual + Trixi.calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + + # Compute contravariant vectors + calc_contravariant_vectors_2d_shell!(contravariant_vectors, + element, + jacobian_matrix, + node_coordinates, + basis) + + # Compute the inverse Jacobian as the norm of the cross product of the covariant vectors + for j in eachnode(basis), i in eachnode(basis) + inverse_jacobian[i, j, element] = 1 / + norm(Trixi.cross(jacobian_matrix[:, 1, i, + j, + element], + jacobian_matrix[:, 2, i, + j, + element])) + end + end + + return nothing +end + +# Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis +function calc_node_coordinates_2d_shell!(node_coordinates, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + basis::LobattoLegendreBasis) + # Hanging nodes will cause holes in the mesh if its polydeg is higher + # than the polydeg of the solver. + @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" + + calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis.nodes) +end + +# Interpolate tree_node_coordinates to each quadrant at the specified nodes +function calc_node_coordinates_2d_shell!(node_coordinates, + mesh::P4estMesh{2}, + nodes::AbstractVector) + # We use `StrideArray`s here since these buffers are used in performance-critical + # places and the additional information passed to the compiler makes them faster + # than native `Array`s. + tmp1 = StrideArray(undef, real(mesh), + StaticInt(size(mesh.tree_node_coordinates, 1)), + Trixi.static_length(nodes), Trixi.static_length(mesh.nodes)) + matrix1 = StrideArray(undef, real(mesh), + Trixi.static_length(nodes), Trixi.static_length(mesh.nodes)) + matrix2 = similar(matrix1) + baryweights_in = Trixi.barycentric_weights(mesh.nodes) + + # Macros from `p4est` + p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL + p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL + p4est_quadrant_len(l) = 1 << (Trixi.P4EST_MAXLEVEL - l) + + trees = Trixi.unsafe_wrap_sc(Trixi.p4est_tree_t, mesh.p4est.trees) + + for tree in eachindex(trees) + offset = trees[tree].quadrants_offset + quadrants = Trixi.unsafe_wrap_sc(Trixi.p4est_quadrant_t, trees[tree].quadrants) + + for i in eachindex(quadrants) + element = offset + i + quad = quadrants[i] + + quad_length = p4est_quadrant_len(quad.level) / p4est_root_len + + nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.x / p4est_root_len) .- 1 + nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.y / p4est_root_len) .- 1 + Trixi.polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, + baryweights_in) + Trixi.polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, + baryweights_in) + + Trixi.multiply_dimensionwise!(view(node_coordinates, :, :, :, element), + matrix1, matrix2, + view(mesh.tree_node_coordinates, :, :, :, + tree), + tmp1) + end + end + + return node_coordinates +end + +function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArray{<:Any, + 5}, + element, + jacobian_matrix, node_coordinates, + basis::LobattoLegendreBasis) + @unpack derivative_matrix = basis + + for j in eachnode(basis), i in eachnode(basis) + for n in 1:3 + # (n, m, l) cyclic + m = (n % 3) + 1 + l = ((n + 1) % 3) + 1 + + contravariant_vectors[n, 1, i, j, element] = (jacobian_matrix[m, 2, i, j, + element] * + node_coordinates[l, i, j, + element] + - + jacobian_matrix[l, 2, i, j, + element] * + node_coordinates[m, i, j, + element]) + + contravariant_vectors[n, 2, i, j, element] = (jacobian_matrix[l, 1, i, j, + element] * + node_coordinates[m, i, j, + element] + - + jacobian_matrix[m, 1, i, j, + element] * + node_coordinates[l, i, j, + element]) + + contravariant_vectors[n, 3, i, j, element] = (jacobian_matrix[m, 1, i, j, + element] * + jacobian_matrix[l, 2, i, j, + element] + - + jacobian_matrix[m, 2, i, j, + element] * + jacobian_matrix[l, 1, i, j, + element]) + end + end + + return contravariant_vectors +end +end # @muladd diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl new file mode 100644 index 0000000..3cd808e --- /dev/null +++ b/src/solvers/solvers.jl @@ -0,0 +1 @@ +include("dgsem_p4est_covariant/containers.jl")