Actual source code: ex14f.F90
1: !
2: !
3: ! Solves a nonlinear system in parallel with a user-defined
4: ! Newton method that uses KSP to solve the linearized Newton systems. This solver
5: ! is a very simplistic inexact Newton method. The intent of this code is to
6: ! demonstrate the repeated solution of linear systems with the same nonzero pattern.
7: !
8: ! This is NOT the recommended approach for solving nonlinear problems with PETSc!
9: ! We urge users to employ the SNES component for solving nonlinear problems whenever
10: ! possible, as it offers many advantages over coding nonlinear solvers independently.
11: !
12: ! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
13: ! domain, using distributed arrays (DMDAs) to partition the parallel grid.
14: !
15: ! The command line options include:
16: ! -par <parameter>, where <parameter> indicates the problem's nonlinearity
17: ! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)
18: ! -mx <xg>, where <xg> = number of grid points in the x-direction
19: ! -my <yg>, where <yg> = number of grid points in the y-direction
20: ! -Nx <npx>, where <npx> = number of processors in the x-direction
21: ! -Ny <npy>, where <npy> = number of processors in the y-direction
22: ! -mf use matrix-free for matrix-vector product
23: !
25: ! ------------------------------------------------------------------------
26: !
27: ! Solid Fuel Ignition (SFI) problem. This problem is modeled by
28: ! the partial differential equation
29: !
30: ! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
31: !
32: ! with boundary conditions
33: !
34: ! u = 0 for x = 0, x = 1, y = 0, y = 1.
35: !
36: ! A finite difference approximation with the usual 5-point stencil
37: ! is used to discretize the boundary value problem to obtain a nonlinear
38: ! system of equations.
39: !
40: ! The SNES version of this problem is: snes/tutorials/ex5f.F
41: !
42: ! -------------------------------------------------------------------------
43: #include <petsc/finclude/petscdmda.h>
44: #include <petsc/finclude/petscksp.h>
45: module ex14fmodule
46: use petscdm
47: use petscdmda
48: use petscksp
49: implicit none
51: Vec localX
52: PetscInt mx, my
53: Mat B
54: DM da
55: contains
56: ! -------------------------------------------------------------------
57: !
58: ! MyMult - user provided matrix multiply
59: !
60: ! Input Parameters:
61: !. X - input vector
62: !
63: ! Output Parameter:
64: !. F - function vector
65: !
66: subroutine MyMult(J, X, F, ierr)
68: Mat J
69: Vec X, F
70: PetscErrorCode ierr
71: !
72: ! Here we use the actual formed matrix B; users would
73: ! instead write their own matrix-vector product routine
74: !
75: PetscCall(MatMult(B, X, F, ierr))
76: end
77: ! -------------------------------------------------------------------
78: !
79: ! FormInitialGuess - Forms initial approximation.
80: !
81: ! Input Parameters:
82: ! X - vector
83: !
84: ! Output Parameter:
85: ! X - vector
86: !
87: subroutine FormInitialGuess(X, ierr)
89: PetscErrorCode ierr
90: Vec X
91: PetscInt i, j, row
92: PetscInt xs, ys, xm
93: PetscInt ym
94: PetscReal one, lambda, temp1, temp, hx, hy
95: PetscScalar, pointer ::xx(:)
97: one = 1.0
98: lambda = 6.0
99: hx = one/(mx - 1)
100: hy = one/(my - 1)
101: temp1 = lambda/(lambda + one)
103: ! Get a pointer to vector data.
104: ! - VecGetArray() returns a pointer to the data array.
105: ! - You MUST call VecRestoreArray() when you no longer need access to
106: ! the array.
107: PetscCall(VecGetArray(X, xx, ierr))
109: ! Get local grid boundaries (for 2-dimensional DMDA):
110: ! xs, ys - starting grid indices (no ghost points)
111: ! xm, ym - widths of local grid (no ghost points)
113: PetscCall(DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
115: ! Compute initial guess over the locally owned part of the grid
117: do j = ys, ys + ym - 1
118: temp = (min(j, my - j - 1))*hy
119: do i = xs, xs + xm - 1
120: row = i - xs + (j - ys)*xm + 1
121: if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
122: xx(row) = 0.0
123: continue
124: end if
125: xx(row) = temp1*sqrt(min((min(i, mx - i - 1))*hx, temp))
126: end do
127: end do
129: ! Restore vector
131: PetscCall(VecRestoreArray(X, xx, ierr))
132: end
134: ! -------------------------------------------------------------------
135: !
136: ! ComputeFunction - Evaluates nonlinear function, F(x).
137: !
138: ! Input Parameters:
139: !. X - input vector
140: !
141: ! Output Parameter:
142: !. F - function vector
143: !
144: subroutine ComputeFunction(X, F, ierr)
146: Vec X, F
147: PetscInt gys, gxm, gym
148: PetscErrorCode ierr
149: PetscInt i, j, row, xs, ys, xm, ym, gxs
150: PetscInt rowf
151: PetscReal two, one, lambda, hx
152: PetscReal hy, hxdhy, hydhx, sc
153: PetscScalar u, uxx, uyy
154: PetscScalar, pointer ::xx(:), ff(:)
156: two = 2.0
157: one = 1.0
158: lambda = 6.0
160: hx = one/(mx - 1)
161: hy = one/(my - 1)
162: sc = hx*hy*lambda
163: hxdhy = hx/hy
164: hydhx = hy/hx
166: ! Scatter ghost points to local vector, using the 2-step process
167: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
168: ! By placing code between these two statements, computations can be
169: ! done while messages are in transition.
170: !
171: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
172: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
174: ! Get pointers to vector data
176: PetscCall(VecGetArrayRead(localX, xx, ierr))
177: PetscCall(VecGetArray(F, ff, ierr))
179: ! Get local grid boundaries
181: PetscCall(DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
182: PetscCall(DMDAGetGhostCorners(da, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
184: ! Compute function over the locally owned part of the grid
185: rowf = 0
186: do j = ys, ys + ym - 1
188: row = (j - gys)*gxm + xs - gxs
189: do i = xs, xs + xm - 1
190: row = row + 1
191: rowf = rowf + 1
193: if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
194: ff(rowf) = xx(row)
195: cycle
196: end if
197: u = xx(row)
198: uxx = (two*u - xx(row - 1) - xx(row + 1))*hydhx
199: uyy = (two*u - xx(row - gxm) - xx(row + gxm))*hxdhy
200: ff(rowf) = uxx + uyy - sc*exp(u)
201: end do
202: end do
204: ! Restore vectors
206: PetscCall(VecRestoreArrayRead(localX, xx, ierr))
207: PetscCall(VecRestoreArray(F, ff, ierr))
208: end
210: ! -------------------------------------------------------------------
211: !
212: ! ComputeJacobian - Evaluates Jacobian matrix.
213: !
214: ! Input Parameters:
215: ! x - input vector
216: !
217: ! Output Parameters:
218: ! jac - Jacobian matrix
219: !
220: ! Notes:
221: ! Due to grid point reordering with DMDAs, we must always work
222: ! with the local grid points, and then transform them to the new
223: ! global numbering with the 'ltog' mapping
224: ! We cannot work directly with the global numbers for the original
225: ! uniprocessor grid!
226: !
227: subroutine ComputeJacobian(X, jac, ierr)
229: Vec X
230: Mat jac
231: PetscErrorCode ierr
232: PetscInt xs, ys, xm, ym
233: PetscInt gxs, gys, gxm, gym
234: PetscInt grow(1), i, j
235: PetscInt row, ione
236: PetscInt col(5), ifive
237: PetscScalar two, one, lambda
238: PetscScalar v(5), hx, hy, hxdhy
239: PetscScalar hydhx, sc
240: ISLocalToGlobalMapping ltogm
241: PetscScalar, pointer ::xx(:)
242: PetscInt, pointer ::ltog(:)
244: ione = 1
245: ifive = 5
246: one = 1.0
247: two = 2.0
248: hx = one/(mx - 1)
249: hy = one/(my - 1)
250: sc = hx*hy
251: hxdhy = hx/hy
252: hydhx = hy/hx
253: lambda = 6.0
255: ! Scatter ghost points to local vector, using the 2-step process
256: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
257: ! By placing code between these two statements, computations can be
258: ! done while messages are in transition.
260: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
261: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
263: ! Get pointer to vector data
265: PetscCall(VecGetArrayRead(localX, xx, ierr))
267: ! Get local grid boundaries
269: PetscCall(DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
270: PetscCall(DMDAGetGhostCorners(da, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
272: ! Get the global node numbers for all local nodes, including ghost points
274: PetscCall(DMGetLocalToGlobalMapping(da, ltogm, ierr))
275: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, ltog, ierr))
277: ! Compute entries for the locally owned part of the Jacobian.
278: ! - Currently, all PETSc parallel matrix formats are partitioned by
279: ! contiguous chunks of rows across the processors. The 'grow'
280: ! parameter computed below specifies the global row number
281: ! corresponding to each local grid point.
282: ! - Each processor needs to insert only elements that it owns
283: ! locally (but any non-local elements will be sent to the
284: ! appropriate processor during matrix assembly).
285: ! - Always specify global row and columns of matrix entries.
286: ! - Here, we set all entries for a particular row at once.
288: do j = ys, ys + ym - 1
289: row = (j - gys)*gxm + xs - gxs
290: do i = xs, xs + xm - 1
291: row = row + 1
292: grow(1) = ltog(row)
293: if (i == 0 .or. j == 0 .or. i == (mx - 1) .or. j == (my - 1)) then
294: PetscCall(MatSetValues(jac, ione, grow, ione, grow, [one], INSERT_VALUES, ierr))
295: cycle
296: end if
297: v(1) = -hxdhy
298: col(1) = ltog(row - gxm)
299: v(2) = -hydhx
300: col(2) = ltog(row - 1)
301: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(xx(row))
302: col(3) = grow(1)
303: v(4) = -hydhx
304: col(4) = ltog(row + 1)
305: v(5) = -hxdhy
306: col(5) = ltog(row + gxm)
307: PetscCall(MatSetValues(jac, ione, grow, ifive, col, v, INSERT_VALUES, ierr))
308: end do
309: end do
311: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, ltog, ierr))
313: ! Assemble matrix, using the 2-step process:
314: ! MatAssemblyBegin(), MatAssemblyEnd().
315: ! By placing code between these two statements, computations can be
316: ! done while messages are in transition.
318: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
319: PetscCall(VecRestoreArrayRead(localX, xx, ierr))
320: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
321: end
322: end module ex14fmodule
324: program main
325: use ex14fmodule
326: implicit none
328: MPIU_Comm comm
329: Vec X, Y, F
330: Mat J
331: KSP ksp
333: PetscInt Nx, Ny, N, ifive, ithree
334: PetscBool flg, nooutput, usemf
336: ! --------------- Data to define nonlinear solver --------------
337: PetscReal rtol, ttol
338: PetscReal fnorm, ynorm, xnorm
339: PetscInt max_nonlin_its, one
340: PetscInt lin_its
341: PetscInt i, m
342: PetscScalar mone
343: PetscErrorCode ierr
345: mone = -1.0
346: rtol = 1.e-8
347: max_nonlin_its = 10
348: one = 1
349: ifive = 5
350: ithree = 3
352: PetscCallA(PetscInitialize(ierr))
353: comm = PETSC_COMM_WORLD
355: ! Initialize problem parameters
357: !
358: mx = 4
359: my = 4
360: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
361: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))
362: N = mx*my
364: nooutput = .false.
365: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-no_output', nooutput, ierr))
367: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
368: ! Create linear solver context
369: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371: PetscCallA(KSPCreate(comm, ksp, ierr))
373: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
374: ! Create vector data structures
375: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
377: !
378: ! Create distributed array (DMDA) to manage parallel grid and vectors
379: !
380: Nx = PETSC_DECIDE
381: Ny = PETSC_DECIDE
382: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-Nx', Nx, flg, ierr))
383: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-Ny', Ny, flg, ierr))
384: PetscCallA(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, mx, my, Nx, Ny, one, one, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr))
385: PetscCallA(DMSetFromOptions(da, ierr))
386: PetscCallA(DMSetUp(da, ierr))
387: !
388: ! Extract global and local vectors from DMDA then duplicate for remaining
389: ! vectors that are the same types
390: !
391: PetscCallA(DMCreateGlobalVector(da, X, ierr))
392: PetscCallA(DMCreateLocalVector(da, localX, ierr))
393: PetscCallA(VecDuplicate(X, F, ierr))
394: PetscCallA(VecDuplicate(X, Y, ierr))
396: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
397: ! Create matrix data structure for Jacobian
398: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
399: !
400: ! Note: For the parallel case, vectors and matrices MUST be partitioned
401: ! accordingly. When using distributed arrays (DMDAs) to create vectors,
402: ! the DMDAs determine the problem partitioning. We must explicitly
403: ! specify the local matrix dimensions upon its creation for compatibility
404: ! with the vector distribution.
405: !
406: ! Note: Here we only approximately preallocate storage space for the
407: ! Jacobian. See the users manual for a discussion of better techniques
408: ! for preallocating matrix memory.
409: !
410: PetscCallA(VecGetLocalSize(X, m, ierr))
411: PetscCallA(MatCreateFromOptions(comm, PETSC_NULL_CHARACTER, one, m, m, N, N, B, ierr))
413: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414: ! if usemf is on then matrix-vector product is done via matrix-free
415: ! approach. Note this is just an example, and not realistic because
416: ! we still use the actual formed matrix, but in reality one would
417: ! provide their own subroutine that would directly do the matrix
418: ! vector product and call MatMult()
419: ! Note: we put B into a module so it will be visible to the
420: ! mymult() routine
421: usemf = .false.
422: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mf', usemf, ierr))
423: if (usemf) then
424: PetscCallA(MatCreateShell(comm, m, m, N, N, PETSC_NULL_INTEGER, J, ierr))
425: PetscCallA(MatShellSetOperation(J, MATOP_MULT, mymult, ierr))
426: else
427: ! If not doing matrix-free then matrix operator, J, and matrix used
428: ! to construct preconditioner, B, are the same
429: J = B
430: end if
432: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
433: ! Customize linear solver set runtime options
434: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
435: !
436: ! Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
437: !
438: PetscCallA(KSPSetFromOptions(ksp, ierr))
440: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
441: ! Evaluate initial guess
442: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
444: PetscCallA(FormInitialGuess(X, ierr))
445: PetscCallA(ComputeFunction(X, F, ierr))
446: PetscCallA(VecNorm(F, NORM_2, fnorm, ierr))
447: ttol = fnorm*rtol
448: if (.not. nooutput) then
449: print *, 'Initial function norm ', fnorm
450: end if
452: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
453: ! Solve nonlinear system with a user-defined method
454: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
456: ! This solver is a very simplistic inexact Newton method, with no
457: ! no damping strategies or bells and whistles. The intent of this code
458: ! is merely to demonstrate the repeated solution with KSP of linear
459: ! systems with the same nonzero structure.
460: !
461: ! This is NOT the recommended approach for solving nonlinear problems
462: ! with PETSc! We urge users to employ the SNES component for solving
463: ! nonlinear problems whenever possible with application codes, as it
464: ! offers many advantages over coding nonlinear solvers independently.
466: do i = 0, max_nonlin_its
468: ! Compute the Jacobian matrix. See the comments in this routine for
469: ! important information about setting the flag mat_flag.
471: PetscCallA(ComputeJacobian(X, B, ierr))
473: ! Solve J Y = F, where J is the Jacobian matrix.
474: ! - First, set the KSP linear operators. Here the matrix that
475: ! defines the linear system also serves as the preconditioning
476: ! matrix.
477: ! - Then solve the Newton system.
479: PetscCallA(KSPSetOperators(ksp, J, B, ierr))
480: PetscCallA(KSPSolve(ksp, F, Y, ierr))
482: ! Compute updated iterate
484: PetscCallA(VecNorm(Y, NORM_2, ynorm, ierr))
485: PetscCallA(VecAYPX(Y, mone, X, ierr))
486: PetscCallA(VecCopy(Y, X, ierr))
487: PetscCallA(VecNorm(X, NORM_2, xnorm, ierr))
488: PetscCallA(KSPGetIterationNumber(ksp, lin_its, ierr))
489: if (.not. nooutput) then
490: print *, 'linear solve iterations = ', lin_its, ' xnorm = ', xnorm, ' ynorm = ', ynorm
491: end if
493: ! Evaluate nonlinear function at new location
495: PetscCallA(ComputeFunction(X, F, ierr))
496: PetscCallA(VecNorm(F, NORM_2, fnorm, ierr))
497: if (.not. nooutput) then
498: print *, 'Iteration ', i + 1, ' function norm', fnorm
499: end if
501: ! Test for convergence
503: if (fnorm <= ttol) then
504: if (.not. nooutput) then
505: print *, 'Converged: function norm ', fnorm, ' tolerance ', ttol
506: end if
507: exit
508: end if
509: end do
511: write (6, 100) i + 1
512: 100 format('Number of SNES iterations =', I2)
514: ! Check if mymult() produces a linear operator
515: if (usemf) then
516: N = 5
517: PetscCallA(MatIsLinear(J, N, flg, ierr))
518: if (.not. flg) then
519: print *, 'IsLinear', flg
520: end if
521: end if
523: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
524: ! Free work space. All PETSc objects should be destroyed when they
525: ! are no longer needed.
526: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
528: PetscCallA(MatDestroy(B, ierr))
529: if (usemf) then
530: PetscCallA(MatDestroy(J, ierr))
531: end if
532: PetscCallA(VecDestroy(localX, ierr))
533: PetscCallA(VecDestroy(X, ierr))
534: PetscCallA(VecDestroy(Y, ierr))
535: PetscCallA(VecDestroy(F, ierr))
536: PetscCallA(KSPDestroy(ksp, ierr))
537: PetscCallA(DMDestroy(da, ierr))
538: PetscCallA(PetscFinalize(ierr))
539: end
541: !/*TEST
542: !
543: ! test:
544: ! args: -no_output -ksp_gmres_cgs_refinement_type refine_always
545: ! output_file: output/ex14_1.out
546: ! requires: !single
547: !
548: !TEST*/