Negative volume handling in Trixi.jl and t8code linkage

[sandro-elsweijer]

Hey y'all,

during our coding week we (the t8code devs and @benegee ) started to update T8code.jl and Trixi.jl to t8code version 4.0.0.
We tried to set up some simulations and we stumbled over an already known problem: Negative volumes.
Because of a bug, an overly specific implementation of the negative volume handling in t8code and our principle to not alter our users data, we decided to discontinue our negative volume handling in t8code some time ago. But now Trixi.jl simulations using our Gmsh reader crash immediately because of said negative volumes. When we disable the negative volume check in Trixi.jl, the simulation also crashes in the first timestep. Presumably because of these negative volumes.

To test this, I altered the t8code_2d_dgsem/elixir_euler_free_stream.jl slightly using this msh file. I used t8code v4.0.0, commit c7635799a08541af48db828eb8815da43aae56a8 of T8code.jl and this commit 84e6f4933d1354852321f4a1d8326686d6883bc0 of Trixi.jl.

using OrdinaryDiffEqLowStorageRK
using Trixi

###############################################################################
# Semidiscretization of the compressible Euler equations.

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_constant

# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
# `StepsizeCallback` (CFL-Condition) and less diffusion.
solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))

###############################################################################
# use a local msh file
mesh_file = ("/home/elsw_sa/workspace/wing/airfoil_0.0.msh")

# disable the mapping
mesh = T8codeMesh(mesh_file, 2; polydeg = 3,
                  initial_refinement_level = 2)

function adapt_callback(forest, ltreeid, scheme, tree_class, lelemntid, elements, is_family,
                        user_data)
    vertex = Vector{Cdouble}(undef, 3)

    Trixi.t8_element_get_vertex_reference_coords(scheme, tree_class, elements[1], 0, vertex)

    level = Trixi.t8_element_get_level(scheme, tree_class, elements[1])

    # TODO: Make this condition more general.
    if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 3
        # return true (refine)
        return 1
    else
        # return false (don't refine)
        return 0
    end
end

# disable adaptation to narrow down the problem
#Trixi.adapt!(mesh, adapt_callback)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition)))

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 100,
                                     save_initial_solution = true,
                                     save_final_solution = true,
                                     solution_variables = cons2prim)

stepsize_callback = StepsizeCallback(cfl = 2.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        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
            ode_default_options()..., callback = callbacks);

Running this gives me negative volume errors and after commenting them out it also crashes:

julia> trixi_include(joinpath("/localdata1/elsw_sa/source/Trixi.jl/run/dev/Trixi/examples", "t8code_2d_dgsem", "elixir_euler_free_stream.jl"))
[...]
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/integrator_interface.jl:618
retcode: DtNaN
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 0.037017738177397674
u: 2-element Vector{Vector{Float64}}:
 [1.0, 0.1, -0.2, 10.0, 1.0, 0.1, -0.2, 10.0, 1.0, 0.1  …  -0.2, 10.0, 1.0, 0.1, -0.2, 10.0, 1.0, 0.1, -0.2, 10.0]
 [0.999990987901522, 0.10001089881356225, -0.19997295892804484, 9.999832607391165, 0.9999999817136326, 0.10000141065568129, -0.20000029072240522, 9.999999979325006, 0.9999988725397543, 0.10000216336901167  …  -0.20000115904847085, 10.000004260436967, 0.999999069208971, 0.10000041388082619, -0.1999977688647179, 9.999986040543485, 1.0000024519825415, 0.09999445270640808, -0.2000262059108428, 10.000101920589348]

But to be honest we do not really understand what exactly negative volumes in Trixi.jl are. Visually, the vtk output of the grid looks fine and also all normals of the grid point in the same direction:

To tackle this we would like to join forces and discuss how to tackle this problem, if we implement this in t8code or Trixi.jl and also we would like to know more about this topic. Some of our developer team is already on holiday and tbh nothing will happen this year. But we want you to know about this and then we can take a look at this next year.

Have a nice Christmas and best wishes for the New Year!
The t8code dev team

[DanielDoehring]

Hm, so If you say the import of a .msh mesh into T8code is fine, then this line

Trixi.jl/src/meshes/t8code_mesh.jl

Line 673 in edde724

cmesh = t8_cmesh_from_msh_file(meshfile_prefix, 0, mpi_comm(), NDIMS, 0, 0)

should be fine, i.e., do not produce negative volumes.

Then you probably need to check

Trixi.jl/src/meshes/t8code_mesh.jl

Lines 527 to 547 in edde724

	function T8codeMesh(cmesh::Ptr{t8_cmesh};
	mapping = nothing, polydeg = 1, RealT = Float64,
	initial_refinement_level = 0)
	@assert (t8_cmesh_get_num_trees(cmesh)>0) "Given `cmesh` does not contain any trees."
	# Infer NDIMS from the geometry of the first tree.
	NDIMS = Int(t8_cmesh_get_dimension(cmesh))
	@assert (NDIMS == 2||NDIMS == 3) "NDIMS should be 2 or 3."
	do_face_ghost = mpi_isparallel()
	scheme = t8_scheme_new_default_cxx()
	forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost,
	mpi_comm())
	# There's no simple and generic way to distinguish boundaries, yet. Name all of them :all.
	boundary_names = fill(:all, 2 * NDIMS, t8_cmesh_get_num_trees(cmesh))
	return T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = polydeg,
	mapping = mapping)
	end

and

Trixi.jl/src/meshes/t8code_mesh.jl

Line 290 in edde724

cmesh = t8_forest_get_cmesh(forest)

and

Trixi.jl/src/meshes/t8code_mesh.jl

Line 306 in edde724

vertices_pointer = t8_cmesh_get_tree_vertices(cmesh, itree - 1)

if something fuzzy happens along the way.

[sandro-elsweijer]

Hey @DanielDoehring,

thank you for the quick reply!
Since the ParaView output seems fine I don't think that this is the reason.
Johannes Markert once implemented a negative volume handling in t8code, which we discontinued until v4.0.0. Now, the same check in Trixi.jl fails.
https://github.com/trixi-framework/Trixi.jl/blob/main/src/meshes/t8code_mesh.jl#L303-L338

We discontinued it, because it caused a bug and it was also overly specific just for 2D meshes in the xy plane. Now we need to handle this, but we do not really get what the requirements from Trixi.jl are.

[benegee]

Looking at https://github.com/trixi-framework/Trixi.jl/blob/main/src/meshes/t8code_mesh.jl#L316-L323, it seems to me the following assumptions are being made:

1. A 2D mesh lies in the x-y-plane in 3D space, i.e. z is always 0.
2. The local coordinate system in each quadrilateral cell is spanned by node2 - node1 (x axis) and node3 - node1 (y axis)
3. The cross product of the local spanning vectors must point in positive z direction. This is supposed to ensure a right-handed local coordinate system.

First of all, we should document these assumption.

The question remains: what makes this example fail?

[DanielDoehring]

I might be missing something here, but the negative volume computation is purely based on the vertices of the cell, i.e., before transformations to higher order occur:

Trixi.jl/src/meshes/t8code_mesh.jl

Lines 306 to 327 in edde724

	vertices_pointer = t8_cmesh_get_tree_vertices(cmesh, itree - 1)
	# Note, `vertices = unsafe_wrap(Array, vertices_pointer, (3, 1 << NDIMS))`
	# sometimes does not work since `vertices_pointer` is not necessarily properly
	# aligned to 8 bytes.
	for icorner in 1:number_of_corners
	vertices[1, icorner] = unsafe_load(vertices_pointer, (icorner - 1) * 3 + 1)
	vertices[2, icorner] = unsafe_load(vertices_pointer, (icorner - 1) * 3 + 2)
	end
	# Check if tree's node ordering is right-handed or print a warning.
	let z = zero(eltype(vertices)), o = one(eltype(vertices))
	u = vertices[:, 2] - vertices[:, 1]
	v = vertices[:, 3] - vertices[:, 1]
	w = [z, z, o]
	# Triple product gives signed volume of spanned parallelepiped.
	vol = dot(cross(u, v), w)
	if vol < z
	error("Discovered negative volumes in `cmesh`: vol = $vol")
	end

So either this volume computation is wrong or the mesh is indeed somehow distorted, right? I would try to pin this down by using other mesh reading tools, or trying to get the exact cell info where things are supposed to go wrong and then check manually.

the answer from @benegee does into the same direction 👍

[benegee]

@sandro-elsweijer's mesh lies in the x-y plane. ✔️

But the normals in the screenshot point in negative z-direction.

In case the computation of the normals in paraview works similar to Trixi's, this would mean the cells in the file do not match with Trixi's assumption of the ordering.

Something like

3---4        2---4
|   |   vs   |   | 
1---2        1---3

[DanielDoehring]

Isn't the standard node order right-handed, i.e., something similar to

4---3
|   |
1---2

?

[benegee]

Yes, that makes sense! Cf. https://trixi-framework.github.io/TrixiDocumentation/stable/meshes/unstructured_quad_mesh/

The problem would then be like

4---3        2---3
|   |   vs   |   | 
1---2        1---4

[sandro-elsweijer]

In t8code we have

3---4
|   |
1---2

, so maybe there is just a conversion step missing

[DanielDoehring]

Hm okay, that might be the reason - and is probably a source of many future evil...

[andrewwinters5000]

so maybe there is just a conversion step missing

I think that is the case. I recall having these ordering issues when creating the file format and read-in for the higher order boundary data on P4estMesh (https://trixi-framework.github.io/TrixiDocumentation/stable/meshes/p4est_mesh/) as the sides are numbered differently from the UnstructuredMesh2D (which is based on the conventions in the book by Kopriva).

[sandro-elsweijer]

Yeah that makes sense. p4est uses the same node ordering as t8code. So if that's the case it should be relatively straightforward to copy the design from the p4est interface.

so maybe there is just a conversion step missing

I think that is the case. I recall having these ordering issues when creating the file format and read-in for the higher order boundary data on P4estMesh (https://trixi-framework.github.io/TrixiDocumentation/stable/meshes/p4est_mesh/) as the sides are numbered differently from the UnstructuredMesh2D (which is based on the conventions in the book by Kopriva).