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/petscdm.h>
     38: #include <petsc/finclude/petscdmda.h>
     39:       PetscInt xs,xe,xm,gxs,gxe,gxm
     40:       PetscInt ys,ye,ym,gys,gye,gym
     41:       PetscInt mx,my
     42:       PetscMPIInt rank,size
     43:       PetscReal lambda
     44:       end module ex5fmodule

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

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

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

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

     76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     77: !  Initialize program
     78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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

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

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

    109: !  Set convergence test routine if desired

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

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

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

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

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

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

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

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

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

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

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

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

    166: !  Set function evaluation routine and vector

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

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

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

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

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

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

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

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

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

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

    248:       ierr = 0

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

    258:       call VecGetArrayF90(X,lx_v,ierr)
    259:       CHKERRQ(ierr)

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

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

    266: !  Restore vector

    268:       call VecRestoreArrayF90(X,lx_v,ierr)
    269:       CHKERRQ(ierr)

    271:       end

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

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

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

    300: !  Set parameters

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

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

    319:       end

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

    341:       DM da

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

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

    355:       xs     = info(DMDA_LOCAL_INFO_XS)+1
    356:       xe     = xs+info(DMDA_LOCAL_INFO_XM)-1
    357:       ys     = info(DMDA_LOCAL_INFO_YS)+1
    358:       ye     = ys+info(DMDA_LOCAL_INFO_YM)-1
    359:       mx     = info(DMDA_LOCAL_INFO_MX)
    360:       my     = info(DMDA_LOCAL_INFO_MY)

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

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

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

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

    388:       end

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

    430:       DM da

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

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

    443: !  Set parameters

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

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

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

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

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

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

    530:       end

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