Actual source code: ex1f.F90

      1: !
      2: !  Description: Uses the Newton method to solve a two-variable system.
      3: !
      4: #include <petsc/finclude/petsc.h>
      5: module ex1fmodule
      6:   use petscvec
      7:   use petscsnesdef
      8:   use petscvec
      9:   use petscmat
     10:   implicit none

     12: contains
     13: !
     14: ! ------------------------------------------------------------------------
     15: !
     16: !  FormFunction - Evaluates nonlinear function, F(x).
     17: !
     18: !  Input Parameters:
     19: !  snes - the SNES context
     20: !  x - input vector
     21: !  dummy - optional user-defined context (not used here)
     22: !
     23: !  Output Parameter:
     24: !  f - function vector
     25: !
     26:   subroutine FormFunction(snes, x, f, dummy, ierr)
     27:     SNES snes
     28:     Vec x, f
     29:     PetscErrorCode, intent(out) :: ierr
     30:     integer dummy(*)

     32: !  Declarations for use with local arrays
     33:     PetscScalar, pointer :: lx_v(:), lf_v(:)

     35: !  Get pointers to vector data.
     36: !    - VecGetArray() returns a pointer to the data array.
     37: !    - You MUST call VecRestoreArray() when you no longer need access to
     38: !      the array.

     40:     PetscCall(VecGetArrayRead(x, lx_v, ierr))
     41:     PetscCall(VecGetArray(f, lf_v, ierr))

     43: !  Compute function

     45:     lf_v(1) = lx_v(1)*lx_v(1) + lx_v(1)*lx_v(2) - 3.0
     46:     lf_v(2) = lx_v(1)*lx_v(2) + lx_v(2)*lx_v(2) - 6.0

     48: !  Restore vectors

     50:     PetscCall(VecRestoreArrayRead(x, lx_v, ierr))
     51:     PetscCall(VecRestoreArray(f, lf_v, ierr))

     53:   end

     55: ! ---------------------------------------------------------------------
     56: !
     57: !  FormJacobian - Evaluates Jacobian matrix.
     58: !
     59: !  Input Parameters:
     60: !  snes - the SNES context
     61: !  x - input vector
     62: !  dummy - optional user-defined context (not used here)
     63: !
     64: !  Output Parameters:
     65: !  A - Jacobian matrix
     66: !  B - optionally different matrix used to construct the preconditioner
     67: !
     68:   subroutine FormJacobian(snes, X, jac, B, dummy, ierr)

     70:     SNES snes
     71:     Vec X
     72:     Mat jac, B
     73:     PetscScalar A(4)
     74:     PetscErrorCode, intent(out) :: ierr
     75:     PetscInt idx(2)
     76:     integer dummy(*)

     78: !  Declarations for use with local arrays

     80:     PetscScalar, pointer :: lx_v(:)

     82: !  Get pointer to vector data

     84:     PetscCall(VecGetArrayRead(x, lx_v, ierr))

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

     92:     idx = [0, 1]
     93:     A = [2.0*lx_v(1) + lx_v(2), lx_v(1), lx_v(2), lx_v(1) + 2.0*lx_v(2)]
     94:     PetscCall(MatSetValues(B, 2_PETSC_INT_KIND, idx, 2_PETSC_INT_KIND, idx, A, INSERT_VALUES, ierr))

     96: !  Restore vector

     98:     PetscCall(VecRestoreArrayRead(x, lx_v, ierr))

    100: !  Assemble matrix

    102:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
    103:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
    104:     if (B /= jac) then
    105:       PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
    106:       PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
    107:     end if

    109:   end

    111:   subroutine MyLineSearch(linesearch, lctx, ierr)

    113:     SNESLineSearch linesearch
    114:     SNES snes
    115:     integer lctx
    116:     Vec x, f, g, y, w
    117:     PetscReal ynorm, gnorm, xnorm
    118:     PetscErrorCode, intent(out) :: ierr

    120:     PetscScalar, parameter :: mone = -1.0

    122:     PetscCall(SNESLineSearchGetSNES(linesearch, snes, ierr))
    123:     PetscCall(SNESLineSearchGetVecs(linesearch, x, f, y, w, g, ierr))
    124:     PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
    125:     PetscCall(VecAXPY(x, mone, y, ierr))
    126:     PetscCall(SNESComputeFunction(snes, x, f, ierr))
    127:     PetscCall(VecNorm(f, NORM_2, gnorm, ierr))
    128:     PetscCall(VecNorm(x, NORM_2, xnorm, ierr))
    129:     PetscCall(VecNorm(y, NORM_2, ynorm, ierr))
    130:     PetscCall(SNESLineSearchSetNorms(linesearch, xnorm, gnorm, ynorm, ierr))
    131:   end
    132: end module

    134: program main
    135:   use petsc
    136:   use ex1fmodule
    137:   implicit none

    139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    140: !                   Variable declarations
    141: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    142: !
    143: !  Variables:
    144: !     snes        - nonlinear solver
    145: !     ksp        - linear solver
    146: !     pc          - preconditioner context
    147: !     ksp         - Krylov subspace method context
    148: !     x, r        - solution, residual vectors
    149: !     J           - Jacobian matrix
    150: !     its         - iterations for convergence
    151: !
    152:   SNES snes
    153:   PC pc
    154:   KSP ksp
    155:   Vec x, r
    156:   Mat J
    157:   SNESLineSearch linesearch
    158:   PetscErrorCode ierr
    159:   PetscInt its
    160:   PetscMPIInt size, rank
    161:   PetscScalar, parameter :: pfive = 0.5
    162:   PetscReal, parameter :: tol = 1.e-4
    163:   PetscBool setls
    164:   PetscReal, pointer :: rhistory(:)
    165:   PetscInt, pointer :: itshistory(:)
    166:   PetscInt nhistory
    167: #if defined(PETSC_USE_LOG)
    168:   PetscViewer viewer
    169: #endif
    170:   double precision threshold, oldthreshold

    172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    173: !                 Beginning of program
    174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    176:   PetscCallA(PetscInitialize(ierr))
    177:   PetscCallA(PetscLogNestedBegin(ierr))
    178:   threshold = 1.0
    179:   PetscCallA(PetscLogSetThreshold(threshold, oldthreshold, ierr))
    180:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
    181:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
    182:   PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'Uniprocessor example')

    184: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -
    185: !  Create nonlinear solver context
    186: ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - -

    188:   PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))

    190:   PetscCallA(SNESSetConvergenceHistory(snes, PETSC_NULL_REAL_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_DECIDE, PETSC_FALSE, ierr))

    192: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    193: !  Create matrix and vector data structures; set corresponding routines
    194: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    196: !  Create vectors for solution and nonlinear function

    198:   PetscCallA(VecCreateSeq(PETSC_COMM_SELF, 2_PETSC_INT_KIND, x, ierr))
    199:   PetscCallA(VecDuplicate(x, r, ierr))

    201: !  Create Jacobian matrix data structure

    203:   PetscCallA(MatCreate(PETSC_COMM_SELF, J, ierr))
    204:   PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, 2_PETSC_INT_KIND, 2_PETSC_INT_KIND, ierr))
    205:   PetscCallA(MatSetFromOptions(J, ierr))
    206:   PetscCallA(MatSetUp(J, ierr))

    208: !  Set function evaluation routine and vector

    210:   PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr))

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

    214:   PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr))

    216: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    217: !  Customize nonlinear solver; set runtime options
    218: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

    224:   PetscCallA(SNESGetKSP(snes, ksp, ierr))
    225:   PetscCallA(KSPGetPC(ksp, pc, ierr))
    226:   PetscCallA(PCSetType(pc, PCNONE, ierr))
    227:   PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, 20_PETSC_INT_KIND, ierr))

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

    235:   PetscCallA(SNESSetFromOptions(snes, ierr))

    237:   PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-setls', setls, ierr))

    239:   if (setls) then
    240:     PetscCallA(SNESGetLineSearch(snes, linesearch, ierr))
    241:     PetscCallA(SNESLineSearchSetType(linesearch, 'shell', ierr))
    242:     PetscCallA(SNESLineSearchShellSetApply(linesearch, MyLineSearch, 0, ierr))
    243:   end if

    245: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    246: !  Evaluate initial guess; then solve nonlinear system
    247: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

    254:   PetscCallA(VecSet(x, pfive, ierr))
    255:   PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))

    257:   PetscCallA(SNESGetConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))
    258:   PetscCallA(SNESRestoreConvergenceHistory(snes, rhistory, itshistory, nhistory, ierr))

    260: ! View solver converged reason; we could instead use the option -snes_converged_reason
    261:   PetscCallA(SNESConvergedReasonView(snes, PETSC_VIEWER_STDOUT_WORLD, ierr))

    263:   PetscCallA(SNESGetIterationNumber(snes, its, ierr))
    264:   if (rank == 0) then
    265:     write (6, 100) its
    266:   end if
    267: 100 format('Number of SNES iterations = ', i5)

    269: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    270: !  Free work space.  All PETSc objects should be destroyed when they
    271: !  are no longer needed.
    272: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    274:   PetscCallA(VecDestroy(x, ierr))
    275:   PetscCallA(VecDestroy(r, ierr))
    276:   PetscCallA(MatDestroy(J, ierr))
    277:   PetscCallA(SNESDestroy(snes, ierr))
    278: #if defined(PETSC_USE_LOG)
    279:   PetscCallA(PetscViewerASCIIOpen(PETSC_COMM_WORLD, 'filename.xml', viewer, ierr))
    280:   PetscCallA(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_XML, ierr))
    281:   PetscCallA(PetscLogView(viewer, ierr))
    282:   PetscCallA(PetscViewerDestroy(viewer, ierr))
    283: #endif
    284:   PetscCallA(PetscFinalize(ierr))
    285: end
    286: !/*TEST
    287: !
    288: !   test:
    289: !      args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
    290: !      requires: !single
    291: !
    292: !TEST*/