Actual source code: ex14f.F90
1: !
2: !
3: ! Solves a nonlinear system in parallel with a user-defined
4: ! Newton method that uses KSP to solve the linearized Newton systems. This solver
5: ! is a very simplistic inexact Newton method. The intent of this code is to
6: ! demonstrate the repeated solution of linear systems 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 (DMDAs) 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: DMDA^using distributed arrays
27: ! Processors: n
28: !T*/
30: ! ------------------------------------------------------------------------
31: !
32: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
33: ! the partial differential equation
34: !
35: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
36: !
37: ! with boundary conditions
38: !
39: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
40: !
41: ! A finite difference approximation with the usual 5-point stencil
42: ! is used to discretize the boundary value problem to obtain a nonlinear
43: ! system of equations.
44: !
45: ! The SNES version of this problem is: snes/tutorials/ex5f.F
46: !
47: ! -------------------------------------------------------------------------
48: module mymoduleex14f
49: #include <petsc/finclude/petscksp.h>
50: use petscdmda
51: use petscksp
52: Vec localX
53: PetscInt mx,my
54: Mat B
55: DM da
56: end module
58: program main
59: use mymoduleex14f
60: implicit none
62: MPI_Comm comm
63: Vec X,Y,F
64: Mat J
65: KSP ksp
67: PetscInt Nx,Ny,N,ifive,ithree
68: PetscBool flg,nooutput,usemf
69: !
70: ! This is the routine to use for matrix-free approach
71: !
72: external mymult
74: ! --------------- Data to define nonlinear solver --------------
75: PetscReal rtol,ttol
76: PetscReal fnorm,ynorm,xnorm
77: PetscInt max_nonlin_its,one
78: PetscInt lin_its
79: PetscInt i,m
80: PetscScalar mone
81: PetscErrorCode ierr
83: mone = -1.0
84: rtol = 1.e-8
85: max_nonlin_its = 10
86: one = 1
87: ifive = 5
88: ithree = 3
90: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
91: if (ierr .ne. 0) then
92: print*,'Unable to initialize PETSc'
93: stop
94: endif
95: comm = PETSC_COMM_WORLD
97: ! Initialize problem parameters
99: !
100: mx = 4
101: my = 4
102: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
103: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
104: N = mx*my
106: nooutput = .false.
107: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_output',nooutput,ierr)
109: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: ! Create linear solver context
111: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: call KSPCreate(comm,ksp,ierr)
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: ! Create vector data structures
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: !
120: ! Create distributed array (DMDA) to manage parallel grid and vectors
121: !
122: Nx = PETSC_DECIDE
123: Ny = PETSC_DECIDE
124: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-Nx',Nx,flg,ierr)
125: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-Ny',Ny,flg,ierr)
126: call DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,mx,my,Nx,Ny,one,one, &
127: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
128: call DMSetFromOptions(da,ierr)
129: call DMSetUp(da,ierr)
130: !
131: ! Extract global and local vectors from DMDA then duplicate for remaining
132: ! vectors that are the same types
133: !
134: call DMCreateGlobalVector(da,X,ierr)
135: call DMCreateLocalVector(da,localX,ierr)
136: call VecDuplicate(X,F,ierr)
137: call VecDuplicate(X,Y,ierr)
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: ! Create matrix data structure for Jacobian
141: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: !
143: ! Note: For the parallel case, vectors and matrices MUST be partitioned
144: ! accordingly. When using distributed arrays (DMDAs) to create vectors,
145: ! the DMDAs determine the problem partitioning. We must explicitly
146: ! specify the local matrix dimensions upon its creation for compatibility
147: ! with the vector distribution.
148: !
149: ! Note: Here we only approximately preallocate storage space for the
150: ! Jacobian. See the users manual for a discussion of better techniques
151: ! for preallocating matrix memory.
152: !
153: call VecGetLocalSize(X,m,ierr)
154: call MatCreateAIJ(comm,m,m,N,N,ifive,PETSC_NULL_INTEGER,ithree,PETSC_NULL_INTEGER,B,ierr)
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! if usemf is on then matrix vector product is done via matrix free
158: ! approach. Note this is just an example, and not realistic because
159: ! we still use the actual formed matrix, but in reality one would
160: ! provide their own subroutine that would directly do the matrix
161: ! vector product and not call MatMult()
162: ! Note: we put B into a module so it will be visible to the
163: ! mymult() routine
164: usemf = .false.
165: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mf',usemf,ierr)
166: if (usemf) then
167: call MatCreateShell(comm,m,m,N,N,PETSC_NULL_INTEGER,J,ierr)
168: call MatShellSetOperation(J,MATOP_MULT,mymult,ierr)
169: else
170: ! If not doing matrix free then matrix operator, J, and matrix used
171: ! to construct preconditioner, B, are the same
172: J = B
173: endif
175: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: ! Customize linear solver set runtime options
177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: !
179: ! Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
180: !
181: call KSPSetFromOptions(ksp,ierr)
183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: ! Evaluate initial guess
185: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: call FormInitialGuess(X,ierr)
188: call ComputeFunction(X,F,ierr)
189: call VecNorm(F,NORM_2,fnorm,ierr)
190: ttol = fnorm*rtol
191: if (.not. nooutput) then
192: print*, 'Initial function norm ',fnorm
193: endif
195: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: ! Solve nonlinear system with a user-defined method
197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: ! This solver is a very simplistic inexact Newton method, with no
200: ! no damping strategies or bells and whistles. The intent of this code
201: ! is merely to demonstrate the repeated solution with KSP of linear
202: ! systems with the same nonzero structure.
203: !
204: ! This is NOT the recommended approach for solving nonlinear problems
205: ! with PETSc! We urge users to employ the SNES component for solving
206: ! nonlinear problems whenever possible with application codes, as it
207: ! offers many advantages over coding nonlinear solvers independently.
209: do 10 i=0,max_nonlin_its
211: ! Compute the Jacobian matrix. See the comments in this routine for
212: ! important information about setting the flag mat_flag.
214: call ComputeJacobian(X,B,ierr)
216: ! Solve J Y = F, where J is the Jacobian matrix.
217: ! - First, set the KSP linear operators. Here the matrix that
218: ! defines the linear system also serves as the preconditioning
219: ! matrix.
220: ! - Then solve the Newton system.
222: call KSPSetOperators(ksp,J,B,ierr)
223: call KSPSolve(ksp,F,Y,ierr)
225: ! Compute updated iterate
227: call VecNorm(Y,NORM_2,ynorm,ierr)
228: call VecAYPX(Y,mone,X,ierr)
229: call VecCopy(Y,X,ierr)
230: call VecNorm(X,NORM_2,xnorm,ierr)
231: call KSPGetIterationNumber(ksp,lin_its,ierr)
232: if (.not. nooutput) then
233: print*,'linear solve iterations = ',lin_its,' xnorm = ',xnorm,' ynorm = ',ynorm
234: endif
236: ! Evaluate nonlinear function at new location
238: call ComputeFunction(X,F,ierr)
239: call VecNorm(F,NORM_2,fnorm,ierr)
240: if (.not. nooutput) then
241: print*, 'Iteration ',i+1,' function norm',fnorm
242: endif
244: ! Test for convergence
246: if (fnorm .le. ttol) then
247: if (.not. nooutput) then
248: print*,'Converged: function norm ',fnorm,' tolerance ',ttol
249: endif
250: goto 20
251: endif
252: 10 continue
253: 20 continue
255: write(6,100) i+1
256: 100 format('Number of SNES iterations =',I2)
258: ! Check if mymult() produces a linear operator
259: if (usemf) then
260: N = 5
261: call MatIsLinear(J,N,flg,ierr)
262: if (.not. flg) then
263: print *, 'IsLinear',flg
264: endif
265: endif
267: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268: ! Free work space. All PETSc objects should be destroyed when they
269: ! are no longer needed.
270: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272: call MatDestroy(B,ierr)
273: if (usemf) then
274: call MatDestroy(J,ierr)
275: endif
276: call VecDestroy(localX,ierr)
277: call VecDestroy(X,ierr)
278: call VecDestroy(Y,ierr)
279: call VecDestroy(F,ierr)
280: call KSPDestroy(ksp,ierr)
281: call DMDestroy(da,ierr)
282: call PetscFinalize(ierr)
283: end
285: ! -------------------------------------------------------------------
286: !
287: ! FormInitialGuess - Forms initial approximation.
288: !
289: ! Input Parameters:
290: ! X - vector
291: !
292: ! Output Parameter:
293: ! X - vector
294: !
295: subroutine FormInitialGuess(X,ierr)
296: use mymoduleex14f
297: implicit none
299: PetscErrorCode ierr
300: PetscOffset idx
301: Vec X
302: PetscInt i,j,row
303: PetscInt xs,ys,xm
304: PetscInt ym
305: PetscReal one,lambda,temp1,temp,hx,hy
306: PetscScalar xx(2)
308: one = 1.0
309: lambda = 6.0
310: hx = one/(mx-1)
311: hy = one/(my-1)
312: temp1 = lambda/(lambda + one)
314: ! Get a pointer to vector data.
315: ! - VecGetArray() returns a pointer to the data array.
316: ! - You MUST call VecRestoreArray() when you no longer need access to
317: ! the array.
318: call VecGetArray(X,xx,idx,ierr)
320: ! Get local grid boundaries (for 2-dimensional DMDA):
321: ! xs, ys - starting grid indices (no ghost points)
322: ! xm, ym - widths of local grid (no ghost points)
324: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
326: ! Compute initial guess over the locally owned part of the grid
328: do 30 j=ys,ys+ym-1
329: temp = (min(j,my-j-1))*hy
330: do 40 i=xs,xs+xm-1
331: row = i - xs + (j - ys)*xm + 1
332: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. j .eq. my-1) then
333: xx(idx+row) = 0.0
334: continue
335: endif
336: xx(idx+row) = temp1*sqrt(min((min(i,mx-i-1))*hx,temp))
337: 40 continue
338: 30 continue
340: ! Restore vector
342: call VecRestoreArray(X,xx,idx,ierr)
343: return
344: end
346: ! -------------------------------------------------------------------
347: !
348: ! ComputeFunction - Evaluates nonlinear function, F(x).
349: !
350: ! Input Parameters:
351: !. X - input vector
352: !
353: ! Output Parameter:
354: !. F - function vector
355: !
356: subroutine ComputeFunction(X,F,ierr)
357: use mymoduleex14f
358: implicit none
360: Vec X,F
361: PetscInt gys,gxm,gym
362: PetscOffset idx,idf
363: PetscErrorCode ierr
364: PetscInt i,j,row,xs,ys,xm,ym,gxs
365: PetscInt rowf
366: PetscReal two,one,lambda,hx
367: PetscReal hy,hxdhy,hydhx,sc
368: PetscScalar u,uxx,uyy,xx(2),ff(2)
370: two = 2.0
371: one = 1.0
372: lambda = 6.0
374: hx = one/(mx-1)
375: hy = one/(my-1)
376: sc = hx*hy*lambda
377: hxdhy = hx/hy
378: hydhx = hy/hx
380: ! Scatter ghost points to local vector, using the 2-step process
381: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
382: ! By placing code between these two statements, computations can be
383: ! done while messages are in transition.
384: !
385: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
386: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
388: ! Get pointers to vector data
390: call VecGetArray(localX,xx,idx,ierr)
391: call VecGetArray(F,ff,idf,ierr)
393: ! Get local grid boundaries
395: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
396: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)
398: ! Compute function over the locally owned part of the grid
399: rowf = 0
400: do 50 j=ys,ys+ym-1
402: row = (j - gys)*gxm + xs - gxs
403: do 60 i=xs,xs+xm-1
404: row = row + 1
405: rowf = rowf + 1
407: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. j .eq. my-1) then
408: ff(idf+rowf) = xx(idx+row)
409: goto 60
410: endif
411: u = xx(idx+row)
412: uxx = (two*u - xx(idx+row-1) - xx(idx+row+1))*hydhx
413: uyy = (two*u - xx(idx+row-gxm) - xx(idx+row+gxm))*hxdhy
414: ff(idf+rowf) = uxx + uyy - sc*exp(u)
415: 60 continue
416: 50 continue
418: ! Restore vectors
420: call VecRestoreArray(localX,xx,idx,ierr)
421: call VecRestoreArray(F,ff,idf,ierr)
422: return
423: end
425: ! -------------------------------------------------------------------
426: !
427: ! ComputeJacobian - Evaluates Jacobian matrix.
428: !
429: ! Input Parameters:
430: ! x - input vector
431: !
432: ! Output Parameters:
433: ! jac - Jacobian matrix
434: ! flag - flag indicating matrix structure
435: !
436: ! Notes:
437: ! Due to grid point reordering with DMDAs, we must always work
438: ! with the local grid points, and then transform them to the new
439: ! global numbering with the 'ltog' mapping
440: ! We cannot work directly with the global numbers for the original
441: ! uniprocessor grid!
442: !
443: subroutine ComputeJacobian(X,jac,ierr)
444: use mymoduleex14f
445: implicit none
447: Vec X
448: Mat jac
449: PetscInt ltog(2)
450: PetscOffset idltog,idx
451: PetscErrorCode ierr
452: PetscInt xs,ys,xm,ym
453: PetscInt gxs,gys,gxm,gym
454: PetscInt grow(1),i,j
455: PetscInt row,ione
456: PetscInt col(5),ifive
457: PetscScalar two,one,lambda
458: PetscScalar v(5),hx,hy,hxdhy
459: PetscScalar hydhx,sc,xx(2)
460: ISLocalToGlobalMapping ltogm
462: ione = 1
463: ifive = 5
464: one = 1.0
465: two = 2.0
466: hx = one/(mx-1)
467: hy = one/(my-1)
468: sc = hx*hy
469: hxdhy = hx/hy
470: hydhx = hy/hx
471: lambda = 6.0
473: ! Scatter ghost points to local vector, using the 2-step process
474: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
475: ! By placing code between these two statements, computations can be
476: ! done while messages are in transition.
478: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
479: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
481: ! Get pointer to vector data
483: call VecGetArray(localX,xx,idx,ierr)
485: ! Get local grid boundaries
487: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
488: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)
490: ! Get the global node numbers for all local nodes, including ghost points
492: call DMGetLocalToGlobalMapping(da,ltogm,ierr)
493: call ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)
495: ! Compute entries for the locally owned part of the Jacobian.
496: ! - Currently, all PETSc parallel matrix formats are partitioned by
497: ! contiguous chunks of rows across the processors. The 'grow'
498: ! parameter computed below specifies the global row number
499: ! corresponding to each local grid point.
500: ! - Each processor needs to insert only elements that it owns
501: ! locally (but any non-local elements will be sent to the
502: ! appropriate processor during matrix assembly).
503: ! - Always specify global row and columns of matrix entries.
504: ! - Here, we set all entries for a particular row at once.
506: do 10 j=ys,ys+ym-1
507: row = (j - gys)*gxm + xs - gxs
508: do 20 i=xs,xs+xm-1
509: row = row + 1
510: grow(1) = ltog(idltog+row)
511: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. (mx-1) .or. j .eq. (my-1)) then
512: call MatSetValues(jac,ione,grow,ione,grow,one,INSERT_VALUES,ierr)
513: go to 20
514: endif
515: v(1) = -hxdhy
516: col(1) = ltog(idltog+row - gxm)
517: v(2) = -hydhx
518: col(2) = ltog(idltog+row - 1)
519: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(xx(idx+row))
520: col(3) = grow(1)
521: v(4) = -hydhx
522: col(4) = ltog(idltog+row + 1)
523: v(5) = -hxdhy
524: col(5) = ltog(idltog+row + gxm)
525: call MatSetValues(jac,ione,grow,ifive,col,v,INSERT_VALUES,ierr)
526: 20 continue
527: 10 continue
529: call ISLocalToGlobalMappingRestoreIndices(ltogm,ltog,idltog,ierr)
531: ! Assemble matrix, using the 2-step process:
532: ! MatAssemblyBegin(), MatAssemblyEnd().
533: ! By placing code between these two statements, computations can be
534: ! done while messages are in transition.
536: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
537: call VecRestoreArray(localX,xx,idx,ierr)
538: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
539: return
540: end
542: ! -------------------------------------------------------------------
543: !
544: ! MyMult - user provided matrix multiply
545: !
546: ! Input Parameters:
547: !. X - input vector
548: !
549: ! Output Parameter:
550: !. F - function vector
551: !
552: subroutine MyMult(J,X,F,ierr)
553: use mymoduleex14f
554: implicit none
556: Mat J
557: Vec X,F
558: PetscErrorCode ierr
559: !
560: ! Here we use the actual formed matrix B; users would
561: ! instead write their own matrix vector product routine
562: !
563: call MatMult(B,X,F,ierr)
564: return
565: end
567: !/*TEST
568: !
569: ! test:
570: ! args: -no_output -ksp_gmres_cgs_refinement_type refine_always
571: ! output_file: output/ex14_1.out
572: ! requires: !single
573: !
574: !TEST*/