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..28deb4964e1 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl @@ -0,0 +1,148 @@ + +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) + +# Initial condition for a Gaussian density profile with constant pressure +# and the velocity of a rotating solid body +function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEquations3D) + # Gaussian density + rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) + # Constant pressure + p = 1.0 + + # Spherical coordinates for the point x + if sign(x[2]) == 0.0 + signy = 1.0 + else + signy = sign(x[2]) + end + # Co-latitude + colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) + # Latitude (auxiliary variable) + lat = -colat + 0.5 * pi + # Longitude + r_xy = sqrt(x[1]^2 + x[2]^2) + if r_xy == 0.0 + phi = pi / 2 + else + phi = signy * acos(x[1] / r_xy) + end + + # Compute the velocity of a rotating solid body + # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system) + v0 = 1.0 # Velocity at the "equator" + alpha = 0.0 #pi / 4 + v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) + v_lat = -v0 * sin(phi) * sin(alpha) + + # Transform to Cartesian coordinate system + v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long + v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long + v3 = sin(colat) * v_lat + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +# Source term function to transform the Euler equations into the linear advection equations with variable advection velocity +function source_terms_convert_to_linear_advection(u, du, x, t, + equations::CompressibleEulerEquations3D) + v1 = u[2] / u[1] + v2 = u[3] / u[1] + v3 = u[4] / u[1] + + s2 = du[1] * v1 - du[2] + s3 = du[1] * v2 - du[3] + s4 = du[1] * v3 - du[4] + s5 = 0.5 * (s2 * v1 + s3 * v2 + s4 * v3) - du[5] + + return SVector(0.0, s2, s3, s4, s5) +end + +# Custom RHS that applies a source term that depends on du (used to convert the 3D Euler equations into the 3D linear advection +# equations with position-dependent advection speed) +function rhs_semi_custom!(du_ode, u_ode, semi, t) + # Compute standard Trixi RHS + Trixi.rhs!(du_ode, u_ode, semi, t) + + # Now apply the custom source term + Trixi.@trixi_timeit Trixi.timer() "custom source term" begin + @unpack solver, equations, cache = semi + @unpack node_coordinates = cache.elements + + # Wrap the solution and RHS + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + + Trixi.@threaded for element in eachelement(solver, cache) + for j in eachnode(solver), i in eachnode(solver) + u_local = Trixi.get_node_vars(u, equations, solver, i, j, element) + du_local = Trixi.get_node_vars(du, equations, solver, i, j, element) + x_local = Trixi.get_node_coords(node_coordinates, equations, solver, + i, j, element) + source = source_terms_convert_to_linear_advection(u_local, du_local, + x_local, t, equations) + Trixi.add_to_node_vars!(du, source, equations, solver, i, j, element) + end + end + end +end + +initial_condition = initial_condition_advection_sphere + +mesh = 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) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +ode = semidiscretize(semi, tspan) + +# Create custom discretization that runs with the custom RHS +ode_semi_custom = ODEProblem(rhs_semi_custom!, + ode.u0, + ode.tspan, + semi) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 10, + save_analysis = true, + extra_analysis_integrals = (Trixi.density,)) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.7) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode_semi_custom, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index 1a509ed92d1..7cb4968b0c6 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -50,7 +50,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/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 2046220aeca..101f981f19e 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -772,6 +772,57 @@ 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 a 2D `P4estMesh` with +`6 * trees_per_face_dimension^2` trees. + +The mesh will have no boundaries. + +# Arguments +- `trees_per_face_dimension::Integer`: the number of trees in the two local dimensions of + each face. +- `radius::Integer`: the radius of the sphere. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. +- `RealT::Type`: the type that should be used for coordinates. +- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +- `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file. +- `p4est_partition_allow_for_coarsening::Bool`: Must be `true` when using AMR to make mesh adaptivity + independent of domain partitioning. Should be `false` for static meshes + to permit more fine-grained partitioning. +""" +function P4estMeshCubedSphere2D(trees_per_face_dimension, radius; + polydeg, RealT = Float64, + initial_refinement_level = 0, unsaved_changes = true, + p4est_partition_allow_for_coarsening = true) + connectivity = connectivity_cubed_sphere_2D(trees_per_face_dimension) + + n_trees = 6 * trees_per_face_dimension^2 + + basis = LobattoLegendreBasis(RealT, polydeg) + nodes = basis.nodes + + tree_node_coordinates = Array{RealT, 4}(undef, 3, + ntuple(_ -> length(nodes), 2)..., + n_trees) + calc_tree_node_coordinates_cubed_sphere_2D!(tree_node_coordinates, nodes, + trees_per_face_dimension, radius) + + p4est = 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 @@ -989,6 +1040,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 @@ -1239,6 +1291,226 @@ 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 @@ -1438,6 +1710,45 @@ 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 7c82a132a0b..007b00dd617 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -41,7 +41,7 @@ struct SemidiscretizationHyperbolic{Mesh, Equations, InitialCondition, 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 f9830d0011c..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} <: 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::Array{RealT, NDIMSP3} # [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 @@ -87,19 +89,29 @@ 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 * NDIMS * nnodes(basis)^NDIMS * + nelements) jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix), - (NDIMS, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., + (ndims_spa, NDIMS, + ntuple(_ -> nnodes(basis), NDIMS)..., nelements)) - _contravariant_vectors = similar(_jacobian_matrix) - contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors), - size(jacobian_matrix)) + _contravariant_vectors = Vector{RealT}(undef, + ndims_spa^2 * 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), @@ -115,12 +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}(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 diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 236d7d24c06..ad6156ac07e 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -13,12 +13,32 @@ 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 + # 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 (cross-product form) + 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) + inverse_jacobian[i, j, element] = 1 / + norm(cross(jacobian_matrix[:, 1, i, j, + element], + jacobian_matrix[:, 2, i, j, + element])) + end + 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_contravariant_vectors!(contravariant_vectors, element, jacobian_matrix) - calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix) + calc_inverse_jacobian!(inverse_jacobian, element, jacobian_matrix) + end end return nothing @@ -43,7 +63,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(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) @@ -84,6 +105,63 @@ 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 +# 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 + }, + element, + jacobian_matrix, node_coordinates, + basis::LobattoLegendreBasis) + @unpack derivative_matrix = basis + + for j in eachnode(basis), i in eachnode(basis) + for n in 1:3 + # (n, m, l) cyclic + m = (n % 3) + 1 + l = ((n + 1) % 3) + 1 + + contravariant_vectors[n, 1, i, j, element] = (jacobian_matrix[m, 2, i, j, + element] * + node_coordinates[l, i, j, + element] + - + jacobian_matrix[l, 2, i, j, + element] * + node_coordinates[m, i, j, + element]) + + contravariant_vectors[n, 2, i, j, element] = (jacobian_matrix[l, 1, i, j, + element] * + node_coordinates[m, i, j, + element] + - + jacobian_matrix[m, 1, i, j, + element] * + node_coordinates[l, i, j, + element]) + + contravariant_vectors[n, 3, i, j, element] = (jacobian_matrix[m, 1, i, j, + element] * + jacobian_matrix[l, 2, i, j, + element] + - + jacobian_matrix[m, 2, i, j, + element] * + jacobian_matrix[l, 1, i, j, + element]) + end + end + + return contravariant_vectors +end + # 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.jl b/src/solvers/dgsem_structured/dg.jl index 5617ae90e3f..94cf8bcd2a1 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -23,9 +23,10 @@ function create_cache(mesh::Union{StructuredMesh, StructuredMeshView}, 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(ndims(contravariant_vectors) - 3))) + static2val(static_size(contravariant_vectors, StaticInt(1))))) end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 51b91eac6be..6ce86844306 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -61,7 +61,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 mesh::Union{StructuredMesh{2}, StructuredMeshView{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. @@ -73,7 +74,6 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 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) @@ -98,6 +98,55 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 return nothing 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 + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + flux1 = flux(u_node, 1, equations) + flux2 = flux(u_node, 2, equations) + flux3 = flux(u_node, 3, equations) + + # Compute the contravariant flux by taking the scalar product of the + # first contravariant vector Ja^1 and the flux vector + Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 + for ii in eachnode(dg) + 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 + + 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]) + end + end + + return nothing +end + @inline function flux_differencing_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index ef5ac2c9de3..dc4660a7ded 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -434,6 +434,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_euler_NACA6412airfoil_mach2.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA6412airfoil_mach2.jl"), l2=[