Actual source code: plate2f.F90
petsc-3.13.6 2020-09-29
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
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: ! Variable declarations
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: !
48: ! Variables:
49: ! (common from plate2f.h):
50: ! Nx, Ny number of processors in x- and y- directions
51: ! mx, my number of grid points in x,y directions
52: ! N global dimension of vector
53: use mymodule
54: implicit none
56: PetscErrorCode ierr ! used to check for functions returning nonzeros
57: Vec x ! solution vector
58: PetscInt m ! number of local elements in vector
59: Tao tao ! Tao solver context
60: Mat H ! Hessian matrix
61: ISLocalToGlobalMapping isltog ! local to global mapping object
62: PetscBool flg
63: PetscInt i1,i3,i7
66: external FormFunctionGradient
67: external FormHessian
68: external MSA_BoundaryConditions
69: external MSA_Plate
70: external MSA_InitialPoint
71: ! Initialize Tao
73: i1=1
74: i3=3
75: i7=7
78: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
79: if (ierr .ne. 0) then
80: print*,'Unable to initialize PETSc'
81: stop
82: endif
84: ! Specify default dimensions of the problem
85: mx = 10
86: my = 10
87: bheight = 0.1
89: ! Check for any command line arguments that override defaults
91: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
92: & '-mx',mx,flg,ierr)
93: CHKERRA(ierr)
94: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
95: & '-my',my,flg,ierr)
96: CHKERRA(ierr)
98: bmx = mx/2
99: bmy = my/2
101: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
102: & '-bmx',bmx,flg,ierr)
103: CHKERRA(ierr)
104: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
105: & '-bmy',bmy,flg,ierr)
106: CHKERRA(ierr)
107: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
108: & '-bheight',bheight,flg,ierr)
109: CHKERRA(ierr)
112: ! Calculate any derived values from parameters
113: N = mx*my
115: ! Let Petsc determine the dimensions of the local vectors
116: Nx = PETSC_DECIDE
117: NY = PETSC_DECIDE
119: ! A two dimensional distributed array will help define this problem, which
120: ! derives from an elliptic PDE on a two-dimensional domain. From the
121: ! distributed array, create the vectors
123: call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, &
124: & DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, &
125: & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
126: & dm,ierr)
127: CHKERRA(ierr)
128: call DMSetFromOptions(dm,ierr)
129: call DMSetUp(dm,ierr)
131: ! Extract global and local vectors from DM; The local vectors are
132: ! used solely as work space for the evaluation of the function,
133: ! gradient, and Hessian. Duplicate for remaining vectors that are
134: ! the same types.
136: call DMCreateGlobalVector(dm,x,ierr)
137: CHKERRA(ierr)
138: call DMCreateLocalVector(dm,localX,ierr)
139: CHKERRA(ierr)
140: call VecDuplicate(localX,localV,ierr)
141: CHKERRA(ierr)
143: ! Create a matrix data structure to store the Hessian.
144: ! Here we (optionally) also associate the local numbering scheme
145: ! with the matrix so that later we can use local indices for matrix
146: ! assembly
148: call VecGetLocalSize(x,m,ierr)
149: CHKERRA(ierr)
150: call MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, &
151: & i3,PETSC_NULL_INTEGER,H,ierr)
152: CHKERRA(ierr)
154: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
155: CHKERRA(ierr)
156: call DMGetLocalToGlobalMapping(dm,isltog,ierr)
157: CHKERRA(ierr)
158: call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr)
159: CHKERRA(ierr)
162: ! The Tao code begins here
163: ! Create TAO solver and set desired solution method.
164: ! This problems uses bounded variables, so the
165: ! method must either be 'tao_tron' or 'tao_blmvm'
167: call TaoCreate(PETSC_COMM_WORLD,tao,ierr)
168: CHKERRA(ierr)
169: call TaoSetType(tao,TAOBLMVM,ierr)
170: CHKERRA(ierr)
172: ! Set minimization function and gradient, hessian evaluation functions
174: call TaoSetObjectiveAndGradientRoutine(tao, &
175: & FormFunctionGradient,0,ierr)
176: CHKERRA(ierr)
178: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
179: & 0, ierr)
180: CHKERRA(ierr)
182: ! Set Variable bounds
183: call MSA_BoundaryConditions(ierr)
184: CHKERRA(ierr)
185: call TaoSetVariableBoundsRoutine(tao,MSA_Plate, &
186: & 0,ierr)
187: CHKERRA(ierr)
189: ! Set the initial solution guess
190: call MSA_InitialPoint(x, ierr)
191: CHKERRA(ierr)
192: call TaoSetInitialVector(tao,x,ierr)
193: CHKERRA(ierr)
195: ! Check for any tao command line options
196: call TaoSetFromOptions(tao,ierr)
197: CHKERRA(ierr)
199: ! Solve the Section 1.5 Writing Application Codes with PETSc
200: call TaoSolve(tao,ierr)
201: CHKERRA(ierr)
203: ! Free TAO data structures
204: call TaoDestroy(tao,ierr)
205: CHKERRA(ierr)
207: ! Free PETSc data structures
208: call VecDestroy(x,ierr)
209: call VecDestroy(Top,ierr)
210: call VecDestroy(Bottom,ierr)
211: call VecDestroy(Left,ierr)
212: call VecDestroy(Right,ierr)
213: call MatDestroy(H,ierr)
214: call VecDestroy(localX,ierr)
215: call VecDestroy(localV,ierr)
216: call DMDestroy(dm,ierr)
218: ! Finalize TAO
220: call PetscFinalize(ierr)
222: end
224: ! ---------------------------------------------------------------------
225: !
226: ! FormFunctionGradient - Evaluates function f(X).
227: !
228: ! Input Parameters:
229: ! tao - the Tao context
230: ! X - the input vector
231: ! dummy - optional user-defined context, as set by TaoSetFunction()
232: ! (not used here)
233: !
234: ! Output Parameters:
235: ! fcn - the newly evaluated function
236: ! G - the gradient vector
237: ! info - error code
238: !
241: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
242: use mymodule
243: implicit none
244:
245: ! Input/output variables
247: Tao tao
248: PetscReal fcn
249: Vec X, G
250: PetscErrorCode ierr
251: PetscInt dummy
253: PetscInt i,j,row
254: PetscInt xs, xm
255: PetscInt gxs, gxm
256: PetscInt ys, ym
257: PetscInt gys, gym
258: PetscReal ft,zero,hx,hy,hydhx,hxdhy
259: PetscReal area,rhx,rhy
260: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
261: PetscReal d4,d5,d6,d7,d8
262: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
263: PetscReal df5dxc,df6dxc
264: PetscReal xc,xl,xr,xt,xb,xlt,xrb
267: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
268: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
269: ! will return an array of doubles referenced by x_array offset by x_index.
270: ! i.e., to reference the kth element of X, use x_array(k + x_index).
271: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
272: PetscReal g_v(0:1),x_v(0:1)
273: PetscReal top_v(0:1),left_v(0:1)
274: PetscReal right_v(0:1),bottom_v(0:1)
275: PetscOffset g_i,left_i,right_i
276: PetscOffset bottom_i,top_i,x_i
278: ft = 0.0
279: zero = 0.0
280: hx = 1.0/real(mx + 1)
281: hy = 1.0/real(my + 1)
282: hydhx = hy/hx
283: hxdhy = hx/hy
284: area = 0.5 * hx * hy
285: rhx = real(mx) + 1.0
286: rhy = real(my) + 1.0
289: ! Get local mesh boundaries
290: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
291: & PETSC_NULL_INTEGER,ierr)
292: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
293: & gxm,gym,PETSC_NULL_INTEGER,ierr)
295: ! Scatter ghost points to local vector
296: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
297: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
299: ! Initialize the vector to zero
300: call VecSet(localV,zero,ierr)
302: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
303: call VecGetArray(localX,x_v,x_i,ierr)
304: call VecGetArray(localV,g_v,g_i,ierr)
305: call VecGetArray(Top,top_v,top_i,ierr)
306: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
307: call VecGetArray(Left,left_v,left_i,ierr)
308: call VecGetArray(Right,right_v,right_i,ierr)
310: ! Compute function over the locally owned part of the mesh
311: do j = ys,ys+ym-1
312: do i = xs,xs+xm-1
313: row = (j-gys)*gxm + (i-gxs)
314: xc = x_v(row+x_i)
315: xt = xc
316: xb = xc
317: xr = xc
318: xl = xc
319: xrb = xc
320: xlt = xc
322: if (i .eq. 0) then !left side
323: xl = left_v(j - ys + 1 + left_i)
324: xlt = left_v(j - ys + 2 + left_i)
325: else
326: xl = x_v(row - 1 + x_i)
327: endif
329: if (j .eq. 0) then !bottom side
330: xb = bottom_v(i - xs + 1 + bottom_i)
331: xrb = bottom_v(i - xs + 2 + bottom_i)
332: else
333: xb = x_v(row - gxm + x_i)
334: endif
336: if (i + 1 .eq. gxs + gxm) then !right side
337: xr = right_v(j - ys + 1 + right_i)
338: xrb = right_v(j - ys + right_i)
339: else
340: xr = x_v(row + 1 + x_i)
341: endif
343: if (j + 1 .eq. gys + gym) then !top side
344: xt = top_v(i - xs + 1 + top_i)
345: xlt = top_v(i - xs + top_i)
346: else
347: xt = x_v(row + gxm + x_i)
348: endif
350: if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
351: xlt = x_v(row - 1 + gxm + x_i)
352: endif
354: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
355: xrb = x_v(row + 1 - gxm + x_i)
356: endif
358: d1 = xc-xl
359: d2 = xc-xr
360: d3 = xc-xt
361: d4 = xc-xb
362: d5 = xr-xrb
363: d6 = xrb-xb
364: d7 = xlt-xl
365: d8 = xt-xlt
367: df1dxc = d1 * hydhx
368: df2dxc = d1 * hydhx + d4 * hxdhy
369: df3dxc = d3 * hxdhy
370: df4dxc = d2 * hydhx + d3 * hxdhy
371: df5dxc = d2 * hydhx
372: df6dxc = d4 * hxdhy
374: d1 = d1 * rhx
375: d2 = d2 * rhx
376: d3 = d3 * rhy
377: d4 = d4 * rhy
378: d5 = d5 * rhy
379: d6 = d6 * rhx
380: d7 = d7 * rhy
381: d8 = d8 * rhx
383: f1 = sqrt(1.0 + d1*d1 + d7*d7)
384: f2 = sqrt(1.0 + d1*d1 + d4*d4)
385: f3 = sqrt(1.0 + d3*d3 + d8*d8)
386: f4 = sqrt(1.0 + d3*d3 + d2*d2)
387: f5 = sqrt(1.0 + d2*d2 + d5*d5)
388: f6 = sqrt(1.0 + d4*d4 + d6*d6)
390: ft = ft + f2 + f4
392: df1dxc = df1dxc / f1
393: df2dxc = df2dxc / f2
394: df3dxc = df3dxc / f3
395: df4dxc = df4dxc / f4
396: df5dxc = df5dxc / f5
397: df6dxc = df6dxc / f6
399: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
400: & df5dxc + df6dxc)
401: enddo
402: enddo
404: ! Compute triangular areas along the border of the domain.
405: if (xs .eq. 0) then ! left side
406: do j=ys,ys+ym-1
407: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
408: & * rhy
409: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
410: & * rhx
411: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
412: enddo
413: endif
416: if (ys .eq. 0) then !bottom side
417: do i=xs,xs+xm-1
418: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
419: & * rhx
420: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
421: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
422: enddo
423: endif
426: if (xs + xm .eq. mx) then ! right side
427: do j=ys,ys+ym-1
428: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
429: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
430: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
431: enddo
432: endif
435: if (ys + ym .eq. my) then
436: do i=xs,xs+xm-1
437: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
438: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
439: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
440: enddo
441: endif
444: if ((ys .eq. 0) .and. (xs .eq. 0)) then
445: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
446: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
447: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
448: endif
450: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
451: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
452: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
453: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
454: endif
456: ft = ft * area
457: call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, &
458: & MPIU_SUM,PETSC_COMM_WORLD,ierr)
461: ! Restore vectors
462: call VecRestoreArray(localX,x_v,x_i,ierr)
463: call VecRestoreArray(localV,g_v,g_i,ierr)
464: call VecRestoreArray(Left,left_v,left_i,ierr)
465: call VecRestoreArray(Top,top_v,top_i,ierr)
466: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
467: call VecRestoreArray(Right,right_v,right_i,ierr)
469: ! Scatter values to global vector
470: call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)
471: call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)
473: call PetscLogFlops(70.0d0*xm*ym,ierr)
475: return
476: end !FormFunctionGradient
482: ! ----------------------------------------------------------------------------
483: !
484: !
485: ! FormHessian - Evaluates Hessian matrix.
486: !
487: ! Input Parameters:
488: !. tao - the Tao context
489: !. X - input vector
490: !. dummy - not used
491: !
492: ! Output Parameters:
493: !. Hessian - Hessian matrix
494: !. Hpc - optionally different preconditioning matrix
495: !. flag - flag indicating matrix structure
496: !
497: ! Notes:
498: ! Due to mesh point reordering with DMs, we must always work
499: ! with the local mesh points, and then transform them to the new
500: ! global numbering with the local-to-global mapping. We cannot work
501: ! directly with the global numbers for the original uniprocessor mesh!
502: !
503: ! MatSetValuesLocal(), using the local ordering (including
504: ! ghost points!)
505: ! - Then set matrix entries using the local ordering
506: ! by calling MatSetValuesLocal()
508: subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
509: use mymodule
510: implicit none
512: Tao tao
513: Vec X
514: Mat Hessian,Hpc
515: PetscInt dummy
516: PetscErrorCode ierr
518: PetscInt i,j,k,row
519: PetscInt xs,xm,gxs,gxm
520: PetscInt ys,ym,gys,gym
521: PetscInt col(0:6)
522: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
523: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
524: PetscReal d4,d5,d6,d7,d8
525: PetscReal xc,xl,xr,xt,xb,xlt,xrb
526: PetscReal hl,hr,ht,hb,hc,htl,hbr
528: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
529: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
530: ! will return an array of doubles referenced by x_array offset by x_index.
531: ! i.e., to reference the kth element of X, use x_array(k + x_index).
532: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
533: PetscReal right_v(0:1),left_v(0:1)
534: PetscReal bottom_v(0:1),top_v(0:1)
535: PetscReal x_v(0:1)
536: PetscOffset x_i,right_i,left_i
537: PetscOffset bottom_i,top_i
538: PetscReal v(0:6)
539: PetscBool assembled
540: PetscInt i1
542: i1=1
544: ! Set various matrix options
545: call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, &
546: & PETSC_TRUE,ierr)
548: ! Get local mesh boundaries
549: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
550: & PETSC_NULL_INTEGER,ierr)
551: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
552: & PETSC_NULL_INTEGER,ierr)
554: ! Scatter ghost points to local vectors
555: call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)
556: call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)
558: ! Get pointers to vector data (see note on Fortran arrays above)
559: call VecGetArray(localX,x_v,x_i,ierr)
560: call VecGetArray(Top,top_v,top_i,ierr)
561: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
562: call VecGetArray(Left,left_v,left_i,ierr)
563: call VecGetArray(Right,right_v,right_i,ierr)
565: ! Initialize matrix entries to zero
566: call MatAssembled(Hessian,assembled,ierr)
567: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr)
570: rhx = real(mx) + 1.0
571: rhy = real(my) + 1.0
572: hx = 1.0/rhx
573: hy = 1.0/rhy
574: hydhx = hy/hx
575: hxdhy = hx/hy
576: ! compute Hessian over the locally owned part of the mesh
578: do i=xs,xs+xm-1
579: do j=ys,ys+ym-1
580: row = (j-gys)*gxm + (i-gxs)
582: xc = x_v(row + x_i)
583: xt = xc
584: xb = xc
585: xr = xc
586: xl = xc
587: xrb = xc
588: xlt = xc
590: if (i .eq. gxs) then ! Left side
591: xl = left_v(left_i + j - ys + 1)
592: xlt = left_v(left_i + j - ys + 2)
593: else
594: xl = x_v(x_i + row -1 )
595: endif
597: if (j .eq. gys) then ! bottom side
598: xb = bottom_v(bottom_i + i - xs + 1)
599: xrb = bottom_v(bottom_i + i - xs + 2)
600: else
601: xb = x_v(x_i + row - gxm)
602: endif
604: if (i+1 .eq. gxs + gxm) then !right side
605: xr = right_v(right_i + j - ys + 1)
606: xrb = right_v(right_i + j - ys)
607: else
608: xr = x_v(x_i + row + 1)
609: endif
611: if (j+1 .eq. gym+gys) then !top side
612: xt = top_v(top_i +i - xs + 1)
613: xlt = top_v(top_i + i - xs)
614: else
615: xt = x_v(x_i + row + gxm)
616: endif
618: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
619: xlt = x_v(x_i + row - 1 + gxm)
620: endif
622: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
623: xrb = x_v(x_i + row + 1 - gxm)
624: endif
626: d1 = (xc-xl)*rhx
627: d2 = (xc-xr)*rhx
628: d3 = (xc-xt)*rhy
629: d4 = (xc-xb)*rhy
630: d5 = (xrb-xr)*rhy
631: d6 = (xrb-xb)*rhx
632: d7 = (xlt-xl)*rhy
633: d8 = (xlt-xt)*rhx
635: f1 = sqrt( 1.0 + d1*d1 + d7*d7)
636: f2 = sqrt( 1.0 + d1*d1 + d4*d4)
637: f3 = sqrt( 1.0 + d3*d3 + d8*d8)
638: f4 = sqrt( 1.0 + d3*d3 + d2*d2)
639: f5 = sqrt( 1.0 + d2*d2 + d5*d5)
640: f6 = sqrt( 1.0 + d4*d4 + d6*d6)
643: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
644: & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
646: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
647: & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
649: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
650: & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
652: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
653: & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
655: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
656: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
658: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
659: & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
660: & hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
661: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
662: & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
663: & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
664: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
666: hl = hl * 0.5
667: hr = hr * 0.5
668: ht = ht * 0.5
669: hb = hb * 0.5
670: hbr = hbr * 0.5
671: htl = htl * 0.5
672: hc = hc * 0.5
674: k = 0
676: if (j .gt. 0) then
677: v(k) = hb
678: col(k) = row - gxm
679: k=k+1
680: endif
682: if ((j .gt. 0) .and. (i .lt. mx-1)) then
683: v(k) = hbr
684: col(k) = row-gxm+1
685: k=k+1
686: endif
688: if (i .gt. 0) then
689: v(k) = hl
690: col(k) = row - 1
691: k = k+1
692: endif
694: v(k) = hc
695: col(k) = row
696: k=k+1
698: if (i .lt. mx-1) then
699: v(k) = hr
700: col(k) = row + 1
701: k=k+1
702: endif
704: if ((i .gt. 0) .and. (j .lt. my-1)) then
705: v(k) = htl
706: col(k) = row + gxm - 1
707: k=k+1
708: endif
710: if (j .lt. my-1) then
711: v(k) = ht
712: col(k) = row + gxm
713: k=k+1
714: endif
716: ! Set matrix values using local numbering, defined earlier in main routine
717: call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, &
718: & ierr)
722: enddo
723: enddo
725: ! restore vectors
726: call VecRestoreArray(localX,x_v,x_i,ierr)
727: call VecRestoreArray(Left,left_v,left_i,ierr)
728: call VecRestoreArray(Right,right_v,right_i,ierr)
729: call VecRestoreArray(Top,top_v,top_i,ierr)
730: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
733: ! Assemble the matrix
734: call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)
735: call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)
737: call PetscLogFlops(199.0d0*xm*ym,ierr)
739: return
740: end
746: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
748: ! ----------------------------------------------------------------------------
749: !
750: !/*
751: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
752: !
753: !
754: !*/
756: subroutine MSA_BoundaryConditions(ierr)
757: use mymodule
758: implicit none
760: PetscErrorCode ierr
761: PetscInt i,j,k,limit,maxits
762: PetscInt xs, xm, gxs, gxm
763: PetscInt ys, ym, gys, gym
764: PetscInt bsize, lsize
765: PetscInt tsize, rsize
766: PetscReal one,two,three,tol
767: PetscReal scl,fnorm,det,xt
768: PetscReal yt,hx,hy,u1,u2,nf1,nf2
769: PetscReal njac11,njac12,njac21,njac22
770: PetscReal b, t, l, r
771: PetscReal boundary_v(0:1)
772: PetscOffset boundary_i
773: logical exitloop
774: PetscBool flg
776: limit=0
777: maxits = 5
778: tol=1e-10
779: b=-0.5
780: t= 0.5
781: l=-0.5
782: r= 0.5
783: xt=0
784: yt=0
785: one=1.0
786: two=2.0
787: three=3.0
790: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
791: & PETSC_NULL_INTEGER,ierr)
792: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, &
793: & gxm,gym,PETSC_NULL_INTEGER,ierr)
795: bsize = xm + 2
796: lsize = ym + 2
797: rsize = ym + 2
798: tsize = xm + 2
801: call VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)
802: call VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)
803: call VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)
804: call VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)
806: hx= (r-l)/(mx+1)
807: hy= (t-b)/(my+1)
809: do j=0,3
811: if (j.eq.0) then
812: yt=b
813: xt=l+hx*xs
814: limit=bsize
815: call VecGetArray(Bottom,boundary_v,boundary_i,ierr)
818: elseif (j.eq.1) then
819: yt=t
820: xt=l+hx*xs
821: limit=tsize
822: call VecGetArray(Top,boundary_v,boundary_i,ierr)
824: elseif (j.eq.2) then
825: yt=b+hy*ys
826: xt=l
827: limit=lsize
828: call VecGetArray(Left,boundary_v,boundary_i,ierr)
830: elseif (j.eq.3) then
831: yt=b+hy*ys
832: xt=r
833: limit=rsize
834: call VecGetArray(Right,boundary_v,boundary_i,ierr)
835: endif
838: do i=0,limit-1
840: u1=xt
841: u2=-yt
842: k = 0
843: exitloop = .false.
844: do while (k .lt. maxits .and. (.not. exitloop) )
846: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
847: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
848: fnorm=sqrt(nf1*nf1+nf2*nf2)
849: if (fnorm .gt. tol) then
850: njac11=one+u2*u2-u1*u1
851: njac12=two*u1*u2
852: njac21=-two*u1*u2
853: njac22=-one - u1*u1 + u2*u2
854: det = njac11*njac22-njac21*njac12
855: u1 = u1-(njac22*nf1-njac12*nf2)/det
856: u2 = u2-(njac11*nf2-njac21*nf1)/det
857: else
858: exitloop = .true.
859: endif
860: k=k+1
861: enddo
863: boundary_v(i + boundary_i) = u1*u1-u2*u2
864: if ((j .eq. 0) .or. (j .eq. 1)) then
865: xt = xt + hx
866: else
867: yt = yt + hy
868: endif
870: enddo
873: if (j.eq.0) then
874: call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr)
875: elseif (j.eq.1) then
876: call VecRestoreArray(Top,boundary_v,boundary_i,ierr)
877: elseif (j.eq.2) then
878: call VecRestoreArray(Left,boundary_v,boundary_i,ierr)
879: elseif (j.eq.3) then
880: call VecRestoreArray(Right,boundary_v,boundary_i,ierr)
881: endif
883: enddo
886: ! Scale the boundary if desired
887: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
888: & '-bottom',scl,flg,ierr)
889: if (flg .eqv. PETSC_TRUE) then
890: call VecScale(Bottom,scl,ierr)
891: endif
893: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
894: & '-top',scl,flg,ierr)
895: if (flg .eqv. PETSC_TRUE) then
896: call VecScale(Top,scl,ierr)
897: endif
899: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
900: & '-right',scl,flg,ierr)
901: if (flg .eqv. PETSC_TRUE) then
902: call VecScale(Right,scl,ierr)
903: endif
905: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
906: & '-left',scl,flg,ierr)
907: if (flg .eqv. PETSC_TRUE) then
908: call VecScale(Left,scl,ierr)
909: endif
912: return
913: end
915: ! ----------------------------------------------------------------------------
916: !
917: !/*
918: ! MSA_Plate - Calculates an obstacle for surface to stretch over
919: !
920: ! Output Parameter:
921: !. xl - lower bound vector
922: !. xu - upper bound vector
923: !
924: !*/
926: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
927: use mymodule
928: implicit none
930: Tao tao
931: Vec xl,xu
932: PetscErrorCode ierr
933: PetscInt i,j,row
934: PetscInt xs, xm, ys, ym
935: PetscReal lb,ub
936: PetscInt dummy
938: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
939: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
940: ! will return an array of doubles referenced by x_array offset by x_index.
941: ! i.e., to reference the kth element of X, use x_array(k + x_index).
942: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
943: PetscReal xl_v(0:1)
944: PetscOffset xl_i
946: lb = PETSC_NINFINITY
947: ub = PETSC_INFINITY
949: if (bmy .lt. 0) bmy = 0
950: if (bmy .gt. my) bmy = my
951: if (bmx .lt. 0) bmx = 0
952: if (bmx .gt. mx) bmx = mx
955: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
956: & PETSC_NULL_INTEGER,ierr)
958: call VecSet(xl,lb,ierr)
959: call VecSet(xu,ub,ierr)
961: call VecGetArray(xl,xl_v,xl_i,ierr)
964: do i=xs,xs+xm-1
966: do j=ys,ys+ym-1
968: row=(j-ys)*xm + (i-xs)
970: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
971: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
972: xl_v(xl_i+row) = bheight
974: endif
976: enddo
977: enddo
980: call VecRestoreArray(xl,xl_v,xl_i,ierr)
982: return
983: end
989: ! ----------------------------------------------------------------------------
990: !
991: !/*
992: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
993: !
994: ! Output Parameter:
995: !. X - vector for initial guess
996: !
997: !*/
999: subroutine MSA_InitialPoint(X, ierr)
1000: use mymodule
1001: implicit none
1003: Vec X
1004: PetscErrorCode ierr
1005: PetscInt start,i,j
1006: PetscInt row
1007: PetscInt xs,xm,gxs,gxm
1008: PetscInt ys,ym,gys,gym
1009: PetscReal zero, np5
1011: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
1012: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
1013: ! will return an array of doubles referenced by x_array offset by x_index.
1014: ! i.e., to reference the kth element of X, use x_array(k + x_index).
1015: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
1016: PetscReal left_v(0:1),right_v(0:1)
1017: PetscReal bottom_v(0:1),top_v(0:1)
1018: PetscReal x_v(0:1)
1019: PetscOffset left_i, right_i, top_i
1020: PetscOffset bottom_i,x_i
1021: PetscBool flg
1022: PetscRandom rctx
1024: zero = 0.0
1025: np5 = -0.5
1027: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
1028: & '-start', start,flg,ierr)
1030: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
1031: call VecSet(X,zero,ierr)
1033: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
1034: call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
1035: do i=0,start-1
1036: call VecSetRandom(X,rctx,ierr)
1037: enddo
1039: call PetscRandomDestroy(rctx,ierr)
1040: call VecShift(X,np5,ierr)
1042: else ! average of boundary conditions
1044: ! Get Local mesh boundaries
1045: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
1046: & PETSC_NULL_INTEGER,ierr)
1047: call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
1048: & PETSC_NULL_INTEGER,ierr)
1052: ! Get pointers to vector data
1053: call VecGetArray(Top,top_v,top_i,ierr)
1054: call VecGetArray(Bottom,bottom_v,bottom_i,ierr)
1055: call VecGetArray(Left,left_v,left_i,ierr)
1056: call VecGetArray(Right,right_v,right_i,ierr)
1058: call VecGetArray(localX,x_v,x_i,ierr)
1060: ! Perform local computations
1061: do j=ys,ys+ym-1
1062: do i=xs,xs+xm-1
1063: row = (j-gys)*gxm + (i-gxs)
1064: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
1065: & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
1066: & (i+1)*left_v(left_i+j-ys+1)/mx + &
1067: & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
1068: enddo
1069: enddo
1071: ! Restore vectors
1072: call VecRestoreArray(localX,x_v,x_i,ierr)
1074: call VecRestoreArray(Left,left_v,left_i,ierr)
1075: call VecRestoreArray(Top,top_v,top_i,ierr)
1076: call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr)
1077: call VecRestoreArray(Right,right_v,right_i,ierr)
1079: call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)
1080: call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)
1082: endif
1084: return
1085: end
1087: !
1088: !/*TEST
1089: !
1090: ! build:
1091: ! requires: !complex
1092: !
1093: ! test:
1094: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1095: ! filter: sort -b
1096: ! filter_output: sort -b
1097: ! requires: !single
1098: !
1099: ! test:
1100: ! suffix: 2
1101: ! nsize: 2
1102: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
1103: ! filter: sort -b
1104: ! filter_output: sort -b
1105: ! requires: !single
1106: !
1107: !TEST*/