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: !
 60: !
 61: ! Macros to make setting/getting  values into vector clearer.
 62: ! The element xx(ib) is the ibth element in the vector indicated by ctx%F
 63: #define xx(ib)  vxx(ixx + (ib))
 64: #define ff(ib)  vff(iff + (ib))
 65: #define F2(ib)  vF2(iF2 + (ib))
 66:       program main
 67:       use ex12fmodule
 68:       implicit none
 69:       type(User) ctx
 70:       PetscMPIInt rank,size
 71:       PetscErrorCode ierr
 72:       PetscInt N,start,end,nn,i
 73:       PetscInt ii,its,i1,i0,i3
 74:       PetscBool  flg
 75:       SNES             snes
 76:       Mat              J
 77:       Vec              x,r,u
 78:       PetscScalar      xp,FF,UU,h
 79:       character*(10)   matrixname
 80:       external         FormJacobian,FormFunction
 81:       external         formmonitor
 82:       type(monctx) :: snesm

 84:       PetscCallA(PetscInitialize(ierr))
 85:       i1 = 1
 86:       i0 = 0
 87:       i3 = 3
 88:       N  = 10
 89:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',N,flg,ierr))
 90:       h = 1.0/real(N-1)
 91:       ctx%N = N
 92:       ctx%comm = PETSC_COMM_WORLD

 94:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 95:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))

 97: ! Set up data structures
 98:       PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,i1,i1,PETSC_NULL_INTEGER,ctx%da,ierr))
 99:       PetscCallA(DMSetFromOptions(ctx%da,ierr))
100:       PetscCallA(DMSetUp(ctx%da,ierr))
101:       PetscCallA(DMCreateGlobalVector(ctx%da,x,ierr))
102:       PetscCallA(DMCreateLocalVector(ctx%da,ctx%xl,ierr))

104:       PetscCallA(PetscObjectSetName(x,'Approximate Solution',ierr))
105:       PetscCallA(VecDuplicate(x,r,ierr))
106:       PetscCallA(VecDuplicate(x,ctx%F,ierr))
107:       PetscCallA(VecDuplicate(x,U,ierr))
108:       PetscCallA(PetscObjectSetName(U,'Exact Solution',ierr))

110:       PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,i3,PETSC_NULL_INTEGER,i0,PETSC_NULL_INTEGER,J,ierr))
111:       PetscCallA(MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE,ierr))
112:       PetscCallA(MatGetType(J,matrixname,ierr))

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

132: ! Create nonlinear solver
133:       PetscCallA(SNESCreate(PETSC_COMM_WORLD,snes,ierr))

135: ! Set various routines and options
136:       PetscCallA(SNESSetFunction(snes,r,FormFunction,ctx,ierr))
137:       PetscCallA(SNESSetJacobian(snes,J,J,FormJacobian,ctx,ierr))

139:       snesm%its = 0
140:       PetscCallA(SNESGetLagJacobian(snes,snesm%lag,ierr))
141:       PetscCallA(SNESMonitorSet(snes,FormMonitor,snesm,PETSC_NULL_FUNCTION,ierr))
142:       PetscCallA(SNESSetFromOptions(snes,ierr))

144: ! Solve nonlinear system
145:       PetscCallA(FormInitialGuess(snes,x,ierr))
146:       PetscCallA(SNESSolve(snes,PETSC_NULL_VEC,x,ierr))
147:       PetscCallA(SNESGetIterationNumber(snes,its,ierr))

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

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

164:       subroutine FormFunction(snes,x,f,ctx,ierr)
165:       use ex12fmodule
166:       implicit none
167:       SNES             snes
168:       Vec              x,f
169:       type(User) ctx
170:       PetscMPIInt  rank,size,zero
171:       PetscInt i,s,n
172:       PetscErrorCode ierr
173:       PetscOffset      ixx,iff,iF2
174:       PetscScalar      h,d,vf2(2),vxx(2),vff(2)

176:       zero = 0
177:       PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr))
178:       PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr))
179:       h     = 1.0/(real(ctx%N) - 1.0)
180:       PetscCall(DMGlobalToLocalBegin(ctx%da,x,INSERT_VALUES,ctx%xl,ierr))
181:       PetscCall(DMGlobalToLocalEnd(ctx%da,x,INSERT_VALUES,ctx%xl,ierr))

183:       PetscCall(VecGetLocalSize(ctx%xl,n,ierr))
184:       if (n .gt. 1000) then
185:         print*, 'Local work array not big enough'
186:         call MPI_Abort(PETSC_COMM_WORLD,zero,ierr)
187:       endif

189: !
190: ! This sets the index ixx so that vxx(ixx+1) is the first local
191: ! element in the vector indicated by ctx%xl.
192: !
193:       PetscCall(VecGetArrayRead(ctx%xl,vxx,ixx,ierr))
194:       PetscCall(VecGetArray(f,vff,iff,ierr))
195:       PetscCall(VecGetArray(ctx%F,vF2,iF2,ierr))

197:       d = h*h

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

222:       do 10 i=1,n-2
223:        ff(i-s+1) = d*(xx(i) - 2.0*xx(i+1) + xx(i+2)) + xx(i+1)*xx(i+1) - F2(i-s+1)
224:  10   continue

226:       if (rank .eq. size-1) then
227:         ff(n-s) = xx(n) - 1.0
228:       endif

230:       PetscCall(VecRestoreArray(f,vff,iff,ierr))
231:       PetscCall(VecRestoreArrayRead(ctx%xl,vxx,ixx,ierr))
232:       PetscCall(VecRestoreArray(ctx%F,vF2,iF2,ierr))
233:       return
234:       end

236: ! --------------------  Form initial approximation -----------------

238:       subroutine FormInitialGuess(snes,x,ierr)
239:       use ex12fmodule
240:       implicit none

242:       PetscErrorCode   ierr
243:       Vec              x
244:       SNES             snes
245:       PetscScalar      five

247:       five = .5
248:       PetscCall(VecSet(x,five,ierr))
249:       return
250:       end

252: ! --------------------  Evaluate Jacobian --------------------

254:       subroutine FormJacobian(snes,x,jac,B,ctx,ierr)
255:       use ex12fmodule
256:       implicit none

258:       SNES             snes
259:       Vec              x
260:       Mat              jac,B
261:       type(User) ctx
262:       PetscInt  ii,istart,iend
263:       PetscInt  i,j,n,end,start,i1
264:       PetscErrorCode ierr
265:       PetscMPIInt rank,size
266:       PetscOffset      ixx
267:       PetscScalar      d,A,h,vxx(2)

269:       i1 = 1
270:       h = 1.0/(real(ctx%N) - 1.0)
271:       d = h*h
272:       PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr))
273:       PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr))

275:       PetscCall(VecGetArrayRead(x,vxx,ixx,ierr))
276:       PetscCall(VecGetOwnershipRange(x,start,end,ierr))
277:       n = end - start

279:       if (rank .eq. 0) then
280:         A = 1.0
281:         PetscCall(MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr))
282:         istart = 1
283:       else
284:         istart = 0
285:       endif
286:       if (rank .eq. size-1) then
287:         i = INT(ctx%N-1)
288:         A = 1.0
289:         PetscCall(MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr))
290:         iend = n-1
291:       else
292:         iend = n
293:       endif
294:       do 10 i=istart,iend-1
295:         ii = i + start
296:         j = start + i - 1
297:         PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
298:         j = start + i + 1
299:         PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
300:         A = -2.0*d + 2.0*xx(i+1)
301:         PetscCall(MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr))
302:  10   continue
303:       PetscCall(VecRestoreArrayRead(x,vxx,ixx,ierr))
304:       PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
305:       PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
306:       return
307:       end

309: !/*TEST
310: !
311: !   test:
312: !      nsize: 2
313: !      args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
314: !      output_file: output/ex12_1.out
315: !
316: !TEST*/