Actual source code: ex12f.F90

  1: !
  2: !
  3: !  This example demonstrates basic use of the SNES Fortran interface.
  4: !
  5: !
  6:         module ex12fmodule
  7: #include <petsc/finclude/petscsnes.h>
  8:         use petscsnes
  9:         type User
 10:           DM  da
 11:           Vec F
 12:           Vec xl
 13:           MPI_Comm comm
 14:           PetscInt N
 15:         end type User
 16:         save
 17:         type monctx
 18:         PetscInt :: its,lag
 19:         end type monctx
 20:       end module

 22: ! ---------------------------------------------------------------------
 23: !  Subroutine FormMonitor
 24: !  This function lets up keep track of the SNES progress at each step
 25: !  In this routine, we determine when the Jacobian is rebuilt with the parameter 'jag'
 26: !
 27: !  Input Parameters:
 28: !    snes    - SNES nonlinear solver context
 29: !    its     - current nonlinear iteration, starting from a call of SNESSolve()
 30: !    norm    - 2-norm of current residual (may be approximate)
 31: !    snesm - monctx designed module (included in Snesmmod)
 32: ! ---------------------------------------------------------------------
 33:       subroutine FormMonitor(snes,its,norm,snesm,ierr)
 34:       use ex12fmodule
 35:       implicit none

 37:       SNES ::           snes
 38:       PetscInt ::       its,one,mone
 39:       PetscScalar ::    norm
 40:       type(monctx) ::   snesm
 41:       PetscErrorCode :: ierr

 43: !      write(6,*) ' '
 44: !      write(6,*) '    its ',its,snesm%its,'lag',
 45: !     &            snesm%lag
 46: !      call flush(6)
 47:       if (mod(snesm%its,snesm%lag).eq.0) then
 48:         one = 1
 49:         PetscCall(SNESSetLagJacobian(snes,one,ierr))  ! build jacobian
 50:       else
 51:         mone = -1
 52:         PetscCall(SNESSetLagJacobian(snes,mone,ierr)) ! do NOT build jacobian
 53:       endif
 54:       snesm%its = snesm%its + 1
 55:       end subroutine FormMonitor

 57: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 58: !  MUST be declared as external.
 59: !

 61:       program main
 62:       use ex12fmodule
 63:       implicit none
 64:       type(User) ctx
 65:       PetscMPIInt rank,size
 66:       PetscErrorCode ierr
 67:       PetscInt N,start,end,nn,i
 68:       PetscInt ii,its,i1,i0,i3
 69:       PetscBool  flg
 70:       SNES             snes
 71:       Mat              J
 72:       Vec              x,r,u
 73:       PetscScalar      xp,FF,UU,h
 74:       character*(10)   matrixname
 75:       external         FormJacobian,FormFunction
 76:       external         formmonitor
 77:       type(monctx) :: snesm

 79:       PetscCallA(PetscInitialize(ierr))
 80:       i1 = 1
 81:       i0 = 0
 82:       i3 = 3
 83:       N  = 10
 84:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',N,flg,ierr))
 85:       h = 1.0/real(N-1)
 86:       ctx%N = N
 87:       ctx%comm = PETSC_COMM_WORLD

 89:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 90:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))

 92: ! Set up data structures
 93:       PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,i1,i1,PETSC_NULL_INTEGER,ctx%da,ierr))
 94:       PetscCallA(DMSetFromOptions(ctx%da,ierr))
 95:       PetscCallA(DMSetUp(ctx%da,ierr))
 96:       PetscCallA(DMCreateGlobalVector(ctx%da,x,ierr))
 97:       PetscCallA(DMCreateLocalVector(ctx%da,ctx%xl,ierr))

 99:       PetscCallA(PetscObjectSetName(x,'Approximate Solution',ierr))
100:       PetscCallA(VecDuplicate(x,r,ierr))
101:       PetscCallA(VecDuplicate(x,ctx%F,ierr))
102:       PetscCallA(VecDuplicate(x,U,ierr))
103:       PetscCallA(PetscObjectSetName(U,'Exact Solution',ierr))

105:       PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,i3,PETSC_NULL_INTEGER,i0,PETSC_NULL_INTEGER,J,ierr))
106:       PetscCallA(MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE,ierr))
107:       PetscCallA(MatGetType(J,matrixname,ierr))

109: ! Store right-hand-side of PDE and exact solution
110:       PetscCallA(VecGetOwnershipRange(x,start,end,ierr))
111:       xp = h*start
112:       nn = end - start
113:       ii = start
114:       do 10, i=0,nn-1
115:         FF = 6.0*xp + (xp+1.e-12)**6.e0
116:         UU = xp*xp*xp
117:         PetscCallA(VecSetValues(ctx%F,i1,ii,FF,INSERT_VALUES,ierr))
118:         PetscCallA(VecSetValues(U,i1,ii,UU,INSERT_VALUES,ierr))
119:         xp = xp + h
120:         ii = ii + 1
121:  10   continue
122:       PetscCallA(VecAssemblyBegin(ctx%F,ierr))
123:       PetscCallA(VecAssemblyEnd(ctx%F,ierr))
124:       PetscCallA(VecAssemblyBegin(U,ierr))
125:       PetscCallA(VecAssemblyEnd(U,ierr))

127: ! Create nonlinear solver
128:       PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))

130: ! Set various routines and options
131:       PetscCallA(SNESSetFunction(snes,r,FormFunction,ctx,ierr))
132:       PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,ctx,ierr))

134:       snesm%its = 0
135:       PetscCallA(SNESGetLagJacobian(snes,snesm%lag,ierr))
136:       PetscCallA(SNESMonitorSet(snes,FormMonitor,snesm,PETSC_NULL_FUNCTION,ierr))
137:       PetscCallA(SNESSetFromOptions(snes,ierr))

139: ! Solve nonlinear system
140:       PetscCallA(FormInitialGuess(snes,x,ierr))
141:       PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
142:       PetscCallA(SNESGetIterationNumber(snes,its,ierr))

144: !  Free work space.  All PETSc objects should be destroyed when they
145: !  are no longer needed.
146:       PetscCallA(VecDestroy(x,ierr))
147:       PetscCallA(VecDestroy(ctx%xl,ierr))
148:       PetscCallA(VecDestroy(r,ierr))
149:       PetscCallA(VecDestroy(U,ierr))
150:       PetscCallA(VecDestroy(ctx%F,ierr))
151:       PetscCallA(MatDestroy(J,ierr))
152:       PetscCallA(SNESDestroy(snes,ierr))
153:       PetscCallA(DMDestroy(ctx%da,ierr))
154:       PetscCallA(PetscFinalize(ierr))
155:       end

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

159:       subroutine FormFunction(snes,x,f,ctx,ierr)
160:       use ex12fmodule
161:       implicit none
162:       SNES             snes
163:       Vec              x,f
164:       type(User) ctx
165:       PetscMPIInt  rank,size,zero
166:       PetscInt i,s,n
167:       PetscErrorCode ierr
168:       PetscScalar      h,d
169:       PetscScalar,pointer :: vf2(:),vxx(:),vff(:)

171:       zero = 0
172:       PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr))
173:       PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr))
174:       h     = 1.0/(real(ctx%N) - 1.0)
175:       PetscCall(DMGlobalToLocalBegin(ctx%da,x,INSERT_VALUES,ctx%xl,ierr))
176:       PetscCall(DMGlobalToLocalEnd(ctx%da,x,INSERT_VALUES,ctx%xl,ierr))

178:       PetscCall(VecGetLocalSize(ctx%xl,n,ierr))
179:       if (n .gt. 1000) then
180:         print*, 'Local work array not big enough'
181:         call MPI_Abort(PETSC_COMM_WORLD,zero,ierr)
182:       endif

184:       PetscCall(VecGetArrayReadF90(ctx%xl,vxx,ierr))
185:       PetscCall(VecGetArrayF90(f,vff,ierr))
186:       PetscCall(VecGetArrayF90(ctx%F,vF2,ierr))

188:       d = h*h

190: !
191: !  Note that the array vxx() was obtained from a ghosted local vector
192: !  ctx%xl while the array vff() was obtained from the non-ghosted parallel
193: !  vector F. This is why there is a need for shift variable s. Since vff()
194: !  does not have locations for the ghost variables we need to index in it
195: !  slightly different then indexing into vxx(). For example on processor
196: !  1 (the second processor)
197: !
198: !        xx(1)        xx(2)             xx(3)             .....
199: !      ^^^^^^^        ^^^^^             ^^^^^
200: !      ghost value   1st local value   2nd local value
201: !
202: !                      ff(1)             ff(2)
203: !                     ^^^^^^^           ^^^^^^^
204: !                    1st local value   2nd local value
205: !
206:        if (rank .eq. 0) then
207:         s = 0
208:         vff(1) = vxx(1)
209:       else
210:         s = 1
211:       endif

213:       do 10 i=1,n-2
214:        vff(i-s+1) = d*(vxx(i) - 2.0*vxx(i+1) + vxx(i+2)) + vxx(i+1)*vxx(i+1) - vF2(i-s+1)
215:  10   continue

217:       if (rank .eq. size-1) then
218:         vff(n-s) = vxx(n) - 1.0
219:       endif

221:       PetscCall(VecRestoreArrayF90(f,vff,ierr))
222:       PetscCall(VecRestoreArrayReadF90(ctx%xl,vxx,ierr))
223:       PetscCall(VecRestoreArrayF90(ctx%F,vF2,ierr))
224:       return
225:       end

227: ! --------------------  Form initial approximation -----------------

229:       subroutine FormInitialGuess(snes,x,ierr)
230:       use ex12fmodule
231:       implicit none

233:       PetscErrorCode   ierr
234:       Vec              x
235:       SNES             snes
236:       PetscScalar      five

238:       five = .5
239:       PetscCall(VecSet(x,five,ierr))
240:       return
241:       end

243: ! --------------------  Evaluate Jacobian --------------------

245:       subroutine FormJacobian(snes,x,jac,B,ctx,ierr)
246:       use ex12fmodule
247:       implicit none

249:       SNES             snes
250:       Vec              x
251:       Mat              jac,B
252:       type(User) ctx
253:       PetscInt  ii,istart,iend
254:       PetscInt  i,j,n,end,start,i1
255:       PetscErrorCode ierr
256:       PetscMPIInt rank,size
257:       PetscScalar      d,A,h
258:       PetscScalar,pointer :: vxx(:)

260:       i1 = 1
261:       h = 1.0/(real(ctx%N) - 1.0)
262:       d = h*h
263:       PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr))
264:       PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr))

266:       PetscCall(VecGetArrayReadF90(x,vxx,ierr))
267:       PetscCall(VecGetOwnershipRange(x,start,end,ierr))
268:       n = end - start

270:       if (rank .eq. 0) then
271:         A = 1.0
272:         PetscCall(MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr))
273:         istart = 1
274:       else
275:         istart = 0
276:       endif
277:       if (rank .eq. size-1) then
278:         i = INT(ctx%N-1)
279:         A = 1.0
280:         PetscCall(MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr))
281:         iend = n-1
282:       else
283:         iend = n
284:       endif
285:       do 10 i=istart,iend-1
286:         ii = i + start
287:         j = start + i - 1
288:         PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
289:         j = start + i + 1
290:         PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
291:         A = -2.0*d + 2.0*vxx(i+1)
292:         PetscCall(MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr))
293:  10   continue
294:       PetscCall(VecRestoreArrayReadF90(x,vxx,ierr))
295:       PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
296:       PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
297:       return
298:       end

300: !/*TEST
301: !
302: !   test:
303: !      nsize: 2
304: !      args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
305: !      output_file: output/ex12_1.out
306: !
307: !TEST*/