Actual source code: plate2.c
1: #include <petscdmda.h>
2: #include <petsctao.h>
4: static char help[] = "This example demonstrates use of the TAO package to \n\
5: solve an unconstrained minimization problem. This example is based on a \n\
6: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain, \n\
7: boundary values along the edges of the domain, and a plate represented by \n\
8: lower boundary conditions, the objective is to find the\n\
9: surface with the minimal area that satisfies the boundary conditions.\n\
10: The command line options are:\n\
11: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
12: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
13: -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
14: -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
15: -bheight <ht>, where <ht> = height of the plate\n\
16: -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
17: for an average of the boundary conditions\n\n";
19: /*
20: User-defined application context - contains data needed by the
21: application-provided call-back routines, FormFunctionGradient(),
22: FormHessian().
23: */
24: typedef struct {
25: /* problem parameters */
26: PetscReal bheight; /* Height of plate under the surface */
27: PetscInt mx, my; /* discretization in x, y directions */
28: PetscInt bmx, bmy; /* Size of plate under the surface */
29: Vec Bottom, Top, Left, Right; /* boundary values */
31: /* Working space */
32: Vec localX, localV; /* ghosted local vector */
33: DM dm; /* distributed array data structure */
34: Mat H;
35: } AppCtx;
37: /* -------- User-defined Routines --------- */
39: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
40: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
41: static PetscErrorCode MSA_Plate(Vec, Vec, void *);
42: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
43: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
45: /* For testing matrix-free submatrices */
46: PetscErrorCode MatrixFreeHessian(Tao, Vec, Mat, Mat, void *);
47: PetscErrorCode MyMatMult(Mat, Vec, Vec);
49: int main(int argc, char **argv)
50: {
51: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
52: PetscInt m, N; /* number of local and global elements in vectors */
53: Vec x, xl, xu; /* solution vector and bounds*/
54: PetscBool flg; /* A return variable when checking for user options */
55: Tao tao; /* Tao solver context */
56: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
57: Mat H_shell; /* to test matrix-free submatrices */
58: AppCtx user; /* user-defined work context */
60: /* Initialize PETSc, TAO */
61: PetscFunctionBeginUser;
62: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
64: /* Specify default dimension of the problem */
65: user.mx = 10;
66: user.my = 10;
67: user.bheight = 0.1;
69: /* Check for any command line arguments that override defaults */
70: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
71: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
72: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bheight", &user.bheight, &flg));
74: user.bmx = user.mx / 2;
75: user.bmy = user.my / 2;
76: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bmx", &user.bmx, &flg));
77: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bmy", &user.bmy, &flg));
79: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n---- Minimum Surface Area With Plate Problem -----\n"));
80: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mx:%" PetscInt_FMT ", my:%" PetscInt_FMT ", bmx:%" PetscInt_FMT ", bmy:%" PetscInt_FMT ", height:%g\n", user.mx, user.my, user.bmx, user.bmy, (double)user.bheight));
82: /* Calculate any derived values from parameters */
83: N = user.mx * user.my;
85: /* Let Petsc determine the dimensions of the local vectors */
86: Nx = PETSC_DECIDE;
87: Ny = PETSC_DECIDE;
89: /*
90: A two dimensional distributed array will help define this problem,
91: which derives from an elliptic PDE on two dimensional domain. From
92: the distributed array, Create the vectors.
93: */
94: PetscCall(DMDACreate2d(MPI_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, user.mx, user.my, Nx, Ny, 1, 1, NULL, NULL, &user.dm));
95: PetscCall(DMSetFromOptions(user.dm));
96: PetscCall(DMSetUp(user.dm));
97: /*
98: Extract global and local vectors from DM; The local vectors are
99: used solely as work space for the evaluation of the function,
100: gradient, and Hessian. Duplicate for remaining vectors that are
101: the same types.
102: */
103: PetscCall(DMCreateGlobalVector(user.dm, &x)); /* Solution */
104: PetscCall(DMCreateLocalVector(user.dm, &user.localX));
105: PetscCall(VecDuplicate(user.localX, &user.localV));
107: PetscCall(VecDuplicate(x, &xl));
108: PetscCall(VecDuplicate(x, &xu));
110: /* The TAO code begins here */
112: /*
113: Create TAO solver and set desired solution method
114: The method must either be TAOTRON or TAOBLMVM
115: If TAOBLMVM is used, then hessian function is not called.
116: */
117: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
118: PetscCall(TaoSetType(tao, TAOBLMVM));
120: /* Set initial solution guess; */
121: PetscCall(MSA_BoundaryConditions(&user));
122: PetscCall(MSA_InitialPoint(&user, x));
123: PetscCall(TaoSetSolution(tao, x));
125: /* Set routines for function, gradient and hessian evaluation */
126: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
128: PetscCall(VecGetLocalSize(x, &m));
129: PetscCall(DMCreateMatrix(user.dm, &user.H));
130: PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
132: PetscCall(DMGetLocalToGlobalMapping(user.dm, &isltog));
133: PetscCall(MatSetLocalToGlobalMapping(user.H, isltog, isltog));
134: PetscCall(PetscOptionsHasName(NULL, NULL, "-matrixfree", &flg));
135: if (flg) {
136: PetscCall(MatCreateShell(PETSC_COMM_WORLD, m, m, N, N, (void *)&user, &H_shell));
137: PetscCall(MatShellSetOperation(H_shell, MATOP_MULT, (void (*)(void))MyMatMult));
138: PetscCall(MatSetOption(H_shell, MAT_SYMMETRIC, PETSC_TRUE));
139: PetscCall(TaoSetHessian(tao, H_shell, H_shell, MatrixFreeHessian, (void *)&user));
140: } else {
141: PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
142: }
144: /* Set Variable bounds */
145: PetscCall(MSA_Plate(xl, xu, (void *)&user));
146: PetscCall(TaoSetVariableBounds(tao, xl, xu));
148: /* Check for any tao command line options */
149: PetscCall(TaoSetFromOptions(tao));
151: /* SOLVE THE APPLICATION */
152: PetscCall(TaoSolve(tao));
154: PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD));
156: /* Free TAO data structures */
157: PetscCall(TaoDestroy(&tao));
159: /* Free PETSc data structures */
160: PetscCall(VecDestroy(&x));
161: PetscCall(VecDestroy(&xl));
162: PetscCall(VecDestroy(&xu));
163: PetscCall(MatDestroy(&user.H));
164: PetscCall(VecDestroy(&user.localX));
165: PetscCall(VecDestroy(&user.localV));
166: PetscCall(VecDestroy(&user.Bottom));
167: PetscCall(VecDestroy(&user.Top));
168: PetscCall(VecDestroy(&user.Left));
169: PetscCall(VecDestroy(&user.Right));
170: PetscCall(DMDestroy(&user.dm));
171: if (flg) PetscCall(MatDestroy(&H_shell));
172: PetscCall(PetscFinalize());
173: return 0;
174: }
176: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
178: Input Parameters:
179: . tao - the Tao context
180: . X - input vector
181: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
183: Output Parameters:
184: . fcn - the function value
185: . G - vector containing the newly evaluated gradient
187: Notes:
188: In this case, we discretize the domain and Create triangles. The
189: surface of each triangle is planar, whose surface area can be easily
190: computed. The total surface area is found by sweeping through the grid
191: and computing the surface area of the two triangles that have their
192: right angle at the grid point. The diagonal line segments on the
193: grid that define the triangles run from top left to lower right.
194: The numbering of points starts at the lower left and runs left to
195: right, then bottom to top.
196: */
197: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
198: {
199: AppCtx *user = (AppCtx *)userCtx;
200: PetscInt i, j, row;
201: PetscInt mx = user->mx, my = user->my;
202: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
203: PetscReal ft = 0;
204: PetscReal zero = 0.0;
205: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
206: PetscReal rhx = mx + 1, rhy = my + 1;
207: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
208: PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
209: PetscReal *g, *x, *left, *right, *bottom, *top;
210: Vec localX = user->localX, localG = user->localV;
212: PetscFunctionBeginUser;
213: /* Get local mesh boundaries */
214: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
215: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
217: /* Scatter ghost points to local vector */
218: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
219: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
221: /* Initialize vector to zero */
222: PetscCall(VecSet(localG, zero));
224: /* Get pointers to vector data */
225: PetscCall(VecGetArray(localX, &x));
226: PetscCall(VecGetArray(localG, &g));
227: PetscCall(VecGetArray(user->Top, &top));
228: PetscCall(VecGetArray(user->Bottom, &bottom));
229: PetscCall(VecGetArray(user->Left, &left));
230: PetscCall(VecGetArray(user->Right, &right));
232: /* Compute function over the locally owned part of the mesh */
233: for (j = ys; j < ys + ym; j++) {
234: for (i = xs; i < xs + xm; i++) {
235: row = (j - gys) * gxm + (i - gxs);
237: xc = x[row];
238: xlt = xrb = xl = xr = xb = xt = xc;
240: if (i == 0) { /* left side */
241: xl = left[j - ys + 1];
242: xlt = left[j - ys + 2];
243: } else {
244: xl = x[row - 1];
245: }
247: if (j == 0) { /* bottom side */
248: xb = bottom[i - xs + 1];
249: xrb = bottom[i - xs + 2];
250: } else {
251: xb = x[row - gxm];
252: }
254: if (i + 1 == gxs + gxm) { /* right side */
255: xr = right[j - ys + 1];
256: xrb = right[j - ys];
257: } else {
258: xr = x[row + 1];
259: }
261: if (j + 1 == gys + gym) { /* top side */
262: xt = top[i - xs + 1];
263: xlt = top[i - xs];
264: } else {
265: xt = x[row + gxm];
266: }
268: if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
269: if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];
271: d1 = (xc - xl);
272: d2 = (xc - xr);
273: d3 = (xc - xt);
274: d4 = (xc - xb);
275: d5 = (xr - xrb);
276: d6 = (xrb - xb);
277: d7 = (xlt - xl);
278: d8 = (xt - xlt);
280: df1dxc = d1 * hydhx;
281: df2dxc = (d1 * hydhx + d4 * hxdhy);
282: df3dxc = d3 * hxdhy;
283: df4dxc = (d2 * hydhx + d3 * hxdhy);
284: df5dxc = d2 * hydhx;
285: df6dxc = d4 * hxdhy;
287: d1 *= rhx;
288: d2 *= rhx;
289: d3 *= rhy;
290: d4 *= rhy;
291: d5 *= rhy;
292: d6 *= rhx;
293: d7 *= rhy;
294: d8 *= rhx;
296: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
297: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
298: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
299: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
300: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
301: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
303: ft = ft + (f2 + f4);
305: df1dxc /= f1;
306: df2dxc /= f2;
307: df3dxc /= f3;
308: df4dxc /= f4;
309: df5dxc /= f5;
310: df6dxc /= f6;
312: g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
313: }
314: }
316: /* Compute triangular areas along the border of the domain. */
317: if (xs == 0) { /* left side */
318: for (j = ys; j < ys + ym; j++) {
319: d3 = (left[j - ys + 1] - left[j - ys + 2]) * rhy;
320: d2 = (left[j - ys + 1] - x[(j - gys) * gxm]) * rhx;
321: ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
322: }
323: }
324: if (ys == 0) { /* bottom side */
325: for (i = xs; i < xs + xm; i++) {
326: d2 = (bottom[i + 1 - xs] - bottom[i - xs + 2]) * rhx;
327: d3 = (bottom[i - xs + 1] - x[i - gxs]) * rhy;
328: ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
329: }
330: }
332: if (xs + xm == mx) { /* right side */
333: for (j = ys; j < ys + ym; j++) {
334: d1 = (x[(j + 1 - gys) * gxm - 1] - right[j - ys + 1]) * rhx;
335: d4 = (right[j - ys] - right[j - ys + 1]) * rhy;
336: ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
337: }
338: }
339: if (ys + ym == my) { /* top side */
340: for (i = xs; i < xs + xm; i++) {
341: d1 = (x[(gym - 1) * gxm + i - gxs] - top[i - xs + 1]) * rhy;
342: d4 = (top[i - xs + 1] - top[i - xs]) * rhx;
343: ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
344: }
345: }
347: if (ys == 0 && xs == 0) {
348: d1 = (left[0] - left[1]) * rhy;
349: d2 = (bottom[0] - bottom[1]) * rhx;
350: ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
351: }
352: if (ys + ym == my && xs + xm == mx) {
353: d1 = (right[ym + 1] - right[ym]) * rhy;
354: d2 = (top[xm + 1] - top[xm]) * rhx;
355: ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
356: }
358: ft = ft * area;
359: PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
361: /* Restore vectors */
362: PetscCall(VecRestoreArray(localX, &x));
363: PetscCall(VecRestoreArray(localG, &g));
364: PetscCall(VecRestoreArray(user->Left, &left));
365: PetscCall(VecRestoreArray(user->Top, &top));
366: PetscCall(VecRestoreArray(user->Bottom, &bottom));
367: PetscCall(VecRestoreArray(user->Right, &right));
369: /* Scatter values to global vector */
370: PetscCall(DMLocalToGlobalBegin(user->dm, localG, INSERT_VALUES, G));
371: PetscCall(DMLocalToGlobalEnd(user->dm, localG, INSERT_VALUES, G));
373: PetscCall(PetscLogFlops(70.0 * xm * ym));
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: /* ------------------------------------------------------------------- */
378: /*
379: FormHessian - Evaluates Hessian matrix.
381: Input Parameters:
382: . tao - the Tao context
383: . x - input vector
384: . ptr - optional user-defined context, as set by TaoSetHessian()
386: Output Parameters:
387: . A - Hessian matrix
388: . B - optionally different preconditioning matrix
390: Notes:
391: Due to mesh point reordering with DMs, we must always work
392: with the local mesh points, and then transform them to the new
393: global numbering with the local-to-global mapping. We cannot work
394: directly with the global numbers for the original uniprocessor mesh!
396: Two methods are available for imposing this transformation
397: when setting matrix entries:
398: (A) MatSetValuesLocal(), using the local ordering (including
399: ghost points!)
400: - Do the following two steps once, before calling TaoSolve()
401: - Use DMGetISLocalToGlobalMapping() to extract the
402: local-to-global map from the DM
403: - Associate this map with the matrix by calling
404: MatSetLocalToGlobalMapping()
405: - Then set matrix entries using the local ordering
406: by calling MatSetValuesLocal()
407: (B) MatSetValues(), using the global ordering
408: - Use DMGetGlobalIndices() to extract the local-to-global map
409: - Then apply this map explicitly yourself
410: - Set matrix entries using the global ordering by calling
411: MatSetValues()
412: Option (A) seems cleaner/easier in many cases, and is the procedure
413: used in this example.
414: */
415: PetscErrorCode FormHessian(Tao tao, Vec X, Mat Hptr, Mat Hessian, void *ptr)
416: {
417: AppCtx *user = (AppCtx *)ptr;
418: PetscInt i, j, k, row;
419: PetscInt mx = user->mx, my = user->my;
420: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym, col[7];
421: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
422: PetscReal rhx = mx + 1, rhy = my + 1;
423: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
424: PetscReal hl, hr, ht, hb, hc, htl, hbr;
425: PetscReal *x, *left, *right, *bottom, *top;
426: PetscReal v[7];
427: Vec localX = user->localX;
428: PetscBool assembled;
430: PetscFunctionBeginUser;
431: /* Set various matrix options */
432: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
434: /* Initialize matrix entries to zero */
435: PetscCall(MatAssembled(Hessian, &assembled));
436: if (assembled) PetscCall(MatZeroEntries(Hessian));
438: /* Get local mesh boundaries */
439: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
440: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
442: /* Scatter ghost points to local vector */
443: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
444: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
446: /* Get pointers to vector data */
447: PetscCall(VecGetArray(localX, &x));
448: PetscCall(VecGetArray(user->Top, &top));
449: PetscCall(VecGetArray(user->Bottom, &bottom));
450: PetscCall(VecGetArray(user->Left, &left));
451: PetscCall(VecGetArray(user->Right, &right));
453: /* Compute Hessian over the locally owned part of the mesh */
455: for (i = xs; i < xs + xm; i++) {
456: for (j = ys; j < ys + ym; j++) {
457: row = (j - gys) * gxm + (i - gxs);
459: xc = x[row];
460: xlt = xrb = xl = xr = xb = xt = xc;
462: /* Left side */
463: if (i == gxs) {
464: xl = left[j - ys + 1];
465: xlt = left[j - ys + 2];
466: } else {
467: xl = x[row - 1];
468: }
470: if (j == gys) {
471: xb = bottom[i - xs + 1];
472: xrb = bottom[i - xs + 2];
473: } else {
474: xb = x[row - gxm];
475: }
477: if (i + 1 == gxs + gxm) {
478: xr = right[j - ys + 1];
479: xrb = right[j - ys];
480: } else {
481: xr = x[row + 1];
482: }
484: if (j + 1 == gys + gym) {
485: xt = top[i - xs + 1];
486: xlt = top[i - xs];
487: } else {
488: xt = x[row + gxm];
489: }
491: if (i > gxs && j + 1 < gys + gym) xlt = x[row - 1 + gxm];
492: if (j > gys && i + 1 < gxs + gxm) xrb = x[row + 1 - gxm];
494: d1 = (xc - xl) * rhx;
495: d2 = (xc - xr) * rhx;
496: d3 = (xc - xt) * rhy;
497: d4 = (xc - xb) * rhy;
498: d5 = (xrb - xr) * rhy;
499: d6 = (xrb - xb) * rhx;
500: d7 = (xlt - xl) * rhy;
501: d8 = (xlt - xt) * rhx;
503: f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
504: f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
505: f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
506: f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
507: f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
508: f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
510: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
511: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
512: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
513: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
515: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
516: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
518: 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) + 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) + hydhx * (1.0 + d3 * d3) - 2 * d2 * d3) / (f4 * f4 * f4);
520: hl *= 0.5;
521: hr *= 0.5;
522: ht *= 0.5;
523: hb *= 0.5;
524: hbr *= 0.5;
525: htl *= 0.5;
526: hc *= 0.5;
528: k = 0;
529: if (j > 0) {
530: v[k] = hb;
531: col[k] = row - gxm;
532: k++;
533: }
535: if (j > 0 && i < mx - 1) {
536: v[k] = hbr;
537: col[k] = row - gxm + 1;
538: k++;
539: }
541: if (i > 0) {
542: v[k] = hl;
543: col[k] = row - 1;
544: k++;
545: }
547: v[k] = hc;
548: col[k] = row;
549: k++;
551: if (i < mx - 1) {
552: v[k] = hr;
553: col[k] = row + 1;
554: k++;
555: }
557: if (i > 0 && j < my - 1) {
558: v[k] = htl;
559: col[k] = row + gxm - 1;
560: k++;
561: }
563: if (j < my - 1) {
564: v[k] = ht;
565: col[k] = row + gxm;
566: k++;
567: }
569: /*
570: Set matrix values using local numbering, which was defined
571: earlier, in the main routine.
572: */
573: PetscCall(MatSetValuesLocal(Hessian, 1, &row, k, col, v, INSERT_VALUES));
574: }
575: }
577: /* Restore vectors */
578: PetscCall(VecRestoreArray(localX, &x));
579: PetscCall(VecRestoreArray(user->Left, &left));
580: PetscCall(VecRestoreArray(user->Top, &top));
581: PetscCall(VecRestoreArray(user->Bottom, &bottom));
582: PetscCall(VecRestoreArray(user->Right, &right));
584: /* Assemble the matrix */
585: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
586: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
588: PetscCall(PetscLogFlops(199.0 * xm * ym));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /* ------------------------------------------------------------------- */
593: /*
594: MSA_BoundaryConditions - Calculates the boundary conditions for
595: the region.
597: Input Parameter:
598: . user - user-defined application context
600: Output Parameter:
601: . user - user-defined application context
602: */
603: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
604: {
605: PetscInt i, j, k, maxits = 5, limit = 0;
606: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
607: PetscInt mx = user->mx, my = user->my;
608: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
609: PetscReal one = 1.0, two = 2.0, three = 3.0, scl = 1.0, tol = 1e-10;
610: PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
611: PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
612: PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
613: PetscReal *boundary;
614: PetscBool flg;
615: Vec Bottom, Top, Right, Left;
617: PetscFunctionBeginUser;
618: /* Get local mesh boundaries */
619: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
620: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
622: bsize = xm + 2;
623: lsize = ym + 2;
624: rsize = ym + 2;
625: tsize = xm + 2;
627: PetscCall(VecCreateMPI(MPI_COMM_WORLD, bsize, PETSC_DECIDE, &Bottom));
628: PetscCall(VecCreateMPI(MPI_COMM_WORLD, tsize, PETSC_DECIDE, &Top));
629: PetscCall(VecCreateMPI(MPI_COMM_WORLD, lsize, PETSC_DECIDE, &Left));
630: PetscCall(VecCreateMPI(MPI_COMM_WORLD, rsize, PETSC_DECIDE, &Right));
632: user->Top = Top;
633: user->Left = Left;
634: user->Bottom = Bottom;
635: user->Right = Right;
637: hx = (r - l) / (mx + 1);
638: hy = (t - b) / (my + 1);
640: for (j = 0; j < 4; j++) {
641: if (j == 0) {
642: yt = b;
643: xt = l + hx * xs;
644: limit = bsize;
645: PetscCall(VecGetArray(Bottom, &boundary));
646: } else if (j == 1) {
647: yt = t;
648: xt = l + hx * xs;
649: limit = tsize;
650: PetscCall(VecGetArray(Top, &boundary));
651: } else if (j == 2) {
652: yt = b + hy * ys;
653: xt = l;
654: limit = lsize;
655: PetscCall(VecGetArray(Left, &boundary));
656: } else if (j == 3) {
657: yt = b + hy * ys;
658: xt = r;
659: limit = rsize;
660: PetscCall(VecGetArray(Right, &boundary));
661: }
663: for (i = 0; i < limit; i++) {
664: u1 = xt;
665: u2 = -yt;
666: for (k = 0; k < maxits; k++) {
667: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
668: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
669: fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
670: if (fnorm <= tol) break;
671: njac11 = one + u2 * u2 - u1 * u1;
672: njac12 = two * u1 * u2;
673: njac21 = -two * u1 * u2;
674: njac22 = -one - u1 * u1 + u2 * u2;
675: det = njac11 * njac22 - njac21 * njac12;
676: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
677: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
678: }
680: boundary[i] = u1 * u1 - u2 * u2;
681: if (j == 0 || j == 1) {
682: xt = xt + hx;
683: } else if (j == 2 || j == 3) {
684: yt = yt + hy;
685: }
686: }
687: if (j == 0) {
688: PetscCall(VecRestoreArray(Bottom, &boundary));
689: } else if (j == 1) {
690: PetscCall(VecRestoreArray(Top, &boundary));
691: } else if (j == 2) {
692: PetscCall(VecRestoreArray(Left, &boundary));
693: } else if (j == 3) {
694: PetscCall(VecRestoreArray(Right, &boundary));
695: }
696: }
698: /* Scale the boundary if desired */
700: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg));
701: if (flg) PetscCall(VecScale(Bottom, scl));
702: PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg));
703: if (flg) PetscCall(VecScale(Top, scl));
704: PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg));
705: if (flg) PetscCall(VecScale(Right, scl));
707: PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg));
708: if (flg) PetscCall(VecScale(Left, scl));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: /* ------------------------------------------------------------------- */
713: /*
714: MSA_Plate - Calculates an obstacle for surface to stretch over.
716: Input Parameter:
717: . user - user-defined application context
719: Output Parameter:
720: . user - user-defined application context
721: */
722: static PetscErrorCode MSA_Plate(Vec XL, Vec XU, void *ctx)
723: {
724: AppCtx *user = (AppCtx *)ctx;
725: PetscInt i, j, row;
726: PetscInt xs, ys, xm, ym;
727: PetscInt mx = user->mx, my = user->my, bmy, bmx;
728: PetscReal t1, t2, t3;
729: PetscReal *xl, lb = PETSC_NINFINITY, ub = PETSC_INFINITY;
730: PetscBool cylinder;
732: PetscFunctionBeginUser;
733: user->bmy = PetscMax(0, user->bmy);
734: user->bmy = PetscMin(my, user->bmy);
735: user->bmx = PetscMax(0, user->bmx);
736: user->bmx = PetscMin(mx, user->bmx);
737: bmy = user->bmy;
738: bmx = user->bmx;
740: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
742: PetscCall(VecSet(XL, lb));
743: PetscCall(VecSet(XU, ub));
745: PetscCall(VecGetArray(XL, &xl));
747: PetscCall(PetscOptionsHasName(NULL, NULL, "-cylinder", &cylinder));
748: /* Compute the optional lower box */
749: if (cylinder) {
750: for (i = xs; i < xs + xm; i++) {
751: for (j = ys; j < ys + ym; j++) {
752: row = (j - ys) * xm + (i - xs);
753: t1 = (2.0 * i - mx) * bmy;
754: t2 = (2.0 * j - my) * bmx;
755: t3 = bmx * bmx * bmy * bmy;
756: if (t1 * t1 + t2 * t2 <= t3) xl[row] = user->bheight;
757: }
758: }
759: } else {
760: /* Compute the optional lower box */
761: for (i = xs; i < xs + xm; i++) {
762: for (j = ys; j < ys + ym; j++) {
763: row = (j - ys) * xm + (i - xs);
764: if (i >= (mx - bmx) / 2 && i < mx - (mx - bmx) / 2 && j >= (my - bmy) / 2 && j < my - (my - bmy) / 2) xl[row] = user->bheight;
765: }
766: }
767: }
768: PetscCall(VecRestoreArray(XL, &xl));
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: /* ------------------------------------------------------------------- */
773: /*
774: MSA_InitialPoint - Calculates the initial guess in one of three ways.
776: Input Parameters:
777: . user - user-defined application context
778: . X - vector for initial guess
780: Output Parameters:
781: . X - newly computed initial guess
782: */
783: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
784: {
785: PetscInt start = -1, i, j;
786: PetscReal zero = 0.0;
787: PetscBool flg;
789: PetscFunctionBeginUser;
790: PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg));
791: if (flg && start == 0) { /* The zero vector is reasonable */
792: PetscCall(VecSet(X, zero));
793: } else if (flg && start > 0) { /* Try a random start between -0.5 and 0.5 */
794: PetscRandom rctx;
795: PetscReal np5 = -0.5;
797: PetscCall(PetscRandomCreate(MPI_COMM_WORLD, &rctx));
798: for (i = 0; i < start; i++) PetscCall(VecSetRandom(X, rctx));
799: PetscCall(PetscRandomDestroy(&rctx));
800: PetscCall(VecShift(X, np5));
802: } else { /* Take an average of the boundary conditions */
804: PetscInt row, xs, xm, gxs, gxm, ys, ym, gys, gym;
805: PetscInt mx = user->mx, my = user->my;
806: PetscReal *x, *left, *right, *bottom, *top;
807: Vec localX = user->localX;
809: /* Get local mesh boundaries */
810: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
811: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
813: /* Get pointers to vector data */
814: PetscCall(VecGetArray(user->Top, &top));
815: PetscCall(VecGetArray(user->Bottom, &bottom));
816: PetscCall(VecGetArray(user->Left, &left));
817: PetscCall(VecGetArray(user->Right, &right));
819: PetscCall(VecGetArray(localX, &x));
820: /* Perform local computations */
821: for (j = ys; j < ys + ym; j++) {
822: for (i = xs; i < xs + xm; i++) {
823: row = (j - gys) * gxm + (i - gxs);
824: x[row] = ((j + 1) * bottom[i - xs + 1] / my + (my - j + 1) * top[i - xs + 1] / (my + 2) + (i + 1) * left[j - ys + 1] / mx + (mx - i + 1) * right[j - ys + 1] / (mx + 2)) / 2.0;
825: }
826: }
828: /* Restore vectors */
829: PetscCall(VecRestoreArray(localX, &x));
831: PetscCall(VecRestoreArray(user->Left, &left));
832: PetscCall(VecRestoreArray(user->Top, &top));
833: PetscCall(VecRestoreArray(user->Bottom, &bottom));
834: PetscCall(VecRestoreArray(user->Right, &right));
836: /* Scatter values into global vector */
837: PetscCall(DMLocalToGlobalBegin(user->dm, localX, INSERT_VALUES, X));
838: PetscCall(DMLocalToGlobalEnd(user->dm, localX, INSERT_VALUES, X));
839: }
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: /* For testing matrix-free submatrices */
844: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
845: {
846: AppCtx *user = (AppCtx *)ptr;
848: PetscFunctionBegin;
849: PetscCall(FormHessian(tao, x, user->H, user->H, ptr));
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }
853: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
854: {
855: void *ptr;
856: AppCtx *user;
858: PetscFunctionBegin;
859: PetscCall(MatShellGetContext(H_shell, &ptr));
860: user = (AppCtx *)ptr;
861: PetscCall(MatMult(user->H, X, Y));
862: PetscFunctionReturn(PETSC_SUCCESS);
863: }
865: /*TEST
867: build:
868: requires: !complex
870: test:
871: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
872: requires: !single
874: test:
875: suffix: 2
876: nsize: 2
877: args: -tao_monitor_short -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
878: requires: !single
880: test:
881: suffix: 3
882: nsize: 3
883: args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
884: requires: !single
886: test:
887: suffix: 4
888: nsize: 3
889: args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
890: requires: !single
892: test:
893: suffix: 5
894: nsize: 3
895: args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -pc_type none -tao_type tron -tao_gatol 1.e-5
896: requires: !single
898: test:
899: suffix: 6
900: nsize: 3
901: args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
902: requires: !single
904: test:
905: suffix: 7
906: nsize: 3
907: args: -tao_monitor_short -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -pc_type none -tao_type gpcg -tao_gatol 1.e-5
908: requires: !single
910: test:
911: suffix: 8
912: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
913: requires: !single
915: test:
916: suffix: 9
917: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
918: requires: !single
920: test:
921: suffix: 10
922: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
923: requires: !single
925: test:
926: suffix: 11
927: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
928: requires: !single
930: test:
931: suffix: 12
932: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
933: requires: !single
935: test:
936: suffix: 13
937: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
938: requires: !single
940: test:
941: suffix: 14
942: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
943: requires: !single
945: test:
946: suffix: 15
947: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
948: requires: !single
950: test:
951: suffix: 16
952: args: -tao_monitor_short -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
953: requires: !single
955: test:
956: suffix: 17
957: args: -tao_monitor_short -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
958: requires: !single
960: test:
961: suffix: 18
962: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
963: requires: !single
965: test:
966: suffix: 19
967: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
968: requires: !single
970: test:
971: suffix: 20
972: args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
973: requires: !single
975: TEST*/