Skip to content

Commit

Permalink
cleanup and remove unnecessary containers
Browse files Browse the repository at this point in the history
  • Loading branch information
@andrewwinters5000
andrewwinters5000 committed Mar 4, 2024
1 parent d28e814 commit f1cabf8
Show file tree
Hide file tree
Showing 3 changed files with 14 additions and 177 deletions.
7 changes: 2 additions & 5 deletions src/equations/numerical_fluxes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -415,7 +415,8 @@ flux vector splitting.
The [`SurfaceIntegralUpwind`](@ref) with a given `splitting` is equivalent to
the [`SurfaceIntegralStrongForm`](@ref) with `FluxUpwind(splitting)`
as numerical flux (up to floating point differences).
as numerical flux (up to floating point differences). Note, that
[`SurfaceIntegralUpwind`](@ref) is only available on [`TreeMesh`](@ref).
!!! warning "Experimental implementation (upwind SBP)"
This is an experimental feature and may change in future releases.
Expand All @@ -431,14 +432,10 @@ end
return fm + fp
end

# TODO: Upwind FD. Figure out a better strategy to compute this, especially
# if we need to pass two versions of the normal direction
@inline function (numflux::FluxUpwind)(u_ll, u_rr,
normal_direction::AbstractVector,
equations::AbstractEquations{2})
@unpack splitting = numflux
# Compute splitting in generic normal direction with specialized
# eigenvalues estimates calculated inside the `splitting` function
f_tilde_m = splitting(u_rr, Val{:minus}(), normal_direction, equations)
f_tilde_p = splitting(u_ll, Val{:plus}() , normal_direction, equations)
return f_tilde_m + f_tilde_p
Expand Down
101 changes: 2 additions & 99 deletions src/solvers/dgsem_unstructured/containers_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,103 +5,6 @@
@muladd begin
#! format: noindent


######################################################################################################
# TODO: FD; where should this live? Want to put it into `solvers/fdsbp_unstructured/containers_2d.jl`
# but then dispatching is not possible to reuse functionality for interfaces/boundaries here
# unless the FDSBP solver files are included before the DG files, which seems strange.
# Container data structure (structure-of-arrays style) for FDSBP upwind solver elements on curved unstructured mesh
struct UpwindElementContainer2D{RealT<:Real, uEltype<:Real}
node_coordinates ::Array{RealT, 4} # [ndims, nnodes, nnodes, nelement]
jacobian_matrix ::Array{RealT, 5} # [ndims, ndims, nnodes, nnodes, nelement]
inverse_jacobian ::Array{RealT, 3} # [nnodes, nnodes, nelement]
contravariant_vectors ::Array{RealT, 5} # [ndims, ndims, nnodes, nnodes, nelement]
contravariant_vectors_plus ::Array{RealT, 5} # [ndims, ndims, nnodes, nnodes, nelement]
contravariant_vectors_minus::Array{RealT, 5} # [ndims, ndims, nnodes, nnodes, nelement]
normal_directions ::Array{RealT, 4} # [ndims, nnodes, local sides, nelement]
normal_directions_plus ::Array{RealT, 4} # [ndims, nnodes, local sides, nelement]
normal_directions_minus ::Array{RealT, 4} # [ndims, nnodes, local sides, nelement]
rotations ::Vector{Int} # [nelement]
surface_flux_values ::Array{uEltype, 4} # [variables, nnodes, local sides, elements]
end


# construct an empty curved element container for the `UpwindOperators` type solver
# to be filled later with geometries in the unstructured mesh constructor
# OBS! Extended version of the `UnstructuredElementContainer2D` with additional arrays to hold
# the contravariant vectors created with the biased derivative operators.
function UpwindElementContainer2D{RealT, uEltype}(capacity::Integer, n_variables, n_nodes) where {RealT<:Real, uEltype<:Real}
nan_RealT = convert(RealT, NaN)
nan_uEltype = convert(uEltype, NaN)

node_coordinates = fill(nan_RealT, (2, n_nodes, n_nodes, capacity))
jacobian_matrix = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
inverse_jacobian = fill(nan_RealT, (n_nodes, n_nodes, capacity))
contravariant_vectors = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
contravariant_vectors_plus = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
contravariant_vectors_minus = fill(nan_RealT, (2, 2, n_nodes, n_nodes, capacity))
normal_directions = fill(nan_RealT, (2, n_nodes, 4, capacity))
normal_directions_plus = fill(nan_RealT, (2, n_nodes, 4, capacity))
normal_directions_minus = fill(nan_RealT, (2, n_nodes, 4, capacity))
rotations = fill(typemin(Int), capacity) # Fill with "nonsense" rotation values
surface_flux_values = fill(nan_uEltype, (n_variables, n_nodes, 4, capacity))

return UpwindElementContainer2D{RealT, uEltype}(node_coordinates,
jacobian_matrix,
inverse_jacobian,
contravariant_vectors,
contravariant_vectors_plus,
contravariant_vectors_minus,
normal_directions,
normal_directions_plus,
normal_directions_minus,
rotations,
surface_flux_values)
end


@inline nelements(elements::UpwindElementContainer2D) = size(elements.surface_flux_values, 4)
@inline eachelement(elements::UpwindElementContainer2D) = Base.OneTo(nelements(elements))

@inline nvariables(elements::UpwindElementContainer2D) = size(elements.surface_flux_values, 1)
@inline nnodes(elements::UpwindElementContainer2D) = size(elements.surface_flux_values, 2)

Base.real(elements::UpwindElementContainer2D) = eltype(elements.node_coordinates)
Base.eltype(elements::UpwindElementContainer2D) = eltype(elements.surface_flux_values)

function init_elements(mesh::UnstructuredMesh2D, equations,
basis::SummationByPartsOperators.UpwindOperators,
RealT, uEltype)
elements = UpwindElementContainer2D{RealT, uEltype}(
mesh.n_elements, nvariables(equations), nnodes(basis))
init_elements!(elements, mesh, basis)
return elements
end

function init_elements!(elements::UpwindElementContainer2D, mesh, basis)
four_corners = zeros(eltype(mesh.corners), 4, 2)

# loop through elements and call the correct constructor based on whether the element is curved
for element in eachelement(elements)
if mesh.element_is_curved[element]
init_element!(elements, element, basis, view(mesh.surface_curves, :, element))
else # straight sided element
for i in 1:4, j in 1:2
# pull the (x,y) values of these corners out of the global corners array
four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
end
init_element!(elements, element, basis, four_corners)
end
end

# Use the mesh element information to determine the rotations, if any,
# of the local coordinate system in each element
# TODO: remove me. unnecessary (and slightly broken)
calc_element_rotations!(elements, mesh)
end

######################################################################################################

# Container data structure (structure-of-arrays style) for DG elements on curved unstructured mesh
struct UnstructuredElementContainer2D{RealT <: Real, uEltype <: Real}
node_coordinates::Array{RealT, 4} # [ndims, nnodes, nnodes, nelement]
Expand Down Expand Up @@ -244,7 +147,7 @@ end
@inline nnodes(interfaces::UnstructuredInterfaceContainer2D) = size(interfaces.u, 3)

function init_interfaces(mesh::UnstructuredMesh2D,
elements::Union{UnstructuredElementContainer2D,UpwindElementContainer2D})
elements::UnstructuredElementContainer2D)
interfaces = UnstructuredInterfaceContainer2D{eltype(elements)}(mesh.n_interfaces,
nvariables(elements),
nnodes(elements))
Expand Down Expand Up @@ -385,7 +288,7 @@ end
end

function init_boundaries(mesh::UnstructuredMesh2D,
elements::Union{UnstructuredElementContainer2D,UpwindElementContainer2D})
elements::UnstructuredElementContainer2D)
boundaries = UnstructuredBoundaryContainer2D{real(elements), eltype(elements)}(mesh.n_boundaries,
nvariables(elements),
nnodes(elements))
Expand Down
83 changes: 10 additions & 73 deletions src/solvers/fdsbp_unstructured/containers_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,21 @@
#! format: noindent

# initialize all the values in the container of a general FD block (either straight sided or curved)
# OBS! Requires the SBP derivative matrix in order to compute metric terms that are free-stream preserving
# OBS! Requires the SBP derivative matrix in order to compute metric terms. For a standard
# SBP operator the matrix is passed directly, whereas for an upwind SBP operator we must
# specifically pass the central differencing matrix.
function init_element!(elements, element, basis::AbstractDerivativeOperator,
corners_or_surface_curves)
calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis),
corners_or_surface_curves)

calc_metric_terms!(elements.jacobian_matrix, element, basis,
elements.node_coordinates)
if basis isa DerivativeOperator
calc_metric_terms!(elements.jacobian_matrix, element, basis,
elements.node_coordinates)
else # basis isa SummationByPartsOperators.UpwindOperators
calc_metric_terms!(elements.jacobian_matrix, element, basis.central,
elements.node_coordinates)
end

calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)

Expand All @@ -29,37 +36,8 @@ function init_element!(elements, element, basis::AbstractDerivativeOperator,
return elements
end

# initialize all the values in the container of a general FD block (either straight sided or curved)
# for a set of upwind finite difference operators
# OBS! (Maybe) Requires the biased derivative matrices in order to compute proper metric terms
# that are free-stream preserving. If this is not necessary, then we do not need this specialized container.
function init_element!(elements, element, basis::SummationByPartsOperators.UpwindOperators, corners_or_surface_curves)

calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), corners_or_surface_curves)

# Create the metric terms and contravariant vectors with the D⁺ operator
calc_metric_terms!(elements.jacobian_matrix, element, basis.plus, elements.node_coordinates)
calc_contravariant_vectors!(elements.contravariant_vectors_plus, element, elements.jacobian_matrix)
calc_normal_directions!(elements.normal_directions_plus, element, elements.jacobian_matrix)

# Create the metric terms and contravariant vectors with the D⁻ operator
calc_metric_terms!(elements.jacobian_matrix, element, basis.minus, elements.node_coordinates)
calc_contravariant_vectors!(elements.contravariant_vectors_minus, element, elements.jacobian_matrix)
calc_normal_directions!(elements.normal_directions_minus, element, elements.jacobian_matrix)

# Create the metric terms, Jacobian information, and normals for analysis
# and the SATs with the central D operator.
calc_metric_terms!(elements.jacobian_matrix, element, basis.central, elements.node_coordinates)
calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)
calc_contravariant_vectors!(elements.contravariant_vectors, element, elements.jacobian_matrix)
calc_normal_directions!(elements.normal_directions, element, elements.jacobian_matrix)

return elements
end

# construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
# applied to the node coordinates.
# TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices
function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator,
node_coordinates)

Expand All @@ -69,9 +47,6 @@ function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeO
# jacobian_matrix[2,1,:,:,:] <- Y_xi
# jacobian_matrix[2,2,:,:,:] <- Y_eta

# TODO: for now the upwind version needs to specify which matrix is used.
# requires more testing / debugging but for now use the central SBP matrix

# Compute the xi derivatives by applying D on the left
# This is basically the same as
# jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element]
Expand Down Expand Up @@ -152,42 +127,4 @@ function calc_normal_directions!(normal_directions, element, jacobian_matrix)

return normal_directions
end

# Compute the rotation, if any, for the local coordinate system on each `element`
# with respect to the standard x-axis. This is necessary because if the local
# element axis is rotated it influences the computation of the +/- directions
# for the flux vector splitting of the upwind scheme.
# Local principle x-axis is computed from the four corners of a particular `element`,
# see page 35 of the "Verdict Library" for details.
# From the local axis the `atan` function is used to determine the value of the local rotation.
#
# - C. J. Stimpson, C. D. Ernst, P. Knupp, P. P. Pébay, and D. Thompson (2007)
# The Verdict Geometric Quality Library
# Technical Report. Sandia National Laboraties
# [SAND2007-1751](https://coreform.com/papers/verdict_quality_library.pdf)
# - `atan(y, x)` function (https://docs.julialang.org/en/v1/base/math/#Base.atan-Tuple%7BNumber%7D)
# TODO: remove me. unnecessary
function calc_element_rotations!(elements, mesh::UnstructuredMesh2D)

for element in 1:size(elements.node_coordinates, 4)
# Pull the corners for the current element
corner_points = mesh.corners[:, mesh.element_node_ids[:, element]]

# Compute the principle x-axis of a right-handed quadrilateral element
local_x_axis = ( (corner_points[:, 2] - corner_points[:, 1])
+ (corner_points[:, 3] - corner_points[:, 4]) )

# Two argument `atan` function retuns the angle in radians in [−pi, pi]
# between the positive x-axis and the point (x, y)
local_angle = atan(local_x_axis[2], local_x_axis[1])

# Possible reference rotations of the local axis for a given quadrilateral element in the mesh
# 0°, 90°, 180°, -90° (or 270°), -180°
reference_angles = [0.0, 0.5*pi, pi, -0.5*pi, -pi]

elements.rotations[element] = argmin( abs.(reference_angles .- local_angle) ) - 1
end

return nothing
end
end # @muladd

0 comments on commit f1cabf8

Please sign in to comment.