mikaem / fenicstools

55 stars 36 forks source link

Simple probe fails #19

Closed jywuce closed 7 years ago

jywuce commented 7 years ago

I have the following simple code test.py for probe:

from dolfin import from numpy import array from fenicstools import

mesh = UnitSquareMesh(10, 10) V = FunctionSpace(mesh, 'CG', 1) u = interpolate(Expression("sin(x[0]pi)cos(x[1]*pi)"), V)

x = array([0.6, 0.5]) probe = Probe(x, V) probe(u)

print "result = ", probe.get_probe_sub(0)

It gives the correct result "0" (almost) if I run it in series: python test.py However, when I run it in parallel: mpirun -np 2 python test.py, it seems that the probe fails.

------------------------------------------------------------------------- DOLFIN encountered an error. If you are not able to resolve this issue *** using the information listed below, you can ask for help at


*** fenics-support@googlegroups.com


Remember to include the error message listed below and, if possible, include a minimal running example to reproduce the error.


------------------------------------------------------------------------- Error: Unable to set probe. Reason: Probe is not found on processor. Where: This error was encountered inside Probe.cpp. *** Process: 0


DOLFIN version: 2016.1.0 Git changeset:
*** -------------------------------------------------------------------------

Very interestingly, if I probe the array [0.5, 0.5], it does give the correct result. Furthermore, I tried to use compute_first_entity_collision() directly in fenics, it also fails. I think it might be some bugs in fenics. Could you please confirm it or fix the bug? Thanks!

MiroK commented 7 years ago

Yes, I can reproduce this. The case [0.5, 0.5] most likely passes because the point lies on the boundary between CPUs. Anyways, I have no issues with entity_collision, in particular the following works as expected (i.e. only one CPU finds the point)

from dolfin import *
from numpy import array
from fenicstools import *

mesh = UnitSquareMesh(10, 10)

x = array([0.6, 0.5])
p = Point(x)

tree = mesh.bounding_box_tree()
foo = tree.compute_first_collision(p) 
print foo

I would like the probe initialization to fail only if the probe is not found on any of the processes, how about it @mikaem ?

jywuce commented 7 years ago

Yes, MiroK (the same on FEniCS Q&A?). Your codes works as expected. I also noticed that in your cpp code, you use tree.compute_first_collision() to probe the point, but it fails. So, where is the bug, in FEnics or FEnicsTools?

mikaem commented 7 years ago

Hi The Probe is not supposed to work in parallel. You need to use Probes. Probes tries to evaluate a Probe on each process, and only when successful is a probe created locally. The error message you see is caught by Probes.Probes()

jywuce commented 7 years ago

Thanks, Mihaem. With the Probes(), now the code works.

from dolfin import from numpy import array from fenicstools import

mesh = UnitSquareMesh(10, 10) V = FunctionSpace(mesh, 'CG', 1) u = interpolate(Expression("sin(x[0]pi)cos(x[1]*pi)"), V)

x = array([[0.25, 0.50], [0.62, 0.34], [0.28, 0.25]]) probes = Probes(x.flatten(), V) probes(u)

print probes.array()

I tried probes.array() as said in wiki, but it also gives "none" for those unprobed process (I used mpirun -np 2).

Process 0: Computed global bounding box tree with 3 boxes. Process 1: Computed global bounding box tree with 3 boxes. None [ 4.27647350e-17 4.44188722e-01 5.29222330e-01]

How can I just pick up the last line

[ 4.27647350e-17 4.44188722e-01 5.29222330e-01]

which are meaningful for me?

jywuce commented 7 years ago

Any one can help?

mikaem commented 7 years ago

If you only want to print on root, then you can do something like

a = probes.array()
if comm.Get_rank() == 0:
    print(a)

probes.array is a function called by all ranks, and all but root gets None as return type.