diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl index f2e14273ae7..2632e80db71 100644 --- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -96,11 +96,9 @@ volume_integral = VolumeIntegralFluxDifferencing(volume_flux) solver = DGSEM(basis, surface_flux, volume_integral) -coordinates_min = (0.0, 0.0) -coordinates_max = (20_000.0, 10_000.0) +coordinates_min = (0.0, -5000.0) +coordinates_max = (20_000.0, 15_000.0) -# Same coordinates as in examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl -# However TreeMesh will generate a 20_000 x 20_000 square domain instead mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 6, n_cells_max = 10_000, diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl index 4c30406680a..874b8845cb1 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -12,8 +12,8 @@ equations_parabolic = LaplaceDiffusion3D(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 = (-1.0, -0.5, -0.25) # minimum coordinates (min(x), min(y), min(z)) -coordinates_max = (0.0, 0.5, 0.25) # maximum coordinates (max(x), max(y), max(z)) +coordinates_min = (-1.0, -0.5, -0.5) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (0.0, 0.5, 0.5) # maximum coordinates (max(x), max(y), max(z)) # Create a uniformly refined mesh with periodic boundaries mesh = TreeMesh(coordinates_min, coordinates_max, diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index f9413e0f1f9..8594b0f5fa8 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -136,10 +136,9 @@ function initial_condition_convergence_test(x, t, rho = ini # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1) - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) prim1 = rho * v1 diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index bac76c9a021..2c6bbdc4a23 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -141,10 +141,9 @@ function initial_condition_convergence_test(x, t, rho = ini # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1) - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) prim1 = rho * v1 @@ -185,10 +184,9 @@ Source terms used for convergence tests in combination with tmp6 = tmp2 + c # Here we compute an arbitrary number of different rhos. (one rho is double the next rho while the sum of all rhos is 1 - du_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - tmp1 + du_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + tmp1 / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) du1 = tmp5 diff --git a/src/equations/ideal_glm_mhd_multicomponent_1d.jl b/src/equations/ideal_glm_mhd_multicomponent_1d.jl index 86a69e9fe10..a5f16bd78d8 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_1d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_1d.jl @@ -89,11 +89,10 @@ function initial_condition_convergence_test(x, t, # domain must be set to [0, 1], γ = 5/3 RealT = eltype(x) - rho = 1 - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + rho = one(RealT) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) v1 = 0 # TODO: sincospi @@ -132,16 +131,16 @@ function initial_condition_weak_blast_wave(x, t, if r > 0.5f0 rho = one(RealT) prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * - (1 - 2) / (1 - - 2^ncomponents(equations)) * - rho + (1 - 2) * rho / + (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) else rho = convert(RealT, 1.1691) prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * - (1 - 2) / (1 - - 2^ncomponents(equations)) * - rho + (1 - 2) * rho / + (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) end v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index c56294b6643..99fe391301d 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -114,11 +114,10 @@ function initial_condition_convergence_test(x, t, alpha = 0.25f0 * convert(RealT, pi) x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) B_perp = convert(RealT, 0.1) * sinpi(2 * x_perp) - rho = 1 - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + rho = one(RealT) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) v1 = -B_perp * sin(alpha) @@ -157,13 +156,12 @@ function initial_condition_weak_blast_wave(x, t, prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ? 2^(i - 1) * (1 - 2) / + (RealT(1) - + 2^ncomponents(equations)) : + 2^(i - 1) * (1 - 2) * + RealT(1.1691) / (1 - - 2^ncomponents(equations)) * - one(RealT) : - 2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - convert(RealT, 1.1691) + 2^ncomponents(equations)) for i in eachcomponent(equations)) v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi diff --git a/src/equations/traffic_flow_lwr_1d.jl b/src/equations/traffic_flow_lwr_1d.jl index c41fbb2809e..89530f70b6b 100644 --- a/src/equations/traffic_flow_lwr_1d.jl +++ b/src/equations/traffic_flow_lwr_1d.jl @@ -17,7 +17,7 @@ the maximum possible speed (e.g. due to speed limits) is $v_{\text{max}}$. For more details see e.g. Section 11.1 of - Randall LeVeque (2002) Finite Volume Methods for Hyperbolic Problems -[DOI: 10.1017/CBO9780511791253]https://doi.org/10.1017/CBO9780511791253 +[DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) """ struct TrafficFlowLWREquations1D{RealT <: Real} <: AbstractTrafficFlowLWREquations{1, 1} v_max::RealT diff --git a/src/meshes/tree_mesh.jl b/src/meshes/tree_mesh.jl index 933bfa62722..458a43203d2 100644 --- a/src/meshes/tree_mesh.jl +++ b/src/meshes/tree_mesh.jl @@ -119,9 +119,12 @@ function TreeMesh(coordinates_min::NTuple{NDIMS, Real}, throw(ArgumentError("`initial_refinement_level` must be a non-negative integer (provided `initial_refinement_level = $initial_refinement_level`)")) end - # Domain length is calculated as the maximum length in any axis direction + # TreeMesh requires equal domain lengths in all dimensions domain_center = @. (coordinates_min + coordinates_max) / 2 - domain_length = maximum(coordinates_max .- coordinates_min) + domain_length = coordinates_max[1] - coordinates_min[1] + if !all(coordinates_max[i] - coordinates_min[i] ≈ domain_length for i in 2:NDIMS) + throw(ArgumentError("The TreeMesh domain must be a hypercube (provided `coordinates_max` .- `coordinates_min` = $(coordinates_max .- coordinates_min))")) + end # TODO: MPI, create nice interface for a parallel tree/mesh if mpi_isparallel() diff --git a/test/test_unit.jl b/test/test_unit.jl index bccdcf8faaa..4e8945eab8a 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -65,6 +65,16 @@ end @timed_testset "TreeMesh" begin @testset "constructors" begin @test TreeMesh{1, Trixi.SerialTree{1}}(1, 5.0, 2.0) isa TreeMesh + + # Invalid domain length check (TreeMesh expects a hypercube) + # 2D + @test_throws ArgumentError TreeMesh((-0.5, 0.0), (1.0, 2.0), + initial_refinement_level = 2, + n_cells_max = 10_000) + # 3D + @test_throws ArgumentError TreeMesh((-0.5, 0.0, -0.2), (1.0, 2.0, 1.5), + initial_refinement_level = 2, + n_cells_max = 10_000) end end