diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index a31e1abaef0..f0a91dfcbc3 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -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. @@ -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 diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl index 97599f7f0a0..f51dd09801b 100644 --- a/src/solvers/dgsem_unstructured/containers_2d.jl +++ b/src/solvers/dgsem_unstructured/containers_2d.jl @@ -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] @@ -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)) @@ -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)) diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl index 73ac0e33974..76facd59394 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -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) @@ -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) @@ -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] @@ -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