Merge branch 'main' into entropy2cons_bug_fix

trixi-framework · Dec 14, 2023 · 2b43473 · 2b43473
2 parents bbbe07b + f191444
commit 2b43473
Show file tree

Hide file tree

Showing 12 changed files with 716 additions and 17 deletions.
diff --git a/Project.toml b/Project.toml
@@ -6,6 +6,7 @@ version = "0.6.5-pre"
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
@@ -52,6 +53,7 @@ TrixiMakieExt = "Makie"
 [compat]
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
+DataStructures = "0.18.15"
 DiffEqCallbacks = "2.25"
 EllipsisNotation = "1.0"
 FillArrays = "0.13.2, 1"

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
@@ -0,0 +1,85 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+"""
+    initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D)
+
+A discontinuous initial condition taken from
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
+    High order entropy stable schemes for the quasi-one-dimensional
+    shallow water and compressible Euler equations
+    [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)   
+"""
+function initial_condition_discontinuity(x, t,
+                                         equations::CompressibleEulerEquationsQuasi1D)
+    rho = (x[1] < 0) ? 3.4718 : 2.0
+    v1 = (x[1] < 0) ? -2.5923 : -3.0
+    p = (x[1] < 0) ? 5.7118 : 2.639
+    a = (x[1] < 0) ? 1.0 : 1.5
+
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+initial_condition = initial_condition_discontinuity
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal)
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_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.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_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 = 100
+
+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.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        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
diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
@@ -0,0 +1,73 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations with a discontinuous nozzle width function.
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+# Setup a truly discontinuous density function and nozzle width for
+# this academic testcase of entropy conservation. The errors from the analysis
+# callback are not important but the entropy error for this test case
+# `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff.
+# Works as intended for TreeMesh1D with `initial_refinement_level=6`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D)
+    v1 = 0.1
+    rho = 2.0 + 0.1 * x[1]
+    p = 3.0
+    a = 2.0 + x[1]
+
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+initial_condition = initial_condition_ec
+
+surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+volume_flux = surface_flux
+solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0,)
+coordinates_max = (1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# 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)
+
+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 = 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
diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
@@ -0,0 +1,60 @@
+using OrdinaryDiffEq
+using Trixi
+using ForwardDiff
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+volume_flux = surface_flux
+solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = -1.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# 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,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+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 = 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
diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl
@@ -67,7 +67,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval = 100,
+save_solution = SaveSolutionCallback(dt = 0.1,
                                      save_initial_solution = true,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)

diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -37,7 +37,8 @@ using SciMLBase: CallbackSet, DiscreteCallback,
 import SciMLBase: get_du, get_tmp_cache, u_modified!,
                   AbstractODEIntegrator, init, step!, check_error,
                   get_proposed_dt, set_proposed_dt!,
-                  terminate!, remake
+                  terminate!, remake, add_tstop!, has_tstop, first_tstop
+
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
@@ -70,6 +71,7 @@ using TriplotBase: TriplotBase
 using TriplotRecipes: DGTriPseudocolor
 @reexport using SimpleUnPack: @unpack
 using SimpleUnPack: @pack!
+using DataStructures: BinaryHeap, FasterForward, extract_all!
 
 # finite difference SBP operators
 using SummationByPartsOperators: AbstractDerivativeOperator,
@@ -139,6 +141,7 @@ export AcousticPerturbationEquations2D,
        CompressibleEulerEquations3D,
        CompressibleEulerMulticomponentEquations1D,
        CompressibleEulerMulticomponentEquations2D,
+       CompressibleEulerEquationsQuasi1D,
        IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D,
        IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D,