Actual source code: ex1f.F90

      1: !
      2: !
      3: !  Description: Uses the Newton method to solve a two-variable system.
      4: !
      5: !!/*T
      6: !  Concepts: SNES^basic uniprocessor example
      7: !  Processors: 1
      8: !T*/

     10:       program main
     11: #include <petsc/finclude/petsc.h>
     12:       use petsc
     13:       implicit none

     15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     16: !                   Variable declarations
     17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     18: !
     19: !  Variables:
     20: !     snes        - nonlinear solver
     21: !     ksp        - linear solver
     22: !     pc          - preconditioner context
     23: !     ksp         - Krylov subspace method context
     24: !     x, r        - solution, residual vectors
     25: !     J           - Jacobian matrix
     26: !     its         - iterations for convergence
     27: !
     28:       SNES     snes
     29:       PC       pc
     30:       KSP      ksp
     31:       Vec      x,r
     32:       Mat      J
     33:       SNESLineSearch linesearch
     34:       PetscErrorCode  ierr
     35:       PetscInt its,i2,i20
     36:       PetscMPIInt size,rank
     37:       PetscScalar   pfive
     38:       PetscReal   tol
     39:       PetscBool   setls
     40: #if defined(PETSC_USE_LOG)
     41:       PetscViewer viewer
     42: #endif
     43:       double precision threshold,oldthreshold

     45: !  Note: Any user-defined Fortran routines (such as FormJacobian)
     46: !  MUST be declared as external.

     48:       external FormFunction, FormJacobian, MyLineSearch

     50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     51: !                   Macro definitions
     52: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     53: !
     54: !  Macros to make clearer the process of setting values in vectors and
     55: !  getting values from vectors.  These vectors are used in the routines
     56: !  FormFunction() and FormJacobian().
     57: !   - The element lx_a(ib) is element ib in the vector x
     58: !
     59: #define lx_a(ib) lx_v(lx_i + (ib))
     60: #define lf_a(ib) lf_v(lf_i + (ib))
     61: !
     62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     63: !                 Beginning of program
     64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     66:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     67:       if (ierr .ne. 0) then
     68:         print*,'Unable to initialize PETSc'
     69:         stop
     70:       endif
     71:       call PetscLogNestedBegin(ierr);CHKERRA(ierr)
     72:       threshold = 1.0
     73:       call PetscLogSetThreshold(threshold,oldthreshold,ierr)
     74: ! dummy test of logging a reduction
     75: #if defined(PETSC_USE_LOG)
     76:       PetscAReduce()
     77: #endif
     78:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
     79:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
     80:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'Uniprocessor example'); endif

     82:       i2  = 2
     83:       i20 = 20
     84: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
     85: !  Create nonlinear solver context
     86: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

     88:       call SNESCreate(PETSC_COMM_WORLD,snes,ierr)

     90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     91: !  Create matrix and vector data structures; set corresponding routines
     92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

     94: !  Create vectors for solution and nonlinear function

     96:       call VecCreateSeq(PETSC_COMM_SELF,i2,x,ierr)
     97:       call VecDuplicate(x,r,ierr)

     99: !  Create Jacobian matrix data structure

    101:       call MatCreate(PETSC_COMM_SELF,J,ierr)
    102:       call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,i2,i2,ierr)
    103:       call MatSetFromOptions(J,ierr)
    104:       call MatSetUp(J,ierr)

    106: !  Set function evaluation routine and vector

    108:       call SNESSetFunction(snes,r,FormFunction,0,ierr)

    110: !  Set Jacobian matrix data structure and Jacobian evaluation routine

    112:       call SNESSetJacobian(snes,J,J,FormJacobian,0,ierr)

    114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    115: !  Customize nonlinear solver; set runtime options
    116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    118: !  Set linear solver defaults for this problem. By extracting the
    119: !  KSP, KSP, and PC contexts from the SNES context, we can then
    120: !  directly call any KSP, KSP, and PC routines to set various options.

    122:       call SNESGetKSP(snes,ksp,ierr)
    123:       call KSPGetPC(ksp,pc,ierr)
    124:       call PCSetType(pc,PCNONE,ierr)
    125:       tol = 1.e-4
    126:       call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,                  &
    127:      &                      PETSC_DEFAULT_REAL,i20,ierr)

    129: !  Set SNES/KSP/KSP/PC runtime options, e.g.,
    130: !      -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
    131: !  These options will override those specified above as long as
    132: !  SNESSetFromOptions() is called _after_ any other customization
    133: !  routines.

    136:       call SNESSetFromOptions(snes,ierr)

    138:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
    139:      &                         '-setls',setls,ierr)

    141:       if (setls) then
    142:         call SNESGetLineSearch(snes, linesearch, ierr)
    143:         call SNESLineSearchSetType(linesearch, 'shell', ierr)
    144:         call SNESLineSearchShellSetUserFunc(linesearch, MyLineSearch,   &
    145:      &                                      0, ierr)
    146:       endif

    148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    149: !  Evaluate initial guess; then solve nonlinear system
    150: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

    157:       pfive = 0.5
    158:       call VecSet(x,pfive,ierr)
    159:       call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)

    161: !  View solver converged reason; we could instead use the option -snes_converged_reason
    162:       call SNESConvergedReasonView(snes,PETSC_VIEWER_STDOUT_WORLD,ierr)

    164:       call SNESGetIterationNumber(snes,its,ierr);
    165:       if (rank .eq. 0) then
    166:          write(6,100) its
    167:       endif
    168:   100 format('Number of SNES iterations = ',i5)

    170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    171: !  Free work space.  All PETSc objects should be destroyed when they
    172: !  are no longer needed.
    173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    175:       call VecDestroy(x,ierr)
    176:       call VecDestroy(r,ierr)
    177:       call MatDestroy(J,ierr)
    178:       call SNESDestroy(snes,ierr)
    179: #if defined(PETSC_USE_LOG)
    180:       call PetscViewerASCIIOpen(PETSC_COMM_WORLD,'filename.xml',viewer,ierr)
    181:       call PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_XML,ierr)
    182:       call PetscLogView(viewer,ierr)
    183:       call PetscViewerDestroy(viewer,ierr)
    184: #endif
    185:       call PetscFinalize(ierr)
    186:       end
    187: !
    188: ! ------------------------------------------------------------------------
    189: !
    190: !  FormFunction - Evaluates nonlinear function, F(x).
    191: !
    192: !  Input Parameters:
    193: !  snes - the SNES context
    194: !  x - input vector
    195: !  dummy - optional user-defined context (not used here)
    196: !
    197: !  Output Parameter:
    198: !  f - function vector
    199: !
    200:       subroutine FormFunction(snes,x,f,dummy,ierr)
    201:       use petscsnes
    202:       implicit none

    204:       SNES     snes
    205:       Vec      x,f
    206:       PetscErrorCode ierr
    207:       integer dummy(*)

    209: !  Declarations for use with local arrays

    211:       PetscScalar  lx_v(2),lf_v(2)
    212:       PetscOffset  lx_i,lf_i

    214: !  Get pointers to vector data.
    215: !    - For default PETSc vectors, VecGetArray() returns a pointer to
    216: !      the data array.  Otherwise, the routine is implementation dependent.
    217: !    - You MUST call VecRestoreArray() when you no longer need access to
    218: !      the array.
    219: !    - Note that the Fortran interface to VecGetArray() differs from the
    220: !      C version.  See the Fortran chapter of the users manual for details.

    222:       call VecGetArrayRead(x,lx_v,lx_i,ierr)
    223:       call VecGetArray(f,lf_v,lf_i,ierr)

    225: !  Compute function

    227:       lf_a(1) = lx_a(1)*lx_a(1)                                         &
    228:      &          + lx_a(1)*lx_a(2) - 3.0
    229:       lf_a(2) = lx_a(1)*lx_a(2)                                         &
    230:      &          + lx_a(2)*lx_a(2) - 6.0

    232: !  Restore vectors

    234:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)
    235:       call VecRestoreArray(f,lf_v,lf_i,ierr)

    237:       return
    238:       end

    240: ! ---------------------------------------------------------------------
    241: !
    242: !  FormJacobian - Evaluates Jacobian matrix.
    243: !
    244: !  Input Parameters:
    245: !  snes - the SNES context
    246: !  x - input vector
    247: !  dummy - optional user-defined context (not used here)
    248: !
    249: !  Output Parameters:
    250: !  A - Jacobian matrix
    251: !  B - optionally different preconditioning matrix
    252: !  flag - flag indicating matrix structure
    253: !
    254:       subroutine FormJacobian(snes,X,jac,B,dummy,ierr)
    255:       use petscsnes
    256:       implicit none

    258:       SNES         snes
    259:       Vec          X
    260:       Mat          jac,B
    261:       PetscScalar  A(4)
    262:       PetscErrorCode ierr
    263:       PetscInt idx(2),i2
    264:       integer dummy(*)

    266: !  Declarations for use with local arrays

    268:       PetscScalar lx_v(2)
    269:       PetscOffset lx_i

    271: !  Get pointer to vector data

    273:       i2 = 2
    274:       call VecGetArrayRead(x,lx_v,lx_i,ierr)

    276: !  Compute Jacobian entries and insert into matrix.
    277: !   - Since this is such a small problem, we set all entries for
    278: !     the matrix at once.
    279: !   - Note that MatSetValues() uses 0-based row and column numbers
    280: !     in Fortran as well as in C (as set here in the array idx).

    282:       idx(1) = 0
    283:       idx(2) = 1
    284:       A(1) = 2.0*lx_a(1) + lx_a(2)
    285:       A(2) = lx_a(1)
    286:       A(3) = lx_a(2)
    287:       A(4) = lx_a(1) + 2.0*lx_a(2)
    288:       call MatSetValues(B,i2,idx,i2,idx,A,INSERT_VALUES,ierr)

    290: !  Restore vector

    292:       call VecRestoreArrayRead(x,lx_v,lx_i,ierr)

    294: !  Assemble matrix

    296:       call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
    297:       call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
    298:       if (B .ne. jac) then
    299:         call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
    300:         call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
    301:       endif

    303:       return
    304:       end

    307:       subroutine MyLineSearch(linesearch, lctx, ierr)
    308:       use petscsnes
    309:       implicit none

    311:       SNESLineSearch    linesearch
    312:       SNES              snes
    313:       integer           lctx
    314:       Vec               x, f,g, y, w
    315:       PetscReal         ynorm,gnorm,xnorm
    316:       PetscBool         flag
    317:       PetscErrorCode    ierr

    319:       PetscScalar       mone

    321:       mone = -1.0
    322:       call SNESLineSearchGetSNES(linesearch, snes, ierr)
    323:       call SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr)
    324:       call VecNorm(y,NORM_2,ynorm,ierr)
    325:       call VecAXPY(x,mone,y,ierr)
    326:       call SNESComputeFunction(snes,x,f,ierr)
    327:       call VecNorm(f,NORM_2,gnorm,ierr)
    328:       call VecNorm(x,NORM_2,xnorm,ierr)
    329:       call VecNorm(y,NORM_2,ynorm,ierr)
    330:       call SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm,      &
    331:      & ierr)
    332:       flag = PETSC_FALSE
    333:       return
    334:       end

    336: !/*TEST
    337: !
    338: !   test:
    339: !      args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
    340: !      requires: !single
    341: !
    342: !TEST*/