Closure orientation changes (WIP)

[connorjward]

Description

This PR contains a first stab at a big refactor of how we compute cell closures and orientations.

This work was motivated by the fact that we have totally distinct code paths for determining closures and orientations for simplices, quads and hexes. This PR is an attempt to unify them and also to make greater use of DMPlex orientations.

One thing in particular that bothered me was that entity_orientations (which contains the orientations of points in the closure relative to some canonical orientation) for simplices and quads was not zero. This doesn't make sense because we cleverly design our closures using the vertex numbering in order to ensure that a canonical ordering is always possible, and hence that entity_orientations should always be zero.

I have figured out that using DMPlexOrientPoint, coupled with a single transformation from plex closure to FIAT closure seems to get the job done. DMPlex does pretty much all the heavy lifting for us in terms of determining orientations so all we have to do is know the transformation from its canonical represent to FIAT's.

Thus far I have only done this for triangles (and lines and points). I haven't yet done the work to get other polytopes transforming correctly.

@ksagiyam I would be particularly interested in getting your feedback on this. I may be barking up totally the wrong tree.

[ksagiyam]

I think the idea is interesting, and some components might be useful. I think, though, we are to abandon global vertex number based approaches as hex requires something more general.

[ksagiyam]

To make this globally consistent, you need to use global section; see

firedrake/firedrake/mesh.py

Line 1092 in d282271

vertex_numbering = self._vertex_numbering.createGlobalSection(plex.getPointSF())

and convert negative global number, corresponding to vertices owned by the other processes, using the -(o + 1) rule.

We decided to order DoFs on each entity relative to the cone of that entity for checkpointing purpose because the ordering of the cone is guaranteed to be preserved in the save-load cycle.
Modifying cones using DMPlexOrientPoint() based on the global vertex numbers is actually not desired as global vertex numbers are not preserved in the save-load cycle (or upon redistribution of the mesh).

I think we can make the global vertex numbers (and actually all the global entity numbers) preserved upon save/load or redistribution, but it has yet to be done. Even if it has been implemented, I think, we should abandon all approaches based on the global vertex numbers as hex requires a more generic approach anyway.

[connorjward]

To make this globally consistent, you need to use global section

Sure. That should be straightforward to add.

global vertex numbers is actually not desired as global vertex numbers are not preserved in the save-load cycle

I didn't realise this. How come checkpointing works currently then? What I have here is a simplification of the prior approach so the same limitations should apply.

I think, we should abandon all approaches based on the global vertex numbers

I disagree with this. I think that we want to have this because it allows us to pack/unpack H(div) and H(curl) spaces without potentially expensive runtime transformations.

[connorjward]

Another contribution of this PR that I want to emphasise is that it unifies the code paths for cell_closure and entity_orientations. I think we can pretty much rely on the orientations from getTransitiveClosure for our own work (with a single transformation from canonical plex to canonical FIAT).

Basically PETSc tells us things like "relative to this cell, this edge is flipped". I believe that, subject to a canonical transformation, this is universal information applicable to both DMPlex and FIAT.

[connorjward]

We decided to order DoFs on each entity relative to the cone of that entity for checkpointing purpose because the ordering of the cone is guaranteed to be preserved in the save-load cycle.
Modifying cones using DMPlexOrientPoint() based on the global vertex numbers is actually not desired as global vertex numbers are not preserved in the save-load cycle (or upon redistribution of the mesh).

I have thought more on this and this makes a lot of sense. The cones should not be arbitrarily modified upon loading the mesh. However, I still think that my approach has value. Could we introduce perhaps a DMCopy of the topology just for tabulating the closure, called, say, closure_dm? The advantage of having such a DM is then we get nicely ordered closures and can preserve the orientable properties of simplices and quads.

[connorjward]

I have had a lengthy discussion with @ksagiyam and we have concluded the following:

• This work is correct and worth doing. Orienting the mesh w.r.t. the global vertex numbering is a valid thing to do even when we have a more general solution.
• For it to work with checkpointing we need to preserve the global numbering when we save the mesh. This should have no impact on the cache locality of the local numbering.
• This constitutes a change in checkpoint files. I think we should have a "firedrake checkpoint file version number" that we can inspect in order to keep backwards compatibility.
• The simplest approach is to modify (orient) the plex immediately upon instantiation. This is only required for a new plex; meshes that have been loaded from checkpoint files are already oriented.

[ksagiyam]

Actually, I/O backward compatibility requires more thoughts. If we load a mesh from an old file, cones will not be ordered in the way you want. If we then reorder the cones using whatever we use to reorder cones, then we will not be able to later load functions onto that mesh as those functions assume old cone orderings.

[connorjward]

Actually, I/O backward compatibility requires more thoughts. If we load a mesh from an old file, cones will not be ordered in the way you want. If we then reorder the cones using whatever we use to reorder cones, then we will not be able to later load functions onto that mesh as those functions assume old cone orderings.

Good point. We didn't discuss preserving saved functions when we discussed backwards compatibility. I think what will need to happen is the following:

1. Load the old mesh and function
2. Orient the old mesh
3. When we orient the old mesh - by looping over all mesh points and determining an orientation - we should stash this orientation somewhere with an IS.
4. We can now transform (permute) the old functions with this IS. The orientations in the IS tell us which permutation to apply.

[ksagiyam]

In theory, I think, we could directly fix it, but it might be cumbersome. We should note that:

• we should be able to load functions any time after we load the mesh.
• if we reorder cones, I think, it will affect RCM reordering, which changes _dm_renumbering, which then changes FunctionSpaces section, and thus anything that depends on that section; e.g.. halo.sf for halo exchange. So just permuting old function`s DoFs will not be enough.

[connorjward]

• we should be able to load functions any time after we load the mesh.

That's fine. We can store some prev_point_orientation property on the mesh and refer to it when we load the function

• if we reorder cones, I think, it will affect RCM reordering, which changes _dm_renumbering, which then changes FunctionSpaces section, and thus anything that depends on that section; e.g.. halo.sf for halo exchange. So just permuting old function`s DoFs will not be enough.

It's very hard to be sure, but I think that this might be OK and "just work". Once we've done the appropriate permutation I imagine that whatever Section and SF we have already set up would work correctly. Definitely something to remain aware of though.