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:       module mymoduleex14f
 44: #include <petsc/finclude/petscksp.h>
 45:       use petscdmda
 46:       use petscksp
 47:       Vec      localX
 48:       PetscInt mx,my
 49:       Mat B
 50:       DM da
 51:       end module

 53:       program main
 54:       use mymoduleex14f
 55:       implicit none

 57:       MPI_Comm comm
 58:       Vec      X,Y,F
 59:       Mat      J
 60:       KSP      ksp

 62:       PetscInt  Nx,Ny,N,ifive,ithree
 63:       PetscBool  flg,nooutput,usemf
 64: !
 65: !      This is the routine to use for matrix-free approach
 66: !
 67:       external mymult

 69: !     --------------- Data to define nonlinear solver --------------
 70:       PetscReal   rtol,ttol
 71:       PetscReal   fnorm,ynorm,xnorm
 72:       PetscInt            max_nonlin_its,one
 73:       PetscInt            lin_its
 74:       PetscInt           i,m
 75:       PetscScalar        mone
 76:       PetscErrorCode ierr

 78:       mone           = -1.0
 79:       rtol           = 1.e-8
 80:       max_nonlin_its = 10
 81:       one            = 1
 82:       ifive          = 5
 83:       ithree         = 3

 85:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 86:       if (ierr .ne. 0) then
 87:         print*,'Unable to initialize PETSc'
 88:         stop
 89:       endif
 90:       comm = PETSC_COMM_WORLD

 92: !  Initialize problem parameters

 94: !
 95:       mx = 4
 96:       my = 4
 97:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
 98:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
 99:       N = mx*my

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

104: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: !     Create linear solver context
106: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

108:       call KSPCreate(comm,ksp,ierr)

110: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: !     Create vector data structures
112: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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

170: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: !     Customize linear solver set runtime options
172: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: !
174: !     Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
175: !
176:        call KSPSetFromOptions(ksp,ierr)

178: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: !     Evaluate initial guess
180: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

182:        call FormInitialGuess(X,ierr)
183:        call ComputeFunction(X,F,ierr)
184:        call VecNorm(F,NORM_2,fnorm,ierr)
185:        ttol = fnorm*rtol
186:        if (.not. nooutput) then
187:          print*, 'Initial function norm ',fnorm
188:        endif

190: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191: !     Solve nonlinear system with a user-defined method
192: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

204:        do 10 i=0,max_nonlin_its

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

209:          call ComputeJacobian(X,B,ierr)

211: !  Solve J Y = F, where J is the Jacobian matrix.
212: !    - First, set the KSP linear operators.  Here the matrix that
213: !      defines the linear system also serves as the preconditioning
214: !      matrix.
215: !    - Then solve the Newton system.

217:          call KSPSetOperators(ksp,J,B,ierr)
218:          call KSPSolve(ksp,F,Y,ierr)

220: !  Compute updated iterate

222:          call VecNorm(Y,NORM_2,ynorm,ierr)
223:          call VecAYPX(Y,mone,X,ierr)
224:          call VecCopy(Y,X,ierr)
225:          call VecNorm(X,NORM_2,xnorm,ierr)
226:          call KSPGetIterationNumber(ksp,lin_its,ierr)
227:          if (.not. nooutput) then
228:            print*,'linear solve iterations = ',lin_its,' xnorm = ',xnorm,' ynorm = ',ynorm
229:          endif

231: !  Evaluate nonlinear function at new location

233:          call ComputeFunction(X,F,ierr)
234:          call VecNorm(F,NORM_2,fnorm,ierr)
235:          if (.not. nooutput) then
236:            print*, 'Iteration ',i+1,' function norm',fnorm
237:          endif

239: !  Test for convergence

241:        if (fnorm .le. ttol) then
242:          if (.not. nooutput) then
243:            print*,'Converged: function norm ',fnorm,' tolerance ',ttol
244:          endif
245:          goto 20
246:        endif
247:  10   continue
248:  20   continue

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

253: !     Check if mymult() produces a linear operator
254:       if (usemf) then
255:          N = 5
256:          call MatIsLinear(J,N,flg,ierr)
257:          if (.not. flg) then
258:             print *, 'IsLinear',flg
259:          endif
260:       endif

262: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
263: !     Free work space.  All PETSc objects should be destroyed when they
264: !     are no longer needed.
265: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

267:        call MatDestroy(B,ierr)
268:        if (usemf) then
269:          call MatDestroy(J,ierr)
270:        endif
271:        call VecDestroy(localX,ierr)
272:        call VecDestroy(X,ierr)
273:        call VecDestroy(Y,ierr)
274:        call VecDestroy(F,ierr)
275:        call KSPDestroy(ksp,ierr)
276:        call DMDestroy(da,ierr)
277:        call PetscFinalize(ierr)
278:        end

280: ! -------------------------------------------------------------------
281: !
282: !   FormInitialGuess - Forms initial approximation.
283: !
284: !   Input Parameters:
285: !   X - vector
286: !
287: !   Output Parameter:
288: !   X - vector
289: !
290:       subroutine FormInitialGuess(X,ierr)
291:       use mymoduleex14f
292:       implicit none

294:       PetscErrorCode    ierr
295:       PetscOffset      idx
296:       Vec       X
297:       PetscInt  i,j,row
298:       PetscInt  xs,ys,xm
299:       PetscInt  ym
300:       PetscReal one,lambda,temp1,temp,hx,hy
301:       PetscScalar      xx(2)

303:       one    = 1.0
304:       lambda = 6.0
305:       hx     = one/(mx-1)
306:       hy     = one/(my-1)
307:       temp1  = lambda/(lambda + one)

309: !  Get a pointer to vector data.
310: !    - VecGetArray() returns a pointer to the data array.
311: !    - You MUST call VecRestoreArray() when you no longer need access to
312: !      the array.
313:        call VecGetArray(X,xx,idx,ierr)

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

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

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

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

335: !     Restore vector

337:        call VecRestoreArray(X,xx,idx,ierr)
338:        return
339:        end

341: ! -------------------------------------------------------------------
342: !
343: !   ComputeFunction - Evaluates nonlinear function, F(x).
344: !
345: !   Input Parameters:
346: !.  X - input vector
347: !
348: !   Output Parameter:
349: !.  F - function vector
350: !
351:       subroutine  ComputeFunction(X,F,ierr)
352:       use mymoduleex14f
353:       implicit none

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

365:       two    = 2.0
366:       one    = 1.0
367:       lambda = 6.0

369:       hx     = one/(mx-1)
370:       hy     = one/(my-1)
371:       sc     = hx*hy*lambda
372:       hxdhy  = hx/hy
373:       hydhx  = hy/hx

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

383: !  Get pointers to vector data

385:       call VecGetArray(localX,xx,idx,ierr)
386:       call VecGetArray(F,ff,idf,ierr)

388: !  Get local grid boundaries

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

393: !  Compute function over the locally owned part of the grid
394:       rowf = 0
395:       do 50 j=ys,ys+ym-1

397:         row  = (j - gys)*gxm + xs - gxs
398:         do 60 i=xs,xs+xm-1
399:           row  = row + 1
400:           rowf = rowf + 1

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

413: !  Restore vectors

415:        call VecRestoreArray(localX,xx,idx,ierr)
416:        call VecRestoreArray(F,ff,idf,ierr)
417:        return
418:        end

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

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

457:       ione   = 1
458:       ifive  = 5
459:       one    = 1.0
460:       two    = 2.0
461:       hx     = one/(mx-1)
462:       hy     = one/(my-1)
463:       sc     = hx*hy
464:       hxdhy  = hx/hy
465:       hydhx  = hy/hx
466:       lambda = 6.0

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

473:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
474:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)

476: !  Get pointer to vector data

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

480: !  Get local grid boundaries

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

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

487:       call DMGetLocalToGlobalMapping(da,ltogm,ierr)
488:       call ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)

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

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

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

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

531:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
532:       call VecRestoreArray(localX,xx,idx,ierr)
533:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
534:       return
535:       end

537: ! -------------------------------------------------------------------
538: !
539: !   MyMult - user provided matrix multiply
540: !
541: !   Input Parameters:
542: !.  X - input vector
543: !
544: !   Output Parameter:
545: !.  F - function vector
546: !
547:       subroutine  MyMult(J,X,F,ierr)
548:       use mymoduleex14f
549:       implicit none

551:       Mat     J
552:       Vec     X,F
553:       PetscErrorCode ierr
554: !
555: !       Here we use the actual formed matrix B; users would
556: !     instead write their own matrix vector product routine
557: !
558:       call MatMult(B,X,F,ierr)
559:       return
560:       end

562: !/*TEST
563: !
564: !   test:
565: !      args: -no_output -ksp_gmres_cgs_refinement_type refine_always
566: !      output_file: output/ex14_1.out
567: !      requires: !single
568: !
569: !TEST*/