Merge pull request #256 from trixi-framework/shock-non-cons-in-taal

Port shock capturing with non-conservative terms into Taal

Author: andrewwinters5000

AUTHORS.md

@@ -24,5 +24,6 @@ are listed in alphabetical order:
 * Gregor Gassner
 * Julia Odenthal
 * Hendrik Ranocha
+* Andrés M. Rueda-Ramírez
 * Michael Schlottke-Lakemper
 * Andrew Winters

examples/2d/elixir_mhd_blast_wave_shockcapturing_amr.jl

@@ -52,13 +52,13 @@ amr_indicator = IndicatorHennemannGassner(semi,
                                           variable=density_pressure)
 amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                       base_level=4,
-                                      max_level =6, max_threshold=0.01)
+                                      max_level =7, max_threshold=0.01)
 amr_callback = AMRCallback(semi, amr_controller,
                            interval=5,
                            adapt_initial_condition=true,
                            adapt_initial_condition_only_refine=true)
 
-stepsize_callback = StepsizeCallback(cfl=0.5) # can probably be increased when shock-capturing is fixed for MHD
+stepsize_callback = StepsizeCallback(cfl=1.0)
 
 save_solution = SaveSolutionCallback(interval=100,
                                      save_initial_solution=true,

examples/2d/elixir_mhd_orszag_tang_shockcapturing_amr.jl

@@ -14,7 +14,7 @@ equations = IdealGlmMhdEquations2D(gamma)
 
 initial_condition = initial_condition_orszag_tang
 
-surface_flux = flux_hll
+surface_flux = flux_lax_friedrichs
 volume_flux  = flux_central
 basis = LobattoLegendreBasis(3)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
@@ -58,7 +58,7 @@ amr_callback = AMRCallback(semi, amr_controller,
                            adapt_initial_condition=true,
                            adapt_initial_condition_only_refine=true)
 
-stepsize_callback = StepsizeCallback(cfl=0.25) # can probably be increased when shock-capturing is fixed for MHD
+stepsize_callback = StepsizeCallback(cfl=1.0)
 
 save_solution = SaveSolutionCallback(interval=100,
                                      save_initial_solution=true,

examples/2d/elixir_mhd_rotor_shockcapturing_amr.jl

@@ -14,7 +14,7 @@ equations = IdealGlmMhdEquations2D(1.4)
 initial_condition = initial_condition_rotor
 
 surface_flux = flux_lax_friedrichs
-volume_flux  = flux_central
+volume_flux  = flux_derigs_etal
 basis = LobattoLegendreBasis(4)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
                                          alpha_max=0.5,
@@ -57,7 +57,7 @@ amr_callback = AMRCallback(semi, amr_controller,
                            adapt_initial_condition=true,
                            adapt_initial_condition_only_refine=true)
 
-stepsize_callback = StepsizeCallback(cfl=0.5) # can probably be increased when shock-capturing is fixed for MHD
+stepsize_callback = StepsizeCallback(cfl=0.95)
 
 save_solution = SaveSolutionCallback(interval=100,
                                      save_initial_solution=true,

examples/2d/parameters_mhd_rotor.toml

@@ -5,7 +5,7 @@ equations = "IdealGlmMhdEquations"
 gamma = 1.4
 initial_condition = "initial_condition_rotor"
 surface_flux = "flux_lax_friedrichs"
-volume_flux = "flux_central"
+volume_flux = "flux_derigs_etal"
 
 amr_indicator = "blast_wave" # for now just for testing
 

src/equations/2d/ideal_glm_mhd.jl

@@ -453,14 +453,25 @@ function flux_hll(u_ll, u_rr, orientation, equation::IdealGlmMhdEquations2D)
   return SVector(f1, f2, f3, f4, f5, f6, f7, f8, f9)
 end
 
-
-# strong form of nonconservative flux on a side, e.g., the Powell term
-#     phi^L 1/2 (B^L+B^R) normal - phi^L B^L normal = phi^L 1/2 (B^R-B^L) normal
-# OBS! 1) "weak" formulation of split DG already includes the contribution -1/2(phi^L B^L normal)
-#         so this routine only adds 1/2(phi^L B^R normal)
-#         analogously for the Galilean nonconservative term
-#      2) this is non-unique along an interface! normal direction is super important
-function noncons_interface_flux(u_left, u_right, orientation, equation::IdealGlmMhdEquations2D)
+"""
+    noncons_interface_flux(u_left, u_right, orientation, mode, equation::IdealGlmMhdEquations2D)
+
+Strong form of non-conservative flux on a surface (Powell and GLM terms)
+```math
+phi^L 1/2 (B^L + B^R)_{normal} - phi^L B^L+{normal} = phi^L 1/2 (B^R - B^L)_{normal}
+```
+!!! note
+    The non-conservative interface flux depends on the discretization. Following "modes" are available:
+    * `:weak`: 'weak' formulation of split DG already includes the contribution
+      ``-1/2 (phi^L B^L_{normal})`` so this mode only adds ``1/2 (phi^L B^R_{normal})``,
+      analogously for the Galilean nonconservative term
+    * `:whole`: This mode adds the whole non-conservative term: phi^L 1/2 (B^R-B^L)
+    * `:inner`: This mode adds the split-form DG volume integral contribution: This is equivalent to
+      ``(2)-(1) - 1/2 (phi^L B^L)``
+!!! warning
+    This is non-unique along an interface! The normal direction is super important.
+"""
+@inline function noncons_interface_flux(u_left, u_right, orientation, mode, equation::IdealGlmMhdEquations2D)
   rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _, B1_ll, B2_ll, B3_ll, psi_ll = u_left
   _, _, _, _, _, B1_rr, B2_rr, _, psi_rr = u_right
 
@@ -471,27 +482,48 @@ function noncons_interface_flux(u_left, u_right, orientation, equation::IdealGlm
   v_dot_B_ll = v1_ll*B1_ll + v2_ll*B2_ll + v3_ll*B3_ll
   # extract magnetic field variable from the right and set the normal velocity
   # Note, both depend upon the orientation and need psi_rr
-  if orientation == 1 # x-direction
-    v_normal = v1_ll
-    B_normal = B1_rr
-  else # y-direction
-    v_normal = v2_ll
-    B_normal = B2_rr
+  if mode==:weak
+    if orientation == 1 # x-direction
+      v_normal = v1_ll
+      B_normal = B1_rr
+    else # y-direction
+      v_normal = v2_ll
+      B_normal = B2_rr
+    end
+    psi_norm = psi_rr
+  elseif mode==:whole
+    if orientation == 1 # x-direction
+      v_normal = v1_ll
+      B_normal = B1_rr - B1_ll
+    else # y-direction
+      v_normal = v2_ll
+      B_normal = B2_rr - B2_ll
+    end
+    psi_norm = psi_rr - psi_ll
+  else #mode==:inner
+    if orientation == 1 # x-direction
+      v_normal = v1_ll
+      B_normal =-B1_ll
+    else # y-direction
+      v_normal = v2_ll
+      B_normal =-B2_ll
+    end
+    psi_norm =-psi_ll
   end
+
   # compute the nonconservative flux: Powell (with B_normal) and Galilean (with v_normal)
   noncons2 = 0.5 * B_normal * B1_ll
   noncons3 = 0.5 * B_normal * B2_ll
   noncons4 = 0.5 * B_normal * B3_ll
-  noncons5 = 0.5 * B_normal * v_dot_B_ll + 0.5 * v_normal * psi_ll * psi_rr
+  noncons5 = 0.5 * B_normal * v_dot_B_ll + 0.5 * v_normal * psi_ll * psi_norm
   noncons6 = 0.5 * B_normal * v1_ll
   noncons7 = 0.5 * B_normal * v2_ll
   noncons8 = 0.5 * B_normal * v3_ll
-  noncons9 = 0.5 * v_normal * psi_rr
+  noncons9 = 0.5 * v_normal * psi_norm
 
   return SVector(0, noncons2, noncons3, noncons4, noncons5, noncons6, noncons7, noncons8, noncons9)
 end
 
-
 # 1) Determine maximum stable time step based on polynomial degree and CFL number
 # 2) Update the GLM cleaning wave speed c_h to be the largest value of the fast
 #    magnetoacoustic over the entire domain (note this routine is called in a loop