From 8cdb93892dead6e33ae581e60fd8199aec4ee19e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 6 Mar 2024 13:24:50 +0100 Subject: [PATCH 01/11] Upwind SBP on curved meshes (#1857) * baseline implementation of the curvilinear USBP for testing * working version of curvilinear upwind solver. Needs significant cleanup of debugging statements and different variants of the rotated flux vector splittings * cleanup of the fdsbp_2d file * clean-up FVS routines in the compressible Euler file * cleanup and remove unnecessary containers * add tests for the new solver * remove extra space * run formatter * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * add specialized calc_metric_terms function for upwind type * revert change to the surface integral function * add reference for curvilinear van Leer splitting * new splitting_drikakis_tsangaris in Cartesian and generalized coordinates * added test for Cartesian splitting_drikakis_tsangaris * run formatter * Update src/equations/compressible_euler_2d.jl * remove orientation_or_normal from Steger-Warming --------- Co-authored-by: Hendrik Ranocha --- .../elixir_euler_free_stream_upwind.jl | 86 ++++++ .../elixir_euler_source_terms_upwind.jl | 87 ++++++ src/Trixi.jl | 3 +- src/equations/compressible_euler_2d.jl | 290 +++++++++++++++++- src/equations/numerical_fluxes.jl | 12 +- src/solvers/dgsem_unstructured/dg_2d.jl | 2 +- src/solvers/fdsbp_tree/fdsbp_2d.jl | 2 +- .../fdsbp_unstructured/containers_2d.jl | 10 +- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 99 +++++- test/test_tree_2d_fdsbp.jl | 26 ++ test/test_unstructured_2d.jl | 70 +++++ 11 files changed, 659 insertions(+), 28 deletions(-) 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..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 From 253c358243ed2a8b6f63c422b53cf09a68188f68 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 6 Mar 2024 13:25:19 +0100 Subject: [PATCH 02/11] set version to v0.7.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 800c7b4c0fa..53b859a422f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.1-pre" +version = "0.7.1" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From ca082c2eb1273611cf38e80c6d2dab04e8f8177f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 6 Mar 2024 13:25:37 +0100 Subject: [PATCH 03/11] set development version to v0.7.2-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 53b859a422f..e8eb7c788ce 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.1" +version = "0.7.2-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 4bf61a0fd3abb21b457ea9d9ee2c19c835d018c3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 6 Mar 2024 14:25:25 +0100 Subject: [PATCH 04/11] force new FFMPEG.jl version for tests (#1858) --- test/Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/Project.toml b/test/Project.toml index a376c2805ea..1a042dab44f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" @@ -15,6 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Aqua = "0.8" CairoMakie = "0.10" Downloads = "1" +FFMPEG = "0.4" ForwardDiff = "0.10.24" LinearAlgebra = "1" MPI = "0.20" From 3bed8285ca2bb3857050bde4fd1408d151cb760d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 7 Mar 2024 10:41:56 +0100 Subject: [PATCH 05/11] Check BCs for periodicity for periodic Tree & Structured meshes (#1860) * Check BCs for periodicity for periodic meshes * default case for periodic bcs * fmt * specialize * error ant fmt * isperiodic TreeMesh * avoid if * shorten * shorten * Make meshes non-periodic * fix rti * shorten dispatch --------- Co-authored-by: Hendrik Ranocha --- .../elixir_advection_coupled.jl | 12 ++-- ...lixir_euler_rayleigh_taylor_instability.jl | 2 +- .../elixir_euler_warm_bubble.jl | 3 +- src/meshes/tree_mesh.jl | 3 + .../semidiscretization_hyperbolic.jl | 70 +++++++++++++++++++ ...semidiscretization_hyperbolic_parabolic.jl | 2 + 6 files changed, 86 insertions(+), 6 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_advection_coupled.jl b/examples/structured_2d_dgsem/elixir_advection_coupled.jl index 43b68f21b03..0002bb8d374 100644 --- a/examples/structured_2d_dgsem/elixir_advection_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_advection_coupled.jl @@ -53,7 +53,8 @@ cells_per_dimension = (8, 8) coordinates_min1 = (-1.0, 0.0) # minimum coordinates (min(x), min(y)) coordinates_max1 = (0.0, 1.0) # maximum coordinates (max(x), max(y)) -mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1) +mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1, + periodicity = false) # Define the coupling functions coupling_function12 = (x, u, equations_other, equations_own) -> u @@ -84,7 +85,8 @@ semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_converg coordinates_min2 = (0.0, 0.0) # minimum coordinates (min(x), min(y)) coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y)) -mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2) +mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2, + periodicity = false) # Define the coupling functions coupling_function21 = (x, u, equations_other, equations_own) -> u @@ -115,7 +117,8 @@ semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_converg coordinates_min3 = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) coordinates_max3 = (0.0, 0.0) # maximum coordinates (max(x), max(y)) -mesh3 = StructuredMesh(cells_per_dimension, coordinates_min3, coordinates_max3) +mesh3 = StructuredMesh(cells_per_dimension, coordinates_min3, coordinates_max3, + periodicity = false) # Define the coupling functions coupling_function34 = (x, u, equations_other, equations_own) -> u @@ -146,7 +149,8 @@ semi3 = SemidiscretizationHyperbolic(mesh3, equations, initial_condition_converg coordinates_min4 = (0.0, -1.0) # minimum coordinates (min(x), min(y)) coordinates_max4 = (1.0, 0.0) # maximum coordinates (max(x), max(y)) -mesh4 = StructuredMesh(cells_per_dimension, coordinates_min4, coordinates_max4) +mesh4 = StructuredMesh(cells_per_dimension, coordinates_min4, coordinates_max4, + periodicity = false) # Define the coupling functions coupling_function43 = (x, u, equations_other, equations_own) -> u diff --git a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl index 6c254e8bd8b..dd0cc339b20 100644 --- a/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl +++ b/examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl @@ -69,7 +69,7 @@ mapping(xi, eta) = SVector(0.25 * 0.5 * (1.0 + xi), 0.5 * (1.0 + eta)) num_elements_per_dimension = 32 cells_per_dimension = (num_elements_per_dimension, num_elements_per_dimension * 4) -mesh = StructuredMesh(cells_per_dimension, mapping) +mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false) initial_condition = initial_condition_rayleigh_taylor_instability boundary_conditions = (x_neg = boundary_condition_slip_wall, diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index 05c09d57530..38b9386e94e 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -100,7 +100,8 @@ coordinates_min = (0.0, 0.0) coordinates_max = (20_000.0, 10_000.0) cells_per_dimension = (64, 32) -mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, false)) semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, source_terms = warm_bubble_setup, diff --git a/src/meshes/tree_mesh.jl b/src/meshes/tree_mesh.jl index 05699d17d16..1092fc54cc1 100644 --- a/src/meshes/tree_mesh.jl +++ b/src/meshes/tree_mesh.jl @@ -228,5 +228,8 @@ function total_volume(mesh::TreeMesh) return mesh.tree.length_level_0^ndims(mesh) end +isperiodic(mesh::TreeMesh) = isperiodic(mesh.tree) +isperiodic(mesh::TreeMesh, dimension) = isperiodic(mesh.tree, dimension) + include("parallel_tree_mesh.jl") end # @muladd diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 7ebd758de37..f61378a7dca 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -72,6 +72,8 @@ function SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache) + check_periodicity_mesh_boundary_conditions(mesh, _boundary_conditions) + SemidiscretizationHyperbolic{typeof(mesh), typeof(equations), typeof(initial_condition), typeof(_boundary_conditions), typeof(source_terms), @@ -210,6 +212,74 @@ function digest_boundary_conditions(boundary_conditions::AbstractArray, mesh, so throw(ArgumentError("Please use a (named) tuple instead of an (abstract) array to supply multiple boundary conditions (to improve performance).")) end +# No checks for these meshes yet available +function check_periodicity_mesh_boundary_conditions(mesh::Union{P4estMesh, + UnstructuredMesh2D, + T8codeMesh, + DGMultiMesh}, + boundary_conditions) +end + +# No actions needed for periodic boundary conditions +function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh, + StructuredMesh}, + boundary_conditions::BoundaryConditionPeriodic) +end + +function check_periodicity_mesh_boundary_conditions_x(mesh, x_neg, x_pos) + if isperiodic(mesh, 1) && + (x_neg != BoundaryConditionPeriodic() || + x_pos != BoundaryConditionPeriodic()) + @error "For periodic mesh non-periodic boundary conditions in x-direction are supplied." + end +end + +function check_periodicity_mesh_boundary_conditions_y(mesh, y_neg, y_pos) + if isperiodic(mesh, 2) && + (y_neg != BoundaryConditionPeriodic() || + y_pos != BoundaryConditionPeriodic()) + @error "For periodic mesh non-periodic boundary conditions in y-direction are supplied." + end +end + +function check_periodicity_mesh_boundary_conditions_z(mesh, z_neg, z_pos) + if isperiodic(mesh, 3) && + (z_neg != BoundaryConditionPeriodic() || + z_pos != BoundaryConditionPeriodic()) + @error "For periodic mesh non-periodic boundary conditions in z-direction are supplied." + end +end + +function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{1}, + StructuredMesh{1}}, + boundary_conditions::Union{NamedTuple, + Tuple}) + check_periodicity_mesh_boundary_conditions_x(mesh, boundary_conditions[1], + boundary_conditions[2]) +end + +function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{2}, + StructuredMesh{2}}, + boundary_conditions::Union{NamedTuple, + Tuple}) + check_periodicity_mesh_boundary_conditions_x(mesh, boundary_conditions[1], + boundary_conditions[2]) + check_periodicity_mesh_boundary_conditions_y(mesh, boundary_conditions[3], + boundary_conditions[4]) +end + +function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{3}, + StructuredMesh{3}}, + boundary_conditions::Union{NamedTuple, + Tuple}) + check_periodicity_mesh_boundary_conditions_x(mesh, boundary_conditions[1], + boundary_conditions[2]) + check_periodicity_mesh_boundary_conditions_y(mesh, boundary_conditions[3], + boundary_conditions[4]) + check_periodicity_mesh_boundary_conditions_z(mesh, boundary_conditions[5], + boundary_conditions[6]) +end + function Base.show(io::IO, semi::SemidiscretizationHyperbolic) @nospecialize semi # reduce precompilation time diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 0f44941390a..57724374acb 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -136,6 +136,8 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo _boundary_conditions_parabolic = digest_boundary_conditions(boundary_conditions_parabolic, mesh, solver, cache) + check_periodicity_mesh_boundary_conditions(mesh, _boundary_conditions) + cache_parabolic = (; create_cache_parabolic(mesh, equations, equations_parabolic, solver, solver_parabolic, RealT, From b8b34ea21cfa22f93daa131c7ef77d8ed8f3f789 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 7 Mar 2024 10:45:17 +0100 Subject: [PATCH 06/11] set version to v0.7.2 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e8eb7c788ce..9a6e8e0d9a2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.2-pre" +version = "0.7.2" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From b6bc7c8cae92cf6877fca90c5ef100771a7781a8 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 7 Mar 2024 10:45:30 +0100 Subject: [PATCH 07/11] set development version to v0.7.3-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9a6e8e0d9a2..6b44af4a3fa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.2" +version = "0.7.3-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From b11bc10d8ac68798ae2e53f9db3d5c4f824af5fa Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Thu, 7 Mar 2024 13:04:19 +0100 Subject: [PATCH 08/11] add warning (#1863) --- src/solvers/dgsem/basis_lobatto_legendre.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index 6a92fd1c066..9e21b88dfa1 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -12,6 +12,7 @@ Create a nodal Lobatto-Legendre basis for polynomials of degree `polydeg`. For the special case `polydeg=0` the DG method reduces to a finite volume method. Therefore, this function sets the center point of the cell as single node. +This exceptional case is currently only supported for TreeMesh! """ struct LobattoLegendreBasis{RealT <: Real, NNODES, VectorT <: AbstractVector{RealT}, From c4bf3df8a4d4e3920b0f774960d129a4af6ee287 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Thu, 7 Mar 2024 13:06:20 +0100 Subject: [PATCH 09/11] Update parallelization.md (#1864) MPI support in T8code came with #1803 --- docs/src/parallelization.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/docs/src/parallelization.md b/docs/src/parallelization.md index f599eb5fafe..2114f30fb87 100644 --- a/docs/src/parallelization.md +++ b/docs/src/parallelization.md @@ -69,8 +69,7 @@ installations. Follow the steps described [here](https://github.com/DLR-AMR/T8co [here](https://github.com/trixi-framework/P4est.jl/blob/main/README.md#installation) for the configuration. The paths that point to `libp4est.so` (and potentially to `libsc.so`) need to be the same for P4est.jl and T8code.jl. This could, e.g., be `libp4est.so` that usually can be found -in `lib/` or `local/lib/` in the installation directory of `t8code`. Note that the `T8codeMesh`, however, -does not support MPI yet. +in `lib/` or `local/lib/` in the installation directory of `t8code`. The preferences for [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) always need to be set, even if you do not want to use `HDF5` from Trixi.jl, see also [issue #1079 in HDF5.jl](https://github.com/JuliaIO/HDF5.jl/issues/1079). To set the preferences for HDF5.jl, follow the instructions described From 1ca37cf2271806d203a832d3d99bfa2c14d3226f Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Fri, 8 Mar 2024 07:48:24 +0100 Subject: [PATCH 10/11] set capacity also when using MPI (#1862) Co-authored-by: Hendrik Ranocha --- src/meshes/mesh_io.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 337e33e6969..28e6efa8c57 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -74,6 +74,7 @@ function save_mesh_file(mesh::TreeMesh, output_directory, timestep, attributes(file)["mesh_type"] = get_name(mesh) attributes(file)["ndims"] = ndims(mesh) attributes(file)["n_cells"] = n_cells + attributes(file)["capacity"] = mesh.tree.capacity attributes(file)["n_leaf_cells"] = count_leaf_cells(mesh.tree) attributes(file)["minimum_level"] = minimum_level(mesh.tree) attributes(file)["maximum_level"] = maximum_level(mesh.tree) From f235619a49dbc8bd7d84a77269558a64b21a155f Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 8 Mar 2024 12:20:45 +0100 Subject: [PATCH 11/11] Mention hyphen/dash caveat for boundary symbols (#1866) * Mention hyphen/dash caveat for boundary symbols * typo * elaborate --- docs/src/meshes/p4est_mesh.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/meshes/p4est_mesh.md b/docs/src/meshes/p4est_mesh.md index 3b35ffcad6f..a14551b3f46 100644 --- a/docs/src/meshes/p4est_mesh.md +++ b/docs/src/meshes/p4est_mesh.md @@ -256,6 +256,8 @@ By doing so, only nodesets with a label present in `boundary_symbols` are treate Other nodesets that could be used for diagnostics are not treated as external boundaries. Note that there is a leading colon `:` compared to the label in the `.inp` mesh file. This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols). +**Important**: In Julia, a symbol _cannot_ contain a hyphen/dash `-`, i.e., `:BC-1` is _not_ a valid symbol. +Keep this in mind when importing boundaries, you might have to convert hyphens/dashes `-` to underscores `_` in the `.inp` mesh file, i.e., `BC_1` instead of `BC-1`. A 2D example for this mesh, which is read-in for an unstructured mesh file created with `gmsh`, is presented in `examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl`.