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: PetscReal, pointer :: g_v(:),x_v(:)
208: PetscReal, pointer :: top_v(:),left_v(:)
209: PetscReal, pointer :: right_v(:),bottom_v(:)
211: ft = 0.0
212: zero = 0.0
213: hx = 1.0/real(mx + 1)
214: hy = 1.0/real(my + 1)
215: hydhx = hy/hx
216: hxdhy = hx/hy
217: area = 0.5 * hx * hy
218: rhx = real(mx) + 1.0
219: rhy = real(my) + 1.0
221: ! Get local mesh boundaries
222: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
223: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
225: ! Scatter ghost points to local vector
226: PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
227: PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))
229: ! Initialize the vector to zero
230: PetscCall(VecSet(localV,zero,ierr))
232: ! Get arrays to vector data (See note above about using VecGetArrayF90 in Fortran)
233: PetscCall(VecGetArrayF90(localX,x_v,ierr))
234: PetscCall(VecGetArrayF90(localV,g_v,ierr))
235: PetscCall(VecGetArrayF90(Top,top_v,ierr))
236: PetscCall(VecGetArrayF90(Bottom,bottom_v,ierr))
237: PetscCall(VecGetArrayF90(Left,left_v,ierr))
238: PetscCall(VecGetArrayF90(Right,right_v,ierr))
240: ! Compute function over the locally owned part of the mesh
241: do j = ys,ys+ym-1
242: do i = xs,xs+xm-1
243: row = (j-gys)*gxm + (i-gxs)
244: xc = x_v(1+row)
245: xt = xc
246: xb = xc
247: xr = xc
248: xl = xc
249: xrb = xc
250: xlt = xc
252: if (i .eq. 0) then !left side
253: xl = left_v(1+j - ys + 1)
254: xlt = left_v(1+j - ys + 2)
255: else
256: xl = x_v(1+row - 1)
257: endif
259: if (j .eq. 0) then !bottom side
260: xb = bottom_v(1+i - xs + 1)
261: xrb = bottom_v(1+i - xs + 2)
262: else
263: xb = x_v(1+row - gxm)
264: endif
266: if (i + 1 .eq. gxs + gxm) then !right side
267: xr = right_v(1+j - ys + 1)
268: xrb = right_v(1+j - ys)
269: else
270: xr = x_v(1+row + 1)
271: endif
273: if (j + 1 .eq. gys + gym) then !top side
274: xt = top_v(1+i - xs + 1)
275: xlt = top_v(1+i - xs)
276: else
277: xt = x_v(1+row + gxm)
278: endif
280: if ((i .gt. gxs) .and. (j + 1 .lt. gys + gym)) then
281: xlt = x_v(1+row - 1 + gxm)
282: endif
284: if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
285: xrb = x_v(1+row + 1 - gxm)
286: endif
288: d1 = xc-xl
289: d2 = xc-xr
290: d3 = xc-xt
291: d4 = xc-xb
292: d5 = xr-xrb
293: d6 = xrb-xb
294: d7 = xlt-xl
295: d8 = xt-xlt
297: df1dxc = d1 * hydhx
298: df2dxc = d1 * hydhx + d4 * hxdhy
299: df3dxc = d3 * hxdhy
300: df4dxc = d2 * hydhx + d3 * hxdhy
301: df5dxc = d2 * hydhx
302: df6dxc = d4 * hxdhy
304: d1 = d1 * rhx
305: d2 = d2 * rhx
306: d3 = d3 * rhy
307: d4 = d4 * rhy
308: d5 = d5 * rhy
309: d6 = d6 * rhx
310: d7 = d7 * rhy
311: d8 = d8 * rhx
313: f1 = sqrt(1.0 + d1*d1 + d7*d7)
314: f2 = sqrt(1.0 + d1*d1 + d4*d4)
315: f3 = sqrt(1.0 + d3*d3 + d8*d8)
316: f4 = sqrt(1.0 + d3*d3 + d2*d2)
317: f5 = sqrt(1.0 + d2*d2 + d5*d5)
318: f6 = sqrt(1.0 + d4*d4 + d6*d6)
320: ft = ft + f2 + f4
322: df1dxc = df1dxc / f1
323: df2dxc = df2dxc / f2
324: df3dxc = df3dxc / f3
325: df4dxc = df4dxc / f4
326: df5dxc = df5dxc / f5
327: df6dxc = df6dxc / f6
329: g_v(1+row) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc)
330: enddo
331: enddo
333: ! Compute triangular areas along the border of the domain.
334: if (xs .eq. 0) then ! left side
335: do j=ys,ys+ym-1
336: d3 = (left_v(1+j-ys+1) - left_v(1+j-ys+2)) * rhy
337: d2 = (left_v(1+j-ys+1) - x_v(1+(j-gys)*gxm)) * rhx
338: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
339: enddo
340: endif
342: if (ys .eq. 0) then !bottom side
343: do i=xs,xs+xm-1
344: d2 = (bottom_v(1+i+1-xs)-bottom_v(1+i-xs+2)) * rhx
345: d3 = (bottom_v(1+i-xs+1)-x_v(1+i-gxs))*rhy
346: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
347: enddo
348: endif
350: if (xs + xm .eq. mx) then ! right side
351: do j=ys,ys+ym-1
352: d1 = (x_v(1+(j+1-gys)*gxm-1)-right_v(1+j-ys+1))*rhx
353: d4 = (right_v(1+j-ys) - right_v(1+j-ys+1))*rhy
354: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
355: enddo
356: endif
358: if (ys + ym .eq. my) then
359: do i=xs,xs+xm-1
360: d1 = (x_v(1+(gym-1)*gxm+i-gxs) - top_v(1+i-xs+1))*rhy
361: d4 = (top_v(1+i-xs+1) - top_v(1+i-xs))*rhx
362: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
363: enddo
364: endif
366: if ((ys .eq. 0) .and. (xs .eq. 0)) then
367: d1 = (left_v(1+0) - left_v(1+1)) * rhy
368: d2 = (bottom_v(1+0)-bottom_v(1+1))*rhx
369: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
370: endif
372: if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
373: d1 = (right_v(1+ym+1) - right_v(1+ym))*rhy
374: d2 = (top_v(1+xm+1) - top_v(1+xm))*rhx
375: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
376: endif
378: ft = ft * area
379: PetscCallMPI(MPI_Allreduce(ft,fcn,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr))
381: ! Restore vectors
382: PetscCall(VecRestoreArrayF90(localX,x_v,ierr))
383: PetscCall(VecRestoreArrayF90(localV,g_v,ierr))
384: PetscCall(VecRestoreArrayF90(Left,left_v,ierr))
385: PetscCall(VecRestoreArrayF90(Top,top_v,ierr))
386: PetscCall(VecRestoreArrayF90(Bottom,bottom_v,ierr))
387: PetscCall(VecRestoreArrayF90(Right,right_v,ierr))
389: ! Scatter values to global vector
390: PetscCall(DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr))
391: PetscCall(DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr))
393: PetscCall(PetscLogFlops(70.0d0*xm*ym,ierr))
395: end !FormFunctionGradient
397: ! ----------------------------------------------------------------------------
398: !
399: !
400: ! FormHessian - Evaluates Hessian matrix.
401: !
402: ! Input Parameters:
403: !. tao - the Tao context
404: !. X - input vector
405: !. dummy - not used
406: !
407: ! Output Parameters:
408: !. Hessian - Hessian matrix
409: !. Hpc - optionally different preconditioning matrix
410: !. flag - flag indicating matrix structure
411: !
412: ! Notes:
413: ! Due to mesh point reordering with DMs, we must always work
414: ! with the local mesh points, and then transform them to the new
415: ! global numbering with the local-to-global mapping. We cannot work
416: ! directly with the global numbers for the original uniprocessor mesh!
417: !
418: ! MatSetValuesLocal(), using the local ordering (including
419: ! ghost points!)
420: ! - Then set matrix entries using the local ordering
421: ! by calling MatSetValuesLocal()
423: subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr)
424: use plate2fmodule
425: implicit none
427: Tao tao
428: Vec X
429: Mat Hessian,Hpc
430: PetscInt dummy
431: PetscErrorCode ierr
433: PetscInt i,j,k,row
434: PetscInt xs,xm,gxs,gxm
435: PetscInt ys,ym,gys,gym
436: PetscInt col(0:6)
437: PetscReal hx,hy,hydhx,hxdhy,rhx,rhy
438: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3
439: PetscReal d4,d5,d6,d7,d8
440: PetscReal xc,xl,xr,xt,xb,xlt,xrb
441: PetscReal hl,hr,ht,hb,hc,htl,hbr
443: PetscReal,pointer :: right_v(:),left_v(:)
444: PetscReal,pointer :: bottom_v(:),top_v(:)
445: PetscReal,pointer :: x_v(:)
446: PetscReal v(0:6)
447: PetscBool assembled
448: PetscInt i1
450: i1=1
452: ! Set various matrix options
453: PetscCall(MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE,ierr))
455: ! Get local mesh boundaries
456: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
457: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
459: ! Scatter ghost points to local vectors
460: PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr))
461: PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr))
463: ! Get pointers to vector data (see note on Fortran arrays above)
464: PetscCall(VecGetArrayF90(localX,x_v,ierr))
465: PetscCall(VecGetArrayF90(Top,top_v,ierr))
466: PetscCall(VecGetArrayF90(Bottom,bottom_v,ierr))
467: PetscCall(VecGetArrayF90(Left,left_v,ierr))
468: PetscCall(VecGetArrayF90(Right,right_v,ierr))
470: ! Initialize matrix entries to zero
471: PetscCall(MatAssembled(Hessian,assembled,ierr))
472: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian,ierr))
474: rhx = real(mx) + 1.0
475: rhy = real(my) + 1.0
476: hx = 1.0/rhx
477: hy = 1.0/rhy
478: hydhx = hy/hx
479: hxdhy = hx/hy
480: ! compute Hessian over the locally owned part of the mesh
482: do i=xs,xs+xm-1
483: do j=ys,ys+ym-1
484: row = (j-gys)*gxm + (i-gxs)
486: xc = x_v(1+row)
487: xt = xc
488: xb = xc
489: xr = xc
490: xl = xc
491: xrb = xc
492: xlt = xc
494: if (i .eq. gxs) then ! Left side
495: xl = left_v(1+j - ys + 1)
496: xlt = left_v(1+j - ys + 2)
497: else
498: xl = x_v(1+row -1)
499: endif
501: if (j .eq. gys) then ! bottom side
502: xb = bottom_v(1+i - xs + 1)
503: xrb = bottom_v(1+i - xs + 2)
504: else
505: xb = x_v(1+row - gxm)
506: endif
508: if (i+1 .eq. gxs + gxm) then !right side
509: xr = right_v(1+j - ys + 1)
510: xrb = right_v(1+j - ys)
511: else
512: xr = x_v(1+row + 1)
513: endif
515: if (j+1 .eq. gym+gys) then !top side
516: xt = top_v(1+i - xs + 1)
517: xlt = top_v(1+i - xs)
518: else
519: xt = x_v(1+row + gxm)
520: endif
522: if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
523: xlt = x_v(1+row - 1 + gxm)
524: endif
526: if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
527: xrb = x_v(1+row + 1 - gxm)
528: endif
530: d1 = (xc-xl)*rhx
531: d2 = (xc-xr)*rhx
532: d3 = (xc-xt)*rhy
533: d4 = (xc-xb)*rhy
534: d5 = (xrb-xr)*rhy
535: d6 = (xrb-xb)*rhx
536: d7 = (xlt-xl)*rhy
537: d8 = (xlt-xt)*rhx
539: f1 = sqrt(1.0 + d1*d1 + d7*d7)
540: f2 = sqrt(1.0 + d1*d1 + d4*d4)
541: f3 = sqrt(1.0 + d3*d3 + d8*d8)
542: f4 = sqrt(1.0 + d3*d3 + d2*d2)
543: f5 = sqrt(1.0 + d2*d2 + d5*d5)
544: f6 = sqrt(1.0 + d4*d4 + d6*d6)
546: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
548: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
550: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
552: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
554: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
555: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
557: 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) + &
558: & 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)+ &
559: & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
561: hl = hl * 0.5
562: hr = hr * 0.5
563: ht = ht * 0.5
564: hb = hb * 0.5
565: hbr = hbr * 0.5
566: htl = htl * 0.5
567: hc = hc * 0.5
569: k = 0
571: if (j .gt. 0) then
572: v(k) = hb
573: col(k) = row - gxm
574: k=k+1
575: endif
577: if ((j .gt. 0) .and. (i .lt. mx-1)) then
578: v(k) = hbr
579: col(k) = row-gxm+1
580: k=k+1
581: endif
583: if (i .gt. 0) then
584: v(k) = hl
585: col(k) = row - 1
586: k = k+1
587: endif
589: v(k) = hc
590: col(k) = row
591: k=k+1
593: if (i .lt. mx-1) then
594: v(k) = hr
595: col(k) = row + 1
596: k=k+1
597: endif
599: if ((i .gt. 0) .and. (j .lt. my-1)) then
600: v(k) = htl
601: col(k) = row + gxm - 1
602: k=k+1
603: endif
605: if (j .lt. my-1) then
606: v(k) = ht
607: col(k) = row + gxm
608: k=k+1
609: endif
611: ! Set matrix values using local numbering, defined earlier in main routine
612: PetscCall(MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES,ierr))
614: enddo
615: enddo
617: ! restore vectors
618: PetscCall(VecRestoreArrayF90(localX,x_v,ierr))
619: PetscCall(VecRestoreArrayF90(Left,left_v,ierr))
620: PetscCall(VecRestoreArrayF90(Right,right_v,ierr))
621: PetscCall(VecRestoreArrayF90(Top,top_v,ierr))
622: PetscCall(VecRestoreArrayF90(Bottom,bottom_v,ierr))
624: ! Assemble the matrix
625: PetscCall(MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr))
626: PetscCall(MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr))
628: PetscCall(PetscLogFlops(199.0d0*xm*ym,ierr))
630: end
632: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
634: ! ----------------------------------------------------------------------------
635: !
636: !/*
637: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
638: !
639: !
640: !*/
642: subroutine MSA_BoundaryConditions(ierr)
643: use plate2fmodule
644: implicit none
646: PetscErrorCode ierr
647: PetscInt i,j,k,limit,maxits
648: PetscInt xs, xm, gxs, gxm
649: PetscInt ys, ym, gys, gym
650: PetscInt bsize, lsize
651: PetscInt tsize, rsize
652: PetscReal one,two,three,tol
653: PetscReal scl,fnorm,det,xt
654: PetscReal yt,hx,hy,u1,u2,nf1,nf2
655: PetscReal njac11,njac12,njac21,njac22
656: PetscReal b, t, l, r
657: PetscReal,pointer :: boundary_v(:)
659: logical exitloop
660: PetscBool flg
662: limit=0
663: maxits = 5
664: tol=1e-10
665: b=-0.5
666: t= 0.5
667: l=-0.5
668: r= 0.5
669: xt=0
670: yt=0
671: one=1.0
672: two=2.0
673: three=3.0
675: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
676: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
678: bsize = xm + 2
679: lsize = ym + 2
680: rsize = ym + 2
681: tsize = xm + 2
683: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr))
684: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr))
685: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr))
686: PetscCall(VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr))
688: hx= (r-l)/(mx+1)
689: hy= (t-b)/(my+1)
691: do j=0,3
693: if (j.eq.0) then
694: yt=b
695: xt=l+hx*xs
696: limit=bsize
697: PetscCall(VecGetArrayF90(Bottom,boundary_v,ierr))
699: elseif (j.eq.1) then
700: yt=t
701: xt=l+hx*xs
702: limit=tsize
703: PetscCall(VecGetArrayF90(Top,boundary_v,ierr))
705: elseif (j.eq.2) then
706: yt=b+hy*ys
707: xt=l
708: limit=lsize
709: PetscCall(VecGetArrayF90(Left,boundary_v,ierr))
711: elseif (j.eq.3) then
712: yt=b+hy*ys
713: xt=r
714: limit=rsize
715: PetscCall(VecGetArrayF90(Right,boundary_v,ierr))
716: endif
718: do i=0,limit-1
720: u1=xt
721: u2=-yt
722: k = 0
723: exitloop = .false.
724: do while (k .lt. maxits .and. (.not. exitloop))
726: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
727: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
728: fnorm=sqrt(nf1*nf1+nf2*nf2)
729: if (fnorm .gt. tol) then
730: njac11=one+u2*u2-u1*u1
731: njac12=two*u1*u2
732: njac21=-two*u1*u2
733: njac22=-one - u1*u1 + u2*u2
734: det = njac11*njac22-njac21*njac12
735: u1 = u1-(njac22*nf1-njac12*nf2)/det
736: u2 = u2-(njac11*nf2-njac21*nf1)/det
737: else
738: exitloop = .true.
739: endif
740: k=k+1
741: enddo
743: boundary_v(1+i) = u1*u1-u2*u2
744: if ((j .eq. 0) .or. (j .eq. 1)) then
745: xt = xt + hx
746: else
747: yt = yt + hy
748: endif
750: enddo
752: if (j.eq.0) then
753: PetscCall(VecRestoreArrayF90(Bottom,boundary_v,ierr))
754: elseif (j.eq.1) then
755: PetscCall(VecRestoreArrayF90(Top,boundary_v,ierr))
756: elseif (j.eq.2) then
757: PetscCall(VecRestoreArrayF90(Left,boundary_v,ierr))
758: elseif (j.eq.3) then
759: PetscCall(VecRestoreArrayF90(Right,boundary_v,ierr))
760: endif
762: enddo
764: ! Scale the boundary if desired
765: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-bottom',scl,flg,ierr))
766: if (flg .eqv. PETSC_TRUE) then
767: PetscCall(VecScale(Bottom,scl,ierr))
768: endif
770: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-top',scl,flg,ierr))
771: if (flg .eqv. PETSC_TRUE) then
772: PetscCall(VecScale(Top,scl,ierr))
773: endif
775: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-right',scl,flg,ierr))
776: if (flg .eqv. PETSC_TRUE) then
777: PetscCall(VecScale(Right,scl,ierr))
778: endif
780: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-left',scl,flg,ierr))
781: if (flg .eqv. PETSC_TRUE) then
782: PetscCall(VecScale(Left,scl,ierr))
783: endif
785: end
787: ! ----------------------------------------------------------------------------
788: !
789: !/*
790: ! MSA_Plate - Calculates an obstacle for surface to stretch over
791: !
792: ! Output Parameter:
793: !. xl - lower bound vector
794: !. xu - upper bound vector
795: !
796: !*/
798: subroutine MSA_Plate(tao,xl,xu,dummy,ierr)
799: use plate2fmodule
800: implicit none
802: Tao tao
803: Vec xl,xu
804: PetscErrorCode ierr
805: PetscInt i,j,row
806: PetscInt xs, xm, ys, ym
807: PetscReal lb,ub
808: PetscInt dummy
809: PetscReal, pointer :: xl_v(:)
811: lb = PETSC_NINFINITY
812: ub = PETSC_INFINITY
814: if (bmy .lt. 0) bmy = 0
815: if (bmy .gt. my) bmy = my
816: if (bmx .lt. 0) bmx = 0
817: if (bmx .gt. mx) bmx = mx
819: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
821: PetscCall(VecSet(xl,lb,ierr))
822: PetscCall(VecSet(xu,ub,ierr))
824: PetscCall(VecGetArrayF90(xl,xl_v,ierr))
826: do i=xs,xs+xm-1
828: do j=ys,ys+ym-1
830: row=(j-ys)*xm + (i-xs)
832: if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
833: & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
834: xl_v(1+row) = bheight
836: endif
838: enddo
839: enddo
841: PetscCall(VecRestoreArrayF90(xl,xl_v,ierr))
843: end
845: ! ----------------------------------------------------------------------------
846: !
847: !/*
848: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
849: !
850: ! Output Parameter:
851: !. X - vector for initial guess
852: !
853: !*/
855: subroutine MSA_InitialPoint(X, ierr)
856: use plate2fmodule
857: implicit none
859: Vec X
860: PetscErrorCode ierr
861: PetscInt start,i,j
862: PetscInt row
863: PetscInt xs,xm,gxs,gxm
864: PetscInt ys,ym,gys,gym
865: PetscReal zero, np5
867: PetscReal,pointer :: left_v(:),right_v(:)
868: PetscReal,pointer :: bottom_v(:),top_v(:)
869: PetscReal,pointer :: x_v(:)
870: PetscBool flg
871: PetscRandom rctx
873: zero = 0.0
874: np5 = -0.5
876: PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-start', start,flg,ierr))
878: if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
879: PetscCall(VecSet(X,zero,ierr))
881: elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
882: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr))
883: do i=0,start-1
884: PetscCall(VecSetRandom(X,rctx,ierr))
885: enddo
887: PetscCall(PetscRandomDestroy(rctx,ierr))
888: PetscCall(VecShift(X,np5,ierr))
890: else ! average of boundary conditions
892: ! Get Local mesh boundaries
893: PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr))
894: PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr))
896: ! Get pointers to vector data
897: PetscCall(VecGetArrayF90(Top,top_v,ierr))
898: PetscCall(VecGetArrayF90(Bottom,bottom_v,ierr))
899: PetscCall(VecGetArrayF90(Left,left_v,ierr))
900: PetscCall(VecGetArrayF90(Right,right_v,ierr))
902: PetscCall(VecGetArrayF90(localX,x_v,ierr))
904: ! Perform local computations
905: do j=ys,ys+ym-1
906: do i=xs,xs+xm-1
907: row = (j-gys)*gxm + (i-gxs)
908: x_v(1+row) = ((j+1)*bottom_v(1+i-xs+1)/my + (my-j+1)*top_v(1+i-xs+1)/(my+2) + &
909: & (i+1)*left_v(1+j-ys+1)/mx + (mx-i+1)*right_v(1+j-ys+1)/(mx+2))*0.5
910: enddo
911: enddo
913: ! Restore vectors
914: PetscCall(VecRestoreArrayF90(localX,x_v,ierr))
916: PetscCall(VecRestoreArrayF90(Left,left_v,ierr))
917: PetscCall(VecRestoreArrayF90(Top,top_v,ierr))
918: PetscCall(VecRestoreArrayF90(Bottom,bottom_v,ierr))
919: PetscCall(VecRestoreArrayF90(Right,right_v,ierr))
921: PetscCall(DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr))
922: PetscCall(DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr))
924: endif
926: end
928: !
929: !/*TEST
930: !
931: ! build:
932: ! requires: !complex
933: !
934: ! test:
935: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
936: ! filter: sort -b
937: ! filter_output: sort -b
938: ! requires: !single
939: !
940: ! test:
941: ! suffix: 2
942: ! nsize: 2
943: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
944: ! filter: sort -b
945: ! filter_output: sort -b
946: ! requires: !single
947: !
948: !TEST*/