Feature: t8code as meshing backend (#1426)

* Initial commit for the new feature using t8code as meshing backend. * Delete t8code_2d_dgsem * Added new examples and tests. Testing updates for T8code.jl. * Worked in the comments. * Fixed spelling. * Update src/auxiliary/auxiliary.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Added whitespace in Unions. * Adapted commented out code block reporting the no. of elements per level. * Added dummy save mesh support for . * Added test . * Added to method signature. * Deleted unnecessary comments. * Removed commented out tests. * Fixed Morton ordering bug in 2D at mortar interfaces. * Disabled `save_solution` callbacks and added more tests. * Added more tests. * Updated code according to the review. * Update src/auxiliary/t8code.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/auxiliary/t8code.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/auxiliary/t8code.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/auxiliary/t8code.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/solvers/dgsem_t8code/containers_2d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update src/meshes/t8code_mesh.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Code cleanup. * Updated to [email protected] * Fixing minor issues. * Fixed typo. * Code cleanup. * Enabled `set_ghost` in examples. * Generalized type info in function signature. * Added namespace qualifier. * Updated comments. * Refactored code and deleted lots of it. * Removed a copy operation. * Fixed some merging issues and formatting. * Fixed spelling. * Fixed spelling and changed assert macro. * Applied automatic formatting. * Backup. * Removed superfluous outer constructor for T8codeMesh. * Added return statement for consistency. * Fixed wrong indentation by autoformatter. * Added comments. * Made sure an exception is thrown. * Changed flags for sc_init for t8code initialization. * Updated formatting. * Workaround for error about calling MPI routines after MPI has been finalized. * Upped to T8code v0.4.1. * Added mpi_finailize_hook for proper memory cleanup. * Added t8code to test_threaded.jl * Added a `save_mesh_file` call in order to satisfy code coverage. * Improved finalizer logic for T8coeMesh. * Refined code. * Restructured to do blocks. * Moved save_mesh_file call to test file. * Fixed spelling error. --------- Co-authored-by: Johannes Markert <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Jul 26, 2023 · fe0e78c · fe0e78c
1 parent 53a826b
commit fe0e78c
Show file tree

Hide file tree

Showing 41 changed files with 2,767 additions and 168 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -69,6 +69,7 @@ jobs:
           - structured
           - p4est_part1
           - p4est_part2
+          - t8code_part1
           - unstructured_dgmulti
           - parabolic
           - paper_self_gravitating_gas_dynamics

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

diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl
@@ -0,0 +1,143 @@
+using OrdinaryDiffEq
+using Trixi
+
+# Define new structs inside a module to allow re-evaluating the file.
+module TrixiExtension
+
+using Trixi
+
+struct IndicatorSolutionIndependent{Cache <: NamedTuple} <: Trixi.AbstractIndicator
+    cache::Cache
+end
+
+function IndicatorSolutionIndependent(semi)
+    basis = semi.solver.basis
+    alpha = Vector{real(basis)}()
+    cache = (; semi.mesh, alpha)
+    return IndicatorSolutionIndependent{typeof(cache)}(cache)
+end
+
+function (indicator::IndicatorSolutionIndependent)(u::AbstractArray{<:Any, 4},
+                                                   mesh, equations, dg, cache;
+                                                   t, kwargs...)
+    mesh = indicator.cache.mesh
+    alpha = indicator.cache.alpha
+    resize!(alpha, nelements(dg, cache))
+
+    # Predict the theoretical center.
+    advection_velocity = (0.2, -0.7)
+    center = t .* advection_velocity
+
+    inner_distance = 1
+    outer_distance = 1.85
+
+    # Iterate over all elements.
+    for element in 1:length(alpha)
+        # Calculate periodic distance between cell and center.
+        # This requires an uncurved mesh!
+        coordinates = SVector(0.5 * (cache.elements.node_coordinates[1, 1, 1, element] +
+                               cache.elements.node_coordinates[1, end, 1, element]),
+                              0.5 * (cache.elements.node_coordinates[2, 1, 1, element] +
+                               cache.elements.node_coordinates[2, 1, end, element]))
+
+        # The geometric shape of the amr should be preserved when the base_level is increased.
+        # This is done by looking at the original coordinates of each cell.
+        cell_coordinates = original_coordinates(coordinates, 5 / 8)
+        cell_distance = periodic_distance_2d(cell_coordinates, center, 10)
+        if cell_distance < (inner_distance + outer_distance) / 2
+            cell_coordinates = original_coordinates(coordinates, 5 / 16)
+            cell_distance = periodic_distance_2d(cell_coordinates, center, 10)
+        end
+
+        # Set alpha according to cells position inside the circles.
+        target_level = (cell_distance < inner_distance) + (cell_distance < outer_distance)
+        alpha[element] = target_level / 2
+    end
+    return alpha
+end
+
+# For periodic domains, distance between two points must take into account
+# periodic extensions of the domain.
+function periodic_distance_2d(coordinates, center, domain_length)
+    dx = coordinates .- center
+    dx_shifted = abs.(dx .% domain_length)
+    dx_periodic = min.(dx_shifted, domain_length .- dx_shifted)
+    return sqrt(sum(dx_periodic .^ 2))
+end
+
+# This takes a cells coordinates and transforms them into the coordinates of a
+# parent-cell it originally refined from.  It does it so that the parent-cell
+# has given cell_length.
+function original_coordinates(coordinates, cell_length)
+    offset = coordinates .% cell_length
+    offset_sign = sign.(offset)
+    border = coordinates - offset
+    center = border + (offset_sign .* cell_length / 2)
+    return center
+end
+
+end # module TrixiExtension
+
+import .TrixiExtension
+
+###############################################################################
+# Semidiscretization of the linear advection equation.
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_gauss
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-5.0, -5.0)
+coordinates_max = (5.0, 5.0)
+
+mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
+
+trees_per_dimension = (1, 1)
+
+mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
+                  mapping = mapping,
+                  initial_refinement_level = 1)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_controller = ControllerThreeLevel(semi,
+                                      TrixiExtension.IndicatorSolutionIndependent(semi),
+                                      base_level = 4,
+                                      med_level = 5, med_threshold = 0.1,
+                                      max_level = 6, max_threshold = 0.6)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        amr_callback, stepsize_callback);
+
+###############################################################################
+# Run the simulation.
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
@@ -0,0 +1,87 @@
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the linear advection equation.
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+initial_condition = initial_condition_gauss
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:all => boundary_condition)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+# Deformed rectangle that looks like a waving flag, lower and upper faces are
+# sinus curves, left and right are vertical lines.
+f1(s) = SVector(-5.0, 5 * s - 5.0)
+f2(s) = SVector(5.0, 5 * s + 5.0)
+f3(s) = SVector(5 * s, -5.0 + 5 * sin(0.5 * pi * s))
+f4(s) = SVector(5 * s, 5.0 + 5 * sin(0.5 * pi * s))
+faces = (f1, f2, f3, f4)
+
+# This creates a mapping that transforms [-1, 1]^2 to the domain with the faces
+# defined above.  It generally doesn't work for meshes loaded from mesh files
+# because these can be meshes of arbitrary domains, but the mesh below is
+# specifically built on the domain [-1, 1]^2.
+Trixi.validate_faces(faces)
+mapping_flag = Trixi.transfinite_mapping(faces)
+
+# Unstructured mesh with 24 cells of the square domain [-1, 1]^n
+mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
+isfile(mesh_file) ||
+    download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
+             mesh_file)
+
+# INP mesh files are only support by p4est. Hence, we
+# create a p4est connecvity object first from which
+# we can create a t8code mesh.
+conn = Trixi.read_inp_p4est(mesh_file, Val(2))
+
+mesh = T8codeMesh{2}(conn, polydeg = 3,
+                     mapping = mapping_flag,
+                     initial_refinement_level = 1)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 1,
+                                      med_level = 2, med_threshold = 0.1,
+                                      max_level = 3, max_threshold = 0.6)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        amr_callback, stepsize_callback)
+
+###############################################################################
+# Run the simulation.
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+summary_callback() # print the timer summary
diff --git a/examples/t8code_2d_dgsem/elixir_advection_basic.jl b/examples/t8code_2d_dgsem/elixir_advection_basic.jl
@@ -0,0 +1,59 @@
+# The same setup as tree_2d_dgsem/elixir_advection_basic.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
+
+trees_per_dimension = (8, 8)
+
+mesh = T8codeMesh(trees_per_dimension, polydeg = 3,
+                  mapping = mapping,
+                  initial_refinement_level = 1)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 1.0));
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()