Actual source code: plate2f.F90
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: module mymodule
29: #include "petsc/finclude/petscdmda.h"
30: #include "petsc/finclude/petsctao.h"
31: use petscdmda
32: use petsctao
34: Vec localX, localV
35: Vec Top, Left
36: Vec Right, Bottom
37: DM dm
38: PetscReal bheight
39: PetscInt bmx, bmy
40: PetscInt mx, my, Nx, Ny, N
41: end module
43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: ! Variable declarations
45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: !
47: ! Variables:
48: ! (common from plate2f.h):
49: ! Nx, Ny number of processors in x- and y- directions
50: ! mx, my number of grid points in x,y directions
51: ! N global dimension of vector
52: use mymodule
53: implicit none
55: PetscErrorCode ierr ! used to check for functions returning nonzeros
56: Vec x ! solution vector
57: PetscInt m ! number of local elements in vector
58: Tao tao ! Tao solver context
59: Mat H ! Hessian matrix
60: ISLocalToGlobalMapping isltog ! local to global mapping object
61: PetscBool flg
62: PetscInt i1,i3,i7
64: external FormFunctionGradient
65: external FormHessian
66: external MSA_BoundaryConditions
67: external MSA_Plate
68: external MSA_InitialPoint
69: ! Initialize Tao
71: i1=1
72: i3=3
73: i7=7
75: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
76: if (ierr .ne. 0) then
77: print*,'Unable to initialize PETSc'
78: stop
79: endif
81: ! Specify default dimensions of the problem
82: mx = 10
83: my = 10
84: bheight = 0.1
86: ! Check for any command line arguments that override defaults
88: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
89: & '-mx',mx,flg,ierr)
90: CHKERRA(ierr)
91: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
92: & '-my',my,flg,ierr)
93: CHKERRA(ierr)
95: bmx = mx/2
96: bmy = my/2
98: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
99: & '-bmx',bmx,flg,ierr)
100: CHKERRA(ierr)
101: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
102: & '-bmy',bmy,flg,ierr)
103: CHKERRA(ierr)
104: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
105: & '-bheight',bheight,flg,ierr)
106: CHKERRA(ierr)
108: ! Calculate any derived values from parameters
109: N = mx*my
111: ! Let Petsc determine the dimensions of the local vectors
112: Nx = PETSC_DECIDE
113: NY = PETSC_DECIDE
115: ! A two dimensional distributed array will help define this problem, which
116: ! derives from an elliptic PDE on a two-dimensional domain. From the
117: ! distributed array, create the vectors
119: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, &
120: & DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
121: & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
122: & dm,ierr)
123: CHKERRA(ierr)
124: call DMSetFromOptions(dm,ierr)
125: call DMSetUp(dm,ierr)
127: ! Extract global and local vectors from DM; The local vectors are
128: ! used solely as work space for the evaluation of the function,
129: ! gradient, and Hessian. Duplicate for remaining vectors that are
130: ! the same types.
132: call DMCreateGlobalVector(dm,x,ierr)
133: CHKERRA(ierr)
134: call DMCreateLocalVector(dm,localX,ierr)
135: CHKERRA(ierr)
136: call VecDuplicate(localX,localV,ierr)
137: CHKERRA(ierr)
139: ! Create a matrix data structure to store the Hessian.
140: ! Here we (optionally) also associate the local numbering scheme
141: ! with the matrix so that later we can use local indices for matrix
142: ! assembly
144: call VecGetLocalSize(x,m,ierr)
145: CHKERRA(ierr)
146: call MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, &
147: & i3,PETSC_NULL_INTEGER,H,ierr)
148: CHKERRA(ierr)
150: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
151: CHKERRA(ierr)
152: call DMGetLocalToGlobalMapping(dm,isltog,ierr)
153: CHKERRA(ierr)
154: call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
155: CHKERRA(ierr)
157: ! The Tao code begins here
158: ! Create TAO solver and set desired solution method.
159: ! This problems uses bounded variables, so the
160: ! method must either be 'tao_tron' or 'tao_blmvm'
162: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
163: CHKERRA(ierr)
164: call TaoSetType(tao,TAOBLMVM,ierr)
165: CHKERRA(ierr)
167: ! Set minimization function and gradient, hessian evaluation functions
169: call TaoSetObjectiveAndGradientRoutine(tao, &
170: & FormFunctionGradient,0,ierr)
171: CHKERRA(ierr)
173: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
174: & 0, ierr)
175: CHKERRA(ierr)
177: ! Set Variable bounds
178: call MSA_BoundaryConditions(ierr)
179: CHKERRA(ierr)
180: call TaoSetVariableBoundsRoutine(tao,MSA_Plate, &
181: & 0,ierr)
182: CHKERRA(ierr)
184: ! Set the initial solution guess
185: call MSA_InitialPoint(x, ierr)
186: CHKERRA(ierr)
187: call TaoSetInitialVector(tao,x,ierr)
188: CHKERRA(ierr)
190: ! Check for any tao command line options
191: call TaoSetFromOptions(tao,ierr)
192: CHKERRA(ierr)
194: ! Solve the application
195: call TaoSolve(tao,ierr)
196: CHKERRA(ierr)
198: ! Free TAO data structures
199: call TaoDestroy(tao,ierr)
200: CHKERRA(ierr)
202: ! Free PETSc data structures
203: call VecDestroy(x,ierr)
204: call VecDestroy(Top,ierr)
205: call VecDestroy(Bottom,ierr)
206: call VecDestroy(Left,ierr)
207: call VecDestroy(Right,ierr)
208: call MatDestroy(H,ierr)
209: call VecDestroy(localX,ierr)
210: call VecDestroy(localV,ierr)
211: call DMDestroy(dm,ierr)
213: ! Finalize TAO
215: call PetscFinalize(ierr)
217: end
219: ! ---------------------------------------------------------------------
220: !
221: ! FormFunctionGradient - Evaluates function f(X).
222: !
223: ! Input Parameters:
224: ! tao - the Tao context
225: ! X - the input vector
226: ! dummy - optional user-defined context, as set by TaoSetFunction()
227: ! (not used here)
228: !
229: ! Output Parameters:
230: ! fcn - the newly evaluated function
231: ! G - the gradient vector
232: ! info - error code
233: !
235: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
236: use mymodule
237: implicit none
239: ! Input/output variables
241: Tao tao
242: PetscReal fcn
243: Vec X, G
244: PetscErrorCode ierr
245: PetscInt dummy
247: PetscInt i,j,row
248: PetscInt xs, xm
249: PetscInt gxs, gxm
250: PetscInt ys, ym
251: PetscInt gys, gym
252: PetscReal ft,zero,hx,hy,hydhx,hxdhy
253: PetscReal area,rhx,rhy
254: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
255: PetscReal d4,d5,d6,d7,d8
256: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
257: PetscReal df5dxc,df6dxc
258: PetscReal xc,xl,xr,xt,xb,xlt,xrb
260: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
261: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
262: ! will return an array of doubles referenced by x_array offset by x_index.
263: ! i.e., to reference the kth element of X, use x_array(k + x_index).
264: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
265: PetscReal g_v(0:1),x_v(0:1)
266: PetscReal top_v(0:1),left_v(0:1)
267: PetscReal right_v(0:1),bottom_v(0:1)
268: PetscOffset g_i,left_i,right_i
269: PetscOffset bottom_i,top_i,x_i
271: ft = 0.0
272: zero = 0.0
273: hx = 1.0/real(mx + 1)
274: hy = 1.0/real(my + 1)
275: hydhx = hy/hx
276: hxdhy = hx/hy
277: area = 0.5 * hx * hy
278: rhx = real(mx) + 1.0
279: rhy = real(my) + 1.0
281: ! Get local mesh boundaries
282: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
283: & PETSC_NULL_INTEGER,ierr)
284: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
285: & gxm,gym,PETSC_NULL_INTEGER,ierr)
287: ! Scatter ghost points to local vector
288: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
289: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
291: ! Initialize the vector to zero
292: call VecSet(localV,zero,ierr)
294: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
295: call VecGetArray(localX,x_v,x_i,ierr)
296: call VecGetArray(localV,g_v,g_i,ierr)
297: call VecGetArray(Top,top_v,top_i,ierr)
298: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
299: call VecGetArray(Left,left_v,left_i,ierr)
300: call VecGetArray(Right,right_v,right_i,ierr)
302: ! Compute function over the locally owned part of the mesh
303: do j = ys,ys+ym-1
304: do i = xs,xs+xm-1
305: row = (j-gys)*gxm + (i-gxs)
306: xc = x_v(row+x_i)
307: xt = xc
308: xb = xc
309: xr = xc
310: xl = xc
311: xrb = xc
312: xlt = xc
314: if (i .eq. 0) then !left side
315: xl = left_v(j - ys + 1 + left_i)
316: xlt = left_v(j - ys + 2 + left_i)
317: else
318: xl = x_v(row - 1 + x_i)
319: endif
321: if (j .eq. 0) then !bottom side
322: xb = bottom_v(i - xs + 1 + bottom_i)
323: xrb = bottom_v(i - xs + 2 + bottom_i)
324: else
325: xb = x_v(row - gxm + x_i)
326: endif
328: if (i + 1 .eq. gxs + gxm) then !right side
329: xr = right_v(j - ys + 1 + right_i)
330: xrb = right_v(j - ys + right_i)
331: else
332: xr = x_v(row + 1 + x_i)
333: endif
335: if (j + 1 .eq. gys + gym) then !top side
336: xt = top_v(i - xs + 1 + top_i)
337: xlt = top_v(i - xs + top_i)
338: else
339: xt = x_v(row + gxm + x_i)
340: endif
342: if ((i .gt. gxs) .and. (j + 1 .lt. gys + gym)) then
343: xlt = x_v(row - 1 + gxm + x_i)
344: endif
346: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
347: xrb = x_v(row + 1 - gxm + x_i)
348: endif
350: d1 = xc-xl
351: d2 = xc-xr
352: d3 = xc-xt
353: d4 = xc-xb
354: d5 = xr-xrb
355: d6 = xrb-xb
356: d7 = xlt-xl
357: d8 = xt-xlt
359: df1dxc = d1 * hydhx
360: df2dxc = d1 * hydhx + d4 * hxdhy
361: df3dxc = d3 * hxdhy
362: df4dxc = d2 * hydhx + d3 * hxdhy
363: df5dxc = d2 * hydhx
364: df6dxc = d4 * hxdhy
366: d1 = d1 * rhx
367: d2 = d2 * rhx
368: d3 = d3 * rhy
369: d4 = d4 * rhy
370: d5 = d5 * rhy
371: d6 = d6 * rhx
372: d7 = d7 * rhy
373: d8 = d8 * rhx
375: f1 = sqrt(1.0 + d1*d1 + d7*d7)
376: f2 = sqrt(1.0 + d1*d1 + d4*d4)
377: f3 = sqrt(1.0 + d3*d3 + d8*d8)
378: f4 = sqrt(1.0 + d3*d3 + d2*d2)
379: f5 = sqrt(1.0 + d2*d2 + d5*d5)
380: f6 = sqrt(1.0 + d4*d4 + d6*d6)
382: ft = ft + f2 + f4
384: df1dxc = df1dxc / f1
385: df2dxc = df2dxc / f2
386: df3dxc = df3dxc / f3
387: df4dxc = df4dxc / f4
388: df5dxc = df5dxc / f5
389: df6dxc = df6dxc / f6
391: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
392: & df5dxc + df6dxc)
393: enddo
394: enddo
396: ! Compute triangular areas along the border of the domain.
397: if (xs .eq. 0) then ! left side
398: do j=ys,ys+ym-1
399: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
400: & * rhy
401: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
402: & * rhx
403: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
404: enddo
405: endif
407: if (ys .eq. 0) then !bottom side
408: do i=xs,xs+xm-1
409: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
410: & * rhx
411: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
412: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
413: enddo
414: endif
416: if (xs + xm .eq. mx) then ! right side
417: do j=ys,ys+ym-1
418: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
419: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
420: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
421: enddo
422: endif
424: if (ys + ym .eq. my) then
425: do i=xs,xs+xm-1
426: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
427: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
428: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
429: enddo
430: endif
432: if ((ys .eq. 0) .and. (xs .eq. 0)) then
433: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
434: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
435: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
436: endif
438: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
439: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
440: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
441: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
442: endif
444: ft = ft * area
445: call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, &
446: & MPIU_SUM,PETSC_COMM_WORLD,ierr)
448: ! Restore vectors
449: call VecRestoreArray(localX,x_v,x_i,ierr)
450: call VecRestoreArray(localV,g_v,g_i,ierr)
451: call VecRestoreArray(Left,left_v,left_i,ierr)
452: call VecRestoreArray(Top,top_v,top_i,ierr)
453: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
454: call VecRestoreArray(Right,right_v,right_i,ierr)
456: ! Scatter values to global vector
457: call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
458: call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)
460: call PetscLogFlops(70.0d0*xm*ym,ierr)
462: return
463: 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: use mymodule
493: implicit none
495: Tao tao
496: Vec X
497: Mat Hessian,Hpc
498: PetscInt dummy
499: PetscErrorCode ierr
501: PetscInt i,j,k,row
502: PetscInt xs,xm,gxs,gxm
503: PetscInt ys,ym,gys,gym
504: PetscInt col(0:6)
505: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
506: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
507: PetscReal d4,d5,d6,d7,d8
508: PetscReal xc,xl,xr,xt,xb,xlt,xrb
509: PetscReal hl,hr,ht,hb,hc,htl,hbr
511: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
512: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
513: ! will return an array of doubles referenced by x_array offset by x_index.
514: ! i.e., to reference the kth element of X, use x_array(k + x_index).
515: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
516: PetscReal right_v(0:1),left_v(0:1)
517: PetscReal bottom_v(0:1),top_v(0:1)
518: PetscReal x_v(0:1)
519: PetscOffset x_i,right_i,left_i
520: PetscOffset bottom_i,top_i
521: PetscReal v(0:6)
522: PetscBool assembled
523: PetscInt i1
525: i1=1
527: ! Set various matrix options
528: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, &
529: & PETSC_TRUE,ierr)
531: ! Get local mesh boundaries
532: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
533: & PETSC_NULL_INTEGER,ierr)
534: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
535: & PETSC_NULL_INTEGER,ierr)
537: ! Scatter ghost points to local vectors
538: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
539: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
541: ! Get pointers to vector data (see note on Fortran arrays above)
542: call VecGetArray(localX,x_v,x_i,ierr)
543: call VecGetArray(Top,top_v,top_i,ierr)
544: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
545: call VecGetArray(Left,left_v,left_i,ierr)
546: call VecGetArray(Right,right_v,right_i,ierr)
548: ! Initialize matrix entries to zero
549: call MatAssembled(Hessian,assembled,ierr)
550: 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)
624: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
625: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
627: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
628: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
630: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
631: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
633: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
634: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
636: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
637: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
639: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
640: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
641: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
642: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
643: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
644: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
645: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
647: hl = hl * 0.5
648: hr = hr * 0.5
649: ht = ht * 0.5
650: hb = hb * 0.5
651: hbr = hbr * 0.5
652: htl = htl * 0.5
653: hc = hc * 0.5
655: k = 0
657: if (j .gt. 0) then
658: v(k) = hb
659: col(k) = row - gxm
660: k=k+1
661: endif
663: if ((j .gt. 0) .and. (i .lt. mx-1)) then
664: v(k) = hbr
665: col(k) = row-gxm+1
666: k=k+1
667: endif
669: if (i .gt. 0) then
670: v(k) = hl
671: col(k) = row - 1
672: k = k+1
673: endif
675: v(k) = hc
676: col(k) = row
677: k=k+1
679: if (i .lt. mx-1) then
680: v(k) = hr
681: col(k) = row + 1
682: k=k+1
683: endif
685: if ((i .gt. 0) .and. (j .lt. my-1)) then
686: v(k) = htl
687: col(k) = row + gxm - 1
688: k=k+1
689: endif
691: if (j .lt. my-1) then
692: v(k) = ht
693: col(k) = row + gxm
694: k=k+1
695: endif
697: ! Set matrix values using local numbering, defined earlier in main routine
698: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
699: & ierr)
701: enddo
702: enddo
704: ! restore vectors
705: call VecRestoreArray(localX,x_v,x_i,ierr)
706: call VecRestoreArray(Left,left_v,left_i,ierr)
707: call VecRestoreArray(Right,right_v,right_i,ierr)
708: call VecRestoreArray(Top,top_v,top_i,ierr)
709: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
711: ! Assemble the matrix
712: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
713: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)
715: call PetscLogFlops(199.0d0*xm*ym,ierr)
717: return
718: end
720: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
722: ! ----------------------------------------------------------------------------
723: !
724: !/*
725: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
726: !
727: !
728: !*/
730: subroutine MSA_BoundaryConditions(ierr)
731: use mymodule
732: implicit none
734: PetscErrorCode ierr
735: PetscInt i,j,k,limit,maxits
736: PetscInt xs, xm, gxs, gxm
737: PetscInt ys, ym, gys, gym
738: PetscInt bsize, lsize
739: PetscInt tsize, rsize
740: PetscReal one,two,three,tol
741: PetscReal scl,fnorm,det,xt
742: PetscReal yt,hx,hy,u1,u2,nf1,nf2
743: PetscReal njac11,njac12,njac21,njac22
744: PetscReal b, t, l, r
745: PetscReal boundary_v(0:1)
746: PetscOffset boundary_i
747: logical exitloop
748: PetscBool flg
750: limit=0
751: maxits = 5
752: tol=1e-10
753: b=-0.5
754: t= 0.5
755: l=-0.5
756: r= 0.5
757: xt=0
758: yt=0
759: one=1.0
760: two=2.0
761: three=3.0
763: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
764: & PETSC_NULL_INTEGER,ierr)
765: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
766: & gxm,gym,PETSC_NULL_INTEGER,ierr)
768: bsize = xm + 2
769: lsize = ym + 2
770: rsize = ym + 2
771: tsize = xm + 2
773: call VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
774: call VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
775: call VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
776: call VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)
778: hx= (r-l)/(mx+1)
779: hy= (t-b)/(my+1)
781: do j=0,3
783: if (j.eq.0) then
784: yt=b
785: xt=l+hx*xs
786: limit=bsize
787: call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
789: elseif (j.eq.1) then
790: yt=t
791: xt=l+hx*xs
792: limit=tsize
793: call VecGetArray(Top,boundary_v,boundary_i,ierr)
795: elseif (j.eq.2) then
796: yt=b+hy*ys
797: xt=l
798: limit=lsize
799: call VecGetArray(Left,boundary_v,boundary_i,ierr)
801: elseif (j.eq.3) then
802: yt=b+hy*ys
803: xt=r
804: limit=rsize
805: call VecGetArray(Right,boundary_v,boundary_i,ierr)
806: endif
808: do i=0,limit-1
810: u1=xt
811: u2=-yt
812: k = 0
813: exitloop = .false.
814: do while (k .lt. maxits .and. (.not. exitloop))
816: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
817: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
818: fnorm=sqrt(nf1*nf1+nf2*nf2)
819: if (fnorm .gt. tol) then
820: njac11=one+u2*u2-u1*u1
821: njac12=two*u1*u2
822: njac21=-two*u1*u2
823: njac22=-one - u1*u1 + u2*u2
824: det = njac11*njac22-njac21*njac12
825: u1 = u1-(njac22*nf1-njac12*nf2)/det
826: u2 = u2-(njac11*nf2-njac21*nf1)/det
827: else
828: exitloop = .true.
829: endif
830: k=k+1
831: enddo
833: boundary_v(i + boundary_i) = u1*u1-u2*u2
834: if ((j .eq. 0) .or. (j .eq. 1)) then
835: xt = xt + hx
836: else
837: yt = yt + hy
838: endif
840: enddo
842: if (j.eq.0) then
843: call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
844: elseif (j.eq.1) then
845: call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
846: elseif (j.eq.2) then
847: call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
848: elseif (j.eq.3) then
849: call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
850: endif
852: enddo
854: ! Scale the boundary if desired
855: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
856: & '-bottom',scl,flg,ierr)
857: if (flg .eqv. PETSC_TRUE) then
858: call VecScale(Bottom,scl,ierr)
859: endif
861: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
862: & '-top',scl,flg,ierr)
863: if (flg .eqv. PETSC_TRUE) then
864: call VecScale(Top,scl,ierr)
865: endif
867: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
868: & '-right',scl,flg,ierr)
869: if (flg .eqv. PETSC_TRUE) then
870: call VecScale(Right,scl,ierr)
871: endif
873: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
874: & '-left',scl,flg,ierr)
875: if (flg .eqv. PETSC_TRUE) then
876: call VecScale(Left,scl,ierr)
877: endif
879: return
880: end
882: ! ----------------------------------------------------------------------------
883: !
884: !/*
885: ! MSA_Plate - Calculates an obstacle for surface to stretch over
886: !
887: ! Output Parameter:
888: !. xl - lower bound vector
889: !. xu - upper bound vector
890: !
891: !*/
893: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
894: use mymodule
895: implicit none
897: Tao tao
898: Vec xl,xu
899: PetscErrorCode ierr
900: PetscInt i,j,row
901: PetscInt xs, xm, ys, ym
902: PetscReal lb,ub
903: PetscInt dummy
905: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
906: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
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: PetscReal xl_v(0:1)
911: PetscOffset xl_i
913: lb = PETSC_NINFINITY
914: ub = PETSC_INFINITY
916: if (bmy .lt. 0) bmy = 0
917: if (bmy .gt. my) bmy = my
918: if (bmx .lt. 0) bmx = 0
919: if (bmx .gt. mx) bmx = mx
921: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
922: & PETSC_NULL_INTEGER,ierr)
924: call VecSet(xl,lb,ierr)
925: call VecSet(xu,ub,ierr)
927: call VecGetArray(xl,xl_v,xl_i,ierr)
929: do i=xs,xs+xm-1
931: do j=ys,ys+ym-1
933: row=(j-ys)*xm + (i-xs)
935: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
936: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
937: xl_v(xl_i+row) = bheight
939: endif
941: enddo
942: enddo
944: call VecRestoreArray(xl,xl_v,xl_i,ierr)
946: return
947: end
949: ! ----------------------------------------------------------------------------
950: !
951: !/*
952: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
953: !
954: ! Output Parameter:
955: !. X - vector for initial guess
956: !
957: !*/
959: subroutine MSA_InitialPoint(X, ierr)
960: use mymodule
961: implicit none
963: Vec X
964: PetscErrorCode ierr
965: PetscInt start,i,j
966: PetscInt row
967: PetscInt xs,xm,gxs,gxm
968: PetscInt ys,ym,gys,gym
969: PetscReal zero, np5
971: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
972: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
973: ! will return an array of doubles referenced by x_array offset by x_index.
974: ! i.e., to reference the kth element of X, use x_array(k + x_index).
975: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
976: PetscReal left_v(0:1),right_v(0:1)
977: PetscReal bottom_v(0:1),top_v(0:1)
978: PetscReal x_v(0:1)
979: PetscOffset left_i, right_i, top_i
980: PetscOffset bottom_i,x_i
981: PetscBool flg
982: PetscRandom rctx
984: zero = 0.0
985: np5 = -0.5
987: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
988: & '-start', start,flg,ierr)
990: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
991: call VecSet(X,zero,ierr)
993: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
994: call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
995: do i=0,start-1
996: call VecSetRandom(X,rctx,ierr)
997: enddo
999: call PetscRandomDestroy(rctx,ierr)
1000: call VecShift(X,np5,ierr)
1002: else ! average of boundary conditions
1004: ! Get Local mesh boundaries
1005: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1006: & PETSC_NULL_INTEGER,ierr)
1007: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1008: & PETSC_NULL_INTEGER,ierr)
1010: ! Get pointers to vector data
1011: call VecGetArray(Top,top_v,top_i,ierr)
1012: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1013: call VecGetArray(Left,left_v,left_i,ierr)
1014: call VecGetArray(Right,right_v,right_i,ierr)
1016: call VecGetArray(localX,x_v,x_i,ierr)
1018: ! Perform local computations
1019: do j=ys,ys+ym-1
1020: do i=xs,xs+xm-1
1021: row = (j-gys)*gxm + (i-gxs)
1022: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1023: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1024: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1025: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1026: enddo
1027: enddo
1029: ! Restore vectors
1030: call VecRestoreArray(localX,x_v,x_i,ierr)
1032: call VecRestoreArray(Left,left_v,left_i,ierr)
1033: call VecRestoreArray(Top,top_v,top_i,ierr)
1034: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1035: call VecRestoreArray(Right,right_v,right_i,ierr)
1037: call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1038: call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)
1040: endif
1042: return
1043: end
1045: !
1046: !/*TEST
1047: !
1048: ! build:
1049: ! requires: !complex
1050: !
1051: ! test:
1052: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1053: ! filter: sort -b
1054: ! filter_output: sort -b
1055: ! requires: !single
1056: !
1057: ! test:
1058: ! suffix: 2
1059: ! nsize: 2
1060: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1061: ! filter: sort -b
1062: ! filter_output: sort -b
1063: ! requires: !single
1064: !
1065: !TEST*/