Skip to content

Closure orientation changes (WIP) #3332

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 6 commits into
base: master
Choose a base branch
from
Draft

Closure orientation changes (WIP) #3332

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
21a4e99
Appears to work...
connorjward Jan 16, 2024
7ab1192
WIP
connorjward Jan 16, 2024
af96718
WIP, think I have something neat
connorjward Jan 17, 2024
b248826
Appears to work for triangles
connorjward Jan 17, 2024
406a828
Cleanup
connorjward Jan 17, 2024
a7c969c
Cleanup
connorjward Jan 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
128 changes: 110 additions & 18 deletions firedrake/cython/dmcommon.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -341,6 +341,7 @@ cdef inline PetscInt _reorder_plex_cone(PETSc.DM dm,
comparison with the associated FIAT cones.
The reordered DMPlex cones are later used to construct cell_closure and to
compute orientations.

"""
if dm.getCellType(p) == PETSc.DM.PolytopeType.POINT:
raise RuntimeError(f"POINT has no cone")
Expand Down Expand Up @@ -392,16 +393,19 @@ cdef inline PetscInt _reorder_plex_cone(PETSc.DM dm,
raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")


cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm,
PetscInt p,
PetscInt *plex_closure,
PetscInt *fiat_closure):
cdef inline PetscInt _reorder_plex_closure(
PETSc.DM dm,
PetscInt p,
PetscInt *plex_closure,
PetscInt *fiat_closure,
PetscInt *orientations,
):
"""Reorder DMPlex closure for FIAT closure.

:arg dm: The DMPlex object
:arg p: The plex point
:arg plex_closure: The original DMPlex closure
:arg fiat_closure: The reorderd closure (output)
:arg fiat_closure: The reordered closure (output)

This function defines fixed rules to reorder DMPlex closures
for FIAT closures.
Expand All @@ -411,27 +415,56 @@ cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm,
Indeed the same FIAT closure can be obtained merely using
_reorder_plex_cone() and ensuring that the cell orientation is 0.
"""
cdef PetscDMPolytopeType ct = dm.getCellType(p)

if dm.getCellType(p) == PETSc.DM.PolytopeType.POINT:
raise RuntimeError(f"POINT has no cone")
elif dm.getCellType(p) == PETSc.DM.PolytopeType.SEGMENT:
# UFCInterval: 0---2---1
#
# PETSc.DM.PolytopeType. 1---0---2
# SEGMENT:
raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")
fiat_closure[0] = plex_closure[2*0]
orientations[0] = plex_closure[2*0+1]

fiat_closure[1] = plex_closure[2*1]
orientations[1] = plex_closure[2*1+1]
elif dm.getCellType(p) == PETSc.DM.PolytopeType.TRIANGLE:
# UFCTriangle: 1
# | \
# ↑ ↘
# 5 3
# | 6 \
# 0---4---2
# ↑ 6 ↘
# 0 → 4 → 2
#
# PETSc.DM.PolytopeType. 4
# TRIANGLE: | \
# 3 1
# | 0 \
# 6---2---5
raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")
# PETSc.DM.PolytopeType. 5
# TRIANGLE: ↑ ↘
# 1 2
# ↑ 0 ↘
# 4 ← 3 ← 6

# TODO: could straightforwardly-ish convert these directly to FIAT orientations
fiat_closure[0] = plex_closure[2*4]
orientations[0] = plex_closure[2*4+1]

fiat_closure[1] = plex_closure[2*5]
orientations[1] = plex_closure[2*5+1]

fiat_closure[2] = plex_closure[2*6]
orientations[2] = plex_closure[2*6+1]

fiat_closure[3] = plex_closure[2*2]
orientations[3] = plex_closure[2*2+1]

# this edge is flipped
fiat_closure[4] = plex_closure[2*3]
orientations[4] = DMPolytopeTypeComposeOrientation(ct, plex_closure[2*3 + 1], -1)

fiat_closure[5] = plex_closure[2 * 1]
orientations[5] = plex_closure[2 * 1+1]

fiat_closure[6] = plex_closure[2 * 0]
orientations[6] = plex_closure[2 * 0+1]

elif dm.getCellType(p) == PETSc.DM.PolytopeType.TETRAHEDRON:
# UFCTetrahedron: 0---9---1---9---0
# \ 12 / \ 13 /
Expand Down Expand Up @@ -523,6 +556,61 @@ cdef inline PetscInt _reorder_plex_closure(PETSc.DM dm,
raise NotImplementedError(f"Not implemented for {dm.getCellType(p)}")


def cell_orientation_from_vertex_arrangement(
cell_type,
vertex_ordering
):
cdef:
PetscDMPolytopeType ct
PetscInt nvert, narr
PetscInt *arr

ct = cell_type
nvert = DMPolytopeTypeGetNumVertices(ct)
narr = DMPolytopeTypeGetNumArrangments(ct)
shift = narr // 2
for o in range(-shift, shift):
arr = DMPolytopeTypeGetVertexArrangment(ct, o)

# compare to vertex_ordering
equal = True
for i in range(nvert):
if arr[i] != vertex_ordering[i]:
equal = False
break
if equal:
return o

raise ValueError(f"No valid orientation found for vertex arrangement {vertex_ordering}")


def _ordering(iterable):
iterable = list(iterable)
return tuple(iterable.index(p) for p in sorted(iterable))


def orient_mesh(PETSc.DM dm, vertex_numbering):
cdef:
PetscDMPolytopeType ct
PetscInt o, nverts

pstart, pend = dm.getChart()
for p in range(pstart, pend):
cell_type = dm.getCellType(p)

# vertices have just a single orientation
if cell_type == PETSc.DM.PolytopeType.POINT:
continue

ct = cell_type
nverts = DMPolytopeTypeGetNumVertices(ct)
verts = dm.getTransitiveClosure(p)[0][-nverts:]
verts_renum = [vertex_numbering.getOffset(v) for v in verts]
vertex_ordering = _ordering(verts_renum)
o = cell_orientation_from_vertex_arrangement(cell_type, vertex_ordering)
DMPlexOrientPoint(dm.dm, p, o)
Comment on lines +592 to +611
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.



@cython.boundscheck(False)
@cython.wraparound(False)
def create_cell_closure(PETSc.DM dm,
Expand All @@ -538,7 +626,7 @@ def create_cell_closure(PETSc.DM dm,
PetscInt c, cStart, cEnd, cell, i
PetscInt closureSize = _closureSize, closureSize1
PetscInt *closure = NULL
PetscInt *fiat_closure = NULL
PetscInt *fiat_closure = NULL, *fiat_orientations = NULL
np.ndarray[PetscInt, ndim=2, mode="c"] cell_closure

get_height_stratum(dm.dm, 0, &cStart, &cEnd)
Expand All @@ -550,16 +638,20 @@ def create_cell_closure(PETSc.DM dm,
raise RuntimeError(f"point 0 and point {c} have different cell types")
restore_transitive_closure(dm.dm, c, PETSC_TRUE, &closureSize1, &closure)
cell_closure = np.empty((cEnd - cStart, closureSize), dtype=IntType)
orientations = np.empty((cEnd - cStart, closureSize), dtype=IntType)
CHKERR(PetscMalloc1(closureSize, &fiat_closure))
CHKERR(PetscMalloc1(closureSize, &fiat_orientations))
for c in range(cStart, cEnd):
CHKERR(PetscSectionGetOffset(cell_numbering.sec, c, &cell))
get_transitive_closure(dm.dm, c, PETSC_TRUE, &closureSize1, &closure)
_reorder_plex_closure(dm, c, closure, fiat_closure)
_reorder_plex_closure(dm, c, closure, fiat_closure, fiat_orientations)
restore_transitive_closure(dm.dm, c, PETSC_TRUE, &closureSize1, &closure)
for i in range(closureSize):
cell_closure[cell, i] = fiat_closure[i]
orientations[cell, i] = fiat_orientations[i]
CHKERR(PetscFree(fiat_closure))
return cell_closure
CHKERR(PetscFree(fiat_orientations))
return cell_closure, orientations


@cython.boundscheck(False)
Expand Down
30 changes: 30 additions & 0 deletions firedrake/cython/petschdr.pxi
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,28 @@ cdef extern from "petsc.h":
PETSC_COMPLEX,
PETSC_DATATYPE_UNKNOWN

cdef extern from * nogil:
ctypedef enum PetscDMPolytopeType "DMPolytopeType":
DM_POLYTOPE_POINT
DM_POLYTOPE_SEGMENT
DM_POLYTOPE_POINT_PRISM_TENSOR
DM_POLYTOPE_TRIANGLE
DM_POLYTOPE_QUADRILATERAL
DM_POLYTOPE_SEG_PRISM_TENSOR
DM_POLYTOPE_TETRAHEDRON
DM_POLYTOPE_HEXAHEDRON
DM_POLYTOPE_TRI_PRISM
DM_POLYTOPE_TRI_PRISM_TENSOR
DM_POLYTOPE_QUAD_PRISM_TENSOR
DM_POLYTOPE_PYRAMID
DM_POLYTOPE_FV_GHOST
DM_POLYTOPE_INTERIOR_GHOST
DM_POLYTOPE_UNKNOWN
DM_POLYTOPE_UNKNOWN_CELL
DM_POLYTOPE_UNKNOWN_FACE
DM_NUM_POLYTOPES


cdef extern from "petscsys.h" nogil:
int PetscMalloc1(PetscInt,void*)
int PetscMalloc2(PetscInt,void*,PetscInt,void*)
Expand Down Expand Up @@ -54,6 +76,8 @@ cdef extern from "petscdmplex.h" nogil:
int DMPlexCreatePointNumbering(PETSc.PetscDM,PETSc.PetscIS*)
int DMPlexLabelComplete(PETSc.PetscDM, PETSc.PetscDMLabel)

int DMPlexOrientPoint(PETSc.PetscDM,PetscInt,PetscInt)

cdef extern from "petscdmlabel.h" nogil:
struct _n_DMLabel
ctypedef _n_DMLabel* DMLabel "DMLabel"
Expand All @@ -72,6 +96,12 @@ cdef extern from "petscdm.h" nogil:
int DMGetLabel(PETSc.PetscDM,char[],DMLabel*)
int DMGetPointSF(PETSc.PetscDM,PETSc.PetscSF*)

PetscInt DMPolytopeTypeGetNumVertices(PetscDMPolytopeType)
PetscInt DMPolytopeTypeComposeOrientation(PetscDMPolytopeType,PetscInt,PetscInt)
PetscInt DMPolytopeTypeGetNumArrangments(PetscDMPolytopeType)
PetscInt* DMPolytopeTypeGetArrangment(PetscDMPolytopeType,PetscInt)
PetscInt* DMPolytopeTypeGetVertexArrangment(PetscDMPolytopeType,PetscInt)

cdef extern from "petscdmswarm.h" nogil:
int DMSwarmGetLocalSize(PETSc.PetscDM,PetscInt*)
int DMSwarmGetCellDM(PETSc.PetscDM, PETSc.PetscDM*)
Expand Down
67 changes: 26 additions & 41 deletions firedrake/mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -553,6 +553,7 @@ def callback(self):
self.sfXC = sfXB.compose(self.sfBC) if self.sfBC else self.sfXB
dmcommon.label_facets(self.topology_dm)
dmcommon.complete_facet_labels(self.topology_dm)

# TODO: Allow users to set distribution name if they want to save
# conceptually the same mesh but with different distributions,
# e.g., those generated by different partitioners.
Expand Down Expand Up @@ -589,6 +590,20 @@ def callback(self):
entity_dofs[-2] = 1
facet_numbering = self.create_section(entity_dofs)
self._facet_ordering = dmcommon.get_facet_ordering(self.topology_dm, facet_numbering)

# Try to give the mesh a globally consistent orientation. This will set
# entity_orientations to all zeros, meaning that no transformations are
# required when packing.
# NOTE: I believe that the quadrilateral numbering algorithm has
# effectively been reproduced in DMPlexOrient which claims to spit
# out an error if the mesh is non-orientable. It appears that all of
# our test meshes *are* orientable at the moment so this is hard to
# verify.
# The difference between having an orientable and a non-orientable mesh
# is that an orientable mesh can be coerced into having a closure ordering
# with a (FIAT) orientation of all zeros.
dmcommon.orient_mesh(self.topology_dm, self._vertex_numbering)

self._callback = callback
self.name = name
# Set/Generate names to be used when checkpointing.
Expand Down Expand Up @@ -1089,55 +1104,25 @@ def cell_closure(self):

Each row contains ordered cell entities for a cell, one row per cell.
"""
plex = self.topology_dm
tdim = plex.getDimension()
return self._closure_and_orientations[0]

# Cell numbering and global vertex numbering
cell_numbering = self._cell_numbering
vertex_numbering = self._vertex_numbering.createGlobalSection(plex.getPointSF())
@utils.cached_property
def entity_orientations(self):
return self._closure_and_orientations[1]

@utils.cached_property
def _closure_and_orientations(self):
cell = self.ufl_cell()
assert tdim == cell.topological_dimension()
if cell.is_simplex():
topology = FIAT.ufc_cell(cell).get_topology()
entity_per_cell = np.zeros(len(topology), dtype=IntType)
for d, ents in topology.items():
entity_per_cell[d] = len(ents)

return dmcommon.closure_ordering(plex, vertex_numbering,
cell_numbering, entity_per_cell)

elif cell.cellname() == "quadrilateral":
# can I replace this with calling DMPlexOrient?
if cell.cellname() == "quadrilateral":
from firedrake_citations import Citations
Citations().register("Homolya2016")
Citations().register("McRae2016")
# Quadrilateral mesh
cell_ranks = dmcommon.get_cell_remote_ranks(plex)

facet_orientations = dmcommon.quadrilateral_facet_orientations(
plex, vertex_numbering, cell_ranks)

cell_orientations = dmcommon.orientations_facet2cell(
plex, vertex_numbering, cell_ranks,
facet_orientations, cell_numbering)

dmcommon.exchange_cell_orientations(plex,
cell_numbering,
cell_orientations)

return dmcommon.quadrilateral_closure_ordering(
plex, vertex_numbering, cell_numbering, cell_orientations)
elif cell.cellname() == "hexahedron":
# TODO: Should change and use create_cell_closure() for all cell types.
topology = FIAT.ufc_cell(cell).get_topology()
closureSize = sum([len(ents) for _, ents in topology.items()])
return dmcommon.create_cell_closure(plex, cell_numbering, closureSize)
else:
raise NotImplementedError("Cell type '%s' not supported." % cell)

@utils.cached_property
def entity_orientations(self):
return dmcommon.entity_orientations(self, self.cell_closure)
topology = FIAT.ufc_cell(cell).get_topology()
closureSize = sum([len(ents) for _, ents in topology.items()])
return dmcommon.create_cell_closure(self.topology_dm, self._cell_numbering, closureSize)

@PETSc.Log.EventDecorator()
def _facets(self, kind):
Expand Down