From 9fa4e27f38c7cfee77bf871dbb3c42de514b751d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 7 Nov 2024 20:57:47 +0100 Subject: [PATCH] Make `TreeMesh` type general (#2129) * Make TreeMesh type general * fmt * get mesh read * work * parallel read in * tests * get things workin before commit tests * Avoid irrational type * Update examples/tree_1d_dgsem/elixir_advection_diffusion.jl * remove unused * try bounding type parameter * another try aqua test pass * Be less breaking * example * restore comments * test vals * demand same type for refinement patches * change vals * avoid breaking changes * warning if non-agreeing coordinates * make examples warning-free * avoid warning * fmt * type general NaN * typo * Consistent constructors * hard-code realT * tests for consistency constructors --------- Co-authored-by: Hendrik Ranocha --- .../elixir_eulergravity_sedov_blast_wave.jl | 4 +- .../elixir_advection_diffusion.jl | 4 +- .../tree_2d_dgsem/elixir_mhd_ec_float32.jl | 62 +++++++++++++ src/auxiliary/precompile.jl | 50 ++++++----- src/meshes/abstract_tree.jl | 26 +++--- src/meshes/mesh_io.jl | 7 +- src/meshes/parallel_tree.jl | 57 +++++++----- src/meshes/serial_tree.jl | 57 +++++++----- src/meshes/tree_mesh.jl | 88 +++++++++++-------- test/test_tree_2d_mhd.jl | 30 +++++++ test/test_unit.jl | 25 ++++-- 11 files changed, 283 insertions(+), 127 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_mhd_ec_float32.jl diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl index bc7ceb97c8b..861619d2344 100644 --- a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl +++ b/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, diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl index b75e3622ac1..e9530692d4b 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl +++ b/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, diff --git a/examples/tree_2d_dgsem/elixir_mhd_ec_float32.jl b/examples/tree_2d_dgsem/elixir_mhd_ec_float32.jl new file mode 100644 index 00000000000..0270ba2ff65 --- /dev/null +++ b/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 diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 4d5399b5ba3..dfda8ece687 100644 --- a/src/auxiliary/precompile.jl +++ b/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 diff --git a/src/meshes/abstract_tree.jl b/src/meshes/abstract_tree.jl index 469189ff50c..243991745e1 100644 --- a/src/meshes/abstract_tree.jl +++ b/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] NaN, NDIMS)) - t.length_level_0 = NaN + t.center_level_0 = SVector(ntuple(_ -> convert(RealT, NaN), NDIMS)) + t.length_level_0 = convert(RealT, NaN) return t 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) +# Constructors accepting a single number as center (as opposed to an array) for 1D +function ParallelTree{1, RealT}(cap::Int, center::RealT, len::RealT, + periodicity = true) where {RealT <: Real} + ParallelTree{1, RealT}(cap, [center], len, periodicity) +end +function ParallelTree{1}(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 @@ -187,7 +198,8 @@ end # Reset range of cells to values that are prone to cause errors as soon as they are used. # # Rationale: If an invalid cell is accidentally used, we want to know it as soon as possible. -function invalidate!(t::ParallelTree, first::Int, last::Int) +function invalidate!(t::ParallelTree{NDIMS, RealT}, + first::Int, last::Int) where {NDIMS, RealT <: Real} @assert first > 0 @assert last <= t.capacity + 1 @@ -196,7 +208,7 @@ function invalidate!(t::ParallelTree, first::Int, last::Int) t.child_ids[:, first:last] .= typemin(Int) t.neighbor_ids[:, first:last] .= typemin(Int) t.levels[first:last] .= typemin(Int) - t.coordinates[:, first:last] .= NaN + t.coordinates[:, first:last] .= convert(RealT, NaN) # `NaN` is of type Float64 t.original_cell_ids[first:last] .= typemin(Int) t.mpi_ranks[first:last] .= typemin(Int) @@ -222,12 +234,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) diff --git a/src/meshes/serial_tree.jl b/src/meshes/serial_tree.jl index 143ac19f6ee..0dea057efbb 100644 --- a/src/meshes/serial_tree.jl +++ b/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 @@ -54,43 +54,54 @@ mutable struct SerialTree{NDIMS} <: AbstractTree{NDIMS} t.child_ids = fill(typemin(Int), 2^NDIMS, capacity + 1) t.neighbor_ids = fill(typemin(Int), 2 * NDIMS, capacity + 1) t.levels = fill(typemin(Int), capacity + 1) - t.coordinates = fill(NaN, NDIMS, capacity + 1) + t.coordinates = fill(convert(RealT, NaN), NDIMS, capacity + 1) # `NaN` is of type Float64 t.original_cell_ids = fill(typemin(Int), capacity + 1) t.capacity = capacity t.length = 0 t.dummy = capacity + 1 - t.center_level_0 = SVector(ntuple(_ -> NaN, NDIMS)) - t.length_level_0 = NaN + t.center_level_0 = SVector(ntuple(_ -> convert(RealT, NaN), NDIMS)) + t.length_level_0 = convert(RealT, NaN) return t end end -# Constructor for passing the dimension as an argument -SerialTree(::Val{NDIMS}, args...) where {NDIMS} = SerialTree{NDIMS}(args...) +# Constructor for passing the dimension as an argument. 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) +# Constructors accepting a single number as center (as opposed to an array) for 1D +function SerialTree{1, RealT}(cap::Int, center::RealT, len::RealT, + periodicity = true) where {RealT <: Real} + SerialTree{1, RealT}(cap, [center], len, periodicity) +end +function SerialTree{1}(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 @@ -170,7 +181,8 @@ end # Reset range of cells to values that are prone to cause errors as soon as they are used. # # Rationale: If an invalid cell is accidentally used, we want to know it as soon as possible. -function invalidate!(t::SerialTree, first::Int, last::Int) +function invalidate!(t::SerialTree{NDIMS, RealT}, + first::Int, last::Int) where {NDIMS, RealT <: Real} @assert first > 0 @assert last <= t.capacity + 1 @@ -179,7 +191,7 @@ function invalidate!(t::SerialTree, first::Int, last::Int) t.child_ids[:, first:last] .= typemin(Int) t.neighbor_ids[:, first:last] .= typemin(Int) t.levels[first:last] .= typemin(Int) - t.coordinates[:, first:last] .= NaN + t.coordinates[:, first:last] .= convert(RealT, NaN) # `NaN` is of type Float64 t.original_cell_ids[first:last] .= typemin(Int) return nothing @@ -203,12 +215,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) diff --git a/src/meshes/tree_mesh.jl b/src/meshes/tree_mesh.jl index 458a43203d2..885e1e3e724 100644 --- a/src/meshes/tree_mesh.jl +++ b/src/meshes/tree_mesh.jl @@ -24,16 +24,19 @@ get_name(mesh::AbstractMesh) = mesh |> typeof |> nameof |> string A Cartesian mesh based on trees of hypercubes to support adaptive mesh refinement. """ -mutable struct TreeMesh{NDIMS, TreeType <: AbstractTree{NDIMS}} <: AbstractMesh{NDIMS} +mutable struct TreeMesh{NDIMS, TreeType <: AbstractTree{NDIMS}, RealT <: Real} <: + AbstractMesh{NDIMS} tree::TreeType current_filename::String unsaved_changes::Bool first_cell_by_rank::OffsetVector{Int, Vector{Int}} n_cells_by_rank::OffsetVector{Int, Vector{Int}} - function TreeMesh{NDIMS, TreeType}(n_cells_max::Integer) where {NDIMS, - TreeType <: - AbstractTree{NDIMS}} + function TreeMesh{NDIMS, TreeType, RealT}(n_cells_max::Integer) where {NDIMS, + TreeType <: + AbstractTree{NDIMS}, + RealT <: + Real} # Create mesh m = new() m.tree = TreeType(n_cells_max) @@ -46,14 +49,13 @@ mutable struct TreeMesh{NDIMS, TreeType <: AbstractTree{NDIMS}} <: AbstractMesh{ end # TODO: Taal refactor, order of important arguments, use of n_cells_max? - # TODO: Taal refactor, allow other RealT for the mesh, not just Float64 - # TODO: Taal refactor, use NTuple instead of domain_center::AbstractArray{Float64} - function TreeMesh{NDIMS, TreeType}(n_cells_max::Integer, - domain_center::AbstractArray{Float64}, - domain_length, - periodicity = true) where {NDIMS, - TreeType <: - AbstractTree{NDIMS}} + function TreeMesh{NDIMS, TreeType, RealT}(n_cells_max::Integer, + domain_center::AbstractArray{RealT}, + domain_length::RealT, + periodicity = true) where {NDIMS, + TreeType <: + AbstractTree{NDIMS}, + RealT <: Real} @assert NDIMS isa Integer && NDIMS > 0 # Create mesh @@ -81,27 +83,29 @@ const ParallelTreeMesh{NDIMS} = TreeMesh{NDIMS, <:ParallelTree{NDIMS}} partition!(mesh::SerialTreeMesh) = nothing # Constructor for passing the dimension and mesh type as an argument -function TreeMesh(::Type{TreeType}, - args...) where {NDIMS, TreeType <: AbstractTree{NDIMS}} - TreeMesh{NDIMS, TreeType}(args...) +function TreeMesh(::Type{TreeType}, args...; + RealT = Float64) where {NDIMS, TreeType <: AbstractTree{NDIMS}} + TreeMesh{NDIMS, TreeType, RealT}(args...) end # Constructor accepting a single number as center (as opposed to an array) for 1D -function TreeMesh{1, TreeType}(n::Int, center::Real, len::Real, - periodicity = true) where {TreeType <: AbstractTree{1}} - # TODO: Taal refactor, allow other RealT for the mesh, not just Float64 - return TreeMesh{1, TreeType}(n, SVector{1, Float64}(center), len, periodicity) +function TreeMesh{1, TreeType, RealT}(n::Int, center::RealT, len::RealT, + periodicity = true) where { + TreeType <: + AbstractTree{1}, + RealT <: Real} + return TreeMesh{1, TreeType, RealT}(n, SVector{1, RealT}(center), len, periodicity) end -function TreeMesh{NDIMS, TreeType}(n_cells_max::Integer, - domain_center::NTuple{NDIMS, Real}, - domain_length::Real, - periodicity = true) where {NDIMS, - TreeType <: - AbstractTree{NDIMS}} - # TODO: Taal refactor, allow other RealT for the mesh, not just Float64 - TreeMesh{NDIMS, TreeType}(n_cells_max, SVector{NDIMS, Float64}(domain_center), - convert(Float64, domain_length), periodicity) +function TreeMesh{NDIMS, TreeType, RealT}(n_cells_max::Integer, + domain_center::NTuple{NDIMS, RealT}, + domain_length::RealT, + periodicity = true) where {NDIMS, + TreeType <: + AbstractTree{NDIMS}, + RealT <: Real} + TreeMesh{NDIMS, TreeType, RealT}(n_cells_max, SVector{NDIMS, RealT}(domain_center), + domain_length, periodicity) end function TreeMesh(coordinates_min::NTuple{NDIMS, Real}, @@ -110,7 +114,8 @@ function TreeMesh(coordinates_min::NTuple{NDIMS, Real}, periodicity = true, initial_refinement_level, refinement_patches = (), - coarsening_patches = ()) where {NDIMS} + coarsening_patches = (), + RealT = Float64) where {NDIMS} # check arguments if !(n_cells_max isa Integer && n_cells_max > 0) throw(ArgumentError("`n_cells_max` must be a positive integer (provided `n_cells_max = $n_cells_max`)")) @@ -119,9 +124,19 @@ 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 + # Check if elements in coordinates_min and coordinates_max are all of type RealT + for i in 1:NDIMS + if !(coordinates_min[i] isa RealT) + @warn "Element $i in `coordinates_min` is not of type $RealT (provided `coordinates_min[$i] = $(coordinates_min[i])`)" + end + if !(coordinates_max[i] isa RealT) + @warn "Element $i in `coordinates_max` is not of type $RealT (provided `coordinates_max[$i] = $(coordinates_max[i])`)" + end + end + # TreeMesh requires equal domain lengths in all dimensions - domain_center = @. (coordinates_min + coordinates_max) / 2 - domain_length = coordinates_max[1] - coordinates_min[1] + domain_center = @. convert(RealT, (coordinates_min + coordinates_max) / 2) + domain_length = convert(RealT, 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 @@ -139,10 +154,10 @@ function TreeMesh(coordinates_min::NTuple{NDIMS, Real}, end # Create mesh - mesh = @trixi_timeit timer() "creation" TreeMesh{NDIMS, TreeType}(n_cells_max, - domain_center, - domain_length, - periodicity) + mesh = @trixi_timeit timer() "creation" TreeMesh{NDIMS, TreeType, RealT}(n_cells_max, + domain_center, + domain_length, + periodicity) # Initialize mesh initialize!(mesh, initial_refinement_level, refinement_patches, coarsening_patches) @@ -184,7 +199,8 @@ function initialize!(mesh::TreeMesh, initial_refinement_level, return nothing end -function TreeMesh(coordinates_min::Real, coordinates_max::Real; kwargs...) +function TreeMesh(coordinates_min::Real, coordinates_max::Real; + kwargs...) TreeMesh((coordinates_min,), (coordinates_max,); kwargs...) end diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index ed6b7f6b527..cf70d0eef5b 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -148,6 +148,36 @@ end end end +@trixi_testset "elixir_mhd_ec_float32.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_float32.jl"), + l2=Float32[0.036355644, + 0.042947736, + 0.04294775, + 0.025748005, + 0.16211236, + 0.017452478, + 0.017452475, + 0.026877576, + 2.0951437f-7], + linf=Float32[0.22100884, + 0.28798944, + 0.2879896, + 0.20858234, + 0.8812654, + 0.09208202, + 0.09208143, + 0.14795381, + 2.0912607f-6]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_mhd_orszag_tang.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"), l2=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index b0925f31cd8..641c7244efe 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -27,6 +27,7 @@ isdir(outdir) && rm(outdir, recursive = true) @timed_testset "SerialTree" begin @testset "constructors" begin @test_nowarn Trixi.SerialTree(Val(1), 10, 0.0, 1.0) + @test_nowarn Trixi.SerialTree{1}(10, 0.0, 1.0) end @testset "helper functions" begin @@ -57,6 +58,7 @@ end @timed_testset "ParallelTree" begin @testset "constructors" begin @test_nowarn Trixi.ParallelTree(Val(1), 10, 0.0, 1.0) + @test_nowarn Trixi.ParallelTree{1}(10, 0.0, 1.0) end @testset "helper functions" begin @@ -68,7 +70,8 @@ end @timed_testset "TreeMesh" begin @testset "constructors" begin - @test TreeMesh{1, Trixi.SerialTree{1}}(1, 5.0, 2.0) isa TreeMesh + @test TreeMesh{1, Trixi.SerialTree{1, Float64}, Float64}(1, 5.0, 2.0) isa + TreeMesh # Invalid domain length check (TreeMesh expects a hypercube) # 2D @@ -89,7 +92,9 @@ end let @test Trixi.mpi_nranks() == 2 - mesh = TreeMesh{2, Trixi.ParallelTree{2}}(30, (0.0, 0.0), 1) + mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(30, + (0.0, 0.0), + 1.0) # Refine twice Trixi.refine!(mesh.tree) Trixi.refine!(mesh.tree) @@ -117,7 +122,9 @@ end let @test Trixi.mpi_nranks() == 3 - mesh = TreeMesh{2, Trixi.ParallelTree{2}}(100, (0.0, 0.0), 1) + mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(100, + (0.0, 0.0), + 1.0) # Refine twice Trixi.refine!(mesh.tree) Trixi.refine!(mesh.tree) @@ -145,7 +152,9 @@ end let @test Trixi.mpi_nranks() == 9 - mesh = TreeMesh{2, Trixi.ParallelTree{2}}(1000, (0.0, 0.0), 1) + mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(1000, + (0.0, 0.0), + 1.0) # Refine twice Trixi.refine!(mesh.tree) Trixi.refine!(mesh.tree) @@ -168,7 +177,9 @@ end let @test Trixi.mpi_nranks() == 3 - mesh = TreeMesh{2, Trixi.ParallelTree{2}}(100, (0.0, 0.0), 1) + mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(100, + (0.0, 0.0), + 1.0) # Refine whole tree Trixi.refine!(mesh.tree) # Refine left leaf @@ -195,7 +206,9 @@ end let @test Trixi.mpi_nranks() == 3 - mesh = TreeMesh{2, Trixi.ParallelTree{2}}(100, (0.0, 0.0), 1) + mesh = TreeMesh{2, Trixi.ParallelTree{2, Float64}, Float64}(100, + (0.0, 0.0), + 1.0) # Only one leaf @test_throws AssertionError("Too many ranks to properly partition the mesh!") Trixi.partition!(mesh)