Feature: DGMultiMesh from t8code forest (prototype)

[jmark]

This PR adds routines to build a DGMultiMesh from a t8code forest. It is a prototype implementation. It is not particular performant, robust and only supports conformal meshes with one element type (no hybrid meshes, no mortar interfaces).

This is a follow-up PR to #2259.

[github-actions]

Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md with its PR number.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

[github-actions]

Remaining comments which cannot be posted as a review comment to avoid GitHub Rate Limit

JuliaFormatter

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 175 in 1aa2762

end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 177 to 178 in 1aa2762

	EToV = zeros(Int, num_elements, num_corners)
	VXYZ = Tuple(Vector{Float64}() for idim = 1:ndims)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 180 to 183 in 1aa2762

	etov = zeros(Int, num_corners)
	min_range = vxyz |> minimum
	max_range = vxyz |> maximum

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 185 to 186 in 1aa2762

	num_bins = 1e9
	inverse_bin_size = num_bins / (max_range - min_range)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 188 in 1aa2762

hashed_vertices = Dict{Tuple{Tuple(UInt128 for idim = 1:ndims)...},Int}()

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 190 to 193 in 1aa2762

	cumulative_index = 0
	for ielem = 1:num_elements
	for ivert = 1:num_corners
	hashed_xyz = Tuple(trunc(UInt128, (vxyz[idim,ivert,ielem] - min_range + 0.5) * inverse_bin_size) for idim in 1:ndims)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 195 to 205 in 1aa2762

	if haskey(hashed_vertices, hashed_xyz)
	index = hashed_vertices[hashed_xyz]
	else
	index = cumulative_index += 1
	hashed_vertices[hashed_xyz] = index
	for idim = 1:ndims
	push!(VXYZ[idim], vxyz[idim,ivert,ielem])
	end
	end
	etov[ivert] = index
	end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 207 to 208 in 1aa2762

	EToV[ielem,:] = t8code2startupdg(eclass, etov)
	end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 210 to 211 in 1aa2762

	return VXYZ, EToV
	end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 213 in 1aa2762

function compute_coordinates(forest::Ptr{t8_forest}, rd::RefElemData, md::MeshData)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 215 to 278 in 1aa2762

	if rd.element_type isa StartUpDG.Tri
	e_rst = [[ 1.0, 0.0, 0.0 ],
	[ 1.0, 1.0, 0.0 ],
	[ 0.0, 0.0, 0.0 ]]
	elseif rd.element_type isa StartUpDG.Quad
	e_rst = [[ 1.0, 0.0, 0.0 ],
	[ 0.0, 1.0, 0.0 ],
	[ 0.0, 0.0, 0.0 ]]
	elseif rd.element_type isa StartUpDG.Tet
	e_rst = [[ 1.0, 0.0, 0.0 ],
	[ 1.0, 0.0, 1.0 ],
	[ 1.0, 1.0, 1.0 ]]
	elseif rd.element_type isa StartUpDG.Wedge
	e_rst = [[ 1.0, 0.0, 0.0 ],
	[ 1.0, 1.0, 0.0 ],
	[ 0.0, 0.0, 1.0 ]]
	elseif rd.element_type isa StartUpDG.Hex
	e_rst = [[ 1.0, 0.0, 0.0 ],
	[ 0.0, 1.0, 0.0 ],
	[ 0.0, 0.0, 1.0 ]]
	else
	@error element_type
	end
	ndims = length(md.xyz)
	xyz = tuple([similar(md.x) for _ = 1:ndims]...)
	num_local_trees = t8_forest_get_num_local_trees(forest)
	current_element = 1
	for itree = 0:num_local_trees-1
	tree_class = t8_forest_get_tree_class(forest, itree)
	eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class)
	eclass = t8_forest_get_eclass(forest, itree)
	num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree)
	for ielement = 0:num_elements_in_tree-1
	element = t8_forest_get_element_in_tree(forest, itree, ielement)
	num_corners = t8_element_num_corners(eclass_scheme, element)
	verts = zeros(3,num_corners)
	for ivert in 1:num_corners
	t8_forest_element_coordinate(forest, itree, element, ivert-1, @view(verts[:,ivert]))
	end
	_, do_perm = correct_negative_volume(eclass, verts)
	if do_perm
	perm = [2,1,3]
	else
	perm = [1,2,3]
	end
	for i = 1:length(rd.r)
	ref_coords = zeros(3)
	out_coords = Vector{Cdouble}(undef,3)
	for iref = 1:ndims
	rst = 0.5*(1.0 + rd.rst[perm[iref]][i])
	ref_coords += rst*e_rst[iref]
	end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 280 in 1aa2762

t8_forest_element_from_ref_coords(forest, itree, element, pointer(ref_coords), 1, pointer(out_coords))

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Lines 282 to 285 in 1aa2762

	for idim = 1:ndims
	xyz[idim][i,current_element] = out_coords[idim]
	end
	end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 287 in 1aa2762

current_element += 1

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 289 in 1aa2762

end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 291 in 1aa2762

return xyz

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/meshes/dgmulti_t8code.jl

Line 304 in 1aa2762

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/solvers/dgmulti/dg.jl

Lines 489 to 518 in 1aa2762

	num_faces_total = num_faces * md.num_elements
	# This function was originally defined as
	# `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`.
	# This results in allocations due to https://github.com/JuliaLang/julia/issues/36313.
	# To avoid allocations, we use Tim Holy's suggestion:
	# https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300.
	reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ())
	u_face_values = reshape_by_face(u_face_values)
	flux_face_values = reshape_by_face(flux_face_values)
	Jf = reshape_by_face(Jf)
	nxyzJ, xyzf = reshape_by_face.(nxyzJ), reshape_by_face.(xyzf) # broadcast over nxyzJ::NTuple{NDIMS,Matrix}
	# loop through boundary faces, which correspond to columns of reshaped u_face_values, ...
	for i in mesh.boundary_faces[boundary_key]
	for i in Base.OneTo(num_pts_per_face)
	face_normal = SVector{NDIMS}(getindex.(nxyzJ, i)) / Jf[i]
	face_coordinates = SVector{NDIMS}(getindex.(xyzf, i))
	flux_face_values[i] = boundary_condition(u_face_values[i],
	face_normal, face_coordinates,
	t,
	surface_flux, equations) *
	Jf[i]
	end
	end
	elseif mesh.boundary_faces_type == :nodes
	for i in mesh.boundary_faces[boundary_key]
	face_normal = SVector{NDIMS}(getindex.(nxyzJ, i)) / Jf[i]
	face_coordinates = SVector{NDIMS}(getindex.(xyzf, i))

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/solvers/dgmulti/dg.jl

Lines 520 to 528 in 1aa2762

	face_normal, face_coordinates,
	t,
	surface_flux, equations)
	flux_face_values[i] = temp_f * Jf[i]
	end
	else
	error("Unknown boundary_faces type.")
	end

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/solvers/dgmulti/types.jl

Line 441 in 1aa2762

function DGMultiMesh(dg::DGMulti, cmesh::Ptr{t8_cmesh}; initial_refinement_level = 0, is_on_boundary::Dict{Symbol, <:Function} = Dict())

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/solvers/dgmulti/types.jl

Line 444 in 1aa2762

forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost,

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/solvers/dgmulti/types.jl

Line 454 in 1aa2762

mesh = DGMultiMesh(dg,GeometricTermsType(Curved(), dg), md, boundary_nodes)

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/src/solvers/dgmulti/types.jl

Line 459 in 1aa2762

---

[JuliaFormatter] _{reported by reviewdog 🐶}

Trixi.jl/test/test_dgmulti_3d.jl

Line 452 in 1aa2762

[codecov]

Codecov Report

Attention: Patch coverage is 21.41280% with 356 lines in your changes missing coverage. Please review.

Project coverage is 93.46%. Comparing base (544a6f0) to head (211ced8).
Report is 83 commits behind head on main.

Files with missing lines	Patch %	Lines
src/meshes/dgmulti_t8code.jl	0.00%	171 Missing ⚠️
.../dgmulti_3d/elixir_euler_baroclinic_instability.jl	0.00%	99 Missing ⚠️
src/meshes/mesh_io.jl	51.79%	27 Missing ⚠️
src/meshes/dgmulti_meshes.jl	19.23%	21 Missing ⚠️
src/callbacks_step/save_solution_dg.jl	72.00%	14 Missing ⚠️
src/solvers/dgmulti/types.jl	0.00%	12 Missing ⚠️
src/solvers/dgmulti/dg.jl	70.83%	7 Missing ⚠️
...amples/dgmulti_3d/elixir_advection_cubed_sphere.jl	0.00%	2 Missing ⚠️
...ples/dgmulti_3d/elixir_advection_prismed_sphere.jl	0.00%	1 Missing ⚠️
src/callbacks_step/save_solution.jl	88.89%	1 Missing ⚠️
... and 1 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2270      +/-   ##
==========================================
- Coverage   96.94%   93.46%   -3.48%     
==========================================
  Files         502      506       +4     
  Lines       41682    42070     +388     
==========================================
- Hits        40405    39318    -1087     
- Misses       1277     2752    +1475     

Flag	Coverage Δ
unittests	93.46% <21.41%> (-3.48%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[tristanmontoya]

The CI hangs (or at least takes an extremely long time) at CI / unstructured_dgmulti - ubuntu-latest. Need to investigate why.