Actual source code: plate2.c

  1: #include <petscdmda.h>
  2: #include <petsctao.h>

  4: static  char help[] =
  5: "This example demonstrates use of the TAO package to \n\
  6: solve an unconstrained minimization problem.  This example is based on a \n\
  7: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain, \n\
  8: boundary values along the edges of the domain, and a plate represented by \n\
  9: lower boundary conditions, the objective is to find the\n\
 10: surface with the minimal area that satisfies the boundary conditions.\n\
 11: The command line options are:\n\
 12:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 13:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 14:   -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
 15:   -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
 16:   -bheight <ht>, where <ht> = height of the plate\n\
 17:   -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
 18:                for an average of the boundary conditions\n\n";

 20: /*T
 21:    Concepts: TAO^Solving a bound constrained minimization problem
 22:    Routines: TaoCreate();
 23:    Routines: TaoSetType(); TaoSetObjectiveAndGradient();
 24:    Routines: TaoSetHessian();
 25:    Routines: TaoSetSolution();
 26:    Routines: TaoSetVariableBounds();
 27:    Routines: TaoSetFromOptions();
 28:    Routines: TaoSolve(); TaoView();
 29:    Routines: TaoDestroy();
 30:    Processors: n
 31: T*/

 33: /*
 34:    User-defined application context - contains data needed by the
 35:    application-provided call-back routines, FormFunctionGradient(),
 36:    FormHessian().
 37: */
 38: typedef struct {
 39:   /* problem parameters */
 40:   PetscReal      bheight;                  /* Height of plate under the surface */
 41:   PetscInt       mx, my;                   /* discretization in x, y directions */
 42:   PetscInt       bmx,bmy;                  /* Size of plate under the surface */
 43:   Vec            Bottom, Top, Left, Right; /* boundary values */

 45:   /* Working space */
 46:   Vec         localX, localV;           /* ghosted local vector */
 47:   DM          dm;                       /* distributed array data structure */
 48:   Mat         H;
 49: } AppCtx;

 51: /* -------- User-defined Routines --------- */

 53: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
 54: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
 55: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
 56: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
 57: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);

 59: /* For testing matrix free submatrices */
 60: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
 61: PetscErrorCode MyMatMult(Mat,Vec,Vec);

 63: int main(int argc, char **argv)
 64: {
 65:   PetscInt               Nx, Ny;               /* number of processors in x- and y- directions */
 66:   PetscInt               m, N;                 /* number of local and global elements in vectors */
 67:   Vec                    x,xl,xu;               /* solution vector  and bounds*/
 68:   PetscBool              flg;                /* A return variable when checking for user options */
 69:   Tao                    tao;                  /* Tao solver context */
 70:   ISLocalToGlobalMapping isltog;   /* local-to-global mapping object */
 71:   Mat                    H_shell;                  /* to test matrix-free submatrices */
 72:   AppCtx                 user;                 /* user-defined work context */

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

 77:   /* Specify default dimension of the problem */
 78:   user.mx = 10; user.my = 10; user.bheight=0.1;

 80:   /* Check for any command line arguments that override defaults */
 81:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
 82:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
 83:   PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);

 85:   user.bmx = user.mx/2; user.bmy = user.my/2;
 86:   PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);
 87:   PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);

 89:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
 90:   PetscPrintf(PETSC_COMM_WORLD,"mx:%D, my:%D, bmx:%D, bmy:%D, height:%g\n",user.mx,user.my,user.bmx,user.bmy,(double)user.bheight);

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

 95:   /* Let Petsc determine the dimensions of the local vectors */
 96:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 98:   /*
 99:      A two dimensional distributed array will help define this problem,
100:      which derives from an elliptic PDE on two dimensional domain.  From
101:      the distributed array, Create the vectors.
102:   */
103:   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);
104:   DMSetFromOptions(user.dm);
105:   DMSetUp(user.dm);
106:   /*
107:      Extract global and local vectors from DM; The local vectors are
108:      used solely as work space for the evaluation of the function,
109:      gradient, and Hessian.  Duplicate for remaining vectors that are
110:      the same types.
111:   */
112:   DMCreateGlobalVector(user.dm,&x); /* Solution */
113:   DMCreateLocalVector(user.dm,&user.localX);
114:   VecDuplicate(user.localX,&user.localV);

116:   VecDuplicate(x,&xl);
117:   VecDuplicate(x,&xu);

119:   /* The TAO code begins here */

121:   /*
122:      Create TAO solver and set desired solution method
123:      The method must either be TAOTRON or TAOBLMVM
124:      If TAOBLMVM is used, then hessian function is not called.
125:   */
126:   TaoCreate(PETSC_COMM_WORLD,&tao);
127:   TaoSetType(tao,TAOBLMVM);

129:   /* Set initial solution guess; */
130:   MSA_BoundaryConditions(&user);
131:   MSA_InitialPoint(&user,x);
132:   TaoSetSolution(tao,x);

134:   /* Set routines for function, gradient and hessian evaluation */
135:   TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*) &user);

137:   VecGetLocalSize(x,&m);
138:   MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
139:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);

141:   DMGetLocalToGlobalMapping(user.dm,&isltog);
142:   MatSetLocalToGlobalMapping(user.H,isltog,isltog);
143:   PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);
144:   if (flg) {
145:       MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
146:       MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
147:       MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
148:       TaoSetHessian(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
149:   } else {
150:       TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user);
151:   }

153:   /* Set Variable bounds */
154:   MSA_Plate(xl,xu,(void*)&user);
155:   TaoSetVariableBounds(tao,xl,xu);

157:   /* Check for any tao command line options */
158:   TaoSetFromOptions(tao);

160:   /* SOLVE THE APPLICATION */
161:   TaoSolve(tao);

163:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

165:   /* Free TAO data structures */
166:   TaoDestroy(&tao);

168:   /* Free PETSc data structures */
169:   VecDestroy(&x);
170:   VecDestroy(&xl);
171:   VecDestroy(&xu);
172:   MatDestroy(&user.H);
173:   VecDestroy(&user.localX);
174:   VecDestroy(&user.localV);
175:   VecDestroy(&user.Bottom);
176:   VecDestroy(&user.Top);
177:   VecDestroy(&user.Left);
178:   VecDestroy(&user.Right);
179:   DMDestroy(&user.dm);
180:   if (flg) {
181:     MatDestroy(&H_shell);
182:   }
183:   PetscFinalize();
184:   return 0;
185: }

187: /*  FormFunctionGradient - Evaluates f(x) and gradient g(x).

189:     Input Parameters:
190: .   tao     - the Tao context
191: .   X      - input vector
192: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()

194:     Output Parameters:
195: .   fcn     - the function value
196: .   G      - vector containing the newly evaluated gradient

198:    Notes:
199:    In this case, we discretize the domain and Create triangles. The
200:    surface of each triangle is planar, whose surface area can be easily
201:    computed.  The total surface area is found by sweeping through the grid
202:    and computing the surface area of the two triangles that have their
203:    right angle at the grid point.  The diagonal line segments on the
204:    grid that define the triangles run from top left to lower right.
205:    The numbering of points starts at the lower left and runs left to
206:    right, then bottom to top.
207: */
208: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
209: {
210:   AppCtx         *user = (AppCtx *) userCtx;
211:   PetscInt       i,j,row;
212:   PetscInt       mx=user->mx, my=user->my;
213:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
214:   PetscReal      ft=0;
215:   PetscReal      zero=0.0;
216:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
217:   PetscReal      rhx=mx+1, rhy=my+1;
218:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
219:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
220:   PetscReal      *g, *x,*left,*right,*bottom,*top;
221:   Vec            localX = user->localX, localG = user->localV;

223:   /* Get local mesh boundaries */
224:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
225:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

227:   /* Scatter ghost points to local vector */
228:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
229:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

231:   /* Initialize vector to zero */
232:   VecSet(localG, zero);

234:   /* Get pointers to vector data */
235:   VecGetArray(localX,&x);
236:   VecGetArray(localG,&g);
237:   VecGetArray(user->Top,&top);
238:   VecGetArray(user->Bottom,&bottom);
239:   VecGetArray(user->Left,&left);
240:   VecGetArray(user->Right,&right);

242:   /* Compute function over the locally owned part of the mesh */
243:   for (j=ys; j<ys+ym; j++) {
244:     for (i=xs; i< xs+xm; i++) {
245:       row=(j-gys)*gxm + (i-gxs);

247:       xc = x[row];
248:       xlt=xrb=xl=xr=xb=xt=xc;

250:       if (i==0) { /* left side */
251:         xl= left[j-ys+1];
252:         xlt = left[j-ys+2];
253:       } else {
254:         xl = x[row-1];
255:       }

257:       if (j==0) { /* bottom side */
258:         xb=bottom[i-xs+1];
259:         xrb = bottom[i-xs+2];
260:       } else {
261:         xb = x[row-gxm];
262:       }

264:       if (i+1 == gxs+gxm) { /* right side */
265:         xr=right[j-ys+1];
266:         xrb = right[j-ys];
267:       } else {
268:         xr = x[row+1];
269:       }

271:       if (j+1==gys+gym) { /* top side */
272:         xt=top[i-xs+1];
273:         xlt = top[i-xs];
274:       }else {
275:         xt = x[row+gxm];
276:       }

278:       if (i>gxs && j+1<gys+gym) {
279:         xlt = x[row-1+gxm];
280:       }
281:       if (j>gys && i+1<gxs+gxm) {
282:         xrb = x[row+1-gxm];
283:       }

285:       d1 = (xc-xl);
286:       d2 = (xc-xr);
287:       d3 = (xc-xt);
288:       d4 = (xc-xb);
289:       d5 = (xr-xrb);
290:       d6 = (xrb-xb);
291:       d7 = (xlt-xl);
292:       d8 = (xt-xlt);

294:       df1dxc = d1*hydhx;
295:       df2dxc = (d1*hydhx + d4*hxdhy);
296:       df3dxc = d3*hxdhy;
297:       df4dxc = (d2*hydhx + d3*hxdhy);
298:       df5dxc = d2*hydhx;
299:       df6dxc = d4*hxdhy;

301:       d1 *= rhx;
302:       d2 *= rhx;
303:       d3 *= rhy;
304:       d4 *= rhy;
305:       d5 *= rhy;
306:       d6 *= rhx;
307:       d7 *= rhy;
308:       d8 *= rhx;

310:       f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
311:       f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
312:       f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
313:       f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
314:       f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
315:       f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);

317:       ft = ft + (f2 + f4);

319:       df1dxc /= f1;
320:       df2dxc /= f2;
321:       df3dxc /= f3;
322:       df4dxc /= f4;
323:       df5dxc /= f5;
324:       df6dxc /= f6;

326:       g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;

328:     }
329:   }

331:   /* Compute triangular areas along the border of the domain. */
332:   if (xs==0) { /* left side */
333:     for (j=ys; j<ys+ym; j++) {
334:       d3=(left[j-ys+1] - left[j-ys+2])*rhy;
335:       d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
336:       ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
337:     }
338:   }
339:   if (ys==0) { /* bottom side */
340:     for (i=xs; i<xs+xm; i++) {
341:       d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
342:       d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
343:       ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
344:     }
345:   }

347:   if (xs+xm==mx) { /* right side */
348:     for (j=ys; j< ys+ym; j++) {
349:       d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
350:       d4=(right[j-ys]-right[j-ys+1])*rhy;
351:       ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
352:     }
353:   }
354:   if (ys+ym==my) { /* top side */
355:     for (i=xs; i<xs+xm; i++) {
356:       d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
357:       d4=(top[i-xs+1] - top[i-xs])*rhx;
358:       ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
359:     }
360:   }

362:   if (ys==0 && xs==0) {
363:     d1=(left[0]-left[1])*rhy;
364:     d2=(bottom[0]-bottom[1])*rhx;
365:     ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
366:   }
367:   if (ys+ym == my && xs+xm == mx) {
368:     d1=(right[ym+1] - right[ym])*rhy;
369:     d2=(top[xm+1] - top[xm])*rhx;
370:     ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
371:   }

373:   ft=ft*area;
374:   MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);

376:   /* Restore vectors */
377:   VecRestoreArray(localX,&x);
378:   VecRestoreArray(localG,&g);
379:   VecRestoreArray(user->Left,&left);
380:   VecRestoreArray(user->Top,&top);
381:   VecRestoreArray(user->Bottom,&bottom);
382:   VecRestoreArray(user->Right,&right);

384:   /* Scatter values to global vector */
385:   DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
386:   DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);

388:   PetscLogFlops(70.0*xm*ym);

390:   return 0;
391: }

393: /* ------------------------------------------------------------------- */
394: /*
395:    FormHessian - Evaluates Hessian matrix.

397:    Input Parameters:
398: .  tao  - the Tao context
399: .  x    - input vector
400: .  ptr  - optional user-defined context, as set by TaoSetHessian()

402:    Output Parameters:
403: .  A    - Hessian matrix
404: .  B    - optionally different preconditioning matrix

406:    Notes:
407:    Due to mesh point reordering with DMs, we must always work
408:    with the local mesh points, and then transform them to the new
409:    global numbering with the local-to-global mapping.  We cannot work
410:    directly with the global numbers for the original uniprocessor mesh!

412:    Two methods are available for imposing this transformation
413:    when setting matrix entries:
414:      (A) MatSetValuesLocal(), using the local ordering (including
415:          ghost points!)
416:          - Do the following two steps once, before calling TaoSolve()
417:            - Use DMGetISLocalToGlobalMapping() to extract the
418:              local-to-global map from the DM
419:            - Associate this map with the matrix by calling
420:              MatSetLocalToGlobalMapping()
421:          - Then set matrix entries using the local ordering
422:            by calling MatSetValuesLocal()
423:      (B) MatSetValues(), using the global ordering
424:          - Use DMGetGlobalIndices() to extract the local-to-global map
425:          - Then apply this map explicitly yourself
426:          - Set matrix entries using the global ordering by calling
427:            MatSetValues()
428:    Option (A) seems cleaner/easier in many cases, and is the procedure
429:    used in this example.
430: */
431: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
432: {
433:   AppCtx         *user = (AppCtx *) ptr;
434:   PetscInt       i,j,k,row;
435:   PetscInt       mx=user->mx, my=user->my;
436:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
437:   PetscReal      hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
438:   PetscReal      rhx=mx+1, rhy=my+1;
439:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
440:   PetscReal      hl,hr,ht,hb,hc,htl,hbr;
441:   PetscReal      *x,*left,*right,*bottom,*top;
442:   PetscReal      v[7];
443:   Vec            localX = user->localX;
444:   PetscBool      assembled;

446:   /* Set various matrix options */
447:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

449:   /* Initialize matrix entries to zero */
450:   MatAssembled(Hessian,&assembled);
451:   if (assembled) MatZeroEntries(Hessian);

453:   /* Get local mesh boundaries */
454:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
455:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

457:   /* Scatter ghost points to local vector */
458:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
459:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

461:   /* Get pointers to vector data */
462:   VecGetArray(localX,&x);
463:   VecGetArray(user->Top,&top);
464:   VecGetArray(user->Bottom,&bottom);
465:   VecGetArray(user->Left,&left);
466:   VecGetArray(user->Right,&right);

468:   /* Compute Hessian over the locally owned part of the mesh */

470:   for (i=xs; i< xs+xm; i++) {

472:     for (j=ys; j<ys+ym; j++) {

474:       row=(j-gys)*gxm + (i-gxs);

476:       xc = x[row];
477:       xlt=xrb=xl=xr=xb=xt=xc;

479:       /* Left side */
480:       if (i==gxs) {
481:         xl= left[j-ys+1];
482:         xlt = left[j-ys+2];
483:       } else {
484:         xl = x[row-1];
485:       }

487:       if (j==gys) {
488:         xb=bottom[i-xs+1];
489:         xrb = bottom[i-xs+2];
490:       } else {
491:         xb = x[row-gxm];
492:       }

494:       if (i+1 == gxs+gxm) {
495:         xr=right[j-ys+1];
496:         xrb = right[j-ys];
497:       } else {
498:         xr = x[row+1];
499:       }

501:       if (j+1==gys+gym) {
502:         xt=top[i-xs+1];
503:         xlt = top[i-xs];
504:       }else {
505:         xt = x[row+gxm];
506:       }

508:       if (i>gxs && j+1<gys+gym) {
509:         xlt = x[row-1+gxm];
510:       }
511:       if (j>gys && i+1<gxs+gxm) {
512:         xrb = x[row+1-gxm];
513:       }

515:       d1 = (xc-xl)*rhx;
516:       d2 = (xc-xr)*rhx;
517:       d3 = (xc-xt)*rhy;
518:       d4 = (xc-xb)*rhy;
519:       d5 = (xrb-xr)*rhy;
520:       d6 = (xrb-xb)*rhx;
521:       d7 = (xlt-xl)*rhy;
522:       d8 = (xlt-xt)*rhx;

524:       f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
525:       f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
526:       f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
527:       f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
528:       f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
529:       f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);

531:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
532:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
533:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
534:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
535:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
536:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
537:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
538:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

543:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
544:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
545:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
546:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

548:       hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5;  hc*=0.5;

550:       k=0;
551:       if (j>0) {
552:         v[k]=hb; col[k]=row - gxm; k++;
553:       }

555:       if (j>0 && i < mx -1) {
556:         v[k]=hbr; col[k]=row - gxm+1; k++;
557:       }

559:       if (i>0) {
560:         v[k]= hl; col[k]=row - 1; k++;
561:       }

563:       v[k]= hc; col[k]=row; k++;

565:       if (i < mx-1) {
566:         v[k]= hr; col[k]=row+1; k++;
567:       }

569:       if (i>0 && j < my-1) {
570:         v[k]= htl; col[k] = row+gxm-1; k++;
571:       }

573:       if (j < my-1) {
574:         v[k]= ht; col[k] = row+gxm; k++;
575:       }

577:       /*
578:          Set matrix values using local numbering, which was defined
579:          earlier, in the main routine.
580:       */
581:       MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);

583:     }
584:   }

586:   /* Restore vectors */
587:   VecRestoreArray(localX,&x);
588:   VecRestoreArray(user->Left,&left);
589:   VecRestoreArray(user->Top,&top);
590:   VecRestoreArray(user->Bottom,&bottom);
591:   VecRestoreArray(user->Right,&right);

593:   /* Assemble the matrix */
594:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
595:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

597:   PetscLogFlops(199.0*xm*ym);
598:   return 0;
599: }

601: /* ------------------------------------------------------------------- */
602: /*
603:    MSA_BoundaryConditions -  Calculates the boundary conditions for
604:    the region.

606:    Input Parameter:
607: .  user - user-defined application context

609:    Output Parameter:
610: .  user - user-defined application context
611: */
612: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
613: {
614:   PetscInt   i,j,k,maxits=5,limit=0;
615:   PetscInt   xs,ys,xm,ym,gxs,gys,gxm,gym;
616:   PetscInt   mx=user->mx,my=user->my;
617:   PetscInt   bsize=0, lsize=0, tsize=0, rsize=0;
618:   PetscReal  one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
619:   PetscReal  fnorm,det,hx,hy,xt=0,yt=0;
620:   PetscReal  u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
621:   PetscReal  b=-0.5, t=0.5, l=-0.5, r=0.5;
622:   PetscReal  *boundary;
623:   PetscBool  flg;
624:   Vec        Bottom,Top,Right,Left;

626:   /* Get local mesh boundaries */
627:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
628:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

630:   bsize=xm+2;
631:   lsize=ym+2;
632:   rsize=ym+2;
633:   tsize=xm+2;

635:   VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
636:   VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
637:   VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
638:   VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);

640:   user->Top=Top;
641:   user->Left=Left;
642:   user->Bottom=Bottom;
643:   user->Right=Right;

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

647:   for (j=0; j<4; j++) {
648:     if (j==0) {
649:       yt=b;
650:       xt=l+hx*xs;
651:       limit=bsize;
652:       VecGetArray(Bottom,&boundary);
653:     } else if (j==1) {
654:       yt=t;
655:       xt=l+hx*xs;
656:       limit=tsize;
657:       VecGetArray(Top,&boundary);
658:     } else if (j==2) {
659:       yt=b+hy*ys;
660:       xt=l;
661:       limit=lsize;
662:       VecGetArray(Left,&boundary);
663:     } else if (j==3) {
664:       yt=b+hy*ys;
665:       xt=r;
666:       limit=rsize;
667:       VecGetArray(Right,&boundary);
668:     }

670:     for (i=0; i<limit; i++) {
671:       u1=xt;
672:       u2=-yt;
673:       for (k=0; k<maxits; k++) {
674:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
675:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
676:         fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
677:         if (fnorm <= tol) break;
678:         njac11=one+u2*u2-u1*u1;
679:         njac12=two*u1*u2;
680:         njac21=-two*u1*u2;
681:         njac22=-one - u1*u1 + u2*u2;
682:         det = njac11*njac22-njac21*njac12;
683:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
684:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
685:       }

687:       boundary[i]=u1*u1-u2*u2;
688:       if (j==0 || j==1) {
689:         xt=xt+hx;
690:       } else if (j==2 || j==3) {
691:         yt=yt+hy;
692:       }
693:     }
694:     if (j==0) {
695:       VecRestoreArray(Bottom,&boundary);
696:     } else if (j==1) {
697:       VecRestoreArray(Top,&boundary);
698:     } else if (j==2) {
699:       VecRestoreArray(Left,&boundary);
700:     } else if (j==3) {
701:       VecRestoreArray(Right,&boundary);
702:     }
703:   }

705:   /* Scale the boundary if desired */

707:   PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
708:   if (flg) {
709:     VecScale(Bottom, scl);
710:   }
711:   PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
712:   if (flg) {
713:     VecScale(Top, scl);
714:   }
715:   PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
716:   if (flg) {
717:     VecScale(Right, scl);
718:   }

720:   PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
721:   if (flg) {
722:     VecScale(Left, scl);
723:   }
724:   return 0;
725: }

727: /* ------------------------------------------------------------------- */
728: /*
729:    MSA_Plate -  Calculates an obstacle for surface to stretch over.

731:    Input Parameter:
732: .  user - user-defined application context

734:    Output Parameter:
735: .  user - user-defined application context
736: */
737: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx)
738: {
739:   AppCtx         *user=(AppCtx *)ctx;
740:   PetscInt       i,j,row;
741:   PetscInt       xs,ys,xm,ym;
742:   PetscInt       mx=user->mx, my=user->my, bmy, bmx;
743:   PetscReal      t1,t2,t3;
744:   PetscReal      *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
745:   PetscBool      cylinder;

747:   user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
748:   user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
749:   bmy=user->bmy; bmx=user->bmx;

751:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);

753:   VecSet(XL, lb);
754:   VecSet(XU, ub);

756:   VecGetArray(XL,&xl);

758:   PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);
759:   /* Compute the optional lower box */
760:   if (cylinder) {
761:     for (i=xs; i< xs+xm; i++) {
762:       for (j=ys; j<ys+ym; j++) {
763:         row=(j-ys)*xm + (i-xs);
764:         t1=(2.0*i-mx)*bmy;
765:         t2=(2.0*j-my)*bmx;
766:         t3=bmx*bmx*bmy*bmy;
767:         if (t1*t1 + t2*t2 <= t3) {
768:           xl[row] = user->bheight;
769:         }
770:       }
771:     }
772:   } else {
773:     /* Compute the optional lower box */
774:     for (i=xs; i< xs+xm; i++) {
775:       for (j=ys; j<ys+ym; j++) {
776:         row=(j-ys)*xm + (i-xs);
777:         if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
778:             j>=(my-bmy)/2 && j<my-(my-bmy)/2) {
779:           xl[row] = user->bheight;
780:         }
781:       }
782:     }
783:   }
784:     VecRestoreArray(XL,&xl);

786:   return 0;
787: }

789: /* ------------------------------------------------------------------- */
790: /*
791:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

793:    Input Parameters:
794: .  user - user-defined application context
795: .  X - vector for initial guess

797:    Output Parameters:
798: .  X - newly computed initial guess
799: */
800: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
801: {
802:   PetscInt       start=-1,i,j;
803:   PetscReal      zero=0.0;
804:   PetscBool      flg;

806:   PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
807:   if (flg && start==0) { /* The zero vector is reasonable */
808:     VecSet(X, zero);
809:   } else if (flg && start>0) { /* Try a random start between -0.5 and 0.5 */
810:     PetscRandom rctx;  PetscReal np5=-0.5;

812:     PetscRandomCreate(MPI_COMM_WORLD,&rctx);
813:     for (i=0; i<start; i++) {
814:       VecSetRandom(X, rctx);
815:     }
816:     PetscRandomDestroy(&rctx);
817:     VecShift(X, np5);

819:   } else { /* Take an average of the boundary conditions */

821:     PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
822:     PetscInt mx=user->mx,my=user->my;
823:     PetscReal *x,*left,*right,*bottom,*top;
824:     Vec    localX = user->localX;

826:     /* Get local mesh boundaries */
827:     DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
828:     DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

830:     /* Get pointers to vector data */
831:     VecGetArray(user->Top,&top);
832:     VecGetArray(user->Bottom,&bottom);
833:     VecGetArray(user->Left,&left);
834:     VecGetArray(user->Right,&right);

836:     VecGetArray(localX,&x);
837:     /* Perform local computations */
838:     for (j=ys; j<ys+ym; j++) {
839:       for (i=xs; i< xs+xm; i++) {
840:         row=(j-gys)*gxm + (i-gxs);
841:         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;
842:       }
843:     }

845:     /* Restore vectors */
846:     VecRestoreArray(localX,&x);

848:     VecRestoreArray(user->Left,&left);
849:     VecRestoreArray(user->Top,&top);
850:     VecRestoreArray(user->Bottom,&bottom);
851:     VecRestoreArray(user->Right,&right);

853:     /* Scatter values into global vector */
854:     DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
855:     DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);

857:   }
858:   return 0;
859: }

861: /* For testing matrix free submatrices */
862: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
863: {
864:   AppCtx         *user = (AppCtx*)ptr;
865:   FormHessian(tao,x,user->H,user->H,ptr);
866:   return 0;
867: }
868: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
869: {
870:   void           *ptr;
871:   AppCtx         *user;
872:   MatShellGetContext(H_shell,&ptr);
873:   user = (AppCtx*)ptr;
874:   MatMult(user->H,X,Y);
875:   return 0;
876: }

878: /*TEST

880:    build:
881:       requires: !complex

883:    test:
884:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
885:       requires: !single

887:    test:
888:       suffix: 2
889:       nsize: 2
890:       args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
891:       requires: !single

893:    test:
894:       suffix: 3
895:       nsize: 3
896:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
897:       requires: !single

899:    test:
900:       suffix: 4
901:       nsize: 3
902:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
903:       requires: !single

905:    test:
906:       suffix: 5
907:       nsize: 3
908:       args: -tao_smonitor -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
909:       requires: !single

911:    test:
912:       suffix: 6
913:       nsize: 3
914:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
915:       requires: !single

917:    test:
918:       suffix: 7
919:       nsize: 3
920:       args: -tao_smonitor -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
921:       requires: !single

923:    test:
924:       suffix: 8
925:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
926:       requires: !single

928:    test:
929:       suffix: 9
930:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
931:       requires: !single

933:    test:
934:       suffix: 10
935:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
936:       requires: !single

938:    test:
939:       suffix: 11
940:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
941:       requires: !single

943:    test:
944:       suffix: 12
945:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
946:       requires: !single

948:    test:
949:       suffix: 13
950:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
951:       requires: !single

953:    test:
954:       suffix: 14
955:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
956:       requires: !single

958:    test:
959:       suffix: 15
960:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
961:       requires: !single

963:    test:
964:      suffix: 16
965:      args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
966:      requires: !single

968:    test:
969:      suffix: 17
970:      args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
971:      requires: !single

973:    test:
974:      suffix: 18
975:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
976:      requires: !single

978:    test:
979:      suffix: 19
980:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
981:      requires: !single

983:    test:
984:      suffix: 20
985:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian

987: TEST*/