Actual source code: plate2f.F
1: ! "$Id$";
2: !
3: ! Program usage: mpirun -np <proc> plate2f [all TAO options]
4: !
5: ! This example demonstrates use of the TAO package to solve a bound constrained
6: ! minimization problem. This example is based on a problem from the
7: ! MINPACK-2 test suite. Given a rectangular 2-D domain and boundary values
8: ! along the edges of the domain, the objective is to find the surface
9: ! with the minimal area that satisfies the boundary conditions.
10: ! The command line options are:
11: ! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
12: ! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
13: ! -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
14: ! -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
15: ! -bheight <ht>, where <ht> = height of the plate
16: !
17: !/*T
18: ! Concepts: TAO - Solving a bound constrained minimization problem
19: ! Routines: TaoInitialize(); TaoFinalize();
20: ! Routines: TaoCreate(); TaoDestroy();
21: ! Routines: TaoAppSetObjectiveAndGradientRoutine();
22: ! Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
23: ! Routines: TaoAppSetInitialSolutionVec(); TaoAppSetVariableBounds();
24: ! Routines: TaoSetOptions();
25: ! Routines: TaoApplicationCreate(); TaoSolve();
26: ! Routines: TaoView(); TaoAppDestroy();
27: ! Processors: n
28: !T*/
32: implicit none
34: #include "plate2f.h"
36: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37: ! Variable declarations
38: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39: !
40: ! Variables:
41: ! (common from plate2f.h):
42: ! Nx, Ny number of processors in x- and y- directions
43: ! mx, my number of grid points in x,y directions
44: ! N global dimension of vector
46: integer info ! used to check for functions returning nonzeros
47: Vec x ! solution vector
48: Vec xl, xu ! lower and upper bounds vectorsp
49: PetscInt m ! number of local elements in vector
50: TAO_SOLVER tao ! TAO_SOLVER solver context
51: TAO_APPLICATION plateapp ! PETSc application
52: Mat H ! Hessian matrix
53: ISLocalToGlobalMapping isltog ! local to global mapping object
54: PetscTruth flg
55: PetscInt i1,i3,i7
58: external FormFunctionGradient
59: external FormHessian
60: external MSA_BoundaryConditions
61: external MSA_Plate
62: external MSA_InitialPoint
63: ! Initialize Tao
65: i1=1
66: i3=3
67: i7=7
70: call PetscInitialize(PETSC_NULL_CHARACTER,info)
71: call TaoInitialize(PETSC_NULL_CHARACTER,info)
72:
74: ! Specify default dimensions of the problem
75: mx = 10
76: my = 10
77: bheight = 0.1
79: ! Check for any command line arguments that override defaults
80:
81: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-mx",mx,flg,info)
82: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-my",my,flg,info)
83:
84: bmx = mx/2
85: bmy = my/2
87: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmx",bmx,flg,info)
88: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmy",bmy,flg,info)
89: call PetscOptionsGetReal(TAO_NULL_CHARACTER,"-bheight",bheight, &
90: & flg,info)
91:
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 DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX, &
105: & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
106: & da,info)
107:
109: ! Extract global and local vectors from DA; The local vectors are
110: ! used solely as work space for the evaluation of the function,
111: ! gradient, and Hessian. Duplicate for remaining vectors that are
112: ! the same types.
114: call DACreateGlobalVector(da,x,info)
115: call DACreateLocalVector(da,localX,info)
116: call VecDuplicate(localX,localV,info)
118: ! Create a matrix data structure to store the Hessian.
119: ! Here we (optionally) also associate the local numbering scheme
120: ! with the matrix so that later we can use local indices for matrix
121: ! assembly
123: call VecGetLocalSize(x,m,info)
124: call MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, &
125: & i3,PETSC_NULL_INTEGER,H,info)
127: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)
128: call DAGetISLocalToGlobalMapping(da,isltog,info)
129: call MatSetLocalToGlobalMapping(H,isltog,info)
130:
132: ! The Tao code begins here
133: ! Create TAO solver and set desired solution method.
134: ! This problems uses bounded variables, so the
135: ! method must either be 'tao_tron' or 'tao_blmvm'
137: call TaoCreate(MPI_COMM_WORLD,'tao_blmvm',tao,info)
138: call TaoApplicationCreate(MPI_COMM_WORLD,plateapp,info)
139:
141: ! Set minimization function and gradient, hessian evaluation functions
143: ! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
144: call TaoAppSetObjectiveAndGradientRo(plateapp, &
145: & FormFunctionGradient,PETSC_NULL_OBJECT,info)
147: call TaoAppSetHessianMat(plateapp,H,H,info)
148: call TaoAppSetHessianRoutine(plateapp,FormHessian, &
149: & PETSC_NULL_OBJECT, info)
150:
152: ! Set Variable bounds
153: call MSA_BoundaryConditions(info)
154: call VecDuplicate(x,xl,info)
155: call VecDuplicate(x,xu,info)
156: call MSA_Plate(xl,xu,info)
157: call TaoAppSetVariableBounds(plateapp,xl,xu,info)
159: ! Set the initial solution guess
160: call MSA_InitialPoint(x, info)
161: call TaoAppSetInitialSolutionVec(plateapp,x,info)
163: ! Check for any tao command line options
164: call TaoSetOptions(plateapp,tao,info)
166: ! Solve the application
167: call TaoSolveApplication(plateapp,tao,info)
169: ! View TAO solver information
170: ! call TaoView(tao,info)
172: ! Free TAO data structures
173: call TaoDestroy(tao,info)
174: call TaoAppDestroy(plateapp,info)
176: ! Free PETSc data structures
177: call VecDestroy(x,info)
178: call VecDestroy(xl,info)
179: call VecDestroy(xu,info)
180: call VecDestroy(Top,info)
181: call VecDestroy(Bottom,info)
182: call VecDestroy(Left,info)
183: call VecDestroy(Right,info)
184: call MatDestroy(H,info)
185: call VecDestroy(localX,info)
186: call VecDestroy(localV,info)
187: call DADestroy(da,info)
189: ! Finalize TAO
191: call TaoFinalize(info)
192: call PetscFinalize(info)
194: end
196: ! ---------------------------------------------------------------------
197: !
198: ! FormFunctionGradient - Evaluates function f(X).
199: !
200: ! Input Parameters:
201: ! tao - the TAO_SOLVER context
202: ! X - the input vector
203: ! dummy - optional user-defined context, as set by TaoSetFunction()
204: ! (not used here)
205: !
206: ! Output Parameters:
207: ! fcn - the newly evaluated function
208: ! G - the gradient vector
209: ! info - error code
210: !
213: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,info)
214: implicit none
216: ! da, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
217: #include "plate2f.h"
218:
219: ! Input/output variables
221: TAO_SOLVER tao
222: PetscScalar fcn
223: Vec X, G
224: integer info
225: PetscInt dummy
226:
227: PetscInt i,j,row
228: PetscInt xs, xm
229: PetscInt gxs, gxm
230: PetscInt ys, ym
231: PetscInt gys, gym
232: PetscScalar ft,zero,hx,hy,hydhx,hxdhy
233: PetscScalar area,rhx,rhy
234: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3
235: PetscScalar d4,d5,d6,d7,d8
236: PetscScalar df1dxc,df2dxc,df3dxc,df4dxc
237: PetscScalar df5dxc,df6dxc
238: PetscScalar xc,xl,xr,xt,xb,xlt,xrb
241: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
242: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
243: ! will return an array of doubles referenced by x_array offset by x_index.
244: ! i.e., to reference the kth element of X, use x_array(k + x_index).
245: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
246: PetscScalar g_v(0:1),x_v(0:1)
247: PetscScalar top_v(0:1),left_v(0:1)
248: PetscScalar right_v(0:1),bottom_v(0:1)
249: PetscOffset g_i,left_i,right_i
250: PetscOffset bottom_i,top_i,x_i
252: ft = 0.0d0
253: zero = 0.0d0
254: hx = 1.0d0/(mx + 1)
255: hy = 1.0d0/(my + 1)
256: hydhx = hy/hx
257: hxdhy = hx/hy
258: area = 0.5d0 * hx * hy
259: rhx = mx + 1.0d0
260: rhy = my + 1.0d0
263: ! Get local mesh boundaries
264: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
265: & PETSC_NULL_INTEGER,info)
266: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
267: & gxm,gym,PETSC_NULL_INTEGER,info)
269: ! Scatter ghost points to local vector
270: call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
271: call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
273: ! Initialize the vector to zero
274: call VecSet(localV,zero,info)
276: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
277: call VecGetArray(localX,x_v,x_i,info)
278: call VecGetArray(localV,g_v,g_i,info)
279: call VecGetArray(Top,top_v,top_i,info)
280: call VecGetArray(Bottom,bottom_v,bottom_i,info)
281: call VecGetArray(Left,left_v,left_i,info)
282: call VecGetArray(Right,right_v,right_i,info)
284: ! Compute function over the locally owned part of the mesh
285: do j = ys,ys+ym-1
286: do i = xs,xs+xm-1
287: row = (j-gys)*gxm + (i-gxs)
288: xc = x_v(row+x_i)
289: xt = xc
290: xb = xc
291: xr = xc
292: xl = xc
293: xrb = xc
294: xlt = xc
296: if (i .eq. 0) then !left side
297: xl = left_v(j - ys + 1 + left_i)
298: xlt = left_v(j - ys + 2 + left_i)
299: else
300: xl = x_v(row - 1 + x_i)
301: endif
303: if (j .eq. 0) then !bottom side
304: xb = bottom_v(i - xs + 1 + bottom_i)
305: xrb = bottom_v(i - xs + 2 + bottom_i)
306: else
307: xb = x_v(row - gxm + x_i)
308: endif
310: if (i + 1 .eq. gxs + gxm) then !right side
311: xr = right_v(j - ys + 1 + right_i)
312: xrb = right_v(j - ys + right_i)
313: else
314: xr = x_v(row + 1 + x_i)
315: endif
317: if (j + 1 .eq. gys + gym) then !top side
318: xt = top_v(i - xs + 1 + top_i)
319: xlt = top_v(i - xs + top_i)
320: else
321: xt = x_v(row + gxm + x_i)
322: endif
324: if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
325: xlt = x_v(row - 1 + gxm + x_i)
326: endif
328: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
329: xrb = x_v(row + 1 - gxm + x_i)
330: endif
332: d1 = xc-xl
333: d2 = xc-xr
334: d3 = xc-xt
335: d4 = xc-xb
336: d5 = xr-xrb
337: d6 = xrb-xb
338: d7 = xlt-xl
339: d8 = xt-xlt
341: df1dxc = d1 * hydhx
342: df2dxc = d1 * hydhx + d4 * hxdhy
343: df3dxc = d3 * hxdhy
344: df4dxc = d2 * hydhx + d3 * hxdhy
345: df5dxc = d2 * hydhx
346: df6dxc = d4 * hxdhy
348: d1 = d1 * rhx
349: d2 = d2 * rhx
350: d3 = d3 * rhy
351: d4 = d4 * rhy
352: d5 = d5 * rhy
353: d6 = d6 * rhx
354: d7 = d7 * rhy
355: d8 = d8 * rhx
357: f1 = sqrt(1.0d0 + d1*d1 + d7*d7)
358: f2 = sqrt(1.0d0 + d1*d1 + d4*d4)
359: f3 = sqrt(1.0d0 + d3*d3 + d8*d8)
360: f4 = sqrt(1.0d0 + d3*d3 + d2*d2)
361: f5 = sqrt(1.0d0 + d2*d2 + d5*d5)
362: f6 = sqrt(1.0d0 + d4*d4 + d6*d6)
364: ft = ft + f2 + f4
366: df1dxc = df1dxc / f1
367: df2dxc = df2dxc / f2
368: df3dxc = df3dxc / f3
369: df4dxc = df4dxc / f4
370: df5dxc = df5dxc / f5
371: df6dxc = df6dxc / f6
373: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
374: & df5dxc + df6dxc)
375: enddo
376: enddo
378: ! Compute triangular areas along the border of the domain.
379: if (xs .eq. 0) then ! left side
380: do j=ys,ys+ym-1
381: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
382: & * rhy
383: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
384: & * rhx
385: ft = ft + sqrt(1.0d0 + d3*d3 + d2*d2)
386: enddo
387: endif
389:
390: if (ys .eq. 0) then !bottom side
391: do i=xs,xs+xm-1
392: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
393: & * rhx
394: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
395: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
396: enddo
397: endif
399:
400: if (xs + xm .eq. mx) then ! right side
401: do j=ys,ys+ym-1
402: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
403: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
404: ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
405: enddo
406: endif
408:
409: if (ys + ym .eq. my) then
410: do i=xs,xs+xm-1
411: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
412: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
413: ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
414: enddo
415: endif
417:
418: if ((ys .eq. 0) .and. (xs .eq. 0)) then
419: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
420: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
421: ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
422: endif
424: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
425: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
426: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
427: ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
428: endif
430: ft = ft * area
431: call MPI_Allreduce(ft,fcn,1,MPI_DOUBLE_PRECISION, &
432: & MPI_SUM,MPI_COMM_WORLD,info)
436: ! Restore vectors
437: call VecRestoreArray(localX,x_v,x_i,info)
438: call VecRestoreArray(localV,g_v,g_i,info)
439: call VecRestoreArray(Left,left_v,left_i,info)
440: call VecRestoreArray(Top,top_v,top_i,info)
441: call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
442: call VecRestoreArray(Right,right_v,right_i,info)
444: ! Scatter values to global vector
445: call DALocalToGlobal(da,localV,INSERT_VALUES,G,info)
447: call PetscLogFlops(70*xm*ym,info)
449: return
450: end !FormFunctionGradient
451:
456: ! ----------------------------------------------------------------------------
457: !
458: !/*
459: ! FormHessian - Evaluates Hessian matrix.
460: !
461: ! Input Parameters:
462: !. tao - the TAO_SOLVER context
463: !. X - input vector
464: !. dummy - not used
465: !
466: ! Output Parameters:
467: !. Hessian - Hessian matrix
468: !. Hpc - optionally different preconditioning matrix
469: !. flag - flag indicating matrix structure
470: !
471: ! Notes:
472: ! Due to mesh point reordering with DAs, we must always work
473: ! with the local mesh points, and then transform them to the new
474: ! global numbering with the local-to-global mapping. We cannot work
475: ! directly with the global numbers for the original uniprocessor mesh!
476: !
477: ! Two methods are available for imposing this transformation
478: ! when setting matrix entries:
479: ! (A) MatSetValuesLocal(), using the local ordering (including
480: ! ghost points!)
481: ! - Do the following two steps once, before calling TaoSolve()
482: ! - Use DAGetISLocalToGlobalMapping() to extract the
483: ! local-to-global map from the DA
484: ! - Associate this map with the matrix by calling
485: ! MatSetLocalToGlobalMapping()
486: ! - Then set matrix entries using the local ordering
487: ! by calling MatSetValuesLocal()
488: ! (B) MatSetValues(), using the global ordering
489: ! - Use DAGetGlobalIndices() to extract the local-to-global map
490: ! - Then apply this map explicitly yourself
491: ! - Set matrix entries using the global ordering by calling
492: ! MatSetValues()
493: ! Option (A) seems cleaner/easier in many cases, and is the procedure
494: ! used in this example.
495: */
496: subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, info)
497: implicit none
499: ! da,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
500: #include "plate2f.h"
501:
502: TAO_SOLVER tao
503: Vec X
504: Mat Hessian,Hpc
505: MatStructure flg
506: PetscInt dummy
507: integer info
509: PetscInt i,j,k,row
510: PetscInt xs,xm,gxs,gxm
511: PetscInt ys,ym,gys,gym
512: PetscInt col(0:6)
513: PetscScalar hx,hy,hydhx,hxdhy,rhx,rhy
514: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3
515: PetscScalar d4,d5,d6,d7,d8
516: PetscScalar xc,xl,xr,xt,xb,xlt,xrb
517: PetscScalar hl,hr,ht,hb,hc,htl,hbr
519: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
520: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
521: ! will return an array of doubles referenced by x_array offset by x_index.
522: ! i.e., to reference the kth element of X, use x_array(k + x_index).
523: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
524: PetscScalar right_v(0:1),left_v(0:1)
525: PetscScalar bottom_v(0:1),top_v(0:1)
526: PetscScalar x_v(0:1)
527: PetscOffset x_i,right_i,left_i
528: PetscOffset bottom_i,top_i
529: PetscScalar v(0:6)
530: PetscTruth assembled
531: PetscInt i1
532:
533: i1=1
534:
535: ! Set various matrix options
536: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,
537: & PETSC_TRUE,info)
539: ! Get local mesh boundaries
540: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
541: & PETSC_NULL_INTEGER,info)
542: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
543: & PETSC_NULL_INTEGER,info)
545: ! Scatter ghost points to local vectors
546: call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
547: call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
549: ! Get pointers to vector data (see note on Fortran arrays above)
550: call VecGetArray(localX,x_v,x_i,info)
551: call VecGetArray(Top,top_v,top_i,info)
552: call VecGetArray(Bottom,bottom_v,bottom_i,info)
553: call VecGetArray(Left,left_v,left_i,info)
554: call VecGetArray(Right,right_v,right_i,info)
556: ! Initialize matrix entries to zero
557: call MatAssembled(Hessian,assembled,info)
558: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,info)
561: rhx = mx + 1.0
562: rhy = my + 1.0
563: hx = 1.0/rhx
564: hy = 1.0/rhy
565: hydhx = hy/hx
566: hxdhy = hx/hy
567: ! compute Hessian over the locally owned part of the mesh
569: do i=xs,xs+xm-1
570: do j=ys,ys+ym-1
571: row = (j-gys)*gxm + (i-gxs)
572:
573: xc = x_v(row + x_i)
574: xt = xc
575: xb = xc
576: xr = xc
577: xl = xc
578: xrb = xc
579: xlt = xc
581: if (i .eq. gxs) then ! Left side
582: xl = left_v(left_i + j - ys + 1)
583: xlt = left_v(left_i + j - ys + 2)
584: else
585: xl = x_v(x_i + row -1 )
586: endif
588: if (j .eq. gys) then ! bottom side
589: xb = bottom_v(bottom_i + i - xs + 1)
590: xrb = bottom_v(bottom_i + i - xs + 2)
591: else
592: xb = x_v(x_i + row - gxm)
593: endif
595: if (i+1 .eq. gxs + gxm) then !right side
596: xr = right_v(right_i + j - ys + 1)
597: xrb = right_v(right_i + j - ys)
598: else
599: xr = x_v(x_i + row + 1)
600: endif
602: if (j+1 .eq. gym+gys) then !top side
603: xt = top_v(top_i +i - xs + 1)
604: xlt = top_v(top_i + i - xs)
605: else
606: xt = x_v(x_i + row + gxm)
607: endif
609: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
610: xlt = x_v(x_i + row - 1 + gxm)
611: endif
612:
613: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
614: xrb = x_v(x_i + row + 1 - gxm)
615: endif
617: d1 = (xc-xl)*rhx
618: d2 = (xc-xr)*rhx
619: d3 = (xc-xt)*rhy
620: d4 = (xc-xb)*rhy
621: d5 = (xrb-xr)*rhy
622: d6 = (xrb-xb)*rhx
623: d7 = (xlt-xl)*rhy
624: d8 = (xlt-xt)*rhx
625:
626: f1 = sqrt( 1.0d0 + d1*d1 + d7*d7)
627: f2 = sqrt( 1.0d0 + d1*d1 + d4*d4)
628: f3 = sqrt( 1.0d0 + d3*d3 + d8*d8)
629: f4 = sqrt( 1.0d0 + d3*d3 + d2*d2)
630: f5 = sqrt( 1.0d0 + d2*d2 + d5*d5)
631: f6 = sqrt( 1.0d0 + d4*d4 + d6*d6)
632:
633:
634: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
635: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
637: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
638: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
640: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
641: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
643: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
644: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
645:
646: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
647: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
648:
649: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
650: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
651: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
652: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
653: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
654: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
655: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
656:
657: hl = hl * 0.5
658: hr = hr * 0.5
659: ht = ht * 0.5
660: hb = hb * 0.5
661: hbr = hbr * 0.5
662: htl = htl * 0.5
663: hc = hc * 0.5
665: k = 0
667: if (j .gt. 0) then
668: v(k) = hb
669: col(k) = row - gxm
670: k=k+1
671: endif
673: if ((j .gt. 0) .and. (i .lt. mx-1)) then
674: v(k) = hbr
675: col(k) = row-gxm+1
676: k=k+1
677: endif
679: if (i .gt. 0) then
680: v(k) = hl
681: col(k) = row - 1
682: k = k+1
683: endif
685: v(k) = hc
686: col(k) = row
687: k=k+1
689: if (i .lt. mx-1) then
690: v(k) = hr
691: col(k) = row + 1
692: k=k+1
693: endif
695: if ((i .gt. 0) .and. (j .lt. my-1)) then
696: v(k) = htl
697: col(k) = row + gxm - 1
698: k=k+1
699: endif
701: if (j .lt. my-1) then
702: v(k) = ht
703: col(k) = row + gxm
704: k=k+1
705: endif
707: ! Set matrix values using local numbering, defined earlier in main routine
708: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
709: & info)
711:
713: enddo
714: enddo
715:
716: ! restore vectors
717: call VecRestoreArray(localX,x_v,x_i,info)
718: call VecRestoreArray(Left,left_v,left_i,info)
719: call VecRestoreArray(Right,right_v,right_i,info)
720: call VecRestoreArray(Top,top_v,top_i,info)
721: call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
724: ! Assemble the matrix
725: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,info)
726: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,info)
728: call PetscLogFlops(199*xm*ym,info)
730: return
731: end
732:
733:
737: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
739: ! ----------------------------------------------------------------------------
740: !
741: !/*
742: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
743: !
744: !
745: !*/
747: subroutine MSA_BoundaryConditions(info)
748: implicit none
750: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
751: #include "plate2f.h"
753: integer info
754: PetscInt i,j,k,limit,maxits
755: PetscInt xs, xm, gxs, gxm
756: PetscInt ys, ym, gys, gym
757: PetscInt bsize, lsize
758: PetscInt tsize, rsize
759: PetscScalar one,two,three,tol
760: PetscScalar scl,fnorm,det,xt
761: PetscScalar yt,hx,hy,u1,u2,nf1,nf2
762: PetscScalar njac11,njac12,njac21,njac22
763: PetscScalar b, t, l, r
764: PetscScalar boundary_v(0:1)
765: PetscOffset boundary_i
766: logical exitloop
767: TaoTruth flg
769: limit=0
770: maxits = 5
771: tol=1e-10
772: b=-0.5d0
773: t= 0.5d0
774: l=-0.5d0
775: r= 0.5d0
776: xt=0
777: yt=0
778: one=1.0d0
779: two=2.0d0
780: three=3.0d0
783: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
784: & PETSC_NULL_INTEGER,info)
785: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
786: & gxm,gym,PETSC_NULL_INTEGER,info)
788: bsize = xm + 2
789: lsize = ym + 2
790: rsize = ym + 2
791: tsize = xm + 2
792:
794: call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,info)
795: call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,info)
796: call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,info)
797: call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,info)
799: hx= (r-l)/(mx+1)
800: hy= (t-b)/(my+1)
802: do j=0,3
803:
804: if (j.eq.0) then
805: yt=b
806: xt=l+hx*xs
807: limit=bsize
808: call VecGetArray(Bottom,boundary_v,boundary_i,info)
809:
811: elseif (j.eq.1) then
812: yt=t
813: xt=l+hx*xs
814: limit=tsize
815: call VecGetArray(Top,boundary_v,boundary_i,info)
817: elseif (j.eq.2) then
818: yt=b+hy*ys
819: xt=l
820: limit=lsize
821: call VecGetArray(Left,boundary_v,boundary_i,info)
823: elseif (j.eq.3) then
824: yt=b+hy*ys
825: xt=r
826: limit=rsize
827: call VecGetArray(Right,boundary_v,boundary_i,info)
828: endif
829:
831: do i=0,limit-1
832:
833: u1=xt
834: u2=-yt
835: k = 0
836: exitloop = .false.
837: do while (k .lt. maxits .and. (.not. exitloop) )
839: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
840: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
841: fnorm=sqrt(nf1*nf1+nf2*nf2)
842: if (fnorm .gt. tol) then
843: njac11=one+u2*u2-u1*u1
844: njac12=two*u1*u2
845: njac21=-two*u1*u2
846: njac22=-one - u1*u1 + u2*u2
847: det = njac11*njac22-njac21*njac12
848: u1 = u1-(njac22*nf1-njac12*nf2)/det
849: u2 = u2-(njac11*nf2-njac21*nf1)/det
850: else
851: exitloop = .true.
852: endif
853: k=k+1
854: enddo
856: boundary_v(i + boundary_i) = u1*u1-u2*u2
857: if ((j .eq. 0) .or. (j .eq. 1)) then
858: xt = xt + hx
859: else
860: yt = yt + hy
861: endif
863: enddo
864:
866: if (j.eq.0) then
867: call VecRestoreArray(Bottom,boundary_v,boundary_i,info)
868: elseif (j.eq.1) then
869: call VecRestoreArray(Top,boundary_v,boundary_i,info)
870: elseif (j.eq.2) then
871: call VecRestoreArray(Left,boundary_v,boundary_i,info)
872: elseif (j.eq.3) then
873: call VecRestoreArray(Right,boundary_v,boundary_i,info)
874: endif
875:
876: enddo
879: ! Scale the boundary if desired
880: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom", &
881: & scl,flg,info)
882: if (flg .eqv. PETSC_TRUE) then
883: call VecScale(scl,Bottom,info)
884: endif
886: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top", &
887: & scl,flg,info)
888: if (flg .eqv. PETSC_TRUE) then
889: call VecScale(scl,Top,info)
890: endif
892: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right", &
893: & scl,flg,info)
894: if (flg .eqv. PETSC_TRUE) then
895: call VecScale(scl,Right,info)
896: endif
898: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left", &
899: & scl,flg,info)
900: if (flg .eqv. PETSC_TRUE) then
901: call VecScale(scl,Left,info)
902: endif
903:
904:
905: return
906: end
908: ! ----------------------------------------------------------------------------
909: !
910: !/*
911: ! MSA_Plate - Calculates an obstacle for surface to stretch over
912: !
913: ! Output Parameter:
914: !. xl - lower bound vector
915: !. xu - upper bound vector
916: !
917: !*/
919: subroutine MSA_Plate(xl,xu,info)
920: implicit none
922: ! mx,my,bmx,bmy,da,bheight defined in plate2f.h
923: #include "plate2f.h"
924: Vec xl,xu
925: integer info
926: PetscInt i,j,row
927: PetscInt xs, xm, ys, ym
928: PetscScalar lb,ub
930: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
931: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
932: ! will return an array of doubles referenced by x_array offset by x_index.
933: ! i.e., to reference the kth element of X, use x_array(k + x_index).
934: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
935: PetscScalar xl_v(0:1)
936: PetscOffset xl_i
939: lb = -1.0d300
940: ub = 1.0d300
942: if (bmy .lt. 0) bmy = 0
943: if (bmy .gt. my) bmy = my
944: if (bmx .lt. 0) bmx = 0
945: if (bmx .gt. mx) bmx = mx
946:
948: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
949: & PETSC_NULL_INTEGER,info)
951: call VecSet(xl,lb,info)
952: call VecSet(xu,ub,info)
954: call VecGetArray(xl,xl_v,xl_i,info)
955:
957: do i=xs,xs+xm-1
959: do j=ys,ys+ym-1
960:
961: row=(j-ys)*xm + (i-xs)
963: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
964: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
965: xl_v(xl_i+row) = bheight
967: endif
969: enddo
970: enddo
973: call VecRestoreArray(xl,xl_v,xl_i,info)
974:
975: return
976: end
980:
981:
982: ! ----------------------------------------------------------------------------
983: !
984: !/*
985: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
986: !
987: ! Output Parameter:
988: !. X - vector for initial guess
989: !
990: !*/
992: subroutine MSA_InitialPoint(X, info)
993: implicit none
995: ! mx,my,localX,da,Top,Left,Bottom,Right defined in plate2f.h
996: #include "plate2f.h"
997: Vec X
998: integer info
999: PetscInt start,i,j
1000: PetscInt row
1001: PetscInt xs,xm,gxs,gxm
1002: PetscInt ys,ym,gys,gym
1003: PetscScalar zero, np5
1005: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1006: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (integer) x_index, info)
1007: ! will return an array of doubles referenced by x_array offset by x_index.
1008: ! i.e., to reference the kth element of X, use x_array(k + x_index).
1009: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1010: PetscScalar left_v(0:1),right_v(0:1)
1011: PetscScalar bottom_v(0:1),top_v(0:1)
1012: PetscScalar x_v(0:1)
1013: PetscOffset left_i, right_i, top_i
1014: PetscOffset bottom_i,x_i
1015: PetscTruth flg
1016: PetscRandom rctx
1017:
1018: zero = 0.0d0
1019: np5 = -0.5d0
1021: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start", &
1022: & start,flg,info)
1024: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
1025: call VecSet(X,zero,info)
1027: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
1028: call PetscRandomCreate(MPI_COMM_WORLD,rctx,info)
1029: do i=0,start-1
1030: call VecSetRandom(X,rctx,info)
1031: enddo
1033: call PetscRandomDestroy(rctx,info)
1034: call VecShift(X,np5,info)
1036: else ! average of boundary conditions
1037:
1038: ! Get Local mesh boundaries
1039: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1040: & PETSC_NULL_INTEGER,info)
1041: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1042: & PETSC_NULL_INTEGER,info)
1046: ! Get pointers to vector data
1047: call VecGetArray(Top,top_v,top_i,info)
1048: call VecGetArray(Bottom,bottom_v,bottom_i,info)
1049: call VecGetArray(Left,left_v,left_i,info)
1050: call VecGetArray(Right,right_v,right_i,info)
1051:
1052: call VecGetArray(localX,x_v,x_i,info)
1053:
1054: ! Perform local computations
1055: do j=ys,ys+ym-1
1056: do i=xs,xs+xm-1
1057: row = (j-gys)*gxm + (i-gxs)
1058: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1059: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1060: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1061: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1062: enddo
1063: enddo
1065: ! Restore vectors
1066: call VecRestoreArray(localX,x_v,x_i,info)
1068: call VecRestoreArray(Left,left_v,left_i,info)
1069: call VecRestoreArray(Top,top_v,top_i,info)
1070: call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
1071: call VecRestoreArray(Right,right_v,right_i,info)
1073: call DALocalToGlobal(da,localX,INSERT_VALUES,X,info)
1075: endif
1077: return
1078: end