WIP: Jin Xin Relaxation

examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability.jl

@@ -1,6 +1,6 @@
 using Trixi, OrdinaryDiffEq
 
-dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
+dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
              surface_integral = SurfaceIntegralWeakForm(FluxLaxFriedrichs()),
              volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
 
@@ -31,7 +31,7 @@ function initial_condition_kelvin_helmholtz_instability(x, t,
 end
 initial_condition = initial_condition_kelvin_helmholtz_instability
 
-cells_per_dimension = (32, 32)
+cells_per_dimension = (64, 64)
 mesh = DGMultiMesh(dg, cells_per_dimension; periodicity = true)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
@@ -40,8 +40,8 @@ tspan = (0.0, 1.0)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
-alive_callback = AliveCallback(alive_interval = 10)
-analysis_interval = 100
+alive_callback = AliveCallback(alive_interval = 100)
+analysis_interval = 1000
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
@@ -50,7 +50,8 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+#sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+sol = solve(ode, SSPRK43(),
             dt = estimate_dt(mesh, dg), save_everystep = false, callback = callbacks);
 
 summary_callback() # print the timer summary

examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl

@@ -11,14 +11,18 @@ initial_condition = initial_condition_convergence_test
 
 surface_flux = flux_hllc
 
-basis = LobattoLegendreBasis(3)
-volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc)
+# basis = LobattoLegendreBasis(3)
+# volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc)
+# solver = DGSEM(basis, surface_flux, volume_integral)
+polydeg = 4
+basis = GaussLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 4)
+volume_integral = VolumeIntegralWeakForm()
 solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinates_min = (0.0, 0.0)
 coordinates_max = (2.0, 2.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 4,
+                initial_refinement_level = 2,
                 n_cells_max = 10_000)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
@@ -32,7 +36,7 @@ ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 100
+analysis_interval = 100 * 10
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
@@ -46,7 +50,7 @@ stepsize_callback = StepsizeCallback(cfl = 0.5)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        # save_solution,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_euler_density_wave.jl

@@ -1,23 +1,52 @@
 
 using OrdinaryDiffEq
 using Trixi
+using LinearAlgebra
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 gamma = 1.4
 equations = CompressibleEulerEquations2D(gamma)
 
-initial_condition = initial_condition_density_wave
-
-solver = DGSEM(polydeg = 5, surface_flux = flux_central)
+function initial_condition_density_wave2(x, t, equations::CompressibleEulerEquations2D)
+    v1 = 0.1
+    v2 = 0.2
+    rho = 1 + 0.98 * sinpi(2 * (x[1] + x[2] - t * (v1 + v2)))
+    rho_v1 = rho * v1
+    rho_v2 = rho * v2
+    p = 20
+    rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2)
+    return SVector(rho, rho_v1, rho_v2, rho_e)
+end
+initial_condition = initial_condition_density_wave2
+
+surface_flux = flux_central
+polydeg = 17
+basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 5)
+volume_integral = VolumeIntegralWeakForm()
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# solver = DGSEM(polydeg = 5, surface_flux = flux_central)
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level = 2,
                 n_cells_max = 30_000)
 
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+# @info "Create semi..."
+# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+# @info "Compute Jacobian..."
+# J = jacobian_ad_forward(semi)
+
+# @info "Compute eigenvalues..."
+# λ = eigvals(J)
+
+# @info "max real part" maximum(real.(λ))
+# # @info "Plot spectrum..."
+# # scatter(real.(λ), imag.(λ), label="central flux")
+# wololo
 
 ###############################################################################
 # ODE solvers, callbacks etc.
@@ -27,7 +56,7 @@ ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 100
+analysis_interval = 100 * 20
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
@@ -41,7 +70,7 @@ stepsize_callback = StepsizeCallback(cfl = 1.6)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        # save_solution,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_euler_ec.jl

@@ -1,16 +1,33 @@
 
 using OrdinaryDiffEq
 using Trixi
-
+using Random: seed!
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 equations = CompressibleEulerEquations2D(1.4)
 
-initial_condition = initial_condition_weak_blast_wave
-
-volume_flux = flux_ranocha
-solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
-               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+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)
+end
+# initial_condition = initial_condition_weak_blast_wave
+initial_condition = initial_condition_random_field
+
+#volume_flux = flux_ranocha
+#solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
+#               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+#
+surface_flux = flux_ranocha
+polydeg = 3
+basis = GaussLegendreBasis(polydeg; polydeg_projection = 1 * polydeg, polydeg_cutoff = 3)
+volume_integral = VolumeIntegralWeakFormProjection()
+#volume_integral = VolumeIntegralWeakForm()
+solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinates_min = (-2.0, -2.0)
 coordinates_max = (2.0, 2.0)
@@ -40,17 +57,23 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 1.0)
+stepsize_callback = StepsizeCallback(cfl = 0.5)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
+                        #save_solution,
                         stepsize_callback)
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+# Create modal filter and filter initial condition
+modal_filter = ModalFilter(solver; polydeg_cutoff = 3,
+                                   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),
             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

@@ -32,15 +32,16 @@ function initial_condition_kelvin_helmholtz_instability(x, t,
 end
 initial_condition = initial_condition_kelvin_helmholtz_instability
 
-surface_flux = flux_lax_friedrichs
-volume_flux = flux_ranocha
+surface_flux = flux_hllc
+volume_flux  = flux_ranocha
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
-                                         alpha_max = 0.002,
+                                         alpha_max = 0.2,
                                          alpha_min = 0.0001,
                                          alpha_smooth = true,
-                                         variable = density_pressure)
+                                         variable = Trixi.density)
+
 volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                  volume_flux_dg = volume_flux,
                                                  volume_flux_fv = surface_flux)
@@ -49,29 +50,29 @@ solver = DGSEM(basis, surface_flux, volume_integral)
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 5,
+                initial_refinement_level = 6,
                 n_cells_max = 100_000)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 3.0)
+tspan = (0.0, 3.7)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 100
+analysis_interval = 1000
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval = 20,
+save_solution = SaveSolutionCallback(interval=1000,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 1.3)
+stepsize_callback = StepsizeCallback(cfl = 0.9)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl

@@ -0,0 +1,96 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
+
+A version of the classical Kelvin-Helmholtz instability based on
+- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
+  A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
+  of the Euler Equations
+  [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
+"""
+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_kelvin_helmholtz_instability
+
+surface_flux = FluxLMARS(1.7)
+# polydeg = 3
+polydeg = 3
+basis = GaussLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3)
+#indicator_sc = IndicatorHennemannGassner(equations, basis,
+#                                         alpha_max = 0.000,
+#                                         alpha_min = 0.0000,
+#                                         alpha_smooth = true,
+#                                         variable = density_pressure)
+#
+#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+#                                                 volume_flux_dg = volume_flux,
+#                                                 volume_flux_fv = surface_flux)
+#volume_integral = VolumeIntegralWeakFormProjection()
+volume_integral = VolumeIntegralWeakForm()
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 100_000)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.7)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, 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)
+
+# stepsize_callback = StepsizeCallback(cfl = 0.5)
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+#stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-4, 5.0e-4),
+#                                                     variables = (Trixi.density, pressure))
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        #save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+#sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+#sol = solve(ode, SSPRK43(stage_limiter!),
+sol = solve(ode, SSPRK43(),
+            dt = 1.0e-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_kelvin_helmholtz_instability_sc_subcell.jl

@@ -32,7 +32,7 @@ function initial_condition_kelvin_helmholtz_instability(x, t,
 end
 initial_condition = initial_condition_kelvin_helmholtz_instability
 
-surface_flux = flux_lax_friedrichs
+surface_flux = flux_hllc
 volume_flux = flux_ranocha
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
@@ -48,7 +48,7 @@ solver = DGSEM(basis, surface_flux, volume_integral)
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 5,
+                initial_refinement_level = 6,
                 n_cells_max = 100_000)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)