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