Probing a function at a point when running in parallel

Question

I understand that the mesh is split up when running in the code parallel. Is there any way to determine whether a particular point resides in which of the split mesh domain belonging to different process?

For example, I have a function v and if I probe v at point (0,0), v(0,0) works fine if I run on a single processor.

However, if I run the code in parallel, I run into error because the point (0,0) only belongs to a mesh domain associated with a particular process.

Is there an easy way to deal with this?

Answer 1

The trick is to only call evaluate on the process which actually contains the point and then distribute the evaluated value to all other processes (in case that you need the information on all other processes as well). Here's a little example:

from dolfin import *
from numpy import infty

mesh = UnitIntervalMesh(100)
V = FunctionSpace(mesh, 'CG', 1)
f = interpolate(Expression("x[0]"), V)

val = -infty
point = Point(0.2)
if mesh.bounding_box_tree().compute_first_entity_collision(point) < mesh.num_cells():
    val = f(point)

val = MPI.max(mpi_comm_world(), val)

print(val)

Note, that this code implies that if the point is located outside the full mesh this will always just return -infty and not yield an error, so take care.

Comments

Careful also with the MPI.sum. The point could live on more than one CPU...

---

Thx, changed it.

Answer 2

Hi

This could be solved using some try-catch+MPI communication. Otherwise the fenicstools package has a nice solution for this. Examples of both below

from dolfin import *
from numpy import array

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, "CG", 1)
u = interpolate(Expression("x[0]"), V)
x = array([0.25, 0.75])

# Simple example, works only for scalar space
try:
    uu = u(x)
except:
    uu = -1e16
    print "Not found on proc ", MPI.rank(mpi_comm_world()) 

ux = MPI.max(mpi_comm_world(), uu)
print "Proc ", MPI.rank(mpi_comm_world()), " ux = ", ux

from fenicstools import Probes
p = Probes(x, V)  # This works for any mixed space
# Probe once
p(u)
# Probe once again (here same result)
p(u)
# Print results of first probings on rank 0
pa = p.array()
if MPI.rank(mpi_comm_world()) == 0:
    print "Probe = ", pa[0]

Comments

Thanks a lot!

---

For me, using the Probes class in a serial job takes about twice as long as doing a simple loop in python (dof=17k).