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 1/7] 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") From eb768d9a39b1562b006b47d64bd85acab4680f92 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 26 Jul 2024 16:51:19 +0200 Subject: [PATCH 2/7] Added element container with PtrArray for performance and modified the structure --- Project.toml | 6 +- src/TrixiAtmo.jl | 2 + .../containers_2d_manifold_in_3d.jl} | 77 +++++++++++++------ src/solvers/dgsem_p4est/dg.jl | 7 ++ src/solvers/solvers.jl | 3 +- 5 files changed, 70 insertions(+), 25 deletions(-) rename src/solvers/{dgsem_p4est_covariant/containers.jl => dgsem_p4est/containers_2d_manifold_in_3d.jl} (71%) create mode 100644 src/solvers/dgsem_p4est/dg.jl diff --git a/Project.toml b/Project.toml index cee0969..aebfcaf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "TrixiAtmo" uuid = "c9ed1054-d170-44a9-8ee2-d5566f5d1389" -authors = ["Benedict Geihe ", "Tristan Montoya ", "Hendrik Ranocha ", "Michael Schlottke-Lakemper "] +authors = ["Benedict Geihe ", "Tristan Montoya ", "Hendrik Ranocha ", "Andrés Rueda-Ramírez ", "Michael Schlottke-Lakemper "] version = "0.1.0-DEV" [deps] @@ -9,6 +9,8 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718" +StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] @@ -17,5 +19,7 @@ MuladdMacro = "0.2.2" Reexport = "1.0" Static = "0.8, 1" StaticArrays = "1" +StaticArrayInterface = "1.5.1" +StrideArrays = "0.1.28" Trixi = "0.7, 0.8" julia = "1.9" diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index c3f068e..ae5b28c 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -12,6 +12,8 @@ using Trixi using MuladdMacro: @muladd using StaticArrays: SVector using Static: True, False +using StrideArrays: PtrArray +using StaticArrayInterface: static_size using LinearAlgebra: norm using Reexport: @reexport @reexport using StaticArrays: SVector diff --git a/src/solvers/dgsem_p4est_covariant/containers.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl similarity index 71% rename from src/solvers/dgsem_p4est_covariant/containers.jl rename to src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl index 4a52f36..1006f44 100644 --- a/src/solvers/dgsem_p4est_covariant/containers.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl @@ -1,10 +1,38 @@ @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, +# New p4est element container that allows the use of a PtrArray for the contravariant_vectors +mutable struct P4estElementContainerPtrArray{NDIMS, RealT <: Real, uEltype <: Real, + NDIMSP1, + NDIMSP2, NDIMSP3, + ContravariantVectors <: + AbstractArray{RealT, NDIMSP3}} <: + Trixi.AbstractContainer + # Physical coordinates at each node + node_coordinates::Array{RealT, NDIMSP2} # [orientation, node_i, node_j, node_k, element] + # Jacobian matrix of the transformation + # [jacobian_i, jacobian_j, node_i, node_j, node_k, element] where jacobian_i is the first index of the Jacobian matrix,... + jacobian_matrix::Array{RealT, NDIMSP3} + # Contravariant vectors, scaled by J, in Kopriva's blue book called Ja^i_n (i index, n dimension) + contravariant_vectors::ContravariantVectors # [dimension, index, node_i, node_j, node_k, element] + # 1/J where J is the Jacobian determinant (determinant of Jacobian matrix) + inverse_jacobian::Array{RealT, NDIMSP1} # [node_i, node_j, node_k, element] + # Buffer for calculated surface flux + surface_flux_values::Array{uEltype, NDIMSP2} # [variable, i, j, direction, element] + + # internal `resize!`able storage + _node_coordinates::Vector{RealT} + _jacobian_matrix::Vector{RealT} + _contravariant_vectors::Vector{RealT} + _inverse_jacobian::Vector{RealT} + _surface_flux_values::Vector{uEltype} +end + +# Create element container and initialize element data. +# This function dispatches on the dimensions of the mesh and the equation (AbstractEquations{3}) +function Trixi.init_elements(mesh::Union{P4estMesh{2, RealT}, + T8codeMesh{2, RealT}}, + equations::AbstractEquations{3}, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} @@ -30,10 +58,10 @@ function Trixi.init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, _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)) + contravariant_vectors = PtrArray(pointer(_contravariant_vectors), + (Trixi.StaticInt(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), @@ -49,24 +77,27 @@ function Trixi.init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, 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) + elements = P4estElementContainerPtrArray{NDIMS, RealT, uEltype, NDIMS + 1, + NDIMS + 2, + NDIMS + 3, typeof(contravariant_vectors)}(node_coordinates, + jacobian_matrix, + contravariant_vectors, + inverse_jacobian, + surface_flux_values, + _node_coordinates, + _jacobian_matrix, + _contravariant_vectors, + _inverse_jacobian, + _surface_flux_values) + + init_elements_2d_manifold_in_3d!(elements, mesh, basis) return elements end -function init_elements_covariant!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) +function init_elements_2d_manifold_in_3d!(elements, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + basis::LobattoLegendreBasis) (; node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian) = elements @@ -118,7 +149,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates, # 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.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)) diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl new file mode 100644 index 0000000..d291bbf --- /dev/null +++ b/src/solvers/dgsem_p4est/dg.jl @@ -0,0 +1,7 @@ +# Extract contravariant vector Ja^i (i = index) as SVector +# This function dispatches on the type of contravariant_vectors +static2val(::Trixi.StaticInt{N}) where {N} = Val{N}() +@inline function get_contravariant_vector(index, contravariant_vectors::PtrArray, indices...) + SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), + static2val(static_size(contravariant_vectors, Trixi.StaticInt(1))))) +end \ No newline at end of file diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index 3cd808e..89c01e1 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -1 +1,2 @@ -include("dgsem_p4est_covariant/containers.jl") +include("dgsem_p4est/dg.jl") +include("dgsem_p4est/containers_2d_manifold_in_3d.jl") \ No newline at end of file From ca44965c311bd1ceb8be95e1a1ef514f178ea135 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 26 Jul 2024 16:56:42 +0200 Subject: [PATCH 3/7] format --- src/solvers/dgsem_p4est/dg.jl | 5 +++-- src/solvers/solvers.jl | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index d291bbf..6757aa5 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -1,7 +1,8 @@ # Extract contravariant vector Ja^i (i = index) as SVector # This function dispatches on the type of contravariant_vectors static2val(::Trixi.StaticInt{N}) where {N} = Val{N}() -@inline function get_contravariant_vector(index, contravariant_vectors::PtrArray, indices...) +@inline function get_contravariant_vector(index, contravariant_vectors::PtrArray, + indices...) SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), static2val(static_size(contravariant_vectors, Trixi.StaticInt(1))))) -end \ No newline at end of file +end diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index 89c01e1..a37105a 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -1,2 +1,2 @@ include("dgsem_p4est/dg.jl") -include("dgsem_p4est/containers_2d_manifold_in_3d.jl") \ No newline at end of file +include("dgsem_p4est/containers_2d_manifold_in_3d.jl") From 8ec3e9c8753d88190a244845de8eed2970671fa1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 26 Jul 2024 17:32:39 +0200 Subject: [PATCH 4/7] Fixed bug in the definition of init_elements and added nelements(). eltype() and ndims() functions for new struct --- .../containers_2d_manifold_in_3d.jl | 27 ++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl index 1006f44..b18708b 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl @@ -28,16 +28,29 @@ mutable struct P4estElementContainerPtrArray{NDIMS, RealT <: Real, uEltype <: Re _surface_flux_values::Vector{uEltype} end +@inline function nelements(elements::P4estElementContainerPtrArray) + 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 + } + uEltype +end + # Create element container and initialize element data. # This function dispatches on the dimensions of the mesh and the equation (AbstractEquations{3}) function Trixi.init_elements(mesh::Union{P4estMesh{2, RealT}, T8codeMesh{2, RealT}}, equations::AbstractEquations{3}, basis, - ::Type{uEltype}) where {NDIMS, RealT <: Real, + ::Type{uEltype}) where {RealT <: Real, uEltype <: Real} nelements = Trixi.ncells(mesh) + NDIMS = 2 #Dimension of the manifold ndims_spa = size(mesh.tree_node_coordinates, 1) _node_coordinates = Vector{RealT}(undef, @@ -148,11 +161,13 @@ function calc_node_coordinates_2d_shell!(node_coordinates, # 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), - Trixi.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)) + tmp1 = Trixi.StrideArray(undef, real(mesh), + Trixi.StaticInt(size(mesh.tree_node_coordinates, 1)), + Trixi.static_length(nodes), + Trixi.static_length(mesh.nodes)) + matrix1 = Trixi.StrideArray(undef, real(mesh), + Trixi.static_length(nodes), + Trixi.static_length(mesh.nodes)) matrix2 = similar(matrix1) baryweights_in = Trixi.barycentric_weights(mesh.nodes) From 4119ad71581aa62bcb856c1fc8977a5d538e9110 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 13 Aug 2024 13:59:57 +0200 Subject: [PATCH 5/7] Apply suggestions --- src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl | 2 +- src/solvers/solvers.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl index b18708b..538b895 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl @@ -28,7 +28,7 @@ mutable struct P4estElementContainerPtrArray{NDIMS, RealT <: Real, uEltype <: Re _surface_flux_values::Vector{uEltype} end -@inline function nelements(elements::P4estElementContainerPtrArray) +@inline function Trixi.nelements(elements::P4estElementContainerPtrArray) size(elements.node_coordinates, ndims(elements) + 2) end @inline Base.ndims(::P4estElementContainerPtrArray{NDIMS}) where {NDIMS} = NDIMS diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index a37105a..b9f211d 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -1,2 +1,3 @@ include("dgsem_p4est/dg.jl") +include("dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl") include("dgsem_p4est/containers_2d_manifold_in_3d.jl") From 3f62b5f12d68af5a92bfb858aa65330da68417bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 13 Aug 2024 14:53:42 +0200 Subject: [PATCH 6/7] Added Cartesian weak-form kernel for 2D manifolds in 3D --- ...xir_euler_spherical_advection_cartesian.jl | 150 ++++++++++++++++++ src/TrixiAtmo.jl | 1 + ...retization_hyperbolic_2d_manifold_in_3d.jl | 32 ++++ .../dg_2d_manifold_in_3d_cartesian.jl | 43 +++++ 4 files changed, 226 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl create mode 100644 src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl create mode 100644 src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl b/examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl new file mode 100644 index 0000000..c9ab9c0 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl @@ -0,0 +1,150 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiAtmo + +############################################################################### +# semidiscretization of the linear advection equation + +equations = CompressibleEulerEquations3D(1.4) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) + +# Initial condition for a Gaussian density profile with constant pressure +# and the velocity of a rotating solid body +function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEquations3D) + # Gaussian density + rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) + # Constant pressure + p = 1.0 + + # Spherical coordinates for the point x + if sign(x[2]) == 0.0 + signy = 1.0 + else + signy = sign(x[2]) + end + # Co-latitude + colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) + # Latitude (auxiliary variable) + lat = -colat + 0.5 * pi + # Longitude + r_xy = sqrt(x[1]^2 + x[2]^2) + if r_xy == 0.0 + phi = pi / 2 + else + phi = signy * acos(x[1] / r_xy) + end + + # Compute the velocity of a rotating solid body + # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system) + v0 = 1.0 # Velocity at the "equator" + alpha = 0.0 #pi / 4 + v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) + v_lat = -v0 * sin(phi) * sin(alpha) + + # Transform to Cartesian coordinate system + v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long + v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long + v3 = sin(colat) * v_lat + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +# Source term function to transform the Euler equations into the linear advection equations with variable advection velocity +function source_terms_convert_to_linear_advection(u, du, x, t, + equations::CompressibleEulerEquations3D) + 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] + s5 = 0.5 * (s2 * v1 + s3 * v2 + s4 * v3) - du[5] + + return SVector(0.0, s2, s3, s4, s5) +end + +# Custom RHS that applies a source term that depends on du (used to convert the 3D Euler equations into the 3D linear advection +# equations with position-dependent advection speed) +function rhs_semi_custom!(du_ode, u_ode, semi, t) + # Compute standard Trixi RHS + Trixi.rhs!(du_ode, u_ode, semi, t) + + # Now apply the custom source term + Trixi.@trixi_timeit Trixi.timer() "custom source term" begin + @unpack solver, equations, cache = semi + @unpack node_coordinates = cache.elements + + # Wrap the solution and RHS + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + + Trixi.@threaded for element in eachelement(solver, cache) + for j in eachnode(solver), i in eachnode(solver) + u_local = Trixi.get_node_vars(u, equations, solver, i, j, element) + du_local = Trixi.get_node_vars(du, equations, solver, i, j, element) + x_local = Trixi.get_node_coords(node_coordinates, equations, solver, + i, j, element) + source = source_terms_convert_to_linear_advection(u_local, du_local, + x_local, t, equations) + Trixi.add_to_node_vars!(du, source, equations, solver, i, j, element) + end + end + end +end + +initial_condition = initial_condition_advection_sphere + +mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, initial_refinement_level = 0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +ode = semidiscretize(semi, tspan) + +# Create custom discretization that runs with the custom RHS +ode_semi_custom = ODEProblem(rhs_semi_custom!, + ode.u0, + ode.tspan, + semi) + +# 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 = 10, + save_analysis = true, + extra_analysis_errors = (:conservation_error, ), + extra_analysis_integrals = (Trixi.density,)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 10, + 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_semi_custom, 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 ae5b28c..3da6d63 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -22,6 +22,7 @@ include("auxiliary/auxiliary.jl") include("equations/equations.jl") include("meshes/meshes.jl") include("solvers/solvers.jl") +include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D diff --git a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl new file mode 100644 index 0000000..c05bd16 --- /dev/null +++ b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl @@ -0,0 +1,32 @@ +# The SemidiscretizationHyperbolic constructor has been modified to remove the assertion that checks +# the compatibility between the mesh dimensionality and the equations' dimensionality. +# Instead, we now directly dispatch using the specific mesh type (P4estMesh{2}) for 2D meshes and +# AbstractEquations{3} for 3D equations. This change is necessary to support the Cartesian implementation +# of a 2D manifold embedded in a 3D space. +function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{2}, + equations::AbstractEquations{3}, + initial_condition, + solver; + source_terms = nothing, + boundary_conditions = boundary_condition_periodic, + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT = real(solver), uEltype = RealT, + initial_cache = NamedTuple()) + cache = (; Trixi.create_cache(mesh, equations, solver, RealT, uEltype)..., + initial_cache...) + _boundary_conditions = Trixi.digest_boundary_conditions(boundary_conditions, mesh, + solver, + cache) + performance_counter = Trixi.PerformanceCounter() + + SemidiscretizationHyperbolic{typeof(mesh), typeof(equations), + typeof(initial_condition), + typeof(_boundary_conditions), typeof(source_terms), + typeof(solver), typeof(cache)}(mesh, equations, + initial_condition, + _boundary_conditions, + source_terms, solver, + cache, + performance_counter) +end diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl new file mode 100644 index 0000000..8035e1a --- /dev/null +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl @@ -0,0 +1,43 @@ +# Weak-form kernel for 3D equations solved in 2D manifolds +@inline function Trixi.weak_form_kernel!(du, u, + element, + mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, + P4estMesh{2}, T8codeMesh{2}}, + nonconservative_terms::False, + equations::AbstractEquations{3}, + dg::DGSEM, cache, alpha = true) + # true * [some floating point value] == [exactly the same floating point value] + # This can (hopefully) be optimized away due to constant propagation. + @unpack derivative_dhat = dg.basis + @unpack contravariant_vectors = cache.elements + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + flux1 = flux(u_node, 1, equations) + flux2 = flux(u_node, 2, equations) + flux3 = flux(u_node, 3, equations) + + # Compute the contravariant flux by taking the scalar product of the + # first contravariant vector Ja^1 and the flux vector + Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 + for ii in eachnode(dg) + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], + contravariant_flux1, equations, dg, ii, j, + element) + end + + # Compute the contravariant flux by taking the scalar product of the + # second contravariant vector Ja^2 and the flux vector + Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 + for jj in eachnode(dg) + Trixi.multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], + contravariant_flux2, equations, dg, i, jj, + element) + end + end + + return nothing +end From bb2b1b2987c537c39c917ddf692b94f8f10f9a54 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 15 Aug 2024 14:02:35 +0200 Subject: [PATCH 7/7] format --- .../elixir_euler_spherical_advection_cartesian.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl b/examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl index c9ab9c0..4a3bb4c 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_spherical_advection_cartesian.jl @@ -99,7 +99,8 @@ end initial_condition = initial_condition_advection_sphere -mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, initial_refinement_level = 0) +mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, + initial_refinement_level = 0) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) @@ -124,7 +125,7 @@ summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results analysis_callback = AnalysisCallback(semi, interval = 10, save_analysis = true, - extra_analysis_errors = (:conservation_error, ), + extra_analysis_errors = (:conservation_error,), extra_analysis_integrals = (Trixi.density,)) # The SaveSolutionCallback allows to save the solution to a file in regular intervals