Actual source code: ex5f.F
petsc-3.3-p7 2013-05-11
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: ! Program usage: mpiexec -n <procs> ex5f [-help] [all PETSc options]
10: !
11: !/*T
12: ! Concepts: SNES^parallel Bratu example
13: ! Concepts: DMDA^using distributed arrays;
14: ! Processors: n
15: !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: implicit none
36: !
37: ! We place common blocks, variable declarations, and other include files
38: ! needed for this code in the single file ex5f.h. We then need to include
39: ! only this file throughout the various routines in this program. See
40: ! additional comments in the file ex5f.h.
41: !
42: #include "ex5f.h"
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Variable declarations
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: !
48: ! Variables:
49: ! snes - nonlinear solver
50: ! x, r - solution, residual vectors
51: ! its - iterations for convergence
52: !
53: ! See additional variable declarations in the file ex5f.h
54: !
55: SNES snes
56: Vec x,r
57: PetscInt its,i1,i4
58: PetscErrorCode ierr
59: PetscReal lambda_max,lambda_min
60: PetscBool flg
63: ! Note: Any user-defined Fortran routines (such as FormJacobianLocal)
64: ! MUST be declared as external.
66: external FormInitialGuess
67: external FormFunctionLocal,FormJacobianLocal
69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: ! Initialize program
71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
74: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
75: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
77: ! Initialize problem parameters
79: i1 = 1
80: i4 = -4
81: lambda_max = 6.81
82: lambda_min = 0.0
83: lambda = 6.0
84: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-par',lambda, &
85: & flg,ierr)
86: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then
87: if (rank .eq. 0) write(6,*) 'Lambda is out of range'
88: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
89: endif
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: ! Create nonlinear solver context
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: ! Create vector data structures; set function evaluation routine
99: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: ! Create distributed array (DMDA) to manage parallel grid and vectors
103: ! This really needs only the star-type stencil, but we use the box
104: ! stencil temporarily.
105: call DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, &
106: & DMDA_BOUNDARY_NONE, &
107: & DMDA_STENCIL_STAR,i4,i4,PETSC_DECIDE,PETSC_DECIDE,i1,i1, &
108: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
110: ! Extract global and local vectors from DMDA; then duplicate for remaining
111: ! vectors that are the same types
113: call DMCreateGlobalVector(da,x,ierr)
114: call VecDuplicate(x,r,ierr)
116: ! Get local grid boundaries (for 2-dimensional DMDA)
118: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER, &
119: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
120: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
121: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
122: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
123: & PETSC_NULL_INTEGER,ierr)
124: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
125: & PETSC_NULL_INTEGER,ierr)
126: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
127: & PETSC_NULL_INTEGER,ierr)
129: ! Here we shift the starting indices up by one so that we can easily
130: ! use the Fortran convention of 1-based indices (rather 0-based indices).
132: xs = xs+1
133: ys = ys+1
134: gxs = gxs+1
135: gys = gys+1
137: ye = ys+ym-1
138: xe = xs+xm-1
139: gye = gys+gym-1
140: gxe = gxs+gxm-1
142: ! Set function evaluation routine and vector
144: call DMDASetLocalFunction(da,FormFunctionLocal,ierr)
145: call DMDASetLocalJacobian(da,FormJacobianLocal,ierr)
146: call SNESSetDM(snes,da,ierr)
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: ! Customize nonlinear solver; set runtime options
150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
154: call SNESSetFromOptions(snes,ierr)
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: ! Evaluate initial guess; then solve nonlinear system.
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Note: The user should initialize the vector, x, with the initial guess
160: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
161: ! to employ an initial guess of zero, the user should explicitly set
162: ! this vector to zero by calling VecSet().
164: call FormInitialGuess(x,ierr)
165: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
166: call SNESGetIterationNumber(snes,its,ierr)
167: if (rank .eq. 0) then
168: write(6,100) its
169: endif
170: 100 format('Number of SNES iterations = ',i5)
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: ! Free work space. All PETSc objects should be destroyed when they
175: ! are no longer needed.
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: call VecDestroy(x,ierr)
179: call VecDestroy(r,ierr)
180: call SNESDestroy(snes,ierr)
181: call DMDestroy(da,ierr)
182: call PetscFinalize(ierr)
183: end
185: ! ---------------------------------------------------------------------
186: !
187: ! FormInitialGuess - Forms initial approximation.
188: !
189: ! Input Parameters:
190: ! X - vector
191: !
192: ! Output Parameter:
193: ! X - vector
194: !
195: ! Notes:
196: ! This routine serves as a wrapper for the lower-level routine
197: ! "ApplicationInitialGuess", where the actual computations are
198: ! done using the standard Fortran style of treating the local
199: ! vector data as a multidimensional array over the local mesh.
200: ! This routine merely handles ghost point scatters and accesses
201: ! the local vector data via VecGetArray() and VecRestoreArray().
202: !
203: subroutine FormInitialGuess(X,ierr)
204: implicit none
206: #include "ex5f.h"
208: ! Input/output variables:
209: Vec X
210: PetscErrorCode ierr
212: ! Declarations for use with local arrays:
213: PetscScalar lx_v(0:1)
214: PetscOffset lx_i
215: Vec localX
217: 0
219: ! Get a pointer to vector data.
220: ! - For default PETSc vectors, VecGetArray() returns a pointer to
221: ! the data array. Otherwise, the routine is implementation dependent.
222: ! - You MUST call VecRestoreArray() when you no longer need access to
223: ! the array.
224: ! - Note that the Fortran interface to VecGetArray() differs from the
225: ! C version. See the users manual for details.
227: call DMGetLocalVector(da,localX,ierr)
228: call VecGetArray(localX,lx_v,lx_i,ierr)
230: ! Compute initial guess over the locally owned part of the grid
232: call InitialGuessLocal(lx_v(lx_i),ierr)
234: ! Restore vector
236: call VecRestoreArray(localX,lx_v,lx_i,ierr)
238: ! Insert values into global vector
240: call DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X,ierr)
241: call DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X,ierr)
242: call DMRestoreLocalVector(da,localX,ierr)
243: return
244: end
246: ! ---------------------------------------------------------------------
247: !
248: ! InitialGuessLocal - Computes initial approximation, called by
249: ! the higher level routine FormInitialGuess().
250: !
251: ! Input Parameter:
252: ! x - local vector data
253: !
254: ! Output Parameters:
255: ! x - local vector data
256: ! ierr - error code
257: !
258: ! Notes:
259: ! This routine uses standard Fortran-style computations over a 2-dim array.
260: !
261: subroutine InitialGuessLocal(x,ierr)
262: implicit none
264: #include "ex5f.h"
266: ! Input/output variables:
267: PetscScalar x(gxs:gxe,gys:gye)
268: PetscErrorCode ierr
270: ! Local variables:
271: PetscInt i,j
272: PetscReal temp1,temp,one,hx,hy
274: ! Set parameters
276: 0
277: one = 1.0
278: hx = one/((mx-1))
279: hy = one/((my-1))
280: temp1 = lambda/(lambda + one)
282: do 20 j=ys,ye
283: temp = (min(j-1,my-j))*hy
284: do 10 i=xs,xe
285: if (i .eq. 1 .or. j .eq. 1 &
286: & .or. i .eq. mx .or. j .eq. my) then
287: x(i,j) = 0.0
288: else
289: x(i,j) = temp1 * &
290: & 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,dummy,ierr)
316: implicit none
318: #include "ex5f.h"
320: ! Input/output variables:
321: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
322: PetscScalar x(gxs:gxe,gys:gye)
323: PetscScalar f(xs:xe,ys:ye)
324: PetscErrorCode ierr
325: PetscObject dummy
327: ! Local variables:
328: PetscScalar two,one,hx,hy
329: PetscScalar hxdhy,hydhx,sc
330: PetscScalar u,uxx,uyy
331: PetscInt i,j
333: xs = info(DMDA_LOCAL_INFO_XS)+1
334: xe = xs+info(DMDA_LOCAL_INFO_XM)-1
335: ys = info(DMDA_LOCAL_INFO_YS)+1
336: ye = ys+info(DMDA_LOCAL_INFO_YM)-1
337: mx = info(DMDA_LOCAL_INFO_MX)
338: my = info(DMDA_LOCAL_INFO_MY)
340: one = 1.0
341: two = 2.0
342: hx = one/(mx-1)
343: hy = one/(my-1)
344: sc = hx*hy*lambda
345: hxdhy = hx/hy
346: hydhx = hy/hx
348: ! Compute function over the locally owned part of the grid
350: do 20 j=ys,ye
351: do 10 i=xs,xe
352: if (i .eq. 1 .or. j .eq. 1 &
353: & .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 &
358: & - x(i-1,j) - x(i+1,j))
359: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
360: f(i,j) = uxx + uyy - sc*exp(u)
361: endif
362: 10 continue
363: 20 continue
365: call PetscLogFlops(11.0d0*ym*xm,ierr)
367: return
368: end
370: ! ---------------------------------------------------------------------
371: !
372: ! FormJacobianLocal - Computes Jacobian matrix, called by
373: ! the higher level routine FormJacobian().
374: !
375: ! Input Parameters:
376: ! x - local vector data
377: !
378: ! Output Parameters:
379: ! jac - Jacobian matrix
380: ! jac_prec - optionally different preconditioning matrix (not used here)
381: ! ierr - error code
382: !
383: ! Notes:
384: ! This routine uses standard Fortran-style computations over a 2-dim array.
385: !
386: ! Notes:
387: ! Due to grid point reordering with DMDAs, we must always work
388: ! with the local grid points, and then transform them to the new
389: ! global numbering with the "ltog" mapping (via DMDAGetGlobalIndices()).
390: ! We cannot work directly with the global numbers for the original
391: ! uniprocessor grid!
392: !
393: ! Two methods are available for imposing this transformation
394: ! when setting matrix entries:
395: ! (A) MatSetValuesLocal(), using the local ordering (including
396: ! ghost points!)
397: ! - Use DMDAGetGlobalIndices() to extract the local-to-global map
398: ! - Associate this map with the matrix by calling
399: ! MatSetLocalToGlobalMapping() once
400: ! - Set matrix entries using the local ordering
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,jac,ctx,ierr)
411: implicit none
413: #include "ex5f.h"
415: ! Input/output variables:
416: PetscScalar x(gxs:gxe,gys:gye)
417: Mat jac
418: PetscErrorCode ierr
419: integer ctx
420: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
421:
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 &
458: & .or. i .eq. mx .or. j .eq. my) then
459: ! Some f90 compilers need 4th arg to be of same type in both calls
460: col(1) = row
461: v(1) = one
462: call MatSetValuesLocal(jac,i1,row,i1,col,v, &
463: & INSERT_VALUES,ierr)
464: ! interior grid points
465: else
466: v(1) = -hxdhy
467: v(2) = -hydhx
468: v(3) = two*(hydhx + hxdhy) &
469: & - sc*lambda*exp(x(i,j))
470: v(4) = -hydhx
471: v(5) = -hxdhy
472: col(1) = row - gxm
473: col(2) = row - 1
474: col(3) = row
475: col(4) = row + 1
476: col(5) = row + gxm
477: call MatSetValuesLocal(jac,i1,row,i5,col,v, &
478: & INSERT_VALUES,ierr)
479: endif
480: 10 continue
481: 20 continue
482: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
483: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
484: return
485: end