Actual source code: ex5f.F90
1: !
2: ! Description: This example solves a nonlinear system in parallel with SNES.
3: ! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
4: ! domain, using distributed arrays (DMDAs) to partition the parallel grid.
5: ! The command line options include:
6: ! -par <param>, where <param> indicates the nonlinearity of the problem
7: ! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)
8: !
9: !
11: !
12: ! --------------------------------------------------------------------------
13: !
14: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
15: ! the partial differential equation
16: !
17: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
18: !
19: ! with boundary conditions
20: !
21: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
22: !
23: ! A finite difference approximation with the usual 5-point stencil
24: ! is used to discretize the boundary value problem to obtain a nonlinear
25: ! system of equations.
26: !
27: ! --------------------------------------------------------------------------
29: program main
30: #include <petsc/finclude/petscsnes.h>
31: use petscdmda
32: use petscsnes
33: implicit none
34: !
35: ! We place common blocks, variable declarations, and other include files
36: ! needed for this code in the single file ex5f.h. We then need to include
37: ! only this file throughout the various routines in this program. See
38: ! additional comments in the file ex5f.h.
39: !
40: #include "ex5f.h"
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Variable declarations
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: !
46: ! Variables:
47: ! snes - nonlinear solver
48: ! x, r - solution, residual vectors
49: ! its - iterations for convergence
50: !
51: ! See additional variable declarations in the file ex5f.h
52: !
53: SNES snes
54: Vec x,r
55: PetscInt its,i1,i4
56: PetscErrorCode ierr
57: PetscReal lambda_max,lambda_min
58: PetscBool flg
59: DM da
61: ! Note: Any user-defined Fortran routines (such as FormJacobianLocal)
62: ! MUST be declared as external.
64: external FormInitialGuess
65: external FormFunctionLocal,FormJacobianLocal
66: external MySNESConverged
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: ! Initialize program
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
73: if (ierr .ne. 0) then
74: print*,'Unable to initialize PETSc'
75: stop
76: endif
77: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
78: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
80: ! Initialize problem parameters
82: i1 = 1
83: i4 = 4
84: lambda_max = 6.81
85: lambda_min = 0.0
86: lambda = 6.0
87: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',lambda,PETSC_NULL_BOOL,ierr)
88: ! this statement is split into multiple-lines to keep lines under 132 char limit - required by 'make check'
89: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then
90: PETSC_ERR_ARG_OUTOFRANGE; SETERRA(PETSC_COMM_WORLD,ierr,'Lambda')
91: endif
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: ! Create nonlinear solver context
95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
99: ! Set convergence test routine if desired
101: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my_snes_convergence',flg,ierr)
102: if (flg) then
103: call SNESSetConvergenceTest(snes,MySNESConverged,0,PETSC_NULL_FUNCTION,ierr)
104: endif
106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: ! Create vector data structures; set function evaluation routine
108: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: ! Create distributed array (DMDA) to manage parallel grid and vectors
112: ! This really needs only the star-type stencil, but we use the box
113: ! stencil temporarily.
114: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE, &
115: & DMDA_STENCIL_STAR,i4,i4,PETSC_DECIDE,PETSC_DECIDE,i1,i1, &
116: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
117: call DMSetFromOptions(da,ierr)
118: call DMSetUp(da,ierr)
120: ! Extract global and local vectors from DMDA; then duplicate for remaining
121: ! vectors that are the same types
123: call DMCreateGlobalVector(da,x,ierr)
124: call VecDuplicate(x,r,ierr)
126: ! Get local grid boundaries (for 2-dimensional DMDA)
128: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER, &
129: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
130: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
131: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
132: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
133: & PETSC_NULL_INTEGER,ierr)
134: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
135: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)
137: ! Here we shift the starting indices up by one so that we can easily
138: ! use the Fortran convention of 1-based indices (rather 0-based indices).
140: xs = xs+1
141: ys = ys+1
142: gxs = gxs+1
143: gys = gys+1
145: ye = ys+ym-1
146: xe = xs+xm-1
147: gye = gys+gym-1
148: gxe = gxs+gxm-1
150: ! Set function evaluation routine and vector
152: call DMDASNESSetFunctionLocal(da,INSERT_VALUES,FormFunctionLocal,da,ierr)
153: call DMDASNESSetJacobianLocal(da,FormJacobianLocal,da,ierr)
154: call SNESSetDM(snes,da,ierr)
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! Customize nonlinear solver; set runtime options
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
162: call SNESSetFromOptions(snes,ierr)
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: ! Evaluate initial guess; then solve nonlinear system.
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: ! Note: The user should initialize the vector, x, with the initial guess
168: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
169: ! to employ an initial guess of zero, the user should explicitly set
170: ! this vector to zero by calling VecSet().
172: call FormInitialGuess(x,ierr)
173: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
174: call SNESGetIterationNumber(snes,its,ierr)
175: if (rank .eq. 0) then
176: write(6,100) its
177: endif
178: 100 format('Number of SNES iterations = ',i5)
180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: ! Free work space. All PETSc objects should be destroyed when they
182: ! are no longer needed.
183: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: call VecDestroy(x,ierr)
186: call VecDestroy(r,ierr)
187: call SNESDestroy(snes,ierr)
188: call DMDestroy(da,ierr)
189: call PetscFinalize(ierr)
190: end
192: ! ---------------------------------------------------------------------
193: !
194: ! FormInitialGuess - Forms initial approximation.
195: !
196: ! Input Parameters:
197: ! X - vector
198: !
199: ! Output Parameter:
200: ! X - vector
201: !
202: ! Notes:
203: ! This routine serves as a wrapper for the lower-level routine
204: ! "ApplicationInitialGuess", where the actual computations are
205: ! done using the standard Fortran style of treating the local
206: ! vector data as a multidimensional array over the local mesh.
207: ! This routine merely handles ghost point scatters and accesses
208: ! the local vector data via VecGetArray() and VecRestoreArray().
209: !
210: subroutine FormInitialGuess(X,ierr)
211: use petscsnes
212: implicit none
214: #include "ex5f.h"
216: ! Input/output variables:
217: Vec X
218: PetscErrorCode ierr
220: ! Declarations for use with local arrays:
221: PetscScalar lx_v(0:1)
222: PetscOffset lx_i
224: 0
226: ! Get a pointer to vector data.
227: ! - For default PETSc vectors, VecGetArray() returns a pointer to
228: ! the data array. Otherwise, the routine is implementation dependent.
229: ! - You MUST call VecRestoreArray() when you no longer need access to
230: ! the array.
231: ! - Note that the Fortran interface to VecGetArray() differs from the
232: ! C version. See the users manual for details.
234: call VecGetArray(X,lx_v,lx_i,ierr)
236: ! Compute initial guess over the locally owned part of the grid
238: call InitialGuessLocal(lx_v(lx_i),ierr)
240: ! Restore vector
242: call VecRestoreArray(X,lx_v,lx_i,ierr)
244: return
245: end
247: ! ---------------------------------------------------------------------
248: !
249: ! InitialGuessLocal - Computes initial approximation, called by
250: ! the higher level routine FormInitialGuess().
251: !
252: ! Input Parameter:
253: ! x - local vector data
254: !
255: ! Output Parameters:
256: ! x - local vector data
257: ! ierr - error code
258: !
259: ! Notes:
260: ! This routine uses standard Fortran-style computations over a 2-dim array.
261: !
262: subroutine InitialGuessLocal(x,ierr)
263: use petscsnes
264: implicit none
266: #include "ex5f.h"
268: ! Input/output variables:
269: PetscScalar x(xs:xe,ys:ye)
270: PetscErrorCode ierr
272: ! Local variables:
273: PetscInt i,j
274: PetscReal temp1,temp,one,hx,hy
276: ! Set parameters
278: 0
279: one = 1.0
280: hx = one/((mx-1))
281: hy = one/((my-1))
282: temp1 = lambda/(lambda + one)
284: do 20 j=ys,ye
285: temp = (min(j-1,my-j))*hy
286: do 10 i=xs,xe
287: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
288: x(i,j) = 0.0
289: else
290: x(i,j) = temp1 * sqrt(min(min(i-1,mx-i)*hx,(temp)))
291: endif
292: 10 continue
293: 20 continue
295: return
296: end
298: ! ---------------------------------------------------------------------
299: !
300: ! FormFunctionLocal - Computes nonlinear function, called by
301: ! the higher level routine FormFunction().
302: !
303: ! Input Parameter:
304: ! x - local vector data
305: !
306: ! Output Parameters:
307: ! f - local vector data, f(x)
308: ! ierr - error code
309: !
310: ! Notes:
311: ! This routine uses standard Fortran-style computations over a 2-dim array.
312: !
313: !
314: subroutine FormFunctionLocal(info,x,f,da,ierr)
315: #include <petsc/finclude/petscdmda.h>
316: use petscsnes
317: implicit none
319: #include "ex5f.h"
320: DM da
322: ! Input/output variables:
323: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
324: PetscScalar x(gxs:gxe,gys:gye)
325: PetscScalar f(xs:xe,ys:ye)
326: PetscErrorCode ierr
328: ! Local variables:
329: PetscScalar two,one,hx,hy
330: PetscScalar hxdhy,hydhx,sc
331: PetscScalar u,uxx,uyy
332: PetscInt i,j
334: xs = info(DMDA_LOCAL_INFO_XS)+1
335: xe = xs+info(DMDA_LOCAL_INFO_XM)-1
336: ys = info(DMDA_LOCAL_INFO_YS)+1
337: ye = ys+info(DMDA_LOCAL_INFO_YM)-1
338: mx = info(DMDA_LOCAL_INFO_MX)
339: my = info(DMDA_LOCAL_INFO_MY)
341: one = 1.0
342: two = 2.0
343: hx = one/(mx-1)
344: hy = one/(my-1)
345: sc = hx*hy*lambda
346: hxdhy = hx/hy
347: hydhx = hy/hx
349: ! Compute function over the locally owned part of the grid
351: do 20 j=ys,ye
352: do 10 i=xs,xe
353: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
354: f(i,j) = x(i,j)
355: else
356: u = x(i,j)
357: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
358: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
359: f(i,j) = uxx + uyy - sc*exp(u)
360: endif
361: 10 continue
362: 20 continue
364: call PetscLogFlops(11.0d0*ym*xm,ierr)
366: return
367: end
369: ! ---------------------------------------------------------------------
370: !
371: ! FormJacobianLocal - Computes Jacobian matrix, called by
372: ! the higher level routine FormJacobian().
373: !
374: ! Input Parameters:
375: ! x - local vector data
376: !
377: ! Output Parameters:
378: ! jac - Jacobian matrix
379: ! jac_prec - optionally different preconditioning matrix (not used here)
380: ! ierr - error code
381: !
382: ! Notes:
383: ! This routine uses standard Fortran-style computations over a 2-dim array.
384: !
385: ! Notes:
386: ! Due to grid point reordering with DMDAs, we must always work
387: ! with the local grid points, and then transform them to the new
388: ! global numbering with the "ltog" mapping
389: ! We cannot work directly with the global numbers for the original
390: ! uniprocessor grid!
391: !
392: ! Two methods are available for imposing this transformation
393: ! when setting matrix entries:
394: ! (A) MatSetValuesLocal(), using the local ordering (including
395: ! ghost points!)
396: ! by calling MatSetValuesLocal()
397: ! (B) MatSetValues(), using the global ordering
398: ! - Use DMDAGetGlobalIndices() to extract the local-to-global map
399: ! - Then apply this map explicitly yourself
400: ! - Set matrix entries using the global ordering by calling
401: ! MatSetValues()
402: ! Option (A) seems cleaner/easier in many cases, and is the procedure
403: ! used in this example.
404: !
405: subroutine FormJacobianLocal(info,x,A,jac,da,ierr)
406: use petscsnes
407: implicit none
409: #include "ex5f.h"
410: DM da
412: ! Input/output variables:
413: PetscScalar x(gxs:gxe,gys:gye)
414: Mat A,jac
415: PetscErrorCode ierr
416: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
418: ! Local variables:
419: PetscInt row,col(5),i,j,i1,i5
420: PetscScalar two,one,hx,hy,v(5)
421: PetscScalar hxdhy,hydhx,sc
423: ! Set parameters
425: i1 = 1
426: i5 = 5
427: one = 1.0
428: two = 2.0
429: hx = one/(mx-1)
430: hy = one/(my-1)
431: sc = hx*hy
432: hxdhy = hx/hy
433: hydhx = hy/hx
435: ! Compute entries for the locally owned part of the Jacobian.
436: ! - Currently, all PETSc parallel matrix formats are partitioned by
437: ! contiguous chunks of rows across the processors.
438: ! - Each processor needs to insert only elements that it owns
439: ! locally (but any non-local elements will be sent to the
440: ! appropriate processor during matrix assembly).
441: ! - Here, we set all entries for a particular row at once.
442: ! - We can set matrix entries either using either
443: ! MatSetValuesLocal() or MatSetValues(), as discussed above.
444: ! - Note that MatSetValues() uses 0-based row and column numbers
445: ! in Fortran as well as in C.
447: do 20 j=ys,ye
448: row = (j - gys)*gxm + xs - gxs - 1
449: do 10 i=xs,xe
450: row = row + 1
451: ! boundary points
452: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
453: ! Some f90 compilers need 4th arg to be of same type in both calls
454: col(1) = row
455: v(1) = one
456: call MatSetValuesLocal(jac,i1,row,i1,col,v,INSERT_VALUES,ierr)
457: ! interior grid points
458: else
459: v(1) = -hxdhy
460: v(2) = -hydhx
461: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
462: v(4) = -hydhx
463: v(5) = -hxdhy
464: col(1) = row - gxm
465: col(2) = row - 1
466: col(3) = row
467: col(4) = row + 1
468: col(5) = row + gxm
469: call MatSetValuesLocal(jac,i1,row,i5,col,v, INSERT_VALUES,ierr)
470: endif
471: 10 continue
472: 20 continue
473: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
474: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
475: if (A .ne. jac) then
476: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
477: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
478: endif
479: return
480: end
482: !
483: ! Simple convergence test based on the infinity norm of the residual being small
484: !
485: subroutine MySNESConverged(snes,it,xnorm,snorm,fnorm,reason,dummy,ierr)
486: use petscsnes
487: implicit none
489: SNES snes
490: PetscInt it,dummy
491: PetscReal xnorm,snorm,fnorm,nrm
492: SNESConvergedReason reason
493: Vec f
494: PetscErrorCode ierr
496: call SNESGetFunction(snes,f,PETSC_NULL_FUNCTION,dummy,ierr)
497: call VecNorm(f,NORM_INFINITY,nrm,ierr)
498: if (nrm .le. 1.e-5) reason = SNES_CONVERGED_FNORM_ABS
500: end
502: !/*TEST
503: !
504: ! build:
505: ! requires: !complex !single
506: !
507: ! test:
508: ! nsize: 4
509: ! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short \
510: ! -ksp_gmres_cgs_refinement_type refine_always
511: !
512: ! test:
513: ! suffix: 2
514: ! nsize: 4
515: ! args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
516: !
517: ! test:
518: ! suffix: 3
519: ! nsize: 3
520: ! args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
521: !
522: ! test:
523: ! suffix: 6
524: ! nsize: 1
525: ! args: -snes_monitor_short -my_snes_convergence
526: !
527: !TEST*/