diff --git a/Project.toml b/Project.toml index bb3bea6aa0..6caf69884c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.8.4-DEV" +version = "0.8.5-DEV" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/README.md b/README.md index 57b99238cc..a76fa26843 100644 --- a/README.md +++ b/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) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 0d6fabcd4a..e551e6e9cb 100644 --- a/benchmark/benchmarks.jl +++ b/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"), diff --git a/docs/src/overview.md b/docs/src/overview.md index 72a341de7e..42e8691142 100644 --- a/docs/src/overview.md +++ b/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). diff --git a/examples/structured_1d_dgsem/elixir_euler_sedov.jl b/examples/structured_1d_dgsem/elixir_euler_sedov.jl index 4ffa100cc7..fc974cb94d 100644 --- a/examples/structured_1d_dgsem/elixir_euler_sedov.jl +++ b/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, diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl index 272dd897ce..bad417d9c8 100644 --- a/src/equations/lattice_boltzmann_2d.jl +++ b/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() diff --git a/src/equations/lattice_boltzmann_3d.jl b/src/equations/lattice_boltzmann_3d.jl index d3eada15f5..bf4c365e0f 100644 --- a/src/equations/lattice_boltzmann_3d.jl +++ b/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` diff --git a/src/equations/maxwell_1d.jl b/src/equations/maxwell_1d.jl index 838e4e7b81..096fe3946c 100644 --- a/src/equations/maxwell_1d.jl +++ b/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 diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 7007bea887..998deb04de 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -71,16 +71,16 @@ A smooth initial condition used for convergence tests in combination with [`source_terms_convergence_test`](@ref) (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). """ - function initial_condition_convergence_test(x, t, equations::ShallowWaterEquations1D) # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] - c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi + RealT = eltype(x) + c = 7 + omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2)) + omega_t = 2 * convert(RealT, pi) H = c + cos(omega_x * x[1]) * cos(omega_t * t) - v = 0.5 - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x[1]) + v = 0.5f0 + b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x[1]) return prim2cons(SVector(H, v, b), equations) end @@ -92,19 +92,19 @@ Source terms used for convergence tests in combination with (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). This manufactured solution source term is specifically designed for the bottom topography function -`b(x) = 2.0 + 0.5 * sin(sqrt(2.0)*pi*x[1])` +`b(x) = 2.0 + 0.5 * sinpi(sqrt(2.0) * x[1])` as defined in [`initial_condition_convergence_test`](@ref). """ - @inline function source_terms_convergence_test(u, x, t, equations::ShallowWaterEquations1D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocity is taken to be constant - c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi - omega_b = sqrt(2.0) * pi - v = 0.5 + RealT = eltype(u) + c = 7 + omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2)) + omega_t = 2 * convert(RealT, pi) + omega_b = sqrt(convert(RealT, 2)) * convert(RealT, pi) + v = 0.5f0 sinX, cosX = sincos(omega_x * x[1]) sinT, cosT = sincos(omega_t * t) @@ -116,12 +116,12 @@ as defined in [`initial_condition_convergence_test`](@ref). H_t = -omega_t * cosX * sinT # bottom topography and its spatial derivative - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x[1]) - b_x = 0.5 * omega_b * cos(omega_b * x[1]) + b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x[1]) + b_x = 0.5f0 * omega_b * cos(omega_b * x[1]) du1 = H_t + v * (H_x - b_x) du2 = v * du1 + equations.gravity * (H - b) * H_x - return SVector(du1, du2, 0.0) + return SVector(du1, du2, 0) end """ @@ -131,13 +131,14 @@ A weak blast wave discontinuity useful for testing, e.g., total energy conservat Note for the shallow water equations to the total energy acts as a mathematical entropy function. """ function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquations1D) - inicenter = 0.7 + RealT = eltype(x) + inicenter = convert(RealT, 0.7) x_norm = x[1] - inicenter r = abs(x_norm) # Calculate primitive variables - H = r > 0.5 ? 3.25 : 4.0 - v = r > 0.5 ? 0.0 : 0.1882 + H = r > 0.5f0 ? 3.25f0 : 4.0f0 + v = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) b = sin(x[1]) # arbitrary continuous function return prim2cons(SVector(H, v, b), equations) @@ -185,12 +186,12 @@ end h, h_v, _ = u v = velocity(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 f1 = h_v f2 = h_v * v + p - return SVector(f1, f2, zero(eltype(u))) + return SVector(f1, f2, 0) end """ @@ -212,10 +213,8 @@ Further details are available in the paper: h_ll = waterheight(u_ll, equations) b_rr = u_rr[3] - z = zero(eltype(u_ll)) - # Bottom gradient nonconservative term: (0, g h b_x, 0) - f = SVector(z, equations.gravity * h_ll * b_rr, z) + f = SVector(0, equations.gravity * h_ll * b_rr, 0) return f end @@ -249,19 +248,17 @@ and for curvilinear 2D case in the paper: h_ll, _, b_ll = u_ll h_rr, _, b_rr = u_rr - h_average = 0.5 * (h_ll + h_rr) + h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll # Includes two parts: # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid # cross-averaging across a discontinuous bottom topography # (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry - z = zero(eltype(u_ll)) - - f = SVector(z, + f = SVector(0, equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump, - z) + 0) return f end @@ -296,15 +293,14 @@ Further details on the hydrostatic reconstruction and its motivation can be foun # Copy the reconstructed water height for easier to read code h_ll_star = u_ll_star[1] - z = zero(eltype(u_ll)) # Includes two parts: # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid # cross-averaging across a discontinuous bottom topography # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry - return SVector(z, + return SVector(0, equations.gravity * h_ll * b_ll + equations.gravity * (h_ll^2 - h_ll_star^2), - z) + 0) end """ @@ -337,10 +333,8 @@ For further details see: # Calculate jump b_jump = b_rr - b_ll - z = zero(eltype(u_ll)) - # Bottom gradient nonconservative term: (0, g h b_x, 0) - f = SVector(z, equations.gravity * h_ll * b_jump, z) + f = SVector(0, equations.gravity * h_ll * b_jump, 0) return f end @@ -367,15 +361,15 @@ Details are available in Eq. (4.1) in the paper: v_rr = velocity(u_rr, equations) # Average each factor of products in flux - h_avg = 0.5 * (h_ll + h_rr) - v_avg = 0.5 * (v_ll + v_rr) - p_avg = 0.25 * equations.gravity * (h_ll^2 + h_rr^2) + h_avg = 0.5f0 * (h_ll + h_rr) + v_avg = 0.5f0 * (v_ll + v_rr) + p_avg = 0.25f0 * equations.gravity * (h_ll^2 + h_rr^2) # Calculate fluxes depending on orientation f1 = h_avg * v_avg f2 = f1 * v_avg + p_avg - return SVector(f1, f2, zero(eltype(u_ll))) + return SVector(f1, f2, 0) end """ @@ -403,14 +397,14 @@ Further details are available in Theorem 1 of the paper: v_rr = velocity(u_rr, equations) # Average each factor of products in flux - v_avg = 0.5 * (v_ll + v_rr) - p_avg = 0.5 * equations.gravity * h_ll * h_rr + v_avg = 0.5f0 * (v_ll + v_rr) + p_avg = 0.5f0 * equations.gravity * h_ll * h_rr # Calculate fluxes depending on orientation - f1 = 0.5 * (h_v_ll + h_v_rr) + f1 = 0.5f0 * (h_v_ll + h_v_rr) f2 = f1 * v_avg + p_avg - return SVector(f1, f2, zero(eltype(u_ll))) + return SVector(f1, f2, 0) end """ @@ -438,8 +432,8 @@ Further details on this hydrostatic reconstruction and its motivation can be fou v1_rr = velocity(u_rr, equations) # Compute the reconstructed water heights - h_ll_star = max(zero(h_ll), h_ll + b_ll - max(b_ll, b_rr)) - h_rr_star = max(zero(h_rr), h_rr + b_rr - max(b_ll, b_rr)) + h_ll_star = max(0, h_ll + b_ll - max(b_ll, b_rr)) + h_rr_star = max(0, h_rr + b_rr - max(b_ll, b_rr)) # Create the conservative variables using the reconstruted water heights u_ll_star = SVector(h_ll_star, h_ll_star * v1_ll, b_ll) @@ -471,8 +465,8 @@ end equations::ShallowWaterEquations1D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) - return SVector(diss[1], diss[2], zero(eltype(u_ll))) + diss = -0.5f0 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], 0) end # Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography @@ -494,7 +488,7 @@ end factor_diss = λ_min * λ_max * inv_λ_max_minus_λ_min diss = u_rr - u_ll return factor_ll * f_ll - factor_rr * f_rr + - factor_diss * SVector(diss[1], diss[2], zero(eltype(u_ll))) + factor_diss * SVector(diss[1], diss[2], 0) end end @@ -581,7 +575,7 @@ end v = velocity(u, equations) - w1 = equations.gravity * (h + b) - 0.5 * v^2 + w1 = equations.gravity * (h + b) - 0.5f0 * v^2 w2 = v return SVector(w1, w2, b) @@ -591,7 +585,7 @@ end @inline function entropy2cons(w, equations::ShallowWaterEquations1D) w1, w2, b = w - h = (w1 + 0.5 * w2^2) / equations.gravity - b + h = (w1 + 0.5f0 * w2^2) / equations.gravity - b h_v = h * w2 return SVector(h, h_v, b) end @@ -612,7 +606,7 @@ end @inline function pressure(u, equations::ShallowWaterEquations1D) h = waterheight(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 return p end @@ -638,7 +632,7 @@ Or equation (9.17) in [this lecture notes](https://metaphor.ethz.ch/x/2019/hs/40 h_rr = waterheight(u_rr, equations) v_rr = velocity(u_rr, equations) - h_roe = 0.5 * (h_ll + h_rr) + h_roe = 0.5f0 * (h_ll + h_rr) c_roe = sqrt(equations.gravity * h_roe) h_ll_sqrt = sqrt(h_ll) @@ -658,7 +652,7 @@ end @inline function energy_total(cons, equations::ShallowWaterEquations1D) h, h_v, b = cons - e = (h_v^2) / (2 * h) + 0.5 * equations.gravity * h^2 + equations.gravity * h * b + e = (h_v^2) / (2 * h) + 0.5f0 * equations.gravity * h^2 + equations.gravity * h * b return e end diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 73fae89a0f..db8f00fc15 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -77,16 +77,18 @@ A smooth initial condition used for convergence tests in combination with """ function initial_condition_convergence_test(x, t, equations::ShallowWaterEquations2D) # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2 - c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi + RealT = eltype(x) + c = 7 + omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2)) + omega_t = 2 * convert(RealT, pi) x1, x2 = x H = c + cos(omega_x * x1) * sin(omega_x * x2) * cos(omega_t * t) - v1 = 0.5 - v2 = 1.5 - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2) + v1 = 0.5f0 + v2 = 1.5f0 + b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x1) + + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x2) return prim2cons(SVector(H, v1, v2, b), equations) end @@ -98,19 +100,20 @@ Source terms used for convergence tests in combination with (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). This manufactured solution source term is specifically designed for the bottom topography function -`b(x,y) = 2 + 0.5 * sin(sqrt(2)*pi*x) + 0.5 * sin(sqrt(2)*pi*y)` +`b(x,y) = 2 + 0.5 * sinpi(sqrt(2) * x) + 0.5 * sinpi(sqrt(2) * y)` as defined in [`initial_condition_convergence_test`](@ref). """ @inline function source_terms_convergence_test(u, x, t, equations::ShallowWaterEquations2D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocities are taken to be constants - c = 7.0 - omega_x = 2.0 * pi * sqrt(2.0) - omega_t = 2.0 * pi - omega_b = sqrt(2.0) * pi - v1 = 0.5 - v2 = 1.5 + RealT = eltype(u) + c = 7 + omega_x = 2 * convert(RealT, pi) * sqrt(convert(RealT, 2)) + omega_t = 2 * convert(RealT, pi) + omega_b = sqrt(convert(RealT, 2)) * convert(RealT, pi) + v1 = 0.5f0 + v2 = 1.5f0 x1, x2 = x @@ -126,15 +129,16 @@ as defined in [`initial_condition_convergence_test`](@ref). H_t = -omega_t * cosX * sinY * sinT # bottom topography and its gradient - b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2) - tmp1 = 0.5 * omega_b + b = 2 + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x1) + + 0.5f0 * sinpi(sqrt(convert(RealT, 2)) * x2) + tmp1 = 0.5f0 * omega_b b_x = tmp1 * cos(omega_b * x1) b_y = tmp1 * cos(omega_b * x2) du1 = H_t + v1 * (H_x - b_x) + v2 * (H_y - b_y) du2 = v1 * du1 + equations.gravity * (H - b) * H_x du3 = v2 * du1 + equations.gravity * (H - b) * H_y - return SVector(du1, du2, du3, 0.0) + return SVector(du1, du2, du3, 0) end """ @@ -145,7 +149,8 @@ Note for the shallow water equations to the total energy acts as a mathematical """ function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquations2D) # Set up polar coordinates - inicenter = SVector(0.7, 0.7) + RealT = eltype(x) + inicenter = SVector(convert(RealT, 0.7), convert(RealT, 0.7)) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] r = sqrt(x_norm^2 + y_norm^2) @@ -153,10 +158,10 @@ function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquation sin_phi, cos_phi = sincos(phi) # Calculate primitive variables - H = r > 0.5 ? 3.25 : 4.0 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - b = 0.0 # by default assume there is no bottom topography + H = r > 0.5f0 ? 3.25f0 : 4.0f0 + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi + b = 0 # by default assume there is no bottom topography return prim2cons(SVector(H, v1, v2, b), equations) end @@ -185,8 +190,8 @@ For details see Section 9.2.5 of the book: # create the "external" boundary solution state u_boundary = SVector(u_inner[1], - u_inner[2] - 2.0 * u_normal * normal[1], - u_inner[3] - 2.0 * u_normal * normal[2], + u_inner[2] - 2 * u_normal * normal[1], + u_inner[3] - 2 * u_normal * normal[2], u_inner[4]) # calculate the boundary flux @@ -228,7 +233,7 @@ end h, h_v1, h_v2, _ = u v1, v2 = velocity(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 if orientation == 1 f1 = h_v1 f2 = h_v1 * v1 + p @@ -238,7 +243,7 @@ end f2 = h_v2 * v1 f3 = h_v2 * v2 + p end - return SVector(f1, f2, f3, zero(eltype(u))) + return SVector(f1, f2, f3, 0) end # Calculate 1D flux for a single point in the normal direction @@ -250,12 +255,12 @@ end v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] h_v_normal = h * v_normal - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 f1 = h_v_normal f2 = h_v_normal * v1 + p * normal_direction[1] f3 = h_v_normal * v2 + p * normal_direction[2] - return SVector(f1, f2, f3, zero(eltype(u))) + return SVector(f1, f2, f3, 0) end """ @@ -286,12 +291,11 @@ Further details are available in the paper: h_ll = waterheight(u_ll, equations) b_rr = u_rr[4] - z = zero(eltype(u_ll)) # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) if orientation == 1 - f = SVector(z, equations.gravity * h_ll * b_rr, z, z) + f = SVector(0, equations.gravity * h_ll * b_rr, 0, 0) else # orientation == 2 - f = SVector(z, z, equations.gravity * h_ll * b_rr, z) + f = SVector(0, 0, equations.gravity * h_ll * b_rr, 0) end return f end @@ -305,10 +309,10 @@ end b_rr = u_rr[4] # Note this routine only uses the `normal_direction_average` and the average of the # bottom topography to get a quadratic split form DG gradient on curved elements - return SVector(zero(eltype(u_ll)), + return SVector(0, normal_direction_average[1] * equations.gravity * h_ll * b_rr, normal_direction_average[2] * equations.gravity * h_ll * b_rr, - zero(eltype(u_ll))) + 0) end """ @@ -349,24 +353,23 @@ and for curvilinear 2D case in the paper: h_ll, _, _, b_ll = u_ll h_rr, _, _, b_rr = u_rr - h_average = 0.5 * (h_ll + h_rr) + h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll # Includes two parts: # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid # cross-averaging across a discontinuous bottom topography # (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry - z = zero(eltype(u_ll)) if orientation == 1 - f = SVector(z, + f = SVector(0, equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump, - z, z) + 0, 0) else # orientation == 2 - f = SVector(z, z, + f = SVector(0, 0, equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump, - z) + 0) end return f @@ -389,14 +392,14 @@ end # (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump` # to handle discontinuous bathymetry - h_average = 0.5 * (h_ll + h_rr) + h_average = 0.5f0 * (h_ll + h_rr) b_jump = b_rr - b_ll f2 += normal_direction_ll[1] * equations.gravity * h_average * b_jump f3 += normal_direction_ll[2] * equations.gravity * h_average * b_jump # First and last equations do not have a nonconservative flux - f1 = f4 = zero(eltype(u_ll)) + f1 = f4 = 0 return SVector(f1, f2, f3, f4) end @@ -426,8 +429,8 @@ Further details for the hydrostatic reconstruction and its motivation can be fou v1_rr, v2_rr = velocity(u_rr, equations) # Compute the reconstructed water heights - h_ll_star = max(zero(h_ll), h_ll + b_ll - max(b_ll, b_rr)) - h_rr_star = max(zero(h_rr), h_rr + b_rr - max(b_ll, b_rr)) + h_ll_star = max(0, h_ll + b_ll - max(b_ll, b_rr)) + h_rr_star = max(0, h_rr + b_rr - max(b_ll, b_rr)) # Create the conservative variables using the reconstruted water heights u_ll_star = SVector(h_ll_star, h_ll_star * v1_ll, h_ll_star * v2_ll, b_ll) @@ -469,21 +472,20 @@ Further details for the hydrostatic reconstruction and its motivation can be fou # Copy the reconstructed water height for easier to read code h_ll_star = u_ll_star[1] - z = zero(eltype(u_ll)) # Includes two parts: # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid # cross-averaging across a discontinuous bottom topography # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry if orientation == 1 - f = SVector(z, + f = SVector(0, equations.gravity * h_ll * b_ll + equations.gravity * (h_ll^2 - h_ll_star^2), - z, z) + 0, 0) else # orientation == 2 - f = SVector(z, z, + f = SVector(0, 0, equations.gravity * h_ll * b_ll + equations.gravity * (h_ll^2 - h_ll_star^2), - z) + 0) end return f @@ -516,7 +518,7 @@ end f3 += normal_direction_ll[2] * equations.gravity * (h_ll^2 - h_ll_star^2) # First and last equations do not have a nonconservative flux - f1 = f4 = zero(eltype(u_ll)) + f1 = f4 = 0 return SVector(f1, f2, f3, f4) end @@ -560,12 +562,11 @@ For further details see: # Calculate jump b_jump = b_rr - b_ll - z = zero(eltype(u_ll)) # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) if orientation == 1 - f = SVector(z, equations.gravity * h_ll * b_jump, z, z) + f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0) else # orientation == 2 - f = SVector(z, z, equations.gravity * h_ll * b_jump, z) + f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0) end return f end @@ -583,10 +584,10 @@ end b_jump = b_rr - b_ll # Note this routine only uses the `normal_direction_average` and the average of the # bottom topography to get a quadratic split form DG gradient on curved elements - return SVector(zero(eltype(u_ll)), + return SVector(0, normal_direction_average[1] * equations.gravity * h_ll * b_jump, normal_direction_average[2] * equations.gravity * h_ll * b_jump, - zero(eltype(u_ll))) + 0) end """ @@ -611,10 +612,10 @@ Details are available in Eq. (4.1) in the paper: v1_rr, v2_rr = velocity(u_rr, equations) # Average each factor of products in flux - h_avg = 0.5 * (h_ll + h_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.25 * equations.gravity * (h_ll^2 + h_rr^2) + h_avg = 0.5f0 * (h_ll + h_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.25f0 * equations.gravity * (h_ll^2 + h_rr^2) # Calculate fluxes depending on orientation if orientation == 1 @@ -627,7 +628,7 @@ Details are available in Eq. (4.1) in the paper: f3 = f1 * v2_avg + p_avg end - return SVector(f1, f2, f3, zero(eltype(u_ll))) + return SVector(f1, f2, f3, 0) end @inline function flux_fjordholm_etal(u_ll, u_rr, normal_direction::AbstractVector, @@ -642,19 +643,19 @@ end v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] # Average each factor of products in flux - h_avg = 0.5 * (h_ll + h_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - h2_avg = 0.5 * (h_ll^2 + h_rr^2) - p_avg = 0.5 * equations.gravity * h2_avg - v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr) + h_avg = 0.5f0 * (h_ll + h_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + h2_avg = 0.5f0 * (h_ll^2 + h_rr^2) + p_avg = 0.5f0 * equations.gravity * h2_avg + v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) # Calculate fluxes depending on normal_direction f1 = h_avg * v_dot_n_avg f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] - return SVector(f1, f2, f3, zero(eltype(u_ll))) + return SVector(f1, f2, f3, 0) end """ @@ -682,22 +683,22 @@ Further details are available in Theorem 1 of the paper: v1_rr, v2_rr = velocity(u_rr, equations) # Average each factor of products in flux - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * equations.gravity * h_ll * h_rr + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * equations.gravity * h_ll * h_rr # Calculate fluxes depending on orientation if orientation == 1 - f1 = 0.5 * (h_v1_ll + h_v1_rr) + f1 = 0.5f0 * (h_v1_ll + h_v1_rr) f2 = f1 * v1_avg + p_avg f3 = f1 * v2_avg else - f1 = 0.5 * (h_v2_ll + h_v2_rr) + f1 = 0.5f0 * (h_v2_ll + h_v2_rr) f2 = f1 * v1_avg f3 = f1 * v2_avg + p_avg end - return SVector(f1, f2, f3, zero(eltype(u_ll))) + return SVector(f1, f2, f3, 0) end @inline function flux_wintermeyer_etal(u_ll, u_rr, normal_direction::AbstractVector, @@ -711,18 +712,18 @@ end v1_rr, v2_rr = velocity(u_rr, equations) # Average each factor of products in flux - h_v1_avg = 0.5 * (h_v1_ll + h_v1_rr) - h_v2_avg = 0.5 * (h_v2_ll + h_v2_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p_avg = 0.5 * equations.gravity * h_ll * h_rr + h_v1_avg = 0.5f0 * (h_v1_ll + h_v1_rr) + h_v2_avg = 0.5f0 * (h_v2_ll + h_v2_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + p_avg = 0.5f0 * equations.gravity * h_ll * h_rr # Calculate fluxes depending on normal_direction f1 = h_v1_avg * normal_direction[1] + h_v2_avg * normal_direction[2] f2 = f1 * v1_avg + p_avg * normal_direction[1] f3 = f1 * v2_avg + p_avg * normal_direction[2] - return SVector(f1, f2, f3, zero(eltype(u_ll))) + return SVector(f1, f2, f3, 0) end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the @@ -771,8 +772,8 @@ end equations::ShallowWaterEquations2D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) - return SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll))) + diss = -0.5f0 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], diss[3], 0) end # Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography @@ -794,7 +795,7 @@ end factor_diss = λ_min * λ_max * inv_λ_max_minus_λ_min diss = u_rr - u_ll return factor_ll * f_ll - factor_rr * f_rr + - factor_diss * SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll))) + factor_diss * SVector(diss[1], diss[2], diss[3], 0) end end @@ -958,7 +959,7 @@ end v1, v2 = velocity(u, equations) v_square = v1^2 + v2^2 - w1 = equations.gravity * (h + b) - 0.5 * v_square + w1 = equations.gravity * (h + b) - 0.5f0 * v_square w2 = v1 w3 = v2 return SVector(w1, w2, w3, b) @@ -968,7 +969,7 @@ end @inline function entropy2cons(w, equations::ShallowWaterEquations2D) w1, w2, w3, b = w - h = (w1 + 0.5 * (w2^2 + w3^2)) / equations.gravity - b + h = (w1 + 0.5f0 * (w2^2 + w3^2)) / equations.gravity - b h_v1 = h * w2 h_v2 = h * w3 return SVector(h, h_v1, h_v2, b) @@ -990,7 +991,7 @@ end @inline function pressure(u, equations::ShallowWaterEquations2D) h = waterheight(u, equations) - p = 0.5 * equations.gravity * h^2 + p = 0.5f0 * equations.gravity * h^2 return p end @@ -1017,7 +1018,7 @@ slides 8 and 9. h_rr = waterheight(u_rr, equations) v1_rr, v2_rr = velocity(u_rr, equations) - h_roe = 0.5 * (h_ll + h_rr) + h_roe = 0.5f0 * (h_ll + h_rr) c_roe = sqrt(equations.gravity * h_roe) h_ll_sqrt = sqrt(h_ll) @@ -1041,7 +1042,7 @@ end norm_ = norm(normal_direction) - h_roe = 0.5 * (h_ll + h_rr) + h_roe = 0.5f0 * (h_ll + h_rr) c_roe = sqrt(equations.gravity * h_roe) * norm_ h_ll_sqrt = sqrt(h_ll) @@ -1064,7 +1065,7 @@ end @inline function energy_total(cons, equations::ShallowWaterEquations2D) h, h_v1, h_v2, b = cons - e = (h_v1^2 + h_v2^2) / (2 * h) + 0.5 * equations.gravity * h^2 + + e = (h_v1^2 + h_v2^2) / (2 * h) + 0.5f0 * equations.gravity * h^2 + equations.gravity * h * b return e end diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index 0862002125..eb90fb53f3 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/src/equations/shallow_water_quasi_1d.jl @@ -72,11 +72,12 @@ function initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D) # generates a manufactured solution. # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] - Omega = sqrt(2) * pi - H = 2.0 + 0.5 * sin(Omega * x[1] - t) - v = 0.25 - b = 0.2 - 0.05 * sin(Omega * x[1]) - a = 1 + 0.1 * cos(Omega * x[1]) + RealT = eltype(x) + Omega = sqrt(convert(RealT, 2)) * convert(RealT, pi) + H = 2 + 0.5f0 * sin(Omega * x[1] - t) + v = 0.25f0 + b = convert(RealT, 0.2) - convert(RealT, 0.05) * sin(Omega * x[1]) + a = 1 + convert(RealT, 0.1) * cos(Omega * x[1]) return prim2cons(SVector(H, v, b, a), equations) end @@ -88,30 +89,31 @@ Source terms used for convergence tests in combination with (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). This manufactured solution source term is specifically designed for the bottom topography function -`b(x) = 0.2 - 0.05 * sin(sqrt(2) * pi *x[1])` and channel width 'a(x)= 1 + 0.1 * cos(sqrt(2) * pi * x[1])' +`b(x) = 0.2 - 0.05 * sinpi(sqrt(2) * x[1])` and channel width 'a(x)= 1 + 0.1 * cospi(sqrt(2) * x[1])' as defined in [`initial_condition_convergence_test`](@ref). """ @inline function source_terms_convergence_test(u, x, t, equations::ShallowWaterEquationsQuasi1D) # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because # this manufactured solution velocity is taken to be constant - Omega = sqrt(2) * pi - H = 2.0 + 0.5 * sin(Omega * x[1] - t) - H_x = 0.5 * cos(Omega * x[1] - t) * Omega - H_t = -0.5 * cos(Omega * x[1] - t) + RealT = eltype(u) + Omega = sqrt(convert(RealT, 2)) * convert(RealT, pi) + H = 2 + 0.5f0 * sin(Omega * x[1] - t) + H_x = 0.5f0 * cos(Omega * x[1] - t) * Omega + H_t = -0.5f0 * cos(Omega * x[1] - t) - v = 0.25 + v = 0.25f0 - b = 0.2 - 0.05 * sin(Omega * x[1]) - b_x = -0.05 * cos(Omega * x[1]) * Omega + b = convert(RealT, 0.2) - convert(RealT, 0.05) * sin(Omega * x[1]) + b_x = -convert(RealT, 0.05) * cos(Omega * x[1]) * Omega - a = 1 + 0.1 * cos(Omega * x[1]) - a_x = -0.1 * sin(Omega * x[1]) * Omega + a = 1 + convert(RealT, 0.1) * cos(Omega * x[1]) + a_x = -convert(RealT, 0.1) * sin(Omega * x[1]) * Omega du1 = a * H_t + v * (a_x * (H - b) + a * (H_x - b_x)) du2 = v * du1 + a * (equations.gravity * (H - b) * H_x) - return SVector(du1, du2, 0.0, 0.0) + return SVector(du1, du2, 0, 0) end # Calculate 1D conservative flux for a single point @@ -123,7 +125,7 @@ end f1 = a_h_v f2 = a_h_v * v - return SVector(f1, f2, zero(eltype(u)), zero(eltype(u))) + return SVector(f1, f2, 0, 0) end """ @@ -153,9 +155,7 @@ Further details are available in the paper: h_ll = waterheight(u_ll, equations) h_rr = waterheight(u_rr, equations) - z = zero(eltype(u_ll)) - - return SVector(z, equations.gravity * a_ll * h_ll * (h_rr + b_rr), z, z) + return SVector(0, equations.gravity * a_ll * h_ll * (h_rr + b_rr), 0, 0) end # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that @@ -200,10 +200,10 @@ Further details are available in the paper: v_ll = velocity(u_ll, equations) v_rr = velocity(u_rr, equations) - f1 = 0.5 * (a_h_v_ll + a_h_v_rr) - f2 = f1 * 0.5 * (v_ll + v_rr) + f1 = 0.5f0 * (a_h_v_ll + a_h_v_rr) + f2 = f1 * 0.5f0 * (v_ll + v_rr) - return SVector(f1, f2, zero(eltype(u_ll)), zero(eltype(u_ll))) + return SVector(f1, f2, 0, 0) end # While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that @@ -240,8 +240,8 @@ end equations::ShallowWaterEquationsQuasi1D) λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations) - diss = -0.5 * λ * (u_rr - u_ll) - return SVector(diss[1], diss[2], zero(eltype(u_ll)), zero(eltype(u_ll))) + diss = -0.5f0 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], 0, 0) end @inline function max_abs_speeds(u, equations::ShallowWaterEquationsQuasi1D) @@ -278,7 +278,7 @@ end h = waterheight(u, equations) v = velocity(u, equations) #entropy variables are the same as ones in standard shallow water equations - w1 = equations.gravity * (h + b) - 0.5 * v^2 + w1 = equations.gravity * (h + b) - 0.5f0 * v^2 w2 = v return SVector(w1, w2, b, a) @@ -306,7 +306,7 @@ end # Calculate total energy for a conservative state `cons` @inline function energy_total(cons, equations::ShallowWaterEquationsQuasi1D) a_h, a_h_v, b, a = cons - e = (a_h_v^2) / (2 * a * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + + e = (a_h_v^2) / (2 * a * a_h) + 0.5f0 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b return e end diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 4a0747d1c0..225b85a059 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -294,7 +294,7 @@ end # TODO: Taal dimension agnostic function calc_volume_integral!(du, u, - mesh::TreeMesh{1}, + mesh::Union{TreeMesh{1}, StructuredMesh{1}}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 6a26b5e95b..4f2c72cd67 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -376,7 +376,9 @@ end # TODO: Taal dimension agnostic function calc_volume_integral!(du, u, - mesh::TreeMesh{2}, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 02ff338e91..2e7e882f5e 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -424,7 +424,8 @@ end # TODO: Taal dimension agnostic function calc_volume_integral!(du, u, - mesh::TreeMesh{3}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DGSEM, cache) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index f4924bdcf7..ef5ac2c9de 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -655,6 +655,35 @@ end @test isapprox(lift, 0.029076443678087403, atol = 1e-13) @test isapprox(drag, 0.13564720009197903, atol = 1e-13) end + +@trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "tree_2d_dgsem"), + "elixir_euler_blast_wave_pure_fv.jl"), + l2=[ + 0.39957047631960346, + 0.21006912294983154, + 0.21006903549932, + 0.6280328163981136, + ], + linf=[ + 2.20417889887697, + 1.5487238480003327, + 1.5486788679247812, + 2.4656795949035857, + ], + tspan=(0.0, 0.5), + mesh=P4estMesh((64, 64), polydeg = 3, + coordinates_min = (-2.0, -2.0), + coordinates_max = (2.0, 2.0))) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index f97696d089..78230e5cf0 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -168,6 +168,30 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_convergence_pure_fv.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "tree_1d_dgsem"), + "elixir_euler_convergence_pure_fv.jl"), + mesh=StructuredMesh(16, (0.0,), (2.0,)), + l2=[ + 0.019355699748523896, + 0.022326984561234497, + 0.02523665947241734, + ], + linf=[ + 0.02895961127645519, + 0.03293442484199227, + 0.04246098278632804, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl index 81d2a7cdd8..5c452444be 100644 --- a/test/test_t8code_3d.jl +++ b/test/test_t8code_3d.jl @@ -282,6 +282,40 @@ mkdir(outdir) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "tree_3d_dgsem"), + "elixir_euler_convergence_pure_fv.jl"), + l2=[ + 0.037182410351406, + 0.032062252638283974, + 0.032062252638283974, + 0.03206225263828395, + 0.12228177813586687, + ], + linf=[ + 0.0693648413632646, + 0.0622101894740843, + 0.06221018947408474, + 0.062210189474084965, + 0.24196451799555962, + ], + mesh=T8codeMesh((4, 4, 4), polydeg = 3, + coordinates_min = (0.0, 0.0, 0.0), + coordinates_max = (2.0, 2.0, 2.0)), + # Remove SaveSolution callback + callbacks=CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_type.jl b/test/test_type.jl index c45a750c07..34fef14fee 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1517,6 +1517,99 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Lattice Boltzmann 2D" begin + for RealT in (Float32, Float64) + equations = @inferred LatticeBoltzmannEquations2D(Ma = RealT(0.1), Re = 1000) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + p = rho = dt = one(RealT) + + surface_flux_function = flux_godunov + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + + for orientation in orientations + for direction in directions + @test eltype(@inferred boundary_condition_noslip_wall(u_inner, + orientation, + direction, x, t, + surface_flux_function, + equations)) == + RealT + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) == + RealT + + @test typeof(@inferred velocity(u, orientation, equations)) == RealT + end + + @test typeof(@inferred density(p, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test eltype(@inferred velocity(u, equations)) == RealT + @test typeof(@inferred pressure(rho, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + + @test eltype(@inferred Trixi.collision_bgk(u, dt, equations)) == RealT + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + end + end + + @timed_testset "Lattice Boltzmann 3D" begin + for RealT in (Float32, Float64) + equations = @inferred LatticeBoltzmannEquations3D(Ma = RealT(0.1), Re = 1000) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT), one(RealT)) + orientations = [1, 2, 3] + p = rho = dt = one(RealT) + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) == + RealT + + @test typeof(@inferred velocity(u, orientation, equations)) == RealT + end + + @test typeof(@inferred density(p, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test eltype(@inferred velocity(u, equations)) == RealT + @test typeof(@inferred pressure(rho, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + + @test eltype(@inferred Trixi.collision_bgk(u, dt, equations)) == RealT + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred energy_kinetic(u, equations)) == RealT + @test typeof(@inferred Trixi.energy_kinetic_nondimensional(u, equations)) == + RealT + end + end + @timed_testset "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) @@ -1759,6 +1852,348 @@ isdir(outdir) && rm(outdir, recursive = true) @test typeof(@inferred energy_total(u, equations)) == RealT end end + + @timed_testset "Maxwell 1D" begin + for RealT in (Float32, Float64) + c = RealT(299_792_458) + equations = @inferred MaxwellEquations1D(c) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = SVector(one(RealT), one(RealT)) + orientation = 1 + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == + RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + end + end + + @timed_testset "Shallow Water 1D" begin + for RealT in (Float32, Float64) + equations = @inferred ShallowWaterEquations1D(gravity_constant = RealT(9.81)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT)) + orientation = 1 + directions = [1, 2] + normal_direction = SVector(one(RealT)) + + surface_flux_function = flux_lax_friedrichs + dissipation = DissipationLocalLaxFriedrichs() + numflux = FluxHLL() + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == RealT + end + + @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, direction, + equations)) == + RealT + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr, + orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr, + orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation, + equations)) == RealT + + @test eltype(eltype(@inferred hydrostatic_reconstruction_audusse_etal(u_ll, + u_rr, + equations))) == + RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred numflux(u_ll, u_rr, orientation, equations)) == RealT + # no matching method + # @test eltype(@inferred numflux(u_ll, u_rr, normal_direction, equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + + @test typeof(@inferred velocity(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test eltype(@inferred entropy2cons(u, equations)) == RealT + @test typeof(@inferred Trixi.waterheight(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred waterheight_pressure(u, equations)) == RealT + + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(u, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred lake_at_rest_error(u, equations)) == RealT + end + end + + @timed_testset "Shallow Water 2D" begin + for RealT in (Float32, Float64) + equations = @inferred ShallowWaterEquations2D(gravity_constant = RealT(9.81)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), + zero(RealT)) + + surface_flux_function = flux_lax_friedrichs + dissipation = DissipationLocalLaxFriedrichs() + numflux = FluxHLL() + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + normal_direction, + x, t, + surface_flux_function, + equations)) == RealT + end + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == RealT + @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, normal_direction, + equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred numflux(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_slip_wall(u_inner, + orientation, + direction, x, t, + surface_flux_function, + equations)) == + RealT + end + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation, + equations))) == + RealT + @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation, + equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred numflux(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred Trixi.calc_wavespeed_roe(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(eltype(@inferred hydrostatic_reconstruction_audusse_etal(u_ll, + u_rr, + equations))) == + RealT + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred velocity(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test eltype(@inferred entropy2cons(u, equations)) == RealT + @test typeof(@inferred Trixi.waterheight(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred waterheight_pressure(u, equations)) == RealT + + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(u, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred lake_at_rest_error(u, equations)) == RealT + end + end + + @timed_testset "Shallow Water Quasi 1D" begin + for RealT in (Float32, Float64) + equations = @inferred ShallowWaterEquationsQuasi1D(gravity_constant = RealT(9.81)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) + orientation = 1 + normal_direction = normal_ll = normal_rr = SVector(one(RealT)) + + dissipation = DissipationLocalLaxFriedrichs() + + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, + normal_direction, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, + normal_rr, + equations)) == RealT + @test eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_chan_etal(u_ll, u_rr, normal_direction, + equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred dissipation(u_ll, u_rr, orientation, equations)) == RealT + @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) == + RealT + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test typeof(@inferred velocity(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred Trixi.waterheight(u, equations)) == RealT + + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred lake_at_rest_error(u, equations)) == RealT + end + end end end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index 09cf506576..aa66ce2556 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -439,14 +439,17 @@ end # Neither Mach number nor velocity set @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = nothing, Re = 1000) # Both Mach number and velocity set - @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, u0 = 1) + @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, + u0 = 1.0) # Neither Reynolds number nor viscosity set @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = nothing) # Both Reynolds number and viscosity set - @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, nu = 1) + @test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, + nu = 1.0) # No non-dimensional values set - @test LatticeBoltzmannEquations2D(Ma = nothing, Re = nothing, u0 = 1, nu = 1) isa + @test LatticeBoltzmannEquations2D(Ma = nothing, Re = nothing, u0 = 1.0, + nu = 1.0) isa LatticeBoltzmannEquations2D end @@ -454,14 +457,17 @@ end # Neither Mach number nor velocity set @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = nothing, Re = 1000) # Both Mach number and velocity set - @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, u0 = 1) + @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, + u0 = 1.0) # Neither Reynolds number nor viscosity set @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = nothing) # Both Reynolds number and viscosity set - @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, nu = 1) + @test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, + nu = 1.0) # No non-dimensional values set - @test LatticeBoltzmannEquations3D(Ma = nothing, Re = nothing, u0 = 1, nu = 1) isa + @test LatticeBoltzmannEquations3D(Ma = nothing, Re = nothing, u0 = 1.0, + nu = 1.0) isa LatticeBoltzmannEquations3D end