Actual source code: plate2f.F
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: TaoInitialize(); TaoFinalize();
18: ! Routines: TaoCreate();
19: ! Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
20: ! Routines: TaoSetHessianRoutine();
21: ! Routines: TaoSetVariableBoundsRoutine();
22: ! Routines: TaoSetInitialVector();
23: ! Routines: TaoSetFromOptions();
24: ! Routines: TaoSolve();
25: ! Routines: TaoDestroy();
26: ! Processors: n
27: !T*/
31: implicit none
33: #include "plate2f.h"
35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: ! Variable declarations
37: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: !
39: ! Variables:
40: ! (common from plate2f.h):
41: ! Nx, Ny number of processors in x- and y- directions
42: ! mx, my number of grid points in x,y directions
43: ! N global dimension of vector
45: PetscErrorCode ierr ! used to check for functions returning nonzeros
46: Vec x ! solution vector
47: Vec xl, xu ! lower and upper bounds vectorsp
48: PetscInt m ! number of local elements in vector
49: TaoSolver tao ! TaoSolver solver context
50: Mat H ! Hessian matrix
51: ISLocalToGlobalMapping isltog ! local to global mapping object
52: PetscBool flg
53: PetscInt i1,i3,i7
56: external FormFunctionGradient
57: external FormHessian
58: external MSA_BoundaryConditions
59: external MSA_Plate
60: external MSA_InitialPoint
61: ! Initialize Tao
63: i1=1
64: i3=3
65: i7=7
68: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
69: CHKERRQ(ierr)
70: call TaoInitialize(PETSC_NULL_CHARACTER,ierr)
71:
73: ! Specify default dimensions of the problem
74: mx = 10
75: my = 10
76: bheight = 0.1
78: ! Check for any command line arguments that override defaults
79:
80: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mx",mx,flg,ierr)
81: CHKERRQ(ierr)
82: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,ierr)
83: CHKERRQ(ierr)
84:
85: bmx = mx/2
86: bmy = my/2
88: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-bmx",bmx,flg,ierr)
89: CHKERRQ(ierr)
90: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-bmy",bmy,flg,ierr)
91: CHKERRQ(ierr)
92: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bheight",bheight, &
93: & flg,ierr)
94: CHKERRQ(ierr)
95:
97: ! Calculate any derived values from parameters
98: N = mx*my
100: ! Let Petsc determine the dimensions of the local vectors
101: Nx = PETSC_DECIDE
102: NY = PETSC_DECIDE
104: ! A two dimensional distributed array will help define this problem, which
105: ! derives from an elliptic PDE on a two-dimensional domain. From the
106: ! distributed array, create the vectors
108: call DMDACreate2d(MPI_COMM_WORLD,DMDA_BOUNDARY_NONE, &
109: & DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
110: & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
111: & dm,ierr)
112: CHKERRQ(ierr)
113:
115: ! Extract global and local vectors from DM; The local vectors are
116: ! used solely as work space for the evaluation of the function,
117: ! gradient, and Hessian. Duplicate for remaining vectors that are
118: ! the same types.
120: call DMCreateGlobalVector(dm,x,ierr)
121: CHKERRQ(ierr)
122: call DMCreateLocalVector(dm,localX,ierr)
123: CHKERRQ(ierr)
124: call VecDuplicate(localX,localV,ierr)
125: CHKERRQ(ierr)
127: ! Create a matrix data structure to store the Hessian.
128: ! Here we (optionally) also associate the local numbering scheme
129: ! with the matrix so that later we can use local indices for matrix
130: ! assembly
132: call VecGetLocalSize(x,m,ierr)
133: CHKERRQ(ierr)
134: call MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, &
135: & i3,PETSC_NULL_INTEGER,H,ierr)
136: CHKERRQ(ierr)
138: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
139: CHKERRQ(ierr)
140: call DMGetLocalToGlobalMapping(dm,isltog,ierr)
141: CHKERRQ(ierr)
142: call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
143: CHKERRQ(ierr)
144:
146: ! The Tao code begins here
147: ! Create TAO solver and set desired solution method.
148: ! This problems uses bounded variables, so the
149: ! method must either be 'tao_tron' or 'tao_blmvm'
151: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
152: CHKERRQ(ierr)
153: call TaoSetType(tao,'tao_blmvm',ierr)
154: CHKERRQ(ierr)
156: ! Set minimization function and gradient, hessian evaluation functions
158: call TaoSetObjectiveAndGradientRoutine(tao, &
159: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
160: CHKERRQ(ierr)
162: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
163: & PETSC_NULL_OBJECT, ierr)
164: CHKERRQ(ierr)
166: ! Set Variable bounds
167: call MSA_BoundaryConditions(ierr)
168: CHKERRQ(ierr)
169: call TaoSetVariableBoundsRoutine(tao,MSA_Plate, &
170: & PETSC_NULL_OBJECT,ierr)
171: CHKERRQ(ierr)
173: ! Set the initial solution guess
174: call MSA_InitialPoint(x, ierr)
175: CHKERRQ(ierr)
176: call TaoSetInitialVector(tao,x,ierr)
177: CHKERRQ(ierr)
179: ! Check for any tao command line options
180: call TaoSetFromOptions(tao,ierr)
181: CHKERRQ(ierr)
183: ! Solve the application
184: call TaoSolve(tao,ierr)
185: CHKERRQ(ierr)
187: ! Free TAO data structures
188: call TaoDestroy(tao,ierr)
189: CHKERRQ(ierr)
191: ! Free PETSc data structures
192: call VecDestroy(x,ierr)
193: call VecDestroy(Top,ierr)
194: call VecDestroy(Bottom,ierr)
195: call VecDestroy(Left,ierr)
196: call VecDestroy(Right,ierr)
197: call MatDestroy(H,ierr)
198: call VecDestroy(localX,ierr)
199: call VecDestroy(localV,ierr)
200: call DMDestroy(dm,ierr)
202: ! Finalize TAO
204: call TaoFinalize(ierr)
205: call PetscFinalize(ierr)
207: end
209: ! ---------------------------------------------------------------------
210: !
211: ! FormFunctionGradient - Evaluates function f(X).
212: !
213: ! Input Parameters:
214: ! tao - the TaoSolver context
215: ! X - the input vector
216: ! dummy - optional user-defined context, as set by TaoSetFunction()
217: ! (not used here)
218: !
219: ! Output Parameters:
220: ! fcn - the newly evaluated function
221: ! G - the gradient vector
222: ! info - error code
223: !
226: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
227: implicit none
229: ! dm, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
230: #include "plate2f.h"
231:
232: ! Input/output variables
234: TaoSolver tao
235: PetscReal fcn
236: Vec X, G
237: PetscErrorCode ierr
238: PetscInt dummy
239:
240: PetscInt i,j,row
241: PetscInt xs, xm
242: PetscInt gxs, gxm
243: PetscInt ys, ym
244: PetscInt gys, gym
245: PetscReal ft,zero,hx,hy,hydhx,hxdhy
246: PetscReal area,rhx,rhy
247: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
248: PetscReal d4,d5,d6,d7,d8
249: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
250: PetscReal df5dxc,df6dxc
251: PetscReal xc,xl,xr,xt,xb,xlt,xrb
254: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
255: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
256: ! will return an array of doubles referenced by x_array offset by x_index.
257: ! i.e., to reference the kth element of X, use x_array(k + x_index).
258: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
259: PetscReal g_v(0:1),x_v(0:1)
260: PetscReal top_v(0:1),left_v(0:1)
261: PetscReal right_v(0:1),bottom_v(0:1)
262: PetscOffset g_i,left_i,right_i
263: PetscOffset bottom_i,top_i,x_i
265: ft = 0.0d0
266: zero = 0.0d0
267: hx = 1.0d0/(mx + 1)
268: hy = 1.0d0/(my + 1)
269: hydhx = hy/hx
270: hxdhy = hx/hy
271: area = 0.5d0 * hx * hy
272: rhx = mx + 1.0d0
273: rhy = my + 1.0d0
276: ! Get local mesh boundaries
277: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
278: & PETSC_NULL_INTEGER,ierr)
279: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
280: & gxm,gym,PETSC_NULL_INTEGER,ierr)
282: ! Scatter ghost points to local vector
283: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
284: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
286: ! Initialize the vector to zero
287: call VecSet(localV,zero,ierr)
289: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
290: call VecGetArray(localX,x_v,x_i,ierr)
291: call VecGetArray(localV,g_v,g_i,ierr)
292: call VecGetArray(Top,top_v,top_i,ierr)
293: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
294: call VecGetArray(Left,left_v,left_i,ierr)
295: call VecGetArray(Right,right_v,right_i,ierr)
297: ! Compute function over the locally owned part of the mesh
298: do j = ys,ys+ym-1
299: do i = xs,xs+xm-1
300: row = (j-gys)*gxm + (i-gxs)
301: xc = x_v(row+x_i)
302: xt = xc
303: xb = xc
304: xr = xc
305: xl = xc
306: xrb = xc
307: xlt = xc
309: if (i .eq. 0) then !left side
310: xl = left_v(j - ys + 1 + left_i)
311: xlt = left_v(j - ys + 2 + left_i)
312: else
313: xl = x_v(row - 1 + x_i)
314: endif
316: if (j .eq. 0) then !bottom side
317: xb = bottom_v(i - xs + 1 + bottom_i)
318: xrb = bottom_v(i - xs + 2 + bottom_i)
319: else
320: xb = x_v(row - gxm + x_i)
321: endif
323: if (i + 1 .eq. gxs + gxm) then !right side
324: xr = right_v(j - ys + 1 + right_i)
325: xrb = right_v(j - ys + right_i)
326: else
327: xr = x_v(row + 1 + x_i)
328: endif
330: if (j + 1 .eq. gys + gym) then !top side
331: xt = top_v(i - xs + 1 + top_i)
332: xlt = top_v(i - xs + top_i)
333: else
334: xt = x_v(row + gxm + x_i)
335: endif
337: if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
338: xlt = x_v(row - 1 + gxm + x_i)
339: endif
341: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
342: xrb = x_v(row + 1 - gxm + x_i)
343: endif
345: d1 = xc-xl
346: d2 = xc-xr
347: d3 = xc-xt
348: d4 = xc-xb
349: d5 = xr-xrb
350: d6 = xrb-xb
351: d7 = xlt-xl
352: d8 = xt-xlt
354: df1dxc = d1 * hydhx
355: df2dxc = d1 * hydhx + d4 * hxdhy
356: df3dxc = d3 * hxdhy
357: df4dxc = d2 * hydhx + d3 * hxdhy
358: df5dxc = d2 * hydhx
359: df6dxc = d4 * hxdhy
361: d1 = d1 * rhx
362: d2 = d2 * rhx
363: d3 = d3 * rhy
364: d4 = d4 * rhy
365: d5 = d5 * rhy
366: d6 = d6 * rhx
367: d7 = d7 * rhy
368: d8 = d8 * rhx
370: f1 = sqrt(1.0d0 + d1*d1 + d7*d7)
371: f2 = sqrt(1.0d0 + d1*d1 + d4*d4)
372: f3 = sqrt(1.0d0 + d3*d3 + d8*d8)
373: f4 = sqrt(1.0d0 + d3*d3 + d2*d2)
374: f5 = sqrt(1.0d0 + d2*d2 + d5*d5)
375: f6 = sqrt(1.0d0 + d4*d4 + d6*d6)
377: ft = ft + f2 + f4
379: df1dxc = df1dxc / f1
380: df2dxc = df2dxc / f2
381: df3dxc = df3dxc / f3
382: df4dxc = df4dxc / f4
383: df5dxc = df5dxc / f5
384: df6dxc = df6dxc / f6
386: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
387: & df5dxc + df6dxc)
388: enddo
389: enddo
391: ! Compute triangular areas along the border of the domain.
392: if (xs .eq. 0) then ! left side
393: do j=ys,ys+ym-1
394: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
395: & * rhy
396: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
397: & * rhx
398: ft = ft + sqrt(1.0d0 + d3*d3 + d2*d2)
399: enddo
400: endif
402:
403: if (ys .eq. 0) then !bottom side
404: do i=xs,xs+xm-1
405: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
406: & * rhx
407: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
408: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
409: enddo
410: endif
412:
413: if (xs + xm .eq. mx) then ! right side
414: do j=ys,ys+ym-1
415: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
416: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
417: ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
418: enddo
419: endif
421:
422: if (ys + ym .eq. my) then
423: do i=xs,xs+xm-1
424: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
425: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
426: ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
427: enddo
428: endif
430:
431: if ((ys .eq. 0) .and. (xs .eq. 0)) then
432: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
433: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
434: ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
435: endif
437: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
438: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
439: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
440: ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
441: endif
443: ft = ft * area
444: call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, &
445: & MPI_SUM,MPI_COMM_WORLD,ierr)
449: ! Restore vectors
450: call VecRestoreArray(localX,x_v,x_i,ierr)
451: call VecRestoreArray(localV,g_v,g_i,ierr)
452: call VecRestoreArray(Left,left_v,left_i,ierr)
453: call VecRestoreArray(Top,top_v,top_i,ierr)
454: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
455: call VecRestoreArray(Right,right_v,right_i,ierr)
457: ! Scatter values to global vector
458: call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
459: call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)
461: call PetscLogFlops(70.0d0*xm*ym,ierr)
463: return
464: end !FormFunctionGradient
465:
470: ! ----------------------------------------------------------------------------
471: !
472: !/*
473: ! FormHessian - Evaluates Hessian matrix.
474: !
475: ! Input Parameters:
476: !. tao - the TaoSolver context
477: !. X - input vector
478: !. dummy - not used
479: !
480: ! Output Parameters:
481: !. Hessian - Hessian matrix
482: !. Hpc - optionally different preconditioning matrix
483: !. flag - flag indicating matrix structure
484: !
485: ! Notes:
486: ! Due to mesh point reordering with DMs, we must always work
487: ! with the local mesh points, and then transform them to the new
488: ! global numbering with the local-to-global mapping. We cannot work
489: ! directly with the global numbers for the original uniprocessor mesh!
490: !
491: ! Two methods are available for imposing this transformation
492: ! when setting matrix entries:
493: ! (A) MatSetValuesLocal(), using the local ordering (including
494: ! ghost points!)
495: ! - Do the following two steps once, before calling TaoSolve()
496: ! - Use DMDAGetISLocalToGlobalMapping() to extract the
497: ! local-to-global map from the DM
498: ! - Associate this map with the matrix by calling
499: ! MatSetLocalToGlobalMapping()
500: ! - Then set matrix entries using the local ordering
501: ! by calling MatSetValuesLocal()
502: ! (B) MatSetValues(), using the global ordering
503: ! - Use DMDAGetGlobalIndices() to extract the local-to-global map
504: ! - Then apply this map explicitly yourself
505: ! - Set matrix entries using the global ordering by calling
506: ! MatSetValues()
507: ! Option (A) seems cleaner/easier in many cases, and is the procedure
508: ! used in this example.
509: */
510: subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, ierr)
511: implicit none
513: ! dm,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
514: #include "plate2f.h"
515:
516: TaoSolver tao
517: Vec X
518: Mat Hessian,Hpc
519: MatStructure flg
520: PetscInt dummy
521: PetscErrorCode ierr
523: PetscInt i,j,k,row
524: PetscInt xs,xm,gxs,gxm
525: PetscInt ys,ym,gys,gym
526: PetscInt col(0:6)
527: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
528: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
529: PetscReal d4,d5,d6,d7,d8
530: PetscReal xc,xl,xr,xt,xb,xlt,xrb
531: PetscReal hl,hr,ht,hb,hc,htl,hbr
533: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
534: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
535: ! will return an array of doubles referenced by x_array offset by x_index.
536: ! i.e., to reference the kth element of X, use x_array(k + x_index).
537: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
538: PetscReal right_v(0:1),left_v(0:1)
539: PetscReal bottom_v(0:1),top_v(0:1)
540: PetscReal x_v(0:1)
541: PetscOffset x_i,right_i,left_i
542: PetscOffset bottom_i,top_i
543: PetscReal v(0:6)
544: PetscBool assembled
545: PetscInt i1
546:
547: i1=1
548:
549: ! Set various matrix options
550: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, &
551: & PETSC_TRUE,ierr)
553: ! Get local mesh boundaries
554: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
555: & PETSC_NULL_INTEGER,ierr)
556: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
557: & PETSC_NULL_INTEGER,ierr)
559: ! Scatter ghost points to local vectors
560: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
561: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
563: ! Get pointers to vector data (see note on Fortran arrays above)
564: call VecGetArray(localX,x_v,x_i,ierr)
565: call VecGetArray(Top,top_v,top_i,ierr)
566: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
567: call VecGetArray(Left,left_v,left_i,ierr)
568: call VecGetArray(Right,right_v,right_i,ierr)
570: ! Initialize matrix entries to zero
571: call MatAssembled(Hessian,assembled,ierr)
572: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)
575: rhx = mx + 1.0
576: rhy = my + 1.0
577: hx = 1.0/rhx
578: hy = 1.0/rhy
579: hydhx = hy/hx
580: hxdhy = hx/hy
581: ! compute Hessian over the locally owned part of the mesh
583: do i=xs,xs+xm-1
584: do j=ys,ys+ym-1
585: row = (j-gys)*gxm + (i-gxs)
586:
587: xc = x_v(row + x_i)
588: xt = xc
589: xb = xc
590: xr = xc
591: xl = xc
592: xrb = xc
593: xlt = xc
595: if (i .eq. gxs) then ! Left side
596: xl = left_v(left_i + j - ys + 1)
597: xlt = left_v(left_i + j - ys + 2)
598: else
599: xl = x_v(x_i + row -1 )
600: endif
602: if (j .eq. gys) then ! bottom side
603: xb = bottom_v(bottom_i + i - xs + 1)
604: xrb = bottom_v(bottom_i + i - xs + 2)
605: else
606: xb = x_v(x_i + row - gxm)
607: endif
609: if (i+1 .eq. gxs + gxm) then !right side
610: xr = right_v(right_i + j - ys + 1)
611: xrb = right_v(right_i + j - ys)
612: else
613: xr = x_v(x_i + row + 1)
614: endif
616: if (j+1 .eq. gym+gys) then !top side
617: xt = top_v(top_i +i - xs + 1)
618: xlt = top_v(top_i + i - xs)
619: else
620: xt = x_v(x_i + row + gxm)
621: endif
623: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
624: xlt = x_v(x_i + row - 1 + gxm)
625: endif
626:
627: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
628: xrb = x_v(x_i + row + 1 - gxm)
629: endif
631: d1 = (xc-xl)*rhx
632: d2 = (xc-xr)*rhx
633: d3 = (xc-xt)*rhy
634: d4 = (xc-xb)*rhy
635: d5 = (xrb-xr)*rhy
636: d6 = (xrb-xb)*rhx
637: d7 = (xlt-xl)*rhy
638: d8 = (xlt-xt)*rhx
639:
640: f1 = sqrt( 1.0d0 + d1*d1 + d7*d7)
641: f2 = sqrt( 1.0d0 + d1*d1 + d4*d4)
642: f3 = sqrt( 1.0d0 + d3*d3 + d8*d8)
643: f4 = sqrt( 1.0d0 + d3*d3 + d2*d2)
644: f5 = sqrt( 1.0d0 + d2*d2 + d5*d5)
645: f6 = sqrt( 1.0d0 + d4*d4 + d6*d6)
646:
647:
648: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
649: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
651: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
652: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
654: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
655: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
657: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
658: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
659:
660: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
661: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
662:
663: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
664: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
665: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
666: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
667: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
668: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
669: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
670:
671: hl = hl * 0.5
672: hr = hr * 0.5
673: ht = ht * 0.5
674: hb = hb * 0.5
675: hbr = hbr * 0.5
676: htl = htl * 0.5
677: hc = hc * 0.5
679: k = 0
681: if (j .gt. 0) then
682: v(k) = hb
683: col(k) = row - gxm
684: k=k+1
685: endif
687: if ((j .gt. 0) .and. (i .lt. mx-1)) then
688: v(k) = hbr
689: col(k) = row-gxm+1
690: k=k+1
691: endif
693: if (i .gt. 0) then
694: v(k) = hl
695: col(k) = row - 1
696: k = k+1
697: endif
699: v(k) = hc
700: col(k) = row
701: k=k+1
703: if (i .lt. mx-1) then
704: v(k) = hr
705: col(k) = row + 1
706: k=k+1
707: endif
709: if ((i .gt. 0) .and. (j .lt. my-1)) then
710: v(k) = htl
711: col(k) = row + gxm - 1
712: k=k+1
713: endif
715: if (j .lt. my-1) then
716: v(k) = ht
717: col(k) = row + gxm
718: k=k+1
719: endif
721: ! Set matrix values using local numbering, defined earlier in main routine
722: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
723: & ierr)
725:
727: enddo
728: enddo
729:
730: ! restore vectors
731: call VecRestoreArray(localX,x_v,x_i,ierr)
732: call VecRestoreArray(Left,left_v,left_i,ierr)
733: call VecRestoreArray(Right,right_v,right_i,ierr)
734: call VecRestoreArray(Top,top_v,top_i,ierr)
735: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
738: ! Assemble the matrix
739: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
740: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)
742: call PetscLogFlops(199.0d0*xm*ym,ierr)
744: return
745: end
746:
747:
751: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
753: ! ----------------------------------------------------------------------------
754: !
755: !/*
756: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
757: !
758: !
759: !*/
761: subroutine MSA_BoundaryConditions(ierr)
762: implicit none
764: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
765: #include "plate2f.h"
767: PetscErrorCode ierr
768: PetscInt i,j,k,limit,maxits
769: PetscInt xs, xm, gxs, gxm
770: PetscInt ys, ym, gys, gym
771: PetscInt bsize, lsize
772: PetscInt tsize, rsize
773: PetscReal one,two,three,tol
774: PetscReal scl,fnorm,det,xt
775: PetscReal yt,hx,hy,u1,u2,nf1,nf2
776: PetscReal njac11,njac12,njac21,njac22
777: PetscReal b, t, l, r
778: PetscReal boundary_v(0:1)
779: PetscOffset boundary_i
780: logical exitloop
781: PetscBool flg
783: limit=0
784: maxits = 5
785: tol=1e-10
786: b=-0.5d0
787: t= 0.5d0
788: l=-0.5d0
789: r= 0.5d0
790: xt=0
791: yt=0
792: one=1.0d0
793: two=2.0d0
794: three=3.0d0
797: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
798: & PETSC_NULL_INTEGER,ierr)
799: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
800: & gxm,gym,PETSC_NULL_INTEGER,ierr)
802: bsize = xm + 2
803: lsize = ym + 2
804: rsize = ym + 2
805: tsize = xm + 2
806:
808: call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
809: call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
810: call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
811: call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)
813: hx= (r-l)/(mx+1)
814: hy= (t-b)/(my+1)
816: do j=0,3
817:
818: if (j.eq.0) then
819: yt=b
820: xt=l+hx*xs
821: limit=bsize
822: call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
823:
825: elseif (j.eq.1) then
826: yt=t
827: xt=l+hx*xs
828: limit=tsize
829: call VecGetArray(Top,boundary_v,boundary_i,ierr)
831: elseif (j.eq.2) then
832: yt=b+hy*ys
833: xt=l
834: limit=lsize
835: call VecGetArray(Left,boundary_v,boundary_i,ierr)
837: elseif (j.eq.3) then
838: yt=b+hy*ys
839: xt=r
840: limit=rsize
841: call VecGetArray(Right,boundary_v,boundary_i,ierr)
842: endif
843:
845: do i=0,limit-1
846:
847: u1=xt
848: u2=-yt
849: k = 0
850: exitloop = .false.
851: do while (k .lt. maxits .and. (.not. exitloop) )
853: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
854: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
855: fnorm=sqrt(nf1*nf1+nf2*nf2)
856: if (fnorm .gt. tol) then
857: njac11=one+u2*u2-u1*u1
858: njac12=two*u1*u2
859: njac21=-two*u1*u2
860: njac22=-one - u1*u1 + u2*u2
861: det = njac11*njac22-njac21*njac12
862: u1 = u1-(njac22*nf1-njac12*nf2)/det
863: u2 = u2-(njac11*nf2-njac21*nf1)/det
864: else
865: exitloop = .true.
866: endif
867: k=k+1
868: enddo
870: boundary_v(i + boundary_i) = u1*u1-u2*u2
871: if ((j .eq. 0) .or. (j .eq. 1)) then
872: xt = xt + hx
873: else
874: yt = yt + hy
875: endif
877: enddo
878:
880: if (j.eq.0) then
881: call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
882: elseif (j.eq.1) then
883: call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
884: elseif (j.eq.2) then
885: call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
886: elseif (j.eq.3) then
887: call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
888: endif
889:
890: enddo
893: ! Scale the boundary if desired
894: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom", &
895: & scl,flg,ierr)
896: if (flg .eqv. PETSC_TRUE) then
897: call VecScale(scl,Bottom,ierr)
898: endif
900: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top", &
901: & scl,flg,ierr)
902: if (flg .eqv. PETSC_TRUE) then
903: call VecScale(scl,Top,ierr)
904: endif
906: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right", &
907: & scl,flg,ierr)
908: if (flg .eqv. PETSC_TRUE) then
909: call VecScale(scl,Right,ierr)
910: endif
912: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left", &
913: & scl,flg,ierr)
914: if (flg .eqv. PETSC_TRUE) then
915: call VecScale(scl,Left,ierr)
916: endif
917:
918:
919: return
920: end
922: ! ----------------------------------------------------------------------------
923: !
924: !/*
925: ! MSA_Plate - Calculates an obstacle for surface to stretch over
926: !
927: ! Output Parameter:
928: !. xl - lower bound vector
929: !. xu - upper bound vector
930: !
931: !*/
933: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
934: implicit none
935: ! mx,my,bmx,bmy,dm,bheight defined in plate2f.h
936: #include "plate2f.h"
937: TaoSolver tao
938: Vec xl,xu
939: PetscErrorCode ierr
940: PetscInt i,j,row
941: PetscInt xs, xm, ys, ym
942: PetscReal lb,ub
943: PetscInt dummy
945: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
946: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
947: ! will return an array of doubles referenced by x_array offset by x_index.
948: ! i.e., to reference the kth element of X, use x_array(k + x_index).
949: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
950: PetscReal xl_v(0:1)
951: PetscOffset xl_i
953: print *,'msa_plate'
955: lb = -1.0d300
956: ub = 1.0d300
958: if (bmy .lt. 0) bmy = 0
959: if (bmy .gt. my) bmy = my
960: if (bmx .lt. 0) bmx = 0
961: if (bmx .gt. mx) bmx = mx
962:
964: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
965: & PETSC_NULL_INTEGER,ierr)
967: call VecSet(xl,lb,ierr)
968: call VecSet(xu,ub,ierr)
970: call VecGetArray(xl,xl_v,xl_i,ierr)
971:
973: do i=xs,xs+xm-1
975: do j=ys,ys+ym-1
976:
977: row=(j-ys)*xm + (i-xs)
979: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
980: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
981: xl_v(xl_i+row) = bheight
983: endif
985: enddo
986: enddo
989: call VecRestoreArray(xl,xl_v,xl_i,ierr)
990:
991: return
992: end
996:
997:
998: ! ----------------------------------------------------------------------------
999: !
1000: !/*
1001: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
1002: !
1003: ! Output Parameter:
1004: !. X - vector for initial guess
1005: !
1006: !*/
1008: subroutine MSA_InitialPoint(X, ierr)
1009: implicit none
1011: ! mx,my,localX,dm,Top,Left,Bottom,Right defined in plate2f.h
1012: #include "plate2f.h"
1013: Vec X
1014: PetscErrorCode ierr
1015: PetscInt start,i,j
1016: PetscInt row
1017: PetscInt xs,xm,gxs,gxm
1018: PetscInt ys,ym,gys,gym
1019: PetscReal zero, np5
1021: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1022: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
1023: ! will return an array of doubles referenced by x_array offset by x_index.
1024: ! i.e., to reference the kth element of X, use x_array(k + x_index).
1025: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1026: PetscReal left_v(0:1),right_v(0:1)
1027: PetscReal bottom_v(0:1),top_v(0:1)
1028: PetscReal x_v(0:1)
1029: PetscOffset left_i, right_i, top_i
1030: PetscOffset bottom_i,x_i
1031: PetscBool flg
1032: PetscRandom rctx
1033:
1034: zero = 0.0d0
1035: np5 = -0.5d0
1037: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start", &
1038: & start,flg,ierr)
1040: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
1041: call VecSet(X,zero,ierr)
1043: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
1044: call PetscRandomCreate(MPI_COMM_WORLD,rctx,ierr)
1045: do i=0,start-1
1046: call VecSetRandom(X,rctx,ierr)
1047: enddo
1049: call PetscRandomDestroy(rctx,ierr)
1050: call VecShift(X,np5,ierr)
1052: else ! average of boundary conditions
1053:
1054: ! Get Local mesh boundaries
1055: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1056: & PETSC_NULL_INTEGER,ierr)
1057: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1058: & PETSC_NULL_INTEGER,ierr)
1062: ! Get pointers to vector data
1063: call VecGetArray(Top,top_v,top_i,ierr)
1064: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1065: call VecGetArray(Left,left_v,left_i,ierr)
1066: call VecGetArray(Right,right_v,right_i,ierr)
1067:
1068: call VecGetArray(localX,x_v,x_i,ierr)
1069:
1070: ! Perform local computations
1071: do j=ys,ys+ym-1
1072: do i=xs,xs+xm-1
1073: row = (j-gys)*gxm + (i-gxs)
1074: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1075: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1076: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1077: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1078: enddo
1079: enddo
1081: ! Restore vectors
1082: call VecRestoreArray(localX,x_v,x_i,ierr)
1084: call VecRestoreArray(Left,left_v,left_i,ierr)
1085: call VecRestoreArray(Top,top_v,top_i,ierr)
1086: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1087: call VecRestoreArray(Right,right_v,right_i,ierr)
1089: call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1090: call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)
1092: endif
1094: return
1095: end