Actual source code: ex5f.F90
1: !
2: ! This example shows how to avoid Fortran line lengths larger than 132 characters.
3: ! It avoids used of certain macros such as PetscCallA() and PetscCheckA() that
4: ! generate very long lines
5: !
6: ! We recommend starting from src/snes/tutorials/ex5f90.F90 instead of this example
7: ! because that does not have the restricted formatting that makes this version
8: ! more difficult to read
9: !
10: ! Description: This example solves a nonlinear system in parallel with SNES.
11: ! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
12: ! domain, using distributed arrays (DMDAs) to partition the parallel grid.
13: ! The command line options include:
14: ! -par <param>, where <param> indicates the nonlinearity of the problem
15: ! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)
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: ! --------------------------------------------------------------------------
33: module ex5fmodule
34: use petscsnes
35: use petscdmda
36: #include <petsc/finclude/petscsnes.h>
37: #include <petsc/finclude/petscdmda.h>
38: PetscInt xs,xe,xm,gxs,gxe,gxm
39: PetscInt ys,ye,ym,gys,gye,gym
40: PetscInt mx,my
41: PetscMPIInt rank,size
42: PetscReal lambda
43: end module ex5fmodule
45: program main
46: use ex5fmodule
47: implicit none
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: ! Variable declarations
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: !
53: ! Variables:
54: ! snes - nonlinear solver
55: ! x, r - solution, residual vectors
56: ! its - iterations for convergence
57: !
58: ! See additional variable declarations in the file ex5f.h
59: !
60: SNES snes
61: Vec x,r
62: PetscInt its,i1,i4
63: PetscErrorCode ierr
64: PetscReal lambda_max,lambda_min
65: PetscBool flg
66: DM da
68: ! Note: Any user-defined Fortran routines (such as FormJacobianLocal)
69: ! MUST be declared as external.
71: external FormInitialGuess
72: external FormFunctionLocal,FormJacobianLocal
73: external MySNESConverged
75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: ! Initialize program
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: call PetscInitialize(ierr)
80: CHKERRA(ierr)
81: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
82: CHKERRMPIA(ierr)
83: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
84: CHKERRMPIA(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: CHKERRA(ierr)
95: ! this statement is split into multiple-lines to keep lines under 132 char limit - required by 'make check'
96: if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then
97: ierr = PETSC_ERR_ARG_OUTOFRANGE;
98: SETERRA(PETSC_COMM_WORLD,ierr,'Lambda')
99: endif
101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: ! Create nonlinear solver context
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
106: CHKERRA(ierr)
108: ! Set convergence test routine if desired
110: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my_snes_convergence',flg,ierr)
111: CHKERRA(ierr)
112: if (flg) then
113: call SNESSetConvergenceTest(snes,MySNESConverged,0,PETSC_NULL_FUNCTION,ierr)
114: CHKERRA(ierr)
115: endif
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: ! Create vector data structures; set function evaluation routine
119: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Create distributed array (DMDA) to manage parallel grid and vectors
123: ! This really needs only the star-type stencil, but we use the box stencil
125: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,i4,i4,PETSC_DECIDE,PETSC_DECIDE, &
126: i1,i1, PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,da,ierr)
127: CHKERRA(ierr)
128: call DMSetFromOptions(da,ierr)
129: CHKERRA(ierr)
130: call DMSetUp(da,ierr)
131: CHKERRA(ierr)
133: ! Extract global and local vectors from DMDA; then duplicate for remaining
134: ! vectors that are the same types
136: call DMCreateGlobalVector(da,x,ierr)
137: CHKERRA(ierr)
138: call VecDuplicate(x,r,ierr)
139: CHKERRA(ierr)
141: ! Get local grid boundaries (for 2-dimensional DMDA)
143: call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
144: PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMBOUNDARYTYPE, &
145: PETSC_NULL_DMBOUNDARYTYPE,PETSC_NULL_DMDASTENCILTYPE,ierr)
146: CHKERRA(ierr)
147: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
148: CHKERRA(ierr)
149: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)
150: CHKERRA(ierr)
152: ! Here we shift the starting indices up by one so that we can easily
153: ! use the Fortran convention of 1-based indices (rather 0-based indices).
155: xs = xs+1
156: ys = ys+1
157: gxs = gxs+1
158: gys = gys+1
160: ye = ys+ym-1
161: xe = xs+xm-1
162: gye = gys+gym-1
163: gxe = gxs+gxm-1
165: ! Set function evaluation routine and vector
167: call DMDASNESSetFunctionLocal(da,INSERT_VALUES,FormFunctionLocal,da,ierr)
168: CHKERRA(ierr)
169: call DMDASNESSetJacobianLocal(da,FormJacobianLocal,da,ierr)
170: CHKERRA(ierr)
171: call SNESSetDM(snes,da,ierr)
172: CHKERRA(ierr)
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: ! Customize nonlinear solver; set runtime options
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: ! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
180: call SNESSetFromOptions(snes,ierr)
181: CHKERRA(ierr)
182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: ! Evaluate initial guess; then solve nonlinear system.
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: ! Note: The user should initialize the vector, x, with the initial guess
187: ! for the nonlinear solver prior to calling SNESSolve(). In particular,
188: ! to employ an initial guess of zero, the user should explicitly set
189: ! this vector to zero by calling VecSet().
191: call FormInitialGuess(x,ierr)
192: CHKERRA(ierr)
193: call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
194: CHKERRA(ierr)
195: call SNESGetIterationNumber(snes,its,ierr)
196: CHKERRA(ierr)
197: if (rank .eq. 0) then
198: write(6,100) its
199: endif
200: 100 format('Number of SNES iterations = ',i5)
202: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: ! Free work space. All PETSc objects should be destroyed when they
204: ! are no longer needed.
205: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207: call VecDestroy(x,ierr)
208: CHKERRA(ierr)
209: call VecDestroy(r,ierr)
210: CHKERRA(ierr)
211: call SNESDestroy(snes,ierr)
212: CHKERRA(ierr)
213: call DMDestroy(da,ierr)
214: CHKERRA(ierr)
215: call PetscFinalize(ierr)
216: CHKERRA(ierr)
217: end
219: ! ---------------------------------------------------------------------
220: !
221: ! FormInitialGuess - Forms initial approximation.
222: !
223: ! Input Parameters:
224: ! X - vector
225: !
226: ! Output Parameter:
227: ! X - vector
228: !
229: ! Notes:
230: ! This routine serves as a wrapper for the lower-level routine
231: ! "ApplicationInitialGuess", where the actual computations are
232: ! done using the standard Fortran style of treating the local
233: ! vector data as a multidimensional array over the local mesh.
234: ! This routine merely handles ghost point scatters and accesses
235: ! the local vector data via VecGetArray() and VecRestoreArray().
236: !
237: subroutine FormInitialGuess(X,ierr)
238: use ex5fmodule
239: implicit none
241: ! Input/output variables:
242: Vec X
243: PetscErrorCode ierr
244: ! Declarations for use with local arrays:
245: PetscScalar, pointer :: lx_v(:)
247: ierr = 0
249: ! Get a pointer to vector data.
250: ! - For default PETSc vectors, VecGetArray() returns a pointer to
251: ! the data array. Otherwise, the routine is implementation dependent.
252: ! - You MUST call VecRestoreArray() when you no longer need access to
253: ! the array.
254: ! - Note that the Fortran interface to VecGetArray() differs from the
255: ! C version. See the users manual for details.
257: call VecGetArray(X,lx_v,ierr)
258: CHKERRQ(ierr)
260: ! Compute initial guess over the locally owned part of the grid
262: call InitialGuessLocal(lx_v,ierr)
263: CHKERRQ(ierr)
265: ! Restore vector
267: call VecRestoreArray(X,lx_v,ierr)
268: CHKERRQ(ierr)
270: end
272: ! ---------------------------------------------------------------------
273: !
274: ! InitialGuessLocal - Computes initial approximation, called by
275: ! the higher level routine FormInitialGuess().
276: !
277: ! Input Parameter:
278: ! x - local vector data
279: !
280: ! Output Parameters:
281: ! x - local vector data
282: ! ierr - error code
283: !
284: ! Notes:
285: ! This routine uses standard Fortran-style computations over a 2-dim array.
286: !
287: subroutine InitialGuessLocal(x,ierr)
288: use ex5fmodule
289: implicit none
291: ! Input/output variables:
292: PetscScalar x(xs:xe,ys:ye)
293: PetscErrorCode ierr
295: ! Local variables:
296: PetscInt i,j
297: PetscReal temp1,temp,one,hx,hy
299: ! Set parameters
301: ierr = 0
302: one = 1.0
303: hx = one/((real(mx)-1))
304: hy = one/((real(my)-1))
305: temp1 = lambda/(lambda + one)
307: do 20 j=ys,ye
308: temp = (real(min(j-1,my-j)))*hy
309: do 10 i=xs,xe
310: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
311: x(i,j) = 0.0
312: else
313: x(i,j) = temp1 * sqrt(min(real(min(i-1,mx-i))*hx,(temp)))
314: endif
315: 10 continue
316: 20 continue
318: end
320: ! ---------------------------------------------------------------------
321: !
322: ! FormFunctionLocal - Computes nonlinear function, called by
323: ! the higher level routine FormFunction().
324: !
325: ! Input Parameter:
326: ! x - local vector data
327: !
328: ! Output Parameters:
329: ! f - local vector data, f(x)
330: ! ierr - error code
331: !
332: ! Notes:
333: ! This routine uses standard Fortran-style computations over a 2-dim array.
334: !
335: !
336: subroutine FormFunctionLocal(info,x,f,da,ierr)
337: use ex5fmodule
338: implicit none
340: DM da
342: ! Input/output variables:
343: DMDALocalInfo info
344: PetscScalar x(gxs:gxe,gys:gye)
345: PetscScalar f(xs:xe,ys:ye)
346: PetscErrorCode ierr
348: ! Local variables:
349: PetscScalar two,one,hx,hy
350: PetscScalar hxdhy,hydhx,sc
351: PetscScalar u,uxx,uyy
352: PetscInt i,j
354: xs = info%XS+1
355: xe = xs+info%XM-1
356: ys = info%YS+1
357: ye = ys+info%YM-1
358: mx = info%MX
359: my = info%MY
361: one = 1.0
362: two = 2.0
363: hx = one/(real(mx)-1)
364: hy = one/(real(my)-1)
365: sc = hx*hy*lambda
366: hxdhy = hx/hy
367: hydhx = hy/hx
369: ! Compute function over the locally owned part of the grid
371: do 20 j=ys,ye
372: do 10 i=xs,xe
373: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
374: f(i,j) = x(i,j)
375: else
376: u = x(i,j)
377: uxx = hydhx * (two*u - x(i-1,j) - x(i+1,j))
378: uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1))
379: f(i,j) = uxx + uyy - sc*exp(u)
380: endif
381: 10 continue
382: 20 continue
384: call PetscLogFlops(11.0d0*ym*xm,ierr)
385: CHKERRQ(ierr)
387: end
389: ! ---------------------------------------------------------------------
390: !
391: ! FormJacobianLocal - Computes Jacobian matrix, called by
392: ! the higher level routine FormJacobian().
393: !
394: ! Input Parameters:
395: ! x - local vector data
396: !
397: ! Output Parameters:
398: ! jac - Jacobian matrix
399: ! jac_prec - optionally different preconditioning matrix (not used here)
400: ! ierr - error code
401: !
402: ! Notes:
403: ! This routine uses standard Fortran-style computations over a 2-dim array.
404: !
405: ! Notes:
406: ! Due to grid point reordering with DMDAs, we must always work
407: ! with the local grid points, and then transform them to the new
408: ! global numbering with the "ltog" mapping
409: ! We cannot work directly with the global numbers for the original
410: ! uniprocessor grid!
411: !
412: ! Two methods are available for imposing this transformation
413: ! when setting matrix entries:
414: ! (A) MatSetValuesLocal(), using the local ordering (including
415: ! ghost points!)
416: ! by calling MatSetValuesLocal()
417: ! (B) MatSetValues(), using the global ordering
418: ! - Use DMDAGetGlobalIndices() to extract the local-to-global map
419: ! - Then apply this map explicitly yourself
420: ! - Set matrix entries using the global ordering by calling
421: ! MatSetValues()
422: ! Option (A) seems cleaner/easier in many cases, and is the procedure
423: ! used in this example.
424: !
425: subroutine FormJacobianLocal(info,x,A,jac,da,ierr)
426: use ex5fmodule
427: implicit none
429: DM da
431: ! Input/output variables:
432: PetscScalar x(gxs:gxe,gys:gye)
433: Mat A,jac
434: PetscErrorCode ierr
435: DMDALocalInfo info
437: ! Local variables:
438: PetscInt row,col(5),i,j,i1,i5
439: PetscScalar two,one,hx,hy,v(5)
440: PetscScalar hxdhy,hydhx,sc
442: ! Set parameters
444: i1 = 1
445: i5 = 5
446: one = 1.0
447: two = 2.0
448: hx = one/(real(mx)-1)
449: hy = one/(real(my)-1)
450: sc = hx*hy
451: hxdhy = hx/hy
452: hydhx = hy/hx
453: ! -Wmaybe-uninitialized
454: v = 0.0
455: col = 0
457: ! Compute entries for the locally owned part of the Jacobian.
458: ! - Currently, all PETSc parallel matrix formats are partitioned by
459: ! contiguous chunks of rows across the processors.
460: ! - Each processor needs to insert only elements that it owns
461: ! locally (but any non-local elements will be sent to the
462: ! appropriate processor during matrix assembly).
463: ! - Here, we set all entries for a particular row at once.
464: ! - We can set matrix entries either using either
465: ! MatSetValuesLocal() or MatSetValues(), as discussed above.
466: ! - Note that MatSetValues() uses 0-based row and column numbers
467: ! in Fortran as well as in C.
469: do 20 j=ys,ye
470: row = (j - gys)*gxm + xs - gxs - 1
471: do 10 i=xs,xe
472: row = row + 1
473: ! boundary points
474: if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
475: ! Some f90 compilers need 4th arg to be of same type in both calls
476: col(1) = row
477: v(1) = one
478: call MatSetValuesLocal(jac,i1,[row],i1,[col],[v],INSERT_VALUES,ierr)
479: CHKERRQ(ierr)
480: ! interior grid points
481: else
482: v(1) = -hxdhy
483: v(2) = -hydhx
484: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i,j))
485: v(4) = -hydhx
486: v(5) = -hxdhy
487: col(1) = row - gxm
488: col(2) = row - 1
489: col(3) = row
490: col(4) = row + 1
491: col(5) = row + gxm
492: call MatSetValuesLocal(jac,i1,[row],i5,[col],[v], INSERT_VALUES,ierr)
493: CHKERRQ(ierr)
494: endif
495: 10 continue
496: 20 continue
497: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
498: CHKERRQ(ierr)
499: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
500: CHKERRQ(ierr)
501: if (A .ne. jac) then
502: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
503: CHKERRQ(ierr)
504: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
505: CHKERRQ(ierr)
506: endif
507: end
509: !
510: ! Simple convergence test based on the infinity norm of the residual being small
511: !
512: subroutine MySNESConverged(snes,it,xnorm,snorm,fnorm,reason,dummy,ierr)
513: use ex5fmodule
514: implicit none
516: SNES snes
517: PetscInt it,dummy
518: PetscReal xnorm,snorm,fnorm,nrm
519: SNESConvergedReason reason
520: Vec f
521: PetscErrorCode ierr
523: call SNESGetFunction(snes,f,PETSC_NULL_FUNCTION,dummy,ierr)
524: CHKERRQ(ierr)
525: call VecNorm(f,NORM_INFINITY,nrm,ierr)
526: CHKERRQ(ierr)
527: if (nrm .le. 1.e-5) reason = SNES_CONVERGED_FNORM_ABS
529: end
531: !/*TEST
532: !
533: ! build:
534: ! requires: !complex !single
535: !
536: ! test:
537: ! nsize: 4
538: ! args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short \
539: ! -ksp_gmres_cgs_refinement_type refine_always
540: !
541: ! test:
542: ! suffix: 2
543: ! nsize: 4
544: ! args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
545: !
546: ! test:
547: ! suffix: 3
548: ! nsize: 3
549: ! args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
550: !
551: ! test:
552: ! suffix: 6
553: ! nsize: 1
554: ! args: -snes_monitor_short -my_snes_convergence
555: !
556: !TEST*/