From 9fa4e27f38c7cfee77bf871dbb3c42de514b751d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 7 Nov 2024 20:57:47 +0100 Subject: [PATCH 1/8] 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) From 806c1b9b76586de7cd15a691ed9e7f65f997fe2d Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Fri, 8 Nov 2024 01:18:43 -1000 Subject: [PATCH 2/8] Type stability of functions from examples tree_1d_dgsem (#2133) * Start * Start again * Add more * Complete * Fix * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * Fix --------- Co-authored-by: Hendrik Ranocha --- .../elixir_advection_diffusion.jl | 14 ++--- .../elixir_burgers_linear_stability.jl | 3 +- .../elixir_burgers_rarefaction.jl | 3 +- .../tree_1d_dgsem/elixir_burgers_shock.jl | 3 +- .../tree_1d_dgsem/elixir_euler_blast_wave.jl | 9 +-- .../tree_1d_dgsem/elixir_euler_positivity.jl | 15 ++--- .../elixir_euler_quasi_1d_discontinuous.jl | 9 +-- .../tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl | 9 +-- .../elixir_euler_sedov_blast_wave.jl | 15 ++--- .../elixir_euler_sedov_blast_wave_pure_fv.jl | 15 ++--- ..._eulermulti_two_interacting_blast_waves.jl | 17 ++--- .../elixir_hypdiff_harmonic_nonperiodic.jl | 9 +-- .../elixir_linearizedeuler_gauss_wall.jl | 2 +- .../elixir_maxwell_E_excitation.jl | 5 +- .../elixir_mhd_briowu_shock_tube.jl | 17 ++--- .../elixir_mhd_ryujones_shock_tube.jl | 15 ++--- .../elixir_mhd_shu_osher_shock_tube.jl | 38 ++++++------ .../elixir_mhd_torrilhon_shock_tube.jl | 17 ++--- .../elixir_mhdmulti_briowu_shock_tube.jl | 35 +++++------ ...lixir_navierstokes_convergence_periodic.jl | 50 ++++++++------- .../elixir_navierstokes_convergence_walls.jl | 62 ++++++++++--------- ...ixir_navierstokes_convergence_walls_amr.jl | 62 ++++++++++--------- .../tree_1d_dgsem/elixir_shallowwater_ec.jl | 15 ++--- ...xir_shallowwater_quasi_1d_discontinuous.jl | 13 ++-- ...xir_shallowwater_quasi_1d_well_balanced.jl | 4 +- .../elixir_shallowwater_shock_capturing.jl | 18 +++--- .../elixir_shallowwater_well_balanced.jl | 8 +-- ..._shallowwater_well_balanced_nonperiodic.jl | 4 +- .../elixir_traffic_flow_lwr_trafficjam.jl | 3 +- 29 files changed, 257 insertions(+), 232 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl index e9530692d4b..9182fd22d0e 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl @@ -22,10 +22,12 @@ mesh = TreeMesh(coordinates_min, coordinates_max, n_cells_max = 30_000, # set maximum capacity of tree data structure periodicity = true) -function x_trans_periodic(x, domain_length = SVector(2 * pi), center = SVector(0.0)) +function x_trans_periodic(x, domain_length = SVector(oftype(x[1], 2 * pi)), + center = SVector(oftype(x[1], 0))) x_normalized = x .- center x_shifted = x_normalized .% domain_length - x_offset = ((x_shifted .< -0.5 * domain_length) - (x_shifted .> 0.5 * domain_length)) .* + x_offset = ((x_shifted .< -0.5f0 * domain_length) - + (x_shifted .> 0.5f0 * domain_length)) .* domain_length return center + x_shifted + x_offset end @@ -37,11 +39,9 @@ function initial_condition_diffusive_convergence_test(x, t, x_trans = x_trans_periodic(x - equation.advection_velocity * t) nu = diffusivity() - c = 0.0 - A = 1.0 - L = 2 - f = 1 / L - omega = 1.0 + c = 0 + A = 1 + omega = 1 scalar = c + A * sin(omega * sum(x_trans)) * exp(-nu * omega^2 * t) return SVector(scalar) end diff --git a/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl b/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl index ae2039edde6..438356144cc 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_linear_stability.jl @@ -8,8 +8,9 @@ using Trixi equations = InviscidBurgersEquation1D() function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D) + RealT = eltype(x) k = 1 - 2 + sinpi(k * (x[1] - 0.7)) |> SVector + SVector(2 + sinpi(k * (x[1] - convert(RealT, 0.7)))) end volume_flux = flux_ec diff --git a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl index 843d8a7aa4c..5e40baf49cc 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl @@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max, # Discontinuous initial condition (Riemann Problem) leading to a rarefaction fan. function initial_condition_rarefaction(x, t, equation::InviscidBurgersEquation1D) - scalar = x[1] < 0.5 ? 0.5 : 1.5 + RealT = eltype(x) + scalar = x[1] < 0.5f0 ? convert(RealT, 0.5f0) : convert(RealT, 1.5f0) return SVector(scalar) end diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl index 09c6cee9f98..b795e8623fa 100644 --- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl +++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl @@ -35,7 +35,8 @@ mesh = TreeMesh(coordinate_min, coordinate_max, # Discontinuous initial condition (Riemann Problem) leading to a shock to test e.g. correct shock speed. function initial_condition_shock(x, t, equation::InviscidBurgersEquation1D) - scalar = x[1] < 0.5 ? 1.5 : 0.5 + RealT = eltype(x) + scalar = x[1] < 0.5f0 ? convert(RealT, 1.5f0) : convert(RealT, 0.5f0) return SVector(scalar) end diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl index 9cba4936d22..5b5eb6fdd33 100644 --- a/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl +++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl @@ -18,7 +18,8 @@ A medium blast wave taken from function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # The following code is equivalent to @@ -28,9 +29,9 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - p = r > 0.5 ? 1.0E-3 : 1.245 + rho = r > 0.5f0 ? one(RealT) : RealT(1.1691) + v1 = r > 0.5f0 ? zero(RealT) : RealT(0.1882) * cos_phi + p = r > 0.5f0 ? RealT(1.0E-3) : RealT(1.245) return prim2cons(SVector(rho, v1, p), equations) end diff --git a/examples/tree_1d_dgsem/elixir_euler_positivity.jl b/examples/tree_1d_dgsem/elixir_euler_positivity.jl index 3fb96fb807b..34a6ad16a43 100644 --- a/examples/tree_1d_dgsem/elixir_euler_positivity.jl +++ b/examples/tree_1d_dgsem/elixir_euler_positivity.jl @@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 - r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup - E = 1.0 - p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0) + p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables - rho = 1.0 - v1 = 0.0 + rho = 1 + v1 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, p), equations) diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl index cc4535be028..0af5b210ee9 100644 --- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl +++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl @@ -18,10 +18,11 @@ A discontinuous initial condition taken from """ function initial_condition_discontinuity(x, t, equations::CompressibleEulerEquationsQuasi1D) - rho = (x[1] < 0) ? 3.4718 : 2.0 - v1 = (x[1] < 0) ? -2.5923 : -3.0 - p = (x[1] < 0) ? 5.7118 : 2.639 - a = (x[1] < 0) ? 1.0 : 1.5 + RealT = eltype(x) + rho = (x[1] < 0) ? RealT(3.4718) : RealT(2.0) + v1 = (x[1] < 0) ? RealT(-2.5923) : RealT(-3.0) + p = (x[1] < 0) ? RealT(5.7118) : RealT(2.639) + a = (x[1] < 0) ? 1.0f0 : 1.5f0 return prim2cons(SVector(rho, v1, p, a), equations) end diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl index ae1b2b24b62..ebd3680cef4 100644 --- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl +++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl @@ -15,10 +15,11 @@ equations = CompressibleEulerEquationsQuasi1D(1.4) # refinement level is changed the initial condition below may need changed as well to # ensure that the discontinuities lie on an element interface. function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D) - v1 = 0.1 - rho = 2.0 + 0.1 * x[1] - p = 3.0 - a = 2.0 + x[1] + RealT = eltype(x) + v1 = convert(RealT, 0.1) + rho = 2 + convert(RealT, 0.1) * x[1] + p = 3 + a = 2 + x[1] return prim2cons(SVector(rho, v1, p, a), equations) end diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl index b67b2bc602e..a89727da147 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave.jl @@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 - r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup - E = 1.0 - p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0) + p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables - rho = 1.0 - v1 = 0.0 + rho = 1 + v1 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, p), equations) diff --git a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl index 8a0241680b9..cf33f7f5b64 100644 --- a/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl +++ b/examples/tree_1d_dgsem/elixir_euler_sedov_blast_wave_pure_fv.jl @@ -15,21 +15,22 @@ The Sedov blast wave setup based on Flash """ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D) # Set up polar coordinates - inicenter = SVector(0.0) + RealT = eltype(x) + inicenter = SVector(0) x_norm = x[1] - inicenter[1] r = abs(x_norm) # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 - r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + r0 = 0.21875f0 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) # r0 = 0.5 # = more reasonable setup - E = 1.0 - p0_inner = 6 * (equations.gamma - 1) * E / (3 * pi * r0) - p0_outer = 1.0e-5 # = true Sedov setup + E = 1 + p0_inner = 6 * (equations.gamma - 1) * E / (3 * convert(RealT, pi) * r0) + p0_outer = convert(RealT, 1.0e-5) # = true Sedov setup # p0_outer = 1.0e-3 # = more reasonable setup # Calculate primitive variables - rho = 1.0 - v1 = 0.0 + rho = 1 + v1 = 0 p = r > r0 ? p0_outer : p0_inner return prim2cons(SVector(rho, v1, p), equations) diff --git a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl index c093420cc7c..3089e0c6809 100644 --- a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl +++ b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl @@ -18,20 +18,21 @@ A multicomponent two interacting blast wave test taken from """ function initial_condition_two_interacting_blast_waves(x, t, equations::CompressibleEulerMulticomponentEquations1D) - rho1 = 0.5 * x[1]^2 - rho2 = 0.5 * (sin(20 * x[1]))^2 + RealT = eltype(x) + rho1 = 0.5f0 * x[1]^2 + rho2 = 0.5f0 * (sin(20 * x[1]))^2 rho3 = 1 - rho1 - rho2 prim_rho = SVector{3, real(equations)}(rho1, rho2, rho3) - v1 = 0.0 + v1 = 0 - if x[1] <= 0.1 - p = 1000 - elseif x[1] < 0.9 - p = 0.01 + if x[1] <= RealT(0.1) + p = convert(RealT, 1000) + elseif x[1] < RealT(0.9) + p = convert(RealT, 0.01) else - p = 100 + p = convert(RealT, 100) end prim_other = SVector{2, real(equations)}(v1, p) diff --git a/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl index b9173ec9f49..dda1009cfa8 100644 --- a/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_hypdiff_harmonic_nonperiodic.jl @@ -19,12 +19,13 @@ A non-priodic harmonic function used in combination with function initial_condition_harmonic_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D) # elliptic equation: -νΔϕ = f - if t == 0.0 - phi = 5.0 - q1 = 0.0 + RealT = eltype(x) + if t == 0 + phi = convert(RealT, 5) + q1 = zero(RealT) else A = 3 - B = exp(1) + B = exp(one(RealT)) phi = A + B * x[1] q1 = B end diff --git a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl index 0884249559a..7a7b83dfe84 100644 --- a/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -22,7 +22,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # that is advected to left with v - c and to the right with v + c. # Correspondigly, the bump splits in half. function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations1D) - v1_prime = 0.0 + v1_prime = 0 rho_prime = p_prime = 2 * exp(-(x[1] - 45)^2 / 25) return SVector(rho_prime, v1_prime, p_prime) end diff --git a/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl b/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl index b5478331a7d..da44be30da2 100644 --- a/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl +++ b/examples/tree_1d_dgsem/elixir_maxwell_E_excitation.jl @@ -20,9 +20,10 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # Excite the electric field which causes a standing wave. # The solution is an undamped exchange between electric and magnetic energy. function initial_condition_E_excitation(x, t, equations::MaxwellEquations1D) + RealT = eltype(x) c = equations.speed_of_light - E = -c * sin(2 * pi * x[1]) - B = 0.0 + E = -c * sinpi(2 * x[1]) + B = 0 return SVector(E, B) end diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl index 4d537508b47..8f33ad77cfc 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl @@ -18,14 +18,15 @@ MHD extension of the Sod shock tube. Taken from Section V of the article """ function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [0, 1], γ = 2, final time = 0.12 - rho = x[1] < 0.5 ? 1.0 : 0.125 - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - p = x[1] < 0.5 ? 1.0 : 0.1 - B1 = 0.75 - B2 = x[1] < 0.5 ? 1.0 : -1.0 - B3 = 0.0 + RealT = eltype(x) + rho = x[1] < 0.5f0 ? 1.0f0 : 0.125f0 + v1 = 0 + v2 = 0 + v3 = 0 + p = x[1] < 0.5f0 ? one(RealT) : RealT(0.1) + B1 = 0.75f0 + B2 = x[1] < 0.5f0 ? 1 : -1 + B3 = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end initial_condition = initial_condition_briowu_shock_tube diff --git a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl index a7d7689a806..fc25ef87df7 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl @@ -23,14 +23,15 @@ present in the one dimensional MHD equations. It is the second test from Section """ function initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [0, 1], γ = 5/3, final time = 0.2 - rho = x[1] <= 0.5 ? 1.08 : 1.0 - v1 = x[1] <= 0.5 ? 1.2 : 0.0 - v2 = x[1] <= 0.5 ? 0.01 : 0.0 - v3 = x[1] <= 0.5 ? 0.5 : 0.0 - p = x[1] <= 0.5 ? 0.95 : 1.0 - inv_sqrt4pi = 1.0 / sqrt(4 * pi) + RealT = eltype(x) + rho = x[1] <= 0.5f0 ? RealT(1.08) : one(RealT) + v1 = x[1] <= 0.5f0 ? RealT(1.2) : zero(RealT) + v2 = x[1] <= 0.5f0 ? RealT(0.01) : zero(RealT) + v3 = x[1] <= 0.5f0 ? 0.5f0 : 0.0f0 + p = x[1] <= 0.5f0 ? RealT(0.95) : one(RealT) + inv_sqrt4pi = 1 / sqrt(4 * convert(RealT, pi)) B1 = 2 * inv_sqrt4pi - B2 = x[1] <= 0.5 ? 3.6 * inv_sqrt4pi : 4.0 * inv_sqrt4pi + B2 = x[1] <= 0.5f0 ? RealT(3.6) * inv_sqrt4pi : 4 * inv_sqrt4pi B3 = B1 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) diff --git a/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl index 689537ebdd4..1ab4c54985b 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl @@ -19,15 +19,16 @@ Taken from Section 4.1 of function initial_condition_shu_osher_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [-5, 5], γ = 5/3, final time = 0.7 # initial shock location is taken to be at x = -4 - x_0 = -4.0 - rho = x[1] <= x_0 ? 3.5 : 1.0 + 0.2 * sin(5.0 * x[1]) - v1 = x[1] <= x_0 ? 5.8846 : 0.0 - v2 = x[1] <= x_0 ? 1.1198 : 0.0 - v3 = 0.0 - p = x[1] <= x_0 ? 42.0267 : 1.0 - B1 = 1.0 - B2 = x[1] <= x_0 ? 3.6359 : 1.0 - B3 = 0.0 + RealT = eltype(x) + x_0 = -4 + rho = x[1] <= x_0 ? RealT(3.5) : 1 + RealT(0.2) * sin(5 * x[1]) + v1 = x[1] <= x_0 ? RealT(5.8846) : zero(RealT) + v2 = x[1] <= x_0 ? RealT(1.1198) : zero(RealT) + v3 = 0 + p = x[1] <= x_0 ? RealT(42.0267) : one(RealT) + B1 = 1 + B2 = x[1] <= x_0 ? RealT(3.6359) : one(RealT) + B3 = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end @@ -45,15 +46,16 @@ function initial_condition_shu_osher_shock_tube_flipped(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [-5, 5], γ = 5/3, final time = 0.7 # initial shock location is taken to be at x = 4 - x_0 = 4.0 - rho = x[1] <= x_0 ? 1.0 + 0.2 * sin(5.0 * x[1]) : 3.5 - v1 = x[1] <= x_0 ? 0.0 : -5.8846 - v2 = x[1] <= x_0 ? 0.0 : -1.1198 - v3 = 0.0 - p = x[1] <= x_0 ? 1.0 : 42.0267 - B1 = 1.0 - B2 = x[1] <= x_0 ? 1.0 : 3.6359 - B3 = 0.0 + RealT = eltype(x) + x_0 = 4 + rho = x[1] <= x_0 ? 1 + RealT(0.2) * sin(5 * x[1]) : RealT(3.5) + v1 = x[1] <= x_0 ? zero(RealT) : RealT(-5.8846) + v2 = x[1] <= x_0 ? zero(RealT) : RealT(-1.1198) + v3 = 0 + p = x[1] <= x_0 ? one(RealT) : RealT(42.0267) + B1 = 1 + B2 = x[1] <= x_0 ? one(RealT) : RealT(3.6359) + B3 = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end diff --git a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl index 3b366c35e0f..c35212b4074 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_torrilhon_shock_tube.jl @@ -17,14 +17,15 @@ Torrilhon's shock tube test case for one dimensional ideal MHD equations. """ function initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D) # domain must be set to [-1, 1.5], γ = 5/3, final time = 0.4 - rho = x[1] <= 0 ? 3.0 : 1.0 - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - p = x[1] <= 0 ? 3.0 : 1.0 - B1 = 1.5 - B2 = x[1] <= 0 ? 1.0 : cos(1.5) - B3 = x[1] <= 0 ? 0.0 : sin(1.5) + RealT = eltype(x) + rho = x[1] <= 0 ? 3 : 1 + v1 = 0 + v2 = 0 + v3 = 0 + p = x[1] <= 0 ? 3 : 1 + B1 = 1.5f0 + B2 = x[1] <= 0 ? one(RealT) : cos(RealT(1.5f0)) + B3 = x[1] <= 0 ? zero(RealT) : sin(RealT(1.5f0)) return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) end initial_condition = initial_condition_torrilhon_shock_tube diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl index 831fa7afedb..999e6c42b29 100644 --- a/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhdmulti_briowu_shock_tube.jl @@ -19,28 +19,27 @@ MHD extension of the Sod shock tube. Taken from Section V of the article function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdMulticomponentEquations1D) # domain must be set to [0, 1], γ = 2, final time = 0.12 - if x[1] < 0.5 - rho = 1.0 - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + RealT = eltype(x) + if x[1] < 0.5f0 + rho = one(RealT) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) else - rho = 0.125 - prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - rho + rho = RealT(0.125) + prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) * + rho / (1 - + 2^ncomponents(equations)) for i in eachcomponent(equations)) end - v1 = 0.0 - v2 = 0.0 - v3 = 0.0 - p = x[1] < 0.5 ? 1.0 : 0.1 - B1 = 0.75 - B2 = x[1] < 0.5 ? 1.0 : -1.0 - B3 = 0.0 + v1 = 0 + v2 = 0 + v3 = 0 + p = x[1] < 0.5f0 ? one(RealT) : RealT(0.1) + B1 = 0.75f0 + B2 = x[1] < 0.5f0 ? 1 : -1 + B3 = 0 prim_other = SVector(v1, v2, v3, p, B1, B2, B3) return prim2cons(vcat(prim_other, prim_rho), equations) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl index eab0840f385..1089cf4bb85 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl @@ -16,12 +16,13 @@ equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu = mu(), # (Simplified version of the 2D) function initial_condition_navier_stokes_convergence_test(x, t, equations) # Amplitude and shift - A = 0.5 - c = 2.0 + RealT = eltype(x) + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x[1] - pi_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t rho = c + A * sin(pi_x) * cos(pi_t) v1 = sin(pi_x) * cos(pi_t) @@ -34,42 +35,43 @@ initial_condition = initial_condition_navier_stokes_convergence_test @inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) # we currently need to hardcode these parameters until we fix the "combined equation" issue # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + RealT = eltype(x) inv_gamma_minus_one = inv(equations.gamma - 1) Pr = prandtl_number() mu_ = mu() # Same settings as in `initial_condition` # Amplitude and shift - A = 0.5 - c = 2.0 + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x[1] - pi_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t # compute the manufactured solution and all necessary derivatives rho = c + A * sin(pi_x) * cos(pi_t) - rho_t = -pi * A * sin(pi_x) * sin(pi_t) - rho_x = pi * A * cos(pi_x) * cos(pi_t) - rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_t) + rho_t = -convert(RealT, pi) * A * sin(pi_x) * sin(pi_t) + rho_x = convert(RealT, pi) * A * cos(pi_x) * cos(pi_t) + rho_xx = -convert(RealT, pi) * convert(RealT, pi) * A * sin(pi_x) * cos(pi_t) v1 = sin(pi_x) * cos(pi_t) - v1_t = -pi * sin(pi_x) * sin(pi_t) - v1_x = pi * cos(pi_x) * cos(pi_t) - v1_xx = -pi * pi * sin(pi_x) * cos(pi_t) + v1_t = -convert(RealT, pi) * sin(pi_x) * sin(pi_t) + v1_x = convert(RealT, pi) * cos(pi_x) * cos(pi_t) + v1_xx = -convert(RealT, pi) * convert(RealT, pi) * sin(pi_x) * cos(pi_t) p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x - E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 - E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + E = p * inv_gamma_minus_one + 0.5f0 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5f0 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5f0 * rho_x * v1^2 + rho * v1 * v1_x # Some convenience constants T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_rho_cubed = 1.0 / (rho^3) + inv_rho_cubed = 1 / (rho^3) # compute the source terms # density equation @@ -77,7 +79,7 @@ initial_condition = initial_condition_navier_stokes_convergence_test # x-momentum equation du2 = (rho_t * v1 + rho * v1_t - + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - + + p_x + rho_x * v1^2 + 2 * rho * v1 * v1_x - # stress tensor from x-direction v1_xx * mu_) @@ -88,8 +90,8 @@ initial_condition = initial_condition_navier_stokes_convergence_test v1_x * v1_x * mu_ - T_const * inv_rho_cubed * (p_xx * rho * rho - - 2.0 * p_x * rho * rho_x + - 2.0 * p * rho_x * rho_x - + 2 * p_x * rho * rho_x + + 2 * p * rho_x * rho_x - p * rho * rho_xx) * mu_) return SVector(du1, du2, du3) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl index 40030d53345..d534b7b189b 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl @@ -31,21 +31,23 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) function initial_condition_navier_stokes_convergence_test(x, t, equations) # Amplitude and shift - A = 0.5 - c = 2.0 + RealT = eltype(x) + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x[1] - pi_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t rho = c + A * cos(pi_x) * cos(pi_t) - v1 = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0))) * cos(pi_t) + v1 = log(x[1] + 2) * (1 - exp(-A * (x[1] - 1))) * cos(pi_t) p = rho^2 return prim2cons(SVector(rho, v1, p), equations) end @inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + RealT = eltype(x) x = x[1] # TODO: parabolic @@ -57,40 +59,40 @@ end # Same settings as in `initial_condition` # Amplitude and shift - A = 0.5 - c = 2.0 + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x - pi_t = pi * t + pi_x = convert(RealT, pi) * x + pi_t = convert(RealT, pi) * t # compute the manufactured solution and all necessary derivatives rho = c + A * cos(pi_x) * cos(pi_t) - rho_t = -pi * A * cos(pi_x) * sin(pi_t) - rho_x = -pi * A * sin(pi_x) * cos(pi_t) - rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t) - - v1 = log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t) - v1_t = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t) - v1_x = (A * log(x + 2.0) * exp(-A * (x - 1.0)) + - (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t) - v1_xx = ((2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0) - - A * A * log(x + 2.0) * exp(-A * (x - 1.0)) - - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t)) + rho_t = -convert(RealT, pi) * A * cos(pi_x) * sin(pi_t) + rho_x = -convert(RealT, pi) * A * sin(pi_x) * cos(pi_t) + rho_xx = -convert(RealT, pi) * convert(RealT, pi) * A * cos(pi_x) * cos(pi_t) + + v1 = log(x + 2) * (1 - exp(-A * (x - 1))) * cos(pi_t) + v1_t = -convert(RealT, pi) * log(x + 2) * (1 - exp(-A * (x - 1))) * sin(pi_t) + v1_x = (A * log(x + 2) * exp(-A * (x - 1)) + + (1 - exp(-A * (x - 1))) / (x + 2)) * cos(pi_t) + v1_xx = ((2 * A * exp(-A * (x - 1)) / (x + 2) - + A * A * log(x + 2) * exp(-A * (x - 1)) - + (1 - exp(-A * (x - 1))) / ((x + 2) * (x + 2))) * cos(pi_t)) p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 - E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + E = p * inv_gamma_minus_one + 0.5f0 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5f0 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5f0 * rho_x * v1^2 + rho * v1 * v1_x # Some convenience constants T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_rho_cubed = 1.0 / (rho^3) + inv_rho_cubed = 1 / (rho^3) # compute the source terms # density equation @@ -98,7 +100,7 @@ end # y-momentum equation du2 = (rho_t * v1 + rho * v1_t - + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - + + p_x + rho_x * v1^2 + 2 * rho * v1 * v1_x - # stress tensor from y-direction v1_xx * mu_) @@ -109,8 +111,8 @@ end v1_x * v1_x * mu_ - T_const * inv_rho_cubed * (p_xx * rho * rho - - 2.0 * p_x * rho * rho_x + - 2.0 * p * rho_x * rho_x - + 2 * p_x * rho * rho_x + + 2 * p * rho_x * rho_x - p * rho * rho_xx) * mu_) return SVector(du1, du2, du3) diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl index e833155a68e..4f9cd6365fa 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl @@ -31,21 +31,23 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) function initial_condition_navier_stokes_convergence_test(x, t, equations) # Amplitude and shift - A = 0.5 - c = 2.0 + RealT = eltype(x) + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x[1] - pi_t = pi * t + pi_x = convert(RealT, pi) * x[1] + pi_t = convert(RealT, pi) * t rho = c + A * cos(pi_x) * cos(pi_t) - v1 = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0))) * cos(pi_t) + v1 = log(x[1] + 2) * (1 - exp(-A * (x[1] - 1))) * cos(pi_t) p = rho^2 return prim2cons(SVector(rho, v1, p), equations) end @inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + RealT = eltype(x) x = x[1] # TODO: parabolic @@ -57,40 +59,40 @@ end # Same settings as in `initial_condition` # Amplitude and shift - A = 0.5 - c = 2.0 + A = 0.5f0 + c = 2 # convenience values for trig. functions - pi_x = pi * x - pi_t = pi * t + pi_x = convert(RealT, pi) * x + pi_t = convert(RealT, pi) * t # compute the manufactured solution and all necessary derivatives rho = c + A * cos(pi_x) * cos(pi_t) - rho_t = -pi * A * cos(pi_x) * sin(pi_t) - rho_x = -pi * A * sin(pi_x) * cos(pi_t) - rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t) - - v1 = log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t) - v1_t = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t) - v1_x = (A * log(x + 2.0) * exp(-A * (x - 1.0)) + - (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t) - v1_xx = ((2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0) - - A * A * log(x + 2.0) * exp(-A * (x - 1.0)) - - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t)) + rho_t = -convert(RealT, pi) * A * cos(pi_x) * sin(pi_t) + rho_x = -convert(RealT, pi) * A * sin(pi_x) * cos(pi_t) + rho_xx = -convert(RealT, pi) * convert(RealT, pi) * A * cos(pi_x) * cos(pi_t) + + v1 = log(x + 2) * (1 - exp(-A * (x - 1))) * cos(pi_t) + v1_t = -convert(RealT, pi) * log(x + 2) * (1 - exp(-A * (x - 1))) * sin(pi_t) + v1_x = (A * log(x + 2) * exp(-A * (x - 1)) + + (1 - exp(-A * (x - 1))) / (x + 2)) * cos(pi_t) + v1_xx = ((2 * A * exp(-A * (x - 1)) / (x + 2) - + A * A * log(x + 2) * exp(-A * (x - 1)) - + (1 - exp(-A * (x - 1))) / ((x + 2) * (x + 2))) * cos(pi_t)) p = rho * rho - p_t = 2.0 * rho * rho_t - p_x = 2.0 * rho * rho_x - p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 - E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + E = p * inv_gamma_minus_one + 0.5f0 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5f0 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5f0 * rho_x * v1^2 + rho * v1 * v1_x # Some convenience constants T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_rho_cubed = 1.0 / (rho^3) + inv_rho_cubed = 1 / (rho^3) # compute the source terms # density equation @@ -98,7 +100,7 @@ end # y-momentum equation du2 = (rho_t * v1 + rho * v1_t - + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x - + + p_x + rho_x * v1^2 + 2 * rho * v1 * v1_x - # stress tensor from y-direction v1_xx * mu_) @@ -109,8 +111,8 @@ end v1_x * v1_x * mu_ - T_const * inv_rho_cubed * (p_xx * rho * rho - - 2.0 * p_x * rho * rho_x + - 2.0 * p * rho_x * rho_x - + 2 * p_x * rho * rho_x + + 2 * p * rho_x * rho_x - p * rho * rho_xx) * mu_) return SVector(du1, du2, du3) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index af0da5d1768..018ce0999eb 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -17,19 +17,20 @@ equations = ShallowWaterEquations1D(gravity_constant = 9.81) # ensure that the discontinuities lie on an element interface. function initial_condition_ec_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D) # Set the background values - H = 4.25 - v = 0.0 + RealT = eltype(x) + H = 4.25f0 + v = zero(RealT) b = sin(x[1]) # arbitrary continuous function # Setup the discontinuous water height and velocity - if x[1] >= 0.125 && x[1] <= 0.25 - H = 5.0 - v = 0.1882 + if x[1] >= 0.125f0 && x[1] <= 0.25f0 + H = 5.0f0 + v = convert(RealT, 0.1882) end # Setup a discontinuous bottom topography - if x[1] >= -0.25 && x[1] <= -0.125 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + if x[1] >= -0.25f0 && x[1] <= -0.125f0 + b = 2 + 0.5f0 * sinpi(2 * x[1]) end return prim2cons(SVector(H, v, b), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl index 76c04759389..0f9998b60a7 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl @@ -8,15 +8,16 @@ using Trixi equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81) function initial_condition_discontinuity(x, t, equations::ShallowWaterEquationsQuasi1D) - H = 2 + 0.1 * exp(-25 * x[1]^2) - v = 0.0 + RealT = eltype(x) + H = 2 + convert(RealT, 0.1) * exp(-25 * x[1]^2) + v = 0 if x[1] > 0 - b = 0.1 - a = 1.0 + b = convert(RealT, 0.1) + a = one(RealT) else - b = 0.0 - a = 1.1 + b = zero(RealT) + a = convert(RealT, 1.1) end return prim2cons(SVector(H, v, b, a), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl index a4f4b0189ba..4a777dc732c 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl @@ -17,11 +17,11 @@ equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81, H0 = 2.0) function initial_condition_discontinuous_well_balancedness(x, t, equations::ShallowWaterEquationsQuasi1D) H = equations.H0 - v = 0.0 + v = 0 # for a periodic domain, this choice of `b` and `a` mimic # discontinuity across the periodic boundary. - b = 0.5 * (x[1] + 1) + b = 0.5f0 * (x[1] + 1) a = 2 + x[1] return prim2cons(SVector(H, v, b, a), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl index 511f33d1101..fba76f6fa99 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl @@ -20,19 +20,19 @@ function initial_condition_stone_throw_discontinuous_bottom(x, t, H = equations.H0 # Discontinuous velocity - v = 0.0 - if x[1] >= -0.75 && x[1] <= 0.0 - v = -1.0 - elseif x[1] >= 0.0 && x[1] <= 0.75 - v = 1.0 + v = 0 + if x[1] >= -0.75f0 && x[1] <= 0 + v = -1 + elseif x[1] >= 0 && x[1] <= 0.75f0 + v = 1 end - b = (1.5 / exp(0.5 * ((x[1] - 1.0)^2)) + - 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + b = (1.5f0 / exp(0.5f0 * ((x[1] - 1)^2)) + + 0.75f0 / exp(0.5f0 * ((x[1] + 1)^2))) # Force a discontinuous bottom topography - if x[1] >= -1.5 && x[1] <= 0.0 - b = 0.5 + if x[1] >= -1.5f0 && x[1] <= 0 + b = 0.5f0 end return prim2cons(SVector(H, v, b), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index 5851530e230..0bb46129ba4 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -19,12 +19,12 @@ function initial_condition_discontinuous_well_balancedness(x, t, equations::ShallowWaterEquations1D) # Set the background values H = equations.H0 - v = 0.0 - b = 0.0 + v = 0 + b = zero(eltype(x)) # Setup a discontinuous bottom topography - if x[1] >= 0.5 && x[1] <= 0.75 - b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + if x[1] >= 0.5f0 && x[1] <= 0.75f0 + b = 2 + 0.5f0 * sinpi(2 * x[1]) end return prim2cons(SVector(H, v, b), equations) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl index 9ed02c0e378..a95b7892c4e 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl @@ -11,9 +11,9 @@ equations = ShallowWaterEquations1D(gravity_constant = 1.0, H0 = 3.0) function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations1D) # Set the background values H = equations.H0 - v = 0.0 + v = 0 - b = (1.5 / exp(0.5 * ((x[1] - 1.0)^2)) + 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + b = (1.5f0 / exp(0.5f0 * ((x[1] - 1)^2)) + 0.75f0 / exp(0.5f0 * ((x[1] + 1)^2))) return prim2cons(SVector(H, v, b), equations) end diff --git a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl index d3a17b513fc..67d3d56c7af 100644 --- a/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl +++ b/examples/tree_1d_dgsem/elixir_traffic_flow_lwr_trafficjam.jl @@ -21,7 +21,8 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # Discontinuous initial condition (Riemann Problem) leading to a shock that moves to the left. # The shock corresponds to the traffic congestion. function initial_condition_traffic_jam(x, t, equation::TrafficFlowLWREquations1D) - scalar = x[1] < 0.0 ? 0.5 : 1.0 + RealT = eltype(x) + scalar = x[1] < 0 ? RealT(0.5f0) : one(RealT) return SVector(scalar) end From a0c14d7088cc2504e765d9293e8717d07f780067 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 8 Nov 2024 17:17:37 +0100 Subject: [PATCH 3/8] Similiarize `Structured1D` to `Tree1D` (#2153) * Similarize Structured1D to Tree1D * Update src/meshes/structured_mesh.jl --------- Co-authored-by: Hendrik Ranocha --- .../elixir_traffic_flow_lwr_greenlight.jl | 4 ++-- src/meshes/structured_mesh.jl | 5 +++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl index e5badf14451..3cb4735e132 100644 --- a/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl +++ b/examples/structured_1d_dgsem/elixir_traffic_flow_lwr_greenlight.jl @@ -8,8 +8,8 @@ equations = TrafficFlowLWREquations1D() solver = DGSEM(polydeg = 3, surface_flux = FluxHLL(min_max_speed_davis)) -coordinates_min = (-1.0,) # minimum coordinate -coordinates_max = (1.0,) # maximum coordinate +coordinates_min = -1.0 # minimum coordinate +coordinates_max = 1.0 # maximum coordinate cells_per_dimension = (64,) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 0e4d9a359ca..ef2ed211f58 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -171,6 +171,11 @@ end function coordinates2mapping(coordinates_min::NTuple{1}, coordinates_max::NTuple{1}) mapping(xi) = linear_interpolate(xi, coordinates_min[1], coordinates_max[1]) end +# Convenience function for 1D: Do not insist on tuples +function coordinates2mapping(coordinates_min::RealT, + coordinates_max::RealT) where {RealT <: Real} + mapping(xi) = linear_interpolate(xi, coordinates_min, coordinates_max) +end function coordinates2mapping(coordinates_min::NTuple{2}, coordinates_max::NTuple{2}) function mapping(xi, eta) From e1950ac92c683b054964458c3a61d18ebf627c57 Mon Sep 17 00:00:00 2001 From: Warisa Roongaraya <81345089+warisa-r@users.noreply.github.com> Date: Fri, 8 Nov 2024 18:37:26 +0100 Subject: [PATCH 4/8] Minor modification on PairedExplicitRK3 to reduce redundant computation (#2152) * remove the redundant recalculation of c in the objective function * put verbose and max_iter on the same line as c --------- Co-authored-by: Hendrik Ranocha Co-authored-by: Daniel Doehring --- ext/TrixiNLsolveExt.jl | 12 +++++------- .../paired_explicit_runge_kutta/methods_PERK3.jl | 2 +- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/ext/TrixiNLsolveExt.jl b/ext/TrixiNLsolveExt.jl index 9247687fdc7..fa188d04c71 100644 --- a/ext/TrixiNLsolveExt.jl +++ b/ext/TrixiNLsolveExt.jl @@ -15,7 +15,7 @@ end using StableRNGs: StableRNG, rand # Use functions that are to be extended and additional symbols that are not exported -using Trixi: Trixi, compute_c_coeffs, @muladd +using Trixi: Trixi, @muladd # By default, Julia/LLVM does not use fused multiply-add operations (FMAs). # Since these FMAs can increase the performance of many numerical algorithms, @@ -29,9 +29,8 @@ using Trixi: Trixi, compute_c_coeffs, @muladd function PairedExplicitRK3_butcher_tableau_objective_function!(c_eq, a_unknown, num_stages, num_stage_evals, - monomial_coeffs, - cS2) - c_ts = compute_c_coeffs(num_stages, cS2) # ts = timestep + monomial_coeffs, c) + c_ts = c # ts = timestep # For explicit methods, a_{1,1} = 0 and a_{2,1} = c_2 (Butcher's condition) a_coeff = [0, c_ts[2], a_unknown...] # Equality constraint array that ensures that the stability polynomial computed from @@ -74,8 +73,7 @@ end # For details, see Proposition 3.2, Equation (3.3) from # Hairer, Wanner: Solving Ordinary Differential Equations 2 function Trixi.solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, monomial_coeffs, - c_s2, c; - verbose, max_iter = 100000) + c; verbose, max_iter = 100000) # Define the objective_function function objective_function!(c_eq, x) @@ -83,7 +81,7 @@ function Trixi.solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, monomial_c num_stages, num_stages, monomial_coeffs, - c_s2) + c) end # RealT is determined as the type of the first element in monomial_coeffs to ensure type consistency diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 384f4b408ca..477748b28c5 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -55,7 +55,7 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, # Butcher array abscissae c to find Butcher matrix A # This function is extended in TrixiNLsolveExt.jl a_unknown = solve_a_butcher_coeffs_unknown!(a_unknown, num_stages, - monomial_coeffs, cS2, c; + monomial_coeffs, c; verbose) end # Fill A-matrix in P-ERK style From 7f6a86018f5e1bcdfd3138769487d55c7e4c05db Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 11 Nov 2024 10:39:25 +0100 Subject: [PATCH 5/8] Make T8CodeMesh fully type general (#2154) --- src/meshes/t8code_mesh.jl | 9 +++++---- src/solvers/dgsem_t8code/containers_2d.jl | 6 +++--- src/solvers/dgsem_t8code/containers_3d.jl | 6 +++--- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index a7d7b0313d6..0d16101257f 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -30,11 +30,12 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: function T8codeMesh{NDIMS}(forest::Ptr{t8_forest}, tree_node_coordinates, nodes, boundary_names, - current_filename) where {NDIMS} + current_filename, + RealT = Float64) where {NDIMS} is_parallel = mpi_isparallel() ? True() : False() - mesh = new{NDIMS, Float64, typeof(is_parallel), NDIMS + 2, length(nodes)}(forest, - is_parallel) + mesh = new{NDIMS, RealT, typeof(is_parallel), NDIMS + 2, length(nodes)}(forest, + is_parallel) mesh.nodes = nodes mesh.boundary_names = boundary_names @@ -272,7 +273,7 @@ function T8codeMesh{NDIMS, RealT}(forest::Ptr{t8_forest}, boundary_names; polyde ntuple(_ -> length(nodes), NDIMS)..., number_of_trees) - reference_coordinates = Vector{Float64}(undef, 3) + reference_coordinates = Vector{RealT}(undef, 3) # Calculate node coordinates of reference mesh. if NDIMS == 2 diff --git a/src/solvers/dgsem_t8code/containers_2d.jl b/src/solvers/dgsem_t8code/containers_2d.jl index ce525bfdf65..4c136441386 100644 --- a/src/solvers/dgsem_t8code/containers_2d.jl +++ b/src/solvers/dgsem_t8code/containers_2d.jl @@ -7,8 +7,8 @@ # Interpolate tree_node_coordinates to each quadrant at the specified nodes. function calc_node_coordinates!(node_coordinates, - mesh::T8codeMesh{2}, - nodes::AbstractVector) + mesh::T8codeMesh{2, RealT}, + nodes::AbstractVector) where {RealT <: Real} # We use `StrideArray`s here since these buffers are used in performance-critical # places and the additional information passed to the compiler makes them faster # than native `Array`s. @@ -37,7 +37,7 @@ function calc_node_coordinates!(node_coordinates, # "integer" length to a float in relation to the unit interval [0,1]. element_length = t8_quad_len(element_level) / t8_quad_root_len - element_coords = Array{Float64}(undef, 3) + element_coords = Array{RealT}(undef, 3) t8_element_vertex_reference_coords(eclass_scheme, element, 0, pointer(element_coords)) diff --git a/src/solvers/dgsem_t8code/containers_3d.jl b/src/solvers/dgsem_t8code/containers_3d.jl index 1375782631a..981dabc397f 100644 --- a/src/solvers/dgsem_t8code/containers_3d.jl +++ b/src/solvers/dgsem_t8code/containers_3d.jl @@ -7,8 +7,8 @@ # Interpolate tree_node_coordinates to each quadrant at the specified nodes function calc_node_coordinates!(node_coordinates, - mesh::T8codeMesh{3}, - nodes::AbstractVector) + mesh::T8codeMesh{3, RealT}, + nodes::AbstractVector) where {RealT <: Real} # We use `StrideArray`s here since these buffers are used in performance-critical # places and the additional information passed to the compiler makes them faster # than native `Array`s. @@ -39,7 +39,7 @@ function calc_node_coordinates!(node_coordinates, # "integer" length to a float in relation to the unit interval [0,1]. element_length = t8_hex_len(element_level) / t8_hex_root_len - element_coords = Vector{Float64}(undef, 3) + element_coords = Vector{RealT}(undef, 3) t8_element_vertex_reference_coords(eclass_scheme, element, 0, pointer(element_coords)) From b9b6f69d8a29579484437e447c871ebf3e790151 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 12 Nov 2024 23:00:56 +0100 Subject: [PATCH 6/8] Restart `HyperbolicParabolic` with `SplitODEProblem` (#2156) * Restart HyperbolicParabolic with SplitODEProblem * cons docstrings * remove split_form flag * cons format * revert doc * docs --- .../elixir_advection_diffusion.jl | 6 +++- .../elixir_advection_diffusion_restart.jl | 30 +++++++++++++++++ src/semidiscretization/semidiscretization.jl | 3 +- ...semidiscretization_hyperbolic_parabolic.jl | 33 +++++++++++++++++++ test/test_parabolic_1d.jl | 15 +++++++++ 5 files changed, 85 insertions(+), 2 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_advection_diffusion_restart.jl diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl index 9182fd22d0e..c557065308f 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion.jl @@ -75,8 +75,12 @@ analysis_callback = AnalysisCallback(semi, interval = 100) # The AliveCallback prints short status information in regular intervals alive_callback = AliveCallback(analysis_interval = 100) +# The SaveRestartCallback allows to save a file from which a Trixi.jl simulation can be restarted +save_restart = SaveRestartCallback(interval = 100, + save_final_restart = true) + # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_restart) ############################################################################### # run the simulation diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_restart.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_restart.jl new file mode 100644 index 00000000000..2b8ad1d5fc3 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_restart.jl @@ -0,0 +1,30 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# create a restart file + +elixir_file = "elixir_advection_diffusion.jl" +trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file)) + +############################################################################### +# initialize the ODE + +restart_file = "restart_000000018.h5" +restart_filename = joinpath("out", restart_file) +tspan = (load_time(restart_filename), 2.0) + +ode = semidiscretize(semi, tspan, restart_filename); + +# Do not save restart files here +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, KenCarp4(autodiff = false), abstol = time_abs_tol, reltol = time_int_tol, + save_everystep = false, callback = callbacks) + +# Print the timer summary +summary_callback() diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index 5f886b350cd..dd5c3c4791d 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -102,7 +102,8 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan; end """ - semidiscretize(semi::AbstractSemidiscretization, tspan, restart_file::AbstractString) + semidiscretize(semi::AbstractSemidiscretization, tspan, + restart_file::AbstractString) Wrap the semidiscretization `semi` as an ODE problem in the time interval `tspan` that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 87204843ed9..1af65995e79 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -295,6 +295,39 @@ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan; return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) end +""" + semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan, + restart_file::AbstractString) + +Wrap the semidiscretization `semi` as a split ODE problem in the time interval `tspan` +that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). +The parabolic right-hand side is the first function of the split ODE problem and +will be used by default by the implicit part of IMEX methods from the +SciML ecosystem. + +The initial condition etc. is taken from the `restart_file`. +""" +function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan, + restart_file::AbstractString; + reset_threads = true) + # Optionally reset Polyester.jl threads. See + # https://github.com/trixi-framework/Trixi.jl/issues/1583 + # https://github.com/JuliaSIMD/Polyester.jl/issues/30 + if reset_threads + Polyester.reset_threads!() + end + + u0_ode = load_restart_file(semi, restart_file) + # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using + # mpi_isparallel() && MPI.Barrier(mpi_comm()) + # See https://github.com/trixi-framework/Trixi.jl/issues/328 + iip = true # is-inplace, i.e., we modify a vector when calling rhs_parabolic!, rhs! + # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the + # first function implicitly and the second one explicitly. Thus, we pass the + # stiffer parabolic function first. + return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) +end + function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) @unpack mesh, equations, boundary_conditions, source_terms, solver, cache = semi diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 062e6363a2f..48b2b9ce4d6 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -28,6 +28,21 @@ isdir(outdir) && rm(outdir, recursive = true) end end +@trixi_testset "TreeMesh1D: elixir_advection_diffusion_restart.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", + "elixir_advection_diffusion_restart.jl"), + l2=[1.0671615777620987e-5], + linf=[3.861509422325993e-5]) + # 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 "TreeMesh1D: elixir_advection_diffusion.jl (AMR)" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_advection_diffusion.jl"), From 3d465cc26832bfabfcf6e6b65f8e686cec6f7636 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 13 Nov 2024 08:45:10 +0100 Subject: [PATCH 7/8] set version to v0.9.5 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 828d1139932..032f511c22c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.9.5-DEV" +version = "0.9.5" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" From c9f07078972b82e92ca4c0c9569cbcd4b70cdb05 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 13 Nov 2024 08:45:24 +0100 Subject: [PATCH 8/8] set development version to v0.9.6-DEV --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 032f511c22c..6266fe97bb5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.9.5" +version = "0.9.6-DEV" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"