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_ARRAY,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_ARRAY,i0,PETSC_NULL_INTEGER_ARRAY,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:       end

226: ! --------------------  Form initial approximation -----------------

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

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

237:       five = .5
238:       PetscCall(VecSet(x,five,ierr))
239:       end

241: ! --------------------  Evaluate Jacobian --------------------

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

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

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

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

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

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