diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b0a2c93db3..4790f93d91 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -69,6 +69,7 @@ jobs: - structured - p4est_part1 - p4est_part2 + - t8code_part1 - unstructured_dgmulti - parabolic - paper_self_gravitating_gas_dynamics diff --git a/Project.toml b/Project.toml index 94c47a35ac..db41031785 100644 --- a/Project.toml +++ b/Project.toml @@ -37,6 +37,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240" +T8code = "d0cc0030-9a40-4274-8435-baadcfd54fa1" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344" TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3" @@ -80,6 +81,7 @@ StaticArrays = "1" StrideArrays = "0.1.18" StructArrays = "0.6" SummationByPartsOperators = "0.5.41" +T8code = "0.4.1" TimerOutputs = "0.5" Triangulate = "2.0" TriplotBase = "0.1" diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl new file mode 100644 index 0000000000..653bab41e2 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -0,0 +1,143 @@ +using OrdinaryDiffEq +using Trixi + +# Define new structs inside a module to allow re-evaluating the file. +module TrixiExtension + +using Trixi + +struct IndicatorSolutionIndependent{Cache <: NamedTuple} <: Trixi.AbstractIndicator + cache::Cache +end + +function IndicatorSolutionIndependent(semi) + basis = semi.solver.basis + alpha = Vector{real(basis)}() + cache = (; semi.mesh, alpha) + return IndicatorSolutionIndependent{typeof(cache)}(cache) +end + +function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4}, + mesh, equations, dg, cache; + t, kwargs...) + mesh = indicator.cache.mesh + alpha = indicator.cache.alpha + resize!(alpha, nelements(dg, cache)) + + # Predict the theoretical center. + advection_velocity = (0.2, -0.7) + center = t .* advection_velocity + + inner_distance = 1 + outer_distance = 1.85 + + # Iterate over all elements. + for element in 1:length(alpha) + # Calculate periodic distance between cell and center. + # This requires an uncurved mesh! + coordinates = SVector(0.5 * (cache.elements.node_coordinates[1, 1, 1, element] + + cache.elements.node_coordinates[1, end, 1, element]), + 0.5 * (cache.elements.node_coordinates[2, 1, 1, element] + + cache.elements.node_coordinates[2, 1, end, element])) + + # The geometric shape of the amr should be preserved when the base_level is increased. + # This is done by looking at the original coordinates of each cell. + cell_coordinates = original_coordinates(coordinates, 5 / 8) + cell_distance = periodic_distance_2d(cell_coordinates, center, 10) + if cell_distance < (inner_distance + outer_distance) / 2 + cell_coordinates = original_coordinates(coordinates, 5 / 16) + cell_distance = periodic_distance_2d(cell_coordinates, center, 10) + end + + # Set alpha according to cells position inside the circles. + target_level = (cell_distance < inner_distance) + (cell_distance < outer_distance) + alpha[element] = target_level / 2 + end + return alpha +end + +# For periodic domains, distance between two points must take into account +# periodic extensions of the domain. +function periodic_distance_2d(coordinates, center, domain_length) + dx = coordinates .- center + dx_shifted = abs.(dx .% domain_length) + dx_periodic = min.(dx_shifted, domain_length .- dx_shifted) + return sqrt(sum(dx_periodic .^ 2)) +end + +# This takes a cells coordinates and transforms them into the coordinates of a +# parent-cell it originally refined from. It does it so that the parent-cell +# has given cell_length. +function original_coordinates(coordinates, cell_length) + offset = coordinates .% cell_length + offset_sign = sign.(offset) + border = coordinates - offset + center = border + (offset_sign .* cell_length / 2) + return center +end + +end # module TrixiExtension + +import .TrixiExtension + +############################################################################### +# Semidiscretization of the linear advection equation. + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +initial_condition = initial_condition_gauss + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-5.0, -5.0) +coordinates_max = (5.0, 5.0) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +trees_per_dimension = (1, 1) + +mesh = T8codeMesh(trees_per_dimension, polydeg = 3, + mapping = mapping, + initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_controller = ControllerThreeLevel(semi, + TrixiExtension.IndicatorSolutionIndependent(semi), + base_level = 4, + med_level = 5, med_threshold = 0.1, + max_level = 6, max_threshold = 0.6) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 1.6) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + amr_callback, stepsize_callback); + +############################################################################### +# Run the simulation. + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # 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/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl new file mode 100644 index 0000000000..adf1d009a5 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -0,0 +1,87 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the linear advection equation. + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +initial_condition = initial_condition_gauss + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +# Deformed rectangle that looks like a waving flag, lower and upper faces are +# sinus curves, left and right are vertical lines. +f1(s) = SVector(-5.0, 5 * s - 5.0) +f2(s) = SVector(5.0, 5 * s + 5.0) +f3(s) = SVector(5 * s, -5.0 + 5 * sin(0.5 * pi * s)) +f4(s) = SVector(5 * s, 5.0 + 5 * sin(0.5 * pi * s)) +faces = (f1, f2, f3, f4) + +# This creates a mapping that transforms [-1, 1]^2 to the domain with the faces +# defined above. It generally doesn't work for meshes loaded from mesh files +# because these can be meshes of arbitrary domains, but the mesh below is +# specifically built on the domain [-1, 1]^2. +Trixi.validate_faces(faces) +mapping_flag = Trixi.transfinite_mapping(faces) + +# Unstructured mesh with 24 cells of the square domain [-1, 1]^n +mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connecvity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(2)) + +mesh = T8codeMesh{2}(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 0.1, + max_level = 3, max_threshold = 0.6) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.7) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + amr_callback, stepsize_callback) + +############################################################################### +# Run the simulation. + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1, # 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/examples/t8code_2d_dgsem/elixir_advection_basic.jl b/examples/t8code_2d_dgsem/elixir_advection_basic.jl new file mode 100644 index 0000000000..efc5122658 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_advection_basic.jl @@ -0,0 +1,59 @@ +# The same setup as tree_2d_dgsem/elixir_advection_basic.jl +# to verify the StructuredMesh implementation against TreeMesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +# 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 = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +trees_per_dimension = (8, 8) + +mesh = T8codeMesh(trees_per_dimension, polydeg = 3, + mapping = mapping, + initial_refinement_level = 1) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0, 1.0)); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl new file mode 100644 index 0000000000..31a8bc9369 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl @@ -0,0 +1,109 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +# Create DG solver with polynomial degree = 4 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs) + +# Deformed rectangle that looks like a waving flag, +# lower and upper faces are sinus curves, left and right are vertical lines. +f1(s) = SVector(-1.0, s - 1.0) +f2(s) = SVector(1.0, s + 1.0) +f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) +f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) + +faces = (f1, f2, f3, f4) +mapping = Trixi.transfinite_mapping(faces) + +# Create P4estMesh with 3 x 2 trees and 6 x 4 elements, +# approximate the geometry with a smaller polydeg for testing. +trees_per_dimension = (3, 2) +mesh = T8codeMesh(trees_per_dimension, polydeg = 3, + mapping = mapping, + initial_refinement_level = 1) + +function adapt_callback(forest, + forest_from, + which_tree, + lelement_id, + ts, + is_family, + num_elements, + elements_ptr)::Cint + vertex = Vector{Cdouble}(undef, 3) + + elements = unsafe_wrap(Array, elements_ptr, num_elements) + + Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) + + level = Trixi.t8_element_level(ts, elements[1]) + + # TODO: Make this condition more general. + if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 4 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.@T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest)!=0); + +# Init new forest. +new_forest_ref = Ref{Trixi.t8_forest_t}() +Trixi.t8_forest_init(new_forest_ref); +new_forest = new_forest_ref[] + +# Check out `examples/t8_step4_partition_balance_ghost.jl` in +# https://github.com/DLR-AMR/T8code.jl for detailed explanations. +let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0 + Trixi.t8_forest_set_user_data(new_forest, C_NULL) + Trixi.t8_forest_set_adapt(new_forest, mesh.forest, + Trixi.@t8_adapt_callback(adapt_callback), recursive) + Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition) + Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening) + Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES) + Trixi.t8_forest_commit(new_forest) +end + +mesh.forest = new_forest + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 0.2 +ode = semidiscretize(semi, (0.0, 0.2)); + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl new file mode 100644 index 0000000000..df9cbc26f6 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl @@ -0,0 +1,81 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the linear advection equation. + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +# 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) + +# Deformed rectangle that looks like a waving flag, +# lower and upper faces are sinus curves, left and right are vertical lines. +f1(s) = SVector(-1.0, s - 1.0) +f2(s) = SVector(1.0, s + 1.0) +f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) +f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) +faces = (f1, f2, f3, f4) + +Trixi.validate_faces(faces) +mapping_flag = Trixi.transfinite_mapping(faces) + +# Unstructured mesh with 24 cells of the square domain [-1, 1]^n. +mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connecvity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(2)) + +mesh = T8codeMesh{2}(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 2) + +# A semidiscretization collects data structures and functions for the spatial discretization. +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 0.2. +tspan = (0.0, 0.2) +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of +# the simulation setup and resets the timers. +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and +# prints the results. +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each +# time step. +stepsize_callback = StepsizeCallback(cfl = 1.4) + +# Create a CallbackSet to collect all callbacks such that they can be passed to +# the ODE solver. +callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback) + +############################################################################### +# Run the simulation. + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks. +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # Solve needs some value here but it will be overwritten by the stepsize_callback. + save_everystep = false, callback = callbacks); + +# Print the timer summary. +summary_callback() diff --git a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl new file mode 100644 index 0000000000..01e0449c67 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl @@ -0,0 +1,122 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the compressible Euler equations. + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_constant + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D +function mapping(xi_, eta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + + y = eta + 3 / 8 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3)) + + x = xi + 3 / 8 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3)) + + return SVector(x, y) +end + +############################################################################### +# Get the uncurved mesh from a file (downloads the file if not available locally) + +# Unstructured mesh with 48 cells of the square domain [-1, 1]^n +mesh_file = joinpath(@__DIR__, "square_unstructured_1.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connecvity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(2)) + +mesh = T8codeMesh{2}(conn, polydeg = 3, + mapping = mapping, + initial_refinement_level = 1) + +function adapt_callback(forest, + forest_from, + which_tree, + lelement_id, + ts, + is_family, + num_elements, + elements_ptr)::Cint + vertex = Vector{Cdouble}(undef, 3) + + elements = unsafe_wrap(Array, elements_ptr, num_elements) + + Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) + + level = Trixi.t8_element_level(ts, elements[1]) + + # TODO: Make this condition more general. + if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 3 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.@T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest)!=0); + +# Init new forest. +new_forest_ref = Ref{Trixi.t8_forest_t}() +Trixi.t8_forest_init(new_forest_ref); +new_forest = new_forest_ref[] + +# Check out `examples/t8_step4_partition_balance_ghost.jl` in +# https://github.com/DLR-AMR/T8code.jl for detailed explanations. +let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0 + Trixi.t8_forest_set_user_data(new_forest, C_NULL) + Trixi.t8_forest_set_adapt(new_forest, mesh.forest, + Trixi.@t8_adapt_callback(adapt_callback), recursive) + Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition) + Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening) + Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES) + Trixi.t8_forest_commit(new_forest) +end + +mesh.forest = new_forest + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition))) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 2.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # 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/examples/t8code_2d_dgsem/elixir_euler_sedov.jl b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl new file mode 100644 index 0000000000..965d794f8d --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl @@ -0,0 +1,97 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the compressible Euler equations. + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +The Sedov blast wave setup based on Flash +- http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 +""" +function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + # Setup based on http://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel/node184.html#SECTION010114000000000000000 + r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) + p0_outer = 1.0e-5 # = true Sedov setup + + # Calculate primitive variables + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_sedov_blast_wave + +# Get the DG approximation space +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +############################################################################### + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +trees_per_dimension = (4, 4) + +mesh = T8codeMesh(trees_per_dimension, polydeg = 4, + mapping = mapping, + initial_refinement_level = 2, periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 300 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # 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/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl new file mode 100644 index 0000000000..55a9063a00 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl @@ -0,0 +1,68 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the compressible Euler equations. + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_weak_blast_wave + +surface_flux = flux_ranocha +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +############################################################################### + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +trees_per_dimension = (4, 4) + +mesh = T8codeMesh(trees_per_dimension, polydeg = 4, + mapping = mapping, + initial_refinement_level = 2, periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # 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/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl new file mode 100644 index 0000000000..21f26d79ba --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl @@ -0,0 +1,122 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_convergence_test + +source_terms = source_terms_convergence_test + +# BCs must be passed as Dict +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +# Deformed rectangle that looks like a waving flag, +# lower and upper faces are sinus curves, left and right are vertical lines. +f1(s) = SVector(-1.0, s - 1.0) +f2(s) = SVector(1.0, s + 1.0) +f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) +f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) +faces = (f1, f2, f3, f4) + +Trixi.validate_faces(faces) +mapping_flag = Trixi.transfinite_mapping(faces) + +# Get the uncurved mesh from a file (downloads the file if not available locally) +# Unstructured mesh with 24 cells of the square domain [-1, 1]^n +mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connecvity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(2)) + +mesh = T8codeMesh{2}(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 1) + +function adapt_callback(forest, + forest_from, + which_tree, + lelement_id, + ts, + is_family, + num_elements, + elements_ptr)::Cint + vertex = Vector{Cdouble}(undef, 3) + + elements = unsafe_wrap(Array, elements_ptr, num_elements) + + Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) + + level = Trixi.t8_element_level(ts, elements[1]) + + # TODO: Make this condition more general. + if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 2 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +@assert(Trixi.t8_forest_is_committed(mesh.forest)!=0); + +# Init new forest. +new_forest_ref = Ref{Trixi.t8_forest_t}() +Trixi.t8_forest_init(new_forest_ref); +new_forest = new_forest_ref[] + +# Check out `examples/t8_step4_partition_balance_ghost.jl` in +# https://github.com/DLR-AMR/T8code.jl for detailed explanations. +let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0 + Trixi.t8_forest_set_user_data(new_forest, C_NULL) + Trixi.t8_forest_set_adapt(new_forest, mesh.forest, + Trixi.@t8_adapt_callback(adapt_callback), recursive) + Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition) + Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening) + Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES) + Trixi.t8_forest_commit(new_forest) +end + +mesh.forest = new_forest + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # 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/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl new file mode 100644 index 0000000000..32649eacff --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl @@ -0,0 +1,77 @@ +using OrdinaryDiffEq +using Trixi + +initial_condition = initial_condition_eoc_test_coupled_euler_gravity + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 2.0 +equations_euler = CompressibleEulerEquations2D(gamma) + +polydeg = 3 +solver_euler = DGSEM(polydeg, flux_hll) + +coordinates_min = (0.0, 0.0) +coordinates_max = (2.0, 2.0) + +trees_per_dimension = (1, 1) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +mesh = T8codeMesh(trees_per_dimension, polydeg = 1, + mapping = mapping, + initial_refinement_level = 2) + +semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition, solver_euler, + source_terms=source_terms_eoc_test_coupled_euler_gravity) + + +############################################################################### +# semidiscretization of the hyperbolic diffusion equations +equations_gravity = HyperbolicDiffusionEquations2D() + +solver_gravity = DGSEM(polydeg, flux_lax_friedrichs) + +semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_condition, solver_gravity, + source_terms=source_terms_harmonic) + + +############################################################################### +# combining both semidiscretizations for Euler + self-gravity +parameters = ParametersEulerGravity(background_density=2.0, # aka rho0 + # rho0 is (ab)used to add a "+8π" term to the source terms + # for the manufactured solution + gravitational_constant=1.0, # aka G + cfl=1.1, + resid_tol=1.0e-10, + n_iterations_max=1000, + timestep_gravity=timestep_gravity_erk52_3Sstar!) + +semi = SemidiscretizationEulerGravity(semi_euler, semi_gravity, parameters) + + +############################################################################### +# ODE solvers, callbacks etc. +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() + +stepsize_callback = StepsizeCallback(cfl=0.8) + +analysis_interval = 100 +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +analysis_callback = AnalysisCallback(semi_euler, interval=analysis_interval, + save_analysis=true) + +callbacks = CallbackSet(summary_callback, stepsize_callback, + analysis_callback, alive_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # 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 +println("Number of gravity subcycles: ", semi.gravity_counter.ncalls_since_readout) diff --git a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl new file mode 100644 index 0000000000..463f916fa2 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl @@ -0,0 +1,60 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the compressible ideal GLM-MHD equations. + +gamma = 5/3 +equations = IdealGlmMhdEquations2D(gamma) + +initial_condition = initial_condition_convergence_test + +# Get the DG approximation space +volume_flux = (flux_central, flux_nonconservative_powell) +solver = DGSEM(polydeg=4, surface_flux=(flux_hll, flux_nonconservative_powell), + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0 , 0.0 ) +coordinates_max = (sqrt(2.0), sqrt(2.0)) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +trees_per_dimension = (8, 8) + +mesh = T8codeMesh(trees_per_dimension, polydeg=3, + mapping=mapping, + initial_refinement_level=0, periodicity=true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +cfl = 0.9 +stepsize_callback = StepsizeCallback(cfl=cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1.0, # 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/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl new file mode 100644 index 0000000000..9a4bd99e44 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl @@ -0,0 +1,134 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible ideal GLM-MHD equations +equations = IdealGlmMhdEquations2D(1.4) + +""" + initial_condition_rotor(x, t, equations::IdealGlmMhdEquations2D) + +The classical MHD rotor test case. Here, the setup is taken from +- Dominik Derigs, Gregor J. Gassner, Stefanie Walch & Andrew R. Winters (2018) + Entropy Stable Finite Volume Approximations for Ideal Magnetohydrodynamics + [doi: 10.1365/s13291-018-0178-9](https://doi.org/10.1365/s13291-018-0178-9) +""" +function initial_condition_rotor(x, t, equations::IdealGlmMhdEquations2D) + # setup taken from Derigs et al. DMV article (2018) + # domain must be [0, 1] x [0, 1], γ = 1.4 + dx = x[1] - 0.5 + dy = x[2] - 0.5 + r = sqrt(dx^2 + dy^2) + f = (0.115 - r) / 0.015 + if r <= 0.1 + rho = 10.0 + v1 = -20.0 * dy + v2 = 20.0 * dx + elseif r >= 0.115 + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + else + rho = 1.0 + 9.0 * f + v1 = -20.0 * f * dy + v2 = 20.0 * f * dx + end + v3 = 0.0 + p = 1.0 + B1 = 5.0 / sqrt(4.0 * pi) + B2 = 0.0 + B3 = 0.0 + psi = 0.0 + return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) +end +initial_condition = initial_condition_rotor + +surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell) +volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +# Affine type mapping to take the [-1,1]^2 domain from the mesh file +# and put it onto the rotor domain [0,1]^2 and then warp it with a mapping +# as described in https://arxiv.org/abs/2012.12040 +function mapping_twist(xi, eta) + y = 0.5 * (eta + 1.0) + + 0.05 * cos(1.5 * pi * (2.0 * xi - 1.0)) * cos(0.5 * pi * (2.0 * eta - 1.0)) + x = 0.5 * (xi + 1.0) + 0.05 * cos(0.5 * pi * (2.0 * xi - 1.0)) * cos(2.0 * pi * y) + return SVector(x, y) +end + +mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connecvity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(2)) + +mesh = T8codeMesh{2}(conn, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 1) + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.15) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, + variable = density_pressure) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 1, + med_level = 3, med_threshold = 0.05, + max_level = 5, max_threshold = 0.1) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +cfl = 0.5 +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # 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/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl new file mode 100644 index 0000000000..c19f440ebc --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -0,0 +1,60 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the shallow water equations. + +equations = ShallowWaterEquations2D(gravity_constant=9.81) + +initial_condition = initial_condition_convergence_test # MMS EOC test + + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg=3, surface_flux=(flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the P4estMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y)) +coordinates_max = (sqrt(2.0), sqrt(2.0)) # maximum coordinates (max(x), max(y)) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +trees_per_dimension = (8, 8) + +mesh = T8codeMesh(trees_per_dimension, polydeg=3, + mapping=mapping, + initial_refinement_level=1) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms=source_terms_convergence_test) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-8, reltol=1.0e-8, + ode_default_options()..., callback=callbacks); +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index b0c872b190..990c33f3c9 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -54,6 +54,7 @@ using Octavian: Octavian, matmul! using Polyester: Polyester, @batch # You know, the cheapest threads you can find... using OffsetArrays: OffsetArray, OffsetVector using P4est +using T8code using Setfield: @set using RecipesBase: RecipesBase using Requires: @require @@ -110,6 +111,7 @@ include("basic_types.jl") include("auxiliary/auxiliary.jl") include("auxiliary/mpi.jl") include("auxiliary/p4est.jl") +include("auxiliary/t8code.jl") include("equations/equations.jl") include("meshes/meshes.jl") include("solvers/solvers.jl") @@ -210,7 +212,7 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, export lake_at_rest_error export ncomponents, eachcomponent -export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh +export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh export DG, DGSEM, LobattoLegendreBasis, @@ -277,6 +279,7 @@ function __init__() init_mpi() init_p4est() + init_t8code() register_error_hints() diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl new file mode 100644 index 0000000000..37cb782bb9 --- /dev/null +++ b/src/auxiliary/t8code.jl @@ -0,0 +1,486 @@ +""" + init_t8code() + +Initialize `t8code` by calling `sc_init`, `p4est_init`, and `t8_init` while +setting the log level to `SC_LP_ERROR`. This function will check if `t8code` +is already initialized and if yes, do nothing, thus it is safe to call it +multiple times. +""" +function init_t8code() + t8code_package_id = t8_get_package_id() + if t8code_package_id >= 0 + return nothing + end + + # Initialize the sc library, has to happen before we initialize t8code. + let catch_signals = 0, print_backtrace = 0, log_handler = C_NULL + T8code.Libt8.sc_init(mpi_comm(), catch_signals, print_backtrace, log_handler, + T8code.Libt8.SC_LP_ERROR) + end + + if T8code.Libt8.p4est_is_initialized() == 0 + # Initialize `p4est` with log level ERROR to prevent a lot of output in AMR simulations + T8code.Libt8.p4est_init(C_NULL, T8code.Libt8.SC_LP_ERROR) + end + + # Initialize t8code with log level ERROR to prevent a lot of output in AMR simulations. + t8_init(T8code.Libt8.SC_LP_ERROR) + + if haskey(ENV, "TRIXI_T8CODE_SC_FINALIZE") + # Normally, `sc_finalize` should always be called during shutdown of an + # application. It checks whether there is still un-freed memory by t8code + # and/or T8code.jl and throws an exception if this is the case. For + # production runs this is not mandatory, but is helpful during + # development. Hence, this option is only activated when environment + # variable TRIXI_T8CODE_SC_FINALIZE exists. + @warn "T8code.jl: sc_finalize will be called during shutdown of Trixi.jl." + MPI.add_finalize_hook!(T8code.Libt8.sc_finalize) + end + + return nothing +end + +function trixi_t8_unref_forest(forest) + t8_forest_unref(Ref(forest)) +end + +function t8_free(ptr) + T8code.Libt8.sc_free(t8_get_package_id(), ptr) +end + +function trixi_t8_count_interfaces(forest) + # Check that forest is a committed, that is valid and usable, forest. + @assert t8_forest_is_committed(forest) != 0 + + # Get the number of local elements of forest. + num_local_elements = t8_forest_get_local_num_elements(forest) + # Get the number of ghost elements of forest. + num_ghost_elements = t8_forest_get_num_ghosts(forest) + # Get the number of trees that have elements of this process. + num_local_trees = t8_forest_get_num_local_trees(forest) + + current_index = t8_locidx_t(0) + + local_num_conform = 0 + local_num_mortars = 0 + local_num_boundary = 0 + + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) + + # Get the number of elements of this tree. + num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) + + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(forest, itree, ielement) + + level = t8_element_level(eclass_scheme, element) + + num_faces = t8_element_num_faces(eclass_scheme, element) + + for iface in 0:(num_faces - 1) + pelement_indices_ref = Ref{Ptr{t8_locidx_t}}() + pneighbor_leafs_ref = Ref{Ptr{Ptr{t8_element}}}() + pneigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() + + dual_faces_ref = Ref{Ptr{Cint}}() + num_neighbors_ref = Ref{Cint}() + + forest_is_balanced = Cint(1) + + t8_forest_leaf_face_neighbors(forest, itree, element, + pneighbor_leafs_ref, iface, dual_faces_ref, + num_neighbors_ref, + pelement_indices_ref, pneigh_scheme_ref, + forest_is_balanced) + + num_neighbors = num_neighbors_ref[] + neighbor_ielements = unsafe_wrap(Array, pelement_indices_ref[], + num_neighbors) + neighbor_leafs = unsafe_wrap(Array, pneighbor_leafs_ref[], num_neighbors) + neighbor_scheme = pneigh_scheme_ref[] + + if num_neighbors > 0 + neighbor_level = t8_element_level(neighbor_scheme, neighbor_leafs[1]) + + # Conforming interface: The second condition ensures we only visit the interface once. + if level == neighbor_level && current_index <= neighbor_ielements[1] + local_num_conform += 1 + elseif level < neighbor_level + local_num_mortars += 1 + end + + else + local_num_boundary += 1 + end + + t8_free(dual_faces_ref[]) + t8_free(pneighbor_leafs_ref[]) + t8_free(pelement_indices_ref[]) + end # for + + current_index += 1 + end # for + end # for + + return (interfaces = local_num_conform, + mortars = local_num_mortars, + boundaries = local_num_boundary) +end + +function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundaries, + boundary_names) + # Check that forest is a committed, that is valid and usable, forest. + @assert t8_forest_is_committed(forest) != 0 + + # Get the number of local elements of forest. + num_local_elements = t8_forest_get_local_num_elements(forest) + # Get the number of ghost elements of forest. + num_ghost_elements = t8_forest_get_num_ghosts(forest) + # Get the number of trees that have elements of this process. + num_local_trees = t8_forest_get_num_local_trees(forest) + + current_index = t8_locidx_t(0) + + local_num_conform = 0 + local_num_mortars = 0 + local_num_boundary = 0 + + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) + + # Get the number of elements of this tree. + num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) + + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(forest, itree, ielement) + + level = t8_element_level(eclass_scheme, element) + + num_faces = t8_element_num_faces(eclass_scheme, element) + + for iface in 0:(num_faces - 1) + + # Compute the `orientation` of the touching faces. + if t8_element_is_root_boundary(eclass_scheme, element, iface) == 1 + cmesh = t8_forest_get_cmesh(forest) + itree_in_cmesh = t8_forest_ltreeid_to_cmesh_ltreeid(forest, itree) + iface_in_tree = t8_element_tree_face(eclass_scheme, element, iface) + orientation_ref = Ref{Cint}() + + t8_cmesh_get_face_neighbor(cmesh, itree_in_cmesh, iface_in_tree, C_NULL, + orientation_ref) + orientation = orientation_ref[] + else + orientation = zero(Cint) + end + + pelement_indices_ref = Ref{Ptr{t8_locidx_t}}() + pneighbor_leafs_ref = Ref{Ptr{Ptr{t8_element}}}() + pneigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() + + dual_faces_ref = Ref{Ptr{Cint}}() + num_neighbors_ref = Ref{Cint}() + + forest_is_balanced = Cint(1) + + t8_forest_leaf_face_neighbors(forest, itree, element, + pneighbor_leafs_ref, iface, dual_faces_ref, + num_neighbors_ref, + pelement_indices_ref, pneigh_scheme_ref, + forest_is_balanced) + + num_neighbors = num_neighbors_ref[] + dual_faces = unsafe_wrap(Array, dual_faces_ref[], num_neighbors) + neighbor_ielements = unsafe_wrap(Array, pelement_indices_ref[], + num_neighbors) + neighbor_leafs = unsafe_wrap(Array, pneighbor_leafs_ref[], num_neighbors) + neighbor_scheme = pneigh_scheme_ref[] + + if num_neighbors > 0 + neighbor_level = t8_element_level(neighbor_scheme, neighbor_leafs[1]) + + # Conforming interface: The second condition ensures we only visit the interface once. + if level == neighbor_level && current_index <= neighbor_ielements[1] + local_num_conform += 1 + + faces = (iface, dual_faces[1]) + interface_id = local_num_conform + + # Write data to interfaces container. + interfaces.neighbor_ids[1, interface_id] = current_index + 1 + interfaces.neighbor_ids[2, interface_id] = neighbor_ielements[1] + 1 + + # Iterate over primary and secondary element. + for side in 1:2 + # Align interface in positive coordinate direction of primary element. + # For orientation == 1, the secondary element needs to be indexed backwards + # relative to the interface. + if side == 1 || orientation == 0 + # Forward indexing + indexing = :i_forward + else + # Backward indexing + indexing = :i_backward + end + + if faces[side] == 0 + # Index face in negative x-direction + interfaces.node_indices[side, interface_id] = (:begin, + indexing) + elseif faces[side] == 1 + # Index face in positive x-direction + interfaces.node_indices[side, interface_id] = (:end, + indexing) + elseif faces[side] == 2 + # Index face in negative y-direction + interfaces.node_indices[side, interface_id] = (indexing, + :begin) + else # faces[side] == 3 + # Index face in positive y-direction + interfaces.node_indices[side, interface_id] = (indexing, + :end) + end + end + + # Non-conforming interface. + elseif level < neighbor_level + local_num_mortars += 1 + + faces = (dual_faces[1], iface) + + mortar_id = local_num_mortars + + # Last entry is the large element. + mortars.neighbor_ids[end, mortar_id] = current_index + 1 + + # First `1:end-1` entries are the smaller elements. + mortars.neighbor_ids[1:(end - 1), mortar_id] .= neighbor_ielements .+ + 1 + + for side in 1:2 + # Align mortar in positive coordinate direction of small side. + # For orientation == 1, the large side needs to be indexed backwards + # relative to the mortar. + if side == 1 || orientation == 0 + # Forward indexing for small side or orientation == 0. + indexing = :i_forward + else + # Backward indexing for large side with reversed orientation. + indexing = :i_backward + # Since the orientation is reversed we have to account for this + # when filling the `neighbor_ids` array. + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + + 1 + end + + if faces[side] == 0 + # Index face in negative x-direction + mortars.node_indices[side, mortar_id] = (:begin, indexing) + elseif faces[side] == 1 + # Index face in positive x-direction + mortars.node_indices[side, mortar_id] = (:end, indexing) + elseif faces[side] == 2 + # Index face in negative y-direction + mortars.node_indices[side, mortar_id] = (indexing, :begin) + else # faces[side] == 3 + # Index face in positive y-direction + mortars.node_indices[side, mortar_id] = (indexing, :end) + end + end + + # else: "level > neighbor_level" is skipped since we visit the mortar interface only once. + end + + # Domain boundary. + else + local_num_boundary += 1 + boundary_id = local_num_boundary + + boundaries.neighbor_ids[boundary_id] = current_index + 1 + + if iface == 0 + # Index face in negative x-direction. + boundaries.node_indices[boundary_id] = (:begin, :i_forward) + elseif iface == 1 + # Index face in positive x-direction. + boundaries.node_indices[boundary_id] = (:end, :i_forward) + elseif iface == 2 + # Index face in negative y-direction. + boundaries.node_indices[boundary_id] = (:i_forward, :begin) + else # iface == 3 + # Index face in positive y-direction. + boundaries.node_indices[boundary_id] = (:i_forward, :end) + end + + # One-based indexing. + boundaries.name[boundary_id] = boundary_names[iface + 1, itree + 1] + end + + t8_free(dual_faces_ref[]) + t8_free(pneighbor_leafs_ref[]) + t8_free(pelement_indices_ref[]) + end # for iface = ... + + current_index += 1 + end # for + end # for + + return (interfaces = local_num_conform, + mortars = local_num_mortars, + boundaries = local_num_boundary) +end + +function trixi_t8_get_local_element_levels(forest) + # Check that forest is a committed, that is valid and usable, forest. + @assert t8_forest_is_committed(forest) != 0 + + levels = Vector{Int}(undef, t8_forest_get_local_num_elements(forest)) + + # Get the number of trees that have elements of this process. + num_local_trees = t8_forest_get_num_local_trees(forest) + + current_index = 0 + + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) + + # Get the number of elements of this tree. + num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) + + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(forest, itree, ielement) + current_index += 1 + levels[current_index] = t8_element_level(eclass_scheme, element) + end # for + end # for + + return levels +end + +# Callback function prototype to decide for refining and coarsening. +# If `is_family` equals 1, the first `num_elements` in elements +# form a family and we decide whether this family should be coarsened +# or only the first element should be refined. +# Otherwise `is_family` must equal zero and we consider the first entry +# of the element array for refinement. +# Entries of the element array beyond the first `num_elements` are undefined. +# \param [in] forest the forest to which the new elements belong +# \param [in] forest_from the forest that is adapted. +# \param [in] which_tree the local tree containing `elements` +# \param [in] lelement_id the local element id in `forest_old` in the tree of the current element +# \param [in] ts the eclass scheme of the tree +# \param [in] is_family if 1, the first `num_elements` entries in `elements` form a family. If 0, they do not. +# \param [in] num_elements the number of entries in `elements` that are defined +# \param [in] elements Pointers to a family or, if `is_family` is zero, +# pointer to one element. +# \return greater zero if the first entry in `elements` should be refined, +# smaller zero if the family `elements` shall be coarsened, +# zero else. +function adapt_callback(forest, + forest_from, + which_tree, + lelement_id, + ts, + is_family, + num_elements, + elements)::Cint + num_levels = t8_forest_get_local_num_elements(forest_from) + + indicator_ptr = Ptr{Int}(t8_forest_get_user_data(forest)) + indicators = unsafe_wrap(Array, indicator_ptr, num_levels) + + offset = t8_forest_get_tree_element_offset(forest_from, which_tree) + + # Only allow coarsening for complete families. + if indicators[offset + lelement_id + 1] < 0 && is_family == 0 + return Cint(0) + end + + return Cint(indicators[offset + lelement_id + 1]) +end + +function trixi_t8_adapt_new(old_forest, indicators) + # Check that forest is a committed, that is valid and usable, forest. + @assert t8_forest_is_committed(old_forest) != 0 + + # Init new forest. + new_forest_ref = Ref{t8_forest_t}() + t8_forest_init(new_forest_ref) + new_forest = new_forest_ref[] + + let set_from = C_NULL, recursive = 0, set_for_coarsening = 0, no_repartition = 0 + t8_forest_set_user_data(new_forest, pointer(indicators)) + t8_forest_set_adapt(new_forest, old_forest, @t8_adapt_callback(adapt_callback), + recursive) + t8_forest_set_balance(new_forest, set_from, no_repartition) + t8_forest_set_partition(new_forest, set_from, set_for_coarsening) + t8_forest_set_ghost(new_forest, 1, T8_GHOST_FACES) # Note: MPI support not available yet so it is a dummy call. + t8_forest_commit(new_forest) + end + + return new_forest +end + +function trixi_t8_get_difference(old_levels, new_levels, num_children) + old_nelems = length(old_levels) + new_nelems = length(new_levels) + + changes = Vector{Int}(undef, old_nelems) + + # Local element indices. + old_index = 1 + new_index = 1 + + while old_index <= old_nelems && new_index <= new_nelems + if old_levels[old_index] < new_levels[new_index] + # Refined. + + changes[old_index] = 1 + + old_index += 1 + new_index += num_children + + elseif old_levels[old_index] > new_levels[new_index] + # Coarsend. + + for child_index in old_index:(old_index + num_children - 1) + changes[child_index] = -1 + end + + old_index += num_children + new_index += 1 + + else + # No changes. + + changes[old_index] = 0 + + old_index += 1 + new_index += 1 + end + end + + return changes +end + +# Coarsen or refine marked cells and rebalance forest. Return a difference between +# old and new mesh. +function trixi_t8_adapt!(mesh, indicators) + old_levels = trixi_t8_get_local_element_levels(mesh.forest) + + forest_cached = trixi_t8_adapt_new(mesh.forest, indicators) + + new_levels = trixi_t8_get_local_element_levels(forest_cached) + + differences = trixi_t8_get_difference(old_levels, new_levels, 2^ndims(mesh)) + + mesh.forest = forest_cached + + return differences +end diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index bef49b4c48..4d80e6e113 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -471,6 +471,65 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::P4estMesh, return has_changed end +function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::SerialT8codeMesh, + equations, dg::DG, cache, semi, + t, iter; + only_refine = false, only_coarsen = false, + passive_args = ()) + has_changed = false + + @unpack controller, adaptor = amr_callback + + u = wrap_array(u_ode, mesh, equations, dg, cache) + indicators = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, + cache, t = t, iter = iter) + + if only_coarsen + indicators[indicators .> 0] .= 0 + end + + if only_refine + indicators[indicators .< 0] .= 0 + end + + @boundscheck begin + @assert axes(indicators)==(Base.OneTo(ncells(mesh)),) ("Indicator array (axes = $(axes(indicators))) and mesh cells (axes = $(Base.OneTo(ncells(mesh)))) have different axes") + end + + @trixi_timeit timer() "adapt" begin + difference = @trixi_timeit timer() "mesh" trixi_t8_adapt!(mesh, indicators) + + @trixi_timeit timer() "solver" adapt!(u_ode, adaptor, mesh, equations, dg, + cache, difference) + end + + # Store whether there were any cells coarsened or refined and perform load balancing. + has_changed = any(difference .!= 0) + + # TODO: T8codeMesh for MPI not implemented yet. + # Check if mesh changed on other processes + # if mpi_isparallel() + # has_changed = MPI.Allreduce!(Ref(has_changed), |, mpi_comm())[] + # end + + if has_changed + # TODO: T8codeMesh for MPI not implemented yet. + # if mpi_isparallel() && amr_callback.dynamic_load_balancing + # @trixi_timeit timer() "dynamic load balancing" begin + # global_first_quadrant = unsafe_wrap(Array, mesh.p4est.global_first_quadrant, mpi_nranks() + 1) + # old_global_first_quadrant = copy(global_first_quadrant) + # partition!(mesh) + # rebalance_solver!(u_ode, mesh, equations, dg, cache, old_global_first_quadrant) + # end + # end + + reinitialize_boundaries!(semi.boundary_conditions, cache) + end + + # Return true if there were any cells coarsened or refined, otherwise false. + return has_changed +end + function reinitialize_boundaries!(boundary_conditions::UnstructuredSortedBoundaryTypes, cache) # Reinitialize boundary types container because boundaries may have changed. @@ -639,6 +698,10 @@ function current_element_levels(mesh::P4estMesh, solver, cache) return current_levels end +function current_element_levels(mesh::T8codeMesh, solver, cache) + return trixi_t8_get_local_element_levels(mesh.forest) +end + # TODO: Taal refactor, merge the two loops of ControllerThreeLevel and IndicatorLöhner etc.? # But that would remove the simplest possibility to write that stuff to a file... # We could of course implement some additional logic and workarounds, but is it worth the effort? diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 400d16347d..1d37dfce03 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -333,9 +333,79 @@ function coarsen_elements!(u::AbstractArray{<:Any, 4}, element_id, end end +# Coarsen and refine elements in the DG solver based on a difference list. +function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, + dg::DGSEM, cache, difference) + + # Return early if there is nothing to do. + if !any(difference .!= 0) + return nothing + end + + # Number of (local) cells/elements. + old_nelems = nelements(dg, cache) + new_nelems = ncells(mesh) + + # Local element indices. + old_index = 1 + new_index = 1 + + # Note: This is true for `quads` only. + T8_CHILDREN = 4 + + # Retain current solution data. + old_u_ode = copy(u_ode) + + GC.@preserve old_u_ode begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + while old_index <= old_nelems && new_index <= new_nelems + if difference[old_index] > 0 # Refine. + + # Refine element and store solution directly in new data structure. + refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg) + + old_index += 1 + new_index += T8_CHILDREN + + elseif difference[old_index] < 0 # Coarsen. + + # If an element is to be removed, sanity check if the following elements + # are also marked - otherwise there would be an error in the way the + # cells/elements are sorted. + @assert all(difference[old_index:(old_index + T8_CHILDREN - 1)] .< 0) "bad cell/element order" + + # Coarsen elements and store solution directly in new data structure. + coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, + dg) + + old_index += T8_CHILDREN + new_index += 1 + + else # No changes. + + # Copy old element data to new element container. + @views u[:, .., new_index] .= old_u[:, .., old_index] + + old_index += 1 + new_index += 1 + end + end # while + end # GC.@preserve old_u_ode + + return nothing +end + # this method is called when an `ControllerThreeLevel` is constructed function create_cache(::Type{ControllerThreeLevel}, - mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DG, cache) + mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, equations, + dg::DG, cache) controller_value = Vector{Int}(undef, nelements(dg, cache)) return (; controller_value) end diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 7c453aab63..fad42b1109 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -534,6 +534,36 @@ function print_amr_information(callbacks, mesh::P4estMesh, solver, cache) return nothing end +# Print level information only if AMR is enabled +function print_amr_information(callbacks, mesh::T8codeMesh, solver, cache) + + # Return early if there is nothing to print + uses_amr(callbacks) || return nothing + + # TODO: Switch to global element levels array when MPI supported or find + # another solution. + levels = trixi_t8_get_local_element_levels(mesh.forest) + + min_level = minimum(levels) + max_level = maximum(levels) + + mpi_println(" minlevel = $min_level") + mpi_println(" maxlevel = $max_level") + + if min_level > 0 + elements_per_level = [count(==(l), levels) for l in 1:max_level] + + for level in max_level:-1:(min_level + 1) + mpi_println(" ├── level $level: " * + @sprintf("% 14d", elements_per_level[level])) + end + mpi_println(" └── level $min_level: " * + @sprintf("% 14d", elements_per_level[min_level])) + end + + return nothing +end + # Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming". function analyze_integrals(analysis_integrals::NTuple{N, Any}, io, du, u, t, semi) where {N} diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 6c74e172e4..4e456f7987 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -31,7 +31,7 @@ end function create_cache_analysis(analyzer, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache, RealT, uEltype) @@ -108,7 +108,7 @@ end function calc_error_norms(func, u, t, analyzer, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, equations, + P4estMesh{2}, T8codeMesh{2}}, equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @unpack node_coordinates, inverse_jacobian = cache.elements @@ -176,7 +176,7 @@ end function integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, equations, + P4estMesh{2}, T8codeMesh{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis @@ -204,7 +204,7 @@ end function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg @@ -215,7 +215,7 @@ end function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, @@ -259,7 +259,8 @@ function analyze(::Val{:l2_divb}, du, u, t, end function analyze(::Val{:l2_divb}, du, u, t, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, + mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) @unpack contravariant_vectors = cache.elements integrate_via_indices(u, mesh, equations, dg, cache, cache, @@ -326,7 +327,8 @@ function analyze(::Val{:linf_divb}, du, u, t, end function analyze(::Val{:linf_divb}, du, u, t, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, + mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) @unpack derivative_matrix, weights = dg.basis @unpack contravariant_vectors = cache.elements diff --git a/src/callbacks_step/save_restart_dg.jl b/src/callbacks_step/save_restart_dg.jl index 5695eb8bed..8db6db2d2b 100644 --- a/src/callbacks_step/save_restart_dg.jl +++ b/src/callbacks_step/save_restart_dg.jl @@ -7,7 +7,8 @@ function save_restart_file(u, time, dt, timestep, mesh::Union{SerialTreeMesh, StructuredMesh, - UnstructuredMesh2D, SerialP4estMesh}, + UnstructuredMesh2D, SerialP4estMesh, + SerialT8codeMesh}, equations, dg::DG, cache, restart_callback) @unpack output_directory = restart_callback diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 6cd4a0ec9c..6d5004ff65 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -7,7 +7,8 @@ function save_solution_file(u, time, dt, timestep, mesh::Union{SerialTreeMesh, StructuredMesh, - UnstructuredMesh2D, SerialP4estMesh}, + UnstructuredMesh2D, SerialP4estMesh, + SerialT8codeMesh}, equations, dg::DG, cache, solution_callback, element_variables = Dict{Symbol, Any}(); system = "") diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 89a2b2b835..673c3ba6aa 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -75,7 +75,9 @@ function max_dt(u, t, mesh::ParallelTreeMesh{2}, return dt end -function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, +function max_dt(u, t, + mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -109,7 +111,9 @@ function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMe return 2 / (nnodes(dg) * max_scaled_speed) end -function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, +function max_dt(u, t, + mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, constant_speed::True, equations, dg::DG, cache) @unpack contravariant_vectors, inverse_jacobian = cache.elements diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index da67fe23e0..b9895e7d45 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -6,7 +6,7 @@ #! format: noindent # Save current mesh with some context information as an HDF5 file. -function save_mesh_file(mesh::Union{TreeMesh, P4estMesh}, output_directory, +function save_mesh_file(mesh::Union{TreeMesh, P4estMesh, T8codeMesh}, output_directory, timestep = 0) save_mesh_file(mesh, output_directory, timestep, mpi_parallel(mesh)) end @@ -220,6 +220,13 @@ function save_mesh_file(mesh::P4estMesh, output_directory, timestep, mpi_paralle return filename end +# TODO: Implement this function as soon as there is support for this in `t8code`. +function save_mesh_file(mesh::T8codeMesh, output_directory, timestep, mpi_parallel) + error("Mesh file output not supported yet for `T8codeMesh`.") + + return joinpath(output_directory, "dummy_mesh.h5") +end + """ load_mesh(restart_file::AbstractString; n_cells_max) diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index 2716aa2007..ed2158b169 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -12,6 +12,7 @@ include("unstructured_mesh.jl") include("face_interpolant.jl") include("transfinite_mappings_3d.jl") include("p4est_mesh.jl") +include("t8code_mesh.jl") include("mesh_io.jl") include("dgmulti_meshes.jl") end # @muladd diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl new file mode 100644 index 0000000000..13edcc2971 --- /dev/null +++ b/src/meshes/t8code_mesh.jl @@ -0,0 +1,345 @@ +""" + T8codeMesh{NDIMS} <: AbstractMesh{NDIMS} + +An unstructured curved mesh based on trees that uses the C library +['t8code'](https://github.com/DLR-AMR/t8code) +to manage trees and mesh refinement. +""" +mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: + AbstractMesh{NDIMS} + cmesh :: Ptr{t8_cmesh} # cpointer to coarse mesh + scheme :: Ptr{t8_eclass_scheme} # cpointer to element scheme + forest :: Ptr{t8_forest} # cpointer to forest + is_parallel :: IsParallel + + # This specifies the geometry interpolation for each tree. + tree_node_coordinates::Array{RealT, NDIMSP2} # [dimension, i, j, k, tree] + + # Stores the quadrature nodes. + nodes::SVector{NNODES, RealT} + + boundary_names :: Array{Symbol, 2} # [face direction, tree] + current_filename :: String + + ninterfaces :: Int + nmortars :: Int + nboundaries :: Int + + function T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, + boundary_names, + current_filename) where {NDIMS} + is_parallel = False() + + mesh = new{NDIMS, Float64, typeof(is_parallel), NDIMS + 2, length(nodes)}(cmesh, + scheme, + forest, + is_parallel) + + mesh.nodes = nodes + mesh.boundary_names = boundary_names + mesh.current_filename = current_filename + mesh.tree_node_coordinates = tree_node_coordinates + + finalizer(mesh) do mesh + # When finalizing `mesh.forest`, `mesh.scheme` and `mesh.cmesh` are + # also cleaned up from within `t8code`. The cleanup code for + # `cmesh` does some MPI calls for deallocating shared memory + # arrays. Due to garbage collection in Julia the order of shutdown + # is not deterministic. The following code might happen after MPI + # is already in finalized state. + # If the environment variable `TRIXI_T8CODE_SC_FINALIZE` is set the + # `finalize_hook` of the MPI module takes care of the cleanup. See + # further down. However, this might cause a pile-up of `mesh` + # objects during long-running sessions. + if !MPI.Finalized() + trixi_t8_unref_forest(mesh.forest) + end + end + + # This finalizer call is only recommended during development and not for + # production runs, especially long-running sessions since a reference to + # the `mesh` object will be kept throughout the lifetime of the session. + # See comments in `init_t8code()` in file `src/auxiliary/t8code.jl` for + # more information. + if haskey(ENV, "TRIXI_T8CODE_SC_FINALIZE") + MPI.add_finalize_hook!() do + trixi_t8_unref_forest(mesh.forest) + end + end + + return mesh + end +end + +const SerialT8codeMesh{NDIMS} = T8codeMesh{NDIMS, <:Real, <:False} +@inline mpi_parallel(mesh::SerialT8codeMesh) = False() + +@inline Base.ndims(::T8codeMesh{NDIMS}) where {NDIMS} = NDIMS +@inline Base.real(::T8codeMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT + +@inline ntrees(mesh::T8codeMesh) = Int(t8_forest_get_num_local_trees(mesh.forest)) +@inline ncells(mesh::T8codeMesh) = Int(t8_forest_get_local_num_elements(mesh.forest)) +@inline ninterfaces(mesh::T8codeMesh) = mesh.ninterfaces +@inline nmortars(mesh::T8codeMesh) = mesh.nmortars +@inline nboundaries(mesh::T8codeMesh) = mesh.nboundaries + +function Base.show(io::IO, mesh::T8codeMesh) + print(io, "T8codeMesh{", ndims(mesh), ", ", real(mesh), "}") +end + +function Base.show(io::IO, ::MIME"text/plain", mesh::T8codeMesh) + if get(io, :compact, false) + show(io, mesh) + else + setup = [ + "#trees" => ntrees(mesh), + "current #cells" => ncells(mesh), + "polydeg" => length(mesh.nodes) - 1, + ] + summary_box(io, + "T8codeMesh{" * string(ndims(mesh)) * ", " * string(real(mesh)) * "}", + setup) + end +end + +""" + T8codeMesh(trees_per_dimension; polydeg, mapping=identity, + RealT=Float64, initial_refinement_level=0, periodicity=true) + +Create a structured potentially curved 'T8codeMesh' of the specified size. + +Non-periodic boundaries will be called ':x_neg', ':x_pos', ':y_neg', ':y_pos', ':z_neg', ':z_pos'. + +# Arguments +- 'trees_per_dimension::NTupleE{NDIMS, Int}': the number of trees in each dimension. +- 'polydeg::Integer': polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. +- 'mapping': a function of 'NDIMS' variables to describe the mapping that transforms + the reference mesh ('[-1, 1]^n') to the physical domain. +- 'RealT::Type': the type that should be used for coordinates. +- 'initial_refinement_level::Integer': refine the mesh uniformly to this level before the simulation starts. +- 'periodicity': either a 'Bool' deciding if all of the boundaries are periodic or an 'NTuple{NDIMS, Bool}' + deciding for each dimension if the boundaries in this dimension are periodic. +""" +function T8codeMesh(trees_per_dimension; polydeg, + mapping = coordinates2mapping((-1.0, -1.0), (1.0, 1.0)), + RealT = Float64, initial_refinement_level = 0, periodicity = true) + NDIMS = length(trees_per_dimension) + + @assert NDIMS == 2 # Only support for NDIMS = 2 yet. + + # Convert periodicity to a Tuple of a Bool for every dimension + if all(periodicity) + # Also catches case where periodicity = true + periodicity = ntuple(_ -> true, NDIMS) + elseif !any(periodicity) + # Also catches case where periodicity = false + periodicity = ntuple(_ -> false, NDIMS) + else + # Default case if periodicity is an iterable + periodicity = Tuple(periodicity) + end + + conn = T8code.Libt8.p4est_connectivity_new_brick(trees_per_dimension..., periodicity...) + do_partition = 0 + cmesh = t8_cmesh_new_from_p4est(conn, mpi_comm(), do_partition) + T8code.Libt8.p4est_connectivity_destroy(conn) + + scheme = t8_scheme_new_default_cxx() + forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) + + basis = LobattoLegendreBasis(RealT, polydeg) + nodes = basis.nodes + + tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, + ntuple(_ -> length(nodes), NDIMS)..., + prod(trees_per_dimension)) + + # Get cell length in reference mesh: Omega_ref = [-1,1]^2. + dx = 2 / trees_per_dimension[1] + dy = 2 / trees_per_dimension[2] + + num_local_trees = t8_cmesh_get_num_local_trees(cmesh) + + # Non-periodic boundaries. + boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) + + for itree in 1:num_local_trees + veptr = t8_cmesh_get_tree_vertices(cmesh, itree - 1) + verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS)) + + # Calculate node coordinates of reference mesh. + cell_x_offset = (verts[1, 1] - 1 / 2 * (trees_per_dimension[1] - 1)) * dx + cell_y_offset = (verts[2, 1] - 1 / 2 * (trees_per_dimension[2] - 1)) * dy + + for j in eachindex(nodes), i in eachindex(nodes) + tree_node_coordinates[:, i, j, itree] .= mapping(cell_x_offset + + dx * nodes[i] / 2, + cell_y_offset + + dy * nodes[j] / 2) + end + + if !periodicity[1] + boundary_names[1, itree] = :x_neg + boundary_names[2, itree] = :x_pos + end + + if !periodicity[2] + boundary_names[3, itree] = :y_neg + boundary_names[4, itree] = :y_pos + end + end + + return T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, + boundary_names, "") +end + +""" + T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}, + mapping=nothing, polydeg=1, RealT=Float64, + initial_refinement_level=0) + +Main mesh constructor for the `T8codeMesh` that imports an unstructured, +conforming mesh from a `t8_cmesh` data structure. + +# Arguments +- `cmesh::Ptr{t8_cmesh}`: Pointer to a cmesh object. +- `mapping`: a function of `NDIMS` variables to describe the mapping that transforms + the imported mesh to the physical domain. Use `nothing` for the identity map. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. + The default of `1` creates an uncurved geometry. Use a higher value if the mapping + will curve the imported uncurved mesh. +- `RealT::Type`: the type that should be used for coordinates. +- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +""" +function T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}; + mapping = nothing, polydeg = 1, RealT = Float64, + initial_refinement_level = 0) where {NDIMS} + @assert NDIMS == 2 # Only support for NDIMS = 2 yet. + + scheme = t8_scheme_new_default_cxx() + forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) + + basis = LobattoLegendreBasis(RealT, polydeg) + nodes = basis.nodes + + num_local_trees = t8_cmesh_get_num_local_trees(cmesh) + + tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, + ntuple(_ -> length(nodes), NDIMS)..., + num_local_trees) + + nodes_in = [-1.0, 1.0] + matrix = polynomial_interpolation_matrix(nodes_in, nodes) + data_in = Array{RealT, 3}(undef, 2, 2, 2) + tmp1 = zeros(RealT, 2, length(nodes), length(nodes_in)) + + for itree in 0:(num_local_trees - 1) + veptr = t8_cmesh_get_tree_vertices(cmesh, itree) + verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS)) + + u = verts[:, 2] - verts[:, 1] + v = verts[:, 3] - verts[:, 1] + w = [0.0, 0.0, 1.0] + + vol = dot(cross(u, v), w) + + if vol < 0.0 + @warn "Discovered negative volumes in `cmesh`: vol = $vol" + end + + # Tree vertices are stored in z-order. + @views data_in[:, 1, 1] .= verts[1:2, 1] + @views data_in[:, 2, 1] .= verts[1:2, 2] + @views data_in[:, 1, 2] .= verts[1:2, 3] + @views data_in[:, 2, 2] .= verts[1:2, 4] + + # Interpolate corner coordinates to specified nodes. + multiply_dimensionwise!(view(tree_node_coordinates, :, :, :, itree + 1), + matrix, matrix, + data_in, + tmp1) + end + + map_node_coordinates!(tree_node_coordinates, mapping) + + # There's no simple and generic way to distinguish boundaries. Name all of them :all. + boundary_names = fill(:all, 2 * NDIMS, num_local_trees) + + return T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, + boundary_names, "") +end + +""" + T8codeMesh{NDIMS}(conn::Ptr{p4est_connectivity}, + mapping=nothing, polydeg=1, RealT=Float64, + initial_refinement_level=0) + +Main mesh constructor for the `T8codeMesh` that imports an unstructured, +conforming mesh from a `p4est_connectivity` data structure. + +# Arguments +- `conn::Ptr{p4est_connectivity}`: Pointer to a P4est connectivity object. +- `mapping`: a function of `NDIMS` variables to describe the mapping that transforms + the imported mesh to the physical domain. Use `nothing` for the identity map. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. + The default of `1` creates an uncurved geometry. Use a higher value if the mapping + will curve the imported uncurved mesh. +- `RealT::Type`: the type that should be used for coordinates. +- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +""" +function T8codeMesh{NDIMS}(conn::Ptr{p4est_connectivity}; kwargs...) where {NDIMS} + @assert NDIMS == 2 # Only support for NDIMS = 2 yet. + + cmesh = t8_cmesh_new_from_p4est(conn, mpi_comm(), 0) + + return T8codeMesh{NDIMS}(cmesh; kwargs...) +end + +""" + T8codeMesh{NDIMS}(meshfile::String; + mapping=nothing, polydeg=1, RealT=Float64, + initial_refinement_level=0) + +Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming +mesh from a Gmsh mesh file (`.msh`). + +# Arguments +- `meshfile::String`: path to a Gmsh mesh file. +- `mapping`: a function of `NDIMS` variables to describe the mapping that transforms + the imported mesh to the physical domain. Use `nothing` for the identity map. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. + The default of `1` creates an uncurved geometry. Use a higher value if the mapping + will curve the imported uncurved mesh. +- `RealT::Type`: the type that should be used for coordinates. +- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +""" +function T8codeMesh{NDIMS}(meshfile::String; kwargs...) where {NDIMS} + @assert NDIMS == 2 # Only support for NDIMS = 2 yet. + + # Prevent `t8code` from crashing Julia if the file doesn't exist. + @assert isfile(meshfile) + + meshfile_prefix, meshfile_suffix = splitext(meshfile) + + cmesh = t8_cmesh_from_msh_file(meshfile_prefix, 0, mpi_comm(), NDIMS, 0, 0) + + return T8codeMesh{NDIMS}(cmesh; kwargs...) +end + +# TODO: Just a placeholder. Will be implemented later when MPI is supported. +function balance!(mesh::T8codeMesh, init_fn = C_NULL) + return nothing +end + +# TODO: Just a placeholder. Will be implemented later when MPI is supported. +function partition!(mesh::T8codeMesh; allow_coarsening = true, weight_fn = C_NULL) + return nothing +end diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 2536cfe0bf..495e0ffc4a 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -363,7 +363,8 @@ function get_element_variables!(element_variables, u, mesh, equations, dg::DG, c dg, cache) end -const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh} +const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, + T8codeMesh} @inline function ndofs(mesh::MeshesDGSEM, dg::DG, cache) nelements(cache.elements) * nnodes(dg)^ndims(mesh) @@ -679,4 +680,5 @@ include("dgsem_tree/dg.jl") include("dgsem_structured/dg.jl") include("dgsem_unstructured/dg.jl") include("dgsem_p4est/dg.jl") +include("dgsem_t8code/dg.jl") end # @muladd diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 2b9c6987d2..0176f5c634 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -81,7 +81,8 @@ function Base.resize!(elements::P4estElementContainer, capacity) end # Create element container and initialize element data -function init_elements(mesh::P4estMesh{NDIMS, RealT}, equations, +function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, RealT}}, + equations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} nelements = ncells(mesh) @@ -165,7 +166,7 @@ function Base.resize!(interfaces::P4estInterfaceContainer, capacity) end # Create interface container and initialize interface data. -function init_interfaces(mesh::P4estMesh, equations, basis, elements) +function init_interfaces(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -240,7 +241,7 @@ function Base.resize!(boundaries::P4estBoundaryContainer, capacity) end # Create interface container and initialize interface data in `elements`. -function init_boundaries(mesh::P4estMesh, equations, basis, elements) +function init_boundaries(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -371,7 +372,7 @@ function Base.resize!(mortars::P4estMortarContainer, capacity) end # Create mortar container and initialize mortar data. -function init_mortars(mesh::P4estMesh, equations, basis, elements) +function init_mortars(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 11747f1f17..236d7d24c0 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -6,7 +6,8 @@ #! format: noindent # Initialize data structures in element container -function init_elements!(elements, mesh::P4estMesh{2}, basis::LobattoLegendreBasis) +function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + basis::LobattoLegendreBasis) @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -25,7 +26,7 @@ end # Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis function calc_node_coordinates!(node_coordinates, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, basis::LobattoLegendreBasis) # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index bc7d9edb6e..97b931fa32 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -7,8 +7,8 @@ # The methods below are specialized on the mortar type # and called from the basic `create_cache` method at the top. -function create_cache(mesh::P4estMesh{2}, equations, mortar_l2::LobattoLegendreMortarL2, - uEltype) +function create_cache(mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, + mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal performance using different types MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, @@ -58,7 +58,7 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack interfaces = cache index_range = eachnode(dg) @@ -114,7 +114,7 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -182,7 +182,7 @@ end # Inlined version of the interface flux computation for conservation laws @inline function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -206,7 +206,7 @@ end # Inlined version of the interface flux computation for equations with conservative and nonconservative terms @inline function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -247,7 +247,7 @@ end end function prolong2boundaries!(cache, u, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack boundaries = cache index_range = eachnode(dg) @@ -276,7 +276,7 @@ function prolong2boundaries!(cache, u, end function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack boundaries = cache @unpack surface_flux_values = cache.elements @@ -312,7 +312,7 @@ end # inlined version of the boundary flux calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, i_index, j_index, @@ -343,7 +343,7 @@ end # inlined version of the boundary flux with nonconservative terms calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, i_index, j_index, @@ -385,7 +385,7 @@ end end function prolong2mortars!(cache, u, - mesh::P4estMesh{2}, equations, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DGSEM) @unpack neighbor_ids, node_indices = cache.mortars @@ -452,7 +452,7 @@ function prolong2mortars!(cache, u, end function calc_mortar_flux!(surface_flux_values, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -511,7 +511,7 @@ end # Inlined version of the mortar flux computation on small elements for conservation laws @inline function calc_mortar_flux!(fstar, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -531,7 +531,7 @@ end # Inlined version of the mortar flux computation on small elements for equations with conservative and # nonconservative terms @inline function calc_mortar_flux!(fstar, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -559,7 +559,8 @@ end end @inline function mortar_fluxes_to_elements!(surface_flux_values, - mesh::P4estMesh{2}, equations, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, mortar, fstar, u_buffer) @unpack neighbor_ids, node_indices = cache.mortars @@ -620,7 +621,7 @@ end end function calc_surface_integral!(du, u, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index c013bf62d9..3e8ce759b3 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -52,7 +52,7 @@ end @inline function weak_form_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -93,8 +93,8 @@ end @inline function flux_differencing_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2} - }, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) @unpack derivative_split = dg.basis @@ -150,8 +150,8 @@ end @inline function flux_differencing_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2} - }, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms::True, equations, volume_flux, dg::DGSEM, cache, alpha = true) @unpack derivative_split = dg.basis @@ -219,7 +219,7 @@ end # [arXiv: 2008.12044v2](https://arxiv.org/pdf/2008.12044) @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @@ -289,7 +289,7 @@ end @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, u::AbstractArray{<:Any, 4}, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @@ -609,9 +609,8 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, end function apply_jacobian!(du, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2 - } - }, + mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, + P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) @unpack inverse_jacobian = cache.elements diff --git a/src/solvers/dgsem_t8code/containers.jl b/src/solvers/dgsem_t8code/containers.jl new file mode 100644 index 0000000000..093feb2985 --- /dev/null +++ b/src/solvers/dgsem_t8code/containers.jl @@ -0,0 +1,60 @@ +function reinitialize_containers!(mesh::T8codeMesh, equations, dg::DGSEM, cache) + # Re-initialize elements container. + @unpack elements = cache + resize!(elements, ncells(mesh)) + init_elements!(elements, mesh, dg.basis) + + count_required_surfaces!(mesh) + + # Resize interfaces container. + @unpack interfaces = cache + resize!(interfaces, mesh.ninterfaces) + + # Resize mortars container. + @unpack mortars = cache + resize!(mortars, mesh.nmortars) + + # Resize boundaries container. + @unpack boundaries = cache + resize!(boundaries, mesh.nboundaries) + + trixi_t8_fill_mesh_info(mesh.forest, elements, interfaces, mortars, boundaries, + mesh.boundary_names) + + return nothing +end + +function count_required_surfaces!(mesh::T8codeMesh) + counts = trixi_t8_count_interfaces(mesh.forest) + + mesh.nmortars = counts.mortars + mesh.ninterfaces = counts.interfaces + mesh.nboundaries = counts.boundaries + + return counts +end + +# Compatibility to `dgsem_p4est/containers.jl`. +function count_required_surfaces(mesh::T8codeMesh) + return (interfaces = mesh.ninterfaces, + mortars = mesh.nmortars, + boundaries = mesh.nboundaries) +end + +# Compatibility to `dgsem_p4est/containers.jl`. +function init_interfaces!(interfaces, mesh::T8codeMesh) + # Already computed. Do nothing. + return nothing +end + +# Compatibility to `dgsem_p4est/containers.jl`. +function init_mortars!(mortars, mesh::T8codeMesh) + # Already computed. Do nothing. + return nothing +end + +# Compatibility to `dgsem_p4est/containers.jl`. +function init_boundaries!(boundaries, mesh::T8codeMesh) + # Already computed. Do nothing. + return nothing +end diff --git a/src/solvers/dgsem_t8code/containers_2d.jl b/src/solvers/dgsem_t8code/containers_2d.jl new file mode 100644 index 0000000000..029e6674af --- /dev/null +++ b/src/solvers/dgsem_t8code/containers_2d.jl @@ -0,0 +1,58 @@ +@muladd begin +#! format: noindent + +# Interpolate tree_node_coordinates to each quadrant at the specified nodes. +function calc_node_coordinates!(node_coordinates, + mesh::T8codeMesh{2}, + nodes::AbstractVector) + # 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. + tmp1 = StrideArray(undef, real(mesh), + StaticInt(2), static_length(nodes), static_length(mesh.nodes)) + matrix1 = StrideArray(undef, real(mesh), + static_length(nodes), static_length(mesh.nodes)) + matrix2 = similar(matrix1) + baryweights_in = barycentric_weights(mesh.nodes) + + num_local_trees = t8_forest_get_num_local_trees(mesh.forest) + + current_index = 0 + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(mesh.forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class) + num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree) + + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) + element_level = t8_element_level(eclass_scheme, element) + + element_length = t8_quad_len(element_level) / t8_quad_root_len + + element_coords = Array{Float64}(undef, 3) + t8_element_vertex_reference_coords(eclass_scheme, element, 0, + pointer(element_coords)) + + nodes_out_x = 2 * + (element_length * 1 / 2 * (nodes .+ 1) .+ element_coords[1]) .- + 1 + nodes_out_y = 2 * + (element_length * 1 / 2 * (nodes .+ 1) .+ element_coords[2]) .- + 1 + + polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, + baryweights_in) + polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, + baryweights_in) + + multiply_dimensionwise!(view(node_coordinates, :, :, :, current_index += 1), + matrix1, matrix2, + view(mesh.tree_node_coordinates, :, :, :, + itree + 1), + tmp1) + end + end + + return node_coordinates +end +end # @muladd diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl new file mode 100644 index 0000000000..16a9d7d35b --- /dev/null +++ b/src/solvers/dgsem_t8code/dg.jl @@ -0,0 +1,31 @@ +@muladd begin +#! format: noindent + +# This method is called when a SemidiscretizationHyperbolic is constructed. +# It constructs the basic `cache` used throughout the simulation to compute +# the RHS etc. +function create_cache(mesh::T8codeMesh, equations::AbstractEquations, dg::DG, ::Any, + ::Type{uEltype}) where {uEltype <: Real} + count_required_surfaces!(mesh) + + elements = init_elements(mesh, equations, dg.basis, uEltype) + interfaces = init_interfaces(mesh, equations, dg.basis, elements) + boundaries = init_boundaries(mesh, equations, dg.basis, elements) + mortars = init_mortars(mesh, equations, dg.basis, elements) + + trixi_t8_fill_mesh_info(mesh.forest, elements, interfaces, mortars, boundaries, + mesh.boundary_names) + + cache = (; elements, interfaces, boundaries, mortars) + + # Add specialized parts of the cache required to compute the volume integral etc. + cache = (; cache..., + create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) + cache = (; cache..., create_cache(mesh, equations, dg.mortar, uEltype)...) + + return cache +end + +include("containers.jl") +include("containers_2d.jl") +end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 6c5e0cee0c..c30d0a8e01 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -37,14 +37,14 @@ end # The methods below are specialized on the volume integral type # and called from the basic `create_cache` method at the top. function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype) NamedTuple() end function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, equations, + P4estMesh{2}, T8codeMesh{2}}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] element_ids_dgfv = Int[] @@ -70,7 +70,7 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMe end function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, equations, + P4estMesh{2}, T8codeMesh{2}}, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG, uEltype) A3dp1_x = Array{uEltype, 3} @@ -92,7 +92,7 @@ end # The methods below are specialized on the mortar type # and called from the basic `create_cache` method at the top. function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}}, + P4estMesh{2}, T8codeMesh{2}}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal performance using different types MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, @@ -110,7 +110,7 @@ end # TODO: Taal discuss/refactor timer, allowing users to pass a custom timer? function rhs!(du, u, t, - mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, + mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -180,7 +180,8 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -226,7 +227,8 @@ end # from the evaluation of the physical fluxes in each Cartesian direction function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) @@ -322,7 +324,8 @@ end # TODO: Taal dimension agnostic function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DGSEM, cache) @@ -381,7 +384,8 @@ end @inline function fv_kernel!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2} + }, nonconservative_terms, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded = cache diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index f7c7854717..2f34e0eb66 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -208,7 +208,8 @@ end end # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha -function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}}, alpha, alpha_tmp, dg, +function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, alpha, + alpha_tmp, dg, cache) # Copy alpha values such that smoothing is indpedenent of the element access order alpha_tmp .= alpha diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 95dec027a8..7b8dafdddd 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -307,14 +307,14 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::Union{UnstructuredMesh2D, P4estMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, T8codeMesh}, equations, surface_integral, dg::DG) @assert isempty(eachboundary(dg, cache)) end # Function barrier for type stability function calc_boundary_flux!(cache, t, boundary_conditions, - mesh::Union{UnstructuredMesh2D, P4estMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, T8codeMesh}, equations, surface_integral, dg::DG) @unpack boundary_condition_types, boundary_indices = boundary_conditions @@ -327,7 +327,8 @@ end # in a type-stable way using "lispy tuple programming". function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any}, BC_indices::NTuple{N, Vector{Int}}, - mesh::Union{UnstructuredMesh2D, P4estMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, + T8codeMesh}, equations, surface_integral, dg::DG) where {N} # Extract the boundary condition type and index vector boundary_condition = first(BCs) @@ -350,7 +351,8 @@ end # terminate the type-stable iteration over tuples function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, - mesh::Union{UnstructuredMesh2D, P4estMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, + T8codeMesh}, equations, surface_integral, dg::DG) nothing end diff --git a/test/runtests.jl b/test/runtests.jl index f76811dddb..1d7eefe1fc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,113 +4,119 @@ using MPI: mpiexec # run tests on Travis CI in parallel const TRIXI_TEST = get(ENV, "TRIXI_TEST", "all") const TRIXI_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3) -const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) +const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time @testset "Trixi.jl tests" begin - # This is placed first since tests error out otherwise if `TRIXI_TEST == "all"`, - # at least on some systems. - @time if TRIXI_TEST == "all" || TRIXI_TEST == "mpi" - # Do a dummy `@test true`: - # If the process errors out the testset would error out as well, - # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 - @test true - - # There are spurious test failures of Trixi.jl with MPI on Windows, see - # https://github.com/trixi-framework/Trixi.jl/issues/901 - # To reduce their impact, we do not test MPI with coverage on Windows. - # This reduces the chance to hit a spurious test failure by one half. - # In addition, it looks like the Linux GitHub runners run out of memory during the 3D tests - # with coverage, so we currently do not test MPI with coverage on Linux. For more details, - # see the discussion at https://github.com/trixi-framework/Trixi.jl/pull/1062#issuecomment-1035901020 - cmd = string(Base.julia_cmd()) - coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", cmd) - if !(coverage && Sys.iswindows()) && !(coverage && Sys.islinux()) - # We provide a `--heap-size-hint` to avoid/reduce out-of-memory errors during CI testing - mpiexec() do cmd - run(`$cmd -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --threads=1 --check-bounds=yes --heap-size-hint=1G $(abspath("test_mpi.jl"))`) - end - end - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "threaded" || TRIXI_TEST == "threaded_legacy" - # Do a dummy `@test true`: - # If the process errors out the testset would error out as well, - # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 - @test true - - run(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none $(abspath("test_threaded.jl"))`) - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part1" - include("test_tree_1d.jl") - include("test_tree_2d_part1.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part2" - include("test_tree_2d_part2.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part3" - include("test_tree_2d_part3.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part4" - include("test_tree_3d_part1.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part5" - include("test_tree_3d_part2.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part6" - include("test_tree_3d_part3.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "structured" - include("test_structured_1d.jl") - include("test_structured_2d.jl") - include("test_structured_3d.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "p4est_part1" - include("test_p4est_2d.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "p4est_part2" - include("test_p4est_3d.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "unstructured_dgmulti" - include("test_unstructured_2d.jl") - include("test_dgmulti_1d.jl") - include("test_dgmulti_2d.jl") - include("test_dgmulti_3d.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "parabolic" - include("test_parabolic_1d.jl") - include("test_parabolic_2d.jl") - include("test_parabolic_3d.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1" - include("test_unit.jl") - include("test_visualization.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part2" - include("test_special_elixirs.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "performance_specializations_part1" - include("test_performance_specializations_2d.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "performance_specializations_part2" - include("test_performance_specializations_3d.jl") - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "paper_self_gravitating_gas_dynamics" - include("test_paper_self_gravitating_gas_dynamics.jl") - end + # This is placed first since tests error out otherwise if `TRIXI_TEST == "all"`, + # at least on some systems. + @time if TRIXI_TEST == "all" || TRIXI_TEST == "mpi" + # Do a dummy `@test true`: + # If the process errors out the testset would error out as well, + # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 + @test true + + # There are spurious test failures of Trixi.jl with MPI on Windows, see + # https://github.com/trixi-framework/Trixi.jl/issues/901 + # To reduce their impact, we do not test MPI with coverage on Windows. + # This reduces the chance to hit a spurious test failure by one half. + # In addition, it looks like the Linux GitHub runners run out of memory during the 3D tests + # with coverage, so we currently do not test MPI with coverage on Linux. For more details, + # see the discussion at https://github.com/trixi-framework/Trixi.jl/pull/1062#issuecomment-1035901020 + cmd = string(Base.julia_cmd()) + coverage = occursin("--code-coverage", cmd) && + !occursin("--code-coverage=none", cmd) + if !(coverage && Sys.iswindows()) && !(coverage && Sys.islinux()) + # We provide a `--heap-size-hint` to avoid/reduce out-of-memory errors during CI testing + mpiexec() do cmd + run(`$cmd -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --threads=1 --check-bounds=yes --heap-size-hint=1G $(abspath("test_mpi.jl"))`) + end + end + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "threaded" || + TRIXI_TEST == "threaded_legacy" + # Do a dummy `@test true`: + # If the process errors out the testset would error out as well, + # cf. https://github.com/JuliaParallel/MPI.jl/pull/391 + @test true + + run(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none $(abspath("test_threaded.jl"))`) + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part1" + include("test_tree_1d.jl") + include("test_tree_2d_part1.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part2" + include("test_tree_2d_part2.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part3" + include("test_tree_2d_part3.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part4" + include("test_tree_3d_part1.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part5" + include("test_tree_3d_part2.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part6" + include("test_tree_3d_part3.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "structured" + include("test_structured_1d.jl") + include("test_structured_2d.jl") + include("test_structured_3d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "p4est_part1" + include("test_p4est_2d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "p4est_part2" + include("test_p4est_3d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "t8code_part1" + include("test_t8code_2d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "unstructured_dgmulti" + include("test_unstructured_2d.jl") + include("test_dgmulti_1d.jl") + include("test_dgmulti_2d.jl") + include("test_dgmulti_3d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "parabolic" + include("test_parabolic_1d.jl") + include("test_parabolic_2d.jl") + include("test_parabolic_3d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1" + include("test_unit.jl") + include("test_visualization.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part2" + include("test_special_elixirs.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "performance_specializations_part1" + include("test_performance_specializations_2d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "performance_specializations_part2" + include("test_performance_specializations_3d.jl") + end + + @time if TRIXI_TEST == "all" || TRIXI_TEST == "paper_self_gravitating_gas_dynamics" + include("test_paper_self_gravitating_gas_dynamics.jl") + end end diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl new file mode 100644 index 0000000000..a424c9df84 --- /dev/null +++ b/test/test_t8code_2d.jl @@ -0,0 +1,182 @@ +module TestExamplesT8codeMesh2D + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = joinpath(examples_dir(), "t8code_2d_dgsem") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) +mkdir(outdir) + +@testset "T8codeMesh2D" begin + + @trixi_testset "test save_mesh_file" begin + @test_throws Exception begin + # Save mesh file support will be added in the future. The following + # lines of code are here for satisfying code coverage. + + # Create dummy mesh. + mesh = T8codeMesh((1, 1), polydeg = 1, + mapping = Trixi.coordinates2mapping((-1.0, -1.0), ( 1.0, 1.0)), + initial_refinement_level = 1) + + # This call throws an error. + Trixi.save_mesh_file(mesh, "dummy") + end + end + + @trixi_testset "elixir_advection_basic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[8.311947673061856e-6], + linf=[6.627000273229378e-5]) + end + + @trixi_testset "elixir_advection_nonconforming_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_nonconforming_flag.jl"), + l2=[3.198940059144588e-5], + linf=[0.00030636069494005547]) + end + + @trixi_testset "elixir_advection_unstructured_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_flag.jl"), + l2=[0.0005379687442422346], + linf=[0.007438525029884735]) + end + + @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_amr_unstructured_flag.jl"), + l2=[0.001993165013217687], + linf=[0.032891018571625796], + coverage_override=(maxiters = 6,)) + end + + @trixi_testset "elixir_advection_amr_solution_independent.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_amr_solution_independent.jl"), + # Expected errors are exactly the same as with StructuredMesh! + l2=[4.949660644033807e-5], + linf=[0.0004867846262313763], + coverage_override=(maxiters = 6,)) + end + + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), + l2=[ + 0.0034516244508588046, + 0.0023420334036925493, + 0.0024261923964557187, + 0.004731710454271893, + ], + linf=[ + 0.04155789011775046, + 0.024772109862748914, + 0.03759938693042297, + 0.08039824959535657, + ]) + end + + @trixi_testset "elixir_euler_free_stream.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), + l2=[ + 2.063350241405049e-15, + 1.8571016296925367e-14, + 3.1769447886391905e-14, + 1.4104095258528071e-14, + ], + linf=[1.9539925233402755e-14, 2e-12, 4.8e-12, 4e-12], + atol=2.0e-12,) + end + + @trixi_testset "elixir_euler_shockcapturing_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_ec.jl"), + l2=[ + 9.53984675e-02, + 1.05633455e-01, + 1.05636158e-01, + 3.50747237e-01, + ], + linf=[ + 2.94357464e-01, + 4.07893014e-01, + 3.97334516e-01, + 1.08142520e+00, + ], + tspan=(0.0, 1.0)) + end + + @trixi_testset "elixir_euler_sedov.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), + l2=[ + 3.76149952e-01, + 2.46970327e-01, + 2.46970327e-01, + 1.28889042e+00, + ], + linf=[ + 1.22139001e+00, + 1.17742626e+00, + 1.17742626e+00, + 6.20638482e+00, + ], + tspan=(0.0, 0.3)) + end + + @trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 9.168126407325352e-5, + 0.0009795410115453788, + 0.002546408320320785, + 3.941189812642317e-6, + ], + linf=[ + 0.0009903782521019089, + 0.0059752684687262025, + 0.010941106525454103, + 1.2129488214718265e-5, + ], + tspan=(0.0, 0.1)) + end + + @trixi_testset "elixir_mhd_alfven_wave.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), + l2=[1.0513414461545583e-5, 1.0517900957166411e-6, + 1.0517900957304043e-6, 1.511816606372376e-6, + 1.0443997728645063e-6, 7.879639064990798e-7, + 7.879639065049896e-7, 1.0628631669056271e-6, + 4.3382328912336153e-7], + linf=[4.255466285174592e-5, 1.0029706745823264e-5, + 1.0029706747467781e-5, 1.2122265939010224e-5, + 5.4791097160444835e-6, 5.18922042269665e-6, + 5.189220422141538e-6, 9.552667261422676e-6, + 1.4237578427628152e-6]) + end + + @trixi_testset "elixir_mhd_rotor.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), + l2=[0.44211360369891683, 0.8805178316216257, 0.8262710688468049, + 0.0, + 0.9616090460973586, 0.10386643568745411, + 0.15403457366543802, 0.0, + 2.8399715649715473e-5], + linf=[10.04369305341599, 17.995640564998403, 9.576041548174265, + 0.0, + 19.429658884314534, 1.3821395681242314, 1.818559351543182, + 0.0, + 0.002261930217575465], + tspan=(0.0, 0.02)) + end +end + +# Clean up afterwards: delete Trixi.jl output directory +@test_nowarn rm(outdir, recursive = true) + +end # module diff --git a/test/test_threaded.jl b/test/test_threaded.jl index 77fa16ad33..9b30836d0e 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -235,6 +235,22 @@ Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive=true) end + @testset "T8codeMesh" begin + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin + @test_trixi_include(joinpath(examples_dir(), "t8code_2d_dgsem", "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), + l2 = [0.0034516244508588046, 0.0023420334036925493, 0.0024261923964557187, 0.004731710454271893], + linf = [0.04155789011775046, 0.024772109862748914, 0.03759938693042297, 0.08039824959535657]) + end + + @trixi_testset "elixir_eulergravity_convergence.jl" begin + @test_trixi_include(joinpath(examples_dir(), "t8code_2d_dgsem", "elixir_eulergravity_convergence.jl"), + l2 = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697], + linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758], + tspan = (0.0, 0.1)) + end + end + + @testset "DGMulti" begin @trixi_testset "elixir_euler_weakform.jl (SBP, EC)" begin @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_euler_weakform.jl"),