Actual source code: ex12f.F

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: !
  2: !
  3: !  This example demonstrates basic use of the SNES Fortran interface.
  4: !
  5: !
  6:         module UserModule
  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: ! ---------------------------------------------------------------------
 24: !  Subroutine FormMonitor
 25: !  This function lets up keep track of the SNES progress at each step
 26: !  In this routine, we determine when the Jacobian is rebuilt with the parameter 'jag'
 27: !
 28: !  Input Parameters:
 29: !    snes    - SNES nonlinear solver context
 30: !    its     - current nonlinear iteration, starting from a call of SNESSolve()
 31: !    norm    - 2-norm of current residual (may be approximate)
 32: !    snesm - monctx designed module (included in Snesmmod)
 33: ! ---------------------------------------------------------------------
 34:       subroutine FormMonitor(snes,its,norm,snesm,ierr)
 35:       use UserModule
 36:       implicit none

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

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

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

 85:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 86:       if (ierr .ne. 0) then
 87:         print*,'Unable to initialize PETSc'
 88:         stop
 89:       endif
 90:       i1 = 1
 91:       i0 = 0
 92:       i3 = 3
 93:       N  = 10
 94:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,        &
 95:      &                        '-n',N,flg,ierr)
 96:       h = 1.0/real(N-1)
 97:       ctx%N = N
 98:       ctx%comm = PETSC_COMM_WORLD

100:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
101:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)

103: ! Set up data structures
104:       call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,i1,i1,            &
105:      &     PETSC_NULL_INTEGER,ctx%da,ierr)
106:       call DMSetFromOptions(ctx%da,ierr)
107:       call DMSetUp(ctx%da,ierr)
108:       call DMCreateGlobalVector(ctx%da,x,ierr)
109:       call DMCreateLocalVector(ctx%da,ctx%xl,ierr)

111:       call PetscObjectSetName(x,'Approximate Solution',ierr)
112:       call VecDuplicate(x,r,ierr)
113:       call VecDuplicate(x,ctx%F,ierr)
114:       call VecDuplicate(x,U,ierr)
115:       call PetscObjectSetName(U,'Exact Solution',ierr)

117:       call MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,           &
118:      &     N,i3,PETSC_NULL_INTEGER,i0,PETSC_NULL_INTEGER,J,ierr)
119:       call MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE,ierr)
120:       call MatGetType(J,matrixname,ierr)

122: ! Store right-hand-side of PDE and exact solution
123:       call VecGetOwnershipRange(x,start,end,ierr)
124:       xp = h*start
125:       nn = end - start
126:       ii = start
127:       do 10, i=0,nn-1
128:         FF = 6.0*xp + (xp+1.e-12)**6.e0
129:         UU = xp*xp*xp
130:         call VecSetValues(ctx%F,i1,ii,FF,INSERT_VALUES,ierr)
131:         call VecSetValues(U,i1,ii,UU,INSERT_VALUES,ierr)
132:         xp = xp + h
133:         ii = ii + 1
134:  10   continue
135:       call VecAssemblyBegin(ctx%F,ierr)
136:       call VecAssemblyEnd(ctx%F,ierr)
137:       call VecAssemblyBegin(U,ierr)
138:       call VecAssemblyEnd(U,ierr)

140: ! Create nonlinear solver
141:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

143: ! Set various routines and options
144:       call SNESSetFunction(snes,r,FormFunction,ctx,ierr)
145:       call SNESSetJacobian(snes,J,J,FormJacobian,ctx,ierr)

147:       snesm%its = 0
148:       call SNESGetLagJacobian(snes,snesm%lag,ierr)
149:       call SNESMonitorSet(snes,FormMonitor,snesm,                        &
150:      &        PETSC_NULL_FUNCTION,ierr)
151:       call SNESSetFromOptions(snes,ierr)

153: ! Solve nonlinear system
154:       call FormInitialGuess(snes,x,ierr)
155:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
156:       call SNESGetIterationNumber(snes,its,ierr);

158: !  Free work space.  All PETSc objects should be destroyed when they
159: !  are no longer needed.
160:       call VecDestroy(x,ierr)
161:       call VecDestroy(ctx%xl,ierr)
162:       call VecDestroy(r,ierr)
163:       call VecDestroy(U,ierr)
164:       call VecDestroy(ctx%F,ierr)
165:       call MatDestroy(J,ierr)
166:       call SNESDestroy(snes,ierr)
167:       call DMDestroy(ctx%da,ierr)
168:       call PetscFinalize(ierr)
169:       end


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

174:       subroutine FormFunction(snes,x,f,ctx,ierr)
175:       use UserModule
176:       implicit none
177:       SNES             snes
178:       Vec              x,f
179:       type(User) ctx
180:       PetscMPIInt  rank,size
181:       PetscInt i,s,n
182:       PetscErrorCode ierr
183:       PetscOffset      ixx,iff,iF2
184:       PetscScalar      h,d,vf2(2),vxx(2),vff(2)

186:       call MPI_Comm_rank(ctx%comm,rank,ierr)
187:       call MPI_Comm_size(ctx%comm,size,ierr)
188:       h     = 1.0/(real(ctx%N) - 1.0)
189:       call DMGlobalToLocalBegin(ctx%da,x,INSERT_VALUES,ctx%xl,ierr)
190:       call DMGlobalToLocalEnd(ctx%da,x,INSERT_VALUES,ctx%xl,ierr)

192:       call VecGetLocalSize(ctx%xl,n,ierr)
193:       if (n .gt. 1000) then
194:         print*, 'Local work array not big enough'
195:         call MPI_Abort(PETSC_COMM_WORLD,0,ierr)
196:       endif

198: !
199: ! This sets the index ixx so that vxx(ixx+1) is the first local
200: ! element in the vector indicated by ctx%xl.
201: !
202:       call VecGetArrayRead(ctx%xl,vxx,ixx,ierr)
203:       call VecGetArray(f,vff,iff,ierr)
204:       call VecGetArray(ctx%F,vF2,iF2,ierr)

206:       d = h*h

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

231:       do 10 i=1,n-2
232:        ff(i-s+1) = d*(xx(i) - 2.0*xx(i+1)                               &
233:      &      + xx(i+2)) + xx(i+1)*xx(i+1)                                &
234:      &      - F2(i-s+1)
235:  10   continue

237:       if (rank .eq. size-1) then
238:         ff(n-s) = xx(n) - 1.0
239:       endif

241:       call VecRestoreArray(f,vff,iff,ierr)
242:       call VecRestoreArrayRead(ctx%xl,vxx,ixx,ierr)
243:       call VecRestoreArray(ctx%F,vF2,iF2,ierr)
244:       return
245:       end

247: ! --------------------  Form initial approximation -----------------

249:       subroutine FormInitialGuess(snes,x,ierr)
250:       use UserModule
251:       implicit none

253:       PetscErrorCode   ierr
254:       Vec              x
255:       SNES             snes
256:       PetscScalar      five

258:       five = .5
259:       call VecSet(x,five,ierr)
260:       return
261:       end

263: ! --------------------  Evaluate Jacobian --------------------

265:       subroutine FormJacobian(snes,x,jac,B,ctx,ierr)
266:       use UserModule
267:       implicit none

269:       SNES             snes
270:       Vec              x
271:       Mat              jac,B
272:       type(User) ctx
273:       PetscInt  ii,istart,iend
274:       PetscInt  i,j,n,end,start,i1
275:       PetscErrorCode ierr
276:       PetscMPIInt rank,size
277:       PetscOffset      ixx
278:       PetscScalar      d,A,h,vxx(2)

280:       i1 = 1
281:       h = 1.0/(real(ctx%N) - 1.0)
282:       d = h*h
283:       call MPI_Comm_rank(ctx%comm,rank,ierr)
284:       call MPI_Comm_size(ctx%comm,size,ierr)

286:       call VecGetArrayRead(x,vxx,ixx,ierr)
287:       call VecGetOwnershipRange(x,start,end,ierr)
288:       n = end - start

290:       if (rank .eq. 0) then
291:         A = 1.0
292:         call MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr)
293:         istart = 1
294:       else
295:         istart = 0
296:       endif
297:       if (rank .eq. size-1) then
298:         i = INT(ctx%N-1)
299:         A = 1.0
300:         call MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr)
301:         iend = n-1
302:       else
303:         iend = n
304:       endif
305:       do 10 i=istart,iend-1
306:         ii = i + start
307:         j = start + i - 1
308:         call MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)
309:         j = start + i + 1
310:         call MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)
311:         A = -2.0*d + 2.0*xx(i+1)
312:         call MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr)
313:  10   continue
314:       call VecRestoreArrayRead(x,vxx,ixx,ierr)
315:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
316:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
317:       return
318:       end

320: !/*TEST
321: !
322: !   test:
323: !      nsize: 2
324: !      args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
325: !      output_file: output/ex12_1.out
326: !
327: !TEST*/