diff --git a/examples/tree_2d_dgsem/convergence_subcell/elixir_mhd_alfven_wave_fluxdifferencing.jl b/examples/tree_2d_dgsem/convergence_subcell/elixir_mhd_alfven_wave_fluxdifferencing.jl new file mode 100644 index 00000000000..bb2df4f0fc6 --- /dev/null +++ b/examples/tree_2d_dgsem/convergence_subcell/elixir_mhd_alfven_wave_fluxdifferencing.jl @@ -0,0 +1,71 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible ideal GLM-MHD equations +gamma = 5 / 3 +equations = IdealGlmMhdEquations2D(gamma) + +initial_condition = initial_condition_convergence_test + +volume_flux = (flux_central, flux_nonconservative_powell2) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell2) + +basis = LobattoLegendreBasis(3) +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(polydeg = 3, + surface_flux = surface_flux, + volume_integral = volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.0)) +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, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (entropy, energy_total, + energy_kinetic, + energy_internal, + energy_magnetic, + cross_helicity)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 10, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +cfl = 0.5 +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(); # + 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/examples/tree_2d_dgsem/convergence_subcell/elixir_mhd_alfven_wave_subcell.jl b/examples/tree_2d_dgsem/convergence_subcell/elixir_mhd_alfven_wave_subcell.jl new file mode 100644 index 00000000000..1c8087811b8 --- /dev/null +++ b/examples/tree_2d_dgsem/convergence_subcell/elixir_mhd_alfven_wave_subcell.jl @@ -0,0 +1,78 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible ideal GLM-MHD equations +gamma = 5 / 3 +equations = IdealGlmMhdEquations2D(gamma) + +initial_condition = initial_condition_convergence_test + +volume_flux = (flux_central, flux_nonconservative_powell2) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell2) + +basis = LobattoLegendreBasis(3) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons=[1], + positivity_correction_factor=0.1) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg=volume_flux, + volume_flux_fv=surface_flux) + +solver = DGSEM(polydeg = 3, + surface_flux = surface_flux, + volume_integral = volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.0)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (entropy, energy_total, + energy_kinetic, + energy_internal, + energy_magnetic, + cross_helicity)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 10, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +cfl = 0.5 +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation +stage_callbacks = (SubcellLimiterIDPCorrection(),) + +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/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 1d8f3543603..b16aa947fac 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -397,8 +397,56 @@ the non-conservative staggered "fluxes" for subcell limiting. See, e.g., return f end """ - flux_nonconservative_powell2(u_ll, orientation::Integer, - equations::IdealGlmMhdEquations2D, + flux_nonconservative_powell2(f, f_sym, 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. +This routine takes the symmetric part of the staggered flux `f_sym`, multiplies it with the +local part, and adds it to `f`. +Even though this routine modifies its arguments and returns nothing, we omit the exclamation mark +in the end to make it compatible with the other non-conservative flux functions +""" +@inline function flux_nonconservative_powell2(f, f_sym, 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 + + if noncons_term == 1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + f[2] += f_sym[2] * B1_ll + f[3] += f_sym[3] * B2_ll + f[4] += f_sym[4] * B3_ll + f[5] += f_sym[5] * v_dot_B_ll + f[6] += f_sym[6] * v1_ll + f[7] += f_sym[7] * v2_ll + f[8] += f_sym[8] * v3_ll + else #noncons_term ==2 + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + if orientation == 1 + v1_ll = rho_v1_ll / rho_ll + f[5] += f_sym[5] * v1_ll * psi_ll + f[9] += f_sym[9] * v1_ll + else #orientation == 2 + v2_ll = rho_v2_ll / rho_ll + f[5] += f_sym[5] * v2_ll * psi_ll + f[9] += f_sym[9] * v2_ll + end + end + return f +end +""" + flux_nonconservative_powell2(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations2D, nonconservative_type::NonConservativeSymmetric, noncons_term::Integer) @@ -455,6 +503,55 @@ the non-conservative staggered "fluxes" for subcell limiting. See, e.g., return f end +""" + flux_nonconservative_powell2(f1, f2, factor1, factor2, u_ll, u_rr, 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. +This routine computes the symmetric part of the staggered flux, multiplies it with factor1 and adds it +to f1, then multiplies it with factor2 and adds it to f2. +Even though this routine modifies its arguments and returns nothing, we omit the exclamation mark +in the end to make it compatible with the other non-conservative flux functions +""" +@inline function flux_nonconservative_powell2(f1, f2, factor1, factor2, 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 + + 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 + for v in 2:8 + f1[v] += factor1 * B1_avg + f2[v] += factor2 * B1_avg + end + 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 + for v in 2:8 + f1[v] += factor1 * B2_avg + f2[v] += factor2 * B2_avg + end + 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 + f1[5] += factor1 * psi_avg + f1[9] += factor1 * psi_avg + f2[5] += factor2 * psi_avg + f2[9] += factor2 * psi_avg + end + + return nothing +end """ flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations2D) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 080f24c6ce5..79919b4b0dd 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -38,18 +38,13 @@ function create_cache(mesh::TreeMesh{2}, 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_nonconservative_temp_threaded, fhat_temp_threaded, - fhat_nonconservative_temp_threaded, phi_threaded) + fhat_nonconservative_temp_threaded) end function calc_volume_integral!(du, u, @@ -82,9 +77,11 @@ 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, - 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 @@ -214,7 +211,7 @@ end volume_flux, dg::DGSEM, element, cache) @unpack weights, derivative_split = dg.basis @unpack flux_temp_threaded, flux_nonconservative_temp_threaded = cache - @unpack fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded = cache + @unpack fhat_temp_threaded, fhat_nonconservative_temp_threaded = cache volume_flux_cons, volume_flux_noncons = volume_flux @@ -223,7 +220,6 @@ 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, @@ -254,14 +250,12 @@ 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 = volume_flux_noncons(u_node, u_node_ii, 1, equations, - NonConservativeSymmetric(), noncons) - multiply_add_to_node_vars!(flux_noncons_temp, 0.5 * derivative_split[i, ii], - flux1_noncons, - equations, dg, noncons, i, j) - multiply_add_to_node_vars!(flux_noncons_temp, 0.5 * derivative_split[ii, i], - flux1_noncons, - equations, dg, noncons, ii, j) + volume_flux_noncons(view(flux_noncons_temp, :, noncons, i, j), + view(flux_noncons_temp, :, noncons, ii, j), + 0.5 * derivative_split[i, ii], + 0.5 * derivative_split[ii, i], + u_node, u_node_ii, 1, equations, + NonConservativeSymmetric(), noncons) end end end @@ -275,17 +269,6 @@ end fhat_temp[:, 1, :] .= zero(eltype(fhat1_L)) fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L)) - # Compute local contribution to non-conservative flux - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - for noncons in 1:nnoncons(equations) - set_node_vars!(phi, - volume_flux_noncons(u_local, 1, equations, - NonConservativeLocal(), noncons), - equations, dg, noncons, i, j) - end - end - for j in eachnode(dg), i in 1:(nnodes(dg) - 1) # Conservative part for v in eachvariable(equations) @@ -295,14 +278,24 @@ end fhat1_R[v, i + 1, j] = value end # Nonconservative part - for noncons in 1:nnoncons(equations), v in eachvariable(equations) - value = fhat_noncons_temp[v, noncons, i, j] + - weights[i] * flux_noncons_temp[v, noncons, i, j] - fhat_noncons_temp[v, noncons, i + 1, j] = value - - fhat1_L[v, i + 1, j] = fhat1_L[v, i + 1, j] + phi[v, noncons, i, j] * value - fhat1_R[v, i + 1, j] = fhat1_R[v, i + 1, j] + - phi[v, noncons, i + 1, j] * value + u_ll = get_node_vars(u, equations, dg, i, j, element) + u_rr = get_node_vars(u, equations, dg, i + 1, j, element) + 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_noncons_temp[v, noncons, + i, + j] + end + volume_flux_noncons(view(fhat1_L, :, i + 1, j), + view(fhat_noncons_temp, :, noncons, i + 1, j), + u_ll, 1, equations, NonConservativeLocal(), noncons) + volume_flux_noncons(view(fhat1_R, :, i + 1, j), + view(fhat_noncons_temp, :, noncons, i + 1, j), + u_rr, 1, equations, NonConservativeLocal(), noncons) end end @@ -321,14 +314,12 @@ 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 = volume_flux_noncons(u_node, u_node_jj, 2, equations, - NonConservativeSymmetric(), noncons) - multiply_add_to_node_vars!(flux_noncons_temp, 0.5 * derivative_split[j, jj], - flux2_noncons, - equations, dg, noncons, i, j) - multiply_add_to_node_vars!(flux_noncons_temp, 0.5 * derivative_split[jj, j], - flux2_noncons, - equations, dg, noncons, i, jj) + volume_flux_noncons(view(flux_noncons_temp, :, noncons, i, j), + view(flux_noncons_temp, :, noncons, i, jj), + 0.5 * derivative_split[j, jj], + 0.5 * derivative_split[jj, j], + u_node, u_node_jj, 2, equations, + NonConservativeSymmetric(), noncons) end end end @@ -342,17 +333,6 @@ end fhat_temp[:, :, 1] .= zero(eltype(fhat1_L)) fhat_noncons_temp[:, :, :, 1] .= zero(eltype(fhat1_L)) - # Compute local contribution to non-conservative flux - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - for noncons in 1:nnoncons(equations) - set_node_vars!(phi, - volume_flux_noncons(u_local, 2, equations, - NonConservativeLocal(), noncons), - equations, dg, noncons, i, j) - end - end - for j in 1:(nnodes(dg) - 1), i in eachnode(dg) # Conservative part for v in eachvariable(equations) @@ -362,14 +342,24 @@ end fhat2_R[v, i, j + 1] = value end # Nonconservative part - for noncons in 1:nnoncons(equations), v in eachvariable(equations) - value = fhat_noncons_temp[v, noncons, i, j] + - weights[j] * flux_noncons_temp[v, noncons, i, j] - fhat_noncons_temp[v, noncons, i, j + 1] = value - - fhat2_L[v, i, j + 1] = fhat2_L[v, i, j + 1] + phi[v, noncons, i, j] * value - fhat2_R[v, i, j + 1] = fhat2_R[v, i, j + 1] + - phi[v, noncons, i, j + 1] * value + u_ll = get_node_vars(u, equations, dg, i, j, element) + u_rr = get_node_vars(u, equations, dg, i, j + 1, element) + 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_noncons_temp[v, noncons, + i, + j] + end + volume_flux_noncons(view(fhat2_L, :, i, j + 1), + view(fhat_noncons_temp, :, noncons, i, j + 1), + u_ll, 2, equations, NonConservativeLocal(), noncons) + volume_flux_noncons(view(fhat2_R, :, i, j + 1), + view(fhat_noncons_temp, :, noncons, i, j + 1), + u_rr, 2, equations, NonConservativeLocal(), noncons) end end return nothing