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 >= lambda_max .or. lambda <= lambda_min) then
97: ierr = PETSC_ERR_ARG_OUTOFRANGE
98: SETERRA(PETSC_COMM_WORLD, ierr, 'Lambda')
99: end if
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: end if
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 == 0) then
198: write (6, 100) its
199: end if
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 == 1 .or. j == 1 .or. i == mx .or. j == 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: end if
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 == 1 .or. j == 1 .or. i == mx .or. j == 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: end if
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 matrix used to construct the preconditioner (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 == 1 .or. j == 1 .or. i == mx .or. j == 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: end if
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 /= 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: end if
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 <= 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*/