Skip to content

Commit

Permalink
Make TreeMesh type general (trixi-framework#2129)
Browse files Browse the repository at this point in the history 
* 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 <[email protected]>
  • Loading branch information
@DanielDoehring @ranocha
DanielDoehring and ranocha authored Nov 7, 2024
1 parent dcf1b58 commit 9fa4e27
Show file tree
Hide file tree
Showing 11 changed files with 283 additions and 127 deletions.
4 changes: 2 additions & 2 deletions examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_sedov_blast_wave.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
4 changes: 2 additions & 2 deletions examples/tree_1d_dgsem/elixir_advection_diffusion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
62 changes: 62 additions & 0 deletions examples/tree_2d_dgsem/elixir_mhd_ec_float32.jl
View file
Original file line number Diff line number Diff line change
@@ -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
50 changes: 26 additions & 24 deletions src/auxiliary/precompile.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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},
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
26 changes: 15 additions & 11 deletions src/meshes/abstract_tree.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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."
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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.
Expand Down
7 changes: 5 additions & 2 deletions src/meshes/mesh_io.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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()
Expand Down
Loading

0 comments on commit 9fa4e27

Please sign in to comment.