Actual source code: ex12f.F

  1: !
  2: !    "$Id: ex12.F,v 1.60 2001/08/07 03:04:13 balay Exp $";
  3: !
  4: !  This example demonstrates basic use of the SNES Fortran interface.
  5: !
  6: !  Note:  The program ex10.f is the same as this example, except that it
  7: !         uses the Fortran .f suffix rather than the .F suffix.
  8: !
  9: !  In this example the application context is a Fortran integer array:
 10: !      ctx(1) = da    - distributed array
 11: !          2  = F     - global vector where the function is stored
 12: !          3  = xl    - local work vector
 13: !          4  = rank  - processor rank
 14: !          5  = size  - number of processors
 15: !          6  = N     - system size
 16: !
 17: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 18: !  MUST be declared as external.
 19: !
 20: !
 21: ! Macros to make setting/getting  values into vector clearer.
 22: ! The element xx(ib) is the ibth element in the vector indicated by ctx(3)
 23: #define xx(ib)  vxx(ixx + (ib))
 24: #define ff(ib)  vff(iff + (ib))
 25: #define F2(ib)  vF2(iF2 + (ib))
 26:       program main
 27:       implicit none

 29:  #include include/finclude/petsc.h
 30:  #include include/finclude/petscvec.h
 31:  #include include/finclude/petscda.h
 32:  #include include/finclude/petscmat.h
 33:  #include include/finclude/petscsnes.h

 35:       PetscFortranAddr ctx(6)
 36:       integer          rank,size,ierr,N,start,end,nn,i,ii,its,flg
 37:       SNES             snes
 38:       Mat              J
 39:       Vec              x,r,u
 40:       PetscScalar      xp,FF,UU,h
 41:       character*(10)   matrixname
 42:       external         FormJacobian,FormFunction

 44:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 45:       N = 10
 46:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',N,flg,ierr)
 47:       h = 1.d0/(N-1.d0)
 48:       ctx(6) = N

 50:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 51:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 52:       ctx(4) = rank
 53:       ctx(5) = size

 55: ! Set up data structures
 56:       call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,N,1,1,            &
 57:      &     PETSC_NULL_INTEGER,ctx(1),ierr)

 59:       call DACreateGlobalVector(ctx(1),x,ierr)
 60:       call DACreateLocalVector(ctx(1),ctx(3),ierr)

 62:       call PetscObjectSetName(x,'Approximate Solution',ierr)
 63:       call VecDuplicate(x,r,ierr)
 64:       call VecDuplicate(x,ctx(2),ierr)
 65:       call VecDuplicate(x,U,ierr)
 66:       call PetscObjectSetName(U,'Exact Solution',ierr)

 68:       call MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N, &
 69:      &     N,3,PETSC_NULL_INTEGER,0,PETSC_NULL_INTEGER,J,ierr)

 71:       call MatGetType(J,matrixname,ierr)

 73: ! Store right-hand-side of PDE and exact solution
 74:       call VecGetOwnershipRange(x,start,end,ierr)
 75:       xp = h*start
 76:       nn = end - start
 77:       ii = start
 78:       do 10, i=0,nn-1
 79:         FF = 6.0*xp + (xp+1.e-12)**6.e0
 80:         UU = xp*xp*xp
 81:         call VecSetValues(ctx(2),1,ii,FF,INSERT_VALUES,ierr)
 82:         call VecSetValues(U,1,ii,UU,INSERT_VALUES,ierr)
 83:         xp = xp + h
 84:         ii = ii + 1
 85:  10   continue
 86:       call VecAssemblyBegin(ctx(2),ierr)
 87:       call VecAssemblyEnd(ctx(2),ierr)
 88:       call VecAssemblyBegin(U,ierr)
 89:       call VecAssemblyEnd(U,ierr)

 91: ! Create nonlinear solver
 92:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

 94: ! Set various routines and options
 95:       call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
 96:       call SNESSetJacobian(snes,J,J,FormJacobian,ctx,ierr)
 97:       call SNESSetFromOptions(snes,ierr)

 99: ! Solve nonlinear system
100:       call FormInitialGuess(snes,x,ierr)
101:       call SNESSolve(snes,x,ierr)
102:       call SNESGetIterationNumber(snes,its,ierr);

104: ! Write results if first processor
105:       if (ctx(4) .eq. 0) then
106:         write(6,100) its
107:       endif
108:   100 format('Number of Newton iterations = ',i5)

110: !  Free work space.  All PETSc objects should be destroyed when they
111: !  are no longer needed.
112:       call VecDestroy(x,ierr)
113:       call VecDestroy(ctx(3),ierr)
114:       call VecDestroy(r,ierr)
115:       call VecDestroy(U,ierr)
116:       call VecDestroy(ctx(2),ierr)
117:       call MatDestroy(J,ierr)
118:       call SNESDestroy(snes,ierr)
119:       call DADestroy(ctx(1),ierr)
120:       call PetscFinalize(ierr)
121:       end


124: ! --------------------  Evaluate Function F(x) ---------------------

126:       subroutine FormFunction(snes,x,f,ctx,ierr)
127:       implicit none
128:       SNES             snes
129:       Vec              x,f
130:       PetscFortranAddr ctx(*)
131:       integer          rank,size,i,s,n,ierr
132:       PetscOffset      ixx,iff,iF2
133:       PetscScalar      h,d,vf2(1),vxx(1),vff(1)
134:  #include include/finclude/petsc.h
135:  #include include/finclude/petscvec.h
136:  #include include/finclude/petscda.h
137:  #include include/finclude/petscmat.h
138:  #include include/finclude/petscsnes.h


141:       rank  = ctx(4)
142:       size  = ctx(5)
143:       h     = 1.d0/(ctx(6) - 1.d0)
144:       call DAGlobalToLocalBegin(ctx(1),x,INSERT_VALUES,ctx(3),ierr)
145:       call DAGlobalToLocalEnd(ctx(1),x,INSERT_VALUES,ctx(3),ierr)

147:       call VecGetLocalSize(ctx(3),n,ierr)
148:       if (n .gt. 1000) then
149:         print*, 'Local work array not big enough'
150:         call MPI_Abort(PETSC_COMM_WORLD,0,ierr)
151:       endif

153: !
154: ! This sets the index ixx so that vxx(ixx+1) is the first local
155: ! element in the vector indicated by ctx(3).
156: !
157:       call VecGetArray(ctx(3),vxx,ixx,ierr)
158:       call VecGetArray(f,vff,iff,ierr)
159:       call VecGetArray(ctx(2),vF2,iF2,ierr)

161:       d = h*h

163: !
164: !  Note that the array vxx() was obtained from a ghosted local vector
165: !  ctx(3) while the array vff() was obtained from the non-ghosted parallel
166: !  vector F. This is why there is a need for shift variable s. Since vff()
167: !  does not have locations for the ghost variables we need to index in it
168: !  slightly different then indexing into vxx(). For example on processor
169: !  1 (the second processor)
170: !
171: !        xx(1)        xx(2)             xx(3)             .....
172: !      ^^^^^^^        ^^^^^             ^^^^^
173: !      ghost value   1st local value   2nd local value
174: !
175: !                      ff(1)             ff(2)
176: !                     ^^^^^^^           ^^^^^^^
177: !                    1st local value   2nd local value
178: !
179:        if (rank .eq. 0) then
180:         s = 0
181:         ff(1) = xx(1)
182:       else
183:         s = 1
184:       endif

186:       do 10 i=1,n-2
187:        ff(i-s+1) = d*(xx(i) - 2.d0*xx(i+1)                              &
188:      &      + xx(i+2)) + xx(i+1)*xx(i+1)                                &
189:      &      - F2(i-s+1)
190:  10   continue

192:       if (rank .eq. size-1) then
193:         ff(n-s) = xx(n) - 1.d0
194:       endif

196:       call VecRestoreArray(f,vff,iff,ierr)
197:       call VecRestoreArray(ctx(3),vxx,ixx,ierr)
198:       call VecRestoreArray(ctx(2),vF2,iF2,ierr)
199:       return
200:       end

202: ! --------------------  Form initial approximation -----------------

204:       subroutine FormInitialGuess(snes,x,ierr)
205:       implicit none
206:  #include include/finclude/petsc.h
207:  #include include/finclude/petscvec.h
208:  #include include/finclude/petscsnes.h
209:       integer          ierr
210:       Vec              x
211:       SNES             snes
212:       PetscScalar      five

214:       five = 5.d-1
215:       call VecSet(five,x,ierr)
216:       return
217:       end

219: ! --------------------  Evaluate Jacobian --------------------

221:       subroutine FormJacobian(snes,x,jac,B,flag,ctx,ierr)
222:       implicit none
223:  #include include/finclude/petsc.h
224:  #include include/finclude/petscvec.h
225:  #include include/finclude/petscda.h
226:  #include include/finclude/petscmat.h
227:  #include include/finclude/petscsnes.h
228:       SNES             snes
229:       Vec              x
230:       Mat              jac,B
231:       PetscFortranAddr ctx(*)
232:       integer          flag,ii,istart
233:       PetscOffset      ixx
234:       integer          iend,i,j,n,rank,size,end,start,ierr
235:       PetscScalar      d,A,h,vxx(1)


238:       h = 1.d0/(ctx(6) - 1.d0)
239:       d = h*h
240:       rank = ctx(4)
241:       size = ctx(5)

243:       call VecGetArray(x,vxx,ixx,ierr)
244:       call VecGetOwnershipRange(x,start,end,ierr)
245:       n = end - start

247:       if (rank .eq. 0) then
248:         A = 1.0
249:         call MatSetValues(jac,1,start,1,start,A,INSERT_VALUES,ierr)
250:         istart = 1
251:       else
252:         istart = 0
253:       endif
254:       if (rank .eq. size-1) then
255:         i = ctx(6)-1
256:         A = 1.0
257:         call MatSetValues(jac,1,i,1,i,A,INSERT_VALUES,ierr)
258:         iend = n-1
259:       else
260:         iend = n
261:       endif
262:       do 10 i=istart,iend-1
263:         ii = i + start
264:         j = start + i - 1
265:         call MatSetValues(jac,1,ii,1,j,d,INSERT_VALUES,ierr)
266:         j = start + i + 1
267:         call MatSetValues(jac,1,ii,1,j,d,INSERT_VALUES,ierr)
268:         A = -2.0*d + 2.0*xx(i+1)
269:         call MatSetValues(jac,1,ii,1,ii,A,INSERT_VALUES,ierr)
270:  10   continue
271:       call VecRestoreArray(x,vxx,ixx,ierr)
272:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
273:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
274:       return
275:       end