Actual source code: ex5f.F90

      1: !
      2: !  This example shows how to avoid Fortran line lengths larger than 132 characters.
      3: !  It avoids used of certain macros such as PetscCallA() and PetscCheckA() that
      4: !  generate very long lines
      5: !
      6: !  We recommend starting from src/snes/tutorials/ex5f90.F90 instead of this example
      7: !  because that does not have the restricted formatting that makes this version
      8: !  more difficult to read
      9: !
     10: !  Description: This example solves a nonlinear system in parallel with SNES.
     11: !  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
     12: !  domain, using distributed arrays (DMDAs) to partition the parallel grid.
     13: !  The command line options include:
     14: !    -par <param>, where <param> indicates the nonlinearity of the problem
     15: !       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
     16: !
     17: !  --------------------------------------------------------------------------
     18: !
     19: !  Solid Fuel Ignition (SFI) problem.  This problem is modeled by
     20: !  the partial differential equation
     21: !
     22: !          -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
     23: !
     24: !  with boundary conditions
     25: !
     26: !           u = 0  for  x = 0, x = 1, y = 0, y = 1.
     27: !
     28: !  A finite difference approximation with the usual 5-point stencil
     29: !  is used to discretize the boundary value problem to obtain a nonlinear
     30: !  system of equations.
     31: !
     32: !  --------------------------------------------------------------------------
     33: module ex5fmodule
     34:   use petscsnes
     35:   use petscdmda
     36: #include <petsc/finclude/petscsnes.h>
     37: #include <petsc/finclude/petscdmda.h>
     38:   PetscInt xs, xe, xm, gxs, gxe, gxm
     39:   PetscInt ys, ye, ym, gys, gye, gym
     40:   PetscInt mx, my
     41:   PetscMPIInt rank, size
     42:   PetscReal lambda
     43: end module ex5fmodule

     45: program main
     46:   use ex5fmodule
     47:   implicit none

     49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     50: !                   Variable declarations
     51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     52: !
     53: !  Variables:
     54: !     snes        - nonlinear solver
     55: !     x, r        - solution, residual vectors
     56: !     its         - iterations for convergence
     57: !
     58: !  See additional variable declarations in the file ex5f.h
     59: !
     60:   SNES snes
     61:   Vec x, r
     62:   PetscInt its, i1, i4
     63:   PetscErrorCode ierr
     64:   PetscReal lambda_max, lambda_min
     65:   PetscBool flg
     66:   DM da

     68: !  Note: Any user-defined Fortran routines (such as FormJacobianLocal)
     69: !  MUST be declared as external.

     71:   external FormInitialGuess
     72:   external FormFunctionLocal, FormJacobianLocal
     73:   external MySNESConverged

     75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76: !  Initialize program
     77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     79:   call PetscInitialize(ierr)
     80:   CHKERRA(ierr)
     81:   call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)
     82:   CHKERRMPIA(ierr)
     83:   call MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr)
     84:   CHKERRMPIA(ierr)
     85: !  Initialize problem parameters

     87:   i1 = 1
     88:   i4 = 4
     89:   lambda_max = 6.81
     90:   lambda_min = 0.0
     91:   lambda = 6.0
     92:   call PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', lambda, PETSC_NULL_BOOL, ierr)
     93:   CHKERRA(ierr)

     95: ! this statement is split into multiple-lines to keep lines under 132 char limit - required by 'make check'
     96:   if (lambda >= lambda_max .or. lambda <= lambda_min) then
     97:     ierr = PETSC_ERR_ARG_OUTOFRANGE
     98:     SETERRA(PETSC_COMM_WORLD, ierr, 'Lambda')
     99:   end if

    101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    102: !  Create nonlinear solver context
    103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    105:   call SNESCreate(PETSC_COMM_WORLD, snes, ierr)
    106:   CHKERRA(ierr)

    108: !  Set convergence test routine if desired

    110:   call PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my_snes_convergence', flg, ierr)
    111:   CHKERRA(ierr)
    112:   if (flg) then
    113:     call SNESSetConvergenceTest(snes, MySNESConverged, 0, PETSC_NULL_FUNCTION, ierr)
    114:     CHKERRA(ierr)
    115:   end if

    117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    118: !  Create vector data structures; set function evaluation routine
    119: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    121: !  Create distributed array (DMDA) to manage parallel grid and vectors

    123: !     This really needs only the star-type stencil, but we use the box stencil

    125:   call DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, i4, i4, PETSC_DECIDE, PETSC_DECIDE, &
    126:                     i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, da, ierr)
    127:   CHKERRA(ierr)
    128:   call DMSetFromOptions(da, ierr)
    129:   CHKERRA(ierr)
    130:   call DMSetUp(da, ierr)
    131:   CHKERRA(ierr)

    133: !  Extract global and local vectors from DMDA; then duplicate for remaining
    134: !  vectors that are the same types

    136:   call DMCreateGlobalVector(da, x, ierr)
    137:   CHKERRA(ierr)
    138:   call VecDuplicate(x, r, ierr)
    139:   CHKERRA(ierr)

    141: !  Get local grid boundaries (for 2-dimensional DMDA)

    143:   call DMDAGetInfo(da, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, &
    144:                    PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, &
    145:                    PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr)
    146:   CHKERRA(ierr)
    147:   call DMDAGetCorners(da, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)
    148:   CHKERRA(ierr)
    149:   call DMDAGetGhostCorners(da, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)
    150:   CHKERRA(ierr)

    152: !  Here we shift the starting indices up by one so that we can easily
    153: !  use the Fortran convention of 1-based indices (rather 0-based indices).

    155:   xs = xs + 1
    156:   ys = ys + 1
    157:   gxs = gxs + 1
    158:   gys = gys + 1

    160:   ye = ys + ym - 1
    161:   xe = xs + xm - 1
    162:   gye = gys + gym - 1
    163:   gxe = gxs + gxm - 1

    165: !  Set function evaluation routine and vector

    167:   call DMDASNESSetFunctionLocal(da, INSERT_VALUES, FormFunctionLocal, da, ierr)
    168:   CHKERRA(ierr)
    169:   call DMDASNESSetJacobianLocal(da, FormJacobianLocal, da, ierr)
    170:   CHKERRA(ierr)
    171:   call SNESSetDM(snes, da, ierr)
    172:   CHKERRA(ierr)

    174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    175: !  Customize nonlinear solver; set runtime options
    176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    178: !  Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)

    180:   call SNESSetFromOptions(snes, ierr)
    181:   CHKERRA(ierr)
    182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    183: !  Evaluate initial guess; then solve nonlinear system.
    184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    186: !  Note: The user should initialize the vector, x, with the initial guess
    187: !  for the nonlinear solver prior to calling SNESSolve().  In particular,
    188: !  to employ an initial guess of zero, the user should explicitly set
    189: !  this vector to zero by calling VecSet().

    191:   call FormInitialGuess(x, ierr)
    192:   CHKERRA(ierr)
    193:   call SNESSolve(snes, PETSC_NULL_VEC, x, ierr)
    194:   CHKERRA(ierr)
    195:   call SNESGetIterationNumber(snes, its, ierr)
    196:   CHKERRA(ierr)
    197:   if (rank == 0) then
    198:     write (6, 100) its
    199:   end if
    200: 100 format('Number of SNES iterations = ', i5)

    202: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    203: !  Free work space.  All PETSc objects should be destroyed when they
    204: !  are no longer needed.
    205: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    207:   call VecDestroy(x, ierr)
    208:   CHKERRA(ierr)
    209:   call VecDestroy(r, ierr)
    210:   CHKERRA(ierr)
    211:   call SNESDestroy(snes, ierr)
    212:   CHKERRA(ierr)
    213:   call DMDestroy(da, ierr)
    214:   CHKERRA(ierr)
    215:   call PetscFinalize(ierr)
    216:   CHKERRA(ierr)
    217: end

    219: ! ---------------------------------------------------------------------
    220: !
    221: !  FormInitialGuess - Forms initial approximation.
    222: !
    223: !  Input Parameters:
    224: !  X - vector
    225: !
    226: !  Output Parameter:
    227: !  X - vector
    228: !
    229: !  Notes:
    230: !  This routine serves as a wrapper for the lower-level routine
    231: !  "ApplicationInitialGuess", where the actual computations are
    232: !  done using the standard Fortran style of treating the local
    233: !  vector data as a multidimensional array over the local mesh.
    234: !  This routine merely handles ghost point scatters and accesses
    235: !  the local vector data via VecGetArray() and VecRestoreArray().
    236: !
    237: subroutine FormInitialGuess(X, ierr)
    238:   use ex5fmodule
    239:   implicit none

    241: !  Input/output variables:
    242:   Vec X
    243:   PetscErrorCode ierr
    244: !  Declarations for use with local arrays:
    245:   PetscScalar, pointer :: lx_v(:)

    247:   ierr = 0

    249: !  Get a pointer to vector data.
    250: !    - For default PETSc vectors, VecGetArray() returns a pointer to
    251: !      the data array.  Otherwise, the routine is implementation dependent.
    252: !    - You MUST call VecRestoreArray() when you no longer need access to
    253: !      the array.
    254: !    - Note that the Fortran interface to VecGetArray() differs from the
    255: !      C version.  See the users manual for details.

    257:   call VecGetArray(X, lx_v, ierr)
    258:   CHKERRQ(ierr)

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

    262:   call InitialGuessLocal(lx_v, ierr)
    263:   CHKERRQ(ierr)

    265: !  Restore vector

    267:   call VecRestoreArray(X, lx_v, ierr)
    268:   CHKERRQ(ierr)

    270: end

    272: ! ---------------------------------------------------------------------
    273: !
    274: !  InitialGuessLocal - Computes initial approximation, called by
    275: !  the higher level routine FormInitialGuess().
    276: !
    277: !  Input Parameter:
    278: !  x - local vector data
    279: !
    280: !  Output Parameters:
    281: !  x - local vector data
    282: !  ierr - error code
    283: !
    284: !  Notes:
    285: !  This routine uses standard Fortran-style computations over a 2-dim array.
    286: !
    287: subroutine InitialGuessLocal(x, ierr)
    288:   use ex5fmodule
    289:   implicit none

    291: !  Input/output variables:
    292:   PetscScalar x(xs:xe, ys:ye)
    293:   PetscErrorCode ierr

    295: !  Local variables:
    296:   PetscInt i, j
    297:   PetscReal temp1, temp, one, hx, hy

    299: !  Set parameters

    301:   ierr = 0
    302:   one = 1.0
    303:   hx = one/((real(mx) - 1))
    304:   hy = one/((real(my) - 1))
    305:   temp1 = lambda/(lambda + one)

    307:   do 20 j = ys, ye
    308:     temp = (real(min(j - 1, my - j)))*hy
    309:     do 10 i = xs, xe
    310:       if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
    311:         x(i, j) = 0.0
    312:       else
    313:         x(i, j) = temp1*sqrt(min(real(min(i - 1, mx - i))*hx, (temp)))
    314:       end if
    315: 10    continue
    316: 20    continue

    318:     end

    320: ! ---------------------------------------------------------------------
    321: !
    322: !  FormFunctionLocal - Computes nonlinear function, called by
    323: !  the higher level routine FormFunction().
    324: !
    325: !  Input Parameter:
    326: !  x - local vector data
    327: !
    328: !  Output Parameters:
    329: !  f - local vector data, f(x)
    330: !  ierr - error code
    331: !
    332: !  Notes:
    333: !  This routine uses standard Fortran-style computations over a 2-dim array.
    334: !
    335: !
    336:     subroutine FormFunctionLocal(info, x, f, da, ierr)
    337:       use ex5fmodule
    338:       implicit none

    340:       DM da

    342: !  Input/output variables:
    343:       DMDALocalInfo info
    344:       PetscScalar x(gxs:gxe, gys:gye)
    345:       PetscScalar f(xs:xe, ys:ye)
    346:       PetscErrorCode ierr

    348: !  Local variables:
    349:       PetscScalar two, one, hx, hy
    350:       PetscScalar hxdhy, hydhx, sc
    351:       PetscScalar u, uxx, uyy
    352:       PetscInt i, j

    354:       xs = info%XS + 1
    355:       xe = xs + info%XM - 1
    356:       ys = info%YS + 1
    357:       ye = ys + info%YM - 1
    358:       mx = info%MX
    359:       my = info%MY

    361:       one = 1.0
    362:       two = 2.0
    363:       hx = one/(real(mx) - 1)
    364:       hy = one/(real(my) - 1)
    365:       sc = hx*hy*lambda
    366:       hxdhy = hx/hy
    367:       hydhx = hy/hx

    369: !  Compute function over the locally owned part of the grid

    371:       do 20 j = ys, ye
    372:         do 10 i = xs, xe
    373:           if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
    374:             f(i, j) = x(i, j)
    375:           else
    376:             u = x(i, j)
    377:             uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j))
    378:             uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1))
    379:             f(i, j) = uxx + uyy - sc*exp(u)
    380:           end if
    381: 10        continue
    382: 20        continue

    384:           call PetscLogFlops(11.0d0*ym*xm, ierr)
    385:           CHKERRQ(ierr)

    387:         end

    389: ! ---------------------------------------------------------------------
    390: !
    391: !  FormJacobianLocal - Computes Jacobian matrix, called by
    392: !  the higher level routine FormJacobian().
    393: !
    394: !  Input Parameters:
    395: !  x        - local vector data
    396: !
    397: !  Output Parameters:
    398: !  jac      - Jacobian matrix
    399: !  jac_prec - optionally different matrix used to construct the preconditioner (not used here)
    400: !  ierr     - error code
    401: !
    402: !  Notes:
    403: !  This routine uses standard Fortran-style computations over a 2-dim array.
    404: !
    405: !  Notes:
    406: !  Due to grid point reordering with DMDAs, we must always work
    407: !  with the local grid points, and then transform them to the new
    408: !  global numbering with the "ltog" mapping
    409: !  We cannot work directly with the global numbers for the original
    410: !  uniprocessor grid!
    411: !
    412: !  Two methods are available for imposing this transformation
    413: !  when setting matrix entries:
    414: !    (A) MatSetValuesLocal(), using the local ordering (including
    415: !        ghost points!)
    416: !          by calling MatSetValuesLocal()
    417: !    (B) MatSetValues(), using the global ordering
    418: !        - Use DMDAGetGlobalIndices() to extract the local-to-global map
    419: !        - Then apply this map explicitly yourself
    420: !        - Set matrix entries using the global ordering by calling
    421: !          MatSetValues()
    422: !  Option (A) seems cleaner/easier in many cases, and is the procedure
    423: !  used in this example.
    424: !
    425:         subroutine FormJacobianLocal(info, x, A, jac, da, ierr)
    426:           use ex5fmodule
    427:           implicit none

    429:           DM da

    431: !  Input/output variables:
    432:           PetscScalar x(gxs:gxe, gys:gye)
    433:           Mat A, jac
    434:           PetscErrorCode ierr
    435:           DMDALocalInfo info

    437: !  Local variables:
    438:           PetscInt row, col(5), i, j, i1, i5
    439:           PetscScalar two, one, hx, hy, v(5)
    440:           PetscScalar hxdhy, hydhx, sc

    442: !  Set parameters

    444:           i1 = 1
    445:           i5 = 5
    446:           one = 1.0
    447:           two = 2.0
    448:           hx = one/(real(mx) - 1)
    449:           hy = one/(real(my) - 1)
    450:           sc = hx*hy
    451:           hxdhy = hx/hy
    452:           hydhx = hy/hx
    453: ! -Wmaybe-uninitialized
    454:           v = 0.0
    455:           col = 0

    457: !  Compute entries for the locally owned part of the Jacobian.
    458: !   - Currently, all PETSc parallel matrix formats are partitioned by
    459: !     contiguous chunks of rows across the processors.
    460: !   - Each processor needs to insert only elements that it owns
    461: !     locally (but any non-local elements will be sent to the
    462: !     appropriate processor during matrix assembly).
    463: !   - Here, we set all entries for a particular row at once.
    464: !   - We can set matrix entries either using either
    465: !     MatSetValuesLocal() or MatSetValues(), as discussed above.
    466: !   - Note that MatSetValues() uses 0-based row and column numbers
    467: !     in Fortran as well as in C.

    469:           do 20 j = ys, ye
    470:             row = (j - gys)*gxm + xs - gxs - 1
    471:             do 10 i = xs, xe
    472:               row = row + 1
    473: !           boundary points
    474:               if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
    475: !       Some f90 compilers need 4th arg to be of same type in both calls
    476:                 col(1) = row
    477:                 v(1) = one
    478:                 call MatSetValuesLocal(jac, i1, [row], i1, [col], [v], INSERT_VALUES, ierr)
    479:                 CHKERRQ(ierr)
    480: !           interior grid points
    481:               else
    482:                 v(1) = -hxdhy
    483:                 v(2) = -hydhx
    484:                 v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i, j))
    485:                 v(4) = -hydhx
    486:                 v(5) = -hxdhy
    487:                 col(1) = row - gxm
    488:                 col(2) = row - 1
    489:                 col(3) = row
    490:                 col(4) = row + 1
    491:                 col(5) = row + gxm
    492:                 call MatSetValuesLocal(jac, i1, [row], i5, [col], [v], INSERT_VALUES, ierr)
    493:                 CHKERRQ(ierr)
    494:               end if
    495: 10            continue
    496: 20            continue
    497:               call MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)
    498:               CHKERRQ(ierr)
    499:               call MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)
    500:               CHKERRQ(ierr)
    501:               if (A /= jac) then
    502:                 call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
    503:                 CHKERRQ(ierr)
    504:                 call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)
    505:                 CHKERRQ(ierr)
    506:               end if
    507:             end

    509: !
    510: !     Simple convergence test based on the infinity norm of the residual being small
    511: !
    512:             subroutine MySNESConverged(snes, it, xnorm, snorm, fnorm, reason, dummy, ierr)
    513:               use ex5fmodule
    514:               implicit none

    516:               SNES snes
    517:               PetscInt it, dummy
    518:               PetscReal xnorm, snorm, fnorm, nrm
    519:               SNESConvergedReason reason
    520:               Vec f
    521:               PetscErrorCode ierr

    523:               call SNESGetFunction(snes, f, PETSC_NULL_FUNCTION, dummy, ierr)
    524:               CHKERRQ(ierr)
    525:               call VecNorm(f, NORM_INFINITY, nrm, ierr)
    526:               CHKERRQ(ierr)
    527:               if (nrm <= 1.e-5) reason = SNES_CONVERGED_FNORM_ABS

    529:             end

    531: !/*TEST
    532: !
    533: !   build:
    534: !      requires: !complex !single
    535: !
    536: !   test:
    537: !      nsize: 4
    538: !      args: -snes_mf -pc_type none -da_processors_x 4 -da_processors_y 1 -snes_monitor_short \
    539: !            -ksp_gmres_cgs_refinement_type refine_always
    540: !
    541: !   test:
    542: !      suffix: 2
    543: !      nsize: 4
    544: !      args: -da_processors_x 2 -da_processors_y 2 -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
    545: !
    546: !   test:
    547: !      suffix: 3
    548: !      nsize: 3
    549: !      args: -snes_fd -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
    550: !
    551: !   test:
    552: !      suffix: 6
    553: !      nsize: 1
    554: !      args: -snes_monitor_short -my_snes_convergence
    555: !
    556: !TEST*/