On 28.10.2013 16:43, Matthew Knepley wrote::
On Mon, Oct 28, 2013 at 6:17 AM, Bernhard Reinhardt
<[email protected]
<mailto:[email protected]>> wrote:
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
aries.
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?
I have modified your example and it is attached. Here are a few points:
1) You ask for the closure of point 4 in your Vec, and get 1.0. This is
correct.
The closure of 4 is {4, 1, 2}, the cell + its vertices. There are
only values on cells,
so you get the value on 4, which is 1.0.
2) The default for preallocation is FEM-style topology, namely that
unknowns are connected to other
unknowns with overlapping support. The support of your unknowns is
only the cell itself, so there
are no connections by default. I put in a call to
DMPlexSetPreallocationCenterDimension(), which
is the dimension of intermediate points. In FEM these are cells,
but in FVM they are vertices (or
faces). I set your topology to FVM, and you should get the
preallocation you want.
Dear Matt,
thank you for your advice. Now I get the right allocation for my
exemplary problem (although I have to admit I´m not entirely convinced
for the right reasons ;-). I expanded the bar from 3 to 5 volumes, so
that the allocation becomes clearer. See attached file.
However, I still don't get what DMPlexMatSetClosure does. What I would
like to do, is set the operator into the matrix which computes the
change of the value at point p. So let's say in my example the Laplacian
operator (1, -2, 1).
E.g. Delta_p7 = A_21 p6 + A_22 p7 + A_23 p8; with A_21=1, A_22 = -2,
A_23 = 1. However, in my example MatClosure(p7) consits only of A_22.
Probably my sieve-layout is not appropriate. I recognize that the
Laplacian consists of 3 elements, but my volumes have only a cone size
of 2. Nevertheless I hoped to be able to set at least A_21 and A_23.
I tried all combinations of dof=1 or dof=0 at the faces or volumes and
PreallocationCenterDimension 0 or 1 without elucidation. I also tried
different mesh layouts without much succes. That is
- - - - - (without any faces)
or
|-||-||-||-|-| (with double faces)
or
_ _ _ _ _
|-||-||-||-|-| (with double faces and an additional "internal" face)
In the best case I get the foreseen allocation, but am stil am not able
to set the right values in the matrix.
This leads me to the question: Is the sieve-layout appropriate
(|-|-|-|-|-)? If yes, what is the envisaged way to set the operator for
point p in the matrix?
Best 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, array_len
PetscInt :: pStart, pEnd, p, cStart, cEnd, c, vStart, vEnd, v, eStart, eEnd, e
PetscInt :: vecsize
PetscSection :: s
PetscInt, pointer :: pEC(:)
PetscInt, pointer :: pES(:)
PetscScalar, allocatable, target, dimension(:) :: p_array
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)
call DMPlexSetDimension(dm, 1*i1, ierr)
!Consider following 1D mesh with periodic bounds. It consists of 5 faces (0,4), and
!5 "volumes" (5,9)
! |--|--|--|--|--
!
! So initalize 10 sieve points in total
!
call DMPlexSetChart(dm, i0, 10*i1, ierr)
! only the volumes have a cone which consists of the face points
do volume = 5, 9 !
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, 5*i1,[0*i1,1*i1],ierr)
call DMPlexSetCone(dm, 6*i1,[1*i1,2*i1],ierr)
call DMPlexSetCone(dm, 7*i1,[2*i1,3*i1],ierr)
call DMPlexSetCone(dm, 8*i1,[3*i1,4*i1],ierr)
call DMPlexSetCone(dm, 9*i1,[4*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
! DMPLexSetSupportSize and DMPlexSetSupport exist, but instead of
! setting each point individually we just 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 faces. However, they
! should reside in a matrix => 0 dof on the faces
!
!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,1,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 contiguous 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)
do p = 0, pEnd-1
call DMPlexGetSupport(dm,p*i1,pES,ierr)
print *, 'point', p, 'support', pES
call DMPlexRestoreCone(dm,p*i1,pES,ierr)
call DMPlexGetCone(dm,p*i1,pEC,ierr)
print *, 'point', p, 'cone', pEC
call DMPlexRestoreCone(dm,p*i1,pEC,ierr)
call DMPlexVecGetClosure(dm,s,globalVec,p, ptr,ierr)
print *, 'point',p,'->closure', ptr
call DMPlexVecRestoreClosure(dm,s, globalVec, p, ptr, ierr)
enddo
! Choose FVM connection topology
call DMPlexSetPreallocationCenterDimension(dm, 0*i1, ierr);
call DMCreateMatrix(dm, A, ierr)
! create vector with length array_len, set with contiguous numbers
! set pointer to array
array_len = 1
allocate(p_array(array_len))
do i=1, array_len
p_array(i) = i*r1
enddo
ptr => p_array
! Set Closure at point 7
call DMPlexMatSetClosure(dm,s, s, A, 7*i1, ptr, INSERT_VALUES, ierr)
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
