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 ex14fmodule
 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 ex14fmodule
 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:       PetscCallA(PetscInitialize(ierr))
 86:       comm = PETSC_COMM_WORLD

 88: !  Initialize problem parameters

 90: !
 91:       mx = 4
 92:       my = 4
 93:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
 94:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
 95:       N = mx*my

 97:       nooutput = .false.
 98:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-no_output',nooutput,ierr))

100: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: !     Create linear solver context
102: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

104:       PetscCallA(KSPCreate(comm,ksp,ierr))

106: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: !     Create vector data structures
108: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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

165: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: !     Customize linear solver set runtime options
167: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: !
169: !     Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
170: !
171:        PetscCallA(KSPSetFromOptions(ksp,ierr))

173: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: !     Evaluate initial guess
175: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

177:        PetscCallA(FormInitialGuess(X,ierr))
178:        PetscCallA(ComputeFunction(X,F,ierr))
179:        PetscCallA(VecNorm(F,NORM_2,fnorm,ierr))
180:        ttol = fnorm*rtol
181:        if (.not. nooutput) then
182:          print*, 'Initial function norm ',fnorm
183:        endif

185: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: !     Solve nonlinear system with a user-defined method
187: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

199:        do 10 i=0,max_nonlin_its

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

204:          PetscCallA(ComputeJacobian(X,B,ierr))

206: !  Solve J Y = F, where J is the Jacobian matrix.
207: !    - First, set the KSP linear operators.  Here the matrix that
208: !      defines the linear system also serves as the preconditioning
209: !      matrix.
210: !    - Then solve the Newton system.

212:          PetscCallA(KSPSetOperators(ksp,J,B,ierr))
213:          PetscCallA(KSPSolve(ksp,F,Y,ierr))

215: !  Compute updated iterate

217:          PetscCallA(VecNorm(Y,NORM_2,ynorm,ierr))
218:          PetscCallA(VecAYPX(Y,mone,X,ierr))
219:          PetscCallA(VecCopy(Y,X,ierr))
220:          PetscCallA(VecNorm(X,NORM_2,xnorm,ierr))
221:          PetscCallA(KSPGetIterationNumber(ksp,lin_its,ierr))
222:          if (.not. nooutput) then
223:            print*,'linear solve iterations = ',lin_its,' xnorm = ',xnorm,' ynorm = ',ynorm
224:          endif

226: !  Evaluate nonlinear function at new location

228:          PetscCallA(ComputeFunction(X,F,ierr))
229:          PetscCallA(VecNorm(F,NORM_2,fnorm,ierr))
230:          if (.not. nooutput) then
231:            print*, 'Iteration ',i+1,' function norm',fnorm
232:          endif

234: !  Test for convergence

236:        if (fnorm .le. ttol) then
237:          if (.not. nooutput) then
238:            print*,'Converged: function norm ',fnorm,' tolerance ',ttol
239:          endif
240:          goto 20
241:        endif
242:  10   continue
243:  20   continue

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

248: !     Check if mymult() produces a linear operator
249:       if (usemf) then
250:          N = 5
251:          PetscCallA(MatIsLinear(J,N,flg,ierr))
252:          if (.not. flg) then
253:             print *, 'IsLinear',flg
254:          endif
255:       endif

257: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: !     Free work space.  All PETSc objects should be destroyed when they
259: !     are no longer needed.
260: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

262:        PetscCallA(MatDestroy(B,ierr))
263:        if (usemf) then
264:          PetscCallA(MatDestroy(J,ierr))
265:        endif
266:        PetscCallA(VecDestroy(localX,ierr))
267:        PetscCallA(VecDestroy(X,ierr))
268:        PetscCallA(VecDestroy(Y,ierr))
269:        PetscCallA(VecDestroy(F,ierr))
270:        PetscCallA(KSPDestroy(ksp,ierr))
271:        PetscCallA(DMDestroy(da,ierr))
272:        PetscCallA(PetscFinalize(ierr))
273:        end

275: ! -------------------------------------------------------------------
276: !
277: !   FormInitialGuess - Forms initial approximation.
278: !
279: !   Input Parameters:
280: !   X - vector
281: !
282: !   Output Parameter:
283: !   X - vector
284: !
285:       subroutine FormInitialGuess(X,ierr)
286:       use ex14fmodule
287:       implicit none

289:       PetscErrorCode    ierr
290:       PetscOffset      idx
291:       Vec       X
292:       PetscInt  i,j,row
293:       PetscInt  xs,ys,xm
294:       PetscInt  ym
295:       PetscReal one,lambda,temp1,temp,hx,hy
296:       PetscScalar      xx(2)

298:       one    = 1.0
299:       lambda = 6.0
300:       hx     = one/(mx-1)
301:       hy     = one/(my-1)
302:       temp1  = lambda/(lambda + one)

304: !  Get a pointer to vector data.
305: !    - VecGetArray() returns a pointer to the data array.
306: !    - You MUST call VecRestoreArray() when you no longer need access to
307: !      the array.
308:        PetscCall(VecGetArray(X,xx,idx,ierr))

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

314:        PetscCall(DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))

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

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

330: !     Restore vector

332:        PetscCall(VecRestoreArray(X,xx,idx,ierr))
333:        return
334:        end

336: ! -------------------------------------------------------------------
337: !
338: !   ComputeFunction - Evaluates nonlinear function, F(x).
339: !
340: !   Input Parameters:
341: !.  X - input vector
342: !
343: !   Output Parameter:
344: !.  F - function vector
345: !
346:       subroutine  ComputeFunction(X,F,ierr)
347:       use ex14fmodule
348:       implicit none

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

360:       two    = 2.0
361:       one    = 1.0
362:       lambda = 6.0

364:       hx     = one/(mx-1)
365:       hy     = one/(my-1)
366:       sc     = hx*hy*lambda
367:       hxdhy  = hx/hy
368:       hydhx  = hy/hx

370: !  Scatter ghost points to local vector, using the 2-step process
371: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
372: !  By placing code between these two statements, computations can be
373: !  done while messages are in transition.
374: !
375:       PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr))
376:       PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr))

378: !  Get pointers to vector data

380:       PetscCall(VecGetArray(localX,xx,idx,ierr))
381:       PetscCall(VecGetArray(F,ff,idf,ierr))

383: !  Get local grid boundaries

385:       PetscCall(DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
386:       PetscCall(DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

388: !  Compute function over the locally owned part of the grid
389:       rowf = 0
390:       do 50 j=ys,ys+ym-1

392:         row  = (j - gys)*gxm + xs - gxs
393:         do 60 i=xs,xs+xm-1
394:           row  = row + 1
395:           rowf = rowf + 1

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

408: !  Restore vectors

410:        PetscCall(VecRestoreArray(localX,xx,idx,ierr))
411:        PetscCall(VecRestoreArray(F,ff,idf,ierr))
412:        return
413:        end

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

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

452:       ione   = 1
453:       ifive  = 5
454:       one    = 1.0
455:       two    = 2.0
456:       hx     = one/(mx-1)
457:       hy     = one/(my-1)
458:       sc     = hx*hy
459:       hxdhy  = hx/hy
460:       hydhx  = hy/hx
461:       lambda = 6.0

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

468:       PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr))
469:       PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr))

471: !  Get pointer to vector data

473:       PetscCall(VecGetArray(localX,xx,idx,ierr))

475: !  Get local grid boundaries

477:       PetscCall(DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
478:       PetscCall(DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

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

482:       PetscCall(DMGetLocalToGlobalMapping(da,ltogm,ierr))
483:       PetscCall(ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr))

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

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

519:       PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm,ltog,idltog,ierr))

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

526:       PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
527:       PetscCall(VecRestoreArray(localX,xx,idx,ierr))
528:       PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
529:       return
530:       end

532: ! -------------------------------------------------------------------
533: !
534: !   MyMult - user provided matrix multiply
535: !
536: !   Input Parameters:
537: !.  X - input vector
538: !
539: !   Output Parameter:
540: !.  F - function vector
541: !
542:       subroutine  MyMult(J,X,F,ierr)
543:       use ex14fmodule
544:       implicit none

546:       Mat     J
547:       Vec     X,F
548:       PetscErrorCode ierr
549: !
550: !       Here we use the actual formed matrix B; users would
551: !     instead write their own matrix vector product routine
552: !
553:       PetscCall(MatMult(B,X,F,ierr))
554:       return
555:       end

557: !/*TEST
558: !
559: !   test:
560: !      args: -no_output -ksp_gmres_cgs_refinement_type refine_always
561: !      output_file: output/ex14_1.out
562: !      requires: !single
563: !
564: !TEST*/