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, intent(out) :: 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, 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: hx = 1.0/real(mx + 1)
83: hy = 1.0/real(my + 1)
84: hydhx = hy/hx
85: hxdhy = hx/hy
86: area = 0.5*hx*hy
87: rhx = real(mx) + 1.0
88: rhy = real(my) + 1.0
90: ! Get local mesh boundaries
91: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
92: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
94: ! Scatter ghost points to local vector
95: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
96: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))
98: ! Initialize the vector to zero
99: PetscCall(VecSet(localV, 0.0_PETSC_REAL_KIND, ierr))
101: ! Get arrays to vector data (See note above about using VecGetArray in Fortran)
102: PetscCall(VecGetArray(localX, x_v, ierr))
103: PetscCall(VecGetArray(localV, g_v, ierr))
104: PetscCall(VecGetArray(Top, top_v, ierr))
105: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
106: PetscCall(VecGetArray(Left, left_v, ierr))
107: PetscCall(VecGetArray(Right, right_v, ierr))
109: ft = 0.0
110: ! Compute function over the locally owned part of the mesh
111: do j = ys, ys + ym - 1
112: do i = xs, xs + xm - 1
113: row = (j - gys)*gxm + (i - gxs)
114: xc = x_v(1 + row)
115: xt = xc
116: xb = xc
117: xr = xc
118: xl = xc
119: xrb = xc
120: xlt = xc
122: if (i == 0) then !left side
123: xl = left_v(1 + j - ys + 1)
124: xlt = left_v(1 + j - ys + 2)
125: else
126: xl = x_v(1 + row - 1)
127: end if
129: if (j == 0) then !bottom side
130: xb = bottom_v(1 + i - xs + 1)
131: xrb = bottom_v(1 + i - xs + 2)
132: else
133: xb = x_v(1 + row - gxm)
134: end if
136: if (i + 1 == gxs + gxm) then !right side
137: xr = right_v(1 + j - ys + 1)
138: xrb = right_v(1 + j - ys)
139: else
140: xr = x_v(1 + row + 1)
141: end if
143: if (j + 1 == gys + gym) then !top side
144: xt = top_v(1 + i - xs + 1)
145: xlt = top_v(1 + i - xs)
146: else
147: xt = x_v(1 + row + gxm)
148: end if
150: if ((i > gxs) .and. (j + 1 < gys + gym)) then
151: xlt = x_v(1 + row - 1 + gxm)
152: end if
154: if ((j > gys) .and. (i + 1 < gxs + gxm)) then
155: xrb = x_v(1 + row + 1 - gxm)
156: end if
158: d1 = xc - xl
159: d2 = xc - xr
160: d3 = xc - xt
161: d4 = xc - xb
162: d5 = xr - xrb
163: d6 = xrb - xb
164: d7 = xlt - xl
165: d8 = xt - xlt
167: df1dxc = d1*hydhx
168: df2dxc = d1*hydhx + d4*hxdhy
169: df3dxc = d3*hxdhy
170: df4dxc = d2*hydhx + d3*hxdhy
171: df5dxc = d2*hydhx
172: df6dxc = d4*hxdhy
174: d1 = d1*rhx
175: d2 = d2*rhx
176: d3 = d3*rhy
177: d4 = d4*rhy
178: d5 = d5*rhy
179: d6 = d6*rhx
180: d7 = d7*rhy
181: d8 = d8*rhx
183: f1 = sqrt(1.0 + d1*d1 + d7*d7)
184: f2 = sqrt(1.0 + d1*d1 + d4*d4)
185: f3 = sqrt(1.0 + d3*d3 + d8*d8)
186: f4 = sqrt(1.0 + d3*d3 + d2*d2)
187: f5 = sqrt(1.0 + d2*d2 + d5*d5)
188: f6 = sqrt(1.0 + d4*d4 + d6*d6)
190: ft = ft + f2 + f4
192: df1dxc = df1dxc/f1
193: df2dxc = df2dxc/f2
194: df3dxc = df3dxc/f3
195: df4dxc = df4dxc/f4
196: df5dxc = df5dxc/f5
197: df6dxc = df6dxc/f6
199: g_v(1 + row) = 0.5*(df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc)
200: end do
201: end do
203: ! Compute triangular areas along the border of the domain.
204: if (xs == 0) then ! left side
205: do j = ys, ys + ym - 1
206: d3 = (left_v(1 + j - ys + 1) - left_v(1 + j - ys + 2))*rhy
207: d2 = (left_v(1 + j - ys + 1) - x_v(1 + (j - gys)*gxm))*rhx
208: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
209: end do
210: end if
212: if (ys == 0) then !bottom side
213: do i = xs, xs + xm - 1
214: d2 = (bottom_v(1 + i + 1 - xs) - bottom_v(1 + i - xs + 2))*rhx
215: d3 = (bottom_v(1 + i - xs + 1) - x_v(1 + i - gxs))*rhy
216: ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
217: end do
218: end if
220: if (xs + xm == mx) then ! right side
221: do j = ys, ys + ym - 1
222: d1 = (x_v(1 + (j + 1 - gys)*gxm - 1) - right_v(1 + j - ys + 1))*rhx
223: d4 = (right_v(1 + j - ys) - right_v(1 + j - ys + 1))*rhy
224: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
225: end do
226: end if
228: if (ys + ym == my) then
229: do i = xs, xs + xm - 1
230: d1 = (x_v(1 + (gym - 1)*gxm + i - gxs) - top_v(1 + i - xs + 1))*rhy
231: d4 = (top_v(1 + i - xs + 1) - top_v(1 + i - xs))*rhx
232: ft = ft + sqrt(1.0 + d1*d1 + d4*d4)
233: end do
234: end if
236: if ((ys == 0) .and. (xs == 0)) then
237: d1 = (left_v(1 + 0) - left_v(1 + 1))*rhy
238: d2 = (bottom_v(1 + 0) - bottom_v(1 + 1))*rhx
239: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
240: end if
242: if ((ys + ym == my) .and. (xs + xm == mx)) then
243: d1 = (right_v(1 + ym + 1) - right_v(1 + ym))*rhy
244: d2 = (top_v(1 + xm + 1) - top_v(1 + xm))*rhx
245: ft = ft + sqrt(1.0 + d1*d1 + d2*d2)
246: end if
248: ft = ft*area
249: PetscCallMPI(MPI_Allreduce(ft, fcn, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD, ierr))
251: ! Restore vectors
252: PetscCall(VecRestoreArray(localX, x_v, ierr))
253: PetscCall(VecRestoreArray(localV, g_v, ierr))
254: PetscCall(VecRestoreArray(Left, left_v, ierr))
255: PetscCall(VecRestoreArray(Top, top_v, ierr))
256: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
257: PetscCall(VecRestoreArray(Right, right_v, ierr))
259: ! Scatter values to global vector
260: PetscCall(DMLocalToGlobalBegin(dm, localV, INSERT_VALUES, G, ierr))
261: PetscCall(DMLocalToGlobalEnd(dm, localV, INSERT_VALUES, G, ierr))
263: PetscCall(PetscLogFlops(70.0_C_DOUBLE*xm*ym, ierr))
265: end !FormFunctionGradient
267: ! ----------------------------------------------------------------------------
268: !
269: !
270: ! FormHessian - Evaluates Hessian matrix.
271: !
272: ! Input Parameters:
273: !. tao - the Tao context
274: !. X - input vector
275: !. dummy - not used
276: !
277: ! Output Parameters:
278: !. Hessian - Hessian matrix
279: !. Hpc - optionally different matrix used to construct the preconditioner
280: !
281: ! Notes:
282: ! Due to mesh point reordering with DMs, we must always work
283: ! with the local mesh points, and then transform them to the new
284: ! global numbering with the local-to-global mapping. We cannot work
285: ! directly with the global numbers for the original uniprocessor mesh!
286: !
287: ! MatSetValuesLocal(), using the local ordering (including
288: ! ghost points!)
289: ! - Then set matrix entries using the local ordering
290: ! by calling MatSetValuesLocal()
292: subroutine FormHessian(ta, X, Hessian, Hpc, dummy, ierr)
294: Tao ta
295: Vec X
296: Mat Hessian, Hpc
297: PetscInt dummy
298: PetscErrorCode, intent(out) :: ierr
300: PetscInt i, j, k, row
301: PetscInt xs, xm, gxs, gxm
302: PetscInt ys, ym, gys, gym
303: PetscInt col(0:6)
304: PetscReal hx, hy, hydhx, hxdhy, rhx, rhy
305: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3
306: PetscReal d4, d5, d6, d7, d8
307: PetscReal xc, xl, xr, xt, xb, xlt, xrb
308: PetscReal hl, hr, ht, hb, hc, htl, hbr
310: PetscReal, pointer :: right_v(:), left_v(:)
311: PetscReal, pointer :: bottom_v(:), top_v(:)
312: PetscReal, pointer :: x_v(:)
313: PetscReal v(0:6)
314: PetscBool assembled
316: ! Set various matrix options
317: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE, ierr))
319: ! Get local mesh boundaries
320: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
321: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
323: ! Scatter ghost points to local vectors
324: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
325: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))
327: ! Get pointers to vector data (see note on Fortran arrays above)
328: PetscCall(VecGetArray(localX, x_v, ierr))
329: PetscCall(VecGetArray(Top, top_v, ierr))
330: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
331: PetscCall(VecGetArray(Left, left_v, ierr))
332: PetscCall(VecGetArray(Right, right_v, ierr))
334: ! Initialize matrix entries to zero
335: PetscCall(MatAssembled(Hessian, assembled, ierr))
336: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian, ierr))
338: rhx = real(mx) + 1.0
339: rhy = real(my) + 1.0
340: hx = 1.0/rhx
341: hy = 1.0/rhy
342: hydhx = hy/hx
343: hxdhy = hx/hy
344: ! compute Hessian over the locally owned part of the mesh
346: do i = xs, xs + xm - 1
347: do j = ys, ys + ym - 1
348: row = (j - gys)*gxm + (i - gxs)
350: xc = x_v(1 + row)
351: xt = xc
352: xb = xc
353: xr = xc
354: xl = xc
355: xrb = xc
356: xlt = xc
358: if (i == gxs) then ! Left side
359: xl = left_v(1 + j - ys + 1)
360: xlt = left_v(1 + j - ys + 2)
361: else
362: xl = x_v(1 + row - 1)
363: end if
365: if (j == gys) then ! bottom side
366: xb = bottom_v(1 + i - xs + 1)
367: xrb = bottom_v(1 + i - xs + 2)
368: else
369: xb = x_v(1 + row - gxm)
370: end if
372: if (i + 1 == gxs + gxm) then !right side
373: xr = right_v(1 + j - ys + 1)
374: xrb = right_v(1 + j - ys)
375: else
376: xr = x_v(1 + row + 1)
377: end if
379: if (j + 1 == gym + gys) then !top side
380: xt = top_v(1 + i - xs + 1)
381: xlt = top_v(1 + i - xs)
382: else
383: xt = x_v(1 + row + gxm)
384: end if
386: if ((i > gxs) .and. (j + 1 < gys + gym)) then
387: xlt = x_v(1 + row - 1 + gxm)
388: end if
390: if ((i + 1 < gxs + gxm) .and. (j > gys)) then
391: xrb = x_v(1 + row + 1 - gxm)
392: end if
394: d1 = (xc - xl)*rhx
395: d2 = (xc - xr)*rhx
396: d3 = (xc - xt)*rhy
397: d4 = (xc - xb)*rhy
398: d5 = (xrb - xr)*rhy
399: d6 = (xrb - xb)*rhx
400: d7 = (xlt - xl)*rhy
401: d8 = (xlt - xt)*rhx
403: f1 = sqrt(1.0 + d1**2 + d7**2)
404: f2 = sqrt(1.0 + d1**2 + d4**2)
405: f3 = sqrt(1.0 + d3**2 + d8**2)
406: f4 = sqrt(1.0 + d3**2 + d2**2)
407: f5 = sqrt(1.0 + d2**2 + d5**2)
408: f6 = sqrt(1.0 + d4**2 + d6**2)
410: hl = (-hydhx*(1.0 + d7**2) + d1*d7)/f1**3 + (-hydhx*(1.0 + d4**2) + d1*d4)/f2**3
412: hr = (-hydhx*(1.0 + d5**2) + d2*d5)/f5**3 + (-hydhx*(1.0 + d3**2) + d2*d3)/f4**3
414: ht = (-hxdhy*(1.0 + d8**2) + d3*d8)/f3**3 + (-hxdhy*(1.0 + d2**2) + d2*d3)/f4**3
416: hb = (-hxdhy*(1.0 + d6**2) + d4*d6)/f6**3 + (-hxdhy*(1.0 + d1**2) + d1*d4)/f2**3
418: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
419: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
421: hc = hydhx*(1.0 + d7**2)/f1**3 + hxdhy*(1.0 + d8**2)/f3**3 + hydhx*(1.0 + d5**2)/f5**3 + hxdhy*(1.0 + d6**2)/f6**3 + &
422: (hxdhy*(1.0 + d1**2) + hydhx*(1.0 + d4**2) - 2*d1*d4)/f2**3 + &
423: (hxdhy*(1.0 + d2**2) + hydhx*(1.0 + d3**2) - 2*d2*d3)/f4**3
425: hl = hl*0.5
426: hr = hr*0.5
427: ht = ht*0.5
428: hb = hb*0.5
429: hbr = hbr*0.5
430: htl = htl*0.5
431: hc = hc*0.5
433: k = 0
435: if (j > 0) then
436: v(k) = hb
437: col(k) = row - gxm
438: k = k + 1
439: end if
441: if ((j > 0) .and. (i < mx - 1)) then
442: v(k) = hbr
443: col(k) = row - gxm + 1
444: k = k + 1
445: end if
447: if (i > 0) then
448: v(k) = hl
449: col(k) = row - 1
450: k = k + 1
451: end if
453: v(k) = hc
454: col(k) = row
455: k = k + 1
457: if (i < mx - 1) then
458: v(k) = hr
459: col(k) = row + 1
460: k = k + 1
461: end if
463: if ((i > 0) .and. (j < my - 1)) then
464: v(k) = htl
465: col(k) = row + gxm - 1
466: k = k + 1
467: end if
469: if (j < my - 1) then
470: v(k) = ht
471: col(k) = row + gxm
472: k = k + 1
473: end if
475: ! Set matrix values using local numbering, defined earlier in main routine
476: PetscCall(MatSetValuesLocal(Hessian, 1_PETSC_INT_KIND, [row], k, col, v, INSERT_VALUES, ierr))
478: end do
479: end do
481: ! restore vectors
482: PetscCall(VecRestoreArray(localX, x_v, ierr))
483: PetscCall(VecRestoreArray(Left, left_v, ierr))
484: PetscCall(VecRestoreArray(Right, right_v, ierr))
485: PetscCall(VecRestoreArray(Top, top_v, ierr))
486: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
488: ! Assemble the matrix
489: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY, ierr))
490: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY, ierr))
492: PetscCall(PetscLogFlops(199.0_C_DOUBLE*xm*ym, ierr))
494: end
496: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H
498: ! ----------------------------------------------------------------------------
499: !
500: !/*
501: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
502: !
503: !
504: !*/
506: subroutine MSA_BoundaryConditions(ierr)
508: PetscErrorCode, intent(out) :: ierr
509: PetscInt, parameter :: maxits = 5
510: PetscInt i, j, k, limit
511: PetscInt xs, xm, gxs, gxm
512: PetscInt ys, ym, gys, gym
513: PetscInt bsize, lsize
514: PetscInt tsize, rsize
515: PetscReal scl, fnorm, det, xt
516: PetscReal yt, hx, hy, u1, u2, nf1, nf2
517: PetscReal njac11, njac12, njac21, njac22
518: PetscReal, parameter :: b = -0.5, t = 0.5, l = -0.5, r = 0.5, tol = 1.e-10
519: PetscReal, pointer :: boundary_v(:)
521: logical exitloop
522: PetscBool flg
524: limit = 0
526: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
527: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
529: bsize = xm + 2
530: lsize = ym + 2
531: rsize = ym + 2
532: tsize = xm + 2
534: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, bsize, PETSC_DECIDE, Bottom, ierr))
535: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, tsize, PETSC_DECIDE, Top, ierr))
536: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, lsize, PETSC_DECIDE, Left, ierr))
537: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, rsize, PETSC_DECIDE, Right, ierr))
539: hx = (r - l)/(mx + 1)
540: hy = (t - b)/(my + 1)
542: xt = 0
543: yt = 0
544: do j = 0, 3
546: if (j == 0) then
547: yt = b
548: xt = l + hx*xs
549: limit = bsize
550: PetscCall(VecGetArray(Bottom, boundary_v, ierr))
552: elseif (j == 1) then
553: yt = t
554: xt = l + hx*xs
555: limit = tsize
556: PetscCall(VecGetArray(Top, boundary_v, ierr))
558: elseif (j == 2) then
559: yt = b + hy*ys
560: xt = l
561: limit = lsize
562: PetscCall(VecGetArray(Left, boundary_v, ierr))
564: elseif (j == 3) then
565: yt = b + hy*ys
566: xt = r
567: limit = rsize
568: PetscCall(VecGetArray(Right, boundary_v, ierr))
569: end if
571: do i = 0, limit - 1
573: u1 = xt
574: u2 = -yt
575: k = 0
576: exitloop = .false.
577: do while (k < maxits .and. (.not. exitloop))
579: nf1 = u1 + u1*u2*u2 - u1*u1*u1/3.0_PETSC_REAL_KIND - xt
580: nf2 = -u2 - u1*u1*u2 + u2*u2*u2/3.0_PETSC_REAL_KIND - yt
581: fnorm = sqrt(nf1*nf1 + nf2*nf2)
582: if (fnorm > tol) then
583: njac11 = 1.0_PETSC_REAL_KIND + u2*u2 - u1*u1
584: njac12 = 2.0_PETSC_REAL_KIND*u1*u2
585: njac21 = -2.0_PETSC_REAL_KIND*u1*u2
586: njac22 = -1.0_PETSC_REAL_KIND - u1*u1 + u2*u2
587: det = njac11*njac22 - njac21*njac12
588: u1 = u1 - (njac22*nf1 - njac12*nf2)/det
589: u2 = u2 - (njac11*nf2 - njac21*nf1)/det
590: else
591: exitloop = .true.
592: end if
593: k = k + 1
594: end do
596: boundary_v(1 + i) = u1*u1 - u2*u2
597: if ((j == 0) .or. (j == 1)) then
598: xt = xt + hx
599: else
600: yt = yt + hy
601: end if
603: end do
605: if (j == 0) then
606: PetscCall(VecRestoreArray(Bottom, boundary_v, ierr))
607: elseif (j == 1) then
608: PetscCall(VecRestoreArray(Top, boundary_v, ierr))
609: elseif (j == 2) then
610: PetscCall(VecRestoreArray(Left, boundary_v, ierr))
611: elseif (j == 3) then
612: PetscCall(VecRestoreArray(Right, boundary_v, ierr))
613: end if
615: end do
617: ! Scale the boundary if desired
618: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bottom', scl, flg, ierr))
619: if (flg .eqv. PETSC_TRUE) then
620: PetscCall(VecScale(Bottom, scl, ierr))
621: end if
623: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-top', scl, flg, ierr))
624: if (flg .eqv. PETSC_TRUE) then
625: PetscCall(VecScale(Top, scl, ierr))
626: end if
628: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-right', scl, flg, ierr))
629: if (flg .eqv. PETSC_TRUE) then
630: PetscCall(VecScale(Right, scl, ierr))
631: end if
633: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-left', scl, flg, ierr))
634: if (flg .eqv. PETSC_TRUE) then
635: PetscCall(VecScale(Left, scl, ierr))
636: end if
638: end
640: ! ----------------------------------------------------------------------------
641: !
642: !/*
643: ! MSA_Plate - Calculates an obstacle for surface to stretch over
644: !
645: ! Output Parameter:
646: !. xl - lower bound vector
647: !. xu - upper bound vector
648: !
649: !*/
651: subroutine MSA_Plate(ta, xl, xu, dummy, ierr)
653: Tao ta
654: Vec xl, xu
655: PetscErrorCode, intent(out) :: ierr
656: PetscInt i, j, row
657: PetscInt xs, xm, ys, ym
658: PetscReal lb, ub
659: PetscInt dummy
660: PetscReal, pointer :: xl_v(:)
662: lb = PETSC_NINFINITY
663: ub = PETSC_INFINITY
665: if (bmy < 0) bmy = 0
666: if (bmy > my) bmy = my
667: if (bmx < 0) bmx = 0
668: if (bmx > mx) bmx = mx
670: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
672: PetscCall(VecSet(xl, lb, ierr))
673: PetscCall(VecSet(xu, ub, ierr))
675: PetscCall(VecGetArray(xl, xl_v, ierr))
677: do i = xs, xs + xm - 1
679: do j = ys, ys + ym - 1
681: row = (j - ys)*xm + (i - xs)
683: if (i >= ((mx - bmx)/2) .and. i < (mx - (mx - bmx)/2) .and. &
684: j >= ((my - bmy)/2) .and. j < (my - (my - bmy)/2)) then
685: xl_v(1 + row) = bheight
687: end if
689: end do
690: end do
692: PetscCall(VecRestoreArray(xl, xl_v, ierr))
694: end
696: ! ----------------------------------------------------------------------------
697: !
698: !/*
699: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
700: !
701: ! Output Parameter:
702: !. X - vector for initial guess
703: !
704: !*/
706: subroutine MSA_InitialPoint(X, ierr)
708: Vec X
709: PetscErrorCode, intent(out) :: ierr
710: PetscInt start, i, j
711: PetscInt row
712: PetscInt xs, xm, gxs, gxm
713: PetscInt ys, ym, gys, gym
714: PetscReal zero, np5
716: PetscReal, pointer :: left_v(:), right_v(:)
717: PetscReal, pointer :: bottom_v(:), top_v(:)
718: PetscReal, pointer :: x_v(:)
719: PetscBool flg
720: PetscRandom rctx
722: zero = 0.0
723: np5 = -0.5
725: PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-start', start, flg, ierr))
727: if ((flg .eqv. PETSC_TRUE) .and. (start == 0)) then ! the zero vector is reasonable
728: PetscCall(VecSet(X, zero, ierr))
730: elseif ((flg .eqv. PETSC_TRUE) .and. (start > 0)) then ! random start -0.5 < xi < 0.5
731: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, rctx, ierr))
732: do i = 0, start - 1
733: PetscCall(VecSetRandom(X, rctx, ierr))
734: end do
736: PetscCall(PetscRandomDestroy(rctx, ierr))
737: PetscCall(VecShift(X, np5, ierr))
739: else ! average of boundary conditions
741: ! Get Local mesh boundaries
742: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
743: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
745: ! Get pointers to vector data
746: PetscCall(VecGetArray(Top, top_v, ierr))
747: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
748: PetscCall(VecGetArray(Left, left_v, ierr))
749: PetscCall(VecGetArray(Right, right_v, ierr))
751: PetscCall(VecGetArray(localX, x_v, ierr))
753: ! Perform local computations
754: do j = ys, ys + ym - 1
755: do i = xs, xs + xm - 1
756: row = (j - gys)*gxm + (i - gxs)
757: x_v(1 + row) = ((j + 1)*bottom_v(1 + i - xs + 1)/my + (my - j + 1)*top_v(1 + i - xs + 1)/(my + 2) + &
758: (i + 1)*left_v(1 + j - ys + 1)/mx + (mx - i + 1)*right_v(1 + j - ys + 1)/(mx + 2))*0.5
759: end do
760: end do
762: ! Restore vectors
763: PetscCall(VecRestoreArray(localX, x_v, ierr))
765: PetscCall(VecRestoreArray(Left, left_v, ierr))
766: PetscCall(VecRestoreArray(Top, top_v, ierr))
767: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
768: PetscCall(VecRestoreArray(Right, right_v, ierr))
770: PetscCall(DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, X, ierr))
771: PetscCall(DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, X, ierr))
773: end if
775: end
777: end module
779: program plate2f
780: use plate2fmodule
781: implicit none
783: PetscErrorCode ierr ! used to check for functions returning nonzeros
784: Vec x ! solution vector
785: PetscInt m ! number of local elements in vector
786: Tao ta ! Tao solver context
787: Mat H ! Hessian matrix
788: ISLocalToGlobalMapping isltog ! local to global mapping object
789: PetscBool flg
791: ! Initialize Tao
792: PetscCallA(PetscInitialize(ierr))
794: ! Specify default dimensions of the problem and
795: ! check for any command line arguments that override defaults
796: mx = 10
797: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr))
798: my = 10
799: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr))
800: bmx = mx/2
801: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmx', bmx, flg, ierr))
802: bmy = my/2
803: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmy', bmy, flg, ierr))
804: bheight = 0.1
805: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bheight', bheight, flg, ierr))
807: ! Calculate any derived values from parameters
808: N = mx*my
810: ! Let PETSc determine the dimensions of the local vectors
811: Nx = PETSC_DECIDE
812: NY = PETSC_DECIDE
814: ! A two dimensional distributed array will help define this problem, which
815: ! derives from an elliptic PDE on a two-dimensional domain. From the
816: ! distributed array, create the vectors
818: PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, Nx, Ny, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
819: PetscCallA(DMSetFromOptions(dm, ierr))
820: PetscCallA(DMSetUp(dm, ierr))
822: ! Extract global and local vectors from DM; The local vectors are
823: ! used solely as work space for the evaluation of the function,
824: ! gradient, and Hessian. Duplicate for remaining vectors that are
825: ! the same types.
827: PetscCallA(DMCreateGlobalVector(dm, x, ierr))
828: PetscCallA(DMCreateLocalVector(dm, localX, ierr))
829: PetscCallA(VecDuplicate(localX, localV, ierr))
831: ! Create a matrix data structure to store the Hessian.
832: ! Here we (optionally) also associate the local numbering scheme
833: ! with the matrix so that later we can use local indices for matrix
834: ! assembly
836: PetscCallA(VecGetLocalSize(x, m, ierr))
837: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, m, m, N, N, 7_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, H, ierr))
839: PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))
840: PetscCallA(DMGetLocalToGlobalMapping(dm, isltog, ierr))
841: PetscCallA(MatSetLocalToGlobalMapping(H, isltog, isltog, ierr))
843: ! The Tao code begins here
844: ! Create TAO solver and set desired solution method.
845: ! This problems uses bounded variables, so the
846: ! method must either be 'tao_tron' or 'tao_blmvm'
848: PetscCallA(TaoCreate(PETSC_COMM_WORLD, ta, ierr))
849: PetscCallA(TaoSetType(ta, TAOBLMVM, ierr))
851: ! Set minimization function and gradient, hessian evaluation functions
853: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
855: PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))
857: ! Set Variable bounds
858: PetscCallA(MSA_BoundaryConditions(ierr))
859: PetscCallA(TaoSetVariableBoundsRoutine(ta, MSA_Plate, 0, ierr))
861: ! Set the initial solution guess
862: PetscCallA(MSA_InitialPoint(x, ierr))
863: PetscCallA(TaoSetSolution(ta, x, ierr))
865: ! Check for any tao command line options
866: PetscCallA(TaoSetFromOptions(ta, ierr))
868: ! Solve the application
869: PetscCallA(TaoSolve(ta, ierr))
871: ! Free TAO data structures
872: PetscCallA(TaoDestroy(ta, ierr))
874: ! Free PETSc data structures
875: PetscCallA(VecDestroy(x, ierr))
876: PetscCallA(VecDestroy(Top, ierr))
877: PetscCallA(VecDestroy(Bottom, ierr))
878: PetscCallA(VecDestroy(Left, ierr))
879: PetscCallA(VecDestroy(Right, ierr))
880: PetscCallA(MatDestroy(H, ierr))
881: PetscCallA(VecDestroy(localX, ierr))
882: PetscCallA(VecDestroy(localV, ierr))
883: PetscCallA(DMDestroy(dm, ierr))
885: ! Finalize TAO
886: PetscCallA(PetscFinalize(ierr))
888: end program plate2f
889: !
890: !/*TEST
891: !
892: ! build:
893: ! requires: !complex
894: !
895: ! test:
896: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
897: ! filter: sort -b
898: ! filter_output: sort -b
899: ! requires: !single
900: !
901: ! test:
902: ! suffix: 2
903: ! nsize: 2
904: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
905: ! filter: sort -b
906: ! filter_output: sort -b
907: ! requires: !single
908: !
909: !TEST*/