Actual source code: eptorsion2f.F90

      1: !  Program usage: mpiexec -n <proc> eptorsion2f [all TAO options]
      2: !
      3: !  Description:  This example demonstrates use of the TAO package to solve
      4: !  unconstrained minimization problems in parallel.  This example is based
      5: !  on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.
      6: !  The command line options are:
      7: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
      8: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
      9: !    -par <param>, where <param> = angle of twist per unit length
     10: !
     11: !
     12: ! ----------------------------------------------------------------------
     13: !
     14: !  Elastic-plastic torsion problem.
     15: !
     16: !  The elastic plastic torsion problem arises from the deconverged
     17: !  of the stress field on an infinitely long cylindrical bar, which is
     18: !  equivalent to the solution of the following problem:
     19: !     min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
     20: !  where C is the torsion angle per unit length.
     21: !
     22: !  The C version of this code is eptorsion2.c
     23: !
     24: ! ----------------------------------------------------------------------

     26:       module eptorsion2fmodule
     27: #include "petsc/finclude/petsctao.h"
     28:       use petscdmda
     29:       use petsctao
     30:       implicit none

     32:       type(tVec)       localX
     33:       type(tDM)        dm
     34:       PetscReal        param
     35:       PetscInt         mx, my
     36:       end module

     38:       use eptorsion2fmodule
     39:       implicit none
     40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     41: !                   Variable declarations
     42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     43: !
     44: !  See additional variable declarations in the file eptorsion2f.h
     45: !
     46:       PetscErrorCode   ierr           ! used to check for functions returning nonzeros
     47:       type(tVec)       x              ! solution vector
     48:       type(tMat)       H              ! hessian matrix
     49:       PetscInt         Nx, Ny         ! number of processes in x- and y- directions
     50:       type(tTao)       tao            ! Tao solver context
     51:       PetscBool        flg
     52:       PetscInt         i1
     53:       PetscInt         dummy

     55: !  Note: Any user-defined Fortran routines (such as FormGradient)
     56: !  MUST be declared as external.

     58:       external FormInitialGuess,FormFunctionGradient,ComputeHessian
     59:       external Monitor,ConvergenceTest

     61:       i1 = 1

     63: !     Initialize TAO, PETSc  contexts
     64:       PetscCallA(PetscInitialize(ierr))

     66: !     Specify default parameters
     67:       param = 5.0
     68:       mx = 10
     69:       my = 10
     70:       Nx = PETSC_DECIDE
     71:       Ny = PETSC_DECIDE

     73: !     Check for any command line arguments that might override defaults
     74:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
     75:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
     76:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-par',param,flg,ierr))

     78: !     Set up distributed array and vectors
     79:       PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,dm,ierr))
     80:       PetscCallA(DMSetFromOptions(dm,ierr))
     81:       PetscCallA(DMSetUp(dm,ierr))

     83: !     Create vectors
     84:       PetscCallA(DMCreateGlobalVector(dm,x,ierr))
     85:       PetscCallA(DMCreateLocalVector(dm,localX,ierr))

     87: !     Create Hessian
     88:       PetscCallA(DMCreateMatrix(dm,H,ierr))
     89:       PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))

     91: !     The TAO code begins here

     93: !     Create TAO solver
     94:       PetscCallA(TaoCreate(PETSC_COMM_WORLD,tao,ierr))
     95:       PetscCallA(TaoSetType(tao,TAOCG,ierr))

     97: !     Set routines for function and gradient evaluation

     99:       PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr))
    100:       PetscCallA(TaoSetHessian(tao,H,H,ComputeHessian,0,ierr))

    102: !     Set initial guess
    103:       PetscCallA(FormInitialGuess(x,ierr))
    104:       PetscCallA(TaoSetSolution(tao,x,ierr))

    106:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testmonitor',flg,ierr))
    107:       if (flg) then
    108:          PetscCallA(TaoMonitorSet(tao,Monitor,dummy,PETSC_NULL_FUNCTION,ierr))
    109:       endif

    111:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-testconvergence',flg, ierr))
    112:       if (flg) then
    113:          PetscCallA(TaoSetConvergenceTest(tao,ConvergenceTest,dummy,ierr))
    114:       endif

    116: !     Check for any TAO command line options
    117:       PetscCallA(TaoSetFromOptions(tao,ierr))

    119: !     SOLVE THE APPLICATION
    120:       PetscCallA(TaoSolve(tao,ierr))

    122: !     Free TAO data structures
    123:       PetscCallA(TaoDestroy(tao,ierr))

    125: !     Free PETSc data structures
    126:       PetscCallA(VecDestroy(x,ierr))
    127:       PetscCallA(VecDestroy(localX,ierr))
    128:       PetscCallA(MatDestroy(H,ierr))
    129:       PetscCallA(DMDestroy(dm,ierr))

    131: !     Finalize TAO and PETSc
    132:       PetscCallA(PetscFinalize(ierr))
    133:       end

    135: ! ---------------------------------------------------------------------
    136: !
    137: !   FormInitialGuess - Computes an initial approximation to the solution.
    138: !
    139: !   Input Parameters:
    140: !   X    - vector
    141: !
    142: !   Output Parameters:
    143: !   X    - vector
    144: !   ierr - error code
    145: !
    146:       subroutine FormInitialGuess(X,ierr)
    147:       use eptorsion2fmodule
    148:       implicit none

    150: !  Input/output variables:
    151:       Vec              X
    152:       PetscErrorCode   ierr

    154: !  Local variables:
    155:       PetscInt         i, j, k, xe, ye
    156:       PetscReal      temp, val, hx, hy
    157:       PetscInt         xs, ys, xm, ym
    158:       PetscInt         gxm, gym, gxs, gys
    159:       PetscInt         i1

    161:       i1 = 1
    162:       hx = 1.0/real(mx + 1)
    163:       hy = 1.0/real(my + 1)

    165: !  Get corner information
    166:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
    167:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

    169: !  Compute initial guess over locally owned part of mesh
    170:       xe = xs+xm
    171:       ye = ys+ym
    172:       do j=ys,ye-1
    173:          temp = min(j+1,my-j)*hy
    174:          do i=xs,xe-1
    175:             k   = (j-gys)*gxm + i-gxs
    176:             val = min((min(i+1,mx-i))*hx,temp)
    177:             PetscCall(VecSetValuesLocal(X,i1,[k],[val],ADD_VALUES,ierr))
    178:          end do
    179:       end do
    180:       PetscCall(VecAssemblyBegin(X,ierr))
    181:       PetscCall(VecAssemblyEnd(X,ierr))
    182:       end

    184: ! ---------------------------------------------------------------------
    185: !
    186: !  FormFunctionGradient - Evaluates gradient G(X).
    187: !
    188: !  Input Parameters:
    189: !  tao   - the Tao context
    190: !  X     - input vector
    191: !  dummy - optional user-defined context (not used here)
    192: !
    193: !  Output Parameters:
    194: !  f     - the function value at X
    195: !  G     - vector containing the newly evaluated gradient
    196: !  ierr  - error code
    197: !
    198: !  Notes:
    199: !  This routine serves as a wrapper for the lower-level routine
    200: !  "ApplicationGradient", where the actual computations are
    201: !  done using the standard Fortran style of treating the local
    202: !  input vector data as an array over the local mesh.
    203: !
    204:       subroutine FormFunctionGradient(tao,X,f,G,dummy,ierr)
    205:       use eptorsion2fmodule
    206:       implicit none

    208: !  Input/output variables:
    209:       type(tTao)       tao
    210:       type(tVec)       X, G
    211:       PetscReal        f
    212:       PetscErrorCode   ierr
    213:       PetscInt         dummy

    215: !  Declarations for use with local array:

    217:       PetscReal, pointer :: lx_v(:)

    219: !  Local variables:
    220:       PetscReal      zero, p5, area, cdiv3
    221:       PetscReal      val, flin, fquad,floc
    222:       PetscReal      v, vb, vl, vr, vt, dvdx
    223:       PetscReal      dvdy, hx, hy
    224:       PetscInt       xe, ye, xsm, ysm
    225:       PetscInt       xep, yep, i, j, k, ind
    226:       PetscInt       xs, ys, xm, ym
    227:       PetscInt       gxs, gys, gxm, gym
    228:       PetscInt       i1

    230:       i1 = 1
    231:       ierr  = 0
    232:       cdiv3 = param/3.0
    233:       zero = 0.0
    234:       p5   = 0.5
    235:       hx = 1.0/real(mx + 1)
    236:       hy = 1.0/real(my + 1)
    237:       fquad = zero
    238:       flin = zero

    240: !  Initialize gradient to zero
    241:       PetscCall(VecSet(G,zero,ierr))

    243: !  Scatter ghost points to local vector
    244:       PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
    245:       PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))

    247: !  Get corner information
    248:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
    249:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

    251: !  Get pointer to vector data.
    252:       PetscCall(VecGetArrayReadF90(localX,lx_v,ierr))

    254: !  Set local loop dimensions
    255:       xe = xs+xm
    256:       ye = ys+ym
    257:       if (xs .eq. 0) then
    258:          xsm = xs-1
    259:       else
    260:          xsm = xs
    261:       endif
    262:       if (ys .eq. 0) then
    263:          ysm = ys-1
    264:       else
    265:          ysm = ys
    266:       endif
    267:       if (xe .eq. mx) then
    268:          xep = xe+1
    269:       else
    270:          xep = xe
    271:       endif
    272:       if (ye .eq. my) then
    273:          yep = ye+1
    274:       else
    275:          yep = ye
    276:       endif

    278: !     Compute local gradient contributions over the lower triangular elements

    280:       do j = ysm, ye-1
    281:          do i = xsm, xe-1
    282:             k  = (j-gys)*gxm + i-gxs
    283:             v  = zero
    284:             vr = zero
    285:             vt = zero
    286:             if (i .ge. 0 .and. j .ge. 0)      v = lx_v(k+1)
    287:             if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(k+2)
    288:             if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(k+1+gxm)
    289:             dvdx = (vr-v)/hx
    290:             dvdy = (vt-v)/hy
    291:             if (i .ne. -1 .and. j .ne. -1) then
    292:                ind = k
    293:                val = - dvdx/hx - dvdy/hy - cdiv3
    294:                PetscCall(VecSetValuesLocal(G,i1,[k],[val],ADD_VALUES,ierr))
    295:             endif
    296:             if (i .ne. mx-1 .and. j .ne. -1) then
    297:                ind = k+1
    298:                val =  dvdx/hx - cdiv3
    299:                PetscCall(VecSetValuesLocal(G,i1,[ind],[val],ADD_VALUES,ierr))
    300:             endif
    301:             if (i .ne. -1 .and. j .ne. my-1) then
    302:               ind = k+gxm
    303:               val = dvdy/hy - cdiv3
    304:               PetscCall(VecSetValuesLocal(G,i1,[ind],[val],ADD_VALUES,ierr))
    305:             endif
    306:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
    307:             flin = flin - cdiv3 * (v+vr+vt)
    308:          end do
    309:       end do

    311: !     Compute local gradient contributions over the upper triangular elements

    313:       do j = ys, yep-1
    314:          do i = xs, xep-1
    315:             k  = (j-gys)*gxm + i-gxs
    316:             vb = zero
    317:             vl = zero
    318:             v  = zero
    319:             if (i .lt. mx .and. j .gt. 0) vb = lx_v(k+1-gxm)
    320:             if (i .gt. 0 .and. j .lt. my) vl = lx_v(k)
    321:             if (i .lt. mx .and. j .lt. my) v = lx_v(1+k)
    322:             dvdx = (v-vl)/hx
    323:             dvdy = (v-vb)/hy
    324:             if (i .ne. mx .and. j .ne. 0) then
    325:                ind = k-gxm
    326:                val = - dvdy/hy - cdiv3
    327:                PetscCall(VecSetValuesLocal(G,i1,[ind],[val],ADD_VALUES,ierr))
    328:             endif
    329:             if (i .ne. 0 .and. j .ne. my) then
    330:                ind = k-1
    331:                val =  - dvdx/hx - cdiv3
    332:                PetscCall(VecSetValuesLocal(G,i1,[ind],[val],ADD_VALUES,ierr))
    333:             endif
    334:             if (i .ne. mx .and. j .ne. my) then
    335:                ind = k
    336:                val =  dvdx/hx + dvdy/hy - cdiv3
    337:                PetscCall(VecSetValuesLocal(G,i1,[ind],[val],ADD_VALUES,ierr))
    338:             endif
    339:             fquad = fquad + dvdx*dvdx + dvdy*dvdy
    340:             flin = flin - cdiv3*(vb + vl + v)
    341:          end do
    342:       end do

    344: !  Restore vector
    345:       PetscCall(VecRestoreArrayReadF90(localX,lx_v,ierr))

    347: !  Assemble gradient vector
    348:       PetscCall(VecAssemblyBegin(G,ierr))
    349:       PetscCall(VecAssemblyEnd(G,ierr))

    351: ! Scale the gradient
    352:       area = p5*hx*hy
    353:       floc = area *(p5*fquad+flin)
    354:       PetscCall(VecScale(G,area,ierr))

    356: !  Sum function contributions from all processes
    357:       PetscCallMPI(MPI_Allreduce(floc,f,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr))
    358:       PetscCall(PetscLogFlops(20.0d0*(ye-ysm)*(xe-xsm)+16.0d0*(xep-xs)*(yep-ys),ierr))
    359:       end

    361:       subroutine ComputeHessian(tao, X, H, Hpre, dummy, ierr)
    362:       use eptorsion2fmodule
    363:       implicit none

    365:       type(tTao)      tao
    366:       type(tVec)      X
    367:       type(tMat)      H,Hpre
    368:       PetscErrorCode  ierr
    369:       PetscInt        dummy

    371:       PetscInt i,j,k
    372:       PetscInt col(0:4),row
    373:       PetscInt xs,xm,gxs,gxm
    374:       PetscInt ys,ym,gys,gym
    375:       PetscReal v(0:4)
    376:       PetscInt i1

    378:       i1 = 1

    380: !     Get local grid boundaries
    381:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
    382:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

    384:       do j=ys,ys+ym-1
    385:          do i=xs,xs+xm-1
    386:             row = (j-gys)*gxm + (i-gxs)

    388:             k = 0
    389:             if (j .gt. gys) then
    390:                v(k) = -1.0
    391:                col(k) = row-gxm
    392:                k = k + 1
    393:             endif

    395:             if (i .gt. gxs) then
    396:                v(k) = -1.0
    397:                col(k) = row - 1
    398:                k = k +1
    399:             endif

    401:             v(k) = 4.0
    402:             col(k) = row
    403:             k = k + 1

    405:             if (i+1 .lt. gxs + gxm) then
    406:                v(k) = -1.0
    407:                col(k) = row + 1
    408:                k = k + 1
    409:             endif

    411:             if (j+1 .lt. gys + gym) then
    412:                v(k) = -1.0
    413:                col(k) = row + gxm
    414:                k = k + 1
    415:             endif

    417:             PetscCall(MatSetValuesLocal(H,i1,[row],k,col,v,INSERT_VALUES,ierr))
    418:          enddo
    419:       enddo

    421: !     Assemble matrix
    422:       PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
    423:       PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))

    425: !     Tell the matrix we will never add a new nonzero location to the
    426: !     matrix.  If we do it will generate an error.

    428:       PetscCall(MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE,ierr))
    429:       PetscCall(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))

    431:       PetscCall(PetscLogFlops(9.0d0*xm*ym + 49.0d0*xm,ierr))

    433:       ierr = 0
    434:       end

    436:       subroutine Monitor(tao, dummy, ierr)
    437:       use eptorsion2fmodule
    438:       implicit none

    440:       type(tTao)        tao
    441:       PetscInt          dummy
    442:       PetscErrorCode    ierr

    444:       PetscInt           its
    445:       PetscReal          f,gnorm,cnorm,xdiff
    446:       TaoConvergedReason reason

    448:       PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
    449:       if (mod(its,5) .ne. 0) then
    450:          PetscCall(PetscPrintf(PETSC_COMM_WORLD,'iteration multiple of 5\n',ierr))
    451:       endif

    453:       ierr = 0

    455:       end

    457:       subroutine ConvergenceTest(tao, dummy, ierr)
    458:       use eptorsion2fmodule
    459:       implicit none

    461:       type(tTao)          tao
    462:       PetscInt           dummy
    463:       PetscErrorCode     ierr

    465:       PetscInt           its
    466:       PetscReal          f,gnorm,cnorm,xdiff
    467:       TaoConvergedReason reason

    469:       PetscCall(TaoGetSolutionStatus(tao,its,f,gnorm,cnorm,xdiff,reason,ierr))
    470:       if (its .eq. 7) then
    471:        PetscCall(TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS,ierr))
    472:       endif

    474:       ierr = 0

    476:       end

    478: !/*TEST
    479: !
    480: !   build:
    481: !      requires: !complex
    482: !
    483: !   test:
    484: !      args: -tao_monitor_short -tao_type nls -tao_gttol 1.e-2
    485: !
    486: !   test:
    487: !      suffix: 2
    488: !      nsize: 2
    489: !      args: -tao_monitor_short -tao_type lmvm -tao_gttol 1.e-2
    490: !
    491: !   test:
    492: !      suffix: 3
    493: !      args: -testmonitor -tao_type lmvm -tao_gttol 1.e-2
    494: !TEST*/