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..2a1956f9d10 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl @@ -0,0 +1,86 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. + +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 + +# 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. +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 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")) + +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..9bd2afa5749 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl @@ -0,0 +1,87 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. + +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 + +# 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. +D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, + derivative_order = 1, + accuracy_order = 4, + xmin = -1.0, xmax = 1.0, + N = 9) + +flux_splitting = splitting_drikakis_tsangaris +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 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 semidiscretization 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/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 f5a632723cf..43f15a3cfb9 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -689,7 +689,9 @@ end orientation::Integer, 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 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 @@ -809,6 +811,174 @@ end return SVector(f1m, f2m, f3m, f4m) end +""" + 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 +[`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. + +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 = (rho_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 + 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 + + 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) + + rhoa_2gamma = 0.5 * rho * a / equations.gamma + 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) +end + +@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 + 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 + + 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) + + rhoa_2gamma = 0.5 * rho * a / equations.gamma + 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) +end + """ FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) @@ -902,10 +1072,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 @@ -913,7 +1083,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 @@ -934,11 +1105,16 @@ 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::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 @@ -1002,11 +1178,57 @@ 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 + + v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + M = v_n / a + 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) +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 + + v_n = normal_direction[1] * v1 + normal_direction[2] * v2 + M = v_n / a + 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) +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)` @@ -1021,10 +1243,12 @@ 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 @@ -1082,6 +1306,48 @@ 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 6794c71a32b..e3e798381ae 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,5 +432,14 @@ end return fm + fp end +@inline function (numflux::FluxUpwind)(u_ll, u_rr, + normal_direction::AbstractVector, + 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) + 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/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_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..f68b1e00f59 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -9,7 +9,7 @@ #! 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. function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves) calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), @@ -29,9 +29,15 @@ 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. -# 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) diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index b459f4c42cc..c35772cdf18 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. @@ -86,6 +84,91 @@ 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 + @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 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 + # 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) + f_minus_plus_element[i, j] = splitting(u_element[i, j], Ja1, equations) + end + + for j in eachnode(dg) + 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 + for j in eachnode(dg), i in eachnode(dg) + # 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) + end + + for i in eachnode(dg) + 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 + # Note! The local side numbering for the unstructured quadrilateral element implementation differs # from the structured TreeMesh or StructuredMesh local side numbering: # @@ -114,8 +197,7 @@ 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_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), @@ -124,8 +206,7 @@ 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) multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), @@ -134,8 +215,7 @@ 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) multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), @@ -144,8 +224,7 @@ 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) multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), 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"), 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