<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sun, Apr 6, 2014 at 9:43 PM, Adrian Croucher <span dir="ltr"><<a href="mailto:a.croucher@auckland.ac.nz" target="_blank">a.croucher@auckland.ac.nz</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF">
<br>
<div>On 05/04/14 02:31, Matthew Knepley
wrote:<br>
</div>
<blockquote type="cite">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote"><br>
<div>Instead of this, use</div>
<div><br>
</div>
<div> call DMPlexGetDefaultSection(cell_dm, section, ierr);
(Move outside of loop)</div>
<div> call PetscSectionGetOffset(section, c, off, ierr);</div>
<div> part(off) =rank</div>
</div>
</div>
</div>
</blockquote>
<br>
OK thanks. The code compiles now.<br>
<br>
However, at runtime I'm getting another error, when it calls
DMSetPointSF():<br>
<br>
[0]PETSC ERROR: #1 DMSetPointSF() line 3204 in
/home/acro018/software/PETSc/code/src/dm/interface/dm.c<br>
[1]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------<br>
[1]PETSC ERROR: Invalid argument<br>
[1]PETSC ERROR: Wrong type of object: Parameter # 1<br>
[1]PETSC ERROR: See
<a href="http://http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for
trouble shooting.<br>
[1]PETSC ERROR: Petsc Development GIT revision:
v3.4.4-5306-g9849faf GIT Date: 2014-04-06 15:35:41 -0500<br>
<br>
This seems to work fine in C, so I'm wondering if I need to do
something different in fortran, or if the interface isn't quite
working right. My test code is below. The offending call to
DMSetPointSF() is in the routine create_partition_vector(). Any tips
much appreciated!<br></div></blockquote><div><br></div><div>1) You do not need those calls since DMClone() now does this automatically.</div><div><br></div><div>2) However, it looks like DMClone() is not functioning correctly for you since cell_dm should be alright.</div>
<div><br></div><div>I will try and run your code today.</div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div text="#000000" bgcolor="#FFFFFF">
It also occurred to me there might just be something wrong with my
mesh, so I also tested it using a mesh created with
DMPLexCreateBoxMesh()- but this gave the same error (and also
indicated that the fortran interface to DMPlexSetRefinementLimit()
is probably also missing in action!).<br></div></blockquote><div><br></div><div>Fixed. Will push it today.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF">
- Adrian<br>
<br>
--<br>
<br>
program testplex<br>
<br>
! Testing PETSc dmplex from fortran, using a line of blocks in
3-D.<br>
! Create the DMPlex by passing in the cell definitions to
DMPlexCreateFromDAG(),<br>
! interpolate the result to form the faces and edges, distribute
over the<br>
! processors, and visualize the partitions in VTK.<br>
<br>
implicit none<br>
<br>
#include <finclude/petsc.h90><br>
<br>
PetscInt, parameter :: num_cells = 20<br>
PetscReal, parameter :: cell_size = 10.<br>
PetscErrorCode :: ierr<br>
DM :: dm, dist_dm, cell_dm<br>
character(5) :: partitioner = "chaco"<br>
PetscInt :: overlap = 0<br>
Vec :: partition<br>
<br>
call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)<br>
<br>
call create_1d_grid(num_cells, cell_size, dm, ierr); CHKERRQ(ierr)<br>
<br>
call DMPlexDistribute(dm, partitioner, overlap, PETSC_NULL_OBJECT,
dist_dm, ierr)<br>
CHKERRQ(ierr)<br>
call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
dm = dist_dm<br>
<br>
call DMView(dm, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr)<br>
<br>
call create_partition_vector(dm, cell_dm, partition, ierr);
CHKERRQ(ierr)<br>
<br>
! Visualize partition vector here...<br>
<br>
call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
call PetscFinalize(ierr); CHKERRQ(ierr)<br>
<br>
end program testplex<br>
<br>
!------------------------------------------------------------------------<br>
<br>
subroutine create_1d_grid(num_cells, cell_size, dm, ierr)<br>
<br>
! Creates grid consisting of a 1D line of cube cells in 3D.<br>
<br>
implicit none<br>
<br>
#include <finclude/petsc.h90><br>
<br>
PetscInt, intent(in) :: num_cells<br>
PetscReal, intent(in):: cell_size<br>
DM, intent(out) :: dm<br>
PetscErrorCode, intent(out) :: ierr<br>
! Locals:<br>
PetscInt, parameter :: dim = 3<br>
PetscMPIInt :: rank<br>
PetscInt :: num_vertices_per_face, num_vertices_per_cell<br>
PetscInt :: num_vertices, total_num_points<br>
PetscInt :: i, i0, i1, l<br>
PetscInt :: depth<br>
real(8) :: x0, dx<br>
PetscInt, allocatable :: num_points(:)<br>
PetscInt, allocatable :: cone_size(:)<br>
PetscInt, allocatable :: cones(:)<br>
PetscInt, allocatable :: cone_orientations(:)<br>
PetscScalar, allocatable :: vertex_coords(:)<br>
DM :: interp_dm<br>
<br>
ierr = 0<br>
call MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr); CHKERRQ(ierr)<br>
<br>
call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr); CHKERRQ(ierr)<br>
call PetscObjectSetName(dm, "testplex", ierr); CHKERRQ(ierr)<br>
call DMPlexSetDimension(dm, dim, ierr); CHKERRQ(ierr)<br>
<br>
if (rank == 0) then<br>
<br>
num_vertices_per_face = 2 ** (dim-1)<br>
num_vertices_per_cell = 2 * num_vertices_per_face<br>
num_vertices = num_vertices_per_face * (num_cells + 1)<br>
l = num_vertices_per_face * dim<br>
total_num_points = num_vertices + num_cells ! # points in DAG<br>
<br>
allocate(vertex_coords(0:num_vertices*dim-1),
cones(0:num_cells*num_vertices_per_cell-1))<br>
<br>
! Set up vertex coordinates<br>
dx = real(cell_size)<br>
<br>
do i = 0, num_cells<br>
i0 = i * l<br>
i1 = (i+1)*l<br>
x0 = real(i * cell_size)<br>
vertex_coords(i0:i1) = \<br>
[x0, 0.d0, 0.d0,\<br>
x0, dx, 0.d0,\<br>
x0, 0.d0, dx,\<br>
x0, dx, dx]<br>
end do<br>
<br>
! Set up cells<br>
do i = 0, num_cells-1<br>
i0 = i * num_vertices_per_cell<br>
i1 = (i+1) * num_vertices_per_cell<br>
cones(i0:i1) = [0, 1, 5, 4, 2, 6, 7, 3] + i *
num_vertices_per_face <br>
end do<br>
cones = cones + num_cells<br>
<br>
num_points = [num_vertices, num_cells]<br>
allocate(cone_size(0: total_num_points-1))<br>
allocate(cone_orientations(0:num_cells*num_vertices_per_cell-1))<br>
cone_size = 0<br>
cone_size(0: num_cells-1) = num_vertices_per_cell<br>
cone_orientations = 0<br>
depth = 1<br>
<br>
call DMPlexCreateFromDAG(dm, depth, num_points, cone_size,
cones, &<br>
cone_orientations, vertex_coords, ierr); CHKERRQ(ierr)<br>
<br>
else ! rank > 0: create empty DMPlex<br>
<br>
num_points = [0, 0, 0]; cone_size = [0]<br>
cones = [0]; cone_orientations = [0]<br>
vertex_coords = [0]<br>
call DMPlexCreateFromDAG(dm, 1, num_points, cone_size, &<br>
cones, cone_orientations, vertex_coords, ierr);
CHKERRQ(ierr)<br>
<br>
end if<br>
<br>
interp_dm = PETSC_NULL_OBJECT<br>
call DMPlexInterpolate(dm, interp_dm, ierr); CHKERRQ(ierr)<br>
call DMPlexCopyCoordinates(dm, interp_dm, ierr); CHKERRQ(ierr)<br>
call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
dm = interp_dm<br>
<br>
deallocate(num_points, cone_size, cones, cone_orientations,
vertex_coords)<br>
<br>
end subroutine create_1d_grid<br>
<br>
!------------------------------------------------------------------------<br>
<br>
subroutine create_partition_vector(dm, cell_dm, partition, ierr)<br>
<br>
! Returns a DM and vector containing partitioning data<br>
<br>
implicit none<br>
<br>
#include <finclude/petsc.h90><br>
<br>
DM, intent(in) :: dm<br>
DM, intent(out) :: cell_dm<br>
Vec, intent(out) :: partition<br>
PetscErrorCode, intent(out) :: ierr<br>
! Locals:<br>
MPI_Comm :: comm<br>
PetscMPIInt :: rank<br>
PetscSF :: sfPoint<br>
PetscSection :: coord_section<br>
Vec :: coordinates<br>
PetscSection :: cell_section<br>
PetscScalar, pointer :: part(:)<br>
PetscInt :: cStart, cEnd, c, off<br>
<br>
call DMGetCoordinateSection(dm, coord_section, ierr);
CHKERRQ(ierr)<br>
call DMGetCoordinatesLocal(dm, coordinates, ierr); CHKERRQ(ierr)<br>
call DMClone(dm, cell_dm, ierr); CHKERRQ(ierr)<br>
call DMGetPointSF(dm, sfPoint, ierr); CHKERRQ(ierr)<br>
call DMSetPointSF(cell_dm, sfPoint, ierr); CHKERRQ(ierr)<br>
call DMSetCoordinateSection(cell_dm, coord_section, ierr);
CHKERRQ(ierr)<br>
call DMSetCoordinatesLocal(cell_dm, coordinates, ierr);
CHKERRQ(ierr)<br>
call PetscObjectGetComm(dm, comm, ierr); CHKERRQ(ierr)<br>
call PetscSectionCreate(comm, cell_section, ierr); CHKERRQ(ierr)<br>
call DMPlexGetHeightStratum(cell_dm, 0, cStart, cEnd, ierr);
CHKERRQ(ierr)<br>
call PetscSectionSetChart(cell_section, cStart, cEnd, ierr);
CHKERRQ(ierr)<br>
do c = cStart, cEnd-1<br>
call PetscSectionSetDof(cell_section, c, 1, ierr);
CHKERRQ(ierr)<br>
end do<br>
call PetscSectionSetUp(cell_section, ierr); CHKERRQ(ierr)<br>
call DMSetDefaultSection(cell_dm, cell_section, ierr);
CHKERRQ(ierr)<br>
call PetscSectionDestroy(cell_section, ierr); CHKERRQ(ierr)<br>
call DMCreateLocalVector(cell_dm, partition, ierr); CHKERRQ(ierr)<br>
call PetscObjectSetName(partition, "partition", ierr);
CHKERRQ(ierr)<br>
call VecGetArrayF90(partition, part, ierr); CHKERRQ(ierr)<br>
call MPI_Comm_rank(comm, rank); CHKERRQ(ierr)<br>
do c = cStart, cEnd-1<br>
call PetscSectionGetOffset(cell_section, c, off, ierr);
CHKERRQ(ierr)<br>
part(off) = rank<br>
end do<br>
call VecRestoreArrayF90(partition, part, ierr); CHKERRQ(ierr)<br>
<br>
end subroutine create_partition_vector<span class="HOEnZb"><font color="#888888"><br>
<br>
<br>
<br>
<pre cols="72">--
Dr Adrian Croucher
Department of Engineering Science
University of Auckland
New Zealand
tel 64-9-373-7599 ext 84611
</pre>
</font></span></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>