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/petscdmda.h>
     45: #include <petsc/finclude/petscksp.h>
     46:   use petscis
     47:   use petscvec
     48:   use petscdm
     49:   use petscdmda
     50:   use petscksp

     52:   Vec localX
     53:   PetscInt mx, my
     54:   Mat B
     55:   DM da
     56: end module

     58: program main
     59:   use ex14fmodule
     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:   PetscCallA(PetscInitialize(ierr))
     91:   comm = PETSC_COMM_WORLD

     93: !  Initialize problem parameters

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

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

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

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

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

    115: !
    116: !  Create distributed array (DMDA) to manage parallel grid and vectors
    117: !
    118:   Nx = PETSC_DECIDE
    119:   Ny = PETSC_DECIDE
    120:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-Nx', Nx, flg, ierr))
    121:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-Ny', Ny, flg, ierr))
    122:   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))
    123:   PetscCallA(DMSetFromOptions(da, ierr))
    124:   PetscCallA(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:   PetscCallA(DMCreateGlobalVector(da, X, ierr))
    130:   PetscCallA(DMCreateLocalVector(da, localX, ierr))
    131:   PetscCallA(VecDuplicate(X, F, ierr))
    132:   PetscCallA(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:   PetscCallA(VecGetLocalSize(X, m, ierr))
    149:   PetscCallA(MatCreateFromOptions(comm, PETSC_NULL_CHARACTER, one, m, m, N, N, 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 call MatMult()
    157: !     Note: we put B into a module so it will be visible to the
    158: !     mymult() routine
    159:   usemf = .false.
    160:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mf', usemf, ierr))
    161:   if (usemf) then
    162:     PetscCallA(MatCreateShell(comm, m, m, N, N, PETSC_NULL_INTEGER, J, ierr))
    163:     PetscCallA(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:   end if

    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:   PetscCallA(KSPSetFromOptions(ksp, ierr))

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

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

    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:     PetscCallA(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:     PetscCallA(KSPSetOperators(ksp, J, B, ierr))
    218:     PetscCallA(KSPSolve(ksp, F, Y, ierr))

    220: !  Compute updated iterate

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

    231: !  Evaluate nonlinear function at new location

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

    239: !  Test for convergence

    241:     if (fnorm <= ttol) then
    242:       if (.not. nooutput) then
    243:         print *, 'Converged: function norm ', fnorm, ' tolerance ', ttol
    244:       end if
    245:       goto 20
    246:     end if
    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:       PetscCallA(MatIsLinear(J, N, flg, ierr))
    257:       if (.not. flg) then
    258:         print *, 'IsLinear', flg
    259:       end if
    260:     end if

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

    267:     PetscCallA(MatDestroy(B, ierr))
    268:     if (usemf) then
    269:       PetscCallA(MatDestroy(J, ierr))
    270:     end if
    271:     PetscCallA(VecDestroy(localX, ierr))
    272:     PetscCallA(VecDestroy(X, ierr))
    273:     PetscCallA(VecDestroy(Y, ierr))
    274:     PetscCallA(VecDestroy(F, ierr))
    275:     PetscCallA(KSPDestroy(ksp, ierr))
    276:     PetscCallA(DMDestroy(da, ierr))
    277:     PetscCallA(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 ex14fmodule
    292:     implicit none

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

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

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

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

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

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

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

    334: !     Restore vector

    336:         PetscCall(VecRestoreArray(X, xx, ierr))
    337:       end

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

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

    363:         two = 2.0
    364:         one = 1.0
    365:         lambda = 6.0

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

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

    381: !  Get pointers to vector data

    383:         PetscCall(VecGetArrayRead(localX, xx, ierr))
    384:         PetscCall(VecGetArray(F, ff, ierr))

    386: !  Get local grid boundaries

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

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

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

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

    411: !  Restore vectors

    413:             PetscCall(VecRestoreArrayRead(localX, xx, ierr))
    414:             PetscCall(VecRestoreArray(F, ff, ierr))
    415:           end

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

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

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

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

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

    472: !  Get pointer to vector data

    474:             PetscCall(VecGetArrayRead(localX, xx, ierr))

    476: !  Get local grid boundaries

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

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

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

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

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

    520:                 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, ltog, ierr))

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

    527:                 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
    528:                 PetscCall(VecRestoreArrayRead(localX, xx, ierr))
    529:                 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
    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:               end

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