From 7107447a14b54ec7d3cf28cd48490a5260304d96 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 21 Mar 2026 10:18:39 +0100 Subject: [PATCH 1/6] start --- .../tree_3d_dgsem/elixir_mhd_convergence.jl | 120 ++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 examples/tree_3d_dgsem/elixir_mhd_convergence.jl diff --git a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl new file mode 100644 index 00000000000..c3944cdf906 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl @@ -0,0 +1,120 @@ +using Trixi +using OrdinaryDiffEqLowStorageRK + +############################################################################### +# semidiscretization of the compressible ideal GLM-MHD equations + +# Fluid parameters +gamma() = 2.0 +equations = IdealGlmMhdEquations3D(gamma()) + +""" + initial_condition_convergence(x, t, equations::IdealGlmMhdEquations3D) + +Manufactured solution for the 3D(only!) compressible ideal GLM-MHD equations. +Proposed in +- Marvin Bohm, Andrew R Winters, Gregor J Gassner, Dominik Derigs, Florian Hindenlang, Joachim Saur (2020): + An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. + Part I: Theory and numerical verification + [10.1016/j.jcp.2018.06.027](https://doi.org/10.1016/j.jcp.2018.06.027) +""" +@inline function initial_condition_convergence(x, t, equations::IdealGlmMhdEquations3D) + h = 0.5 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + 2 + + u_1 = h + u_2 = h + u_3 = h + u_4 = 0.0 + u_5 = 2 * h^2 + h + + u_6 = h + u_7 = -h + u_8 = 0.0 + u_9 = 0.0 + + return SVector(u_1, u_2, u_3, u_4, u_5, u_6, u_7, u_8, u_9) +end + +""" + source_terms_convergence(x, t, equations::IdealGlmMhdEquations3D) + +Manufactured solution for the 3D(only!) compressible ideal GLM-MHD equations. +Proposed in +- Marvin Bohm, Andrew R Winters, Gregor J Gassner, Dominik Derigs, Florian Hindenlang, Joachim Saur (2020): + An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. + Part I: Theory and numerical verification + [10.1016/j.jcp.2018.06.027](https://doi.org/10.1016/j.jcp.2018.06.027) + +For the version without parabolic terms see the implementation in FLUXO: +https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/mhd/equation.f90#L1539-L1554 +""" +function source_terms_convergence(u, x, t, equations::IdealGlmMhdEquations3D) + h = 0.5 * sinpi(2 * (x[1] + x[2] + x[3]- t)) + 2 + h_x = pi * cospi(2 * (x[1] + x[2] +x[3] - t)) + h_xx = -2 * pi^2 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + + s_1 = h_x + s_2 = h_x + 4 * h * h_x + s_3 = h_x + 4 * h * h_x + s_4 = 4 * h * h_x + s_5 = h_x + 12 * h * h_x + s_6 = h_x + s_7 = -h_x + s_8 = 0.0 + s_9 = 0.0 + + return SVector(s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8, s_9) +end + +surface_flux = (flux_hll, flux_nonconservative_powell) +volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (1.0, 1.0, 1.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + periodicity = true, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, + initial_condition_convergence, solver; + source_terms = source_terms_convergence, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 50 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +cfl = 1.8 +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + ode_default_options()..., callback = callbacks); \ No newline at end of file From f6b39e79c9e0cf2d0fa8d43b3797772ee93db4ff Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sat, 21 Mar 2026 10:22:20 +0100 Subject: [PATCH 2/6] fmt --- .../tree_3d_dgsem/elixir_mhd_convergence.jl | 7 ++--- test/test_tree_3d_mhd.jl | 29 +++++++++++++++++++ 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl index c3944cdf906..9b437bc0f62 100644 --- a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl @@ -4,7 +4,6 @@ using OrdinaryDiffEqLowStorageRK ############################################################################### # semidiscretization of the compressible ideal GLM-MHD equations -# Fluid parameters gamma() = 2.0 equations = IdealGlmMhdEquations3D(gamma()) @@ -49,8 +48,8 @@ For the version without parabolic terms see the implementation in FLUXO: https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/mhd/equation.f90#L1539-L1554 """ function source_terms_convergence(u, x, t, equations::IdealGlmMhdEquations3D) - h = 0.5 * sinpi(2 * (x[1] + x[2] + x[3]- t)) + 2 - h_x = pi * cospi(2 * (x[1] + x[2] +x[3] - t)) + h = 0.5 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + 2 + h_x = pi * cospi(2 * (x[1] + x[2] + x[3] - t)) h_xx = -2 * pi^2 * sinpi(2 * (x[1] + x[2] + x[3] - t)) s_1 = h_x @@ -117,4 +116,4 @@ callbacks = CallbackSet(summary_callback, sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false); dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - ode_default_options()..., callback = callbacks); \ No newline at end of file + ode_default_options()..., callback = callbacks); diff --git a/test/test_tree_3d_mhd.jl b/test/test_tree_3d_mhd.jl index 2f1dc252409..b2d42eb8b6f 100644 --- a/test/test_tree_3d_mhd.jl +++ b/test/test_tree_3d_mhd.jl @@ -221,6 +221,35 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_mhd_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_convergence.jl"), + l2=[ + 0.009193403522877426, + 0.013723993903671483, + 0.013723993903671291, + 0.00748865999127907, + 0.02205355095416697, + 0.009989774083000539, + 0.009989774083000584, + 0.004353727273126007, + 0.0013769046942994955 + ], + linf=[ + 0.027349064091453767, + 0.04489621477501449, + 0.04489621477501515, + 0.02699974461840463, + 0.15656498361977533, + 0.028682970290997645, + 0.028682970290998533, + 0.0124215159258882, + 0.007286608218242374 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"), l2=[ From 85b97278271dcd5a13e8bd677a58a2771fd1b2a6 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 21 Mar 2026 10:34:35 +0100 Subject: [PATCH 3/6] Update examples/tree_3d_dgsem/elixir_mhd_convergence.jl --- examples/tree_3d_dgsem/elixir_mhd_convergence.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl index 9b437bc0f62..27636971901 100644 --- a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl @@ -50,7 +50,6 @@ https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a848 function source_terms_convergence(u, x, t, equations::IdealGlmMhdEquations3D) h = 0.5 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + 2 h_x = pi * cospi(2 * (x[1] + x[2] + x[3] - t)) - h_xx = -2 * pi^2 * sinpi(2 * (x[1] + x[2] + x[3] - t)) s_1 = h_x s_2 = h_x + 4 * h * h_x From d241e5759eb9131dddfdf4fa6b7dc4a05ea8380d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Mar 2026 08:27:30 +0100 Subject: [PATCH 4/6] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/tree_3d_dgsem/elixir_mhd_convergence.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl index 27636971901..4bac968a2fe 100644 --- a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl @@ -18,18 +18,18 @@ Proposed in [10.1016/j.jcp.2018.06.027](https://doi.org/10.1016/j.jcp.2018.06.027) """ @inline function initial_condition_convergence(x, t, equations::IdealGlmMhdEquations3D) - h = 0.5 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + 2 + h = 0.5f0 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + 2 u_1 = h u_2 = h u_3 = h - u_4 = 0.0 + u_4 = 0 u_5 = 2 * h^2 + h u_6 = h u_7 = -h - u_8 = 0.0 - u_9 = 0.0 + u_8 = 0 + u_9 = 0 return SVector(u_1, u_2, u_3, u_4, u_5, u_6, u_7, u_8, u_9) end @@ -48,7 +48,7 @@ For the version without parabolic terms see the implementation in FLUXO: https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/mhd/equation.f90#L1539-L1554 """ function source_terms_convergence(u, x, t, equations::IdealGlmMhdEquations3D) - h = 0.5 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + 2 + h = 0.5f0 * sinpi(2 * (x[1] + x[2] + x[3] - t)) + 2 h_x = pi * cospi(2 * (x[1] + x[2] + x[3] - t)) s_1 = h_x @@ -58,8 +58,8 @@ function source_terms_convergence(u, x, t, equations::IdealGlmMhdEquations3D) s_5 = h_x + 12 * h * h_x s_6 = h_x s_7 = -h_x - s_8 = 0.0 - s_9 = 0.0 + s_8 = 0 + s_9 = 0 return SVector(s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8, s_9) end From 3cc4dc548fcf65c78154f6085583804dd551c820 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Mar 2026 08:28:05 +0100 Subject: [PATCH 5/6] Update examples/tree_3d_dgsem/elixir_mhd_convergence.jl --- examples/tree_3d_dgsem/elixir_mhd_convergence.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl index 4bac968a2fe..deb56dcea55 100644 --- a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl @@ -4,7 +4,7 @@ using OrdinaryDiffEqLowStorageRK ############################################################################### # semidiscretization of the compressible ideal GLM-MHD equations -gamma() = 2.0 +gamma() = 2.0 # required to make solution below valid! equations = IdealGlmMhdEquations3D(gamma()) """ From 757a77f33748d1a08a0de8ca2700f20f619aeaa6 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 23 Mar 2026 08:28:19 +0100 Subject: [PATCH 6/6] Update examples/tree_3d_dgsem/elixir_mhd_convergence.jl --- examples/tree_3d_dgsem/elixir_mhd_convergence.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl index deb56dcea55..c92141fd2f1 100644 --- a/examples/tree_3d_dgsem/elixir_mhd_convergence.jl +++ b/examples/tree_3d_dgsem/elixir_mhd_convergence.jl @@ -4,7 +4,7 @@ using OrdinaryDiffEqLowStorageRK ############################################################################### # semidiscretization of the compressible ideal GLM-MHD equations -gamma() = 2.0 # required to make solution below valid! +gamma() = 2.0 # required to make solution below working! equations = IdealGlmMhdEquations3D(gamma()) """