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: !
16: module mymodule
17: #include "petsc/finclude/petscdmda.h"
18: #include "petsc/finclude/petsctao.h"
19: use petscdmda
20: use petsctao
22: Vec localX, localV
23: Vec Top, Left
24: Vec Right, Bottom
25: DM dm
26: PetscReal bheight
27: PetscInt bmx, bmy
28: PetscInt mx, my, Nx, Ny, N
29: end module
31: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32: ! Variable declarations
33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
34: !
35: ! Variables:
36: ! (common from plate2f.h):
37: ! Nx, Ny number of processors in x- and y- directions
38: ! mx, my number of grid points in x,y directions
39: ! N global dimension of vector
40: use mymodule
41: implicit none
43: PetscErrorCode ierr ! used to check for functions returning nonzeros
44: Vec x ! solution vector
45: PetscInt m ! number of local elements in vector
46: Tao tao ! Tao solver context
47: Mat H ! Hessian matrix
48: ISLocalToGlobalMapping isltog ! local to global mapping object
49: PetscBool flg
50: PetscInt i1,i3,i7
52: external FormFunctionGradient
53: external FormHessian
54: external MSA_BoundaryConditions
55: external MSA_Plate
56: external MSA_InitialPoint
57: ! Initialize Tao
59: i1=1
60: i3=3
61: i7=7
63: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
64: if (ierr .ne. 0) then
65: print*,'Unable to initialize PETSc'
66: stop
67: endif
69: ! Specify default dimensions of the problem
70: mx = 10
71: my = 10
72: bheight = 0.1
74: ! Check for any command line arguments that override defaults
76: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
77: & '-mx',mx,flg,ierr)
78: CHKERRA(ierr)
79: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
80: & '-my',my,flg,ierr)
81: CHKERRA(ierr)
83: bmx = mx/2
84: bmy = my/2
86: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
87: & '-bmx',bmx,flg,ierr)
88: CHKERRA(ierr)
89: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
90: & '-bmy',bmy,flg,ierr)
91: CHKERRA(ierr)
92: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
93: & '-bheight',bheight,flg,ierr)
94: 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(PETSC_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(PETSC_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)
145: ! The Tao code begins here
146: ! Create TAO solver and set desired solution method.
147: ! This problems uses bounded variables, so the
148: ! method must either be 'tao_tron' or 'tao_blmvm'
150: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
151: CHKERRA(ierr)
152: call TaoSetType(tao,TAOBLMVM,ierr)
153: CHKERRA(ierr)
155: ! Set minimization function and gradient, hessian evaluation functions
157: call TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC, &
158: & FormFunctionGradient,0,ierr)
159: CHKERRA(ierr)
161: call TaoSetHessian(tao,H,H,FormHessian, &
162: & 0, ierr)
163: CHKERRA(ierr)
165: ! Set Variable bounds
166: call MSA_BoundaryConditions(ierr)
167: CHKERRA(ierr)
168: call TaoSetVariableBoundsRoutine(tao,MSA_Plate, &
169: & 0,ierr)
170: CHKERRA(ierr)
172: ! Set the initial solution guess
173: call MSA_InitialPoint(x, ierr)
174: CHKERRA(ierr)
175: call TaoSetSolution(tao,x,ierr)
176: CHKERRA(ierr)
178: ! Check for any tao command line options
179: call TaoSetFromOptions(tao,ierr)
180: CHKERRA(ierr)
182: ! Solve the application
183: call TaoSolve(tao,ierr)
184: CHKERRA(ierr)
186: ! Free TAO data structures
187: call TaoDestroy(tao,ierr)
188: CHKERRA(ierr)
190: ! Free PETSc data structures
191: call VecDestroy(x,ierr)
192: call VecDestroy(Top,ierr)
193: call VecDestroy(Bottom,ierr)
194: call VecDestroy(Left,ierr)
195: call VecDestroy(Right,ierr)
196: call MatDestroy(H,ierr)
197: call VecDestroy(localX,ierr)
198: call VecDestroy(localV,ierr)
199: call DMDestroy(dm,ierr)
201: ! Finalize TAO
203: call PetscFinalize(ierr)
205: end
207: ! ---------------------------------------------------------------------
208: !
209: ! FormFunctionGradient - Evaluates function f(X).
210: !
211: ! Input Parameters:
212: ! tao - the Tao context
213: ! X - the input vector
214: ! dummy - optional user-defined context, as set by TaoSetFunction()
215: ! (not used here)
216: !
217: ! Output Parameters:
218: ! fcn - the newly evaluated function
219: ! G - the gradient vector
220: ! info - error code
221: !
223: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
224: use mymodule
225: implicit none
227: ! Input/output variables
229: Tao tao
230: PetscReal fcn
231: Vec X, G
232: PetscErrorCode ierr
233: PetscInt dummy
235: PetscInt i,j,row
236: PetscInt xs, xm
237: PetscInt gxs, gxm
238: PetscInt ys, ym
239: PetscInt gys, gym
240: PetscReal ft,zero,hx,hy,hydhx,hxdhy
241: PetscReal area,rhx,rhy
242: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
243: PetscReal d4,d5,d6,d7,d8
244: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
245: PetscReal df5dxc,df6dxc
246: PetscReal xc,xl,xr,xt,xb,xlt,xrb
248: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
249: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
250: ! will return an array of doubles referenced by x_array offset by x_index.
251: ! i.e., to reference the kth element of X, use x_array(k + x_index).
252: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
253: PetscReal g_v(0:1),x_v(0:1)
254: PetscReal top_v(0:1),left_v(0:1)
255: PetscReal right_v(0:1),bottom_v(0:1)
256: PetscOffset g_i,left_i,right_i
257: PetscOffset bottom_i,top_i,x_i
259: ft = 0.0
260: zero = 0.0
261: hx = 1.0/real(mx + 1)
262: hy = 1.0/real(my + 1)
263: hydhx = hy/hx
264: hxdhy = hx/hy
265: area = 0.5 * hx * hy
266: rhx = real(mx) + 1.0
267: rhy = real(my) + 1.0
269: ! Get local mesh boundaries
270: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
271: & PETSC_NULL_INTEGER,ierr)
272: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
273: & gxm,gym,PETSC_NULL_INTEGER,ierr)
275: ! Scatter ghost points to local vector
276: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
277: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
279: ! Initialize the vector to zero
280: call VecSet(localV,zero,ierr)
282: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
283: call VecGetArray(localX,x_v,x_i,ierr)
284: call VecGetArray(localV,g_v,g_i,ierr)
285: call VecGetArray(Top,top_v,top_i,ierr)
286: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
287: call VecGetArray(Left,left_v,left_i,ierr)
288: call VecGetArray(Right,right_v,right_i,ierr)
290: ! Compute function over the locally owned part of the mesh
291: do j = ys,ys+ym-1
292: do i = xs,xs+xm-1
293: row = (j-gys)*gxm + (i-gxs)
294: xc = x_v(row+x_i)
295: xt = xc
296: xb = xc
297: xr = xc
298: xl = xc
299: xrb = xc
300: xlt = xc
302: if (i .eq. 0) then !left side
303: xl = left_v(j - ys + 1 + left_i)
304: xlt = left_v(j - ys + 2 + left_i)
305: else
306: xl = x_v(row - 1 + x_i)
307: endif
309: if (j .eq. 0) then !bottom side
310: xb = bottom_v(i - xs + 1 + bottom_i)
311: xrb = bottom_v(i - xs + 2 + bottom_i)
312: else
313: xb = x_v(row - gxm + x_i)
314: endif
316: if (i + 1 .eq. gxs + gxm) then !right side
317: xr = right_v(j - ys + 1 + right_i)
318: xrb = right_v(j - ys + right_i)
319: else
320: xr = x_v(row + 1 + x_i)
321: endif
323: if (j + 1 .eq. gys + gym) then !top side
324: xt = top_v(i - xs + 1 + top_i)
325: xlt = top_v(i - xs + top_i)
326: else
327: xt = x_v(row + gxm + x_i)
328: endif
330: if ((i .gt. gxs) .and. (j + 1 .lt. gys + gym)) then
331: xlt = x_v(row - 1 + gxm + x_i)
332: endif
334: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
335: xrb = x_v(row + 1 - gxm + x_i)
336: endif
338: d1 = xc-xl
339: d2 = xc-xr
340: d3 = xc-xt
341: d4 = xc-xb
342: d5 = xr-xrb
343: d6 = xrb-xb
344: d7 = xlt-xl
345: d8 = xt-xlt
347: df1dxc = d1 * hydhx
348: df2dxc = d1 * hydhx + d4 * hxdhy
349: df3dxc = d3 * hxdhy
350: df4dxc = d2 * hydhx + d3 * hxdhy
351: df5dxc = d2 * hydhx
352: df6dxc = d4 * hxdhy
354: d1 = d1 * rhx
355: d2 = d2 * rhx
356: d3 = d3 * rhy
357: d4 = d4 * rhy
358: d5 = d5 * rhy
359: d6 = d6 * rhx
360: d7 = d7 * rhy
361: d8 = d8 * rhx
363: f1 = sqrt(1.0 + d1*d1 + d7*d7)
364: f2 = sqrt(1.0 + d1*d1 + d4*d4)
365: f3 = sqrt(1.0 + d3*d3 + d8*d8)
366: f4 = sqrt(1.0 + d3*d3 + d2*d2)
367: f5 = sqrt(1.0 + d2*d2 + d5*d5)
368: f6 = sqrt(1.0 + d4*d4 + d6*d6)
370: ft = ft + f2 + f4
372: df1dxc = df1dxc / f1
373: df2dxc = df2dxc / f2
374: df3dxc = df3dxc / f3
375: df4dxc = df4dxc / f4
376: df5dxc = df5dxc / f5
377: df6dxc = df6dxc / f6
379: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
380: & df5dxc + df6dxc)
381: enddo
382: enddo
384: ! Compute triangular areas along the border of the domain.
385: if (xs .eq. 0) then ! left side
386: do j=ys,ys+ym-1
387: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
388: & * rhy
389: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
390: & * rhx
391: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
392: enddo
393: endif
395: if (ys .eq. 0) then !bottom side
396: do i=xs,xs+xm-1
397: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
398: & * rhx
399: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
400: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
401: enddo
402: endif
404: if (xs + xm .eq. mx) then ! right side
405: do j=ys,ys+ym-1
406: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
407: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
408: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
409: enddo
410: endif
412: if (ys + ym .eq. my) then
413: do i=xs,xs+xm-1
414: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
415: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
416: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
417: enddo
418: endif
420: if ((ys .eq. 0) .and. (xs .eq. 0)) then
421: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
422: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
423: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
424: endif
426: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
427: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
428: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
429: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
430: endif
432: ft = ft * area
433: call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, &
434: & MPIU_SUM,PETSC_COMM_WORLD,ierr)
436: ! Restore vectors
437: call VecRestoreArray(localX,x_v,x_i,ierr)
438: call VecRestoreArray(localV,g_v,g_i,ierr)
439: call VecRestoreArray(Left,left_v,left_i,ierr)
440: call VecRestoreArray(Top,top_v,top_i,ierr)
441: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
442: call VecRestoreArray(Right,right_v,right_i,ierr)
444: ! Scatter values to global vector
445: call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
446: call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)
448: call PetscLogFlops(70.0d0*xm*ym,ierr)
450: return
451: end !FormFunctionGradient
453: ! ----------------------------------------------------------------------------
454: !
455: !
456: ! FormHessian - Evaluates Hessian matrix.
457: !
458: ! Input Parameters:
459: !. tao - the Tao context
460: !. X - input vector
461: !. dummy - not used
462: !
463: ! Output Parameters:
464: !. Hessian - Hessian matrix
465: !. Hpc - optionally different preconditioning matrix
466: !. flag - flag indicating matrix structure
467: !
468: ! Notes:
469: ! Due to mesh point reordering with DMs, we must always work
470: ! with the local mesh points, and then transform them to the new
471: ! global numbering with the local-to-global mapping. We cannot work
472: ! directly with the global numbers for the original uniprocessor mesh!
473: !
474: ! MatSetValuesLocal(), using the local ordering (including
475: ! ghost points!)
476: ! - Then set matrix entries using the local ordering
477: ! by calling MatSetValuesLocal()
479: subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
480: use mymodule
481: implicit none
483: Tao tao
484: Vec X
485: Mat Hessian,Hpc
486: PetscInt dummy
487: PetscErrorCode ierr
489: PetscInt i,j,k,row
490: PetscInt xs,xm,gxs,gxm
491: PetscInt ys,ym,gys,gym
492: PetscInt col(0:6)
493: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
494: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
495: PetscReal d4,d5,d6,d7,d8
496: PetscReal xc,xl,xr,xt,xb,xlt,xrb
497: PetscReal hl,hr,ht,hb,hc,htl,hbr
499: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
500: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
501: ! will return an array of doubles referenced by x_array offset by x_index.
502: ! i.e., to reference the kth element of X, use x_array(k + x_index).
503: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
504: PetscReal right_v(0:1),left_v(0:1)
505: PetscReal bottom_v(0:1),top_v(0:1)
506: PetscReal x_v(0:1)
507: PetscOffset x_i,right_i,left_i
508: PetscOffset bottom_i,top_i
509: PetscReal v(0:6)
510: PetscBool assembled
511: PetscInt i1
513: i1=1
515: ! Set various matrix options
516: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, &
517: & PETSC_TRUE,ierr)
519: ! Get local mesh boundaries
520: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
521: & PETSC_NULL_INTEGER,ierr)
522: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
523: & PETSC_NULL_INTEGER,ierr)
525: ! Scatter ghost points to local vectors
526: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
527: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
529: ! Get pointers to vector data (see note on Fortran arrays above)
530: call VecGetArray(localX,x_v,x_i,ierr)
531: call VecGetArray(Top,top_v,top_i,ierr)
532: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
533: call VecGetArray(Left,left_v,left_i,ierr)
534: call VecGetArray(Right,right_v,right_i,ierr)
536: ! Initialize matrix entries to zero
537: call MatAssembled(Hessian,assembled,ierr)
538: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)
540: rhx = real(mx) + 1.0
541: rhy = real(my) + 1.0
542: hx = 1.0/rhx
543: hy = 1.0/rhy
544: hydhx = hy/hx
545: hxdhy = hx/hy
546: ! compute Hessian over the locally owned part of the mesh
548: do i=xs,xs+xm-1
549: do j=ys,ys+ym-1
550: row = (j-gys)*gxm + (i-gxs)
552: xc = x_v(row + x_i)
553: xt = xc
554: xb = xc
555: xr = xc
556: xl = xc
557: xrb = xc
558: xlt = xc
560: if (i .eq. gxs) then ! Left side
561: xl = left_v(left_i + j - ys + 1)
562: xlt = left_v(left_i + j - ys + 2)
563: else
564: xl = x_v(x_i + row -1)
565: endif
567: if (j .eq. gys) then ! bottom side
568: xb = bottom_v(bottom_i + i - xs + 1)
569: xrb = bottom_v(bottom_i + i - xs + 2)
570: else
571: xb = x_v(x_i + row - gxm)
572: endif
574: if (i+1 .eq. gxs + gxm) then !right side
575: xr = right_v(right_i + j - ys + 1)
576: xrb = right_v(right_i + j - ys)
577: else
578: xr = x_v(x_i + row + 1)
579: endif
581: if (j+1 .eq. gym+gys) then !top side
582: xt = top_v(top_i +i - xs + 1)
583: xlt = top_v(top_i + i - xs)
584: else
585: xt = x_v(x_i + row + gxm)
586: endif
588: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
589: xlt = x_v(x_i + row - 1 + gxm)
590: endif
592: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
593: xrb = x_v(x_i + row + 1 - gxm)
594: endif
596: d1 = (xc-xl)*rhx
597: d2 = (xc-xr)*rhx
598: d3 = (xc-xt)*rhy
599: d4 = (xc-xb)*rhy
600: d5 = (xrb-xr)*rhy
601: d6 = (xrb-xb)*rhx
602: d7 = (xlt-xl)*rhy
603: d8 = (xlt-xt)*rhx
605: f1 = sqrt(1.0 + d1*d1 + d7*d7)
606: f2 = sqrt(1.0 + d1*d1 + d4*d4)
607: f3 = sqrt(1.0 + d3*d3 + d8*d8)
608: f4 = sqrt(1.0 + d3*d3 + d2*d2)
609: f5 = sqrt(1.0 + d2*d2 + d5*d5)
610: f6 = sqrt(1.0 + d4*d4 + d6*d6)
612: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
613: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
615: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
616: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
618: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
619: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
621: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
622: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
624: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
625: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
627: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
628: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
629: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
630: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
631: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
632: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
633: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
635: hl = hl * 0.5
636: hr = hr * 0.5
637: ht = ht * 0.5
638: hb = hb * 0.5
639: hbr = hbr * 0.5
640: htl = htl * 0.5
641: hc = hc * 0.5
643: k = 0
645: if (j .gt. 0) then
646: v(k) = hb
647: col(k) = row - gxm
648: k=k+1
649: endif
651: if ((j .gt. 0) .and. (i .lt. mx-1)) then
652: v(k) = hbr
653: col(k) = row-gxm+1
654: k=k+1
655: endif
657: if (i .gt. 0) then
658: v(k) = hl
659: col(k) = row - 1
660: k = k+1
661: endif
663: v(k) = hc
664: col(k) = row
665: k=k+1
667: if (i .lt. mx-1) then
668: v(k) = hr
669: col(k) = row + 1
670: k=k+1
671: endif
673: if ((i .gt. 0) .and. (j .lt. my-1)) then
674: v(k) = htl
675: col(k) = row + gxm - 1
676: k=k+1
677: endif
679: if (j .lt. my-1) then
680: v(k) = ht
681: col(k) = row + gxm
682: k=k+1
683: endif
685: ! Set matrix values using local numbering, defined earlier in main routine
686: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
687: & ierr)
689: enddo
690: enddo
692: ! restore vectors
693: call VecRestoreArray(localX,x_v,x_i,ierr)
694: call VecRestoreArray(Left,left_v,left_i,ierr)
695: call VecRestoreArray(Right,right_v,right_i,ierr)
696: call VecRestoreArray(Top,top_v,top_i,ierr)
697: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
699: ! Assemble the matrix
700: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
701: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)
703: call PetscLogFlops(199.0d0*xm*ym,ierr)
705: return
706: end
708: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
710: ! ----------------------------------------------------------------------------
711: !
712: !/*
713: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
714: !
715: !
716: !*/
718: subroutine MSA_BoundaryConditions(ierr)
719: use mymodule
720: implicit none
722: PetscErrorCode ierr
723: PetscInt i,j,k,limit,maxits
724: PetscInt xs, xm, gxs, gxm
725: PetscInt ys, ym, gys, gym
726: PetscInt bsize, lsize
727: PetscInt tsize, rsize
728: PetscReal one,two,three,tol
729: PetscReal scl,fnorm,det,xt
730: PetscReal yt,hx,hy,u1,u2,nf1,nf2
731: PetscReal njac11,njac12,njac21,njac22
732: PetscReal b, t, l, r
733: PetscReal boundary_v(0:1)
734: PetscOffset boundary_i
735: logical exitloop
736: PetscBool flg
738: limit=0
739: maxits = 5
740: tol=1e-10
741: b=-0.5
742: t= 0.5
743: l=-0.5
744: r= 0.5
745: xt=0
746: yt=0
747: one=1.0
748: two=2.0
749: three=3.0
751: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
752: & PETSC_NULL_INTEGER,ierr)
753: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
754: & gxm,gym,PETSC_NULL_INTEGER,ierr)
756: bsize = xm + 2
757: lsize = ym + 2
758: rsize = ym + 2
759: tsize = xm + 2
761: call VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
762: call VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
763: call VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
764: call VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)
766: hx= (r-l)/(mx+1)
767: hy= (t-b)/(my+1)
769: do j=0,3
771: if (j.eq.0) then
772: yt=b
773: xt=l+hx*xs
774: limit=bsize
775: call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
777: elseif (j.eq.1) then
778: yt=t
779: xt=l+hx*xs
780: limit=tsize
781: call VecGetArray(Top,boundary_v,boundary_i,ierr)
783: elseif (j.eq.2) then
784: yt=b+hy*ys
785: xt=l
786: limit=lsize
787: call VecGetArray(Left,boundary_v,boundary_i,ierr)
789: elseif (j.eq.3) then
790: yt=b+hy*ys
791: xt=r
792: limit=rsize
793: call VecGetArray(Right,boundary_v,boundary_i,ierr)
794: endif
796: do i=0,limit-1
798: u1=xt
799: u2=-yt
800: k = 0
801: exitloop = .false.
802: do while (k .lt. maxits .and. (.not. exitloop))
804: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
805: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
806: fnorm=sqrt(nf1*nf1+nf2*nf2)
807: if (fnorm .gt. tol) then
808: njac11=one+u2*u2-u1*u1
809: njac12=two*u1*u2
810: njac21=-two*u1*u2
811: njac22=-one - u1*u1 + u2*u2
812: det = njac11*njac22-njac21*njac12
813: u1 = u1-(njac22*nf1-njac12*nf2)/det
814: u2 = u2-(njac11*nf2-njac21*nf1)/det
815: else
816: exitloop = .true.
817: endif
818: k=k+1
819: enddo
821: boundary_v(i + boundary_i) = u1*u1-u2*u2
822: if ((j .eq. 0) .or. (j .eq. 1)) then
823: xt = xt + hx
824: else
825: yt = yt + hy
826: endif
828: enddo
830: if (j.eq.0) then
831: call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
832: elseif (j.eq.1) then
833: call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
834: elseif (j.eq.2) then
835: call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
836: elseif (j.eq.3) then
837: call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
838: endif
840: enddo
842: ! Scale the boundary if desired
843: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
844: & '-bottom',scl,flg,ierr)
845: if (flg .eqv. PETSC_TRUE) then
846: call VecScale(Bottom,scl,ierr)
847: endif
849: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
850: & '-top',scl,flg,ierr)
851: if (flg .eqv. PETSC_TRUE) then
852: call VecScale(Top,scl,ierr)
853: endif
855: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
856: & '-right',scl,flg,ierr)
857: if (flg .eqv. PETSC_TRUE) then
858: call VecScale(Right,scl,ierr)
859: endif
861: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
862: & '-left',scl,flg,ierr)
863: if (flg .eqv. PETSC_TRUE) then
864: call VecScale(Left,scl,ierr)
865: endif
867: return
868: end
870: ! ----------------------------------------------------------------------------
871: !
872: !/*
873: ! MSA_Plate - Calculates an obstacle for surface to stretch over
874: !
875: ! Output Parameter:
876: !. xl - lower bound vector
877: !. xu - upper bound vector
878: !
879: !*/
881: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
882: use mymodule
883: implicit none
885: Tao tao
886: Vec xl,xu
887: PetscErrorCode ierr
888: PetscInt i,j,row
889: PetscInt xs, xm, ys, ym
890: PetscReal lb,ub
891: PetscInt dummy
893: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
894: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
895: ! will return an array of doubles referenced by x_array offset by x_index.
896: ! i.e., to reference the kth element of X, use x_array(k + x_index).
897: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
898: PetscReal xl_v(0:1)
899: PetscOffset xl_i
901: lb = PETSC_NINFINITY
902: ub = PETSC_INFINITY
904: if (bmy .lt. 0) bmy = 0
905: if (bmy .gt. my) bmy = my
906: if (bmx .lt. 0) bmx = 0
907: if (bmx .gt. mx) bmx = mx
909: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
910: & PETSC_NULL_INTEGER,ierr)
912: call VecSet(xl,lb,ierr)
913: call VecSet(xu,ub,ierr)
915: call VecGetArray(xl,xl_v,xl_i,ierr)
917: do i=xs,xs+xm-1
919: do j=ys,ys+ym-1
921: row=(j-ys)*xm + (i-xs)
923: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
924: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
925: xl_v(xl_i+row) = bheight
927: endif
929: enddo
930: enddo
932: call VecRestoreArray(xl,xl_v,xl_i,ierr)
934: return
935: end
937: ! ----------------------------------------------------------------------------
938: !
939: !/*
940: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
941: !
942: ! Output Parameter:
943: !. X - vector for initial guess
944: !
945: !*/
947: subroutine MSA_InitialPoint(X, ierr)
948: use mymodule
949: implicit none
951: Vec X
952: PetscErrorCode ierr
953: PetscInt start,i,j
954: PetscInt row
955: PetscInt xs,xm,gxs,gxm
956: PetscInt ys,ym,gys,gym
957: PetscReal zero, np5
959: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
960: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
961: ! will return an array of doubles referenced by x_array offset by x_index.
962: ! i.e., to reference the kth element of X, use x_array(k + x_index).
963: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
964: PetscReal left_v(0:1),right_v(0:1)
965: PetscReal bottom_v(0:1),top_v(0:1)
966: PetscReal x_v(0:1)
967: PetscOffset left_i, right_i, top_i
968: PetscOffset bottom_i,x_i
969: PetscBool flg
970: PetscRandom rctx
972: zero = 0.0
973: np5 = -0.5
975: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
976: & '-start', start,flg,ierr)
978: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
979: call VecSet(X,zero,ierr)
981: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
982: call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
983: do i=0,start-1
984: call VecSetRandom(X,rctx,ierr)
985: enddo
987: call PetscRandomDestroy(rctx,ierr)
988: call VecShift(X,np5,ierr)
990: else ! average of boundary conditions
992: ! Get Local mesh boundaries
993: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
994: & PETSC_NULL_INTEGER,ierr)
995: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
996: & PETSC_NULL_INTEGER,ierr)
998: ! Get pointers to vector data
999: call VecGetArray(Top,top_v,top_i,ierr)
1000: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1001: call VecGetArray(Left,left_v,left_i,ierr)
1002: call VecGetArray(Right,right_v,right_i,ierr)
1004: call VecGetArray(localX,x_v,x_i,ierr)
1006: ! Perform local computations
1007: do j=ys,ys+ym-1
1008: do i=xs,xs+xm-1
1009: row = (j-gys)*gxm + (i-gxs)
1010: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1011: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1012: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1013: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1014: enddo
1015: enddo
1017: ! Restore vectors
1018: call VecRestoreArray(localX,x_v,x_i,ierr)
1020: call VecRestoreArray(Left,left_v,left_i,ierr)
1021: call VecRestoreArray(Top,top_v,top_i,ierr)
1022: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1023: call VecRestoreArray(Right,right_v,right_i,ierr)
1025: call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1026: call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)
1028: endif
1030: return
1031: end
1033: !
1034: !/*TEST
1035: !
1036: ! build:
1037: ! requires: !complex
1038: !
1039: ! test:
1040: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1041: ! filter: sort -b
1042: ! filter_output: sort -b
1043: ! requires: !single
1044: !
1045: ! test:
1046: ! suffix: 2
1047: ! nsize: 2
1048: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1049: ! filter: sort -b
1050: ! filter_output: sort -b
1051: ! requires: !single
1052: !
1053: !TEST*/