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