More flexibility & consistency (across dimensions) to compute_coefficients! for discontinuous ICs

docs/literate/src/files/adding_nonconservative_equation.jl

@@ -145,7 +145,7 @@ plot(sol)
 # above.
 
 error_1 = analysis_callback(sol).l2 |> first
-@test isapprox(error_1, 0.00029609575838969394) #src
+@test isapprox(error_1, 0.0002961027497380263) #src
 # Next, we increase the grid resolution by one refinement level and run the
 # simulation again.
 

examples/p4est_3d_dgsem/elixir_euler_sedov.jl

@@ -6,15 +6,28 @@ using Trixi
 
 equations = CompressibleEulerEquations3D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for the Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes are at `[-1, 1]` are shifted to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionMediumSedovBlast <: DiscontinuousInitialCondition end
+
 """
-    initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D)
+    (initial_condition_medium_sedov_blast_wave::InitialConditionMediumSedovBlast)(x, t, 
+                                                                                  equations::CompressibleEulerEquations3D)
 
 The Sedov blast wave setup based on Flash
 - https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
 with smaller strength of the initial discontinuity.
 """
-function initial_condition_medium_sedov_blast_wave(x, t,
-                                                   equations::CompressibleEulerEquations3D)
+function (initial_condition_medium_sedov_blast_wave::InitialConditionMediumSedovBlast)(x, t,
+                                                                                       equations::CompressibleEulerEquations3D)
     # Set up polar coordinates
     inicenter = SVector(0.0, 0.0, 0.0)
     x_norm = x[1] - inicenter[1]
@@ -37,8 +50,9 @@ function initial_condition_medium_sedov_blast_wave(x, t,
 
     return prim2cons(SVector(rho, v1, v2, v3, p), equations)
 end
-
-initial_condition = initial_condition_medium_sedov_blast_wave
+# Note calling the constructor of the struct: `InitialConditionMediumSedovBlast()` instead of 
+# `initial_condition_medium_sedov_blast_wave` !
+const initial_condition = InitialConditionMediumSedovBlast()
 
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha

examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl

@@ -6,8 +6,21 @@ using Trixi
 
 advection_velocity = 1.0
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionComposite <: DiscontinuousInitialCondition end
+
 """
-    initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D)
+    (initial_condition_composite::InitialConditionComposite)(x, t, 
+                                                             equations::LinearScalarAdvectionEquation1D)
 
 Wave form that is a combination of a Gaussian pulse, a square wave, a triangle wave,
 and half an ellipse with periodic boundary conditions.
@@ -16,7 +29,8 @@ Slight simplification from
   Efficient Implementation of Weighted ENO Schemes
   [DOI: 10.1006/jcph.1996.0130](https://doi.org/10.1006/jcph.1996.0130)
 """
-function initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D)
+function (initial_condition_composite::InitialConditionComposite)(x, t,
+                                                                  equations::LinearScalarAdvectionEquation1D)
     xmin, xmax = -1.0, 1.0 # Only works if the domain is [-1.0,1.0]
     x_trans = x[1] - t
     L = xmax - xmin
@@ -42,8 +56,8 @@ function initial_condition_composite(x, t, equations::LinearScalarAdvectionEquat
 
     return SVector(value)
 end
-
-initial_condition = initial_condition_composite
+# Note calling the constructor of the struct: `InitialConditionComposite()` instead of `initial_condition_composite` !
+const initial_condition = InitialConditionComposite()
 
 equations = LinearScalarAdvectionEquation1D(advection_velocity)
 

examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl

@@ -14,23 +14,39 @@ cells_per_dimension = (64,)
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
                       periodicity = false)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionGreenlight <: DiscontinuousInitialCondition end
+
 # Example inspired from http://www.clawpack.org/riemann_book/html/Traffic_flow.html#Example:-green-light
 # Green light that at x = 0 which switches at t = 0 from red to green.
 # To the left there are cars bumper to bumper, to the right there are no cars.
-function initial_condition_greenlight(x, t, equation::TrafficFlowLWREquations1D)
+function (initial_condition_greenlight::InitialConditionGreenlight)(x, t,
+                                                                    equation::TrafficFlowLWREquations1D)
     RealT = eltype(x)
     scalar = x[1] < 0 ? one(RealT) : zero(RealT)
 
     return SVector(scalar)
 end
+# Note calling the constructor of the struct: `InitialConditionGreenlight()` instead of
+# `initial_condition_greenlight` !
+const initial_condition = InitialConditionGreenlight()
 
 ###############################################################################
 # Specify non-periodic boundary conditions
 
 # Assume that there are always cars waiting at the left
 function inflow(x, t, equations::TrafficFlowLWREquations1D)
     # -1.0 = coordinates_min
-    return initial_condition_greenlight(-1.0, t, equations)
+    return initial_condition(-1.0, t, equations)
 end
 boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
 
@@ -47,8 +63,6 @@ end
 boundary_conditions = (x_neg = boundary_condition_inflow,
                        x_pos = boundary_condition_outflow)
 
-initial_condition = initial_condition_greenlight
-
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_conditions)
 

examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl

@@ -32,18 +32,34 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
                 n_cells_max = 10_000,
                 periodicity = false)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionRarefaction <: DiscontinuousInitialCondition end
+
 # Discontinuous initial condition (Riemann Problem) leading to a rarefaction fan.
-function initial_condition_rarefaction(x, t, equation::InviscidBurgersEquation1D)
+function (initial_condition_rarefaction::InitialConditionRarefaction)(x, t,
+                                                                      equation::InviscidBurgersEquation1D)
     RealT = eltype(x)
     scalar = x[1] < 0.5f0 ? convert(RealT, 0.5f0) : convert(RealT, 1.5f0)
 
     return SVector(scalar)
 end
+# Note calling the constructor of the struct: `InitialConditionRarefaction()` instead of 
+# `initial_condition_rarefaction` !
+const initial_condition = InitialConditionRarefaction()
 
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_rarefaction)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,
@@ -57,8 +73,6 @@ end
 boundary_conditions = (x_neg = boundary_condition_inflow,
                        x_pos = boundary_condition_outflow)
 
-initial_condition = initial_condition_rarefaction
-
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_conditions)
 

examples/tree_1d_dgsem/elixir_burgers_shock.jl

@@ -32,18 +32,33 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
                 n_cells_max = 10_000,
                 periodicity = false)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionShock <: DiscontinuousInitialCondition end
+
 # Discontinuous initial condition (Riemann Problem) leading to a shock to test e.g. correct shock speed.
-function initial_condition_shock(x, t, equation::InviscidBurgersEquation1D)
+function (initial_condition_shock::InitialConditionShock)(x, t,
+                                                          equation::InviscidBurgersEquation1D)
     RealT = eltype(x)
     scalar = x[1] < 0.5f0 ? convert(RealT, 1.5f0) : convert(RealT, 0.5f0)
 
     return SVector(scalar)
 end
+# Note calling the constructor of the struct: `InitialConditionShock()` instead of `initial_condition_shock` !
+const initial_condition = InitialConditionShock()
 
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_shock)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,
@@ -57,8 +72,6 @@ end
 boundary_conditions = (x_neg = boundary_condition_inflow,
                        x_pos = boundary_condition_outflow)
 
-initial_condition = initial_condition_shock
-
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_conditions)
 

examples/tree_1d_dgsem/elixir_euler_blast_wave.jl

@@ -6,15 +6,29 @@ using Trixi
 
 equations = CompressibleEulerEquations1D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionBlastWave <: DiscontinuousInitialCondition end
+
 """
-    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+    (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                              equations::CompressibleEulerEquations1D)
 
 A medium blast wave taken from
 - Sebastian Hennemann, Gregor J. Gassner (2020)
   A provably entropy stable subcell shock capturing approach for high order split form DG
   [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
 """
-function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+function (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                                   equations::CompressibleEulerEquations1D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
     RealT = eltype(x)
@@ -34,7 +48,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation
 
     return prim2cons(SVector(rho, v1, p), equations)
 end
-initial_condition = initial_condition_blast_wave
+# Note calling the constructor of the struct: `InitialConditionBlastWave()` instead of
+# `initial_condition_blast_wave` !
+const initial_condition = InitialConditionBlastWave()
 
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha

examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl

@@ -6,15 +6,29 @@ using Trixi
 
 equations = CompressibleEulerEquations1D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionBlastWave <: DiscontinuousInitialCondition end
+
 """
-    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+    (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                              equations::CompressibleEulerEquations1D)
 
 A medium blast wave taken from
 - Sebastian Hennemann, Gregor J. Gassner (2020)
   A provably entropy stable subcell shock capturing approach for high order split form DG
   [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
 """
-function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
+function (initial_condition_blast_wave::InitialConditionBlastWave)(x, t,
+                                                                   equations::CompressibleEulerEquations1D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
     RealT = eltype(x)
@@ -34,7 +48,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation
 
     return prim2cons(SVector(rho, v1, p), equations)
 end
-initial_condition = initial_condition_blast_wave
+# Note calling the constructor of the struct: `InitialConditionBlastWave()` instead of
+# `initial_condition_blast_wave` !
+const initial_condition = InitialConditionBlastWave()
 
 # Note: We do not need to use the shock-capturing methodology here,
 # in contrast to the standard `euler_blast_wave.jl` example.

examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl

@@ -7,17 +7,30 @@ using Trixi
 
 equations = CompressibleEulerEquationsQuasi1D(1.4)
 
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousInitialCondition` for more information) which comes with a specialized 
+# initialization routine suited for Riemann problems.
+# In short, if a discontinuity is right at an interface, the boundary nodes (which are at the same location)
+# on that interface will be initialized with the left and right state of the discontinuity, i.e., 
+#                         { u_1, if element = left element and x_{element}^{(n)} = x_jump
+# u(x_jump, t, element) = {
+#                         { u_2, if element = right element and x_{element}^{(1)} = x_jump
+# This is realized by shifting the outer DG nodes inwards, i.e., on reference element
+# the outer nodes at `[-1, 1]` are shifted inwards to `[-1 + ε, 1 - ε]` with machine precision `ε`.
+struct InitialConditionDiscontinuity <: DiscontinuousInitialCondition end
+
 """
-    initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D)
+    (initial_condition_discontinuity::InitialConditionDiscontinuity)(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)
+function (initial_condition_discontinuity::InitialConditionDiscontinuity)(x, t,
+                                                                          equations::CompressibleEulerEquationsQuasi1D)
     RealT = eltype(x)
     rho = (x[1] < 0) ? RealT(3.4718) : RealT(2.0)
     v1 = (x[1] < 0) ? RealT(-2.5923) : RealT(-3.0)
@@ -26,8 +39,9 @@ function initial_condition_discontinuity(x, t,
 
     return prim2cons(SVector(rho, v1, p, a), equations)
 end
-
-initial_condition = initial_condition_discontinuity
+# Note calling the constructor of the struct: `InitialConditionDiscontinuity()` instead of
+# `initial_condition_discontinuity` !
+const initial_condition = InitialConditionDiscontinuity()
 
 surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal)
 volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)