Actual source code: plate2f.F90
1: ! Program usage: mpiexec -n <proc> plate2f [all TAO options]
2: !
3: ! This example demonstrates use of the TAO package to solve a bound constrained
4: ! minimization problem. This example is based on a problem from the
5: ! MINPACK-2 test suite. Given a rectangular 2-D domain and boundary values
6: ! along the edges of the domain, the objective is to find the surface
7: ! with the minimal area that satisfies the boundary conditions.
8: ! The command line options are:
9: ! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
10: ! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
11: ! -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
12: ! -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
13: ! -bheight <ht>, where <ht> = height of the plate
14: !
15: #include "petsc/finclude/petscdmda.h"
16: #include "petsc/finclude/petsctao.h"
17: module plate2fmodule
18: use petscdmda
19: use petsctao
21: implicit none
22: Vec localX, localV
23: Vec Top, Left
24: Vec Right, Bottom
25: DM dm
26: PetscReal bheight
27: PetscInt bmx, bmy
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: ! Variable declarations
30: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: !
32: ! Variables:
33: ! (common from plate2f.h):
34: ! Nx, Ny number of processors in x- and y- directions
35: ! mx, my number of grid points in x,y directions
36: ! N global dimension of vector
37: PetscInt mx, my, Nx, Ny, N
38: contains
40: ! ---------------------------------------------------------------------
41: !
42: ! FormFunctionGradient - Evaluates function f(X).
43: !
44: ! Input Parameters:
45: ! tao - the Tao context
46: ! X - the input vector
47: ! dummy - optional user-defined context, as set by TaoSetFunction()
48: ! (not used here)
49: !
50: ! Output Parameters:
51: ! fcn - the newly evaluated function
52: ! G - the gradient vector
53: ! info - error code
54: !
56: subroutine FormFunctionGradient(ta, X, fcn, G, dummy, ierr)
57: ! Input/output variables
59: Tao ta
60: PetscReal fcn
61: Vec X, G
62: PetscErrorCode ierr
63: PetscInt dummy
65: PetscInt i, j, row
66: PetscInt xs, xm
67: PetscInt gxs, gxm
68: PetscInt ys, ym
69: PetscInt gys, gym
70: PetscReal ft, zero, hx, hy, hydhx, hxdhy
71: PetscReal area, rhx, rhy
72: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3
73: PetscReal d4, d5, d6, d7, d8
74: PetscReal df1dxc, df2dxc, df3dxc, df4dxc
75: PetscReal df5dxc, df6dxc
76: PetscReal xc, xl, xr, xt, xb, xlt, xrb
78: PetscReal, pointer :: g_v(:), x_v(:)
79: PetscReal, pointer :: top_v(:), left_v(:)
80: PetscReal, pointer :: right_v(:), bottom_v(:)
82: ft = 0.0
83: zero = 0.0
84: hx = 1.0/real(mx + 1)
85: hy = 1.0/real(my + 1)
86: hydhx = hy/hx
87: hxdhy = hx/hy
88: area = 0.5*hx*hy
89: rhx = real(mx) + 1.0
90: rhy = real(my) + 1.0
92: ! Get local mesh boundaries
93: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
94: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
96: ! Scatter ghost points to local vector
97: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
98: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))
100: ! Initialize the vector to zero
101: PetscCall(VecSet(localV, zero, ierr))
103: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
104: PetscCall(VecGetArray(localX, x_v, ierr))
105: PetscCall(VecGetArray(localV, g_v, ierr))
106: PetscCall(VecGetArray(Top, top_v, ierr))
107: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
108: PetscCall(VecGetArray(Left, left_v, ierr))
109: PetscCall(VecGetArray(Right, right_v, ierr))
111: ! Compute function over the locally owned part of the mesh
112: do j = ys, ys + ym - 1
113: do i = xs, xs + xm - 1
114: row = (j - gys)*gxm + (i - gxs)
115: xc = x_v(1 + row)
116: xt = xc
117: xb = xc
118: xr = xc
119: xl = xc
120: xrb = xc
121: xlt = xc
123: if (i == 0) then !left side
124: xl = left_v(1 + j - ys + 1)
125: xlt = left_v(1 + j - ys + 2)
126: else
127: xl = x_v(1 + row - 1)
128: end if
130: if (j == 0) then !bottom side
131: xb = bottom_v(1 + i - xs + 1)
132: xrb = bottom_v(1 + i - xs + 2)
133: else
134: xb = x_v(1 + row - gxm)
135: end if
137: if (i + 1 == gxs + gxm) then !right side
138: xr = right_v(1 + j - ys + 1)
139: xrb = right_v(1 + j - ys)
140: else
141: xr = x_v(1 + row + 1)
142: end if
144: if (j + 1 == gys + gym) then !top side
145: xt = top_v(1 + i - xs + 1)
146: xlt = top_v(1 + i - xs)
147: else
148: xt = x_v(1 + row + gxm)
149: end if
151: if ((i > gxs) .and. (j + 1 < gys + gym)) then
152: xlt = x_v(1 + row - 1 + gxm)
153: end if
155: if ((j > gys) .and. (i + 1 < gxs + gxm)) then
156: xrb = x_v(1 + row + 1 - gxm)
157: end if
159: d1 = xc - xl
160: d2 = xc - xr
161: d3 = xc - xt
162: d4 = xc - xb
163: d5 = xr - xrb
164: d6 = xrb - xb
165: d7 = xlt - xl
166: d8 = xt - xlt
168: df1dxc = d1*hydhx
169: df2dxc = d1*hydhx + d4*hxdhy
170: df3dxc = d3*hxdhy
171: df4dxc = d2*hydhx + d3*hxdhy
172: df5dxc = d2*hydhx
173: df6dxc = d4*hxdhy
175: d1 = d1*rhx
176: d2 = d2*rhx
177: d3 = d3*rhy
178: d4 = d4*rhy
179: d5 = d5*rhy
180: d6 = d6*rhx
181: d7 = d7*rhy
182: d8 = d8*rhx
184: f1 = sqrt(1.0 + d1*d1 + d7*d7)
185: f2 = sqrt(1.0 + d1*d1 + d4*d4)
186: f3 = sqrt(1.0 + d3*d3 + d8*d8)
187: f4 = sqrt(1.0 + d3*d3 + d2*d2)
188: f5 = sqrt(1.0 + d2*d2 + d5*d5)
189: f6 = sqrt(1.0 + d4*d4 + d6*d6)
191: ft = ft + f2 + f4
193: df1dxc = df1dxc/f1
194: df2dxc = df2dxc/f2
195: df3dxc = df3dxc/f3
196: df4dxc = df4dxc/f4
197: df5dxc = df5dxc/f5
198: df6dxc = df6dxc/f6
200: g_v(1 + row) = 0.5*(df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc)
201: end do
202: end do
204: ! Compute triangular areas along the border of the domain.
205: if (xs == 0) then ! left side
206: do j = ys, ys + ym - 1
207: d3 = (left_v(1 + j - ys + 1) - left_v(1 + j - ys + 2))*rhy
208: d2 = (left_v(1 + j - ys + 1) - x_v(1 + (j - gys)*gxm))*rhx
209: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
210: end do
211: end if
213: if (ys == 0) then !bottom side
214: do i = xs, xs + xm - 1
215: d2 = (bottom_v(1 + i + 1 - xs) - bottom_v(1 + i - xs + 2))*rhx
216: d3 = (bottom_v(1 + i - xs + 1) - x_v(1 + i - gxs))*rhy
217: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
218: end do
219: end if
221: if (xs + xm == mx) then ! right side
222: do j = ys, ys + ym - 1
223: d1 = (x_v(1 + (j + 1 - gys)*gxm - 1) - right_v(1 + j - ys + 1))*rhx
224: d4 = (right_v(1 + j - ys) - right_v(1 + j - ys + 1))*rhy
225: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
226: end do
227: end if
229: if (ys + ym == my) then
230: do i = xs, xs + xm - 1
231: d1 = (x_v(1 + (gym - 1)*gxm + i - gxs) - top_v(1 + i - xs + 1))*rhy
232: d4 = (top_v(1 + i - xs + 1) - top_v(1 + i - xs))*rhx
233: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
234: end do
235: end if
237: if ((ys == 0) .and. (xs == 0)) then
238: d1 = (left_v(1 + 0) - left_v(1 + 1))*rhy
239: d2 = (bottom_v(1 + 0) - bottom_v(1 + 1))*rhx
240: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
241: end if
243: if ((ys + ym == my) .and. (xs + xm == mx)) then
244: d1 = (right_v(1 + ym + 1) - right_v(1 + ym))*rhy
245: d2 = (top_v(1 + xm + 1) - top_v(1 + xm))*rhx
246: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
247: end if
249: ft = ft*area
250: PetscCallMPI(MPI_Allreduce(ft, fcn, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD, ierr))
252: ! Restore vectors
253: PetscCall(VecRestoreArray(localX, x_v, ierr))
254: PetscCall(VecRestoreArray(localV, g_v, ierr))
255: PetscCall(VecRestoreArray(Left, left_v, ierr))
256: PetscCall(VecRestoreArray(Top, top_v, ierr))
257: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
258: PetscCall(VecRestoreArray(Right, right_v, ierr))
260: ! Scatter values to global vector
261: PetscCall(DMLocalToGlobalBegin(dm, localV, INSERT_VALUES, G, ierr))
262: PetscCall(DMLocalToGlobalEnd(dm, localV, INSERT_VALUES, G, ierr))
264: PetscCall(PetscLogFlops(70.0d0*xm*ym, ierr))
266: end !FormFunctionGradient
268: ! ----------------------------------------------------------------------------
269: !
270: !
271: ! FormHessian - Evaluates Hessian matrix.
272: !
273: ! Input Parameters:
274: !. tao - the Tao context
275: !. X - input vector
276: !. dummy - not used
277: !
278: ! Output Parameters:
279: !. Hessian - Hessian matrix
280: !. Hpc - optionally different matrix used to construct the preconditioner
281: !
282: ! Notes:
283: ! Due to mesh point reordering with DMs, we must always work
284: ! with the local mesh points, and then transform them to the new
285: ! global numbering with the local-to-global mapping. We cannot work
286: ! directly with the global numbers for the original uniprocessor mesh!
287: !
288: ! MatSetValuesLocal(), using the local ordering (including
289: ! ghost points!)
290: ! - Then set matrix entries using the local ordering
291: ! by calling MatSetValuesLocal()
293: subroutine FormHessian(ta, X, Hessian, Hpc, dummy, ierr)
295: Tao ta
296: Vec X
297: Mat Hessian, Hpc
298: PetscInt dummy
299: PetscErrorCode ierr
301: PetscInt i, j, k, row
302: PetscInt xs, xm, gxs, gxm
303: PetscInt ys, ym, gys, gym
304: PetscInt col(0:6)
305: PetscReal hx, hy, hydhx, hxdhy, rhx, rhy
306: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3
307: PetscReal d4, d5, d6, d7, d8
308: PetscReal xc, xl, xr, xt, xb, xlt, xrb
309: PetscReal hl, hr, ht, hb, hc, htl, hbr
311: PetscReal, pointer :: right_v(:), left_v(:)
312: PetscReal, pointer :: bottom_v(:), top_v(:)
313: PetscReal, pointer :: x_v(:)
314: PetscReal v(0:6)
315: PetscBool assembled
316: PetscInt i1
318: i1 = 1
320: ! Set various matrix options
321: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE, ierr))
323: ! Get local mesh boundaries
324: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
325: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
327: ! Scatter ghost points to local vectors
328: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
329: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))
331: ! Get pointers to vector data (see note on Fortran arrays above)
332: PetscCall(VecGetArray(localX, x_v, ierr))
333: PetscCall(VecGetArray(Top, top_v, ierr))
334: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
335: PetscCall(VecGetArray(Left, left_v, ierr))
336: PetscCall(VecGetArray(Right, right_v, ierr))
338: ! Initialize matrix entries to zero
339: PetscCall(MatAssembled(Hessian, assembled, ierr))
340: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian, ierr))
342: rhx = real(mx) + 1.0
343: rhy = real(my) + 1.0
344: hx = 1.0/rhx
345: hy = 1.0/rhy
346: hydhx = hy/hx
347: hxdhy = hx/hy
348: ! compute Hessian over the locally owned part of the mesh
350: do i = xs, xs + xm - 1
351: do j = ys, ys + ym - 1
352: row = (j - gys)*gxm + (i - gxs)
354: xc = x_v(1 + row)
355: xt = xc
356: xb = xc
357: xr = xc
358: xl = xc
359: xrb = xc
360: xlt = xc
362: if (i == gxs) then ! Left side
363: xl = left_v(1 + j - ys + 1)
364: xlt = left_v(1 + j - ys + 2)
365: else
366: xl = x_v(1 + row - 1)
367: end if
369: if (j == gys) then ! bottom side
370: xb = bottom_v(1 + i - xs + 1)
371: xrb = bottom_v(1 + i - xs + 2)
372: else
373: xb = x_v(1 + row - gxm)
374: end if
376: if (i + 1 == gxs + gxm) then !right side
377: xr = right_v(1 + j - ys + 1)
378: xrb = right_v(1 + j - ys)
379: else
380: xr = x_v(1 + row + 1)
381: end if
383: if (j + 1 == gym + gys) then !top side
384: xt = top_v(1 + i - xs + 1)
385: xlt = top_v(1 + i - xs)
386: else
387: xt = x_v(1 + row + gxm)
388: end if
390: if ((i > gxs) .and. (j + 1 < gys + gym)) then
391: xlt = x_v(1 + row - 1 + gxm)
392: end if
394: if ((i + 1 < gxs + gxm) .and. (j > gys)) then
395: xrb = x_v(1 + row + 1 - gxm)
396: end if
398: d1 = (xc - xl)*rhx
399: d2 = (xc - xr)*rhx
400: d3 = (xc - xt)*rhy
401: d4 = (xc - xb)*rhy
402: d5 = (xrb - xr)*rhy
403: d6 = (xrb - xb)*rhx
404: d7 = (xlt - xl)*rhy
405: d8 = (xlt - xt)*rhx
407: f1 = sqrt(1.0 + d1*d1 + d7*d7)
408: f2 = sqrt(1.0 + d1*d1 + d4*d4)
409: f3 = sqrt(1.0 + d3*d3 + d8*d8)
410: f4 = sqrt(1.0 + d3*d3 + d2*d2)
411: f5 = sqrt(1.0 + d2*d2 + d5*d5)
412: f6 = sqrt(1.0 + d4*d4 + d6*d6)
414: hl = (-hydhx*(1.0 + d7*d7) + d1*d7)/(f1*f1*f1) + (-hydhx*(1.0 + d4*d4) + d1*d4)/(f2*f2*f2)
416: hr = (-hydhx*(1.0 + d5*d5) + d2*d5)/(f5*f5*f5) + (-hydhx*(1.0 + d3*d3) + d2*d3)/(f4*f4*f4)
418: ht = (-hxdhy*(1.0 + d8*d8) + d3*d8)/(f3*f3*f3) + (-hxdhy*(1.0 + d2*d2) + d2*d3)/(f4*f4*f4)
420: hb = (-hxdhy*(1.0 + d6*d6) + d4*d6)/(f6*f6*f6) + (-hxdhy*(1.0 + d1*d1) + d1*d4)/(f2*f2*f2)
422: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
423: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
425: 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) + &
426: & 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) + &
427: & hydhx*(1.0 + d3*d3) - 2*d2*d3)/(f4*f4*f4)
429: hl = hl*0.5
430: hr = hr*0.5
431: ht = ht*0.5
432: hb = hb*0.5
433: hbr = hbr*0.5
434: htl = htl*0.5
435: hc = hc*0.5
437: k = 0
439: if (j > 0) then
440: v(k) = hb
441: col(k) = row - gxm
442: k = k + 1
443: end if
445: if ((j > 0) .and. (i < mx - 1)) then
446: v(k) = hbr
447: col(k) = row - gxm + 1
448: k = k + 1
449: end if
451: if (i > 0) then
452: v(k) = hl
453: col(k) = row - 1
454: k = k + 1
455: end if
457: v(k) = hc
458: col(k) = row
459: k = k + 1
461: if (i < mx - 1) then
462: v(k) = hr
463: col(k) = row + 1
464: k = k + 1
465: end if
467: if ((i > 0) .and. (j < my - 1)) then
468: v(k) = htl
469: col(k) = row + gxm - 1
470: k = k + 1
471: end if
473: if (j < my - 1) then
474: v(k) = ht
475: col(k) = row + gxm
476: k = k + 1
477: end if
479: ! Set matrix values using local numbering, defined earlier in main routine
480: PetscCall(MatSetValuesLocal(Hessian, i1, [row], k, col, v, INSERT_VALUES, ierr))
482: end do
483: end do
485: ! restore vectors
486: PetscCall(VecRestoreArray(localX, x_v, ierr))
487: PetscCall(VecRestoreArray(Left, left_v, ierr))
488: PetscCall(VecRestoreArray(Right, right_v, ierr))
489: PetscCall(VecRestoreArray(Top, top_v, ierr))
490: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
492: ! Assemble the matrix
493: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY, ierr))
494: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY, ierr))
496: PetscCall(PetscLogFlops(199.0d0*xm*ym, ierr))
498: end
500: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
502: ! ----------------------------------------------------------------------------
503: !
504: !/*
505: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
506: !
507: !
508: !*/
510: subroutine MSA_BoundaryConditions(ierr)
512: PetscErrorCode ierr
513: PetscInt i, j, k, limit, maxits
514: PetscInt xs, xm, gxs, gxm
515: PetscInt ys, ym, gys, gym
516: PetscInt bsize, lsize
517: PetscInt tsize, rsize
518: PetscReal one, two, three, tol
519: PetscReal scl, fnorm, det, xt
520: PetscReal yt, hx, hy, u1, u2, nf1, nf2
521: PetscReal njac11, njac12, njac21, njac22
522: PetscReal b, t, l, r
523: PetscReal, pointer :: boundary_v(:)
525: logical exitloop
526: PetscBool flg
528: limit = 0
529: maxits = 5
530: tol = 1e-10
531: b = -0.5
532: t = 0.5
533: l = -0.5
534: r = 0.5
535: xt = 0
536: yt = 0
537: one = 1.0
538: two = 2.0
539: three = 3.0
541: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
542: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
544: bsize = xm + 2
545: lsize = ym + 2
546: rsize = ym + 2
547: tsize = xm + 2
549: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, bsize, PETSC_DECIDE, Bottom, ierr))
550: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, tsize, PETSC_DECIDE, Top, ierr))
551: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, lsize, PETSC_DECIDE, Left, ierr))
552: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, rsize, PETSC_DECIDE, Right, ierr))
554: hx = (r - l)/(mx + 1)
555: hy = (t - b)/(my + 1)
557: do j = 0, 3
559: if (j == 0) then
560: yt = b
561: xt = l + hx*xs
562: limit = bsize
563: PetscCall(VecGetArray(Bottom, boundary_v, ierr))
565: elseif (j == 1) then
566: yt = t
567: xt = l + hx*xs
568: limit = tsize
569: PetscCall(VecGetArray(Top, boundary_v, ierr))
571: elseif (j == 2) then
572: yt = b + hy*ys
573: xt = l
574: limit = lsize
575: PetscCall(VecGetArray(Left, boundary_v, ierr))
577: elseif (j == 3) then
578: yt = b + hy*ys
579: xt = r
580: limit = rsize
581: PetscCall(VecGetArray(Right, boundary_v, ierr))
582: end if
584: do i = 0, limit - 1
586: u1 = xt
587: u2 = -yt
588: k = 0
589: exitloop = .false.
590: do while (k < maxits .and. (.not. exitloop))
592: nf1 = u1 + u1*u2*u2 - u1*u1*u1/three - xt
593: nf2 = -u2 - u1*u1*u2 + u2*u2*u2/three - yt
594: fnorm = sqrt(nf1*nf1 + nf2*nf2)
595: if (fnorm > tol) then
596: njac11 = one + u2*u2 - u1*u1
597: njac12 = two*u1*u2
598: njac21 = -two*u1*u2
599: njac22 = -one - u1*u1 + u2*u2
600: det = njac11*njac22 - njac21*njac12
601: u1 = u1 - (njac22*nf1 - njac12*nf2)/det
602: u2 = u2 - (njac11*nf2 - njac21*nf1)/det
603: else
604: exitloop = .true.
605: end if
606: k = k + 1
607: end do
609: boundary_v(1 + i) = u1*u1 - u2*u2
610: if ((j == 0) .or. (j == 1)) then
611: xt = xt + hx
612: else
613: yt = yt + hy
614: end if
616: end do
618: if (j == 0) then
619: PetscCall(VecRestoreArray(Bottom, boundary_v, ierr))
620: elseif (j == 1) then
621: PetscCall(VecRestoreArray(Top, boundary_v, ierr))
622: elseif (j == 2) then
623: PetscCall(VecRestoreArray(Left, boundary_v, ierr))
624: elseif (j == 3) then
625: PetscCall(VecRestoreArray(Right, boundary_v, ierr))
626: end if
628: end do
630: ! Scale the boundary if desired
631: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bottom', scl, flg, ierr))
632: if (flg .eqv. PETSC_TRUE) then
633: PetscCall(VecScale(Bottom, scl, ierr))
634: end if
636: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-top', scl, flg, ierr))
637: if (flg .eqv. PETSC_TRUE) then
638: PetscCall(VecScale(Top, scl, ierr))
639: end if
641: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-right', scl, flg, ierr))
642: if (flg .eqv. PETSC_TRUE) then
643: PetscCall(VecScale(Right, scl, ierr))
644: end if
646: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-left', scl, flg, ierr))
647: if (flg .eqv. PETSC_TRUE) then
648: PetscCall(VecScale(Left, scl, ierr))
649: end if
651: end
653: ! ----------------------------------------------------------------------------
654: !
655: !/*
656: ! MSA_Plate - Calculates an obstacle for surface to stretch over
657: !
658: ! Output Parameter:
659: !. xl - lower bound vector
660: !. xu - upper bound vector
661: !
662: !*/
664: subroutine MSA_Plate(ta, xl, xu, dummy, ierr)
666: Tao ta
667: Vec xl, xu
668: PetscErrorCode ierr
669: PetscInt i, j, row
670: PetscInt xs, xm, ys, ym
671: PetscReal lb, ub
672: PetscInt dummy
673: PetscReal, pointer :: xl_v(:)
675: lb = PETSC_NINFINITY
676: ub = PETSC_INFINITY
678: if (bmy < 0) bmy = 0
679: if (bmy > my) bmy = my
680: if (bmx < 0) bmx = 0
681: if (bmx > mx) bmx = mx
683: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
685: PetscCall(VecSet(xl, lb, ierr))
686: PetscCall(VecSet(xu, ub, ierr))
688: PetscCall(VecGetArray(xl, xl_v, ierr))
690: do i = xs, xs + xm - 1
692: do j = ys, ys + ym - 1
694: row = (j - ys)*xm + (i - xs)
696: if (i >= ((mx - bmx)/2) .and. i < (mx - (mx - bmx)/2) .and. &
697: & j >= ((my - bmy)/2) .and. j < (my - (my - bmy)/2)) then
698: xl_v(1 + row) = bheight
700: end if
702: end do
703: end do
705: PetscCall(VecRestoreArray(xl, xl_v, ierr))
707: end
709: ! ----------------------------------------------------------------------------
710: !
711: !/*
712: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
713: !
714: ! Output Parameter:
715: !. X - vector for initial guess
716: !
717: !*/
719: subroutine MSA_InitialPoint(X, ierr)
721: Vec X
722: PetscErrorCode ierr
723: PetscInt start, i, j
724: PetscInt row
725: PetscInt xs, xm, gxs, gxm
726: PetscInt ys, ym, gys, gym
727: PetscReal zero, np5
729: PetscReal, pointer :: left_v(:), right_v(:)
730: PetscReal, pointer :: bottom_v(:), top_v(:)
731: PetscReal, pointer :: x_v(:)
732: PetscBool flg
733: PetscRandom rctx
735: zero = 0.0
736: np5 = -0.5
738: PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-start', start, flg, ierr))
740: if ((flg .eqv. PETSC_TRUE) .and. (start == 0)) then ! the zero vector is reasonable
741: PetscCall(VecSet(X, zero, ierr))
743: elseif ((flg .eqv. PETSC_TRUE) .and. (start > 0)) then ! random start -0.5 < xi < 0.5
744: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, rctx, ierr))
745: do i = 0, start - 1
746: PetscCall(VecSetRandom(X, rctx, ierr))
747: end do
749: PetscCall(PetscRandomDestroy(rctx, ierr))
750: PetscCall(VecShift(X, np5, ierr))
752: else ! average of boundary conditions
754: ! Get Local mesh boundaries
755: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
756: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
758: ! Get pointers to vector data
759: PetscCall(VecGetArray(Top, top_v, ierr))
760: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
761: PetscCall(VecGetArray(Left, left_v, ierr))
762: PetscCall(VecGetArray(Right, right_v, ierr))
764: PetscCall(VecGetArray(localX, x_v, ierr))
766: ! Perform local computations
767: do j = ys, ys + ym - 1
768: do i = xs, xs + xm - 1
769: row = (j - gys)*gxm + (i - gxs)
770: x_v(1 + row) = ((j + 1)*bottom_v(1 + i - xs + 1)/my + (my - j + 1)*top_v(1 + i - xs + 1)/(my + 2) + &
771: & (i + 1)*left_v(1 + j - ys + 1)/mx + (mx - i + 1)*right_v(1 + j - ys + 1)/(mx + 2))*0.5
772: end do
773: end do
775: ! Restore vectors
776: PetscCall(VecRestoreArray(localX, x_v, ierr))
778: PetscCall(VecRestoreArray(Left, left_v, ierr))
779: PetscCall(VecRestoreArray(Top, top_v, ierr))
780: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
781: PetscCall(VecRestoreArray(Right, right_v, ierr))
783: PetscCall(DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, X, ierr))
784: PetscCall(DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, X, ierr))
786: end if
788: end
790: end module
792: program plate2f
793: use plate2fmodule
794: implicit none
796: PetscErrorCode ierr ! used to check for functions returning nonzeros
797: Vec x ! solution vector
798: PetscInt m ! number of local elements in vector
799: Tao ta ! Tao solver context
800: Mat H ! Hessian matrix
801: ISLocalToGlobalMapping isltog ! local to global mapping object
802: PetscBool flg
803: PetscInt i1, i3, i7
805: ! Initialize Tao
807: i1 = 1
808: i3 = 3
809: i7 = 7
811: PetscCallA(PetscInitialize(ierr))
813: ! Specify default dimensions of the problem
814: mx = 10
815: my = 10
816: bheight = 0.1
818: ! Check for any command line arguments that override defaults
820: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
821: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))
823: bmx = mx/2
824: bmy = my/2
826: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmx', bmx, flg, ierr))
827: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmy', bmy, flg, ierr))
828: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bheight', bheight, flg, ierr))
830: ! Calculate any derived values from parameters
831: N = mx*my
833: ! Let PETSc determine the dimensions of the local vectors
834: Nx = PETSC_DECIDE
835: NY = PETSC_DECIDE
837: ! A two dimensional distributed array will help define this problem, which
838: ! derives from an elliptic PDE on a two-dimensional domain. From the
839: ! distributed array, create the vectors
841: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, Nx, Ny, i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
842: PetscCallA(DMSetFromOptions(dm, ierr))
843: PetscCallA(DMSetUp(dm, ierr))
845: ! Extract global and local vectors from DM; The local vectors are
846: ! used solely as work space for the evaluation of the function,
847: ! gradient, and Hessian. Duplicate for remaining vectors that are
848: ! the same types.
850: PetscCallA(DMCreateGlobalVector(dm, x, ierr))
851: PetscCallA(DMCreateLocalVector(dm, localX, ierr))
852: PetscCallA(VecDuplicate(localX, localV, ierr))
854: ! Create a matrix data structure to store the Hessian.
855: ! Here we (optionally) also associate the local numbering scheme
856: ! with the matrix so that later we can use local indices for matrix
857: ! assembly
859: PetscCallA(VecGetLocalSize(x, m, ierr))
860: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, m, m, N, N, i7, PETSC_NULL_INTEGER_ARRAY, i3, PETSC_NULL_INTEGER_ARRAY, H, ierr))
862: PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))
863: PetscCallA(DMGetLocalToGlobalMapping(dm, isltog, ierr))
864: PetscCallA(MatSetLocalToGlobalMapping(H, isltog, isltog, ierr))
866: ! The Tao code begins here
867: ! Create TAO solver and set desired solution method.
868: ! This problems uses bounded variables, so the
869: ! method must either be 'tao_tron' or 'tao_blmvm'
871: PetscCallA(TaoCreate(PETSC_COMM_WORLD, ta, ierr))
872: PetscCallA(TaoSetType(ta, TAOBLMVM, ierr))
874: ! Set minimization function and gradient, hessian evaluation functions
876: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
878: PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))
880: ! Set Variable bounds
881: PetscCallA(MSA_BoundaryConditions(ierr))
882: PetscCallA(TaoSetVariableBoundsRoutine(ta, MSA_Plate, 0, ierr))
884: ! Set the initial solution guess
885: PetscCallA(MSA_InitialPoint(x, ierr))
886: PetscCallA(TaoSetSolution(ta, x, ierr))
888: ! Check for any tao command line options
889: PetscCallA(TaoSetFromOptions(ta, ierr))
891: ! Solve the application
892: PetscCallA(TaoSolve(ta, ierr))
894: ! Free TAO data structures
895: PetscCallA(TaoDestroy(ta, ierr))
897: ! Free PETSc data structures
898: PetscCallA(VecDestroy(x, ierr))
899: PetscCallA(VecDestroy(Top, ierr))
900: PetscCallA(VecDestroy(Bottom, ierr))
901: PetscCallA(VecDestroy(Left, ierr))
902: PetscCallA(VecDestroy(Right, ierr))
903: PetscCallA(MatDestroy(H, ierr))
904: PetscCallA(VecDestroy(localX, ierr))
905: PetscCallA(VecDestroy(localV, ierr))
906: PetscCallA(DMDestroy(dm, ierr))
908: ! Finalize TAO
910: PetscCallA(PetscFinalize(ierr))
912: end program plate2f
913: !
914: !/*TEST
915: !
916: ! build:
917: ! requires: !complex
918: !
919: ! test:
920: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
921: ! filter: sort -b
922: ! filter_output: sort -b
923: ! requires: !single
924: !
925: ! test:
926: ! suffix: 2
927: ! nsize: 2
928: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
929: ! filter: sort -b
930: ! filter_output: sort -b
931: ! requires: !single
932: !
933: !TEST*/