Actual source code: plate2f.F
petsc-3.7.5 2017-01-01
1: ! Program usage: mpiexec -n <proc> plate2f [all TAO options]
2: !
3: ! This example demonstrates use of the TAO package to solve a bound constrained
4: ! minimization problem. This example is based on a problem from the
5: ! MINPACK-2 test suite. Given a rectangular 2-D domain and boundary values
6: ! along the edges of the domain, the objective is to find the surface
7: ! with the minimal area that satisfies the boundary conditions.
8: ! The command line options are:
9: ! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
10: ! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
11: ! -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
12: ! -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
13: ! -bheight <ht>, where <ht> = height of the plate
14: !
15: !/*T
16: ! Concepts: TAO^Solving a bound constrained minimization problem
17: ! Routines: TaoCreate();
18: ! Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
19: ! Routines: TaoSetHessianRoutine();
20: ! Routines: TaoSetVariableBoundsRoutine();
21: ! Routines: TaoSetInitialVector();
22: ! Routines: TaoSetFromOptions();
23: ! Routines: TaoSolve();
24: ! Routines: TaoDestroy();
25: ! Processors: n
26: !T*/
30: implicit none
32: #include "plate2f.h"
34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: ! Variable declarations
36: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37: !
38: ! Variables:
39: ! (common from plate2f.h):
40: ! Nx, Ny number of processors in x- and y- directions
41: ! mx, my number of grid points in x,y directions
42: ! N global dimension of vector
44: PetscErrorCode ierr ! used to check for functions returning nonzeros
45: Vec x ! solution vector
46: PetscInt m ! number of local elements in vector
47: Tao tao ! Tao solver context
48: Mat H ! Hessian matrix
49: ISLocalToGlobalMapping isltog ! local to global mapping object
50: PetscBool flg
51: PetscInt i1,i3,i7
54: external FormFunctionGradient
55: external FormHessian
56: external MSA_BoundaryConditions
57: external MSA_Plate
58: external MSA_InitialPoint
59: ! Initialize Tao
61: i1=1
62: i3=3
63: i7=7
66: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
68: ! Specify default dimensions of the problem
69: mx = 10
70: my = 10
71: bheight = 0.1
73: ! Check for any command line arguments that override defaults
75: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
76: & '-mx',mx,flg,ierr)
77: CHKERRQ(ierr)
78: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
79: & '-my',my,flg,ierr)
80: CHKERRQ(ierr)
82: bmx = mx/2
83: bmy = my/2
85: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
86: & '-bmx',bmx,flg,ierr)
87: CHKERRQ(ierr)
88: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
89: & '-bmy',bmy,flg,ierr)
90: CHKERRQ(ierr)
91: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
92: & '-bheight',bheight,flg,ierr)
93: CHKERRQ(ierr)
96: ! Calculate any derived values from parameters
97: N = mx*my
99: ! Let Petsc determine the dimensions of the local vectors
100: Nx = PETSC_DECIDE
101: NY = PETSC_DECIDE
103: ! A two dimensional distributed array will help define this problem, which
104: ! derives from an elliptic PDE on a two-dimensional domain. From the
105: ! distributed array, create the vectors
107: call DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE, &
108: & DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
109: & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
110: & dm,ierr)
111: CHKERRQ(ierr)
114: ! Extract global and local vectors from DM; The local vectors are
115: ! used solely as work space for the evaluation of the function,
116: ! gradient, and Hessian. Duplicate for remaining vectors that are
117: ! the same types.
119: call DMCreateGlobalVector(dm,x,ierr)
120: CHKERRQ(ierr)
121: call DMCreateLocalVector(dm,localX,ierr)
122: CHKERRQ(ierr)
123: call VecDuplicate(localX,localV,ierr)
124: CHKERRQ(ierr)
126: ! Create a matrix data structure to store the Hessian.
127: ! Here we (optionally) also associate the local numbering scheme
128: ! with the matrix so that later we can use local indices for matrix
129: ! assembly
131: call VecGetLocalSize(x,m,ierr)
132: CHKERRQ(ierr)
133: call MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, &
134: & i3,PETSC_NULL_INTEGER,H,ierr)
135: CHKERRQ(ierr)
137: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
138: CHKERRQ(ierr)
139: call DMGetLocalToGlobalMapping(dm,isltog,ierr)
140: CHKERRQ(ierr)
141: call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
142: CHKERRQ(ierr)
145: ! The Tao code begins here
146: ! Create TAO solver and set desired solution method.
147: ! This problems uses bounded variables, so the
148: ! method must either be 'tao_tron' or 'tao_blmvm'
150: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
151: CHKERRQ(ierr)
152: call TaoSetType(tao,TAOBLMVM,ierr)
153: CHKERRQ(ierr)
155: ! Set minimization function and gradient, hessian evaluation functions
157: call TaoSetObjectiveAndGradientRoutine(tao, &
158: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
159: CHKERRQ(ierr)
161: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
162: & PETSC_NULL_OBJECT, ierr)
163: CHKERRQ(ierr)
165: ! Set Variable bounds
166: call MSA_BoundaryConditions(ierr)
167: CHKERRQ(ierr)
168: call TaoSetVariableBoundsRoutine(tao,MSA_Plate, &
169: & PETSC_NULL_OBJECT,ierr)
170: CHKERRQ(ierr)
172: ! Set the initial solution guess
173: call MSA_InitialPoint(x, ierr)
174: CHKERRQ(ierr)
175: call TaoSetInitialVector(tao,x,ierr)
176: CHKERRQ(ierr)
178: ! Check for any tao command line options
179: call TaoSetFromOptions(tao,ierr)
180: CHKERRQ(ierr)
182: ! Solve the application
183: call TaoSolve(tao,ierr)
184: CHKERRQ(ierr)
186: ! Free TAO data structures
187: call TaoDestroy(tao,ierr)
188: CHKERRQ(ierr)
190: ! Free PETSc data structures
191: call VecDestroy(x,ierr)
192: call VecDestroy(Top,ierr)
193: call VecDestroy(Bottom,ierr)
194: call VecDestroy(Left,ierr)
195: call VecDestroy(Right,ierr)
196: call MatDestroy(H,ierr)
197: call VecDestroy(localX,ierr)
198: call VecDestroy(localV,ierr)
199: call DMDestroy(dm,ierr)
201: ! Finalize TAO
203: call PetscFinalize(ierr)
205: end
207: ! ---------------------------------------------------------------------
208: !
209: ! FormFunctionGradient - Evaluates function f(X).
210: !
211: ! Input Parameters:
212: ! tao - the Tao context
213: ! X - the input vector
214: ! dummy - optional user-defined context, as set by TaoSetFunction()
215: ! (not used here)
216: !
217: ! Output Parameters:
218: ! fcn - the newly evaluated function
219: ! G - the gradient vector
220: ! info - error code
221: !
224: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
225: implicit none
227: ! dm, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
228: #include "plate2f.h"
230: ! Input/output variables
232: Tao tao
233: PetscReal fcn
234: Vec X, G
235: PetscErrorCode ierr
236: PetscInt dummy
238: PetscInt i,j,row
239: PetscInt xs, xm
240: PetscInt gxs, gxm
241: PetscInt ys, ym
242: PetscInt gys, gym
243: PetscReal ft,zero,hx,hy,hydhx,hxdhy
244: PetscReal area,rhx,rhy
245: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
246: PetscReal d4,d5,d6,d7,d8
247: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
248: PetscReal df5dxc,df6dxc
249: PetscReal xc,xl,xr,xt,xb,xlt,xrb
252: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
253: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
254: ! will return an array of doubles referenced by x_array offset by x_index.
255: ! i.e., to reference the kth element of X, use x_array(k + x_index).
256: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
257: PetscReal g_v(0:1),x_v(0:1)
258: PetscReal top_v(0:1),left_v(0:1)
259: PetscReal right_v(0:1),bottom_v(0:1)
260: PetscOffset g_i,left_i,right_i
261: PetscOffset bottom_i,top_i,x_i
263: ft = 0.0
264: zero = 0.0
265: hx = 1.0/(mx + 1)
266: hy = 1.0/(my + 1)
267: hydhx = hy/hx
268: hxdhy = hx/hy
269: area = 0.5 * hx * hy
270: rhx = mx + 1.0
271: rhy = my + 1.0
274: ! Get local mesh boundaries
275: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
276: & PETSC_NULL_INTEGER,ierr)
277: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
278: & gxm,gym,PETSC_NULL_INTEGER,ierr)
280: ! Scatter ghost points to local vector
281: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
282: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
284: ! Initialize the vector to zero
285: call VecSet(localV,zero,ierr)
287: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
288: call VecGetArray(localX,x_v,x_i,ierr)
289: call VecGetArray(localV,g_v,g_i,ierr)
290: call VecGetArray(Top,top_v,top_i,ierr)
291: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
292: call VecGetArray(Left,left_v,left_i,ierr)
293: call VecGetArray(Right,right_v,right_i,ierr)
295: ! Compute function over the locally owned part of the mesh
296: do j = ys,ys+ym-1
297: do i = xs,xs+xm-1
298: row = (j-gys)*gxm + (i-gxs)
299: xc = x_v(row+x_i)
300: xt = xc
301: xb = xc
302: xr = xc
303: xl = xc
304: xrb = xc
305: xlt = xc
307: if (i .eq. 0) then !left side
308: xl = left_v(j - ys + 1 + left_i)
309: xlt = left_v(j - ys + 2 + left_i)
310: else
311: xl = x_v(row - 1 + x_i)
312: endif
314: if (j .eq. 0) then !bottom side
315: xb = bottom_v(i - xs + 1 + bottom_i)
316: xrb = bottom_v(i - xs + 2 + bottom_i)
317: else
318: xb = x_v(row - gxm + x_i)
319: endif
321: if (i + 1 .eq. gxs + gxm) then !right side
322: xr = right_v(j - ys + 1 + right_i)
323: xrb = right_v(j - ys + right_i)
324: else
325: xr = x_v(row + 1 + x_i)
326: endif
328: if (j + 1 .eq. gys + gym) then !top side
329: xt = top_v(i - xs + 1 + top_i)
330: xlt = top_v(i - xs + top_i)
331: else
332: xt = x_v(row + gxm + x_i)
333: endif
335: if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
336: xlt = x_v(row - 1 + gxm + x_i)
337: endif
339: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
340: xrb = x_v(row + 1 - gxm + x_i)
341: endif
343: d1 = xc-xl
344: d2 = xc-xr
345: d3 = xc-xt
346: d4 = xc-xb
347: d5 = xr-xrb
348: d6 = xrb-xb
349: d7 = xlt-xl
350: d8 = xt-xlt
352: df1dxc = d1 * hydhx
353: df2dxc = d1 * hydhx + d4 * hxdhy
354: df3dxc = d3 * hxdhy
355: df4dxc = d2 * hydhx + d3 * hxdhy
356: df5dxc = d2 * hydhx
357: df6dxc = d4 * hxdhy
359: d1 = d1 * rhx
360: d2 = d2 * rhx
361: d3 = d3 * rhy
362: d4 = d4 * rhy
363: d5 = d5 * rhy
364: d6 = d6 * rhx
365: d7 = d7 * rhy
366: d8 = d8 * rhx
368: f1 = sqrt(1.0 + d1*d1 + d7*d7)
369: f2 = sqrt(1.0 + d1*d1 + d4*d4)
370: f3 = sqrt(1.0 + d3*d3 + d8*d8)
371: f4 = sqrt(1.0 + d3*d3 + d2*d2)
372: f5 = sqrt(1.0 + d2*d2 + d5*d5)
373: f6 = sqrt(1.0 + d4*d4 + d6*d6)
375: ft = ft + f2 + f4
377: df1dxc = df1dxc / f1
378: df2dxc = df2dxc / f2
379: df3dxc = df3dxc / f3
380: df4dxc = df4dxc / f4
381: df5dxc = df5dxc / f5
382: df6dxc = df6dxc / f6
384: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
385: & df5dxc + df6dxc)
386: enddo
387: enddo
389: ! Compute triangular areas along the border of the domain.
390: if (xs .eq. 0) then ! left side
391: do j=ys,ys+ym-1
392: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
393: & * rhy
394: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
395: & * rhx
396: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
397: enddo
398: endif
401: if (ys .eq. 0) then !bottom side
402: do i=xs,xs+xm-1
403: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
404: & * rhx
405: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
406: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
407: enddo
408: endif
411: if (xs + xm .eq. mx) then ! right side
412: do j=ys,ys+ym-1
413: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
414: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
415: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
416: enddo
417: endif
420: if (ys + ym .eq. my) then
421: do i=xs,xs+xm-1
422: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
423: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
424: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
425: enddo
426: endif
429: if ((ys .eq. 0) .and. (xs .eq. 0)) then
430: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
431: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
432: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
433: endif
435: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
436: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
437: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
438: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
439: endif
441: ft = ft * area
442: call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, &
443: & MPIU_SUM,MPI_COMM_WORLD,ierr)
446: ! Restore vectors
447: call VecRestoreArray(localX,x_v,x_i,ierr)
448: call VecRestoreArray(localV,g_v,g_i,ierr)
449: call VecRestoreArray(Left,left_v,left_i,ierr)
450: call VecRestoreArray(Top,top_v,top_i,ierr)
451: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
452: call VecRestoreArray(Right,right_v,right_i,ierr)
454: ! Scatter values to global vector
455: call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
456: call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)
458: call PetscLogFlops(70.0d0*xm*ym,ierr)
460: return
461: end !FormFunctionGradient
467: ! ----------------------------------------------------------------------------
468: !
469: !
470: ! FormHessian - Evaluates Hessian matrix.
471: !
472: ! Input Parameters:
473: !. tao - the Tao context
474: !. X - input vector
475: !. dummy - not used
476: !
477: ! Output Parameters:
478: !. Hessian - Hessian matrix
479: !. Hpc - optionally different preconditioning matrix
480: !. flag - flag indicating matrix structure
481: !
482: ! Notes:
483: ! Due to mesh point reordering with DMs, we must always work
484: ! with the local mesh points, and then transform them to the new
485: ! global numbering with the local-to-global mapping. We cannot work
486: ! directly with the global numbers for the original uniprocessor mesh!
487: !
488: ! MatSetValuesLocal(), using the local ordering (including
489: ! ghost points!)
490: ! - Then set matrix entries using the local ordering
491: ! by calling MatSetValuesLocal()
493: subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
494: implicit none
496: ! dm,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
497: #include "plate2f.h"
499: Tao tao
500: Vec X
501: Mat Hessian,Hpc
502: PetscInt dummy
503: PetscErrorCode ierr
505: PetscInt i,j,k,row
506: PetscInt xs,xm,gxs,gxm
507: PetscInt ys,ym,gys,gym
508: PetscInt col(0:6)
509: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
510: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
511: PetscReal d4,d5,d6,d7,d8
512: PetscReal xc,xl,xr,xt,xb,xlt,xrb
513: PetscReal hl,hr,ht,hb,hc,htl,hbr
515: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
516: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
517: ! will return an array of doubles referenced by x_array offset by x_index.
518: ! i.e., to reference the kth element of X, use x_array(k + x_index).
519: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
520: PetscReal right_v(0:1),left_v(0:1)
521: PetscReal bottom_v(0:1),top_v(0:1)
522: PetscReal x_v(0:1)
523: PetscOffset x_i,right_i,left_i
524: PetscOffset bottom_i,top_i
525: PetscReal v(0:6)
526: PetscBool assembled
527: PetscInt i1
529: i1=1
531: ! Set various matrix options
532: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, &
533: & PETSC_TRUE,ierr)
535: ! Get local mesh boundaries
536: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
537: & PETSC_NULL_INTEGER,ierr)
538: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
539: & PETSC_NULL_INTEGER,ierr)
541: ! Scatter ghost points to local vectors
542: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
543: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
545: ! Get pointers to vector data (see note on Fortran arrays above)
546: call VecGetArray(localX,x_v,x_i,ierr)
547: call VecGetArray(Top,top_v,top_i,ierr)
548: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
549: call VecGetArray(Left,left_v,left_i,ierr)
550: call VecGetArray(Right,right_v,right_i,ierr)
552: ! Initialize matrix entries to zero
553: call MatAssembled(Hessian,assembled,ierr)
554: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)
557: rhx = mx + 1.0
558: rhy = my + 1.0
559: hx = 1.0/rhx
560: hy = 1.0/rhy
561: hydhx = hy/hx
562: hxdhy = hx/hy
563: ! compute Hessian over the locally owned part of the mesh
565: do i=xs,xs+xm-1
566: do j=ys,ys+ym-1
567: row = (j-gys)*gxm + (i-gxs)
569: xc = x_v(row + x_i)
570: xt = xc
571: xb = xc
572: xr = xc
573: xl = xc
574: xrb = xc
575: xlt = xc
577: if (i .eq. gxs) then ! Left side
578: xl = left_v(left_i + j - ys + 1)
579: xlt = left_v(left_i + j - ys + 2)
580: else
581: xl = x_v(x_i + row -1 )
582: endif
584: if (j .eq. gys) then ! bottom side
585: xb = bottom_v(bottom_i + i - xs + 1)
586: xrb = bottom_v(bottom_i + i - xs + 2)
587: else
588: xb = x_v(x_i + row - gxm)
589: endif
591: if (i+1 .eq. gxs + gxm) then !right side
592: xr = right_v(right_i + j - ys + 1)
593: xrb = right_v(right_i + j - ys)
594: else
595: xr = x_v(x_i + row + 1)
596: endif
598: if (j+1 .eq. gym+gys) then !top side
599: xt = top_v(top_i +i - xs + 1)
600: xlt = top_v(top_i + i - xs)
601: else
602: xt = x_v(x_i + row + gxm)
603: endif
605: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
606: xlt = x_v(x_i + row - 1 + gxm)
607: endif
609: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
610: xrb = x_v(x_i + row + 1 - gxm)
611: endif
613: d1 = (xc-xl)*rhx
614: d2 = (xc-xr)*rhx
615: d3 = (xc-xt)*rhy
616: d4 = (xc-xb)*rhy
617: d5 = (xrb-xr)*rhy
618: d6 = (xrb-xb)*rhx
619: d7 = (xlt-xl)*rhy
620: d8 = (xlt-xt)*rhx
622: f1 = sqrt( 1.0 + d1*d1 + d7*d7)
623: f2 = sqrt( 1.0 + d1*d1 + d4*d4)
624: f3 = sqrt( 1.0 + d3*d3 + d8*d8)
625: f4 = sqrt( 1.0 + d3*d3 + d2*d2)
626: f5 = sqrt( 1.0 + d2*d2 + d5*d5)
627: f6 = sqrt( 1.0 + d4*d4 + d6*d6)
630: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
631: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
633: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
634: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
636: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
637: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
639: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
640: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
642: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
643: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
645: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
646: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
647: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
648: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
649: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
650: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
651: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
653: hl = hl * 0.5
654: hr = hr * 0.5
655: ht = ht * 0.5
656: hb = hb * 0.5
657: hbr = hbr * 0.5
658: htl = htl * 0.5
659: hc = hc * 0.5
661: k = 0
663: if (j .gt. 0) then
664: v(k) = hb
665: col(k) = row - gxm
666: k=k+1
667: endif
669: if ((j .gt. 0) .and. (i .lt. mx-1)) then
670: v(k) = hbr
671: col(k) = row-gxm+1
672: k=k+1
673: endif
675: if (i .gt. 0) then
676: v(k) = hl
677: col(k) = row - 1
678: k = k+1
679: endif
681: v(k) = hc
682: col(k) = row
683: k=k+1
685: if (i .lt. mx-1) then
686: v(k) = hr
687: col(k) = row + 1
688: k=k+1
689: endif
691: if ((i .gt. 0) .and. (j .lt. my-1)) then
692: v(k) = htl
693: col(k) = row + gxm - 1
694: k=k+1
695: endif
697: if (j .lt. my-1) then
698: v(k) = ht
699: col(k) = row + gxm
700: k=k+1
701: endif
703: ! Set matrix values using local numbering, defined earlier in main routine
704: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
705: & ierr)
709: enddo
710: enddo
712: ! restore vectors
713: call VecRestoreArray(localX,x_v,x_i,ierr)
714: call VecRestoreArray(Left,left_v,left_i,ierr)
715: call VecRestoreArray(Right,right_v,right_i,ierr)
716: call VecRestoreArray(Top,top_v,top_i,ierr)
717: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
720: ! Assemble the matrix
721: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
722: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)
724: call PetscLogFlops(199.0d0*xm*ym,ierr)
726: return
727: end
733: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
735: ! ----------------------------------------------------------------------------
736: !
737: !/*
738: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
739: !
740: !
741: !*/
743: subroutine MSA_BoundaryConditions(ierr)
744: implicit none
746: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
747: #include "plate2f.h"
749: PetscErrorCode ierr
750: PetscInt i,j,k,limit,maxits
751: PetscInt xs, xm, gxs, gxm
752: PetscInt ys, ym, gys, gym
753: PetscInt bsize, lsize
754: PetscInt tsize, rsize
755: PetscReal one,two,three,tol
756: PetscReal scl,fnorm,det,xt
757: PetscReal yt,hx,hy,u1,u2,nf1,nf2
758: PetscReal njac11,njac12,njac21,njac22
759: PetscReal b, t, l, r
760: PetscReal boundary_v(0:1)
761: PetscOffset boundary_i
762: logical exitloop
763: PetscBool flg
765: limit=0
766: maxits = 5
767: tol=1e-10
768: b=-0.5
769: t= 0.5
770: l=-0.5
771: r= 0.5
772: xt=0
773: yt=0
774: one=1.0
775: two=2.0
776: three=3.0
779: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
780: & PETSC_NULL_INTEGER,ierr)
781: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
782: & gxm,gym,PETSC_NULL_INTEGER,ierr)
784: bsize = xm + 2
785: lsize = ym + 2
786: rsize = ym + 2
787: tsize = xm + 2
790: call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
791: call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
792: call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
793: call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)
795: hx= (r-l)/(mx+1)
796: hy= (t-b)/(my+1)
798: do j=0,3
800: if (j.eq.0) then
801: yt=b
802: xt=l+hx*xs
803: limit=bsize
804: call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
807: elseif (j.eq.1) then
808: yt=t
809: xt=l+hx*xs
810: limit=tsize
811: call VecGetArray(Top,boundary_v,boundary_i,ierr)
813: elseif (j.eq.2) then
814: yt=b+hy*ys
815: xt=l
816: limit=lsize
817: call VecGetArray(Left,boundary_v,boundary_i,ierr)
819: elseif (j.eq.3) then
820: yt=b+hy*ys
821: xt=r
822: limit=rsize
823: call VecGetArray(Right,boundary_v,boundary_i,ierr)
824: endif
827: do i=0,limit-1
829: u1=xt
830: u2=-yt
831: k = 0
832: exitloop = .false.
833: do while (k .lt. maxits .and. (.not. exitloop) )
835: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
836: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
837: fnorm=sqrt(nf1*nf1+nf2*nf2)
838: if (fnorm .gt. tol) then
839: njac11=one+u2*u2-u1*u1
840: njac12=two*u1*u2
841: njac21=-two*u1*u2
842: njac22=-one - u1*u1 + u2*u2
843: det = njac11*njac22-njac21*njac12
844: u1 = u1-(njac22*nf1-njac12*nf2)/det
845: u2 = u2-(njac11*nf2-njac21*nf1)/det
846: else
847: exitloop = .true.
848: endif
849: k=k+1
850: enddo
852: boundary_v(i + boundary_i) = u1*u1-u2*u2
853: if ((j .eq. 0) .or. (j .eq. 1)) then
854: xt = xt + hx
855: else
856: yt = yt + hy
857: endif
859: enddo
862: if (j.eq.0) then
863: call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
864: elseif (j.eq.1) then
865: call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
866: elseif (j.eq.2) then
867: call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
868: elseif (j.eq.3) then
869: call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
870: endif
872: enddo
875: ! Scale the boundary if desired
876: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
877: & '-bottom',scl,flg,ierr)
878: if (flg .eqv. PETSC_TRUE) then
879: call VecScale(scl,Bottom,ierr)
880: endif
882: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
883: & '-top',scl,flg,ierr)
884: if (flg .eqv. PETSC_TRUE) then
885: call VecScale(scl,Top,ierr)
886: endif
888: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
889: & '-right',scl,flg,ierr)
890: if (flg .eqv. PETSC_TRUE) then
891: call VecScale(scl,Right,ierr)
892: endif
894: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
895: & '-left',scl,flg,ierr)
896: if (flg .eqv. PETSC_TRUE) then
897: call VecScale(scl,Left,ierr)
898: endif
901: return
902: end
904: ! ----------------------------------------------------------------------------
905: !
906: !/*
907: ! MSA_Plate - Calculates an obstacle for surface to stretch over
908: !
909: ! Output Parameter:
910: !. xl - lower bound vector
911: !. xu - upper bound vector
912: !
913: !*/
915: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
916: implicit none
917: ! mx,my,bmx,bmy,dm,bheight defined in plate2f.h
918: #include "plate2f.h"
919: Tao tao
920: Vec xl,xu
921: PetscErrorCode ierr
922: PetscInt i,j,row
923: PetscInt xs, xm, ys, ym
924: PetscReal lb,ub
925: PetscInt dummy
927: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
928: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
929: ! will return an array of doubles referenced by x_array offset by x_index.
930: ! i.e., to reference the kth element of X, use x_array(k + x_index).
931: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
932: PetscReal xl_v(0:1)
933: PetscOffset xl_i
935: print *,'msa_plate'
937: lb = PETSC_NINFINITY
938: ub = PETSC_INFINITY
940: if (bmy .lt. 0) bmy = 0
941: if (bmy .gt. my) bmy = my
942: if (bmx .lt. 0) bmx = 0
943: if (bmx .gt. mx) bmx = mx
946: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
947: & PETSC_NULL_INTEGER,ierr)
949: call VecSet(xl,lb,ierr)
950: call VecSet(xu,ub,ierr)
952: call VecGetArray(xl,xl_v,xl_i,ierr)
955: do i=xs,xs+xm-1
957: do j=ys,ys+ym-1
959: row=(j-ys)*xm + (i-xs)
961: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
962: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
963: xl_v(xl_i+row) = bheight
965: endif
967: enddo
968: enddo
971: call VecRestoreArray(xl,xl_v,xl_i,ierr)
973: return
974: end
980: ! ----------------------------------------------------------------------------
981: !
982: !/*
983: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
984: !
985: ! Output Parameter:
986: !. X - vector for initial guess
987: !
988: !*/
990: subroutine MSA_InitialPoint(X, ierr)
991: implicit none
993: ! mx,my,localX,dm,Top,Left,Bottom,Right defined in plate2f.h
994: #include "plate2f.h"
995: Vec X
996: PetscErrorCode ierr
997: PetscInt start,i,j
998: PetscInt row
999: PetscInt xs,xm,gxs,gxm
1000: PetscInt ys,ym,gys,gym
1001: PetscReal zero, np5
1003: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1004: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
1005: ! will return an array of doubles referenced by x_array offset by x_index.
1006: ! i.e., to reference the kth element of X, use x_array(k + x_index).
1007: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1008: PetscReal left_v(0:1),right_v(0:1)
1009: PetscReal bottom_v(0:1),top_v(0:1)
1010: PetscReal x_v(0:1)
1011: PetscOffset left_i, right_i, top_i
1012: PetscOffset bottom_i,x_i
1013: PetscBool flg
1014: PetscRandom rctx
1016: zero = 0.0
1017: np5 = -0.5
1019: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
1020: & '-start', start,flg,ierr)
1022: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
1023: call VecSet(X,zero,ierr)
1025: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
1026: call PetscRandomCreate(MPI_COMM_WORLD,rctx,ierr)
1027: do i=0,start-1
1028: call VecSetRandom(X,rctx,ierr)
1029: enddo
1031: call PetscRandomDestroy(rctx,ierr)
1032: call VecShift(X,np5,ierr)
1034: else ! average of boundary conditions
1036: ! Get Local mesh boundaries
1037: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1038: & PETSC_NULL_INTEGER,ierr)
1039: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1040: & PETSC_NULL_INTEGER,ierr)
1044: ! Get pointers to vector data
1045: call VecGetArray(Top,top_v,top_i,ierr)
1046: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1047: call VecGetArray(Left,left_v,left_i,ierr)
1048: call VecGetArray(Right,right_v,right_i,ierr)
1050: call VecGetArray(localX,x_v,x_i,ierr)
1052: ! Perform local computations
1053: do j=ys,ys+ym-1
1054: do i=xs,xs+xm-1
1055: row = (j-gys)*gxm + (i-gxs)
1056: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1057: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1058: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1059: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1060: enddo
1061: enddo
1063: ! Restore vectors
1064: call VecRestoreArray(localX,x_v,x_i,ierr)
1066: call VecRestoreArray(Left,left_v,left_i,ierr)
1067: call VecRestoreArray(Top,top_v,top_i,ierr)
1068: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1069: call VecRestoreArray(Right,right_v,right_i,ierr)
1071: call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1072: call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)
1074: endif
1076: return
1077: end