Actual source code: plate2f.F90
petsc-3.8.4 2018-03-24
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*/
28: #include "plate2f.h"
30: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: ! Variable declarations
32: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
33: !
34: ! Variables:
35: ! (common from plate2f.h):
36: ! Nx, Ny number of processors in x- and y- directions
37: ! mx, my number of grid points in x,y directions
38: ! N global dimension of vector
40: PetscErrorCode ierr ! used to check for functions returning nonzeros
41: Vec x ! solution vector
42: PetscInt m ! number of local elements in vector
43: Tao tao ! Tao solver context
44: Mat H ! Hessian matrix
45: ISLocalToGlobalMapping isltog ! local to global mapping object
46: PetscBool flg
47: PetscInt i1,i3,i7
50: external FormFunctionGradient
51: external FormHessian
52: external MSA_BoundaryConditions
53: external MSA_Plate
54: external MSA_InitialPoint
55: ! Initialize Tao
57: i1=1
58: i3=3
59: i7=7
62: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
63: if (ierr .ne. 0) then
64: print*,'Unable to initialize PETSc'
65: stop
66: endif
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_OPTIONS,PETSC_NULL_CHARACTER, &
76: & '-mx',mx,flg,ierr)
77: CHKERRA(ierr)
78: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
79: & '-my',my,flg,ierr)
80: CHKERRA(ierr)
82: bmx = mx/2
83: bmy = my/2
85: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
86: & '-bmx',bmx,flg,ierr)
87: CHKERRA(ierr)
88: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
89: & '-bmy',bmy,flg,ierr)
90: CHKERRA(ierr)
91: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
92: & '-bheight',bheight,flg,ierr)
93: CHKERRA(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: CHKERRA(ierr)
112: call DMSetFromOptions(dm,ierr)
113: call DMSetUp(dm,ierr)
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: CHKERRA(ierr)
122: call DMCreateLocalVector(dm,localX,ierr)
123: CHKERRA(ierr)
124: call VecDuplicate(localX,localV,ierr)
125: CHKERRA(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: CHKERRA(ierr)
134: call MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, &
135: & i3,PETSC_NULL_INTEGER,H,ierr)
136: CHKERRA(ierr)
138: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
139: CHKERRA(ierr)
140: call DMGetLocalToGlobalMapping(dm,isltog,ierr)
141: CHKERRA(ierr)
142: call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
143: CHKERRA(ierr)
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: CHKERRA(ierr)
153: call TaoSetType(tao,TAOBLMVM,ierr)
154: CHKERRA(ierr)
156: ! Set minimization function and gradient, hessian evaluation functions
158: call TaoSetObjectiveAndGradientRoutine(tao, &
159: & FormFunctionGradient,0,ierr)
160: CHKERRA(ierr)
162: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
163: & 0, ierr)
164: CHKERRA(ierr)
166: ! Set Variable bounds
167: call MSA_BoundaryConditions(ierr)
168: CHKERRA(ierr)
169: call TaoSetVariableBoundsRoutine(tao,MSA_Plate, &
170: & 0,ierr)
171: CHKERRA(ierr)
173: ! Set the initial solution guess
174: call MSA_InitialPoint(x, ierr)
175: CHKERRA(ierr)
176: call TaoSetInitialVector(tao,x,ierr)
177: CHKERRA(ierr)
179: ! Check for any tao command line options
180: call TaoSetFromOptions(tao,ierr)
181: CHKERRA(ierr)
183: ! Solve the application
184: call TaoSolve(tao,ierr)
185: CHKERRA(ierr)
187: ! Free TAO data structures
188: call TaoDestroy(tao,ierr)
189: CHKERRA(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 PetscFinalize(ierr)
206: end
208: ! ---------------------------------------------------------------------
209: !
210: ! FormFunctionGradient - Evaluates function f(X).
211: !
212: ! Input Parameters:
213: ! tao - the Tao context
214: ! X - the input vector
215: ! dummy - optional user-defined context, as set by TaoSetFunction()
216: ! (not used here)
217: !
218: ! Output Parameters:
219: ! fcn - the newly evaluated function
220: ! G - the gradient vector
221: ! info - error code
222: !
225: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
226: #include "plate2f.h"
228: ! Input/output variables
230: Tao tao
231: PetscReal fcn
232: Vec X, G
233: PetscErrorCode ierr
234: PetscInt dummy
236: PetscInt i,j,row
237: PetscInt xs, xm
238: PetscInt gxs, gxm
239: PetscInt ys, ym
240: PetscInt gys, gym
241: PetscReal ft,zero,hx,hy,hydhx,hxdhy
242: PetscReal area,rhx,rhy
243: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
244: PetscReal d4,d5,d6,d7,d8
245: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
246: PetscReal df5dxc,df6dxc
247: PetscReal xc,xl,xr,xt,xb,xlt,xrb
250: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
251: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
252: ! will return an array of doubles referenced by x_array offset by x_index.
253: ! i.e., to reference the kth element of X, use x_array(k + x_index).
254: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
255: PetscReal g_v(0:1),x_v(0:1)
256: PetscReal top_v(0:1),left_v(0:1)
257: PetscReal right_v(0:1),bottom_v(0:1)
258: PetscOffset g_i,left_i,right_i
259: PetscOffset bottom_i,top_i,x_i
261: ft = 0.0
262: zero = 0.0
263: hx = 1.0/real(mx + 1)
264: hy = 1.0/real(my + 1)
265: hydhx = hy/hx
266: hxdhy = hx/hy
267: area = 0.5 * hx * hy
268: rhx = real(mx) + 1.0
269: rhy = real(my) + 1.0
272: ! Get local mesh boundaries
273: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
274: & PETSC_NULL_INTEGER,ierr)
275: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
276: & gxm,gym,PETSC_NULL_INTEGER,ierr)
278: ! Scatter ghost points to local vector
279: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
280: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
282: ! Initialize the vector to zero
283: call VecSet(localV,zero,ierr)
285: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
286: call VecGetArray(localX,x_v,x_i,ierr)
287: call VecGetArray(localV,g_v,g_i,ierr)
288: call VecGetArray(Top,top_v,top_i,ierr)
289: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
290: call VecGetArray(Left,left_v,left_i,ierr)
291: call VecGetArray(Right,right_v,right_i,ierr)
293: ! Compute function over the locally owned part of the mesh
294: do j = ys,ys+ym-1
295: do i = xs,xs+xm-1
296: row = (j-gys)*gxm + (i-gxs)
297: xc = x_v(row+x_i)
298: xt = xc
299: xb = xc
300: xr = xc
301: xl = xc
302: xrb = xc
303: xlt = xc
305: if (i .eq. 0) then !left side
306: xl = left_v(j - ys + 1 + left_i)
307: xlt = left_v(j - ys + 2 + left_i)
308: else
309: xl = x_v(row - 1 + x_i)
310: endif
312: if (j .eq. 0) then !bottom side
313: xb = bottom_v(i - xs + 1 + bottom_i)
314: xrb = bottom_v(i - xs + 2 + bottom_i)
315: else
316: xb = x_v(row - gxm + x_i)
317: endif
319: if (i + 1 .eq. gxs + gxm) then !right side
320: xr = right_v(j - ys + 1 + right_i)
321: xrb = right_v(j - ys + right_i)
322: else
323: xr = x_v(row + 1 + x_i)
324: endif
326: if (j + 1 .eq. gys + gym) then !top side
327: xt = top_v(i - xs + 1 + top_i)
328: xlt = top_v(i - xs + top_i)
329: else
330: xt = x_v(row + gxm + x_i)
331: endif
333: if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
334: xlt = x_v(row - 1 + gxm + x_i)
335: endif
337: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
338: xrb = x_v(row + 1 - gxm + x_i)
339: endif
341: d1 = xc-xl
342: d2 = xc-xr
343: d3 = xc-xt
344: d4 = xc-xb
345: d5 = xr-xrb
346: d6 = xrb-xb
347: d7 = xlt-xl
348: d8 = xt-xlt
350: df1dxc = d1 * hydhx
351: df2dxc = d1 * hydhx + d4 * hxdhy
352: df3dxc = d3 * hxdhy
353: df4dxc = d2 * hydhx + d3 * hxdhy
354: df5dxc = d2 * hydhx
355: df6dxc = d4 * hxdhy
357: d1 = d1 * rhx
358: d2 = d2 * rhx
359: d3 = d3 * rhy
360: d4 = d4 * rhy
361: d5 = d5 * rhy
362: d6 = d6 * rhx
363: d7 = d7 * rhy
364: d8 = d8 * rhx
366: f1 = sqrt(1.0 + d1*d1 + d7*d7)
367: f2 = sqrt(1.0 + d1*d1 + d4*d4)
368: f3 = sqrt(1.0 + d3*d3 + d8*d8)
369: f4 = sqrt(1.0 + d3*d3 + d2*d2)
370: f5 = sqrt(1.0 + d2*d2 + d5*d5)
371: f6 = sqrt(1.0 + d4*d4 + d6*d6)
373: ft = ft + f2 + f4
375: df1dxc = df1dxc / f1
376: df2dxc = df2dxc / f2
377: df3dxc = df3dxc / f3
378: df4dxc = df4dxc / f4
379: df5dxc = df5dxc / f5
380: df6dxc = df6dxc / f6
382: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
383: & df5dxc + df6dxc)
384: enddo
385: enddo
387: ! Compute triangular areas along the border of the domain.
388: if (xs .eq. 0) then ! left side
389: do j=ys,ys+ym-1
390: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
391: & * rhy
392: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
393: & * rhx
394: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
395: enddo
396: endif
399: if (ys .eq. 0) then !bottom side
400: do i=xs,xs+xm-1
401: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
402: & * rhx
403: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
404: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
405: enddo
406: endif
409: if (xs + xm .eq. mx) then ! right side
410: do j=ys,ys+ym-1
411: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
412: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
413: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
414: enddo
415: endif
418: if (ys + ym .eq. my) then
419: do i=xs,xs+xm-1
420: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
421: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
422: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
423: enddo
424: endif
427: if ((ys .eq. 0) .and. (xs .eq. 0)) then
428: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
429: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
430: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
431: endif
433: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
434: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
435: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
436: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
437: endif
439: ft = ft * area
440: call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, &
441: & MPIU_SUM,MPI_COMM_WORLD,ierr)
444: ! Restore vectors
445: call VecRestoreArray(localX,x_v,x_i,ierr)
446: call VecRestoreArray(localV,g_v,g_i,ierr)
447: call VecRestoreArray(Left,left_v,left_i,ierr)
448: call VecRestoreArray(Top,top_v,top_i,ierr)
449: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
450: call VecRestoreArray(Right,right_v,right_i,ierr)
452: ! Scatter values to global vector
453: call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
454: call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)
456: call PetscLogFlops(70.0d0*xm*ym,ierr)
458: return
459: end !FormFunctionGradient
465: ! ----------------------------------------------------------------------------
466: !
467: !
468: ! FormHessian - Evaluates Hessian matrix.
469: !
470: ! Input Parameters:
471: !. tao - the Tao context
472: !. X - input vector
473: !. dummy - not used
474: !
475: ! Output Parameters:
476: !. Hessian - Hessian matrix
477: !. Hpc - optionally different preconditioning matrix
478: !. flag - flag indicating matrix structure
479: !
480: ! Notes:
481: ! Due to mesh point reordering with DMs, we must always work
482: ! with the local mesh points, and then transform them to the new
483: ! global numbering with the local-to-global mapping. We cannot work
484: ! directly with the global numbers for the original uniprocessor mesh!
485: !
486: ! MatSetValuesLocal(), using the local ordering (including
487: ! ghost points!)
488: ! - Then set matrix entries using the local ordering
489: ! by calling MatSetValuesLocal()
491: subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
492: #include "plate2f.h"
494: Tao tao
495: Vec X
496: Mat Hessian,Hpc
497: PetscInt dummy
498: PetscErrorCode ierr
500: PetscInt i,j,k,row
501: PetscInt xs,xm,gxs,gxm
502: PetscInt ys,ym,gys,gym
503: PetscInt col(0:6)
504: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
505: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
506: PetscReal d4,d5,d6,d7,d8
507: PetscReal xc,xl,xr,xt,xb,xlt,xrb
508: PetscReal hl,hr,ht,hb,hc,htl,hbr
510: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
511: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
512: ! will return an array of doubles referenced by x_array offset by x_index.
513: ! i.e., to reference the kth element of X, use x_array(k + x_index).
514: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
515: PetscReal right_v(0:1),left_v(0:1)
516: PetscReal bottom_v(0:1),top_v(0:1)
517: PetscReal x_v(0:1)
518: PetscOffset x_i,right_i,left_i
519: PetscOffset bottom_i,top_i
520: PetscReal v(0:6)
521: PetscBool assembled
522: PetscInt i1
524: i1=1
526: ! Set various matrix options
527: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, &
528: & PETSC_TRUE,ierr)
530: ! Get local mesh boundaries
531: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
532: & PETSC_NULL_INTEGER,ierr)
533: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
534: & PETSC_NULL_INTEGER,ierr)
536: ! Scatter ghost points to local vectors
537: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
538: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
540: ! Get pointers to vector data (see note on Fortran arrays above)
541: call VecGetArray(localX,x_v,x_i,ierr)
542: call VecGetArray(Top,top_v,top_i,ierr)
543: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
544: call VecGetArray(Left,left_v,left_i,ierr)
545: call VecGetArray(Right,right_v,right_i,ierr)
547: ! Initialize matrix entries to zero
548: call MatAssembled(Hessian,assembled,ierr)
549: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)
552: rhx = real(mx + 1.0)
553: rhy = real(my + 1.0)
554: hx = 1.0/rhx
555: hy = 1.0/rhy
556: hydhx = hy/hx
557: hxdhy = hx/hy
558: ! compute Hessian over the locally owned part of the mesh
560: do i=xs,xs+xm-1
561: do j=ys,ys+ym-1
562: row = (j-gys)*gxm + (i-gxs)
564: xc = x_v(row + x_i)
565: xt = xc
566: xb = xc
567: xr = xc
568: xl = xc
569: xrb = xc
570: xlt = xc
572: if (i .eq. gxs) then ! Left side
573: xl = left_v(left_i + j - ys + 1)
574: xlt = left_v(left_i + j - ys + 2)
575: else
576: xl = x_v(x_i + row -1 )
577: endif
579: if (j .eq. gys) then ! bottom side
580: xb = bottom_v(bottom_i + i - xs + 1)
581: xrb = bottom_v(bottom_i + i - xs + 2)
582: else
583: xb = x_v(x_i + row - gxm)
584: endif
586: if (i+1 .eq. gxs + gxm) then !right side
587: xr = right_v(right_i + j - ys + 1)
588: xrb = right_v(right_i + j - ys)
589: else
590: xr = x_v(x_i + row + 1)
591: endif
593: if (j+1 .eq. gym+gys) then !top side
594: xt = top_v(top_i +i - xs + 1)
595: xlt = top_v(top_i + i - xs)
596: else
597: xt = x_v(x_i + row + gxm)
598: endif
600: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
601: xlt = x_v(x_i + row - 1 + gxm)
602: endif
604: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
605: xrb = x_v(x_i + row + 1 - gxm)
606: endif
608: d1 = (xc-xl)*rhx
609: d2 = (xc-xr)*rhx
610: d3 = (xc-xt)*rhy
611: d4 = (xc-xb)*rhy
612: d5 = (xrb-xr)*rhy
613: d6 = (xrb-xb)*rhx
614: d7 = (xlt-xl)*rhy
615: d8 = (xlt-xt)*rhx
617: f1 = sqrt( 1.0 + d1*d1 + d7*d7)
618: f2 = sqrt( 1.0 + d1*d1 + d4*d4)
619: f3 = sqrt( 1.0 + d3*d3 + d8*d8)
620: f4 = sqrt( 1.0 + d3*d3 + d2*d2)
621: f5 = sqrt( 1.0 + d2*d2 + d5*d5)
622: f6 = sqrt( 1.0 + d4*d4 + d6*d6)
625: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
626: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
628: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
629: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
631: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
632: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
634: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
635: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
637: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
638: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
640: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
641: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
642: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
643: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
644: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
645: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
646: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
648: hl = hl * 0.5
649: hr = hr * 0.5
650: ht = ht * 0.5
651: hb = hb * 0.5
652: hbr = hbr * 0.5
653: htl = htl * 0.5
654: hc = hc * 0.5
656: k = 0
658: if (j .gt. 0) then
659: v(k) = hb
660: col(k) = row - gxm
661: k=k+1
662: endif
664: if ((j .gt. 0) .and. (i .lt. mx-1)) then
665: v(k) = hbr
666: col(k) = row-gxm+1
667: k=k+1
668: endif
670: if (i .gt. 0) then
671: v(k) = hl
672: col(k) = row - 1
673: k = k+1
674: endif
676: v(k) = hc
677: col(k) = row
678: k=k+1
680: if (i .lt. mx-1) then
681: v(k) = hr
682: col(k) = row + 1
683: k=k+1
684: endif
686: if ((i .gt. 0) .and. (j .lt. my-1)) then
687: v(k) = htl
688: col(k) = row + gxm - 1
689: k=k+1
690: endif
692: if (j .lt. my-1) then
693: v(k) = ht
694: col(k) = row + gxm
695: k=k+1
696: endif
698: ! Set matrix values using local numbering, defined earlier in main routine
699: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
700: & ierr)
704: enddo
705: enddo
707: ! restore vectors
708: call VecRestoreArray(localX,x_v,x_i,ierr)
709: call VecRestoreArray(Left,left_v,left_i,ierr)
710: call VecRestoreArray(Right,right_v,right_i,ierr)
711: call VecRestoreArray(Top,top_v,top_i,ierr)
712: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
715: ! Assemble the matrix
716: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
717: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)
719: call PetscLogFlops(199.0d0*xm*ym,ierr)
721: return
722: end
728: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
730: ! ----------------------------------------------------------------------------
731: !
732: !/*
733: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
734: !
735: !
736: !*/
738: subroutine MSA_BoundaryConditions(ierr)
739: #include "plate2f.h"
741: PetscErrorCode ierr
742: PetscInt i,j,k,limit,maxits
743: PetscInt xs, xm, gxs, gxm
744: PetscInt ys, ym, gys, gym
745: PetscInt bsize, lsize
746: PetscInt tsize, rsize
747: PetscReal one,two,three,tol
748: PetscReal scl,fnorm,det,xt
749: PetscReal yt,hx,hy,u1,u2,nf1,nf2
750: PetscReal njac11,njac12,njac21,njac22
751: PetscReal b, t, l, r
752: PetscReal boundary_v(0:1)
753: PetscOffset boundary_i
754: logical exitloop
755: PetscBool flg
757: limit=0
758: maxits = 5
759: tol=1e-10
760: b=-0.5
761: t= 0.5
762: l=-0.5
763: r= 0.5
764: xt=0
765: yt=0
766: one=1.0
767: two=2.0
768: three=3.0
771: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
772: & PETSC_NULL_INTEGER,ierr)
773: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
774: & gxm,gym,PETSC_NULL_INTEGER,ierr)
776: bsize = xm + 2
777: lsize = ym + 2
778: rsize = ym + 2
779: tsize = xm + 2
782: call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
783: call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
784: call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
785: call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)
787: hx= (r-l)/(mx+1)
788: hy= (t-b)/(my+1)
790: do j=0,3
792: if (j.eq.0) then
793: yt=b
794: xt=l+hx*xs
795: limit=bsize
796: call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
799: elseif (j.eq.1) then
800: yt=t
801: xt=l+hx*xs
802: limit=tsize
803: call VecGetArray(Top,boundary_v,boundary_i,ierr)
805: elseif (j.eq.2) then
806: yt=b+hy*ys
807: xt=l
808: limit=lsize
809: call VecGetArray(Left,boundary_v,boundary_i,ierr)
811: elseif (j.eq.3) then
812: yt=b+hy*ys
813: xt=r
814: limit=rsize
815: call VecGetArray(Right,boundary_v,boundary_i,ierr)
816: endif
819: do i=0,limit-1
821: u1=xt
822: u2=-yt
823: k = 0
824: exitloop = .false.
825: do while (k .lt. maxits .and. (.not. exitloop) )
827: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
828: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
829: fnorm=sqrt(nf1*nf1+nf2*nf2)
830: if (fnorm .gt. tol) then
831: njac11=one+u2*u2-u1*u1
832: njac12=two*u1*u2
833: njac21=-two*u1*u2
834: njac22=-one - u1*u1 + u2*u2
835: det = njac11*njac22-njac21*njac12
836: u1 = u1-(njac22*nf1-njac12*nf2)/det
837: u2 = u2-(njac11*nf2-njac21*nf1)/det
838: else
839: exitloop = .true.
840: endif
841: k=k+1
842: enddo
844: boundary_v(i + boundary_i) = u1*u1-u2*u2
845: if ((j .eq. 0) .or. (j .eq. 1)) then
846: xt = xt + hx
847: else
848: yt = yt + hy
849: endif
851: enddo
854: if (j.eq.0) then
855: call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
856: elseif (j.eq.1) then
857: call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
858: elseif (j.eq.2) then
859: call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
860: elseif (j.eq.3) then
861: call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
862: endif
864: enddo
867: ! Scale the boundary if desired
868: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
869: & '-bottom',scl,flg,ierr)
870: if (flg .eqv. PETSC_TRUE) then
871: call VecScale(Bottom,scl,ierr)
872: endif
874: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
875: & '-top',scl,flg,ierr)
876: if (flg .eqv. PETSC_TRUE) then
877: call VecScale(Top,scl,ierr)
878: endif
880: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
881: & '-right',scl,flg,ierr)
882: if (flg .eqv. PETSC_TRUE) then
883: call VecScale(Right,scl,ierr)
884: endif
886: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
887: & '-left',scl,flg,ierr)
888: if (flg .eqv. PETSC_TRUE) then
889: call VecScale(Left,scl,ierr)
890: endif
893: return
894: end
896: ! ----------------------------------------------------------------------------
897: !
898: !/*
899: ! MSA_Plate - Calculates an obstacle for surface to stretch over
900: !
901: ! Output Parameter:
902: !. xl - lower bound vector
903: !. xu - upper bound vector
904: !
905: !*/
907: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
908: #include "plate2f.h"
910: Tao tao
911: Vec xl,xu
912: PetscErrorCode ierr
913: PetscInt i,j,row
914: PetscInt xs, xm, ys, ym
915: PetscReal lb,ub
916: PetscInt dummy
918: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
919: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
920: ! will return an array of doubles referenced by x_array offset by x_index.
921: ! i.e., to reference the kth element of X, use x_array(k + x_index).
922: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
923: PetscReal xl_v(0:1)
924: PetscOffset xl_i
926: print *,'msa_plate'
928: lb = PETSC_NINFINITY
929: ub = PETSC_INFINITY
931: if (bmy .lt. 0) bmy = 0
932: if (bmy .gt. my) bmy = my
933: if (bmx .lt. 0) bmx = 0
934: if (bmx .gt. mx) bmx = mx
937: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
938: & PETSC_NULL_INTEGER,ierr)
940: call VecSet(xl,lb,ierr)
941: call VecSet(xu,ub,ierr)
943: call VecGetArray(xl,xl_v,xl_i,ierr)
946: do i=xs,xs+xm-1
948: do j=ys,ys+ym-1
950: row=(j-ys)*xm + (i-xs)
952: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
953: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
954: xl_v(xl_i+row) = bheight
956: endif
958: enddo
959: enddo
962: call VecRestoreArray(xl,xl_v,xl_i,ierr)
964: return
965: end
971: ! ----------------------------------------------------------------------------
972: !
973: !/*
974: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
975: !
976: ! Output Parameter:
977: !. X - vector for initial guess
978: !
979: !*/
981: subroutine MSA_InitialPoint(X, ierr)
982: #include "plate2f.h"
984: Vec X
985: PetscErrorCode ierr
986: PetscInt start,i,j
987: PetscInt row
988: PetscInt xs,xm,gxs,gxm
989: PetscInt ys,ym,gys,gym
990: PetscReal zero, np5
992: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
993: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
994: ! will return an array of doubles referenced by x_array offset by x_index.
995: ! i.e., to reference the kth element of X, use x_array(k + x_index).
996: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
997: PetscReal left_v(0:1),right_v(0:1)
998: PetscReal bottom_v(0:1),top_v(0:1)
999: PetscReal x_v(0:1)
1000: PetscOffset left_i, right_i, top_i
1001: PetscOffset bottom_i,x_i
1002: PetscBool flg
1003: PetscRandom rctx
1005: zero = 0.0
1006: np5 = -0.5
1008: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
1009: & '-start', start,flg,ierr)
1011: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
1012: call VecSet(X,zero,ierr)
1014: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
1015: call PetscRandomCreate(MPI_COMM_WORLD,rctx,ierr)
1016: do i=0,start-1
1017: call VecSetRandom(X,rctx,ierr)
1018: enddo
1020: call PetscRandomDestroy(rctx,ierr)
1021: call VecShift(X,np5,ierr)
1023: else ! average of boundary conditions
1025: ! Get Local mesh boundaries
1026: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1027: & PETSC_NULL_INTEGER,ierr)
1028: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1029: & PETSC_NULL_INTEGER,ierr)
1033: ! Get pointers to vector data
1034: call VecGetArray(Top,top_v,top_i,ierr)
1035: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1036: call VecGetArray(Left,left_v,left_i,ierr)
1037: call VecGetArray(Right,right_v,right_i,ierr)
1039: call VecGetArray(localX,x_v,x_i,ierr)
1041: ! Perform local computations
1042: do j=ys,ys+ym-1
1043: do i=xs,xs+xm-1
1044: row = (j-gys)*gxm + (i-gxs)
1045: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1046: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1047: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1048: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1049: enddo
1050: enddo
1052: ! Restore vectors
1053: call VecRestoreArray(localX,x_v,x_i,ierr)
1055: call VecRestoreArray(Left,left_v,left_i,ierr)
1056: call VecRestoreArray(Top,top_v,top_i,ierr)
1057: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1058: call VecRestoreArray(Right,right_v,right_i,ierr)
1060: call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1061: call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)
1063: endif
1065: return
1066: end