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              tao           ! Tao solver context
47:       Mat              H             ! Hessian matrix
48:       ISLocalToGlobalMapping isltog  ! local to global mapping object
49:       PetscBool        flg
50:       PetscInt         i1,i3,i7

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,PETSC_NULL_INTEGER,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,i3,PETSC_NULL_INTEGER,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,tao,ierr))
124:       PetscCallA(TaoSetType(tao,TAOBLMVM,ierr))

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

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

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

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

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

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

146: ! Free TAO data structures
147:       PetscCallA(TaoDestroy(tao,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: !

183:       use plate2fmodule
184:       implicit none

186: ! Input/output variables

188:       Tao        tao
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 VecGetArrayF90 in Fortran)
233:       PetscCall(VecGetArrayF90(localX,x_v,ierr))
234:       PetscCall(VecGetArrayF90(localV,g_v,ierr))
235:       PetscCall(VecGetArrayF90(Top,top_v,ierr))
236:       PetscCall(VecGetArrayF90(Bottom,bottom_v,ierr))
237:       PetscCall(VecGetArrayF90(Left,left_v,ierr))
238:       PetscCall(VecGetArrayF90(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 .eq. 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:             endif

259:             if (j .eq. 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:             endif

266:             if (i + 1 .eq. 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:             endif

273:             if (j + 1 .eq. 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:             endif

280:             if ((i .gt. gxs) .and. (j + 1 .lt. gys + gym)) then
281:                xlt = x_v(1+row - 1 + gxm)
282:             endif

284:             if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
285:                xrb = x_v(1+row + 1 - gxm)
286:             endif

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:          enddo
331:       enddo

333: ! Compute triangular areas along the border of the domain.
334:       if (xs .eq. 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:          enddo
340:       endif

342:       if (ys .eq. 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:          enddo
348:       endif

350:       if (xs + xm .eq. 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:          enddo
356:       endif

358:       if (ys + ym .eq. 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:          enddo
364:       endif

366:       if ((ys .eq. 0) .and. (xs .eq. 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:       endif

372:       if ((ys + ym .eq. my) .and. (xs + xm .eq. 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:       endif

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

381: ! Restore vectors
382:       PetscCall(VecRestoreArrayF90(localX,x_v,ierr))
383:       PetscCall(VecRestoreArrayF90(localV,g_v,ierr))
384:       PetscCall(VecRestoreArrayF90(Left,left_v,ierr))
385:       PetscCall(VecRestoreArrayF90(Top,top_v,ierr))
386:       PetscCall(VecRestoreArrayF90(Bottom,bottom_v,ierr))
387:       PetscCall(VecRestoreArrayF90(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))

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 preconditioning matrix
410: !.  flag - flag indicating matrix structure
411: !
412: !   Notes:
413: !   Due to mesh point reordering with DMs, we must always work
414: !   with the local mesh points, and then transform them to the new
415: !   global numbering with the local-to-global mapping.  We cannot work
416: !   directly with the global numbers for the original uniprocessor mesh!
417: !
418: !      MatSetValuesLocal(), using the local ordering (including
419: !         ghost points!)
420: !         - Then set matrix entries using the local ordering
421: !           by calling MatSetValuesLocal()

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

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

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

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

450:       i1=1

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

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

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

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

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

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

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

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

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

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

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

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

522:             if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
523:                xlt = x_v(1+row - 1 + gxm)
524:             endif

526:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
527:                xrb = x_v(1+row + 1 - gxm)
528:             endif

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

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

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

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

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

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

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

557:             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) +                      &
558:      &           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)+   &
559:      &           hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)

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

569:             k = 0

571:             if (j .gt. 0) then
572:                v(k) = hb
573:                col(k) = row - gxm
574:                k=k+1
575:             endif

577:             if ((j .gt. 0) .and. (i .lt. mx-1)) then
578:                v(k) = hbr
579:                col(k) = row-gxm+1
580:                k=k+1
581:             endif

583:             if (i .gt. 0) then
584:                v(k) = hl
585:                col(k) = row - 1
586:                k = k+1
587:             endif

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

593:             if (i .lt. mx-1) then
594:                v(k) = hr
595:                col(k) = row + 1
596:                k=k+1
597:             endif

599:             if ((i .gt. 0) .and. (j .lt. my-1)) then
600:                v(k) = htl
601:                col(k) = row + gxm - 1
602:                k=k+1
603:             endif

605:             if (j .lt. my-1) then
606:                v(k) = ht
607:                col(k) = row + gxm
608:                k=k+1
609:             endif

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

614:          enddo
615:       enddo

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

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

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

630:       end

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

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

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

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

659:       logical exitloop
660:       PetscBool flg

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

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

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

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

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

691:       do j=0,3

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

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

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

711:          elseif (j.eq.3) then
712:             yt=b+hy*ys
713:             xt=r
714:             limit=rsize
715:             PetscCall(VecGetArrayF90(Right,boundary_v,ierr))
716:          endif

718:          do i=0,limit-1

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

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

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

750:          enddo

752:          if (j.eq.0) then
753:             PetscCall(VecRestoreArrayF90(Bottom,boundary_v,ierr))
754:          elseif (j.eq.1) then
755:             PetscCall(VecRestoreArrayF90(Top,boundary_v,ierr))
756:          elseif (j.eq.2) then
757:             PetscCall(VecRestoreArrayF90(Left,boundary_v,ierr))
758:          elseif (j.eq.3) then
759:             PetscCall(VecRestoreArrayF90(Right,boundary_v,ierr))
760:          endif

762:       enddo

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

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

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

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

785:       end

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

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

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

811:       lb = PETSC_NINFINITY
812:       ub = PETSC_INFINITY

814:       if (bmy .lt. 0) bmy = 0
815:       if (bmy .gt. my) bmy = my
816:       if (bmx .lt. 0) bmx = 0
817:       if (bmx .gt. mx) bmx = mx

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

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

824:       PetscCall(VecGetArrayF90(xl,xl_v,ierr))

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

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

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

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

836:             endif

838:          enddo
839:       enddo

841:       PetscCall(VecRestoreArrayF90(xl,xl_v,ierr))

843:       end

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

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

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

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

873:       zero = 0.0
874:       np5 = -0.5

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

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

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

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

890:       else   ! average of boundary conditions

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

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

902:          PetscCall(VecGetArrayF90(localX,x_v,ierr))

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

913: !        Restore vectors
914:          PetscCall(VecRestoreArrayF90(localX,x_v,ierr))

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

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

924:       endif

926:       end

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