Problem with evaluation at a point in parallel

1

Hello
I have a nonlinear problem and I want to evaluate an unknown at a point. Here is my minimal work:

from dolfin import *

mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), "0file.h5", "r")
hdf.read(mesh, "/mesh", False)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

hdf.read(subdomains, "/subdomains")
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
hdf.read(boundaries, "/boundaries")

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)

# Defining the mixed function space
W_elem = MixedElement([Element1, Element1 ,Element1])

W = FunctionSpace(mesh, W_elem)

z = Function(W)
dz=TrialFunction(W)

alpha, beta,  gamma = split(z)

(v_1, v_2,v_3) = TestFunctions(W)

# Time variables
dt = 0.5
t = 0
T = 2

F = ....

while t <= T:

    bcs = [...]
    J = derivative(F, z, dz)

    problem = NonlinearVariationalProblem(F, z, bcs, J)
    solver = NonlinearVariationalSolver(problem)
    solver.solve()

    (alpha, beta,gamma) = z.split(True)

    print(alpha(0.0245, 0.02))
    print(beta(0.0245, 0.02))

    t += dt

It works well in serial without any error. The problem is when I want to do it in parallel using MPI:
mpirun -n 4 python TEST.py

I get this error:

*** Error:   Unable to evaluate function at point.
*** Reason:  The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.

I found some similar reported issues like here and here but they did not work. It is also suspected as a bug as reported here.

Any help would be appreciated!

2

The mesh is partitioned between processes when running in parallel. You need to ensure that the point at which you wish to evaluate the FE function is contained in the mesh partition on your current process, otherwise extrapolation is required.

However, if you would like every process to see the same point evaluation of a function, you can do something along these lines which I’ve sketched out below

from dolfin import *
import numpy as np

def evaluate_function(u, x):
    comm = u.function_space().mesh().mpi_comm()
    if comm.size == 1:
        return u(*x)

    # Find whether the point lies on the partition of the mesh local
    # to this process, and evaulate u(x)
    cell, distance = mesh.bounding_box_tree().compute_closest_entity(Point(*x))
    u_eval = u(*x) if distance < DOLFIN_EPS else None

    # Gather the results on process 0
    comm = mesh.mpi_comm()
    computed_u = comm.gather(u_eval, root=0)

    # Verify the results on process 0 to ensure we see the same value
    # on a process boundary
    if comm.rank == 0:
        global_u_evals = np.array([y for y in computed_u if y is not None], dtype=np.double)
        assert np.all(np.abs(global_u_evals[0] - global_u_evals) < 1e-9)
    
        computed_u = global_u_evals[0]
    else:
        computed_u = None

    # Broadcast the verified result to all processes
    computed_u = comm.bcast(computed_u, root=0)

    return computed_u


mesh = UnitSquareMesh(32, 32)
f = interpolate(Expression("sin(pi*x[0])*sin(pi*x[1])", degree=4),
                FunctionSpace(mesh, "CG", 1))

x = (0.5, 0.5)
f_x = evaluate_function(f, x)
info("f(%f, %f) = %.3e" % (x[0], x[1], f_x))
3 Likes
3

Thanks for your answer. I tested the example you provided. It works fine when I use 1 process by running: mpirun -n 1 python example.py
But when I want to use multiple processes by running:
mpirun -n 4 python example.py

I get this error:

AttributeError: 'petsc4py.PETSc.Comm' object has no attribute 'gather'

4

You’re getting this error because the version of DOLFIN you are using has not been configured with mpi4py.

If you’re using a FEniCS docker container, this is available in the 2019 stable release.

If you’ve compiled from source, you can simply install mpi4py, e.g: pip3 install mpi4py and recompile the python layer of DOLFIN.

If neither of this appeal to you, you can use the PETSc communicator as you’ve shown in your error message, and change the syntax of my code example accordingly.

5

I resolved this error by changing the comm.size to the number of processes that I want to use. For example if I want to run the code with 4 processes (mpirun -n 4 python example.py) I should do:

def evaluate_function(u, x):
    comm = u.function_space().mesh().mpi_comm()

    #Here is where I should change the number of processes 
    if comm.size == 4:
        return u(*x)

Now MPI seems to work on my machine but the very first issue still exists by this error:

RuntimeError:

*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate function at point.
*** Reason:  The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 3
*** 
*** DOLFIN version: 2017.2.0

And the same error appears on other processes as well. I could not figure it out. Could you please look into it again?
Thanks

6

Hi, consider

from mpi4py import MPI as pyMPI
import numpy as np


def mpi4py_comm(comm):
    '''Get mpi4py communicator'''
    try:
        return comm.tompi4py()
    except AttributeError:
        return comm

    
def peval(f, x):
    '''Parallel synced eval'''
    try:
        yloc = f(x)
    except RuntimeError:
        yloc = np.inf*np.ones(f.value_shape())

    comm = mpi4py_comm(f.function_space().mesh().mpi_comm())
    yglob = np.zeros_like(yloc)
    comm.Allreduce(yloc, yglob, op=pyMPI.MIN)

    return yglob

# -------------------------------------------------------------------

if __name__ == '__main__':
    from dolfin import (UnitSquareMesh, Function, FunctionSpace, interpolate,
                        Expression, VectorFunctionSpace)

    mesh = UnitSquareMesh(32, 32)
    V = VectorFunctionSpace(mesh, 'CG', 1)
    f_ = Expression(('x[0]+x[1]', 'x[0]'), degree=1)
    f = interpolate(f_, V)

    x = np.array([0.5, 0.5])
    assert(np.linalg.norm(f_(x) - peval(f, x)) < 1E-13)
1 Like
7

Thank you both (Mirok and Nate) for your clarifications. I tried what Mirok suggested and it worked.

8

I have a quick question. The method suggested here by @MiroK works perfectly on my machine (Ubuntu and FEniCS 2019.1). However I recently faced with an issue on another machine (with the same environment). When running the code by:

mpirun -np 4 python3 CODE.py

It returns this error:

ValueError: mismatch in send count 8 and receive count 1

It is kind of strange because it works fine on my laptop but returns this error on another computer. Does anybody know where the problem is?
Thanks!

9

I had a similar error when trying @MiroK 's solution:

Traceback (most recent call last):
  File "mwe2.py", line 72, in <module>
    peval(my_model.T.T, x=(0,0))
  File "mwe2.py", line 68, in peval
    comm.Allreduce(yloc, yglob, op=pyMPI.MIN)
  File "mpi4py/MPI/Comm.pyx", line 713, in mpi4py.MPI.Comm.Allreduce
  File "mpi4py/MPI/msgbuffer.pxi", line 699, in mpi4py.MPI._p_msg_cco.for_allreduce
ValueError: mismatch in send and receive MPI datatypes