From a6ba38e63c540fa2342df89412de511a42782c86 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 6 Feb 2024 22:10:51 +0100 Subject: [PATCH 01/17] baseline implementation of the curvilinear USBP for testing --- src/equations/compressible_euler_2d.jl | 74 ++++++- src/equations/numerical_fluxes.jl | 13 ++ src/solvers/fdsbp_tree/fdsbp_2d.jl | 2 +- .../fdsbp_unstructured/containers_2d.jl | 11 +- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 203 +++++++++++++++++- 5 files changed, 294 insertions(+), 9 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index f5a632723cf..89d61ecdf0a 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -706,10 +706,10 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() With Application to Finite Difference Methods [NASA Technical Memorandum](https://ntrs.nasa.gov/api/citations/19790020779/downloads/19790020779.pdf) """ -@inline function splitting_steger_warming(u, orientation::Integer, +@inline function splitting_steger_warming(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) - fm = splitting_steger_warming(u, Val{:minus}(), orientation, equations) - fp = splitting_steger_warming(u, Val{:plus}(), orientation, equations) + fm = splitting_steger_warming(u, Val{:minus}(), orientation_or_normal_direction, equations) + fp = splitting_steger_warming(u, Val{:plus}(), orientation_or_normal_direction, equations) return fm, fp end @@ -809,6 +809,74 @@ end return SVector(f1m, f2m, f3m, f4m) end +@inline function splitting_steger_warming(u, ::Val{:plus}, + normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + a = sqrt(equations.gamma * p / rho) + + s_hat = norm(normal_direction) # sqrt(normal_direction[1]^2 + normal_direction[2]^2) + normal_vector = normal_direction / s_hat + v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + + lambda1 = normal_direction[1] * v1 + normal_direction[2] * v2 + lambda2 = lambda1 + a * s_hat + lambda3 = lambda1 - a * s_hat + + lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) + lambda2_p = positive_part(lambda2) + lambda3_p = positive_part(lambda3) + + alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p + + rho_2gamma = 0.5 * rho / equations.gamma + f1p = rho_2gamma * alpha_p + f2p = rho_2gamma * (alpha_p * v1 + a * normal_vector[1] * (lambda2_p - lambda3_p)) + f3p = rho_2gamma * (alpha_p * v2 + a * normal_vector[2] * (lambda2_p - lambda3_p)) + f4p = rho_2gamma * + (alpha_p * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_p - lambda3_p) + + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) + + return SVector(f1p, f2p, f3p, f4p) +end + +@inline function splitting_steger_warming(u, ::Val{:minus}, + normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + a = sqrt(equations.gamma * p / rho) + + s_hat = norm(normal_direction) # sqrt(normal_direction[1]^2 + normal_direction[2]^2) + normal_vector = normal_direction / s_hat + v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + + lambda1 = normal_direction[1] * v1 + normal_direction[2] * v2 + lambda2 = lambda1 + a * s_hat + lambda3 = lambda1 - a * s_hat + + lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) + lambda2_m = negative_part(lambda2) + lambda3_m = negative_part(lambda3) + + alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m + + rho_2gamma = 0.5 * rho / equations.gamma + f1m = rho_2gamma * alpha_m + f2m = rho_2gamma * (alpha_m * v1 + a * normal_vector[1] * (lambda2_m - lambda3_m)) + f3m = rho_2gamma * (alpha_m * v2 + a * normal_vector[2] * (lambda2_m - lambda3_m)) + f4m = rho_2gamma * + (alpha_m * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_m - lambda3_m) + + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) + + return SVector(f1m, f2m, f3m, f4m) +end + """ FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 44d523b6e89..32762b37d14 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -444,5 +444,18 @@ end return fm + fp end +# TODO: Upwind FD. Figure out a better strategy to compute this +# (From an old PR and seemed to not work correctly!) +@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 +end + Base.show(io::IO, f::FluxUpwind) = print(io, "FluxUpwind(", f.splitting, ")") end # @muladd diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl index 09d18cecd75..36afbbc022f 100644 --- a/src/solvers/fdsbp_tree/fdsbp_2d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_2d.jl @@ -19,7 +19,7 @@ function create_cache(mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, equations, return (; f_threaded) end -function create_cache(mesh::TreeMesh{2}, equations, +function create_cache(mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, equations, volume_integral::VolumeIntegralUpwind, dg, uEltype) u_node = SVector{nvariables(equations), uEltype}(ntuple(_ -> zero(uEltype), Val{nvariables(equations)}())) diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl index 3857c2d8a20..d83d4d37557 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -41,29 +41,32 @@ 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] # but uses only matrix-vector products instead of a matrix-matrix product. for j in eachnode(D_SBP) - mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP, + mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP.central, view(node_coordinates, 1, :, j, element)) end # jacobian_matrix[2, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[2, :, :, element] for j in eachnode(D_SBP) - mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP, + mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP.central, view(node_coordinates, 2, :, j, element)) end # Compute the eta derivatives by applying transpose of D on the right # jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(D_SBP)' for i in eachnode(D_SBP) - mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP, + mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP.central, view(node_coordinates, 1, i, :, element)) end # jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(D_SBP)' for i in eachnode(D_SBP) - mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP, + mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP.central, view(node_coordinates, 2, i, :, element)) end diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index b459f4c42cc..4ee2c873641 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -86,6 +86,148 @@ end return nothing end +# 2D volume integral contributions for `VolumeIntegralUpwind`. +# Note that the plus / minus notation of the operators does not refer to the +# upwind / downwind directions of the fluxes. +# Instead, the plus / minus refers to the direction of the biasing within +# the finite difference stencils. Thus, the D^- operator acts on the positive +# part of the flux splitting f^+ and the D^+ operator acts on the negative part +# of the flux splitting f^-. +function calc_volume_integral!(du, u, + mesh::UnstructuredMesh2D, + nonconservative_terms::False, equations, + volume_integral::VolumeIntegralUpwind, + dg::FDSBP, cache) + # Assume that + # dg.basis isa SummationByPartsOperators.UpwindOperators + D_minus = dg.basis.minus # Upwind SBP D^- derivative operator + D_plus = dg.basis.plus # Upwind SBP D^+ derivative operator + # D_minus = dg.basis.central # Central SBP D derivative operator + # D_plus = dg.basis.central # Central SBP D derivative operator + @unpack f_minus_plus_threaded, f_minus_threaded, f_plus_threaded = cache + @unpack splitting = volume_integral + @unpack contravariant_vectors = cache.elements + + # SBP operators from SummationByPartsOperators.jl implement the basic interface + # of matrix-vector multiplication. Thus, we pass an "array of structures", + # packing all variables per node in an `SVector`. + if nvariables(equations) == 1 + # `reinterpret(reshape, ...)` removes the leading dimension only if more + # than one variable is used. + u_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, + du), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + else + u_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u) + du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, + du) + end + + # Use the tensor product structure to compute the discrete derivatives of + # the fluxes line-by-line and add them to `du` for each element. + @threaded for element in eachelement(dg, cache) + # f_minus_plus_element wraps the storage provided by f_minus_element and + # f_plus_element such that we can use a single plain broadcasting below. + # f_minus_element and f_plus_element are updated in broadcasting calls + # of the form `@. f_minus_plus_element = ...`. + f_minus_plus_element = f_minus_plus_threaded[Threads.threadid()] + f_minus_element = f_minus_threaded[Threads.threadid()] + f_plus_element = f_plus_threaded[Threads.threadid()] + u_element = view(u_vectors, :, :, element) + + # x direction + # @. f_minus_plus_element = splitting(u_element, 1, equations) + for j in eachnode(dg) + for i in eachnode(dg) + Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja1, equations) + end + mul!(view(du_vectors, :, j, element), D_minus, view(f_plus_element, :, j), + one(eltype(du)), one(eltype(du))) + mul!(view(du_vectors, :, j, element), D_plus, view(f_minus_element, :, j), + one(eltype(du)), one(eltype(du))) + end + + # y direction + # @. f_minus_plus_element = splitting(u_element, 2, equations) + for i in eachnode(dg) + for j in eachnode(dg) + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja2, equations) + end + mul!(view(du_vectors, i, :, element), D_minus, view(f_plus_element, i, :), + one(eltype(du)), one(eltype(du))) + mul!(view(du_vectors, i, :, element), D_plus, view(f_minus_element, i, :), + one(eltype(du)), one(eltype(du))) + end + end + + return nothing +end + +# TODO: Write a better comment here. Specialized computation of the interface +# values using flux vector splitting in an arbitrary direction +function calc_interface_flux!(surface_flux_values, + mesh::UnstructuredMesh2D, + nonconservative_terms::False, equations, + surface_integral::SurfaceIntegralUpwind, + dg::FDSBP, cache) + @unpack splitting = surface_integral + @unpack u, start_index, index_increment, element_ids, element_side_ids = cache.interfaces + @unpack normal_directions = cache.elements + + @threaded for interface in eachinterface(dg, cache) + # Get neighboring elements + primary_element = element_ids[1, interface] + secondary_element = element_ids[2, interface] + + # Get the local side id on which to compute the flux + primary_side = element_side_ids[1, interface] + secondary_side = element_side_ids[2, interface] + + # initial index for the coordinate system on the secondary element + secondary_index = start_index[interface] + + # loop through the primary element coordinate system and compute the interface coupling + for primary_index in eachnode(dg) + # Pull the primary and secondary states from the boundary u values + u_ll = get_one_sided_surface_node_vars(u, equations, dg, 1, primary_index, + interface) + u_rr = get_one_sided_surface_node_vars(u, equations, dg, 2, secondary_index, + interface) + + # pull the outward pointing (normal) directional vector + # Note! this assumes a conforming approximation, more must be done in terms of the normals + # for hanging nodes and other non-conforming approximation spaces + outward_direction = get_surface_normal(normal_directions, primary_index, + primary_side, + primary_element) + + # Compute the upwind coupling terms where right-traveling + # information comes from the left and left-traveling information + # comes from the right + flux_minus_rr = splitting(u_rr, Val{:minus}(), outward_direction, equations) + flux_plus_ll = splitting(u_ll, Val{:plus}(), outward_direction, equations) + + # Copy upwind fluxes back to primary/secondary element storage + # Note the sign change for the normal flux in the secondary element! + # TODO: do we need this sign flip on the neighbour? + # OBS! if we switch the sign flip here then the `SurfaceIntegralUpwind` also needs adjusted + for v in eachvariable(equations) + surface_flux_values[v, primary_index, primary_side, primary_element] = flux_minus_rr[v] + surface_flux_values[v, secondary_index, secondary_side, secondary_element] = -flux_plus_ll[v] + end + + # increment the index of the coordinate system in the secondary element + secondary_index += index_increment[interface] + end + end + + return nothing +end + # Note! The local side numbering for the unstructured quadrilateral element implementation differs # from the structured TreeMesh or StructuredMesh local side numbering: # @@ -104,7 +246,7 @@ end # surface contributions are added. function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, equations, surface_integral::SurfaceIntegralStrongForm, - dg::DG, cache) + dg::FDSBP, cache) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack normal_directions, surface_flux_values = cache.elements @@ -156,6 +298,65 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, return nothing end +# Implementation of fully upwind SATs. The surface flux values are pre-computed +# in the specialized `calc_interface_flux` routine. These SATs are still of +# a strong form penalty type, except that the interior flux at a particular +# side of the element are computed in the upwind direction. +function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, + equations, surface_integral::SurfaceIntegralUpwind, + dg::FDSBP, cache) + inv_weight_left = inv(left_boundary_weight(dg.basis)) + inv_weight_right = inv(right_boundary_weight(dg.basis)) + @unpack normal_directions, surface_flux_values = cache.elements + @unpack splitting = surface_integral + + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + # surface at -x + u_node = get_node_vars(u, equations, dg, 1, l, element) + # compute internal flux in normal direction on side 4 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, + element) + f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) + multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), + equations, dg, 1, l, element) + + # surface at +x + u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) + # compute internal flux in normal direction on side 2 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, + element) + f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) + multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), + equations, dg, nnodes(dg), l, element) + + # surface at -y + u_node = get_node_vars(u, equations, dg, l, 1, element) + # compute internal flux in normal direction on side 1 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, + element) + f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) + multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), + equations, dg, l, 1, element) + + # surface at +y + u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) + # compute internal flux in normal direction on side 3 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, + element) + f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) + multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), + equations, dg, l, nnodes(dg), element) + end + end + + return nothing +end + # AnalysisCallback function integrate_via_indices(func::Func, u, mesh::UnstructuredMesh2D, equations, From 3aa5ea4c2d747dcef7fc59a351d1c6792e86b642 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 1 Mar 2024 07:40:25 +0100 Subject: [PATCH 02/17] working version of curvilinear upwind solver. Needs significant cleanup of debugging statements and different variants of the rotated flux vector splittings --- src/equations/compressible_euler_2d.jl | 328 ++++++++++++-- src/equations/numerical_fluxes.jl | 4 +- .../dgsem_unstructured/containers_2d.jl | 101 ++++- src/solvers/dgsem_unstructured/dg_2d.jl | 2 +- .../fdsbp_unstructured/containers_2d.jl | 74 ++- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 427 ++++++++++++------ 6 files changed, 752 insertions(+), 184 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 89d61ecdf0a..8935cb12eee 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -818,29 +818,63 @@ end p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) - s_hat = norm(normal_direction) # sqrt(normal_direction[1]^2 + normal_direction[2]^2) - normal_vector = normal_direction / s_hat - v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # ## equivalent to rotate u, do splitting and rotate back + # s_hat = norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - lambda1 = normal_direction[1] * v1 + normal_direction[2] * v2 - lambda2 = lambda1 + a * s_hat - lambda3 = lambda1 - a * s_hat + # lambda1 = v_n + # lambda2 = lambda1 + a * s_hat + # lambda3 = lambda1 - a * s_hat + + # lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) + # lambda2_p = positive_part(lambda2) + # lambda3_p = positive_part(lambda3) + + # alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p + + # rho_2gamma = 0.5 * rho / equations.gamma + # f1p = rho_2gamma * alpha_p + # f2p = rho_2gamma * (alpha_p * v1 + a * normal_vector[1] * (lambda2_p - lambda3_p)) + # f3p = rho_2gamma * (alpha_p * v2 + a * normal_vector[2] * (lambda2_p - lambda3_p)) + # f4p = rho_2gamma * + # (alpha_p * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_p - lambda3_p) + # + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) + + # return SVector(f1p, f2p, f3p, f4p) + + ## alternative scaling from the paper of Drikakis and Tsangris + # s_hat = norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + lambda1 = v_n + a + lambda2 = v_n - a lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) lambda2_p = positive_part(lambda2) - lambda3_p = positive_part(lambda3) - alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p + ## original form from the paper + # H = (rho_e + p) / rho + # f1p = 0.5 * rho * (lambda1_p + lambda2_p) + # f2p = (0.5 * rho * (v1 + normal_direction[1] * a / equations.gamma) * lambda1_p + # + 0.5 * rho * (v1 - normal_direction[1] * a / equations.gamma) * lambda2_p) + # f3p = (0.5 * rho * (v2 + normal_direction[2] * a / equations.gamma) * lambda1_p + # + 0.5 * rho * (v2 - normal_direction[2] * a / equations.gamma) * lambda2_p) + # f4p = f1p * H + + # slightly rewritten form that is cleaner + rhoa_2gamma = 0.5 * rho * a / equations.gamma + H = (rho_e + p) / rho - rho_2gamma = 0.5 * rho / equations.gamma - f1p = rho_2gamma * alpha_p - f2p = rho_2gamma * (alpha_p * v1 + a * normal_vector[1] * (lambda2_p - lambda3_p)) - f3p = rho_2gamma * (alpha_p * v2 + a * normal_vector[2] * (lambda2_p - lambda3_p)) - f4p = rho_2gamma * - (alpha_p * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_p - lambda3_p) - + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) + f1p = 0.5 * rho * (lambda1_p + lambda2_p) + f2p = f1p * v1 + rhoa_2gamma * normal_direction[1] * (lambda1_p - lambda2_p) + f3p = f1p * v2 + rhoa_2gamma * normal_direction[2] * (lambda1_p - lambda2_p) + f4p = f1p * H - return SVector(f1p, f2p, f3p, f4p) + return SVector(f1p, f2p, f3p, f4p)# * s_hat end @inline function splitting_steger_warming(u, ::Val{:minus}, @@ -852,29 +886,63 @@ end p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) - s_hat = norm(normal_direction) # sqrt(normal_direction[1]^2 + normal_direction[2]^2) - normal_vector = normal_direction / s_hat - v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # # equivalent to rotate u, do splitting and rotate back + # s_hat = norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + # lambda1 = v_n + # lambda2 = lambda1 + a * s_hat + # lambda3 = lambda1 - a * s_hat - lambda1 = normal_direction[1] * v1 + normal_direction[2] * v2 - lambda2 = lambda1 + a * s_hat - lambda3 = lambda1 - a * s_hat + # lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) + # lambda2_m = negative_part(lambda2) + # lambda3_m = negative_part(lambda3) + + # alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m + + # rho_2gamma = 0.5 * rho / equations.gamma + # f1m = rho_2gamma * alpha_m + # f2m = rho_2gamma * (alpha_m * v1 + a * normal_vector[1] * (lambda2_m - lambda3_m)) + # f3m = rho_2gamma * (alpha_m * v2 + a * normal_vector[2] * (lambda2_m - lambda3_m)) + # f4m = rho_2gamma * + # (alpha_m * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_m - lambda3_m) + # + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) + + # return SVector(f1m, f2m, f3m, f4m) + + ## alternative scaling from the paper of Drikakis and Tsangris + # s_hat = norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + lambda1 = v_n + a + lambda2 = v_n - a lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) lambda2_m = negative_part(lambda2) - lambda3_m = negative_part(lambda3) - alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m + ## original version from the paper + # H = (rho_e + p) / rho + # f1m = 0.5 * rho * (lambda1_m + lambda2_m) + # f2m = (0.5 * rho * (v1 + normal_direction[1] * a / equations.gamma) * lambda1_m + # + 0.5 * rho * (v1 - normal_direction[1] * a / equations.gamma) * lambda2_m) + # f3m = (0.5 * rho * (v2 + normal_direction[2] * a / equations.gamma) * lambda1_m + # + 0.5 * rho * (v2 - normal_direction[2] * a / equations.gamma) * lambda2_m) + # f4m = f1m * H + + # slightly rewritten form that is cleaner + rhoa_2gamma = 0.5 * rho * a / equations.gamma + H = (rho_e + p) / rho - rho_2gamma = 0.5 * rho / equations.gamma - f1m = rho_2gamma * alpha_m - f2m = rho_2gamma * (alpha_m * v1 + a * normal_vector[1] * (lambda2_m - lambda3_m)) - f3m = rho_2gamma * (alpha_m * v2 + a * normal_vector[2] * (lambda2_m - lambda3_m)) - f4m = rho_2gamma * - (alpha_m * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_m - lambda3_m) - + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) + f1m = 0.5 * rho * (lambda1_m + lambda2_m) + f2m = f1m * v1 + rhoa_2gamma * normal_direction[1] * (lambda1_m - lambda2_m) + f3m = f1m * v2 + rhoa_2gamma * normal_direction[2] * (lambda1_m - lambda2_m) + f4m = f1m * H - return SVector(f1m, f2m, f3m, f4m) + return SVector(f1m, f2m, f3m, f4m)# * s_hat end """ @@ -1003,10 +1071,10 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() High-Order Polynomial Expansions (HOPE) for Flux-Vector Splitting [NASA Technical Memorandum](https://ntrs.nasa.gov/citations/19910016425) """ -@inline function splitting_vanleer_haenel(u, orientation::Integer, +@inline function splitting_vanleer_haenel(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) - fm = splitting_vanleer_haenel(u, Val{:minus}(), orientation, equations) - fp = splitting_vanleer_haenel(u, Val{:plus}(), orientation, equations) + fm = splitting_vanleer_haenel(u, Val{:minus}(), orientation_or_normal_direction, equations) + fp = splitting_vanleer_haenel(u, Val{:plus}(), orientation_or_normal_direction, equations) return fm, fp end @@ -1070,6 +1138,150 @@ end return SVector(f1m, f2m, f3m, f4m) end +@inline function splitting_vanleer_haenel(u, ::Val{:plus}, + normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + + a = sqrt(equations.gamma * p / rho) + H = (rho_e + p) / rho + + ## equivalent to the rotate u, calc flux, and then backrotate + # s_hat = norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + M = v_n / a + p_plus = 0.5 * (1 + equations.gamma * M) * p # HOPE style + # p_plus = 0.5 * (1 + M) * p # Liou-Steffan style + # p_plus = 0.25 * (M + 1)^2 * (2 - M) * p # van Leer style + + f1p = 0.25 * rho * a * (M + 1)^2 + f2p = f1p * v1 + normal_direction[1] * p_plus + f3p = f1p * v2 + normal_direction[2] * p_plus + f4p = f1p * H + + return SVector(f1p, f2p, f3p, f4p)# * s_hat + + # ## yet another form taken from Palmer (https://doi.org/10.2514/3.26178) + # s_hat = 1.0 # norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + # g = equations.gamma + # f1p = s_hat * 0.25 * rho * (v_n + a)^2 / a + # f2p = f1p * (v1 + normal_vector[1] / g * (2 * a - v_n)) + # f3p = f1p * (v2 + normal_vector[2] / g * (2 * a - v_n)) + # f4p = f1p * (0.5 * (v1^2 + v2^2) + (2 * a^2 + 2 * (g - 1) * v_n * a - (g - 1) * v_n^2)/(g^2 - 1)) + + # return SVector(f1p, f2p, f3p, f4p) + + # ## this is actually Zha and Bilgen but just put it in here for now + # # s_hat = norm(normal_direction) + # # normal_vector = normal_direction / s_hat + # # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + # p_plus = 0.5 * p * (1 + v_n / a) + # phi_plus = 0.5 * (v_n + a) + + # f1p = max(v_n, zero(eltype(u))) * rho + # f2p = max(v_n, zero(eltype(u))) * rho_v1 + normal_direction[1] * p_plus + # f3p = max(v_n, zero(eltype(u))) * rho_v2 + normal_direction[2] * p_plus + # f4p = max(v_n, zero(eltype(u))) * rho_e + p * phi_plus + + # return SVector(f1p, f2p, f3p, f4p) + + # ## alternative scaling from the paper of Drikakis and Tsangris + # # s_hat = norm(normal_direction) + # # normal_vector = normal_direction / s_hat + # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + # f1p = 0.25 * rho * a * (v_n / a + 1)^2 + # f2p = f1p * (v1 - normal_direction[1] / equations.gamma * (v_n - 2 * a)) + # f3p = f1p * (v2 - normal_direction[2] / equations.gamma * (v_n - 2 * a)) + # f4p = f1p * H + + # return SVector(f1p, f2p, f3p, f4p)# * s_hat +end + +@inline function splitting_vanleer_haenel(u, ::Val{:minus}, + normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + + a = sqrt(equations.gamma * p / rho) + H = (rho_e + p) / rho + + ## equivalent to the rotate u, calc flux, and then backrotate strategy + # s_hat = norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + M = v_n / a + p_minus = 0.5 * (1 - equations.gamma * M) * p # HOPE style + # p_minus = 0.5 * (1 - M) * p # Liou-Steffan style + # p_minus = 0.25 * (M - 1)^2 * (2 + M) * p # van Leer style + + f1m = -0.25 * rho * a * (M - 1)^2 + f2m = f1m * v1 + normal_direction[1] * p_minus + f3m = f1m * v2 + normal_direction[2] * p_minus + f4m = f1m * H + + return SVector(f1m, f2m, f3m, f4m)# * s_hat + + # ## yet another form taken from Palmer (https://doi.org/10.2514/3.26178) + # s_hat = 1.0 # norm(normal_direction) + # normal_vector = normal_direction / s_hat + # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + # g = equations.gamma + # f1m = -s_hat * 0.25 * rho * (v_n - a)^2 / a + # f2m = f1m * (v1 + normal_vector[1] / g * (-2 * a - v_n)) + # f3m = f1m * (v2 + normal_vector[2] / g * (-2 * a - v_n)) + # f4m = f1m * (0.5 * (v1^2 + v2^2) + (2 * a^2 - 2 * (g - 1) * v_n * a - (g - 1) * v_n^2)/(g^2 - 1)) + + # return SVector(f1m, f2m, f3m, f4m) + + # ## this is actually Zha and Bilgen but just put it in here for now + # # s_hat = norm(normal_direction) + # # normal_vector = normal_direction / s_hat + # # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 + # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + # p_minus = 0.5 * p * (1 - v_n / a) + # phi_minus = 0.5 * (v_n - a) + + # f1m = min(v_n, zero(eltype(u))) * rho + # f2m = min(v_n, zero(eltype(u))) * rho_v1 + normal_direction[1] * p_minus + # f3m = min(v_n, zero(eltype(u))) * rho_v2 + normal_direction[2] * p_minus + # f4m = min(v_n, zero(eltype(u))) * rho_e + p * phi_minus + + # return SVector(f1m, f2m, f3m, f4m) + + # ## alternative scaling from the paper of Drikakis and Tsangris + # # s_hat = norm(normal_direction) + # # normal_vector = normal_direction / s_hat + # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + + # f1m = -0.25 * rho * a * (v_n / a - 1)^2 + # f2m = f1m * (v1 - normal_direction[1] / equations.gamma * (v_n + 2 * a)) + # f3m = f1m * (v2 - normal_direction[2] / equations.gamma * (v_n + 2 * a)) + # f4m = f1m * H + + # return SVector(f1m, f2m, f3m, f4m)# * s_hat +end + """ splitting_lax_friedrichs(u, orientation::Integer, equations::CompressibleEulerEquations2D) @@ -1089,10 +1301,10 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() !!! warning "Experimental implementation (upwind SBP)" This is an experimental feature and may change in future releases. """ -@inline function splitting_lax_friedrichs(u, orientation::Integer, +@inline function splitting_lax_friedrichs(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) - fm = splitting_lax_friedrichs(u, Val{:minus}(), orientation, equations) - fp = splitting_lax_friedrichs(u, Val{:plus}(), orientation, equations) + fm = splitting_lax_friedrichs(u, Val{:minus}(), orientation_or_normal_direction, equations) + fp = splitting_lax_friedrichs(u, Val{:plus}(), orientation_or_normal_direction, equations) return fm, fp end @@ -1150,6 +1362,46 @@ end return SVector(f1m, f2m, f3m, f4m) end +@inline function splitting_lax_friedrichs(u, ::Val{:plus}, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + rho_e = last(u) + rho, v1, v2, p = cons2prim(u, equations) + + a = sqrt(equations.gamma * p / rho) + H = (rho_e + p) / rho + lambda = 0.5 * (sqrt(v1^2 + v2^2) + a) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + + f1p = 0.5 * rho_v_normal + lambda * u[1] + f2p = 0.5 * rho_v_normal * v1 + 0.5 * p * normal_direction[1] + lambda * u[2] + f3p = 0.5 * rho_v_normal * v2 + 0.5 * p * normal_direction[2] + lambda * u[3] + f4p = 0.5 * rho_v_normal * H + lambda * u[4] + + return SVector(f1p, f2p, f3p, f4p) +end + +@inline function splitting_lax_friedrichs(u, ::Val{:minus}, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + rho_e = last(u) + rho, v1, v2, p = cons2prim(u, equations) + + a = sqrt(equations.gamma * p / rho) + H = (rho_e + p) / rho + lambda = 0.5 * (sqrt(v1^2 + v2^2) + a) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + + f1m = 0.5 * rho_v_normal - lambda * u[1] + f2m = 0.5 * rho_v_normal * v1 + 0.5 * p * normal_direction[1] - lambda * u[2] + f3m = 0.5 * rho_v_normal * v2 + 0.5 * p * normal_direction[2] - lambda * u[3] + f4m = 0.5 * rho_v_normal * H - lambda * u[4] + + return SVector(f1m, f2m, f3m, f4m) +end + # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the # maximum velocity magnitude plus the maximum speed of sound @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 32762b37d14..49f3f3eef6d 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -444,8 +444,8 @@ end return fm + fp end -# TODO: Upwind FD. Figure out a better strategy to compute this -# (From an old PR and seemed to not work correctly!) +# 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}) diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl index f51dd09801b..97599f7f0a0 100644 --- a/src/solvers/dgsem_unstructured/containers_2d.jl +++ b/src/solvers/dgsem_unstructured/containers_2d.jl @@ -5,6 +5,103 @@ @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] @@ -147,7 +244,7 @@ end @inline nnodes(interfaces::UnstructuredInterfaceContainer2D) = size(interfaces.u, 3) function init_interfaces(mesh::UnstructuredMesh2D, - elements::UnstructuredElementContainer2D) + elements::Union{UnstructuredElementContainer2D,UpwindElementContainer2D}) interfaces = UnstructuredInterfaceContainer2D{eltype(elements)}(mesh.n_interfaces, nvariables(elements), nnodes(elements)) @@ -288,7 +385,7 @@ end end function init_boundaries(mesh::UnstructuredMesh2D, - elements::UnstructuredElementContainer2D) + elements::Union{UnstructuredElementContainer2D,UpwindElementContainer2D}) boundaries = UnstructuredBoundaryContainer2D{real(elements), eltype(elements)}(mesh.n_boundaries, nvariables(elements), nnodes(elements)) diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index b12a96c4c31..988e995d6b7 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -77,7 +77,7 @@ function rhs!(du, u, t, end # Apply Jacobian from mapping to reference element - # Note! this routine is reused from dg_curved/dg_2d.jl + # Note! this routine is reused from dgsem_structured/dg_2d.jl @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache) # Calculate source terms diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl index d83d4d37557..73ac0e33974 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -29,6 +29,34 @@ 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 @@ -49,24 +77,24 @@ function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeO # jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element] # but uses only matrix-vector products instead of a matrix-matrix product. for j in eachnode(D_SBP) - mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP.central, + mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP, view(node_coordinates, 1, :, j, element)) end # jacobian_matrix[2, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[2, :, :, element] for j in eachnode(D_SBP) - mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP.central, + mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP, view(node_coordinates, 2, :, j, element)) end # Compute the eta derivatives by applying transpose of D on the right # jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(D_SBP)' for i in eachnode(D_SBP) - mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP.central, + mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP, view(node_coordinates, 1, i, :, element)) end # jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(D_SBP)' for i in eachnode(D_SBP) - mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP.central, + mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP, view(node_coordinates, 2, i, :, element)) end @@ -124,4 +152,42 @@ 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 diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index 4ee2c873641..c24b056228b 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -98,15 +98,91 @@ function calc_volume_integral!(du, u, nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, dg::FDSBP, cache) + + # ###################### + # # This is an attempt at the encapsulated version of the USBP operators. The standard way + # # of doing this from Åland cannot be done because it is NOT true that + # # ftilde^\pm = Ja_1 f1^\pm + Ja_2 f2^\pm + # # Thus, we explicitly separate the conservative and nonconserative forms of the mapped + # # fluxes. Earlier experiments with the nonconservative form had issues with rotated elements + # # which I suspect will still be the case here. This still appears to be an issue. + # # The other issue is that the conservative form is not FSP and this already spoils any hope + # # that this strategy would be FSP as well. + # ###################### + + # @unpack f_minus_plus_threaded, f_minus_threaded, f_plus_threaded = cache + # @unpack splitting = volume_integral + # @unpack contravariant_vectors, contravariant_vectors_plus, contravariant_vectors_minus = cache.elements + + # D_minus = Matrix(dg.basis.minus) # Upwind SBP D^- derivative operator + # D_plus = Matrix(dg.basis.plus) # Upwind SBP D^+ derivative operator + + # # Use the tensor product structure to compute the discrete derivatives of + # # the fluxes line-by-line and add them to `du` for each element. + # @threaded for element in eachelement(dg, cache) + # for j in eachnode(dg), i in eachnode(dg) + # u_node = get_node_vars(u, equations, dg, i, j, element) + + # # Grab the contravariant vectors Ja^1 computed with D^+ and D^- + # Ja1_plus = get_contravariant_vector(1, contravariant_vectors, i, j, element) + # Ja1_minus = get_contravariant_vector(1, contravariant_vectors, i, j, element) + + # # Grab the contravariant vectors Ja^2 computed with D^+ and D^- + # Ja2_plus = get_contravariant_vector(2, contravariant_vectors, i, j, element) + # Ja2_minus = get_contravariant_vector(2, contravariant_vectors, i, j, element) + + # # specialized FVS terms that use different combinations of the biased metric terms + # ftilde1_minus = splitting(u_node, Val{:minus}(), Ja1_plus, equations) + # ftilde1_plus = splitting(u_node, Val{:plus}(), Ja1_minus, equations) + # f1_minus, f1_plus = splitting(u_node, 1, equations) + + # ftilde2_minus = splitting(u_node, Val{:minus}(), Ja2_plus, equations) + # ftilde2_plus = splitting(u_node, Val{:plus}(), Ja2_minus, equations) + # f2_minus, f2_plus = splitting(u_node, 2, equations) + + # # xi-direction + # for ii in eachnode(dg) + # # Term: 1/2 D^+ ftilde1_minus + # multiply_add_to_node_vars!(du, 0.5*D_plus[ii, i], ftilde1_minus, equations, dg, ii, j, element) + # # Term: 1/2 D^- ftilde1_plus + # multiply_add_to_node_vars!(du, 0.5*D_minus[ii, i], ftilde1_plus, equations, dg, ii, j, element) + # # Term: 1/2 Ja11^+ D^+ f1_minus + # multiply_add_to_node_vars!(du, 0.5*Ja1_plus[1]*D_plus[ii, i], f1_minus, equations, dg, ii, j, element) + # # Term: 1/2 Ja11^- D^- f1_plus + # multiply_add_to_node_vars!(du, 0.5*Ja1_minus[1]*D_minus[ii, i], f1_plus, equations, dg, ii, j, element) + # # Term: 1/2 Ja12^+ D^+ f2_minus + # multiply_add_to_node_vars!(du, 0.5*Ja1_plus[2]*D_plus[ii, i], f2_minus, equations, dg, ii, j, element) + # # Term: 1/2 Ja12^- D^- f2_plus + # multiply_add_to_node_vars!(du, 0.5*Ja1_minus[2]*D_minus[ii, i], f2_plus, equations, dg, ii, j, element) + # end + + # # eta-direction + # for jj in eachnode(dg) + # # Term: 1/2 D^+ ftilde2_minus + # multiply_add_to_node_vars!(du, 0.5*D_plus[jj, j], ftilde2_minus, equations, dg, i, jj, element) + # # Term: 1/2 D^- ftilde2_plus + # multiply_add_to_node_vars!(du, 0.5*D_minus[jj, j], ftilde2_plus, equations, dg, i, jj, element) + # # Term: 1/2 Ja21^+ D^+ f1_minus + # multiply_add_to_node_vars!(du, 0.5*Ja2_plus[1]*D_plus[jj, j], f1_minus, equations, dg, i, jj, element) + # # Term: 1/2 Ja21^+ D^- f1_plus + # multiply_add_to_node_vars!(du, 0.5*Ja2_minus[1]*D_minus[jj, j], f1_plus, equations, dg, i, jj, element) + # # Term: 1/2 Ja22^+ D^+ f2_minus + # multiply_add_to_node_vars!(du, 0.5*Ja2_plus[2]*D_plus[jj, j], f2_minus, equations, dg, i, jj, element) + # # Term: 1/2 Ja22^+ D^- f2_plus + # multiply_add_to_node_vars!(du, 0.5*Ja2_minus[2]*D_minus[jj, j], f2_plus, equations, dg, i, jj, element) + # end + # end + # end + + # VERSION THAT IS MORE OR LESS WORKING # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator D_plus = dg.basis.plus # Upwind SBP D^+ derivative operator - # D_minus = dg.basis.central # Central SBP D derivative operator - # D_plus = dg.basis.central # Central SBP D derivative operator + D_central = dg.basis.central # Central operator 0.5 (D^+ + D^-) @unpack f_minus_plus_threaded, f_minus_threaded, f_plus_threaded = cache @unpack splitting = volume_integral - @unpack contravariant_vectors = cache.elements + @unpack contravariant_vectors, contravariant_vectors_plus, contravariant_vectors_minus = cache.elements # SBP operators from SummationByPartsOperators.jl implement the basic interface # of matrix-vector multiplication. Thus, we pass an "array of structures", @@ -139,95 +215,162 @@ function calc_volume_integral!(du, u, # x direction # @. f_minus_plus_element = splitting(u_element, 1, equations) + for j in eachnode(dg), i in eachnode(dg) + # contravariant vector and flux computed with central D matrix + # OBS! converges for MMS on flipped mesh but not FSP + Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja1, equations) + + # # For testing, use biased metric terms and rotationally invariant splitting + # # OBS! converges for MMS on flipped mesh but not FSP + # Ja1 = get_contravariant_vector(1, contravariant_vectors_plus, i, j, element) + # f_minus_element[i, j] = splitting(u_element[i, j], Val{:minus}(), Ja1, equations) + + # Ja1 = get_contravariant_vector(1, contravariant_vectors_minus, i, j, element) + # f_plus_element[i, j] = splitting(u_element[i, j], Val{:plus}(), Ja1, equations) + + # # Testing, the f_minus stores the actually flux and f_plus stores the difference + # # between f^- and f^+ + # Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + # f_minus_element[i, j] = flux(u_element[i, j], Ja1, equations) + # f_plus_element[i, j] = (splitting(u_element[i, j], Val{:minus}(), Ja1, equations) + # - splitting(u_element[i, j], Val{:plus}(), Ja1, equations)) + end + for j in eachnode(dg) - for i in eachnode(dg) - Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) - f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja1, equations) - end + # first instinct upwind differences mul!(view(du_vectors, :, j, element), D_minus, view(f_plus_element, :, j), one(eltype(du)), one(eltype(du))) mul!(view(du_vectors, :, j, element), D_plus, view(f_minus_element, :, j), one(eltype(du)), one(eltype(du))) + + # # use a combination of the central and diffusion matrices that comes from, e.g., + # # D⁺ f⁻ + D⁻ f⁺ = D⁺ f⁻ + D⁻ (f - f⁻) + # # = D⁺ (0.5 f + f⁻ - 0.5 f) + D⁻ (0.5 f - f⁻ + 0.5 f) + # # = 0.5 (D⁺ + D⁻) f + 0.5 (D⁺ - D⁻) (2 f⁻ - f) + # # = Dᶜ f + 0.5 (D⁺ - D⁻) (f⁻ - f⁺) + # mul!(view(du_vectors, :, j, element), dg.basis.central, view(f_minus_element, :, j), + # one(eltype(du)), one(eltype(du))) + # mul!(view(du_vectors, :, j, element), D_plus, view(f_plus_element, :, j), + # 0.5*one(eltype(du)), one(eltype(du))) + # mul!(view(du_vectors, :, j, element), D_minus, view(f_plus_element, :, j), + # -0.5*one(eltype(du)), one(eltype(du))) end # y direction # @. f_minus_plus_element = splitting(u_element, 2, equations) + for j in eachnode(dg), i in eachnode(dg) + # contravariant vector and flux computed with central D matrix (not FSP) + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja2, equations) + + # # For testing, use biased metric terms and rotationally invariant splitting + # Ja2 = get_contravariant_vector(2, contravariant_vectors_plus, i, j, element) + # f_minus_element[i, j] = splitting(u_element[i, j], Val{:minus}(), Ja2, equations) + + # Ja2 = get_contravariant_vector(2, contravariant_vectors_minus, i, j, element) + # f_plus_element[i, j] = splitting(u_element[i, j], Val{:plus}(), Ja2, equations) + + # # Testing, the f_minus stores the full flux and f_plus stores the difference + # # between f^- and f^+ + # Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + # f_minus_element[i, j] = flux(u_element[i, j], Ja2, equations) + # f_plus_element[i, j] = (splitting(u_element[i, j], Val{:minus}(), Ja2, equations) + # - splitting(u_element[i, j], Val{:plus}(), Ja2, equations)) + end + for i in eachnode(dg) - for j in eachnode(dg) - Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) - f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja2, equations) - end + # first instinct upwind differences mul!(view(du_vectors, i, :, element), D_minus, view(f_plus_element, i, :), one(eltype(du)), one(eltype(du))) mul!(view(du_vectors, i, :, element), D_plus, view(f_minus_element, i, :), one(eltype(du)), one(eltype(du))) - end - end - - return nothing -end -# TODO: Write a better comment here. Specialized computation of the interface -# values using flux vector splitting in an arbitrary direction -function calc_interface_flux!(surface_flux_values, - mesh::UnstructuredMesh2D, - nonconservative_terms::False, equations, - surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache) - @unpack splitting = surface_integral - @unpack u, start_index, index_increment, element_ids, element_side_ids = cache.interfaces - @unpack normal_directions = cache.elements - - @threaded for interface in eachinterface(dg, cache) - # Get neighboring elements - primary_element = element_ids[1, interface] - secondary_element = element_ids[2, interface] - - # Get the local side id on which to compute the flux - primary_side = element_side_ids[1, interface] - secondary_side = element_side_ids[2, interface] - - # initial index for the coordinate system on the secondary element - secondary_index = start_index[interface] - - # loop through the primary element coordinate system and compute the interface coupling - for primary_index in eachnode(dg) - # Pull the primary and secondary states from the boundary u values - u_ll = get_one_sided_surface_node_vars(u, equations, dg, 1, primary_index, - interface) - u_rr = get_one_sided_surface_node_vars(u, equations, dg, 2, secondary_index, - interface) - - # pull the outward pointing (normal) directional vector - # Note! this assumes a conforming approximation, more must be done in terms of the normals - # for hanging nodes and other non-conforming approximation spaces - outward_direction = get_surface_normal(normal_directions, primary_index, - primary_side, - primary_element) - - # Compute the upwind coupling terms where right-traveling - # information comes from the left and left-traveling information - # comes from the right - flux_minus_rr = splitting(u_rr, Val{:minus}(), outward_direction, equations) - flux_plus_ll = splitting(u_ll, Val{:plus}(), outward_direction, equations) - - # Copy upwind fluxes back to primary/secondary element storage - # Note the sign change for the normal flux in the secondary element! - # TODO: do we need this sign flip on the neighbour? - # OBS! if we switch the sign flip here then the `SurfaceIntegralUpwind` also needs adjusted - for v in eachvariable(equations) - surface_flux_values[v, primary_index, primary_side, primary_element] = flux_minus_rr[v] - surface_flux_values[v, secondary_index, secondary_side, secondary_element] = -flux_plus_ll[v] - end - - # increment the index of the coordinate system in the secondary element - secondary_index += index_increment[interface] + # # use a combination of the central and diffusion matrices + # mul!(view(du_vectors, i, :, element), dg.basis.central, view(f_minus_element, i, :), + # one(eltype(du)), one(eltype(du))) + # mul!(view(du_vectors, i, :, element), D_plus, view(f_plus_element, i, :), + # 0.5*one(eltype(du)), one(eltype(du))) + # mul!(view(du_vectors, i, :, element), D_minus, view(f_plus_element, i, :), + # -0.5*one(eltype(du)), one(eltype(du))) end end return nothing end +# # TODO: Write a better comment here. Specialized computation of the interface +# # values using flux vector splitting in an arbitrary direction +# # TODO: OBS! `SurfaceIntegralUpwind` is buggy and possibly unnecessary. +# function calc_interface_flux!(surface_flux_values, +# mesh::UnstructuredMesh2D, +# nonconservative_terms::False, equations, +# surface_integral::SurfaceIntegralUpwind, +# dg::FDSBP, cache) +# @unpack splitting = surface_integral +# @unpack u, start_index, index_increment, element_ids, element_side_ids = cache.interfaces +# @unpack normal_directions, normal_directions_plus, normal_directions_minus = cache.elements + +# @threaded for interface in eachinterface(dg, cache) +# # Get neighboring elements +# primary_element = element_ids[1, interface] +# secondary_element = element_ids[2, interface] + +# # Get the local side id on which to compute the flux +# primary_side = element_side_ids[1, interface] +# secondary_side = element_side_ids[2, interface] + +# # initial index for the coordinate system on the secondary element +# secondary_index = start_index[interface] + +# # loop through the primary element coordinate system and compute the interface coupling +# for primary_index in eachnode(dg) +# # Pull the primary and secondary states from the boundary u values +# u_ll = get_one_sided_surface_node_vars(u, equations, dg, 1, primary_index, +# interface) +# u_rr = get_one_sided_surface_node_vars(u, equations, dg, 2, secondary_index, +# interface) + +# # pull the outward pointing (normal) directional vectors +# # Note! this assumes a conforming approximation, more must be done in terms of the normals +# # for hanging nodes and other non-conforming approximation spaces +# outward_direction = get_surface_normal(normal_directions, primary_index, +# primary_side, +# primary_element) +# # outward_direction_plus = get_surface_normal(normal_directions_plus, +# # primary_index, +# # primary_side, +# # primary_element) +# # outward_direction_minus = get_surface_normal(normal_directions_minus, +# # primary_index, +# # primary_side, +# # primary_element) + +# # Compute the upwind coupling terms where right-traveling +# # information comes from the left and left-traveling information +# # comes from the right +# flux_minus_rr = splitting(u_rr, Val{:minus}(), outward_direction, equations) +# flux_plus_ll = splitting(u_ll, Val{:plus}(), outward_direction, equations) +# # flux_minus_rr = splitting(u_rr, Val{:minus}(), outward_direction_plus, equations) +# # flux_plus_ll = splitting(u_ll, Val{:plus}(), outward_direction_minus, equations) + +# # Copy upwind fluxes back to primary/secondary element storage +# # Note the sign change for the normal flux in the secondary element! +# # TODO: do we need this sign flip on the neighbour? +# # OBS! if we switch the sign flip here then the `SurfaceIntegralUpwind` also needs adjusted +# for v in eachvariable(equations) +# surface_flux_values[v, primary_index, primary_side, primary_element] = flux_minus_rr[v] +# surface_flux_values[v, secondary_index, secondary_side, secondary_element] = -flux_plus_ll[v] +# end + +# # increment the index of the coordinate system in the secondary element +# secondary_index += index_increment[interface] +# end +# end + +# return nothing +# end + # Note! The local side numbering for the unstructured quadrilateral element implementation differs # from the structured TreeMesh or StructuredMesh local side numbering: # @@ -256,8 +399,10 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at -x u_node = get_node_vars(u, equations, dg, 1, l, element) # compute internal flux in normal direction on side 4 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, - element) + # outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, + # element) + # TODO: can probably replace the above with the simpler (same goes with the other instances below) + outward_direction = get_surface_normal(normal_directions, l, 4, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), @@ -266,8 +411,9 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at +x u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) # compute internal flux in normal direction on side 2 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, - element) + # outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, + # element) + outward_direction = get_surface_normal(normal_directions, l, 2, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), @@ -276,8 +422,9 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at -y u_node = get_node_vars(u, equations, dg, l, 1, element) # compute internal flux in normal direction on side 1 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, - element) + # outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, + # element) + outward_direction = get_surface_normal(normal_directions, l, 1, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), @@ -286,8 +433,9 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at +y u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) # compute internal flux in normal direction on side 3 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, - element) + # outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, + # element) + outward_direction = get_surface_normal(normal_directions, l, 3, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), @@ -298,64 +446,69 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, return nothing end -# Implementation of fully upwind SATs. The surface flux values are pre-computed -# in the specialized `calc_interface_flux` routine. These SATs are still of -# a strong form penalty type, except that the interior flux at a particular -# side of the element are computed in the upwind direction. -function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, - equations, surface_integral::SurfaceIntegralUpwind, - dg::FDSBP, cache) - inv_weight_left = inv(left_boundary_weight(dg.basis)) - inv_weight_right = inv(right_boundary_weight(dg.basis)) - @unpack normal_directions, surface_flux_values = cache.elements - @unpack splitting = surface_integral - - @threaded for element in eachelement(dg, cache) - for l in eachnode(dg) - # surface at -x - u_node = get_node_vars(u, equations, dg, 1, l, element) - # compute internal flux in normal direction on side 4 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, - element) - f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) - f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) - multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), - equations, dg, 1, l, element) - - # surface at +x - u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) - # compute internal flux in normal direction on side 2 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, - element) - f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) - f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) - multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), - equations, dg, nnodes(dg), l, element) - - # surface at -y - u_node = get_node_vars(u, equations, dg, l, 1, element) - # compute internal flux in normal direction on side 1 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, - element) - f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) - f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) - multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), - equations, dg, l, 1, element) - - # surface at +y - u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) - # compute internal flux in normal direction on side 3 - outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, - element) - f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) - f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) - multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), - equations, dg, l, nnodes(dg), element) - end - end - - return nothing -end +# # Implementation of fully upwind SATs. The surface flux values are pre-computed +# # in the specialized `calc_interface_flux` routine. These SATs are still of +# # a strong form penalty type, except that the interior flux at a particular +# # side of the element are computed in the upwind direction. +# # TODO: OBS! `SurfaceIntegralUpwind` is buggy and possibly unnecessary. +# function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, +# equations, surface_integral::SurfaceIntegralUpwind, +# dg::FDSBP, cache) +# inv_weight_left = inv(left_boundary_weight(dg.basis)) +# inv_weight_right = inv(right_boundary_weight(dg.basis)) +# @unpack normal_directions, normal_directions_plus, normal_directions_minus, surface_flux_values = cache.elements +# @unpack splitting = surface_integral + +# @threaded for element in eachelement(dg, cache) +# for l in eachnode(dg) +# # surface at -x +# u_node = get_node_vars(u, equations, dg, 1, l, element) +# # compute internal flux in normal direction on side 4 +# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, +# # element) +# outward_direction = get_surface_normal(normal_directions, l, 4, element) +# f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) +# f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) +# multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), +# equations, dg, 1, l, element) + +# # surface at +x +# u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) +# # compute internal flux in normal direction on side 2 +# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, +# # element) +# outward_direction = get_surface_normal(normal_directions, l, 2, element) +# f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) +# f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) +# multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), +# equations, dg, nnodes(dg), l, element) + +# # surface at -y +# u_node = get_node_vars(u, equations, dg, l, 1, element) +# # compute internal flux in normal direction on side 1 +# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, +# # element) +# outward_direction = get_surface_normal(normal_directions, l, 1, element) +# f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) +# f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) +# multiply_add_to_node_vars!(du, inv_weight_left, -(f_num - f_node), +# equations, dg, l, 1, element) + +# # surface at +y +# u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) +# # compute internal flux in normal direction on side 3 +# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, +# # element) +# outward_direction = get_surface_normal(normal_directions, l, 3, element) +# f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) +# f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) +# multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), +# equations, dg, l, nnodes(dg), element) +# end +# end + +# return nothing +# end # AnalysisCallback function integrate_via_indices(func::Func, u, From 4324ec6f25b34c3932082fed480b8b697ee9501e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 4 Mar 2024 12:33:09 +0100 Subject: [PATCH 03/17] cleanup of the fdsbp_2d file --- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 284 +-------------------- 1 file changed, 3 insertions(+), 281 deletions(-) diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index c24b056228b..948c020a1fb 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -25,8 +25,6 @@ function create_cache(mesh::UnstructuredMesh2D, equations, dg::FDSBP, RealT, uEl return cache end -# TODO: FD; Upwind versions of surface / volume integral - # 2D volume integral contributions for `VolumeIntegralStrongForm` # OBS! This is the standard (not de-aliased) form of the volume integral. # So it is not provably stable for variable coefficients due to the the metric terms. @@ -98,91 +96,13 @@ function calc_volume_integral!(du, u, nonconservative_terms::False, equations, volume_integral::VolumeIntegralUpwind, dg::FDSBP, cache) - - # ###################### - # # This is an attempt at the encapsulated version of the USBP operators. The standard way - # # of doing this from Åland cannot be done because it is NOT true that - # # ftilde^\pm = Ja_1 f1^\pm + Ja_2 f2^\pm - # # Thus, we explicitly separate the conservative and nonconserative forms of the mapped - # # fluxes. Earlier experiments with the nonconservative form had issues with rotated elements - # # which I suspect will still be the case here. This still appears to be an issue. - # # The other issue is that the conservative form is not FSP and this already spoils any hope - # # that this strategy would be FSP as well. - # ###################### - - # @unpack f_minus_plus_threaded, f_minus_threaded, f_plus_threaded = cache - # @unpack splitting = volume_integral - # @unpack contravariant_vectors, contravariant_vectors_plus, contravariant_vectors_minus = cache.elements - - # D_minus = Matrix(dg.basis.minus) # Upwind SBP D^- derivative operator - # D_plus = Matrix(dg.basis.plus) # Upwind SBP D^+ derivative operator - - # # Use the tensor product structure to compute the discrete derivatives of - # # the fluxes line-by-line and add them to `du` for each element. - # @threaded for element in eachelement(dg, cache) - # for j in eachnode(dg), i in eachnode(dg) - # u_node = get_node_vars(u, equations, dg, i, j, element) - - # # Grab the contravariant vectors Ja^1 computed with D^+ and D^- - # Ja1_plus = get_contravariant_vector(1, contravariant_vectors, i, j, element) - # Ja1_minus = get_contravariant_vector(1, contravariant_vectors, i, j, element) - - # # Grab the contravariant vectors Ja^2 computed with D^+ and D^- - # Ja2_plus = get_contravariant_vector(2, contravariant_vectors, i, j, element) - # Ja2_minus = get_contravariant_vector(2, contravariant_vectors, i, j, element) - - # # specialized FVS terms that use different combinations of the biased metric terms - # ftilde1_minus = splitting(u_node, Val{:minus}(), Ja1_plus, equations) - # ftilde1_plus = splitting(u_node, Val{:plus}(), Ja1_minus, equations) - # f1_minus, f1_plus = splitting(u_node, 1, equations) - - # ftilde2_minus = splitting(u_node, Val{:minus}(), Ja2_plus, equations) - # ftilde2_plus = splitting(u_node, Val{:plus}(), Ja2_minus, equations) - # f2_minus, f2_plus = splitting(u_node, 2, equations) - - # # xi-direction - # for ii in eachnode(dg) - # # Term: 1/2 D^+ ftilde1_minus - # multiply_add_to_node_vars!(du, 0.5*D_plus[ii, i], ftilde1_minus, equations, dg, ii, j, element) - # # Term: 1/2 D^- ftilde1_plus - # multiply_add_to_node_vars!(du, 0.5*D_minus[ii, i], ftilde1_plus, equations, dg, ii, j, element) - # # Term: 1/2 Ja11^+ D^+ f1_minus - # multiply_add_to_node_vars!(du, 0.5*Ja1_plus[1]*D_plus[ii, i], f1_minus, equations, dg, ii, j, element) - # # Term: 1/2 Ja11^- D^- f1_plus - # multiply_add_to_node_vars!(du, 0.5*Ja1_minus[1]*D_minus[ii, i], f1_plus, equations, dg, ii, j, element) - # # Term: 1/2 Ja12^+ D^+ f2_minus - # multiply_add_to_node_vars!(du, 0.5*Ja1_plus[2]*D_plus[ii, i], f2_minus, equations, dg, ii, j, element) - # # Term: 1/2 Ja12^- D^- f2_plus - # multiply_add_to_node_vars!(du, 0.5*Ja1_minus[2]*D_minus[ii, i], f2_plus, equations, dg, ii, j, element) - # end - - # # eta-direction - # for jj in eachnode(dg) - # # Term: 1/2 D^+ ftilde2_minus - # multiply_add_to_node_vars!(du, 0.5*D_plus[jj, j], ftilde2_minus, equations, dg, i, jj, element) - # # Term: 1/2 D^- ftilde2_plus - # multiply_add_to_node_vars!(du, 0.5*D_minus[jj, j], ftilde2_plus, equations, dg, i, jj, element) - # # Term: 1/2 Ja21^+ D^+ f1_minus - # multiply_add_to_node_vars!(du, 0.5*Ja2_plus[1]*D_plus[jj, j], f1_minus, equations, dg, i, jj, element) - # # Term: 1/2 Ja21^+ D^- f1_plus - # multiply_add_to_node_vars!(du, 0.5*Ja2_minus[1]*D_minus[jj, j], f1_plus, equations, dg, i, jj, element) - # # Term: 1/2 Ja22^+ D^+ f2_minus - # multiply_add_to_node_vars!(du, 0.5*Ja2_plus[2]*D_plus[jj, j], f2_minus, equations, dg, i, jj, element) - # # Term: 1/2 Ja22^+ D^- f2_plus - # multiply_add_to_node_vars!(du, 0.5*Ja2_minus[2]*D_minus[jj, j], f2_plus, equations, dg, i, jj, element) - # end - # end - # end - - # VERSION THAT IS MORE OR LESS WORKING # Assume that # dg.basis isa SummationByPartsOperators.UpwindOperators D_minus = dg.basis.minus # Upwind SBP D^- derivative operator D_plus = dg.basis.plus # Upwind SBP D^+ derivative operator - D_central = dg.basis.central # Central operator 0.5 (D^+ + D^-) @unpack f_minus_plus_threaded, f_minus_threaded, f_plus_threaded = cache @unpack splitting = volume_integral - @unpack contravariant_vectors, contravariant_vectors_plus, contravariant_vectors_minus = cache.elements + @unpack contravariant_vectors = cache.elements # SBP operators from SummationByPartsOperators.jl implement the basic interface # of matrix-vector multiplication. Thus, we pass an "array of structures", @@ -216,161 +136,36 @@ function calc_volume_integral!(du, u, # x direction # @. f_minus_plus_element = splitting(u_element, 1, equations) for j in eachnode(dg), i in eachnode(dg) - # contravariant vector and flux computed with central D matrix - # OBS! converges for MMS on flipped mesh but not FSP + # contravariant vectors computed with central D matrix Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja1, equations) - - # # For testing, use biased metric terms and rotationally invariant splitting - # # OBS! converges for MMS on flipped mesh but not FSP - # Ja1 = get_contravariant_vector(1, contravariant_vectors_plus, i, j, element) - # f_minus_element[i, j] = splitting(u_element[i, j], Val{:minus}(), Ja1, equations) - - # Ja1 = get_contravariant_vector(1, contravariant_vectors_minus, i, j, element) - # f_plus_element[i, j] = splitting(u_element[i, j], Val{:plus}(), Ja1, equations) - - # # Testing, the f_minus stores the actually flux and f_plus stores the difference - # # between f^- and f^+ - # Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) - # f_minus_element[i, j] = flux(u_element[i, j], Ja1, equations) - # f_plus_element[i, j] = (splitting(u_element[i, j], Val{:minus}(), Ja1, equations) - # - splitting(u_element[i, j], Val{:plus}(), Ja1, equations)) end for j in eachnode(dg) - # first instinct upwind differences mul!(view(du_vectors, :, j, element), D_minus, view(f_plus_element, :, j), one(eltype(du)), one(eltype(du))) mul!(view(du_vectors, :, j, element), D_plus, view(f_minus_element, :, j), one(eltype(du)), one(eltype(du))) - - # # use a combination of the central and diffusion matrices that comes from, e.g., - # # D⁺ f⁻ + D⁻ f⁺ = D⁺ f⁻ + D⁻ (f - f⁻) - # # = D⁺ (0.5 f + f⁻ - 0.5 f) + D⁻ (0.5 f - f⁻ + 0.5 f) - # # = 0.5 (D⁺ + D⁻) f + 0.5 (D⁺ - D⁻) (2 f⁻ - f) - # # = Dᶜ f + 0.5 (D⁺ - D⁻) (f⁻ - f⁺) - # mul!(view(du_vectors, :, j, element), dg.basis.central, view(f_minus_element, :, j), - # one(eltype(du)), one(eltype(du))) - # mul!(view(du_vectors, :, j, element), D_plus, view(f_plus_element, :, j), - # 0.5*one(eltype(du)), one(eltype(du))) - # mul!(view(du_vectors, :, j, element), D_minus, view(f_plus_element, :, j), - # -0.5*one(eltype(du)), one(eltype(du))) end # y direction - # @. f_minus_plus_element = splitting(u_element, 2, equations) for j in eachnode(dg), i in eachnode(dg) - # contravariant vector and flux computed with central D matrix (not FSP) + # contravariant vectors computed with central D matrix Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja2, equations) - - # # For testing, use biased metric terms and rotationally invariant splitting - # Ja2 = get_contravariant_vector(2, contravariant_vectors_plus, i, j, element) - # f_minus_element[i, j] = splitting(u_element[i, j], Val{:minus}(), Ja2, equations) - - # Ja2 = get_contravariant_vector(2, contravariant_vectors_minus, i, j, element) - # f_plus_element[i, j] = splitting(u_element[i, j], Val{:plus}(), Ja2, equations) - - # # Testing, the f_minus stores the full flux and f_plus stores the difference - # # between f^- and f^+ - # Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) - # f_minus_element[i, j] = flux(u_element[i, j], Ja2, equations) - # f_plus_element[i, j] = (splitting(u_element[i, j], Val{:minus}(), Ja2, equations) - # - splitting(u_element[i, j], Val{:plus}(), Ja2, equations)) end for i in eachnode(dg) - # first instinct upwind differences mul!(view(du_vectors, i, :, element), D_minus, view(f_plus_element, i, :), one(eltype(du)), one(eltype(du))) mul!(view(du_vectors, i, :, element), D_plus, view(f_minus_element, i, :), one(eltype(du)), one(eltype(du))) - - # # use a combination of the central and diffusion matrices - # mul!(view(du_vectors, i, :, element), dg.basis.central, view(f_minus_element, i, :), - # one(eltype(du)), one(eltype(du))) - # mul!(view(du_vectors, i, :, element), D_plus, view(f_plus_element, i, :), - # 0.5*one(eltype(du)), one(eltype(du))) - # mul!(view(du_vectors, i, :, element), D_minus, view(f_plus_element, i, :), - # -0.5*one(eltype(du)), one(eltype(du))) end end return nothing end -# # TODO: Write a better comment here. Specialized computation of the interface -# # values using flux vector splitting in an arbitrary direction -# # TODO: OBS! `SurfaceIntegralUpwind` is buggy and possibly unnecessary. -# function calc_interface_flux!(surface_flux_values, -# mesh::UnstructuredMesh2D, -# nonconservative_terms::False, equations, -# surface_integral::SurfaceIntegralUpwind, -# dg::FDSBP, cache) -# @unpack splitting = surface_integral -# @unpack u, start_index, index_increment, element_ids, element_side_ids = cache.interfaces -# @unpack normal_directions, normal_directions_plus, normal_directions_minus = cache.elements - -# @threaded for interface in eachinterface(dg, cache) -# # Get neighboring elements -# primary_element = element_ids[1, interface] -# secondary_element = element_ids[2, interface] - -# # Get the local side id on which to compute the flux -# primary_side = element_side_ids[1, interface] -# secondary_side = element_side_ids[2, interface] - -# # initial index for the coordinate system on the secondary element -# secondary_index = start_index[interface] - -# # loop through the primary element coordinate system and compute the interface coupling -# for primary_index in eachnode(dg) -# # Pull the primary and secondary states from the boundary u values -# u_ll = get_one_sided_surface_node_vars(u, equations, dg, 1, primary_index, -# interface) -# u_rr = get_one_sided_surface_node_vars(u, equations, dg, 2, secondary_index, -# interface) - -# # pull the outward pointing (normal) directional vectors -# # Note! this assumes a conforming approximation, more must be done in terms of the normals -# # for hanging nodes and other non-conforming approximation spaces -# outward_direction = get_surface_normal(normal_directions, primary_index, -# primary_side, -# primary_element) -# # outward_direction_plus = get_surface_normal(normal_directions_plus, -# # primary_index, -# # primary_side, -# # primary_element) -# # outward_direction_minus = get_surface_normal(normal_directions_minus, -# # primary_index, -# # primary_side, -# # primary_element) - -# # Compute the upwind coupling terms where right-traveling -# # information comes from the left and left-traveling information -# # comes from the right -# flux_minus_rr = splitting(u_rr, Val{:minus}(), outward_direction, equations) -# flux_plus_ll = splitting(u_ll, Val{:plus}(), outward_direction, equations) -# # flux_minus_rr = splitting(u_rr, Val{:minus}(), outward_direction_plus, equations) -# # flux_plus_ll = splitting(u_ll, Val{:plus}(), outward_direction_minus, equations) - -# # Copy upwind fluxes back to primary/secondary element storage -# # Note the sign change for the normal flux in the secondary element! -# # TODO: do we need this sign flip on the neighbour? -# # OBS! if we switch the sign flip here then the `SurfaceIntegralUpwind` also needs adjusted -# for v in eachvariable(equations) -# surface_flux_values[v, primary_index, primary_side, primary_element] = flux_minus_rr[v] -# surface_flux_values[v, secondary_index, secondary_side, secondary_element] = -flux_plus_ll[v] -# end - -# # increment the index of the coordinate system in the secondary element -# secondary_index += index_increment[interface] -# end -# end - -# return nothing -# end - # Note! The local side numbering for the unstructured quadrilateral element implementation differs # from the structured TreeMesh or StructuredMesh local side numbering: # @@ -399,9 +194,6 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at -x u_node = get_node_vars(u, equations, dg, 1, l, element) # compute internal flux in normal direction on side 4 - # outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, - # element) - # TODO: can probably replace the above with the simpler (same goes with the other instances below) outward_direction = get_surface_normal(normal_directions, l, 4, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) @@ -411,8 +203,6 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at +x u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) # compute internal flux in normal direction on side 2 - # outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, - # element) outward_direction = get_surface_normal(normal_directions, l, 2, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) @@ -422,8 +212,6 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at -y u_node = get_node_vars(u, equations, dg, l, 1, element) # compute internal flux in normal direction on side 1 - # outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, - # element) outward_direction = get_surface_normal(normal_directions, l, 1, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) @@ -433,8 +221,6 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # surface at +y u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) # compute internal flux in normal direction on side 3 - # outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, - # element) outward_direction = get_surface_normal(normal_directions, l, 3, element) f_node = flux(u_node, outward_direction, equations) f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) @@ -446,70 +232,6 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, return nothing end -# # Implementation of fully upwind SATs. The surface flux values are pre-computed -# # in the specialized `calc_interface_flux` routine. These SATs are still of -# # a strong form penalty type, except that the interior flux at a particular -# # side of the element are computed in the upwind direction. -# # TODO: OBS! `SurfaceIntegralUpwind` is buggy and possibly unnecessary. -# function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, -# equations, surface_integral::SurfaceIntegralUpwind, -# dg::FDSBP, cache) -# inv_weight_left = inv(left_boundary_weight(dg.basis)) -# inv_weight_right = inv(right_boundary_weight(dg.basis)) -# @unpack normal_directions, normal_directions_plus, normal_directions_minus, surface_flux_values = cache.elements -# @unpack splitting = surface_integral - -# @threaded for element in eachelement(dg, cache) -# for l in eachnode(dg) -# # surface at -x -# u_node = get_node_vars(u, equations, dg, 1, l, element) -# # compute internal flux in normal direction on side 4 -# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, -# # element) -# outward_direction = get_surface_normal(normal_directions, l, 4, element) -# f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) -# f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) -# multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), -# equations, dg, 1, l, element) - -# # surface at +x -# u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) -# # compute internal flux in normal direction on side 2 -# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, -# # element) -# outward_direction = get_surface_normal(normal_directions, l, 2, element) -# f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) -# f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) -# multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), -# equations, dg, nnodes(dg), l, element) - -# # surface at -y -# u_node = get_node_vars(u, equations, dg, l, 1, element) -# # compute internal flux in normal direction on side 1 -# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, -# # element) -# outward_direction = get_surface_normal(normal_directions, l, 1, element) -# f_node = splitting(u_node, Val{:plus}(), outward_direction, equations) -# f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) -# multiply_add_to_node_vars!(du, inv_weight_left, -(f_num - f_node), -# equations, dg, l, 1, element) - -# # surface at +y -# u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) -# # compute internal flux in normal direction on side 3 -# # outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, -# # element) -# outward_direction = get_surface_normal(normal_directions, l, 3, element) -# f_node = splitting(u_node, Val{:minus}(), outward_direction, equations) -# f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) -# multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), -# equations, dg, l, nnodes(dg), element) -# end -# end - -# return nothing -# end - # AnalysisCallback function integrate_via_indices(func::Func, u, mesh::UnstructuredMesh2D, equations, From d28e814b115e0c0670da1d7d192be4e7f33aa6f4 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 4 Mar 2024 14:14:42 +0100 Subject: [PATCH 04/17] clean-up FVS routines in the compressible Euler file --- src/equations/compressible_euler_2d.jl | 216 +++---------------------- 1 file changed, 21 insertions(+), 195 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 8935cb12eee..69a7f777f49 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -683,13 +683,15 @@ end end """ - splitting_steger_warming(u, orientation::Integer, + splitting_steger_warming(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) splitting_steger_warming(u, which::Union{Val{:minus}, Val{:plus}} - orientation::Integer, + orientation_or_normal_direction, equations::CompressibleEulerEquations2D) -Splitting of the compressible Euler flux of Steger and Warming. +Splitting of the compressible Euler flux of Steger and Warming. For +curvilinear coordinates to compute in the `normal_direction` we use +a modified version of this splitting proposed by Drikakis and Tsangris. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the @@ -705,6 +707,10 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() Flux Vector Splitting of the Inviscid Gasdynamic Equations With Application to Finite Difference Methods [NASA Technical Memorandum](https://ntrs.nasa.gov/api/citations/19790020779/downloads/19790020779.pdf) +- D. Drikakis and S. Tsangaris (1993) + On the solution of the compressible Navier-Stokes equations using + improved flux vector splitting methods + [DOI: 10.1016/0307-904X(93)90054-K](https://doi.org/10.1016/0307-904X(93)90054-K) """ @inline function splitting_steger_warming(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) @@ -817,37 +823,8 @@ end v2 = rho_v2 / rho p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) + H = (rho_e + p) / rho - # ## equivalent to rotate u, do splitting and rotate back - # s_hat = norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 - # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # lambda1 = v_n - # lambda2 = lambda1 + a * s_hat - # lambda3 = lambda1 - a * s_hat - - # lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) - # lambda2_p = positive_part(lambda2) - # lambda3_p = positive_part(lambda3) - - # alpha_p = 2 * (equations.gamma - 1) * lambda1_p + lambda2_p + lambda3_p - - # rho_2gamma = 0.5 * rho / equations.gamma - # f1p = rho_2gamma * alpha_p - # f2p = rho_2gamma * (alpha_p * v1 + a * normal_vector[1] * (lambda2_p - lambda3_p)) - # f3p = rho_2gamma * (alpha_p * v2 + a * normal_vector[2] * (lambda2_p - lambda3_p)) - # f4p = rho_2gamma * - # (alpha_p * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_p - lambda3_p) - # + a^2 * (lambda2_p + lambda3_p) * equations.inv_gamma_minus_one) - - # return SVector(f1p, f2p, f3p, f4p) - - ## alternative scaling from the paper of Drikakis and Tsangris - # s_hat = norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 v_n = normal_direction[1] * v1 + normal_direction[2] * v2 lambda1 = v_n + a @@ -856,25 +833,13 @@ end lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) lambda2_p = positive_part(lambda2) - ## original form from the paper - # H = (rho_e + p) / rho - # f1p = 0.5 * rho * (lambda1_p + lambda2_p) - # f2p = (0.5 * rho * (v1 + normal_direction[1] * a / equations.gamma) * lambda1_p - # + 0.5 * rho * (v1 - normal_direction[1] * a / equations.gamma) * lambda2_p) - # f3p = (0.5 * rho * (v2 + normal_direction[2] * a / equations.gamma) * lambda1_p - # + 0.5 * rho * (v2 - normal_direction[2] * a / equations.gamma) * lambda2_p) - # f4p = f1p * H - - # slightly rewritten form that is cleaner rhoa_2gamma = 0.5 * rho * a / equations.gamma - H = (rho_e + p) / rho - f1p = 0.5 * rho * (lambda1_p + lambda2_p) f2p = f1p * v1 + rhoa_2gamma * normal_direction[1] * (lambda1_p - lambda2_p) f3p = f1p * v2 + rhoa_2gamma * normal_direction[2] * (lambda1_p - lambda2_p) f4p = f1p * H - return SVector(f1p, f2p, f3p, f4p)# * s_hat + return SVector(f1p, f2p, f3p, f4p) end @inline function splitting_steger_warming(u, ::Val{:minus}, @@ -885,37 +850,8 @@ end v2 = rho_v2 / rho p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) + H = (rho_e + p) / rho - # # equivalent to rotate u, do splitting and rotate back - # s_hat = norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 - # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # lambda1 = v_n - # lambda2 = lambda1 + a * s_hat - # lambda3 = lambda1 - a * s_hat - - # lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) - # lambda2_m = negative_part(lambda2) - # lambda3_m = negative_part(lambda3) - - # alpha_m = 2 * (equations.gamma - 1) * lambda1_m + lambda2_m + lambda3_m - - # rho_2gamma = 0.5 * rho / equations.gamma - # f1m = rho_2gamma * alpha_m - # f2m = rho_2gamma * (alpha_m * v1 + a * normal_vector[1] * (lambda2_m - lambda3_m)) - # f3m = rho_2gamma * (alpha_m * v2 + a * normal_vector[2] * (lambda2_m - lambda3_m)) - # f4m = rho_2gamma * - # (alpha_m * 0.5 * (v1^2 + v2^2) + a * v_n * (lambda2_m - lambda3_m) - # + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one) - - # return SVector(f1m, f2m, f3m, f4m) - - ## alternative scaling from the paper of Drikakis and Tsangris - # s_hat = norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 v_n = normal_direction[1] * v1 + normal_direction[2] * v2 lambda1 = v_n + a @@ -924,25 +860,13 @@ end lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) lambda2_m = negative_part(lambda2) - ## original version from the paper - # H = (rho_e + p) / rho - # f1m = 0.5 * rho * (lambda1_m + lambda2_m) - # f2m = (0.5 * rho * (v1 + normal_direction[1] * a / equations.gamma) * lambda1_m - # + 0.5 * rho * (v1 - normal_direction[1] * a / equations.gamma) * lambda2_m) - # f3m = (0.5 * rho * (v2 + normal_direction[2] * a / equations.gamma) * lambda1_m - # + 0.5 * rho * (v2 - normal_direction[2] * a / equations.gamma) * lambda2_m) - # f4m = f1m * H - - # slightly rewritten form that is cleaner rhoa_2gamma = 0.5 * rho * a / equations.gamma - H = (rho_e + p) / rho - f1m = 0.5 * rho * (lambda1_m + lambda2_m) f2m = f1m * v1 + rhoa_2gamma * normal_direction[1] * (lambda1_m - lambda2_m) f3m = f1m * v2 + rhoa_2gamma * normal_direction[2] * (lambda1_m - lambda2_m) f4m = f1m * H - return SVector(f1m, f2m, f3m, f4m)# * s_hat + return SVector(f1m, f2m, f3m, f4m) end """ @@ -1038,10 +962,10 @@ end end """ - splitting_vanleer_haenel(u, orientation::Integer, + splitting_vanleer_haenel(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) splitting_vanleer_haenel(u, which::Union{Val{:minus}, Val{:plus}} - orientation::Integer, + orientation_or_normal_direction, equations::CompressibleEulerEquations2D) Splitting of the compressible Euler flux from van Leer. This splitting further @@ -1149,65 +1073,16 @@ end a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho - ## equivalent to the rotate u, calc flux, and then backrotate - # s_hat = norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - M = v_n / a - p_plus = 0.5 * (1 + equations.gamma * M) * p # HOPE style - # p_plus = 0.5 * (1 + M) * p # Liou-Steffan style - # p_plus = 0.25 * (M + 1)^2 * (2 - M) * p # van Leer style + p_plus = 0.5 * (1 + equations.gamma * M) * p f1p = 0.25 * rho * a * (M + 1)^2 f2p = f1p * v1 + normal_direction[1] * p_plus f3p = f1p * v2 + normal_direction[2] * p_plus f4p = f1p * H - return SVector(f1p, f2p, f3p, f4p)# * s_hat - - # ## yet another form taken from Palmer (https://doi.org/10.2514/3.26178) - # s_hat = 1.0 # norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 - # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # g = equations.gamma - # f1p = s_hat * 0.25 * rho * (v_n + a)^2 / a - # f2p = f1p * (v1 + normal_vector[1] / g * (2 * a - v_n)) - # f3p = f1p * (v2 + normal_vector[2] / g * (2 * a - v_n)) - # f4p = f1p * (0.5 * (v1^2 + v2^2) + (2 * a^2 + 2 * (g - 1) * v_n * a - (g - 1) * v_n^2)/(g^2 - 1)) - - # return SVector(f1p, f2p, f3p, f4p) - - # ## this is actually Zha and Bilgen but just put it in here for now - # # s_hat = norm(normal_direction) - # # normal_vector = normal_direction / s_hat - # # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 - # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # p_plus = 0.5 * p * (1 + v_n / a) - # phi_plus = 0.5 * (v_n + a) - - # f1p = max(v_n, zero(eltype(u))) * rho - # f2p = max(v_n, zero(eltype(u))) * rho_v1 + normal_direction[1] * p_plus - # f3p = max(v_n, zero(eltype(u))) * rho_v2 + normal_direction[2] * p_plus - # f4p = max(v_n, zero(eltype(u))) * rho_e + p * phi_plus - - # return SVector(f1p, f2p, f3p, f4p) - - # ## alternative scaling from the paper of Drikakis and Tsangris - # # s_hat = norm(normal_direction) - # # normal_vector = normal_direction / s_hat - # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # f1p = 0.25 * rho * a * (v_n / a + 1)^2 - # f2p = f1p * (v1 - normal_direction[1] / equations.gamma * (v_n - 2 * a)) - # f3p = f1p * (v2 - normal_direction[2] / equations.gamma * (v_n - 2 * a)) - # f4p = f1p * H - - # return SVector(f1p, f2p, f3p, f4p)# * s_hat + return SVector(f1p, f2p, f3p, f4p) end @inline function splitting_vanleer_haenel(u, ::Val{:minus}, @@ -1221,72 +1096,23 @@ end a = sqrt(equations.gamma * p / rho) H = (rho_e + p) / rho - ## equivalent to the rotate u, calc flux, and then backrotate strategy - # s_hat = norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - M = v_n / a - p_minus = 0.5 * (1 - equations.gamma * M) * p # HOPE style - # p_minus = 0.5 * (1 - M) * p # Liou-Steffan style - # p_minus = 0.25 * (M - 1)^2 * (2 + M) * p # van Leer style + p_minus = 0.5 * (1 - equations.gamma * M) * p f1m = -0.25 * rho * a * (M - 1)^2 f2m = f1m * v1 + normal_direction[1] * p_minus f3m = f1m * v2 + normal_direction[2] * p_minus f4m = f1m * H - return SVector(f1m, f2m, f3m, f4m)# * s_hat - - # ## yet another form taken from Palmer (https://doi.org/10.2514/3.26178) - # s_hat = 1.0 # norm(normal_direction) - # normal_vector = normal_direction / s_hat - # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 - # # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # g = equations.gamma - # f1m = -s_hat * 0.25 * rho * (v_n - a)^2 / a - # f2m = f1m * (v1 + normal_vector[1] / g * (-2 * a - v_n)) - # f3m = f1m * (v2 + normal_vector[2] / g * (-2 * a - v_n)) - # f4m = f1m * (0.5 * (v1^2 + v2^2) + (2 * a^2 - 2 * (g - 1) * v_n * a - (g - 1) * v_n^2)/(g^2 - 1)) - - # return SVector(f1m, f2m, f3m, f4m) - - # ## this is actually Zha and Bilgen but just put it in here for now - # # s_hat = norm(normal_direction) - # # normal_vector = normal_direction / s_hat - # # v_n = normal_vector[1] * v1 + normal_vector[2] * v2 - # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # p_minus = 0.5 * p * (1 - v_n / a) - # phi_minus = 0.5 * (v_n - a) - - # f1m = min(v_n, zero(eltype(u))) * rho - # f2m = min(v_n, zero(eltype(u))) * rho_v1 + normal_direction[1] * p_minus - # f3m = min(v_n, zero(eltype(u))) * rho_v2 + normal_direction[2] * p_minus - # f4m = min(v_n, zero(eltype(u))) * rho_e + p * phi_minus - - # return SVector(f1m, f2m, f3m, f4m) - - # ## alternative scaling from the paper of Drikakis and Tsangris - # # s_hat = norm(normal_direction) - # # normal_vector = normal_direction / s_hat - # v_n = normal_direction[1] * v1 + normal_direction[2] * v2 - - # f1m = -0.25 * rho * a * (v_n / a - 1)^2 - # f2m = f1m * (v1 - normal_direction[1] / equations.gamma * (v_n + 2 * a)) - # f3m = f1m * (v2 - normal_direction[2] / equations.gamma * (v_n + 2 * a)) - # f4m = f1m * H - - # return SVector(f1m, f2m, f3m, f4m)# * s_hat + return SVector(f1m, f2m, f3m, f4m) end """ - splitting_lax_friedrichs(u, orientation::Integer, + splitting_lax_friedrichs(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) splitting_lax_friedrichs(u, which::Union{Val{:minus}, Val{:plus}} - orientation::Integer, + orientation_or_normal_direction, equations::CompressibleEulerEquations2D) Naive local Lax-Friedrichs style flux splitting of the form `f⁺ = 0.5 (f + λ u)` From f1cabf83d36f9a03168cbf8a0ee63b0ea3f83919 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 4 Mar 2024 14:56:00 +0100 Subject: [PATCH 05/17] cleanup and remove unnecessary containers --- src/equations/numerical_fluxes.jl | 7 +- .../dgsem_unstructured/containers_2d.jl | 101 +----------------- .../fdsbp_unstructured/containers_2d.jl | 83 ++------------ 3 files changed, 14 insertions(+), 177 deletions(-) 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 From a8cf6916357dc5baa2fbefcb721d5f2de0e6c096 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 4 Mar 2024 15:52:21 +0100 Subject: [PATCH 06/17] add tests for the new solver --- .../elixir_euler_free_stream_upwind.jl | 82 ++++++++++++++++++ .../elixir_euler_source_terms_upwind.jl | 83 +++++++++++++++++++ test/test_unstructured_2d.jl | 70 ++++++++++++++++ 3 files changed, 235 insertions(+) create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl new file mode 100644 index 00000000000..351f1dc7df1 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl @@ -0,0 +1,82 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_constant + +# Boundary conditions for free-stream preservation test +boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) + +boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream, + :cone1 => boundary_condition_free_stream, + :cone2 => boundary_condition_free_stream, + :iceCream => boundary_condition_free_stream) + +############################################################################### +# Get the Upwind FDSBP approximation space + +# Note, one must set `xmin=-1` and `xmax=1` due to the reuse +# of interpolation routines from `calc_node_coordinates!` to create +# the physical coordinates in the mappings. +D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, + derivative_order=1, + accuracy_order=8, + xmin=-1.0, xmax=1.0, + N=17) + +flux_splitting = splitting_vanleer_haenel +solver = FDSBP(D_upw, + surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)), + volume_integral = VolumeIntegralUpwind(flux_splitting)) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) + +# Mesh with second order boundary polynomials requires an upwind SBP operator +# with (at least) 4th order boundary closure to guarantee the approximation is free-stream +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/ec9a345f09199ebe471d35d5c1e4e08f/raw/15975943d8642e42f8292235314b6f1b30aa860d/mesh_inner_outer_boundaries.mesh", + joinpath(@__DIR__, "mesh_inner_outer_boundaries.mesh")) + +mesh = UnstructuredMesh2D(mesh_file) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + save_solution, + alive_callback) + +############################################################################### +# run the simulation + +# set small tolerances for the free-stream preservation test +sol = solve(ode, SSPRK43(), abstol = 1.0e-12, reltol = 1.0e-12, + save_everystep = false, callback = callbacks) + +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl new file mode 100644 index 00000000000..250dccc4cbf --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl @@ -0,0 +1,83 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_convergence_test + +source_term = source_terms_convergence_test + +boundary_condition_eoc = BoundaryConditionDirichlet(initial_condition) + +boundary_conditions = Dict( :Top => boundary_condition_eoc, + :Bottom => boundary_condition_eoc, + :Right => boundary_condition_eoc, + :Left => boundary_condition_eoc ) + +############################################################################### +# Get the Upwind FDSBP approximation space + +# Note, one must set `xmin=-1` and `xmax=1` due to the reuse +# of interpolation routines from `calc_node_coordinates!` to create +# the physical coordinates in the mappings. +D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, + derivative_order=1, + accuracy_order=4, + xmin=-1.0, xmax=1.0, + N=9) + +flux_splitting = splitting_steger_warming +solver = FDSBP(D_upw, + surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)), + volume_integral = VolumeIntegralUpwind(flux_splitting)) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) + +# Mesh with first order boundary polynomials requires an upwind SBP operator +# with (at least) 2nd order boundary closure to guarantee the approximation is free-stream +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a4f4743008bf3233957a9ea6ac7a62e0/raw/8b36cc6649153fe0a5723b200368a210a1d74eaf/mesh_refined_box.mesh", + joinpath(@__DIR__, "mesh_refined_box.mesh")) + +mesh = UnstructuredMesh2D(mesh_file) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_term, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + save_solution, + alive_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(), abstol=1.0e-6, reltol=1.0e-6, + save_everystep = false, callback = callbacks) + +summary_callback() # print the timer summary diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 87d677e1623..8a62dcaec3c 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -610,6 +610,76 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "FDSBP (upwind): elixir_euler_source_terms_upwind.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_source_terms_upwind.jl"), + l2=[4.085391175504837e-5, + 7.19179253772227e-5, + 7.191792537723135e-5, + 0.00021775241532855398], + linf=[0.0004054489124620808, + 0.0006164432358217731, + 0.0006164432358186644, + 0.001363103391379461], + tspan=(0.0, 0.05), + atol=1.0e-10) + # 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 + +@trixi_testset "FDSBP (upwind): elixir_euler_source_terms_upwind.jl with LF splitting" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_source_terms_upwind.jl"), + l2=[3.8300267071890586e-5, + 5.295846741663533e-5, + 5.295846741663526e-5, + 0.00017564759295593478], + linf=[0.00018810716496542312, + 0.0003794187430412599, + 0.0003794187430412599, + 0.0009632958510650269], + tspan=(0.0, 0.025), + flux_splitting=splitting_lax_friedrichs, + atol=1.0e-10) + # 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 + +@trixi_testset "FDSBP (upwind): elixir_euler_free_stream_upwind.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_free_stream_upwind.jl"), + l2=[3.2114065566681054e-14, + 2.132488788134846e-14, + 2.106144937311659e-14, + 8.609642264224197e-13], + linf=[3.354871935812298e-11, + 7.006478730531285e-12, + 1.148153794261475e-11, + 9.041265514042607e-10], + tspan=(0.0, 0.05), + atol=1.0e-10) + # 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 # Clean up afterwards: delete Trixi.jl output directory From a622fb69206997d9fe26b54f52757e90775794d2 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 5 Mar 2024 06:11:23 +0100 Subject: [PATCH 07/17] remove extra space --- src/equations/numerical_fluxes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index f0a91dfcbc3..e3e798381ae 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -437,7 +437,7 @@ end equations::AbstractEquations{2}) @unpack splitting = numflux f_tilde_m = splitting(u_rr, Val{:minus}(), normal_direction, equations) - f_tilde_p = splitting(u_ll, Val{:plus}() , normal_direction, equations) + f_tilde_p = splitting(u_ll, Val{:plus}(), normal_direction, equations) return f_tilde_m + f_tilde_p end From 376e84d055dbc9daa707716cf1018ca08605268d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 5 Mar 2024 06:31:35 +0100 Subject: [PATCH 08/17] run formatter --- .../elixir_euler_free_stream_upwind.jl | 8 +++---- .../elixir_euler_source_terms_upwind.jl | 18 +++++++------- src/equations/compressible_euler_2d.jl | 24 ++++++++++++------- 3 files changed, 29 insertions(+), 21 deletions(-) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl index 351f1dc7df1..4f391f81f89 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl @@ -24,10 +24,10 @@ boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream, # of interpolation routines from `calc_node_coordinates!` to create # the physical coordinates in the mappings. D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, - derivative_order=1, - accuracy_order=8, - xmin=-1.0, xmax=1.0, - N=17) + derivative_order = 1, + accuracy_order = 8, + xmin = -1.0, xmax = 1.0, + N = 17) flux_splitting = splitting_vanleer_haenel solver = FDSBP(D_upw, diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl index 250dccc4cbf..8a4cc70e764 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl @@ -13,10 +13,10 @@ source_term = source_terms_convergence_test boundary_condition_eoc = BoundaryConditionDirichlet(initial_condition) -boundary_conditions = Dict( :Top => boundary_condition_eoc, - :Bottom => boundary_condition_eoc, - :Right => boundary_condition_eoc, - :Left => boundary_condition_eoc ) +boundary_conditions = Dict(:Top => boundary_condition_eoc, + :Bottom => boundary_condition_eoc, + :Right => boundary_condition_eoc, + :Left => boundary_condition_eoc) ############################################################################### # Get the Upwind FDSBP approximation space @@ -25,10 +25,10 @@ boundary_conditions = Dict( :Top => boundary_condition_eoc, # of interpolation routines from `calc_node_coordinates!` to create # the physical coordinates in the mappings. D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, - derivative_order=1, - accuracy_order=4, - xmin=-1.0, xmax=1.0, - N=9) + derivative_order = 1, + accuracy_order = 4, + xmin = -1.0, xmax = 1.0, + N = 9) flux_splitting = splitting_steger_warming solver = FDSBP(D_upw, @@ -77,7 +77,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, SSPRK43(), abstol=1.0e-6, reltol=1.0e-6, +sol = solve(ode, SSPRK43(), abstol = 1.0e-6, reltol = 1.0e-6, save_everystep = false, callback = callbacks) summary_callback() # print the timer summary diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 69a7f777f49..616f133c737 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -714,8 +714,10 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() """ @inline function splitting_steger_warming(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) - fm = splitting_steger_warming(u, Val{:minus}(), orientation_or_normal_direction, equations) - fp = splitting_steger_warming(u, Val{:plus}(), orientation_or_normal_direction, equations) + fm = splitting_steger_warming(u, Val{:minus}(), orientation_or_normal_direction, + equations) + fp = splitting_steger_warming(u, Val{:plus}(), orientation_or_normal_direction, + equations) return fm, fp end @@ -997,8 +999,10 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() """ @inline function splitting_vanleer_haenel(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) - fm = splitting_vanleer_haenel(u, Val{:minus}(), orientation_or_normal_direction, equations) - fp = splitting_vanleer_haenel(u, Val{:plus}(), orientation_or_normal_direction, equations) + fm = splitting_vanleer_haenel(u, Val{:minus}(), orientation_or_normal_direction, + equations) + fp = splitting_vanleer_haenel(u, Val{:plus}(), orientation_or_normal_direction, + equations) return fm, fp end @@ -1129,8 +1133,10 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() """ @inline function splitting_lax_friedrichs(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) - fm = splitting_lax_friedrichs(u, Val{:minus}(), orientation_or_normal_direction, equations) - fp = splitting_lax_friedrichs(u, Val{:plus}(), orientation_or_normal_direction, equations) + fm = splitting_lax_friedrichs(u, Val{:minus}(), orientation_or_normal_direction, + equations) + fp = splitting_lax_friedrichs(u, Val{:plus}(), orientation_or_normal_direction, + equations) return fm, fp end @@ -1188,7 +1194,8 @@ end return SVector(f1m, f2m, f3m, f4m) end -@inline function splitting_lax_friedrichs(u, ::Val{:plus}, normal_direction::AbstractVector, +@inline function splitting_lax_friedrichs(u, ::Val{:plus}, + normal_direction::AbstractVector, equations::CompressibleEulerEquations2D) rho_e = last(u) rho, v1, v2, p = cons2prim(u, equations) @@ -1208,7 +1215,8 @@ end return SVector(f1p, f2p, f3p, f4p) end -@inline function splitting_lax_friedrichs(u, ::Val{:minus}, normal_direction::AbstractVector, +@inline function splitting_lax_friedrichs(u, ::Val{:minus}, + normal_direction::AbstractVector, equations::CompressibleEulerEquations2D) rho_e = last(u) rho, v1, v2, p = cons2prim(u, equations) From 40f6547a2fc227ce4c10d15da9a03790b430f0bb Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 5 Mar 2024 21:40:42 +0100 Subject: [PATCH 09/17] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- .../elixir_euler_free_stream_upwind.jl | 8 ++++++-- .../elixir_euler_source_terms_upwind.jl | 10 +++++++--- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 11 +++++++---- 3 files changed, 20 insertions(+), 9 deletions(-) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl index 4f391f81f89..2a1956f9d10 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl @@ -1,3 +1,5 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. using OrdinaryDiffEq using Trixi @@ -20,6 +22,7 @@ boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream, ############################################################################### # Get the Upwind FDSBP approximation space +# TODO: FDSBP # Note, one must set `xmin=-1` and `xmax=1` due to the reuse # of interpolation routines from `calc_node_coordinates!` to create # the physical coordinates in the mappings. @@ -37,8 +40,9 @@ solver = FDSBP(D_upw, ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) -# Mesh with second order boundary polynomials requires an upwind SBP operator -# with (at least) 4th order boundary closure to guarantee the approximation is free-stream +# Mesh with second-order boundary polynomials requires an upwind SBP operator +# with (at least) 4th order boundary closure to guarantee the approximation is +# free-stream preserving mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/ec9a345f09199ebe471d35d5c1e4e08f/raw/15975943d8642e42f8292235314b6f1b30aa860d/mesh_inner_outer_boundaries.mesh", joinpath(@__DIR__, "mesh_inner_outer_boundaries.mesh")) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl index 8a4cc70e764..3bb2c4f1cc9 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl @@ -1,3 +1,5 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. using OrdinaryDiffEq using Trixi @@ -21,6 +23,7 @@ boundary_conditions = Dict(:Top => boundary_condition_eoc, ############################################################################### # Get the Upwind FDSBP approximation space +# TODO: FDSBP # Note, one must set `xmin=-1` and `xmax=1` due to the reuse # of interpolation routines from `calc_node_coordinates!` to create # the physical coordinates in the mappings. @@ -38,15 +41,16 @@ solver = FDSBP(D_upw, ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) -# Mesh with first order boundary polynomials requires an upwind SBP operator -# with (at least) 2nd order boundary closure to guarantee the approximation is free-stream +# Mesh with first-order boundary polynomials requires an upwind SBP operator +# with (at least) 2nd order boundary closure to guarantee the approximation is +# free-stream preserving mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a4f4743008bf3233957a9ea6ac7a62e0/raw/8b36cc6649153fe0a5723b200368a210a1d74eaf/mesh_refined_box.mesh", joinpath(@__DIR__, "mesh_refined_box.mesh")) mesh = UnstructuredMesh2D(mesh_file) ############################################################################### -# create the semi discretization object +# create the semidiscretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_term, diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index 948c020a1fb..f0aca6ec690 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -125,16 +125,19 @@ function calc_volume_integral!(du, u, # the fluxes line-by-line and add them to `du` for each element. @threaded for element in eachelement(dg, cache) # f_minus_plus_element wraps the storage provided by f_minus_element and - # f_plus_element such that we can use a single plain broadcasting below. - # f_minus_element and f_plus_element are updated in broadcasting calls - # of the form `@. f_minus_plus_element = ...`. + # f_plus_element such that we can use a single assignment below. + # f_minus_element and f_plus_element are updated whenever we update + # `f_minus_plus_element[i, j] = ...` below. f_minus_plus_element = f_minus_plus_threaded[Threads.threadid()] f_minus_element = f_minus_threaded[Threads.threadid()] f_plus_element = f_plus_threaded[Threads.threadid()] u_element = view(u_vectors, :, :, element) # x direction - # @. f_minus_plus_element = splitting(u_element, 1, equations) + # We use flux vector splittings in the directions of the contravariant + # basis vectors. Thus, we do not use a broadcasting operation like + # @. f_minus_plus_element = splitting(u_element, 1, equations) + # in the Cartesian case but loop over all nodes. for j in eachnode(dg), i in eachnode(dg) # contravariant vectors computed with central D matrix Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) From b836c21248899f08992dfe8525aca03b82a1f1fe Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 08:29:39 +0100 Subject: [PATCH 10/17] add specialized calc_metric_terms function for upwind type --- .../fdsbp_unstructured/containers_2d.jl | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl index 76facd59394..f68b1e00f59 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -9,21 +9,14 @@ #! 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. For a standard -# SBP operator the matrix is passed directly, whereas for an upwind SBP operator we must -# specifically pass the central differencing matrix. +# OBS! Requires the SBP derivative matrix in order to compute metric terms. 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) - 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_metric_terms!(elements.jacobian_matrix, element, basis, + elements.node_coordinates) calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix) @@ -36,6 +29,13 @@ function init_element!(elements, element, basis::AbstractDerivativeOperator, return elements end +# Specialization to pass the central differencing matrix from an upwind SBP operator +function calc_metric_terms!(jacobian_matrix, element, + D_SBP::SummationByPartsOperators.UpwindOperators, + node_coordinates) + calc_metric_terms!(jacobian_matrix, element, D_SBP.central, node_coordinates) +end + # construct the metric terms for a FDSBP element "block". Directly use the derivative matrix # applied to the node coordinates. function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator, From b2a4a2de4668d60372a4e9dafa66d00d64d4e18c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 08:30:37 +0100 Subject: [PATCH 11/17] revert change to the surface integral function --- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index f0aca6ec690..c35772cdf18 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -187,7 +187,7 @@ end # surface contributions are added. function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, equations, surface_integral::SurfaceIntegralStrongForm, - dg::FDSBP, cache) + dg::DG, cache) inv_weight_left = inv(left_boundary_weight(dg.basis)) inv_weight_right = inv(right_boundary_weight(dg.basis)) @unpack normal_directions, surface_flux_values = cache.elements From 5f6238e5eedb1c2da9b8fa5a6d739b34e74fe9c2 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 08:37:13 +0100 Subject: [PATCH 12/17] add reference for curvilinear van Leer splitting --- src/equations/compressible_euler_2d.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 616f133c737..557ea3dcc69 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -975,7 +975,8 @@ contains a reformulation due to Hänel et al. where the energy flux uses the enthalpy. The pressure splitting is independent from the splitting of the convective terms. As such there are many pressure splittings suggested across the literature. We implement the 'p4' variant suggested by Liou and Steffen as -it proved the most robust in practice. +it proved the most robust in practice. For details on the curvilinear variant +of this flux vector splitting see Anderson et al. Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the @@ -996,6 +997,9 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() - Meng-Sing Liou and Chris J. Steffen, Jr. (1991) High-Order Polynomial Expansions (HOPE) for Flux-Vector Splitting [NASA Technical Memorandum](https://ntrs.nasa.gov/citations/19910016425) +- W. Kyle Anderson, James L. Thomas, and Bram van Leer (1986) + Comparison of Finite Volume Flux Vector Splittings for the Euler Equations + [DOI: 10.2514/3.9465](https://doi.org/10.2514/3.9465) """ @inline function splitting_vanleer_haenel(u, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) From 49c87c57a72c3f4037c23f1f9c298558a96c37e0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 09:42:02 +0100 Subject: [PATCH 13/17] new splitting_drikakis_tsangaris in Cartesian and generalized coordinates --- .../elixir_euler_source_terms_upwind.jl | 2 +- src/Trixi.jl | 3 +- src/equations/compressible_euler_2d.jl | 139 ++++++++++++++++-- 3 files changed, 127 insertions(+), 17 deletions(-) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl index 3bb2c4f1cc9..9bd2afa5749 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl @@ -33,7 +33,7 @@ D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, xmin = -1.0, xmax = 1.0, N = 9) -flux_splitting = splitting_steger_warming +flux_splitting = splitting_drikakis_tsangaris solver = FDSBP(D_upw, surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)), volume_integral = VolumeIntegralUpwind(flux_splitting)) diff --git a/src/Trixi.jl b/src/Trixi.jl index 5f8cd9cae8e..da7359999c5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -191,7 +191,8 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, FluxUpwind export splitting_steger_warming, splitting_vanleer_haenel, - splitting_coirier_vanleer, splitting_lax_friedrichs + splitting_coirier_vanleer, splitting_lax_friedrichs, + splitting_drikakis_tsangaris export initial_condition_constant, initial_condition_gauss, diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 557ea3dcc69..6457967fc8a 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -683,15 +683,15 @@ end end """ - splitting_steger_warming(u, orientation_or_normal_direction, + splitting_steger_warming(u, orientation::Integer, equations::CompressibleEulerEquations2D) splitting_steger_warming(u, which::Union{Val{:minus}, Val{:plus}} - orientation_or_normal_direction, + orientation::Integer, equations::CompressibleEulerEquations2D) Splitting of the compressible Euler flux of Steger and Warming. For -curvilinear coordinates to compute in the `normal_direction` we use -a modified version of this splitting proposed by Drikakis and Tsangris. +curvilinear coordinates use the improved Steger-Warming-type splitting +[`splitting_drikakis_tsangaris`](@ref). Returns a tuple of the fluxes "minus" (associated with waves going into the negative axis direction) and "plus" (associated with waves going into the @@ -707,12 +707,8 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() Flux Vector Splitting of the Inviscid Gasdynamic Equations With Application to Finite Difference Methods [NASA Technical Memorandum](https://ntrs.nasa.gov/api/citations/19790020779/downloads/19790020779.pdf) -- D. Drikakis and S. Tsangaris (1993) - On the solution of the compressible Navier-Stokes equations using - improved flux vector splitting methods - [DOI: 10.1016/0307-904X(93)90054-K](https://doi.org/10.1016/0307-904X(93)90054-K) """ -@inline function splitting_steger_warming(u, orientation_or_normal_direction, +@inline function splitting_steger_warming(u, orientation::Integer, equations::CompressibleEulerEquations2D) fm = splitting_steger_warming(u, Val{:minus}(), orientation_or_normal_direction, equations) @@ -817,9 +813,122 @@ end return SVector(f1m, f2m, f3m, f4m) end -@inline function splitting_steger_warming(u, ::Val{:plus}, - normal_direction::AbstractVector, - equations::CompressibleEulerEquations2D) +""" + splitting_drikakis_tsangaris(u, orientation_or_normal_direction, + equations::CompressibleEulerEquations2D) + splitting_drikakis_tsangaris(u, which::Union{Val{:minus}, Val{:plus}} + orientation_or_normal_direction, + equations::CompressibleEulerEquations2D) + +Improved variant of the Steger-Warming flux vector splitting for +generalized coordinates. This splitting also reformulates the energy +flux as in Hänel et al. to obtain conservation of the total temperature +for inviscid flows. + +Returns a tuple of the fluxes "minus" (associated with waves going into the +negative axis direction) and "plus" (associated with waves going into the +positive axis direction). If only one of the fluxes is required, use the +function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}()`. + +!!! warning "Experimental implementation (upwind SBP)" + This is an experimental feature and may change in future releases. + +## References + +- D. Drikakis and S. Tsangaris (1993) + On the solution of the compressible Navier-Stokes equations using + improved flux vector splitting methods + [DOI: 10.1016/0307-904X(93)90054-K](https://doi.org/10.1016/0307-904X(93)90054-K) +- D. Hänel, R. Schwane and G. Seider (1987) + On the accuracy of upwind schemes for the solution of the Navier-Stokes equations + [DOI: 10.2514/6.1987-1105](https://doi.org/10.2514/6.1987-1105) +""" +@inline function splitting_drikakis_tsangaris(u, orientation_or_normal_direction, + equations::CompressibleEulerEquations2D) + fm = splitting_drikakis_tsangaris(u, Val{:minus}(), orientation_or_normal_direction, + equations) + fp = splitting_drikakis_tsangaris(u, Val{:plus}(), orientation_or_normal_direction, + equations) + return fm, fp +end + +@inline function splitting_drikakis_tsangaris(u, ::Val{:plus}, orientation::Integer, + equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + a = sqrt(equations.gamma * p / rho) + H = (eho_e + p) / rho + + if orientation == 1 + lambda1 = v1 + a + lambda2 = v1 - a + + lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) + lambda2_p = positive_part(lambda2) + + rhoa_2gamma = 0.5 * rho * a / equations.gamma + f1p = 0.5 * rho * (lambda1_p + lambda2_p) + f2p = f1p * v1 + rhoa_2gamma * (lambda1_p - lambda2_p) + f3p = f1p * v2 + f4p = f1p * H + else # orientation == 2 + lambda1 = v2 + a + lambda2 = v2 - a + + lambda1_p = positive_part(lambda1) # Same as (lambda_i + abs(lambda_i)) / 2, but faster :) + lambda2_p = positive_part(lambda2) + + rhoa_2gamma = 0.5 * rho * a / equations.gamma + f1p = 0.5 * rho * (lambda1_p + lambda2_p) + f2p = f1p * v1 + f3p = f1p * v2 + rhoa_2gamma * (lambda1_p - lambda2_p) + f4p = f1p * H + end + return SVector(f1p, f2p, f3p, f4p) +end + +@inline function splitting_drikakis_tsangaris(u, ::Val{:minus}, orientation::Integer, + equations::CompressibleEulerEquations2D) + rho, rho_v1, rho_v2, rho_e = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + a = sqrt(equations.gamma * p / rho) + H = (rho_e + p) / rho + + if orientation == 1 + lambda1 = v1 + a + lambda2 = v1 - a + + lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) + lambda2_m = negative_part(lambda2) + + rhoa_2gamma = 0.5 * rho * a / equations.gamma + f1m = 0.5 * rho * (lambda1_m + lambda2_m) + f2m = f1m * v1 + rhoa_2gamma * (lambda1_m - lambda2_m) + f3m = f1m * v2 + f4m = f1m * H + else # orientation == 2 + lambda1 = v2 + a + lambda2 = v2 - a + + lambda1_m = negative_part(lambda1) # Same as (lambda_i - abs(lambda_i)) / 2, but faster :) + lambda2_m = negative_part(lambda2) + + rhoa_2gamma = 0.5 * rho * a / equations.gamma + f1m = 0.5 * rho * (lambda1_m + lambda2_m) + f2m = f1m * v1 + f3m = f1m * v2 + rhoa_2gamma * (lambda1_m - lambda2_m) + f4m = f1m * H + end + return SVector(f1m, f2m, f3m, f4m) +end + +@inline function splitting_drikakis_tsangaris(u, ::Val{:plus}, + normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho @@ -844,9 +953,9 @@ end return SVector(f1p, f2p, f3p, f4p) end -@inline function splitting_steger_warming(u, ::Val{:minus}, - normal_direction::AbstractVector, - equations::CompressibleEulerEquations2D) +@inline function splitting_drikakis_tsangaris(u, ::Val{:minus}, + normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u v1 = rho_v1 / rho v2 = rho_v2 / rho From 22fbda840f68bb73e1f56dca3866eecf718b9614 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 10:28:26 +0100 Subject: [PATCH 14/17] added test for Cartesian splitting_drikakis_tsangaris --- src/equations/compressible_euler_2d.jl | 2 +- test/test_tree_2d_fdsbp.jl | 26 ++++++++++++++++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 6457967fc8a..892eacb5d21 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -859,7 +859,7 @@ end v2 = rho_v2 / rho p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) a = sqrt(equations.gamma * p / rho) - H = (eho_e + p) / rho + H = (rho_e + p) / rho if orientation == 1 lambda1 = v1 + a diff --git a/test/test_tree_2d_fdsbp.jl b/test/test_tree_2d_fdsbp.jl index c0844ee5dba..d477cab0563 100644 --- a/test/test_tree_2d_fdsbp.jl +++ b/test/test_tree_2d_fdsbp.jl @@ -102,6 +102,32 @@ end end end + @trixi_testset "elixir_euler_convergence.jl with Drikakis-Tsangaris splitting" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence.jl"), + l2=[ + 1.708838999643608e-6, + 1.7437997854485807e-6, + 1.7437997854741082e-6, + 5.457223460116349e-6, + ], + linf=[ + 9.796504911285808e-6, + 9.614745899888533e-6, + 9.614745899444443e-6, + 4.02610718399643e-5, + ], + tspan=(0.0, 0.1), flux_splitting=splitting_drikakis_tsangaris) + + # 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 + @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"), From 59233928bb9deae8b9b0d5c7e5fe8086b53f84ef Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 10:31:35 +0100 Subject: [PATCH 15/17] run formatter --- src/equations/compressible_euler_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 892eacb5d21..4d6aa862445 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -844,7 +844,7 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() [DOI: 10.2514/6.1987-1105](https://doi.org/10.2514/6.1987-1105) """ @inline function splitting_drikakis_tsangaris(u, orientation_or_normal_direction, - equations::CompressibleEulerEquations2D) + equations::CompressibleEulerEquations2D) fm = splitting_drikakis_tsangaris(u, Val{:minus}(), orientation_or_normal_direction, equations) fp = splitting_drikakis_tsangaris(u, Val{:plus}(), orientation_or_normal_direction, From 43fd5b18236fd6fe78ae451bfc22909232704b86 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 6 Mar 2024 11:04:55 +0100 Subject: [PATCH 16/17] Update src/equations/compressible_euler_2d.jl --- src/equations/compressible_euler_2d.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 4d6aa862445..84740ba429a 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -820,8 +820,9 @@ end orientation_or_normal_direction, equations::CompressibleEulerEquations2D) -Improved variant of the Steger-Warming flux vector splitting for -generalized coordinates. This splitting also reformulates the energy +Improved variant of the Steger-Warming flux vector splitting +[`splitting_steger_warming`](@ref) for generalized coordinates. +This splitting also reformulates the energy flux as in Hänel et al. to obtain conservation of the total temperature for inviscid flows. From c11af1165130bde5aa02e60b17e644195d35a4e0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 11:28:43 +0100 Subject: [PATCH 17/17] remove orientation_or_normal from Steger-Warming --- src/equations/compressible_euler_2d.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 84740ba429a..43f15a3cfb9 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -710,10 +710,8 @@ function signature with argument `which` set to `Val{:minus}()` or `Val{:plus}() """ @inline function splitting_steger_warming(u, orientation::Integer, equations::CompressibleEulerEquations2D) - fm = splitting_steger_warming(u, Val{:minus}(), orientation_or_normal_direction, - equations) - fp = splitting_steger_warming(u, Val{:plus}(), orientation_or_normal_direction, - equations) + fm = splitting_steger_warming(u, Val{:minus}(), orientation, equations) + fp = splitting_steger_warming(u, Val{:plus}(), orientation, equations) return fm, fp end