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*/