Merge branch 'main' into jc/p4est_parabolic_BCs

NEWS.md

@@ -9,6 +9,7 @@ for human readability.
 #### Added
 
 - Experimental support for 3D parabolic diffusion terms has been added.
+- Capability to set truly discontinuous initial conditions in 1D.
 
 #### Changed
 

docs/literate/src/files/adding_nonconservative_equation.jl

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

examples/dgmulti_3d/elixir_advection_tensor_wedge.jl

@@ -0,0 +1,56 @@
+using OrdinaryDiffEq
+using Trixi
+using LinearAlgebra
+
+###############################################################################
+equations = LinearScalarAdvectionEquation3D(1.0, 1.0, 1.0)
+
+initial_condition = initial_condition_convergence_test
+
+# Define the polynomial degrees for the polynoms of the triangular base and the line
+# of the tensor-prism
+tensor_polydeg = (3, 4)
+
+dg = DGMulti(element_type = Wedge(),
+             approximation_type = Polynomial(),
+             surface_flux = flux_lax_friedrichs,
+             polydeg = tensor_polydeg)
+
+
+cells_per_dimension = (8, 8, 8)
+mesh = DGMultiMesh(dg, 
+                   cells_per_dimension,
+                   coordinates_min = (-1.0, -1.0, -1.0), 
+                   coordinates_max = (1.0, 1.0, 1.0),
+                   periodicity = true)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                boundary_conditions=boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0,
+            save_everystep=false, callback=callbacks);
+
+summary_callback() # print the timer summary

examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl

@@ -9,7 +9,9 @@ advection_velocity = 1.0
 """
     initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D)
 
-Slight simplification of
+Wave form that is a combination of a Gaussian pulse, a square wave, a triangle wave,
+and half an ellipse with periodic boundary conditions.
+Slight simplification from
 - Jiang, Shu (1996)
   Efficient Implementation of Weighted ENO Schemes
   [DOI: 10.1006/jcph.1996.0130](https://doi.org/10.1006/jcph.1996.0130)
@@ -60,7 +62,7 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
 solver = DGSEM(basis, surface_flux, volume_integral)
 
 # Create curved mesh
-cells_per_dimension = (125,)
+cells_per_dimension = (120,)
 coordinates_min = (-1.0,) # minimum coordinate
 coordinates_max = (1.0,)  # maximum coordinate
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,

examples/tree_1d_dgsem/elixir_burgers_shock.jl

@@ -21,7 +21,7 @@ surface_flux = flux_lax_friedrichs
 volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                  volume_flux_dg=surface_flux,
                                                  volume_flux_fv=surface_flux)
-                                                 
+
 solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinate_min = 0.0
@@ -59,7 +59,7 @@ 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,
@@ -79,7 +79,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval=analysis_interval)
 
-stepsize_callback = StepsizeCallback(cfl=0.9)
+stepsize_callback = StepsizeCallback(cfl=0.85)
 
 
 callbacks = CallbackSet(summary_callback,

examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl

@@ -88,7 +88,7 @@ ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 100
+analysis_interval = 1000
 
 analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
 

examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl

@@ -94,7 +94,7 @@ amr_callback = AMRCallback(semi, amr_controller,
                            adapt_initial_condition=true,
                            adapt_initial_condition_only_refine=true)
 
-stepsize_callback = StepsizeCallback(cfl=0.8)
+stepsize_callback = StepsizeCallback(cfl=0.65)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl

@@ -4,8 +4,8 @@ using Trixi
 
 ###############################################################################
 # semidiscretization of the ideal MHD equations
-equations = IdealGlmMhdMulticomponentEquations1D(gammas           = (2.0, 2.0, 2.0),
-                                                 gas_constants    = (2.0, 2.0, 2.0))
+equations = IdealGlmMhdMulticomponentEquations1D(gammas        = (2.0, 2.0, 2.0),
+                                                 gas_constants = (2.0, 2.0, 2.0))
 
 initial_condition = initial_condition_weak_blast_wave
 

examples/tree_1d_dgsem/elixir_mhdmulti_es.jl

@@ -4,8 +4,8 @@ using Trixi
 
 ###############################################################################
 # semidiscretization of the ideal MHD equations
-equations = IdealGlmMhdMulticomponentEquations1D(gammas           = (2.0, 2.0, 2.0),
-                                                 gas_constants    = (2.0, 2.0, 2.0))
+equations = IdealGlmMhdMulticomponentEquations1D(gammas        = (2.0, 2.0, 2.0),
+                                                 gas_constants = (2.0, 2.0, 2.0))
 
 initial_condition = initial_condition_weak_blast_wave
 

examples/tree_1d_dgsem/elixir_shallowwater_ec.jl

@@ -8,9 +8,34 @@ using Trixi
 
 equations = ShallowWaterEquations1D(gravity_constant=9.81)
 
-# Note, this initial condition is used to compute errors in the analysis callback but the initialization is
-# overwritten by `initial_condition_ec_discontinuous_bottom` below.
-initial_condition = initial_condition_weak_blast_wave
+# Initial condition with a truly discontinuous water height, velocity, and bottom
+# topography function as an academic testcase for entropy conservation.
+# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should
+# be around machine roundoff.
+# Works as intended for TreeMesh1D with `initial_refinement_level=4`. 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_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D)
+  # Set the background values
+  H = 4.25
+  v = 0.0
+  b = sin(x[1]) # arbitrary continuous function
+
+  # Setup the discontinuous water height and velocity
+  if x[1] >= 0.125 && x[1] <= 0.25
+    H = 5.0
+    v = 0.1882
+  end
+
+  # Setup a discontinuous bottom topography
+  if x[1] >= -0.25 && x[1] <= -0.125
+    b = 2.0 + 0.5 * sin(2.0 * pi * x[1])
+  end
+
+  return prim2cons(SVector(H, v, b), equations)
+end
+
+initial_condition = initial_condition_ec_discontinuous_bottom
 
 ###############################################################################
 # Get the DG approximation space
@@ -37,46 +62,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 tspan = (0.0, 2.0)
 ode = semidiscretize(semi, tspan)
 
-###############################################################################
-# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.
-
-# alternative version of the initial conditinon used to setup a truly discontinuous
-# bottom topography function and initial condition for this academic testcase of entropy conservation.
-# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff
-# In contrast to the usual signature of initial conditions, this one get passed the
-# `element_id` explicitly. In particular, this initial conditions works as intended
-# only for the TreeMesh1D with `initial_refinement_level=4`.
-function initial_condition_ec_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D)
-  # Set the background values
-  H = 4.25
-  v = 0.0
-  b = sin(x[1]) # arbitrary continuous function
-
-  # setup the discontinuous water height and velocity
-  if element_id == 10
-    H = 5.0
-    v = 0.1882
-  end
-
-  # Setup a discontinuous bottom topography using the element id number
-  if element_id == 7
-    b = 2.0 + 0.5 * sin(2.0 * pi * x[1])
-  end
-
-  return prim2cons(SVector(H, v, b), equations)
-end
-
-# point to the data we want to augment
-u = Trixi.wrap_array(ode.u0, semi)
-# reset the initial condition
-for element in eachelement(semi.solver, semi.cache)
-  for i in eachnode(semi.solver)
-    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element)
-    u_node = initial_condition_ec_discontinuous_bottom(x_node, first(tspan), element, equations)
-    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
-  end
-end
-
 ###############################################################################
 # Callbacks
 

examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl

@@ -7,24 +7,37 @@ using Trixi
 
 equations = ShallowWaterEquations1D(gravity_constant=9.812, H0=1.75)
 
-function initial_condition_stone_throw(x, t, equations::ShallowWaterEquations1D)
-    # Set up polar coordinates
-    inicenter = 0.15
-    x_norm = x[1] - inicenter[1]
-    r = abs(x_norm)
+# Initial condition with a truly discontinuous velocity and bottom topography.
+# Works as intended for TreeMesh1D with `initial_refinement_level=3`. 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_stone_throw_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D)
+
+  # Calculate primitive variables
+
+  # flat lake
+  H = equations.H0
+
+  # Discontinuous velocity
+  v = 0.0
+  if x[1] >= -0.75 && x[1] <= 0.0
+      v = -1.0
+  elseif x[1] >= 0.0 && x[1] <= 0.75
+      v = 1.0
+  end
 
-    # Calculate primitive variables
-    H = equations.H0
-    # v = 0.0 # for well-balanced test
-    v = r < 0.6 ? 1.75 : 0.0 # for stone throw
+  b = (  1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) )
+     + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) )
 
-    b = (  1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) )
-       + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) )
+  # Force a discontinuous bottom topography
+  if x[1] >= -1.5 && x[1] <= 0.0
+    b = 0.5
+  end
 
-    return prim2cons(SVector(H, v, b), equations)
+  return prim2cons(SVector(H, v, b), equations)
 end
 
-initial_condition = initial_condition_stone_throw
+initial_condition = initial_condition_stone_throw_discontinuous_bottom
 
 boundary_condition = boundary_condition_slip_wall
 
@@ -62,49 +75,13 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_condition)
 
 ###############################################################################
-# ODE solvers, callbacks etc.
+# ODE solver
 
 tspan = (0.0, 3.0)
 ode = semidiscretize(semi, tspan)
 
-# Hack in a discontinuous bottom for a more interesting test
-function initial_condition_stone_throw_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D)
-
-    inicenter = 0.15
-    x_norm = x[1] - inicenter[1]
-    r = abs(x_norm)
-
-    # Calculate primitive variables
-    H = equations.H0 # flat lake
-    # Discontinuous velocity set via element id number
-    v = 0.0
-    if element_id == 4
-        v = -1.0
-    elseif element_id == 5
-        v = 1.0
-    end
-
-    b = (  1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) )
-       + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) )
-
-    # Setup a discontinuous bottom topography using the element id number
-    if element_id == 3 || element_id == 4
-      b = 0.5
-    end
-
-    return prim2cons(SVector(H, v, b), equations)
-end
-
-# point to the data we want to augment
-u = Trixi.wrap_array(ode.u0, semi)
-# reset the initial condition
-for element in eachelement(semi.solver, semi.cache)
-  for i in eachnode(semi.solver)
-    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element)
-    u_node = initial_condition_stone_throw_discontinuous_bottom(x_node, first(tspan), element, equations)
-    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
-  end
-end
+###############################################################################
+# Callbacks
 
 summary_callback = SummaryCallback()