Skip to content

Commit

Permalink
Feature t8code: Extending to 3D (#1535)
Browse files Browse the repository at this point in the history
* 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.

* Initial commit.

* Fxinig minor bugs.

* Fixed minor typo.

* Added first 3d example and fixed segfault.

* Added many 3D examples and tests.

* Backup.

* Fixed merging issues.

* Adding more tests.

* Fixed some merging issues and formatting.

* Fixed spelling.

* Fixed spelling and changed assert macro.

* Applied automatic formatting.

* 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.

* Made sc_finalize optional.

* Fixed spelling.

* Cleaned up examples.

* Updated and cleaned t8code solver codes.

* Updated tests for t8code 3D code.

* Fixed spelling.

* Update elixir_euler_source_terms_nonconforming_unstructured_curved.jl

* Update elixir_euler_source_terms_nonconforming_unstructured_curved.jl

* Fixed indentation.

* Update src/solvers/dgsem_structured/dg_3d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/solvers/dgsem_t8code/containers_3d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/callbacks_step/amr_dg3d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update examples/t8code_3d_dgsem/elixir_euler_ec.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/solvers/dgsem_structured/dg_3d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/meshes/t8code_mesh.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/callbacks_step/analysis_dg3d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update examples/t8code_3d_dgsem/elixir_euler_free_stream.jl

Co-authored-by: Andrew Winters <[email protected]>

* Removed NDIMS from T8codeMesh construction in case of p4est/p8est connectivity input.

* Aligned T8codeMesh constructur with other mesh constructors.

* Update examples/t8code_3d_dgsem/elixir_euler_sedov.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update examples/t8code_3d_dgsem/elixir_euler_sedov.jl

Co-authored-by: Andrew Winters <[email protected]>

* Cleanup up.

* Added @allocated test.

* Fixed formatting.

* Applied formatter.

* Code cleanup.

* Removed unused member variable.

* Apply suggestions from code review

Co-authored-by: Daniel Doehring <[email protected]>

* suggestions from review

* fix format (strange?)

* Added comments to help interpreting the source code.

* Update src/callbacks_step/amr_dg3d.jl

Co-authored-by: Benedict <[email protected]>

* Adhered to unified mesh constructor calling scheme.

* Applied formatter.

* Switched to Float64 instead of Cdouble.

* Update src/meshes/t8code_mesh.jl

Co-authored-by: Daniel Doehring <[email protected]>

* Refactored negative volume check.

* Applied formatter.

* Fixed typo resp. bug.

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>

* add missing allocation checks

* Some refactoring.

* Deleted msh file.

* Fixed a bug.

* Code cleanup.

* Ignore gmsh files.

* Removed adapt! from global namespace.

* Added documentation.

* Added @test_warn to test.

* Applied formatter.

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Joshua Lampert <[email protected]>

* Turned @warn to @info.

* Code cleanup and added @deprecated routines in order to avoid breaking release.

* Applied formatter.

* Added formatter pragmas to avoid ugly formatting.

---------

Co-authored-by: Johannes Markert <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Andrew Winters <[email protected]>
Co-authored-by: Benedict <[email protected]>
Co-authored-by: Daniel Doehring <[email protected]>
Co-authored-by: Benedict Geihe <[email protected]>
Co-authored-by: Joshua Lampert <[email protected]>
  • Loading branch information
8 people authored Jan 19, 2024
1 parent 4bb74f8 commit 1946f9d
Show file tree
Hide file tree
Showing 43 changed files with 2,102 additions and 335 deletions.
1 change: 1 addition & 0 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ jobs:
- p4est_part1
- p4est_part2
- t8code_part1
- t8code_part2
- unstructured_dgmulti
- parabolic
- paper_self_gravitating_gas_dynamics
Expand Down
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
*.mesh
*.bson
*.inp
*.msh
**/Manifest.toml
out*/
docs/build
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 1 addition & 3 deletions examples/t8code_2d_dgsem/elixir_advection_basic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
52 changes: 15 additions & 37 deletions examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
45 changes: 9 additions & 36 deletions examples/t8code_2d_dgsem/elixir_euler_free_stream.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)))
Expand Down
4 changes: 1 addition & 3 deletions examples/t8code_2d_dgsem/elixir_euler_sedov.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 1 addition & 3 deletions examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand Down
4 changes: 1 addition & 3 deletions examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
5 changes: 2 additions & 3 deletions examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions examples/t8code_2d_dgsem/elixir_mhd_rotor.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
4 changes: 1 addition & 3 deletions examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 1946f9d

Please sign in to comment.