hi Thanks Matt for the explanation about this.
I have been trying a test which does the following: 1) read in DMPlex from file 2) distribute it, with overlap = 1, using DMPlexDistribute()3) create FVM cell and face geometry vectors using DMPlexComputeGeometryFVM()
4) re-distribute, again with overlap = 1, using DMPlexDistribute()5) distribute the cell and face geometry vectors using DMPlexDistributeField()
Steps 4) and 5) should do essentially nothing, because the mesh has already been distributed (but in my actual non-test code, there is additional stuff between steps 3) and 4) where dual porosity cells are added to the DM).
So I expect the cell and face geometry vectors to be essentially unchanged from the redistribution. And the redistribution SF (from the second distribution) should be just an identity mapping on the cell and face points (except for the overlap ghost points).
This is true for the cells, but not the faces. I've attached the example code and mesh. It is a simple mesh with 10 cells in a horizontal line, each cell 50x50x50 m.
If I run on 2 processes, there are 5 cells (points 0 - 4) on each rank, with centroids at 25, 75, 125, 175 and 225 m on rank 0, and 275, 325, 375, 425 and 475 m on rank 1. The internal faces are the points 36, 42, 47 and 52 on rank 0, and 34, 37, 42, 47 and 52 on rank 1. On rank 0 these should have centroids at 50, 100, 150 and 200 m respectively; on rank 1 they should be at 250, 300, 350 and 400 m. This is true before redistribution.
After redistribution, the cells centroids are still correct, and the face data on rank 1 are OK, but the face data on rank 0 are all wrong.
If you look at the redistribution SF the entries for the rank 0 face data are 36 <- (0,40), 42 <- (0,46), 47 <- (0,51), 52 <- (0,56), instead of the expected 36 <- (0,36), 42 <- (0,42), 47 <- (0,47), 52 <- (0,52). The SF for the rank 1 faces is OK.
If you change the overlap from 1 to 0, it works as expected. So it looks to me like something isn't quite right with the SF for faces when there is overlap. On rank 0 all the entries seem to be shifted up by 4.
I know you originally recommended using overlap = 0 for the initial distribution and only adding overlap for the redistribution. But then Stefano indicated that it should work with overlap now. And it would simplify my code if I could use overlap for the initial distribution (because if dual porosity cells are not being used, then there is no second redistribution).
Is this a bug or is there something I'm doing wrong? - Adrian On 23/06/19 4:39 PM, Matthew Knepley wrote:
On Fri, Jun 21, 2019 at 12:49 AM Adrian Croucher via petsc-users <[email protected] <mailto:[email protected]>> wrote:I have been trying to get this FVM geometry data re-distribution to work using DMPlexDistributeField(). It seems to be working OK for the cell geometry data (cell volumes and centroids). But it is making a mess of the face geometry data (face normals and centroids). Should I even expect DMPlexDistributeField() to work for redistributing a vector of data defined on mesh faces? Or is there some reason I haven't thought of, which means that will never work?Sorry this took a long time. The place I was in France did not have internet.Here is how this stuff works:1) You start with a Section and local vector. The Section describes layout of data in the local vector by mapping mesh points to {# dof, offset}. For a small example, suppose I had two triangles sharing an edge on a sequential mesh for 2 procs.The mesh points would be [0, 1]: Cells [2, 5]: Vertices, where 3,4 are shared [6, 10]: Edges, where 8 is shared A Section for face normals would then look like Process 0 [0, 5]: {0, 0} Meaning no variables lie on cells or vertices 6: {2, 0} One vector per face 7: {2, 2} 8: {2, 4} 9: {2, 6} 10 {2, 8} Process 1 emptyThe vector would just have the face normal values in the canonical order. You can usePetscSectionView() to check that yours looks similar.2) Now we add a PetscSF describing the redistribution. An SF is a map from a given set of integers (leaves) to pairs (int, rank) called roots. Many leaves can point to one root. To begin,we provide an SF mapping mesh points to the new distribution Process 0 0 -> {0, 0} 1 -> {2, 0} 2 -> {3, 0} 3 -> {4, 0} 4 -> {6, 0} 5 -> {7, 0} 6 -> {8, 0} Process 1 0 -> {1, 0} 1 -> {4, 0} 2 -> {3, 0} 3 -> {5, 0} 4 -> {8, 0} 5 -> {9, 0} 6 -> {10, 0}This tells Petsc how to move the points to the parallel distribution. You can use PetscSFView() to check that your point SF is similar. We first use the SF to move the Section data to make a new parallel Section,Process 0 [0, 3]: {0, 0} 4: {2, 0} 5: {2, 2} 6: {2, 4} Process 1 [0, 3]: {0, 0} 4: {2, 0} 5: {2, 2} 6: {2, 4}and then we create a new SF that maps dofs, instead of points, using this new SectionProcess 0: 0 -> {0, 0} First face normal 1 -> {1, 0} 2 -> {2, 0} Second face normal 3 -> {3, 0} 4 -> {4, 0} Third face normal 5 -> {5, 0} Process 1: 0 -> {4, 0} First face normal 1 -> {5, 0} 2 -> {6, 0} Second face normal 3 -> {7, 0} 4 -> {8, 0} Third face normal 5 -> {9, 0}and then this moves the face normal data to the new layout just by broadcasting.You can View the structures you get back to see what Petsc thought it should do.Does this make sense? Thanks, Matt The only example I could find anywhere of DMPlexDistributeField() being used is in DMPlexDistributeCoordinates() so I've been basing what I'm doing on that. (I think I have answered my own questions below by experiment- 1) local vectors should work (they do in DMPlexDistributeCoordinates); 2) probably doesn't matter; 3) yes.) - Adrian On 6/06/19 1:42 PM, Adrian Croucher wrote: > hi > > I have some questions about using the DMPlexDistributeField() > function. I have finite volume mesh geometry data stored in two local > vectors created using DMPlexComputeGeometryFVM(), and I need to > redistribute these after calling DMPlexDistribute() to redistribute my > mesh. (I need the geometry data before redistribution, so I can't just > wait until after redistribution to create them.) > > So I figured DMPlexDistributeField() looks like the thing to use for > that, using the SF that comes out of DMPlexDistribute(). > > 1) Does DMPlexDistributeField() work on local vectors or do they have > to be global ones? > > 2) It takes a 'dm' parameter and the documentation says this is "The > DMPlex object", but is that the original DM (before redistribution) or > the redistributed one, or does it not matter? It looks like it only > uses the DM to get the vector type. > > 3) It looks like you need to manually create the newSection and newVec > output parameters before passing them in to this routine, is that > correct? > > - Adrian >-- Dr Adrian CroucherSenior Research Fellow Department of Engineering Science University of Auckland, New Zealand email: [email protected] <mailto:[email protected]> tel: +64 (0)9 923 4611 --What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.-- Norbert Wienerhttps://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-- Dr Adrian Croucher Senior Research Fellow Department of Engineering Science University of Auckland, New Zealand email: [email protected] tel: +64 (0)9 923 4611
program distgeom
! Redistribution of face geometry- minimal example
#include <petsc/finclude/petsc.h>
use petsc
implicit none
PetscInt, parameter :: overlap = 1
PetscMPIInt :: rank
DM :: dm, dist_dm, redist_dm
PetscErrorCode :: ierr
character(40) :: filename = "gminc_1d.exo"
PetscSF :: dist_sf, redist_sf
Vec :: cell_geom, face_geom
call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
call MPI_comm_rank(PETSC_COMM_WORLD, rank, ierr)
call DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, PETSC_TRUE, dm, ierr)
call DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE, ierr)
call DMPlexDistribute(dm, overlap, dist_sf, dist_dm, ierr)
if (dist_dm .ne. PETSC_NULL_DM) then
call DMDestroy(dm, ierr)
dm = dist_dm
end if
call DMPlexComputeGeometryFVM(dm, cell_geom, face_geom, ierr)
call dm_view_geometry(dm, cell_geom, face_geom)
! redistribute:
call DMPlexDistribute(dm, overlap, redist_sf, redist_dm, ierr)
call redistribute_vec(cell_geom, redist_sf, redist_dm)
call redistribute_vec(face_geom, redist_sf, redist_dm)
call PetscSFView(redist_sf, PETSC_VIEWER_STDOUT_WORLD, ierr)
if (rank == 0) then
write(*,'(/a)') 'redistributed:'
end if
call dm_view_geometry(redist_dm, cell_geom, face_geom)
call VecDestroy(cell_geom, ierr)
call VecDestroy(face_geom, ierr)
call DMDestroy(redist_dm, ierr)
call DMDestroy(dm, ierr)
call PetscSFDestroy(dist_sf, ierr)
call PetscSFDestroy(redist_sf, ierr)
call PetscFinalize(ierr)
contains
subroutine redistribute_vec(v, sf, dist_dm)
Vec, intent(in out) :: v
PetscSF, intent(in) :: sf
DM, intent(in) :: dist_dm
! Locals:
DM :: dm, dist_v_dm
PetscSection :: section, dist_section
Vec :: dist_v
PetscErrorCode :: ierr
call VecGetDM(v, dm, ierr)
call DMGetSection(dm, section, ierr)
call VecCreate(PETSC_COMM_SELF, dist_v, ierr)
call PetscSectionCreate(PETSC_COMM_WORLD, dist_section, ierr)
call DMPlexDistributeField(dm, sf, section, &
v, dist_section, dist_v, ierr)
call DMClone(dist_dm, dist_v_dm, ierr)
call DMSetSection(dist_v_dm, dist_section, ierr)
call VecSetDM(dist_v, dist_v_dm, ierr)
call PetscSectionDestroy(dist_section, ierr)
call VecDestroy(v, ierr)
v = dist_v
end subroutine redistribute_vec
subroutine dm_view_geometry(dm, cell_geom, face_geom)
!! Views FVM geometry on DM
DM, intent(in) :: dm
Vec, intent(in) :: cell_geom, face_geom
! Locals:
DM :: cell_dm, face_dm
PetscInt :: c, start_cell, end_cell
PetscInt :: f, start_face, end_face
PetscInt :: offset
PetscSection :: cell_section, face_section
PetscReal, pointer :: cell_array(:), face_array(:)
PetscReal, pointer :: volume, centroid(:), anormal(:)
PetscInt, pointer :: cells(:)
PetscMPIInt :: rank
PetscErrorCode :: ierr
call mpi_comm_rank(PETSC_COMM_WORLD, rank, ierr)
! cells:
call DMPlexGetHeightStratum(dm, 0, start_cell, end_cell, ierr)
call VecGetDM(cell_geom, cell_dm, ierr)
call DMGetSection(cell_dm, cell_section, ierr)
call VecGetArrayReadF90(cell_geom, cell_array, ierr)
do c = start_cell, end_cell - 1
call PetscSectionGetOffset(cell_section, c, offset, ierr)
offset = offset + 1
centroid => cell_array(offset: offset + 2)
volume => cell_array(offset + 3)
write (*, '(a, i1, a, i2, a, 3(f6.1, 1x), a, e8.3)') &
'rank: ', rank, ' c: ', c, &
' centroid: ', centroid(1), centroid(2), centroid(3), &
' vol: ', volume
end do
call VecRestoreArrayReadF90(cell_geom, cell_array, ierr)
call mpi_barrier(PETSC_COMM_WORLD, ierr)
! faces:
call DMPlexGetHeightStratum(dm, 1, start_face, end_face, ierr)
call VecGetDM(face_geom, face_dm, ierr)
call DMGetSection(face_dm, face_section, ierr)
call VecGetArrayReadF90(face_geom, face_array, ierr)
do f = start_face, end_face - 1
call DMPlexGetSupport(dm, f, cells, ierr)
if (cells(1) /= cells(2)) then
call PetscSectionGetOffset(face_section, f, offset, ierr)
offset = offset + 1
anormal => face_array(offset: offset + 2)
centroid => face_array(offset + 3: offset + 5)
write (*, '(a, i1, a, i2, a, i2, 1x, i2, a, 3(f6.1, 1x))') &
'rank: ', rank, ' f: ', f, &
' cells: ', cells(1), cells(2), &
' centroid: ', centroid(1), centroid(2), centroid(3)
end if
end do
call VecRestoreArrayReadF90(cell_geom, cell_array, ierr)
end subroutine dm_view_geometry
end program distgeom
gminc_1d.exo
Description: Binary data
