- edited description
Global_index() of vertices wrong if more than one processor and shared_facet ghost cells used
Dear all, I am experiencing bugged behavior if I access the global index of vertices with shared_facet ghost cells on more than one processor.
Please see attached minimal example: I am counting how many facets are connected to each vertex by first iterating over the facets and then incrementing a counter based on the global_index() of each vertex associated with the facet I am currently processing.
I attach two pictures of the results I get: Of course the parallel partitioning is clearly visible, which is as expected. But if you study the results with ghost cells closely, one can spot a vertex that is not attached to any facet! It is located near but not on the partition boundary in the top of the picture.
To me, either the connectivity facets -> vertices is bugged or the global_index() of said vertices. I have tested with facets.normals() and those normals seem fine, so I think the global_index() routine is bugged if ghost layers + parallel happens at the same time.
The results were created with mpiexec -np 2 python IndexBug.py
Comments (37)
-
reporter -
- changed milestone to 1.6
-
assigned issue to
- changed version to dev
- marked as critical
-
The following code seems to be OK for me:
parameters["ghost_mode"] = "shared_facet" mesh = UnitSquareMesh(100,100) mesh.init(0, 1) fmax = 0 fmin = 100 for v in vertices(mesh, "regular"): nfacets = v.num_entities(1) fmax = max(fmax, nfacets) fmin = min(fmin, nfacets) print fmax, fmin
or how about:
parameters["ghost_mode"] = "shared_facet" mesh = UnitSquareMesh(100,100) dim = mesh.topology().dim() Touched = numpy.zeros(mesh.topology().size_global(0)) for MyFace in facets(mesh, "regular"): for Vert in vertices(MyFace): VertIndex = Vert.global_index() Touched[VertIndex] += 1.0 minval = 10 for i in range(len(Touched)): sum_proc = MPI.sum(mesh.mpi_comm(), Touched[i]) minval = min(minval, sum_proc) print minval
I agree that your code produces a strange zero, but I can't quite understand what you are doing - I'm also not sure that dof_to_vertex etc. works well in parallel...
-
reporter Sorry for posting this twice, my internet connection is fairly bad right now.
My code is a very much stripped down minimal example of what I actually need: A FEniCS function that contains the vertex normal on boundary vertices and a vector zero inside the domain, in parallel with ghost layers. Hence, I need to count boundary vertices for averaging the facet normals to boundary vertices.
I can confirm that the FEniCS mesh.move(V) functionality also does not work for vertices for which my minimal example above gives zero answer when used with my vertex normal CG 1 VectorFunction. It appears that there are vertices, which are not reachable when iterating over facets...
A vertex normal (rather than a facet normal) would be a great feature for doing shape optimization and deforming geometries with FEniCS: http://www.mathematik.uni-wuerzburg.de/~schmidt/femorph
-
reporter The dof_to_vertex() has never failed me in parallel so far, with or without ghost layers enabled, it always maps from local to local indices.
-
dof_to_vertex_map
is well-defined in parallel although I have no idea whether and how it interacts withghost_mode
. Only pitfall is thatdof_to_vertex_map
did not map unowned DOFs at the time the function has been written (because of global DOF numbering at that time) and I'm not sure whether the function has been rewritten to contain also unowned portion of DOFs with the switch to local DOF numbering by Garth. -
I'd suggest that @Epoxid tries to prove that global indices are incorrect by a more obvious way. By definition: global index is the same on all processes for each entity, rank-wise unique and < global number of entities. Is some of those requirements violated and easily observable?
-
reporter Hmm... I will try to help where I can. But can we also include the dof_to_vertex map in this?
-
Well, it'd would be good if you were able to prove the claim without touching any dofmap at all. But if there's no other way...
-
reporter Oh, I am not strictly claiming that it's the global index, I am just saying, that under the circumstances described above, the chain facet -> vertex -> dof -> function seems to fail occasionally...
-
reporter Returning to the two alternatives Chris Richardson has kindly provided, without checking them in all detail, they don't produce the same answer for me when more than two processors are used. The first example at some point starts to report a minimum of fmin = 3, but the second version always reports 2 for minval, independent of the number of processors...
-
What about
for MyFace in facets(mesh, "regular"): if True: #if MyFace.exterior() == True: for Vert in vertices(MyFace): VertIndex = Vert.global_index() Touched[VertIndex] += 1.0 for MyFace in facets(mesh, "ghost"): if True: #if MyFace.exterior() == True: for Vert in vertices(MyFace): VertIndex = Vert.global_index() Touched[VertIndex] += 1.0
This solves the problem. Although I'm not sure that I understand how a ghosted facet is defined. Also notice vertex iterator in inner loop. This is probably also exclusively non-ghosted.
-
@Epoxid - the results seem correct to me. In the first case, as you increase the process count, some processes will no longer have the corners (with two facets per vertex). In the second case, all vertices are being summed by MPI, so it is always the same.
@blechta - using "all" should be the same as "ghost" + "regular". A ghost facet is one which only belongs to a ghost cell (and is not shared with a regular cell).
-
Ok, let's make it more simple and avoid dofmaps at all:
from dolfin import * parameters["ghost_mode"] = "shared_facet" mesh = UnitSquareMesh(100, 100) # This loops all the vertices local_indices = set() global_indices = set() for f in facets(mesh, "all"): for v in vertices(f): local_indices.add(v.index()) global_indices.add(v.global_index()) assert len(local_indices) == len(global_indices) assert len(local_indices) == mesh.num_vertices() # This loops a subset local_indices = set() global_indices = set() for f in facets(mesh, "regular"): for v in vertices(f): local_indices.add(v.index()) global_indices.add(v.global_index()) assert len(local_indices) == len(global_indices) assert len(local_indices) < mesh.num_vertices() or MPI.size(mesh.mpi_comm()) == 1 assert len(local_indices) == mesh.topology().ghost_offset(0)
@chris_richardson, can you confirm the second loop (vertices of regular facets) does what intended? Should it really keep some vertices unvisited?
BTW, there should be a check for
opt=="all"
inMeshEntityIterator(const Mesh& mesh, std::size_t dim, std::string opt)
otherwise it keeps space for uncatched typos, likefor c in cells(mesh, 'regullar'):
. -
Lets' make
std::string opt
anenum
. -
Then I suggest using new-style enums, see example. For instance
enum class EntityType : int32_t { all = 0, regular = 1, ghost = 2 };
which is wrapped by SWIG automatically as
EntityType_all
,EntityType_regular
,EntityType_ghost
. (I'm not sure whether this is what IPython users are dying for. Can somebody test thatvoid list_timings(TimingClear clear, std::set<TimingType> type);
works with all the IPython magic?) -
Yes, I guess it can be an enum.
The following code illustrates that all vertices are visited, even when using "regular" facets - at least, it works for me. The vertices of "regular" facets provide full coverage of all vertices of the global mesh. Maybe there is some confusion about
num_vertices()
which I think also includes the ghosts.global_indices = set() for f in facets(mesh, "regular"): for v in vertices(f): global_indices.add(v.global_index()) rmin = 100000 for i in range(mesh.size_global(0)): r = MPI.sum(mesh.mpi_comm(), int(i in global_indices)) rmin = min(r, rmin) print rmin
-
You're right. Adding line
assert len(local_indices) == mesh.topology().ghost_offset(0)
makes it clear.The question is whether the vertex coverage by (regular facets vertices) is the same as ownership of respective vertex nodes in dofmap. This would explain Stephan's problem. (I'm currently to busy to investigate this.)
-
reporter Dear all, thank you very much for the support! On a side node from a user: Might I suggest to make this difference between mesh.num_vertices() and mesh.size_global(0) more clear in the documentation?
http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/cpp/mesh/Mesh.html#dolfin.cpp.mesh.Mesh.num_vertices and http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/cpp/mesh/Mesh.html#dolfin.cpp.mesh.Mesh.size_global
To me as a user, reading both documentations, I would expect those to do the same, which is returning what "mesh.topology().ghost_offset(0)" actually returns...
-
Yes, there is always room to improve documentation.
num_vertices()
is the same issize(0)
which is the number of local verticessize_global(0)
is the number of global vertices across all processes.topology().ghost_offset(0)
is the number of local, non-ghosted vertices. -
reporter Dear all, thank you so much for looking into this. In the meantime, I studied the behavior of the dof_to_vertex_map with ghost facets in parallel. I think the dof_to_vertex_map returns ghost vertices instead of regular vertices, which caused the original problem? Anyway... please have a look at this minimal example:
from dolfin import * parameters["ghost_mode"] = "shared_facet" mesh = UnitSquareMesh(100, 100) mesh.init() NSpaceS = FunctionSpace(mesh, "CG", 1) d2v = dof_to_vertex_map(NSpaceS) N = Function(NSpaceS) dmin, dmax = N.vector().local_range() for dof in range(dmax-dmin): vert = Vertex(mesh, d2v[dof]) index = vert.index() if index > mesh.topology().ghost_offset(0): print "I touched a ghost vertex"
As a user, I would expect the dof_to_vertex_map to always give me regular or physical vertices and not ghost ones...
-
The ghost vertices will have dofs associated with them, just as shared vertices do in a non-ghosted mesh. If you don't need them, you can just ignore them, I guess.
-
reporter Ah! Great! Thanks for the answer! I can very much understand having these dofs in the DG context. But in the CG context, somehow all of these different local dofs need to be synchronised into one value.
I suppose this is done in the N.vector().apply("") step. Maybe the problem is rooted there? unfortunately, I could not find the apply function in the online documentation...
-
reporter Thank you again for the great discussion. Please have a look at this little program. To me, this demonstrates that the local vertex associated with the local dof number zero is never touched when looping through the local mesh...
from dolfin import * parameters["ghost_mode"] = "shared_facet" mesh = UnitSquareMesh(100, 100) mesh.init() NSpaceS = FunctionSpace(mesh, "CG", 1) d2v = dof_to_vertex_map(NSpaceS) VertexOfDofZero = d2v[0] print "Vertex of dof zero", VertexOfDofZero Touched = False for MyFace in facets(mesh, "regular"): if True: #if MyFace.exterior() == True: for Vert in vertices(MyFace): lIndex = Vert.index() if lIndex == VertexOfDofZero: Touched = True if Touched == False: print "I have not touched the vertex for dof 0"
I have executed the program a couple of times with different processors and it always reports that the vertex for dof 0 is not being processed when looping though the mesh. It is processed as soon as I turn off ghost.
-
What you say makes sense. If you loop over "regular" vertices, there will be some ghost dofs which you do not touch. Change "regular" to "all" and you will reach all dofs.
-
reporter Thank you for the great support and feedback! I took heart of your suggestions and now I loop over "all" facets instead of the "regular" facets. But I am still seeing paradox behavior. Please have a look at the program below:
1) I am now looping over "all" facets and mark their vertices if the facet is exterior by the vertex global index.
2) A CG 1 function is created and the result is used as the vector() of the function
3) I loop over the function dofs and check if all marked vertices are really on the exterior and all unmarked vertices are really in the interior. This is done independently by checking their coordinates.
The remarkable thing: The check in 3) reports all is okay: All interior vertices are marked, all exterior vertices are unmarked and vice versa.
But if I open the output file in Paraview, atleast on my system, the file has exactly one pair of wrongly marked vertices: One interior is 1, one exterior is 0. It is as if these two are mixed-up...
from dolfin import * import numpy parameters["ghost_mode"] = "shared_facet" def VertOnSurface(vert): if abs(vert.x(0)) < DOLFIN_EPS: return True if abs(vert.x(0)-1.0) < DOLFIN_EPS: return True if abs(vert.x(1)) < DOLFIN_EPS: return True if abs(vert.x(1)-1.0) < DOLFIN_EPS: return True return False mesh = UnitSquareMesh(100, 100) mesh.init() dim = mesh.topology().dim() Touched = numpy.zeros(mesh.topology().size_global(0)) for MyFace in facets(mesh, "all"): if MyFace.exterior() == True: for Vert in vertices(MyFace): VertIndex = Vert.global_index() Touched[VertIndex] = 1.0 NSpaceS = FunctionSpace(mesh, "CG", 1) d2v = dof_to_vertex_map(NSpaceS) N = Function(NSpaceS) dmin, dmax = N.vector().local_range() LocValues = numpy.zeros(dmax-dmin) for dof in range(dmax-dmin): Vert = Vertex(mesh, d2v[dof]) gIndex = Vert.global_index() value = Touched[gIndex] LocValues[dof] = value N.vector().set_local(LocValues) N.vector().apply("") #test if all marked vertices are really on the surface: for dof in range(dmax-dmin): Vert = Vertex(mesh, d2v[dof]) if abs(N.vector()[dof] - 1.0) < DOLFIN_EPS and VertOnSurface(Vert) == False: print "An interior vertex was marked!" if abs(N.vector()[dof]) < DOLFIN_EPS and VertOnSurface(Vert) == True: print "An exterior vertex was not marked!" File("./output/N.pvd", "compressed") << N
-
reporter Here is the result of the above program when I run it on two processors in parallel
-
It could be the vtk output which is wrong. Have you tried xdmf.
-
reporter I am very sure that this will also happen with xdmf.
Unfortunately, I cannot test this right now, because the precompiled FEniCS binaries available here (http://fenicsproject.org/download/osx_details.html) still crash for me when writing xdmf in parallel.
So, I kindly have to ask you to change one line to verify this.
-
reporter I now had the opportunity to test this further with the following minimal example that works entirely without global indices or connectivity:
from dolfin import * import numpy parameters["ghost_mode"] = "shared_facet" mesh = UnitSquareMesh(100,100) NSpaceS = FunctionSpace(mesh, "CG", 1) d2v = dof_to_vertex_map(NSpaceS) N = Function(NSpaceS) dmin, dmax = N.vector().local_range() LocValues = numpy.zeros(dmax-dmin) for dof in range(dmax-dmin): Vert = Vertex(mesh, d2v[dof]) LocValues[dof] = Vert.x(0) N.vector().set_local(LocValues) N.vector().apply("") File("./output/X1.pvd", "compressed") << N XDMFFile(mpi_comm_world(), "./output/X1.xdmf") << N
This little demo is just supposed to turn the x-coordinate of the (x,y) points into a CG 1 function. And it always fails to me if more than one processor is used, no matter if I use xdmf or pvd, OS X or debian. If you look closely, you will notice there is always one node that is not correctly colored. On my OS X machine, the odd vertex is on the lower edge, on my debian machine, it is exactly the lower right corner...
-
Looks like
dof_to_vertex_map
is not right with ghosts, then.mesh = UnitSquareMesh(3, 3) Q = FunctionSpace(mesh, "CG", 1) d2vi = [] for dof in dof_to_vertex_map(Q): d2vi.append(Vertex(mesh, dof).global_index()) vi = [] for v in vertices(mesh, "all"): vi.append(v.global_index()) print sorted(d2vi), sorted(vi)
One dof is repeated, and one missed out in the
dof_to_vertex_map
. -
I reminf,
dof_to_vertex_map
does not map unowned because dofs where originally globally indexed. With "ghosts" it should be similar story. Let me look into this later this week. -
Filed as
#501. -
reporter I am glad we narrowed this down! Do you have any idea of a work-around? This bug is pretty much a show stopper for my projects with FEniCS...
-
- changed status to invalid
Looks like a problem with dof_to_vertex map, see Issue
#501 -
Now fixed in master.
-
- removed milestone
Removing milestone: 1.6 (automated comment)
- Log in to comment