Actual source code: ex5f.F
petsc-3.7.5 2017-01-01
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*/
15: !
16: ! --------------------------------------------------------------------------
17: !
18: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
19: ! the partial differential equation
20: !
21: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
22: !
23: ! with boundary conditions
24: !
25: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
26: !
27: ! A finite difference approximation with the usual 5-point stencil
28: ! is used to discretize the boundary value problem to obtain a nonlinear
29: ! system of equations.
30: !
31: ! --------------------------------------------------------------------------
33: program main
34: implicit none
35: !
36: ! We place common blocks, variable declarations, and other include files
37: ! needed for this code in the single file ex5f.h. We then need to include
38: ! only this file throughout the various routines in this program. See
39: ! additional comments in the file ex5f.h.
40: !
41: #include "ex5f.h"
43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: ! Variable declarations
45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: !
47: ! Variables:
48: ! snes - nonlinear solver
49: ! x, r - solution, residual vectors
50: ! its - iterations for convergence
51: !
52: ! See additional variable declarations in the file ex5f.h
53: !
54: SNES snes
55: Vec x,r
56: PetscInt its,i1,i4
57: PetscErrorCode ierr
58: PetscReal lambda_max,lambda_min
59: PetscBool flg
62: ! Note: Any user-defined Fortran routines (such as FormJacobianLocal)
63: ! MUST be declared as external.
65: external FormInitialGuess
66: external FormFunctionLocal,FormJacobianLocal
68: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: ! Initialize program
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
73: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
74: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
76: ! Initialize problem parameters
78: i1 = 1
79: i4 = -4
80: lambda_max = 6.81
81: lambda_min = 0.0
82: lambda = 6.0
83: call PetscOptionsGetReal(PETSC_NULL_OBJECT, &
84: & PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr)
85: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then
86: if (rank .eq. 0) write(6,*) 'Lambda is out of range'
87: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
88: endif
90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: ! Create nonlinear solver context
92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: ! Create vector data structures; set function evaluation routine
98: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: ! Create distributed array (DMDA) to manage parallel grid and vectors
102: ! This really needs only the star-type stencil, but we use the box
103: ! stencil temporarily.
104: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, &
105: & DM_BOUNDARY_NONE, &
106: & DMDA_STENCIL_STAR,i4,i4,PETSC_DECIDE,PETSC_DECIDE,i1,i1, &
107: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
109: ! Extract global and local vectors from DMDA; then duplicate for remaining
110: ! vectors that are the same types
112: call DMCreateGlobalVector(da,x,ierr)
113: call VecDuplicate(x,r,ierr)
115: ! Get local grid boundaries (for 2-dimensional DMDA)
117: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER, &
118: & PETSC_NULL_INTEGER,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,ierr)
123: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
124: & PETSC_NULL_INTEGER,ierr)
125: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
126: & PETSC_NULL_INTEGER,ierr)
128: ! Here we shift the starting indices up by one so that we can easily
129: ! use the Fortran convention of 1-based indices (rather 0-based indices).
131: xs = xs+1
132: ys = ys+1
133: gxs = gxs+1
134: gys = gys+1
136: ye = ys+ym-1
137: xe = xs+xm-1
138: gye = gys+gym-1
139: gxe = gxs+gxm-1
141: ! Set function evaluation routine and vector
143: call DMDASNESSetFunctionLocal(da,INSERT_VALUES,FormFunctionLocal, &
144: & PETSC_NULL_OBJECT,ierr)
145: call DMDASNESSetJacobianLocal(da,FormJacobianLocal, &
146: & PETSC_NULL_OBJECT,ierr)
147: call SNESSetDM(snes,da,ierr)
149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: ! Customize nonlinear solver; set runtime options
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
155: call SNESSetFromOptions(snes,ierr)
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! Evaluate initial guess; then solve nonlinear system.
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: ! Note: The user should initialize the vector, x, with the initial guess
161: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
162: ! to employ an initial guess of zero, the user should explicitly set
163: ! this vector to zero by calling VecSet().
165: call FormInitialGuess(x,ierr)
166: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
167: call SNESGetIterationNumber(snes,its,ierr)
168: if (rank .eq. 0) then
169: write(6,100) its
170: endif
171: 100 format('Number of SNES iterations = ',i5)
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: ! Free work space. All PETSc objects should be destroyed when they
176: ! are no longer needed.
177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: call VecDestroy(x,ierr)
180: call VecDestroy(r,ierr)
181: call SNESDestroy(snes,ierr)
182: call DMDestroy(da,ierr)
183: call PetscFinalize(ierr)
184: end
186: ! ---------------------------------------------------------------------
187: !
188: ! FormInitialGuess - Forms initial approximation.
189: !
190: ! Input Parameters:
191: ! X - vector
192: !
193: ! Output Parameter:
194: ! X - vector
195: !
196: ! Notes:
197: ! This routine serves as a wrapper for the lower-level routine
198: ! "ApplicationInitialGuess", where the actual computations are
199: ! done using the standard Fortran style of treating the local
200: ! vector data as a multidimensional array over the local mesh.
201: ! This routine merely handles ghost point scatters and accesses
202: ! the local vector data via VecGetArray() and VecRestoreArray().
203: !
204: subroutine FormInitialGuess(X,ierr)
205: implicit none
207: #include "ex5f.h"
209: ! Input/output variables:
210: Vec X
211: PetscErrorCode ierr
213: ! Declarations for use with local arrays:
214: PetscScalar lx_v(0:1)
215: PetscOffset lx_i
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 VecGetArray(X,lx_v,lx_i,ierr)
229: ! Compute initial guess over the locally owned part of the grid
231: call InitialGuessLocal(lx_v(lx_i),ierr)
233: ! Restore vector
235: call VecRestoreArray(X,lx_v,lx_i,ierr)
237: return
238: end
240: ! ---------------------------------------------------------------------
241: !
242: ! InitialGuessLocal - Computes initial approximation, called by
243: ! the higher level routine FormInitialGuess().
244: !
245: ! Input Parameter:
246: ! x - local vector data
247: !
248: ! Output Parameters:
249: ! x - local vector data
250: ! ierr - error code
251: !
252: ! Notes:
253: ! This routine uses standard Fortran-style computations over a 2-dim array.
254: !
255: subroutine InitialGuessLocal(x,ierr)
256: implicit none
258: #include "ex5f.h"
260: ! Input/output variables:
261: PetscScalar x(xs:xe,ys:ye)
262: PetscErrorCode ierr
264: ! Local variables:
265: PetscInt i,j
266: PetscReal temp1,temp,one,hx,hy
268: ! Set parameters
270: 0
271: one = 1.0
272: hx = one/((mx-1))
273: hy = one/((my-1))
274: temp1 = lambda/(lambda + one)
276: do 20 j=ys,ye
277: temp = (min(j-1,my-j))*hy
278: do 10 i=xs,xe
279: if (i .eq. 1 .or. j .eq. 1 &
280: & .or. i .eq. mx .or. j .eq. my) then
281: x(i,j) = 0.0
282: else
283: x(i,j) = temp1 * &
284: & sqrt(min(min(i-1,mx-i)*hx,(temp)))
285: endif
286: 10 continue
287: 20 continue
289: return
290: end
292: ! ---------------------------------------------------------------------
293: !
294: ! FormFunctionLocal - Computes nonlinear function, called by
295: ! the higher level routine FormFunction().
296: !
297: ! Input Parameter:
298: ! x - local vector data
299: !
300: ! Output Parameters:
301: ! f - local vector data, f(x)
302: ! ierr - error code
303: !
304: ! Notes:
305: ! This routine uses standard Fortran-style computations over a 2-dim array.
306: !
307: !
308: subroutine FormFunctionLocal(info,x,f,dummy,ierr)
310: implicit none
312: #include "ex5f.h"
314: ! Input/output variables:
315: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
316: PetscScalar x(gxs:gxe,gys:gye)
317: PetscScalar f(xs:xe,ys:ye)
318: PetscErrorCode ierr
319: PetscObject dummy
321: ! Local variables:
322: PetscScalar two,one,hx,hy
323: PetscScalar hxdhy,hydhx,sc
324: PetscScalar u,uxx,uyy
325: PetscInt i,j
327: xs = info(DMDA_LOCAL_INFO_XS)+1
328: xe = xs+info(DMDA_LOCAL_INFO_XM)-1
329: ys = info(DMDA_LOCAL_INFO_YS)+1
330: ye = ys+info(DMDA_LOCAL_INFO_YM)-1
331: mx = info(DMDA_LOCAL_INFO_MX)
332: my = info(DMDA_LOCAL_INFO_MY)
334: one = 1.0
335: two = 2.0
336: hx = one/(mx-1)
337: hy = one/(my-1)
338: sc = hx*hy*lambda
339: hxdhy = hx/hy
340: hydhx = hy/hx
342: ! Compute function over the locally owned part of the grid
344: do 20 j=ys,ye
345: do 10 i=xs,xe
346: if (i .eq. 1 .or. j .eq. 1 &
347: & .or. i .eq. mx .or. j .eq. my) then
348: f(i,j) = x(i,j)
349: else
350: u = x(i,j)
351: uxx = hydhx * (two*u &
352: & - x(i-1,j) - x(i+1,j))
353: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
354: f(i,j) = uxx + uyy - sc*exp(u)
355: endif
356: 10 continue
357: 20 continue
359: call PetscLogFlops(11.0d0*ym*xm,ierr)
361: return
362: end
364: ! ---------------------------------------------------------------------
365: !
366: ! FormJacobianLocal - Computes Jacobian matrix, called by
367: ! the higher level routine FormJacobian().
368: !
369: ! Input Parameters:
370: ! x - local vector data
371: !
372: ! Output Parameters:
373: ! jac - Jacobian matrix
374: ! jac_prec - optionally different preconditioning matrix (not used here)
375: ! ierr - error code
376: !
377: ! Notes:
378: ! This routine uses standard Fortran-style computations over a 2-dim array.
379: !
380: ! Notes:
381: ! Due to grid point reordering with DMDAs, we must always work
382: ! with the local grid points, and then transform them to the new
383: ! global numbering with the "ltog" mapping
384: ! We cannot work directly with the global numbers for the original
385: ! uniprocessor grid!
386: !
387: ! Two methods are available for imposing this transformation
388: ! when setting matrix entries:
389: ! (A) MatSetValuesLocal(), using the local ordering (including
390: ! ghost points!)
391: ! by calling MatSetValuesLocal()
392: ! (B) MatSetValues(), using the global ordering
393: ! - Use DMDAGetGlobalIndices() to extract the local-to-global map
394: ! - Then apply this map explicitly yourself
395: ! - Set matrix entries using the global ordering by calling
396: ! MatSetValues()
397: ! Option (A) seems cleaner/easier in many cases, and is the procedure
398: ! used in this example.
399: !
400: subroutine FormJacobianLocal(info,x,A,jac,ctx,ierr)
401: implicit none
403: #include "ex5f.h"
405: ! Input/output variables:
406: PetscScalar x(gxs:gxe,gys:gye)
407: Mat A,jac
408: PetscErrorCode ierr
409: integer ctx
410: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
413: ! Local variables:
414: PetscInt row,col(5),i,j,i1,i5
415: PetscScalar two,one,hx,hy,v(5)
416: PetscScalar hxdhy,hydhx,sc
418: ! Set parameters
420: i1 = 1
421: i5 = 5
422: one = 1.0
423: two = 2.0
424: hx = one/(mx-1)
425: hy = one/(my-1)
426: sc = hx*hy
427: hxdhy = hx/hy
428: hydhx = hy/hx
430: ! Compute entries for the locally owned part of the Jacobian.
431: ! - Currently, all PETSc parallel matrix formats are partitioned by
432: ! contiguous chunks of rows across the processors.
433: ! - Each processor needs to insert only elements that it owns
434: ! locally (but any non-local elements will be sent to the
435: ! appropriate processor during matrix assembly).
436: ! - Here, we set all entries for a particular row at once.
437: ! - We can set matrix entries either using either
438: ! MatSetValuesLocal() or MatSetValues(), as discussed above.
439: ! - Note that MatSetValues() uses 0-based row and column numbers
440: ! in Fortran as well as in C.
442: do 20 j=ys,ye
443: row = (j - gys)*gxm + xs - gxs - 1
444: do 10 i=xs,xe
445: row = row + 1
446: ! boundary points
447: if (i .eq. 1 .or. j .eq. 1 &
448: & .or. i .eq. mx .or. j .eq. my) then
449: ! Some f90 compilers need 4th arg to be of same type in both calls
450: col(1) = row
451: v(1) = one
452: call MatSetValuesLocal(jac,i1,row,i1,col,v, &
453: & INSERT_VALUES,ierr)
454: ! interior grid points
455: else
456: v(1) = -hxdhy
457: v(2) = -hydhx
458: v(3) = two*(hydhx + hxdhy) &
459: & - sc*lambda*exp(x(i,j))
460: v(4) = -hydhx
461: v(5) = -hxdhy
462: col(1) = row - gxm
463: col(2) = row - 1
464: col(3) = row
465: col(4) = row + 1
466: col(5) = row + gxm
467: call MatSetValuesLocal(jac,i1,row,i5,col,v, &
468: & INSERT_VALUES,ierr)
469: endif
470: 10 continue
471: 20 continue
472: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
473: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
474: if (A .ne. jac) then
475: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
476: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
477: endif
478: return
479: end