Actual source code: plate2f.F
petsc-3.6.1 2015-08-06
1: ! Program usage: mpirun -np <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: Vec xl, xu ! lower and upper bounds vectorsp
47: PetscInt m ! number of local elements in vector
48: Tao tao ! Tao solver context
49: Mat H ! Hessian matrix
50: ISLocalToGlobalMapping isltog ! local to global mapping object
51: PetscBool flg
52: PetscInt i1,i3,i7
55: external FormFunctionGradient
56: external FormHessian
57: external MSA_BoundaryConditions
58: external MSA_Plate
59: external MSA_InitialPoint
60: ! Initialize Tao
62: i1=1
63: i3=3
64: i7=7
67: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
69: ! Specify default dimensions of the problem
70: mx = 10
71: my = 10
72: bheight = 0.1
74: ! Check for any command line arguments that override defaults
76: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mx",mx,flg,ierr)
77: CHKERRQ(ierr)
78: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,ierr)
79: CHKERRQ(ierr)
81: bmx = mx/2
82: bmy = my/2
84: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-bmx",bmx,flg,ierr)
85: CHKERRQ(ierr)
86: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-bmy",bmy,flg,ierr)
87: CHKERRQ(ierr)
88: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bheight",bheight, &
89: & flg,ierr)
90: CHKERRQ(ierr)
93: ! Calculate any derived values from parameters
94: N = mx*my
96: ! Let Petsc determine the dimensions of the local vectors
97: Nx = PETSC_DECIDE
98: NY = PETSC_DECIDE
100: ! A two dimensional distributed array will help define this problem, which
101: ! derives from an elliptic PDE on a two-dimensional domain. From the
102: ! distributed array, create the vectors
104: call DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE, &
105: & DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
106: & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
107: & dm,ierr)
108: CHKERRQ(ierr)
111: ! Extract global and local vectors from DM; The local vectors are
112: ! used solely as work space for the evaluation of the function,
113: ! gradient, and Hessian. Duplicate for remaining vectors that are
114: ! the same types.
116: call DMCreateGlobalVector(dm,x,ierr)
117: CHKERRQ(ierr)
118: call DMCreateLocalVector(dm,localX,ierr)
119: CHKERRQ(ierr)
120: call VecDuplicate(localX,localV,ierr)
121: CHKERRQ(ierr)
123: ! Create a matrix data structure to store the Hessian.
124: ! Here we (optionally) also associate the local numbering scheme
125: ! with the matrix so that later we can use local indices for matrix
126: ! assembly
128: call VecGetLocalSize(x,m,ierr)
129: CHKERRQ(ierr)
130: call MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, &
131: & i3,PETSC_NULL_INTEGER,H,ierr)
132: CHKERRQ(ierr)
134: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
135: CHKERRQ(ierr)
136: call DMGetLocalToGlobalMapping(dm,isltog,ierr)
137: CHKERRQ(ierr)
138: call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
139: CHKERRQ(ierr)
142: ! The Tao code begins here
143: ! Create TAO solver and set desired solution method.
144: ! This problems uses bounded variables, so the
145: ! method must either be 'tao_tron' or 'tao_blmvm'
147: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
148: CHKERRQ(ierr)
149: call TaoSetType(tao,TAOBLMVM,ierr)
150: CHKERRQ(ierr)
152: ! Set minimization function and gradient, hessian evaluation functions
154: call TaoSetObjectiveAndGradientRoutine(tao, &
155: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
156: CHKERRQ(ierr)
158: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
159: & PETSC_NULL_OBJECT, ierr)
160: CHKERRQ(ierr)
162: ! Set Variable bounds
163: call MSA_BoundaryConditions(ierr)
164: CHKERRQ(ierr)
165: call TaoSetVariableBoundsRoutine(tao,MSA_Plate, &
166: & PETSC_NULL_OBJECT,ierr)
167: CHKERRQ(ierr)
169: ! Set the initial solution guess
170: call MSA_InitialPoint(x, ierr)
171: CHKERRQ(ierr)
172: call TaoSetInitialVector(tao,x,ierr)
173: CHKERRQ(ierr)
175: ! Check for any tao command line options
176: call TaoSetFromOptions(tao,ierr)
177: CHKERRQ(ierr)
179: ! Solve the application
180: call TaoSolve(tao,ierr)
181: CHKERRQ(ierr)
183: ! Free TAO data structures
184: call TaoDestroy(tao,ierr)
185: CHKERRQ(ierr)
187: ! Free PETSc data structures
188: call VecDestroy(x,ierr)
189: call VecDestroy(Top,ierr)
190: call VecDestroy(Bottom,ierr)
191: call VecDestroy(Left,ierr)
192: call VecDestroy(Right,ierr)
193: call MatDestroy(H,ierr)
194: call VecDestroy(localX,ierr)
195: call VecDestroy(localV,ierr)
196: call DMDestroy(dm,ierr)
198: ! Finalize TAO
200: call PetscFinalize(ierr)
202: end
204: ! ---------------------------------------------------------------------
205: !
206: ! FormFunctionGradient - Evaluates function f(X).
207: !
208: ! Input Parameters:
209: ! tao - the Tao context
210: ! X - the input vector
211: ! dummy - optional user-defined context, as set by TaoSetFunction()
212: ! (not used here)
213: !
214: ! Output Parameters:
215: ! fcn - the newly evaluated function
216: ! G - the gradient vector
217: ! info - error code
218: !
221: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
222: implicit none
224: ! dm, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
225: #include "plate2f.h"
227: ! Input/output variables
229: Tao tao
230: PetscReal fcn
231: Vec X, G
232: PetscErrorCode ierr
233: PetscInt dummy
235: PetscInt i,j,row
236: PetscInt xs, xm
237: PetscInt gxs, gxm
238: PetscInt ys, ym
239: PetscInt gys, gym
240: PetscReal ft,zero,hx,hy,hydhx,hxdhy
241: PetscReal area,rhx,rhy
242: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
243: PetscReal d4,d5,d6,d7,d8
244: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
245: PetscReal df5dxc,df6dxc
246: PetscReal xc,xl,xr,xt,xb,xlt,xrb
249: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
250: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
251: ! will return an array of doubles referenced by x_array offset by x_index.
252: ! i.e., to reference the kth element of X, use x_array(k + x_index).
253: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
254: PetscReal g_v(0:1),x_v(0:1)
255: PetscReal top_v(0:1),left_v(0:1)
256: PetscReal right_v(0:1),bottom_v(0:1)
257: PetscOffset g_i,left_i,right_i
258: PetscOffset bottom_i,top_i,x_i
260: ft = 0.0
261: zero = 0.0
262: hx = 1.0/(mx + 1)
263: hy = 1.0/(my + 1)
264: hydhx = hy/hx
265: hxdhy = hx/hy
266: area = 0.5 * hx * hy
267: rhx = mx + 1.0
268: rhy = my + 1.0
271: ! Get local mesh boundaries
272: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
273: & PETSC_NULL_INTEGER,ierr)
274: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
275: & gxm,gym,PETSC_NULL_INTEGER,ierr)
277: ! Scatter ghost points to local vector
278: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
279: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
281: ! Initialize the vector to zero
282: call VecSet(localV,zero,ierr)
284: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
285: call VecGetArray(localX,x_v,x_i,ierr)
286: call VecGetArray(localV,g_v,g_i,ierr)
287: call VecGetArray(Top,top_v,top_i,ierr)
288: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
289: call VecGetArray(Left,left_v,left_i,ierr)
290: call VecGetArray(Right,right_v,right_i,ierr)
292: ! Compute function over the locally owned part of the mesh
293: do j = ys,ys+ym-1
294: do i = xs,xs+xm-1
295: row = (j-gys)*gxm + (i-gxs)
296: xc = x_v(row+x_i)
297: xt = xc
298: xb = xc
299: xr = xc
300: xl = xc
301: xrb = xc
302: xlt = xc
304: if (i .eq. 0) then !left side
305: xl = left_v(j - ys + 1 + left_i)
306: xlt = left_v(j - ys + 2 + left_i)
307: else
308: xl = x_v(row - 1 + x_i)
309: endif
311: if (j .eq. 0) then !bottom side
312: xb = bottom_v(i - xs + 1 + bottom_i)
313: xrb = bottom_v(i - xs + 2 + bottom_i)
314: else
315: xb = x_v(row - gxm + x_i)
316: endif
318: if (i + 1 .eq. gxs + gxm) then !right side
319: xr = right_v(j - ys + 1 + right_i)
320: xrb = right_v(j - ys + right_i)
321: else
322: xr = x_v(row + 1 + x_i)
323: endif
325: if (j + 1 .eq. gys + gym) then !top side
326: xt = top_v(i - xs + 1 + top_i)
327: xlt = top_v(i - xs + top_i)
328: else
329: xt = x_v(row + gxm + x_i)
330: endif
332: if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
333: xlt = x_v(row - 1 + gxm + x_i)
334: endif
336: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
337: xrb = x_v(row + 1 - gxm + x_i)
338: endif
340: d1 = xc-xl
341: d2 = xc-xr
342: d3 = xc-xt
343: d4 = xc-xb
344: d5 = xr-xrb
345: d6 = xrb-xb
346: d7 = xlt-xl
347: d8 = xt-xlt
349: df1dxc = d1 * hydhx
350: df2dxc = d1 * hydhx + d4 * hxdhy
351: df3dxc = d3 * hxdhy
352: df4dxc = d2 * hydhx + d3 * hxdhy
353: df5dxc = d2 * hydhx
354: df6dxc = d4 * hxdhy
356: d1 = d1 * rhx
357: d2 = d2 * rhx
358: d3 = d3 * rhy
359: d4 = d4 * rhy
360: d5 = d5 * rhy
361: d6 = d6 * rhx
362: d7 = d7 * rhy
363: d8 = d8 * rhx
365: f1 = sqrt(1.0 + d1*d1 + d7*d7)
366: f2 = sqrt(1.0 + d1*d1 + d4*d4)
367: f3 = sqrt(1.0 + d3*d3 + d8*d8)
368: f4 = sqrt(1.0 + d3*d3 + d2*d2)
369: f5 = sqrt(1.0 + d2*d2 + d5*d5)
370: f6 = sqrt(1.0 + d4*d4 + d6*d6)
372: ft = ft + f2 + f4
374: df1dxc = df1dxc / f1
375: df2dxc = df2dxc / f2
376: df3dxc = df3dxc / f3
377: df4dxc = df4dxc / f4
378: df5dxc = df5dxc / f5
379: df6dxc = df6dxc / f6
381: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
382: & df5dxc + df6dxc)
383: enddo
384: enddo
386: ! Compute triangular areas along the border of the domain.
387: if (xs .eq. 0) then ! left side
388: do j=ys,ys+ym-1
389: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
390: & * rhy
391: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
392: & * rhx
393: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
394: enddo
395: endif
398: if (ys .eq. 0) then !bottom side
399: do i=xs,xs+xm-1
400: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
401: & * rhx
402: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
403: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
404: enddo
405: endif
408: if (xs + xm .eq. mx) then ! right side
409: do j=ys,ys+ym-1
410: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
411: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
412: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
413: enddo
414: endif
417: if (ys + ym .eq. my) then
418: do i=xs,xs+xm-1
419: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
420: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
421: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
422: enddo
423: endif
426: if ((ys .eq. 0) .and. (xs .eq. 0)) then
427: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
428: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
429: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
430: endif
432: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
433: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
434: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
435: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
436: endif
438: ft = ft * area
439: call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, &
440: & MPIU_SUM,MPI_COMM_WORLD,ierr)
443: ! Restore vectors
444: call VecRestoreArray(localX,x_v,x_i,ierr)
445: call VecRestoreArray(localV,g_v,g_i,ierr)
446: call VecRestoreArray(Left,left_v,left_i,ierr)
447: call VecRestoreArray(Top,top_v,top_i,ierr)
448: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
449: call VecRestoreArray(Right,right_v,right_i,ierr)
451: ! Scatter values to global vector
452: call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
453: call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)
455: call PetscLogFlops(70.0d0*xm*ym,ierr)
457: return
458: end !FormFunctionGradient
464: ! ----------------------------------------------------------------------------
465: !
466: !
467: ! FormHessian - Evaluates Hessian matrix.
468: !
469: ! Input Parameters:
470: !. tao - the Tao context
471: !. X - input vector
472: !. dummy - not used
473: !
474: ! Output Parameters:
475: !. Hessian - Hessian matrix
476: !. Hpc - optionally different preconditioning matrix
477: !. flag - flag indicating matrix structure
478: !
479: ! Notes:
480: ! Due to mesh point reordering with DMs, we must always work
481: ! with the local mesh points, and then transform them to the new
482: ! global numbering with the local-to-global mapping. We cannot work
483: ! directly with the global numbers for the original uniprocessor mesh!
484: !
485: ! MatSetValuesLocal(), using the local ordering (including
486: ! ghost points!)
487: ! - Then set matrix entries using the local ordering
488: ! by calling MatSetValuesLocal()
490: subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
491: implicit none
493: ! dm,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
494: #include "plate2f.h"
496: Tao tao
497: Vec X
498: Mat Hessian,Hpc
499: PetscInt dummy
500: PetscErrorCode ierr
502: PetscInt i,j,k,row
503: PetscInt xs,xm,gxs,gxm
504: PetscInt ys,ym,gys,gym
505: PetscInt col(0:6)
506: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
507: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
508: PetscReal d4,d5,d6,d7,d8
509: PetscReal xc,xl,xr,xt,xb,xlt,xrb
510: PetscReal hl,hr,ht,hb,hc,htl,hbr
512: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
513: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
514: ! will return an array of doubles referenced by x_array offset by x_index.
515: ! i.e., to reference the kth element of X, use x_array(k + x_index).
516: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
517: PetscReal right_v(0:1),left_v(0:1)
518: PetscReal bottom_v(0:1),top_v(0:1)
519: PetscReal x_v(0:1)
520: PetscOffset x_i,right_i,left_i
521: PetscOffset bottom_i,top_i
522: PetscReal v(0:6)
523: PetscBool assembled
524: PetscInt i1
526: i1=1
528: ! Set various matrix options
529: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, &
530: & PETSC_TRUE,ierr)
532: ! Get local mesh boundaries
533: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
534: & PETSC_NULL_INTEGER,ierr)
535: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
536: & PETSC_NULL_INTEGER,ierr)
538: ! Scatter ghost points to local vectors
539: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
540: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
542: ! Get pointers to vector data (see note on Fortran arrays above)
543: call VecGetArray(localX,x_v,x_i,ierr)
544: call VecGetArray(Top,top_v,top_i,ierr)
545: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
546: call VecGetArray(Left,left_v,left_i,ierr)
547: call VecGetArray(Right,right_v,right_i,ierr)
549: ! Initialize matrix entries to zero
550: call MatAssembled(Hessian,assembled,ierr)
551: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)
554: rhx = mx + 1.0
555: rhy = my + 1.0
556: hx = 1.0/rhx
557: hy = 1.0/rhy
558: hydhx = hy/hx
559: hxdhy = hx/hy
560: ! compute Hessian over the locally owned part of the mesh
562: do i=xs,xs+xm-1
563: do j=ys,ys+ym-1
564: row = (j-gys)*gxm + (i-gxs)
566: xc = x_v(row + x_i)
567: xt = xc
568: xb = xc
569: xr = xc
570: xl = xc
571: xrb = xc
572: xlt = xc
574: if (i .eq. gxs) then ! Left side
575: xl = left_v(left_i + j - ys + 1)
576: xlt = left_v(left_i + j - ys + 2)
577: else
578: xl = x_v(x_i + row -1 )
579: endif
581: if (j .eq. gys) then ! bottom side
582: xb = bottom_v(bottom_i + i - xs + 1)
583: xrb = bottom_v(bottom_i + i - xs + 2)
584: else
585: xb = x_v(x_i + row - gxm)
586: endif
588: if (i+1 .eq. gxs + gxm) then !right side
589: xr = right_v(right_i + j - ys + 1)
590: xrb = right_v(right_i + j - ys)
591: else
592: xr = x_v(x_i + row + 1)
593: endif
595: if (j+1 .eq. gym+gys) then !top side
596: xt = top_v(top_i +i - xs + 1)
597: xlt = top_v(top_i + i - xs)
598: else
599: xt = x_v(x_i + row + gxm)
600: endif
602: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
603: xlt = x_v(x_i + row - 1 + gxm)
604: endif
606: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
607: xrb = x_v(x_i + row + 1 - gxm)
608: endif
610: d1 = (xc-xl)*rhx
611: d2 = (xc-xr)*rhx
612: d3 = (xc-xt)*rhy
613: d4 = (xc-xb)*rhy
614: d5 = (xrb-xr)*rhy
615: d6 = (xrb-xb)*rhx
616: d7 = (xlt-xl)*rhy
617: d8 = (xlt-xt)*rhx
619: f1 = sqrt( 1.0 + d1*d1 + d7*d7)
620: f2 = sqrt( 1.0 + d1*d1 + d4*d4)
621: f3 = sqrt( 1.0 + d3*d3 + d8*d8)
622: f4 = sqrt( 1.0 + d3*d3 + d2*d2)
623: f5 = sqrt( 1.0 + d2*d2 + d5*d5)
624: f6 = sqrt( 1.0 + d4*d4 + d6*d6)
627: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
628: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
630: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
631: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
633: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
634: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
636: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
637: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
639: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
640: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
642: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
643: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
644: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
645: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
646: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
647: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
648: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
650: hl = hl * 0.5
651: hr = hr * 0.5
652: ht = ht * 0.5
653: hb = hb * 0.5
654: hbr = hbr * 0.5
655: htl = htl * 0.5
656: hc = hc * 0.5
658: k = 0
660: if (j .gt. 0) then
661: v(k) = hb
662: col(k) = row - gxm
663: k=k+1
664: endif
666: if ((j .gt. 0) .and. (i .lt. mx-1)) then
667: v(k) = hbr
668: col(k) = row-gxm+1
669: k=k+1
670: endif
672: if (i .gt. 0) then
673: v(k) = hl
674: col(k) = row - 1
675: k = k+1
676: endif
678: v(k) = hc
679: col(k) = row
680: k=k+1
682: if (i .lt. mx-1) then
683: v(k) = hr
684: col(k) = row + 1
685: k=k+1
686: endif
688: if ((i .gt. 0) .and. (j .lt. my-1)) then
689: v(k) = htl
690: col(k) = row + gxm - 1
691: k=k+1
692: endif
694: if (j .lt. my-1) then
695: v(k) = ht
696: col(k) = row + gxm
697: k=k+1
698: endif
700: ! Set matrix values using local numbering, defined earlier in main routine
701: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
702: & ierr)
706: enddo
707: enddo
709: ! restore vectors
710: call VecRestoreArray(localX,x_v,x_i,ierr)
711: call VecRestoreArray(Left,left_v,left_i,ierr)
712: call VecRestoreArray(Right,right_v,right_i,ierr)
713: call VecRestoreArray(Top,top_v,top_i,ierr)
714: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
717: ! Assemble the matrix
718: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
719: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)
721: call PetscLogFlops(199.0d0*xm*ym,ierr)
723: return
724: end
730: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
732: ! ----------------------------------------------------------------------------
733: !
734: !/*
735: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
736: !
737: !
738: !*/
740: subroutine MSA_BoundaryConditions(ierr)
741: implicit none
743: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
744: #include "plate2f.h"
746: PetscErrorCode ierr
747: PetscInt i,j,k,limit,maxits
748: PetscInt xs, xm, gxs, gxm
749: PetscInt ys, ym, gys, gym
750: PetscInt bsize, lsize
751: PetscInt tsize, rsize
752: PetscReal one,two,three,tol
753: PetscReal scl,fnorm,det,xt
754: PetscReal yt,hx,hy,u1,u2,nf1,nf2
755: PetscReal njac11,njac12,njac21,njac22
756: PetscReal b, t, l, r
757: PetscReal boundary_v(0:1)
758: PetscOffset boundary_i
759: logical exitloop
760: PetscBool flg
762: limit=0
763: maxits = 5
764: tol=1e-10
765: b=-0.5
766: t= 0.5
767: l=-0.5
768: r= 0.5
769: xt=0
770: yt=0
771: one=1.0
772: two=2.0
773: three=3.0
776: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
777: & PETSC_NULL_INTEGER,ierr)
778: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
779: & gxm,gym,PETSC_NULL_INTEGER,ierr)
781: bsize = xm + 2
782: lsize = ym + 2
783: rsize = ym + 2
784: tsize = xm + 2
787: call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
788: call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
789: call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
790: call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)
792: hx= (r-l)/(mx+1)
793: hy= (t-b)/(my+1)
795: do j=0,3
797: if (j.eq.0) then
798: yt=b
799: xt=l+hx*xs
800: limit=bsize
801: call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
804: elseif (j.eq.1) then
805: yt=t
806: xt=l+hx*xs
807: limit=tsize
808: call VecGetArray(Top,boundary_v,boundary_i,ierr)
810: elseif (j.eq.2) then
811: yt=b+hy*ys
812: xt=l
813: limit=lsize
814: call VecGetArray(Left,boundary_v,boundary_i,ierr)
816: elseif (j.eq.3) then
817: yt=b+hy*ys
818: xt=r
819: limit=rsize
820: call VecGetArray(Right,boundary_v,boundary_i,ierr)
821: endif
824: do i=0,limit-1
826: u1=xt
827: u2=-yt
828: k = 0
829: exitloop = .false.
830: do while (k .lt. maxits .and. (.not. exitloop) )
832: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
833: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
834: fnorm=sqrt(nf1*nf1+nf2*nf2)
835: if (fnorm .gt. tol) then
836: njac11=one+u2*u2-u1*u1
837: njac12=two*u1*u2
838: njac21=-two*u1*u2
839: njac22=-one - u1*u1 + u2*u2
840: det = njac11*njac22-njac21*njac12
841: u1 = u1-(njac22*nf1-njac12*nf2)/det
842: u2 = u2-(njac11*nf2-njac21*nf1)/det
843: else
844: exitloop = .true.
845: endif
846: k=k+1
847: enddo
849: boundary_v(i + boundary_i) = u1*u1-u2*u2
850: if ((j .eq. 0) .or. (j .eq. 1)) then
851: xt = xt + hx
852: else
853: yt = yt + hy
854: endif
856: enddo
859: if (j.eq.0) then
860: call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
861: elseif (j.eq.1) then
862: call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
863: elseif (j.eq.2) then
864: call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
865: elseif (j.eq.3) then
866: call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
867: endif
869: enddo
872: ! Scale the boundary if desired
873: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom", &
874: & scl,flg,ierr)
875: if (flg .eqv. PETSC_TRUE) then
876: call VecScale(scl,Bottom,ierr)
877: endif
879: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top", &
880: & scl,flg,ierr)
881: if (flg .eqv. PETSC_TRUE) then
882: call VecScale(scl,Top,ierr)
883: endif
885: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right", &
886: & scl,flg,ierr)
887: if (flg .eqv. PETSC_TRUE) then
888: call VecScale(scl,Right,ierr)
889: endif
891: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left", &
892: & scl,flg,ierr)
893: if (flg .eqv. PETSC_TRUE) then
894: call VecScale(scl,Left,ierr)
895: endif
898: return
899: end
901: ! ----------------------------------------------------------------------------
902: !
903: !/*
904: ! MSA_Plate - Calculates an obstacle for surface to stretch over
905: !
906: ! Output Parameter:
907: !. xl - lower bound vector
908: !. xu - upper bound vector
909: !
910: !*/
912: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
913: implicit none
914: ! mx,my,bmx,bmy,dm,bheight defined in plate2f.h
915: #include "plate2f.h"
916: Tao tao
917: Vec xl,xu
918: PetscErrorCode ierr
919: PetscInt i,j,row
920: PetscInt xs, xm, ys, ym
921: PetscReal lb,ub
922: PetscInt dummy
924: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
925: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
926: ! will return an array of doubles referenced by x_array offset by x_index.
927: ! i.e., to reference the kth element of X, use x_array(k + x_index).
928: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
929: PetscReal xl_v(0:1)
930: PetscOffset xl_i
932: print *,'msa_plate'
934: lb = PETSC_NINFINITY
935: ub = PETSC_INFINITY
937: if (bmy .lt. 0) bmy = 0
938: if (bmy .gt. my) bmy = my
939: if (bmx .lt. 0) bmx = 0
940: if (bmx .gt. mx) bmx = mx
943: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
944: & PETSC_NULL_INTEGER,ierr)
946: call VecSet(xl,lb,ierr)
947: call VecSet(xu,ub,ierr)
949: call VecGetArray(xl,xl_v,xl_i,ierr)
952: do i=xs,xs+xm-1
954: do j=ys,ys+ym-1
956: row=(j-ys)*xm + (i-xs)
958: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
959: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
960: xl_v(xl_i+row) = bheight
962: endif
964: enddo
965: enddo
968: call VecRestoreArray(xl,xl_v,xl_i,ierr)
970: return
971: end
977: ! ----------------------------------------------------------------------------
978: !
979: !/*
980: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
981: !
982: ! Output Parameter:
983: !. X - vector for initial guess
984: !
985: !*/
987: subroutine MSA_InitialPoint(X, ierr)
988: implicit none
990: ! mx,my,localX,dm,Top,Left,Bottom,Right defined in plate2f.h
991: #include "plate2f.h"
992: Vec X
993: PetscErrorCode ierr
994: PetscInt start,i,j
995: PetscInt row
996: PetscInt xs,xm,gxs,gxm
997: PetscInt ys,ym,gys,gym
998: PetscReal zero, np5
1000: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1001: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
1002: ! will return an array of doubles referenced by x_array offset by x_index.
1003: ! i.e., to reference the kth element of X, use x_array(k + x_index).
1004: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1005: PetscReal left_v(0:1),right_v(0:1)
1006: PetscReal bottom_v(0:1),top_v(0:1)
1007: PetscReal x_v(0:1)
1008: PetscOffset left_i, right_i, top_i
1009: PetscOffset bottom_i,x_i
1010: PetscBool flg
1011: PetscRandom rctx
1013: zero = 0.0
1014: np5 = -0.5
1016: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start", &
1017: & start,flg,ierr)
1019: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
1020: call VecSet(X,zero,ierr)
1022: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
1023: call PetscRandomCreate(MPI_COMM_WORLD,rctx,ierr)
1024: do i=0,start-1
1025: call VecSetRandom(X,rctx,ierr)
1026: enddo
1028: call PetscRandomDestroy(rctx,ierr)
1029: call VecShift(X,np5,ierr)
1031: else ! average of boundary conditions
1033: ! Get Local mesh boundaries
1034: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1035: & PETSC_NULL_INTEGER,ierr)
1036: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1037: & PETSC_NULL_INTEGER,ierr)
1041: ! Get pointers to vector data
1042: call VecGetArray(Top,top_v,top_i,ierr)
1043: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1044: call VecGetArray(Left,left_v,left_i,ierr)
1045: call VecGetArray(Right,right_v,right_i,ierr)
1047: call VecGetArray(localX,x_v,x_i,ierr)
1049: ! Perform local computations
1050: do j=ys,ys+ym-1
1051: do i=xs,xs+xm-1
1052: row = (j-gys)*gxm + (i-gxs)
1053: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1054: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1055: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1056: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1057: enddo
1058: enddo
1060: ! Restore vectors
1061: call VecRestoreArray(localX,x_v,x_i,ierr)
1063: call VecRestoreArray(Left,left_v,left_i,ierr)
1064: call VecRestoreArray(Top,top_v,top_i,ierr)
1065: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1066: call VecRestoreArray(Right,right_v,right_i,ierr)
1068: call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1069: call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)
1071: endif
1073: return
1074: end