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