diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2c0ea798b49..f287cc5feb2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -70,6 +70,7 @@ jobs: - p4est_part1 - p4est_part2 - t8code_part1 + - t8code_part2 - unstructured_dgmulti - parabolic - paper_self_gravitating_gas_dynamics diff --git a/.gitignore b/.gitignore index 3132b9af38b..b4f1cf6bb47 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ *.mesh *.bson *.inp +*.msh **/Manifest.toml out*/ docs/build diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl index 653bab41e2d..0589e76a6a9 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -93,12 +93,10 @@ 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, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 1) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index adf1d009a59..f285d24fc6c 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -41,9 +41,9 @@ isfile(mesh_file) || # 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) +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 1) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) diff --git a/examples/t8code_2d_dgsem/elixir_advection_basic.jl b/examples/t8code_2d_dgsem/elixir_advection_basic.jl index efc51226586..26ced0970fe 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_basic.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_basic.jl @@ -16,12 +16,10 @@ 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, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 1) # A semidiscretization collects data structures and functions for the spatial discretization diff --git a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl index 31a8bc93697..a39f3a7e195 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl @@ -20,31 +20,28 @@ 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, +# Create T8codeMesh 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)) +# Note: This is actually a `p4est_quadrant_t` which is much bigger than the +# following struct. But we only need the first three fields for our purpose. +struct t8_dquad_t + x::Int32 + y::Int32 + level::Int8 + # [...] # See `p4est.h` in `p4est` for more info. +end - level = Trixi.t8_element_level(ts, elements[1]) +# Refine quadrants of each tree at lower left edge to level 4. +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dquad_t}(elements[1])) - # TODO: Make this condition more general. - if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 4 + if el.x == 0 && el.y == 0 && el.level < 4 # return true (refine) return 1 else @@ -53,26 +50,7 @@ function adapt_callback(forest, 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 +Trixi.adapt!(mesh, adapt_callback) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, diff --git a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl index df9cbc26f6e..5ba1ab15489 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl @@ -38,9 +38,9 @@ isfile(mesh_file) || # 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) +mesh = T8codeMesh(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, diff --git a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl index 01e0449c67e..37d15f38566 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl @@ -40,25 +40,17 @@ isfile(mesh_file) || # 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) +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping, + initial_refinement_level = 1) - elements = unsafe_wrap(Array, elements_ptr, num_elements) +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + vertex = Vector{Cdouble}(undef, 3) - Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) + Trixi.t8_element_vertex_reference_coords(eclass_scheme, elements[1], 0, vertex) - level = Trixi.t8_element_level(ts, elements[1]) + level = Trixi.t8_element_level(eclass_scheme, elements[1]) # TODO: Make this condition more general. if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 3 @@ -70,26 +62,7 @@ function adapt_callback(forest, 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 +Trixi.adapt!(mesh, adapt_callback) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition))) diff --git a/examples/t8code_2d_dgsem/elixir_euler_sedov.jl b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl index 965d794f8dc..82770a4050b 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl @@ -58,12 +58,10 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, 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, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 2, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl index 55a9063a001..9ebbd1d28c4 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl @@ -29,12 +29,10 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, 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, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 2, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) 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 index 21f26d79ba8..bcc1abc560e 100644 --- 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 @@ -40,25 +40,17 @@ isfile(mesh_file) || # 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) +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 1) - elements = unsafe_wrap(Array, elements_ptr, num_elements) +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + vertex = Vector{Cdouble}(undef, 3) - Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) + Trixi.t8_element_vertex_reference_coords(eclass_scheme, elements[1], 0, pointer(vertex)) - level = Trixi.t8_element_level(ts, elements[1]) + level = Trixi.t8_element_level(eclass_scheme, elements[1]) # TODO: Make this condition more general. if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 2 @@ -70,26 +62,7 @@ function adapt_callback(forest, 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 +Trixi.adapt!(mesh, adapt_callback) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms, diff --git a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl index 6d6bb27e0c3..98a9a5521a9 100644 --- a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl +++ b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl @@ -16,10 +16,8 @@ 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, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 2) semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition, diff --git a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl index 1e2362a123c..e184cb3fd05 100644 --- a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl +++ b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl @@ -11,18 +11,17 @@ 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_hlle, 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, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 0, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl index 9a4bd99e444..adb154948fb 100644 --- a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl +++ b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl @@ -78,9 +78,9 @@ isfile(mesh_file) || # 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) +mesh = T8codeMesh(conn, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 1) boundary_condition = BoundaryConditionDirichlet(initial_condition) boundary_conditions = Dict(:all => boundary_condition) diff --git a/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl index b2d5097036f..3610639d554 100644 --- a/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl +++ b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -22,12 +22,10 @@ solver = DGSEM(polydeg = 3, 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, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 1) # A semidiscretization collects data structures and functions for the spatial discretization diff --git a/examples/t8code_3d_dgsem/elixir_advection_amr.jl b/examples/t8code_3d_dgsem/elixir_advection_amr.jl new file mode 100644 index 00000000000..5a4b2218d57 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_amr.jl @@ -0,0 +1,66 @@ +# The same setup as tree_3d_dgsem/elixir_advection_amr.jl +# to verify the T8codeMesh implementation against TreeMesh. + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +initial_condition = initial_condition_gauss +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-5.0, -5.0, -5.0) +coordinates_max = (5.0, 5.0, 5.0) +trees_per_dimension = (1, 1, 1) + +# Note that it is not necessary to use mesh polydeg lower than the solver polydeg +# on a Cartesian mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +mesh = T8codeMesh(trees_per_dimension, polydeg = 1, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 4) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.3) +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 = 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.2) + +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_3d_dgsem/elixir_advection_amr_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl new file mode 100644 index 00000000000..617736afbdd --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl @@ -0,0 +1,105 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# Solver with polydeg=4 to ensure free stream preservation (FSP) on non-conforming meshes. +# The polydeg of the solver must be at least twice as big as the polydeg of the mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs) + +initial_condition = initial_condition_gauss +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +boundary_conditions = Dict(:all => boundary_condition) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. The original mapping applied to this unstructured mesh +# causes some Jacobians to be negative, which makes the mesh invalid. +function mapping(xi, eta, zeta) + # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries + # xi = 1.5 * xi_ + 1.5 + # eta = 1.5 * eta_ + 1.5 + # zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 4 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 4 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 4 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + # Transform the weird deformed cube to be approximately the size of [-5,5]^3 to match IC + return SVector(5 * x, 5 * y, 5 * z) +end + +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 2, + mapping = mapping, + initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 8.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 = 1.2) + +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_3d_dgsem/elixir_advection_basic.jl b/examples/t8code_3d_dgsem/elixir_advection_basic.jl new file mode 100644 index 00000000000..f49462035aa --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_basic.jl @@ -0,0 +1,59 @@ +# The same setup as tree_3d_dgsem/elixir_advection_basic.jl +# to verify the T8codeMesh implementation against TreeMesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(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, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# Create P4estMesh with 8 x 8 x 8 elements (note `refinement_level=1`) +trees_per_dimension = (4, 4, 4) +mesh = T8codeMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + 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 +tspan = (0.0, 1.0) +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.2) + +# 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_3d_dgsem/elixir_advection_nonconforming.jl b/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl new file mode 100644 index 00000000000..8d7a48370f5 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl @@ -0,0 +1,85 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the linear advection equation. + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(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, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) +trees_per_dimension = (1, 1, 1) + +# Note that it is not necessary to use mesh polydeg lower than the solver polydeg +# on a Cartesian mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +mesh = T8codeMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 2) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +# Refine bottom left quadrant of each second tree to level 2 +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if iseven(convert(Int, ltreeid)) && el.x == 0 && el.y == 0 && el.z == 0 && + el.level < 3 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +# 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 +tspan = (0.0, 1.0) +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.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_3d_dgsem/elixir_advection_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl new file mode 100644 index 00000000000..df358435c9a --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl @@ -0,0 +1,98 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(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) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. The original mapping applied to this unstructured mesh +# causes some Jacobians to be negative, which makes the mesh invalid. +function mapping(xi, eta, zeta) + # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries + # xi = 1.5 * xi_ + 1.5 + # eta = 1.5 * eta_ + 1.5 + # zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 6 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 6 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 6 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + return SVector(x, y, z) +end + +# Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping, + 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.1 +ode = semidiscretize(semi, (0.0, 0.1)); + +# 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_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 1.2) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + 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_3d_dgsem/elixir_euler_ec.jl b/examples/t8code_3d_dgsem/elixir_euler_ec.jl new file mode 100644 index 00000000000..07745c3ac56 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_ec.jl @@ -0,0 +1,92 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(5 / 3) + +initial_condition = initial_condition_weak_blast_wave + +boundary_conditions = Dict(:all => boundary_condition_slip_wall) + +# Get the DG approximation space + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 5, surface_flux = flux_ranocha, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# Get the curved quad mesh from a file + +# Mapping as described in https://arxiv.org/abs/2012.12040 +function mapping(xi_, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 3 / 8 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 3 / 8 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 3 / 8 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + return SVector(x, y, z) +end + +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 5, + mapping = mapping, + initial_refinement_level = 0) + +# Create the semidiscretization object. +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# 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_3d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl new file mode 100644 index 00000000000..e135d464810 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl @@ -0,0 +1,118 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_constant + +boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition)) + +# Solver with polydeg=4 to ensure free stream preservation (FSP) on non-conforming meshes. +# The polydeg of the solver must be at least twice as big as the polydeg of the mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. This can yield problematic geometries if the unrefined mesh +# is not fine enough. +function mapping(xi_, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 6 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 6 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 6 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + return SVector(x, y, z) +end + +# Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 2, + mapping = mapping, + initial_refinement_level = 0) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +# Refine bottom left quadrant of each second tree to level 2 +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if iseven(convert(Int, ltreeid)) && el.x == 0 && el.y == 0 && el.z == 0 && + el.level < 2 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + 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 = 1.2) + +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_3d_dgsem/elixir_euler_free_stream_extruded.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl new file mode 100644 index 00000000000..d129b59826e --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl @@ -0,0 +1,106 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_constant + +boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition)) + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D. +# This particular mesh is unstructured in the yz-plane, but extruded in x-direction. +# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded +# in x-direction to ensure free stream preservation on a non-conforming mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +function mapping(xi, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + z = zeta + + 1 / 6 * (cos(1.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + y = eta + 1 / 6 * (cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(2 * pi * (2 * z - 3) / 3)) + + return SVector(xi, y, z) +end + +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_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(3)) + +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping, + initial_refinement_level = 0) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +# Refine quadrants in y-direction of each tree at one edge to level 2 +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if convert(Int, ltreeid) < 4 && el.x == 0 && el.y == 0 && el.level < 2 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + 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 = 1.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), #maxiters=1, + 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_3d_dgsem/elixir_euler_sedov.jl b/examples/t8code_3d_dgsem/elixir_euler_sedov.jl new file mode 100644 index 00000000000..618b170b661 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_sedov.jl @@ -0,0 +1,97 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +""" + initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) + +The Sedov blast wave setup based on Flash +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 +with smaller strength of the initial discontinuity. +""" +function initial_condition_medium_sedov_blast_wave(x, t, + equations::CompressibleEulerEquations3D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 + r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) + p0_outer = 1.0e-3 + + # Calculate primitive variables + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +initial_condition = initial_condition_medium_sedov_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 5 +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, -1.0) +coordinates_max = (1.0, 1.0, 1.0) + +trees_per_dimension = (4, 4, 4) +mesh = T8codeMesh(trees_per_dimension, + polydeg = 4, initial_refinement_level = 0, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = true) + +# create the semi discretization object +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 = 100 +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_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl new file mode 100644 index 00000000000..d4664522bea --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl @@ -0,0 +1,120 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +# Solver with polydeg=4 to ensure free stream preservation (FSP) on non-conforming meshes. +# The polydeg of the solver must be at least twice as big as the polydeg of the mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. The original mapping applied to this unstructured mesh +# causes some Jacobians to be negative, which makes the mesh invalid. +function mapping(xi, eta, zeta) + # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries + # xi = 1.5 * xi_ + 1.5 + # eta = 1.5 * eta_ + 1.5 + # zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 6 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 6 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 6 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + # Transform the weird deformed cube to be approximately the cube [0,2]^3 + return SVector(x + 1, y + 1, z + 1) +end + +# Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_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(3)) + +# Mesh polydeg of 2 (half the solver polydeg) to ensure FSP (see above). +mesh = T8codeMesh(conn, polydeg = 2, + mapping = mapping, + initial_refinement_level = 0) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if el.x == 0 && el.y == 0 && el.z == 0 && el.level < 2 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.045) +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.6) + +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_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl new file mode 100644 index 00000000000..7cb03bb312d --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -0,0 +1,62 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:x_neg => boundary_condition, + :x_pos => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :z_neg => boundary_condition, + :z_pos => boundary_condition) + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (2.0, 2.0, 2.0) + +trees_per_dimension = (2, 2, 2) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +mesh = T8codeMesh(trees_per_dimension, polydeg = 1, + mapping = mapping, + periodicity = false, initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.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.6) + +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/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl index bd781b21c1e..db01476bb86 100644 --- a/src/auxiliary/t8code.jl +++ b/src/auxiliary/t8code.jl @@ -35,7 +35,7 @@ function init_t8code() # 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." + @info "T8code.jl: `sc_finalize` will be called during shutdown of Trixi.jl." MPI.add_finalize_hook!(T8code.Libt8.sc_finalize) end else @@ -116,7 +116,6 @@ function trixi_t8_count_interfaces(forest) elseif level < neighbor_level local_num_mortars += 1 end - else local_num_boundary += 1 end @@ -219,38 +218,9 @@ function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundari 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 - + # Save interfaces.node_indices dimension specific in containers_3d.jl. + init_interface_node_indices!(interfaces, faces, orientation, + interface_id) # Non-conforming interface. elseif level < neighbor_level local_num_mortars += 1 @@ -262,42 +232,13 @@ function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundari # 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 + # Fill in the `mortars.neighbor_ids` array and reorder if necessary. + init_mortar_neighbor_ids!(mortars, faces[2], faces[1], + orientation, neighbor_ielements, + mortar_id) + + # Fill in the `mortars.node_indices` array. + init_mortar_node_indices!(mortars, faces, orientation, mortar_id) # else: "level > neighbor_level" is skipped since we visit the mortar interface only once. end @@ -309,19 +250,7 @@ function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundari 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 + init_boundary_node_indices!(boundaries, iface, boundary_id) # One-based indexing. boundaries.name[boundary_id] = boundary_names[iface + 1, itree + 1] @@ -420,13 +349,15 @@ function trixi_t8_adapt_new(old_forest, indicators) 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 + let set_from = C_NULL, recursive = 0, set_for_coarsening = 0, no_repartition = 0, + do_ghost = 1 + 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_set_ghost(new_forest, do_ghost, T8_GHOST_FACES) # Note: MPI support not available yet so it is a dummy call. t8_forest_commit(new_forest) end diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 98e531295b7..b816bc06e65 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -396,7 +396,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, old_index = 1 new_index = 1 - # Note: This is true for `quads` only. + # Note: This is true for `quads`. T8_CHILDREN = 4 # Retain current solution data. diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index c8abe6fdb05..392cbba9e28 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -304,9 +304,84 @@ end # this method is called when an `ControllerThreeLevel` is constructed function create_cache(::Type{ControllerThreeLevel}, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, dg::DG, cache) controller_value = Vector{Int}(undef, nelements(dg, cache)) return (; controller_value) end + +# Coarsen and refine elements in the DG solver based on a difference list. +function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, 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 only true for `hexs`. + T8_CHILDREN = 8 + + # 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) * ndofs(mesh, dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + + 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, + u_tmp1, u_tmp2) + + 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, u_tmp1, u_tmp2) + + 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 end # @muladd diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 81d0795a159..27e8a2b722f 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -35,7 +35,9 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{3}, return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2) end -function create_cache_analysis(analyzer, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, +function create_cache_analysis(analyzer, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DG, cache, RealT, uEltype) @@ -118,7 +120,7 @@ function calc_error_norms(func, u, t, analyzer, end function calc_error_norms(func, u, t, analyzer, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @@ -190,7 +192,8 @@ function integrate_via_indices(func::Func, u, end function integrate_via_indices(func::Func, u, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis @@ -218,7 +221,8 @@ function integrate_via_indices(func::Func, u, end function integrate(func::Func, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DG, cache; normalize = true) where {Func} integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, k, element, equations, dg @@ -248,7 +252,8 @@ function integrate(func::Func, u, end function analyze(::typeof(entropy_timederivative), du, u, t, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DG, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, @@ -277,7 +282,7 @@ function analyze(::Val{:l2_divb}, du, u, t, end function analyze(::Val{:l2_divb}, du, u, t, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations::IdealGlmMhdEquations3D, dg::DGSEM, cache) @unpack contravariant_vectors = cache.elements @@ -333,7 +338,7 @@ function analyze(::Val{:linf_divb}, du, u, t, end function analyze(::Val{:linf_divb}, du, u, t, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations::IdealGlmMhdEquations3D, dg::DGSEM, cache) @unpack derivative_matrix, weights = dg.basis diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl index c9ab7c478a8..822ab2f87ec 100644 --- a/src/callbacks_step/stepsize_dg3d.jl +++ b/src/callbacks_step/stepsize_dg3d.jl @@ -44,7 +44,7 @@ function max_dt(u, t, mesh::TreeMesh{3}, return 2 / (nnodes(dg) * max_scaled_speed) end -function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, +function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, 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 @@ -82,7 +82,7 @@ function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, return 2 / (nnodes(dg) * max_scaled_speed) end -function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, +function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, constant_speed::True, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 13edcc29711..c9665a22af9 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -115,19 +115,49 @@ Non-periodic boundaries will be called ':x_neg', ':x_pos', ':y_neg', ':y_pos', ' - '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. +- `mapping`: a function of `NDIMS` variables to describe the mapping that transforms + the reference mesh (`[-1, 1]^n`) to the physical domain. + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. +- `faces::NTuple{2*NDIMS}`: a tuple of `2 * NDIMS` functions that describe the faces of the domain. + Each function must take `NDIMS-1` arguments. + `faces[1]` describes the face onto which the face in negative x-direction + of the unit hypercube is mapped. The face in positive x-direction of + the unit hypercube will be mapped onto the face described by `faces[2]`. + `faces[3:4]` describe the faces in positive and negative y-direction respectively + (in 2D and 3D). + `faces[5:6]` describe the faces in positive and negative z-direction respectively (in 3D). + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. +- `coordinates_min`: vector or tuple of the coordinates of the corner in the negative direction of each dimension + to create a rectangular mesh. + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. +- `coordinates_max`: vector or tuple of the coordinates of the corner in the positive direction of each dimension + to create a rectangular mesh. + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. - '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) +function T8codeMesh(trees_per_dimension; polydeg = 1, + mapping = nothing, faces = nothing, coordinates_min = nothing, + coordinates_max = nothing, + RealT = Float64, initial_refinement_level = 0, + periodicity = true) + @assert ((coordinates_min === nothing)===(coordinates_max === nothing)) "Either both or none of coordinates_min and coordinates_max must be specified" + + @assert count(i -> i !== nothing, + (mapping, faces, coordinates_min))==1 "Exactly one of mapping, faces and coordinates_min/max must be specified" + + # Extract mapping + if faces !== nothing + validate_faces(faces) + mapping = transfinite_mapping(faces) + elseif coordinates_min !== nothing + mapping = coordinates2mapping(coordinates_min, coordinates_max) + end - @assert NDIMS == 2 # Only support for NDIMS = 2 yet. + NDIMS = length(trees_per_dimension) + @assert (NDIMS == 2||NDIMS == 3) "NDIMS should be 2 or 3." # Convert periodicity to a Tuple of a Bool for every dimension if all(periodicity) @@ -141,10 +171,18 @@ function T8codeMesh(trees_per_dimension; polydeg, 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) + if NDIMS == 2 + conn = T8code.Libt8.p4est_connectivity_new_brick(trees_per_dimension..., + periodicity...) + cmesh = t8_cmesh_new_from_p4est(conn, mpi_comm(), do_partition) + T8code.Libt8.p4est_connectivity_destroy(conn) + elseif NDIMS == 3 + conn = T8code.Libt8.p8est_connectivity_new_brick(trees_per_dimension..., + periodicity...) + cmesh = t8_cmesh_new_from_p8est(conn, mpi_comm(), do_partition) + T8code.Libt8.p8est_connectivity_destroy(conn) + end scheme = t8_scheme_new_default_cxx() forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) @@ -156,28 +194,48 @@ function T8codeMesh(trees_per_dimension; polydeg, 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] + # Get cell length in reference mesh: Omega_ref = [-1,1]^NDIMS. + dx = [2 / n for n in trees_per_dimension] num_local_trees = t8_cmesh_get_num_local_trees(cmesh) # Non-periodic boundaries. boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) + if mapping === nothing + mapping_ = coordinates2mapping(ntuple(_ -> -1.0, NDIMS), ntuple(_ -> 1.0, NDIMS)) + else + mapping_ = mapping + end + 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) + if NDIMS == 2 + cell_x_offset = (verts[1, 1] - 0.5 * (trees_per_dimension[1] - 1)) * dx[1] + cell_y_offset = (verts[2, 1] - 0.5 * (trees_per_dimension[2] - 1)) * dx[2] + + for j in eachindex(nodes), i in eachindex(nodes) + tree_node_coordinates[:, i, j, itree] .= mapping_(cell_x_offset + + dx[1] * nodes[i] / 2, + cell_y_offset + + dx[2] * nodes[j] / 2) + end + elseif NDIMS == 3 + cell_x_offset = (verts[1, 1] - 0.5 * (trees_per_dimension[1] - 1)) * dx[1] + cell_y_offset = (verts[2, 1] - 0.5 * (trees_per_dimension[2] - 1)) * dx[2] + cell_z_offset = (verts[3, 1] - 0.5 * (trees_per_dimension[3] - 1)) * dx[3] + + for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) + tree_node_coordinates[:, i, j, k, itree] .= mapping_(cell_x_offset + + dx[1] * nodes[i] / 2, + cell_y_offset + + dx[2] * nodes[j] / 2, + cell_z_offset + + dx[3] * nodes[k] / 2) + end end if !periodicity[1] @@ -189,6 +247,13 @@ function T8codeMesh(trees_per_dimension; polydeg, boundary_names[3, itree] = :y_neg boundary_names[4, itree] = :y_pos end + + if NDIMS > 2 + if !periodicity[3] + boundary_names[5, itree] = :z_neg + boundary_names[6, itree] = :z_pos + end + end end return T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, @@ -196,9 +261,9 @@ function T8codeMesh(trees_per_dimension; polydeg, end """ - T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}, - mapping=nothing, polydeg=1, RealT=Float64, - initial_refinement_level=0) + T8codeMesh(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. @@ -215,10 +280,15 @@ conforming mesh from a `t8_cmesh` data structure. - `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. +function T8codeMesh(cmesh::Ptr{t8_cmesh}; + mapping = nothing, polydeg = 1, RealT = Float64, + initial_refinement_level = 0) + @assert (t8_cmesh_get_num_trees(cmesh)>0) "Given `cmesh` does not contain any trees." + + # Infer NDIMS from the geometry of the first tree. + NDIMS = Int(t8_geom_get_dimension(t8_cmesh_get_tree_geometry(cmesh, 0))) + + @assert (NDIMS == 2||NDIMS == 3) "NDIMS should be 2 or 3." scheme = t8_scheme_new_default_cxx() forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) @@ -234,34 +304,84 @@ function T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}; 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)) + if NDIMS == 2 + data_in = Array{RealT, 3}(undef, 2, 2, 2) + tmp1 = zeros(RealT, 2, length(nodes), length(nodes_in)) + verts = zeros(3, 4) - u = verts[:, 2] - verts[:, 1] - v = verts[:, 3] - verts[:, 1] - w = [0.0, 0.0, 1.0] + for itree in 0:(num_local_trees - 1) + veptr = t8_cmesh_get_tree_vertices(cmesh, itree) - vol = dot(cross(u, v), w) + # Note, `verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS))` + # sometimes does not work since `veptr` is not necessarily properly + # aligned to 8 bytes. + for icorner in 1:4 + verts[1, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 1) + verts[2, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 2) + end - if vol < 0.0 - @warn "Discovered negative volumes in `cmesh`: vol = $vol" + # Check if tree's node ordering is right-handed or print a warning. + let z = zero(eltype(verts)), o = one(eltype(verts)) + u = verts[:, 2] - verts[:, 1] + v = verts[:, 3] - verts[:, 1] + w = [z, z, o] + + # Triple product gives signed volume of spanned parallelepiped. + vol = dot(cross(u, v), w) + + if vol < z + @warn "Discovered negative volumes in `cmesh`: vol = $vol" + end + 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 - # 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) + elseif NDIMS == 3 + data_in = Array{RealT, 4}(undef, 3, 2, 2, 2) + tmp1 = zeros(RealT, 3, length(nodes), length(nodes_in), length(nodes_in)) + verts = zeros(3, 8) + + for itree in 0:(num_local_trees - 1) + veptr = t8_cmesh_get_tree_vertices(cmesh, itree) + + # Note, `verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS))` + # sometimes does not work since `veptr` is not necessarily properly + # aligned to 8 bytes. + for icorner in 1:8 + verts[1, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 1) + verts[2, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 2) + verts[3, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 3) + end + + # Tree vertices are stored in z-order. + @views data_in[:, 1, 1, 1] .= verts[1:3, 1] + @views data_in[:, 2, 1, 1] .= verts[1:3, 2] + @views data_in[:, 1, 2, 1] .= verts[1:3, 3] + @views data_in[:, 2, 2, 1] .= verts[1:3, 4] + + @views data_in[:, 1, 1, 2] .= verts[1:3, 5] + @views data_in[:, 2, 1, 2] .= verts[1:3, 6] + @views data_in[:, 1, 2, 2] .= verts[1:3, 7] + @views data_in[:, 2, 2, 2] .= verts[1:3, 8] + + # Interpolate corner coordinates to specified nodes. + multiply_dimensionwise!(view(tree_node_coordinates, :, :, :, :, itree + 1), + matrix, matrix, matrix, + data_in, + tmp1) + end end map_node_coordinates!(tree_node_coordinates, mapping) @@ -274,9 +394,7 @@ function T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}; end """ - T8codeMesh{NDIMS}(conn::Ptr{p4est_connectivity}, - mapping=nothing, polydeg=1, RealT=Float64, - initial_refinement_level=0) + T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...) Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming mesh from a `p4est_connectivity` data structure. @@ -293,24 +411,45 @@ conforming mesh from a `p4est_connectivity` data structure. - `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. - +function T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...) cmesh = t8_cmesh_new_from_p4est(conn, mpi_comm(), 0) - return T8codeMesh{NDIMS}(cmesh; kwargs...) + return T8codeMesh(cmesh; kwargs...) +end + +""" + T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) + +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(conn::Ptr{p8est_connectivity}; kwargs...) + cmesh = t8_cmesh_new_from_p8est(conn, mpi_comm(), 0) + + return T8codeMesh(cmesh; kwargs...) end """ - T8codeMesh{NDIMS}(meshfile::String; - mapping=nothing, polydeg=1, RealT=Float64, - initial_refinement_level=0) + T8codeMesh{NDIMS}(meshfile::String; kwargs...) 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. +- `ndims`: Mesh file dimension: `2` or `3`. - `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. @@ -321,17 +460,130 @@ mesh from a Gmsh mesh file (`.msh`). - `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. +function T8codeMesh(meshfile::String, ndims; kwargs...) # 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) + cmesh = t8_cmesh_from_msh_file(meshfile_prefix, 0, mpi_comm(), ndims, 0, 0) + + return T8codeMesh(cmesh; kwargs...) +end + +struct adapt_callback_passthrough + adapt_callback::Function + user_data::Any +end - return T8codeMesh{NDIMS}(cmesh; kwargs...) +# 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_wrapper(forest, + forest_from, + which_tree, + lelement_id, + ts, + is_family, + num_elements, + elements_ptr)::Cint + passthrough = unsafe_pointer_to_objref(t8_forest_get_user_data(forest))[] + + elements = unsafe_wrap(Array, elements_ptr, num_elements) + + return passthrough.adapt_callback(forest_from, which_tree, ts, lelement_id, elements, + Bool(is_family), passthrough.user_data) +end + +""" + Trixi.adapt!(mesh::T8codeMesh, adapt_callback; kwargs...) + +Adapt a `T8codeMesh` according to a user-defined `adapt_callback`. + +# Arguments +- `mesh::T8codeMesh`: Initialized mesh object. +- `adapt_callback`: A user-defined callback which tells the adaption routines + if an element should be refined, coarsened or stay unchanged. + + The expected callback signature is as follows: + + `adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, user_data)` + # Arguments + - `forest`: Pointer to the analyzed forest. + - `ltreeid`: Local index of the current tree where the analyzed elements are part of. + - `eclass_scheme`: Element class of `elements`. + - `lelemntid`: Local index of the first element in `elements`. + - `elements`: Array of elements. If consecutive elements form a family + they are passed together, otherwise `elements` consists of just one element. + - `is_family`: Boolean signifying if `elements` represents a family or not. + - `user_data`: Void pointer to some arbitrary user data. Default value is `C_NULL`. + # Returns + -1 : Coarsen family of elements. + 0 : Stay unchanged. + 1 : Refine element. + +- `kwargs`: + - `recursive = true`: Adapt the forest recursively. If true the caller must ensure that the callback + returns 0 for every analyzed element at some point to stop the recursion. + - `balance = true`: Make sure the adapted forest is 2^(NDIMS-1):1 balanced. + - `partition = true`: Partition the forest to redistribute elements evenly among MPI ranks. + - `ghost = true`: Create a ghost layer for MPI data exchange. + - `user_data = C_NULL`: Pointer to some arbitrary user-defined data. +""" +function adapt!(mesh::T8codeMesh, adapt_callback; recursive = true, balance = true, + partition = true, ghost = true, user_data = C_NULL) + # Check that forest is a committed, that is valid and usable, forest. + @assert t8_forest_is_committed(mesh.forest) != 0 + + # Init new forest. + new_forest_ref = Ref{t8_forest_t}() + 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, set_for_coarsening = 0, no_repartition = !partition + t8_forest_set_user_data(new_forest, + pointer_from_objref(Ref(adapt_callback_passthrough(adapt_callback, + user_data)))) + t8_forest_set_adapt(new_forest, mesh.forest, + @t8_adapt_callback(adapt_callback_wrapper), + recursive) + if balance + t8_forest_set_balance(new_forest, set_from, no_repartition) + end + + if partition + t8_forest_set_partition(new_forest, set_from, set_for_coarsening) + end + + t8_forest_set_ghost(new_forest, ghost, T8_GHOST_FACES) # Note: MPI support not available yet so it is a dummy call. + + # The old forest is destroyed here. + # Call `t8_forest_ref(Ref(mesh.forest))` to keep it. + t8_forest_commit(new_forest) + end + + mesh.forest = new_forest + + return nothing end # TODO: Just a placeholder. Will be implemented later when MPI is supported. @@ -343,3 +595,10 @@ end function partition!(mesh::T8codeMesh; allow_coarsening = true, weight_fn = C_NULL) return nothing end + +#! format: off +@deprecate T8codeMesh{2}(conn::Ptr{p4est_connectivity}; kwargs...) T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...) +@deprecate T8codeMesh{3}(conn::Ptr{p8est_connectivity}; kwargs...) T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) +@deprecate T8codeMesh{2}(meshfile::String; kwargs...) T8codeMesh(meshfile::String, 2; kwargs...) +@deprecate T8codeMesh{3}(meshfile::String; kwargs...) T8codeMesh(meshfile::String, 3; kwargs...) +#! format: on diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index e9994fe4569..7e383924ba7 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -6,7 +6,8 @@ #! format: noindent # Initialize data structures in element container -function init_elements!(elements, mesh::P4estMesh{3}, basis::LobattoLegendreBasis) +function init_elements!(elements, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + basis::LobattoLegendreBasis) @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -26,7 +27,7 @@ end # Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis function calc_node_coordinates!(node_coordinates, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, 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_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index 4c0845ba9af..5b3c5ae5ca8 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.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{3}, equations, mortar_l2::LobattoLegendreMortarL2, - uEltype) +function create_cache(mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, + mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal compare performance of different types fstar_threaded = [Array{uEltype, 4}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2), 4) @@ -88,7 +88,7 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral, dg::DG) @unpack interfaces = cache index_range = eachnode(dg) @@ -163,7 +163,7 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -244,7 +244,7 @@ end # Inlined function for interface flux computation for conservative flux terms @inline function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -271,7 +271,7 @@ end # Inlined function for interface flux computation for flux + nonconservative terms @inline function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -314,7 +314,7 @@ end end function prolong2boundaries!(cache, u, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral, dg::DG) @unpack boundaries = cache index_range = eachnode(dg) @@ -355,7 +355,7 @@ function prolong2boundaries!(cache, u, end function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral, dg::DG) @unpack boundaries = cache @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements @@ -417,7 +417,7 @@ function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, end function prolong2mortars!(cache, u, - mesh::P4estMesh{3}, equations, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DGSEM) @unpack fstar_tmp_threaded = cache @@ -521,7 +521,7 @@ function prolong2mortars!(cache, u, end function calc_mortar_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -595,7 +595,7 @@ end # Inlined version of the mortar flux computation on small elements for conservation fluxes @inline function calc_mortar_flux!(fstar, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -616,7 +616,7 @@ end # Inlined version of the mortar flux computation on small elements for conservation fluxes # with nonconservative terms @inline function calc_mortar_flux!(fstar, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -643,7 +643,8 @@ end end @inline function mortar_fluxes_to_elements!(surface_flux_values, - mesh::P4estMesh{3}, equations, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, mortar, fstar, u_buffer, fstar_tmp) @@ -727,7 +728,7 @@ end end function calc_surface_integral!(du, u, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index cdb085e9008..1df9f408895 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -58,7 +58,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 =# @inline function weak_form_kernel!(du, u, element, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -115,7 +116,8 @@ end # the physical fluxes in each Cartesian direction @inline function flux_differencing_kernel!(du, u, element, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -189,7 +191,8 @@ end @inline function flux_differencing_kernel!(du, u, element, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::True, equations, volume_flux, dg::DGSEM, cache, alpha = true) @unpack derivative_split = dg.basis @@ -274,7 +277,8 @@ end # [arXiv: 2008.12044v2](https://arxiv.org/pdf/2008.12044) @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::False, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @@ -369,7 +373,8 @@ end # # Calculate the finite volume fluxes inside curvilinear elements (**with non-conservative terms**). @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::True, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @@ -783,7 +788,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, end function apply_jacobian!(du, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, dg::DG, cache) @threaded for element in eachelement(dg, cache) for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) diff --git a/src/solvers/dgsem_t8code/containers_2d.jl b/src/solvers/dgsem_t8code/containers_2d.jl index 029e6674afb..bf77826a34b 100644 --- a/src/solvers/dgsem_t8code/containers_2d.jl +++ b/src/solvers/dgsem_t8code/containers_2d.jl @@ -1,3 +1,7 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. @muladd begin #! format: noindent @@ -27,6 +31,9 @@ function calc_node_coordinates!(node_coordinates, element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) element_level = t8_element_level(eclass_scheme, element) + # Note, `t8_quad_len` is encoded as an integer (Morton encoding) in + # relation to `t8_quad_root_len`. This line transforms the + # "integer" length to a float in relation to the unit interval [0,1]. element_length = t8_quad_len(element_level) / t8_quad_root_len element_coords = Array{Float64}(undef, 3) @@ -55,4 +62,16 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end + +function init_mortar_neighbor_ids!(mortars::P4estMortarContainer{2}, my_face, + other_face, orientation, neighbor_ielements, + mortar_id) + if orientation == 0 + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1 + else + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1 + end +end end # @muladd diff --git a/src/solvers/dgsem_t8code/containers_3d.jl b/src/solvers/dgsem_t8code/containers_3d.jl new file mode 100644 index 00000000000..f2d54ff07da --- /dev/null +++ b/src/solvers/dgsem_t8code/containers_3d.jl @@ -0,0 +1,235 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# Interpolate tree_node_coordinates to each quadrant at the specified nodes +function calc_node_coordinates!(node_coordinates, + mesh::T8codeMesh{3}, + 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(3), static_length(nodes), static_length(mesh.nodes), + static_length(mesh.nodes)) + matrix1 = StrideArray(undef, real(mesh), + static_length(nodes), static_length(mesh.nodes)) + matrix2 = similar(matrix1) + matrix3 = 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) + + # Note, `t8_hex_len` is encoded as an integer (Morton encoding) in + # relation to `t8_hex_root_len`. This line transforms the + # "integer" length to a float in relation to the unit interval [0,1]. + element_length = t8_hex_len(element_level) / t8_hex_root_len + + element_coords = Vector{Float64}(undef, 3) + t8_element_vertex_reference_coords(eclass_scheme, element, 0, + pointer(element_coords)) + + nodes_out_x = (2 * + (element_length * 0.5 * (nodes .+ 1) .+ element_coords[1]) .- + 1) + nodes_out_y = (2 * + (element_length * 0.5 * (nodes .+ 1) .+ element_coords[2]) .- + 1) + nodes_out_z = (2 * + (element_length * 0.5 * (nodes .+ 1) .+ element_coords[3]) .- + 1) + + polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, + baryweights_in) + polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, + baryweights_in) + polynomial_interpolation_matrix!(matrix3, mesh.nodes, nodes_out_z, + baryweights_in) + + multiply_dimensionwise!(view(node_coordinates, :, :, :, :, + current_index += 1), + matrix1, matrix2, matrix3, + view(mesh.tree_node_coordinates, :, :, :, :, + itree + 1), + tmp1) + end + end + + return node_coordinates +end + +# This routine was copied and adapted from `src/dgsem_p4est/containers_3d.jl`: `orientation_to_indices_p4est`. +function init_mortar_neighbor_ids!(mortars::P4estMortarContainer{3}, my_face, + other_face, orientation, neighbor_ielements, + mortar_id) + # my_face and other_face are the face directions (zero-based) + # of "my side" and "other side" respectively. + # Face corner 0 of the face with the lower face direction connects to a corner of the other face. + # The number of this corner is the orientation code in `p4est`. + lower = my_face <= other_face + + # x_pos, y_neg, and z_pos are the directions in which the face has right-handed coordinates + # when looked at from the outside. + my_right_handed = my_face in (1, 2, 5) + other_right_handed = other_face in (1, 2, 5) + + # If both or none are right-handed when looked at from the outside, they will have different + # orientations when looked at from the same side of the interface. + flipped = my_right_handed == other_right_handed + + # In the following illustrations, the face corner numbering of `p4est` is shown. + # ξ and η are the local coordinates of the respective face. + # We're looking at both faces from the same side of the interface, so that "other side" + # (in the illustrations on the left) has right-handed coordinates. + if !flipped + if orientation == 0 + # Corner 0 of other side matches corner 0 of my side + # 2┌──────┐3 2┌──────┐3 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 0└──────┘1 + # η η + # ↑ ↑ + # │ │ + # └───> ξ └───> ξ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[4] + 1 + + elseif ((lower && orientation == 2) # Corner 0 of my side matches corner 2 of other side + || + (!lower && orientation == 1)) # Corner 0 of other side matches corner 1 of my side + # 2┌──────┐3 0┌──────┐2 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 1└──────┘3 + # η ┌───> η + # ↑ │ + # │ ↓ + # └───> ξ ξ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[3] + 1 + + elseif ((lower && orientation == 1) # Corner 0 of my side matches corner 1 of other side + || + (!lower && orientation == 2)) # Corner 0 of other side matches corner 2 of my side + # 2┌──────┐3 3┌──────┐1 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 2└──────┘0 + # η ξ + # ↑ ↑ + # │ │ + # └───> ξ η <───┘ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[2] + 1 + + else # orientation == 3 + # Corner 0 of my side matches corner 3 of other side and + # corner 0 of other side matches corner 3 of my side. + # 2┌──────┐3 1┌──────┐0 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 3└──────┘2 + # η ξ <───┐ + # ↑ │ + # │ ↓ + # └───> ξ η + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[1] + 1 + end + else # flipped + if orientation == 0 + # Corner 0 of other side matches corner 0 of my side + # 2┌──────┐3 1┌──────┐3 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 0└──────┘2 + # η ξ + # ↑ ↑ + # │ │ + # └───> ξ └───> η + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[4] + 1 + + elseif orientation == 2 + # Corner 0 of my side matches corner 2 of other side and + # corner 0 of other side matches corner 2 of my side. + # 2┌──────┐3 0┌──────┐1 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 2└──────┘3 + # η ┌───> ξ + # ↑ │ + # │ ↓ + # └───> ξ η + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[2] + 1 + + elseif orientation == 1 + # Corner 0 of my side matches corner 1 of other side and + # corner 0 of other side matches corner 1 of my side. + # 2┌──────┐3 3┌──────┐2 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 1└──────┘0 + # η η + # ↑ ↑ + # │ │ + # └───> ξ ξ <───┘ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[3] + 1 + + else # orientation == 3 + # Corner 0 of my side matches corner 3 of other side and + # corner 0 of other side matches corner 3 of my side. + # 2┌──────┐3 2┌──────┐0 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 3└──────┘1 + # η η <───┐ + # ↑ │ + # │ ↓ + # └───> ξ ξ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[1] + 1 + end + end +end +end # @muladd diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index 16a9d7d35b1..6e9660c917d 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -28,4 +28,5 @@ end include("containers.jl") include("containers_2d.jl") +include("containers_3d.jl") end # @muladd diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 0955dc38655..02ff338e912 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -36,13 +36,15 @@ 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{3}, StructuredMesh{3}, P4estMesh{3}}, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype) NamedTuple() end -function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -79,8 +81,8 @@ function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, fstar2_L_threaded, fstar2_R_threaded, fstar3_L_threaded, fstar3_R_threaded) end -function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, - equations, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG, uEltype) A4dp1_x = Array{uEltype, 4} @@ -112,7 +114,8 @@ 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{3}, StructuredMesh{3}, P4estMesh{3}}, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal compare performance of different types A3d = Array{uEltype, 3} @@ -140,7 +143,7 @@ end # TODO: Taal discuss/refactor timer, allowing users to pass a custom timer? function rhs!(du, u, t, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, + mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -209,8 +212,8 @@ function rhs!(du, u, t, end function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, - P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -264,8 +267,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 end function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, - P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) @@ -378,8 +381,8 @@ end # TODO: Taal dimension agnostic function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, - P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DGSEM, cache) @@ -437,7 +440,8 @@ function calc_volume_integral!(du, u, end @inline function fv_kernel!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, 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, fstar3_L_threaded, fstar3_R_threaded = cache diff --git a/src/solvers/dgsem_tree/indicators_3d.jl b/src/solvers/dgsem_tree/indicators_3d.jl index 40362889397..a11a8e06e4b 100644 --- a/src/solvers/dgsem_tree/indicators_3d.jl +++ b/src/solvers/dgsem_tree/indicators_3d.jl @@ -101,7 +101,8 @@ end alpha[element] = min(alpha_max, alpha_element) end -function apply_smoothing!(mesh::Union{TreeMesh{3}, P4estMesh{3}}, alpha, alpha_tmp, dg, +function apply_smoothing!(mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, alpha, + alpha_tmp, dg, cache) # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha diff --git a/test/runtests.jl b/test/runtests.jl index 7e195fe7402..49f0977bb70 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -89,6 +89,10 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) include("test_t8code_2d.jl") end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "t8code_part2" + include("test_t8code_3d.jl") + end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "unstructured_dgmulti" include("test_unstructured_2d.jl") include("test_dgmulti_1d.jl") diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index b3e19471323..ab95e068d02 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -30,7 +30,21 @@ mkdir(outdir) end end +@trixi_testset "test check_for_negative_volumes" begin + @test_warn "Discovered negative volumes" begin + # Unstructured mesh with six cells which have left-handed node ordering. + mesh_file = joinpath(EXAMPLES_DIR, "rectangle_with_negative_volumes.msh") + isfile(mesh_file) || + download("https://gist.githubusercontent.com/jmark/bfe0d45f8e369298d6cc637733819013/raw/cecf86edecc736e8b3e06e354c494b2052d41f7a/rectangle_with_negative_volumes.msh", + mesh_file) + + # This call should throw a warning about negative volumes detected. + mesh = T8codeMesh(mesh_file, 2) + end +end + @trixi_testset "elixir_advection_basic.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), # Expected errors are exactly the same as with TreeMesh! l2=[8.311947673061856e-6], @@ -46,6 +60,7 @@ end end @trixi_testset "elixir_advection_nonconforming_flag.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming_flag.jl"), l2=[3.198940059144588e-5], @@ -61,6 +76,7 @@ end end @trixi_testset "elixir_advection_unstructured_flag.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_flag.jl"), l2=[0.0005379687442422346], linf=[0.007438525029884735]) @@ -91,6 +107,7 @@ end end @trixi_testset "elixir_advection_amr_solution_independent.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_solution_independent.jl"), # Expected errors are exactly the same as with StructuredMesh! @@ -108,6 +125,7 @@ end end @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), l2=[ @@ -133,6 +151,7 @@ end end @trixi_testset "elixir_euler_free_stream.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), l2=[ 2.063350241405049e-15, @@ -153,6 +172,7 @@ end end @trixi_testset "elixir_euler_shockcapturing_ec.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_ec.jl"), l2=[ 9.53984675e-02, @@ -178,6 +198,8 @@ end end @trixi_testset "elixir_euler_sedov.jl" begin + # This test is identical to the one in `test_p4est_2d.jl` besides minor + # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[ 3.76149952e-01, @@ -203,6 +225,7 @@ end end @trixi_testset "elixir_shallowwater_source_terms.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ 9.168126407325352e-5, @@ -228,6 +251,7 @@ end end @trixi_testset "elixir_mhd_alfven_wave.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), l2=[1.0513414461545583e-5, 1.0517900957166411e-6, 1.0517900957304043e-6, 1.511816606372376e-6, @@ -250,6 +274,8 @@ end end @trixi_testset "elixir_mhd_rotor.jl" begin + # This test is identical to the one in `test_p4est_2d.jl` besides minor + # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), l2=[0.44211360369891683, 0.8805178316216257, 0.8262710688468049, 0.0, diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl new file mode 100644 index 00000000000..4232cf04094 --- /dev/null +++ b/test/test_t8code_3d.jl @@ -0,0 +1,279 @@ +module TestExamplesT8codeMesh3D + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = joinpath(examples_dir(), "t8code_3d_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 "T8codeMesh3D" begin + # This test is identical to the one in `test_p4est_3d.jl`. + @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=[0.00016263963870641478], + linf=[0.0014537194925779984]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_advection_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_unstructured_curved.jl"), + l2=[0.0004750004258546538], + linf=[0.026527551737137167]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_advection_nonconforming.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"), + l2=[0.00253595715323843], + linf=[0.016486952252155795]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl` besides minor + # deviations from the expected error norms. + @trixi_testset "elixir_advection_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[1.1302812803902801e-5], + linf=[0.0007889950196294793], + coverage_override=(maxiters = 6, initial_refinement_level = 1, + base_level = 1, med_level = 2, max_level = 3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl` besides minor + # deviations from the expected error norms. + @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_amr_unstructured_curved.jl"), + l2=[2.0556575425846923e-5], + linf=[0.00105682693484822], + tspan=(0.0, 1.0), + coverage_override=(maxiters = 6, initial_refinement_level = 0, + base_level = 0, med_level = 1, max_level = 2)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonconforming_unstructured_curved.jl"), + l2=[ + 4.070355207909268e-5, + 4.4993257426833716e-5, + 5.10588457841744e-5, + 5.102840924036687e-5, + 0.00019986264001630542, + ], + linf=[ + 0.0016987332417202072, + 0.003622956808262634, + 0.002029576258317789, + 0.0024206977281964193, + 0.008526972236273522, + ], + tspan=(0.0, 0.01)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonperiodic.jl"), + l2=[ + 0.0015106060984283647, + 0.0014733349038567685, + 0.00147333490385685, + 0.001473334903856929, + 0.0028149479453087093, + ], + linf=[ + 0.008070806335238156, + 0.009007245083113125, + 0.009007245083121784, + 0.009007245083102688, + 0.01562861968368434, + ], + tspan=(0.0, 1.0)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_free_stream.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), + l2=[ + 5.162664597942288e-15, + 1.941857343642486e-14, + 2.0232366394187278e-14, + 2.3381518645408552e-14, + 7.083114561232324e-14, + ], + linf=[ + 7.269740365245525e-13, + 3.289868377720495e-12, + 4.440087186807773e-12, + 3.8686831516088205e-12, + 9.412914891981927e-12, + ], + tspan=(0.0, 0.03)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_free_stream_extruded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"), + l2=[ + 8.444868392439035e-16, + 4.889826056731442e-15, + 2.2921260987087585e-15, + 4.268460455702414e-15, + 1.1356712092620279e-14, + ], + linf=[ + 7.749356711883593e-14, + 2.8792246364872653e-13, + 1.1121659149182506e-13, + 3.3228975127030935e-13, + 9.592326932761353e-13, + ], + tspan=(0.0, 0.1)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), + l2=[ + 0.010380390326164493, + 0.006192950051354618, + 0.005970674274073704, + 0.005965831290564327, + 0.02628875593094754, + ], + linf=[ + 0.3326911600075694, + 0.2824952141320467, + 0.41401037398065543, + 0.45574161423218573, + 0.8099577682187109, + ], + tspan=(0.0, 0.2), + coverage_override=(polydeg = 3,)) # Prevent long compile time in CI + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + # This test is identical to the one in `test_p4est_3d.jl` besides minor + # deviations in the expected error norms. + @trixi_testset "elixir_euler_sedov.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), + l2=[ + 7.82070951e-02, + 4.33260474e-02, + 4.33260474e-02, + 4.33260474e-02, + 3.75260911e-01, + ], + linf=[ + 7.45329845e-01, + 3.21754792e-01, + 3.21754792e-01, + 3.21754792e-01, + 4.76151527e+00, + ], + tspan=(0.0, 0.3), + coverage_override=(polydeg = 3,)) # Prevent long compile time in CI + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + +# Clean up afterwards: delete Trixi.jl output directory +@test_nowarn rm(outdir, recursive = true) + +end # module