Actual source code: ex5f.F
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: ! J - Jacobian matrix
52: ! its - iterations for convergence
53: !
54: ! See additional variable declarations in the file ex5f.h
55: !
56: SNES snes
57: Vec x,r
58: Mat J,A
59: PetscInt its,i1,i4
60: PetscErrorCode ierr
61: PetscReal lambda_max,lambda_min
62: PetscBool flg
63: #if defined(PETSC_HAVE_ADIFOR)
64: ISColoring coloring
65: PetscBool adifor_jacobian,adiformf_jacobian
66: #endif
68: ! Note: Any user-defined Fortran routines (such as FormJacobianLocal)
69: ! MUST be declared as external.
71: external FormInitialGuess
72: external FormFunctionLocal,FormJacobianLocal
73: #if defined(PETSC_HAVE_ADIFOR)
74: external g_FormFunctionLocal,m_FormFunctionLocal
75: #endif
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: ! Initialize program
79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
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_CHARACTER,'-par',lambda, &
93: & flg,ierr)
94: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then
95: if (rank .eq. 0) write(6,*) 'Lambda is out of range'
96: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
97: endif
99: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: ! Create nonlinear solver context
101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: ! Create vector data structures; set function evaluation routine
107: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: ! Create distributed array (DMDA) to manage parallel grid and vectors
111: ! This really needs only the star-type stencil, but we use the box
112: ! stencil temporarily.
113: call DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, &
114: & DMDA_BOUNDARY_NONE, &
115: & DMDA_STENCIL_STAR,i4,i4,PETSC_DECIDE,PETSC_DECIDE,i1,i1, &
116: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
118: ! Extract global and local vectors from DMDA; then duplicate for remaining
119: ! vectors that are the same types
121: call DMCreateGlobalVector(da,x,ierr)
122: call VecDuplicate(x,r,ierr)
124: ! Get local grid boundaries (for 2-dimensional DMDA)
126: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER, &
127: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
128: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
129: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
130: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
131: & PETSC_NULL_INTEGER,ierr)
132: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
133: & PETSC_NULL_INTEGER,ierr)
134: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
135: & 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 DMDASetLocalFunction(da,FormFunctionLocal,ierr)
153: call DMDASetLocalJacobian(da,FormJacobianLocal,ierr)
154: #if defined(PETSC_HAVE_ADIFOR)
155: call DMDASetLocalAdiforFunction(da, &
156: & g_FormFunctionLocal,ierr)
157: #endif
158: call SNESSetDM(snes,da,ierr)
159: call SNESSetFunction(snes,r,SNESDAFormFunction,da,ierr)
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: ! Create matrix data structure; set Jacobian evaluation routine
162: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: ! Set Jacobian matrix data structure and default Jacobian evaluation
165: ! routine. User can override with:
166: ! -snes_mf : matrix-free Newton-Krylov method with no preconditioning
167: ! (unless user explicitly sets preconditioner)
168: ! -snes_mf_operator : form preconditioning matrix as set by the user,
169: ! but use matrix-free approx for Jacobian-vector
170: ! products within Newton-Krylov method
171: !
173: call DMGetMatrix(da,MATAIJ,J,ierr)
174: #if defined(PETSC_HAVE_ADIFOR)
175: call PetscOptionsGetBool(PETSC_NULL_CHARACTER &
176: & ,'-adiformf_jacobian', &
177: & adiformf_jacobian,PETSC_NULL_INTEGER,ierr)
178: if (adiformf_jacobian .eq. 1) then
179: call DMDASetLocalAdiforMFFunction(da, &
180: & m_FormFunctionLocal,ierr)
181: call MatRegisterDAAD(ierr)
182: call MatCreateDAAD(da,A,ierr)
183: call MatDAADSetSNES(A,snes,ierr)
184: else
185: A = J
186: endif
187: #else
188: A = J
189: #endif
191: call SNESSetJacobian(snes,A,J,SNESDAComputeJacobian, &
192: & da,ierr)
193: #if defined(PETSC_HAVE_ADIFOR)
194: call PetscOptionsGetBool(PETSC_NULL_CHARACTER &
195: & ,'-adifor_jacobian', &
196: & adifor_jacobian,PETSC_NULL_INTEGER,ierr)
197: if (adifor_jacobian .eq. 1) then
198: call DMGetColoring(da,IS_COLORING_GHOSTED, &
199: & coloring,MATAIJ,ierr)
200: call MatSetColoring(J,coloring,ierr)
201: call SNESSetJacobian(snes,A,J,SNESDAComputeJacobianWithAdifor, &
202: & da,ierr)
203: call ISColoringDestroy(coloring,ierr)
204: endif
205: #endif
208: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209: ! Customize nonlinear solver; set runtime options
210: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
214: call SNESSetFromOptions(snes,ierr)
215: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: ! Evaluate initial guess; then solve nonlinear system.
217: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219: ! Note: The user should initialize the vector, x, with the initial guess
220: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
221: ! to employ an initial guess of zero, the user should explicitly set
222: ! this vector to zero by calling VecSet().
224: call FormInitialGuess(x,ierr)
225: call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr)
226: call SNESGetIterationNumber(snes,its,ierr)
227: if (rank .eq. 0) then
228: write(6,100) its
229: endif
230: 100 format('Number of Newton iterations = ',i5)
233: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234: ! Free work space. All PETSc objects should be destroyed when they
235: ! are no longer needed.
236: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238: if (A .ne. J) call MatDestroy(A,ierr)
239: call MatDestroy(J,ierr)
240: call VecDestroy(x,ierr)
241: call VecDestroy(r,ierr)
242: call SNESDestroy(snes,ierr)
243: call DMDestroy(da,ierr)
244: call PetscFinalize(ierr)
245: end
247: ! ---------------------------------------------------------------------
248: !
249: ! FormInitialGuess - Forms initial approximation.
250: !
251: ! Input Parameters:
252: ! X - vector
253: !
254: ! Output Parameter:
255: ! X - vector
256: !
257: ! Notes:
258: ! This routine serves as a wrapper for the lower-level routine
259: ! "ApplicationInitialGuess", where the actual computations are
260: ! done using the standard Fortran style of treating the local
261: ! vector data as a multidimensional array over the local mesh.
262: ! This routine merely handles ghost point scatters and accesses
263: ! the local vector data via VecGetArray() and VecRestoreArray().
264: !
265: subroutine FormInitialGuess(X,ierr)
266: implicit none
268: #include "ex5f.h"
270: ! Input/output variables:
271: Vec X
272: PetscErrorCode ierr
274: ! Declarations for use with local arrays:
275: PetscScalar lx_v(0:1)
276: PetscOffset lx_i
277: Vec localX
279: 0
281: ! Get a pointer to vector data.
282: ! - For default PETSc vectors, VecGetArray() returns a pointer to
283: ! the data array. Otherwise, the routine is implementation dependent.
284: ! - You MUST call VecRestoreArray() when you no longer need access to
285: ! the array.
286: ! - Note that the Fortran interface to VecGetArray() differs from the
287: ! C version. See the users manual for details.
289: call DMGetLocalVector(da,localX,ierr)
290: call VecGetArray(localX,lx_v,lx_i,ierr)
292: ! Compute initial guess over the locally owned part of the grid
294: call InitialGuessLocal(lx_v(lx_i),ierr)
296: ! Restore vector
298: call VecRestoreArray(localX,lx_v,lx_i,ierr)
300: ! Insert values into global vector
302: call DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X,ierr)
303: call DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X,ierr)
304: call DMRestoreLocalVector(da,localX,ierr)
305: return
306: end
308: ! ---------------------------------------------------------------------
309: !
310: ! InitialGuessLocal - Computes initial approximation, called by
311: ! the higher level routine FormInitialGuess().
312: !
313: ! Input Parameter:
314: ! x - local vector data
315: !
316: ! Output Parameters:
317: ! x - local vector data
318: ! ierr - error code
319: !
320: ! Notes:
321: ! This routine uses standard Fortran-style computations over a 2-dim array.
322: !
323: subroutine InitialGuessLocal(x,ierr)
324: implicit none
326: #include "ex5f.h"
328: ! Input/output variables:
329: PetscScalar x(gxs:gxe,gys:gye)
330: PetscErrorCode ierr
332: ! Local variables:
333: PetscInt i,j
334: PetscReal temp1,temp,one,hx,hy
336: ! Set parameters
338: 0
339: one = 1.0
340: hx = one/((mx-1))
341: hy = one/((my-1))
342: temp1 = lambda/(lambda + one)
344: do 20 j=ys,ye
345: temp = (min(j-1,my-j))*hy
346: do 10 i=xs,xe
347: if (i .eq. 1 .or. j .eq. 1 &
348: & .or. i .eq. mx .or. j .eq. my) then
349: x(i,j) = 0.0
350: else
351: x(i,j) = temp1 * &
352: & sqrt(min(min(i-1,mx-i)*hx,(temp)))
353: endif
354: 10 continue
355: 20 continue
357: return
358: end
360: ! ---------------------------------------------------------------------
361: !
362: ! FormFunctionLocal - Computes nonlinear function, called by
363: ! the higher level routine FormFunction().
364: !
365: ! Input Parameter:
366: ! x - local vector data
367: !
368: ! Output Parameters:
369: ! f - local vector data, f(x)
370: ! ierr - error code
371: !
372: ! Notes:
373: ! This routine uses standard Fortran-style computations over a 2-dim array.
374: !
375: ! Process adifor: FormFunctionLocal
376: !
377: subroutine FormFunctionLocal(info,x,f,dummy,ierr)
379: implicit none
381: #include "ex5f.h"
383: ! Input/output variables:
384: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
385: PetscScalar x(gxs:gxe,gys:gye)
386: PetscScalar f(xs:xe,ys:ye)
387: PetscErrorCode ierr
388: PetscObject dummy
390: ! Local variables:
391: PetscScalar two,one,hx,hy
392: PetscScalar hxdhy,hydhx,sc
393: PetscScalar u,uxx,uyy
394: PetscInt i,j
396: xs = info(DMDA_LOCAL_INFO_XS)+1
397: xe = xs+info(DMDA_LOCAL_INFO_XM)-1
398: ys = info(DMDA_LOCAL_INFO_YS)+1
399: ye = ys+info(DMDA_LOCAL_INFO_YM)-1
400: mx = info(DMDA_LOCAL_INFO_MX)
401: my = info(DMDA_LOCAL_INFO_MY)
403: one = 1.0
404: two = 2.0
405: hx = one/(mx-1)
406: hy = one/(my-1)
407: sc = hx*hy*lambda
408: hxdhy = hx/hy
409: hydhx = hy/hx
411: ! Compute function over the locally owned part of the grid
413: do 20 j=ys,ye
414: do 10 i=xs,xe
415: if (i .eq. 1 .or. j .eq. 1 &
416: & .or. i .eq. mx .or. j .eq. my) then
417: f(i,j) = x(i,j)
418: else
419: u = x(i,j)
420: uxx = hydhx * (two*u &
421: & - x(i-1,j) - x(i+1,j))
422: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
423: f(i,j) = uxx + uyy - sc*exp(u)
424: endif
425: 10 continue
426: 20 continue
428: call PetscLogFlops(11.0d0*ym*xm,ierr)
430: return
431: end
433: ! ---------------------------------------------------------------------
434: !
435: ! FormJacobianLocal - Computes Jacobian matrix, called by
436: ! the higher level routine FormJacobian().
437: !
438: ! Input Parameters:
439: ! x - local vector data
440: !
441: ! Output Parameters:
442: ! jac - Jacobian matrix
443: ! jac_prec - optionally different preconditioning matrix (not used here)
444: ! ierr - error code
445: !
446: ! Notes:
447: ! This routine uses standard Fortran-style computations over a 2-dim array.
448: !
449: ! Notes:
450: ! Due to grid point reordering with DMDAs, we must always work
451: ! with the local grid points, and then transform them to the new
452: ! global numbering with the "ltog" mapping (via DMDAGetGlobalIndices()).
453: ! We cannot work directly with the global numbers for the original
454: ! uniprocessor grid!
455: !
456: ! Two methods are available for imposing this transformation
457: ! when setting matrix entries:
458: ! (A) MatSetValuesLocal(), using the local ordering (including
459: ! ghost points!)
460: ! - Use DMDAGetGlobalIndices() to extract the local-to-global map
461: ! - Associate this map with the matrix by calling
462: ! MatSetLocalToGlobalMapping() once
463: ! - Set matrix entries using the local ordering
464: ! by calling MatSetValuesLocal()
465: ! (B) MatSetValues(), using the global ordering
466: ! - Use DMDAGetGlobalIndices() to extract the local-to-global map
467: ! - Then apply this map explicitly yourself
468: ! - Set matrix entries using the global ordering by calling
469: ! MatSetValues()
470: ! Option (A) seems cleaner/easier in many cases, and is the procedure
471: ! used in this example.
472: !
473: subroutine FormJacobianLocal(info,x,jac,ctx,ierr)
474: implicit none
476: #include "ex5f.h"
478: ! Input/output variables:
479: PetscScalar x(gxs:gxe,gys:gye)
480: Mat jac
481: PetscErrorCode ierr
482: integer ctx
483: DMDALocalInfo info(DMDA_LOCAL_INFO_SIZE)
484:
486: ! Local variables:
487: PetscInt row,col(5),i,j,i1,i5
488: PetscScalar two,one,hx,hy,v(5)
489: PetscScalar hxdhy,hydhx,sc
491: ! Set parameters
493: i1 = 1
494: i5 = 5
495: one = 1.0
496: two = 2.0
497: hx = one/(mx-1)
498: hy = one/(my-1)
499: sc = hx*hy
500: hxdhy = hx/hy
501: hydhx = hy/hx
503: ! Compute entries for the locally owned part of the Jacobian.
504: ! - Currently, all PETSc parallel matrix formats are partitioned by
505: ! contiguous chunks of rows across the processors.
506: ! - Each processor needs to insert only elements that it owns
507: ! locally (but any non-local elements will be sent to the
508: ! appropriate processor during matrix assembly).
509: ! - Here, we set all entries for a particular row at once.
510: ! - We can set matrix entries either using either
511: ! MatSetValuesLocal() or MatSetValues(), as discussed above.
512: ! - Note that MatSetValues() uses 0-based row and column numbers
513: ! in Fortran as well as in C.
515: do 20 j=ys,ye
516: row = (j - gys)*gxm + xs - gxs - 1
517: do 10 i=xs,xe
518: row = row + 1
519: ! boundary points
520: if (i .eq. 1 .or. j .eq. 1 &
521: & .or. i .eq. mx .or. j .eq. my) then
522: ! Some f90 compilers need 4th arg to be of same type in both calls
523: col(1) = row
524: v(1) = one
525: call MatSetValuesLocal(jac,i1,row,i1,col,v, &
526: & INSERT_VALUES,ierr)
527: ! interior grid points
528: else
529: v(1) = -hxdhy
530: v(2) = -hydhx
531: v(3) = two*(hydhx + hxdhy) &
532: & - sc*lambda*exp(x(i,j))
533: v(4) = -hydhx
534: v(5) = -hxdhy
535: col(1) = row - gxm
536: col(2) = row - 1
537: col(3) = row
538: col(4) = row + 1
539: col(5) = row + gxm
540: call MatSetValuesLocal(jac,i1,row,i5,col,v, &
541: & INSERT_VALUES,ierr)
542: endif
543: 10 continue
544: 20 continue
545: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
546: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
547: return
548: end