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: PetscScalar, parameter :: two = 2.0, one = 1.0, mone = -1.0, lambda = 6.0
56: contains
57: ! -------------------------------------------------------------------
58: !
59: ! MyMult - user provided matrix multiply
60: !
61: ! Input Parameters:
62: !. X - input vector
63: !
64: ! Output Parameter:
65: !. F - function vector
66: !
67: subroutine MyMult(J, X, F, ierr)
69: Mat J
70: Vec X, F
71: PetscErrorCode ierr
72: !
73: ! Here we use the actual formed matrix B; users would
74: ! instead write their own matrix-vector product routine
75: !
76: PetscCall(MatMult(B, X, F, ierr))
77: end
78: ! -------------------------------------------------------------------
79: !
80: ! FormInitialGuess - Forms initial approximation.
81: !
82: ! Input Parameters:
83: ! X - vector
84: !
85: ! Output Parameter:
86: ! X - vector
87: !
88: subroutine FormInitialGuess(X, ierr)
90: PetscErrorCode, intent(out) :: ierr
91: Vec X
92: PetscInt i, j, row
93: PetscInt xs, ys, xm
94: PetscInt ym
95: PetscReal temp1, temp, hx, hy
96: PetscScalar, pointer ::xx(:)
98: hx = real(one, PETSC_REAL_KIND)/(mx - 1)
99: hy = real(one, PETSC_REAL_KIND)/(my - 1)
100: temp1 = real(lambda/(lambda + one), PETSC_REAL_KIND)
102: ! Get a pointer to vector data.
103: ! - VecGetArray() returns a pointer to the data array.
104: ! - You MUST call VecRestoreArray() when you no longer need access to
105: ! the array.
106: PetscCall(VecGetArray(X, xx, ierr))
108: ! Get local grid boundaries (for 2-dimensional DMDA):
109: ! xs, ys - starting grid indices (no ghost points)
110: ! xm, ym - widths of local grid (no ghost points)
112: PetscCall(DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
114: ! Compute initial guess over the locally owned part of the grid
116: do j = ys, ys + ym - 1
117: temp = (min(j, my - j - 1))*hy
118: do i = xs, xs + xm - 1
119: row = i - xs + (j - ys)*xm + 1
120: if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
121: xx(row) = 0.0
122: continue
123: end if
124: xx(row) = temp1*sqrt(min((min(i, mx - i - 1))*hx, temp))
125: end do
126: end do
128: ! Restore vector
130: PetscCall(VecRestoreArray(X, xx, ierr))
131: end
133: ! -------------------------------------------------------------------
134: !
135: ! ComputeFunction - Evaluates nonlinear function, F(x).
136: !
137: ! Input Parameters:
138: !. X - input vector
139: !
140: ! Output Parameter:
141: !. F - function vector
142: !
143: subroutine ComputeFunction(X, F, ierr)
145: Vec X, F
146: PetscInt gys, gxm, gym
147: PetscErrorCode, intent(out) :: ierr
148: PetscInt i, j, row, xs, ys, xm, ym, gxs
149: PetscInt rowf
150: PetscReal hx, hy, hxdhy, hydhx, sc
151: PetscScalar u, uxx, uyy
152: PetscScalar, pointer ::xx(:), ff(:)
154: hx = real(one, PETSC_REAL_KIND)/(mx - 1)
155: hy = real(one, PETSC_REAL_KIND)/(my - 1)
156: sc = hx*hy*real(lambda, PETSC_REAL_KIND)
157: hxdhy = hx/hy
158: hydhx = hy/hx
160: ! Scatter ghost points to local vector, using the 2-step process
161: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
162: ! By placing code between these two statements, computations can be
163: ! done while messages are in transition.
164: !
165: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
166: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
168: ! Get pointers to vector data
170: PetscCall(VecGetArrayRead(localX, xx, ierr))
171: PetscCall(VecGetArray(F, ff, ierr))
173: ! Get local grid boundaries
175: PetscCall(DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
176: PetscCall(DMDAGetGhostCorners(da, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
178: ! Compute function over the locally owned part of the grid
179: rowf = 0
180: do j = ys, ys + ym - 1
182: row = (j - gys)*gxm + xs - gxs
183: do i = xs, xs + xm - 1
184: row = row + 1
185: rowf = rowf + 1
187: if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
188: ff(rowf) = xx(row)
189: cycle
190: end if
191: u = xx(row)
192: uxx = (two*u - xx(row - 1) - xx(row + 1))*hydhx
193: uyy = (two*u - xx(row - gxm) - xx(row + gxm))*hxdhy
194: ff(rowf) = uxx + uyy - sc*exp(u)
195: end do
196: end do
198: ! Restore vectors
200: PetscCall(VecRestoreArrayRead(localX, xx, ierr))
201: PetscCall(VecRestoreArray(F, ff, ierr))
202: end
204: ! -------------------------------------------------------------------
205: !
206: ! ComputeJacobian - Evaluates Jacobian matrix.
207: !
208: ! Input Parameters:
209: ! x - input vector
210: !
211: ! Output Parameters:
212: ! jac - Jacobian matrix
213: !
214: ! Notes:
215: ! Due to grid point reordering with DMDAs, we must always work
216: ! with the local grid points, and then transform them to the new
217: ! global numbering with the 'ltog' mapping
218: ! We cannot work directly with the global numbers for the original
219: ! uniprocessor grid!
220: !
221: subroutine ComputeJacobian(X, jac, ierr)
223: Vec X
224: Mat jac
225: PetscErrorCode, intent(out) :: ierr
226: PetscInt xs, ys, xm, ym
227: PetscInt gxs, gys, gxm, gym
228: PetscInt grow(1), i, j
229: PetscInt row
230: PetscInt col(5)
231: PetscScalar v(5), hx, hy, hxdhy
232: PetscScalar hydhx, sc
233: ISLocalToGlobalMapping ltogm
234: PetscScalar, pointer ::xx(:)
235: PetscInt, pointer ::ltog(:)
237: hx = one/(mx - 1)
238: hy = one/(my - 1)
239: sc = hx*hy
240: hxdhy = hx/hy
241: hydhx = hy/hx
243: ! Scatter ghost points to local vector, using the 2-step process
244: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
245: ! By placing code between these two statements, computations can be
246: ! done while messages are in transition.
248: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX, ierr))
249: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX, ierr))
251: ! Get pointer to vector data
253: PetscCall(VecGetArrayRead(localX, xx, ierr))
255: ! Get local grid boundaries
257: PetscCall(DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
258: PetscCall(DMDAGetGhostCorners(da, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
260: ! Get the global node numbers for all local nodes, including ghost points
262: PetscCall(DMGetLocalToGlobalMapping(da, ltogm, ierr))
263: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, ltog, ierr))
265: ! Compute entries for the locally owned part of the Jacobian.
266: ! - Currently, all PETSc parallel matrix formats are partitioned by
267: ! contiguous chunks of rows across the processors. The 'grow'
268: ! parameter computed below specifies the global row number
269: ! corresponding to each local grid point.
270: ! - Each processor needs to insert only elements that it owns
271: ! locally (but any non-local elements will be sent to the
272: ! appropriate processor during matrix assembly).
273: ! - Always specify global row and columns of matrix entries.
274: ! - Here, we set all entries for a particular row at once.
276: do j = ys, ys + ym - 1
277: row = (j - gys)*gxm + xs - gxs
278: do i = xs, xs + xm - 1
279: row = row + 1
280: grow(1) = ltog(row)
281: if (i == 0 .or. j == 0 .or. i == (mx - 1) .or. j == (my - 1)) then
282: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, grow, 1_PETSC_INT_KIND, grow, [one], INSERT_VALUES, ierr))
283: cycle
284: end if
285: v(1) = -hxdhy
286: col(1) = ltog(row - gxm)
287: v(2) = -hydhx
288: col(2) = ltog(row - 1)
289: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(xx(row))
290: col(3) = grow(1)
291: v(4) = -hydhx
292: col(4) = ltog(row + 1)
293: v(5) = -hxdhy
294: col(5) = ltog(row + gxm)
295: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, grow, 5_PETSC_INT_KIND, col, v, INSERT_VALUES, ierr))
296: end do
297: end do
299: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, ltog, ierr))
301: ! Assemble matrix, using the 2-step process:
302: ! MatAssemblyBegin(), MatAssemblyEnd().
303: ! By placing code between these two statements, computations can be
304: ! done while messages are in transition.
306: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
307: PetscCall(VecRestoreArrayRead(localX, xx, ierr))
308: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
309: end
310: end module ex14fmodule
312: program main
313: use ex14fmodule
314: implicit none
316: MPIU_Comm comm
317: Vec X, Y, F
318: Mat J
319: KSP ksp
321: PetscInt Nx, Ny, N
322: PetscBool flg, nooutput, usemf
324: ! --------------- Data to define nonlinear solver --------------
325: PetscReal, parameter :: rtol = 1.e-8
326: PetscReal fnorm, ynorm, xnorm, ttol
327: PetscInt, parameter :: max_nonlin_its = 10
328: PetscInt lin_its
329: PetscInt i, m
330: PetscErrorCode ierr
332: PetscCallA(PetscInitialize(ierr))
333: comm = PETSC_COMM_WORLD
335: ! Initialize problem parameters
337: !
338: mx = 4
339: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
340: my = 4
341: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))
342: N = mx*my
344: nooutput = .false.
345: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-no_output', nooutput, ierr))
347: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: ! Create linear solver context
349: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351: PetscCallA(KSPCreate(comm, ksp, ierr))
353: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
354: ! Create vector data structures
355: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357: !
358: ! Create distributed array (DMDA) to manage parallel grid and vectors
359: !
360: Nx = PETSC_DECIDE
361: Ny = PETSC_DECIDE
362: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-Nx', Nx, flg, ierr))
363: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-Ny', Ny, flg, ierr))
364: PetscCallA(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, mx, my, Nx, Ny, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr))
365: PetscCallA(DMSetFromOptions(da, ierr))
366: PetscCallA(DMSetUp(da, ierr))
367: !
368: ! Extract global and local vectors from DMDA then duplicate for remaining
369: ! vectors that are the same types
370: !
371: PetscCallA(DMCreateGlobalVector(da, X, ierr))
372: PetscCallA(DMCreateLocalVector(da, localX, ierr))
373: PetscCallA(VecDuplicate(X, F, ierr))
374: PetscCallA(VecDuplicate(X, Y, ierr))
376: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
377: ! Create matrix data structure for Jacobian
378: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
379: !
380: ! Note: For the parallel case, vectors and matrices MUST be partitioned
381: ! accordingly. When using distributed arrays (DMDAs) to create vectors,
382: ! the DMDAs determine the problem partitioning. We must explicitly
383: ! specify the local matrix dimensions upon its creation for compatibility
384: ! with the vector distribution.
385: !
386: ! Note: Here we only approximately preallocate storage space for the
387: ! Jacobian. See the users manual for a discussion of better techniques
388: ! for preallocating matrix memory.
389: !
390: PetscCallA(VecGetLocalSize(X, m, ierr))
391: PetscCallA(MatCreateFromOptions(comm, PETSC_NULL_CHARACTER, 1_PETSC_INT_KIND, m, m, N, N, B, ierr))
393: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
394: ! if usemf is on then matrix-vector product is done via matrix-free
395: ! approach. Note this is just an example, and not realistic because
396: ! we still use the actual formed matrix, but in reality one would
397: ! provide their own subroutine that would directly do the matrix
398: ! vector product and call MatMult()
399: ! Note: we put B into a module so it will be visible to the
400: ! mymult() routine
401: usemf = .false.
402: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mf', usemf, ierr))
403: if (usemf) then
404: PetscCallA(MatCreateShell(comm, m, m, N, N, PETSC_NULL_INTEGER, J, ierr))
405: PetscCallA(MatShellSetOperation(J, MATOP_MULT, mymult, ierr))
406: else
407: ! If not doing matrix-free then matrix operator, J, and matrix used
408: ! to construct preconditioner, B, are the same
409: J = B
410: end if
412: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413: ! Customize linear solver set runtime options
414: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
415: !
416: ! Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
417: !
418: PetscCallA(KSPSetFromOptions(ksp, ierr))
420: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
421: ! Evaluate initial guess
422: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
424: PetscCallA(FormInitialGuess(X, ierr))
425: PetscCallA(ComputeFunction(X, F, ierr))
426: PetscCallA(VecNorm(F, NORM_2, fnorm, ierr))
427: ttol = fnorm*rtol
428: if (.not. nooutput) print *, 'Initial function norm ', fnorm
430: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
431: ! Solve nonlinear system with a user-defined method
432: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434: ! This solver is a very simplistic inexact Newton method, with no
435: ! no damping strategies or bells and whistles. The intent of this code
436: ! is merely to demonstrate the repeated solution with KSP of linear
437: ! systems with the same nonzero structure.
438: !
439: ! This is NOT the recommended approach for solving nonlinear problems
440: ! with PETSc! We urge users to employ the SNES component for solving
441: ! nonlinear problems whenever possible with application codes, as it
442: ! offers many advantages over coding nonlinear solvers independently.
444: do i = 0, max_nonlin_its
446: ! Compute the Jacobian matrix. See the comments in this routine for
447: ! important information about setting the flag mat_flag.
449: PetscCallA(ComputeJacobian(X, B, ierr))
451: ! Solve J Y = F, where J is the Jacobian matrix.
452: ! - First, set the KSP linear operators. Here the matrix that
453: ! defines the linear system also serves as the preconditioning
454: ! matrix.
455: ! - Then solve the Newton system.
457: PetscCallA(KSPSetOperators(ksp, J, B, ierr))
458: PetscCallA(KSPSolve(ksp, F, Y, ierr))
460: ! Compute updated iterate
462: PetscCallA(VecNorm(Y, NORM_2, ynorm, ierr))
463: PetscCallA(VecAYPX(Y, mone, X, ierr))
464: PetscCallA(VecCopy(Y, X, ierr))
465: PetscCallA(VecNorm(X, NORM_2, xnorm, ierr))
466: PetscCallA(KSPGetIterationNumber(ksp, lin_its, ierr))
467: if (.not. nooutput) print *, 'linear solve iterations = ', lin_its, ' xnorm = ', xnorm, ' ynorm = ', ynorm
469: ! Evaluate nonlinear function at new location
471: PetscCallA(ComputeFunction(X, F, ierr))
472: PetscCallA(VecNorm(F, NORM_2, fnorm, ierr))
473: if (.not. nooutput) print *, 'Iteration ', i + 1, ' function norm', fnorm
475: ! Test for convergence
477: if (fnorm <= ttol) then
478: if (.not. nooutput) print *, 'Converged: function norm ', fnorm, ' tolerance ', ttol
479: exit
480: end if
481: end do
483: write (6, 100) i + 1
484: 100 format('Number of SNES iterations =', I2)
486: ! Check if mymult() produces a linear operator
487: if (usemf) then
488: N = 5
489: PetscCallA(MatIsLinear(J, N, flg, ierr))
490: if (.not. flg) print *, 'IsLinear', flg
491: end if
493: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
494: ! Free work space. All PETSc objects should be destroyed when they
495: ! are no longer needed.
496: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
498: PetscCallA(MatDestroy(B, ierr))
499: if (usemf) then
500: PetscCallA(MatDestroy(J, ierr))
501: end if
502: PetscCallA(VecDestroy(localX, ierr))
503: PetscCallA(VecDestroy(X, ierr))
504: PetscCallA(VecDestroy(Y, ierr))
505: PetscCallA(VecDestroy(F, ierr))
506: PetscCallA(KSPDestroy(ksp, ierr))
507: PetscCallA(DMDestroy(da, ierr))
508: PetscCallA(PetscFinalize(ierr))
509: end
511: !/*TEST
512: !
513: ! test:
514: ! args: -no_output -ksp_gmres_cgs_refinement_type refine_always
515: ! output_file: output/ex14_1.out
516: ! requires: !single
517: !
518: !TEST*/