Actual source code: minsurf1.c

  1: #include <petsctao.h>

  3: static char help[] = "This example demonstrates use of the TAO package to\n\
  4: solve an unconstrained system of equations.  This example is based on a\n\
  5: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain and\n\
  6: boundary values along the edges of the domain, the objective is to find the\n\
  7: surface with the minimal area that satisfies the boundary conditions.\n\
  8: This application solves this problem using complimentarity -- We are actually\n\
  9: solving the system  (grad f)_i >= 0, if x_i == l_i \n\
 10:                     (grad f)_i = 0, if l_i < x_i < u_i \n\
 11:                     (grad f)_i <= 0, if x_i == u_i  \n\
 12: where f is the function to be minimized. \n\
 13: \n\
 14: The command line options are:\n\
 15:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 16:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 17:   -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \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:   PetscInt   mx, my;
 26:   PetscReal *bottom, *top, *left, *right;
 27: } AppCtx;

 29: /* -------- User-defined Routines --------- */

 31: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
 32: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
 33: PetscErrorCode        FormConstraints(Tao, Vec, Vec, void *);
 34: PetscErrorCode        FormJacobian(Tao, Vec, Mat, Mat, void *);

 36: int main(int argc, char **argv)
 37: {
 38:   Vec         x;                    /* solution vector */
 39:   Vec         c;                    /* Constraints function vector */
 40:   Vec         xl, xu;               /* Bounds on the variables */
 41:   PetscBool   flg;                  /* A return variable when checking for user options */
 42:   Tao         tao;                  /* TAO solver context */
 43:   Mat         J;                    /* Jacobian matrix */
 44:   PetscInt    N;                    /* Number of elements in vector */
 45:   PetscScalar lb = PETSC_NINFINITY; /* lower bound constant */
 46:   PetscScalar ub = PETSC_INFINITY;  /* upper bound constant */
 47:   AppCtx      user;                 /* user-defined work context */

 49:   /* Initialize PETSc, TAO */
 51:   PetscInitialize(&argc, &argv, (char *)0, help);

 53:   /* Specify default dimension of the problem */
 54:   user.mx = 4;
 55:   user.my = 4;

 57:   /* Check for any command line arguments that override defaults */
 58:   PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg);
 59:   PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg);

 61:   /* Calculate any derived values from parameters */
 62:   N = user.mx * user.my;

 64:   PetscPrintf(PETSC_COMM_SELF, "\n---- Minimum Surface Area Problem -----\n");
 65:   PetscPrintf(PETSC_COMM_SELF, "mx:%" PetscInt_FMT ", my:%" PetscInt_FMT "\n", user.mx, user.my);

 67:   /* Create appropriate vectors and matrices */
 68:   VecCreateSeq(MPI_COMM_SELF, N, &x);
 69:   VecDuplicate(x, &c);
 70:   MatCreateSeqAIJ(MPI_COMM_SELF, N, N, 7, NULL, &J);

 72:   /* The TAO code begins here */

 74:   /* Create TAO solver and set desired solution method */
 75:   TaoCreate(PETSC_COMM_SELF, &tao);
 76:   TaoSetType(tao, TAOSSILS);

 78:   /* Set data structure */
 79:   TaoSetSolution(tao, x);

 81:   /*  Set routines for constraints function and Jacobian evaluation */
 82:   TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user);
 83:   TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user);

 85:   /* Set the variable bounds */
 86:   MSA_BoundaryConditions(&user);

 88:   /* Set initial solution guess */
 89:   MSA_InitialPoint(&user, x);

 91:   /* Set Bounds on variables */
 92:   VecDuplicate(x, &xl);
 93:   VecDuplicate(x, &xu);
 94:   VecSet(xl, lb);
 95:   VecSet(xu, ub);
 96:   TaoSetVariableBounds(tao, xl, xu);

 98:   /* Check for any tao command line options */
 99:   TaoSetFromOptions(tao);

101:   /* Solve the application */
102:   TaoSolve(tao);

104:   /* Free Tao data structures */
105:   TaoDestroy(&tao);

107:   /* Free PETSc data structures */
108:   VecDestroy(&x);
109:   VecDestroy(&xl);
110:   VecDestroy(&xu);
111:   VecDestroy(&c);
112:   MatDestroy(&J);

114:   /* Free user-created data structures */
115:   PetscFree(user.bottom);
116:   PetscFree(user.top);
117:   PetscFree(user.left);
118:   PetscFree(user.right);

120:   PetscFinalize();
121:   return 0;
122: }

124: /* -------------------------------------------------------------------- */

126: /*  FormConstraints - Evaluates gradient of f.

128:     Input Parameters:
129: .   tao  - the TAO_APPLICATION context
130: .   X    - input vector
131: .   ptr  - optional user-defined context, as set by TaoSetConstraintsRoutine()

133:     Output Parameters:
134: .   G - vector containing the newly evaluated gradient
135: */
136: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec G, void *ptr)
137: {
138:   AppCtx      *user = (AppCtx *)ptr;
139:   PetscInt     i, j, row;
140:   PetscInt     mx = user->mx, my = user->my;
141:   PetscReal    hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
142:   PetscReal    f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
143:   PetscReal    df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
144:   PetscScalar  zero = 0.0;
145:   PetscScalar *g, *x;

147:   /* Initialize vector to zero */
148:   VecSet(G, zero);

150:   /* Get pointers to vector data */
151:   VecGetArray(X, &x);
152:   VecGetArray(G, &g);

154:   /* Compute function over the locally owned part of the mesh */
155:   for (j = 0; j < my; j++) {
156:     for (i = 0; i < mx; i++) {
157:       row = j * mx + i;

159:       xc  = x[row];
160:       xlt = xrb = xl = xr = xb = xt = xc;

162:       if (i == 0) { /* left side */
163:         xl  = user->left[j + 1];
164:         xlt = user->left[j + 2];
165:       } else {
166:         xl = x[row - 1];
167:       }

169:       if (j == 0) { /* bottom side */
170:         xb  = user->bottom[i + 1];
171:         xrb = user->bottom[i + 2];
172:       } else {
173:         xb = x[row - mx];
174:       }

176:       if (i + 1 == mx) { /* right side */
177:         xr  = user->right[j + 1];
178:         xrb = user->right[j];
179:       } else {
180:         xr = x[row + 1];
181:       }

183:       if (j + 1 == 0 + my) { /* top side */
184:         xt  = user->top[i + 1];
185:         xlt = user->top[i];
186:       } else {
187:         xt = x[row + mx];
188:       }

190:       if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
191:       if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];

193:       d1 = (xc - xl);
194:       d2 = (xc - xr);
195:       d3 = (xc - xt);
196:       d4 = (xc - xb);
197:       d5 = (xr - xrb);
198:       d6 = (xrb - xb);
199:       d7 = (xlt - xl);
200:       d8 = (xt - xlt);

202:       df1dxc = d1 * hydhx;
203:       df2dxc = (d1 * hydhx + d4 * hxdhy);
204:       df3dxc = d3 * hxdhy;
205:       df4dxc = (d2 * hydhx + d3 * hxdhy);
206:       df5dxc = d2 * hydhx;
207:       df6dxc = d4 * hxdhy;

209:       d1 /= hx;
210:       d2 /= hx;
211:       d3 /= hy;
212:       d4 /= hy;
213:       d5 /= hy;
214:       d6 /= hx;
215:       d7 /= hy;
216:       d8 /= hx;

218:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
219:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
220:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
221:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
222:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
223:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

225:       df1dxc /= f1;
226:       df2dxc /= f2;
227:       df3dxc /= f3;
228:       df4dxc /= f4;
229:       df5dxc /= f5;
230:       df6dxc /= f6;

232:       g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
233:     }
234:   }

236:   /* Restore vectors */
237:   VecRestoreArray(X, &x);
238:   VecRestoreArray(G, &g);
239:   PetscLogFlops(67 * mx * my);
240:   return 0;
241: }

243: /* ------------------------------------------------------------------- */
244: /*
245:    FormJacobian - Evaluates Jacobian matrix.

247:    Input Parameters:
248: .  tao  - the TAO_APPLICATION context
249: .  X    - input vector
250: .  ptr  - optional user-defined context, as set by TaoSetJacobian()

252:    Output Parameters:
253: .  tH    - Jacobian matrix

255: */
256: PetscErrorCode FormJacobian(Tao tao, Vec X, Mat H, Mat tHPre, void *ptr)
257: {
258:   AppCtx            *user = (AppCtx *)ptr;
259:   PetscInt           i, j, k, row;
260:   PetscInt           mx = user->mx, my = user->my;
261:   PetscInt           col[7];
262:   PetscReal          hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
263:   PetscReal          f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
264:   PetscReal          hl, hr, ht, hb, hc, htl, hbr;
265:   const PetscScalar *x;
266:   PetscScalar        v[7];
267:   PetscBool          assembled;

269:   /* Set various matrix options */
270:   MatSetOption(H, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
271:   MatAssembled(H, &assembled);
272:   if (assembled) MatZeroEntries(H);

274:   /* Get pointers to vector data */
275:   VecGetArrayRead(X, &x);

277:   /* Compute Jacobian over the locally owned part of the mesh */
278:   for (i = 0; i < mx; i++) {
279:     for (j = 0; j < my; j++) {
280:       row = j * mx + i;

282:       xc  = x[row];
283:       xlt = xrb = xl = xr = xb = xt = xc;

285:       /* Left side */
286:       if (i == 0) {
287:         xl  = user->left[j + 1];
288:         xlt = user->left[j + 2];
289:       } else {
290:         xl = x[row - 1];
291:       }

293:       if (j == 0) {
294:         xb  = user->bottom[i + 1];
295:         xrb = user->bottom[i + 2];
296:       } else {
297:         xb = x[row - mx];
298:       }

300:       if (i + 1 == mx) {
301:         xr  = user->right[j + 1];
302:         xrb = user->right[j];
303:       } else {
304:         xr = x[row + 1];
305:       }

307:       if (j + 1 == my) {
308:         xt  = user->top[i + 1];
309:         xlt = user->top[i];
310:       } else {
311:         xt = x[row + mx];
312:       }

314:       if (i > 0 && j + 1 < my) xlt = x[row - 1 + mx];
315:       if (j > 0 && i + 1 < mx) xrb = x[row + 1 - mx];

317:       d1 = (xc - xl) / hx;
318:       d2 = (xc - xr) / hx;
319:       d3 = (xc - xt) / hy;
320:       d4 = (xc - xb) / hy;
321:       d5 = (xrb - xr) / hy;
322:       d6 = (xrb - xb) / hx;
323:       d7 = (xlt - xl) / hy;
324:       d8 = (xlt - xt) / hx;

326:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
327:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
328:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
329:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
330:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
331:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

333:       hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
334:       hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
335:       ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
336:       hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);

338:       hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
339:       htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);

341:       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);

343:       hl /= 2.0;
344:       hr /= 2.0;
345:       ht /= 2.0;
346:       hb /= 2.0;
347:       hbr /= 2.0;
348:       htl /= 2.0;
349:       hc /= 2.0;

351:       k = 0;
352:       if (j > 0) {
353:         v[k]   = hb;
354:         col[k] = row - mx;
355:         k++;
356:       }

358:       if (j > 0 && i < mx - 1) {
359:         v[k]   = hbr;
360:         col[k] = row - mx + 1;
361:         k++;
362:       }

364:       if (i > 0) {
365:         v[k]   = hl;
366:         col[k] = row - 1;
367:         k++;
368:       }

370:       v[k]   = hc;
371:       col[k] = row;
372:       k++;

374:       if (i < mx - 1) {
375:         v[k]   = hr;
376:         col[k] = row + 1;
377:         k++;
378:       }

380:       if (i > 0 && j < my - 1) {
381:         v[k]   = htl;
382:         col[k] = row + mx - 1;
383:         k++;
384:       }

386:       if (j < my - 1) {
387:         v[k]   = ht;
388:         col[k] = row + mx;
389:         k++;
390:       }

392:       /*
393:          Set matrix values using local numbering, which was defined
394:          earlier, in the main routine.
395:       */
396:       MatSetValues(H, 1, &row, k, col, v, INSERT_VALUES);
397:     }
398:   }

400:   /* Restore vectors */
401:   VecRestoreArrayRead(X, &x);

403:   /* Assemble the matrix */
404:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
405:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
406:   PetscLogFlops(199 * mx * my);
407:   return 0;
408: }

410: /* ------------------------------------------------------------------- */
411: /*
412:    MSA_BoundaryConditions -  Calculates the boundary conditions for
413:    the region.

415:    Input Parameter:
416: .  user - user-defined application context

418:    Output Parameter:
419: .  user - user-defined application context
420: */
421: static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
422: {
423:   PetscInt   i, j, k, limit = 0, maxits = 5;
424:   PetscInt   mx = user->mx, my = user->my;
425:   PetscInt   bsize = 0, lsize = 0, tsize = 0, rsize = 0;
426:   PetscReal  one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
427:   PetscReal  fnorm, det, hx, hy, xt = 0, yt = 0;
428:   PetscReal  u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
429:   PetscReal  b = -0.5, t = 0.5, l = -0.5, r = 0.5;
430:   PetscReal *boundary;

432:   bsize = mx + 2;
433:   lsize = my + 2;
434:   rsize = my + 2;
435:   tsize = mx + 2;

437:   PetscMalloc1(bsize, &user->bottom);
438:   PetscMalloc1(tsize, &user->top);
439:   PetscMalloc1(lsize, &user->left);
440:   PetscMalloc1(rsize, &user->right);

442:   hx = (r - l) / (mx + 1);
443:   hy = (t - b) / (my + 1);

445:   for (j = 0; j < 4; j++) {
446:     if (j == 0) {
447:       yt       = b;
448:       xt       = l;
449:       limit    = bsize;
450:       boundary = user->bottom;
451:     } else if (j == 1) {
452:       yt       = t;
453:       xt       = l;
454:       limit    = tsize;
455:       boundary = user->top;
456:     } else if (j == 2) {
457:       yt       = b;
458:       xt       = l;
459:       limit    = lsize;
460:       boundary = user->left;
461:     } else { /* if  (j==3) */
462:       yt       = b;
463:       xt       = r;
464:       limit    = rsize;
465:       boundary = user->right;
466:     }

468:     for (i = 0; i < limit; i++) {
469:       u1 = xt;
470:       u2 = -yt;
471:       for (k = 0; k < maxits; k++) {
472:         nf1   = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
473:         nf2   = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
474:         fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
475:         if (fnorm <= tol) break;
476:         njac11 = one + u2 * u2 - u1 * u1;
477:         njac12 = two * u1 * u2;
478:         njac21 = -two * u1 * u2;
479:         njac22 = -one - u1 * u1 + u2 * u2;
480:         det    = njac11 * njac22 - njac21 * njac12;
481:         u1     = u1 - (njac22 * nf1 - njac12 * nf2) / det;
482:         u2     = u2 - (njac11 * nf2 - njac21 * nf1) / det;
483:       }

485:       boundary[i] = u1 * u1 - u2 * u2;
486:       if (j == 0 || j == 1) {
487:         xt = xt + hx;
488:       } else { /* if (j==2 || j==3) */
489:         yt = yt + hy;
490:       }
491:     }
492:   }
493:   return 0;
494: }

496: /* ------------------------------------------------------------------- */
497: /*
498:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

500:    Input Parameters:
501: .  user - user-defined application context
502: .  X - vector for initial guess

504:    Output Parameters:
505: .  X - newly computed initial guess
506: */
507: static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
508: {
509:   PetscInt    start = -1, i, j;
510:   PetscScalar zero  = 0.0;
511:   PetscBool   flg;

513:   PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg);

515:   if (flg && start == 0) { /* The zero vector is reasonable */
516:     VecSet(X, zero);
517:   } else { /* Take an average of the boundary conditions */
518:     PetscInt     row;
519:     PetscInt     mx = user->mx, my = user->my;
520:     PetscScalar *x;

522:     /* Get pointers to vector data */
523:     VecGetArray(X, &x);

525:     /* Perform local computations */
526:     for (j = 0; j < my; j++) {
527:       for (i = 0; i < mx; i++) {
528:         row    = (j)*mx + (i);
529:         x[row] = (((j + 1) * user->bottom[i + 1] + (my - j + 1) * user->top[i + 1]) / (my + 2) + ((i + 1) * user->left[j + 1] + (mx - i + 1) * user->right[j + 1]) / (mx + 2)) / 2.0;
530:       }
531:     }

533:     /* Restore vectors */
534:     VecRestoreArray(X, &x);
535:   }
536:   return 0;
537: }

539: /*TEST

541:    build:
542:       requires: !complex

544:    test:
545:       args: -tao_monitor -tao_view -tao_type ssils -tao_gttol 1.e-5
546:       requires: !single

548:    test:
549:       suffix: 2
550:       args: -tao_monitor -tao_view -tao_type ssfls -tao_gttol 1.e-5

552: TEST*/