Actual source code: ex14f.F
1: !
2: !
3: ! Solves a nonlinear system in parallel with a user-defined
4: ! Newton method that uses KSP to solve the linearized Newton sytems. This solver
5: ! is a very simplistic inexact Newton method. The intent of this code is to
6: ! demonstrate the repeated solution of linear sytems with the same nonzero pattern.
7: !
8: ! This is NOT the recommended approach for solving nonlinear problems with PETSc!
9: ! We urge users to employ the SNES component for solving nonlinear problems whenever
10: ! possible, as it offers many advantages over coding nonlinear solvers independently.
11: !
12: ! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
13: ! domain, using distributed arrays (DAs) to partition the parallel grid.
14: !
15: ! The command line options include:
16: ! -par <parameter>, where <parameter> indicates the problem's nonlinearity
17: ! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)
18: ! -mx <xg>, where <xg> = number of grid points in the x-direction
19: ! -my <yg>, where <yg> = number of grid points in the y-direction
20: ! -Nx <npx>, where <npx> = number of processors in the x-direction
21: ! -Ny <npy>, where <npy> = number of processors in the y-direction
22: ! -mf use matrix free for matrix vector product
23: !
24: !/*T
25: ! Concepts: KSP^writing a user-defined nonlinear solver
26: ! Concepts: DA^using distributed arrays
27: ! Processors: n
28: !T*/
29: ! ------------------------------------------------------------------------
30: !
31: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
32: ! the partial differential equation
33: !
34: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
35: !
36: ! with boundary conditions
37: !
38: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
39: !
40: ! A finite difference approximation with the usual 5-point stencil
41: ! is used to discretize the boundary value problem to obtain a nonlinear
42: ! system of equations.
43: !
44: ! The SNES version of this problem is: snes/examples/tutorials/ex5f.F
45: !
46: ! -------------------------------------------------------------------------
48: program main
49: implicit none
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: ! Include files
53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: !
55: ! petsc.h - base PETSc routines petscvec.h - vectors
56: ! petscsys.h - system routines petscmat.h - matrices
57: ! petscis.h - index sets petscksp.h - Krylov subspace methods
58: ! petscviewer.h - viewers petscpc.h - preconditioners
60: #include finclude/petsc.h
61: #include finclude/petscis.h
62: #include finclude/petscvec.h
63: #include finclude/petscmat.h
64: #include finclude/petscpc.h
65: #include finclude/petscksp.h
66: #include finclude/petscda.h
68: MPI_Comm comm
69: Vec X,Y,F,localX,localF
70: Mat J,B
71: DA da
72: KSP ksp
74: PetscInt Nx,Ny,N,mx,my,ifive,ithree
75: PetscTruth flg,nooutput,usemf
76: common /mycommon/ mx,my,B,localX,localF,da
77: !
78: !
79: ! This is the routine to use for matrix-free approach
80: !
81: external mymult
83: ! --------------- Data to define nonlinear solver --------------
84: double precision rtol,ttol
85: double precision fnorm,ynorm,xnorm
86: PetscInt max_nonlin_its,one
87: PetscInt lin_its
88: PetscInt i,m
89: PetscScalar mone
90: PetscErrorCode ierr
92: mone = -1.d0
93: rtol = 1.d-8
94: max_nonlin_its = 10
95: one = 1
96: ifive = 5
97: ithree = 3
99: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
100: comm = PETSC_COMM_WORLD
102: ! Initialize problem parameters
104: !
105: mx = 4
106: my = 4
107: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
108: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
109: N = mx*my
111: nooutput = .false.
112: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-no_output', &
113: & nooutput,ierr)
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: ! Create linear solver context
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: call KSPCreate(comm,ksp,ierr)
121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: ! Create vector data structures
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: !
126: ! Create distributed array (DA) to manage parallel grid and vectors
127: !
128: Nx = PETSC_DECIDE
129: Ny = PETSC_DECIDE
130: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-Nx',Nx,flg,ierr)
131: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-Ny',Ny,flg,ierr)
132: call DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,mx, &
133: & my,Nx,Ny,one,one,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
134: & da,ierr)
136: !
137: ! Extract global and local vectors from DA then duplicate for remaining
138: ! vectors that are the same types
139: !
140: call DACreateGlobalVector(da,X,ierr)
141: call DACreateLocalVector(da,localX,ierr)
142: call VecDuplicate(X,F,ierr)
143: call VecDuplicate(X,Y,ierr)
144: call VecDuplicate(localX,localF,ierr)
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: ! Create matrix data structure for Jacobian
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: !
151: ! Note: For the parallel case, vectors and matrices MUST be partitioned
152: ! accordingly. When using distributed arrays (DAs) to create vectors,
153: ! the DAs determine the problem partitioning. We must explicitly
154: ! specify the local matrix dimensions upon its creation for compatibility
155: ! with the vector distribution.
156: !
157: ! Note: Here we only approximately preallocate storage space for the
158: ! Jacobian. See the users manual for a discussion of better techniques
159: ! for preallocating matrix memory.
160: !
161: call VecGetLocalSize(X,m,ierr)
162: call MatCreateMPIAIJ(comm,m,m,N,N,ifive,PETSC_NULL_INTEGER,ithree, &
163: & PETSC_NULL_INTEGER,B,ierr)
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: ! if usemf is on then matrix vector product is done via matrix free
167: ! approach. Note this is just an example, and not realistic because
168: ! we still use the actual formed matrix, but in reality one would
169: ! provide their own subroutine that would directly do the matrix
170: ! vector product and not call MatMult()
171: ! Note: we put B into a common block so it will be visible to the
172: ! mymult() routine
173: usemf = .false.
174: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-mf',usemf,ierr)
175: if (usemf) then
176: call MatCreateShell(comm,m,m,N,N,PETSC_NULL_INTEGER,J,ierr)
177: call MatShellSetOperation(J,MATOP_MULT,mymult,ierr)
178: else
179: ! If not doing matrix free then matrix operator, J, and matrix used
180: ! to construct preconditioner, B, are the same
181: J = B
182: endif
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: ! Customize linear solver set runtime options
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: !
188: ! Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
189: !
190: call KSPSetFromOptions(ksp,ierr)
192: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: ! Evaluate initial guess
194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: call FormInitialGuess(X,ierr)
197: call ComputeFunction(X,F,ierr)
198: call VecNorm(F,NORM_2,fnorm,ierr)
199: ttol = fnorm*rtol
200: if (.not. nooutput) then
201: print*, 'Initial function norm ',fnorm
202: endif
204: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: ! Solve nonlinear system with a user-defined method
206: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: ! This solver is a very simplistic inexact Newton method, with no
209: ! no damping strategies or bells and whistles. The intent of this code
210: ! is merely to demonstrate the repeated solution with KSP of linear
211: ! sytems with the same nonzero structure.
212: !
213: ! This is NOT the recommended approach for solving nonlinear problems
214: ! with PETSc! We urge users to employ the SNES component for solving
215: ! nonlinear problems whenever possible with application codes, as it
216: ! offers many advantages over coding nonlinear solvers independently.
218: do 10 i=0,max_nonlin_its
220: ! Compute the Jacobian matrix. See the comments in this routine for
221: ! important information about setting the flag mat_flag.
223: call ComputeJacobian(X,B,ierr)
225: ! Solve J Y = F, where J is the Jacobian matrix.
226: ! - First, set the KSP linear operators. Here the matrix that
227: ! defines the linear system also serves as the preconditioning
228: ! matrix.
229: ! - Then solve the Newton system.
231: call KSPSetOperators(ksp,J,B,SAME_NONZERO_PATTERN,ierr)
232: call KSPSolve(ksp,F,Y,ierr)
234: ! Compute updated iterate
236: call VecNorm(Y,NORM_2,ynorm,ierr)
237: call VecAYPX(Y,mone,X,ierr)
238: call VecCopy(Y,X,ierr)
239: call VecNorm(X,NORM_2,xnorm,ierr)
240: call KSPGetIterationNumber(ksp,lin_its,ierr)
241: if (.not. nooutput) then
242: print*,'linear solve iterations = ',lin_its,' xnorm = ', &
243: & xnorm,' ynorm = ',ynorm
244: endif
246: ! Evaluate nonlinear function at new location
248: call ComputeFunction(X,F,ierr)
249: call VecNorm(F,NORM_2,fnorm,ierr)
250: if (.not. nooutput) then
251: print*, 'Iteration ',i+1,' function norm',fnorm
252: endif
254: ! Test for convergence
256: if (fnorm .le. ttol) then
257: if (.not. nooutput) then
258: print*,'Converged: function norm ',fnorm,' tolerance ',ttol
259: endif
260: goto 20
261: endif
262: 10 continue
263: 20 continue
265: write(6,100) i+1
266: 100 format('Number of Newton iterations =',I2)
268: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: ! Free work space. All PETSc objects should be destroyed when they
270: ! are no longer needed.
271: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
273: call MatDestroy(B,ierr)
274: if (usemf) then
275: call MatDestroy(J,ierr)
276: endif
277: call VecDestroy(localX,ierr)
278: call VecDestroy(X,ierr)
279: call VecDestroy(Y,ierr)
280: call VecDestroy(localF,ierr)
281: call VecDestroy(F,ierr)
282: call KSPDestroy(ksp,ierr)
283: call DADestroy(da,ierr)
284: call PetscFinalize(ierr)
285: end
287: ! -------------------------------------------------------------------
288: !
289: ! FormInitialGuess - Forms initial approximation.
290: !
291: ! Input Parameters:
292: ! X - vector
293: !
294: ! Output Parameter:
295: ! X - vector
296: !
297: subroutine FormInitialGuess(X,ierr)
298: implicit none
300: ! petsc.h - base PETSc routines petscvec.h - vectors
301: ! petscsys.h - system routines petscmat.h - matrices
302: ! petscis.h - index sets petscksp.h - Krylov subspace methods
303: ! petscviewer.h - viewers petscpc.h - preconditioners
305: #include finclude/petsc.h
306: #include finclude/petscis.h
307: #include finclude/petscvec.h
308: #include finclude/petscmat.h
309: #include finclude/petscpc.h
310: #include finclude/petscksp.h
311: #include finclude/petscda.h
312: PetscErrorCode ierr
313: PetscOffset idx
314: Vec X,localX,localF
315: PetscInt i,j,row,mx
316: PetscInt my, xs,ys,xm
317: PetscInt ym,gxm,gym,gxs,gys
318: double precision one,lambda,temp1,temp,hx,hy
319: PetscScalar xx(1)
320: DA da
321: Mat B
322: common /mycommon/ mx,my,B,localX,localF,da
323:
324: one = 1.d0
325: lambda = 6.d0
326: hx = one/(mx-1)
327: hy = one/(my-1)
328: temp1 = lambda/(lambda + one)
330: ! Get a pointer to vector data.
331: ! - VecGetArray() returns a pointer to the data array.
332: ! - You MUST call VecRestoreArray() when you no longer need access to
333: ! the array.
334: call VecGetArray(localX,xx,idx,ierr)
336: ! Get local grid boundaries (for 2-dimensional DA):
337: ! xs, ys - starting grid indices (no ghost points)
338: ! xm, ym - widths of local grid (no ghost points)
339: ! gxs, gys - starting grid indices (including ghost points)
340: ! gxm, gym - widths of local grid (including ghost points)
342: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
343: & PETSC_NULL_INTEGER,ierr)
344: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
345: & PETSC_NULL_INTEGER,ierr)
347: ! Compute initial guess over the locally owned part of the grid
349: do 30 j=ys,ys+ym-1
350: temp = (min(j,my-j-1))*hy
351: do 40 i=xs,xs+xm-1
352: row = i - gxs + (j - gys)*gxm + 1
353: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. &
354: & j .eq. my-1) then
355: xx(idx+row) = 0.d0
356: continue
357: endif
358: xx(idx+row) = temp1*sqrt(min((min(i,mx-i-1))*hx,temp))
359: 40 continue
360: 30 continue
362: ! Restore vector
364: call VecRestoreArray(localX,xx,idx,ierr)
366: ! Insert values into global vector
368: call DALocalToGlobal(da,localX,INSERT_VALUES,X,ierr)
369: return
370: end
372: ! -------------------------------------------------------------------
373: !
374: ! ComputeFunction - Evaluates nonlinear function, F(x).
375: !
376: ! Input Parameters:
377: !. X - input vector
378: !
379: ! Output Parameter:
380: !. F - function vector
381: !
382: subroutine ComputeFunction(X,F,ierr)
383: implicit none
385: ! petsc.h - base PETSc routines petscvec.h - vectors
386: ! petscsys.h - system routines petscmat.h - matrices
387: ! petscis.h - index sets petscksp.h - Krylov subspace methods
388: ! petscviewer.h - viewers petscpc.h - preconditioners
390: #include finclude/petsc.h
391: #include finclude/petscis.h
392: #include finclude/petscvec.h
393: #include finclude/petscmat.h
394: #include finclude/petscpc.h
395: #include finclude/petscksp.h
396: #include finclude/petscda.h
398: Vec X,F,localX,localF
399: PetscInt gys,gxm,gym
400: PetscOffset idx,idf
401: PetscErrorCode ierr
402: PetscInt i,j,row,mx,my,xs,ys,xm,ym,gxs
403: double precision two,one,lambda,hx
404: double precision hy,hxdhy,hydhx,sc
405: PetscScalar u,uxx,uyy,xx(1),ff(1)
406: DA da
407: Mat B
408: common /mycommon/ mx,my,B,localX,localF,da
410: two = 2.d0
411: one = 1.d0
412: lambda = 6.d0
414: hx = one/(mx-1)
415: hy = one/(my-1)
416: sc = hx*hy*lambda
417: hxdhy = hx/hy
418: hydhx = hy/hx
420: ! Scatter ghost points to local vector, using the 2-step process
421: ! DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
422: ! By placing code between these two statements, computations can be
423: ! done while messages are in transition.
424: !
425: call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
426: call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
428: ! Get pointers to vector data
430: call VecGetArray(localX,xx,idx,ierr)
431: call VecGetArray(localF,ff,idf,ierr)
433: ! Get local grid boundaries
435: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
436: & PETSC_NULL_INTEGER,ierr)
437: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
438: & PETSC_NULL_INTEGER,ierr)
440: ! Compute function over the locally owned part of the grid
442: do 50 j=ys,ys+ym-1
444: row = (j - gys)*gxm + xs - gxs
445: do 60 i=xs,xs+xm-1
446: row = row + 1
448: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. &
449: & j .eq. my-1) then
450: ff(idf+row) = xx(idx+row)
451: goto 60
452: endif
453: u = xx(idx+row)
454: uxx = (two*u - xx(idx+row-1) - xx(idx+row+1))*hydhx
455: uyy = (two*u - xx(idx+row-gxm) - xx(idx+row+gxm))*hxdhy
456: ff(idf+row) = uxx + uyy - sc*exp(u)
457: 60 continue
458: 50 continue
460: ! Restore vectors
462: call VecRestoreArray(localX,xx,idx,ierr)
463: call VecRestoreArray(localF,ff,idf,ierr)
465: ! Insert values into global vector
467: call DALocalToGlobal(da,localF,INSERT_VALUES,F,ierr)
468: return
469: end
471: ! -------------------------------------------------------------------
472: !
473: ! ComputeJacobian - Evaluates Jacobian matrix.
474: !
475: ! Input Parameters:
476: ! x - input vector
477: !
478: ! Output Parameters:
479: ! jac - Jacobian matrix
480: ! flag - flag indicating matrix structure
481: !
482: ! Notes:
483: ! Due to grid point reordering with DAs, we must always work
484: ! with the local grid points, and then transform them to the new
485: ! global numbering with the 'ltog' mapping (via DAGetGlobalIndices()).
486: ! We cannot work directly with the global numbers for the original
487: ! uniprocessor grid!
488: !
489: subroutine ComputeJacobian(X,jac,ierr)
490: implicit none
492: ! petsc.h - base PETSc routines petscvec.h - vectors
493: ! petscsys.h - system routines petscmat.h - matrices
494: ! petscis.h - index sets petscksp.h - Krylov subspace methods
495: ! petscviewer.h - viewers petscpc.h - preconditioners
497: #include finclude/petsc.h
498: #include finclude/petscis.h
499: #include finclude/petscvec.h
500: #include finclude/petscmat.h
501: #include finclude/petscpc.h
502: #include finclude/petscksp.h
503: #include finclude/petscda.h
505: Vec X
506: Mat jac
507: Vec localX,localF
508: DA da
509: PetscInt ltog(1)
510: PetscOffset idltog,idx
511: PetscErrorCode ierr
512: PetscInt nloc,xs,ys,xm,ym
513: PetscInt gxs,gys,gxm,gym
514: PetscInt grow(1),i,j
515: PetscInt row,mx,my,ione
516: PetscInt col(5),ifive
517: PetscScalar two,one,lambda
518: PetscScalar v(5),hx,hy,hxdhy
519: PetscScalar hydhx,sc,xx(1)
520: Mat B
521: common /mycommon/ mx,my,B,localX,localF,da
523: ione = 1
524: ifive = 5
525: one = 1.d0
526: two = 2.d0
527: hx = one/(mx-1)
528: hy = one/(my-1)
529: sc = hx*hy
530: hxdhy = hx/hy
531: hydhx = hy/hx
532: lambda = 6.d0
534: ! Scatter ghost points to local vector, using the 2-step process
535: ! DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
536: ! By placing code between these two statements, computations can be
537: ! done while messages are in transition.
539: call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
540: call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
542: ! Get pointer to vector data
544: call VecGetArray(localX,xx,idx,ierr)
546: ! Get local grid boundaries
548: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
549: & PETSC_NULL_INTEGER,ierr)
550: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
551: & PETSC_NULL_INTEGER,ierr)
553: ! Get the global node numbers for all local nodes, including ghost points
555: call DAGetGlobalIndices(da,nloc,ltog,idltog,ierr)
557: ! Compute entries for the locally owned part of the Jacobian.
558: ! - Currently, all PETSc parallel matrix formats are partitioned by
559: ! contiguous chunks of rows across the processors. The 'grow'
560: ! parameter computed below specifies the global row number
561: ! corresponding to each local grid point.
562: ! - Each processor needs to insert only elements that it owns
563: ! locally (but any non-local elements will be sent to the
564: ! appropriate processor during matrix assembly).
565: ! - Always specify global row and columns of matrix entries.
566: ! - Here, we set all entries for a particular row at once.
568: do 10 j=ys,ys+ym-1
569: row = (j - gys)*gxm + xs - gxs
570: do 20 i=xs,xs+xm-1
571: row = row + 1
572: grow(1) = ltog(idltog+row)
573: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. (mx-1) .or. &
574: & j .eq. (my-1)) then
575: call MatSetValues(jac,ione,grow,ione,grow,one, &
576: & INSERT_VALUES,ierr)
577: go to 20
578: endif
579: v(1) = -hxdhy
580: col(1) = ltog(idltog+row - gxm)
581: v(2) = -hydhx
582: col(2) = ltog(idltog+row - 1)
583: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(xx(idx+row))
584: col(3) = grow(1)
585: v(4) = -hydhx
586: col(4) = ltog(idltog+row + 1)
587: v(5) = -hxdhy
588: col(5) = ltog(idltog+row + gxm)
589: call MatSetValues(jac,ione,grow,ifive,col,v,INSERT_VALUES, &
590: & ierr)
591: 20 continue
592: 10 continue
594: ! Assemble matrix, using the 2-step process:
595: ! MatAssemblyBegin(), MatAssemblyEnd().
596: ! By placing code between these two statements, computations can be
597: ! done while messages are in transition.
599: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
600: call VecRestoreArray(localX,xx,idx,ierr)
601: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
602: return
603: end
606: ! -------------------------------------------------------------------
607: !
608: ! MyMult - user provided matrix multiply
609: !
610: ! Input Parameters:
611: !. X - input vector
612: !
613: ! Output Parameter:
614: !. F - function vector
615: !
616: subroutine MyMult(J,X,F,ierr)
617: implicit none
618: Mat J,B
619: Vec X,F
620: PetscErrorCode ierr
621: PetscInt mx,my
622: DA da
623: Vec localX,localF
625: common /mycommon/ mx,my,B,localX,localF,da
626: !
627: ! Here we use the actual formed matrix B; users would
628: ! instead write their own matrix vector product routine
629: !
630: call MatMult(B,X,F,ierr)
631: return
632: end