Merge branch 'main' into Ind3Level_Comb_Unified

trixi-framework · Oct 16, 2024 · 5ce7271 · 5ce7271
2 parents 87ab1a6 + cd6f796
commit 5ce7271
Show file tree

Hide file tree

Showing 9 changed files with 48 additions and 43 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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