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