Merge branch 'main' into PERK-documentation

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.9.5-DEV"
+version = "0.9.6-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl

@@ -99,8 +99,8 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                  volume_flux_fv = surface_flux)
 solver_euler = DGSEM(basis, surface_flux, volume_integral)
 
-coordinates_min = (-4, -4)
-coordinates_max = (4, 4)
+coordinates_min = (-4.0, -4.0)
+coordinates_max = (4.0, 4.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level = 2,
                 n_cells_max = 100_000,

examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl

@@ -8,8 +8,8 @@ equations = TrafficFlowLWREquations1D()
 
 solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis))
 
-coordinates_min = (-1.0,) # minimum coordinate
-coordinates_max = (1.0,) # maximum coordinate
+coordinates_min = -1.0 # minimum coordinate
+coordinates_max = 1.0 # maximum coordinate
 cells_per_dimension = (64,)
 
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,

examples/tree_1d_dgsem/elixir_advection_diffusion.jl

@@ -13,19 +13,21 @@ equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
-coordinates_min = -pi # minimum coordinate
-coordinates_max = pi # maximum coordinate
+coordinates_min = -convert(Float64, pi) # minimum coordinate
+coordinates_max = convert(Float64, pi) # maximum coordinate
 
 # Create a uniformly refined mesh with periodic boundaries
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level = 4,
                 n_cells_max = 30_000, # set maximum capacity of tree data structure
                 periodicity = true)
 
-function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0))
+function x_trans_periodic(x, domain_length = SVector(oftype(x[1], 2 * pi)),
+                          center = SVector(oftype(x[1], 0)))
     x_normalized = x .- center
     x_shifted = x_normalized .% domain_length
-    x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .*
+    x_offset = ((x_shifted .< -0.5f0 * domain_length) -
+                (x_shifted .> 0.5f0 * domain_length)) .*
                domain_length
     return center + x_shifted + x_offset
 end
@@ -37,11 +39,9 @@ function initial_condition_diffusive_convergence_test(x, t,
     x_trans = x_trans_periodic(x - equation.advection_velocity * t)
 
     nu = diffusivity()
-    c = 0.0
-    A = 1.0
-    L = 2
-    f = 1 / L
-    omega = 1.0
+    c = 0
+    A = 1
+    omega = 1
     scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t)
     return SVector(scalar)
 end
@@ -75,8 +75,12 @@ analysis_callback = AnalysisCallback(semi, interval = 100)
 # The AliveCallback prints short status information in regular intervals
 alive_callback = AliveCallback(analysis_interval = 100)
 
+# The SaveRestartCallback allows to save a file from which a Trixi.jl simulation can be restarted
+save_restart = SaveRestartCallback(interval = 100,
+                                   save_final_restart = true)
+
 # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
-callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_restart)
 
 ###############################################################################
 # run the simulation

examples/tree_1d_dgsem/elixir_advection_diffusion_restart.jl

@@ -0,0 +1,30 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# create a restart file
+
+elixir_file = "elixir_advection_diffusion.jl"
+trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file))
+
+###############################################################################
+# initialize the ODE
+
+restart_file = "restart_000000018.h5"
+restart_filename = joinpath("out", restart_file)
+tspan = (load_time(restart_filename), 2.0)
+
+ode = semidiscretize(semi, tspan, restart_filename);
+
+# Do not save restart files here
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, KenCarp4(autodiff = false), abstol = time_abs_tol, reltol = time_int_tol,
+            save_everystep = false, callback = callbacks)
+
+# Print the timer summary
+summary_callback()

examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl

@@ -8,8 +8,9 @@ using Trixi
 equations = InviscidBurgersEquation1D()
 
 function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D)
+    RealT = eltype(x)
     k = 1
-    2 + sinpi(k * (x[1] - 0.7)) |> SVector
+    SVector(2 + sinpi(k * (x[1] - convert(RealT, 0.7))))
 end
 
 volume_flux = flux_ec

examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl

@@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
 
 # Discontinuous initial condition (Riemann Problem) leading to a rarefaction fan.
 function initial_condition_rarefaction(x, t, equation::InviscidBurgersEquation1D)
-    scalar = x[1] < 0.5 ? 0.5 : 1.5
+    RealT = eltype(x)
+    scalar = x[1] < 0.5f0 ? convert(RealT, 0.5f0) : convert(RealT, 1.5f0)
 
     return SVector(scalar)
 end

examples/tree_1d_dgsem/elixir_burgers_shock.jl

@@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max,
 
 # Discontinuous initial condition (Riemann Problem) leading to a shock to test e.g. correct shock speed.
 function initial_condition_shock(x, t, equation::InviscidBurgersEquation1D)
-    scalar = x[1] < 0.5 ? 1.5 : 0.5
+    RealT = eltype(x)
+    scalar = x[1] < 0.5f0 ? convert(RealT, 1.5f0) : convert(RealT, 0.5f0)
 
     return SVector(scalar)
 end

examples/tree_1d_dgsem/elixir_euler_blast_wave.jl

@@ -18,7 +18,8 @@ A medium blast wave taken from
 function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
     # The following code is equivalent to
@@ -28,9 +29,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation
     cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)
 
     # Calculate primitive variables
-    rho = r > 0.5 ? 1.0 : 1.1691
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
-    p = r > 0.5 ? 1.0E-3 : 1.245
+    rho = r > 0.5f0 ? one(RealT) : RealT(1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi
+    p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245)
 
     return prim2cons(SVector(rho, v1, p), equations)
 end

examples/tree_1d_dgsem/elixir_euler_positivity.jl

@@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash
 """
 function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
 
     # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
-    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
     # r0 = 0.5 # = more reasonable setup
-    E = 1.0
-    p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
-    p0_outer = 1.0e-5 # = true Sedov setup
+    E = 1
+    p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0)
+    p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup
     # p0_outer = 1.0e-3 # = more reasonable setup
 
     # Calculate primitive variables
-    rho = 1.0
-    v1 = 0.0
+    rho = 1
+    v1 = 0
     p = r > r0 ? p0_outer : p0_inner
 
     return prim2cons(SVector(rho, v1, p), equations)

examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl

@@ -18,10 +18,11 @@ A discontinuous initial condition taken from
 """
 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
+    RealT = eltype(x)
+    rho = (x[1] < 0) ? RealT(3.4718) : RealT(2.0)
+    v1 = (x[1] < 0) ? RealT(-2.5923) : RealT(-3.0)
+    p = (x[1] < 0) ? RealT(5.7118) : RealT(2.639)
+    a = (x[1] < 0) ? 1.0f0 : 1.5f0
 
     return prim2cons(SVector(rho, v1, p, a), equations)
 end

examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl

@@ -15,10 +15,11 @@ equations = CompressibleEulerEquationsQuasi1D(1.4)
 # 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]
+    RealT = eltype(x)
+    v1 = convert(RealT, 0.1)
+    rho = 2 + convert(RealT, 0.1) * x[1]
+    p = 3
+    a = 2 + x[1]
 
     return prim2cons(SVector(rho, v1, p, a), equations)
 end

examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl

@@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash
 """
 function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
 
     # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
-    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
     # r0 = 0.5 # = more reasonable setup
-    E = 1.0
-    p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
-    p0_outer = 1.0e-5 # = true Sedov setup
+    E = 1
+    p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0)
+    p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup
     # p0_outer = 1.0e-3 # = more reasonable setup
 
     # Calculate primitive variables
-    rho = 1.0
-    v1 = 0.0
+    rho = 1
+    v1 = 0
     p = r > r0 ? p0_outer : p0_inner
 
     return prim2cons(SVector(rho, v1, p), equations)

examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl

@@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash
 """
 function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
     # Set up polar coordinates
-    inicenter = SVector(0.0)
+    RealT = eltype(x)
+    inicenter = SVector(0)
     x_norm = x[1] - inicenter[1]
     r = abs(x_norm)
 
     # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
-    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
     # r0 = 0.5 # = more reasonable setup
-    E = 1.0
-    p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0)
-    p0_outer = 1.0e-5 # = true Sedov setup
+    E = 1
+    p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0)
+    p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup
     # p0_outer = 1.0e-3 # = more reasonable setup
 
     # Calculate primitive variables
-    rho = 1.0
-    v1 = 0.0
+    rho = 1
+    v1 = 0
     p = r > r0 ? p0_outer : p0_inner
 
     return prim2cons(SVector(rho, v1, p), equations)

examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl

@@ -18,20 +18,21 @@ A multicomponent two interacting blast wave test taken from
 """
 function initial_condition_two_interacting_blast_waves(x, t,
                                                        equations::CompressibleEulerMulticomponentEquations1D)
-    rho1 = 0.5 * x[1]^2
-    rho2 = 0.5 * (sin(20 * x[1]))^2
+    RealT = eltype(x)
+    rho1 = 0.5f0 * x[1]^2
+    rho2 = 0.5f0 * (sin(20 * x[1]))^2
     rho3 = 1 - rho1 - rho2
 
     prim_rho = SVector{3, real(equations)}(rho1, rho2, rho3)
 
-    v1 = 0.0
+    v1 = 0
 
-    if x[1] <= 0.1
-        p = 1000
-    elseif x[1] < 0.9
-        p = 0.01
+    if x[1] <= RealT(0.1)
+        p = convert(RealT, 1000)
+    elseif x[1] < RealT(0.9)
+        p = convert(RealT, 0.01)
     else
-        p = 100
+        p = convert(RealT, 100)
     end
 
     prim_other = SVector{2, real(equations)}(v1, p)

examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl

@@ -19,12 +19,13 @@ A non-priodic harmonic function used in combination with
 function initial_condition_harmonic_nonperiodic(x, t,
                                                 equations::HyperbolicDiffusionEquations1D)
     # elliptic equation: -νΔϕ = f
-    if t == 0.0
-        phi = 5.0
-        q1 = 0.0
+    RealT = eltype(x)
+    if t == 0
+        phi = convert(RealT, 5)
+        q1 = zero(RealT)
     else
         A = 3
-        B = exp(1)
+        B = exp(one(RealT))
         phi = A + B * x[1]
         q1 = B
     end

examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl

@@ -22,7 +22,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 # that is advected to left with v - c and to the right with v + c.
 # Correspondigly, the bump splits in half.
 function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations1D)
-    v1_prime = 0.0
+    v1_prime = 0
     rho_prime = p_prime = 2 * exp(-(x[1] - 45)^2 / 25)
     return SVector(rho_prime, v1_prime, p_prime)
 end

examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl

@@ -20,9 +20,10 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 # Excite the electric field which causes a standing wave.
 # The solution is an undamped exchange between electric and magnetic energy.
 function initial_condition_E_excitation(x, t, equations::MaxwellEquations1D)
+    RealT = eltype(x)
     c = equations.speed_of_light
-    E = -c * sin(2 * pi * x[1])
-    B = 0.0
+    E = -c * sinpi(2 * x[1])
+    B = 0
 
     return SVector(E, B)
 end

examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl

@@ -18,14 +18,15 @@ MHD extension of the Sod shock tube. Taken from Section V of the article
 """
 function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
     # domain must be set to [0, 1], γ = 2, final time = 0.12
-    rho = x[1] < 0.5 ? 1.0 : 0.125
-    v1 = 0.0
-    v2 = 0.0
-    v3 = 0.0
-    p = x[1] < 0.5 ? 1.0 : 0.1
-    B1 = 0.75
-    B2 = x[1] < 0.5 ? 1.0 : -1.0
-    B3 = 0.0
+    RealT = eltype(x)
+    rho = x[1] < 0.5f0 ? 1.0f0 : 0.125f0
+    v1 = 0
+    v2 = 0
+    v3 = 0
+    p = x[1] < 0.5f0 ? one(RealT) : RealT(0.1)
+    B1 = 0.75f0
+    B2 = x[1] < 0.5f0 ? 1 : -1
+    B3 = 0
     return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
 end
 initial_condition = initial_condition_briowu_shock_tube