<div dir="ltr">Dear all, <div><br></div><div>I haven't been on this mailing list since december, so first of all, happy new year to everyone !</div><div><br></div><div>I am facing another challenge with the DMPlex structure : although I am doing cell-centered finite volume (using PetscFV), I also need to do some manipulations and store information that are vertex-based. </div><div>The problem I have is when I run the code in parallel : when I visualize the nodal data using a vtk viewer, it shows NaN at all the vertices in close proximity to the MPI partition line.</div><div>Here is how I proceed to work on the vertex-based data :</div><div><br></div><div>First comes the creation of the DM that I'll be using for the cell-centered FVM</div><div><br></div><div><font face="tahoma, sans-serif">call DMPlexCreateFromFile(comm, name, PETSC_TRUE, dm, ierr)</font></div><div><font face="tahoma, sans-serif">call DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_TRUE, ierr)</font></div><div><font face="tahoma, sans-serif">call DMPlexDistribute(dm, 1, PETSC_NULL_SF, dmDist, ierr)</font></div><div><font face="tahoma, sans-serif">...</font></div><div><font face="tahoma, sans-serif">dm = dmDist</font></div><div><font face="tahoma, sans-serif">...</font></div><div><font face="tahoma, sans-serif">call DMConstructGhostCells(dm, ..., ..., dmGhost, ierr)</font></div><div><font face="tahoma, sans-serif">...</font></div><div><font face="tahoma, sans-serif">dm = dmGhost</font></div><div><font face="tahoma, sans-serif">...</font></div><div><br></div><div>I set both boolean to PETSC_TRUE in the SetBasicAdjacency so even in MPI the neighbouring relationships are conserved for all strata. Is that true ? Like, does doing this actually enable to get all the transitive closure of a given vertex, even if part of it is in another MPI block ? </div><div><br></div><div>Then I need a vertex-based DM, so I figure it's pretty much what's done in CreateMassMatrix in TS ex11.c.</div><div><br></div><div><font face="tahoma, sans-serif">call DMClone(dm, dmNodal, ierr)</font></div><div><font face="tahoma, sans-serif">call DMGetCoordinateSection(dm, coordSection, ierr); CHKERRA(ierr)<br>call DMGetCoordinatesLocal(dm, coordinates, ierr); CHKERRA(ierr)<br>call DMSetCoordinateSection(dmNodal, PETSC_DETERMINE, coordSection, ierr);<br>call PetscSectionCreate(PETSC_COMM_WORLD, sectionNodal, ierr); CHKERRA(ierr)<br>call DMPlexGetDepthStratum(dm, 0, vstart, vend, ierr); CHKERRA(ierr)<br>call PetscSectionSetChart(sectionNodal, vstart, vend, ierr); CHKERRA(ierr)<br>do inode = vstart, vend-1<br> ! 3 dof per node, ndim = 3<br>call PetscSectionSetDof(sectionNodal, inode, ndim, ierr); CHKERRA(ierr)<br>end do<br>call PetscSectionSetUp(sectionNodal, ierr); CHKERRA(ierr)<br>call DMSetLocalSection(dmNodal, sectionNodal, ierr); CHKERRA(ierr)<br>call PetscSectionDestroy(sectionNodal, ierr); CHKERRA(ierr)</font><br></div><div><br></div><div>From what I understood, this new dmNodal basically has all the characteristics of the ori dm, except it's primary layer is vertex-based with 3 DOF per vertex.</div><div><br></div><div>Then I get a local vector and work on it :</div><div><br></div><div><font face="tahoma, sans-serif">call DMGetLocalVector(dmNodal, nodalvec, ierr)</font></div><div><font face="tahoma, sans-serif">call VecGetArrayF90(nodalvec, nodalvecarray, ierr)</font></div><div><font face="tahoma, sans-serif">call DMPlexGetDepthStratum(dm, 0, vstart, vend, ierr)</font></div><div><font face="tahoma, sans-serif">call DMPlexGetDepthStratum(dm, 2, fstart, fend, ierr)<br></font></div><div><font face="tahoma, sans-serif">! For each node ...</font></div><div><font face="tahoma, sans-serif">do inode = vstart, vend-1</font></div><div><font face="tahoma, sans-serif">A = 0 ! matrix, ndim x ndim</font></div><div><font face="tahoma, sans-serif">b = 0 ! vector, ndim</font></div><div><font face="tahoma, sans-serif">call DMPlexGetTransitiveClosure(dmNodal, inode, PETSC_FALSE, vertexclosure, ierr);<br></font></div><div><font face="tahoma, sans-serif">! For each face attached to that node</font></div><div><font face="tahoma, sans-serif">do iface = 1, size(vertexclosure)</font></div><div><font face="tahoma, sans-serif">if ((vertexclosure(iface) < fstart) .or. (vertexclosure(iface) >= fend)) cycle<br></font></div><div><font face="tahoma, sans-serif">A = A + ...</font></div><div><font face="tahoma, sans-serif">b = b + ...</font></div><div><font face="tahoma, sans-serif">end do</font></div><div><font face="tahoma, sans-serif">nodalvecarray(1+(inode-vstart)*ndim:3+(inode-vstart)*ndim) = matmul(A, b)</font></div><div><font face="tahoma, sans-serif">...</font></div><div><font face="tahoma, sans-serif">end do</font></div><div><br></div><div>Then I want to visualize it :</div><div><br></div><div><font face="tahoma, sans-serif">call VecRestoreArrayF90(nodalvec, nodalvecarray, ierr)</font></div><div><font face="tahoma, sans-serif">call DMCreateGlobalVector(dmNodal, globalnodalvec, ierr)</font></div><div><font face="tahoma, sans-serif">call DMLocalToGlobal(dmNodal, nodalvec, INSERT_VALUES, globalnodalvec, ierr)</font></div><div><font face="tahoma, sans-serif">call PetscViewerCreate(PETSC_COMM_WORLD, vtkViewer, ierr); CHKERRA(ierr)<br>call PetscViewerSetType(vtkViewer, PETSCVIEWERVTK, ierr); CHKERRA(ierr)<br>call PetscViewerFileSetName(vtkViewer, "toto.vtk", ierr); CHKERRA(ierr)<br>call VecView(globalnodalvec, vtkViewer, ierr); CHKERRA(ierr)<br>call PetscViewerDestroy(vtkViewer, ierr); CHKERRA(ierr)</font><br></div><div><br></div><div>In the toto.vtk file, I get vertex-based values as expected, and the whole shebang works in sequential. Now when I run it in parallel, I get vertex-based values inside the MPI blocks, but in the neighbourhood of the join between the MPI blocks, all the vertex-based values are NaN (see attachment).</div><div><br></div><div>All in all, I think I do not understand fully the consequences of the DMPlexDistribute with useCone=true & useClosure=true ... but I cannot figure out where the NaN are coming from : is it because the vertices do not actually know all their closure through the MPI join, or is it some Global / Local error I made ?</div><div><br></div><div>Thank you very much in advance for your support !!</div><div><br clear="all"><div><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><div><div><div><div><div>Thibault</div></div></div></div></div></div></div></div></div></div></div></div></div>