Shu-Osher testcase non-uniform Structured 1D mesh with FV O2 (#2651)

* Shu-Osher testcase non-uniform Structured 1D mesh with FV O2

* Apply suggestions from code review

* Update examples/structured_1d_dgsem/elixir_euler_shu_osher_nonuniform_fvO2.jl

* Update examples/structured_1d_dgsem/elixir_euler_shu_osher_nonuniform_fvO2.jl

* controller

* comment

* toml

* compat

* compat

* compat

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

Author: DanielDoehring

examples/structured_1d_dgsem/elixir_euler_shu_osher_nonuniform_fvO2.jl

@@ -0,0 +1,93 @@
+using OrdinaryDiffEqSSPRK
+using OrdinaryDiffEqCore: PIDController
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+# Shu-Osher initial condition for 1D compressible Euler equations
+# Example 8 from Shu, Osher (1989).
+# https://doi.org/10.1016/0021-9991(89)90222-2
+function initial_condition_shu_osher(x, t, equations::CompressibleEulerEquations1D)
+    x0 = -4
+
+    rho_left = 27 / 7
+    v_left = 4 * sqrt(35) / 9
+    p_left = 31 / 3
+
+    v_right = 0.0
+    p_right = 1.0
+
+    rho = ifelse(x[1] > x0, 1 + 1 / 5 * sin(5 * x[1]), rho_left)
+    v = ifelse(x[1] > x0, v_right, v_left)
+    p = ifelse(x[1] > x0, p_right, p_left)
+
+    return prim2cons(SVector(rho, v, p), equations)
+end
+
+initial_condition = initial_condition_shu_osher
+
+surface_flux = flux_hllc
+polydeg = 3 # governs in this case only the number of FV subcells per DG cell
+basis = LobattoLegendreBasis(polydeg)
+volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
+                                                      volume_flux_fv = surface_flux,
+                                                      reconstruction_mode = reconstruction_O2_inner,
+                                                      slope_limiter = monotonized_central)
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+function refined_mapping(xi)
+    fine_fraction = 0.8
+    a_fine = 4.0
+    x_max = 5.0
+
+    # Compute the value at the boundary for continuity
+    x_boundary = a_fine * fine_fraction
+    # Slope for the outer region
+    a_outer = (x_max - x_boundary) / (1 - fine_fraction)
+    if abs(xi) <= fine_fraction
+        x = a_fine * xi
+    else
+        x = x_boundary + a_outer * (abs(xi) - fine_fraction)
+        x *= sign(xi)
+    end
+    return x
+end
+
+cells_per_dimension = (128,)
+mesh = StructuredMesh(cells_per_dimension, refined_mapping, periodicity = false)
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = boundary_condition_default(mesh, boundary_condition)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+# Solve ODE with optimized timestep controller from https://doi.org/10.1007/s42967-021-00159-w
+sol = solve(ode, SSPRK43();
+            controller = PIDController(0.55, -0.27, 0.05),
+            abstol = 1e-4, reltol = 1e-4,
+            callback = callbacks, ode_default_options()...)

test/Project.toml

@@ -18,6 +18,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqFeagin = "101fe9f7-ebb6-4678-b671-3a81e7194747"
 OrdinaryDiffEqHighOrderRK = "d28bc4f8-55e1-4f49-af69-84c1a99f0f58"
 OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
@@ -56,6 +57,7 @@ LinearSolve = "2.36.1, 3"
 MPI = "0.20.6"
 NLsolve = "4.5.1"
 OrdinaryDiffEqBDF = "1.1"
+OrdinaryDiffEqCore = "1.6"
 OrdinaryDiffEqFeagin = "1"
 OrdinaryDiffEqHighOrderRK = "1.1"
 OrdinaryDiffEqLowOrderRK = "1.2"

test/test_structured_1d.jl

@@ -139,6 +139,18 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_euler_shu_osher_nonuniform_fvO2.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_shu_osher_nonuniform_fvO2.jl"),
+                        abstol=1e-11, reltol=1e-11,
+                        l2=[0.5000856244205995, 1.7958065770030094, 6.50409305116869],
+                        linf=[2.981082161129505, 10.10855083991654, 36.26666086388062],
+                        tspan=(0.0, 0.1))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_euler_source_terms.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!