first draft IMEX JinXin

examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl

@@ -1,4 +1,3 @@
-
 using OrdinaryDiffEq
 using Trixi
 

examples/tree_1d_dgsem/elixir_jin_xin_burgers.jl

@@ -0,0 +1,56 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the (inviscid) Burgers' equation
+epsilon_relaxation = 1.0e-4
+a1 = 9.0
+
+equations_base = InviscidBurgersEquation1D()
+velocities = (SVector(a1),)
+equations = JinXinEquations(equations_base, epsilon_relaxation, velocities)
+function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D)
+    k = 1
+    u = 2 + sinpi(k * (x[1] - 0.7))
+    return prim2cons(SVector(u), equations)
+end
+
+basis = LobattoLegendreBasis(3; polydeg_projection = 12, polydeg_cutoff = 3)
+solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakForm())
+
+coordinates_min = -1.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000)
+initial_condition = Trixi.InitialConditionJinXin(initial_condition_linear_stability)
+
+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)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = Trixi.solve(ode, Trixi.SimpleIMEX(),
+                  dt = 1.0e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks, maxiters = 1e7);
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_jin_xin_euler_density_wave.jl

@@ -0,0 +1,68 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+epsilon_relaxation = 1.0e-6
+
+equations_base = CompressibleEulerEquations1D(1.4)
+velocities = (SVector(a1, a2, a3),)
+equations = JinXinEquations(equations_base, epsilon_relaxation, velocities)
+
+function initial_condition_density_wave(x, t, equations::CompressibleEulerEquations1D)
+    v1 = 0.1
+    rho = 1 + 0.98 * sinpi(2 * (x[1] - t * v1))
+    p = 20
+    return prim2cons(SVector(rho, v1, p), equations)
+end
+
+initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave)
+polydeg = 3
+#basis = LobattoLegendreBasis(polydeg; polydeg_projection = 0)
+basis = LobattoLegendreBasis(polydeg; polydeg_projection = 6)
+
+volume_integral = VolumeIntegralWeakForm()
+#solver = DGSEM(basis, Trixi.flux_upwind,VolumeIntegralWeakForm())
+solver = DGSEM(basis, Trixi.flux_upwind)
+
+coordinates_min = -1.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                n_cells_max = 30_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 = 2000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = Trixi.solve(ode, Trixi.SimpleIMEX(),
+                  dt = 1e-3, # 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

examples/tree_2d_dgsem/elixir_euler_ec.jl

@@ -8,12 +8,12 @@ equations = CompressibleEulerEquations2D(1.4)
 
 seed!(1)
 function initial_condition_random_field(x, t, equations::CompressibleEulerEquations2D)
-amplitude = 1.5
-rho = 2 + amplitude * rand() 
-v1 = -3.1 + amplitude * rand()
-v2 = 1.3 + amplitude * rand() 
-p = 7.54 + amplitude * rand()
-return prim2cons(SVector(rho, v1, v2, p), equations)
+    amplitude = 1.5
+    rho = 2 + amplitude * rand()
+    v1 = -3.1 + amplitude * rand()
+    v2 = 1.3 + amplitude * rand()
+    p = 7.54 + amplitude * rand()
+    return prim2cons(SVector(rho, v1, v2, p), equations)
 end
 # initial_condition = initial_condition_weak_blast_wave
 initial_condition = initial_condition_random_field
@@ -69,11 +69,13 @@ callbacks = CallbackSet(summary_callback,
 
 # Create modal filter and filter initial condition
 modal_filter = ModalFilter(solver; polydeg_cutoff = 3,
-                                   cons2filter = cons2prim, filter2cons = prim2cons)
+                           cons2filter = cons2prim, filter2cons = prim2cons)
 modal_filter(ode.u0, semi)
 
 # sol = solve(ode, CarpenterKennedy2N54(; williamson_condition = false),
-sol = solve(ode, CarpenterKennedy2N54(; stage_limiter! = modal_filter, williamson_condition = false),
+sol = solve(ode,
+            CarpenterKennedy2N54(; stage_limiter! = modal_filter,
+                                 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

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl

@@ -33,7 +33,7 @@ end
 initial_condition = initial_condition_kelvin_helmholtz_instability
 
 surface_flux = flux_hllc
-volume_flux  = flux_ranocha
+volume_flux = flux_ranocha
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
@@ -67,7 +67,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval=1000,
+save_solution = SaveSolutionCallback(interval = 1000,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl

@@ -69,7 +69,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval=1000,
+save_solution = SaveSolutionCallback(interval = 1000,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)

examples/tree_2d_dgsem/elixir_jin_xin_euler.jl

@@ -1,4 +1,3 @@
-
 using OrdinaryDiffEq
 using Trixi
 
@@ -13,29 +12,31 @@ equations_base = CompressibleEulerEquations2D(1.4)
 velocities = (SVector(a1, a2, a3, a4), SVector(b1, b2, b3, b4))
 equations = JinXinEquations(equations_base, epsilon_relaxation, velocities)
 
-function initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
-  # change discontinuity to tanh
-  # typical resolution 128^2, 256^2
-  # domain size is [-1,+1]^2
-  slope = 15
-  amplitude = 0.02
-  B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
-  rho = 0.5 + 0.75 * B
-  v1 = 0.5 * (B - 1)
-  v2 = 0.1 * sin(2 * pi * x[1])
-  p = 1.0
-  return prim2cons(SVector(rho, v1, v2, p), equations)
+function initial_condition_kelvin_helmholtz_instability(x, t,
+                                                        equations::CompressibleEulerEquations2D)
+    # change discontinuity to tanh
+    # typical resolution 128^2, 256^2
+    # domain size is [-1,+1]^2
+    slope = 15
+    amplitude = 0.02
+    B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
+    rho = 0.5 + 0.75 * B
+    v1 = 0.5 * (B - 1)
+    v2 = 0.1 * sin(2 * pi * x[1])
+    p = 1.0
+    return prim2cons(SVector(rho, v1, v2, p), equations)
 end
 
 #initial_condition = initial_condition_constant
 initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability)
 #initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave)
 polydeg = 3
 polydeg_cutoff = 3
-basis = GaussLegendreBasis(polydeg; polydeg_projection = polydeg, polydeg_cutoff = polydeg_cutoff)
-# solver = DGSEM(basis, Trixi.flux_upwind)
+#basis = GaussLegendreBasis(polydeg; polydeg_projection = polydeg, polydeg_cutoff = polydeg_cutoff)
+basis = LobattoLegendreBasis(polydeg; polydeg_projection = 6)
+solver = DGSEM(basis, Trixi.flux_upwind)
 #solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakFormProjection())
-solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakForm())
+#solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakForm())
 #solver = DGSEM(polydeg = 3, surface_flux = Trixi.flux_upwind)
 
 #surface_flux = Trixi.flux_upwind
@@ -50,49 +51,50 @@ solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakForm())
 #solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinates_min = (-1.0, -1.0)
-coordinates_max = ( 1.0,  1.0)
+coordinates_max = (1.0, 1.0)
 
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level=6,
-                n_cells_max=1_000_000)
+                initial_refinement_level = 6,
+                n_cells_max = 1_000_000)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)#,source_terms=source_terms_JinXin_Relaxation)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 3.7)
+tspan = (0.0, 3.6)
 #tspan = (0.0, 1.0)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
 analysis_interval = 1000
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
-alive_callback = AliveCallback(analysis_interval=analysis_interval)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval=1000,
-                                     save_initial_solution=true,
-                                     save_final_solution=true,
-                                     solution_variables=cons2prim)
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl=0.5)
+stepsize_callback = StepsizeCallback(cfl = 0.5)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,stepsize_callback)
+                        save_solution, stepsize_callback)
 
 ###############################################################################
 # run the simulation
-stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
-                                                     variables=(Trixi.density, pressure))
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6),
+                                                     variables = (Trixi.density, pressure))
 
 #sol = solve(ode, CarpenterKennedy2N54(stage_limiter!,williamson_condition=false),
-sol = solve(ode, SSPRK43(stage_limiter!),
-#sol = solve(ode, SSPRK33(stage_limiter!),
-#sol = solve(ode, RDPK3SpFSAL49(),
-#sol = solve(ode, AutoTsit5(Rosenbrock23()),
-            dt=1.0e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks,maxiters=1e7);
+#sol = solve(ode, SSPRK43(stage_limiter!),
+sol = Trixi.solve(ode, Trixi.SimpleIMEX(),
+                  #sol = solve(ode, SSPRK33(stage_limiter!),
+                  #sol = solve(ode, RDPK3SpFSAL49(),
+                  #sol = solve(ode, AutoTsit5(Rosenbrock23()),
+                  dt = 1.0e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks, maxiters = 1e7);
 summary_callback() # print the timer summary