Actual source code: ex12f.F
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,zero
181: PetscInt i,s,n
182: PetscErrorCode ierr
183: PetscOffset ixx,iff,iF2
184: PetscScalar h,d,vf2(2),vxx(2),vff(2)
186: zero = 0
187: call MPI_Comm_rank(ctx%comm,rank,ierr)
188: call MPI_Comm_size(ctx%comm,size,ierr)
189: h = 1.0/(real(ctx%N) - 1.0)
190: call DMGlobalToLocalBegin(ctx%da,x,INSERT_VALUES,ctx%xl,ierr)
191: call DMGlobalToLocalEnd(ctx%da,x,INSERT_VALUES,ctx%xl,ierr)
193: call VecGetLocalSize(ctx%xl,n,ierr)
194: if (n .gt. 1000) then
195: print*, 'Local work array not big enough'
196: call MPI_Abort(PETSC_COMM_WORLD,zero,ierr)
197: endif
199: !
200: ! This sets the index ixx so that vxx(ixx+1) is the first local
201: ! element in the vector indicated by ctx%xl.
202: !
203: call VecGetArrayRead(ctx%xl,vxx,ixx,ierr)
204: call VecGetArray(f,vff,iff,ierr)
205: call VecGetArray(ctx%F,vF2,iF2,ierr)
207: d = h*h
209: !
210: ! Note that the array vxx() was obtained from a ghosted local vector
211: ! ctx%xl while the array vff() was obtained from the non-ghosted parallel
212: ! vector F. This is why there is a need for shift variable s. Since vff()
213: ! does not have locations for the ghost variables we need to index in it
214: ! slightly different then indexing into vxx(). For example on processor
215: ! 1 (the second processor)
216: !
217: ! xx(1) xx(2) xx(3) .....
218: ! ^^^^^^^ ^^^^^ ^^^^^
219: ! ghost value 1st local value 2nd local value
220: !
221: ! ff(1) ff(2)
222: ! ^^^^^^^ ^^^^^^^
223: ! 1st local value 2nd local value
224: !
225: if (rank .eq. 0) then
226: s = 0
227: ff(1) = xx(1)
228: else
229: s = 1
230: endif
232: do 10 i=1,n-2
233: ff(i-s+1) = d*(xx(i) - 2.0*xx(i+1) &
234: & + xx(i+2)) + xx(i+1)*xx(i+1) &
235: & - F2(i-s+1)
236: 10 continue
238: if (rank .eq. size-1) then
239: ff(n-s) = xx(n) - 1.0
240: endif
242: call VecRestoreArray(f,vff,iff,ierr)
243: call VecRestoreArrayRead(ctx%xl,vxx,ixx,ierr)
244: call VecRestoreArray(ctx%F,vF2,iF2,ierr)
245: return
246: end
248: ! -------------------- Form initial approximation -----------------
250: subroutine FormInitialGuess(snes,x,ierr)
251: use UserModule
252: implicit none
254: PetscErrorCode ierr
255: Vec x
256: SNES snes
257: PetscScalar five
259: five = .5
260: call VecSet(x,five,ierr)
261: return
262: end
264: ! -------------------- Evaluate Jacobian --------------------
266: subroutine FormJacobian(snes,x,jac,B,ctx,ierr)
267: use UserModule
268: implicit none
270: SNES snes
271: Vec x
272: Mat jac,B
273: type(User) ctx
274: PetscInt ii,istart,iend
275: PetscInt i,j,n,end,start,i1
276: PetscErrorCode ierr
277: PetscMPIInt rank,size
278: PetscOffset ixx
279: PetscScalar d,A,h,vxx(2)
281: i1 = 1
282: h = 1.0/(real(ctx%N) - 1.0)
283: d = h*h
284: call MPI_Comm_rank(ctx%comm,rank,ierr)
285: call MPI_Comm_size(ctx%comm,size,ierr)
287: call VecGetArrayRead(x,vxx,ixx,ierr)
288: call VecGetOwnershipRange(x,start,end,ierr)
289: n = end - start
291: if (rank .eq. 0) then
292: A = 1.0
293: call MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr)
294: istart = 1
295: else
296: istart = 0
297: endif
298: if (rank .eq. size-1) then
299: i = INT(ctx%N-1)
300: A = 1.0
301: call MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr)
302: iend = n-1
303: else
304: iend = n
305: endif
306: do 10 i=istart,iend-1
307: ii = i + start
308: j = start + i - 1
309: call MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)
310: j = start + i + 1
311: call MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr)
312: A = -2.0*d + 2.0*xx(i+1)
313: call MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr)
314: 10 continue
315: call VecRestoreArrayRead(x,vxx,ixx,ierr)
316: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
317: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
318: return
319: end
321: !/*TEST
322: !
323: ! test:
324: ! nsize: 2
325: ! args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
326: ! output_file: output/ex12_1.out
327: !
328: !TEST*/