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,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,i0,PETSC_NULL_INTEGER,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: return
225: end
227: ! -------------------- Form initial approximation -----------------
229: subroutine FormInitialGuess(snes,x,ierr)
230: use ex12fmodule
231: implicit none
233: PetscErrorCode ierr
234: Vec x
235: SNES snes
236: PetscScalar five
238: five = .5
239: PetscCall(VecSet(x,five,ierr))
240: return
241: end
243: ! -------------------- Evaluate Jacobian --------------------
245: subroutine FormJacobian(snes,x,jac,B,ctx,ierr)
246: use ex12fmodule
247: implicit none
249: SNES snes
250: Vec x
251: Mat jac,B
252: type(User) ctx
253: PetscInt ii,istart,iend
254: PetscInt i,j,n,end,start,i1
255: PetscErrorCode ierr
256: PetscMPIInt rank,size
257: PetscScalar d,A,h
258: PetscScalar,pointer :: vxx(:)
260: i1 = 1
261: h = 1.0/(real(ctx%N) - 1.0)
262: d = h*h
263: PetscCallMPI(MPI_Comm_rank(ctx%comm,rank,ierr))
264: PetscCallMPI(MPI_Comm_size(ctx%comm,size,ierr))
266: PetscCall(VecGetArrayReadF90(x,vxx,ierr))
267: PetscCall(VecGetOwnershipRange(x,start,end,ierr))
268: n = end - start
270: if (rank .eq. 0) then
271: A = 1.0
272: PetscCall(MatSetValues(jac,i1,start,i1,start,A,INSERT_VALUES,ierr))
273: istart = 1
274: else
275: istart = 0
276: endif
277: if (rank .eq. size-1) then
278: i = INT(ctx%N-1)
279: A = 1.0
280: PetscCall(MatSetValues(jac,i1,i,i1,i,A,INSERT_VALUES,ierr))
281: iend = n-1
282: else
283: iend = n
284: endif
285: do 10 i=istart,iend-1
286: ii = i + start
287: j = start + i - 1
288: PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
289: j = start + i + 1
290: PetscCall(MatSetValues(jac,i1,ii,i1,j,d,INSERT_VALUES,ierr))
291: A = -2.0*d + 2.0*vxx(i+1)
292: PetscCall(MatSetValues(jac,i1,ii,i1,ii,A,INSERT_VALUES,ierr))
293: 10 continue
294: PetscCall(VecRestoreArrayReadF90(x,vxx,ierr))
295: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
296: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
297: return
298: end
300: !/*TEST
301: !
302: ! test:
303: ! nsize: 2
304: ! args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
305: ! output_file: output/ex12_1.out
306: !
307: !TEST*/