Actual source code: minsurf2.c
1: /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */
3: /*
4: Include "petsctao.h" so we can use TAO solvers.
5: petscdmda.h for distributed array
6: */
7: #include <petsctao.h>
8: #include <petscdmda.h>
10: static char help[] = "This example demonstrates use of the TAO package to \n\
11: solve an unconstrained minimization problem. This example is based on a \n\
12: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
13: boundary values along the edges of the domain, the objective is to find the\n\
14: surface with the minimal area that satisfies the boundary conditions.\n\
15: The command line options are:\n\
16: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
17: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
18: -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
19: for an average of the boundary conditions\n\n";
21: /*
22: User-defined application context - contains data needed by the
23: application-provided call-back routines, FormFunctionGradient()
24: and FormHessian().
25: */
26: typedef struct {
27: PetscInt mx, my; /* discretization in x, y directions */
28: PetscReal *bottom, *top, *left, *right; /* boundary values */
29: DM dm; /* distributed array data structure */
30: Mat H; /* Hessian */
31: } AppCtx;
33: /* -------- User-defined Routines --------- */
35: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
36: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
37: PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
38: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
39: PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
40: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
41: PetscErrorCode My_Monitor(Tao, void *);
43: int main(int argc, char **argv)
44: {
45: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
46: Vec x; /* solution, gradient vectors */
47: PetscBool flg, viewmat; /* flags */
48: PetscBool fddefault, fdcoloring; /* flags */
49: Tao tao; /* TAO solver context */
50: AppCtx user; /* user-defined work context */
51: ISColoring iscoloring;
52: MatFDColoring matfdcoloring;
54: /* Initialize TAO */
55: PetscFunctionBeginUser;
56: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
58: /* Specify dimension of the problem */
59: user.mx = 10;
60: user.my = 10;
62: /* Check for any command line arguments that override defaults */
63: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
64: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
66: PetscCall(PetscPrintf(MPI_COMM_WORLD, "\n---- Minimum Surface Area Problem -----\n"));
67: PetscCall(PetscPrintf(MPI_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my));
69: /* Let PETSc determine the vector distribution */
70: Nx = PETSC_DECIDE;
71: Ny = PETSC_DECIDE;
73: /* Create distributed array (DM) to manage parallel grid and vectors */
74: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, user.mx, user.my, Nx, Ny, 1, 1, NULL, NULL, &user.dm));
75: PetscCall(DMSetFromOptions(user.dm));
76: PetscCall(DMSetUp(user.dm));
78: /* Create TAO solver and set desired solution method.*/
79: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
80: PetscCall(TaoSetType(tao, TAOCG));
82: /*
83: Extract global vector from DA for the vector of variables -- PETSC routine
84: Compute the initial solution -- application specific, see below
85: Set this vector for use by TAO -- TAO routine
86: */
87: PetscCall(DMCreateGlobalVector(user.dm, &x));
88: PetscCall(MSA_BoundaryConditions(&user));
89: PetscCall(MSA_InitialPoint(&user, x));
90: PetscCall(TaoSetSolution(tao, x));
92: /*
93: Initialize the Application context for use in function evaluations -- application specific, see below.
94: Set routines for function and gradient evaluation
95: */
96: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
98: /*
99: Given the command line arguments, calculate the hessian with either the user-
100: provided function FormHessian, or the default finite-difference driven Hessian
101: functions
102: */
103: PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fddefault", &fddefault));
104: PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fdcoloring", &fdcoloring));
106: /*
107: Create a matrix data structure to store the Hessian and set
108: the Hessian evaluation routine.
109: Set the matrix structure to be used for Hessian evaluations
110: */
111: PetscCall(DMCreateMatrix(user.dm, &user.H));
112: PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
114: if (fdcoloring) {
115: PetscCall(DMCreateColoring(user.dm, IS_COLORING_GLOBAL, &iscoloring));
116: PetscCall(MatFDColoringCreate(user.H, iscoloring, &matfdcoloring));
117: PetscCall(ISColoringDestroy(&iscoloring));
118: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormGradient, (void *)&user));
119: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
120: PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessianColor, (void *)matfdcoloring));
121: } else if (fddefault) {
122: PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessian, (void *)NULL));
123: } else {
124: PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
125: }
127: /*
128: If my_monitor option is in command line, then use the user-provided
129: monitoring function
130: */
131: PetscCall(PetscOptionsHasName(NULL, NULL, "-my_monitor", &viewmat));
132: if (viewmat) PetscCall(TaoSetMonitor(tao, My_Monitor, NULL, NULL));
134: /* Check for any tao command line options */
135: PetscCall(TaoSetFromOptions(tao));
137: /* SOLVE THE APPLICATION */
138: PetscCall(TaoSolve(tao));
140: PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD));
142: /* Free TAO data structures */
143: PetscCall(TaoDestroy(&tao));
145: /* Free PETSc data structures */
146: PetscCall(VecDestroy(&x));
147: PetscCall(MatDestroy(&user.H));
148: if (fdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
149: PetscCall(PetscFree(user.bottom));
150: PetscCall(PetscFree(user.top));
151: PetscCall(PetscFree(user.left));
152: PetscCall(PetscFree(user.right));
153: PetscCall(DMDestroy(&user.dm));
154: PetscCall(PetscFinalize());
155: return 0;
156: }
158: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *userCtx)
159: {
160: PetscReal fcn;
162: PetscFunctionBegin;
163: PetscCall(FormFunctionGradient(tao, X, &fcn, G, userCtx));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /* -------------------------------------------------------------------- */
168: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
170: Input Parameters:
171: . tao - the Tao context
172: . XX - input vector
173: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
175: Output Parameters:
176: . fcn - the newly evaluated function
177: . GG - vector containing the newly evaluated gradient
178: */
179: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
180: {
181: AppCtx *user = (AppCtx *)userCtx;
182: PetscInt i, j;
183: PetscInt mx = user->mx, my = user->my;
184: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
185: PetscReal ft = 0.0;
186: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
187: PetscReal rhx = mx + 1, rhy = my + 1;
188: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
189: PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
190: PetscReal **g, **x;
191: Vec localX;
193: PetscFunctionBegin;
194: /* Get local mesh boundaries */
195: PetscCall(DMGetLocalVector(user->dm, &localX));
196: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
197: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
199: /* Scatter ghost points to local vector */
200: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
201: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
203: /* Get pointers to vector data */
204: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
205: PetscCall(DMDAVecGetArray(user->dm, G, (void **)&g));
207: /* Compute function and gradient over the locally owned part of the mesh */
208: for (j = ys; j < ys + ym; j++) {
209: for (i = xs; i < xs + xm; i++) {
210: xc = x[j][i];
211: xlt = xrb = xl = xr = xb = xt = xc;
213: if (i == 0) { /* left side */
214: xl = user->left[j - ys + 1];
215: xlt = user->left[j - ys + 2];
216: } else {
217: xl = x[j][i - 1];
218: }
220: if (j == 0) { /* bottom side */
221: xb = user->bottom[i - xs + 1];
222: xrb = user->bottom[i - xs + 2];
223: } else {
224: xb = x[j - 1][i];
225: }
227: if (i + 1 == gxs + gxm) { /* right side */
228: xr = user->right[j - ys + 1];
229: xrb = user->right[j - ys];
230: } else {
231: xr = x[j][i + 1];
232: }
234: if (j + 1 == gys + gym) { /* top side */
235: xt = user->top[i - xs + 1];
236: xlt = user->top[i - xs];
237: } else {
238: xt = x[j + 1][i];
239: }
241: if (i > gxs && j + 1 < gys + gym) xlt = x[j + 1][i - 1];
242: if (j > gys && i + 1 < gxs + gxm) xrb = x[j - 1][i + 1];
244: d1 = (xc - xl);
245: d2 = (xc - xr);
246: d3 = (xc - xt);
247: d4 = (xc - xb);
248: d5 = (xr - xrb);
249: d6 = (xrb - xb);
250: d7 = (xlt - xl);
251: d8 = (xt - xlt);
253: df1dxc = d1 * hydhx;
254: df2dxc = (d1 * hydhx + d4 * hxdhy);
255: df3dxc = d3 * hxdhy;
256: df4dxc = (d2 * hydhx + d3 * hxdhy);
257: df5dxc = d2 * hydhx;
258: df6dxc = d4 * hxdhy;
260: d1 *= rhx;
261: d2 *= rhx;
262: d3 *= rhy;
263: d4 *= rhy;
264: d5 *= rhy;
265: d6 *= rhx;
266: d7 *= rhy;
267: d8 *= rhx;
269: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
270: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
271: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
272: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
273: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
274: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
276: ft = ft + (f2 + f4);
278: df1dxc /= f1;
279: df2dxc /= f2;
280: df3dxc /= f3;
281: df4dxc /= f4;
282: df5dxc /= f5;
283: df6dxc /= f6;
285: g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
286: }
287: }
289: /* Compute triangular areas along the border of the domain. */
290: if (xs == 0) { /* left side */
291: for (j = ys; j < ys + ym; j++) {
292: d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
293: d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
294: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
295: }
296: }
297: if (ys == 0) { /* bottom side */
298: for (i = xs; i < xs + xm; i++) {
299: d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
300: d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
301: ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
302: }
303: }
305: if (xs + xm == mx) { /* right side */
306: for (j = ys; j < ys + ym; j++) {
307: d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
308: d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
309: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
310: }
311: }
312: if (ys + ym == my) { /* top side */
313: for (i = xs; i < xs + xm; i++) {
314: d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
315: d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
316: ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
317: }
318: }
320: if (ys == 0 && xs == 0) {
321: d1 = (user->left[0] - user->left[1]) / hy;
322: d2 = (user->bottom[0] - user->bottom[1]) * rhx;
323: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
324: }
325: if (ys + ym == my && xs + xm == mx) {
326: d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
327: d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
328: ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
329: }
331: ft = ft * area;
332: PetscCall(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
334: /* Restore vectors */
335: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
336: PetscCall(DMDAVecRestoreArray(user->dm, G, (void **)&g));
338: /* Scatter values to global vector */
339: PetscCall(DMRestoreLocalVector(user->dm, &localX));
340: PetscCall(PetscLogFlops(67.0 * xm * ym));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: /* ------------------------------------------------------------------- */
345: /*
346: FormHessian - Evaluates Hessian matrix.
348: Input Parameters:
349: . tao - the Tao context
350: . x - input vector
351: . ptr - optional user-defined context, as set by TaoSetHessian()
353: Output Parameters:
354: . H - Hessian matrix
355: . Hpre - optionally different preconditioning matrix
356: . flg - flag indicating matrix structure
358: */
359: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
360: {
361: AppCtx *user = (AppCtx *)ptr;
363: PetscFunctionBegin;
364: /* Evaluate the Hessian entries*/
365: PetscCall(QuadraticH(user, X, H));
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: /* ------------------------------------------------------------------- */
370: /*
371: QuadraticH - Evaluates Hessian matrix.
373: Input Parameters:
374: . user - user-defined context, as set by TaoSetHessian()
375: . X - input vector
377: Output Parameter:
378: . H - Hessian matrix
379: */
380: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
381: {
382: PetscInt i, j, k;
383: PetscInt mx = user->mx, my = user->my;
384: PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
385: PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
386: PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
387: PetscReal hl, hr, ht, hb, hc, htl, hbr;
388: PetscReal **x, v[7];
389: MatStencil col[7], row;
390: Vec localX;
391: PetscBool assembled;
393: PetscFunctionBegin;
394: /* Get local mesh boundaries */
395: PetscCall(DMGetLocalVector(user->dm, &localX));
397: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
398: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
400: /* Scatter ghost points to local vector */
401: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
402: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
404: /* Get pointers to vector data */
405: PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
407: /* Initialize matrix entries to zero */
408: PetscCall(MatAssembled(Hessian, &assembled));
409: if (assembled) PetscCall(MatZeroEntries(Hessian));
411: /* Set various matrix options */
412: PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
414: /* Compute Hessian over the locally owned part of the mesh */
416: for (j = ys; j < ys + ym; j++) {
417: for (i = xs; i < xs + xm; i++) {
418: xc = x[j][i];
419: xlt = xrb = xl = xr = xb = xt = xc;
421: /* Left side */
422: if (i == 0) {
423: xl = user->left[j - ys + 1];
424: xlt = user->left[j - ys + 2];
425: } else {
426: xl = x[j][i - 1];
427: }
429: if (j == 0) {
430: xb = user->bottom[i - xs + 1];
431: xrb = user->bottom[i - xs + 2];
432: } else {
433: xb = x[j - 1][i];
434: }
436: if (i + 1 == mx) {
437: xr = user->right[j - ys + 1];
438: xrb = user->right[j - ys];
439: } else {
440: xr = x[j][i + 1];
441: }
443: if (j + 1 == my) {
444: xt = user->top[i - xs + 1];
445: xlt = user->top[i - xs];
446: } else {
447: xt = x[j + 1][i];
448: }
450: if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
451: if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
453: d1 = (xc - xl) / hx;
454: d2 = (xc - xr) / hx;
455: d3 = (xc - xt) / hy;
456: d4 = (xc - xb) / hy;
457: d5 = (xrb - xr) / hy;
458: d6 = (xrb - xb) / hx;
459: d7 = (xlt - xl) / hy;
460: d8 = (xlt - xt) / hx;
462: f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
463: f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
464: f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
465: f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
466: f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
467: f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
469: hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
470: hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
471: ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
472: hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
474: hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
475: htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
477: 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);
479: hl /= 2.0;
480: hr /= 2.0;
481: ht /= 2.0;
482: hb /= 2.0;
483: hbr /= 2.0;
484: htl /= 2.0;
485: hc /= 2.0;
487: row.j = j;
488: row.i = i;
489: k = 0;
490: if (j > 0) {
491: v[k] = hb;
492: col[k].j = j - 1;
493: col[k].i = i;
494: k++;
495: }
497: if (j > 0 && i < mx - 1) {
498: v[k] = hbr;
499: col[k].j = j - 1;
500: col[k].i = i + 1;
501: k++;
502: }
504: if (i > 0) {
505: v[k] = hl;
506: col[k].j = j;
507: col[k].i = i - 1;
508: k++;
509: }
511: v[k] = hc;
512: col[k].j = j;
513: col[k].i = i;
514: k++;
516: if (i < mx - 1) {
517: v[k] = hr;
518: col[k].j = j;
519: col[k].i = i + 1;
520: k++;
521: }
523: if (i > 0 && j < my - 1) {
524: v[k] = htl;
525: col[k].j = j + 1;
526: col[k].i = i - 1;
527: k++;
528: }
530: if (j < my - 1) {
531: v[k] = ht;
532: col[k].j = j + 1;
533: col[k].i = i;
534: k++;
535: }
537: /*
538: Set matrix values using local numbering, which was defined
539: earlier, in the main routine.
540: */
541: PetscCall(MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES));
542: }
543: }
545: PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
546: PetscCall(DMRestoreLocalVector(user->dm, &localX));
548: PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
549: PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
551: PetscCall(PetscLogFlops(199.0 * xm * ym));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /* ------------------------------------------------------------------- */
556: /*
557: MSA_BoundaryConditions - Calculates the boundary conditions for
558: the region.
560: Input Parameter:
561: . user - user-defined application context
563: Output Parameter:
564: . user - user-defined application context
565: */
566: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
567: {
568: PetscInt i, j, k, limit = 0, maxits = 5;
569: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
570: PetscInt mx = user->mx, my = user->my;
571: PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
572: PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
573: PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
574: PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
575: PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
576: PetscReal *boundary;
577: PetscBool flg;
579: PetscFunctionBegin;
580: /* Get local mesh boundaries */
581: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
582: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
584: bsize = xm + 2;
585: lsize = ym + 2;
586: rsize = ym + 2;
587: tsize = xm + 2;
589: PetscCall(PetscMalloc1(bsize, &user->bottom));
590: PetscCall(PetscMalloc1(tsize, &user->top));
591: PetscCall(PetscMalloc1(lsize, &user->left));
592: PetscCall(PetscMalloc1(rsize, &user->right));
594: hx = (r - l) / (mx + 1);
595: hy = (t - b) / (my + 1);
597: for (j = 0; j < 4; j++) {
598: if (j == 0) {
599: yt = b;
600: xt = l + hx * xs;
601: limit = bsize;
602: boundary = user->bottom;
603: } else if (j == 1) {
604: yt = t;
605: xt = l + hx * xs;
606: limit = tsize;
607: boundary = user->top;
608: } else if (j == 2) {
609: yt = b + hy * ys;
610: xt = l;
611: limit = lsize;
612: boundary = user->left;
613: } else { /* if (j==3) */
614: yt = b + hy * ys;
615: xt = r;
616: limit = rsize;
617: boundary = user->right;
618: }
620: for (i = 0; i < limit; i++) {
621: u1 = xt;
622: u2 = -yt;
623: for (k = 0; k < maxits; k++) {
624: nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
625: nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
626: fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2);
627: if (fnorm <= tol) break;
628: njac11 = one + u2 * u2 - u1 * u1;
629: njac12 = two * u1 * u2;
630: njac21 = -two * u1 * u2;
631: njac22 = -one - u1 * u1 + u2 * u2;
632: det = njac11 * njac22 - njac21 * njac12;
633: u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
634: u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
635: }
637: boundary[i] = u1 * u1 - u2 * u2;
638: if (j == 0 || j == 1) {
639: xt = xt + hx;
640: } else { /* if (j==2 || j==3) */
641: yt = yt + hy;
642: }
643: }
644: }
646: /* Scale the boundary if desired */
647: if (1 == 1) {
648: PetscReal scl = 1.0;
650: PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg));
651: if (flg) {
652: for (i = 0; i < bsize; i++) user->bottom[i] *= scl;
653: }
655: PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg));
656: if (flg) {
657: for (i = 0; i < tsize; i++) user->top[i] *= scl;
658: }
660: PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg));
661: if (flg) {
662: for (i = 0; i < rsize; i++) user->right[i] *= scl;
663: }
665: PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg));
666: if (flg) {
667: for (i = 0; i < lsize; i++) user->left[i] *= scl;
668: }
669: }
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: /* ------------------------------------------------------------------- */
674: /*
675: MSA_InitialPoint - Calculates the initial guess in one of three ways.
677: Input Parameters:
678: . user - user-defined application context
679: . X - vector for initial guess
681: Output Parameters:
682: . X - newly computed initial guess
683: */
684: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
685: {
686: PetscInt start2 = -1, i, j;
687: PetscReal start1 = 0;
688: PetscBool flg1, flg2;
690: PetscFunctionBegin;
691: PetscCall(PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1));
692: PetscCall(PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2));
694: if (flg1) { /* The zero vector is reasonable */
696: PetscCall(VecSet(X, start1));
698: } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */
699: PetscRandom rctx;
700: PetscReal np5 = -0.5;
702: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
703: PetscCall(VecSetRandom(X, rctx));
704: PetscCall(PetscRandomDestroy(&rctx));
705: PetscCall(VecShift(X, np5));
707: } else { /* Take an average of the boundary conditions */
708: PetscInt xs, xm, ys, ym;
709: PetscInt mx = user->mx, my = user->my;
710: PetscReal **x;
712: /* Get local mesh boundaries */
713: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
715: /* Get pointers to vector data */
716: PetscCall(DMDAVecGetArray(user->dm, X, (void **)&x));
718: /* Perform local computations */
719: for (j = ys; j < ys + ym; j++) {
720: for (i = xs; i < xs + xm; i++) x[j][i] = (((j + 1) * user->bottom[i - xs + 1] + (my - j + 1) * user->top[i - xs + 1]) / (my + 2) + ((i + 1) * user->left[j - ys + 1] + (mx - i + 1) * user->right[j - ys + 1]) / (mx + 2)) / 2.0;
721: }
722: PetscCall(DMDAVecRestoreArray(user->dm, X, (void **)&x));
723: PetscCall(PetscLogFlops(9.0 * xm * ym));
724: }
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: /*-----------------------------------------------------------------------*/
729: PetscErrorCode My_Monitor(Tao tao, void *ctx)
730: {
731: Vec X;
733: PetscFunctionBegin;
734: PetscCall(TaoGetSolution(tao, &X));
735: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: /*TEST
741: build:
742: requires: !complex
744: test:
745: args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
746: requires: !single
748: test:
749: suffix: 2
750: nsize: 2
751: args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
752: filter: grep -v "nls ksp"
753: requires: !single
755: test:
756: suffix: 3
757: nsize: 3
758: args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
759: requires: !single
761: test:
762: suffix: 5
763: nsize: 2
764: args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
765: requires: !single
767: TEST*/