Make TreeMesh type general

examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl

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

examples/tree_1d_dgsem/elixir_advection_diffusion.jl

@@ -13,8 +13,8 @@ equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 
-coordinates_min = -pi # minimum coordinate
-coordinates_max = pi # maximum coordinate
+coordinates_min = -convert(Float64, pi) # minimum coordinate
+coordinates_max = convert(Float64, pi) # maximum coordinate
 
 # Create a uniformly refined mesh with periodic boundaries
 mesh = TreeMesh(coordinates_min, coordinates_max,

examples/tree_2d_dgsem/elixir_mhd_ec_float32.jl

@@ -0,0 +1,62 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations2D(1.4f0)
+
+initial_condition = initial_condition_weak_blast_wave
+
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+solver = DGSEM(polydeg = 3, RealT = Float32,
+               surface_flux = (flux_hindenlang_gassner, flux_nonconservative_powell),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2.0f0, -2.0f0)
+coordinates_max = (2.0f0, 2.0f0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                RealT = Float32)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0f0, 0.4f0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+cfl = 1.0f0
+stepsize_callback = StepsizeCallback(cfl = cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5f0, cfl = cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0f0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

src/auxiliary/precompile.jl

@@ -226,8 +226,8 @@ function _precompile_manual_()
     for RealT in (Int, Float64)
         @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),
                                       NamedTuple{(:initial_refinement_level, :n_cells_max),
-                                                 Tuple{Int, Int}}, Type{TreeMesh}, RealT,
-                                      RealT})
+                                                 Tuple{Int, Int}}, Type{TreeMesh},
+                                      RealT, RealT})
         @assert Base.precompile(Tuple{Core.kwftype(typeof(Trixi.Type)),
                                       NamedTuple{(:initial_refinement_level, :n_cells_max),
                                                  Tuple{Int, Int}}, Type{TreeMesh},
@@ -244,10 +244,12 @@ function _precompile_manual_()
     end
     for TreeType in (SerialTree, ParallelTree), NDIMS in 1:3
         @assert Base.precompile(Tuple{typeof(Trixi.initialize!),
-                                      TreeMesh{NDIMS, TreeType{NDIMS}}, Int, Tuple{},
+                                      TreeMesh{NDIMS, TreeType{NDIMS}, Float64}, Int,
+                                      Tuple{},
                                       Tuple{}})
         @assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
-                                      TreeMesh{NDIMS, TreeType{NDIMS}}, String, Int})
+                                      TreeMesh{NDIMS, TreeType{NDIMS}, Float64}, String,
+                                      Int})
     end
 
     # Constructors: linear advection
@@ -396,58 +398,58 @@ function _precompile_manual_()
 
         # 1D, serial
         @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
-                                      TreeMesh{1, Trixi.SerialTree{1}},
+                                      TreeMesh{1, Trixi.SerialTree{1}, RealT},
                                       Trixi.ElementContainer1D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
-                                      TreeMesh{1, Trixi.SerialTree{1}},
+                                      TreeMesh{1, Trixi.SerialTree{1}, RealT},
                                       Trixi.ElementContainer1D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
-                                      TreeMesh{1, Trixi.SerialTree{1}}, String})
+                                      TreeMesh{1, Trixi.SerialTree{1}, RealT}, String})
 
         # 2D, serial
         @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
-                                      TreeMesh{2, Trixi.SerialTree{2}},
+                                      TreeMesh{2, Trixi.SerialTree{2}, RealT},
                                       Trixi.ElementContainer2D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
-                                      TreeMesh{2, Trixi.SerialTree{2}},
+                                      TreeMesh{2, Trixi.SerialTree{2}, RealT},
                                       Trixi.ElementContainer2D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.init_mortars), Array{Int, 1},
-                                      TreeMesh{2, Trixi.SerialTree{2}},
+                                      TreeMesh{2, Trixi.SerialTree{2}, RealT},
                                       Trixi.ElementContainer2D{RealT, uEltype},
                                       mortar_type})
         @assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
-                                      TreeMesh{2, Trixi.SerialTree{2}}, String})
+                                      TreeMesh{2, Trixi.SerialTree{2}, RealT}, String})
 
         # 2D, parallel
         @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
-                                      TreeMesh{2, Trixi.ParallelTree{2}},
+                                      TreeMesh{2, Trixi.ParallelTree{2}, RealT},
                                       Trixi.ElementContainer2D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
-                                      TreeMesh{2, Trixi.ParallelTree{2}},
+                                      TreeMesh{2, Trixi.ParallelTree{2}, RealT},
                                       Trixi.ElementContainer2D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.init_mortars), Array{Int, 1},
-                                      TreeMesh{2, Trixi.ParallelTree{2}},
+                                      TreeMesh{2, Trixi.ParallelTree{2}, RealT},
                                       Trixi.ElementContainer2D{RealT, uEltype},
                                       mortar_type})
         @assert Base.precompile(Tuple{typeof(Trixi.init_mpi_interfaces), Array{Int, 1},
-                                      TreeMesh{2, Trixi.ParallelTree{2}},
+                                      TreeMesh{2, Trixi.ParallelTree{2}, RealT},
                                       Trixi.ElementContainer2D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
-                                      TreeMesh{2, Trixi.ParallelTree{2}}, String})
+                                      TreeMesh{2, Trixi.ParallelTree{2}, RealT}, String})
 
         # 3D, serial
         @assert Base.precompile(Tuple{typeof(Trixi.init_boundaries), Array{Int, 1},
-                                      TreeMesh{3, Trixi.SerialTree{3}},
+                                      TreeMesh{3, Trixi.SerialTree{3}, RealT},
                                       Trixi.ElementContainer3D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.init_interfaces), Array{Int, 1},
-                                      TreeMesh{3, Trixi.SerialTree{3}},
+                                      TreeMesh{3, Trixi.SerialTree{3}, RealT},
                                       Trixi.ElementContainer3D{RealT, uEltype}})
         @assert Base.precompile(Tuple{typeof(Trixi.init_mortars), Array{Int, 1},
-                                      TreeMesh{3, Trixi.SerialTree{3}},
+                                      TreeMesh{3, Trixi.SerialTree{3}, RealT},
                                       Trixi.ElementContainer3D{RealT, uEltype},
                                       mortar_type})
         @assert Base.precompile(Tuple{typeof(Trixi.save_mesh_file),
-                                      TreeMesh{3, Trixi.SerialTree{3}}, String})
+                                      TreeMesh{3, Trixi.SerialTree{3}, RealT}, String})
     end
 
     # various `show` methods
@@ -456,16 +458,16 @@ function _precompile_manual_()
         for NDIMS in 1:3
             # serial
             @assert Base.precompile(Tuple{typeof(show), Base.TTY,
-                                          TreeMesh{NDIMS, Trixi.SerialTree{NDIMS}}})
+                                          TreeMesh{NDIMS, Trixi.SerialTree{NDIMS}, RealT}})
             @assert Base.precompile(Tuple{typeof(show), IOContext{Base.TTY},
                                           MIME"text/plain",
-                                          TreeMesh{NDIMS, Trixi.SerialTree{NDIMS}}})
+                                          TreeMesh{NDIMS, Trixi.SerialTree{NDIMS}, RealT}})
             # parallel
             @assert Base.precompile(Tuple{typeof(show), Base.TTY,
-                                          TreeMesh{NDIMS, Trixi.ParallelTree{NDIMS}}})
+                                          TreeMesh{NDIMS, Trixi.ParallelTree{NDIMS}, RealT}})
             @assert Base.precompile(Tuple{typeof(show), IOContext{Base.TTY},
                                           MIME"text/plain",
-                                          TreeMesh{NDIMS, Trixi.ParallelTree{NDIMS}}})
+                                          TreeMesh{NDIMS, Trixi.ParallelTree{NDIMS}, RealT}})
         end
 
         # equations

src/meshes/abstract_tree.jl

@@ -388,7 +388,8 @@ function refine!(t::AbstractTree, cell_ids,
 end
 
 # Refine all leaf cells with coordinates in a given rectangular box
-function refine_box!(t::AbstractTree{NDIMS}, coordinates_min,
+function refine_box!(t::AbstractTree{NDIMS},
+                     coordinates_min,
                      coordinates_max) where {NDIMS}
     for dim in 1:NDIMS
         @assert coordinates_min[dim]<coordinates_max[dim] "Minimum coordinates are not minimum."
@@ -404,10 +405,11 @@ function refine_box!(t::AbstractTree{NDIMS}, coordinates_min,
     refine!(t, cells)
 end
 
-# Convenience method for 1D
-function refine_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real)
-    return refine_box!(t, [convert(Float64, coordinates_min)],
-                       [convert(Float64, coordinates_max)])
+# Convenience method for 1D (arguments are no arrays)
+function refine_box!(t::AbstractTree{1},
+                     coordinates_min::Real,
+                     coordinates_max::Real)
+    return refine_box!(t, [coordinates_min], [coordinates_max])
 end
 
 # Refine all leaf cells with coordinates in a given sphere
@@ -617,8 +619,9 @@ end
 coarsen!(t::AbstractTree, cell_id::Int) = coarsen!(t::AbstractTree, [cell_id])
 
 # Coarsen all viable parent cells with coordinates in a given rectangular box
-function coarsen_box!(t::AbstractTree{NDIMS}, coordinates_min::AbstractArray{Float64},
-                      coordinates_max::AbstractArray{Float64}) where {NDIMS}
+function coarsen_box!(t::AbstractTree{NDIMS},
+                      coordinates_min,
+                      coordinates_max) where {NDIMS}
     for dim in 1:NDIMS
         @assert coordinates_min[dim]<coordinates_max[dim] "Minimum coordinates are not minimum."
     end
@@ -642,10 +645,11 @@ function coarsen_box!(t::AbstractTree{NDIMS}, coordinates_min::AbstractArray{Flo
     coarsen!(t, parents)
 end
 
-# Convenience method for 1D
-function coarsen_box!(t::AbstractTree{1}, coordinates_min::Real, coordinates_max::Real)
-    return coarsen_box!(t, [convert(Float64, coordinates_min)],
-                        [convert(Float64, coordinates_max)])
+# Convenience method for 1D (arguments are no arrays)
+function coarsen_box!(t::AbstractTree{1},
+                      coordinates_min::Real,
+                      coordinates_max::Real)
+    return coarsen_box!(t, [coordinates_min], [coordinates_max])
 end
 
 # Return coordinates of a child cell based on its relative position to the parent.

src/meshes/mesh_io.jl

@@ -335,7 +335,8 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
         capacity = h5open(mesh_file, "r") do file
             return read(attributes(file)["capacity"])
         end
-        mesh = TreeMesh(SerialTree{ndims}, max(n_cells_max, capacity))
+        mesh = TreeMesh(SerialTree{ndims, RealT}, max(n_cells_max, capacity),
+                        RealT = RealT)
         load_mesh!(mesh, mesh_file)
     elseif mesh_type in ("StructuredMesh", "StructuredMeshView")
         size_, mapping_as_string = h5open(mesh_file, "r") do file
@@ -484,7 +485,9 @@ function load_mesh_parallel(mesh_file::AbstractString; n_cells_max, RealT)
             capacity = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
         end
 
-        mesh = TreeMesh(ParallelTree{ndims_}, max(n_cells, n_cells_max, capacity))
+        mesh = TreeMesh(ParallelTree{ndims_, RealT},
+                        max(n_cells, n_cells_max, capacity),
+                        RealT = RealT)
         load_mesh!(mesh, mesh_file)
     elseif mesh_type == "P4estMesh"
         if mpi_isroot()