Skip to content

Commit

Permalink
Removed hard-coded computation of rhs for the transformation of Euler…
Browse files Browse the repository at this point in the history 
… into advection with variable coefficients, and improved formatting
  • Loading branch information
@amrueda
amrueda committed Nov 21, 2023
1 parent 1eb17ac commit e180aa4
Show file tree
Hide file tree
Showing 5 changed files with 96 additions and 55 deletions.
18 changes: 16 additions & 2 deletions examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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)
Expand Down
43 changes: 24 additions & 19 deletions src/meshes/p4est_mesh.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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
Expand Down
21 changes: 14 additions & 7 deletions src/solvers/dgsem_p4est/containers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down
23 changes: 12 additions & 11 deletions src/solvers/dgsem_structured/dg_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -120,20 +120,21 @@ 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)
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])
du[v, i, j, element] += (Ja31 * flux1[v] + Ja32 * flux2[v] +
Ja33 * flux3[v])
end
end
end

return nothing
end

Expand Down
46 changes: 30 additions & 16 deletions src/solvers/dgsem_tree/dg_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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

0 comments on commit e180aa4

Please sign in to comment.