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