Hi there,
I am trying to familiarize myself with unstructured grids. For this I
started with the example given in the petsc-manual (mesh of 2 triangles)
and tried to modify it according to my needs.
However, the example seems to be incomplete. At least I was not able to
create a matrix without applying a DMPlexSetDimension.
In the end I decided to set up an an example as simple as possible:
1D transport of a scalar in a bar with periodic boundaries.
Lets say I have a bar consisting of 3 volumes and 3 faces.
You could write this down as a vector of length 3 containing the scalar
in the volumes and a 3x3 matrix containing the transport coefficients on
the faces.
I have already tried many different ways of setting this up with plex
but always failed. Attached you find the version that makes most sense
to me. It gives me vectors of length 3 and a 3x3 matrix. However, the
matrix contains no writable entries and I think the returned closures
are not correct.
Is there a simple example with good documentation of the individual
steps like in the manual which is complete and working?
Kind regards
Bernhard
program main
implicit none
#include "finclude/petsc.h90"
PetscInt :: N
PetscErrorCode :: ierr
PetscInt, parameter :: i1 = 1, i0 = 0
PetscScalar,parameter :: r1 = 1.0, r0 = 0.0
PetscBool :: flg
Vec :: globalVec
Mat :: A
DM :: dm
integer :: i, t, myrank, comm_size, loc_vec_size, volume
PetscInt :: pStart, pEnd, p, cStart, cEnd, c, vStart, vEnd, v, eStart, eEnd, e
PetscInt :: vecsize
PetscSection :: s
real(8), dimension(:), pointer :: ptr
! Setup PETSC
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
call MPI_Comm_rank(PETSC_COMM_WORLD,myrank,ierr)
call MPI_Comm_size(PETSC_COMM_WORLD,comm_size,ierr)
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)
call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr)
!Consider following 1D mesh with periodic bounds. It consits of 3 boundaries (0,2), and
!3 "volumes" (3,5)
! |--|--|--
!
! So initalize 6 sieve points in total
!
call DMPlexSetChart(dm, i0, 6*i1, ierr)
! only the volumes have a cone which consits of the boundary points
do volume = 3, 5 !
call DMPlexSetConeSize(dm, volume*i1, 2*i1, ierr)
enddo
call DMSetUp(dm, ierr)
! tell plex the cone relations. periodic conditions show in last line
call DMPlexSetCone(dm, 3*i1,[0*i1,1*i1],ierr)
call DMPlexSetCone(dm, 4*i1,[1*i1,2*i1],ierr)
call DMPlexSetCone(dm, 5*i1,[2*i1,0*i1],ierr)
! Now setup the supporting relationships for each sieve point
! That is reverse the arrows in the Hasse diagram and from which points
! p can now be reached (hope this description is correct)
! DMPLexSetSupportSize and DMPlexSetSupport exist, but instead of
! setting each point individually we jus call DMPlexSymmetrize
call DMPlexSymmetrize(dm,ierr)
!In order to support efficient queries, we also want to construct fast
!search structures, indices, for the different
!types of points, which is done using
call DMPlexStratify(dm,ierr)
! This will give us the number of points in each level of the Hasse diag
call DMPlexGetChart(dm, pStart, pEnd, ierr)
call DMPlexGetHeightStratum(dm, 0, cStart, cEnd, ierr)
call DMPlexGetHeightStratum(dm, 1, eStart, eEnd, ierr)
call DMPlexGetDepthStratum(dm, 0, vStart, vEnd, ierr)
call PetscSectionCreate(PETSC_COMM_WORLD, s, ierr)
call PetscSectionSetChart(s, pStart, pEnd, ierr)
print *, 'pstart', pstart, 'pEnd', pEnd
print *, 'cStart', cStart, 'cEnd', cEnd
print *, 'eStart', eStart, 'eEnd', eEnd
print *, 'vStart', vStart, 'vEnd', vEnd
! Set the degrees of freedom
! Assume we want to calc. the transport of a scalar through the 1D mesh.
! So we have 1 dof on the "volumes"
! We also need transport coefficients on the boundaries. However, they
! should reside in a matrix => 0 dof on the boundaries
!
!do p = pStart, pEnd-1
!call PetscSectionSetDof(s,p,1,ierr)
!enddo
do c = cStart, cEnd-1
call PetscSectionSetDof(s,c,1,ierr)
enddo
!do v = vStart, vEnd-1
!call PetscSectionSetDof(s,v,1,ierr)
!enddo
!do e = eStart, eEnd-1
!call PetscSectionSetDof(s,e,0,ierr)
!enddo
call PetscSectionSetUp(s, ierr)
! Now lets get vectors!
call DMSetDefaultSection(dm, s, ierr)
call DMGetGlobalVector(dm, globalVec,ierr)
call VecGetSize(globalVec,vecsize, ierr)
print *, 'size of globalVec ', vecsize
!Set the vector with coninuous numbers to see if we can retrieve
!the right values using the plex routines
do c = 0, vecsize-1
call VecSetValue(globalVec, c, c*r1, INSERT_VALUES, ierr)
enddo
call VecAssemblyBegin(globalVec,ierr)
call VecAssemblyEnd(globalVec,ierr)
call VecView(globalVec, PETSC_VIEWER_STDOUT_SELF, ierr)
! Why does is the closure of sieve point 4 not (1 , 2)?
call DMPlexVecGetClosure(dm,s,globalVec,4*i1, ptr,ierr)
print *, 'closure', ptr
! Why do I need to set a dimension (else create Matrix will fail)?
! What are the consequences? It seems the number of Dimension does not
! have any effect. What would be the right number?
!
call DMPlexSetDimension(dm,1*i1, ierr)
call DMCreateMatrix(dm, MATAIJ, A, ierr)
! I guess know we would use DMPlexMatSetClosure to set the transport
! coefficients. However in the current setup the matrix contains 0
! writable elements :(
call MatAssemblyBegin(A,ierr)
call MatAssemblyEnd(A,ierr)
call MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr)
call MatView(A, PETSC_VIEWER_DRAW_WORLD, ierr)
read *,ierr
call DMDestroy(dm,ierr)
call PETSCFINALIZE(ierr)
end program main
