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_2d_dgsem/elixir_euler_riemannproblem_quadrants_amr.jl

@@ -0,0 +1,128 @@
+using OrdinaryDiffEqSSPRK
+using Trixi
+
+###############################################################################
+## Semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Variant of the 4-quadrant Riemann problem considered in 
+# - Carsten W. Schulz-Rinne:
+#   Classification of the Riemann Problem for Two-Dimensional Gas Dynamics
+#   https://doi.org/10.1137/0524006
+# and 
+# - Carsten W. Schulz-Rinne, James P. Collins, and Harland M. Glaz
+#   Numerical Solution of the Riemann Problem for Two-Dimensional Gas Dynamics
+#   https://doi.org/10.1137/0914082
+
+# Specify the initial condition as a discontinuous initial condition (see docstring of 
+# `DiscontinuousFunction` 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 InitialConditionRP <: DiscontinuousFunction end
+
+function (initial_condition_rp::InitialConditionRP)(x_, t,
+                                                    equations::CompressibleEulerEquations2D)
+    x, y = x_[1], x_[2]
+
+    if x >= 0.5 && y >= 0.5
+        rho, v1, v2, p = (0.5313, 0.0, 0.0, 0.4)
+    elseif x < 0.5 && y >= 0.5
+        rho, v1, v2, p = (1.0, 0.7276, 0.0, 1.0)
+    elseif x < 0.5 && y < 0.5
+        rho, v1, v2, p = (0.8, 0.0, 0.0, 1.0)
+    elseif x >= 0.5 && y < 0.5
+        rho, v1, v2, p = (1.0, 0.0, 0.7276, 1.0)
+    end
+
+    prim = SVector(rho, v1, v2, p)
+    return prim2cons(prim, equations)
+end
+# Note calling the constructor of the struct: `InitialConditionRP()` instead of `initial_condition_rp` !
+initial_condition = InitialConditionRP()
+
+# Extend domain by specifying free outflow
+@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                            surface_flux_function,
+                                            equations::CompressibleEulerEquations2D)
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_outflow,
+                           :x_pos => boundary_condition_outflow,
+                           :y_neg => boundary_condition_outflow,
+                           :y_pos => boundary_condition_outflow)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+trees_per_dimension = (4, 4) # Set base discretization relatively coarse to allow large cells
+mesh = P4estMesh(trees_per_dimension, polydeg = 1,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, false),
+                 initial_refinement_level = 2)
+
+surface_flux = flux_hllc
+volume_flux = flux_ranocha
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+shock_indicator = IndicatorHennemannGassner(equations, basis,
+                                            alpha_max = 0.5,
+                                            alpha_min = 0.001,
+                                            alpha_smooth = true,
+                                            variable = Trixi.density)
+
+volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+## ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.25)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.8)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 3, med_threshold = 0.02,
+                                      max_level = 6, max_threshold = 0.04)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 10,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+## Run the simulation
+
+sol = solve(ode, SSPRK54();
+            dt = 1.0, save_everystep = false, callback = callbacks);

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 
+# `DiscontinuousFunction` 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 <: DiscontinuousFunction 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` !
+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 
+# `DiscontinuousFunction` 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 <: DiscontinuousFunction 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` !
+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 
+# `DiscontinuousFunction` 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 <: DiscontinuousFunction 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` !
+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 
+# `DiscontinuousFunction` 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 <: DiscontinuousFunction 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` !
+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 
+# `DiscontinuousFunction` 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 <: DiscontinuousFunction 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` !
+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)