added another (more difficult) shock tube test from Manuel Torrilhon

examples/1d/elixir_mhd_torrilhon_shock_tube.jl

@@ -0,0 +1,77 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal MHD equations
+gamma = 5/3
+equations = IdealGlmMhdEquations1D(gamma)
+
+initial_condition = initial_condition_torrilhon_shock_tube
+
+boundary_conditions = boundary_condition_torrilhon_shock_tube
+
+surface_flux = flux_lax_friedrichs
+volume_flux  = flux_central
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max=0.5,
+                                         alpha_min=0.001,
+                                         alpha_smooth=true,
+                                         variable=pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg=volume_flux,
+                                                 volume_flux_fv=surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = -1.0
+coordinates_max =  1.5
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=7,
+                n_cells_max=10_000,
+                periodicity=false)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_errors=(:l2_error_primitive,
+                                                            :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=:primitive)
+
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart, save_solution,
+                        stepsize_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
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -81,11 +81,11 @@ export initial_condition_constant,
        initial_condition_blob,
        initial_condition_orszag_tang,
        initial_condition_rotor,
-       initial_condition_briowu_shock_tube
+       initial_condition_briowu_shock_tube, initial_condition_torrilhon_shock_tube
 
 export boundary_condition_periodic,
        boundary_condition_gauss,
-       boundary_condition_briowu_shock_tube
+       boundary_condition_briowu_shock_tube, boundary_condition_torrilhon_shock_tube
 
 export initial_condition_convergence_test, source_terms_convergence_test, boundary_condition_convergence_test
 export initial_condition_harmonic_nonperiodic, source_terms_harmonic, boundary_condition_harmonic_nonperiodic

src/equations/ideal_glm_mhd_1d.jl

@@ -99,6 +99,52 @@ function boundary_condition_briowu_shock_tube(u_inner, orientation, direction, x
   return flux
 end
 
+
+"""
+    initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
+
+Torrilhon's shock tube test case for one dimensional ideal MHD equations.
+- Torrilhon (2003)
+  Uniqueness conditions for Riemann problems of ideal magnetohydrodynamics
+  [DOI: 10.1017/S0022377803002186](https://doi.org/10.1017/S0022377803002186)
+
+"""
+function initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
+  # domain must be set to [-1, 1.5], γ = 5/3, final time = 0.4
+  rho = x[1] <= 0 ? 3.0 : 1.0
+  v1 = 0.0
+  v2 = 0.0
+  v3 = 0.0
+  p = x[1] <= 0 ? 3.0 : 1.0
+  B1 = 1.5
+  B2 = x[1] <= 0 ? 1.0 : cos(1.5)
+  B3 = x[1] <= 0 ? 0.0 : sin(1.5)
+  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
+end
+
+"""
+    boundary_condition_torrilhon_shock_tube(u_inner, orientation, direction, x, t,
+                                            surface_flux_function,
+                                            equations::IdealGlmMhdEquations1D)
+
+Boundary conditions used for the Torrilhon shock tube in combination with
+[`initial_condition_torrilhon_shock_tube`](@ref).
+"""
+function boundary_condition_torrilhon_shock_tube(u_inner, orientation, direction, x, t,
+                                             surface_flux_function,
+                                             equations::IdealGlmMhdEquations1D)
+  u_boundary = initial_condition_torrilhon_shock_tube(x, t, equations)
+  # Calculate boundary flux
+  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
+    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+  end
+
+  return flux
+end
+
+
 # Calculate 1D flux in for a single point
 @inline function calcflux(u, orientation, equations::IdealGlmMhdEquations1D)
   rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u

test/test_examples_1d_mhd.jl

@@ -35,6 +35,12 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "1d")
       l2   = [0.17764301067932906, 0.19693621875378622, 0.3635136528288642, 0.0, 0.3757321708837591, 8.593007507325741e-16, 0.36473438378159656, 0.0],
       linf = [0.5601530250396535, 0.43867368105486537, 1.0960903616351099, 0.0, 1.0551794137886303, 4.107825191113079e-15, 1.5374410890043144, 0.0])
   end
+
+  @testset "elixir_mhd_torrilhon_shock_tube.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"),
+      l2   = [0.45702524713047815, 0.4792453322653316, 0.3406561762729464, 0.4478037841848423, 0.9204646319222438, 1.3216517820475193e-16, 0.2889772516628848, 0.2552071770333965],
+      linf = [1.221507578527215, 0.8927062606782392, 0.8536666880170606, 0.9739151339925485, 1.6602343391033632, 2.220446049250313e-16, 0.686743631371773, 0.6555428379163893])
+  end
 end
 
 end # module