Merge branch 'main' into jc/StartUpDG_v1.1.0

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.8.4-DEV"
+version = "0.8.5-DEV"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

README.md

@@ -4,7 +4,7 @@
 [![Docs-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://trixi-framework.github.io/Trixi.jl/dev)
 [![Slack](https://img.shields.io/badge/chat-slack-e01e5a)](https://join.slack.com/t/trixi-framework/shared_invite/zt-sgkc6ppw-6OXJqZAD5SPjBYqLd8MU~g)
 [![Youtube](https://img.shields.io/youtube/channel/views/UCpd92vU2HjjTPup-AIN0pkg?style=social)](https://www.youtube.com/@trixi-framework)
-[![Build Status](https://github.com/trixi-framework/Trixi.jl/workflows/CI/badge.svg)](https://github.com/trixi-framework/Trixi.jl/actions?query=workflow%3ACI)
+[![Build Status](https://github.com/trixi-framework/Trixi.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/trixi-framework/Trixi.jl/actions?query=workflow%3ACI)
 [![Codecov](https://codecov.io/gh/trixi-framework/Trixi.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/trixi-framework/Trixi.jl)
 [![Coveralls](https://coveralls.io/repos/github/trixi-framework/Trixi.jl/badge.svg?branch=main)](https://coveralls.io/github/trixi-framework/Trixi.jl?branch=main)
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)

benchmark/benchmarks.jl

@@ -10,7 +10,13 @@ using Trixi
 
 const SUITE = BenchmarkGroup()
 
-for elixir in [joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_extended.jl"),
+for elixir in [# 1D
+               joinpath(examples_dir(), "structured_1d_dgsem", "elixir_euler_sedov.jl"),
+               joinpath(examples_dir(), "tree_1d_dgsem", "elixir_mhd_ec.jl"), 
+               joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls_amr.jl"),
+               joinpath(examples_dir(), "tree_1d_dgsem", "elixir_shallowwater_well_balanced_nonperiodic.jl"),
+               # 2D
+               joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_extended.jl"),
                joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_amr_nonperiodic.jl"),
                joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_ec.jl"),
                joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_vortex_mortar.jl"),
@@ -27,6 +33,7 @@ for elixir in [joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_exten
                joinpath(@__DIR__, "elixir_2d_euler_vortex_structured.jl"),
                joinpath(@__DIR__, "elixir_2d_euler_vortex_unstructured.jl"),
                joinpath(@__DIR__, "elixir_2d_euler_vortex_p4est.jl"),
+               # 3D
                joinpath(examples_dir(), "tree_3d_dgsem", "elixir_advection_extended.jl"),
                joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_ec.jl"),
                joinpath(examples_dir(), "tree_3d_dgsem", "elixir_euler_mortar.jl"),

docs/src/overview.md

@@ -60,15 +60,15 @@ different features on different mesh types.
 | Flux differencing                                            |          ✅         |             ✅            |               ✅              |          ✅          |               ✅            | [`VolumeIntegralFluxDifferencing`](@ref)
 | Shock capturing                                              |          ✅         |             ✅            |               ✅              |          ✅          |               ❌            | [`VolumeIntegralShockCapturingHG`](@ref)
 | Nonconservative equations                                    |          ✅         |             ✅            |               ✅              |          ✅          |               ✅            | e.g., GLM MHD or shallow water equations
-| Parabolic terms                                              |          ✅         |             ✅            |               ❌              |          ✅          |               ✅            | e.g., [`CompressibleNavierStokesDiffusion2D`](@ref)
+| Parabolic terms                                              |          ✅         |             ❌            |               ❌              |          ✅          |               ✅            | e.g., [`CompressibleNavierStokesDiffusion2D`](@ref)
 
 ᵃ: quad = quadrilateral, hex = hexahedron
 
-Note that except for [`TreeMesh`](@ref) all meshes are of *curvilinear* type, 
-which means that a (unit) vector normal to the interface (`normal_direction`) needs to be supplied to the 
+Note that except for [`TreeMesh`](@ref) all meshes are of *curvilinear* type,
+which means that a (unit) vector normal to the interface (`normal_direction`) needs to be supplied to the
 numerical flux function.
-You can check the [reference](https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/) if a certain 
-numerical flux is implemented with a `normal_direction` 
+You can check the [reference](https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/) if a certain
+numerical flux is implemented with a `normal_direction`
 or if currently only the *Cartesian* version (for [`TreeMesh`](@ref)) exists.
 In this case, you can still use this flux on curvilinear meshes by rotating it, see [`FluxRotated`](@ref).
 

examples/structured_1d_dgsem/elixir_euler_sedov.jl

@@ -38,7 +38,8 @@ initial_condition = initial_condition_sedov_blast_wave
 
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha
-basis = LobattoLegendreBasis(3)
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
 shock_indicator_variable = density_pressure
 indicator_sc = IndicatorHennemannGassner(equations, basis,
                                          alpha_max = 1.0,

src/equations/lattice_boltzmann_2d.jl

@@ -101,12 +101,16 @@ function LatticeBoltzmannEquations2D(; Ma, Re, collision_op = collision_bgk,
     # The relation between the isothermal speed of sound `c_s` and the mean thermal molecular velocity
     # `c` depends on the used phase space discretization, and is valid for D2Q9 (and others). For
     # details, see, e.g., [3] in the docstring above.
-    c_s = c / sqrt(3)
+    # c_s = c / sqrt(3) 
 
     # Calculate missing quantities
     if isnothing(Ma)
+        RealT = eltype(u0)
+        c_s = c / sqrt(convert(RealT, 3))
         Ma = u0 / c_s
     elseif isnothing(u0)
+        RealT = eltype(Ma)
+        c_s = c / sqrt(convert(RealT, 3))
         u0 = Ma * c_s
     end
     if isnothing(Re)
@@ -119,9 +123,10 @@ function LatticeBoltzmannEquations2D(; Ma, Re, collision_op = collision_bgk,
     Ma, Re, c, L, rho0, u0, nu = promote(Ma, Re, c, L, rho0, u0, nu)
 
     # Source for weights and speeds: [4] in the docstring above
-    weights = SVector(1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36, 1 / 36, 4 / 9)
-    v_alpha1 = SVector(c, 0, -c, 0, c, -c, -c, c, 0)
-    v_alpha2 = SVector(0, c, 0, -c, c, c, -c, -c, 0)
+    weights = SVector{9, RealT}(1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36,
+                                1 / 36, 4 / 9)
+    v_alpha1 = SVector{9, RealT}(c, 0, -c, 0, c, -c, -c, c, 0)
+    v_alpha2 = SVector{9, RealT}(0, c, 0, -c, c, c, -c, -c, 0)
 
     LatticeBoltzmannEquations2D(c, c_s, rho0, Ma, u0, Re, L, nu,
                                 weights, v_alpha1, v_alpha2,
@@ -154,7 +159,9 @@ A constant initial condition to test free-stream preservation.
 """
 function initial_condition_constant(x, t, equations::LatticeBoltzmannEquations2D)
     @unpack u0 = equations
-    rho = pi
+
+    RealT = eltype(x)
+    rho = convert(RealT, pi)
     v1 = u0
     v2 = u0
 
@@ -253,7 +260,7 @@ end
     else
         v_alpha = equations.v_alpha2
     end
-    return 0.5 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
+    return 0.5f0 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
 end
 
 """
@@ -358,7 +365,7 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
 @inline function collision_bgk(u, dt, equations::LatticeBoltzmannEquations2D)
     @unpack c_s, nu = equations
     tau = nu / (c_s^2 * dt)
-    return -(u - equilibrium_distribution(u, equations)) / (tau + 1 / 2)
+    return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
 end
 
 @inline have_constant_speed(::LatticeBoltzmannEquations2D) = True()

src/equations/lattice_boltzmann_3d.jl

@@ -141,12 +141,16 @@ function LatticeBoltzmannEquations3D(; Ma, Re, collision_op = collision_bgk,
     # The relation between the isothermal speed of sound `c_s` and the mean thermal molecular velocity
     # `c` depends on the used phase space discretization, and is valid for D3Q27 (and others). For
     # details, see, e.g., [3] in the docstring above.
-    c_s = c / sqrt(3)
+    # c_s = c / sqrt(3)
 
     # Calculate missing quantities
     if isnothing(Ma)
+        RealT = eltype(u0)
+        c_s = c / sqrt(convert(RealT, 3))
         Ma = u0 / c_s
     elseif isnothing(u0)
+        RealT = eltype(Ma)
+        c_s = c / sqrt(convert(RealT, 3))
         u0 = Ma * c_s
     end
     if isnothing(Re)
@@ -159,21 +163,24 @@ function LatticeBoltzmannEquations3D(; Ma, Re, collision_op = collision_bgk,
     Ma, Re, c, L, rho0, u0, nu = promote(Ma, Re, c, L, rho0, u0, nu)
 
     # Source for weights and speeds: [4] in docstring above
-    weights = SVector(2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 1 / 54, 1 / 54,
-                      1 / 54,
-                      1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54,
-                      1 / 54,
-                      1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216,
-                      1 / 216, 8 / 27)
-    v_alpha1 = SVector(c, -c, 0, 0, 0, 0, c, -c, c,
-                       -c, 0, 0, c, -c, c, -c, 0, 0,
-                       c, -c, c, -c, c, -c, -c, c, 0)
-    v_alpha2 = SVector(0, 0, c, -c, 0, 0, c, -c, 0,
-                       0, c, -c, -c, c, 0, 0, c, -c,
-                       c, -c, c, -c, -c, c, c, -c, 0)
-    v_alpha3 = SVector(0, 0, 0, 0, c, -c, 0, 0, c,
-                       -c, c, -c, 0, 0, -c, c, -c, c,
-                       c, -c, -c, c, c, -c, c, -c, 0)
+    weights = SVector{27, RealT}(2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 1 / 54,
+                                 1 / 54,
+                                 1 / 54,
+                                 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54,
+                                 1 / 54,
+                                 1 / 54,
+                                 1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216,
+                                 1 / 216,
+                                 1 / 216, 8 / 27)
+    v_alpha1 = SVector{27, RealT}(c, -c, 0, 0, 0, 0, c, -c, c,
+                                  -c, 0, 0, c, -c, c, -c, 0, 0,
+                                  c, -c, c, -c, c, -c, -c, c, 0)
+    v_alpha2 = SVector{27, RealT}(0, 0, c, -c, 0, 0, c, -c, 0,
+                                  0, c, -c, -c, c, 0, 0, c, -c,
+                                  c, -c, c, -c, -c, c, c, -c, 0)
+    v_alpha3 = SVector{27, RealT}(0, 0, 0, 0, c, -c, 0, 0, c,
+                                  -c, c, -c, 0, 0, -c, c, -c, c,
+                                  c, -c, -c, c, c, -c, c, -c, 0)
 
     LatticeBoltzmannEquations3D(c, c_s, rho0, Ma, u0, Re, L, nu,
                                 weights, v_alpha1, v_alpha2, v_alpha3,
@@ -206,7 +213,9 @@ A constant initial condition to test free-stream preservation.
 """
 function initial_condition_constant(x, t, equations::LatticeBoltzmannEquations3D)
     @unpack u0 = equations
-    rho = pi
+
+    RealT = eltype(x)
+    rho = convert(RealT, pi)
     v1 = u0
     v2 = u0
     v3 = u0
@@ -243,7 +252,7 @@ end
     else # z-direction
         v_alpha = equations.v_alpha3
     end
-    return 0.5 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
+    return 0.5f0 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
 end
 
 """
@@ -369,7 +378,7 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
 @inline function collision_bgk(u, dt, equations::LatticeBoltzmannEquations3D)
     @unpack c_s, nu = equations
     tau = nu / (c_s^2 * dt)
-    return -(u - equilibrium_distribution(u, equations)) / (tau + 1 / 2)
+    return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
 end
 
 @inline have_constant_speed(::LatticeBoltzmannEquations3D) = True()
@@ -391,7 +400,7 @@ end
     rho = density(u, equations)
     v1, v2, v3 = velocity(u, equations)
 
-    return 0.5 * (v1^2 + v2^2 + v3^2) / rho / equations.rho0
+    return 0.5f0 * (v1^2 + v2^2 + v3^2) / rho / equations.rho0
 end
 
 # Calculate nondimensionalized kinetic energy for a conservative state `u`

src/equations/maxwell_1d.jl

@@ -58,7 +58,7 @@ function initial_condition_convergence_test(x, t, equations::MaxwellEquations1D)
     c = equations.speed_of_light
     char_pos = c * t + x[1]
 
-    sin_char_pos = sin(2 * pi * char_pos)
+    sin_char_pos = sinpi(2 * char_pos)
 
     E = -c * sin_char_pos
     B = sin_char_pos