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 plate2fmodule
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 plate2fmodule
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: PetscCallA(PetscInitialize(ierr))
65: ! Specify default dimensions of the problem
66: mx = 10
67: my = 10
68: bheight = 0.1
70: ! Check for any command line arguments that override defaults
72: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr))
73: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-my',my,flg,ierr))
75: bmx = mx/2
76: bmy = my/2
78: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-bmx',bmx,flg,ierr))
79: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-bmy',bmy,flg,ierr))
80: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-bheight',bheight,flg,ierr))
82: ! Calculate any derived values from parameters
83: N = mx*my
85: ! Let Petsc determine the dimensions of the local vectors
86: Nx = PETSC_DECIDE
87: NY = PETSC_DECIDE
89: ! A two dimensional distributed array will help define this problem, which
90: ! derives from an elliptic PDE on a two-dimensional domain. From the
91: ! distributed array, create the vectors
93: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,dm,ierr))
94: PetscCallA(DMSetFromOptions(dm,ierr))
95: PetscCallA(DMSetUp(dm,ierr))
97: ! Extract global and local vectors from DM; The local vectors are
98: ! used solely as work space for the evaluation of the function,
99: ! gradient, and Hessian. Duplicate for remaining vectors that are
100: ! the same types.
102: PetscCallA(DMCreateGlobalVector(dm,x,ierr))
103: PetscCallA(DMCreateLocalVector(dm,localX,ierr))
104: PetscCallA(VecDuplicate(localX,localV,ierr))
106: ! Create a matrix data structure to store the Hessian.
107: ! Here we (optionally) also associate the local numbering scheme
108: ! with the matrix so that later we can use local indices for matrix
109: ! assembly
111: PetscCallA(VecGetLocalSize(x,m,ierr))
112: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER,i3,PETSC_NULL_INTEGER,H,ierr))
114: PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))
115: PetscCallA(DMGetLocalToGlobalMapping(dm,isltog,ierr))
116: PetscCallA(MatSetLocalToGlobalMapping(H,isltog,isltog,ierr))
118: ! The Tao code begins here
119: ! Create TAO solver and set desired solution method.
120: ! This problems uses bounded variables, so the
121: ! method must either be 'tao_tron' or 'tao_blmvm'
123: PetscCallA(TaoCreate(PETSC_COMM_WORLD,tao,ierr))
124: PetscCallA(TaoSetType(tao,TAOBLMVM,ierr))
126: ! Set minimization function and gradient, hessian evaluation functions
128: PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr))
130: PetscCallA(TaoSetHessian(tao,H,H,FormHessian,0, ierr))
132: ! Set Variable bounds
133: PetscCallA(MSA_BoundaryConditions(ierr))
134: PetscCallA(TaoSetVariableBoundsRoutine(tao,MSA_Plate,0,ierr))
136: ! Set the initial solution guess
137: PetscCallA(MSA_InitialPoint(x, ierr))
138: PetscCallA(TaoSetSolution(tao,x,ierr))
140: ! Check for any tao command line options
141: PetscCallA(TaoSetFromOptions(tao,ierr))
143: ! Solve the application
144: PetscCallA(TaoSolve(tao,ierr))
146: ! Free TAO data structures
147: PetscCallA(TaoDestroy(tao,ierr))
149: ! Free PETSc data structures
150: PetscCallA(VecDestroy(x,ierr))
151: PetscCallA(VecDestroy(Top,ierr))
152: PetscCallA(VecDestroy(Bottom,ierr))
153: PetscCallA(VecDestroy(Left,ierr))
154: PetscCallA(VecDestroy(Right,ierr))
155: PetscCallA(MatDestroy(H,ierr))
156: PetscCallA(VecDestroy(localX,ierr))
157: PetscCallA(VecDestroy(localV,ierr))
158: PetscCallA(DMDestroy(dm,ierr))
160: ! Finalize TAO
162: PetscCallA(PetscFinalize(ierr))
164: end
166: ! ---------------------------------------------------------------------
167: !
168: ! FormFunctionGradient - Evaluates function f(X).
169: !
170: ! Input Parameters:
171: ! tao - the Tao context
172: ! X - the input vector
173: ! dummy - optional user-defined context, as set by TaoSetFunction()
174: ! (not used here)
175: !
176: ! Output Parameters:
177: ! fcn - the newly evaluated function
178: ! G - the gradient vector
179: ! info - error code
180: !
182: subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr)
183: use plate2fmodule
184: implicit none
186: ! Input/output variables
188: Tao tao
189: PetscReal fcn
190: Vec X, G
191: PetscErrorCode ierr
192: PetscInt dummy
194: PetscInt i,j,row
195: PetscInt xs, xm
196: PetscInt gxs, gxm
197: PetscInt ys, ym
198: PetscInt gys, gym
199: PetscReal ft,zero,hx,hy,hydhx,hxdhy
200: PetscReal area,rhx,rhy
201: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
202: PetscReal d4,d5,d6,d7,d8
203: PetscReal df1dxc,df2dxc,df3dxc,df4dxc
204: PetscReal df5dxc,df6dxc
205: PetscReal xc,xl,xr,xt,xb,xlt,xrb
207: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
208: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
209: ! will return an array of doubles referenced by x_array offset by x_index.
210: ! i.e., to reference the kth element of X, use x_array(k + x_index).
211: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
212: PetscReal g_v(0:1),x_v(0:1)
213: PetscReal top_v(0:1),left_v(0:1)
214: PetscReal right_v(0:1),bottom_v(0:1)
215: PetscOffset g_i,left_i,right_i
216: PetscOffset bottom_i,top_i,x_i
218: ft = 0.0
219: zero = 0.0
220: hx = 1.0/real(mx + 1)
221: hy = 1.0/real(my + 1)
222: hydhx = hy/hx
223: hxdhy = hx/hy
224: area = 0.5 * hx * hy
225: rhx = real(mx) + 1.0
226: rhy = real(my) + 1.0
228: ! Get local mesh boundaries
229: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
230: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
232: ! Scatter ghost points to local vector
233: PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
234: PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))
236: ! Initialize the vector to zero
237: PetscCall(VecSet(localV,zero,ierr))
239: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
240: PetscCall(VecGetArray(localX,x_v,x_i,ierr))
241: PetscCall(VecGetArray(localV,g_v,g_i,ierr))
242: PetscCall(VecGetArray(Top,top_v,top_i,ierr))
243: PetscCall(VecGetArray(Bottom,bottom_v,bottom_i,ierr))
244: PetscCall(VecGetArray(Left,left_v,left_i,ierr))
245: PetscCall(VecGetArray(Right,right_v,right_i,ierr))
247: ! Compute function over the locally owned part of the mesh
248: do j = ys,ys+ym-1
249: do i = xs,xs+xm-1
250: row = (j-gys)*gxm + (i-gxs)
251: xc = x_v(row+x_i)
252: xt = xc
253: xb = xc
254: xr = xc
255: xl = xc
256: xrb = xc
257: xlt = xc
259: if (i .eq. 0) then !left side
260: xl = left_v(j - ys + 1 + left_i)
261: xlt = left_v(j - ys + 2 + left_i)
262: else
263: xl = x_v(row - 1 + x_i)
264: endif
266: if (j .eq. 0) then !bottom side
267: xb = bottom_v(i - xs + 1 + bottom_i)
268: xrb = bottom_v(i - xs + 2 + bottom_i)
269: else
270: xb = x_v(row - gxm + x_i)
271: endif
273: if (i + 1 .eq. gxs + gxm) then !right side
274: xr = right_v(j - ys + 1 + right_i)
275: xrb = right_v(j - ys + right_i)
276: else
277: xr = x_v(row + 1 + x_i)
278: endif
280: if (j + 1 .eq. gys + gym) then !top side
281: xt = top_v(i - xs + 1 + top_i)
282: xlt = top_v(i - xs + top_i)
283: else
284: xt = x_v(row + gxm + x_i)
285: endif
287: if ((i .gt. gxs) .and. (j + 1 .lt. gys + gym)) then
288: xlt = x_v(row - 1 + gxm + x_i)
289: endif
291: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
292: xrb = x_v(row + 1 - gxm + x_i)
293: endif
295: d1 = xc-xl
296: d2 = xc-xr
297: d3 = xc-xt
298: d4 = xc-xb
299: d5 = xr-xrb
300: d6 = xrb-xb
301: d7 = xlt-xl
302: d8 = xt-xlt
304: df1dxc = d1 * hydhx
305: df2dxc = d1 * hydhx + d4 * hxdhy
306: df3dxc = d3 * hxdhy
307: df4dxc = d2 * hydhx + d3 * hxdhy
308: df5dxc = d2 * hydhx
309: df6dxc = d4 * hxdhy
311: d1 = d1 * rhx
312: d2 = d2 * rhx
313: d3 = d3 * rhy
314: d4 = d4 * rhy
315: d5 = d5 * rhy
316: d6 = d6 * rhx
317: d7 = d7 * rhy
318: d8 = d8 * rhx
320: f1 = sqrt(1.0 + d1*d1 + d7*d7)
321: f2 = sqrt(1.0 + d1*d1 + d4*d4)
322: f3 = sqrt(1.0 + d3*d3 + d8*d8)
323: f4 = sqrt(1.0 + d3*d3 + d2*d2)
324: f5 = sqrt(1.0 + d2*d2 + d5*d5)
325: f6 = sqrt(1.0 + d4*d4 + d6*d6)
327: ft = ft + f2 + f4
329: df1dxc = df1dxc / f1
330: df2dxc = df2dxc / f2
331: df3dxc = df3dxc / f3
332: df4dxc = df4dxc / f4
333: df5dxc = df5dxc / f5
334: df6dxc = df6dxc / f6
336: g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc)
337: enddo
338: enddo
340: ! Compute triangular areas along the border of the domain.
341: if (xs .eq. 0) then ! left side
342: do j=ys,ys+ym-1
343: d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) * rhy
344: d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) * rhx
345: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
346: enddo
347: endif
349: if (ys .eq. 0) then !bottom side
350: do i=xs,xs+xm-1
351: d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) * rhx
352: d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
353: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
354: enddo
355: endif
357: if (xs + xm .eq. mx) then ! right side
358: do j=ys,ys+ym-1
359: d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
360: d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
361: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
362: enddo
363: endif
365: if (ys + ym .eq. my) then
366: do i=xs,xs+xm-1
367: d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
368: d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
369: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
370: enddo
371: endif
373: if ((ys .eq. 0) .and. (xs .eq. 0)) then
374: d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
375: d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
376: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
377: endif
379: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
380: d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
381: d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
382: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
383: endif
385: ft = ft * area
386: PetscCallMPI(MPI_Allreduce(ft,fcn,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr))
388: ! Restore vectors
389: PetscCall(VecRestoreArray(localX,x_v,x_i,ierr))
390: PetscCall(VecRestoreArray(localV,g_v,g_i,ierr))
391: PetscCall(VecRestoreArray(Left,left_v,left_i,ierr))
392: PetscCall(VecRestoreArray(Top,top_v,top_i,ierr))
393: PetscCall(VecRestoreArray(Bottom,bottom_v,bottom_i,ierr))
394: PetscCall(VecRestoreArray(Right,right_v,right_i,ierr))
396: ! Scatter values to global vector
397: PetscCall(DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr))
398: PetscCall(DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr))
400: PetscCall(PetscLogFlops(70.0d0*xm*ym,ierr))
402: return
403: end !FormFunctionGradient
405: ! ----------------------------------------------------------------------------
406: !
407: !
408: ! FormHessian - Evaluates Hessian matrix.
409: !
410: ! Input Parameters:
411: !. tao - the Tao context
412: !. X - input vector
413: !. dummy - not used
414: !
415: ! Output Parameters:
416: !. Hessian - Hessian matrix
417: !. Hpc - optionally different preconditioning matrix
418: !. flag - flag indicating matrix structure
419: !
420: ! Notes:
421: ! Due to mesh point reordering with DMs, we must always work
422: ! with the local mesh points, and then transform them to the new
423: ! global numbering with the local-to-global mapping. We cannot work
424: ! directly with the global numbers for the original uniprocessor mesh!
425: !
426: ! MatSetValuesLocal(), using the local ordering (including
427: ! ghost points!)
428: ! - Then set matrix entries using the local ordering
429: ! by calling MatSetValuesLocal()
431: subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
432: use plate2fmodule
433: implicit none
435: Tao tao
436: Vec X
437: Mat Hessian,Hpc
438: PetscInt dummy
439: PetscErrorCode ierr
441: PetscInt i,j,k,row
442: PetscInt xs,xm,gxs,gxm
443: PetscInt ys,ym,gys,gym
444: PetscInt col(0:6)
445: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
446: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
447: PetscReal d4,d5,d6,d7,d8
448: PetscReal xc,xl,xr,xt,xb,xlt,xrb
449: PetscReal hl,hr,ht,hb,hc,htl,hbr
451: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
452: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
453: ! will return an array of doubles referenced by x_array offset by x_index.
454: ! i.e., to reference the kth element of X, use x_array(k + x_index).
455: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
456: PetscReal right_v(0:1),left_v(0:1)
457: PetscReal bottom_v(0:1),top_v(0:1)
458: PetscReal x_v(0:1)
459: PetscOffset x_i,right_i,left_i
460: PetscOffset bottom_i,top_i
461: PetscReal v(0:6)
462: PetscBool assembled
463: PetscInt i1
465: i1=1
467: ! Set various matrix options
468: PetscCall(MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE,ierr))
470: ! Get local mesh boundaries
471: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
472: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
474: ! Scatter ghost points to local vectors
475: PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
476: PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))
478: ! Get pointers to vector data (see note on Fortran arrays above)
479: PetscCall(VecGetArray(localX,x_v,x_i,ierr))
480: PetscCall(VecGetArray(Top,top_v,top_i,ierr))
481: PetscCall(VecGetArray(Bottom,bottom_v,bottom_i,ierr))
482: PetscCall(VecGetArray(Left,left_v,left_i,ierr))
483: PetscCall(VecGetArray(Right,right_v,right_i,ierr))
485: ! Initialize matrix entries to zero
486: PetscCall(MatAssembled(Hessian,assembled,ierr))
487: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian,ierr))
489: rhx = real(mx) + 1.0
490: rhy = real(my) + 1.0
491: hx = 1.0/rhx
492: hy = 1.0/rhy
493: hydhx = hy/hx
494: hxdhy = hx/hy
495: ! compute Hessian over the locally owned part of the mesh
497: do i=xs,xs+xm-1
498: do j=ys,ys+ym-1
499: row = (j-gys)*gxm + (i-gxs)
501: xc = x_v(row + x_i)
502: xt = xc
503: xb = xc
504: xr = xc
505: xl = xc
506: xrb = xc
507: xlt = xc
509: if (i .eq. gxs) then ! Left side
510: xl = left_v(left_i + j - ys + 1)
511: xlt = left_v(left_i + j - ys + 2)
512: else
513: xl = x_v(x_i + row -1)
514: endif
516: if (j .eq. gys) then ! bottom side
517: xb = bottom_v(bottom_i + i - xs + 1)
518: xrb = bottom_v(bottom_i + i - xs + 2)
519: else
520: xb = x_v(x_i + row - gxm)
521: endif
523: if (i+1 .eq. gxs + gxm) then !right side
524: xr = right_v(right_i + j - ys + 1)
525: xrb = right_v(right_i + j - ys)
526: else
527: xr = x_v(x_i + row + 1)
528: endif
530: if (j+1 .eq. gym+gys) then !top side
531: xt = top_v(top_i +i - xs + 1)
532: xlt = top_v(top_i + i - xs)
533: else
534: xt = x_v(x_i + row + gxm)
535: endif
537: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
538: xlt = x_v(x_i + row - 1 + gxm)
539: endif
541: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
542: xrb = x_v(x_i + row + 1 - gxm)
543: endif
545: d1 = (xc-xl)*rhx
546: d2 = (xc-xr)*rhx
547: d3 = (xc-xt)*rhy
548: d4 = (xc-xb)*rhy
549: d5 = (xrb-xr)*rhy
550: d6 = (xrb-xb)*rhx
551: d7 = (xlt-xl)*rhy
552: d8 = (xlt-xt)*rhx
554: f1 = sqrt(1.0 + d1*d1 + d7*d7)
555: f2 = sqrt(1.0 + d1*d1 + d4*d4)
556: f3 = sqrt(1.0 + d3*d3 + d8*d8)
557: f4 = sqrt(1.0 + d3*d3 + d2*d2)
558: f5 = sqrt(1.0 + d2*d2 + d5*d5)
559: f6 = sqrt(1.0 + d4*d4 + d6*d6)
561: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
563: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
565: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
567: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
569: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
570: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
572: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
573: & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
574: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
576: hl = hl * 0.5
577: hr = hr * 0.5
578: ht = ht * 0.5
579: hb = hb * 0.5
580: hbr = hbr * 0.5
581: htl = htl * 0.5
582: hc = hc * 0.5
584: k = 0
586: if (j .gt. 0) then
587: v(k) = hb
588: col(k) = row - gxm
589: k=k+1
590: endif
592: if ((j .gt. 0) .and. (i .lt. mx-1)) then
593: v(k) = hbr
594: col(k) = row-gxm+1
595: k=k+1
596: endif
598: if (i .gt. 0) then
599: v(k) = hl
600: col(k) = row - 1
601: k = k+1
602: endif
604: v(k) = hc
605: col(k) = row
606: k=k+1
608: if (i .lt. mx-1) then
609: v(k) = hr
610: col(k) = row + 1
611: k=k+1
612: endif
614: if ((i .gt. 0) .and. (j .lt. my-1)) then
615: v(k) = htl
616: col(k) = row + gxm - 1
617: k=k+1
618: endif
620: if (j .lt. my-1) then
621: v(k) = ht
622: col(k) = row + gxm
623: k=k+1
624: endif
626: ! Set matrix values using local numbering, defined earlier in main routine
627: PetscCall(MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES,ierr))
629: enddo
630: enddo
632: ! restore vectors
633: PetscCall(VecRestoreArray(localX,x_v,x_i,ierr))
634: PetscCall(VecRestoreArray(Left,left_v,left_i,ierr))
635: PetscCall(VecRestoreArray(Right,right_v,right_i,ierr))
636: PetscCall(VecRestoreArray(Top,top_v,top_i,ierr))
637: PetscCall(VecRestoreArray(Bottom,bottom_v,bottom_i,ierr))
639: ! Assemble the matrix
640: PetscCall(MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr))
641: PetscCall(MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr))
643: PetscCall(PetscLogFlops(199.0d0*xm*ym,ierr))
645: return
646: end
648: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
650: ! ----------------------------------------------------------------------------
651: !
652: !/*
653: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
654: !
655: !
656: !*/
658: subroutine MSA_BoundaryConditions(ierr)
659: use plate2fmodule
660: implicit none
662: PetscErrorCode ierr
663: PetscInt i,j,k,limit,maxits
664: PetscInt xs, xm, gxs, gxm
665: PetscInt ys, ym, gys, gym
666: PetscInt bsize, lsize
667: PetscInt tsize, rsize
668: PetscReal one,two,three,tol
669: PetscReal scl,fnorm,det,xt
670: PetscReal yt,hx,hy,u1,u2,nf1,nf2
671: PetscReal njac11,njac12,njac21,njac22
672: PetscReal b, t, l, r
673: PetscReal boundary_v(0:1)
674: PetscOffset boundary_i
675: logical exitloop
676: PetscBool flg
678: limit=0
679: maxits = 5
680: tol=1e-10
681: b=-0.5
682: t= 0.5
683: l=-0.5
684: r= 0.5
685: xt=0
686: yt=0
687: one=1.0
688: two=2.0
689: three=3.0
691: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
692: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
694: bsize = xm + 2
695: lsize = ym + 2
696: rsize = ym + 2
697: tsize = xm + 2
699: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr))
700: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr))
701: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr))
702: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr))
704: hx= (r-l)/(mx+1)
705: hy= (t-b)/(my+1)
707: do j=0,3
709: if (j.eq.0) then
710: yt=b
711: xt=l+hx*xs
712: limit=bsize
713: PetscCall(VecGetArray(Bottom,boundary_v,boundary_i,ierr))
715: elseif (j.eq.1) then
716: yt=t
717: xt=l+hx*xs
718: limit=tsize
719: PetscCall(VecGetArray(Top,boundary_v,boundary_i,ierr))
721: elseif (j.eq.2) then
722: yt=b+hy*ys
723: xt=l
724: limit=lsize
725: PetscCall(VecGetArray(Left,boundary_v,boundary_i,ierr))
727: elseif (j.eq.3) then
728: yt=b+hy*ys
729: xt=r
730: limit=rsize
731: PetscCall(VecGetArray(Right,boundary_v,boundary_i,ierr))
732: endif
734: do i=0,limit-1
736: u1=xt
737: u2=-yt
738: k = 0
739: exitloop = .false.
740: do while (k .lt. maxits .and. (.not. exitloop))
742: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
743: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
744: fnorm=sqrt(nf1*nf1+nf2*nf2)
745: if (fnorm .gt. tol) then
746: njac11=one+u2*u2-u1*u1
747: njac12=two*u1*u2
748: njac21=-two*u1*u2
749: njac22=-one - u1*u1 + u2*u2
750: det = njac11*njac22-njac21*njac12
751: u1 = u1-(njac22*nf1-njac12*nf2)/det
752: u2 = u2-(njac11*nf2-njac21*nf1)/det
753: else
754: exitloop = .true.
755: endif
756: k=k+1
757: enddo
759: boundary_v(i + boundary_i) = u1*u1-u2*u2
760: if ((j .eq. 0) .or. (j .eq. 1)) then
761: xt = xt + hx
762: else
763: yt = yt + hy
764: endif
766: enddo
768: if (j.eq.0) then
769: PetscCall(VecRestoreArray(Bottom,boundary_v,boundary_i,ierr))
770: elseif (j.eq.1) then
771: PetscCall(VecRestoreArray(Top,boundary_v,boundary_i,ierr))
772: elseif (j.eq.2) then
773: PetscCall(VecRestoreArray(Left,boundary_v,boundary_i,ierr))
774: elseif (j.eq.3) then
775: PetscCall(VecRestoreArray(Right,boundary_v,boundary_i,ierr))
776: endif
778: enddo
780: ! Scale the boundary if desired
781: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-bottom',scl,flg,ierr))
782: if (flg .eqv. PETSC_TRUE) then
783: PetscCall(VecScale(Bottom,scl,ierr))
784: endif
786: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-top',scl,flg,ierr))
787: if (flg .eqv. PETSC_TRUE) then
788: PetscCall(VecScale(Top,scl,ierr))
789: endif
791: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-right',scl,flg,ierr))
792: if (flg .eqv. PETSC_TRUE) then
793: PetscCall(VecScale(Right,scl,ierr))
794: endif
796: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-left',scl,flg,ierr))
797: if (flg .eqv. PETSC_TRUE) then
798: PetscCall(VecScale(Left,scl,ierr))
799: endif
801: return
802: end
804: ! ----------------------------------------------------------------------------
805: !
806: !/*
807: ! MSA_Plate - Calculates an obstacle for surface to stretch over
808: !
809: ! Output Parameter:
810: !. xl - lower bound vector
811: !. xu - upper bound vector
812: !
813: !*/
815: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
816: use plate2fmodule
817: implicit none
819: Tao tao
820: Vec xl,xu
821: PetscErrorCode ierr
822: PetscInt i,j,row
823: PetscInt xs, xm, ys, ym
824: PetscReal lb,ub
825: PetscInt dummy
827: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
828: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
829: ! will return an array of doubles referenced by x_array offset by x_index.
830: ! i.e., to reference the kth element of X, use x_array(k + x_index).
831: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
832: PetscReal xl_v(0:1)
833: PetscOffset xl_i
835: lb = PETSC_NINFINITY
836: ub = PETSC_INFINITY
838: if (bmy .lt. 0) bmy = 0
839: if (bmy .gt. my) bmy = my
840: if (bmx .lt. 0) bmx = 0
841: if (bmx .gt. mx) bmx = mx
843: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
845: PetscCall(VecSet(xl,lb,ierr))
846: PetscCall(VecSet(xu,ub,ierr))
848: PetscCall(VecGetArray(xl,xl_v,xl_i,ierr))
850: do i=xs,xs+xm-1
852: do j=ys,ys+ym-1
854: row=(j-ys)*xm + (i-xs)
856: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
857: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
858: xl_v(xl_i+row) = bheight
860: endif
862: enddo
863: enddo
865: PetscCall(VecRestoreArray(xl,xl_v,xl_i,ierr))
867: return
868: end
870: ! ----------------------------------------------------------------------------
871: !
872: !/*
873: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
874: !
875: ! Output Parameter:
876: !. X - vector for initial guess
877: !
878: !*/
880: subroutine MSA_InitialPoint(X, ierr)
881: use plate2fmodule
882: implicit none
884: Vec X
885: PetscErrorCode ierr
886: PetscInt start,i,j
887: PetscInt row
888: PetscInt xs,xm,gxs,gxm
889: PetscInt ys,ym,gys,gym
890: PetscReal zero, np5
892: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
893: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr)
894: ! will return an array of doubles referenced by x_array offset by x_index.
895: ! i.e., to reference the kth element of X, use x_array(k + x_index).
896: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
897: PetscReal left_v(0:1),right_v(0:1)
898: PetscReal bottom_v(0:1),top_v(0:1)
899: PetscReal x_v(0:1)
900: PetscOffset left_i, right_i, top_i
901: PetscOffset bottom_i,x_i
902: PetscBool flg
903: PetscRandom rctx
905: zero = 0.0
906: np5 = -0.5
908: PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-start', start,flg,ierr))
910: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
911: PetscCall(VecSet(X,zero,ierr))
913: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
914: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr))
915: do i=0,start-1
916: PetscCall(VecSetRandom(X,rctx,ierr))
917: enddo
919: PetscCall(PetscRandomDestroy(rctx,ierr))
920: PetscCall(VecShift(X,np5,ierr))
922: else ! average of boundary conditions
924: ! Get Local mesh boundaries
925: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
926: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
928: ! Get pointers to vector data
929: PetscCall(VecGetArray(Top,top_v,top_i,ierr))
930: PetscCall(VecGetArray(Bottom,bottom_v,bottom_i,ierr))
931: PetscCall(VecGetArray(Left,left_v,left_i,ierr))
932: PetscCall(VecGetArray(Right,right_v,right_i,ierr))
934: PetscCall(VecGetArray(localX,x_v,x_i,ierr))
936: ! Perform local computations
937: do j=ys,ys+ym-1
938: do i=xs,xs+xm-1
939: row = (j-gys)*gxm + (i-gxs)
940: x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
941: & (i+1)*left_v(left_i+j-ys+1)/mx + (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
942: enddo
943: enddo
945: ! Restore vectors
946: PetscCall(VecRestoreArray(localX,x_v,x_i,ierr))
948: PetscCall(VecRestoreArray(Left,left_v,left_i,ierr))
949: PetscCall(VecRestoreArray(Top,top_v,top_i,ierr))
950: PetscCall(VecRestoreArray(Bottom,bottom_v,bottom_i,ierr))
951: PetscCall(VecRestoreArray(Right,right_v,right_i,ierr))
953: PetscCall(DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr))
954: PetscCall(DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr))
956: endif
958: return
959: end
961: !
962: !/*TEST
963: !
964: ! build:
965: ! requires: !complex
966: !
967: ! test:
968: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
969: ! filter: sort -b
970: ! filter_output: sort -b
971: ! requires: !single
972: !
973: ! test:
974: ! suffix: 2
975: ! nsize: 2
976: ! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
977: ! filter: sort -b
978: ! filter_output: sort -b
979: ! requires: !single
980: !
981: !TEST*/