Open jrmaddison opened 1 year ago
A serial example, where evaluation at an out-of-bounds point returns a value.
from firedrake import *
mesh = UnitTriangleMesh()
space = FunctionSpace(mesh, "Lagrange", 1)
X = SpatialCoordinate(mesh)
print(interpolate(X, VectorFunctionSpace(mesh, "Lagrange", 1)).dat.data_ro)
F = Function(space, name="F")
F.interpolate(X[1])
print(F((0.9, 0.9))) # Does not raise an error
I think this originates from #2662. Since libspatialindex indexes items by their bounding box, the update finds elements whose bounding box contains the evaluation point, even if the point lies within the partition for a different process.
This is clearly a bug, and it's quite possibly a result of the FIAT change you point to, however that change is not wrong in and of itself. The change among other things switches to default tolerances for being outside cells which are relative to cell size. This is the right thing to do in order to deal with discretisation error at domain boundaries (and you can override it by user argument if you feel the need). The changes you found are also needed to prevent points getting lost in the boundaries between cells, which has been an ongoing problem.
I think what is going on here is that @ReubenHill has fixed the parallel implications of this via a voting algorithm for VertexOnlyMesh, but this has not been applied to the legacy at
functionality. I think the right answer is probably to modify at
so that it just creates at throw-away VertexOnlyMesh. Can you verify if the same parallel issue occurs using a one-point VertexOnlyMesh?
Can you verify if the same parallel issue occurs using a one-point VertexOnlyMesh?
There is no error when using a VertexOnlyMesh
, the following works.
from firedrake import *
import numpy as np
N_x, N_y, N_z = (2, 2, 2)
mesh = UnitCubeMesh(N_x, N_y, N_z)
vmesh = VertexOnlyMesh(mesh, np.array([[0.2, 0.3, 0.4]]))
X = SpatialCoordinate(mesh)
space = FunctionSpace(mesh, "Lagrange", 1)
vspace = FunctionSpace(vmesh, "Discontinuous Lagrange", 0)
F = Function(space, name="F")
F.interpolate(exp(X[0]) * exp(X[1]))
v = Function(vspace, name="value")
v.interpolate(F)
print(v.dat.data_ro)
This is a real issue - see discussion #2813. If you have a point exactly on the boundary you get point duplication and lose uniqueness. This was introduced in #2662 and requires #2773 (the voting algorithm) to be complete. I will work on this ASAP - in the meantime we can either revert the offending commits or note that this needs fixing.
I don't know how this slipped between my tests but it did! 🤦
@jrmaddison give ReubenHill/voting-algorithm
(PR #2773 ) a go - it should have fixed VertexOnlyMesh
's point duplication behaviour
@ReubenHill Yes, my code works on that branch.
Fixed for VertexOnlyMesh
by #2773
But .at
still demonstrates this behaviour
The following example
leads to a failure when run on 4 processes
There is no error if the expression is changed to
suggesting that one or more processes is extrapolating outside of its partition.