From 99c17828996e95c7de9309efb67a5f5f0bd26056 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 5 Oct 2023 17:47:09 +0200 Subject: [PATCH 01/16] Added first (ugly) implementation of subcell limiting for non-conservative systems -> A working version of this implementation is added for the GLM-MHD system -> The flux-differencing formula requires non-conservative terms of the form (local * symmetric)... I modified equations/ideal_glm_mhd_2d.jl and solvers/dgsem_tree/dg_2d.jl to make it work -> In this first implementation, we only use the Powell term and deactivate the GLM term --- .../elixir_mhd_shockcapturing_subcell.jl | 108 ++++++++ .../subcell_limiter_idp_correction_2d.jl | 10 +- src/equations/ideal_glm_mhd_2d.jl | 124 ++++++++- src/solvers/dgsem_tree/containers_2d.jl | 62 +++-- src/solvers/dgsem_tree/dg_2d.jl | 16 +- .../dgsem_tree/dg_2d_subcell_limiters.jl | 245 +++++++++++++++--- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 10 +- src/time_integration/methods_SSP.jl | 5 + 8 files changed, 508 insertions(+), 72 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..37294973315 --- /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 taken 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) +""" +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 = 10.0 + + rhomin = 0.5 + rhomax = 1.0 + 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) +volume_flux = (flux_derigs_etal, flux_nonconservative_powell) #central +basis = LobattoLegendreBasis(3) + +#volume_integral=VolumeIntegralFluxDifferencing(volume_flux) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons=[1], + positivity_correction_factor=0.8) +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 = 30 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=10) #analysis_interval + +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +cfl = 0.3 +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(),) + +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/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/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 8fef1ee22c9..82d4844b166 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -23,12 +23,17 @@ mutable struct IdealGlmMhdEquations2D{RealT <: Real} <: end end +struct NonConservativeLocal end +struct NonConservativeSymmetric end + function IdealGlmMhdEquations2D(gamma; initial_c_h = convert(typeof(gamma), NaN)) # Use `promote` to ensure that `gamma` and `initial_c_h` have the same type IdealGlmMhdEquations2D(promote(gamma, initial_c_h)...) end have_nonconservative_terms(::IdealGlmMhdEquations2D) = True() +nnoncons(::IdealGlmMhdEquations2D) = 2 + function varnames(::typeof(cons2cons), ::IdealGlmMhdEquations2D) ("rho", "rho_v1", "rho_v2", "rho_v3", "rho_e", "B1", "B2", "B3", "psi") end @@ -128,10 +133,10 @@ end f4 = rho_v1 * v3 - B1 * B3 f5 = (kin_en + equations.gamma * p_over_gamma_minus_one + 2 * mag_en) * v1 - B1 * (v1 * B1 + v2 * B2 + v3 * B3) + equations.c_h * psi * B1 - f6 = equations.c_h * psi + f6 = 0#equations.c_h * psi f7 = v1 * B2 - v2 * B1 f8 = v1 * B3 - v3 * B1 - f9 = equations.c_h * B1 + f9 = 0#equations.c_h * B1 else #if orientation == 2 f1 = rho_v2 f2 = rho_v2 * v1 - B2 * B1 @@ -140,9 +145,9 @@ end f5 = (kin_en + equations.gamma * p_over_gamma_minus_one + 2 * mag_en) * v2 - B2 * (v1 * B1 + v2 * B2 + v3 * B3) + equations.c_h * psi * B2 f6 = v2 * B1 - v1 * B2 - f7 = equations.c_h * psi + f7 = 0#equations.c_h * psi f8 = v2 * B3 - v3 * B2 - f9 = equations.c_h * B2 + f9 = 0#equations.c_h * B2 end return SVector(f1, f2, f3, f4, f5, f6, f7, f8, f9) @@ -206,7 +211,7 @@ terms. equations. Part I: Theory and numerical verification [DOI: 10.1016/j.jcp.2018.06.027](https://doi.org/10.1016/j.jcp.2018.06.027) """ -@inline function flux_nonconservative_powell(u_ll, u_rr, orientation::Integer, +#= @inline function flux_nonconservative_powell(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 @@ -240,6 +245,115 @@ terms. v2_ll * psi_rr) end + return f +end =# +@inline function flux_nonconservative_powell(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 = 0.5 * (psi_ll + psi_rr) + if orientation == 1 + B1_avg = 0.5 * (B1_ll + B1_rr) + 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_rr, + v1_ll * B1_avg, + v2_ll * B1_avg, + v3_ll * B1_avg, + 0)#v1_ll * psi_avg) + else # orientation == 2 + B2_avg = 0.5 * (B2_ll + B2_rr) + 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_rr, + v1_ll * B2_avg, + v2_ll * B2_avg, + v3_ll * B2_avg, + 0)#v2_ll * psi_avg) + end + + return f +end +""" + +""" +@inline function flux_nonconservative_powell(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeLocal) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + + 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}) + f = SVector(0, + B1_ll, + B2_ll, + B3_ll, + v_dot_B_ll, # The term (v1_ll * psi_ll) is missing because we need to define several non-conservative terms + v1_ll, + v2_ll, + v3_ll, + 0)#v1_ll) + + return f +end +""" + +""" +@inline function flux_nonconservative_powell(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeSymmetric) + 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 = 0.5 * (psi_ll + psi_rr) + if orientation == 1 + B1_avg = 0.5 * (B1_ll + B1_rr) + f = SVector(0, + B1_avg, + B1_avg, + B1_avg, + B1_avg, # The term (psi_avg) is missing because we need to define several non-conservative terms + B1_avg, + B1_avg, + B1_avg, + 0)#psi_avg) + else # orientation == 2 + B2_avg = 0.5 * (B2_ll + B2_rr) + f = SVector(0, + B2_avg, + B2_avg, + B2_avg, + B2_avg, # The term (psi_avg) is missing because we need to define several non-conservative terms + B2_avg, + B2_avg, + B2_avg, + 0)#psi_avg) + end + return f end diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 9148b936312..576732b9219 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, + _antidiffusive_flux1_L = fill(nan_uEltype, n_variables * (n_nodes + 1) * n_nodes * capacity) - antidiffusive_flux1 = unsafe_wrap(Array, pointer(_antidiffusive_flux1), + 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 = fill(nan_uEltype, + _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 = unsafe_wrap(Array, pointer(_antidiffusive_flux2), + antidiffusive_flux2_R = unsafe_wrap(Array, pointer(_antidiffusive_flux2_R), (n_variables, n_nodes, n_nodes + 1, capacity)) - return ContainerAntidiffusiveFlux2D{uEltype}(antidiffusive_flux1, - antidiffusive_flux2, - _antidiffusive_flux1, - _antidiffusive_flux2) + 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,14 +1322,22 @@ function Base.resize!(fluxes::ContainerAntidiffusiveFlux2D, capacity) n_nodes = nnodes(fluxes) n_variables = nvariables(fluxes) - @unpack _antidiffusive_flux1, _antidiffusive_flux2 = fluxes + @unpack _antidiffusive_flux1_L, _antidiffusive_flux2_L, _antidiffusive_flux1_R, _antidiffusive_flux2_R = fluxes - resize!(_antidiffusive_flux1, n_variables * (n_nodes + 1) * n_nodes * capacity) - fluxes.antidiffusive_flux1 = unsafe_wrap(Array, pointer(_antidiffusive_flux1), + 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_flux2, n_variables * n_nodes * (n_nodes + 1) * capacity) - fluxes.antidiffusive_flux2 = unsafe_wrap(Array, pointer(_antidiffusive_flux2), + 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)) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index c30d0a8e01a..679b6b33f25 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -316,7 +316,7 @@ end end # The factor 0.5 cancels the factor 2 in the flux differencing form - multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, + multiply_add_to_node_vars!(du, alpha, integral_contribution, equations, dg, i, j, element) end end @@ -493,8 +493,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - f1_L = f1 + 0.5 * nonconservative_flux(u_ll, u_rr, 1, equations) - f1_R = f1 + 0.5 * nonconservative_flux(u_rr, u_ll, 1, equations) + f1_L = f1 + nonconservative_flux(u_ll, u_rr, 1, equations) + f1_R = f1 + nonconservative_flux(u_rr, u_ll, 1, equations) # Copy to temporary storage set_node_vars!(fstar1_L, f1_L, equations, dg, i, j) @@ -519,8 +519,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - f2_L = f2 + 0.5 * nonconservative_flux(u_ll, u_rr, 2, equations) - f2_R = f2 + 0.5 * nonconservative_flux(u_rr, u_ll, 2, equations) + f2_L = f2 + nonconservative_flux(u_ll, u_rr, 2, equations) + f2_R = f2 + nonconservative_flux(u_rr, u_ll, 2, equations) # Copy to temporary storage set_node_vars!(fstar2_L, f2_L, equations, dg, i, j) @@ -626,10 +626,10 @@ function calc_interface_flux!(surface_flux_values, # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, i, left_direction, left_id] = flux[v] + - 0.5 * + #0.5 * noncons_left[v] surface_flux_values[v, i, right_direction, right_id] = flux[v] + - 0.5 * + #0.5 * noncons_right[v] end end @@ -779,7 +779,7 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A # Copy flux to left and right element storage for v in eachvariable(equations) surface_flux_values[v, i, direction, neighbor] = flux[v] + - 0.5 * noncons_flux[v] + noncons_flux[v] end end 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..57053f97e02 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -15,19 +15,24 @@ function create_cache(mesh::TreeMesh{2}, equations, A3dp1_y = Array{uEltype, 3} A3d = Array{uEltype, 3} - fhat1_threaded = A3dp1_x[A3dp1_x(undef, nvariables(equations), nnodes(dg) + 1, + fhat1_L_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), + 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()] - + flux_temp_nonconservative_threaded = A3d[A3d(undef, nvariables(equations), nnodes(dg), nnodes(dg)) #nnoncons(equations) + 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) + return (; cache..., antidiffusive_fluxes, + fhat1_L_threaded, fhat2_L_threaded, fhat1_R_threaded, fhat2_R_threaded, + flux_temp_threaded, flux_temp_nonconservative_threaded) end function calc_volume_integral!(du, u, @@ -47,18 +52,20 @@ 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, + 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 @@ -72,7 +79,9 @@ end nonconservative_terms, equations, volume_flux_fv, dg, element, cache) # antidiffusive flux - calcflux_antidiffusive!(fhat1, fhat2, fstar1_L, fstar2_L, u, mesh, + 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) @@ -93,7 +102,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) @@ -101,7 +110,7 @@ end @unpack flux_temp_threaded = cache flux_temp = flux_temp_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, @@ -132,11 +141,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 +167,211 @@ 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!`. +@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_temp_nonconservative_threaded = cache + + volume_flux_cons, volume_flux_noncons = volume_flux + + flux_temp = flux_temp_threaded[Threads.threadid()] + flux_temp_noncons = flux_temp_nonconservative_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_temp_noncons .= zero(eltype(flux_temp_noncons)) + + 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) + flux1_noncons = volume_flux_noncons(u_node, u_node_ii, 1, equations, NonConservativeSymmetric()) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[i, ii], flux1_noncons, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[ii, i], flux1_noncons, + equations, dg, ii, j) + 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 = zero(MVector{nvariables(equations), eltype(fhat1_L)}) + fhat_noncons_temp = zero(MVector{nvariables(equations), eltype(fhat1_L)}) + + for j in eachnode(dg) + fhat_temp .= zero(eltype(fhat1_L)) + fhat_noncons_temp .= zero(eltype(fhat1_L)) + for i in 1:(nnodes(dg) - 1) + # Get the local contribution to the nonconservative flux + u_node_L = get_node_vars(u, equations, dg, i, j, element) + phi_L = volume_flux_noncons(u_node_L, 1, equations, NonConservativeLocal()) + + u_node_R = get_node_vars(u, equations, dg, i + 1, j, element) + phi_R = volume_flux_noncons(u_node_R, 1, equations, NonConservativeLocal()) + for v in eachvariable(equations) + fhat_temp[v] = fhat_temp[v] + weights[i] * flux_temp[v, i, j] + fhat_noncons_temp[v] = fhat_noncons_temp[v] + weights[i] * flux_temp_noncons[v, i, j] + + fhat1_L[v, i + 1, j] = fhat_temp[v] + phi_L[v] * fhat_noncons_temp[v] + fhat1_R[v, i + 1, j] = fhat_temp[v] + phi_R[v] * fhat_noncons_temp[v] + end + end + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + flux_temp_noncons .= zero(eltype(flux_temp_noncons)) + + 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) + flux2_noncons = volume_flux_noncons(u_node, u_node_jj, 2, equations, NonConservativeSymmetric()) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[j, jj], flux2_noncons, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[jj, j], flux2_noncons, + equations, dg, i, jj) + 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)) + + for i in eachnode(dg) + fhat_temp .= zero(eltype(fhat1_L)) + fhat_noncons_temp .= zero(eltype(fhat1_L)) + for j in 1:(nnodes(dg) - 1) + # Get the local contribution to the nonconservative flux + u_node_L = get_node_vars(u, equations, dg, i, j, element) + phi_L = volume_flux_noncons(u_node_L, 2, equations, NonConservativeLocal()) + + u_node_R = get_node_vars(u, equations, dg, i, j + 1, element) + phi_R = volume_flux_noncons(u_node_R, 2, equations, NonConservativeLocal()) + for v in eachvariable(equations) + fhat_temp[v] = fhat_temp[v] + weights[j] * flux_temp[v, i, j] + fhat_noncons_temp[v] = fhat_noncons_temp[v] + weights[j] * flux_temp_noncons[v, i, j] + + fhat2_L[v, i, j + 1] = fhat_temp[v] + phi_L[v] * fhat_noncons_temp[v] + fhat2_R[v, i, j + 1] = fhat_temp[v] + phi_R[v] * fhat_noncons_temp[v] + end + end end 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::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` 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 09ab84ed11a..07736c6bfa3 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -59,7 +59,7 @@ end @inline function idp_positivity!(alpha, limiter, u, dt, semi, variable, index) mesh, equations, dg, cache = mesh_equations_solver_cache(semi) - @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 inverse_weights = dg.basis @unpack positivity_correction_factor = limiter @@ -88,13 +88,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 a0ed889968a..fed3fb70dde 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -213,6 +213,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 From 67e379ff17af6a84cbddfeb63b712528434bedf1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Fri, 6 Oct 2023 15:03:56 +0200 Subject: [PATCH 02/16] Modified non-conservative fluxes to revert src/solvers/dgsem_tree/dg_2d.jl back to its original state --- src/equations/ideal_glm_mhd_2d.jl | 12 ++++++------ src/solvers/dgsem_tree/dg_2d.jl | 16 ++++++++-------- src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl | 6 ++++-- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 82d4844b166..c032a778c3d 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -259,9 +259,9 @@ end =# # 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 = 0.5 * (psi_ll + psi_rr) + psi_avg = (psi_ll + psi_rr) #* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 if orientation == 1 - B1_avg = 0.5 * (B1_ll + B1_rr) + B1_avg = (B1_ll + B1_rr) #* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 f = SVector(0, B1_ll * B1_avg, B2_ll * B1_avg, @@ -272,7 +272,7 @@ end =# v3_ll * B1_avg, 0)#v1_ll * psi_avg) else # orientation == 2 - B2_avg = 0.5 * (B2_ll + B2_rr) + B2_avg = (B2_ll + B2_rr) #* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 f = SVector(0, B1_ll * B2_avg, B2_ll * B2_avg, @@ -329,9 +329,9 @@ end # 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 = 0.5 * (psi_ll + psi_rr) + psi_avg = (psi_ll + psi_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 if orientation == 1 - B1_avg = 0.5 * (B1_ll + B1_rr) + B1_avg = (B1_ll + B1_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 f = SVector(0, B1_avg, B1_avg, @@ -342,7 +342,7 @@ end B1_avg, 0)#psi_avg) else # orientation == 2 - B2_avg = 0.5 * (B2_ll + B2_rr) + B2_avg = (B2_ll + B2_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 f = SVector(0, B2_avg, B2_avg, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 679b6b33f25..c30d0a8e01a 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -316,7 +316,7 @@ end end # The factor 0.5 cancels the factor 2 in the flux differencing form - multiply_add_to_node_vars!(du, alpha, integral_contribution, equations, + multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, dg, i, j, element) end end @@ -493,8 +493,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - f1_L = f1 + nonconservative_flux(u_ll, u_rr, 1, equations) - f1_R = f1 + nonconservative_flux(u_rr, u_ll, 1, equations) + f1_L = f1 + 0.5 * nonconservative_flux(u_ll, u_rr, 1, equations) + f1_R = f1 + 0.5 * nonconservative_flux(u_rr, u_ll, 1, equations) # Copy to temporary storage set_node_vars!(fstar1_L, f1_L, equations, dg, i, j) @@ -519,8 +519,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - f2_L = f2 + nonconservative_flux(u_ll, u_rr, 2, equations) - f2_R = f2 + nonconservative_flux(u_rr, u_ll, 2, equations) + f2_L = f2 + 0.5 * nonconservative_flux(u_ll, u_rr, 2, equations) + f2_R = f2 + 0.5 * nonconservative_flux(u_rr, u_ll, 2, equations) # Copy to temporary storage set_node_vars!(fstar2_L, f2_L, equations, dg, i, j) @@ -626,10 +626,10 @@ function calc_interface_flux!(surface_flux_values, # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, i, left_direction, left_id] = flux[v] + - #0.5 * + 0.5 * noncons_left[v] surface_flux_values[v, i, right_direction, right_id] = flux[v] + - #0.5 * + 0.5 * noncons_right[v] end end @@ -779,7 +779,7 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A # Copy flux to left and right element storage for v in eachvariable(equations) surface_flux_values[v, i, direction, neighbor] = flux[v] + - noncons_flux[v] + 0.5 * noncons_flux[v] end end end diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 57053f97e02..7c8bdff4b0e 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -222,7 +222,8 @@ end equations, dg, i, j) multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], flux1, equations, dg, ii, j) - flux1_noncons = volume_flux_noncons(u_node, u_node_ii, 1, equations, NonConservativeSymmetric()) + # We multiply by 0.5 because that is done in other parts of Trixi + flux1_noncons = 0.5 * volume_flux_noncons(u_node, u_node_ii, 1, equations, NonConservativeSymmetric()) multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[i, ii], flux1_noncons, equations, dg, i, j) multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[ii, i], flux1_noncons, @@ -272,7 +273,8 @@ end equations, dg, i, j) multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], flux2, equations, dg, i, jj) - flux2_noncons = volume_flux_noncons(u_node, u_node_jj, 2, equations, NonConservativeSymmetric()) + # We multiply by 0.5 because that is done in other parts of Trixi + flux2_noncons = 0.5 * volume_flux_noncons(u_node, u_node_jj, 2, equations, NonConservativeSymmetric()) multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[j, jj], flux2_noncons, equations, dg, i, j) multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[jj, j], flux2_noncons, From 12c6c1d56dc9780839a92e6c8adfcd64a37da89e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 10 Oct 2023 11:51:12 +0200 Subject: [PATCH 03/16] Subcell limiting: Added the possibility to use multiple nonconservative terms --- .../elixir_mhd_shockcapturing_subcell.jl | 17 +- src/Trixi.jl | 2 +- src/equations/equations.jl | 6 + src/equations/ideal_glm_mhd_2d.jl | 165 +++++++++--------- .../dgsem_tree/dg_2d_subcell_limiters.jl | 91 ++++++---- 5 files changed, 153 insertions(+), 128 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl index 37294973315..db362b74cdb 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -19,11 +19,11 @@ 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 = 10.0 - - rhomin = 0.5 + pmin = 1.0 rhomax = 1.0 + rhomin = 0.01 if r <= 0.09 p = pmax rho = rhomax @@ -46,13 +46,12 @@ end initial_condition = initial_condition_blast_wave surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell) -volume_flux = (flux_derigs_etal, flux_nonconservative_powell) #central +volume_flux = (flux_derigs_etal, flux_nonconservative_powell) basis = LobattoLegendreBasis(3) -#volume_integral=VolumeIntegralFluxDifferencing(volume_flux) limiter_idp = SubcellLimiterIDP(equations, basis; positivity_variables_cons=[1], - positivity_correction_factor=0.8) + positivity_correction_factor=0.5) volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg=volume_flux, volume_flux_fv=surface_flux) @@ -76,10 +75,10 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 30 +analysis_interval = 100 analysis_callback = AnalysisCallback(semi, interval=analysis_interval) -alive_callback = AliveCallback(analysis_interval=10) #analysis_interval +alive_callback = AliveCallback(analysis_interval=analysis_interval) save_solution = SaveSolutionCallback(interval=100, save_initial_solution=true, @@ -102,7 +101,7 @@ callbacks = CallbackSet(summary_callback, # run the simulation stage_callbacks = (SubcellLimiterIDPCorrection(),) -sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks=stage_callbacks); # +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 b65d03e7975..b4485217411 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -60,7 +60,7 @@ using RecipesBase: RecipesBase using Requires: @require using Static: Static, One, True, False @reexport using StaticArrays: SVector -using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix +using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix, MMatrix using StrideArrays: PtrArray, StrideArray, StaticInt @reexport using StructArrays: StructArrays, StructArray using TimerOutputs: TimerOutputs, @notimeit, TimerOutput, print_timer, reset_timer! diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 9bae563d85f..a941f750a68 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -220,6 +220,12 @@ 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() +""" + nnoncons(equations) +Number of nonconservative terms for a particular equation. The default is 0 and +it must be defined for each nonconservative equation independently. +""" +nnoncons(::AbstractEquations) = 0 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 c032a778c3d..9a26503397d 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -133,10 +133,10 @@ end f4 = rho_v1 * v3 - B1 * B3 f5 = (kin_en + equations.gamma * p_over_gamma_minus_one + 2 * mag_en) * v1 - B1 * (v1 * B1 + v2 * B2 + v3 * B3) + equations.c_h * psi * B1 - f6 = 0#equations.c_h * psi + f6 = equations.c_h * psi f7 = v1 * B2 - v2 * B1 f8 = v1 * B3 - v3 * B1 - f9 = 0#equations.c_h * B1 + f9 = equations.c_h * B1 else #if orientation == 2 f1 = rho_v2 f2 = rho_v2 * v1 - B2 * B1 @@ -145,9 +145,9 @@ end f5 = (kin_en + equations.gamma * p_over_gamma_minus_one + 2 * mag_en) * v2 - B2 * (v1 * B1 + v2 * B2 + v3 * B3) + equations.c_h * psi * B2 f6 = v2 * B1 - v1 * B2 - f7 = 0#equations.c_h * psi + f7 = equations.c_h * psi f8 = v2 * B3 - v3 * B2 - f9 = 0#equations.c_h * B2 + f9 = equations.c_h * B2 end return SVector(f1, f2, f3, f4, f5, f6, f7, f8, f9) @@ -211,42 +211,6 @@ terms. equations. Part I: Theory and numerical verification [DOI: 10.1016/j.jcp.2018.06.027](https://doi.org/10.1016/j.jcp.2018.06.027) """ -#= @inline function flux_nonconservative_powell(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}) - if orientation == 1 - f = SVector(0, - B1_ll * B1_rr, - B2_ll * B1_rr, - B3_ll * B1_rr, - v_dot_B_ll * B1_rr + v1_ll * psi_ll * psi_rr, - v1_ll * B1_rr, - v2_ll * B1_rr, - v3_ll * B1_rr, - v1_ll * psi_rr) - else # orientation == 2 - f = SVector(0, - B1_ll * B2_rr, - B2_ll * B2_rr, - B3_ll * B2_rr, - v_dot_B_ll * B2_rr + v2_ll * psi_ll * psi_rr, - v1_ll * B2_rr, - v2_ll * B2_rr, - v3_ll * B2_rr, - v2_ll * psi_rr) - end - - return f -end =# @inline function flux_nonconservative_powell(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 @@ -266,22 +230,22 @@ end =# B1_ll * B1_avg, B2_ll * B1_avg, B3_ll * B1_avg, - v_dot_B_ll * B1_avg, # + v1_ll * psi_ll * psi_rr, + v_dot_B_ll * B1_avg + v1_ll * psi_ll * psi_avg, v1_ll * B1_avg, v2_ll * B1_avg, v3_ll * B1_avg, - 0)#v1_ll * psi_avg) + v1_ll * psi_avg) else # orientation == 2 B2_avg = (B2_ll + B2_rr) #* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 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_rr, + v_dot_B_ll * B2_avg + v2_ll * psi_ll * psi_avg, v1_ll * B2_avg, v2_ll * B2_avg, v3_ll * B2_avg, - 0)#v2_ll * psi_avg) + v2_ll * psi_avg) end return f @@ -291,7 +255,8 @@ end """ @inline function flux_nonconservative_powell(u_ll, orientation::Integer, equations::IdealGlmMhdEquations2D, - nonconservative_type::NonConservativeLocal) + nonconservative_type::NonConservativeLocal, + noncons_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 v1_ll = rho_v1_ll / rho_ll @@ -299,18 +264,41 @@ end 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}) - f = SVector(0, - B1_ll, - B2_ll, - B3_ll, - v_dot_B_ll, # The term (v1_ll * psi_ll) is missing because we need to define several non-conservative terms - v1_ll, - v2_ll, - v3_ll, - 0)#v1_ll) - + if noncons_term ==1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + f = SVector(0, + B1_ll, + B2_ll, + B3_ll, + v_dot_B_ll, + v1_ll, + v2_ll, + v3_ll, + 0) + else #noncons_term ==2 + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + if orientation == 1 + f = SVector(0, + 0, + 0, + 0, + v1_ll * psi_ll, + 0, + 0, + 0, + v1_ll) + else #orientation == 2 + f = SVector(0, + 0, + 0, + 0, + v2_ll * psi_ll, + 0, + 0, + 0, + v2_ll) + end + end return f end """ @@ -318,7 +306,8 @@ end """ @inline function flux_nonconservative_powell(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D, - nonconservative_type::NonConservativeSymmetric) + nonconservative_type::NonConservativeSymmetric, + noncons_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 @@ -327,31 +316,43 @@ end 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 # We remove the 0.5 because the flux is always multiplied by 0.5 - if orientation == 1 - B1_avg = (B1_ll + B1_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 - f = SVector(0, - B1_avg, - B1_avg, - B1_avg, - B1_avg, # The term (psi_avg) is missing because we need to define several non-conservative terms - B1_avg, - B1_avg, - B1_avg, - 0)#psi_avg) - else # orientation == 2 - B2_avg = (B2_ll + B2_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5 + if noncons_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 # We remove the 0.5 because the flux is always multiplied by 0.5 + 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 # We remove the 0.5 because the flux is always multiplied by 0.5 + f = SVector(0, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + 0) + end + else #noncons_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 # We remove the 0.5 because the flux is always multiplied by 0.5 f = SVector(0, - B2_avg, - B2_avg, - B2_avg, - B2_avg, # The term (psi_avg) is missing because we need to define several non-conservative terms - B2_avg, - B2_avg, - B2_avg, - 0)#psi_avg) + 0, + 0, + 0, + psi_avg, + 0, + 0, + 0, + psi_avg) end return f diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 7c8bdff4b0e..f658f00fcf1 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -14,6 +14,7 @@ function create_cache(mesh::TreeMesh{2}, equations, A3dp1_x = Array{uEltype, 3} A3dp1_y = Array{uEltype, 3} A3d = Array{uEltype, 3} + A4d = Array{uEltype, 4} fhat1_L_threaded = A3dp1_x[A3dp1_x(undef, nvariables(equations), nnodes(dg) + 1, nnodes(dg)) for _ in 1:Threads.nthreads()] @@ -25,8 +26,8 @@ function create_cache(mesh::TreeMesh{2}, equations, 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()] - flux_temp_nonconservative_threaded = A3d[A3d(undef, nvariables(equations), nnodes(dg), nnodes(dg)) #nnoncons(equations) - for _ in 1:Threads.nthreads()] + flux_temp_nonconservative_threaded = A4d[A4d(undef, nvariables(equations), nnoncons(equations), + nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] antidiffusive_fluxes = Trixi.ContainerAntidiffusiveFlux2D{uEltype}(0, nvariables(equations), nnodes(dg)) @@ -222,12 +223,14 @@ end equations, dg, i, j) multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], flux1, equations, dg, ii, j) - # We multiply by 0.5 because that is done in other parts of Trixi - flux1_noncons = 0.5 * volume_flux_noncons(u_node, u_node_ii, 1, equations, NonConservativeSymmetric()) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[i, ii], flux1_noncons, - equations, dg, i, j) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[ii, i], flux1_noncons, - equations, dg, ii, j) + for noncons in 1:nnoncons(equations) + # We multiply by 0.5 because that is done in other parts of Trixi + flux1_noncons = 0.5 * volume_flux_noncons(u_node, u_node_ii, 1, equations, NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[i, ii], flux1_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[ii, i], flux1_noncons, + equations, dg, noncons, ii, j) + end end end @@ -238,24 +241,31 @@ end fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) fhat_temp = zero(MVector{nvariables(equations), eltype(fhat1_L)}) - fhat_noncons_temp = zero(MVector{nvariables(equations), eltype(fhat1_L)}) + fhat_noncons_temp = zero(MMatrix{nvariables(equations), nnoncons(equations), eltype(fhat1_L)}) for j in eachnode(dg) fhat_temp .= zero(eltype(fhat1_L)) fhat_noncons_temp .= zero(eltype(fhat1_L)) for i in 1:(nnodes(dg) - 1) - # Get the local contribution to the nonconservative flux + # Conservative part + for v in eachvariable(equations) + fhat_temp[v] += weights[i] * flux_temp[v, i, j] + fhat1_L[v, i + 1, j] = fhat_temp[v] + fhat1_R[v, i + 1, j] = fhat_temp[v] + end + # Nonconservative part u_node_L = get_node_vars(u, equations, dg, i, j, element) - phi_L = volume_flux_noncons(u_node_L, 1, equations, NonConservativeLocal()) - u_node_R = get_node_vars(u, equations, dg, i + 1, j, element) - phi_R = volume_flux_noncons(u_node_R, 1, equations, NonConservativeLocal()) - for v in eachvariable(equations) - fhat_temp[v] = fhat_temp[v] + weights[i] * flux_temp[v, i, j] - fhat_noncons_temp[v] = fhat_noncons_temp[v] + weights[i] * flux_temp_noncons[v, i, j] - - fhat1_L[v, i + 1, j] = fhat_temp[v] + phi_L[v] * fhat_noncons_temp[v] - fhat1_R[v, i + 1, j] = fhat_temp[v] + phi_R[v] * fhat_noncons_temp[v] + for noncons in 1:nnoncons(equations) + # Get the local contribution to the nonconservative flux + phi_L = volume_flux_noncons(u_node_L, 1, equations, NonConservativeLocal(), noncons) + phi_R = volume_flux_noncons(u_node_R, 1, equations, NonConservativeLocal(), noncons) + for v in eachvariable(equations) + fhat_noncons_temp[v, noncons] += weights[i] * flux_temp_noncons[v, noncons, i, j] + + fhat1_L[v, i + 1, j] += phi_L[v] * fhat_noncons_temp[v, noncons] + fhat1_R[v, i + 1, j] += phi_R[v] * fhat_noncons_temp[v, noncons] + end end end end @@ -273,12 +283,14 @@ end equations, dg, i, j) multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], flux2, equations, dg, i, jj) - # We multiply by 0.5 because that is done in other parts of Trixi - flux2_noncons = 0.5 * volume_flux_noncons(u_node, u_node_jj, 2, equations, NonConservativeSymmetric()) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[j, jj], flux2_noncons, - equations, dg, i, j) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[jj, j], flux2_noncons, - equations, dg, i, jj) + for noncons in 1:nnoncons(equations) + # We multiply by 0.5 because that is done in other parts of Trixi + flux2_noncons = 0.5 * volume_flux_noncons(u_node, u_node_jj, 2, equations, NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[j, jj], flux2_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[jj, j], flux2_noncons, + equations, dg, noncons, i, jj) + end end end @@ -291,19 +303,26 @@ end for i in eachnode(dg) fhat_temp .= zero(eltype(fhat1_L)) fhat_noncons_temp .= zero(eltype(fhat1_L)) - for j in 1:(nnodes(dg) - 1) - # Get the local contribution to the nonconservative flux + for j in 1:(nnodes(dg) - 1) + # Conservative part + for v in eachvariable(equations) + fhat_temp[v] += weights[j] * flux_temp[v, i, j] + fhat2_L[v, i, j + 1] = fhat_temp[v] + fhat2_R[v, i, j + 1] = fhat_temp[v] + end + # Nonconservative part u_node_L = get_node_vars(u, equations, dg, i, j, element) - phi_L = volume_flux_noncons(u_node_L, 2, equations, NonConservativeLocal()) - u_node_R = get_node_vars(u, equations, dg, i, j + 1, element) - phi_R = volume_flux_noncons(u_node_R, 2, equations, NonConservativeLocal()) - for v in eachvariable(equations) - fhat_temp[v] = fhat_temp[v] + weights[j] * flux_temp[v, i, j] - fhat_noncons_temp[v] = fhat_noncons_temp[v] + weights[j] * flux_temp_noncons[v, i, j] - - fhat2_L[v, i, j + 1] = fhat_temp[v] + phi_L[v] * fhat_noncons_temp[v] - fhat2_R[v, i, j + 1] = fhat_temp[v] + phi_R[v] * fhat_noncons_temp[v] + for noncons in 1:nnoncons(equations) + # Get the local contribution to the nonconservative flux + phi_L = volume_flux_noncons(u_node_L, 2, equations, NonConservativeLocal(), noncons) + phi_R = volume_flux_noncons(u_node_R, 2, equations, NonConservativeLocal(), noncons) + for v in eachvariable(equations) + fhat_noncons_temp[v, noncons] += weights[j] * flux_temp_noncons[v, noncons, i, j] + + fhat2_L[v, i, j + 1] += phi_L[v] * fhat_noncons_temp[v, noncons] + fhat2_R[v, i, j + 1] += phi_R[v] * fhat_noncons_temp[v, noncons] + end end end end From 7017914a1ff465c6d75095397d0e375c8e2e6f1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 10 Oct 2023 12:09:30 +0200 Subject: [PATCH 04/16] Added some comments and improved formatting --- .../elixir_mhd_shockcapturing_subcell.jl | 5 +- src/equations/ideal_glm_mhd_2d.jl | 28 ++++- src/solvers/dgsem_tree/containers_2d.jl | 32 +++--- .../dgsem_tree/dg_2d_subcell_limiters.jl | 108 ++++++++++++------ 4 files changed, 111 insertions(+), 62 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl index db362b74cdb..31855449050 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -10,10 +10,11 @@ equations = IdealGlmMhdEquations2D(1.4) """ initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) -An MHD blast wave taken from +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) @@ -85,7 +86,7 @@ save_solution = SaveSolutionCallback(interval=100, save_final_solution=true, solution_variables=cons2prim) -cfl = 0.3 +cfl = 0.5 stepsize_callback = StepsizeCallback(cfl=cfl) glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 9a26503397d..f9290539e72 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -251,10 +251,18 @@ terms. return f end """ - + flux_nonconservative_powell(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeLocal, + noncons_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. """ @inline function flux_nonconservative_powell(u_ll, orientation::Integer, - equations::IdealGlmMhdEquations2D, + equations::IdealGlmMhdEquations2D, nonconservative_type::NonConservativeLocal, noncons_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 @@ -264,7 +272,7 @@ end v3_ll = rho_v3_ll / rho_ll v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll - if noncons_term ==1 + if noncons_term == 1 # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) f = SVector(0, B1_ll, @@ -302,10 +310,18 @@ end return f end """ - + flux_nonconservative_powell(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeSymmetric, + noncons_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. """ @inline function flux_nonconservative_powell(u_ll, u_rr, orientation::Integer, - equations::IdealGlmMhdEquations2D, + equations::IdealGlmMhdEquations2D, nonconservative_type::NonConservativeSymmetric, noncons_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 @@ -316,7 +332,7 @@ end v3_ll = rho_v3_ll / rho_ll v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll - if noncons_term ==1 + if noncons_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 # We remove the 0.5 because the flux is always multiplied by 0.5 diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl index 576732b9219..0784a834033 100644 --- a/src/solvers/dgsem_tree/containers_2d.jl +++ b/src/solvers/dgsem_tree/containers_2d.jl @@ -1283,22 +1283,22 @@ function ContainerAntidiffusiveFlux2D{uEltype}(capacity::Integer, n_variables, # Initialize fields with defaults _antidiffusive_flux1_L = fill(nan_uEltype, - n_variables * (n_nodes + 1) * n_nodes * capacity) + 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)) + (n_variables, n_nodes + 1, n_nodes, capacity)) _antidiffusive_flux1_R = fill(nan_uEltype, - n_variables * (n_nodes + 1) * n_nodes * capacity) + 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)) + (n_variables, n_nodes + 1, n_nodes, capacity)) _antidiffusive_flux2_L = fill(nan_uEltype, - n_variables * n_nodes * (n_nodes + 1) * capacity) + 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)) + (n_variables, n_nodes, n_nodes + 1, capacity)) _antidiffusive_flux2_R = fill(nan_uEltype, - n_variables * n_nodes * (n_nodes + 1) * capacity) + 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)) + (n_variables, n_nodes, n_nodes + 1, capacity)) return ContainerAntidiffusiveFlux2D{uEltype}(antidiffusive_flux1_L, antidiffusive_flux1_R, @@ -1326,20 +1326,20 @@ function Base.resize!(fluxes::ContainerAntidiffusiveFlux2D, capacity) 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)) + (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)) + (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)) + (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)) + (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 f658f00fcf1..4af22c2f0e3 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -17,21 +17,23 @@ function create_cache(mesh::TreeMesh{2}, equations, A4d = Array{uEltype, 4} fhat1_L_threaded = A3dp1_x[A3dp1_x(undef, nvariables(equations), nnodes(dg) + 1, - nnodes(dg)) for _ in 1:Threads.nthreads()] + 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()] + 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()] + 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()] + 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()] - flux_temp_nonconservative_threaded = A4d[A4d(undef, nvariables(equations), nnoncons(equations), - nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + flux_temp_nonconservative_threaded = A4d[A4d(undef, nvariables(equations), + nnoncons(equations), + nnodes(dg), nnodes(dg)) + for _ in 1:Threads.nthreads()] antidiffusive_fluxes = Trixi.ContainerAntidiffusiveFlux2D{uEltype}(0, nvariables(equations), nnodes(dg)) - return (; cache..., antidiffusive_fluxes, + return (; cache..., antidiffusive_fluxes, fhat1_L_threaded, fhat2_L_threaded, fhat1_R_threaded, fhat2_R_threaded, flux_temp_threaded, flux_temp_nonconservative_threaded) end @@ -80,8 +82,8 @@ end nonconservative_terms, equations, volume_flux_fv, dg, element, cache) # antidiffusive flux - calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, - fstar1_L, fstar1_R, fstar2_L, fstar2_R, + 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) @@ -111,7 +113,7 @@ end @unpack flux_temp_threaded = cache flux_temp = flux_temp_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, @@ -184,6 +186,13 @@ end # (**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, @@ -225,11 +234,15 @@ end equations, dg, ii, j) for noncons in 1:nnoncons(equations) # We multiply by 0.5 because that is done in other parts of Trixi - flux1_noncons = 0.5 * volume_flux_noncons(u_node, u_node_ii, 1, equations, NonConservativeSymmetric(), noncons) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[i, ii], flux1_noncons, - equations, dg, noncons, i, j) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[ii, i], flux1_noncons, - equations, dg, noncons, ii, j) + flux1_noncons = 0.5 * + volume_flux_noncons(u_node, u_node_ii, 1, equations, + NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[i, ii], + flux1_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[ii, i], + flux1_noncons, + equations, dg, noncons, ii, j) end end end @@ -241,7 +254,8 @@ end fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) fhat_temp = zero(MVector{nvariables(equations), eltype(fhat1_L)}) - fhat_noncons_temp = zero(MMatrix{nvariables(equations), nnoncons(equations), eltype(fhat1_L)}) + fhat_noncons_temp = zero(MMatrix{nvariables(equations), nnoncons(equations), + eltype(fhat1_L)}) for j in eachnode(dg) fhat_temp .= zero(eltype(fhat1_L)) @@ -258,11 +272,14 @@ end u_node_R = get_node_vars(u, equations, dg, i + 1, j, element) for noncons in 1:nnoncons(equations) # Get the local contribution to the nonconservative flux - phi_L = volume_flux_noncons(u_node_L, 1, equations, NonConservativeLocal(), noncons) - phi_R = volume_flux_noncons(u_node_R, 1, equations, NonConservativeLocal(), noncons) + phi_L = volume_flux_noncons(u_node_L, 1, equations, + NonConservativeLocal(), noncons) + phi_R = volume_flux_noncons(u_node_R, 1, equations, + NonConservativeLocal(), noncons) for v in eachvariable(equations) - fhat_noncons_temp[v, noncons] += weights[i] * flux_temp_noncons[v, noncons, i, j] - + fhat_noncons_temp[v, noncons] += weights[i] * + flux_temp_noncons[v, noncons, i, j] + fhat1_L[v, i + 1, j] += phi_L[v] * fhat_noncons_temp[v, noncons] fhat1_R[v, i + 1, j] += phi_R[v] * fhat_noncons_temp[v, noncons] end @@ -285,11 +302,15 @@ end equations, dg, i, jj) for noncons in 1:nnoncons(equations) # We multiply by 0.5 because that is done in other parts of Trixi - flux2_noncons = 0.5 * volume_flux_noncons(u_node, u_node_jj, 2, equations, NonConservativeSymmetric(), noncons) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[j, jj], flux2_noncons, - equations, dg, noncons, i, j) - multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[jj, j], flux2_noncons, - equations, dg, noncons, i, jj) + flux2_noncons = 0.5 * + volume_flux_noncons(u_node, u_node_jj, 2, equations, + NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[j, jj], + flux2_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_temp_noncons, derivative_split[jj, j], + flux2_noncons, + equations, dg, noncons, i, jj) end end end @@ -315,10 +336,13 @@ end u_node_R = get_node_vars(u, equations, dg, i, j + 1, element) for noncons in 1:nnoncons(equations) # Get the local contribution to the nonconservative flux - phi_L = volume_flux_noncons(u_node_L, 2, equations, NonConservativeLocal(), noncons) - phi_R = volume_flux_noncons(u_node_R, 2, equations, NonConservativeLocal(), noncons) + phi_L = volume_flux_noncons(u_node_L, 2, equations, + NonConservativeLocal(), noncons) + phi_R = volume_flux_noncons(u_node_R, 2, equations, + NonConservativeLocal(), noncons) for v in eachvariable(equations) - fhat_noncons_temp[v, noncons] += weights[j] * flux_temp_noncons[v, noncons, i, j] + fhat_noncons_temp[v, noncons] += weights[j] * + flux_temp_noncons[v, noncons, i, j] fhat2_L[v, i, j + 1] += phi_L[v] * fhat_noncons_temp[v, noncons] fhat2_R[v, i, j + 1] += phi_R[v] * fhat_noncons_temp[v, noncons] @@ -331,8 +355,8 @@ end 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, +@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) @@ -340,14 +364,18 @@ end 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] + 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] + 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 @@ -373,14 +401,18 @@ end 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] = fhat1_R[v, i, j] - fstar1_R[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_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] + 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 From 5119d97c85bea9e0b89c49cda6d32c5365663768 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 10 Oct 2023 12:25:30 +0200 Subject: [PATCH 05/16] Added expected results for MHD subcell limiting test --- test/test_tree_2d_mhd.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index 3e104da3e91..da9d34ef74a 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -73,6 +73,12 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") # of this IC to have negative density/pressure values, crashing the simulation. coverage_override = (maxiters=9,)) end + + @trixi_testset "elixir_mhd_shockcapturing_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"), + l2 = [1.72786992e-01, 6.33650587e-02, 6.86872255e-02, 0.00000000e+00, 4.31337885e+00, 1.67036008e-01, 1.06316839e-01, 0.00000000e+00, 4.67098356e-03], + linf = [9.80256401e-01, 9.20713091e-01, 5.55986508e-01, 0.00000000e+00, 2.49132899e+01, 6.39041960e-01, 6.08144670e-01, 0.00000000e+00, 5.83546136e-02]) + end end end # module From ad6177e64cda4ff04b457757ee43aef48565662b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 10 Oct 2023 16:52:03 +0200 Subject: [PATCH 06/16] Restored old Powell source term and created a new function name for the modified Powell source term. This was needed due to incompatibility on non-conforming meshes. --- .../elixir_mhd_shockcapturing_subcell.jl | 4 +- src/Trixi.jl | 2 +- src/equations/ideal_glm_mhd_2d.jl | 141 ++++++++++++------ 3 files changed, 100 insertions(+), 47 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl index 31855449050..d7ef23332fe 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl @@ -46,8 +46,8 @@ function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D) end initial_condition = initial_condition_blast_wave -surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell) -volume_flux = (flux_derigs_etal, flux_nonconservative_powell) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell2) +volume_flux = (flux_derigs_etal, flux_nonconservative_powell2) basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; diff --git a/src/Trixi.jl b/src/Trixi.jl index b4485217411..1346d8ff7af 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -161,7 +161,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_powell2 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/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index f9290539e72..76cfc7b4335 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -221,6 +221,95 @@ terms. 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}) + if orientation == 1 + f = SVector(0, + B1_ll * B1_rr, + B2_ll * B1_rr, + B3_ll * B1_rr, + v_dot_B_ll * B1_rr + v1_ll * psi_ll * psi_rr, + v1_ll * B1_rr, + v2_ll * B1_rr, + v3_ll * B1_B1_rravg, + v1_ll * psi_rr) + else # orientation == 2 + f = SVector(0, + B1_ll * B2_rr, + B2_ll * B2_rr, + B3_ll * B2_rr, + v_dot_B_ll * B2_rr + v2_ll * psi_ll * psi_rr, + v1_ll * B2_rr, + v2_ll * B2_rr, + v3_ll * B2_rr, + v2_ll * psi_rr) + end + + return f +end + +@inline function flux_nonconservative_powell(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + 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 + + # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector + # at the same node location) while `B_dot_n_rr` uses the averaged normal + # direction. The reason for this is that `v_dot_n_ll` depends only on the left + # state and multiplies some gradient while `B_dot_n_rr` is used to compute + # the divergence of B. + v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] + B_dot_n_rr = B1_rr * normal_direction_average[1] + + B2_rr * normal_direction_average[2] + + # 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}) + f = SVector(0, + B1_ll * B_dot_n_rr, + B2_ll * B_dot_n_rr, + B3_ll * B_dot_n_rr, + v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr, + v1_ll * B_dot_n_rr, + v2_ll * B_dot_n_rr, + v3_ll * B_dot_n_rr, + v_dot_n_ll * psi_rr) + + return f +end + +""" + flux_nonconservative_powell2(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(!). +## 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_powell2(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 # We remove the 0.5 because the flux is always multiplied by 0.5 @@ -261,10 +350,10 @@ 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. """ -@inline function flux_nonconservative_powell(u_ll, orientation::Integer, - equations::IdealGlmMhdEquations2D, - nonconservative_type::NonConservativeLocal, - noncons_term::Integer) +@inline function flux_nonconservative_powell2(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeLocal, + noncons_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 v1_ll = rho_v1_ll / rho_ll @@ -320,10 +409,10 @@ 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. """ -@inline function flux_nonconservative_powell(u_ll, u_rr, orientation::Integer, - equations::IdealGlmMhdEquations2D, - nonconservative_type::NonConservativeSymmetric, - noncons_term::Integer) +@inline function flux_nonconservative_powell2(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeSymmetric, + noncons_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 @@ -374,42 +463,6 @@ the non-conservative staggered "fluxes" for subcell limiting. See, e.g., return f end -@inline function flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - 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 - - # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector - # at the same node location) while `B_dot_n_rr` uses the averaged normal - # direction. The reason for this is that `v_dot_n_ll` depends only on the left - # state and multiplies some gradient while `B_dot_n_rr` is used to compute - # the divergence of B. - v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] - B_dot_n_rr = B1_rr * normal_direction_average[1] + - B2_rr * normal_direction_average[2] - - # 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}) - f = SVector(0, - B1_ll * B_dot_n_rr, - B2_ll * B_dot_n_rr, - B3_ll * B_dot_n_rr, - v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr, - v1_ll * B_dot_n_rr, - v2_ll * B_dot_n_rr, - v3_ll * B_dot_n_rr, - v_dot_n_ll * psi_rr) - - return f -end - """ flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations2D) From 91580938ef79ce9783d4bfb17aff034ecb04d537 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 10 Oct 2023 16:54:00 +0200 Subject: [PATCH 07/16] Fixed bug --- src/equations/ideal_glm_mhd_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 76cfc7b4335..fcd337d6467 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -231,7 +231,7 @@ terms. v_dot_B_ll * B1_rr + v1_ll * psi_ll * psi_rr, v1_ll * B1_rr, v2_ll * B1_rr, - v3_ll * B1_B1_rravg, + v3_ll * B1_rr, v1_ll * psi_rr) else # orientation == 2 f = SVector(0, From 3dcccc40d97857cb29a1d96aca97b3277409824c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 10 Oct 2023 16:58:53 +0200 Subject: [PATCH 08/16] Fixed bug --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 1346d8ff7af..14687cbb992 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -161,7 +161,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_powell2 + flux_nonconservative_powell, flux_nonconservative_powell2, 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, From adb930e08f2bf5902f24ee20c3cefd8e54a28f44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 11 Oct 2023 11:30:20 +0200 Subject: [PATCH 09/16] Moved new multiple-dispatch structs for to equations.jl --- src/equations/equations.jl | 13 +++++++++++++ src/equations/ideal_glm_mhd_2d.jl | 3 --- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a941f750a68..ac232a8bb32 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -207,7 +207,20 @@ where `x` specifies the coordinates, `t` is the current time, and `equation` is 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) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index fcd337d6467..09b1788ed8f 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -23,9 +23,6 @@ mutable struct IdealGlmMhdEquations2D{RealT <: Real} <: end end -struct NonConservativeLocal end -struct NonConservativeSymmetric end - function IdealGlmMhdEquations2D(gamma; initial_c_h = convert(typeof(gamma), NaN)) # Use `promote` to ensure that `gamma` and `initial_c_h` have the same type IdealGlmMhdEquations2D(promote(gamma, initial_c_h)...) From b0e3aad73c8162287f5049e7be000ee9c742418b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 11 Oct 2023 12:48:41 +0200 Subject: [PATCH 10/16] Improved allocations --- .../dgsem_tree/dg_2d_subcell_limiters.jl | 80 ++++++++++++------- 1 file changed, 52 insertions(+), 28 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 4af22c2f0e3..adcb06837bf 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -30,12 +30,20 @@ function create_cache(mesh::TreeMesh{2}, equations, nnoncons(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()] + fhat_nonconservative_temp_threaded = A4d[A4d(undef, nvariables(equations), + nnoncons(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_L_threaded, fhat2_L_threaded, fhat1_R_threaded, fhat2_R_threaded, - flux_temp_threaded, flux_temp_nonconservative_threaded) + flux_temp_threaded, flux_temp_nonconservative_threaded, fhat_temp_threaded, + fhat_nonconservative_temp_threaded) end function calc_volume_integral!(du, u, @@ -198,13 +206,16 @@ end equations, volume_flux, dg::DGSEM, element, cache) @unpack weights, derivative_split = dg.basis - @unpack flux_temp_threaded, flux_temp_nonconservative_threaded = cache + @unpack flux_temp_threaded, flux_temp_nonconservative_threaded, fhat_temp_threaded, fhat_nonconservative_temp_threaded = cache volume_flux_cons, volume_flux_noncons = volume_flux flux_temp = flux_temp_threaded[Threads.threadid()] flux_temp_noncons = flux_temp_nonconservative_threaded[Threads.threadid()] + fhat_temp = fhat_temp_threaded[Threads.threadid()] + fhat_noncons_temp = fhat_nonconservative_temp_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, @@ -253,19 +264,17 @@ end fhat1_R[:, 1, :] .= zero(eltype(fhat1_R)) fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) - fhat_temp = zero(MVector{nvariables(equations), eltype(fhat1_L)}) - fhat_noncons_temp = zero(MMatrix{nvariables(equations), nnoncons(equations), - eltype(fhat1_L)}) + fhat_temp .= zero(eltype(fhat1_L)) + fhat_noncons_temp .= zero(eltype(fhat1_L)) - for j in eachnode(dg) - fhat_temp .= zero(eltype(fhat1_L)) - fhat_noncons_temp .= zero(eltype(fhat1_L)) + @trixi_timeit timer() "subcell volume integral (x)" for j in eachnode(dg) for i in 1:(nnodes(dg) - 1) # Conservative part for v in eachvariable(equations) - fhat_temp[v] += weights[i] * flux_temp[v, i, j] - fhat1_L[v, i + 1, j] = fhat_temp[v] - fhat1_R[v, i + 1, j] = fhat_temp[v] + fhat_temp[v, i + 1, j] = fhat_temp[v, i, j] + + weights[i] * flux_temp[v, i, j] + fhat1_L[v, i + 1, j] = fhat_temp[v, i + 1, j] + fhat1_R[v, i + 1, j] = fhat_temp[v, i + 1, j] end # Nonconservative part u_node_L = get_node_vars(u, equations, dg, i, j, element) @@ -277,11 +286,18 @@ end phi_R = volume_flux_noncons(u_node_R, 1, equations, NonConservativeLocal(), noncons) for v in eachvariable(equations) - fhat_noncons_temp[v, noncons] += weights[i] * - flux_temp_noncons[v, noncons, i, j] - - fhat1_L[v, i + 1, j] += phi_L[v] * fhat_noncons_temp[v, noncons] - fhat1_R[v, i + 1, j] += phi_R[v] * fhat_noncons_temp[v, noncons] + fhat_noncons_temp[v, noncons, i + 1, j] = fhat_noncons_temp[v, + noncons, + i, j] + + weights[i] * + flux_temp_noncons[v, + noncons, + i, j] + + fhat1_L[v, i + 1, j] += phi_L[v] * + fhat_noncons_temp[v, noncons, i + 1, j] + fhat1_R[v, i + 1, j] += phi_R[v] * + fhat_noncons_temp[v, noncons, i + 1, j] end end end @@ -321,15 +337,17 @@ end fhat2_R[:, :, 1] .= zero(eltype(fhat2_R)) fhat2_R[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_R)) - for i in eachnode(dg) - fhat_temp .= zero(eltype(fhat1_L)) - fhat_noncons_temp .= zero(eltype(fhat1_L)) + fhat_temp .= zero(eltype(fhat1_L)) + fhat_noncons_temp .= zero(eltype(fhat1_L)) + + @trixi_timeit timer() "subcell volume integral (y)" for i in eachnode(dg) for j in 1:(nnodes(dg) - 1) # Conservative part for v in eachvariable(equations) - fhat_temp[v] += weights[j] * flux_temp[v, i, j] - fhat2_L[v, i, j + 1] = fhat_temp[v] - fhat2_R[v, i, j + 1] = fhat_temp[v] + fhat_temp[v, i, j + 1] = fhat_temp[v, i, j] + + weights[j] * flux_temp[v, i, j] + fhat2_L[v, i, j + 1] = fhat_temp[v, i, j + 1] + fhat2_R[v, i, j + 1] = fhat_temp[v, i, j + 1] end # Nonconservative part u_node_L = get_node_vars(u, equations, dg, i, j, element) @@ -341,16 +359,22 @@ end phi_R = volume_flux_noncons(u_node_R, 2, equations, NonConservativeLocal(), noncons) for v in eachvariable(equations) - fhat_noncons_temp[v, noncons] += weights[j] * - flux_temp_noncons[v, noncons, i, j] - - fhat2_L[v, i, j + 1] += phi_L[v] * fhat_noncons_temp[v, noncons] - fhat2_R[v, i, j + 1] += phi_R[v] * fhat_noncons_temp[v, noncons] + fhat_noncons_temp[v, noncons, i, j + 1] = fhat_noncons_temp[v, + noncons, + i, j] + + weights[j] * + flux_temp_noncons[v, + noncons, + i, j] + + fhat2_L[v, i, j + 1] += phi_L[v] * + fhat_noncons_temp[v, noncons, i, j + 1] + fhat2_R[v, i, j + 1] += phi_R[v] * + fhat_noncons_temp[v, noncons, i, j + 1] end end end end - return nothing end From e0a3fdfcb577683ebcf57eafce65539f67cab8d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 11 Oct 2023 16:55:30 +0200 Subject: [PATCH 11/16] Avoid double computation of local part of non-conservative flux --- .../dgsem_tree/dg_2d_subcell_limiters.jl | 52 ++++++++++++------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index adcb06837bf..da5b438e424 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -37,13 +37,19 @@ function create_cache(mesh::TreeMesh{2}, equations, nnoncons(equations), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + + phi_threaded = A4d[A4d(undef, nvariables(equations), + nnoncons(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_L_threaded, fhat2_L_threaded, fhat1_R_threaded, fhat2_R_threaded, flux_temp_threaded, flux_temp_nonconservative_threaded, fhat_temp_threaded, - fhat_nonconservative_temp_threaded) + fhat_nonconservative_temp_threaded, phi_threaded) end function calc_volume_integral!(du, u, @@ -206,7 +212,8 @@ end equations, volume_flux, dg::DGSEM, element, cache) @unpack weights, derivative_split = dg.basis - @unpack flux_temp_threaded, flux_temp_nonconservative_threaded, fhat_temp_threaded, fhat_nonconservative_temp_threaded = cache + @unpack flux_temp_threaded, flux_temp_nonconservative_threaded = cache + @unpack fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded = cache volume_flux_cons, volume_flux_noncons = volume_flux @@ -215,6 +222,7 @@ end 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, @@ -267,6 +275,15 @@ end fhat_temp .= zero(eltype(fhat1_L)) fhat_noncons_temp .= zero(eltype(fhat1_L)) + # Compute local contribution to non-conservative flux + for j in eachnode(dg), i in eachnode(dg), noncons in 1:nnoncons(equations) + u_local = get_node_vars(u, equations, dg, i, j, element) + set_node_vars!(phi, + volume_flux_noncons(u_local, 1, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j) + end + @trixi_timeit timer() "subcell volume integral (x)" for j in eachnode(dg) for i in 1:(nnodes(dg) - 1) # Conservative part @@ -277,14 +294,7 @@ end fhat1_R[v, i + 1, j] = fhat_temp[v, i + 1, j] end # Nonconservative part - u_node_L = get_node_vars(u, equations, dg, i, j, element) - u_node_R = get_node_vars(u, equations, dg, i + 1, j, element) for noncons in 1:nnoncons(equations) - # Get the local contribution to the nonconservative flux - phi_L = volume_flux_noncons(u_node_L, 1, equations, - NonConservativeLocal(), noncons) - phi_R = volume_flux_noncons(u_node_R, 1, equations, - NonConservativeLocal(), noncons) for v in eachvariable(equations) fhat_noncons_temp[v, noncons, i + 1, j] = fhat_noncons_temp[v, noncons, @@ -294,9 +304,9 @@ end noncons, i, j] - fhat1_L[v, i + 1, j] += phi_L[v] * + fhat1_L[v, i + 1, j] += phi[v, noncons, i, j] * fhat_noncons_temp[v, noncons, i + 1, j] - fhat1_R[v, i + 1, j] += phi_R[v] * + fhat1_R[v, i + 1, j] += phi[v, noncons, i + 1, j] * fhat_noncons_temp[v, noncons, i + 1, j] end end @@ -340,6 +350,15 @@ end fhat_temp .= zero(eltype(fhat1_L)) fhat_noncons_temp .= zero(eltype(fhat1_L)) + # Compute local contribution to non-conservative flux + for j in eachnode(dg), i in eachnode(dg), noncons in 1:nnoncons(equations) + u_local = get_node_vars(u, equations, dg, i, j, element) + set_node_vars!(phi, + volume_flux_noncons(u_local, 2, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j) + end + @trixi_timeit timer() "subcell volume integral (y)" for i in eachnode(dg) for j in 1:(nnodes(dg) - 1) # Conservative part @@ -350,14 +369,7 @@ end fhat2_R[v, i, j + 1] = fhat_temp[v, i, j + 1] end # Nonconservative part - u_node_L = get_node_vars(u, equations, dg, i, j, element) - u_node_R = get_node_vars(u, equations, dg, i, j + 1, element) for noncons in 1:nnoncons(equations) - # Get the local contribution to the nonconservative flux - phi_L = volume_flux_noncons(u_node_L, 2, equations, - NonConservativeLocal(), noncons) - phi_R = volume_flux_noncons(u_node_R, 2, equations, - NonConservativeLocal(), noncons) for v in eachvariable(equations) fhat_noncons_temp[v, noncons, i, j + 1] = fhat_noncons_temp[v, noncons, @@ -367,9 +379,9 @@ end noncons, i, j] - fhat2_L[v, i, j + 1] += phi_L[v] * + fhat2_L[v, i, j + 1] += phi[v, noncons, i, j] * fhat_noncons_temp[v, noncons, i, j + 1] - fhat2_R[v, i, j + 1] += phi_R[v] * + fhat2_R[v, i, j + 1] += phi[v, noncons, i, j + 1] * fhat_noncons_temp[v, noncons, i, j + 1] end end From 86884e81e7cfbe18d6f7390fc64d8d14f21bd574 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Wed, 11 Oct 2023 17:15:20 +0200 Subject: [PATCH 12/16] Improved subcell volume integral in y direction and formatting --- .../dgsem_tree/dg_2d_subcell_limiters.jl | 100 ++++++++---------- 1 file changed, 44 insertions(+), 56 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index da5b438e424..3a9d898bc41 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -272,8 +272,8 @@ end fhat1_R[:, 1, :] .= zero(eltype(fhat1_R)) fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) - fhat_temp .= zero(eltype(fhat1_L)) - fhat_noncons_temp .= zero(eltype(fhat1_L)) + 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), noncons in 1:nnoncons(equations) @@ -284,32 +284,26 @@ end equations, dg, noncons, i, j) end - @trixi_timeit timer() "subcell volume integral (x)" for j in eachnode(dg) - for i in 1:(nnodes(dg) - 1) - # Conservative part - for v in eachvariable(equations) - fhat_temp[v, i + 1, j] = fhat_temp[v, i, j] + - weights[i] * flux_temp[v, i, j] - fhat1_L[v, i + 1, j] = fhat_temp[v, i + 1, j] - fhat1_R[v, i + 1, j] = fhat_temp[v, i + 1, j] - end - # Nonconservative part - for noncons in 1:nnoncons(equations) - for v in eachvariable(equations) - fhat_noncons_temp[v, noncons, i + 1, j] = fhat_noncons_temp[v, - noncons, - i, j] + - weights[i] * - flux_temp_noncons[v, - noncons, - i, j] - - fhat1_L[v, i + 1, j] += phi[v, noncons, i, j] * - fhat_noncons_temp[v, noncons, i + 1, j] - fhat1_R[v, i + 1, j] += phi[v, noncons, i + 1, j] * - fhat_noncons_temp[v, noncons, i + 1, j] - end - end + for j in eachnode(dg), i in 1:(nnodes(dg) - 1) + # Conservative part + for v in eachvariable(equations) + fhat_temp[v, i + 1, j] = fhat_temp[v, i, j] + + weights[i] * flux_temp[v, i, j] + fhat1_L[v, i + 1, j] = fhat_temp[v, i + 1, j] + fhat1_R[v, i + 1, j] = fhat_temp[v, i + 1, j] + end + # Nonconservative part + for noncons in 1:nnoncons(equations), v in eachvariable(equations) + fhat_noncons_temp[v, noncons, i + 1, j] = fhat_noncons_temp[v, noncons, i, + j] + + weights[i] * + flux_temp_noncons[v, noncons, i, + j] + + fhat1_L[v, i + 1, j] += phi[v, noncons, i, j] * + fhat_noncons_temp[v, noncons, i + 1, j] + fhat1_R[v, i + 1, j] += phi[v, noncons, i + 1, j] * + fhat_noncons_temp[v, noncons, i + 1, j] end end @@ -347,8 +341,8 @@ end fhat2_R[:, :, 1] .= zero(eltype(fhat2_R)) fhat2_R[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_R)) - fhat_temp .= zero(eltype(fhat1_L)) - fhat_noncons_temp .= zero(eltype(fhat1_L)) + 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), noncons in 1:nnoncons(equations) @@ -359,32 +353,26 @@ end equations, dg, noncons, i, j) end - @trixi_timeit timer() "subcell volume integral (y)" for i in eachnode(dg) - for j in 1:(nnodes(dg) - 1) - # Conservative part - for v in eachvariable(equations) - fhat_temp[v, i, j + 1] = fhat_temp[v, i, j] + - weights[j] * flux_temp[v, i, j] - fhat2_L[v, i, j + 1] = fhat_temp[v, i, j + 1] - fhat2_R[v, i, j + 1] = fhat_temp[v, i, j + 1] - end - # Nonconservative part - for noncons in 1:nnoncons(equations) - for v in eachvariable(equations) - fhat_noncons_temp[v, noncons, i, j + 1] = fhat_noncons_temp[v, - noncons, - i, j] + - weights[j] * - flux_temp_noncons[v, - noncons, - i, j] - - fhat2_L[v, i, j + 1] += phi[v, noncons, i, j] * - fhat_noncons_temp[v, noncons, i, j + 1] - fhat2_R[v, i, j + 1] += phi[v, noncons, i, j + 1] * - fhat_noncons_temp[v, noncons, i, j + 1] - end - end + for j in 1:(nnodes(dg) - 1), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + fhat_temp[v, i, j + 1] = fhat_temp[v, i, j] + + weights[j] * flux_temp[v, i, j] + fhat2_L[v, i, j + 1] = fhat_temp[v, i, j + 1] + fhat2_R[v, i, j + 1] = fhat_temp[v, i, j + 1] + end + # Nonconservative part + for noncons in 1:nnoncons(equations), v in eachvariable(equations) + fhat_noncons_temp[v, noncons, i, j + 1] = fhat_noncons_temp[v, noncons, i, + j] + + weights[j] * + flux_temp_noncons[v, noncons, i, + j] + + fhat2_L[v, i, j + 1] += phi[v, noncons, i, j] * + fhat_noncons_temp[v, noncons, i, j + 1] + fhat2_R[v, i, j + 1] += phi[v, noncons, i, j + 1] * + fhat_noncons_temp[v, noncons, i, j + 1] end end return nothing From 0d61cfdb64e7eb97046f87898ff6c0beecb7b51c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 12 Oct 2023 10:22:10 +0200 Subject: [PATCH 13/16] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/equations/equations.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index ac232a8bb32..63151f5dff7 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -208,17 +208,19 @@ struct BoundaryConditionNeumann{B} boundary_normal_flux_function::B end """ - NonConservativeLocal + 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. +When the argument `nonconservative_type` is of type `NonConservativeLocal`, +the function returns the local part of the non-conservative term. """ struct NonConservativeLocal end """ - NonConservativeSymmetric + 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. +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 From b0caf8eedd81ee1fe3558535560f1e08c365a45e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 12 Oct 2023 10:25:10 +0200 Subject: [PATCH 14/16] Added timers and corrected docstrings --- src/equations/ideal_glm_mhd_2d.jl | 4 ++-- src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 09b1788ed8f..45aeba3ac2a 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -337,7 +337,7 @@ results on non-conforming meshes(!). return f end """ - flux_nonconservative_powell(u_ll, orientation::Integer, + flux_nonconservative_powell2(u_ll, orientation::Integer, equations::IdealGlmMhdEquations2D, nonconservative_type::NonConservativeLocal, noncons_term::Integer) @@ -396,7 +396,7 @@ the non-conservative staggered "fluxes" for subcell limiting. See, e.g., return f end """ - flux_nonconservative_powell(u_ll, orientation::Integer, + flux_nonconservative_powell2(u_ll, orientation::Integer, equations::IdealGlmMhdEquations2D, nonconservative_type::NonConservativeSymmetric, noncons_term::Integer) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 3a9d898bc41..17354cd27ba 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -82,7 +82,7 @@ end 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, + @trixi_timeit timer() "calcflux_fhat!" 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 @@ -92,11 +92,11 @@ end fstar2_L = fstar2_L_threaded[Threads.threadid()] 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, + @trixi_timeit timer() "calcflux_fv!" calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, mesh, nonconservative_terms, equations, volume_flux_fv, dg, element, cache) # antidiffusive flux - calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, + @trixi_timeit timer() "calcflux_antidiffusive!" 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, @@ -214,7 +214,7 @@ end @unpack weights, derivative_split = dg.basis @unpack flux_temp_threaded, flux_temp_nonconservative_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()] From 6dae80a214fbbae1f0506043f9caadc60557f6cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 12 Oct 2023 10:32:09 +0200 Subject: [PATCH 15/16] Improved formatting --- .../dgsem_tree/dg_2d_subcell_limiters.jl | 28 +++++++++++-------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 17354cd27ba..3731f322eee 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -82,9 +82,11 @@ end fhat1_R = fhat1_R_threaded[Threads.threadid()] fhat2_L = fhat2_L_threaded[Threads.threadid()] fhat2_R = fhat2_R_threaded[Threads.threadid()] - @trixi_timeit timer() "calcflux_fhat!" calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, mesh, - nonconservative_terms, equations, volume_flux_dg, dg, element, cache) - + @trixi_timeit timer() "calcflux_fhat!" begin + calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, mesh, + nonconservative_terms, equations, volume_flux_dg, dg, element, + cache) + end # low-order FV fluxes @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded = cache @@ -92,15 +94,19 @@ end fstar2_L = fstar2_L_threaded[Threads.threadid()] fstar1_R = fstar1_R_threaded[Threads.threadid()] fstar2_R = fstar2_R_threaded[Threads.threadid()] - @trixi_timeit timer() "calcflux_fv!" calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, mesh, - nonconservative_terms, equations, volume_flux_fv, dg, element, cache) + @trixi_timeit timer() "calcflux_fv!" begin + calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, mesh, + nonconservative_terms, equations, volume_flux_fv, dg, element, + cache) + end # antidiffusive flux - @trixi_timeit timer() "calcflux_antidiffusive!" 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) + @trixi_timeit timer() "calcflux_antidiffusive!" begin + 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) + end # Calculate volume integral contribution of low-order FV flux for j in eachnode(dg), i in eachnode(dg) @@ -214,7 +220,7 @@ end @unpack weights, derivative_split = dg.basis @unpack flux_temp_threaded, flux_temp_nonconservative_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()] From 0c2e24e4d0e7f4a36a6cc412e666301f3fc7bd75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 12 Oct 2023 10:50:59 +0200 Subject: [PATCH 16/16] Reduced testing time --- test/test_tree_2d_mhd.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index da9d34ef74a..c3911a741c4 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -76,8 +76,9 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") @trixi_testset "elixir_mhd_shockcapturing_subcell.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"), - l2 = [1.72786992e-01, 6.33650587e-02, 6.86872255e-02, 0.00000000e+00, 4.31337885e+00, 1.67036008e-01, 1.06316839e-01, 0.00000000e+00, 4.67098356e-03], - linf = [9.80256401e-01, 9.20713091e-01, 5.55986508e-01, 0.00000000e+00, 2.49132899e+01, 6.39041960e-01, 6.08144670e-01, 0.00000000e+00, 5.83546136e-02]) + 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)) end end