Actual source code: plate2f.F

petsc-master 2016-09-25
Report Typos and Errors
  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: !
 15: !/*T
 16: !   Concepts: TAO^Solving a bound constrained minimization problem
 17: !   Routines: TaoCreate();
 18: !   Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
 19: !   Routines: TaoSetHessianRoutine();
 20: !   Routines: TaoSetVariableBoundsRoutine();
 21: !   Routines: TaoSetInitialVector();
 22: !   Routines: TaoSetFromOptions();
 23: !   Routines: TaoSolve();
 24: !   Routines: TaoDestroy();
 25: !   Processors: n
 26: !T*/



 30:       implicit none

 32: #include "plate2f.h"

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

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


 54:       external FormFunctionGradient
 55:       external FormHessian
 56:       external MSA_BoundaryConditions
 57:       external MSA_Plate
 58:       external MSA_InitialPoint
 59: ! Initialize Tao

 61:       i1=1
 62:       i3=3
 63:       i7=7


 66:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 67:       if (ierr .ne. 0) then
 68:         print*,'Unable to initialize PETSc'
 69:         stop
 70:       endif

 72: ! Specify default dimensions of the problem
 73:       mx = 10
 74:       my = 10
 75:       bheight = 0.1

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

 79:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
 80:      &                        '-mx',mx,flg,ierr)
 81:       CHKERRQ(ierr)
 82:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
 83:      &                        '-my',my,flg,ierr)
 84:       CHKERRQ(ierr)

 86:       bmx = mx/2
 87:       bmy = my/2

 89:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
 90:      &                        '-bmx',bmx,flg,ierr)
 91:       CHKERRQ(ierr)
 92:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,     &
 93:      &                        '-bmy',bmy,flg,ierr)
 94:       CHKERRQ(ierr)
 95:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,    &
 96:      &                         '-bheight',bheight,flg,ierr)
 97:       CHKERRQ(ierr)


100: ! Calculate any derived values from parameters
101:       N = mx*my

103: ! Let Petsc determine the dimensions of the local vectors
104:       Nx = PETSC_DECIDE
105:       NY = PETSC_DECIDE

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

111:       call DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,                    &
112:      &     DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,                              &
113:      &     mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,           &
114:      &     dm,ierr)
115:       CHKERRQ(ierr)
116:       call DMSetFromOptions(dm,ierr)
117:       call DMSetUp(dm,ierr)

119: ! Extract global and local vectors from DM; The local vectors are
120: ! used solely as work space for the evaluation of the function,
121: ! gradient, and Hessian.  Duplicate for remaining vectors that are
122: ! the same types.

124:       call DMCreateGlobalVector(dm,x,ierr)
125:       CHKERRQ(ierr)
126:       call DMCreateLocalVector(dm,localX,ierr)
127:       CHKERRQ(ierr)
128:       call VecDuplicate(localX,localV,ierr)
129:       CHKERRQ(ierr)

131: ! Create a matrix data structure to store the Hessian.
132: ! Here we (optionally) also associate the local numbering scheme
133: ! with the matrix so that later we can use local indices for matrix
134: ! assembly

136:       call VecGetLocalSize(x,m,ierr)
137:       CHKERRQ(ierr)
138:       call MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER,   &
139:      &     i3,PETSC_NULL_INTEGER,H,ierr)
140:       CHKERRQ(ierr)

142:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
143:       CHKERRQ(ierr)
144:       call DMGetLocalToGlobalMapping(dm,isltog,ierr)
145:       CHKERRQ(ierr)
146:       call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
147:       CHKERRQ(ierr)


150: ! The Tao code begins here
151: ! Create TAO solver and set desired solution method.
152: ! This problems uses bounded variables, so the
153: ! method must either be 'tao_tron' or 'tao_blmvm'

155:       call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
156:       CHKERRQ(ierr)
157:       call TaoSetType(tao,TAOBLMVM,ierr)
158:       CHKERRQ(ierr)

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

162:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
163:      &     FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
164:       CHKERRQ(ierr)

166:       call TaoSetHessianRoutine(tao,H,H,FormHessian,                    &
167:      &     PETSC_NULL_OBJECT, ierr)
168:       CHKERRQ(ierr)

170: ! Set Variable bounds
171:       call MSA_BoundaryConditions(ierr)
172:       CHKERRQ(ierr)
173:       call TaoSetVariableBoundsRoutine(tao,MSA_Plate,                   &
174:      &     PETSC_NULL_OBJECT,ierr)
175:       CHKERRQ(ierr)

177: ! Set the initial solution guess
178:       call MSA_InitialPoint(x, ierr)
179:       CHKERRQ(ierr)
180:       call TaoSetInitialVector(tao,x,ierr)
181:       CHKERRQ(ierr)

183: ! Check for any tao command line options
184:       call TaoSetFromOptions(tao,ierr)
185:       CHKERRQ(ierr)

187: ! Solve the application
188:       call TaoSolve(tao,ierr)
189:       CHKERRQ(ierr)

191: ! Free TAO data structures
192:       call TaoDestroy(tao,ierr)
193:       CHKERRQ(ierr)

195: ! Free PETSc data structures
196:       call VecDestroy(x,ierr)
197:       call VecDestroy(Top,ierr)
198:       call VecDestroy(Bottom,ierr)
199:       call VecDestroy(Left,ierr)
200:       call VecDestroy(Right,ierr)
201:       call MatDestroy(H,ierr)
202:       call VecDestroy(localX,ierr)
203:       call VecDestroy(localV,ierr)
204:       call DMDestroy(dm,ierr)

206: ! Finalize TAO

208:       call PetscFinalize(ierr)

210:       end

212: ! ---------------------------------------------------------------------
213: !
214: !  FormFunctionGradient - Evaluates function f(X).
215: !
216: !  Input Parameters:
217: !  tao   - the Tao context
218: !  X     - the input vector
219: !  dummy - optional user-defined context, as set by TaoSetFunction()
220: !          (not used here)
221: !
222: !  Output Parameters:
223: !  fcn     - the newly evaluated function
224: !  G       - the gradient vector
225: !  info  - error code
226: !


229:       subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
230:       implicit none

232: ! dm, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
233: #include "plate2f.h"

235: ! Input/output variables

237:       Tao        tao
238:       PetscReal      fcn
239:       Vec              X, G
240:       PetscErrorCode   ierr
241:       PetscInt         dummy

243:       PetscInt         i,j,row
244:       PetscInt         xs, xm
245:       PetscInt         gxs, gxm
246:       PetscInt         ys, ym
247:       PetscInt         gys, gym
248:       PetscReal      ft,zero,hx,hy,hydhx,hxdhy
249:       PetscReal      area,rhx,rhy
250:       PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3
251:       PetscReal      d4,d5,d6,d7,d8
252:       PetscReal      df1dxc,df2dxc,df3dxc,df4dxc
253:       PetscReal      df5dxc,df6dxc
254:       PetscReal      xc,xl,xr,xt,xb,xlt,xrb


257: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
258: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
259: ! will return an array of doubles referenced by x_array offset by x_index.
260: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
261: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
262:       PetscReal      g_v(0:1),x_v(0:1)
263:       PetscReal      top_v(0:1),left_v(0:1)
264:       PetscReal      right_v(0:1),bottom_v(0:1)
265:       PetscOffset      g_i,left_i,right_i
266:       PetscOffset      bottom_i,top_i,x_i

268:       ft = 0.0
269:       zero = 0.0
270:       hx = 1.0/(mx + 1)
271:       hy = 1.0/(my + 1)
272:       hydhx = hy/hx
273:       hxdhy = hx/hy
274:       area = 0.5 * hx * hy
275:       rhx = mx + 1.0
276:       rhy = my + 1.0


279: ! Get local mesh boundaries
280:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
281:      &                  PETSC_NULL_INTEGER,ierr)
282:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,             &
283:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

285: ! Scatter ghost points to local vector
286:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
287:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

289: ! Initialize the vector to zero
290:       call VecSet(localV,zero,ierr)

292: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
293:       call VecGetArray(localX,x_v,x_i,ierr)
294:       call VecGetArray(localV,g_v,g_i,ierr)
295:       call VecGetArray(Top,top_v,top_i,ierr)
296:       call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
297:       call VecGetArray(Left,left_v,left_i,ierr)
298:       call VecGetArray(Right,right_v,right_i,ierr)

300: ! Compute function over the locally owned part of the mesh
301:       do j = ys,ys+ym-1
302:          do i = xs,xs+xm-1
303:             row = (j-gys)*gxm + (i-gxs)
304:             xc = x_v(row+x_i)
305:             xt = xc
306:             xb = xc
307:             xr = xc
308:             xl = xc
309:             xrb = xc
310:             xlt = xc

312:             if (i .eq. 0) then !left side
313:                xl = left_v(j - ys + 1 + left_i)
314:                xlt = left_v(j - ys + 2 + left_i)
315:             else
316:                xl = x_v(row - 1 + x_i)
317:             endif

319:             if (j .eq. 0) then !bottom side
320:                xb = bottom_v(i - xs + 1 + bottom_i)
321:                xrb = bottom_v(i - xs + 2 + bottom_i)
322:             else
323:                xb = x_v(row - gxm + x_i)
324:             endif

326:             if (i + 1 .eq. gxs + gxm) then !right side
327:                xr = right_v(j - ys + 1 + right_i)
328:                xrb = right_v(j - ys + right_i)
329:             else
330:                xr = x_v(row + 1 + x_i)
331:             endif

333:             if (j + 1 .eq. gys + gym) then !top side
334:                xt = top_v(i - xs + 1 + top_i)
335:                xlt = top_v(i - xs + top_i)
336:             else
337:                xt = x_v(row + gxm + x_i)
338:             endif

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

344:             if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
345:                xrb = x_v(row + 1 - gxm + x_i)
346:             endif

348:             d1 = xc-xl
349:             d2 = xc-xr
350:             d3 = xc-xt
351:             d4 = xc-xb
352:             d5 = xr-xrb
353:             d6 = xrb-xb
354:             d7 = xlt-xl
355:             d8 = xt-xlt

357:             df1dxc = d1 * hydhx
358:             df2dxc = d1 * hydhx + d4 * hxdhy
359:             df3dxc = d3 * hxdhy
360:             df4dxc = d2 * hydhx + d3 * hxdhy
361:             df5dxc = d2 * hydhx
362:             df6dxc = d4 * hxdhy

364:             d1 = d1 * rhx
365:             d2 = d2 * rhx
366:             d3 = d3 * rhy
367:             d4 = d4 * rhy
368:             d5 = d5 * rhy
369:             d6 = d6 * rhx
370:             d7 = d7 * rhy
371:             d8 = d8 * rhx

373:             f1 = sqrt(1.0 + d1*d1 + d7*d7)
374:             f2 = sqrt(1.0 + d1*d1 + d4*d4)
375:             f3 = sqrt(1.0 + d3*d3 + d8*d8)
376:             f4 = sqrt(1.0 + d3*d3 + d2*d2)
377:             f5 = sqrt(1.0 + d2*d2 + d5*d5)
378:             f6 = sqrt(1.0 + d4*d4 + d6*d6)

380:             ft = ft + f2 + f4

382:             df1dxc = df1dxc / f1
383:             df2dxc = df2dxc / f2
384:             df3dxc = df3dxc / f3
385:             df4dxc = df4dxc / f4
386:             df5dxc = df5dxc / f5
387:             df6dxc = df6dxc / f6

389:             g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc +  &
390:      &                              df5dxc + df6dxc)
391:          enddo
392:       enddo

394: ! Compute triangular areas along the border of the domain.
395:       if (xs .eq. 0) then  ! left side
396:          do j=ys,ys+ym-1
397:             d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i))         &
398:      &                 * rhy
399:             d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i))        &
400:      &                 * rhx
401:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
402:          enddo
403:       endif


406:       if (ys .eq. 0) then !bottom side
407:          do i=xs,xs+xm-1
408:             d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i))    &
409:      &                    * rhx
410:             d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
411:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
412:          enddo
413:       endif


416:       if (xs + xm .eq. mx) then ! right side
417:          do j=ys,ys+ym-1
418:             d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
419:             d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
420:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
421:          enddo
422:       endif


425:       if (ys + ym .eq. my) then
426:          do i=xs,xs+xm-1
427:             d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
428:             d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
429:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
430:          enddo
431:       endif


434:       if ((ys .eq. 0) .and. (xs .eq. 0)) then
435:          d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
436:          d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
437:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
438:       endif

440:       if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
441:          d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
442:          d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
443:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
444:       endif

446:       ft = ft * area
447:       call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR,                            &
448:      &             MPIU_SUM,MPI_COMM_WORLD,ierr)


451: ! Restore vectors
452:       call VecRestoreArray(localX,x_v,x_i,ierr)
453:       call VecRestoreArray(localV,g_v,g_i,ierr)
454:       call VecRestoreArray(Left,left_v,left_i,ierr)
455:       call VecRestoreArray(Top,top_v,top_i,ierr)
456:       call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
457:       call VecRestoreArray(Right,right_v,right_i,ierr)

459: ! Scatter values to global vector
460:       call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
461:       call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)

463:       call PetscLogFlops(70.0d0*xm*ym,ierr)

465:       return
466:       end  !FormFunctionGradient





472: ! ----------------------------------------------------------------------------
473: !
474: !
475: !   FormHessian - Evaluates Hessian matrix.
476: !
477: !   Input Parameters:
478: !.  tao  - the Tao context
479: !.  X    - input vector
480: !.  dummy  - not used
481: !
482: !   Output Parameters:
483: !.  Hessian    - Hessian matrix
484: !.  Hpc    - optionally different preconditioning matrix
485: !.  flag - flag indicating matrix structure
486: !
487: !   Notes:
488: !   Due to mesh point reordering with DMs, we must always work
489: !   with the local mesh points, and then transform them to the new
490: !   global numbering with the local-to-global mapping.  We cannot work
491: !   directly with the global numbers for the original uniprocessor mesh!
492: !
493: !      MatSetValuesLocal(), using the local ordering (including
494: !         ghost points!)
495: !         - Then set matrix entries using the local ordering
496: !           by calling MatSetValuesLocal()

498:       subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
499:       implicit none

501: ! dm,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
502: #include "plate2f.h"

504:       Tao     tao
505:       Vec            X
506:       Mat            Hessian,Hpc
507:       PetscInt       dummy
508:       PetscErrorCode ierr

510:       PetscInt       i,j,k,row
511:       PetscInt       xs,xm,gxs,gxm
512:       PetscInt       ys,ym,gys,gym
513:       PetscInt       col(0:6)
514:       PetscReal    hx,hy,hydhx,hxdhy,rhx,rhy
515:       PetscReal    f1,f2,f3,f4,f5,f6,d1,d2,d3
516:       PetscReal    d4,d5,d6,d7,d8
517:       PetscReal    xc,xl,xr,xt,xb,xlt,xrb
518:       PetscReal    hl,hr,ht,hb,hc,htl,hbr

520: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
521: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
522: ! will return an array of doubles referenced by x_array offset by x_index.
523: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
524: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
525:       PetscReal   right_v(0:1),left_v(0:1)
526:       PetscReal   bottom_v(0:1),top_v(0:1)
527:       PetscReal   x_v(0:1)
528:       PetscOffset   x_i,right_i,left_i
529:       PetscOffset   bottom_i,top_i
530:       PetscReal   v(0:6)
531:       PetscBool     assembled
532:       PetscInt      i1

534:       i1=1

536: ! Set various matrix options
537:       call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,              &
538:      &                  PETSC_TRUE,ierr)

540: ! Get local mesh boundaries
541:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
542:      &                  PETSC_NULL_INTEGER,ierr)
543:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
544:      &                       PETSC_NULL_INTEGER,ierr)

546: ! Scatter ghost points to local vectors
547:       call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
548:       call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)

550: ! Get pointers to vector data (see note on Fortran arrays above)
551:       call VecGetArray(localX,x_v,x_i,ierr)
552:       call VecGetArray(Top,top_v,top_i,ierr)
553:       call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
554:       call VecGetArray(Left,left_v,left_i,ierr)
555:       call VecGetArray(Right,right_v,right_i,ierr)

557: ! Initialize matrix entries to zero
558:       call MatAssembled(Hessian,assembled,ierr)
559:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)


562:       rhx = mx + 1.0
563:       rhy = my + 1.0
564:       hx = 1.0/rhx
565:       hy = 1.0/rhy
566:       hydhx = hy/hx
567:       hxdhy = hx/hy
568: ! compute Hessian over the locally owned part of the mesh

570:       do  i=xs,xs+xm-1
571:          do  j=ys,ys+ym-1
572:             row = (j-gys)*gxm + (i-gxs)

574:             xc = x_v(row + x_i)
575:             xt = xc
576:             xb = xc
577:             xr = xc
578:             xl = xc
579:             xrb = xc
580:             xlt = xc

582:             if (i .eq. gxs) then   ! Left side
583:                xl = left_v(left_i + j - ys + 1)
584:                xlt = left_v(left_i + j - ys + 2)
585:             else
586:                xl = x_v(x_i + row -1 )
587:             endif

589:             if (j .eq. gys) then ! bottom side
590:                xb = bottom_v(bottom_i + i - xs + 1)
591:                xrb = bottom_v(bottom_i + i - xs + 2)
592:             else
593:                xb = x_v(x_i + row - gxm)
594:             endif

596:             if (i+1 .eq. gxs + gxm) then !right side
597:                xr = right_v(right_i + j - ys + 1)
598:                xrb = right_v(right_i + j - ys)
599:             else
600:                xr = x_v(x_i + row + 1)
601:             endif

603:             if (j+1 .eq. gym+gys) then !top side
604:                xt = top_v(top_i +i - xs + 1)
605:                xlt = top_v(top_i + i - xs)
606:             else
607:                xt = x_v(x_i + row + gxm)
608:             endif

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

614:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
615:                xrb = x_v(x_i + row + 1 - gxm)
616:             endif

618:             d1 = (xc-xl)*rhx
619:             d2 = (xc-xr)*rhx
620:             d3 = (xc-xt)*rhy
621:             d4 = (xc-xb)*rhy
622:             d5 = (xrb-xr)*rhy
623:             d6 = (xrb-xb)*rhx
624:             d7 = (xlt-xl)*rhy
625:             d8 = (xlt-xt)*rhx

627:             f1 = sqrt( 1.0 + d1*d1 + d7*d7)
628:             f2 = sqrt( 1.0 + d1*d1 + d4*d4)
629:             f3 = sqrt( 1.0 + d3*d3 + d8*d8)
630:             f4 = sqrt( 1.0 + d3*d3 + d2*d2)
631:             f5 = sqrt( 1.0 + d2*d2 + d5*d5)
632:             f6 = sqrt( 1.0 + d4*d4 + d6*d6)


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

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

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

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

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

650:             hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) +                         &
651:      &              hxdhy*(1.0+d8*d8)/(f3*f3*f3) +                      &
652:      &              hydhx*(1.0+d5*d5)/(f5*f5*f5) +                      &
653:      &              hxdhy*(1.0+d6*d6)/(f6*f6*f6) +                      &
654:      &              (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-               &
655:      &              2*d1*d4)/(f2*f2*f2) +  (hxdhy*(1.0+d2*d2)+          &
656:      &              hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)

658:             hl = hl * 0.5
659:             hr = hr * 0.5
660:             ht = ht * 0.5
661:             hb = hb * 0.5
662:             hbr = hbr * 0.5
663:             htl = htl * 0.5
664:             hc = hc * 0.5

666:             k = 0

668:             if (j .gt. 0) then
669:                v(k) = hb
670:                col(k) = row - gxm
671:                k=k+1
672:             endif

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

680:             if (i .gt. 0) then
681:                v(k) = hl
682:                col(k) = row - 1
683:                k = k+1
684:             endif

686:             v(k) = hc
687:             col(k) = row
688:             k=k+1

690:             if (i .lt. mx-1) then
691:                v(k) = hr
692:                col(k) = row + 1
693:                k=k+1
694:             endif

696:             if ((i .gt. 0) .and. (j .lt. my-1)) then
697:                v(k) = htl
698:                col(k) = row + gxm - 1
699:                k=k+1
700:             endif

702:             if (j .lt. my-1) then
703:                v(k) = ht
704:                col(k) = row + gxm
705:                k=k+1
706:             endif

708: ! Set matrix values using local numbering, defined earlier in main routine
709:             call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES,      &
710:      &                              ierr)



714:          enddo
715:       enddo

717: ! restore vectors
718:       call VecRestoreArray(localX,x_v,x_i,ierr)
719:       call VecRestoreArray(Left,left_v,left_i,ierr)
720:       call VecRestoreArray(Right,right_v,right_i,ierr)
721:       call VecRestoreArray(Top,top_v,top_i,ierr)
722:       call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)


725: ! Assemble the matrix
726:       call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
727:       call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)

729:       call PetscLogFlops(199.0d0*xm*ym,ierr)

731:       return
732:       end





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

740: ! ----------------------------------------------------------------------------
741: !
742: !/*
743: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
744: !
745: !
746: !*/

748:       subroutine MSA_BoundaryConditions(ierr)
749:       implicit none

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

754:       PetscErrorCode   ierr
755:       PetscInt         i,j,k,limit,maxits
756:       PetscInt         xs, xm, gxs, gxm
757:       PetscInt         ys, ym, gys, gym
758:       PetscInt         bsize, lsize
759:       PetscInt         tsize, rsize
760:       PetscReal      one,two,three,tol
761:       PetscReal      scl,fnorm,det,xt
762:       PetscReal      yt,hx,hy,u1,u2,nf1,nf2
763:       PetscReal      njac11,njac12,njac21,njac22
764:       PetscReal      b, t, l, r
765:       PetscReal      boundary_v(0:1)
766:       PetscOffset      boundary_i
767:       logical exitloop
768:       PetscBool flg

770:       limit=0
771:       maxits = 5
772:       tol=1e-10
773:       b=-0.5
774:       t= 0.5
775:       l=-0.5
776:       r= 0.5
777:       xt=0
778:       yt=0
779:       one=1.0
780:       two=2.0
781:       three=3.0


784:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
785:      &                  PETSC_NULL_INTEGER,ierr)
786:       call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,              &
787:      &                       gxm,gym,PETSC_NULL_INTEGER,ierr)

789:       bsize = xm + 2
790:       lsize = ym + 2
791:       rsize = ym + 2
792:       tsize = xm + 2


795:       call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
796:       call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
797:       call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
798:       call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)

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

803:       do j=0,3

805:          if (j.eq.0) then
806:             yt=b
807:             xt=l+hx*xs
808:             limit=bsize
809:             call VecGetArray(Bottom,boundary_v,boundary_i,ierr)


812:          elseif (j.eq.1) then
813:             yt=t
814:             xt=l+hx*xs
815:             limit=tsize
816:             call VecGetArray(Top,boundary_v,boundary_i,ierr)

818:          elseif (j.eq.2) then
819:             yt=b+hy*ys
820:             xt=l
821:             limit=lsize
822:             call VecGetArray(Left,boundary_v,boundary_i,ierr)

824:          elseif (j.eq.3) then
825:             yt=b+hy*ys
826:             xt=r
827:             limit=rsize
828:             call VecGetArray(Right,boundary_v,boundary_i,ierr)
829:          endif


832:          do i=0,limit-1

834:             u1=xt
835:             u2=-yt
836:             k = 0
837:             exitloop = .false.
838:             do while (k .lt. maxits .and. (.not. exitloop) )

840:                nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
841:                nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
842:                fnorm=sqrt(nf1*nf1+nf2*nf2)
843:                if (fnorm .gt. tol) then
844:                   njac11=one+u2*u2-u1*u1
845:                   njac12=two*u1*u2
846:                   njac21=-two*u1*u2
847:                   njac22=-one - u1*u1 + u2*u2
848:                   det = njac11*njac22-njac21*njac12
849:                   u1 = u1-(njac22*nf1-njac12*nf2)/det
850:                   u2 = u2-(njac11*nf2-njac21*nf1)/det
851:                else
852:                   exitloop = .true.
853:                endif
854:                k=k+1
855:             enddo

857:             boundary_v(i + boundary_i) = u1*u1-u2*u2
858:             if ((j .eq. 0) .or. (j .eq. 1)) then
859:                xt = xt + hx
860:             else
861:                yt = yt + hy
862:             endif

864:          enddo


867:          if (j.eq.0) then
868:             call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
869:          elseif (j.eq.1) then
870:             call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
871:          elseif (j.eq.2) then
872:             call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
873:          elseif (j.eq.3) then
874:             call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
875:          endif

877:       enddo


880: ! Scale the boundary if desired
881:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,          &
882:      &                         '-bottom',scl,flg,ierr)
883:       if (flg .eqv. PETSC_TRUE) then
884:          call VecScale(scl,Bottom,ierr)
885:       endif

887:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,          &
888:      &                         '-top',scl,flg,ierr)
889:       if (flg .eqv. PETSC_TRUE) then
890:          call VecScale(scl,Top,ierr)
891:       endif

893:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,          &
894:      &                         '-right',scl,flg,ierr)
895:       if (flg .eqv. PETSC_TRUE) then
896:          call VecScale(scl,Right,ierr)
897:       endif

899:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,          &
900:      &                         '-left',scl,flg,ierr)
901:       if (flg .eqv. PETSC_TRUE) then
902:          call VecScale(scl,Left,ierr)
903:       endif


906:       return
907:       end

909: ! ----------------------------------------------------------------------------
910: !
911: !/*
912: !     MSA_Plate - Calculates an obstacle for surface to stretch over
913: !
914: !     Output Parameter:
915: !.    xl - lower bound vector
916: !.    xu - upper bound vector
917: !
918: !*/

920:       subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
921:       implicit none
922: ! mx,my,bmx,bmy,dm,bheight defined in plate2f.h
923: #include "plate2f.h"
924:       Tao        tao
925:       Vec              xl,xu
926:       PetscErrorCode   ierr
927:       PetscInt         i,j,row
928:       PetscInt         xs, xm, ys, ym
929:       PetscReal      lb,ub
930:       PetscInt         dummy

932: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
933: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
934: ! will return an array of doubles referenced by x_array offset by x_index.
935: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
936: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
937:       PetscReal      xl_v(0:1)
938:       PetscOffset      xl_i

940:       print *,'msa_plate'

942:       lb = PETSC_NINFINITY
943:       ub = PETSC_INFINITY

945:       if (bmy .lt. 0) bmy = 0
946:       if (bmy .gt. my) bmy = my
947:       if (bmx .lt. 0) bmx = 0
948:       if (bmx .gt. mx) bmx = mx


951:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,               &
952:      &             PETSC_NULL_INTEGER,ierr)

954:       call VecSet(xl,lb,ierr)
955:       call VecSet(xu,ub,ierr)

957:       call VecGetArray(xl,xl_v,xl_i,ierr)


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

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

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

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

970:             endif

972:          enddo
973:       enddo


976:       call VecRestoreArray(xl,xl_v,xl_i,ierr)

978:       return
979:       end





985: ! ----------------------------------------------------------------------------
986: !
987: !/*
988: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
989: !
990: !     Output Parameter:
991: !.    X - vector for initial guess
992: !
993: !*/

995:       subroutine MSA_InitialPoint(X, ierr)
996:       implicit none

998: ! mx,my,localX,dm,Top,Left,Bottom,Right defined in plate2f.h
999: #include "plate2f.h"
1000:       Vec               X
1001:       PetscErrorCode    ierr
1002:       PetscInt          start,i,j
1003:       PetscInt          row
1004:       PetscInt          xs,xm,gxs,gxm
1005:       PetscInt          ys,ym,gys,gym
1006:       PetscReal       zero, np5

1008: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1009: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
1010: ! will return an array of doubles referenced by x_array offset by x_index.
1011: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
1012: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1013:       PetscReal   left_v(0:1),right_v(0:1)
1014:       PetscReal   bottom_v(0:1),top_v(0:1)
1015:       PetscReal   x_v(0:1)
1016:       PetscOffset   left_i, right_i, top_i
1017:       PetscOffset   bottom_i,x_i
1018:       PetscBool     flg
1019:       PetscRandom   rctx

1021:       zero = 0.0
1022:       np5 = -0.5

1024:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,        &
1025:      &                        '-start', start,flg,ierr)

1027:       if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then  ! the zero vector is reasonable
1028:          call VecSet(X,zero,ierr)

1030:       elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then  ! random start -0.5 < xi < 0.5
1031:          call PetscRandomCreate(MPI_COMM_WORLD,rctx,ierr)
1032:          do i=0,start-1
1033:             call VecSetRandom(X,rctx,ierr)
1034:          enddo

1036:          call PetscRandomDestroy(rctx,ierr)
1037:          call VecShift(X,np5,ierr)

1039:       else   ! average of boundary conditions

1041: !        Get Local mesh boundaries
1042:          call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,             &
1043:      &                     PETSC_NULL_INTEGER,ierr)
1044:          call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,    &
1045:      &                     PETSC_NULL_INTEGER,ierr)



1049: !        Get pointers to vector data
1050:          call VecGetArray(Top,top_v,top_i,ierr)
1051:          call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1052:          call VecGetArray(Left,left_v,left_i,ierr)
1053:          call VecGetArray(Right,right_v,right_i,ierr)

1055:          call VecGetArray(localX,x_v,x_i,ierr)

1057: !        Perform local computations
1058:          do  j=ys,ys+ym-1
1059:             do i=xs,xs+xm-1
1060:                row = (j-gys)*gxm  + (i-gxs)
1061:                x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my        &
1062:      &             + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) +                  &
1063:      &              (i+1)*left_v(left_i+j-ys+1)/mx       +                  &
1064:      &              (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1065:             enddo
1066:          enddo

1068: !        Restore vectors
1069:          call VecRestoreArray(localX,x_v,x_i,ierr)

1071:          call VecRestoreArray(Left,left_v,left_i,ierr)
1072:          call VecRestoreArray(Top,top_v,top_i,ierr)
1073:          call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1074:          call VecRestoreArray(Right,right_v,right_i,ierr)

1076:          call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1077:          call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)

1079:       endif

1081:       return
1082:       end