VertexOnlyMesh Improvements Plan

[ReubenHill]

The work plan for improvements to VertexOnlyMesh was outlined here #2223.

However, the recent merge of #2437 has made using setting redundant == PETSC_TRUE in DMSwarmSetPointCoordinates no longer viable since #2437 switches to using mesh.locate_cell for point location.

First item on old plan

The issue with using mesh.locate_cell was outlined in #2224: namely that you get point duplication for points on cell boundaries. The plan had been to update DMSwarmSetPointCoordinates to use redundant == PETSC_TRUE by default (meaning that rank 0 would broadcast its point list to all other ranks) such that this would not happen. From #2224:

Notably it seems that DMSwarmSetPointCoordinates doesn't suffer from this problem so it would be worth checking how the implementation there works.

It seems this isn't correct. Consider test_parallel_vom_create.py:

from firedrake import *
m = UnitSquareMesh(1,1)
coords = [[0.5, 0.5], [0.9, 0.9]]
vm = VertexOnlyMesh(m, coords)
print(f'Rank {COMM_WORLD.rank}, coordinates {vm.coordinates.dat.data_ro}')

Before and after the merging of #2437 running this would give you point duplication:

~ mpiexec -n 3 python test_parallel_vom_create.py
Rank 0, coordinates [[0.5 0.5]
 [0.9 0.9]]
Rank 1, coordinates []
Rank 2, coordinates [[0.5 0.5]]

Prior to merging #2437 this diff

diff --git a/firedrake/mesh.py b/firedrake/mesh.py
index ab8563da..dba3b23d 100644
--- a/firedrake/mesh.py
+++ b/firedrake/mesh.py
@@ -2283,7 +2283,7 @@ def _pic_swarm_in_plex(plex, coords, fields=[]):
     # coordinates associated with them. The DMPlex cell id associated
     # with each PIC in the DMSwarm is accessed with the `DMSwarm_cellid`
     # field.
-    swarm.setPointCoordinates(coords, redundant=False, mode=PETSc.InsertMode.INSERT_VALUES)
+    swarm.setPointCoordinates(coords, redundant=True, mode=PETSc.InsertMode.INSERT_VALUES)
 
     # Remove PICs which have been placed into ghost cells of a distributed DMPlex
     dmcommon.remove_ghosts_pic(swarm, plex)

ought to have avoided point duplication (if my comment on #2224 is correct) but, confusingly, it still has the issue

~ mpiexec -n 3 python test_parallel_vom_create.py
Rank 0, coordinates [[0.5 0.5]
 [0.9 0.9]]
Rank 1, coordinates []
Rank 2, coordinates [[0.5 0.5]]

So what's going on here?
Do I need to do something I've forgotten about to make my diff have an effect?
Was I just wrong when I said it avoided the issue?

What to do next instead

It's been a while since I wrote the issues and my work plan due to my extended leave of absence.

Any suggestions on how to proceed?

I will tidy up the various issues above following discussion here.

[ksagiyam]

As far as I understand, DMSwarmSetPointCoordinates (followed by DMLocatePoints) does not have a mechanism to remove such point duplication.

Maybe we have enough PETSc functions exposed to implement the voting algorithm in Firedrake.

[ReubenHill]

Looks like my comment

Notably it seems that DMSwarmSetPointCoordinates doesn't suffer from this problem so it would be worth checking how the implementation there works.

was wrong then. Thanks Koki

[ReubenHill]

This is in the process of changing - see https://gitlab.com/petsc/petsc/-/merge_requests/5354 and https://gitlab.com/petsc/petsc/-/merge_requests/5266

[jwallwork23]

Also https://gitlab.com/petsc/petsc/-/merge_requests/5345

[dham]

I think this is all basically fine and nothing we've done really gets in the way.

1. We can move to requiring broadcast from rank 0. That doesn't really change anything else from where we are now. The issues with duplicating or potentially losing points on process ownership boundaries are there now and wouldn't be fixed by this.
2. Given that we're assuming that broadcasting the entire set of points is an OK thing to do, we can fix the process boundary implementation by a simple vote. What we do is change a bit the procedure for dropping vertices we don't own. To do this we would:
   a. modify the current code so that we attempt to find cells in the halo as well as the core (easy).
   b. delay actually dropping verticess until after we vote (also easy). We write into number of vertices long integer array the process number we think owns the process. We write -1 if this rank doesn't know. We can then do a collective max over this array to get unambiguous ownership.
   c. now drop the vertices we don't own.

This is all orthogonal to the tolerance issue, which can happen on cell boundaries interior to a process. We should fix that too.

I think we should also change the default behaviour for point loss to raising an error. We shouldn't silently drop vertices unless the user has explicitly indicated that this is OK by them.

[ReubenHill]

This is great.

API Thoughts

For anyone using VertexOnlyMesh at the moment we are proposing default API changes.

There are, I think, 3 variations

1. Broadcasting from rank 0:

vm = VertexOnlyMesh(parent_mesh, coords, redundant=True,  point_migration=None)

2. Each rank may have different points and points on process boundaries can get duplicated or lost (current behaviour)

vm = VertexOnlyMesh(parent_mesh, coords, redundant=False, point_migration=False)

3. Migrating points to correct rank where each rank has different points (i.e. coming up with rank candidates as outlined in [meta] VertexOnlyMesh parallel creation improvements #2223)

vm = VertexOnlyMesh(parent_mesh, coords, redundant=False, point_migration=True)

I think in most cases people will be using redundant lists of points so (1.) should be the default.
If redundant=False then we should default to case (3.) (point_migration=True) but give NotImplementedError('Point migration is not yet supported, set point_migration to False') for now forcing it to be explicitly set by the user.

Error Behaviour

We already have 3 error modes to consider which are already implemented missing_points_behaviour=None, missing_points_behaviour='warn' and missing_points_behaviour='error' which were added in PR #2187.

The existing implementation of this error requires all ranks to advertise the same list of points but the check scales with number of vertices and ranks (see

firedrake/firedrake/mesh.py

Lines 2194 to 2199 in e6e778a

	# check all ranks have the same vertexcoords so that check is valid
	# NOTE this operation scales with number of vertices and ranks
	_, allequal = mesh.comm.allreduce((vertexcoords, True), op=op)
	op.Free()
	if not allequal:
	raise ValueError("Cannot check for missing points if different vertices on each MPI rank!")

)
With the new proposed implementation I think there must be a way of avoiding this scaling tracking if no process is found for any point but I can't quite work out what it would be. This should also allow us to support case (3.) above.

Thoughts?

[wence-]

To determine if any points were dropped, you need to intersect over the processes.

For the case where everyone starts with the same list of vertices, this is trivial given the voting strategy above. Any point whose allocated rank is -1 after voting was not found anywhere.

If you start with distributed vertices (each rank has its own list of points), then you can compute the intersection as follows:

Determine the set of points I didn't find locally as $M_\text{j} := {p_i}$, we have one such set $M_j$ for each rank. Now just compute the intersection over all the ranks of the $\bigcap_{j} M_j$. If this set is not empty, then any points in it were found on no process.

This same strategy also works if you implement point migration for distributed vertices. You just insert a round where you distribute the $M_j$ to those ranks that might have the points, let them search and then you get some new $M_j^\prime$ which you intersect.

[ReubenHill]

Re draft API, we discussed this in the firedrake meeting and came to the conclusions

• DH: Remove point_migration keyword argument as this is premature optimisation.
• DH: redundant=True should also be strongly encouraged since, if you have so many points as to need it, the locate-cell functionality will kill performance.

[ReubenHill]

@dham you mentioned in #2452 that switching to a voting algorithm would allow for points in mesh halos but still ensure uniqueness. How so?

[ReubenHill]

You don't need to. A vertexonlymesh just has some points that are split between processes arbitrarily but in a way which guarantees that we can evaluate at those points using data from the parent mesh that is available locally. So when you interpolate into a space defined on the vertexonlymesh you have all the data you need locally.

That was a helpful explanation, thank you.

Assuming that no one can think of any cases where having points in halos is a problem (maybe particle methods?) then I'm satisfied that we can have points in the halos as long as those points are unique (as decided by something like a voting algorithm). We can add point migration later. The PETSc SF for the VertexOnlyMesh will need to be updated - something along the lines of making sure the leaves are correct I think.

[ReubenHill]

When we do get to doing point migration, can anyone think why DMSwarmMigrate might not work out of the box if we relabel the DMSWARM with process ownership as necessary?

[ksagiyam]

Not sure about DMSwarmMigrate().

After knowing how many particles exist in each cell, I think we can recast the DMSwarm as PetscSection (with appropriate permutation set as in

firedrake/firedrake/cython/dmcommon.pyx

Line 948 in 93e1c64

CHKERR(PetscSectionSetPermutation(section.sec, renumbering.iset))

) representing the non-uniform DoF layout over the cell stratum. Once this is done, I think we can use PetscSFCreateSectionSF() with the parent plex's pointSF to get swarm's pointSF.

Then I can directly use this PetscSection to checkpoint a VertexOnlyMesh and the fields defined on it.

[ReubenHill]

recast the DMSwarm as PetscSection

What do you mean by this? All I think I know about PetscSection is that we use them to hold certain bits of function space information. A DMSwarm, here taking a similar role to a DMPlex in other meshes, seems to me to be a very different data type

[dham]

@ksagiyam wants to attach the swarm to the plex of the parent mesh in order to checkpoint it. He observes that you can think of a swarm as a function space over the parent mesh in which there are varying numbers of degrees of freedom in each cell. This effectively makes a way to store the swarm in the plex. It's probably not what you want in memory when you're computing (because a swarm is the right abstraction) but it's a feasible way of checkpointing and loading a swarm and the data on it.

[colinjcotter]

It is for computing values for points that depend on values on points that belong to other processors. It seems like you don't need to have points in a parent mesh halo, because another processor will compute them.

…

________________________________ From: Reuben W. Nixon-Hill ***@***.***> Sent: 21 June 2022 14:45 To: firedrakeproject/firedrake ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [firedrakeproject/firedrake] VertexOnlyMesh Improvements Plan (Discussion #2450) Can you remind me what the point of the halo is on a normal mesh? Voting algorithm aside, why should I need to have points in a parent mesh halo? — Reply to this email directly, view it on GitHub<#2450 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4RS243YT3V2B5MQCC3VQHBPZANCNFSM5WZEJG4A>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[danshapero]

Are there possible improvements to single processor performance from using a different spatial data structure for point location? I'm not familiar with how a lot of the parallelism works under the hood, so take this with a grain of salt, but there could be a lot of value to be gained for less effort from putting bounding boxes for each triangle into an R-tree or quadtree.

[wence-]

we do that already