Feature: t8code as meshing backend

Project.toml

@@ -37,6 +37,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
+T8code = "d0cc0030-9a40-4274-8435-baadcfd54fa1"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3"
@@ -73,6 +74,7 @@ StaticArrays = "1"
 StrideArrays = "0.1.18"
 StructArrays = "0.6"
 SummationByPartsOperators = "0.5.25"
+T8code = "0.2.0"
 TimerOutputs = "0.5"
 Triangulate = "2.0"
 TriplotBase = "0.1"

examples/t8code_2d_dgsem/elixir_advection_basic.jl

@@ -0,0 +1,66 @@
+# The same setup as tree_2d_dgsem/elixir_advection_basic.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+
+mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
+
+trees_per_dimension = (8, 8)
+
+# Create T8codeMesh with 8 x 8 trees and 16 x 16 elements
+mesh = T8codeMesh(trees_per_dimension, polydeg=3,
+                  mapping=mapping,
+                  initial_refinement_level=1)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 1.0));
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+# save_solution = SaveSolutionCallback(interval=100,
+#                                     solution_variables=cons2prim)
+
+# 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, save_solution, stepsize_callback)
+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()

examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl

@@ -0,0 +1,115 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 4 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg=4, surface_flux=flux_lax_friedrichs)
+
+# Deformed rectangle that looks like a waving flag,
+# lower and upper faces are sinus curves, left and right are vertical lines.
+f1(s) = SVector(-1.0, s - 1.0)
+f2(s) = SVector( 1.0, s + 1.0)
+f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
+f4(s) = SVector(s,  1.0 + sin(0.5 * pi * s))
+
+faces = (f1, f2, f3, f4)
+mapping = Trixi.transfinite_mapping(faces)
+
+# Create P4estMesh with 3 x 2 trees and 6 x 4 elements,
+# approximate the geometry with a smaller polydeg for testing.
+trees_per_dimension = (3, 2)
+mesh = T8codeMesh(trees_per_dimension, polydeg=3,
+                 mapping=mapping,
+                 initial_refinement_level=1)
+
+function adapt_callback(forest,
+                        forest_from,
+                        which_tree,
+                        lelement_id,
+                        ts,
+                        is_family, 
+                        num_elements,
+                        elements_ptr) :: Cint
+
+  vertex = Vector{Cdouble}(undef,3)
+
+  elements = unsafe_wrap(Array, elements_ptr, num_elements)
+
+  Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex))
+
+  level = Trixi.t8_element_level(ts,elements[1])
+
+  # TODO: Make this condition more general.
+  if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 4
+    # return true (refine)
+    return 1
+  else
+    # return false (don't refine)
+    return 0
+  end
+end
+
+@Trixi.T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest) != 0);
+
+# Init new forest.
+new_forest_ref = Ref{Trixi.t8_forest_t}()
+Trixi.t8_forest_init(new_forest_ref);
+new_forest = new_forest_ref[]
+
+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)
+  # t8_forest_set_ghost(new_forest, 1, T8_GHOST_FACES);
+  Trixi.t8_forest_commit(new_forest)
+end
+
+mesh.forest = new_forest
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 0.2
+ode = semidiscretize(semi, (0.0, 0.2));
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+# save_solution = SaveSolutionCallback(interval=100,
+#                                      solution_variables=cons2prim)
+
+# 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, save_solution, stepsize_callback)
+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()

examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl

@@ -0,0 +1,87 @@
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the linear advection equation.
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(
+  :all => boundary_condition
+)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux.
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+
+# Deformed rectangle that looks like a waving flag,
+# lower and upper faces are sinus curves, left and right are vertical lines.
+f1(s) = SVector(-1.0, s - 1.0)
+f2(s) = SVector( 1.0, s + 1.0)
+f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
+f4(s) = SVector(s,  1.0 + sin(0.5 * pi * s))
+faces = (f1, f2, f3, f4)
+
+Trixi.validate_faces(faces)
+mapping_flag = Trixi.transfinite_mapping(faces)
+
+# Unstructured mesh with 24 cells of the square domain [-1, 1]^n.
+mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
+isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
+                              mesh_file)
+
+# INP mesh files are only support by p4est. Hence, we
+# create a p4est connecvity object first from which
+# we can create a t8code mesh.
+conn = Trixi.read_inp_p4est(mesh_file,Val(2))
+
+mesh = T8codeMesh{2}(conn, polydeg=3,
+                 mapping=mapping_flag,
+                 initial_refinement_level=2)
+
+# A semidiscretization collects data structures and functions for the spatial discretization.
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 0.2.
+tspan = (0.0, 0.2)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of
+# the simulation setup and resets the timers.
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and
+# prints the results.
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular
+# intervals.
+# save_solution = SaveSolutionCallback(interval=100,
+#                                      solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each
+# time step.
+stepsize_callback = StepsizeCallback(cfl=1.4)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to
+# the ODE solver.
+# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+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()