Actual source code: plate2f.F90

      1: !  Program usage: mpiexec -n <proc> plate2f [all TAO options]
      2: !
      3: !  This example demonstrates use of the TAO package to solve a bound constrained
      4: !  minimization problem.  This example is based on a problem from the
      5: !  MINPACK-2 test suite.  Given a rectangular 2-D domain and boundary values
      6: !  along the edges of the domain, the objective is to find the surface
      7: !  with the minimal area that satisfies the boundary conditions.
      8: !  The command line options are:
      9: !    -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
     10: !    -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
     11: !    -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
     12: !    -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
     13: !    -bheight <ht>, where <ht> = height of the plate
     14: !

     16: module plate2fmodule
     17: #include "petsc/finclude/petscdmda.h"
     18: #include "petsc/finclude/petsctao.h"
     19:   use petscdmda
     20:   use petsctao

     22:   Vec localX, localV
     23:   Vec Top, Left
     24:   Vec Right, Bottom
     25:   DM dm
     26:   PetscReal bheight
     27:   PetscInt bmx, bmy
     28:   PetscInt mx, my, Nx, Ny, N
     29: end module

     31: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     32: !                   Variable declarations
     33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     34: !
     35: !  Variables:
     36: !    (common from plate2f.h):
     37: !    Nx, Ny           number of processors in x- and y- directions
     38: !    mx, my           number of grid points in x,y directions
     39: !    N    global dimension of vector
     40: use plate2fmodule
     41: implicit none

     43: PetscErrorCode ierr          ! used to check for functions returning nonzeros
     44: Vec x             ! solution vector
     45: PetscInt m             ! number of local elements in vector
     46: Tao ta           ! Tao solver context
     47: Mat H             ! Hessian matrix
     48: ISLocalToGlobalMapping isltog  ! local to global mapping object
     49: PetscBool flg
     50: PetscInt i1, i3, i7

     52: external FormFunctionGradient
     53: external FormHessian
     54: external MSA_BoundaryConditions
     55: external MSA_Plate
     56: external MSA_InitialPoint
     57: ! Initialize Tao

     59: i1 = 1
     60: i3 = 3
     61: i7 = 7

     63: PetscCallA(PetscInitialize(ierr))

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

     70: ! Check for any command line arguments that override defaults

     72: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
     73: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))

     75: bmx = mx/2
     76: bmy = my/2

     78: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmx', bmx, flg, ierr))
     79: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmy', bmy, flg, ierr))
     80: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bheight', bheight, flg, ierr))

     82: ! Calculate any derived values from parameters
     83: N = mx*my

     85: ! Let PETSc determine the dimensions of the local vectors
     86: Nx = PETSC_DECIDE
     87: NY = PETSC_DECIDE

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

     93: 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))
     94: PetscCallA(DMSetFromOptions(dm, ierr))
     95: PetscCallA(DMSetUp(dm, ierr))

     97: ! Extract global and local vectors from DM; The local vectors are
     98: ! used solely as work space for the evaluation of the function,
     99: ! gradient, and Hessian.  Duplicate for remaining vectors that are
    100: ! the same types.

    102: PetscCallA(DMCreateGlobalVector(dm, x, ierr))
    103: PetscCallA(DMCreateLocalVector(dm, localX, ierr))
    104: PetscCallA(VecDuplicate(localX, localV, ierr))

    106: ! Create a matrix data structure to store the Hessian.
    107: ! Here we (optionally) also associate the local numbering scheme
    108: ! with the matrix so that later we can use local indices for matrix
    109: ! assembly

    111: PetscCallA(VecGetLocalSize(x, m, ierr))
    112: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, m, m, N, N, i7, PETSC_NULL_INTEGER_ARRAY, i3, PETSC_NULL_INTEGER_ARRAY, H, ierr))

    114: PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))
    115: PetscCallA(DMGetLocalToGlobalMapping(dm, isltog, ierr))
    116: PetscCallA(MatSetLocalToGlobalMapping(H, isltog, isltog, ierr))

    118: ! The Tao code begins here
    119: ! Create TAO solver and set desired solution method.
    120: ! This problems uses bounded variables, so the
    121: ! method must either be 'tao_tron' or 'tao_blmvm'

    123: PetscCallA(TaoCreate(PETSC_COMM_WORLD, ta, ierr))
    124: PetscCallA(TaoSetType(ta, TAOBLMVM, ierr))

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

    128: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))

    130: PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))

    132: ! Set Variable bounds
    133: PetscCallA(MSA_BoundaryConditions(ierr))
    134: PetscCallA(TaoSetVariableBoundsRoutine(ta, MSA_Plate, 0, ierr))

    136: ! Set the initial solution guess
    137: PetscCallA(MSA_InitialPoint(x, ierr))
    138: PetscCallA(TaoSetSolution(ta, x, ierr))

    140: ! Check for any tao command line options
    141: PetscCallA(TaoSetFromOptions(ta, ierr))

    143: ! Solve the application
    144: PetscCallA(TaoSolve(ta, ierr))

    146: ! Free TAO data structures
    147: PetscCallA(TaoDestroy(ta, ierr))

    149: ! Free PETSc data structures
    150: PetscCallA(VecDestroy(x, ierr))
    151: PetscCallA(VecDestroy(Top, ierr))
    152: PetscCallA(VecDestroy(Bottom, ierr))
    153: PetscCallA(VecDestroy(Left, ierr))
    154: PetscCallA(VecDestroy(Right, ierr))
    155: PetscCallA(MatDestroy(H, ierr))
    156: PetscCallA(VecDestroy(localX, ierr))
    157: PetscCallA(VecDestroy(localV, ierr))
    158: PetscCallA(DMDestroy(dm, ierr))

    160: ! Finalize TAO

    162: PetscCallA(PetscFinalize(ierr))

    164: end

    166: ! ---------------------------------------------------------------------
    167: !
    168: !  FormFunctionGradient - Evaluates function f(X).
    169: !
    170: !  Input Parameters:
    171: !  tao   - the Tao context
    172: !  X     - the input vector
    173: !  dummy - optional user-defined context, as set by TaoSetFunction()
    174: !          (not used here)
    175: !
    176: !  Output Parameters:
    177: !  fcn     - the newly evaluated function
    178: !  G       - the gradient vector
    179: !  info  - error code
    180: !

    182: subroutine FormFunctionGradient(ta, X, fcn, G, dummy, ierr)
    183:   use plate2fmodule
    184:   implicit none

    186: ! Input/output variables

    188:   Tao ta
    189:   PetscReal fcn
    190:   Vec X, G
    191:   PetscErrorCode ierr
    192:   PetscInt dummy

    194:   PetscInt i, j, row
    195:   PetscInt xs, xm
    196:   PetscInt gxs, gxm
    197:   PetscInt ys, ym
    198:   PetscInt gys, gym
    199:   PetscReal ft, zero, hx, hy, hydhx, hxdhy
    200:   PetscReal area, rhx, rhy
    201:   PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3
    202:   PetscReal d4, d5, d6, d7, d8
    203:   PetscReal df1dxc, df2dxc, df3dxc, df4dxc
    204:   PetscReal df5dxc, df6dxc
    205:   PetscReal xc, xl, xr, xt, xb, xlt, xrb

    207:   PetscReal, pointer :: g_v(:), x_v(:)
    208:   PetscReal, pointer :: top_v(:), left_v(:)
    209:   PetscReal, pointer :: right_v(:), bottom_v(:)

    211:   ft = 0.0
    212:   zero = 0.0
    213:   hx = 1.0/real(mx + 1)
    214:   hy = 1.0/real(my + 1)
    215:   hydhx = hy/hx
    216:   hxdhy = hx/hy
    217:   area = 0.5*hx*hy
    218:   rhx = real(mx) + 1.0
    219:   rhy = real(my) + 1.0

    221: ! Get local mesh boundaries
    222:   PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
    223:   PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))

    225: ! Scatter ghost points to local vector
    226:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
    227:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))

    229: ! Initialize the vector to zero
    230:   PetscCall(VecSet(localV, zero, ierr))

    232: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
    233:   PetscCall(VecGetArray(localX, x_v, ierr))
    234:   PetscCall(VecGetArray(localV, g_v, ierr))
    235:   PetscCall(VecGetArray(Top, top_v, ierr))
    236:   PetscCall(VecGetArray(Bottom, bottom_v, ierr))
    237:   PetscCall(VecGetArray(Left, left_v, ierr))
    238:   PetscCall(VecGetArray(Right, right_v, ierr))

    240: ! Compute function over the locally owned part of the mesh
    241:   do j = ys, ys + ym - 1
    242:     do i = xs, xs + xm - 1
    243:       row = (j - gys)*gxm + (i - gxs)
    244:       xc = x_v(1 + row)
    245:       xt = xc
    246:       xb = xc
    247:       xr = xc
    248:       xl = xc
    249:       xrb = xc
    250:       xlt = xc

    252:       if (i == 0) then !left side
    253:         xl = left_v(1 + j - ys + 1)
    254:         xlt = left_v(1 + j - ys + 2)
    255:       else
    256:         xl = x_v(1 + row - 1)
    257:       end if

    259:       if (j == 0) then !bottom side
    260:         xb = bottom_v(1 + i - xs + 1)
    261:         xrb = bottom_v(1 + i - xs + 2)
    262:       else
    263:         xb = x_v(1 + row - gxm)
    264:       end if

    266:       if (i + 1 == gxs + gxm) then !right side
    267:         xr = right_v(1 + j - ys + 1)
    268:         xrb = right_v(1 + j - ys)
    269:       else
    270:         xr = x_v(1 + row + 1)
    271:       end if

    273:       if (j + 1 == gys + gym) then !top side
    274:         xt = top_v(1 + i - xs + 1)
    275:         xlt = top_v(1 + i - xs)
    276:       else
    277:         xt = x_v(1 + row + gxm)
    278:       end if

    280:       if ((i > gxs) .and. (j + 1 < gys + gym)) then
    281:         xlt = x_v(1 + row - 1 + gxm)
    282:       end if

    284:       if ((j > gys) .and. (i + 1 < gxs + gxm)) then
    285:         xrb = x_v(1 + row + 1 - gxm)
    286:       end if

    288:       d1 = xc - xl
    289:       d2 = xc - xr
    290:       d3 = xc - xt
    291:       d4 = xc - xb
    292:       d5 = xr - xrb
    293:       d6 = xrb - xb
    294:       d7 = xlt - xl
    295:       d8 = xt - xlt

    297:       df1dxc = d1*hydhx
    298:       df2dxc = d1*hydhx + d4*hxdhy
    299:       df3dxc = d3*hxdhy
    300:       df4dxc = d2*hydhx + d3*hxdhy
    301:       df5dxc = d2*hydhx
    302:       df6dxc = d4*hxdhy

    304:       d1 = d1*rhx
    305:       d2 = d2*rhx
    306:       d3 = d3*rhy
    307:       d4 = d4*rhy
    308:       d5 = d5*rhy
    309:       d6 = d6*rhx
    310:       d7 = d7*rhy
    311:       d8 = d8*rhx

    313:       f1 = sqrt(1.0 + d1*d1 + d7*d7)
    314:       f2 = sqrt(1.0 + d1*d1 + d4*d4)
    315:       f3 = sqrt(1.0 + d3*d3 + d8*d8)
    316:       f4 = sqrt(1.0 + d3*d3 + d2*d2)
    317:       f5 = sqrt(1.0 + d2*d2 + d5*d5)
    318:       f6 = sqrt(1.0 + d4*d4 + d6*d6)

    320:       ft = ft + f2 + f4

    322:       df1dxc = df1dxc/f1
    323:       df2dxc = df2dxc/f2
    324:       df3dxc = df3dxc/f3
    325:       df4dxc = df4dxc/f4
    326:       df5dxc = df5dxc/f5
    327:       df6dxc = df6dxc/f6

    329:       g_v(1 + row) = 0.5*(df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc)
    330:     end do
    331:   end do

    333: ! Compute triangular areas along the border of the domain.
    334:   if (xs == 0) then  ! left side
    335:     do j = ys, ys + ym - 1
    336:       d3 = (left_v(1 + j - ys + 1) - left_v(1 + j - ys + 2))*rhy
    337:       d2 = (left_v(1 + j - ys + 1) - x_v(1 + (j - gys)*gxm))*rhx
    338:       ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
    339:     end do
    340:   end if

    342:   if (ys == 0) then !bottom side
    343:     do i = xs, xs + xm - 1
    344:       d2 = (bottom_v(1 + i + 1 - xs) - bottom_v(1 + i - xs + 2))*rhx
    345:       d3 = (bottom_v(1 + i - xs + 1) - x_v(1 + i - gxs))*rhy
    346:       ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
    347:     end do
    348:   end if

    350:   if (xs + xm == mx) then ! right side
    351:     do j = ys, ys + ym - 1
    352:       d1 = (x_v(1 + (j + 1 - gys)*gxm - 1) - right_v(1 + j - ys + 1))*rhx
    353:       d4 = (right_v(1 + j - ys) - right_v(1 + j - ys + 1))*rhy
    354:       ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
    355:     end do
    356:   end if

    358:   if (ys + ym == my) then
    359:     do i = xs, xs + xm - 1
    360:       d1 = (x_v(1 + (gym - 1)*gxm + i - gxs) - top_v(1 + i - xs + 1))*rhy
    361:       d4 = (top_v(1 + i - xs + 1) - top_v(1 + i - xs))*rhx
    362:       ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
    363:     end do
    364:   end if

    366:   if ((ys == 0) .and. (xs == 0)) then
    367:     d1 = (left_v(1 + 0) - left_v(1 + 1))*rhy
    368:     d2 = (bottom_v(1 + 0) - bottom_v(1 + 1))*rhx
    369:     ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
    370:   end if

    372:   if ((ys + ym == my) .and. (xs + xm == mx)) then
    373:     d1 = (right_v(1 + ym + 1) - right_v(1 + ym))*rhy
    374:     d2 = (top_v(1 + xm + 1) - top_v(1 + xm))*rhx
    375:     ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
    376:   end if

    378:   ft = ft*area
    379:   PetscCallMPI(MPI_Allreduce(ft, fcn, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD, ierr))

    381: ! Restore vectors
    382:   PetscCall(VecRestoreArray(localX, x_v, ierr))
    383:   PetscCall(VecRestoreArray(localV, g_v, ierr))
    384:   PetscCall(VecRestoreArray(Left, left_v, ierr))
    385:   PetscCall(VecRestoreArray(Top, top_v, ierr))
    386:   PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
    387:   PetscCall(VecRestoreArray(Right, right_v, ierr))

    389: ! Scatter values to global vector
    390:   PetscCall(DMLocalToGlobalBegin(dm, localV, INSERT_VALUES, G, ierr))
    391:   PetscCall(DMLocalToGlobalEnd(dm, localV, INSERT_VALUES, G, ierr))

    393:   PetscCall(PetscLogFlops(70.0d0*xm*ym, ierr))

    395: end  !FormFunctionGradient

    397: ! ----------------------------------------------------------------------------
    398: !
    399: !
    400: !   FormHessian - Evaluates Hessian matrix.
    401: !
    402: !   Input Parameters:
    403: !.  tao  - the Tao context
    404: !.  X    - input vector
    405: !.  dummy  - not used
    406: !
    407: !   Output Parameters:
    408: !.  Hessian    - Hessian matrix
    409: !.  Hpc    - optionally different matrix used to construct the preconditioner
    410: !
    411: !   Notes:
    412: !   Due to mesh point reordering with DMs, we must always work
    413: !   with the local mesh points, and then transform them to the new
    414: !   global numbering with the local-to-global mapping.  We cannot work
    415: !   directly with the global numbers for the original uniprocessor mesh!
    416: !
    417: !      MatSetValuesLocal(), using the local ordering (including
    418: !         ghost points!)
    419: !         - Then set matrix entries using the local ordering
    420: !           by calling MatSetValuesLocal()

    422: subroutine FormHessian(ta, X, Hessian, Hpc, dummy, ierr)
    423:   use plate2fmodule
    424:   implicit none

    426:   Tao ta
    427:   Vec X
    428:   Mat Hessian, Hpc
    429:   PetscInt dummy
    430:   PetscErrorCode ierr

    432:   PetscInt i, j, k, row
    433:   PetscInt xs, xm, gxs, gxm
    434:   PetscInt ys, ym, gys, gym
    435:   PetscInt col(0:6)
    436:   PetscReal hx, hy, hydhx, hxdhy, rhx, rhy
    437:   PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3
    438:   PetscReal d4, d5, d6, d7, d8
    439:   PetscReal xc, xl, xr, xt, xb, xlt, xrb
    440:   PetscReal hl, hr, ht, hb, hc, htl, hbr

    442:   PetscReal, pointer ::  right_v(:), left_v(:)
    443:   PetscReal, pointer ::  bottom_v(:), top_v(:)
    444:   PetscReal, pointer ::  x_v(:)
    445:   PetscReal v(0:6)
    446:   PetscBool assembled
    447:   PetscInt i1

    449:   i1 = 1

    451: ! Set various matrix options
    452:   PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE, ierr))

    454: ! Get local mesh boundaries
    455:   PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
    456:   PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))

    458: ! Scatter ghost points to local vectors
    459:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
    460:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))

    462: ! Get pointers to vector data (see note on Fortran arrays above)
    463:   PetscCall(VecGetArray(localX, x_v, ierr))
    464:   PetscCall(VecGetArray(Top, top_v, ierr))
    465:   PetscCall(VecGetArray(Bottom, bottom_v, ierr))
    466:   PetscCall(VecGetArray(Left, left_v, ierr))
    467:   PetscCall(VecGetArray(Right, right_v, ierr))

    469: ! Initialize matrix entries to zero
    470:   PetscCall(MatAssembled(Hessian, assembled, ierr))
    471:   if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian, ierr))

    473:   rhx = real(mx) + 1.0
    474:   rhy = real(my) + 1.0
    475:   hx = 1.0/rhx
    476:   hy = 1.0/rhy
    477:   hydhx = hy/hx
    478:   hxdhy = hx/hy
    479: ! compute Hessian over the locally owned part of the mesh

    481:   do i = xs, xs + xm - 1
    482:     do j = ys, ys + ym - 1
    483:       row = (j - gys)*gxm + (i - gxs)

    485:       xc = x_v(1 + row)
    486:       xt = xc
    487:       xb = xc
    488:       xr = xc
    489:       xl = xc
    490:       xrb = xc
    491:       xlt = xc

    493:       if (i == gxs) then   ! Left side
    494:         xl = left_v(1 + j - ys + 1)
    495:         xlt = left_v(1 + j - ys + 2)
    496:       else
    497:         xl = x_v(1 + row - 1)
    498:       end if

    500:       if (j == gys) then ! bottom side
    501:         xb = bottom_v(1 + i - xs + 1)
    502:         xrb = bottom_v(1 + i - xs + 2)
    503:       else
    504:         xb = x_v(1 + row - gxm)
    505:       end if

    507:       if (i + 1 == gxs + gxm) then !right side
    508:         xr = right_v(1 + j - ys + 1)
    509:         xrb = right_v(1 + j - ys)
    510:       else
    511:         xr = x_v(1 + row + 1)
    512:       end if

    514:       if (j + 1 == gym + gys) then !top side
    515:         xt = top_v(1 + i - xs + 1)
    516:         xlt = top_v(1 + i - xs)
    517:       else
    518:         xt = x_v(1 + row + gxm)
    519:       end if

    521:       if ((i > gxs) .and. (j + 1 < gys + gym)) then
    522:         xlt = x_v(1 + row - 1 + gxm)
    523:       end if

    525:       if ((i + 1 < gxs + gxm) .and. (j > gys)) then
    526:         xrb = x_v(1 + row + 1 - gxm)
    527:       end if

    529:       d1 = (xc - xl)*rhx
    530:       d2 = (xc - xr)*rhx
    531:       d3 = (xc - xt)*rhy
    532:       d4 = (xc - xb)*rhy
    533:       d5 = (xrb - xr)*rhy
    534:       d6 = (xrb - xb)*rhx
    535:       d7 = (xlt - xl)*rhy
    536:       d8 = (xlt - xt)*rhx

    538:       f1 = sqrt(1.0 + d1*d1 + d7*d7)
    539:       f2 = sqrt(1.0 + d1*d1 + d4*d4)
    540:       f3 = sqrt(1.0 + d3*d3 + d8*d8)
    541:       f4 = sqrt(1.0 + d3*d3 + d2*d2)
    542:       f5 = sqrt(1.0 + d2*d2 + d5*d5)
    543:       f6 = sqrt(1.0 + d4*d4 + d6*d6)

    545:       hl = (-hydhx*(1.0 + d7*d7) + d1*d7)/(f1*f1*f1) + (-hydhx*(1.0 + d4*d4) + d1*d4)/(f2*f2*f2)

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

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

    551:       hb = (-hxdhy*(1.0 + d6*d6) + d4*d6)/(f6*f6*f6) + (-hxdhy*(1.0 + d1*d1) + d1*d4)/(f2*f2*f2)

    553:       hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
    554:       htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)

    556:       hc = hydhx*(1.0 + d7*d7)/(f1*f1*f1) + hxdhy*(1.0 + d8*d8)/(f3*f3*f3) + hydhx*(1.0 + d5*d5)/(f5*f5*f5) +                      &
    557: &           hxdhy*(1.0 + d6*d6)/(f6*f6*f6) + (hxdhy*(1.0 + d1*d1) + hydhx*(1.0 + d4*d4) - 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0 + d2*d2) +   &
    558: &           hydhx*(1.0 + d3*d3) - 2*d2*d3)/(f4*f4*f4)

    560:       hl = hl*0.5
    561:       hr = hr*0.5
    562:       ht = ht*0.5
    563:       hb = hb*0.5
    564:       hbr = hbr*0.5
    565:       htl = htl*0.5
    566:       hc = hc*0.5

    568:       k = 0

    570:       if (j > 0) then
    571:         v(k) = hb
    572:         col(k) = row - gxm
    573:         k = k + 1
    574:       end if

    576:       if ((j > 0) .and. (i < mx - 1)) then
    577:         v(k) = hbr
    578:         col(k) = row - gxm + 1
    579:         k = k + 1
    580:       end if

    582:       if (i > 0) then
    583:         v(k) = hl
    584:         col(k) = row - 1
    585:         k = k + 1
    586:       end if

    588:       v(k) = hc
    589:       col(k) = row
    590:       k = k + 1

    592:       if (i < mx - 1) then
    593:         v(k) = hr
    594:         col(k) = row + 1
    595:         k = k + 1
    596:       end if

    598:       if ((i > 0) .and. (j < my - 1)) then
    599:         v(k) = htl
    600:         col(k) = row + gxm - 1
    601:         k = k + 1
    602:       end if

    604:       if (j < my - 1) then
    605:         v(k) = ht
    606:         col(k) = row + gxm
    607:         k = k + 1
    608:       end if

    610: ! Set matrix values using local numbering, defined earlier in main routine
    611:       PetscCall(MatSetValuesLocal(Hessian, i1, [row], k, col, v, INSERT_VALUES, ierr))

    613:     end do
    614:   end do

    616: ! restore vectors
    617:   PetscCall(VecRestoreArray(localX, x_v, ierr))
    618:   PetscCall(VecRestoreArray(Left, left_v, ierr))
    619:   PetscCall(VecRestoreArray(Right, right_v, ierr))
    620:   PetscCall(VecRestoreArray(Top, top_v, ierr))
    621:   PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))

    623: ! Assemble the matrix
    624:   PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY, ierr))
    625:   PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY, ierr))

    627:   PetscCall(PetscLogFlops(199.0d0*xm*ym, ierr))

    629: end

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

    633: ! ----------------------------------------------------------------------------
    634: !
    635: !/*
    636: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
    637: !
    638: !
    639: !*/

    641: subroutine MSA_BoundaryConditions(ierr)
    642:   use plate2fmodule
    643:   implicit none

    645:   PetscErrorCode ierr
    646:   PetscInt i, j, k, limit, maxits
    647:   PetscInt xs, xm, gxs, gxm
    648:   PetscInt ys, ym, gys, gym
    649:   PetscInt bsize, lsize
    650:   PetscInt tsize, rsize
    651:   PetscReal one, two, three, tol
    652:   PetscReal scl, fnorm, det, xt
    653:   PetscReal yt, hx, hy, u1, u2, nf1, nf2
    654:   PetscReal njac11, njac12, njac21, njac22
    655:   PetscReal b, t, l, r
    656:   PetscReal, pointer :: boundary_v(:)

    658:   logical exitloop
    659:   PetscBool flg

    661:   limit = 0
    662:   maxits = 5
    663:   tol = 1e-10
    664:   b = -0.5
    665:   t = 0.5
    666:   l = -0.5
    667:   r = 0.5
    668:   xt = 0
    669:   yt = 0
    670:   one = 1.0
    671:   two = 2.0
    672:   three = 3.0

    674:   PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
    675:   PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))

    677:   bsize = xm + 2
    678:   lsize = ym + 2
    679:   rsize = ym + 2
    680:   tsize = xm + 2

    682:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, bsize, PETSC_DECIDE, Bottom, ierr))
    683:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, tsize, PETSC_DECIDE, Top, ierr))
    684:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, lsize, PETSC_DECIDE, Left, ierr))
    685:   PetscCall(VecCreateMPI(PETSC_COMM_WORLD, rsize, PETSC_DECIDE, Right, ierr))

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

    690:   do j = 0, 3

    692:     if (j == 0) then
    693:       yt = b
    694:       xt = l + hx*xs
    695:       limit = bsize
    696:       PetscCall(VecGetArray(Bottom, boundary_v, ierr))

    698:     elseif (j == 1) then
    699:       yt = t
    700:       xt = l + hx*xs
    701:       limit = tsize
    702:       PetscCall(VecGetArray(Top, boundary_v, ierr))

    704:     elseif (j == 2) then
    705:       yt = b + hy*ys
    706:       xt = l
    707:       limit = lsize
    708:       PetscCall(VecGetArray(Left, boundary_v, ierr))

    710:     elseif (j == 3) then
    711:       yt = b + hy*ys
    712:       xt = r
    713:       limit = rsize
    714:       PetscCall(VecGetArray(Right, boundary_v, ierr))
    715:     end if

    717:     do i = 0, limit - 1

    719:       u1 = xt
    720:       u2 = -yt
    721:       k = 0
    722:       exitloop = .false.
    723:       do while (k < maxits .and. (.not. exitloop))

    725:         nf1 = u1 + u1*u2*u2 - u1*u1*u1/three - xt
    726:         nf2 = -u2 - u1*u1*u2 + u2*u2*u2/three - yt
    727:         fnorm = sqrt(nf1*nf1 + nf2*nf2)
    728:         if (fnorm > tol) then
    729:           njac11 = one + u2*u2 - u1*u1
    730:           njac12 = two*u1*u2
    731:           njac21 = -two*u1*u2
    732:           njac22 = -one - u1*u1 + u2*u2
    733:           det = njac11*njac22 - njac21*njac12
    734:           u1 = u1 - (njac22*nf1 - njac12*nf2)/det
    735:           u2 = u2 - (njac11*nf2 - njac21*nf1)/det
    736:         else
    737:           exitloop = .true.
    738:         end if
    739:         k = k + 1
    740:       end do

    742:       boundary_v(1 + i) = u1*u1 - u2*u2
    743:       if ((j == 0) .or. (j == 1)) then
    744:         xt = xt + hx
    745:       else
    746:         yt = yt + hy
    747:       end if

    749:     end do

    751:     if (j == 0) then
    752:       PetscCall(VecRestoreArray(Bottom, boundary_v, ierr))
    753:     elseif (j == 1) then
    754:       PetscCall(VecRestoreArray(Top, boundary_v, ierr))
    755:     elseif (j == 2) then
    756:       PetscCall(VecRestoreArray(Left, boundary_v, ierr))
    757:     elseif (j == 3) then
    758:       PetscCall(VecRestoreArray(Right, boundary_v, ierr))
    759:     end if

    761:   end do

    763: ! Scale the boundary if desired
    764:   PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bottom', scl, flg, ierr))
    765:   if (flg .eqv. PETSC_TRUE) then
    766:     PetscCall(VecScale(Bottom, scl, ierr))
    767:   end if

    769:   PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-top', scl, flg, ierr))
    770:   if (flg .eqv. PETSC_TRUE) then
    771:     PetscCall(VecScale(Top, scl, ierr))
    772:   end if

    774:   PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-right', scl, flg, ierr))
    775:   if (flg .eqv. PETSC_TRUE) then
    776:     PetscCall(VecScale(Right, scl, ierr))
    777:   end if

    779:   PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-left', scl, flg, ierr))
    780:   if (flg .eqv. PETSC_TRUE) then
    781:     PetscCall(VecScale(Left, scl, ierr))
    782:   end if

    784: end

    786: ! ----------------------------------------------------------------------------
    787: !
    788: !/*
    789: !     MSA_Plate - Calculates an obstacle for surface to stretch over
    790: !
    791: !     Output Parameter:
    792: !.    xl - lower bound vector
    793: !.    xu - upper bound vector
    794: !
    795: !*/

    797: subroutine MSA_Plate(ta, xl, xu, dummy, ierr)
    798:   use plate2fmodule
    799:   implicit none

    801:   Tao ta
    802:   Vec xl, xu
    803:   PetscErrorCode ierr
    804:   PetscInt i, j, row
    805:   PetscInt xs, xm, ys, ym
    806:   PetscReal lb, ub
    807:   PetscInt dummy
    808:   PetscReal, pointer :: xl_v(:)

    810:   lb = PETSC_NINFINITY
    811:   ub = PETSC_INFINITY

    813:   if (bmy < 0) bmy = 0
    814:   if (bmy > my) bmy = my
    815:   if (bmx < 0) bmx = 0
    816:   if (bmx > mx) bmx = mx

    818:   PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))

    820:   PetscCall(VecSet(xl, lb, ierr))
    821:   PetscCall(VecSet(xu, ub, ierr))

    823:   PetscCall(VecGetArray(xl, xl_v, ierr))

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

    827:     do j = ys, ys + ym - 1

    829:       row = (j - ys)*xm + (i - xs)

    831:       if (i >= ((mx - bmx)/2) .and. i < (mx - (mx - bmx)/2) .and.           &
    832: &          j >= ((my - bmy)/2) .and. j < (my - (my - bmy)/2)) then
    833:         xl_v(1 + row) = bheight

    835:       end if

    837:     end do
    838:   end do

    840:   PetscCall(VecRestoreArray(xl, xl_v, ierr))

    842: end

    844: ! ----------------------------------------------------------------------------
    845: !
    846: !/*
    847: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
    848: !
    849: !     Output Parameter:
    850: !.    X - vector for initial guess
    851: !
    852: !*/

    854: subroutine MSA_InitialPoint(X, ierr)
    855:   use plate2fmodule
    856:   implicit none

    858:   Vec X
    859:   PetscErrorCode ierr
    860:   PetscInt start, i, j
    861:   PetscInt row
    862:   PetscInt xs, xm, gxs, gxm
    863:   PetscInt ys, ym, gys, gym
    864:   PetscReal zero, np5

    866:   PetscReal, pointer :: left_v(:), right_v(:)
    867:   PetscReal, pointer :: bottom_v(:), top_v(:)
    868:   PetscReal, pointer :: x_v(:)
    869:   PetscBool flg
    870:   PetscRandom rctx

    872:   zero = 0.0
    873:   np5 = -0.5

    875:   PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-start', start, flg, ierr))

    877:   if ((flg .eqv. PETSC_TRUE) .and. (start == 0)) then  ! the zero vector is reasonable
    878:     PetscCall(VecSet(X, zero, ierr))

    880:   elseif ((flg .eqv. PETSC_TRUE) .and. (start > 0)) then  ! random start -0.5 < xi < 0.5
    881:     PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, rctx, ierr))
    882:     do i = 0, start - 1
    883:       PetscCall(VecSetRandom(X, rctx, ierr))
    884:     end do

    886:     PetscCall(PetscRandomDestroy(rctx, ierr))
    887:     PetscCall(VecShift(X, np5, ierr))

    889:   else   ! average of boundary conditions

    891: !        Get Local mesh boundaries
    892:     PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
    893:     PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))

    895: !        Get pointers to vector data
    896:     PetscCall(VecGetArray(Top, top_v, ierr))
    897:     PetscCall(VecGetArray(Bottom, bottom_v, ierr))
    898:     PetscCall(VecGetArray(Left, left_v, ierr))
    899:     PetscCall(VecGetArray(Right, right_v, ierr))

    901:     PetscCall(VecGetArray(localX, x_v, ierr))

    903: !        Perform local computations
    904:     do j = ys, ys + ym - 1
    905:       do i = xs, xs + xm - 1
    906:         row = (j - gys)*gxm + (i - gxs)
    907:         x_v(1 + row) = ((j + 1)*bottom_v(1 + i - xs + 1)/my + (my - j + 1)*top_v(1 + i - xs + 1)/(my + 2) +                  &
    908: &                          (i + 1)*left_v(1 + j - ys + 1)/mx + (mx - i + 1)*right_v(1 + j - ys + 1)/(mx + 2))*0.5
    909:       end do
    910:     end do

    912: !        Restore vectors
    913:     PetscCall(VecRestoreArray(localX, x_v, ierr))

    915:     PetscCall(VecRestoreArray(Left, left_v, ierr))
    916:     PetscCall(VecRestoreArray(Top, top_v, ierr))
    917:     PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
    918:     PetscCall(VecRestoreArray(Right, right_v, ierr))

    920:     PetscCall(DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, X, ierr))
    921:     PetscCall(DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, X, ierr))

    923:   end if

    925: end

    927: !
    928: !/*TEST
    929: !
    930: !   build:
    931: !      requires: !complex
    932: !
    933: !   test:
    934: !      args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
    935: !      filter: sort -b
    936: !      filter_output: sort -b
    937: !      requires: !single
    938: !
    939: !   test:
    940: !      suffix: 2
    941: !      nsize: 2
    942: !      args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
    943: !      filter: sort -b
    944: !      filter_output: sort -b
    945: !      requires: !single
    946: !
    947: !TEST*/