VertexOnlyMesh in parallel

[jrmaddison]

Is that a way to map between the coordinates used to instantiate a VertexOnlyMesh, and the elements of a Dat defined for function on the mesh?

For example with the construction

# 'u' a function defined on function space 'space'

vmesh = VertexOnlyMesh(space.mesh(), X_coords)
vspace = FunctionSpace(vmesh, "Discontinuous Lagrange", 0)

interp = Interpolator(TestFunction(space), vspace)

p = Function(vspace, name="p")
interp.interpolate(u, output=p)

the length of p.dat.data_ro can vary between different processes, and the global sum of the lengths may be longer than the number of points in X_coords.

Is there a way to identify which for which i and j p.dat.data_ro[i] corresponds to X_coords[j, :]?

[ReubenHill]

Not at the moment, which I'm aware is an issue particularly now that broadcasting the coordinates from rank 0 is the default behaviour. For now, assuming none of your coordinates are the same, you can manually look at the coordinates using vmesh.coordinates.data.data_ro_with_halos on each rank and match them to the data you have. At present I empty the halos of points as a way of ensuring uniqueness, which I know isn't ideal, but you ought to find that the total length of vmesh.coordinates.data.data_ro_with_halos does correspond to the length of X_coords - if not, you must be losing points. You can check this by doing

vmesh = VertexOnlyMesh(space.mesh(), X_coords, missing_points_behaviour='error')

to throw an error if points are lost. This will soon become the default behaviour.

I'm currently in the process of allowing points in the halos and adding point ids which you can query to deal with the case where coordinates may match.

Of course the ideal scenario would be to supply a global list of point values which matches the global list of point coordinates and have them automatically assigned. This could use the PETSc SF which I need to set up for having points in the halos but I'm not sure what the API should be. Perhaps we could use Function.assign? Suggestions are welcome! @connorjward you're good at this sort of thing.

[connorjward]

Of course the ideal scenario would be to supply a global list of point values which matches the global list of point coordinates and have them automatically assigned. This could use the PETSc SF which I need to set up for having points in the halos but I'm not sure what the API should be. Perhaps we could use Function.assign? Suggestions are welcome! @connorjward you're good at this sort of thing.

I suppose that assign makes sense here. It does seem strange to supply a global list of values rather than local but it's consistent with how VertexOnlyMesh is set up. I'm happy to discuss this with you more at some point.

This does seem like a bit of a tangent though. It seems like the answer to the question asked here is to inspect vmesh.coordinates.data.data_ro_with_halos.

[jrmaddison]

Thanks, I'll work from vmesh.coordinates.data.data_ro_with_halos.

you ought to find that the total length of vmesh.coordinates.data.data_ro_with_halos does correspond to the length of X_coords

It isn't shorter, but can be longer, e.g. as in

from firedrake import *

import numpy as np

X_coords = np.array([[0.4, 0.2, 0.3]], dtype=np.float64)

mesh = UnitCubeMesh(10, 10, 10)
vmesh = VertexOnlyMesh(mesh, X_coords)
vspace = FunctionSpace(vmesh, "Discontinuous Lagrange", 0)

F = Function(vspace, name="F")

comm = vmesh.comm
for r in range(comm.size):
    comm.barrier()
    if r == comm.rank:
        print(r, F.dat.data_ro.shape, flush=True)
    comm.barrier()

which displays

0 (0,)
1 (1,)
2 (0,)
3 (1,)

[ReubenHill]

Of course the ideal scenario would be to supply a global list of point values which matches the global list of point coordinates and have them automatically assigned. This could use the PETSc SF which I need to set up for having points in the halos but I'm not sure what the API should be. Perhaps we could use Function.assign? Suggestions are welcome! @connorjward you're good at this sort of thing.

I suppose that assign makes sense here. It does seem strange to supply a global list of values rather than local but it's consistent with how VertexOnlyMesh is set up. I'm happy to discuss this with you more at some point.

This does seem like a bit of a tangent though. It seems like the answer to the question asked here is to inspect vmesh.coordinates.data.data_ro_with_halos.

Well I had originally thought that everyone would supply a local list of coordinates which you'd then have the local list of values to match - you can still do this using vmesh = VertexOnlyMesh(space.mesh(), X_coords, redundant = False) such that you don't broadcast the coordinates from rank 0. I changed the default to be redundant = True a few months ago.

[ReubenHill]

Thanks, I'll work from vmesh.coordinates.data.data_ro_with_halos.

you ought to find that the total length of vmesh.coordinates.data.data_ro_with_halos does correspond to the length of X_coords

It isn't shorter, but can be longer, e.g. as in

from firedrake import *

import numpy as np

X_coords = np.array([[0.4, 0.2, 0.3]], dtype=np.float64)

mesh = UnitCubeMesh(10, 10, 10)
vmesh = VertexOnlyMesh(mesh, X_coords)
vspace = FunctionSpace(vmesh, "Discontinuous Lagrange", 0)

F = Function(vspace, name="F")

comm = vmesh.comm
for r in range(comm.size):
    comm.barrier()
    if r == comm.rank:
        print(r, F.dat.data_ro.shape, flush=True)
    comm.barrier()

which displays

0 (0,)
1 (1,)
2 (0,)
3 (1,)

Argh - so what's happening here is we have a point that's exactly on the cell boundary and both ranks think they own it. I thought I had tests to guard against this but clearly not. This is a bug - we're getting point coordinate duplication - which will have been introduced by this PR and which will be fixed by my new voting algorithm (which I haven't finished yet). Really sorry! I'll note this on your issue and get it fixed ASAP.