Actual source code: plate2f.F
1: ! "$Id: plate2f.F 1.47 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
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: TaoSetApplication(); 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: integer 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
57: external FormFunctionGradient
58: external FormHessian
59: external MSA_BoundaryConditions
60: external MSA_Plate
61: external MSA_InitialPoint
62: ! Initialize Tao
63: call PetscInitialize(PETSC_NULL_CHARACTER,info)
64: call TaoInitialize(PETSC_NULL_CHARACTER,info)
65:
67: ! Specify default dimensions of the problem
68: mx = 10
69: my = 10
70: bheight = 0.1
72: ! Check for any command line arguments that override defaults
73:
74: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-mx",mx,flg,info)
75: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-my",my,flg,info)
76:
77: bmx = mx/2
78: bmy = my/2
80: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmx",bmx,flg,info)
81: call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmy",bmy,flg,info)
82: call PetscOptionsGetReal(TAO_NULL_CHARACTER,"-bheight",bheight, &
83: & flg,info)
84:
86: ! Calculate any derived values from parameters
87: N = mx*my
89: ! Let Petsc determine the dimensions of the local vectors
90: Nx = PETSC_DECIDE
91: NY = PETSC_DECIDE
93: ! A two dimensional distributed array will help define this problem, which
94: ! derives from an elliptic PDE on a two-dimensional domain. From the
95: ! distributed array, create the vectors
97: call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX, &
98: & mx,my,Nx,Ny,1,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
99: & da,info)
100:
102: ! Extract global and local vectors from DA; The local vectors are
103: ! used solely as work space for the evaluation of the function,
104: ! gradient, and Hessian. Duplicate for remaining vectors that are
105: ! the same types.
107: call DACreateGlobalVector(da,x,info)
108: call DACreateLocalVector(da,localX,info)
109: call VecDuplicate(localX,localV,info)
111: ! Create a matrix data structure to store the Hessian.
112: ! Here we (optionally) also associate the local numbering scheme
113: ! with the matrix so that later we can use local indices for matrix
114: ! assembly
116: call VecGetLocalSize(x,m,info)
117: call MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL_INTEGER, &
118: & 3,PETSC_NULL_INTEGER,H,info)
120: call MatSetOption(H,MAT_SYMMETRIC,info)
121: call DAGetISLocalToGlobalMapping(da,isltog,info)
122: call MatSetLocalToGlobalMapping(H,isltog,info)
123:
125: ! The Tao code begins here
126: ! Create TAO solver and set desired solution method.
127: ! This problems uses bounded variables, so the
128: ! method must either be 'tao_tron' or 'tao_blmvm'
130: call TaoCreate(MPI_COMM_WORLD,'tao_blmvm',tao,info)
131: call TaoApplicationCreate(MPI_COMM_WORLD,plateapp,info)
132:
134: ! Set minimization function and gradient, hessian evaluation functions
136: ! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
137: call TaoAppSetObjectiveAndGradientRo(plateapp, &
138: & FormFunctionGradient,PETSC_NULL_OBJECT,info)
140: call TaoAppSetHessianMat(plateapp,H,H,info)
141: call TaoAppSetHessianRoutine(plateapp,FormHessian, &
142: & PETSC_NULL_OBJECT, info)
143:
145: ! Set Variable bounds
146: call MSA_BoundaryConditions(info)
147: call VecDuplicate(x,xl,info)
148: call VecDuplicate(x,xu,info)
149: call MSA_Plate(xl,xu,info)
150: call TaoAppSetVariableBounds(plateapp,xl,xu,info)
152: ! Set the initial solution guess
153: call MSA_InitialPoint(x, info)
154: call TaoAppSetInitialSolutionVec(plateapp,x,info)
156: ! Now that the PETSc application is all set, attach to tao context
157: call TaoSetApplication(tao,plateapp,info)
159: ! Check for any tao command line options
160: call TaoSetOptions(plateapp,tao,info)
162: ! Solve the application
163: call TaoSolveApplication(plateapp,tao,info)
165: ! View TAO solver information
166: ! call TaoView(tao,info)
168: ! Free TAO data structures
169: call TaoDestroy(tao,info)
170: call TaoAppDestroy(plateapp,info)
172: ! Free PETSc data structures
173: call VecDestroy(x,info)
174: call VecDestroy(xl,info)
175: call VecDestroy(xu,info)
176: call VecDestroy(Top,info)
177: call VecDestroy(Bottom,info)
178: call VecDestroy(Left,info)
179: call VecDestroy(Right,info)
180: call MatDestroy(H,info)
181: call VecDestroy(localX,info)
182: call VecDestroy(localV,info)
183: call DADestroy(da,info)
185: ! Finalize TAO
187: call TaoFinalize(info)
188: call PetscFinalize(info)
190: end
192: ! ---------------------------------------------------------------------
193: !
194: ! FormFunctionGradient - Evaluates function f(X).
195: !
196: ! Input Parameters:
197: ! tao - the TAO_SOLVER context
198: ! X - the input vector
199: ! dummy - optional user-defined context, as set by TaoSetFunction()
200: ! (not used here)
201: !
202: ! Output Parameters:
203: ! fcn - the newly evaluated function
204: ! G - the gradient vector
205: ! info - error code
206: !
209: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,info)
210: implicit none
212: ! da, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
213: #include "plate2f.h"
214:
215: ! Input/output variables
217: TAO_SOLVER tao
218: PetscScalar fcn
219: Vec X, G
220: integer dummy, info
221:
222: integer i,j,row
223: integer xs, xm, gxs, gxm, ys, ym, gys, gym
224: PetscScalar ft,zero,hx,hy,hydhx,hxdhy,area,rhx,rhy
225: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8
226: PetscScalar df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc
227: PetscScalar xc,xl,xr,xt,xb,xlt,xrb
230: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
231: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
232: ! will return an array of doubles referenced by x_array offset by x_index.
233: ! i.e., to reference the kth element of X, use x_array(k + x_index).
234: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
235: PetscScalar g_v(0:1),x_v(0:1), top_v(0:1)
236: PetscScalar left_v(0:1), right_v(0:1),bottom_v(0:1)
237: PetscOffset g_i,left_i,right_i,bottom_i,top_i
238: PetscOffset x_i
240: ft = 0.0d0
241: zero = 0.0d0
242: hx = 1.0d0/(mx + 1)
243: hy = 1.0d0/(my + 1)
244: hydhx = hy/hx
245: hxdhy = hx/hy
246: area = 0.5d0 * hx * hy
247: rhx = mx + 1.0d0
248: rhy = my + 1.0d0
251: ! Get local mesh boundaries
252: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
253: & PETSC_NULL_INTEGER,info)
254: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
255: & gxm,gym,PETSC_NULL_INTEGER,info)
257: ! Scatter ghost points to local vector
258: call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
259: call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
261: ! Initialize the vector to zero
262: call VecSet(localV,zero,info)
264: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
265: call VecGetArray(localX,x_v,x_i,info)
266: call VecGetArray(localV,g_v,g_i,info)
267: call VecGetArray(Top,top_v,top_i,info)
268: call VecGetArray(Bottom,bottom_v,bottom_i,info)
269: call VecGetArray(Left,left_v,left_i,info)
270: call VecGetArray(Right,right_v,right_i,info)
272: ! Compute function over the locally owned part of the mesh
273: do j = ys,ys+ym-1
274: do i = xs,xs+xm-1
275: row = (j-gys)*gxm + (i-gxs)
276: xc = x_v(row+x_i)
277: xt = xc
278: xb = xc
279: xr = xc
280: xl = xc
281: xrb = xc
282: xlt = xc
284: if (i .eq. 0) then !left side
285: xl = left_v(j - ys + 1 + left_i)
286: xlt = left_v(j - ys + 2 + left_i)
287: else
288: xl = x_v(row - 1 + x_i)
289: endif
291: if (j .eq. 0) then !bottom side
292: xb = bottom_v(i - xs + 1 + bottom_i)
293: xrb = bottom_v(i - xs + 2 + bottom_i)
294: else
295: xb = x_v(row - gxm + x_i)
296: endif
298: if (i + 1 .eq. gxs + gxm) then !right side
299: xr = right_v(j - ys + 1 + right_i)
300: xrb = right_v(j - ys + right_i)
301: else
302: xr = x_v(row + 1 + x_i)
303: endif
305: if (j + 1 .eq. gys + gym) then !top side
306: xt = top_v(i - xs + 1 + top_i)
307: xlt = top_v(i - xs + top_i)
308: else
309: xt = x_v(row + gxm + x_i)
310: endif
312: if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
313: xlt = x_v(row - 1 + gxm + x_i)
314: endif
316: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
317: xrb = x_v(row + 1 - gxm + x_i)
318: endif
320: d1 = xc-xl
321: d2 = xc-xr
322: d3 = xc-xt
323: d4 = xc-xb
324: d5 = xr-xrb
325: d6 = xrb-xb
326: d7 = xlt-xl
327: d8 = xt-xlt
329: df1dxc = d1 * hydhx
330: df2dxc = d1 * hydhx + d4 * hxdhy
331: df3dxc = d3 * hxdhy
332: df4dxc = d2 * hydhx + d3 * hxdhy
333: df5dxc = d2 * hydhx
334: df6dxc = d4 * hxdhy
336: d1 = d1 * rhx
337: d2 = d2 * rhx
338: d3 = d3 * rhy
339: d4 = d4 * rhy
340: d5 = d5 * rhy
341: d6 = d6 * rhx
342: d7 = d7 * rhy
343: d8 = d8 * rhx
345: f1 = sqrt(1.0d0 + d1*d1 + d7*d7)
346: f2 = sqrt(1.0d0 + d1*d1 + d4*d4)
347: f3 = sqrt(1.0d0 + d3*d3 + d8*d8)
348: f4 = sqrt(1.0d0 + d3*d3 + d2*d2)
349: f5 = sqrt(1.0d0 + d2*d2 + d5*d5)
350: f6 = sqrt(1.0d0 + d4*d4 + d6*d6)
352: ft = ft + f2 + f4
354: df1dxc = df1dxc / f1
355: df2dxc = df2dxc / f2
356: df3dxc = df3dxc / f3
357: df4dxc = df4dxc / f4
358: df5dxc = df5dxc / f5
359: df6dxc = df6dxc / f6
361: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
362: & df5dxc + df6dxc)
363: enddo
364: enddo
366: ! Compute triangular areas along the border of the domain.
367: if (xs .eq. 0) then ! left side
368: do j=ys,ys+ym-1
369: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
370: & * rhy
371: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
372: & * rhx
373: ft = ft + sqrt(1.0d0 + d3*d3 + d2*d2)
374: enddo
375: endif
377:
378: if (ys .eq. 0) then !bottom side
379: do i=xs,xs+xm-1
380: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
381: & * rhx
382: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
383: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
384: enddo
385: endif
387:
388: if (xs + xm .eq. mx) then ! right side
389: do j=ys,ys+ym-1
390: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
391: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
392: ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
393: enddo
394: endif
396:
397: if (ys + ym .eq. my) then
398: do i=xs,xs+xm-1
399: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
400: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
401: ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
402: enddo
403: endif
405:
406: if ((ys .eq. 0) .and. (xs .eq. 0)) then
407: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
408: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
409: ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
410: endif
412: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
413: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
414: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
415: ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
416: endif
418: ft = ft * area
419: call MPI_Allreduce(ft,fcn,1,MPI_DOUBLE_PRECISION, &
420: & MPI_SUM,MPI_COMM_WORLD,info)
424: ! Restore vectors
425: call VecRestoreArray(localX,x_v,x_i,info)
426: call VecRestoreArray(localV,g_v,g_i,info)
427: call VecRestoreArray(Left,left_v,left_i,info)
428: call VecRestoreArray(Top,top_v,top_i,info)
429: call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
430: call VecRestoreArray(Right,right_v,right_i,info)
432: ! Scatter values to global vector
433: call DALocalToGlobal(da,localV,INSERT_VALUES,G,info)
435: call PetscLogFlops(70*xm*ym,info)
437: return
438: end !FormFunctionGradient
439:
444: ! ----------------------------------------------------------------------------
445: !
446: !/*
447: ! FormHessian - Evaluates Hessian matrix.
448: !
449: ! Input Parameters:
450: !. tao - the TAO_SOLVER context
451: !. X - input vector
452: !. dummy - not used
453: !
454: ! Output Parameters:
455: !. Hessian - Hessian matrix
456: !. Hpc - optionally different preconditioning matrix
457: !. flag - flag indicating matrix structure
458: !
459: ! Notes:
460: ! Due to mesh point reordering with DAs, we must always work
461: ! with the local mesh points, and then transform them to the new
462: ! global numbering with the local-to-global mapping. We cannot work
463: ! directly with the global numbers for the original uniprocessor mesh!
464: !
465: ! Two methods are available for imposing this transformation
466: ! when setting matrix entries:
467: ! (A) MatSetValuesLocal(), using the local ordering (including
468: ! ghost points!)
469: ! - Do the following two steps once, before calling TaoSolve()
470: ! - Use DAGetISLocalToGlobalMapping() to extract the
471: ! local-to-global map from the DA
472: ! - Associate this map with the matrix by calling
473: ! MatSetLocalToGlobalMapping()
474: ! - Then set matrix entries using the local ordering
475: ! by calling MatSetValuesLocal()
476: ! (B) MatSetValues(), using the global ordering
477: ! - Use DAGetGlobalIndices() to extract the local-to-global map
478: ! - Then apply this map explicitly yourself
479: ! - Set matrix entries using the global ordering by calling
480: ! MatSetValues()
481: ! Option (A) seems cleaner/easier in many cases, and is the procedure
482: ! used in this example.
483: */
484: subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, info)
485: implicit none
487: ! da,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
488: #include "plate2f.h"
489:
490: TAO_SOLVER tao
491: Vec X
492: Mat Hessian,Hpc
493: MatStructure flg
494: integer dummy,info
496: integer i,j,k,row
497: integer xs,xm,gxs,gxm,ys,ym,gys,gym,col(0:6)
498: PetscScalar hx,hy,hydhx,hxdhy,rhx,rhy
499: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8
500: PetscScalar xc,xl,xr,xt,xb,xlt,xrb
501: PetscScalar hl,hr,ht,hb,hc,htl,hbr
503: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
504: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
505: ! will return an array of doubles referenced by x_array offset by x_index.
506: ! i.e., to reference the kth element of X, use x_array(k + x_index).
507: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
508: PetscScalar right_v(0:1),left_v(0:1),bottom_v(0:1),top_v(0:1)
509: PetscScalar x_v(0:1)
510: PetscOffset x_i,right_i,left_i,bottom_i,top_i
511: PetscScalar v(0:6)
512: PetscTruth assembled
513:
514: ! Set various matrix options
515: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,info)
516: call MatSetOption(Hessian,MAT_COLUMNS_SORTED,info)
517: call MatSetOption(Hessian,MAT_ROWS_SORTED,info)
519: ! Get local mesh boundaries
520: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
521: & PETSC_NULL_INTEGER,info)
522: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
523: & PETSC_NULL_INTEGER,info)
525: ! Scatter ghost points to local vectors
526: call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
527: call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
529: ! Get pointers to vector data (see note on Fortran arrays above)
530: call VecGetArray(localX,x_v,x_i,info)
531: call VecGetArray(Top,top_v,top_i,info)
532: call VecGetArray(Bottom,bottom_v,bottom_i,info)
533: call VecGetArray(Left,left_v,left_i,info)
534: call VecGetArray(Right,right_v,right_i,info)
536: ! Initialize matrix entries to zero
537: call MatAssembled(Hessian,assembled,info)
538: if (assembled .eq. PETSC_TRUE) call MatZeroEntries(Hessian,info)
541: rhx = mx + 1.0
542: rhy = my + 1.0
543: hx = 1.0/rhx
544: hy = 1.0/rhy
545: hydhx = hy/hx
546: hxdhy = hx/hy
547: ! compute Hessian over the locally owned part of the mesh
549: do i=xs,xs+xm-1
550: do j=ys,ys+ym-1
551: row = (j-gys)*gxm + (i-gxs)
552:
553: xc = x_v(row + x_i)
554: xt = xc
555: xb = xc
556: xr = xc
557: xl = xc
558: xrb = xc
559: xlt = xc
561: if (i .eq. gxs) then ! Left side
562: xl = left_v(left_i + j - ys + 1)
563: xlt = left_v(left_i + j - ys + 2)
564: else
565: xl = x_v(x_i + row -1 )
566: endif
568: if (j .eq. gys) then ! bottom side
569: xb = bottom_v(bottom_i + i - xs + 1)
570: xrb = bottom_v(bottom_i + i - xs + 2)
571: else
572: xb = x_v(x_i + row - gxm)
573: endif
575: if (i+1 .eq. gxs + gxm) then !right side
576: xr = right_v(right_i + j - ys + 1)
577: xrb = right_v(right_i + j - ys)
578: else
579: xr = x_v(x_i + row + 1)
580: endif
582: if (j+1 .eq. gym+gys) then !top side
583: xt = top_v(top_i +i - xs + 1)
584: xlt = top_v(top_i + i - xs)
585: else
586: xt = x_v(x_i + row + gxm)
587: endif
589: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
590: xlt = x_v(x_i + row - 1 + gxm)
591: endif
592:
593: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
594: xrb = x_v(x_i + row + 1 - gxm)
595: endif
597: d1 = (xc-xl)*rhx
598: d2 = (xc-xr)*rhx
599: d3 = (xc-xt)*rhy
600: d4 = (xc-xb)*rhy
601: d5 = (xrb-xr)*rhy
602: d6 = (xrb-xb)*rhx
603: d7 = (xlt-xl)*rhy
604: d8 = (xlt-xt)*rhx
605:
606: f1 = sqrt( 1.0d0 + d1*d1 + d7*d7)
607: f2 = sqrt( 1.0d0 + d1*d1 + d4*d4)
608: f3 = sqrt( 1.0d0 + d3*d3 + d8*d8)
609: f4 = sqrt( 1.0d0 + d3*d3 + d2*d2)
610: f5 = sqrt( 1.0d0 + d2*d2 + d5*d5)
611: f6 = sqrt( 1.0d0 + d4*d4 + d6*d6)
612:
613:
614: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
615: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
617: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
618: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
620: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
621: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
623: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
624: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
625:
626: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
627: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
628:
629: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
630: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
631: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
632: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
633: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
634: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
635: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
636:
637: hl = hl * 0.5
638: hr = hr * 0.5
639: ht = ht * 0.5
640: hb = hb * 0.5
641: hbr = hbr * 0.5
642: htl = htl * 0.5
643: hc = hc * 0.5
645: k = 0
647: if (j .gt. 0) then
648: v(k) = hb
649: col(k) = row - gxm
650: k=k+1
651: endif
653: if ((j .gt. 0) .and. (i .lt. mx-1)) then
654: v(k) = hbr
655: col(k) = row-gxm+1
656: k=k+1
657: endif
659: if (i .gt. 0) then
660: v(k) = hl
661: col(k) = row - 1
662: k = k+1
663: endif
665: v(k) = hc
666: col(k) = row
667: k=k+1
669: if (i .lt. mx-1) then
670: v(k) = hr
671: col(k) = row + 1
672: k=k+1
673: endif
675: if ((i .gt. 0) .and. (j .lt. my-1)) then
676: v(k) = htl
677: col(k) = row + gxm - 1
678: k=k+1
679: endif
681: if (j .lt. my-1) then
682: v(k) = ht
683: col(k) = row + gxm
684: k=k+1
685: endif
687: ! Set matrix values using local numbering, defined earlier in main routine
688: call MatSetValuesLocal(Hessian,1,row,k,col,v,INSERT_VALUES, &
689: & info)
691:
693: enddo
694: enddo
695:
696: ! restore vectors
697: call VecRestoreArray(localX,x_v,x_i,info)
698: call VecRestoreArray(Left,left_v,left_i,info)
699: call VecRestoreArray(Right,right_v,right_i,info)
700: call VecRestoreArray(Top,top_v,top_i,info)
701: call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
704: ! Assemble the matrix
705: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,info)
706: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,info)
708: call PetscLogFlops(199*xm*ym,info)
710: return
711: end
712:
713:
717: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
719: ! ----------------------------------------------------------------------------
720: !
721: !/*
722: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
723: !
724: !
725: !*/
727: subroutine MSA_BoundaryConditions(info)
728: implicit none
730: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
731: #include "plate2f.h"
733: integer i,j,k,limit,info,maxits
734: integer xs, xm, gxs, gxm, ys, ym, gys, gym
735: integer bsize, lsize, tsize, rsize
736: PetscScalar one, two, three, tol
737: PetscScalar scl,fnorm,det,xt,yt,hx,hy
738: PetscScalar u1,u2,nf1,nf2,njac11,njac12,njac21,njac22
739: PetscScalar b, t, l, r
740: PetscScalar boundary_v(0:1)
741: PetscOffset boundary_i
742: logical exitloop
743: TaoTruth flg
745: limit=0
746: maxits = 5
747: tol=1e-10
748: b=-0.5d0
749: t= 0.5d0
750: l=-0.5d0
751: r= 0.5d0
752: xt=0
753: yt=0
754: one=1.0d0
755: two=2.0d0
756: three=3.0d0
759: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
760: & PETSC_NULL_INTEGER,info)
761: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
762: & gxm,gym,PETSC_NULL_INTEGER,info)
764: bsize = xm + 2
765: lsize = ym + 2
766: rsize = ym + 2
767: tsize = xm + 2
768:
770: call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,info)
771: call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,info)
772: call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,info)
773: call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,info)
775: hx= (r-l)/(mx+1)
776: hy= (t-b)/(my+1)
778: do j=0,3
779:
780: if (j.eq.0) then
781: yt=b
782: xt=l+hx*xs
783: limit=bsize
784: call VecGetArray(Bottom,boundary_v,boundary_i,info)
785:
787: elseif (j.eq.1) then
788: yt=t
789: xt=l+hx*xs
790: limit=tsize
791: call VecGetArray(Top,boundary_v,boundary_i,info)
793: elseif (j.eq.2) then
794: yt=b+hy*ys
795: xt=l
796: limit=lsize
797: call VecGetArray(Left,boundary_v,boundary_i,info)
799: elseif (j.eq.3) then
800: yt=b+hy*ys
801: xt=r
802: limit=rsize
803: call VecGetArray(Right,boundary_v,boundary_i,info)
804: endif
805:
807: do i=0,limit-1
808:
809: u1=xt
810: u2=-yt
811: k = 0
812: exitloop = .false.
813: do while (k .lt. maxits .and. (.not. exitloop) )
815: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
816: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
817: fnorm=sqrt(nf1*nf1+nf2*nf2)
818: if (fnorm .gt. tol) then
819: njac11=one+u2*u2-u1*u1
820: njac12=two*u1*u2
821: njac21=-two*u1*u2
822: njac22=-one - u1*u1 + u2*u2
823: det = njac11*njac22-njac21*njac12
824: u1 = u1-(njac22*nf1-njac12*nf2)/det
825: u2 = u2-(njac11*nf2-njac21*nf1)/det
826: else
827: exitloop = .true.
828: endif
829: k=k+1
830: enddo
832: boundary_v(i + boundary_i) = u1*u1-u2*u2
833: if ((j .eq. 0) .or. (j .eq. 1)) then
834: xt = xt + hx
835: else
836: yt = yt + hy
837: endif
839: enddo
840:
842: if (j.eq.0) then
843: call VecRestoreArray(Bottom,boundary_v,boundary_i,info)
844: elseif (j.eq.1) then
845: call VecRestoreArray(Top,boundary_v,boundary_i,info)
846: elseif (j.eq.2) then
847: call VecRestoreArray(Left,boundary_v,boundary_i,info)
848: elseif (j.eq.3) then
849: call VecRestoreArray(Right,boundary_v,boundary_i,info)
850: endif
851:
852: enddo
855: ! Scale the boundary if desired
856: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom", &
857: & scl,flg,info)
858: if (flg .eq. PETSC_TRUE) then
859: call VecScale(scl,Bottom,info)
860: endif
862: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top", &
863: & scl,flg,info)
864: if (flg .eq. PETSC_TRUE) then
865: call VecScale(scl,Top,info)
866: endif
868: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right", &
869: & scl,flg,info)
870: if (flg .eq. PETSC_TRUE) then
871: call VecScale(scl,Right,info)
872: endif
874: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left", &
875: & scl,flg,info)
876: if (flg .eq. PETSC_TRUE) then
877: call VecScale(scl,Left,info)
878: endif
879:
880:
881: return
882: end
884: ! ----------------------------------------------------------------------------
885: !
886: !/*
887: ! MSA_Plate - Calculates an obstacle for surface to stretch over
888: !
889: ! Output Parameter:
890: !. xl - lower bound vector
891: !. xu - upper bound vector
892: !
893: !*/
895: subroutine MSA_Plate(xl,xu,info)
896: implicit none
898: ! mx,my,bmx,bmy,da,bheight defined in plate2f.h
899: #include "plate2f.h"
900: Vec xl,xu
901: integer i,j,row,info
902: integer xs, xm, ys, ym
903: PetscScalar lb,ub
905: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
906: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
907: ! will return an array of doubles referenced by x_array offset by x_index.
908: ! i.e., to reference the kth element of X, use x_array(k + x_index).
909: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
910: PetscScalar xl_v(0:1)
911: PetscOffset xl_i
914: lb = -1.0d300
915: ub = 1.0d300
917: if (bmy .lt. 0) bmy = 0
918: if (bmy .gt. my) bmy = my
919: if (bmx .lt. 0) bmx = 0
920: if (bmx .gt. mx) bmx = mx
921:
923: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
924: & PETSC_NULL_INTEGER,info)
926: call VecSet(xl,lb,info)
927: call VecSet(xu,ub,info)
929: call VecGetArray(xl,xl_v,xl_i,info)
930:
932: do i=xs,xs+xm-1
934: do j=ys,ys+ym-1
935:
936: row=(j-ys)*xm + (i-xs)
938: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
939: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
940: xl_v(xl_i+row) = bheight
942: endif
944: enddo
945: enddo
948: call VecRestoreArray(xl,xl_v,xl_i,info)
949:
950: return
951: end
955:
956:
957: ! ----------------------------------------------------------------------------
958: !
959: !/*
960: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
961: !
962: ! Output Parameter:
963: !. X - vector for initial guess
964: !
965: !*/
967: subroutine MSA_InitialPoint(X, info)
968: implicit none
970: ! mx,my,localX,da,Top,Left,Bottom,Right defined in plate2f.h
971: #include "plate2f.h"
972: Vec X
974: integer start,i,j,info
975: integer row,xs,xm,gxs,gxm,ys,ym,gys,gym
976: PetscScalar zero, np5
978: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
979: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (integer) x_index, info)
980: ! will return an array of doubles referenced by x_array offset by x_index.
981: ! i.e., to reference the kth element of X, use x_array(k + x_index).
982: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
983: PetscScalar left_v(0:1), right_v(0:1), top_v(0:1)
984: PetscScalar x_v(0:1), bottom_v(0:1)
985: PetscOffset left_i, right_i, top_i, bottom_i, x_i
986: PetscTruth flg
987: PetscRandom rctx
988:
989: zero = 0.0d0
990: np5 = -0.5d0
992: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start", &
993: & start,flg,info)
995: if ((flg .eq. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
996: call VecSet(X,zero,info)
998: elseif ((flg .eq. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
999: call PetscRandomCreate(MPI_COMM_WORLD,RANDOM_DEFAULT,rctx,info)
1000: do i=0,start-1
1001: call VecSetRandom(X,rctx,info)
1002: enddo
1004: call PetscRandomDestroy(rctx,info)
1005: call VecShift(X,np5,info)
1007: else ! average of boundary conditions
1008:
1009: ! Get Local mesh boundaries
1010: call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1011: & PETSC_NULL_INTEGER,info)
1012: call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1013: & PETSC_NULL_INTEGER,info)
1017: ! Get pointers to vector data
1018: call VecGetArray(Top,top_v,top_i,info)
1019: call VecGetArray(Bottom,bottom_v,bottom_i,info)
1020: call VecGetArray(Left,left_v,left_i,info)
1021: call VecGetArray(Right,right_v,right_i,info)
1022:
1023: call VecGetArray(localX,x_v,x_i,info)
1024:
1025: ! Perform local computations
1026: do j=ys,ys+ym-1
1027: do i=xs,xs+xm-1
1028: row = (j-gys)*gxm + (i-gxs)
1029: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1030: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1031: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1032: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1033: enddo
1034: enddo
1036: ! Restore vectors
1037: call VecRestoreArray(localX,x_v,x_i,info)
1039: call VecRestoreArray(Left,left_v,left_i,info)
1040: call VecRestoreArray(Top,top_v,top_i,info)
1041: call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
1042: call VecRestoreArray(Right,right_v,right_i,info)
1044: call DALocalToGlobal(da,localX,INSERT_VALUES,X,info)
1046: endif
1048: return
1049: end