From 61c33b0af6a7a49ed11258e8f230b471b05c6ed8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 31 Oct 2023 09:14:11 +0100 Subject: [PATCH 1/5] Implement subcell limiting for non-conservative systems (#1670) Co-authored-by: Hendrik Ranocha Co-authored-by: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Co-authored-by: Michael Schlottke-Lakemper --- .../elixir_mhd_shockcapturing_subcell.jl | 108 ++++++ src/Trixi.jl | 2 +- .../subcell_limiter_idp_correction_2d.jl | 10 +- src/equations/equations.jl | 26 ++ src/equations/ideal_glm_mhd_2d.jl | 190 ++++++++++ src/solvers/dgsem_tree/containers_2d.jl | 84 +++-- .../dgsem_tree/dg_2d_subcell_limiters.jl | 343 ++++++++++++++++-- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 10 +- src/time_integration/methods_SSP.jl | 5 + test/test_tree_2d_mhd.jl | 31 ++ 10 files changed, 733 insertions(+), 76 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl new file mode 100644 index 00000000000..f40da6676c2 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -0,0 +1,108 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible ideal GLM-MHD equations + +equations = IdealGlmMhdEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) + +An MHD blast wave modified from: +- Dominik Derigs, Gregor J. Gassner, Stefanie Walch & Andrew R. Winters (2018) + Entropy Stable Finite Volume Approximations for Ideal Magnetohydrodynamics + [doi: 10.1365/s13291-018-0178-9](https://doi.org/10.1365/s13291-018-0178-9) +This setup needs a positivity limiter for the density. +""" +function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) + # setup taken from Derigs et al. DMV article (2018) + # domain must be [-0.5, 0.5] x [-0.5, 0.5], γ = 1.4 + r = sqrt(x[1]^2 + x[2]^2) + + pmax = 10.0 + pmin = 1.0 + rhomax = 1.0 + rhomin = 0.01 + if r <= 0.09 + p = pmax + rho = rhomax + elseif r >= 0.1 + p = pmin + rho = rhomin + else + p = pmin + (0.1 - r) * (pmax - pmin) / 0.01 + rho = rhomin + (0.1 - r) * (rhomax - rhomin) / 0.01 + end + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + B1 = 1.0/sqrt(4.0*pi) + B2 = 0.0 + B3 = 0.0 + psi = 0.0 + return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) +end +initial_condition = initial_condition_blast_wave + +surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell_local_symmetric) +volume_flux = (flux_derigs_etal, flux_nonconservative_powell_local_symmetric) +basis = LobattoLegendreBasis(3) + +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons=[1], + positivity_correction_factor=0.5) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg=volume_flux, + volume_flux_fv=surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-0.5, -0.5) +coordinates_max = ( 0.5, 0.5) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=10_000) + + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.1) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +cfl = 0.5 +stepsize_callback = StepsizeCallback(cfl=cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks); + dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep=false, callback=callbacks); +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index 5cb3cf0a9fe..97d518d5b78 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -162,7 +162,7 @@ export GradientVariablesPrimitive, GradientVariablesEntropy export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, - flux_nonconservative_powell, + flux_nonconservative_powell, flux_nonconservative_powell_local_symmetric, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index f6b91444578..6f1723e2a98 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -7,7 +7,7 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache) @unpack inverse_weights = dg.basis - @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes + @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients @threaded for element in eachelement(dg, cache) @@ -17,16 +17,16 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache) for j in eachnode(dg), i in eachnode(dg) # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1} alpha_flux1 = (1 - alpha1[i, j, element]) * - get_node_vars(antidiffusive_flux1, equations, dg, i, j, + get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, element) alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) * - get_node_vars(antidiffusive_flux1, equations, dg, i + 1, + get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j, element) alpha_flux2 = (1 - alpha2[i, j, element]) * - get_node_vars(antidiffusive_flux2, equations, dg, i, j, + get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, element) alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) * - get_node_vars(antidiffusive_flux2, equations, dg, i, + get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1, element) for v in eachvariable(equations) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 3142dcc2765..0e77b92e045 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -208,6 +208,24 @@ struct BoundaryConditionNeumann{B} boundary_normal_flux_function::B end +""" + NonConservativeLocal() + +Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". +When the argument `nonconservative_type` is of type `NonConservativeLocal`, +the function returns the local part of the non-conservative term. +""" +struct NonConservativeLocal end + +""" + NonConservativeSymmetric() + +Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". +When the argument `nonconservative_type` is of type `NonConservativeSymmetric`, +the function returns the symmetric part of the non-conservative term. +""" +struct NonConservativeSymmetric end + # set sensible default values that may be overwritten by specific equations """ have_nonconservative_terms(equations) @@ -220,6 +238,14 @@ example of equations with nonconservative terms. The return value will be `True()` or `False()` to allow dispatching on the return type. """ have_nonconservative_terms(::AbstractEquations) = False() +""" + n_nonconservative_terms(equations) + +Number of nonconservative terms in the form local * symmetric for a particular equation. +This function needs to be specialized only if equations with nonconservative terms are +combined with certain solvers (e.g., subcell limiting). +""" +function n_nonconservative_terms end have_constant_speed(::AbstractEquations) = False() default_analysis_errors(::AbstractEquations) = (:l2_error, :linf_error) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 8fef1ee22c9..e8de0cedde1 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -29,6 +29,8 @@ function IdealGlmMhdEquations2D(gamma; initial_c_h = convert(typeof(gamma), NaN) end have_nonconservative_terms(::IdealGlmMhdEquations2D) = True() +n_nonconservative_terms(::IdealGlmMhdEquations2D) = 2 + function varnames(::typeof(cons2cons), ::IdealGlmMhdEquations2D) ("rho", "rho_v1", "rho_v2", "rho_v3", "rho_e", "B1", "B2", "B3", "psi") end @@ -279,6 +281,194 @@ end return f end +""" + flux_nonconservative_powell_local_symmetric(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdEquations2D) + +Non-symmetric two-point flux discretizing the nonconservative (source) term of +Powell and the Galilean nonconservative term associated with the GLM multiplier +of the [`IdealGlmMhdEquations2D`](@ref). + +This implementation uses a non-conservative term that can be written as the product +of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm +et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different +results on non-conforming meshes(!). + +The two other flux functions with the same name return either the local +or symmetric portion of the non-conservative flux based on the type of the +nonconservative_type argument, employing multiple dispatch. They are used to +compute the subcell fluxes in dg_2d_subcell_limiters.jl. + +## References +- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts + Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +""" +@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdEquations2D) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + psi_avg = (psi_ll + psi_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + if orientation == 1 + B1_avg = (B1_ll + B1_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B1_ll * B1_avg, + B2_ll * B1_avg, + B3_ll * B1_avg, + v_dot_B_ll * B1_avg + v1_ll * psi_ll * psi_avg, + v1_ll * B1_avg, + v2_ll * B1_avg, + v3_ll * B1_avg, + v1_ll * psi_avg) + else # orientation == 2 + B2_avg = (B2_ll + B2_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B1_ll * B2_avg, + B2_ll * B2_avg, + B3_ll * B2_avg, + v_dot_B_ll * B2_avg + v2_ll * psi_ll * psi_avg, + v1_ll * B2_avg, + v2_ll * B2_avg, + v3_ll * B2_avg, + v2_ll * psi_avg) + end + + return f +end + +""" + flux_nonconservative_powell_local_symmetric(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeLocal, + nonconservative_term::Integer) + +Local part of the Powell and GLM non-conservative terms. Needed for the calculation of +the non-conservative staggered "fluxes" for subcell limiting. See, e.g., +- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts + Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +This function is used to compute the subcell fluxes in dg_2d_subcell_limiters.jl. +""" +@inline function flux_nonconservative_powell_local_symmetric(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeLocal, + nonconservative_term::Integer) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + + if nonconservative_term == 1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + f = SVector(0, + B1_ll, + B2_ll, + B3_ll, + v_dot_B_ll, + v1_ll, + v2_ll, + v3_ll, + 0) + else #nonconservative_term ==2 + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + if orientation == 1 + v1_ll = rho_v1_ll / rho_ll + f = SVector(0, + 0, + 0, + 0, + v1_ll * psi_ll, + 0, + 0, + 0, + v1_ll) + else #orientation == 2 + v2_ll = rho_v2_ll / rho_ll + f = SVector(0, + 0, + 0, + 0, + v2_ll * psi_ll, + 0, + 0, + 0, + v2_ll) + end + end + return f +end + +""" + flux_nonconservative_powell_local_symmetric(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeSymmetric, + nonconservative_term::Integer) + +Symmetric part of the Powell and GLM non-conservative terms. Needed for the calculation of +the non-conservative staggered "fluxes" for subcell limiting. See, e.g., +- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts + Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +This function is used to compute the subcell fluxes in dg_2d_subcell_limiters.jl. +""" +@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeSymmetric, + nonconservative_term::Integer) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr + + if nonconservative_term == 1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + if orientation == 1 + B1_avg = (B1_ll + B1_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + 0) + else # orientation == 2 + B2_avg = (B2_ll + B2_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + 0) + end + else #nonconservative_term == 2 + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + psi_avg = (psi_ll + psi_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + 0, + 0, + 0, + psi_avg, + 0, + 0, + 0, + psi_avg) + end + + return f +end + """ flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations2D) diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 9e9fe88c15b..4bfbddead9a 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -1266,11 +1266,15 @@ end # | # (i, j-1) mutable struct ContainerAntidiffusiveFlux2D{uEltype <: Real} - antidiffusive_flux1::Array{uEltype, 4} # [variables, i, j, elements] - antidiffusive_flux2::Array{uEltype, 4} # [variables, i, j, elements] + antidiffusive_flux1_L::Array{uEltype, 4} # [variables, i, j, elements] + antidiffusive_flux1_R::Array{uEltype, 4} # [variables, i, j, elements] + antidiffusive_flux2_L::Array{uEltype, 4} # [variables, i, j, elements] + antidiffusive_flux2_R::Array{uEltype, 4} # [variables, i, j, elements] # internal `resize!`able storage - _antidiffusive_flux1::Vector{uEltype} - _antidiffusive_flux2::Vector{uEltype} + _antidiffusive_flux1_L::Vector{uEltype} + _antidiffusive_flux1_R::Vector{uEltype} + _antidiffusive_flux2_L::Vector{uEltype} + _antidiffusive_flux2_R::Vector{uEltype} end function ContainerAntidiffusiveFlux2D{uEltype}(capacity::Integer, n_variables, @@ -1278,24 +1282,36 @@ function ContainerAntidiffusiveFlux2D{uEltype}(capacity::Integer, n_variables, nan_uEltype = convert(uEltype, NaN) # Initialize fields with defaults - _antidiffusive_flux1 = fill(nan_uEltype, - n_variables * (n_nodes + 1) * n_nodes * capacity) - antidiffusive_flux1 = unsafe_wrap(Array, pointer(_antidiffusive_flux1), - (n_variables, n_nodes + 1, n_nodes, capacity)) - - _antidiffusive_flux2 = fill(nan_uEltype, - n_variables * n_nodes * (n_nodes + 1) * capacity) - antidiffusive_flux2 = unsafe_wrap(Array, pointer(_antidiffusive_flux2), - (n_variables, n_nodes, n_nodes + 1, capacity)) - - return ContainerAntidiffusiveFlux2D{uEltype}(antidiffusive_flux1, - antidiffusive_flux2, - _antidiffusive_flux1, - _antidiffusive_flux2) + _antidiffusive_flux1_L = fill(nan_uEltype, + n_variables * (n_nodes + 1) * n_nodes * capacity) + antidiffusive_flux1_L = unsafe_wrap(Array, pointer(_antidiffusive_flux1_L), + (n_variables, n_nodes + 1, n_nodes, capacity)) + _antidiffusive_flux1_R = fill(nan_uEltype, + n_variables * (n_nodes + 1) * n_nodes * capacity) + antidiffusive_flux1_R = unsafe_wrap(Array, pointer(_antidiffusive_flux1_R), + (n_variables, n_nodes + 1, n_nodes, capacity)) + + _antidiffusive_flux2_L = fill(nan_uEltype, + n_variables * n_nodes * (n_nodes + 1) * capacity) + antidiffusive_flux2_L = unsafe_wrap(Array, pointer(_antidiffusive_flux2_L), + (n_variables, n_nodes, n_nodes + 1, capacity)) + _antidiffusive_flux2_R = fill(nan_uEltype, + n_variables * n_nodes * (n_nodes + 1) * capacity) + antidiffusive_flux2_R = unsafe_wrap(Array, pointer(_antidiffusive_flux2_R), + (n_variables, n_nodes, n_nodes + 1, capacity)) + + return ContainerAntidiffusiveFlux2D{uEltype}(antidiffusive_flux1_L, + antidiffusive_flux1_R, + antidiffusive_flux2_L, + antidiffusive_flux2_R, + _antidiffusive_flux1_L, + _antidiffusive_flux1_R, + _antidiffusive_flux2_L, + _antidiffusive_flux2_R) end -nvariables(fluxes::ContainerAntidiffusiveFlux2D) = size(fluxes.antidiffusive_flux1, 1) -nnodes(fluxes::ContainerAntidiffusiveFlux2D) = size(fluxes.antidiffusive_flux1, 3) +nvariables(fluxes::ContainerAntidiffusiveFlux2D) = size(fluxes.antidiffusive_flux1_L, 1) +nnodes(fluxes::ContainerAntidiffusiveFlux2D) = size(fluxes.antidiffusive_flux1_L, 3) # Only one-dimensional `Array`s are `resize!`able in Julia. # Hence, we use `Vector`s as internal storage and `resize!` @@ -1306,16 +1322,24 @@ function Base.resize!(fluxes::ContainerAntidiffusiveFlux2D, capacity) n_nodes = nnodes(fluxes) n_variables = nvariables(fluxes) - @unpack _antidiffusive_flux1, _antidiffusive_flux2 = fluxes - - resize!(_antidiffusive_flux1, n_variables * (n_nodes + 1) * n_nodes * capacity) - fluxes.antidiffusive_flux1 = unsafe_wrap(Array, pointer(_antidiffusive_flux1), - (n_variables, n_nodes + 1, n_nodes, - capacity)) - resize!(_antidiffusive_flux2, n_variables * n_nodes * (n_nodes + 1) * capacity) - fluxes.antidiffusive_flux2 = unsafe_wrap(Array, pointer(_antidiffusive_flux2), - (n_variables, n_nodes, n_nodes + 1, - capacity)) + @unpack _antidiffusive_flux1_L, _antidiffusive_flux2_L, _antidiffusive_flux1_R, _antidiffusive_flux2_R = fluxes + + resize!(_antidiffusive_flux1_L, n_variables * (n_nodes + 1) * n_nodes * capacity) + fluxes.antidiffusive_flux1_L = unsafe_wrap(Array, pointer(_antidiffusive_flux1_L), + (n_variables, n_nodes + 1, n_nodes, + capacity)) + resize!(_antidiffusive_flux1_R, n_variables * (n_nodes + 1) * n_nodes * capacity) + fluxes.antidiffusive_flux1_R = unsafe_wrap(Array, pointer(_antidiffusive_flux1_R), + (n_variables, n_nodes + 1, n_nodes, + capacity)) + resize!(_antidiffusive_flux2_L, n_variables * n_nodes * (n_nodes + 1) * capacity) + fluxes.antidiffusive_flux2_L = unsafe_wrap(Array, pointer(_antidiffusive_flux2_L), + (n_variables, n_nodes, n_nodes + 1, + capacity)) + resize!(_antidiffusive_flux2_R, n_variables * n_nodes * (n_nodes + 1) * capacity) + fluxes.antidiffusive_flux2_R = unsafe_wrap(Array, pointer(_antidiffusive_flux2_R), + (n_variables, n_nodes, n_nodes + 1, + capacity)) return nothing end diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 70ff346740d..97843db7743 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -14,20 +14,45 @@ function create_cache(mesh::TreeMesh{2}, equations, A3dp1_x = Array{uEltype, 3} A3dp1_y = Array{uEltype, 3} A3d = Array{uEltype, 3} - - fhat1_threaded = A3dp1_x[A3dp1_x(undef, nvariables(equations), nnodes(dg) + 1, - nnodes(dg)) for _ in 1:Threads.nthreads()] - fhat2_threaded = A3dp1_y[A3dp1_y(undef, nvariables(equations), nnodes(dg), - nnodes(dg) + 1) for _ in 1:Threads.nthreads()] + A4d = Array{uEltype, 4} + + fhat1_L_threaded = A3dp1_x[A3dp1_x(undef, nvariables(equations), nnodes(dg) + 1, + nnodes(dg)) for _ in 1:Threads.nthreads()] + fhat2_L_threaded = A3dp1_y[A3dp1_y(undef, nvariables(equations), nnodes(dg), + nnodes(dg) + 1) for _ in 1:Threads.nthreads()] + fhat1_R_threaded = A3dp1_x[A3dp1_x(undef, nvariables(equations), nnodes(dg) + 1, + nnodes(dg)) for _ in 1:Threads.nthreads()] + fhat2_R_threaded = A3dp1_y[A3dp1_y(undef, nvariables(equations), nnodes(dg), + nnodes(dg) + 1) for _ in 1:Threads.nthreads()] flux_temp_threaded = A3d[A3d(undef, nvariables(equations), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] - + fhat_temp_threaded = A3d[A3d(undef, nvariables(equations), nnodes(dg), + nnodes(dg)) + for _ in 1:Threads.nthreads()] antidiffusive_fluxes = Trixi.ContainerAntidiffusiveFlux2D{uEltype}(0, nvariables(equations), nnodes(dg)) - return (; cache..., antidiffusive_fluxes, fhat1_threaded, fhat2_threaded, - flux_temp_threaded) + if have_nonconservative_terms(equations) == true + flux_nonconservative_temp_threaded = A4d[A4d(undef, nvariables(equations), + n_nonconservative_terms(equations), + nnodes(dg), nnodes(dg)) + for _ in 1:Threads.nthreads()] + fhat_nonconservative_temp_threaded = A4d[A4d(undef, nvariables(equations), + n_nonconservative_terms(equations), + nnodes(dg), nnodes(dg)) + for _ in 1:Threads.nthreads()] + phi_threaded = A4d[A4d(undef, nvariables(equations), + n_nonconservative_terms(equations), + nnodes(dg), nnodes(dg)) + for _ in 1:Threads.nthreads()] + cache = (; cache..., flux_nonconservative_temp_threaded, + fhat_nonconservative_temp_threaded, phi_threaded) + end + + return (; cache..., antidiffusive_fluxes, + fhat1_L_threaded, fhat2_L_threaded, fhat1_R_threaded, fhat2_R_threaded, + flux_temp_threaded, fhat_temp_threaded) end function calc_volume_integral!(du, u, @@ -47,19 +72,22 @@ end @inline function subcell_limiting_kernel!(du, u, element, mesh::TreeMesh{2}, - nonconservative_terms::False, equations, + nonconservative_terms, equations, volume_integral, limiter::SubcellLimiterIDP, dg::DGSEM, cache) @unpack inverse_weights = dg.basis @unpack volume_flux_dg, volume_flux_fv = volume_integral # high-order DG fluxes - @unpack fhat1_threaded, fhat2_threaded = cache + @unpack fhat1_L_threaded, fhat1_R_threaded, fhat2_L_threaded, fhat2_R_threaded = cache - fhat1 = fhat1_threaded[Threads.threadid()] - fhat2 = fhat2_threaded[Threads.threadid()] - calcflux_fhat!(fhat1, fhat2, u, mesh, - nonconservative_terms, equations, volume_flux_dg, dg, element, cache) + fhat1_L = fhat1_L_threaded[Threads.threadid()] + fhat1_R = fhat1_R_threaded[Threads.threadid()] + fhat2_L = fhat2_L_threaded[Threads.threadid()] + fhat2_R = fhat2_R_threaded[Threads.threadid()] + calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, mesh, + nonconservative_terms, equations, volume_flux_dg, dg, element, + cache) # low-order FV fluxes @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded = cache @@ -69,12 +97,14 @@ end fstar1_R = fstar1_R_threaded[Threads.threadid()] fstar2_R = fstar2_R_threaded[Threads.threadid()] calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, mesh, - nonconservative_terms, equations, volume_flux_fv, dg, element, cache) + nonconservative_terms, equations, volume_flux_fv, dg, element, + cache) # antidiffusive flux - calcflux_antidiffusive!(fhat1, fhat2, fstar1_L, fstar2_L, u, mesh, - nonconservative_terms, equations, limiter, dg, element, - cache) + calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, + fstar1_L, fstar1_R, fstar2_L, fstar2_R, + u, mesh, nonconservative_terms, equations, limiter, dg, + element, cache) # Calculate volume integral contribution of low-order FV flux for j in eachnode(dg), i in eachnode(dg) @@ -93,7 +123,7 @@ end # (**without non-conservative terms**). # # See also `flux_differencing_kernel!`. -@inline function calcflux_fhat!(fhat1, fhat2, u, +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, mesh::TreeMesh{2}, nonconservative_terms::False, equations, volume_flux, dg::DGSEM, element, cache) @@ -132,11 +162,14 @@ end end # FV-form flux `fhat` in x direction - fhat1[:, 1, :] .= zero(eltype(fhat1)) - fhat1[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1)) + fhat1_L[:, 1, :] .= zero(eltype(fhat1_L)) + fhat1_L[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_L)) + fhat1_R[:, 1, :] .= zero(eltype(fhat1_R)) + fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) for j in eachnode(dg), i in 1:(nnodes(dg) - 1), v in eachvariable(equations) - fhat1[v, i + 1, j] = fhat1[v, i, j] + weights[i] * flux_temp[v, i, j] + fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j] + fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j] end # Split form volume flux in orientation 2: y direction @@ -155,38 +188,278 @@ end end # FV-form flux `fhat` in y direction - fhat2[:, :, 1] .= zero(eltype(fhat2)) - fhat2[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2)) + fhat2_L[:, :, 1] .= zero(eltype(fhat2_L)) + fhat2_L[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_L)) + fhat2_R[:, :, 1] .= zero(eltype(fhat2_R)) + fhat2_R[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_R)) for j in 1:(nnodes(dg) - 1), i in eachnode(dg), v in eachvariable(equations) - fhat2[v, i, j + 1] = fhat2[v, i, j] + weights[j] * flux_temp[v, i, j] + fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j] + fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1] + end + + return nothing +end + +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**with non-conservative terms**). +# +# See also `flux_differencing_kernel!`. +# +# The calculation of the non-conservative staggered "fluxes" requires non-conservative +# terms that can be written as a product of local and a symmetric contributions. See, e.g., +# +# - Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts +# Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +# +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, + mesh::TreeMesh{2}, nonconservative_terms::True, + equations, + volume_flux, dg::DGSEM, element, cache) + @unpack weights, derivative_split = dg.basis + @unpack flux_temp_threaded, flux_nonconservative_temp_threaded = cache + @unpack fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded = cache + + volume_flux_cons, volume_flux_noncons = volume_flux + + flux_temp = flux_temp_threaded[Threads.threadid()] + flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()] + + fhat_temp = fhat_temp_threaded[Threads.threadid()] + fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()] + phi = phi_threaded[Threads.threadid()] + + # The FV-form fluxes are calculated in a recursive manner, i.e.: + # fhat_(0,1) = w_0 * FVol_0, + # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j, for j=1,...,N-1, + # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i). + + # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated + # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing` + # and saved in in `flux_temp`. + + # Split form volume flux in orientation 1: x direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + # All diagonal entries of `derivative_split` are zero. Thus, we can skip + # the computation of the diagonal terms. In addition, we use the symmetry + # of `volume_flux_cons` and `volume_flux_noncons` to save half of the possible two-point flux + # computations. + for ii in (i + 1):nnodes(dg) + u_node_ii = get_node_vars(u, equations, dg, ii, j, element) + flux1 = volume_flux_cons(u_node, u_node_ii, 1, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], flux1, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], flux1, + equations, dg, ii, j) + for noncons in 1:n_nonconservative_terms(equations) + # We multiply by 0.5 because that is done in other parts of Trixi + flux1_noncons = volume_flux_noncons(u_node, u_node_ii, 1, equations, + NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[i, ii], + flux1_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[ii, i], + flux1_noncons, + equations, dg, noncons, ii, j) + end + end end + # FV-form flux `fhat` in x direction + fhat1_L[:, 1, :] .= zero(eltype(fhat1_L)) + fhat1_L[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_L)) + fhat1_R[:, 1, :] .= zero(eltype(fhat1_R)) + fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) + + fhat_temp[:, 1, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + for noncons in 1:n_nonconservative_terms(equations) + set_node_vars!(phi, + volume_flux_noncons(u_local, 1, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j) + end + end + + for j in eachnode(dg), i in 1:(nnodes(dg) - 1) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j] + weights[i] * flux_temp[v, i, j] + fhat_temp[v, i + 1, j] = value + fhat1_L[v, i + 1, j] = value + fhat1_R[v, i + 1, j] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(equations), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j] + + weights[i] * flux_noncons_temp[v, noncons, i, j] + fhat_noncons_temp[v, noncons, i + 1, j] = value + + fhat1_L[v, i + 1, j] = fhat1_L[v, i + 1, j] + phi[v, noncons, i, j] * value + fhat1_R[v, i + 1, j] = fhat1_R[v, i + 1, j] + + phi[v, noncons, i + 1, j] * value + end + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + for jj in (j + 1):nnodes(dg) + u_node_jj = get_node_vars(u, equations, dg, i, jj, element) + flux2 = volume_flux_cons(u_node, u_node_jj, 2, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], flux2, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], flux2, + equations, dg, i, jj) + for noncons in 1:n_nonconservative_terms(equations) + # We multiply by 0.5 because that is done in other parts of Trixi + flux2_noncons = volume_flux_noncons(u_node, u_node_jj, 2, equations, + NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[j, jj], + flux2_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[jj, j], + flux2_noncons, + equations, dg, noncons, i, jj) + end + end + end + + # FV-form flux `fhat` in y direction + fhat2_L[:, :, 1] .= zero(eltype(fhat2_L)) + fhat2_L[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_L)) + fhat2_R[:, :, 1] .= zero(eltype(fhat2_R)) + fhat2_R[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_R)) + + fhat_temp[:, :, 1] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, :, 1] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + for noncons in 1:n_nonconservative_terms(equations) + set_node_vars!(phi, + volume_flux_noncons(u_local, 2, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j) + end + end + + for j in 1:(nnodes(dg) - 1), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j] + weights[j] * flux_temp[v, i, j] + fhat_temp[v, i, j + 1] = value + fhat2_L[v, i, j + 1] = value + fhat2_R[v, i, j + 1] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(equations), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j] + + weights[j] * flux_noncons_temp[v, noncons, i, j] + fhat_noncons_temp[v, noncons, i, j + 1] = value + + fhat2_L[v, i, j + 1] = fhat2_L[v, i, j + 1] + phi[v, noncons, i, j] * value + fhat2_R[v, i, j + 1] = fhat2_R[v, i, j + 1] + + phi[v, noncons, i, j + 1] * value + end + end + + return nothing +end + +# Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems. +@inline function calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, + fstar1_L, fstar1_R, fstar2_L, fstar2_R, + u, mesh, + nonconservative_terms::False, equations, + limiter::SubcellLimiterIDP, dg, element, cache) + @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes + + for j in eachnode(dg), i in 2:nnodes(dg) + for v in eachvariable(equations) + antidiffusive_flux1_L[v, i, j, element] = fhat1_L[v, i, j] - + fstar1_L[v, i, j] + antidiffusive_flux1_R[v, i, j, element] = antidiffusive_flux1_L[v, i, j, + element] + end + end + for j in 2:nnodes(dg), i in eachnode(dg) + for v in eachvariable(equations) + antidiffusive_flux2_L[v, i, j, element] = fhat2_L[v, i, j] - + fstar2_L[v, i, j] + antidiffusive_flux2_R[v, i, j, element] = antidiffusive_flux2_L[v, i, j, + element] + end + end + + antidiffusive_flux1_L[:, 1, :, element] .= zero(eltype(antidiffusive_flux1_L)) + antidiffusive_flux1_L[:, nnodes(dg) + 1, :, element] .= zero(eltype(antidiffusive_flux1_L)) + antidiffusive_flux1_R[:, 1, :, element] .= zero(eltype(antidiffusive_flux1_R)) + antidiffusive_flux1_R[:, nnodes(dg) + 1, :, element] .= zero(eltype(antidiffusive_flux1_R)) + + antidiffusive_flux2_L[:, :, 1, element] .= zero(eltype(antidiffusive_flux2_L)) + antidiffusive_flux2_L[:, :, nnodes(dg) + 1, element] .= zero(eltype(antidiffusive_flux2_L)) + antidiffusive_flux2_R[:, :, 1, element] .= zero(eltype(antidiffusive_flux2_R)) + antidiffusive_flux2_R[:, :, nnodes(dg) + 1, element] .= zero(eltype(antidiffusive_flux2_R)) + return nothing end -# Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar`. -@inline function calcflux_antidiffusive!(fhat1, fhat2, fstar1, fstar2, u, mesh, - nonconservative_terms, equations, +# Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems. +@inline function calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, + fstar1_L, fstar1_R, fstar2_L, fstar2_R, + u, mesh, + nonconservative_terms::True, equations, limiter::SubcellLimiterIDP, dg, element, cache) - @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes + @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes for j in eachnode(dg), i in 2:nnodes(dg) for v in eachvariable(equations) - antidiffusive_flux1[v, i, j, element] = fhat1[v, i, j] - fstar1[v, i, j] + antidiffusive_flux1_L[v, i, j, element] = fhat1_L[v, i, j] - + fstar1_L[v, i, j] + antidiffusive_flux1_R[v, i, j, element] = fhat1_R[v, i, j] - + fstar1_R[v, i, j] end end for j in 2:nnodes(dg), i in eachnode(dg) for v in eachvariable(equations) - antidiffusive_flux2[v, i, j, element] = fhat2[v, i, j] - fstar2[v, i, j] + antidiffusive_flux2_L[v, i, j, element] = fhat2_L[v, i, j] - + fstar2_L[v, i, j] + antidiffusive_flux2_R[v, i, j, element] = fhat2_R[v, i, j] - + fstar2_R[v, i, j] end end - antidiffusive_flux1[:, 1, :, element] .= zero(eltype(antidiffusive_flux1)) - antidiffusive_flux1[:, nnodes(dg) + 1, :, element] .= zero(eltype(antidiffusive_flux1)) + antidiffusive_flux1_L[:, 1, :, element] .= zero(eltype(antidiffusive_flux1_L)) + antidiffusive_flux1_L[:, nnodes(dg) + 1, :, element] .= zero(eltype(antidiffusive_flux1_L)) + antidiffusive_flux1_R[:, 1, :, element] .= zero(eltype(antidiffusive_flux1_R)) + antidiffusive_flux1_R[:, nnodes(dg) + 1, :, element] .= zero(eltype(antidiffusive_flux1_R)) - antidiffusive_flux2[:, :, 1, element] .= zero(eltype(antidiffusive_flux2)) - antidiffusive_flux2[:, :, nnodes(dg) + 1, element] .= zero(eltype(antidiffusive_flux2)) + antidiffusive_flux2_L[:, :, 1, element] .= zero(eltype(antidiffusive_flux2_L)) + antidiffusive_flux2_L[:, :, nnodes(dg) + 1, element] .= zero(eltype(antidiffusive_flux2_L)) + antidiffusive_flux2_R[:, :, 1, element] .= zero(eltype(antidiffusive_flux2_R)) + antidiffusive_flux2_R[:, :, nnodes(dg) + 1, element] .= zero(eltype(antidiffusive_flux2_R)) return nothing end diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 5e00ab4e903..0a72b79ea3f 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -63,7 +63,7 @@ end @inline function idp_positivity!(alpha, limiter, u, dt, semi, variable) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) - (; antidiffusive_flux1, antidiffusive_flux2) = cache.antidiffusive_fluxes + (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes (; inverse_weights) = dg.basis (; positivity_correction_factor) = limiter @@ -91,13 +91,13 @@ end # Calculate Pm # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here. val_flux1_local = inverse_weights[i] * - antidiffusive_flux1[variable, i, j, element] + antidiffusive_flux1_R[variable, i, j, element] val_flux1_local_ip1 = -inverse_weights[i] * - antidiffusive_flux1[variable, i + 1, j, element] + antidiffusive_flux1_L[variable, i + 1, j, element] val_flux2_local = inverse_weights[j] * - antidiffusive_flux2[variable, i, j, element] + antidiffusive_flux2_R[variable, i, j, element] val_flux2_local_jp1 = -inverse_weights[j] * - antidiffusive_flux2[variable, i, j + 1, element] + antidiffusive_flux2_L[variable, i, j + 1, element] Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) + min(0, val_flux2_local) + min(0, val_flux2_local_jp1) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 733f89c2158..dbb9e51121b 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -218,6 +218,11 @@ function set_proposed_dt!(integrator::SimpleIntegratorSSP, dt) integrator.dt = dt end +# used by adaptive timestepping algorithms in DiffEq +function get_proposed_dt(integrator::SimpleIntegratorSSP) + return integrator.dt +end + # stop the time integration function terminate!(integrator::SimpleIntegratorSSP) integrator.finalstep = true diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index af264561027..bd6a95bba50 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -328,6 +328,37 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_mhd_shockcapturing_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"), + l2=[2.9974425783503109e-02, + 7.2849646345685956e-02, + 7.2488477174662239e-02, + 0.0000000000000000e+00, + 1.2507971380965512e+00, + 1.8929505145499678e-02, + 1.2218606317164420e-02, + 0.0000000000000000e+00, + 3.0154796910479838e-03], + linf=[3.2147382412340830e-01, + 1.3709471664007811e+00, + 1.3465154685288383e+00, + 0.0000000000000000e+00, + 1.6051257523415284e+01, + 3.0564266749926644e-01, + 2.3908016329805595e-01, + 0.0000000000000000e+00, + 1.3711262178549158e-01], + tspan=(0.0, 0.003)) + # 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)) < 15000 + end +end end end # module From 28500ea8965bd8ecf54ce35166d2ff6bf35ab212 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 31 Oct 2023 16:03:26 +0100 Subject: [PATCH 2/5] Fix allocs (#1695) * Fix allocs * remove unnecessary code * rerun fmt * format --- src/solvers/dgmulti/dg.jl | 21 +++++---------------- 1 file changed, 5 insertions(+), 16 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 182a486dce5..695260f4b9b 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -465,24 +465,13 @@ function calc_boundary_flux!(cache, t, boundary_conditions::BoundaryConditionPer nothing end -# "lispy tuple programming" instead of for loop for type stability function calc_boundary_flux!(cache, t, boundary_conditions, mesh, have_nonconservative_terms, equations, dg::DGMulti) - - # peel off first boundary condition - calc_single_boundary_flux!(cache, t, first(boundary_conditions), - first(keys(boundary_conditions)), - mesh, have_nonconservative_terms, equations, dg) - - # recurse on the remainder of the boundary conditions - calc_boundary_flux!(cache, t, Base.tail(boundary_conditions), - mesh, have_nonconservative_terms, equations, dg) -end - -# terminate recursion -function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple{(), Tuple{}}, - mesh, have_nonconservative_terms, equations, dg::DGMulti) - nothing + for (key, value) in zip(keys(boundary_conditions), boundary_conditions) + calc_single_boundary_flux!(cache, t, value, + key, + mesh, have_nonconservative_terms, equations, dg) + end end function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, mesh, From abd97b17ed34dd8269bf07e482826e347629811e Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 1 Nov 2023 15:40:37 +0100 Subject: [PATCH 3/5] Allocation tests dgmulti 2d (#1698) --- test/test_dgmulti_2d.jl | 265 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 264 insertions(+), 1 deletion(-) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 8fd00df72ea..892c8ed37f0 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -30,6 +30,14 @@ isdir(outdir) && rm(outdir, recursive = true) 0.002062399194485476, 0.004897700392503701, ]) + # 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_weakform.jl (SBP)" begin @@ -49,6 +57,14 @@ end 0.013593978324845324, 0.03270995869587523, ]) + # 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_weakform.jl (Quadrilateral elements)" begin @@ -68,6 +84,14 @@ end 0.0013568139808290969, 0.0032249020004324613, ]) + # 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_weakform.jl (EC) " begin @@ -88,6 +112,14 @@ end 0.01396529873508534, 0.04227683691807904, ]) + # 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_weakform.jl (SBP, EC)" begin @@ -109,6 +141,14 @@ end 0.05321027922608157, 0.13392025411844033, ]) + # 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_weakform.jl (Quadrilateral elements, SBP, EC)" begin @@ -131,6 +171,14 @@ end 0.012674028479251254, 0.02210545278615017, ]) + # 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_bilinear.jl (Bilinear quadrilateral elements, SBP, flux differencing)" begin @@ -147,6 +195,14 @@ end 6.874188052830021e-5, 0.0001912435192696904, ]) + # 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_curved.jl (Quadrilateral elements, SBP, flux differencing)" begin @@ -163,6 +219,14 @@ end 0.00010003778085621029, 0.00036426282101720275, ]) + # 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_curved.jl (Quadrilateral elements, GaussSBP, flux differencing)" begin @@ -181,6 +245,14 @@ end 3.679582619220412e-5, ], rtol=2 * sqrt(eps())) + # 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_curved.jl (Triangular elements, Polynomial, weak formulation)" begin @@ -199,6 +271,14 @@ end 4.032272526011127e-5, 0.00012013725458537294, ]) + # 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_hohqmesh.jl (Quadrilateral elements, SBP, flux differencing)" begin @@ -215,6 +295,14 @@ end 0.0028721569841545502, 0.011125365074589944, ]) + # 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_weakform.jl (convergence)" begin @@ -228,6 +316,14 @@ end 4.128314378397532, 4.081366752807379, ], rtol = 0.05) + # 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_weakform_periodic.jl" begin @@ -245,6 +341,14 @@ end 0.0019373508504538783, 0.004742686826709086, ]) + # 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_triangulate_pkg_mesh.jl" begin @@ -261,6 +365,14 @@ end 2.6180448559287584e-5, 5.5752932611508044e-5, ]) + # 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 @@ -280,6 +392,14 @@ end 0.1235468309508845, 0.26911424973147735, ]) + # 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 (Quadrilateral elements, GaussSBP)" begin @@ -300,6 +420,14 @@ end 0.12344140473946856, 0.26978428189564774, ]) + # 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_rayleigh_taylor_instability.jl" begin @@ -318,6 +446,14 @@ end 0.040307056293369226, 0.0852365428206836, ]) + # 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_brown_minion_vortex.jl" begin @@ -335,6 +471,14 @@ end 0.021957154972646383, 0.33773439650806303, ]) + # 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_shockcapturing.jl" begin @@ -352,6 +496,14 @@ end 0.1668441768697189, 0.8572572782118661, ]) + # 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_shockcapturing_curved.jl" begin @@ -369,6 +521,14 @@ end 0.16930148707515688, 0.8587706761131937, ]) + # 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_weakform.jl (FD SBP)" begin @@ -395,6 +555,14 @@ end 0.0016243096040242655, 0.004447503691245913, ]) + # 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_weakform.jl (FD SBP, EC)" begin @@ -423,6 +591,14 @@ end 0.003452046522624652, 0.007677153211004928, ]) + # 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_fdsbp_periodic.jl" begin @@ -439,6 +615,14 @@ end 3.9885950273710336e-6, 8.848583042286862e-6, ]) + # 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_fdsbp_periodic.jl (arbitrary reference domain)" begin @@ -456,6 +640,14 @@ end 3.988595024928543e-6, 8.84858303740188e-6, ]) + # 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_fdsbp_periodic.jl (arbitrary reference and physical domains)" begin @@ -479,6 +671,14 @@ end 0.5244468027020814, 1.2210130256735705, ]) + # 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_fdsbp_periodic.jl (CGSEM)" begin @@ -503,6 +703,14 @@ end 5.067812786885284e-5, 9.887976803746312e-5, ]) + # 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_mhd_weak_blast_wave.jl (Quad)" begin @@ -518,6 +726,14 @@ end 0.7450328339751362, 0.06357382685763713, 0.0635738268576378, 0.1058830287485999, 0.005740591170062146]) + # 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_mhd_weak_blast_wave.jl (Tri)" begin @@ -533,6 +749,14 @@ end 0.6906308908961734, 0.05349939593012487, 0.05349939593013042, 0.08830587480616725, 0.0029551359803035027]) + # 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_mhd_weak_blast_wave_SBP.jl (Quad)" begin @@ -601,6 +825,14 @@ end 0.0, 0.0002992631411931389, ]) + # 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_shallowwater_source_terms.jl (Quad, SBP)" begin @@ -619,6 +851,14 @@ end 0.12088392994348407, 9.325873406851315e-15, ]) + # 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_shallowwater_source_terms.jl (Tri, SBP)" begin @@ -637,6 +877,14 @@ end 0.5674183379872275, 1.1546319456101628e-14, ]) + # 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_shallowwater_source_terms.jl (Tri, Polynomial)" begin @@ -657,6 +905,14 @@ end 0.06345885709961152, 3.3989933098554914e-5, ]) + # 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_shallowwater_source_terms.jl (Quad, Polynomial)" begin @@ -679,9 +935,16 @@ end 0.010364675200833062, 1.021405182655144e-14, ]) + # 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 @test_nowarn isdir(outdir) && rm(outdir, recursive = true) From c438d34331a69fdbe9bea831d52b9b518949c446 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 1 Nov 2023 15:41:31 +0100 Subject: [PATCH 4/5] HLLE CEE 2D3D NonCartesian Meshes (#1692) * HLLE CEE 2D3D NonCartesian Meshes * format * hlle via hll * format test * format test * format * do not export hlle * Correct test vals * test values CI * Update src/equations/compressible_euler_2d.jl Co-authored-by: Hendrik Ranocha * Update src/equations/compressible_euler_1d.jl Co-authored-by: Hendrik Ranocha * Update src/equations/compressible_euler_2d.jl Co-authored-by: Hendrik Ranocha * Update src/equations/compressible_euler_3d.jl Co-authored-by: Hendrik Ranocha * Update src/equations/compressible_euler_3d.jl Co-authored-by: Hendrik Ranocha * apply suggestions * additional sentence * Fix typo * typos * correct test vals --------- Co-authored-by: Hendrik Ranocha --- src/equations/compressible_euler_1d.jl | 36 +-------- src/equations/compressible_euler_2d.jl | 92 +++++++++++++++-------- src/equations/compressible_euler_3d.jl | 100 ++++++++++++++++--------- src/equations/numerical_fluxes.jl | 8 ++ test/test_p4est_2d.jl | 26 +++++++ test/test_p4est_3d.jl | 28 +++++++ test/test_tree_1d_euler.jl | 16 ++++ test/test_tree_2d_euler.jl | 26 +++++++ test/test_tree_3d_euler.jl | 28 +++++++ test/test_unit.jl | 18 ++--- 10 files changed, 265 insertions(+), 113 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 9204989e8be..05c38ce791d 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -808,7 +808,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, end """ - flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D) + min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations1D) Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations. Special estimates of the signal velocites and linearization of the Riemann problem developed @@ -825,8 +825,8 @@ Compactly summarized: Numerical methods for conservation laws and related equations. [Link](https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf) """ -function flux_hlle(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerEquations1D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquations1D) # Calculate primitive variables, enthalpy and speed of sound rho_ll, v_ll, p_ll = cons2prim(u_ll, equations) rho_rr, v_rr, p_rr = cons2prim(u_rr, equations) @@ -858,35 +858,7 @@ function flux_hlle(u_ll, u_rr, orientation::Integer, SsL = min(v_roe - c_roe, v_ll - beta * c_ll, zero(v_roe)) SsR = max(v_roe + c_roe, v_rr + beta * c_rr, zero(v_roe)) - if SsL >= 0.0 && SsR > 0.0 - # Positive supersonic speed - f_ll = flux(u_ll, orientation, equations) - - f1 = f_ll[1] - f2 = f_ll[2] - f3 = f_ll[3] - elseif SsR <= 0.0 && SsL < 0.0 - # Negative supersonic speed - f_rr = flux(u_rr, orientation, equations) - - f1 = f_rr[1] - f2 = f_rr[2] - f3 = f_rr[3] - else - # Subsonic case - # Compute left and right fluxes - f_ll = flux(u_ll, orientation, equations) - f_rr = flux(u_rr, orientation, equations) - - f1 = (SsR * f_ll[1] - SsL * f_rr[1] + SsL * SsR * (u_rr[1] - u_ll[1])) / - (SsR - SsL) - f2 = (SsR * f_ll[2] - SsL * f_rr[2] + SsL * SsR * (u_rr[2] - u_ll[2])) / - (SsR - SsL) - f3 = (SsR * f_ll[3] - SsL * f_rr[3] + SsL * SsR * (u_rr[3] - u_ll[3])) / - (SsR - SsL) - end - - return SVector(f1, f2, f3) + return SsL, SsR end @inline function max_abs_speeds(u, equations::CompressibleEulerEquations1D) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 27b92f41953..6c8f2e1e848 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1253,7 +1253,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, end """ - flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) + min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations. Special estimates of the signal velocites and linearization of the Riemann problem developed @@ -1267,8 +1267,8 @@ of the numerical flux. On Godunov-type methods near low densities. [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3) """ -function flux_hlle(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerEquations2D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquations2D) # Calculate primitive variables, enthalpy and speed of sound rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) @@ -1306,39 +1306,65 @@ function flux_hlle(u_ll, u_rr, orientation::Integer, SsR = max(v2_roe + c_roe, v2_rr + beta * c_rr, zero(v2_roe)) end - if SsL >= 0.0 && SsR > 0.0 - # Positive supersonic speed - f_ll = flux(u_ll, orientation, equations) + return SsL, SsR +end - f1 = f_ll[1] - f2 = f_ll[2] - f3 = f_ll[3] - f4 = f_ll[4] - elseif SsR <= 0.0 && SsL < 0.0 - # Negative supersonic speed - f_rr = flux(u_rr, orientation, equations) +""" + min_max_speed_einfeldt(u_ll, u_rr, normal_direction, equations::CompressibleEulerEquations2D) - f1 = f_rr[1] - f2 = f_rr[2] - f3 = f_rr[3] - f4 = f_rr[4] - else - # Subsonic case - # Compute left and right fluxes - f_ll = flux(u_ll, orientation, equations) - f_rr = flux(u_rr, orientation, equations) - - f1 = (SsR * f_ll[1] - SsL * f_rr[1] + SsL * SsR * (u_rr[1] - u_ll[1])) / - (SsR - SsL) - f2 = (SsR * f_ll[2] - SsL * f_rr[2] + SsL * SsR * (u_rr[2] - u_ll[2])) / - (SsR - SsL) - f3 = (SsR * f_ll[3] - SsL * f_rr[3] + SsL * SsR * (u_rr[3] - u_ll[3])) / - (SsR - SsL) - f4 = (SsR * f_ll[4] - SsL * f_rr[4] + SsL * SsR * (u_rr[4] - u_ll[4])) / - (SsR - SsL) - end +Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations. +Special estimates of the signal velocites and linearization of the Riemann problem developed +by Einfeldt to ensure that the internal energy and density remain positive during the computation +of the numerical flux. - return SVector(f1, f2, f3, f4) +- Bernd Einfeldt (1988) + On Godunov-type methods for gas dynamics. + [DOI: 10.1137/0725021](https://doi.org/10.1137/0725021) +- Bernd Einfeldt, Claus-Dieter Munz, Philip L. Roe and Björn Sjögreen (1991) + On Godunov-type methods near low densities. + [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3) +""" +@inline function min_max_speed_einfeldt(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + # Calculate primitive variables, enthalpy and speed of sound + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + + # `u_ll[4]` is total energy `rho_e_ll` on the left + H_ll = (u_ll[4] + p_ll) / rho_ll + c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_ + + # `u_rr[4]` is total energy `rho_e_rr` on the right + H_rr = (u_rr[4] + p_rr) / rho_rr + c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + inv_sum_sqrt_rho = inv(sqrt_rho_ll + sqrt_rho_rr) + + v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) * inv_sum_sqrt_rho + v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) * inv_sum_sqrt_rho + v_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + v_roe_mag = (v1_roe * normal_direction[1])^2 + (v2_roe * normal_direction[2])^2 + + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) * norm_ + + # Compute convenience constant for positivity preservation, see + # https://doi.org/10.1016/0021-9991(91)90211-3 + beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + + # Estimate the edges of the Riemann fan (with positivity conservation) + SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(v_roe)) + SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(v_roe)) + + return SsL, SsR end @inline function max_abs_speeds(u, equations::CompressibleEulerEquations2D) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 7f25bde31fd..fc56f58025b 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1322,7 +1322,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, end """ - flux_hlle(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D) + min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D) Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations. Special estimates of the signal velocites and linearization of the Riemann problem developed @@ -1336,8 +1336,8 @@ of the numerical flux. On Godunov-type methods near low densities. [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3) """ -function flux_hlle(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerEquations3D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquations3D) # Calculate primitive variables, enthalpy and speed of sound rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations) rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations) @@ -1379,43 +1379,69 @@ function flux_hlle(u_ll, u_rr, orientation::Integer, SsR = max(v3_roe + c_roe, v3_rr + beta * c_rr, zero(v3_roe)) end - if SsL >= 0.0 && SsR > 0.0 - # Positive supersonic speed - f_ll = flux(u_ll, orientation, equations) + return SsL, SsR +end - f1 = f_ll[1] - f2 = f_ll[2] - f3 = f_ll[3] - f4 = f_ll[4] - f5 = f_ll[5] - elseif SsR <= 0.0 && SsL < 0.0 - # Negative supersonic speed - f_rr = flux(u_rr, orientation, equations) +""" + min_max_speed_einfeldt(u_ll, u_rr, normal_direction, equations::CompressibleEulerEquations3D) - f1 = f_rr[1] - f2 = f_rr[2] - f3 = f_rr[3] - f4 = f_rr[4] - f5 = f_rr[5] - else - # Subsonic case - # Compute left and right fluxes - f_ll = flux(u_ll, orientation, equations) - f_rr = flux(u_rr, orientation, equations) - - f1 = (SsR * f_ll[1] - SsL * f_rr[1] + SsL * SsR * (u_rr[1] - u_ll[1])) / - (SsR - SsL) - f2 = (SsR * f_ll[2] - SsL * f_rr[2] + SsL * SsR * (u_rr[2] - u_ll[2])) / - (SsR - SsL) - f3 = (SsR * f_ll[3] - SsL * f_rr[3] + SsL * SsR * (u_rr[3] - u_ll[3])) / - (SsR - SsL) - f4 = (SsR * f_ll[4] - SsL * f_rr[4] + SsL * SsR * (u_rr[4] - u_ll[4])) / - (SsR - SsL) - f5 = (SsR * f_ll[5] - SsL * f_rr[5] + SsL * SsR * (u_rr[5] - u_ll[5])) / - (SsR - SsL) - end +Computes the HLLE (Harten-Lax-van Leer-Einfeldt) flux for the compressible Euler equations. +Special estimates of the signal velocites and linearization of the Riemann problem developed +by Einfeldt to ensure that the internal energy and density remain positive during the computation +of the numerical flux. - return SVector(f1, f2, f3, f4, f5) +- Bernd Einfeldt (1988) + On Godunov-type methods for gas dynamics. + [DOI: 10.1137/0725021](https://doi.org/10.1137/0725021) +- Bernd Einfeldt, Claus-Dieter Munz, Philip L. Roe and Björn Sjögreen (1991) + On Godunov-type methods near low densities. + [DOI: 10.1016/0021-9991(91)90211-3](https://doi.org/10.1016/0021-9991(91)90211-3) +""" +@inline function min_max_speed_einfeldt(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations3D) + # Calculate primitive variables, enthalpy and speed of sound + rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations) + + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + norm_ = norm(normal_direction) + + # `u_ll[5]` is total energy `rho_e_ll` on the left + H_ll = (u_ll[5] + p_ll) / rho_ll + c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_ + + # `u_rr[5]` is total energy `rho_e_rr` on the right + H_rr = (u_rr[5] + p_rr) / rho_rr + c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + inv_sum_sqrt_rho = inv(sqrt_rho_ll + sqrt_rho_rr) + + v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) * inv_sum_sqrt_rho + v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) * inv_sum_sqrt_rho + v3_roe = (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr) * inv_sum_sqrt_rho + v_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + + v3_roe * normal_direction[3] + v_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2 + + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) * norm_ + + # Compute convenience constant for positivity preservation, see + # https://doi.org/10.1016/0021-9991(91)90211-3 + beta = sqrt(0.5 * (equations.gamma - 1) / equations.gamma) + + # Estimate the edges of the Riemann fan (with positivity conservation) + SsL = min(v_roe - c_roe, v_dot_n_ll - beta * c_ll, zero(v_roe)) + SsR = max(v_roe + c_roe, v_dot_n_rr + beta * c_rr, zero(v_roe)) + + return SsL, SsR end @inline function max_abs_speeds(u, equations::CompressibleEulerEquations3D) diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 87010275f2c..71782644b17 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -304,6 +304,14 @@ See [`FluxHLL`](@ref). """ const flux_hll = FluxHLL() +""" + flux_hlle + +See [`min_max_speed_einfeldt`](@ref). +This is a [`FluxHLL`](@ref)-type two-wave solver with special estimates of the wave speeds. +""" +const flux_hlle = FluxHLL(min_max_speed_einfeldt) + # TODO: TrixiShallowWater: move the chen_noelle flux structure to the new package # An empty version of the `min_max_speed_chen_noelle` function is declared here diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index cf7250b246b..546e5bff8a6 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -200,6 +200,32 @@ end end end +@trixi_testset "elixir_euler_sedov.jl (HLLE)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), + l2=[ + 0.411541263324004, + 0.2558929632770186, + 0.2558929632770193, + 1.298715766843915, + ], + linf=[ + 1.3457201726152221, + 1.3138961427140758, + 1.313896142714079, + 6.293305112638921, + ], + surface_flux=flux_hlle, + tspan=(0.0, 0.3)) + # 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_blast_wave_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"), l2=[ diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 63077deb436..346c61c7448 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -287,6 +287,34 @@ end end end +@trixi_testset "elixir_euler_sedov.jl (HLLE)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), + l2=[ + 0.09946224487902565, + 0.04863386374672001, + 0.048633863746720116, + 0.04863386374672032, + 0.3751015774232693, + ], + linf=[ + 0.789241521871487, + 0.42046970270100276, + 0.42046970270100276, + 0.4204697027010028, + 4.730877375538398, + ], + tspan=(0.0, 0.3), + surface_flux=flux_hlle) + # 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_source_terms_nonconforming_earth.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_earth.jl"), diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index f01124509cb..62afc5baee3 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -281,6 +281,22 @@ end end end +@trixi_testset "elixir_euler_sedov_blast_wave.jl (HLLE)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"), + l2=[0.6442208390304879, 0.508817280068289, 0.9482809853033687], + linf=[3.007059066482486, 2.4678899558345506, 2.3952311739389787], + tspan=(0.0, 0.5), + surface_flux=flux_hlle) + # 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_sedov_blast_wave_pure_fv.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_pure_fv.jl"), diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index af7c5d8324c..db36cb7d79f 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -422,6 +422,32 @@ end end end +@trixi_testset "elixir_euler_sedov_blast_wave.jl (HLLE)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"), + l2=[ + 0.35267161504176747, + 0.17218309138797958, + 0.17218307467125854, + 0.6236143054619037, + ], + linf=[ + 2.77484045816607, + 1.8281111268370718, + 1.8281110470490887, + 6.24263735888126, + ], + tspan=(0.0, 0.5), + surface_flux=flux_hlle) + # 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_sedov_blast_wave_neuralnetwork_perssonperaire.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl"), diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl index d96dcff3977..02e657e001a 100644 --- a/test/test_tree_3d_euler.jl +++ b/test/test_tree_3d_euler.jl @@ -473,6 +473,34 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_sedov_blast_wave.jl (HLLE)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"), + l2=[ + 0.0007871241159752619, + 0.0037168004033428146, + 0.0037168004033428094, + 0.0037168004033428514, + 0.011119869089205635, + ], + linf=[ + 0.13982864363612468, + 0.786004687738243, + 0.786004687738243, + 0.7860046877382431, + 1.7082524045150382, + ], + tspan=(0.0, 0.01), + surface_flux=flux_hlle) + # 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 end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index 5422db7a66d..a73dfab5504 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -911,7 +911,7 @@ end SVector(-1.2, 0.3)] for normal_direction in normal_directions - @test flux_hll(u, u, normal_direction, equations) ≈ + @test flux_hlle(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) end @@ -930,18 +930,15 @@ end SVector(-1.2, 0.3, 1.4)] for normal_direction in normal_directions - @test flux_hll(u, u, normal_direction, equations) ≈ + @test flux_hlle(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) end end @timed_testset "Consistency check for HLLE flux: SWE" begin - # Test HLL flux with min_max_speed_einfeldt - flux_hll = FluxHLL(min_max_speed_einfeldt) - equations = ShallowWaterEquations1D(gravity_constant = 9.81) u = SVector(1, 0.5, 0.0) - @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations) + @test flux_hlle(u, u, 1, equations) ≈ flux(u, 1, equations) equations = ShallowWaterEquations2D(gravity_constant = 9.81) normal_directions = [SVector(1.0, 0.0), @@ -953,18 +950,17 @@ end u = SVector(1, 0.5, 0.5, 0.0) for orientation in orientations - @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations) + @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations) end for normal_direction in normal_directions - @test flux_hll(u, u, normal_direction, equations) ≈ + @test flux_hlle(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) end end @timed_testset "Consistency check for HLLE flux: MHD" begin - # Test HLL flux with min_max_speed_einfeldt - flux_hll = FluxHLL(min_max_speed_naive) + # Note: min_max_speed_naive for MHD is essentially min_max_speed_einfeldt equations = IdealGlmMhdEquations1D(1.4) u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1), @@ -1242,7 +1238,7 @@ end u = SVector(1, 0.5, 0.5, 0.0) fluxes = [flux_central, flux_fjordholm_etal, flux_wintermeyer_etal, - flux_hll, FluxHLL(min_max_speed_davis), FluxHLL(min_max_speed_einfeldt)] + flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] end @timed_testset "IdealGlmMhdEquations2D" begin From 08a2e09fe41460c9f15c77cb084e8609c5d39e55 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 2 Nov 2023 06:10:43 +0100 Subject: [PATCH 5/5] Bump crate-ci/typos from 1.16.15 to 1.16.21 (#1700) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.16.15 to 1.16.21. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.16.15...v1.16.21) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index eae6d8e0be9..fb71fe45af9 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.16.15 + uses: crate-ci/typos@v1.16.21