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,64 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations2D(1.4f0)
+
+RealT = Float32
+
+initial_condition = initial_condition_weak_blast_wave
+
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+solver = DGSEM(polydeg = 3, RealT = RealT,
+               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 = RealT)
+
+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

@@ -257,7 +257,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
@@ -381,7 +382,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()

src/meshes/parallel_tree.jl

@@ -25,24 +25,24 @@
 # function, which is required for implementing level-wise refinement in a sane
 # way. Also, depth-first ordering *might* not be guaranteed during
 # refinement/coarsening operations.
-mutable struct ParallelTree{NDIMS} <: AbstractTree{NDIMS}
+mutable struct ParallelTree{NDIMS, RealT <: Real} <: AbstractTree{NDIMS}
     parent_ids::Vector{Int}
     child_ids::Matrix{Int}
     neighbor_ids::Matrix{Int}
     levels::Vector{Int}
-    coordinates::Matrix{Float64}
+    coordinates::Matrix{RealT}
     original_cell_ids::Vector{Int}
     mpi_ranks::Vector{Int}
 
     capacity::Int
     length::Int
     dummy::Int
 
-    center_level_0::SVector{NDIMS, Float64}
-    length_level_0::Float64
+    center_level_0::SVector{NDIMS, RealT}
+    length_level_0::RealT
     periodicity::NTuple{NDIMS, Bool}
 
-    function ParallelTree{NDIMS}(capacity::Integer) where {NDIMS}
+    function ParallelTree{NDIMS, RealT}(capacity::Integer) where {NDIMS, RealT <: Real}
         # Verify that NDIMS is an integer
         @assert NDIMS isa Integer
 
@@ -70,29 +70,37 @@ mutable struct ParallelTree{NDIMS} <: AbstractTree{NDIMS}
     end
 end
 
-# Constructor for passing the dimension as an argument
-ParallelTree(::Val{NDIMS}, args...) where {NDIMS} = ParallelTree{NDIMS}(args...)
+# Constructor for passing the dimension as an argument. Default datatype: Float64
+ParallelTree(::Val{NDIMS}, args...) where {NDIMS} = ParallelTree{NDIMS, Float64}(args...)
 
 # Create and initialize tree
-function ParallelTree{NDIMS}(capacity::Int, center::AbstractArray{Float64},
-                             length::Real, periodicity = true) where {NDIMS}
+function ParallelTree{NDIMS, RealT}(capacity::Int, center::AbstractArray{RealT},
+                                    length::RealT,
+                                    periodicity = true) where {NDIMS, RealT <: Real}
     # Create instance
-    t = ParallelTree{NDIMS}(capacity)
+    t = ParallelTree{NDIMS, RealT}(capacity)
 
     # Initialize root cell
     init!(t, center, length, periodicity)
 
     return t
 end
 
+function ParallelTree{NDIMS}(capacity::Int, center::AbstractArray{RealT},
+                             length::RealT,
+                             periodicity = true) where {NDIMS, RealT <: Real}
+    ParallelTree{NDIMS, RealT}(capacity, center, length, periodicity)
+end
+
 # Constructor accepting a single number as center (as opposed to an array) for 1D
-function ParallelTree{1}(cap::Int, center::Real, len::Real, periodicity = true)
-    ParallelTree{1}(cap, [convert(Float64, center)], len, periodicity)
+function ParallelTree{1, RealT}(cap::Int, center::RealT, len::RealT,
+                                periodicity = true) where {RealT <: Real}
+    ParallelTree{1, RealT}(cap, [center], len, periodicity)
 end
 
 # Clear tree with deleting data structures, store center and length, and create root cell
-function init!(t::ParallelTree, center::AbstractArray{Float64}, length::Real,
-               periodicity = true)
+function init!(t::ParallelTree, center::AbstractArray{RealT}, length::RealT,
+               periodicity = true) where {RealT}
     clear!(t)
 
     # Set domain information
@@ -222,12 +230,13 @@ function raw_copy!(target::ParallelTree, source::ParallelTree, first::Int, last:
 end
 
 # Reset data structures by recreating all internal storage containers and invalidating all elements
-function reset_data_structures!(t::ParallelTree{NDIMS}) where {NDIMS}
+function reset_data_structures!(t::ParallelTree{NDIMS, RealT}) where {NDIMS,
+                                                                      RealT <: Real}
     t.parent_ids = Vector{Int}(undef, t.capacity + 1)
     t.child_ids = Matrix{Int}(undef, 2^NDIMS, t.capacity + 1)
     t.neighbor_ids = Matrix{Int}(undef, 2 * NDIMS, t.capacity + 1)
     t.levels = Vector{Int}(undef, t.capacity + 1)
-    t.coordinates = Matrix{Float64}(undef, NDIMS, t.capacity + 1)
+    t.coordinates = Matrix{RealT}(undef, NDIMS, t.capacity + 1)
     t.original_cell_ids = Vector{Int}(undef, t.capacity + 1)
     t.mpi_ranks = Vector{Int}(undef, t.capacity + 1)