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

 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,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

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

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: !

182:       subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
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: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
208: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
209: ! will return an array of doubles referenced by x_array offset by x_index.
210: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
211: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
212:       PetscReal      g_v(0:1),x_v(0:1)
213:       PetscReal      top_v(0:1),left_v(0:1)
214:       PetscReal      right_v(0:1),bottom_v(0:1)
215:       PetscOffset      g_i,left_i,right_i
216:       PetscOffset      bottom_i,top_i,x_i

218:       ft = 0.0
219:       zero = 0.0
220:       hx = 1.0/real(mx + 1)
221:       hy = 1.0/real(my + 1)
222:       hydhx = hy/hx
223:       hxdhy = hx/hy
224:       area = 0.5 * hx * hy
225:       rhx = real(mx) + 1.0
226:       rhy = real(my) + 1.0

228: ! Get local mesh boundaries
229:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
230:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

232: ! Scatter ghost points to local vector
233:       PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
234:       PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))

236: ! Initialize the vector to zero
237:       PetscCall(VecSet(localV,zero,ierr))

239: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
240:       PetscCall(VecGetArray(localX,x_v,x_i,ierr))
241:       PetscCall(VecGetArray(localV,g_v,g_i,ierr))
242:       PetscCall(VecGetArray(Top,top_v,top_i,ierr))
243:       PetscCall(VecGetArray(Bottom,bottom_v,bottom_i,ierr))
244:       PetscCall(VecGetArray(Left,left_v,left_i,ierr))
245:       PetscCall(VecGetArray(Right,right_v,right_i,ierr))

247: ! Compute function over the locally owned part of the mesh
248:       do j = ys,ys+ym-1
249:          do i = xs,xs+xm-1
250:             row = (j-gys)*gxm + (i-gxs)
251:             xc = x_v(row+x_i)
252:             xt = xc
253:             xb = xc
254:             xr = xc
255:             xl = xc
256:             xrb = xc
257:             xlt = xc

259:             if (i .eq. 0) then !left side
260:                xl = left_v(j - ys + 1 + left_i)
261:                xlt = left_v(j - ys + 2 + left_i)
262:             else
263:                xl = x_v(row - 1 + x_i)
264:             endif

266:             if (j .eq. 0) then !bottom side
267:                xb = bottom_v(i - xs + 1 + bottom_i)
268:                xrb = bottom_v(i - xs + 2 + bottom_i)
269:             else
270:                xb = x_v(row - gxm + x_i)
271:             endif

273:             if (i + 1 .eq. gxs + gxm) then !right side
274:                xr = right_v(j - ys + 1 + right_i)
275:                xrb = right_v(j - ys + right_i)
276:             else
277:                xr = x_v(row + 1 + x_i)
278:             endif

280:             if (j + 1 .eq. gys + gym) then !top side
281:                xt = top_v(i - xs + 1 + top_i)
282:                xlt = top_v(i - xs + top_i)
283:             else
284:                xt = x_v(row + gxm + x_i)
285:             endif

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

291:             if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
292:                xrb = x_v(row + 1 - gxm + x_i)
293:             endif

295:             d1 = xc-xl
296:             d2 = xc-xr
297:             d3 = xc-xt
298:             d4 = xc-xb
299:             d5 = xr-xrb
300:             d6 = xrb-xb
301:             d7 = xlt-xl
302:             d8 = xt-xlt

304:             df1dxc = d1 * hydhx
305:             df2dxc = d1 * hydhx + d4 * hxdhy
306:             df3dxc = d3 * hxdhy
307:             df4dxc = d2 * hydhx + d3 * hxdhy
308:             df5dxc = d2 * hydhx
309:             df6dxc = d4 * hxdhy

311:             d1 = d1 * rhx
312:             d2 = d2 * rhx
313:             d3 = d3 * rhy
314:             d4 = d4 * rhy
315:             d5 = d5 * rhy
316:             d6 = d6 * rhx
317:             d7 = d7 * rhy
318:             d8 = d8 * rhx

320:             f1 = sqrt(1.0 + d1*d1 + d7*d7)
321:             f2 = sqrt(1.0 + d1*d1 + d4*d4)
322:             f3 = sqrt(1.0 + d3*d3 + d8*d8)
323:             f4 = sqrt(1.0 + d3*d3 + d2*d2)
324:             f5 = sqrt(1.0 + d2*d2 + d5*d5)
325:             f6 = sqrt(1.0 + d4*d4 + d6*d6)

327:             ft = ft + f2 + f4

329:             df1dxc = df1dxc / f1
330:             df2dxc = df2dxc / f2
331:             df3dxc = df3dxc / f3
332:             df4dxc = df4dxc / f4
333:             df5dxc = df5dxc / f5
334:             df6dxc = df6dxc / f6

336:             g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc)
337:          enddo
338:       enddo

340: ! Compute triangular areas along the border of the domain.
341:       if (xs .eq. 0) then  ! left side
342:          do j=ys,ys+ym-1
343:             d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) * rhy
344:             d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) * rhx
345:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
346:          enddo
347:       endif

349:       if (ys .eq. 0) then !bottom side
350:          do i=xs,xs+xm-1
351:             d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) * rhx
352:             d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
353:             ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
354:          enddo
355:       endif

357:       if (xs + xm .eq. mx) then ! right side
358:          do j=ys,ys+ym-1
359:             d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
360:             d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
361:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
362:          enddo
363:       endif

365:       if (ys + ym .eq. my) then
366:          do i=xs,xs+xm-1
367:             d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
368:             d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
369:             ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
370:          enddo
371:       endif

373:       if ((ys .eq. 0) .and. (xs .eq. 0)) then
374:          d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
375:          d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
376:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
377:       endif

379:       if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
380:          d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
381:          d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
382:          ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
383:       endif

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

388: ! Restore vectors
389:       PetscCall(VecRestoreArray(localX,x_v,x_i,ierr))
390:       PetscCall(VecRestoreArray(localV,g_v,g_i,ierr))
391:       PetscCall(VecRestoreArray(Left,left_v,left_i,ierr))
392:       PetscCall(VecRestoreArray(Top,top_v,top_i,ierr))
393:       PetscCall(VecRestoreArray(Bottom,bottom_v,bottom_i,ierr))
394:       PetscCall(VecRestoreArray(Right,right_v,right_i,ierr))

396: ! Scatter values to global vector
397:       PetscCall(DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr))
398:       PetscCall(DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr))

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

402:       return
403:       end  !FormFunctionGradient

405: ! ----------------------------------------------------------------------------
406: !
407: !
408: !   FormHessian - Evaluates Hessian matrix.
409: !
410: !   Input Parameters:
411: !.  tao  - the Tao context
412: !.  X    - input vector
413: !.  dummy  - not used
414: !
415: !   Output Parameters:
416: !.  Hessian    - Hessian matrix
417: !.  Hpc    - optionally different preconditioning matrix
418: !.  flag - flag indicating matrix structure
419: !
420: !   Notes:
421: !   Due to mesh point reordering with DMs, we must always work
422: !   with the local mesh points, and then transform them to the new
423: !   global numbering with the local-to-global mapping.  We cannot work
424: !   directly with the global numbers for the original uniprocessor mesh!
425: !
426: !      MatSetValuesLocal(), using the local ordering (including
427: !         ghost points!)
428: !         - Then set matrix entries using the local ordering
429: !           by calling MatSetValuesLocal()

431:       subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
432:       use plate2fmodule
433:       implicit none

435:       Tao     tao
436:       Vec            X
437:       Mat            Hessian,Hpc
438:       PetscInt       dummy
439:       PetscErrorCode ierr

441:       PetscInt       i,j,k,row
442:       PetscInt       xs,xm,gxs,gxm
443:       PetscInt       ys,ym,gys,gym
444:       PetscInt       col(0:6)
445:       PetscReal    hx,hy,hydhx,hxdhy,rhx,rhy
446:       PetscReal    f1,f2,f3,f4,f5,f6,d1,d2,d3
447:       PetscReal    d4,d5,d6,d7,d8
448:       PetscReal    xc,xl,xr,xt,xb,xlt,xrb
449:       PetscReal    hl,hr,ht,hb,hc,htl,hbr

451: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
452: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
453: ! will return an array of doubles referenced by x_array offset by x_index.
454: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
455: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
456:       PetscReal   right_v(0:1),left_v(0:1)
457:       PetscReal   bottom_v(0:1),top_v(0:1)
458:       PetscReal   x_v(0:1)
459:       PetscOffset   x_i,right_i,left_i
460:       PetscOffset   bottom_i,top_i
461:       PetscReal   v(0:6)
462:       PetscBool     assembled
463:       PetscInt      i1

465:       i1=1

467: ! Set various matrix options
468:       PetscCall(MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE,ierr))

470: ! Get local mesh boundaries
471:       PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
472:       PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

474: ! Scatter ghost points to local vectors
475:       PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
476:       PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))

478: ! Get pointers to vector data (see note on Fortran arrays above)
479:       PetscCall(VecGetArray(localX,x_v,x_i,ierr))
480:       PetscCall(VecGetArray(Top,top_v,top_i,ierr))
481:       PetscCall(VecGetArray(Bottom,bottom_v,bottom_i,ierr))
482:       PetscCall(VecGetArray(Left,left_v,left_i,ierr))
483:       PetscCall(VecGetArray(Right,right_v,right_i,ierr))

485: ! Initialize matrix entries to zero
486:       PetscCall(MatAssembled(Hessian,assembled,ierr))
487:       if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian,ierr))

489:       rhx = real(mx) + 1.0
490:       rhy = real(my) + 1.0
491:       hx = 1.0/rhx
492:       hy = 1.0/rhy
493:       hydhx = hy/hx
494:       hxdhy = hx/hy
495: ! compute Hessian over the locally owned part of the mesh

497:       do  i=xs,xs+xm-1
498:          do  j=ys,ys+ym-1
499:             row = (j-gys)*gxm + (i-gxs)

501:             xc = x_v(row + x_i)
502:             xt = xc
503:             xb = xc
504:             xr = xc
505:             xl = xc
506:             xrb = xc
507:             xlt = xc

509:             if (i .eq. gxs) then   ! Left side
510:                xl = left_v(left_i + j - ys + 1)
511:                xlt = left_v(left_i + j - ys + 2)
512:             else
513:                xl = x_v(x_i + row -1)
514:             endif

516:             if (j .eq. gys) then ! bottom side
517:                xb = bottom_v(bottom_i + i - xs + 1)
518:                xrb = bottom_v(bottom_i + i - xs + 2)
519:             else
520:                xb = x_v(x_i + row - gxm)
521:             endif

523:             if (i+1 .eq. gxs + gxm) then !right side
524:                xr = right_v(right_i + j - ys + 1)
525:                xrb = right_v(right_i + j - ys)
526:             else
527:                xr = x_v(x_i + row + 1)
528:             endif

530:             if (j+1 .eq. gym+gys) then !top side
531:                xt = top_v(top_i +i - xs + 1)
532:                xlt = top_v(top_i + i - xs)
533:             else
534:                xt = x_v(x_i + row + gxm)
535:             endif

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

541:             if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
542:                xrb = x_v(x_i + row + 1 - gxm)
543:             endif

545:             d1 = (xc-xl)*rhx
546:             d2 = (xc-xr)*rhx
547:             d3 = (xc-xt)*rhy
548:             d4 = (xc-xb)*rhy
549:             d5 = (xrb-xr)*rhy
550:             d6 = (xrb-xb)*rhx
551:             d7 = (xlt-xl)*rhy
552:             d8 = (xlt-xt)*rhx

554:             f1 = sqrt(1.0 + d1*d1 + d7*d7)
555:             f2 = sqrt(1.0 + d1*d1 + d4*d4)
556:             f3 = sqrt(1.0 + d3*d3 + d8*d8)
557:             f4 = sqrt(1.0 + d3*d3 + d2*d2)
558:             f5 = sqrt(1.0 + d2*d2 + d5*d5)
559:             f6 = sqrt(1.0 + d4*d4 + d6*d6)

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

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

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

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

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

572:             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) +                      &
573:      &           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)+   &
574:      &           hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)

576:             hl = hl * 0.5
577:             hr = hr * 0.5
578:             ht = ht * 0.5
579:             hb = hb * 0.5
580:             hbr = hbr * 0.5
581:             htl = htl * 0.5
582:             hc = hc * 0.5

584:             k = 0

586:             if (j .gt. 0) then
587:                v(k) = hb
588:                col(k) = row - gxm
589:                k=k+1
590:             endif

592:             if ((j .gt. 0) .and. (i .lt. mx-1)) then
593:                v(k) = hbr
594:                col(k) = row-gxm+1
595:                k=k+1
596:             endif

598:             if (i .gt. 0) then
599:                v(k) = hl
600:                col(k) = row - 1
601:                k = k+1
602:             endif

604:             v(k) = hc
605:             col(k) = row
606:             k=k+1

608:             if (i .lt. mx-1) then
609:                v(k) = hr
610:                col(k) = row + 1
611:                k=k+1
612:             endif

614:             if ((i .gt. 0) .and. (j .lt. my-1)) then
615:                v(k) = htl
616:                col(k) = row + gxm - 1
617:                k=k+1
618:             endif

620:             if (j .lt. my-1) then
621:                v(k) = ht
622:                col(k) = row + gxm
623:                k=k+1
624:             endif

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

629:          enddo
630:       enddo

632: ! restore vectors
633:       PetscCall(VecRestoreArray(localX,x_v,x_i,ierr))
634:       PetscCall(VecRestoreArray(Left,left_v,left_i,ierr))
635:       PetscCall(VecRestoreArray(Right,right_v,right_i,ierr))
636:       PetscCall(VecRestoreArray(Top,top_v,top_i,ierr))
637:       PetscCall(VecRestoreArray(Bottom,bottom_v,bottom_i,ierr))

639: ! Assemble the matrix
640:       PetscCall(MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr))
641:       PetscCall(MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr))

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

645:       return
646:       end

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

650: ! ----------------------------------------------------------------------------
651: !
652: !/*
653: !     MSA_BoundaryConditions - calculates the boundary conditions for the region
654: !
655: !
656: !*/

658:       subroutine MSA_BoundaryConditions(ierr)
659:       use plate2fmodule
660:       implicit none

662:       PetscErrorCode   ierr
663:       PetscInt         i,j,k,limit,maxits
664:       PetscInt         xs, xm, gxs, gxm
665:       PetscInt         ys, ym, gys, gym
666:       PetscInt         bsize, lsize
667:       PetscInt         tsize, rsize
668:       PetscReal      one,two,three,tol
669:       PetscReal      scl,fnorm,det,xt
670:       PetscReal      yt,hx,hy,u1,u2,nf1,nf2
671:       PetscReal      njac11,njac12,njac21,njac22
672:       PetscReal      b, t, l, r
673:       PetscReal      boundary_v(0:1)
674:       PetscOffset      boundary_i
675:       logical exitloop
676:       PetscBool flg

678:       limit=0
679:       maxits = 5
680:       tol=1e-10
681:       b=-0.5
682:       t= 0.5
683:       l=-0.5
684:       r= 0.5
685:       xt=0
686:       yt=0
687:       one=1.0
688:       two=2.0
689:       three=3.0

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

694:       bsize = xm + 2
695:       lsize = ym + 2
696:       rsize = ym + 2
697:       tsize = xm + 2

699:       PetscCall(VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr))
700:       PetscCall(VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr))
701:       PetscCall(VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr))
702:       PetscCall(VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr))

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

707:       do j=0,3

709:          if (j.eq.0) then
710:             yt=b
711:             xt=l+hx*xs
712:             limit=bsize
713:             PetscCall(VecGetArray(Bottom,boundary_v,boundary_i,ierr))

715:          elseif (j.eq.1) then
716:             yt=t
717:             xt=l+hx*xs
718:             limit=tsize
719:             PetscCall(VecGetArray(Top,boundary_v,boundary_i,ierr))

721:          elseif (j.eq.2) then
722:             yt=b+hy*ys
723:             xt=l
724:             limit=lsize
725:             PetscCall(VecGetArray(Left,boundary_v,boundary_i,ierr))

727:          elseif (j.eq.3) then
728:             yt=b+hy*ys
729:             xt=r
730:             limit=rsize
731:             PetscCall(VecGetArray(Right,boundary_v,boundary_i,ierr))
732:          endif

734:          do i=0,limit-1

736:             u1=xt
737:             u2=-yt
738:             k = 0
739:             exitloop = .false.
740:             do while (k .lt. maxits .and. (.not. exitloop))

742:                nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
743:                nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
744:                fnorm=sqrt(nf1*nf1+nf2*nf2)
745:                if (fnorm .gt. tol) then
746:                   njac11=one+u2*u2-u1*u1
747:                   njac12=two*u1*u2
748:                   njac21=-two*u1*u2
749:                   njac22=-one - u1*u1 + u2*u2
750:                   det = njac11*njac22-njac21*njac12
751:                   u1 = u1-(njac22*nf1-njac12*nf2)/det
752:                   u2 = u2-(njac11*nf2-njac21*nf1)/det
753:                else
754:                   exitloop = .true.
755:                endif
756:                k=k+1
757:             enddo

759:             boundary_v(i + boundary_i) = u1*u1-u2*u2
760:             if ((j .eq. 0) .or. (j .eq. 1)) then
761:                xt = xt + hx
762:             else
763:                yt = yt + hy
764:             endif

766:          enddo

768:          if (j.eq.0) then
769:             PetscCall(VecRestoreArray(Bottom,boundary_v,boundary_i,ierr))
770:          elseif (j.eq.1) then
771:             PetscCall(VecRestoreArray(Top,boundary_v,boundary_i,ierr))
772:          elseif (j.eq.2) then
773:             PetscCall(VecRestoreArray(Left,boundary_v,boundary_i,ierr))
774:          elseif (j.eq.3) then
775:             PetscCall(VecRestoreArray(Right,boundary_v,boundary_i,ierr))
776:          endif

778:       enddo

780: ! Scale the boundary if desired
781:       PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-bottom',scl,flg,ierr))
782:       if (flg .eqv. PETSC_TRUE) then
783:          PetscCall(VecScale(Bottom,scl,ierr))
784:       endif

786:       PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-top',scl,flg,ierr))
787:       if (flg .eqv. PETSC_TRUE) then
788:          PetscCall(VecScale(Top,scl,ierr))
789:       endif

791:       PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-right',scl,flg,ierr))
792:       if (flg .eqv. PETSC_TRUE) then
793:          PetscCall(VecScale(Right,scl,ierr))
794:       endif

796:       PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-left',scl,flg,ierr))
797:       if (flg .eqv. PETSC_TRUE) then
798:          PetscCall(VecScale(Left,scl,ierr))
799:       endif

801:       return
802:       end

804: ! ----------------------------------------------------------------------------
805: !
806: !/*
807: !     MSA_Plate - Calculates an obstacle for surface to stretch over
808: !
809: !     Output Parameter:
810: !.    xl - lower bound vector
811: !.    xu - upper bound vector
812: !
813: !*/

815:       subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
816:       use plate2fmodule
817:       implicit none

819:       Tao        tao
820:       Vec              xl,xu
821:       PetscErrorCode   ierr
822:       PetscInt         i,j,row
823:       PetscInt         xs, xm, ys, ym
824:       PetscReal      lb,ub
825:       PetscInt         dummy

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

835:       lb = PETSC_NINFINITY
836:       ub = PETSC_INFINITY

838:       if (bmy .lt. 0) bmy = 0
839:       if (bmy .gt. my) bmy = my
840:       if (bmx .lt. 0) bmx = 0
841:       if (bmx .gt. mx) bmx = mx

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

845:       PetscCall(VecSet(xl,lb,ierr))
846:       PetscCall(VecSet(xu,ub,ierr))

848:       PetscCall(VecGetArray(xl,xl_v,xl_i,ierr))

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

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

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

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

860:             endif

862:          enddo
863:       enddo

865:       PetscCall(VecRestoreArray(xl,xl_v,xl_i,ierr))

867:       return
868:       end

870: ! ----------------------------------------------------------------------------
871: !
872: !/*
873: !     MSA_InitialPoint - Calculates an obstacle for surface to stretch over
874: !
875: !     Output Parameter:
876: !.    X - vector for initial guess
877: !
878: !*/

880:       subroutine MSA_InitialPoint(X, ierr)
881:       use plate2fmodule
882:       implicit none

884:       Vec               X
885:       PetscErrorCode    ierr
886:       PetscInt          start,i,j
887:       PetscInt          row
888:       PetscInt          xs,xm,gxs,gxm
889:       PetscInt          ys,ym,gys,gym
890:       PetscReal       zero, np5

892: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
893: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
894: ! will return an array of doubles referenced by x_array offset by x_index.
895: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
896: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
897:       PetscReal   left_v(0:1),right_v(0:1)
898:       PetscReal   bottom_v(0:1),top_v(0:1)
899:       PetscReal   x_v(0:1)
900:       PetscOffset   left_i, right_i, top_i
901:       PetscOffset   bottom_i,x_i
902:       PetscBool     flg
903:       PetscRandom   rctx

905:       zero = 0.0
906:       np5 = -0.5

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

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

913:       elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then  ! random start -0.5 < xi < 0.5
914:          PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr))
915:          do i=0,start-1
916:             PetscCall(VecSetRandom(X,rctx,ierr))
917:          enddo

919:          PetscCall(PetscRandomDestroy(rctx,ierr))
920:          PetscCall(VecShift(X,np5,ierr))

922:       else   ! average of boundary conditions

924: !        Get Local mesh boundaries
925:          PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
926:          PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))

928: !        Get pointers to vector data
929:          PetscCall(VecGetArray(Top,top_v,top_i,ierr))
930:          PetscCall(VecGetArray(Bottom,bottom_v,bottom_i,ierr))
931:          PetscCall(VecGetArray(Left,left_v,left_i,ierr))
932:          PetscCall(VecGetArray(Right,right_v,right_i,ierr))

934:          PetscCall(VecGetArray(localX,x_v,x_i,ierr))

936: !        Perform local computations
937:          do  j=ys,ys+ym-1
938:             do i=xs,xs+xm-1
939:                row = (j-gys)*gxm  + (i-gxs)
940:                x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) +                  &
941:      &                          (i+1)*left_v(left_i+j-ys+1)/mx + (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
942:             enddo
943:          enddo

945: !        Restore vectors
946:          PetscCall(VecRestoreArray(localX,x_v,x_i,ierr))

948:          PetscCall(VecRestoreArray(Left,left_v,left_i,ierr))
949:          PetscCall(VecRestoreArray(Top,top_v,top_i,ierr))
950:          PetscCall(VecRestoreArray(Bottom,bottom_v,bottom_i,ierr))
951:          PetscCall(VecRestoreArray(Right,right_v,right_i,ierr))

953:          PetscCall(DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr))
954:          PetscCall(DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr))

956:       endif

958:       return
959:       end

961: !
962: !/*TEST
963: !
964: !   build:
965: !      requires: !complex
966: !
967: !   test:
968: !      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
969: !      filter: sort -b
970: !      filter_output: sort -b
971: !      requires: !single
972: !
973: !   test:
974: !      suffix: 2
975: !      nsize: 2
976: !      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
977: !      filter: sort -b
978: !      filter_output: sort -b
979: !      requires: !single
980: !
981: !TEST*/