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 ta ! 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_ARRAY, PETSC_NULL_INTEGER_ARRAY, 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_ARRAY, i3, PETSC_NULL_INTEGER_ARRAY, 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, ta, ierr))
124: PetscCallA(TaoSetType(ta, TAOBLMVM, ierr))
126: ! Set minimization function and gradient, hessian evaluation functions
128: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
130: PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))
132: ! Set Variable bounds
133: PetscCallA(MSA_BoundaryConditions(ierr))
134: PetscCallA(TaoSetVariableBoundsRoutine(ta, MSA_Plate, 0, ierr))
136: ! Set the initial solution guess
137: PetscCallA(MSA_InitialPoint(x, ierr))
138: PetscCallA(TaoSetSolution(ta, x, ierr))
140: ! Check for any tao command line options
141: PetscCallA(TaoSetFromOptions(ta, ierr))
143: ! Solve the application
144: PetscCallA(TaoSolve(ta, ierr))
146: ! Free TAO data structures
147: PetscCallA(TaoDestroy(ta, 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(ta, X, fcn, G, dummy, ierr)
183: use plate2fmodule
184: implicit none
186: ! Input/output variables
188: Tao ta
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 VecGetArray in Fortran)
233: PetscCall(VecGetArray(localX, x_v, ierr))
234: PetscCall(VecGetArray(localV, g_v, ierr))
235: PetscCall(VecGetArray(Top, top_v, ierr))
236: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
237: PetscCall(VecGetArray(Left, left_v, ierr))
238: PetscCall(VecGetArray(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 == 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: end if
259: if (j == 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: end if
266: if (i + 1 == 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: end if
273: if (j + 1 == 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: end if
280: if ((i > gxs) .and. (j + 1 < gys + gym)) then
281: xlt = x_v(1 + row - 1 + gxm)
282: end if
284: if ((j > gys) .and. (i + 1 < gxs + gxm)) then
285: xrb = x_v(1 + row + 1 - gxm)
286: end if
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: end do
331: end do
333: ! Compute triangular areas along the border of the domain.
334: if (xs == 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: end do
340: end if
342: if (ys == 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: end do
348: end if
350: if (xs + xm == 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: end do
356: end if
358: if (ys + ym == 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: end do
364: end if
366: if ((ys == 0) .and. (xs == 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: end if
372: if ((ys + ym == my) .and. (xs + xm == 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: end if
378: ft = ft*area
379: PetscCallMPI(MPI_Allreduce(ft, fcn, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD, ierr))
381: ! Restore vectors
382: PetscCall(VecRestoreArray(localX, x_v, ierr))
383: PetscCall(VecRestoreArray(localV, g_v, ierr))
384: PetscCall(VecRestoreArray(Left, left_v, ierr))
385: PetscCall(VecRestoreArray(Top, top_v, ierr))
386: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
387: PetscCall(VecRestoreArray(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 matrix used to construct the preconditioner
410: !
411: ! Notes:
412: ! Due to mesh point reordering with DMs, we must always work
413: ! with the local mesh points, and then transform them to the new
414: ! global numbering with the local-to-global mapping. We cannot work
415: ! directly with the global numbers for the original uniprocessor mesh!
416: !
417: ! MatSetValuesLocal(), using the local ordering (including
418: ! ghost points!)
419: ! - Then set matrix entries using the local ordering
420: ! by calling MatSetValuesLocal()
422: subroutine FormHessian(ta, X, Hessian, Hpc, dummy, ierr)
423: use plate2fmodule
424: implicit none
426: Tao ta
427: Vec X
428: Mat Hessian, Hpc
429: PetscInt dummy
430: PetscErrorCode ierr
432: PetscInt i, j, k, row
433: PetscInt xs, xm, gxs, gxm
434: PetscInt ys, ym, gys, gym
435: PetscInt col(0:6)
436: PetscReal hx, hy, hydhx, hxdhy, rhx, rhy
437: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3
438: PetscReal d4, d5, d6, d7, d8
439: PetscReal xc, xl, xr, xt, xb, xlt, xrb
440: PetscReal hl, hr, ht, hb, hc, htl, hbr
442: PetscReal, pointer :: right_v(:), left_v(:)
443: PetscReal, pointer :: bottom_v(:), top_v(:)
444: PetscReal, pointer :: x_v(:)
445: PetscReal v(0:6)
446: PetscBool assembled
447: PetscInt i1
449: i1 = 1
451: ! Set various matrix options
452: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE, ierr))
454: ! Get local mesh boundaries
455: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
456: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
458: ! Scatter ghost points to local vectors
459: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr))
460: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr))
462: ! Get pointers to vector data (see note on Fortran arrays above)
463: PetscCall(VecGetArray(localX, x_v, ierr))
464: PetscCall(VecGetArray(Top, top_v, ierr))
465: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
466: PetscCall(VecGetArray(Left, left_v, ierr))
467: PetscCall(VecGetArray(Right, right_v, ierr))
469: ! Initialize matrix entries to zero
470: PetscCall(MatAssembled(Hessian, assembled, ierr))
471: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian, ierr))
473: rhx = real(mx) + 1.0
474: rhy = real(my) + 1.0
475: hx = 1.0/rhx
476: hy = 1.0/rhy
477: hydhx = hy/hx
478: hxdhy = hx/hy
479: ! compute Hessian over the locally owned part of the mesh
481: do i = xs, xs + xm - 1
482: do j = ys, ys + ym - 1
483: row = (j - gys)*gxm + (i - gxs)
485: xc = x_v(1 + row)
486: xt = xc
487: xb = xc
488: xr = xc
489: xl = xc
490: xrb = xc
491: xlt = xc
493: if (i == gxs) then ! Left side
494: xl = left_v(1 + j - ys + 1)
495: xlt = left_v(1 + j - ys + 2)
496: else
497: xl = x_v(1 + row - 1)
498: end if
500: if (j == gys) then ! bottom side
501: xb = bottom_v(1 + i - xs + 1)
502: xrb = bottom_v(1 + i - xs + 2)
503: else
504: xb = x_v(1 + row - gxm)
505: end if
507: if (i + 1 == gxs + gxm) then !right side
508: xr = right_v(1 + j - ys + 1)
509: xrb = right_v(1 + j - ys)
510: else
511: xr = x_v(1 + row + 1)
512: end if
514: if (j + 1 == gym + gys) then !top side
515: xt = top_v(1 + i - xs + 1)
516: xlt = top_v(1 + i - xs)
517: else
518: xt = x_v(1 + row + gxm)
519: end if
521: if ((i > gxs) .and. (j + 1 < gys + gym)) then
522: xlt = x_v(1 + row - 1 + gxm)
523: end if
525: if ((i + 1 < gxs + gxm) .and. (j > gys)) then
526: xrb = x_v(1 + row + 1 - gxm)
527: end if
529: d1 = (xc - xl)*rhx
530: d2 = (xc - xr)*rhx
531: d3 = (xc - xt)*rhy
532: d4 = (xc - xb)*rhy
533: d5 = (xrb - xr)*rhy
534: d6 = (xrb - xb)*rhx
535: d7 = (xlt - xl)*rhy
536: d8 = (xlt - xt)*rhx
538: f1 = sqrt(1.0 + d1*d1 + d7*d7)
539: f2 = sqrt(1.0 + d1*d1 + d4*d4)
540: f3 = sqrt(1.0 + d3*d3 + d8*d8)
541: f4 = sqrt(1.0 + d3*d3 + d2*d2)
542: f5 = sqrt(1.0 + d2*d2 + d5*d5)
543: f6 = sqrt(1.0 + d4*d4 + d6*d6)
545: hl = (-hydhx*(1.0 + d7*d7) + d1*d7)/(f1*f1*f1) + (-hydhx*(1.0 + d4*d4) + d1*d4)/(f2*f2*f2)
547: hr = (-hydhx*(1.0 + d5*d5) + d2*d5)/(f5*f5*f5) + (-hydhx*(1.0 + d3*d3) + d2*d3)/(f4*f4*f4)
549: ht = (-hxdhy*(1.0 + d8*d8) + d3*d8)/(f3*f3*f3) + (-hxdhy*(1.0 + d2*d2) + d2*d3)/(f4*f4*f4)
551: hb = (-hxdhy*(1.0 + d6*d6) + d4*d6)/(f6*f6*f6) + (-hxdhy*(1.0 + d1*d1) + d1*d4)/(f2*f2*f2)
553: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
554: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
556: 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) + &
557: & 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) + &
558: & hydhx*(1.0 + d3*d3) - 2*d2*d3)/(f4*f4*f4)
560: hl = hl*0.5
561: hr = hr*0.5
562: ht = ht*0.5
563: hb = hb*0.5
564: hbr = hbr*0.5
565: htl = htl*0.5
566: hc = hc*0.5
568: k = 0
570: if (j > 0) then
571: v(k) = hb
572: col(k) = row - gxm
573: k = k + 1
574: end if
576: if ((j > 0) .and. (i < mx - 1)) then
577: v(k) = hbr
578: col(k) = row - gxm + 1
579: k = k + 1
580: end if
582: if (i > 0) then
583: v(k) = hl
584: col(k) = row - 1
585: k = k + 1
586: end if
588: v(k) = hc
589: col(k) = row
590: k = k + 1
592: if (i < mx - 1) then
593: v(k) = hr
594: col(k) = row + 1
595: k = k + 1
596: end if
598: if ((i > 0) .and. (j < my - 1)) then
599: v(k) = htl
600: col(k) = row + gxm - 1
601: k = k + 1
602: end if
604: if (j < my - 1) then
605: v(k) = ht
606: col(k) = row + gxm
607: k = k + 1
608: end if
610: ! Set matrix values using local numbering, defined earlier in main routine
611: PetscCall(MatSetValuesLocal(Hessian, i1, [row], k, col, v, INSERT_VALUES, ierr))
613: end do
614: end do
616: ! restore vectors
617: PetscCall(VecRestoreArray(localX, x_v, ierr))
618: PetscCall(VecRestoreArray(Left, left_v, ierr))
619: PetscCall(VecRestoreArray(Right, right_v, ierr))
620: PetscCall(VecRestoreArray(Top, top_v, ierr))
621: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
623: ! Assemble the matrix
624: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY, ierr))
625: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY, ierr))
627: PetscCall(PetscLogFlops(199.0d0*xm*ym, ierr))
629: end
631: ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
633: ! ----------------------------------------------------------------------------
634: !
635: !/*
636: ! MSA_BoundaryConditions - calculates the boundary conditions for the region
637: !
638: !
639: !*/
641: subroutine MSA_BoundaryConditions(ierr)
642: use plate2fmodule
643: implicit none
645: PetscErrorCode ierr
646: PetscInt i, j, k, limit, maxits
647: PetscInt xs, xm, gxs, gxm
648: PetscInt ys, ym, gys, gym
649: PetscInt bsize, lsize
650: PetscInt tsize, rsize
651: PetscReal one, two, three, tol
652: PetscReal scl, fnorm, det, xt
653: PetscReal yt, hx, hy, u1, u2, nf1, nf2
654: PetscReal njac11, njac12, njac21, njac22
655: PetscReal b, t, l, r
656: PetscReal, pointer :: boundary_v(:)
658: logical exitloop
659: PetscBool flg
661: limit = 0
662: maxits = 5
663: tol = 1e-10
664: b = -0.5
665: t = 0.5
666: l = -0.5
667: r = 0.5
668: xt = 0
669: yt = 0
670: one = 1.0
671: two = 2.0
672: three = 3.0
674: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
675: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
677: bsize = xm + 2
678: lsize = ym + 2
679: rsize = ym + 2
680: tsize = xm + 2
682: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, bsize, PETSC_DECIDE, Bottom, ierr))
683: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, tsize, PETSC_DECIDE, Top, ierr))
684: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, lsize, PETSC_DECIDE, Left, ierr))
685: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, rsize, PETSC_DECIDE, Right, ierr))
687: hx = (r - l)/(mx + 1)
688: hy = (t - b)/(my + 1)
690: do j = 0, 3
692: if (j == 0) then
693: yt = b
694: xt = l + hx*xs
695: limit = bsize
696: PetscCall(VecGetArray(Bottom, boundary_v, ierr))
698: elseif (j == 1) then
699: yt = t
700: xt = l + hx*xs
701: limit = tsize
702: PetscCall(VecGetArray(Top, boundary_v, ierr))
704: elseif (j == 2) then
705: yt = b + hy*ys
706: xt = l
707: limit = lsize
708: PetscCall(VecGetArray(Left, boundary_v, ierr))
710: elseif (j == 3) then
711: yt = b + hy*ys
712: xt = r
713: limit = rsize
714: PetscCall(VecGetArray(Right, boundary_v, ierr))
715: end if
717: do i = 0, limit - 1
719: u1 = xt
720: u2 = -yt
721: k = 0
722: exitloop = .false.
723: do while (k < maxits .and. (.not. exitloop))
725: nf1 = u1 + u1*u2*u2 - u1*u1*u1/three - xt
726: nf2 = -u2 - u1*u1*u2 + u2*u2*u2/three - yt
727: fnorm = sqrt(nf1*nf1 + nf2*nf2)
728: if (fnorm > tol) then
729: njac11 = one + u2*u2 - u1*u1
730: njac12 = two*u1*u2
731: njac21 = -two*u1*u2
732: njac22 = -one - u1*u1 + u2*u2
733: det = njac11*njac22 - njac21*njac12
734: u1 = u1 - (njac22*nf1 - njac12*nf2)/det
735: u2 = u2 - (njac11*nf2 - njac21*nf1)/det
736: else
737: exitloop = .true.
738: end if
739: k = k + 1
740: end do
742: boundary_v(1 + i) = u1*u1 - u2*u2
743: if ((j == 0) .or. (j == 1)) then
744: xt = xt + hx
745: else
746: yt = yt + hy
747: end if
749: end do
751: if (j == 0) then
752: PetscCall(VecRestoreArray(Bottom, boundary_v, ierr))
753: elseif (j == 1) then
754: PetscCall(VecRestoreArray(Top, boundary_v, ierr))
755: elseif (j == 2) then
756: PetscCall(VecRestoreArray(Left, boundary_v, ierr))
757: elseif (j == 3) then
758: PetscCall(VecRestoreArray(Right, boundary_v, ierr))
759: end if
761: end do
763: ! Scale the boundary if desired
764: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bottom', scl, flg, ierr))
765: if (flg .eqv. PETSC_TRUE) then
766: PetscCall(VecScale(Bottom, scl, ierr))
767: end if
769: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-top', scl, flg, ierr))
770: if (flg .eqv. PETSC_TRUE) then
771: PetscCall(VecScale(Top, scl, ierr))
772: end if
774: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-right', scl, flg, ierr))
775: if (flg .eqv. PETSC_TRUE) then
776: PetscCall(VecScale(Right, scl, ierr))
777: end if
779: PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-left', scl, flg, ierr))
780: if (flg .eqv. PETSC_TRUE) then
781: PetscCall(VecScale(Left, scl, ierr))
782: end if
784: end
786: ! ----------------------------------------------------------------------------
787: !
788: !/*
789: ! MSA_Plate - Calculates an obstacle for surface to stretch over
790: !
791: ! Output Parameter:
792: !. xl - lower bound vector
793: !. xu - upper bound vector
794: !
795: !*/
797: subroutine MSA_Plate(ta, xl, xu, dummy, ierr)
798: use plate2fmodule
799: implicit none
801: Tao ta
802: Vec xl, xu
803: PetscErrorCode ierr
804: PetscInt i, j, row
805: PetscInt xs, xm, ys, ym
806: PetscReal lb, ub
807: PetscInt dummy
808: PetscReal, pointer :: xl_v(:)
810: lb = PETSC_NINFINITY
811: ub = PETSC_INFINITY
813: if (bmy < 0) bmy = 0
814: if (bmy > my) bmy = my
815: if (bmx < 0) bmx = 0
816: if (bmx > mx) bmx = mx
818: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
820: PetscCall(VecSet(xl, lb, ierr))
821: PetscCall(VecSet(xu, ub, ierr))
823: PetscCall(VecGetArray(xl, xl_v, ierr))
825: do i = xs, xs + xm - 1
827: do j = ys, ys + ym - 1
829: row = (j - ys)*xm + (i - xs)
831: if (i >= ((mx - bmx)/2) .and. i < (mx - (mx - bmx)/2) .and. &
832: & j >= ((my - bmy)/2) .and. j < (my - (my - bmy)/2)) then
833: xl_v(1 + row) = bheight
835: end if
837: end do
838: end do
840: PetscCall(VecRestoreArray(xl, xl_v, ierr))
842: end
844: ! ----------------------------------------------------------------------------
845: !
846: !/*
847: ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
848: !
849: ! Output Parameter:
850: !. X - vector for initial guess
851: !
852: !*/
854: subroutine MSA_InitialPoint(X, ierr)
855: use plate2fmodule
856: implicit none
858: Vec X
859: PetscErrorCode ierr
860: PetscInt start, i, j
861: PetscInt row
862: PetscInt xs, xm, gxs, gxm
863: PetscInt ys, ym, gys, gym
864: PetscReal zero, np5
866: PetscReal, pointer :: left_v(:), right_v(:)
867: PetscReal, pointer :: bottom_v(:), top_v(:)
868: PetscReal, pointer :: x_v(:)
869: PetscBool flg
870: PetscRandom rctx
872: zero = 0.0
873: np5 = -0.5
875: PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-start', start, flg, ierr))
877: if ((flg .eqv. PETSC_TRUE) .and. (start == 0)) then ! the zero vector is reasonable
878: PetscCall(VecSet(X, zero, ierr))
880: elseif ((flg .eqv. PETSC_TRUE) .and. (start > 0)) then ! random start -0.5 < xi < 0.5
881: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, rctx, ierr))
882: do i = 0, start - 1
883: PetscCall(VecSetRandom(X, rctx, ierr))
884: end do
886: PetscCall(PetscRandomDestroy(rctx, ierr))
887: PetscCall(VecShift(X, np5, ierr))
889: else ! average of boundary conditions
891: ! Get Local mesh boundaries
892: PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
893: PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr))
895: ! Get pointers to vector data
896: PetscCall(VecGetArray(Top, top_v, ierr))
897: PetscCall(VecGetArray(Bottom, bottom_v, ierr))
898: PetscCall(VecGetArray(Left, left_v, ierr))
899: PetscCall(VecGetArray(Right, right_v, ierr))
901: PetscCall(VecGetArray(localX, x_v, ierr))
903: ! Perform local computations
904: do j = ys, ys + ym - 1
905: do i = xs, xs + xm - 1
906: row = (j - gys)*gxm + (i - gxs)
907: x_v(1 + row) = ((j + 1)*bottom_v(1 + i - xs + 1)/my + (my - j + 1)*top_v(1 + i - xs + 1)/(my + 2) + &
908: & (i + 1)*left_v(1 + j - ys + 1)/mx + (mx - i + 1)*right_v(1 + j - ys + 1)/(mx + 2))*0.5
909: end do
910: end do
912: ! Restore vectors
913: PetscCall(VecRestoreArray(localX, x_v, ierr))
915: PetscCall(VecRestoreArray(Left, left_v, ierr))
916: PetscCall(VecRestoreArray(Top, top_v, ierr))
917: PetscCall(VecRestoreArray(Bottom, bottom_v, ierr))
918: PetscCall(VecRestoreArray(Right, right_v, ierr))
920: PetscCall(DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, X, ierr))
921: PetscCall(DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, X, ierr))
923: end if
925: end
927: !
928: !/*TEST
929: !
930: ! build:
931: ! requires: !complex
932: !
933: ! test:
934: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
935: ! filter: sort -b
936: ! filter_output: sort -b
937: ! requires: !single
938: !
939: ! test:
940: ! suffix: 2
941: ! nsize: 2
942: ! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4
943: ! filter: sort -b
944: ! filter_output: sort -b
945: ! requires: !single
946: !
947: !TEST*/