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: !
 24: !!/*T
 25: !   Concepts: KSP^writing a user-defined nonlinear solver
 26: !   Concepts: DMDA^using distributed arrays
 27: !   Processors: n
 28: !T*/

 30: !  ------------------------------------------------------------------------
 31: !
 32: !    Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 33: !    the partial differential equation
 34: !
 35: !            -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 36: !
 37: !    with boundary conditions
 38: !
 39: !             u = 0  for  x = 0, x = 1, y = 0, y = 1.
 40: !
 41: !    A finite difference approximation with the usual 5-point stencil
 42: !    is used to discretize the boundary value problem to obtain a nonlinear
 43: !    system of equations.
 44: !
 45: !    The SNES version of this problem is:  snes/tutorials/ex5f.F
 46: !
 47: !  -------------------------------------------------------------------------
 48:       module mymoduleex14f
 49: #include <petsc/finclude/petscksp.h>
 50:       use petscdmda
 51:       use petscksp
 52:       Vec      localX
 53:       PetscInt mx,my
 54:       Mat B
 55:       DM da
 56:       end module

 58:       program main
 59:       use mymoduleex14f
 60:       implicit none

 62:       MPI_Comm comm
 63:       Vec      X,Y,F
 64:       Mat      J
 65:       KSP      ksp

 67:       PetscInt  Nx,Ny,N,ifive,ithree
 68:       PetscBool  flg,nooutput,usemf
 69: !
 70: !      This is the routine to use for matrix-free approach
 71: !
 72:       external mymult

 74: !     --------------- Data to define nonlinear solver --------------
 75:       PetscReal   rtol,ttol
 76:       PetscReal   fnorm,ynorm,xnorm
 77:       PetscInt            max_nonlin_its,one
 78:       PetscInt            lin_its
 79:       PetscInt           i,m
 80:       PetscScalar        mone
 81:       PetscErrorCode ierr

 83:       mone           = -1.0
 84:       rtol           = 1.e-8
 85:       max_nonlin_its = 10
 86:       one            = 1
 87:       ifive          = 5
 88:       ithree         = 3

 90:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 91:       if (ierr .ne. 0) then
 92:         print*,'Unable to initialize PETSc'
 93:         stop
 94:       endif
 95:       comm = PETSC_COMM_WORLD

 97: !  Initialize problem parameters

 99: !
100:       mx = 4
101:       my = 4
102:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
103:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
104:       N = mx*my

106:       nooutput = .false.
107:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_output',nooutput,ierr)

109: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: !     Create linear solver context
111: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

113:       call KSPCreate(comm,ksp,ierr)

115: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: !     Create vector data structures
117: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

119: !
120: !  Create distributed array (DMDA) to manage parallel grid and vectors
121: !
122:       Nx = PETSC_DECIDE
123:       Ny = PETSC_DECIDE
124:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-Nx',Nx,flg,ierr)
125:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-Ny',Ny,flg,ierr)
126:       call DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,mx,my,Nx,Ny,one,one,          &
127:      &                  PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
128:       call DMSetFromOptions(da,ierr)
129:       call DMSetUp(da,ierr)
130: !
131: !  Extract global and local vectors from DMDA then duplicate for remaining
132: !  vectors that are the same types
133: !
134:        call DMCreateGlobalVector(da,X,ierr)
135:        call DMCreateLocalVector(da,localX,ierr)
136:        call VecDuplicate(X,F,ierr)
137:        call VecDuplicate(X,Y,ierr)

139: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: !     Create matrix data structure for Jacobian
141: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: !
143: !     Note:  For the parallel case, vectors and matrices MUST be partitioned
144: !     accordingly.  When using distributed arrays (DMDAs) to create vectors,
145: !     the DMDAs determine the problem partitioning.  We must explicitly
146: !     specify the local matrix dimensions upon its creation for compatibility
147: !     with the vector distribution.
148: !
149: !     Note: Here we only approximately preallocate storage space for the
150: !     Jacobian.  See the users manual for a discussion of better techniques
151: !     for preallocating matrix memory.
152: !
153:       call VecGetLocalSize(X,m,ierr)
154:       call MatCreateAIJ(comm,m,m,N,N,ifive,PETSC_NULL_INTEGER,ithree,PETSC_NULL_INTEGER,B,ierr)

156: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: !     if usemf is on then matrix vector product is done via matrix free
158: !     approach. Note this is just an example, and not realistic because
159: !     we still use the actual formed matrix, but in reality one would
160: !     provide their own subroutine that would directly do the matrix
161: !     vector product and not call MatMult()
162: !     Note: we put B into a module so it will be visible to the
163: !     mymult() routine
164:       usemf = .false.
165:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mf',usemf,ierr)
166:       if (usemf) then
167:          call MatCreateShell(comm,m,m,N,N,PETSC_NULL_INTEGER,J,ierr)
168:          call MatShellSetOperation(J,MATOP_MULT,mymult,ierr)
169:       else
170: !        If not doing matrix free then matrix operator, J,  and matrix used
171: !        to construct preconditioner, B, are the same
172:         J = B
173:       endif

175: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: !     Customize linear solver set runtime options
177: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: !
179: !     Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
180: !
181:        call KSPSetFromOptions(ksp,ierr)

183: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: !     Evaluate initial guess
185: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

187:        call FormInitialGuess(X,ierr)
188:        call ComputeFunction(X,F,ierr)
189:        call VecNorm(F,NORM_2,fnorm,ierr)
190:        ttol = fnorm*rtol
191:        if (.not. nooutput) then
192:          print*, 'Initial function norm ',fnorm
193:        endif

195: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: !     Solve nonlinear system with a user-defined method
197: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

199: !  This solver is a very simplistic inexact Newton method, with no
200: !  no damping strategies or bells and whistles. The intent of this code
201: !  is merely to demonstrate the repeated solution with KSP of linear
202: !  systems with the same nonzero structure.
203: !
204: !  This is NOT the recommended approach for solving nonlinear problems
205: !  with PETSc!  We urge users to employ the SNES component for solving
206: !  nonlinear problems whenever possible with application codes, as it
207: !  offers many advantages over coding nonlinear solvers independently.

209:        do 10 i=0,max_nonlin_its

211: !  Compute the Jacobian matrix.  See the comments in this routine for
212: !  important information about setting the flag mat_flag.

214:          call ComputeJacobian(X,B,ierr)

216: !  Solve J Y = F, where J is the Jacobian matrix.
217: !    - First, set the KSP linear operators.  Here the matrix that
218: !      defines the linear system also serves as the preconditioning
219: !      matrix.
220: !    - Then solve the Newton system.

222:          call KSPSetOperators(ksp,J,B,ierr)
223:          call KSPSolve(ksp,F,Y,ierr)

225: !  Compute updated iterate

227:          call VecNorm(Y,NORM_2,ynorm,ierr)
228:          call VecAYPX(Y,mone,X,ierr)
229:          call VecCopy(Y,X,ierr)
230:          call VecNorm(X,NORM_2,xnorm,ierr)
231:          call KSPGetIterationNumber(ksp,lin_its,ierr)
232:          if (.not. nooutput) then
233:            print*,'linear solve iterations = ',lin_its,' xnorm = ',xnorm,' ynorm = ',ynorm
234:          endif

236: !  Evaluate nonlinear function at new location

238:          call ComputeFunction(X,F,ierr)
239:          call VecNorm(F,NORM_2,fnorm,ierr)
240:          if (.not. nooutput) then
241:            print*, 'Iteration ',i+1,' function norm',fnorm
242:          endif

244: !  Test for convergence

246:        if (fnorm .le. ttol) then
247:          if (.not. nooutput) then
248:            print*,'Converged: function norm ',fnorm,' tolerance ',ttol
249:          endif
250:          goto 20
251:        endif
252:  10   continue
253:  20   continue

255:       write(6,100) i+1
256:  100  format('Number of SNES iterations =',I2)

258: !     Check if mymult() produces a linear operator
259:       if (usemf) then
260:          N = 5
261:          call MatIsLinear(J,N,flg,ierr)
262:          if (.not. flg) then
263:             print *, 'IsLinear',flg
264:          endif
265:       endif

267: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268: !     Free work space.  All PETSc objects should be destroyed when they
269: !     are no longer needed.
270: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

272:        call MatDestroy(B,ierr)
273:        if (usemf) then
274:          call MatDestroy(J,ierr)
275:        endif
276:        call VecDestroy(localX,ierr)
277:        call VecDestroy(X,ierr)
278:        call VecDestroy(Y,ierr)
279:        call VecDestroy(F,ierr)
280:        call KSPDestroy(ksp,ierr)
281:        call DMDestroy(da,ierr)
282:        call PetscFinalize(ierr)
283:        end

285: ! -------------------------------------------------------------------
286: !
287: !   FormInitialGuess - Forms initial approximation.
288: !
289: !   Input Parameters:
290: !   X - vector
291: !
292: !   Output Parameter:
293: !   X - vector
294: !
295:       subroutine FormInitialGuess(X,ierr)
296:       use mymoduleex14f
297:       implicit none

299:       PetscErrorCode    ierr
300:       PetscOffset      idx
301:       Vec       X
302:       PetscInt  i,j,row
303:       PetscInt  xs,ys,xm
304:       PetscInt  ym
305:       PetscReal one,lambda,temp1,temp,hx,hy
306:       PetscScalar      xx(2)

308:       one    = 1.0
309:       lambda = 6.0
310:       hx     = one/(mx-1)
311:       hy     = one/(my-1)
312:       temp1  = lambda/(lambda + one)

314: !  Get a pointer to vector data.
315: !    - VecGetArray() returns a pointer to the data array.
316: !    - You MUST call VecRestoreArray() when you no longer need access to
317: !      the array.
318:        call VecGetArray(X,xx,idx,ierr)

320: !  Get local grid boundaries (for 2-dimensional DMDA):
321: !    xs, ys   - starting grid indices (no ghost points)
322: !    xm, ym   - widths of local grid (no ghost points)

324:        call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)

326: !  Compute initial guess over the locally owned part of the grid

328:       do 30 j=ys,ys+ym-1
329:         temp = (min(j,my-j-1))*hy
330:         do 40 i=xs,xs+xm-1
331:           row = i - xs + (j - ys)*xm + 1
332:           if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. j .eq. my-1) then
333:             xx(idx+row) = 0.0
334:             continue
335:           endif
336:           xx(idx+row) = temp1*sqrt(min((min(i,mx-i-1))*hx,temp))
337:  40     continue
338:  30   continue

340: !     Restore vector

342:        call VecRestoreArray(X,xx,idx,ierr)
343:        return
344:        end

346: ! -------------------------------------------------------------------
347: !
348: !   ComputeFunction - Evaluates nonlinear function, F(x).
349: !
350: !   Input Parameters:
351: !.  X - input vector
352: !
353: !   Output Parameter:
354: !.  F - function vector
355: !
356:       subroutine  ComputeFunction(X,F,ierr)
357:       use mymoduleex14f
358:       implicit none

360:       Vec              X,F
361:       PetscInt         gys,gxm,gym
362:       PetscOffset      idx,idf
363:       PetscErrorCode ierr
364:       PetscInt i,j,row,xs,ys,xm,ym,gxs
365:       PetscInt rowf
366:       PetscReal two,one,lambda,hx
367:       PetscReal hy,hxdhy,hydhx,sc
368:       PetscScalar      u,uxx,uyy,xx(2),ff(2)

370:       two    = 2.0
371:       one    = 1.0
372:       lambda = 6.0

374:       hx     = one/(mx-1)
375:       hy     = one/(my-1)
376:       sc     = hx*hy*lambda
377:       hxdhy  = hx/hy
378:       hydhx  = hy/hx

380: !  Scatter ghost points to local vector, using the 2-step process
381: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
382: !  By placing code between these two statements, computations can be
383: !  done while messages are in transition.
384: !
385:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
386:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)

388: !  Get pointers to vector data

390:       call VecGetArray(localX,xx,idx,ierr)
391:       call VecGetArray(F,ff,idf,ierr)

393: !  Get local grid boundaries

395:       call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
396:       call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)

398: !  Compute function over the locally owned part of the grid
399:       rowf = 0
400:       do 50 j=ys,ys+ym-1

402:         row  = (j - gys)*gxm + xs - gxs
403:         do 60 i=xs,xs+xm-1
404:           row  = row + 1
405:           rowf = rowf + 1

407:           if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. j .eq. my-1) then
408:             ff(idf+rowf) = xx(idx+row)
409:             goto 60
410:           endif
411:           u   = xx(idx+row)
412:           uxx = (two*u - xx(idx+row-1) - xx(idx+row+1))*hydhx
413:           uyy = (two*u - xx(idx+row-gxm) - xx(idx+row+gxm))*hxdhy
414:           ff(idf+rowf) = uxx + uyy - sc*exp(u)
415:  60     continue
416:  50   continue

418: !  Restore vectors

420:        call VecRestoreArray(localX,xx,idx,ierr)
421:        call VecRestoreArray(F,ff,idf,ierr)
422:        return
423:        end

425: ! -------------------------------------------------------------------
426: !
427: !   ComputeJacobian - Evaluates Jacobian matrix.
428: !
429: !   Input Parameters:
430: !   x - input vector
431: !
432: !   Output Parameters:
433: !   jac - Jacobian matrix
434: !   flag - flag indicating matrix structure
435: !
436: !   Notes:
437: !   Due to grid point reordering with DMDAs, we must always work
438: !   with the local grid points, and then transform them to the new
439: !   global numbering with the 'ltog' mapping
440: !   We cannot work directly with the global numbers for the original
441: !   uniprocessor grid!
442: !
443:       subroutine ComputeJacobian(X,jac,ierr)
444:       use mymoduleex14f
445:       implicit none

447:       Vec         X
448:       Mat         jac
449:       PetscInt     ltog(2)
450:       PetscOffset idltog,idx
451:       PetscErrorCode ierr
452:       PetscInt xs,ys,xm,ym
453:       PetscInt gxs,gys,gxm,gym
454:       PetscInt grow(1),i,j
455:       PetscInt row,ione
456:       PetscInt col(5),ifive
457:       PetscScalar two,one,lambda
458:       PetscScalar v(5),hx,hy,hxdhy
459:       PetscScalar hydhx,sc,xx(2)
460:       ISLocalToGlobalMapping ltogm

462:       ione   = 1
463:       ifive  = 5
464:       one    = 1.0
465:       two    = 2.0
466:       hx     = one/(mx-1)
467:       hy     = one/(my-1)
468:       sc     = hx*hy
469:       hxdhy  = hx/hy
470:       hydhx  = hy/hx
471:       lambda = 6.0

473: !  Scatter ghost points to local vector, using the 2-step process
474: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
475: !  By placing code between these two statements, computations can be
476: !  done while messages are in transition.

478:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
479:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)

481: !  Get pointer to vector data

483:       call VecGetArray(localX,xx,idx,ierr)

485: !  Get local grid boundaries

487:       call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
488:       call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)

490: !  Get the global node numbers for all local nodes, including ghost points

492:       call DMGetLocalToGlobalMapping(da,ltogm,ierr)
493:       call ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)

495: !  Compute entries for the locally owned part of the Jacobian.
496: !   - Currently, all PETSc parallel matrix formats are partitioned by
497: !     contiguous chunks of rows across the processors. The 'grow'
498: !     parameter computed below specifies the global row number
499: !     corresponding to each local grid point.
500: !   - Each processor needs to insert only elements that it owns
501: !     locally (but any non-local elements will be sent to the
502: !     appropriate processor during matrix assembly).
503: !   - Always specify global row and columns of matrix entries.
504: !   - Here, we set all entries for a particular row at once.

506:       do 10 j=ys,ys+ym-1
507:         row = (j - gys)*gxm + xs - gxs
508:         do 20 i=xs,xs+xm-1
509:           row = row + 1
510:           grow(1) = ltog(idltog+row)
511:           if (i .eq. 0 .or. j .eq. 0 .or. i .eq. (mx-1) .or. j .eq. (my-1)) then
512:              call MatSetValues(jac,ione,grow,ione,grow,one,INSERT_VALUES,ierr)
513:              go to 20
514:           endif
515:           v(1)   = -hxdhy
516:           col(1) = ltog(idltog+row - gxm)
517:           v(2)   = -hydhx
518:           col(2) = ltog(idltog+row - 1)
519:           v(3)   = two*(hydhx + hxdhy) - sc*lambda*exp(xx(idx+row))
520:           col(3) = grow(1)
521:           v(4)   = -hydhx
522:           col(4) = ltog(idltog+row + 1)
523:           v(5)   = -hxdhy
524:           col(5) = ltog(idltog+row + gxm)
525:           call MatSetValues(jac,ione,grow,ifive,col,v,INSERT_VALUES,ierr)
526:  20     continue
527:  10   continue

529:       call ISLocalToGlobalMappingRestoreIndices(ltogm,ltog,idltog,ierr)

531: !  Assemble matrix, using the 2-step process:
532: !    MatAssemblyBegin(), MatAssemblyEnd().
533: !  By placing code between these two statements, computations can be
534: !  done while messages are in transition.

536:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
537:       call VecRestoreArray(localX,xx,idx,ierr)
538:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
539:       return
540:       end

542: ! -------------------------------------------------------------------
543: !
544: !   MyMult - user provided matrix multiply
545: !
546: !   Input Parameters:
547: !.  X - input vector
548: !
549: !   Output Parameter:
550: !.  F - function vector
551: !
552:       subroutine  MyMult(J,X,F,ierr)
553:       use mymoduleex14f
554:       implicit none

556:       Mat     J
557:       Vec     X,F
558:       PetscErrorCode ierr
559: !
560: !       Here we use the actual formed matrix B; users would
561: !     instead write their own matrix vector product routine
562: !
563:       call MatMult(B,X,F,ierr)
564:       return
565:       end

567: !/*TEST
568: !
569: !   test:
570: !      args: -no_output -ksp_gmres_cgs_refinement_type refine_always
571: !      output_file: output/ex14_1.out
572: !      requires: !single
573: !
574: !TEST*/