Actual source code: plate2f.F

      1: ! "$Id: plate2f.F 1.47 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
      2: !
      3: !  Program usage: mpirun -np <proc> plate2f [all TAO options]
      4: !
      5: !  This example demonstrates use of the TAO package to solve a bound constrained
      6: !  minimization problem.  This example is based on a problem from the
      7: !  MINPACK-2 test suite.  Given a rectangular 2-D domain and boundary values
      8: !  along the edges of the domain, the objective is to find the surface
      9: !  with the minimal area that satisfies the boundary conditions.
     10: !  The command line options are:
     11: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
     12: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
     13: !    -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
     14: !    -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
     15: !    -bheight <ht>, where <ht> = height of the plate
     16: !
     17: !/*T
     18: !   Concepts: TAO - Solving a bound constrained minimization problem
     19: !   Routines: TaoInitialize(); TaoFinalize();
     20: !   Routines: TaoCreate(); TaoDestroy();
     21: !   Routines: TaoAppSetObjectiveAndGradientRoutine();
     22: !   Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
     23: !   Routines: TaoAppSetInitialSolutionVec(); TaoAppSetVariableBounds();
     24: !   Routines: TaoSetApplication(); TaoSetOptions();
     25: !   Routines: TaoApplicationCreate(); TaoSolve();
     26: !   Routines: TaoView(); TaoAppDestroy();
     27: !   Processors: n
     28: !T*/

     32:       implicit none

     34: #include "plate2f.h"

     36: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     37: !                   Variable declarations
     38: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     39: !
     40: !  Variables:
     41: !    (common from plate2f.h):
     42: !    Nx, Ny           number of processors in x- and y- directions
     43: !    mx, my           number of grid points in x,y directions
     44: !    N    global dimension of vector

     46:       integer          info          ! used to check for functions returning nonzeros
     47:       Vec              x             ! solution vector
     48:       Vec              xl, xu        ! lower and upper bounds vectorsp
     49:       integer          m             ! number of local elements in vector
     50:       TAO_SOLVER       tao           ! TAO_SOLVER solver context
     51:       TAO_APPLICATION  plateapp      ! PETSc application
     52:       Mat              H             ! Hessian matrix
     53:       ISLocalToGlobalMapping isltog  ! local to global mapping object
     54:       PetscTruth       flg

     57:       external FormFunctionGradient
     58:       external FormHessian
     59:       external MSA_BoundaryConditions
     60:       external MSA_Plate
     61:       external MSA_InitialPoint
     62: ! Initialize Tao
     63:       call PetscInitialize(PETSC_NULL_CHARACTER,info)
     64:       call TaoInitialize(PETSC_NULL_CHARACTER,info)
     65:

     67: ! Specify default dimensions of the problem
     68:       mx = 10
     69:       my = 10
     70:       bheight = 0.1

     72: ! Check for any command line arguments that override defaults
     73:
     74:       call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-mx",mx,flg,info)
     75:       call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-my",my,flg,info)
     76:
     77:       bmx = mx/2
     78:       bmy = my/2

     80:       call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmx",bmx,flg,info)
     81:       call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmy",bmy,flg,info)
     82:       call PetscOptionsGetReal(TAO_NULL_CHARACTER,"-bheight",bheight,   &
     83:      &      flg,info)
     84:

     86: ! Calculate any derived values from parameters
     87:       N = mx*my

     89: ! Let Petsc determine the dimensions of the local vectors
     90:       Nx = PETSC_DECIDE
     91:       NY = PETSC_DECIDE

     93: ! A two dimensional distributed array will help define this problem, which
     94: ! derives from an elliptic PDE on a two-dimensional domain.  From the
     95: ! distributed array, create the vectors

     97:       call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,         &
     98:      &     mx,my,Nx,Ny,1,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,           &
     99:      &     da,info)
    100:

    102: ! Extract global and local vectors from DA; The local vectors are
    103: ! used solely as work space for the evaluation of the function,
    104: ! gradient, and Hessian.  Duplicate for remaining vectors that are
    105: ! the same types.

    107:       call DACreateGlobalVector(da,x,info)
    108:       call DACreateLocalVector(da,localX,info)
    109:       call VecDuplicate(localX,localV,info)

    111: ! Create a matrix data structure to store the Hessian.
    112: ! Here we (optionally) also associate the local numbering scheme
    113: ! with the matrix so that later we can use local indices for matrix
    114: ! assembly

    116:       call VecGetLocalSize(x,m,info)
    117:       call MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL_INTEGER,  &
    118:      &     3,PETSC_NULL_INTEGER,H,info)

    120:       call MatSetOption(H,MAT_SYMMETRIC,info)
    121:       call DAGetISLocalToGlobalMapping(da,isltog,info)
    122:       call MatSetLocalToGlobalMapping(H,isltog,info)
    123:

    125: ! The Tao code begins here
    126: ! Create TAO solver and set desired solution method.
    127: ! This problems uses bounded variables, so the
    128: ! method must either be 'tao_tron' or 'tao_blmvm'

    130:       call TaoCreate(MPI_COMM_WORLD,'tao_blmvm',tao,info)
    131:       call TaoApplicationCreate(MPI_COMM_WORLD,plateapp,info)
    132:

    134: !     Set minimization function and gradient, hessian evaluation functions

    136: !     TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
    137:       call TaoAppSetObjectiveAndGradientRo(plateapp,                    &
    138:      &     FormFunctionGradient,PETSC_NULL_OBJECT,info)

    140:       call TaoAppSetHessianMat(plateapp,H,H,info)
    141:       call TaoAppSetHessianRoutine(plateapp,FormHessian,                 &
    142:      &     PETSC_NULL_OBJECT, info)
    143:

    145: ! Set Variable bounds
    146:       call MSA_BoundaryConditions(info)
    147:       call VecDuplicate(x,xl,info)
    148:       call VecDuplicate(x,xu,info)
    149:       call MSA_Plate(xl,xu,info)
    150:       call TaoAppSetVariableBounds(plateapp,xl,xu,info)

    152: ! Set the initial solution guess
    153:       call MSA_InitialPoint(x, info)
    154:       call TaoAppSetInitialSolutionVec(plateapp,x,info)

    156: ! Now that the PETSc application is all set, attach to tao context
    157:       call TaoSetApplication(tao,plateapp,info)

    159: ! Check for any tao command line options
    160:       call TaoSetOptions(plateapp,tao,info)

    162: ! Solve the application
    163:       call TaoSolveApplication(plateapp,tao,info)

    165: !     View TAO solver information
    166: !      call TaoView(tao,info)

    168: ! Free TAO data structures
    169:       call TaoDestroy(tao,info)
    170:       call TaoAppDestroy(plateapp,info)

    172: ! Free PETSc data structures
    173:       call VecDestroy(x,info)
    174:       call VecDestroy(xl,info)
    175:       call VecDestroy(xu,info)
    176:       call VecDestroy(Top,info)
    177:       call VecDestroy(Bottom,info)
    178:       call VecDestroy(Left,info)
    179:       call VecDestroy(Right,info)
    180:       call MatDestroy(H,info)
    181:       call VecDestroy(localX,info)
    182:       call VecDestroy(localV,info)
    183:       call DADestroy(da,info)

    185: ! Finalize TAO

    187:       call TaoFinalize(info)
    188:       call PetscFinalize(info)

    190:       end

    192: ! ---------------------------------------------------------------------
    193: !
    194: !  FormFunctionGradient - Evaluates function f(X).
    195: !
    196: !  Input Parameters:
    197: !  tao   - the TAO_SOLVER context
    198: !  X     - the input vector
    199: !  dummy - optional user-defined context, as set by TaoSetFunction()
    200: !          (not used here)
    201: !
    202: !  Output Parameters:
    203: !  fcn     - the newly evaluated function
    204: !  G       - the gradient vector
    205: !  info  - error code
    206: !

    209:       subroutine FormFunctionGradient(tao,X,fcn,G,dummy,info)
    210:       implicit none

    212: ! da, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
    213: #include "plate2f.h"
    214:
    215: ! Input/output variables

    217:       TAO_SOLVER       tao
    218:       PetscScalar fcn
    219:       Vec              X, G
    220:       integer          dummy, info
    221:
    222:       integer          i,j,row
    223:       integer          xs, xm, gxs, gxm, ys, ym, gys, gym
    224:       PetscScalar      ft,zero,hx,hy,hydhx,hxdhy,area,rhx,rhy
    225:       PetscScalar      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8
    226:       PetscScalar      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc
    227:       PetscScalar      xc,xl,xr,xt,xb,xlt,xrb

    230: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    231: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
    232: ! will return an array of doubles referenced by x_array offset by x_index.
    233: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    234: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    235:       PetscScalar          g_v(0:1),x_v(0:1), top_v(0:1)
    236:       PetscScalar          left_v(0:1), right_v(0:1),bottom_v(0:1)
    237:       PetscOffset          g_i,left_i,right_i,bottom_i,top_i
    238:       PetscOffset          x_i

    240:       ft = 0.0d0
    241:       zero = 0.0d0
    242:       hx = 1.0d0/(mx + 1)
    243:       hy = 1.0d0/(my + 1)
    244:       hydhx = hy/hx
    245:       hxdhy = hx/hy
    246:       area = 0.5d0 * hx * hy
    247:       rhx = mx + 1.0d0
    248:       rhy = my + 1.0d0

    251: ! Get local mesh boundaries
    252:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
    253:      &                  PETSC_NULL_INTEGER,info)
    254:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,             &
    255:      &                       gxm,gym,PETSC_NULL_INTEGER,info)

    257: ! Scatter ghost points to local vector
    258:       call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
    259:       call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)

    261: ! Initialize the vector to zero
    262:       call VecSet(localV,zero,info)

    264: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
    265:       call VecGetArray(localX,x_v,x_i,info)
    266:       call VecGetArray(localV,g_v,g_i,info)
    267:       call VecGetArray(Top,top_v,top_i,info)
    268:       call VecGetArray(Bottom,bottom_v,bottom_i,info)
    269:       call VecGetArray(Left,left_v,left_i,info)
    270:       call VecGetArray(Right,right_v,right_i,info)

    272: ! Compute function over the locally owned part of the mesh
    273:       do j = ys,ys+ym-1
    274:          do i = xs,xs+xm-1
    275:             row = (j-gys)*gxm + (i-gxs)
    276:             xc = x_v(row+x_i)
    277:             xt = xc
    278:             xb = xc
    279:             xr = xc
    280:             xl = xc
    281:             xrb = xc
    282:             xlt = xc

    284:             if (i .eq. 0) then !left side
    285:                xl = left_v(j - ys + 1 + left_i)
    286:                xlt = left_v(j - ys + 2 + left_i)
    287:             else
    288:                xl = x_v(row - 1 + x_i)
    289:             endif

    291:             if (j .eq. 0) then !bottom side
    292:                xb = bottom_v(i - xs + 1 + bottom_i)
    293:                xrb = bottom_v(i - xs + 2 + bottom_i)
    294:             else
    295:                xb = x_v(row - gxm + x_i)
    296:             endif

    298:             if (i + 1 .eq. gxs + gxm) then !right side
    299:                xr = right_v(j - ys + 1 + right_i)
    300:                xrb = right_v(j - ys + right_i)
    301:             else
    302:                xr = x_v(row + 1 + x_i)
    303:             endif

    305:             if (j + 1 .eq. gys + gym) then !top side
    306:                xt = top_v(i - xs + 1 + top_i)
    307:                xlt = top_v(i - xs + top_i)
    308:             else
    309:                xt = x_v(row + gxm + x_i)
    310:             endif

    312:             if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
    313:                xlt = x_v(row - 1 + gxm + x_i)
    314:             endif

    316:             if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
    317:                xrb = x_v(row + 1 - gxm + x_i)
    318:             endif

    320:             d1 = xc-xl
    321:             d2 = xc-xr
    322:             d3 = xc-xt
    323:             d4 = xc-xb
    324:             d5 = xr-xrb
    325:             d6 = xrb-xb
    326:             d7 = xlt-xl
    327:             d8 = xt-xlt

    329:             df1dxc = d1 * hydhx
    330:             df2dxc = d1 * hydhx + d4 * hxdhy
    331:             df3dxc = d3 * hxdhy
    332:             df4dxc = d2 * hydhx + d3 * hxdhy
    333:             df5dxc = d2 * hydhx
    334:             df6dxc = d4 * hxdhy

    336:             d1 = d1 * rhx
    337:             d2 = d2 * rhx
    338:             d3 = d3 * rhy
    339:             d4 = d4 * rhy
    340:             d5 = d5 * rhy
    341:             d6 = d6 * rhx
    342:             d7 = d7 * rhy
    343:             d8 = d8 * rhx

    345:             f1 = sqrt(1.0d0 + d1*d1 + d7*d7)
    346:             f2 = sqrt(1.0d0 + d1*d1 + d4*d4)
    347:             f3 = sqrt(1.0d0 + d3*d3 + d8*d8)
    348:             f4 = sqrt(1.0d0 + d3*d3 + d2*d2)
    349:             f5 = sqrt(1.0d0 + d2*d2 + d5*d5)
    350:             f6 = sqrt(1.0d0 + d4*d4 + d6*d6)

    352:             ft = ft + f2 + f4

    354:             df1dxc = df1dxc / f1
    355:             df2dxc = df2dxc / f2
    356:             df3dxc = df3dxc / f3
    357:             df4dxc = df4dxc / f4
    358:             df5dxc = df5dxc / f5
    359:             df6dxc = df6dxc / f6

    361:             g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc +  &
    362:      &                              df5dxc + df6dxc)
    363:          enddo
    364:       enddo

    366: ! Compute triangular areas along the border of the domain.
    367:       if (xs .eq. 0) then  ! left side
    368:          do j=ys,ys+ym-1
    369:             d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i))         &
    370:      &                 * rhy
    371:             d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i))        &
    372:      &                 * rhx
    373:             ft = ft + sqrt(1.0d0 + d3*d3 + d2*d2)
    374:          enddo
    375:       endif

    377:
    378:       if (ys .eq. 0) then !bottom side
    379:          do i=xs,xs+xm-1
    380:             d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i))    &
    381:      &                    * rhx
    382:             d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
    383:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
    384:          enddo
    385:       endif

    387:
    388:       if (xs + xm .eq. mx) then ! right side
    389:          do j=ys,ys+ym-1
    390:             d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
    391:             d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
    392:             ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
    393:          enddo
    394:       endif

    396:
    397:       if (ys + ym .eq. my) then
    398:          do i=xs,xs+xm-1
    399:             d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
    400:             d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
    401:             ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
    402:          enddo
    403:       endif

    405:
    406:       if ((ys .eq. 0) .and. (xs .eq. 0)) then
    407:          d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
    408:          d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
    409:          ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
    410:       endif

    412:       if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
    413:          d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
    414:          d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
    415:          ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
    416:       endif

    418:       ft = ft * area
    419:       call MPI_Allreduce(ft,fcn,1,MPI_DOUBLE_PRECISION,                  &
    420:      &             MPI_SUM,MPI_COMM_WORLD,info)

    424: ! Restore vectors
    425:       call VecRestoreArray(localX,x_v,x_i,info)
    426:       call VecRestoreArray(localV,g_v,g_i,info)
    427:       call VecRestoreArray(Left,left_v,left_i,info)
    428:       call VecRestoreArray(Top,top_v,top_i,info)
    429:       call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
    430:       call VecRestoreArray(Right,right_v,right_i,info)

    432: ! Scatter values to global vector
    433:       call DALocalToGlobal(da,localV,INSERT_VALUES,G,info)

    435:       call PetscLogFlops(70*xm*ym,info)

    437:       return
    438:       end  !FormFunctionGradient
    439:

    444: ! ----------------------------------------------------------------------------
    445: !
    446: !/*
    447: !   FormHessian - Evaluates Hessian matrix.
    448: !
    449: !   Input Parameters:
    450: !.  tao  - the TAO_SOLVER context
    451: !.  X    - input vector
    452: !.  dummy  - not used
    453: !
    454: !   Output Parameters:
    455: !.  Hessian    - Hessian matrix
    456: !.  Hpc    - optionally different preconditioning matrix
    457: !.  flag - flag indicating matrix structure
    458: !
    459: !   Notes:
    460: !   Due to mesh point reordering with DAs, we must always work
    461: !   with the local mesh points, and then transform them to the new
    462: !   global numbering with the local-to-global mapping.  We cannot work
    463: !   directly with the global numbers for the original uniprocessor mesh!
    464: !
    465: !   Two methods are available for imposing this transformation
    466: !   when setting matrix entries:
    467: !     (A) MatSetValuesLocal(), using the local ordering (including
    468: !         ghost points!)
    469: !         - Do the following two steps once, before calling TaoSolve()
    470: !           - Use DAGetISLocalToGlobalMapping() to extract the
    471: !             local-to-global map from the DA
    472: !           - Associate this map with the matrix by calling
    473: !             MatSetLocalToGlobalMapping()
    474: !         - Then set matrix entries using the local ordering
    475: !           by calling MatSetValuesLocal()
    476: !     (B) MatSetValues(), using the global ordering
    477: !         - Use DAGetGlobalIndices() to extract the local-to-global map
    478: !         - Then apply this map explicitly yourself
    479: !         - Set matrix entries using the global ordering by calling
    480: !           MatSetValues()
    481: !   Option (A) seems cleaner/easier in many cases, and is the procedure
    482: !   used in this example.
    483: */
    484:       subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, info)
    485:       implicit none

    487: ! da,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
    488: #include "plate2f.h"
    489:
    490:       TAO_SOLVER       tao
    491:       Vec              X
    492:       Mat              Hessian,Hpc
    493:       MatStructure     flg
    494:       integer          dummy,info

    496:       integer          i,j,k,row
    497:       integer          xs,xm,gxs,gxm,ys,ym,gys,gym,col(0:6)
    498:       PetscScalar      hx,hy,hydhx,hxdhy,rhx,rhy
    499:       PetscScalar      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8
    500:       PetscScalar      xc,xl,xr,xt,xb,xlt,xrb
    501:       PetscScalar      hl,hr,ht,hb,hc,htl,hbr

    503: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    504: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
    505: ! will return an array of doubles referenced by x_array offset by x_index.
    506: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    507: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    508:       PetscScalar      right_v(0:1),left_v(0:1),bottom_v(0:1),top_v(0:1)
    509:       PetscScalar      x_v(0:1)
    510:       PetscOffset      x_i,right_i,left_i,bottom_i,top_i
    511:       PetscScalar      v(0:6)
    512:       PetscTruth       assembled
    513:
    514: ! Set various matrix options
    515:       call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,info)
    516:       call MatSetOption(Hessian,MAT_COLUMNS_SORTED,info)
    517:       call MatSetOption(Hessian,MAT_ROWS_SORTED,info)

    519: ! Get local mesh boundaries
    520:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
    521:      &                  PETSC_NULL_INTEGER,info)
    522:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
    523:      &                       PETSC_NULL_INTEGER,info)

    525: ! Scatter ghost points to local vectors
    526:       call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
    527:       call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)

    529: ! Get pointers to vector data (see note on Fortran arrays above)
    530:       call VecGetArray(localX,x_v,x_i,info)
    531:       call VecGetArray(Top,top_v,top_i,info)
    532:       call VecGetArray(Bottom,bottom_v,bottom_i,info)
    533:       call VecGetArray(Left,left_v,left_i,info)
    534:       call VecGetArray(Right,right_v,right_i,info)

    536: ! Initialize matrix entries to zero
    537:       call MatAssembled(Hessian,assembled,info)
    538:       if (assembled .eq. PETSC_TRUE) call MatZeroEntries(Hessian,info)

    541:       rhx = mx + 1.0
    542:       rhy = my + 1.0
    543:       hx = 1.0/rhx
    544:       hy = 1.0/rhy
    545:       hydhx = hy/hx
    546:       hxdhy = hx/hy
    547: ! compute Hessian over the locally owned part of the mesh

    549:       do  i=xs,xs+xm-1
    550:          do  j=ys,ys+ym-1
    551:             row = (j-gys)*gxm + (i-gxs)
    552:
    553:             xc = x_v(row + x_i)
    554:             xt = xc
    555:             xb = xc
    556:             xr = xc
    557:             xl = xc
    558:             xrb = xc
    559:             xlt = xc

    561:             if (i .eq. gxs) then   ! Left side
    562:                xl = left_v(left_i + j - ys + 1)
    563:                xlt = left_v(left_i + j - ys + 2)
    564:             else
    565:                xl = x_v(x_i + row -1 )
    566:             endif

    568:             if (j .eq. gys) then ! bottom side
    569:                xb = bottom_v(bottom_i + i - xs + 1)
    570:                xrb = bottom_v(bottom_i + i - xs + 2)
    571:             else
    572:                xb = x_v(x_i + row - gxm)
    573:             endif

    575:             if (i+1 .eq. gxs + gxm) then !right side
    576:                xr = right_v(right_i + j - ys + 1)
    577:                xrb = right_v(right_i + j - ys)
    578:             else
    579:                xr = x_v(x_i + row + 1)
    580:             endif

    582:             if (j+1 .eq. gym+gys) then !top side
    583:                xt = top_v(top_i +i - xs + 1)
    584:                xlt = top_v(top_i + i - xs)
    585:             else
    586:                xt = x_v(x_i + row + gxm)
    587:             endif

    589:             if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
    590:                xlt = x_v(x_i + row - 1 + gxm)
    591:             endif
    592:
    593:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
    594:                xrb = x_v(x_i + row + 1 - gxm)
    595:             endif

    597:             d1 = (xc-xl)*rhx
    598:             d2 = (xc-xr)*rhx
    599:             d3 = (xc-xt)*rhy
    600:             d4 = (xc-xb)*rhy
    601:             d5 = (xrb-xr)*rhy
    602:             d6 = (xrb-xb)*rhx
    603:             d7 = (xlt-xl)*rhy
    604:             d8 = (xlt-xt)*rhx
    605:
    606:             f1 = sqrt( 1.0d0 + d1*d1 + d7*d7)
    607:             f2 = sqrt( 1.0d0 + d1*d1 + d4*d4)
    608:             f3 = sqrt( 1.0d0 + d3*d3 + d8*d8)
    609:             f4 = sqrt( 1.0d0 + d3*d3 + d2*d2)
    610:             f5 = sqrt( 1.0d0 + d2*d2 + d5*d5)
    611:             f6 = sqrt( 1.0d0 + d4*d4 + d6*d6)
    612:
    613:
    614:             hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+                 &
    615:      &              (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)

    617:             hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+                 &
    618:      &            (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)

    620:             ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+                 &
    621:      &                (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)

    623:             hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+                 &
    624:      &              (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
    625:
    626:             hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
    627:             htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
    628:
    629:             hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) +                         &
    630:      &              hxdhy*(1.0+d8*d8)/(f3*f3*f3) +                      &
    631:      &              hydhx*(1.0+d5*d5)/(f5*f5*f5) +                      &
    632:      &              hxdhy*(1.0+d6*d6)/(f6*f6*f6) +                      &
    633:      &              (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-               &
    634:      &              2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+          &
    635:      &              hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
    636:
    637:             hl = hl * 0.5
    638:             hr = hr * 0.5
    639:             ht = ht * 0.5
    640:             hb = hb * 0.5
    641:             hbr = hbr * 0.5
    642:             htl = htl * 0.5
    643:             hc = hc * 0.5

    645:             k = 0

    647:             if (j .gt. 0) then
    648:                v(k) = hb
    649:                col(k) = row - gxm
    650:                k=k+1
    651:             endif

    653:             if ((j .gt. 0) .and. (i .lt. mx-1)) then
    654:                v(k) = hbr
    655:                col(k) = row-gxm+1
    656:                k=k+1
    657:             endif

    659:             if (i .gt. 0) then
    660:                v(k) = hl
    661:                col(k) = row - 1
    662:                k = k+1
    663:             endif

    665:             v(k) = hc
    666:             col(k) = row
    667:             k=k+1

    669:             if (i .lt. mx-1) then
    670:                v(k) = hr
    671:                col(k) = row + 1
    672:                k=k+1
    673:             endif

    675:             if ((i .gt. 0) .and. (j .lt. my-1)) then
    676:                v(k) = htl
    677:                col(k) = row + gxm - 1
    678:                k=k+1
    679:             endif

    681:             if (j .lt. my-1) then
    682:                v(k) = ht
    683:                col(k) = row + gxm
    684:                k=k+1
    685:             endif

    687: ! Set matrix values using local numbering, defined earlier in main routine
    688:             call MatSetValuesLocal(Hessian,1,row,k,col,v,INSERT_VALUES,       &
    689:      &                              info)

    691:

    693:          enddo
    694:       enddo
    695:
    696: ! restore vectors
    697:       call VecRestoreArray(localX,x_v,x_i,info)
    698:       call VecRestoreArray(Left,left_v,left_i,info)
    699:       call VecRestoreArray(Right,right_v,right_i,info)
    700:       call VecRestoreArray(Top,top_v,top_i,info)
    701:       call VecRestoreArray(Bottom,bottom_v,bottom_i,info)

    704: ! Assemble the matrix
    705:       call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,info)
    706:       call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,info)

    708:       call PetscLogFlops(199*xm*ym,info)

    710:       return
    711:       end
    712:
    713:

    717: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h

    719: ! ----------------------------------------------------------------------------
    720: !
    721: !/*
    722: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
    723: !
    724: !
    725: !*/

    727:       subroutine MSA_BoundaryConditions(info)
    728:       implicit none

    730: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
    731: #include "plate2f.h"

    733:       integer i,j,k,limit,info,maxits
    734:       integer          xs, xm, gxs, gxm, ys, ym, gys, gym
    735:       integer bsize, lsize, tsize, rsize
    736:       PetscScalar           one, two, three, tol
    737:       PetscScalar           scl,fnorm,det,xt,yt,hx,hy
    738:       PetscScalar           u1,u2,nf1,nf2,njac11,njac12,njac21,njac22
    739:       PetscScalar           b, t, l, r
    740:       PetscScalar           boundary_v(0:1)
    741:       PetscOffset           boundary_i
    742:       logical exitloop
    743:       TaoTruth flg

    745:       limit=0
    746:       maxits = 5
    747:       tol=1e-10
    748:       b=-0.5d0
    749:       t= 0.5d0
    750:       l=-0.5d0
    751:       r= 0.5d0
    752:       xt=0
    753:       yt=0
    754:       one=1.0d0
    755:       two=2.0d0
    756:       three=3.0d0

    759:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    760:      &                  PETSC_NULL_INTEGER,info)
    761:       call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,              &
    762:      &                       gxm,gym,PETSC_NULL_INTEGER,info)

    764:       bsize = xm + 2
    765:       lsize = ym + 2
    766:       rsize = ym + 2
    767:       tsize = xm + 2
    768:

    770:       call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,info)
    771:       call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,info)
    772:       call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,info)
    773:       call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,info)

    775:       hx= (r-l)/(mx+1)
    776:       hy= (t-b)/(my+1)

    778:       do j=0,3
    779:
    780:          if (j.eq.0) then
    781:             yt=b
    782:             xt=l+hx*xs
    783:             limit=bsize
    784:             call VecGetArray(Bottom,boundary_v,boundary_i,info)
    785:

    787:          elseif (j.eq.1) then
    788:             yt=t
    789:             xt=l+hx*xs
    790:             limit=tsize
    791:             call VecGetArray(Top,boundary_v,boundary_i,info)

    793:          elseif (j.eq.2) then
    794:             yt=b+hy*ys
    795:             xt=l
    796:             limit=lsize
    797:             call VecGetArray(Left,boundary_v,boundary_i,info)

    799:          elseif (j.eq.3) then
    800:             yt=b+hy*ys
    801:             xt=r
    802:             limit=rsize
    803:             call VecGetArray(Right,boundary_v,boundary_i,info)
    804:          endif
    805:

    807:          do i=0,limit-1
    808:
    809:             u1=xt
    810:             u2=-yt
    811:             k = 0
    812:             exitloop = .false.
    813:             do while (k .lt. maxits .and. (.not. exitloop) )

    815:                nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
    816:                nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
    817:                fnorm=sqrt(nf1*nf1+nf2*nf2)
    818:                if (fnorm .gt. tol) then
    819:                   njac11=one+u2*u2-u1*u1
    820:                   njac12=two*u1*u2
    821:                   njac21=-two*u1*u2
    822:                   njac22=-one - u1*u1 + u2*u2
    823:                   det = njac11*njac22-njac21*njac12
    824:                   u1 = u1-(njac22*nf1-njac12*nf2)/det
    825:                   u2 = u2-(njac11*nf2-njac21*nf1)/det
    826:                else
    827:                   exitloop = .true.
    828:                endif
    829:                k=k+1
    830:             enddo

    832:             boundary_v(i + boundary_i) = u1*u1-u2*u2
    833:             if ((j .eq. 0) .or. (j .eq. 1)) then
    834:                xt = xt + hx
    835:             else
    836:                yt = yt + hy
    837:             endif

    839:          enddo
    840:

    842:          if (j.eq.0) then
    843:             call VecRestoreArray(Bottom,boundary_v,boundary_i,info)
    844:          elseif (j.eq.1) then
    845:             call VecRestoreArray(Top,boundary_v,boundary_i,info)
    846:          elseif (j.eq.2) then
    847:             call VecRestoreArray(Left,boundary_v,boundary_i,info)
    848:          elseif (j.eq.3) then
    849:             call VecRestoreArray(Right,boundary_v,boundary_i,info)
    850:          endif
    851:
    852:       enddo

    855: ! Scale the boundary if desired
    856:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom",            &
    857:      &                         scl,flg,info)
    858:       if (flg .eq. PETSC_TRUE) then
    859:          call VecScale(scl,Bottom,info)
    860:       endif

    862:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top",               &
    863:      &                         scl,flg,info)
    864:       if (flg .eq. PETSC_TRUE) then
    865:          call VecScale(scl,Top,info)
    866:       endif

    868:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right",             &
    869:      &                         scl,flg,info)
    870:       if (flg .eq. PETSC_TRUE) then
    871:          call VecScale(scl,Right,info)
    872:       endif

    874:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left",              &
    875:      &                         scl,flg,info)
    876:       if (flg .eq. PETSC_TRUE) then
    877:          call VecScale(scl,Left,info)
    878:       endif
    879:
    880:
    881:       return
    882:       end

    884: ! ----------------------------------------------------------------------------
    885: !
    886: !/*
    887: !     MSA_Plate - Calculates an obstacle for surface to stretch over
    888: !
    889: !     Output Parameter:
    890: !.    xl - lower bound vector
    891: !.    xu - upper bound vector
    892: !
    893: !*/

    895:       subroutine MSA_Plate(xl,xu,info)
    896:       implicit none

    898: ! mx,my,bmx,bmy,da,bheight defined in plate2f.h
    899: #include "plate2f.h"
    900:       Vec              xl,xu
    901:       integer          i,j,row,info
    902:       integer          xs, xm, ys, ym
    903:       PetscScalar      lb,ub

    905: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    906: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
    907: ! will return an array of doubles referenced by x_array offset by x_index.
    908: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    909: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    910:       PetscScalar      xl_v(0:1)
    911:       PetscOffset      xl_i

    914:       lb = -1.0d300
    915:       ub = 1.0d300

    917:       if (bmy .lt. 0) bmy = 0
    918:       if (bmy .gt. my) bmy = my
    919:       if (bmx .lt. 0) bmx = 0
    920:       if (bmx .gt. mx) bmx = mx
    921:

    923:       call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
    924:      &             PETSC_NULL_INTEGER,info)

    926:       call VecSet(xl,lb,info)
    927:       call VecSet(xu,ub,info)

    929:       call VecGetArray(xl,xl_v,xl_i,info)
    930:

    932:       do i=xs,xs+xm-1

    934:          do j=ys,ys+ym-1
    935:
    936:             row=(j-ys)*xm + (i-xs)

    938:             if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and.           &
    939:      &          j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
    940:                xl_v(xl_i+row) = bheight

    942:             endif

    944:          enddo
    945:       enddo

    948:       call VecRestoreArray(xl,xl_v,xl_i,info)
    949:
    950:       return
    951:       end

    955:
    956:
    957: ! ----------------------------------------------------------------------------
    958: !
    959: !/*
    960: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
    961: !
    962: !     Output Parameter:
    963: !.    X - vector for initial guess
    964: !
    965: !*/

    967:       subroutine MSA_InitialPoint(X, info)
    968:       implicit none

    970: ! mx,my,localX,da,Top,Left,Bottom,Right defined in plate2f.h
    971: #include "plate2f.h"
    972:       Vec               X

    974:       integer           start,i,j,info
    975:       integer           row,xs,xm,gxs,gxm,ys,ym,gys,gym
    976:       PetscScalar       zero, np5

    978: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    979: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (integer) x_index, info)
    980: ! will return an array of doubles referenced by x_array offset by x_index.
    981: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    982: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    983:       PetscScalar       left_v(0:1), right_v(0:1), top_v(0:1)
    984:       PetscScalar       x_v(0:1), bottom_v(0:1)
    985:       PetscOffset       left_i, right_i, top_i, bottom_i, x_i
    986:       PetscTruth        flg
    987:       PetscRandom       rctx
    988:
    989:       zero = 0.0d0
    990:       np5 = -0.5d0

    992:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start",            &
    993:      &                        start,flg,info)

    995:       if ((flg .eq. PETSC_TRUE) .and. (start .eq. 0)) then  ! the zero vector is reasonable
    996:          call VecSet(X,zero,info)

    998:       elseif ((flg .eq. PETSC_TRUE) .and. (start .gt. 0)) then  ! random start -0.5 < xi < 0.5
    999:          call PetscRandomCreate(MPI_COMM_WORLD,RANDOM_DEFAULT,rctx,info)
    1000:          do i=0,start-1
    1001:             call VecSetRandom(X,rctx,info)
    1002:          enddo

    1004:          call PetscRandomDestroy(rctx,info)
    1005:          call VecShift(X,np5,info)

    1007:       else   ! average of boundary conditions
    1008:
    1009: !        Get Local mesh boundaries
    1010:          call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,             &
    1011:      &                     PETSC_NULL_INTEGER,info)
    1012:          call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,    &
    1013:      &                     PETSC_NULL_INTEGER,info)

    1017: !        Get pointers to vector data
    1018:          call VecGetArray(Top,top_v,top_i,info)
    1019:          call VecGetArray(Bottom,bottom_v,bottom_i,info)
    1020:          call VecGetArray(Left,left_v,left_i,info)
    1021:          call VecGetArray(Right,right_v,right_i,info)
    1022:
    1023:          call VecGetArray(localX,x_v,x_i,info)
    1024:
    1025: !        Perform local computations
    1026:          do  j=ys,ys+ym-1
    1027:             do i=xs,xs+xm-1
    1028:                row = (j-gys)*gxm  + (i-gxs)
    1029:                x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my        &
    1030:      &             + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) +                  &
    1031:      &              (i+1)*left_v(left_i+j-ys+1)/mx       +                  &
    1032:      &              (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
    1033:             enddo
    1034:          enddo

    1036: !        Restore vectors
    1037:          call VecRestoreArray(localX,x_v,x_i,info)

    1039:          call VecRestoreArray(Left,left_v,left_i,info)
    1040:          call VecRestoreArray(Top,top_v,top_i,info)
    1041:          call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
    1042:          call VecRestoreArray(Right,right_v,right_i,info)

    1044:          call DALocalToGlobal(da,localX,INSERT_VALUES,X,info)

    1046:       endif

    1048:       return
    1049:       end