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