From 02a35e56fa9e27300d5d1265e637d24f61f6c360 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 14 Nov 2023 18:37:59 +0100 Subject: [PATCH 01/19] First version of the p4est solver on a spherical shell --- src/meshes/p4est_mesh.jl | 306 ++++++++++++++++++ .../semidiscretization_hyperbolic.jl | 2 +- src/solvers/dgsem_p4est/containers.jl | 10 +- src/solvers/dgsem_p4est/containers_2d.jl | 100 +++++- src/solvers/dgsem_structured/dg.jl | 2 +- 5 files changed, 409 insertions(+), 11 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 60db285e04f..c884f8761e4 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -532,6 +532,56 @@ function P4estMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, th p4est_partition_allow_for_coarsening) end +""" + 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 `P4estMesh` with +`6 * trees_per_face_dimension^2 * layers` trees. + +The mesh will have two boundaries, `:inside` and `:outside`. + +# Arguments +- `trees_per_face_dimension::Integer`: the number of trees in the first two local dimensions of + each face. +- `radius::Integer`: the inner 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 = 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 + # Create a new p4est_connectivity that represents a structured rectangle. # Similar to p4est_connectivity_new_brick, but doesn't use Morton order. # This order makes `calc_tree_node_coordinates!` below and the calculation @@ -749,6 +799,7 @@ function connectivity_structured(n_cells_x, n_cells_y, n_cells_z, periodicity) return connectivity end +# Connectivity of a 3D cubed sphere function connectivity_cubed_sphere(trees_per_face_dimension, layers) n_cells_x = n_cells_y = trees_per_face_dimension n_cells_z = layers @@ -999,6 +1050,225 @@ function connectivity_cubed_sphere(trees_per_face_dimension, layers) return connectivity 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 = C_NULL + tree_to_vertex = C_NULL + + tree_to_tree = Array{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 = 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(p4est_topidx_t, 1) + + corner_to_tree = C_NULL + corner_to_corner = C_NULL + + connectivity = 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 p4est_connectivity_is_valid(connectivity) == 1 + + return connectivity +end + # Calculate physical coordinates of each node of a structured mesh. # This function assumes a structured mesh with trees in row order. # 2D version @@ -1198,6 +1468,42 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, end 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] .= cubed_sphere_mapping(x_offset + + dx / 2 * + nodes[i], + y_offset + + dy / 2 * + nodes[j], + z_offset, + radius, + 0, + direction) + end + end + end +end + # Map the computational coordinates xi, eta, zeta to the specified side of a cubed sphere # with the specified inner radius and thickness. function cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction) diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 50b2c21c14e..dfc63dab2c5 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -40,7 +40,7 @@ struct SemidiscretizationHyperbolic{Mesh, Equations, InitialCondition, BoundaryConditions, SourceTerms, Solver, Cache} - @assert ndims(mesh) == ndims(equations) + #@assert ndims(mesh) == ndims(equations) performance_counter = PerformanceCounter() diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 0176f5c6346..1d19b346ff3 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -87,14 +87,16 @@ function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, Re ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} nelements = ncells(mesh) - _node_coordinates = Vector{RealT}(undef, NDIMS * nnodes(basis)^NDIMS * nelements) + ndims_spa = size(mesh.tree_node_coordinates,1) + + _node_coordinates = Vector{RealT}(undef, ndims_spa * nnodes(basis)^NDIMS * nelements) node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), - (NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., + (ndims_spa, ntuple(_ -> nnodes(basis), NDIMS)..., nelements)) - _jacobian_matrix = Vector{RealT}(undef, NDIMS^2 * nnodes(basis)^NDIMS * nelements) + _jacobian_matrix = Vector{RealT}(undef, ndims_spa^2 * nnodes(basis)^NDIMS * nelements) jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix), - (NDIMS, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., + (ndims_spa, ndims_spa, ntuple(_ -> nnodes(basis), NDIMS)..., nelements)) _contravariant_vectors = similar(_jacobian_matrix) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 236d7d24c06..c163c64e160 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -13,17 +13,107 @@ function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, calc_node_coordinates!(node_coordinates, mesh, basis) - for element in 1:ncells(mesh) - calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + if size(node_coordinates,1) == 3 # 2D cubed sphere + for element in 1:ncells(mesh) + calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, node_coordinates, basis) - calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix) + #calc_contravariant_vectors_cubed_sphere!(contravariant_vectors, element, jacobian_matrix) - calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix) + calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix) + + #calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix) + calc_inverse_jacobian_cubed_sphere!(inverse_jacobian, element, jacobian_matrix, basis) + end + else + for element in 1:ncells(mesh) + calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + + calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix) + + calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix) + end end + + return nothing end +""" + calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, + node_coordinates::AbstractArray{<:Any, 4}, + basis::LobattoLegendreBasis) +Compute Jacobian matrix for cubed sphere. We compute the Jacobian components in ξ and η +direction as usual, and then compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. + +* Giraldo, F. X., Hesthaven, J. S., & Warburton, T. (2002). Nodal high-order discontinuous + Galerkin methods for the spherical shallow water equations. Journal of Computational Physics, 181(2), 499-525. +""" +function calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, + node_coordinates::AbstractArray{<:Any, 4}, + basis::LobattoLegendreBasis) + # Compute 2D Jacobian matrix as usual + calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + + # Compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. + for j in indices((jacobian_matrix, node_coordinates), (4, 3)), + i in indices((jacobian_matrix, node_coordinates), (3, 2)) + + x = node_coordinates[1, i, j, element] + y = node_coordinates[2, i, j, element] + z = node_coordinates[3, i, j, element] + theta = acos(z / sqrt(x^2 + y^2 + z^2)) + phi = sign(y) * acos(x / sqrt(x^2 + y^2)) + + jacobian_matrix[1, 3, i, j, element] = sin(theta) * cos(phi) + jacobian_matrix[2, 3, i, j, element] = sin(theta) * sin(phi) + jacobian_matrix[3, 3, i, j, element] = cos(theta) + #= jacobian_matrix[1, 3, i, j, element] = theta + jacobian_matrix[2, 3, i, j, element] = phi + jacobian_matrix[3, 3, i, j, element] = sqrt(sum(node_coordinates[:, i, j, element].^2)) =# + end +end + +# Calculate inverse Jacobian for the cubed sphere in 2D (determinant of Jacobian matrix of the mapping) in each node +function calc_inverse_jacobian_cubed_sphere!(inverse_jacobian::AbstractArray{<:Any, 3}, element, + jacobian_matrix, basis) + @turbo for j in eachnode(basis), i in eachnode(basis) + # Calculate Determinant by using Sarrus formula (about 100 times faster than LinearAlgebra.det()) + inverse_jacobian[i, j, element] = inv(jacobian_matrix[1, 1, i, j, + element] * + jacobian_matrix[2, 2, i, j, + element] * + jacobian_matrix[3, 3, i, j, element] + + jacobian_matrix[1, 2, i, j, + element] * + jacobian_matrix[2, 3, i, j, + element] * + jacobian_matrix[3, 1, i, j, element] + + jacobian_matrix[1, 3, i, j, + element] * + jacobian_matrix[2, 1, i, j, + element] * + jacobian_matrix[3, 2, i, j, element] - + jacobian_matrix[3, 1, i, j, + element] * + jacobian_matrix[2, 2, i, j, + element] * + jacobian_matrix[1, 3, i, j, element] - + jacobian_matrix[3, 2, i, j, + element] * + jacobian_matrix[2, 3, i, j, + element] * + jacobian_matrix[1, 1, i, j, element] - + jacobian_matrix[3, 3, i, j, + element] * + jacobian_matrix[2, 1, i, j, + element] * + jacobian_matrix[1, 2, i, j, element]) + end + + return inverse_jacobian +end + # Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis function calc_node_coordinates!(node_coordinates, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, @@ -43,7 +133,7 @@ function calc_node_coordinates!(node_coordinates, # places and the additional information passed to the compiler makes them faster # than native `Array`s. tmp1 = StrideArray(undef, real(mesh), - StaticInt(2), static_length(nodes), static_length(mesh.nodes)) + StaticInt(size(mesh.tree_node_coordinates,1)), static_length(nodes), static_length(mesh.nodes)) matrix1 = StrideArray(undef, real(mesh), static_length(nodes), static_length(mesh.nodes)) matrix2 = similar(matrix1) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 5cf4c4ef78c..88372693122 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -24,7 +24,7 @@ end # Extract contravariant vector Ja^i (i = index) as SVector @inline function get_contravariant_vector(index, contravariant_vectors, indices...) SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), - Val(ndims(contravariant_vectors) - 3))) + Val(size(contravariant_vectors, 1)))) end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, From 7be12aac476ada2f6faa439a65f8a8430d8d2db0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 15 Nov 2023 15:24:40 +0100 Subject: [PATCH 02/19] New strategy to compute metric terms by inheriting them from a 3D cubed sphere --- src/solvers/dgsem_p4est/containers_2d.jl | 192 +++++++++++++---------- src/solvers/dgsem_structured/dg_2d.jl | 59 ++++--- 2 files changed, 152 insertions(+), 99 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index c163c64e160..36af64e5937 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -13,16 +13,49 @@ function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, calc_node_coordinates!(node_coordinates, mesh, basis) - if size(node_coordinates,1) == 3 # 2D cubed sphere - for element in 1:ncells(mesh) - calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, node_coordinates, basis) - - #calc_contravariant_vectors_cubed_sphere!(contravariant_vectors, element, jacobian_matrix) - - calc_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix) + if size(node_coordinates,1) == 3 + # This is a 2D cubed sphere. + # We compute the metric terms by creating a 3D cubed sphere mesh with one layer and computing + # the metric terms there... We then transfer those metric terms to our spherical shell. + # TODO: Is there a more efficient way to do this?!? + mesh2 = P4estMeshCubedSphere(Int(sqrt(size(mesh.tree_node_coordinates,4)/6)), # trees_per_face_dimension read from the 2D mesh + 1, # One layer + sqrt(sum(node_coordinates[:,1,1,1].^2)), # Sphere radius as read from node 1 + 2.0, # thickness of spherical shell (2.0 seems to provide the right scaling for the jacobian_matrix. TODO: Check if this is correct.) + polydeg = size(basis.nodes, 1) - 1, # Same as current basis + initial_refinement_level = 0) #TODO: check if this needs to be changed + + # Allocate storage for 3D metric terms + # TODO: Do this without causing too many allocations! + + nelements = ncells(mesh2) + ndims_spa = size(mesh2.tree_node_coordinates,1) + + node_coordinates2 = zeros(ndims_spa, nnodes(basis), nnodes(basis), nnodes(basis), nelements) + jacobian_matrix2 = zeros(ndims_spa, ndims_spa, nnodes(basis), nnodes(basis), nnodes(basis), nelements) + contravariant_vectors2 = zeros(ndims_spa, ndims_spa, nnodes(basis), nnodes(basis), nnodes(basis), nelements) + #inverse_jacobian2 = zeros(nnodes(basis), nnodes(basis), nnodes(basis), nelements) + + # Compute 3D metric terms + calc_node_coordinates!(node_coordinates2, mesh2, basis) + for element in 1:ncells(mesh2) + calc_jacobian_matrix!(jacobian_matrix2, element, node_coordinates2, basis) + calc_contravariant_vectors!(contravariant_vectors2, element, jacobian_matrix2, + node_coordinates2, basis) + # We don't need to compute the inverse Jacobian, as we will compute the 2D inverse + # Jacobian from the 3D jacobian_matrix + end - #calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix) - calc_inverse_jacobian_cubed_sphere!(inverse_jacobian, element, jacobian_matrix, basis) + # Transfer 3D metric terms to 2D + for element in 1:ncells(mesh) + for j in eachnode(basis), i in eachnode(basis) + # We take the contravariant vectors from node k=1 + contravariant_vectors[:, :, i, j, element] = contravariant_vectors2[:, :, i, j, 1, element] + # We compute the inverse Jacobian using a cross product + inverse_jacobian[i, j, element] = 1 / norm(cross(jacobian_matrix2[:,1,i,j,1,element], jacobian_matrix2[:, 2, i, j, 1, element])) + # We don't really need jacobian_matrix here.. So we do nothing. + jacobian_matrix[:, :, i, j, element] = jacobian_matrix2[:, :, i, j, 1, element] + end end else for element in 1:ncells(mesh) @@ -39,80 +72,77 @@ function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, return nothing end -""" - calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, - node_coordinates::AbstractArray{<:Any, 4}, - basis::LobattoLegendreBasis) -Compute Jacobian matrix for cubed sphere. We compute the Jacobian components in ξ and η -direction as usual, and then compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. - -* Giraldo, F. X., Hesthaven, J. S., & Warburton, T. (2002). Nodal high-order discontinuous - Galerkin methods for the spherical shallow water equations. Journal of Computational Physics, 181(2), 499-525. -""" -function calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, - node_coordinates::AbstractArray{<:Any, 4}, - basis::LobattoLegendreBasis) - # Compute 2D Jacobian matrix as usual - calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) - - # Compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. - for j in indices((jacobian_matrix, node_coordinates), (4, 3)), - i in indices((jacobian_matrix, node_coordinates), (3, 2)) - - x = node_coordinates[1, i, j, element] - y = node_coordinates[2, i, j, element] - z = node_coordinates[3, i, j, element] - theta = acos(z / sqrt(x^2 + y^2 + z^2)) - phi = sign(y) * acos(x / sqrt(x^2 + y^2)) - - jacobian_matrix[1, 3, i, j, element] = sin(theta) * cos(phi) - jacobian_matrix[2, 3, i, j, element] = sin(theta) * sin(phi) - jacobian_matrix[3, 3, i, j, element] = cos(theta) - #= jacobian_matrix[1, 3, i, j, element] = theta - jacobian_matrix[2, 3, i, j, element] = phi - jacobian_matrix[3, 3, i, j, element] = sqrt(sum(node_coordinates[:, i, j, element].^2)) =# - end -end +#""" +# calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, +# node_coordinates::AbstractArray{<:Any, 4}, +# basis::LobattoLegendreBasis) +#Compute Jacobian matrix for cubed sphere. We compute the Jacobian components in ξ and η +#direction as usual, and then compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. +# +#* Giraldo, F. X., Hesthaven, J. S., & Warburton, T. (2002). Nodal high-order discontinuous +# Galerkin methods for the spherical shallow water equations. Journal of Computational Physics, 181(2), 499-525. +#""" +#function calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, +# node_coordinates::AbstractArray{<:Any, 4}, +# basis::LobattoLegendreBasis) +# # Compute 2D Jacobian matrix as usual +# calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) +# +# # Compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. +# for j in indices((jacobian_matrix, node_coordinates), (4, 3)), +# i in indices((jacobian_matrix, node_coordinates), (3, 2)) +# +# x = node_coordinates[1, i, j, element] +# y = node_coordinates[2, i, j, element] +# z = node_coordinates[3, i, j, element] +# theta = acos(z / sqrt(x^2 + y^2 + z^2)) +# phi = sign(y) * acos(x / sqrt(x^2 + y^2)) +# +# jacobian_matrix[1, 3, i, j, element] = sin(theta) * cos(phi) +# jacobian_matrix[2, 3, i, j, element] = sin(theta) * sin(phi) +# jacobian_matrix[3, 3, i, j, element] = cos(theta) +# end +#end # Calculate inverse Jacobian for the cubed sphere in 2D (determinant of Jacobian matrix of the mapping) in each node -function calc_inverse_jacobian_cubed_sphere!(inverse_jacobian::AbstractArray{<:Any, 3}, element, - jacobian_matrix, basis) - @turbo for j in eachnode(basis), i in eachnode(basis) - # Calculate Determinant by using Sarrus formula (about 100 times faster than LinearAlgebra.det()) - inverse_jacobian[i, j, element] = inv(jacobian_matrix[1, 1, i, j, - element] * - jacobian_matrix[2, 2, i, j, - element] * - jacobian_matrix[3, 3, i, j, element] + - jacobian_matrix[1, 2, i, j, - element] * - jacobian_matrix[2, 3, i, j, - element] * - jacobian_matrix[3, 1, i, j, element] + - jacobian_matrix[1, 3, i, j, - element] * - jacobian_matrix[2, 1, i, j, - element] * - jacobian_matrix[3, 2, i, j, element] - - jacobian_matrix[3, 1, i, j, - element] * - jacobian_matrix[2, 2, i, j, - element] * - jacobian_matrix[1, 3, i, j, element] - - jacobian_matrix[3, 2, i, j, - element] * - jacobian_matrix[2, 3, i, j, - element] * - jacobian_matrix[1, 1, i, j, element] - - jacobian_matrix[3, 3, i, j, - element] * - jacobian_matrix[2, 1, i, j, - element] * - jacobian_matrix[1, 2, i, j, element]) - end - - return inverse_jacobian -end +#function calc_inverse_jacobian_cubed_sphere!(inverse_jacobian::AbstractArray{<:Any, 3}, element, +# jacobian_matrix, basis) +# @turbo for j in eachnode(basis), i in eachnode(basis) +# # Calculate Determinant by using Sarrus formula (about 100 times faster than LinearAlgebra.det()) +# inverse_jacobian[i, j, element] = inv(jacobian_matrix[1, 1, i, j, +# element] * +# jacobian_matrix[2, 2, i, j, +# element] * +# jacobian_matrix[3, 3, i, j, element] + +# jacobian_matrix[1, 2, i, j, +# element] * +# jacobian_matrix[2, 3, i, j, +# element] * +# jacobian_matrix[3, 1, i, j, element] + +# jacobian_matrix[1, 3, i, j, +# element] * +# jacobian_matrix[2, 1, i, j, +# element] * +# jacobian_matrix[3, 2, i, j, element] - +# jacobian_matrix[3, 1, i, j, +# element] * +# jacobian_matrix[2, 2, i, j, +# element] * +# jacobian_matrix[1, 3, i, j, element] - +# jacobian_matrix[3, 2, i, j, +# element] * +# jacobian_matrix[2, 3, i, j, +# element] * +# jacobian_matrix[1, 1, i, j, element] - +# jacobian_matrix[3, 3, i, j, +# element] * +# jacobian_matrix[2, 1, i, j, +# element] * +# jacobian_matrix[1, 2, i, j, element]) +# end +# +# return inverse_jacobian +#end # Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis function calc_node_coordinates!(node_coordinates, diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 3e8ce759b30..6a6db220bd0 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -66,27 +66,50 @@ end flux1 = flux(u_node, 1, equations) flux2 = flux(u_node, 2, equations) - # Compute the contravariant flux by taking the scalar product of the - # first contravariant vector Ja^1 and the flux vector - Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) - contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 - for ii in eachnode(dg) - multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], - contravariant_flux1, equations, dg, ii, j, - element) - end + if size(contravariant_vectors,1) == 2 + # Compute the contravariant flux by taking the scalar product of the + # first contravariant vector Ja^1 and the flux vector + Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + for ii in eachnode(dg) + 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 = get_contravariant_vector(2, contravariant_vectors, i, j, element) - contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 - for jj in eachnode(dg) - multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], - contravariant_flux2, equations, dg, i, jj, - element) + # Compute the contravariant flux by taking the scalar product of the + # second contravariant vector Ja^2 and the flux vector + Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + for jj in eachnode(dg) + multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], + contravariant_flux2, equations, dg, i, jj, + element) + end + else #size(contravariant_vectors,1) == 3 + 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) + 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) + multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], + contravariant_flux2, equations, dg, i, jj, + element) + end end end - return nothing end From 27f1a811eceee12febf21622ddce2a668784a95d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 15 Nov 2023 15:38:42 +0100 Subject: [PATCH 03/19] Added elixir for advection test --- .../elixir_advection_cubed_sphere.jl | 80 +++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl diff --git a/examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl b/examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl new file mode 100644 index 00000000000..334d1061b69 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl @@ -0,0 +1,80 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +#= advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) =# + +advection_velocity = (0.0, 0.0, 0.0) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +######## +# For diffusion +#= diffusivity() = 5.0e-2 +equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations) =# +######## + +# 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 = initial_condition_convergence_test +initial_condition = initial_condition_constant + +mesh = Trixi.P4estMeshCubedSphere2D(5, 0.5, polydeg=polydeg, initial_refinement_level=0) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +####### +# For diffusion +#semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver) +####### + +# compute area of the sphere +area = zero(Float64) +for element in 1:Trixi.ncells(mesh) + for j in 1:polydeg+1, i in 1:polydeg+1 + global area += solver.basis.weights[i] * solver.basis.weights[j] / semi.cache.elements.inverse_jacobian[i, j, element] + end +end +println("Area of sphere: ", area) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +# 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=1) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval=1, + solution_variables=cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl=1.2) + +# 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() \ No newline at end of file From d23b56cd34e197978c8492c907376c576c86073c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 17 Nov 2023 16:30:53 +0100 Subject: [PATCH 04/19] Discovered a problem with allocations and implemented a temporary fix --- src/solvers/dgsem_structured/dg.jl | 1 + src/solvers/dgsem_structured/dg_2d.jl | 38 +++++++++++++++++++-------- 2 files changed, 28 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 88372693122..b9c5b70d8ea 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -22,6 +22,7 @@ function create_cache(mesh::StructuredMesh, equations::AbstractEquations, dg::DG end # Extract contravariant vector Ja^i (i = index) as SVector +# TODO: Here, size() causes a lot of allocations! Find an alternative to improve performance @inline function get_contravariant_vector(index, contravariant_vectors, indices...) SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), Val(size(contravariant_vectors, 1)))) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 6a6db220bd0..5c9ef46ec62 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -60,16 +60,17 @@ end @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) + if size(contravariant_vectors,1) == 2 + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) - if size(contravariant_vectors,1) == 2 + flux1 = flux(u_node, 1, equations) + flux2 = flux(u_node, 2, equations) # Compute the contravariant flux by taking the scalar product of the # first contravariant vector Ja^1 and the flux vector - Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + #Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + Ja11 = contravariant_vectors[1, 1, i, j, element] + Ja12 = contravariant_vectors[2, 1, i, j, element] contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 for ii in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], @@ -79,19 +80,30 @@ end # Compute the contravariant flux by taking the scalar product of the # second contravariant vector Ja^2 and the flux vector - Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + #Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + Ja21 = contravariant_vectors[1, 2, i, j, element] + Ja22 = contravariant_vectors[2, 2, i, j, element] contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 for jj in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], contravariant_flux2, equations, dg, i, jj, element) end - else #size(contravariant_vectors,1) == 3 + end + else #size(contravariant_vectors,1) == 3 + 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) + #Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + Ja11 = contravariant_vectors[1, 1, i, j, element] + Ja12 = contravariant_vectors[2, 1, i, j, element] + Ja13 = contravariant_vectors[3, 1, i, j, element] contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 for ii in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], @@ -101,7 +113,10 @@ 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) + #Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + Ja21 = contravariant_vectors[1, 2, i, j, element] + Ja22 = contravariant_vectors[2, 2, i, j, element] + Ja23 = contravariant_vectors[3, 2, i, j, element] contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 for jj in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], @@ -110,6 +125,7 @@ end end end end + return nothing end From 9338d51bf4f90fcd5182103e5c5201a48d5a1456 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 17 Nov 2023 17:02:05 +0100 Subject: [PATCH 05/19] Implemented Giraldo's metric terms and cleaned up a bit --- src/solvers/dgsem_p4est/containers.jl | 9 +- src/solvers/dgsem_p4est/containers_2d.jl | 115 ++++++++++++++--------- src/solvers/dgsem_structured/dg_2d.jl | 8 ++ 3 files changed, 85 insertions(+), 47 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 1d19b346ff3..43e9b3b32ed 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -94,14 +94,15 @@ function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, Re (ndims_spa, ntuple(_ -> nnodes(basis), NDIMS)..., nelements)) - _jacobian_matrix = Vector{RealT}(undef, ndims_spa^2 * nnodes(basis)^NDIMS * nelements) + _jacobian_matrix = Vector{RealT}(undef, ndims_spa * NDIMS * nnodes(basis)^NDIMS * nelements) jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix), - (ndims_spa, ndims_spa, ntuple(_ -> nnodes(basis), NDIMS)..., + (ndims_spa, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements)) - _contravariant_vectors = similar(_jacobian_matrix) + _contravariant_vectors = Vector{RealT}(undef, ndims_spa^2 * nnodes(basis)^NDIMS * nelements) contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors), - size(jacobian_matrix)) + (ndims_spa, ndims_spa, ntuple(_ -> nnodes(basis), NDIMS)..., + nelements)) _inverse_jacobian = Vector{RealT}(undef, nnodes(basis)^NDIMS * nelements) inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian), diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 36af64e5937..25a847a925c 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -13,48 +13,22 @@ function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, calc_node_coordinates!(node_coordinates, mesh, basis) - if size(node_coordinates,1) == 3 - # This is a 2D cubed sphere. - # We compute the metric terms by creating a 3D cubed sphere mesh with one layer and computing - # the metric terms there... We then transfer those metric terms to our spherical shell. - # TODO: Is there a more efficient way to do this?!? - mesh2 = P4estMeshCubedSphere(Int(sqrt(size(mesh.tree_node_coordinates,4)/6)), # trees_per_face_dimension read from the 2D mesh - 1, # One layer - sqrt(sum(node_coordinates[:,1,1,1].^2)), # Sphere radius as read from node 1 - 2.0, # thickness of spherical shell (2.0 seems to provide the right scaling for the jacobian_matrix. TODO: Check if this is correct.) - polydeg = size(basis.nodes, 1) - 1, # Same as current basis - initial_refinement_level = 0) #TODO: check if this needs to be changed - - # Allocate storage for 3D metric terms - # TODO: Do this without causing too many allocations! - - nelements = ncells(mesh2) - ndims_spa = size(mesh2.tree_node_coordinates,1) - - node_coordinates2 = zeros(ndims_spa, nnodes(basis), nnodes(basis), nnodes(basis), nelements) - jacobian_matrix2 = zeros(ndims_spa, ndims_spa, nnodes(basis), nnodes(basis), nnodes(basis), nelements) - contravariant_vectors2 = zeros(ndims_spa, ndims_spa, nnodes(basis), nnodes(basis), nnodes(basis), nelements) - #inverse_jacobian2 = zeros(nnodes(basis), nnodes(basis), nnodes(basis), nelements) - - # Compute 3D metric terms - calc_node_coordinates!(node_coordinates2, mesh2, basis) - for element in 1:ncells(mesh2) - calc_jacobian_matrix!(jacobian_matrix2, element, node_coordinates2, basis) - calc_contravariant_vectors!(contravariant_vectors2, element, jacobian_matrix2, - node_coordinates2, basis) - # We don't need to compute the inverse Jacobian, as we will compute the 2D inverse - # Jacobian from the 3D jacobian_matrix - end - - # Transfer 3D metric terms to 2D + if size(node_coordinates, 1) == 3 + # The mesh is a spherical shell for element in 1:ncells(mesh) + # Compute Jacobian matrix as usual + calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) + # Compute contravariant vectors with Giraldo's formula + calc_contravariant_vectors_cubed_sphere!(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) - # We take the contravariant vectors from node k=1 - contravariant_vectors[:, :, i, j, element] = contravariant_vectors2[:, :, i, j, 1, element] - # We compute the inverse Jacobian using a cross product - inverse_jacobian[i, j, element] = 1 / norm(cross(jacobian_matrix2[:,1,i,j,1,element], jacobian_matrix2[:, 2, i, j, 1, element])) - # We don't really need jacobian_matrix here.. So we do nothing. - jacobian_matrix[:, :, i, j, element] = jacobian_matrix2[:, :, i, j, 1, element] + inverse_jacobian[i, j, element] = 1 / + norm(cross(jacobian_matrix[:, 1, i, j, + element], + jacobian_matrix[:, 2, i, j, + element])) end end else @@ -67,8 +41,6 @@ function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, end end - - return nothing end @@ -163,7 +135,8 @@ function calc_node_coordinates!(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)), static_length(nodes), static_length(mesh.nodes)) + StaticInt(size(mesh.tree_node_coordinates, 1)), + static_length(nodes), static_length(mesh.nodes)) matrix1 = StrideArray(undef, real(mesh), static_length(nodes), static_length(mesh.nodes)) matrix2 = similar(matrix1) @@ -204,6 +177,62 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +# Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, +# using eq (12) of : +# Giraldo, F. X. (2001). A spectral element shallow water model on spherical geodesic grids. +# International Journal for Numerical Methods in Fluids, 35(8), 869-901. https://doi.org/10.1002/1097-0363(20010430)35:8<869::AID-FLD116>3.0.CO;2-S +function calc_contravariant_vectors_cubed_sphere!(contravariant_vectors::AbstractArray{ + <:Any, + 5 + }, + element, + jacobian_matrix, node_coordinates, + basis::LobattoLegendreBasis) + @unpack derivative_matrix = basis + + # The general form is + # Jaⁱₙ = 0.5 * ( ∇ × (Xₘ ∇ Xₗ - Xₗ ∇ Xₘ) )ᵢ where (n, m, l) cyclic and ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ + 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 + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{2}, faces, orientation, interface_id) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 5c9ef46ec62..dfcdfee793b 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -123,6 +123,14 @@ end contravariant_flux2, equations, dg, i, jj, element) end + + #Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element) + Ja31 = contravariant_vectors[1, 3, i, j, element] + Ja32 = contravariant_vectors[2, 3, i, j, element] + Ja33 = contravariant_vectors[3, 3, i, j, element] + for v in eachvariable(equations) + du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] + Ja33 * flux3[v]) + end end end From b4d0becd130310c8a697a1ed16854a306979e677 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 17 Nov 2023 17:19:19 +0100 Subject: [PATCH 06/19] Added elixir to run convergence and FSP test with the Euler equations -> We manually (hard-coded) transform the Euler equations into linear advection --- .../elixir_euler_cubed_sphere_shell.jl | 89 +++++++++++++++++++ src/solvers/dgsem_tree/dg_2d.jl | 17 +++- 2 files changed, 105 insertions(+), 1 deletion(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl new file mode 100644 index 00000000000..3502f4a067a --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -0,0 +1,89 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# 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) + +function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEquations3D) + # Constant velocity and pressure + rho = 1.0 + #if x[2] < 0 + rho = 1.0 + exp(-20*(x[1]^2 + x[3]^2)) + #end + + p = 1.0 + + colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) + lat = -colat + 0.5 * pi # Latitude, not co-latitude! + phi = sign(x[2]) * acos(x[1] / sqrt(x[1]^2 + x[2]^2)) + + v0 = 1.0 + alpha = 0.0 + v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) + v_lat = -v0 * sin(phi) * sin(alpha) + + 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 + +initial_condition = initial_condition_advection_sphere + +mesh = Trixi.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) + +# compute area of the sphere to test +area = zero(Float64) +for element in 1:Trixi.ncells(mesh) + for j in 1:polydeg+1, i in 1:polydeg+1 + global area += solver.basis.weights[i] * solver.basis.weights[j] / semi.cache.elements.inverse_jacobian[i, j, element] + end +end +println("Area of sphere: ", area) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +ode = semidiscretize(semi, tspan) + +# 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) + +# 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.8) + +# 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() \ No newline at end of file diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index c30d0a8e01a..fd88e77aa3f 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -174,7 +174,22 @@ function rhs!(du, u, t, @trixi_timeit timer() "source terms" begin calc_sources!(du, u, t, source_terms, equations, dg, cache) end - + + # Quick fix (TODO: Change) + # For the spherical shell model, transform the Euler equations into linear advection + if size(cache.elements.node_coordinates,1) == 3 && typeof(equations) == CompressibleEulerEquations3D{Float64} + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + v1 = u[2, i, j, element] / u[1, i, j, element] + v2 = u[3, i, j, element] / u[1, i, j, element] + v3 = u[4, i, j, element] / u[1, i, j, element] + du[2, i, j, element] = du[1, i, j, element] * v1 + du[3, i, j, element] = du[1, i, j, element] * v2 + du[4, i, j, element] = du[1, i, j, element] * v3 + du[5, i, j, element] = du[2, i, j, element] * v1 + du[3, i, j, element] * v2 + du[4, i, j, element] * v3 + end + end + end return nothing end From 1eb17ac7cfc8e577201310593633d44bc53a0a5d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 10:03:33 +0100 Subject: [PATCH 07/19] Fixed bugs in 2D cubed sphere elixir and in the transformation of 3D Euler into linear advection --- .../elixir_euler_cubed_sphere_shell.jl | 29 +++++++++++++------ src/solvers/dgsem_tree/dg_2d.jl | 2 +- 2 files changed, 21 insertions(+), 10 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl index 3502f4a067a..7225af72938 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -19,19 +19,28 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq #end p = 1.0 - + if sign(x[2]) == 0.0 + signy = 1.0 + else + signy = sign(x[2]) + end colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) - lat = -colat + 0.5 * pi # Latitude, not co-latitude! - phi = sign(x[2]) * acos(x[1] / sqrt(x[1]^2 + x[2]^2)) + lat = -colat + 0.5 * pi # Latitude, not co-latitude! + 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 v0 = 1.0 - alpha = 0.0 + alpha = pi / 4 v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) v_lat = -v0 * sin(phi) * sin(alpha) - 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 + 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 @@ -64,14 +73,16 @@ ode = semidiscretize(semi, tspan) 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=10, + save_analysis=true, + 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.8) +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) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index fd88e77aa3f..52967505387 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -186,7 +186,7 @@ function rhs!(du, u, t, du[2, i, j, element] = du[1, i, j, element] * v1 du[3, i, j, element] = du[1, i, j, element] * v2 du[4, i, j, element] = du[1, i, j, element] * v3 - du[5, i, j, element] = du[2, i, j, element] * v1 + du[3, i, j, element] * v2 + du[4, i, j, element] * v3 + du[5, i, j, element] = 0.5 * (du[2, i, j, element] * v1 + du[3, i, j, element] * v2 + du[4, i, j, element] * v3) end end end From e180aa43653bd561927f4d45452eb608a70adb39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 13:22:57 +0100 Subject: [PATCH 08/19] Removed hard-coded computation of rhs for the transformation of Euler into advection with variable coefficients, and improved formatting --- .../elixir_euler_cubed_sphere_shell.jl | 18 +++++++- src/meshes/p4est_mesh.jl | 43 +++++++++-------- src/solvers/dgsem_p4est/containers.jl | 21 ++++++--- src/solvers/dgsem_structured/dg_2d.jl | 23 +++++----- src/solvers/dgsem_tree/dg_2d.jl | 46 ++++++++++++------- 5 files changed, 96 insertions(+), 55 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl index 7225af72938..62b0204dce4 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -34,7 +34,7 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq end v0 = 1.0 - alpha = pi / 4 + 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) @@ -45,12 +45,26 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq 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 + initial_condition = initial_condition_advection_sphere mesh = Trixi.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) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convert_to_linear_advection) # compute area of the sphere to test area = zero(Float64) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index c884f8761e4..5a6e410e01d 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -538,15 +538,15 @@ end initial_refinement_level=0, unsaved_changes=true, p4est_partition_allow_for_coarsening=true) -Build a "Cubed Sphere" mesh as `P4estMesh` with -`6 * trees_per_face_dimension^2 * layers` trees. +Build a "Cubed Sphere" mesh as a 2D `P4estMesh` with +`6 * trees_per_face_dimension^2` trees. -The mesh will have two boundaries, `:inside` and `:outside`. +The mesh will have no boundaries. # Arguments -- `trees_per_face_dimension::Integer`: the number of trees in the first two local dimensions of +- `trees_per_face_dimension::Integer`: the number of trees in the two local dimensions of each face. -- `radius::Integer`: the inner radius of the sphere. +- `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. @@ -571,7 +571,8 @@ function P4estMeshCubedSphere2D(trees_per_face_dimension, radius; tree_node_coordinates = Array{RealT, 4}(undef, 3, ntuple(_ -> length(nodes), 2)..., n_trees) - calc_tree_node_coordinates_cubed_sphere_2D!(tree_node_coordinates, nodes, trees_per_face_dimension, radius) + calc_tree_node_coordinates_cubed_sphere_2D!(tree_node_coordinates, nodes, + trees_per_face_dimension, radius) p4est = new_p4est(connectivity, initial_refinement_level) @@ -1054,14 +1055,15 @@ end 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)) + 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 = C_NULL @@ -1469,8 +1471,11 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, end # Calculate physical coordinates of each node of a 2D cubed sphere mesh. -function calc_tree_node_coordinates_cubed_sphere_2D!(node_coordinates::AbstractArray{<:Any, 4}, - nodes, trees_per_face_dimension, radius) +function calc_tree_node_coordinates_cubed_sphere_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)) @@ -1490,15 +1495,15 @@ function calc_tree_node_coordinates_cubed_sphere_2D!(node_coordinates::AbstractA for j in eachindex(nodes), i in eachindex(nodes) # node_coordinates are the mapped reference node coordinates node_coordinates[:, i, j, tree] .= cubed_sphere_mapping(x_offset + - dx / 2 * - nodes[i], - y_offset + - dy / 2 * - nodes[j], - z_offset, - radius, - 0, - direction) + dx / 2 * + nodes[i], + y_offset + + dy / 2 * + nodes[j], + z_offset, + radius, + 0, + direction) end end end diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 43e9b3b32ed..36163293344 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -87,22 +87,29 @@ function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, Re ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} nelements = ncells(mesh) - ndims_spa = size(mesh.tree_node_coordinates,1) + ndims_spa = size(mesh.tree_node_coordinates, 1) - _node_coordinates = Vector{RealT}(undef, ndims_spa * nnodes(basis)^NDIMS * nelements) + _node_coordinates = Vector{RealT}(undef, + ndims_spa * nnodes(basis)^NDIMS * nelements) node_coordinates = 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 = Vector{RealT}(undef, + ndims_spa * NDIMS * nnodes(basis)^NDIMS * + nelements) jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix), - (ndims_spa, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., + (ndims_spa, NDIMS, + ntuple(_ -> nnodes(basis), NDIMS)..., nelements)) - _contravariant_vectors = Vector{RealT}(undef, ndims_spa^2 * nnodes(basis)^NDIMS * nelements) + _contravariant_vectors = Vector{RealT}(undef, + ndims_spa^2 * nnodes(basis)^NDIMS * + nelements) contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors), - (ndims_spa, ndims_spa, ntuple(_ -> nnodes(basis), NDIMS)..., - nelements)) + (ndims_spa, ndims_spa, + ntuple(_ -> nnodes(basis), NDIMS)..., + nelements)) _inverse_jacobian = Vector{RealT}(undef, nnodes(basis)^NDIMS * nelements) inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian), diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index dfcdfee793b..80bfd00d0e1 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -60,7 +60,7 @@ end @unpack derivative_dhat = dg.basis @unpack contravariant_vectors = cache.elements - if size(contravariant_vectors,1) == 2 + if size(contravariant_vectors, 1) == 2 for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, element) @@ -74,8 +74,8 @@ end contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 for ii in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], - contravariant_flux1, equations, dg, ii, j, - element) + contravariant_flux1, equations, dg, ii, j, + element) end # Compute the contravariant flux by taking the scalar product of the @@ -86,8 +86,8 @@ end contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 for jj in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], - contravariant_flux2, equations, dg, i, jj, - element) + contravariant_flux2, equations, dg, i, jj, + element) end end else #size(contravariant_vectors,1) == 3 @@ -107,8 +107,8 @@ end contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 for ii in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], - contravariant_flux1, equations, dg, ii, j, - element) + contravariant_flux1, equations, dg, ii, j, + element) end # Compute the contravariant flux by taking the scalar product of the @@ -120,8 +120,8 @@ end contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 for jj in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], - contravariant_flux2, equations, dg, i, jj, - element) + contravariant_flux2, equations, dg, i, jj, + element) end #Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element) @@ -129,11 +129,12 @@ end Ja32 = contravariant_vectors[2, 3, i, j, element] Ja33 = contravariant_vectors[3, 3, i, j, element] for v in eachvariable(equations) - du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] + Ja33 * flux3[v]) + du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] + + Ja33 * flux3[v]) end end end - + return nothing end diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 52967505387..8b185e8dc7b 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -172,24 +172,13 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) - end - - # Quick fix (TODO: Change) - # For the spherical shell model, transform the Euler equations into linear advection - if size(cache.elements.node_coordinates,1) == 3 && typeof(equations) == CompressibleEulerEquations3D{Float64} - @threaded for element in eachelement(dg, cache) - for j in eachnode(dg), i in eachnode(dg) - v1 = u[2, i, j, element] / u[1, i, j, element] - v2 = u[3, i, j, element] / u[1, i, j, element] - v3 = u[4, i, j, element] / u[1, i, j, element] - du[2, i, j, element] = du[1, i, j, element] * v1 - du[3, i, j, element] = du[1, i, j, element] * v2 - du[4, i, j, element] = du[1, i, j, element] * v3 - du[5, i, j, element] = 0.5 * (du[2, i, j, element] * v1 + du[3, i, j, element] * v2 + du[4, i, j, element] * v3) - end + if size(cache.elements.node_coordinates,1) == 3 + calc_extended_sources!(du, u, t, source_terms, equations, dg, cache) + else + calc_sources!(du, u, t, source_terms, equations, dg, cache) end end + return nothing end @@ -1173,4 +1162,29 @@ function calc_sources!(du, u, t, source_terms, return nothing end + +function calc_extended_sources!(du, u, t, source_terms::Nothing, + equations::AbstractEquations{3}, dg::DG, cache) + return nothing +end + +# Source term that depends on du for 3D equations +# TODO: Is there a better way to do this? +function calc_extended_sources!(du, u, t, source_terms, + equations::AbstractEquations{3}, dg::DG, cache) + @unpack node_coordinates = cache.elements + + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + du_local = get_node_vars(du, equations, dg, i, j, element) + x_local = get_node_coords(node_coordinates, equations, dg, + i, j, element) + source = source_terms(u_local, du_local, x_local, t, equations) + add_to_node_vars!(du, source, equations, dg, i, j, element) + end + end + + return nothing +end end # @muladd From 2d6fc2989c2352d018a3572c3b8a5916ea3a2b4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 13:33:57 +0100 Subject: [PATCH 09/19] Removed unused code --- .../elixir_advection_cubed_sphere.jl | 80 ------------------ src/solvers/dgsem_p4est/containers_2d.jl | 81 ++----------------- 2 files changed, 5 insertions(+), 156 deletions(-) delete mode 100644 examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl diff --git a/examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl b/examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl deleted file mode 100644 index 334d1061b69..00000000000 --- a/examples/p4est_2d_dgsem/elixir_advection_cubed_sphere.jl +++ /dev/null @@ -1,80 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the linear advection equation - -#= advection_velocity = (0.2, -0.7) -equations = LinearScalarAdvectionEquation2D(advection_velocity) =# - -advection_velocity = (0.0, 0.0, 0.0) -equations = LinearScalarAdvectionEquation3D(advection_velocity) - -######## -# For diffusion -#= diffusivity() = 5.0e-2 -equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations) =# -######## - -# 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 = initial_condition_convergence_test -initial_condition = initial_condition_constant - -mesh = Trixi.P4estMeshCubedSphere2D(5, 0.5, polydeg=polydeg, initial_refinement_level=0) - -# A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -####### -# For diffusion -#semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver) -####### - -# compute area of the sphere -area = zero(Float64) -for element in 1:Trixi.ncells(mesh) - for j in 1:polydeg+1, i in 1:polydeg+1 - global area += solver.basis.weights[i] * solver.basis.weights[j] / semi.cache.elements.inverse_jacobian[i, j, element] - end -end -println("Area of sphere: ", area) - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span from 0.0 to 1.0 -tspan = (0.0, 1.0) -ode = semidiscretize(semi, tspan) - -# 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=1) - -# The SaveSolutionCallback allows to save the solution to a file in regular intervals -save_solution = SaveSolutionCallback(interval=1, - solution_variables=cons2prim) - -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl=1.2) - -# 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() \ No newline at end of file diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 25a847a925c..ad6156ac07e 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -18,7 +18,7 @@ function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, for element in 1:ncells(mesh) # Compute Jacobian matrix as usual calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) - # Compute contravariant vectors with Giraldo's formula + # Compute contravariant vectors with Giraldo's formula (cross-product form) calc_contravariant_vectors_cubed_sphere!(contravariant_vectors, element, jacobian_matrix, node_coordinates, basis) @@ -44,78 +44,6 @@ function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, return nothing end -#""" -# calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, -# node_coordinates::AbstractArray{<:Any, 4}, -# basis::LobattoLegendreBasis) -#Compute Jacobian matrix for cubed sphere. We compute the Jacobian components in ξ and η -#direction as usual, and then compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. -# -#* Giraldo, F. X., Hesthaven, J. S., & Warburton, T. (2002). Nodal high-order discontinuous -# Galerkin methods for the spherical shallow water equations. Journal of Computational Physics, 181(2), 499-525. -#""" -#function calc_jacobian_matrix_cubed_sphere!(jacobian_matrix, element, -# node_coordinates::AbstractArray{<:Any, 4}, -# basis::LobattoLegendreBasis) -# # Compute 2D Jacobian matrix as usual -# calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates, basis) -# -# # Compute third component (dx⃗/dζ) analytically as (dx⃗/dr). See, e.g. -# for j in indices((jacobian_matrix, node_coordinates), (4, 3)), -# i in indices((jacobian_matrix, node_coordinates), (3, 2)) -# -# x = node_coordinates[1, i, j, element] -# y = node_coordinates[2, i, j, element] -# z = node_coordinates[3, i, j, element] -# theta = acos(z / sqrt(x^2 + y^2 + z^2)) -# phi = sign(y) * acos(x / sqrt(x^2 + y^2)) -# -# jacobian_matrix[1, 3, i, j, element] = sin(theta) * cos(phi) -# jacobian_matrix[2, 3, i, j, element] = sin(theta) * sin(phi) -# jacobian_matrix[3, 3, i, j, element] = cos(theta) -# end -#end - -# Calculate inverse Jacobian for the cubed sphere in 2D (determinant of Jacobian matrix of the mapping) in each node -#function calc_inverse_jacobian_cubed_sphere!(inverse_jacobian::AbstractArray{<:Any, 3}, element, -# jacobian_matrix, basis) -# @turbo for j in eachnode(basis), i in eachnode(basis) -# # Calculate Determinant by using Sarrus formula (about 100 times faster than LinearAlgebra.det()) -# inverse_jacobian[i, j, element] = inv(jacobian_matrix[1, 1, i, j, -# element] * -# jacobian_matrix[2, 2, i, j, -# element] * -# jacobian_matrix[3, 3, i, j, element] + -# jacobian_matrix[1, 2, i, j, -# element] * -# jacobian_matrix[2, 3, i, j, -# element] * -# jacobian_matrix[3, 1, i, j, element] + -# jacobian_matrix[1, 3, i, j, -# element] * -# jacobian_matrix[2, 1, i, j, -# element] * -# jacobian_matrix[3, 2, i, j, element] - -# jacobian_matrix[3, 1, i, j, -# element] * -# jacobian_matrix[2, 2, i, j, -# element] * -# jacobian_matrix[1, 3, i, j, element] - -# jacobian_matrix[3, 2, i, j, -# element] * -# jacobian_matrix[2, 3, i, j, -# element] * -# jacobian_matrix[1, 1, i, j, element] - -# jacobian_matrix[3, 3, i, j, -# element] * -# jacobian_matrix[2, 1, i, j, -# element] * -# jacobian_matrix[1, 2, i, j, element]) -# end -# -# return inverse_jacobian -#end - # Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis function calc_node_coordinates!(node_coordinates, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, @@ -180,7 +108,10 @@ end # Calculate contravariant vectors, multiplied by the Jacobian determinant J of the transformation mapping, # using eq (12) of : # Giraldo, F. X. (2001). A spectral element shallow water model on spherical geodesic grids. -# International Journal for Numerical Methods in Fluids, 35(8), 869-901. https://doi.org/10.1002/1097-0363(20010430)35:8<869::AID-FLD116>3.0.CO;2-S +# International Journal for Numerical Methods in Fluids, 35(8), 869-901. +# https://doi.org/10.1002/1097-0363(20010430)35:8<869::AID-FLD116>3.0.CO;2-S +# This is nothing but the cross-product form, but we end up with three contravariant vectors +# because there are three space dimensions. function calc_contravariant_vectors_cubed_sphere!(contravariant_vectors::AbstractArray{ <:Any, 5 @@ -190,8 +121,6 @@ function calc_contravariant_vectors_cubed_sphere!(contravariant_vectors::Abstrac basis::LobattoLegendreBasis) @unpack derivative_matrix = basis - # The general form is - # Jaⁱₙ = 0.5 * ( ∇ × (Xₘ ∇ Xₗ - Xₗ ∇ Xₘ) )ᵢ where (n, m, l) cyclic and ∇ = (∂/∂ξ, ∂/∂η, ∂/∂ζ)ᵀ for j in eachnode(basis), i in eachnode(basis) for n in 1:3 # (n, m, l) cyclic From 3715ea8d35fcde5ab65c0cd8d3329720a084bbb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 13:46:06 +0100 Subject: [PATCH 10/19] Format elixir --- .../elixir_euler_cubed_sphere_shell.jl | 47 ++++++++++--------- 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl index 62b0204dce4..cd4a7d30cb3 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -9,13 +9,13 @@ 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) +solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEquations3D) # Constant velocity and pressure rho = 1.0 #if x[2] < 0 - rho = 1.0 + exp(-20*(x[1]^2 + x[3]^2)) + rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) #end p = 1.0 @@ -28,9 +28,9 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq lat = -colat + 0.5 * pi # Latitude, not co-latitude! r_xy = sqrt(x[1]^2 + x[2]^2) if r_xy == 0.0 - phi = pi/2 + phi = pi / 2 else - phi = signy * acos(x[1] / r_xy) + phi = signy * acos(x[1] / r_xy) end v0 = 1.0 @@ -46,11 +46,12 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq 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) +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] @@ -61,16 +62,18 @@ end initial_condition = initial_condition_advection_sphere -mesh = Trixi.P4estMeshCubedSphere2D(5, 1.0, polydeg=polydeg, initial_refinement_level=0) +mesh = Trixi.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, source_terms = source_terms_convert_to_linear_advection) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convert_to_linear_advection) # compute area of the sphere to test area = zero(Float64) for element in 1:Trixi.ncells(mesh) - for j in 1:polydeg+1, i in 1:polydeg+1 - global area += solver.basis.weights[i] * solver.basis.weights[j] / semi.cache.elements.inverse_jacobian[i, j, element] + for j in 1:(polydeg + 1), i in 1:(polydeg + 1) + global area += solver.basis.weights[i] * solver.basis.weights[j] / + semi.cache.elements.inverse_jacobian[i, j, element] end end println("Area of sphere: ", area) @@ -87,28 +90,28 @@ ode = semidiscretize(semi, tspan) 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_integrals = (Trixi.density, )) +analysis_callback = AnalysisCallback(semi, interval = 10, + save_analysis = true, + 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) +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) +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) - +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); +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() \ No newline at end of file +summary_callback() From b30fe7dc615191c72ee599f2da7a7570975cc61e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 14:12:45 +0100 Subject: [PATCH 11/19] Added test. Allocations are failing due to function `get_contravariant_vector` --- .../elixir_euler_cubed_sphere_shell.jl | 31 +++++++++---------- test/test_p4est_2d.jl | 24 ++++++++++++++ 2 files changed, 38 insertions(+), 17 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl index cd4a7d30cb3..83035397039 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -11,21 +11,25 @@ equations = CompressibleEulerEquations3D(1.4) 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) - # Constant velocity and pressure - rho = 1.0 - #if x[2] < 0 + # Gaussian density rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) - #end - + # 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)) - lat = -colat + 0.5 * pi # Latitude, not co-latitude! + # 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 @@ -33,11 +37,14 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq phi = signy * acos(x[1] / r_xy) end - v0 = 1.0 + # 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 @@ -68,16 +75,6 @@ mesh = Trixi.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, initial_refinemen semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convert_to_linear_advection) -# compute area of the sphere to test -area = zero(Float64) -for element in 1:Trixi.ncells(mesh) - for j in 1:(polydeg + 1), i in 1:(polydeg + 1) - global area += solver.basis.weights[i] * solver.basis.weights[j] / - semi.cache.elements.inverse_jacobian[i, j, element] - end -end -println("Area of sphere: ", area) - ############################################################################### # ODE solvers, callbacks etc. diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index db34aecc168..453c79fa3c5 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -365,6 +365,30 @@ end end end +@trixi_testset "elixir_euler_cubed_sphere_shell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_cubed_sphere_shell.jl"), + l2=[ + 0.005386774714265701, 0.0053034549371091975, + 0.0009984808604061753, 0.0, 0.0011482642682649322, + ], + linf=[ + 0.05757879060005222, 0.05734888774521418, + 0.0099936366211697, 0.0, 0.021317676978318545, + ], + tspan=(0.0, 0.01), + skip_coverage=true) + if @isdefined sol # Skipped in coverage run + # 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 Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + @trixi_testset "elixir_eulergravity_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"), l2=[ From 2942c796dd7d7542ccadb540849a12ea6582e781 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 14:15:21 +0100 Subject: [PATCH 12/19] Format --- examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl index 83035397039..c5e7f6645b0 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -18,7 +18,7 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq 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 @@ -28,7 +28,7 @@ function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEq # Co-latitude colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) # Latitude (auxiliary variable) - lat = -colat + 0.5 * pi + lat = -colat + 0.5 * pi # Longitude r_xy = sqrt(x[1]^2 + x[2]^2) if r_xy == 0.0 From b89e333af4fd1b9c84a3413511c6939e4f7d8458 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 14:18:37 +0100 Subject: [PATCH 13/19] Format --- src/solvers/dgsem_tree/dg_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 9173bd7fdb8..f1eb927e397 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -172,13 +172,13 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - if size(cache.elements.node_coordinates,1) == 3 + if size(cache.elements.node_coordinates, 1) == 3 calc_extended_sources!(du, u, t, source_terms, equations, dg, cache) else calc_sources!(du, u, t, source_terms, equations, dg, cache) end end - + return nothing end From a18fcd623ff5d51b4dc808a319c6ba109cd20bcc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 5 Dec 2023 07:16:19 +0100 Subject: [PATCH 14/19] Apply custom source term (that depends on du) with a custom RHS --- .../elixir_euler_cubed_sphere_shell.jl | 40 +++++++++++++++++-- src/solvers/dgsem_tree/dg_2d.jl | 31 +------------- 2 files changed, 38 insertions(+), 33 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl index c5e7f6645b0..28deb4964e1 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -67,13 +67,41 @@ function source_terms_convert_to_linear_advection(u, du, x, t, 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 = Trixi.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, - source_terms = source_terms_convert_to_linear_advection) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) ############################################################################### # ODE solvers, callbacks etc. @@ -82,6 +110,12 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, 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() @@ -106,7 +140,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +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); diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 6f078431142..547ed352ef3 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -172,11 +172,7 @@ function rhs!(du, u, t, # Calculate source terms @trixi_timeit timer() "source terms" begin - if size(cache.elements.node_coordinates, 1) == 3 - calc_extended_sources!(du, u, t, source_terms, equations, dg, cache) - else - calc_sources!(du, u, t, source_terms, equations, dg, cache) - end + calc_sources!(du, u, t, source_terms, equations, dg, cache) end return nothing @@ -1169,29 +1165,4 @@ function calc_sources!(du, u, t, source_terms, return nothing end - -function calc_extended_sources!(du, u, t, source_terms::Nothing, - equations::AbstractEquations{3}, dg::DG, cache) - return nothing -end - -# Source term that depends on du for 3D equations -# TODO: Is there a better way to do this? -function calc_extended_sources!(du, u, t, source_terms, - equations::AbstractEquations{3}, dg::DG, cache) - @unpack node_coordinates = cache.elements - - @threaded for element in eachelement(dg, cache) - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - du_local = get_node_vars(du, equations, dg, i, j, element) - x_local = get_node_coords(node_coordinates, equations, dg, - i, j, element) - source = source_terms(u_local, du_local, x_local, t, equations) - add_to_node_vars!(du, source, equations, dg, i, j, element) - end - end - - return nothing -end end # @muladd From d2b5294e538531e233db39c64c0c8e121e3a9d15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 5 Dec 2023 07:37:59 +0100 Subject: [PATCH 15/19] Using the number of dimensions of the equations to dispatch the volume integral of p4est 2D --- src/solvers/dgsem_structured/dg_2d.jl | 146 ++++++++++++++------------ 1 file changed, 80 insertions(+), 66 deletions(-) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 6d058af9ee2..d2bdffeeb8a 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -60,85 +60,99 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 element, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, - nonconservative_terms::False, equations, + nonconservative_terms::False, + equations::AbstractEquations{2}, 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 - if size(contravariant_vectors, 1) == 2 - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, element) + 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) - # Compute the contravariant flux by taking the scalar product of the - # first contravariant vector Ja^1 and the flux vector - #Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) - Ja11 = contravariant_vectors[1, 1, i, j, element] - Ja12 = contravariant_vectors[2, 1, i, j, element] - contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 - for ii in eachnode(dg) - multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], - contravariant_flux1, equations, dg, ii, j, - element) - end + flux1 = flux(u_node, 1, equations) + flux2 = flux(u_node, 2, equations) + # Compute the contravariant flux by taking the scalar product of the + # first contravariant vector Ja^1 and the flux vector + #Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + Ja11 = contravariant_vectors[1, 1, i, j, element] + Ja12 = contravariant_vectors[2, 1, i, j, element] + contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + for ii in eachnode(dg) + 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 = get_contravariant_vector(2, contravariant_vectors, i, j, element) - Ja21 = contravariant_vectors[1, 2, i, j, element] - Ja22 = contravariant_vectors[2, 2, i, j, element] - contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 - for jj in eachnode(dg) - multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], - contravariant_flux2, equations, dg, i, jj, - element) - end + # Compute the contravariant flux by taking the scalar product of the + # second contravariant vector Ja^2 and the flux vector + #Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + Ja21 = contravariant_vectors[1, 2, i, j, element] + Ja22 = contravariant_vectors[2, 2, i, j, element] + contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + for jj in eachnode(dg) + multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], + contravariant_flux2, equations, dg, i, jj, + element) end - else #size(contravariant_vectors,1) == 3 - for j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, element) + end - flux1 = flux(u_node, 1, equations) - flux2 = flux(u_node, 2, equations) - flux3 = flux(u_node, 3, equations) + return nothing +end - # 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) - Ja11 = contravariant_vectors[1, 1, i, j, element] - Ja12 = contravariant_vectors[2, 1, i, j, element] - Ja13 = contravariant_vectors[3, 1, i, j, element] - contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 - for ii in eachnode(dg) - multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], - contravariant_flux1, equations, dg, ii, j, - element) - end +@inline function 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 - # 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) - Ja21 = contravariant_vectors[1, 2, i, j, element] - Ja22 = contravariant_vectors[2, 2, i, j, element] - Ja23 = contravariant_vectors[3, 2, i, j, element] - contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 - for jj in eachnode(dg) - multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], - contravariant_flux2, equations, dg, i, jj, - element) - end + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) - #Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element) - Ja31 = contravariant_vectors[1, 3, i, j, element] - Ja32 = contravariant_vectors[2, 3, i, j, element] - Ja33 = contravariant_vectors[3, 3, i, j, element] - for v in eachvariable(equations) - du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] + - Ja33 * flux3[v]) - end + 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) + Ja11 = contravariant_vectors[1, 1, i, j, element] + Ja12 = contravariant_vectors[2, 1, i, j, element] + Ja13 = contravariant_vectors[3, 1, i, j, element] + contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 + for ii in eachnode(dg) + 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) + Ja21 = contravariant_vectors[1, 2, i, j, element] + Ja22 = contravariant_vectors[2, 2, i, j, element] + Ja23 = contravariant_vectors[3, 2, i, j, element] + contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 + for jj in eachnode(dg) + multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], + contravariant_flux2, equations, dg, i, jj, + element) + end + + #Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element) + Ja31 = contravariant_vectors[1, 3, i, j, element] + Ja32 = contravariant_vectors[2, 3, i, j, element] + Ja33 = contravariant_vectors[3, 3, i, j, element] + for v in eachvariable(equations) + du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] + + Ja33 * flux3[v]) end end From 245b20d93f684c355ae4e363a87cbb9e9039832c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 6 Dec 2023 15:39:51 +0100 Subject: [PATCH 16/19] Use PtrArray for contravariant_vectors --- src/Trixi.jl | 2 +- src/solvers/dgsem_p4est/containers.jl | 10 +++++----- src/solvers/dgsem_structured/dg.jl | 4 ++-- src/solvers/dgsem_structured/dg_2d.jl | 23 +++++------------------ 4 files changed, 13 insertions(+), 26 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index b8110cf5bdd..bf7e1d0534a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -48,7 +48,7 @@ using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace using IfElse: ifelse using LinearMaps: LinearMap using LoopVectorization: LoopVectorization, @turbo, indices -using StaticArrayInterface: static_length # used by LoopVectorization +using StaticArrayInterface: static_length, static_size # used by LoopVectorization using MuladdMacro: @muladd using Octavian: Octavian, matmul! using Polyester: Polyester, @batch # You know, the cheapest threads you can find... diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index b093cf5c018..b8f38b57cda 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -13,7 +13,7 @@ mutable struct P4estElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIM # [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::Array{RealT, NDIMSP3} # [dimension, index, node_i, node_j, node_k, element] + contravariant_vectors::PtrArray{RealT, NDIMSP3} # [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 @@ -106,10 +106,10 @@ function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, Re _contravariant_vectors = Vector{RealT}(undef, ndims_spa^2 * nnodes(basis)^NDIMS * nelements) - contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors), - (ndims_spa, ndims_spa, - ntuple(_ -> nnodes(basis), NDIMS)..., - nelements)) + contravariant_vectors = PtrArray(pointer(_contravariant_vectors), + (StaticInt(ndims_spa), ndims_spa, + ntuple(_ -> nnodes(basis), NDIMS)..., + nelements)) _inverse_jacobian = Vector{RealT}(undef, nnodes(basis)^NDIMS * nelements) inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian), diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index b9c5b70d8ea..7e920adfd50 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -22,12 +22,12 @@ function create_cache(mesh::StructuredMesh, equations::AbstractEquations, dg::DG end # Extract contravariant vector Ja^i (i = index) as SVector -# TODO: Here, size() causes a lot of allocations! Find an alternative to improve performance @inline function get_contravariant_vector(index, contravariant_vectors, indices...) SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), - Val(size(contravariant_vectors, 1)))) + Val(static_size(contravariant_vectors, 1)))) end + @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionPeriodic, diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index d2bdffeeb8a..625b2186a54 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -75,9 +75,7 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 flux2 = flux(u_node, 2, equations) # Compute the contravariant flux by taking the scalar product of the # first contravariant vector Ja^1 and the flux vector - #Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) - Ja11 = contravariant_vectors[1, 1, i, j, element] - Ja12 = contravariant_vectors[2, 1, i, j, element] + Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 for ii in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], @@ -87,9 +85,7 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 # Compute the contravariant flux by taking the scalar product of the # second contravariant vector Ja^2 and the flux vector - #Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element) - Ja21 = contravariant_vectors[1, 2, i, j, element] - Ja22 = contravariant_vectors[2, 2, i, j, element] + Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element) contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 for jj in eachnode(dg) multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], @@ -122,10 +118,7 @@ end # 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) - Ja11 = contravariant_vectors[1, 1, i, j, element] - Ja12 = contravariant_vectors[2, 1, i, j, element] - Ja13 = contravariant_vectors[3, 1, i, j, element] + 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) multiply_add_to_node_vars!(du, alpha * derivative_dhat[ii, i], @@ -135,10 +128,7 @@ 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) - Ja21 = contravariant_vectors[1, 2, i, j, element] - Ja22 = contravariant_vectors[2, 2, i, j, element] - Ja23 = contravariant_vectors[3, 2, i, j, element] + 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) multiply_add_to_node_vars!(du, alpha * derivative_dhat[jj, j], @@ -146,10 +136,7 @@ end element) end - #Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element) - Ja31 = contravariant_vectors[1, 3, i, j, element] - Ja32 = contravariant_vectors[2, 3, i, j, element] - Ja33 = contravariant_vectors[3, 3, i, j, element] + Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, element) for v in eachvariable(equations) du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] + Ja33 * flux3[v]) From 32b6836652612320f671da43188ace035714f501 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 7 Dec 2023 14:50:54 +0100 Subject: [PATCH 17/19] Added function static2val --- src/solvers/dgsem_structured/dg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 7e920adfd50..47eecc7bf8f 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -22,12 +22,12 @@ function create_cache(mesh::StructuredMesh, equations::AbstractEquations, dg::DG end # Extract contravariant vector Ja^i (i = index) as SVector +static2val(::StaticInt{N}) where {N} = Val{N}() @inline function get_contravariant_vector(index, contravariant_vectors, indices...) SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), - Val(static_size(contravariant_vectors, 1)))) + static2val(static_size(contravariant_vectors, StaticInt(1))))) end - @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionPeriodic, From d959d77feda318438abd4806dff5859718e415a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 8 Dec 2023 10:51:45 +0100 Subject: [PATCH 18/19] Initialize contravariant_vectors with full type information --- src/solvers/dgsem_p4est/containers.jl | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index b8f38b57cda..ede056c0a51 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -6,14 +6,14 @@ #! format: noindent mutable struct P4estElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, - NDIMSP2, NDIMSP3} <: AbstractContainer + NDIMSP2, NDIMSP3, R, S, X, O} <: 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::PtrArray{RealT, NDIMSP3} # [dimension, index, node_i, node_j, node_k, element] + contravariant_vectors::PtrArray{RealT, NDIMSP3, R, S, X, O} # [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 @@ -125,12 +125,20 @@ function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, Re NDIMS * 2, nelements)) elements = 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) + NDIMS + 3, ntuple(@inline(i->i), NDIMS + 3), + Tuple{StaticInt{ndims_spa}, + ntuple(i -> Int64, NDIMS + 2)... + }, NTuple{NDIMS + 3, Nothing}, + NTuple{NDIMS + 3, StaticInt{1}}}(node_coordinates, + jacobian_matrix, + contravariant_vectors, + inverse_jacobian, + surface_flux_values, + _node_coordinates, + _jacobian_matrix, + _contravariant_vectors, + _inverse_jacobian, + _surface_flux_values) init_elements!(elements, mesh, basis) return elements From 500db958c8fad5777ac9e3e7e605510484a6b147 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 8 Dec 2023 11:08:59 +0100 Subject: [PATCH 19/19] Improved initialization of P4estElementContainer --- src/solvers/dgsem_p4est/containers.jl | 30 +++++++++++++-------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index ede056c0a51..c8cd69c66b5 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -6,14 +6,16 @@ #! format: noindent mutable struct P4estElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, - NDIMSP2, NDIMSP3, R, S, X, O} <: AbstractContainer + NDIMSP2, NDIMSP3, + ContravariantVectors <: + AbstractArray{RealT, NDIMSP3}} <: 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::PtrArray{RealT, NDIMSP3, R, S, X, O} # [dimension, index, node_i, node_j, node_k, element] + 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 @@ -125,20 +127,16 @@ function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, Re NDIMS * 2, nelements)) elements = P4estElementContainer{NDIMS, RealT, uEltype, NDIMS + 1, NDIMS + 2, - NDIMS + 3, ntuple(@inline(i->i), NDIMS + 3), - Tuple{StaticInt{ndims_spa}, - ntuple(i -> Int64, NDIMS + 2)... - }, NTuple{NDIMS + 3, Nothing}, - NTuple{NDIMS + 3, StaticInt{1}}}(node_coordinates, - jacobian_matrix, - contravariant_vectors, - inverse_jacobian, - surface_flux_values, - _node_coordinates, - _jacobian_matrix, - _contravariant_vectors, - _inverse_jacobian, - _surface_flux_values) + 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!(elements, mesh, basis) return elements