Actual source code: ex12f.F
petsc-3.13.4 2020-08-01
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*/