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(-pi, Float64) # minimum coordinate
+coordinates_max = convert(pi, Float64) # maximum coordinate
 
 # Create a uniformly refined mesh with periodic boundaries
 mesh = TreeMesh(coordinates_min, coordinates_max,

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/mesh_io.jl

@@ -257,7 +257,7 @@ 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}, RealT, max(n_cells_max, capacity))
         load_mesh!(mesh, mesh_file)
     elseif mesh_type in ("StructuredMesh", "StructuredMeshView")
         size_, mapping_as_string = h5open(mesh_file, "r") do file
@@ -381,7 +381,8 @@ 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}, RealT,
+                        max(n_cells, n_cells_max, capacity))
         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,40 @@ 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 and datatype as an argument
+ParallelTree(::Val{NDIMS}, RealT::DataType, args...) where {NDIMS} = ParallelTree{NDIMS,
+                                                                                  RealT}(args...)
+# 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 +233,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)
 

src/meshes/serial_tree.jl

@@ -25,23 +25,23 @@
 # function, which is required for implementing level-wise refinement in a sane
 # way. Also, depth-first ordering *might* not by guaranteed during
 # refinement/coarsening operations.
-mutable struct SerialTree{NDIMS} <: AbstractTree{NDIMS}
+mutable struct SerialTree{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}
 
     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 SerialTree{NDIMS}(capacity::Integer) where {NDIMS}
+    function SerialTree{NDIMS, RealT}(capacity::Integer) where {NDIMS, RealT <: Real}
         # Verify that NDIMS is an integer
         @assert NDIMS isa Integer
 
@@ -68,29 +68,40 @@ mutable struct SerialTree{NDIMS} <: AbstractTree{NDIMS}
     end
 end
 
-# Constructor for passing the dimension as an argument
-SerialTree(::Val{NDIMS}, args...) where {NDIMS} = SerialTree{NDIMS}(args...)
+# Constructor for passing the dimension and datatype as an argument
+SerialTree(::Val{NDIMS}, RealT::DataType, args...) where {NDIMS} = SerialTree{NDIMS,
+                                                                              RealT}(args...)
+# Default datatype: Float64
+SerialTree(::Val{NDIMS}, args...) where {NDIMS} = SerialTree{NDIMS, Float64}(args...)
 
 # Create and initialize tree
-function SerialTree{NDIMS}(capacity::Int, center::AbstractArray{Float64},
-                           length::Real, periodicity = true) where {NDIMS}
+function SerialTree{NDIMS, RealT}(capacity::Int, center::AbstractArray{RealT},
+                                  length::RealT,
+                                  periodicity = true) where {NDIMS, RealT <: Real}
     # Create instance
-    t = SerialTree{NDIMS}(capacity)
+    t = SerialTree{NDIMS, RealT}(capacity)
 
     # Initialize root cell
     init!(t, center, length, periodicity)
 
     return t
 end
 
+function SerialTree{NDIMS}(capacity::Int, center::AbstractArray{RealT},
+                           length::RealT,
+                           periodicity = true) where {NDIMS, RealT <: Real}
+    t = SerialTree{NDIMS, RealT}(capacity, center, length, periodicity)
+end
+
 # Constructor accepting a single number as center (as opposed to an array) for 1D
-function SerialTree{1}(cap::Int, center::Real, len::Real, periodicity = true)
-    SerialTree{1}(cap, [convert(Float64, center)], len, periodicity)
+function SerialTree{1, RealT}(cap::Int, center::RealT, len::RealT,
+                              periodicity = true) where {RealT <: Real}
+    SerialTree{1, RealT}(cap, [center], len, periodicity)
 end
 
 # Clear tree with deleting data structures, store center and length, and create root cell
-function init!(t::SerialTree, center::AbstractArray{Float64}, length::Real,
-               periodicity = true)
+function init!(t::SerialTree, center::AbstractArray{RealT}, length::RealT,
+               periodicity = true) where {RealT}
     clear!(t)
 
     # Set domain information
@@ -203,12 +214,13 @@ function raw_copy!(target::SerialTree, source::SerialTree, first::Int, last::Int
 end
 
 # Reset data structures by recreating all internal storage containers and invalidating all elements
-function reset_data_structures!(t::SerialTree{NDIMS}) where {NDIMS}
+function reset_data_structures!(t::SerialTree{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)
 
     invalidate!(t, 1, capacity(t) + 1)