Feature t8code: Extending to 3D

.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

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)

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)

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

examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl

@@ -20,13 +20,23 @@ 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)
 
+# 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_dhex_t
+    x::Int32
+    y::Int32
+    level::Int8
+    # [...] # See `p4est.h` in `p4est` for more info.
+end
+
+# Refine quadrants of each tree at lower left edge to level 4.
 function adapt_callback(forest,
                         forest_from,
                         which_tree,
@@ -35,16 +45,10 @@ function adapt_callback(forest,
                         is_family,
                         num_elements,
                         elements_ptr)::Cint
-    vertex = Vector{Cdouble}(undef, 3)
-
     elements = unsafe_wrap(Array, elements_ptr, num_elements)
+    el = unsafe_load(Ptr{t8_dhex_t}(elements[1]))
 
-    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
+    if el.x == 0 && el.y == 0 && el.level < 4
         # return true (refine)
         return 1
     else
@@ -53,7 +57,7 @@ function adapt_callback(forest,
     end
 end
 
-Trixi.@T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest)!=0);
+@assert(Trixi.t8_forest_is_committed(mesh.forest)!=0);
 
 # Init new forest.
 new_forest_ref = Ref{Trixi.t8_forest_t}()

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,

examples/t8code_2d_dgsem/elixir_euler_free_stream.jl

@@ -40,9 +40,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,
-                     initial_refinement_level = 1)
+mesh = T8codeMesh(conn, polydeg = 3,
+                  mapping = mapping,
+                  initial_refinement_level = 1)
 
 function adapt_callback(forest,
                         forest_from,

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)

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)

examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl

@@ -40,9 +40,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)
 
 function adapt_callback(forest,
                         forest_from,

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,

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)

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)

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

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