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[] =
 11: "This example demonstrates use of the TAO package to \n\
 12: solve an unconstrained minimization problem.  This example is based on a \n\
 13: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain and \n\
 14: boundary values along the edges of the domain, the objective is to find the\n\
 15: surface with the minimal area that satisfies the boundary conditions.\n\
 16: The command line options are:\n\
 17:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 18:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 19:   -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
 20:                for an average of the boundary conditions\n\n";

 22: /*T
 23:    Concepts: TAO^Solving an unconstrained minimization problem
 24:    Routines: TaoCreate(); TaoSetType();
 25:    Routines: TaoSetSolution();
 26:    Routines: TaoSetObjectiveAndGradient();
 27:    Routines: TaoSetHessian(); TaoSetFromOptions();
 28:    Routines: TaoSetMonitor();
 29:    Routines: TaoSolve(); TaoView();
 30:    Routines: TaoDestroy();
 31:    Processors: n
 32: T*/

 34: /*
 35:    User-defined application context - contains data needed by the
 36:    application-provided call-back routines, FormFunctionGradient()
 37:    and FormHessian().
 38: */
 39: typedef struct {
 40:   PetscInt    mx, my;                 /* discretization in x, y directions */
 41:   PetscReal   *bottom, *top, *left, *right;             /* boundary values */
 42:   DM          dm;                      /* distributed array data structure */
 43:   Mat         H;                       /* Hessian */
 44: } AppCtx;

 46: /* -------- User-defined Routines --------- */

 48: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
 49: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
 50: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
 51: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*);
 52: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
 53: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
 54: PetscErrorCode My_Monitor(Tao, void *);

 56: int main(int argc, char **argv)
 57: {
 58:   PetscInt           Nx, Ny;              /* number of processors in x- and y- directions */
 59:   Vec                x;                   /* solution, gradient vectors */
 60:   PetscBool          flg, viewmat;        /* flags */
 61:   PetscBool          fddefault, fdcoloring;   /* flags */
 62:   Tao                tao;                 /* TAO solver context */
 63:   AppCtx             user;                /* user-defined work context */
 64:   ISColoring         iscoloring;
 65:   MatFDColoring      matfdcoloring;

 67:   /* Initialize TAO */
 68:   PetscInitialize(&argc, &argv,(char *)0,help);

 70:   /* Specify dimension of the problem */
 71:   user.mx = 10; user.my = 10;

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

 77:   PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
 78:   PetscPrintf(MPI_COMM_WORLD,"mx: %D     my: %D   \n\n",user.mx,user.my);

 80:   /* Let PETSc determine the vector distribution */
 81:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 83:   /* Create distributed array (DM) to manage parallel grid and vectors  */
 84:   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);
 85:   DMSetFromOptions(user.dm);
 86:   DMSetUp(user.dm);

 88:   /* Create TAO solver and set desired solution method.*/
 89:   TaoCreate(PETSC_COMM_WORLD,&tao);
 90:   TaoSetType(tao,TAOCG);

 92:   /*
 93:      Extract global vector from DA for the vector of variables --  PETSC routine
 94:      Compute the initial solution                              --  application specific, see below
 95:      Set this vector for use by TAO                            --  TAO routine
 96:   */
 97:   DMCreateGlobalVector(user.dm,&x);
 98:   MSA_BoundaryConditions(&user);
 99:   MSA_InitialPoint(&user,x);
100:   TaoSetSolution(tao,x);

102:   /*
103:      Initialize the Application context for use in function evaluations  --  application specific, see below.
104:      Set routines for function and gradient evaluation
105:   */
106:   TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);

108:   /*
109:      Given the command line arguments, calculate the hessian with either the user-
110:      provided function FormHessian, or the default finite-difference driven Hessian
111:      functions
112:   */
113:   PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault);
114:   PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);

116:   /*
117:      Create a matrix data structure to store the Hessian and set
118:      the Hessian evalution routine.
119:      Set the matrix structure to be used for Hessian evalutions
120:   */
121:   DMCreateMatrix(user.dm,&user.H);
122:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);

124:   if (fdcoloring) {
125:     DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
126:     MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
127:     ISColoringDestroy(&iscoloring);
128:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
129:     MatFDColoringSetFromOptions(matfdcoloring);
130:     TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
131:   } else if (fddefault) {
132:     TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
133:   } else {
134:     TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user);
135:   }

137:   /*
138:      If my_monitor option is in command line, then use the user-provided
139:      monitoring function
140:   */
141:   PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
142:   if (viewmat) {
143:     TaoSetMonitor(tao,My_Monitor,NULL,NULL);
144:   }

146:   /* Check for any tao command line options */
147:   TaoSetFromOptions(tao);

149:   /* SOLVE THE APPLICATION */
150:   TaoSolve(tao);

152:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

154:   /* Free TAO data structures */
155:   TaoDestroy(&tao);

157:   /* Free PETSc data structures */
158:   VecDestroy(&x);
159:   MatDestroy(&user.H);
160:   if (fdcoloring) {
161:     MatFDColoringDestroy(&matfdcoloring);
162:   }
163:   PetscFree(user.bottom);
164:   PetscFree(user.top);
165:   PetscFree(user.left);
166:   PetscFree(user.right);
167:   DMDestroy(&user.dm);
168:   PetscFinalize();
169:   return 0;
170: }

172: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
173: {
174:   PetscReal      fcn;

176:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
177:   return 0;
178: }

180: /* -------------------------------------------------------------------- */
181: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

183:     Input Parameters:
184: .   tao     - the Tao context
185: .   XX      - input vector
186: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()

188:     Output Parameters:
189: .   fcn     - the newly evaluated function
190: .   GG       - vector containing the newly evaluated gradient
191: */
192: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx)
193: {
194:   AppCtx         *user = (AppCtx *) userCtx;
195:   PetscInt       i,j;
196:   PetscInt       mx=user->mx, my=user->my;
197:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
198:   PetscReal      ft=0.0;
199:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
200:   PetscReal      rhx=mx+1, rhy=my+1;
201:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
202:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
203:   PetscReal      **g, **x;
204:   Vec            localX;

206:   /* Get local mesh boundaries */
207:   DMGetLocalVector(user->dm,&localX);
208:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
209:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

211:   /* Scatter ghost points to local vector */
212:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
213:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

215:   /* Get pointers to vector data */
216:   DMDAVecGetArray(user->dm,localX,(void**)&x);
217:   DMDAVecGetArray(user->dm,G,(void**)&g);

219:   /* Compute function and gradient over the locally owned part of the mesh */
220:   for (j=ys; j<ys+ym; j++) {
221:     for (i=xs; i< xs+xm; i++) {

223:       xc = x[j][i];
224:       xlt=xrb=xl=xr=xb=xt=xc;

226:       if (i==0) { /* left side */
227:         xl= user->left[j-ys+1];
228:         xlt = user->left[j-ys+2];
229:       } else {
230:         xl = x[j][i-1];
231:       }

233:       if (j==0) { /* bottom side */
234:         xb=user->bottom[i-xs+1];
235:         xrb =user->bottom[i-xs+2];
236:       } else {
237:         xb = x[j-1][i];
238:       }

240:       if (i+1 == gxs+gxm) { /* right side */
241:         xr=user->right[j-ys+1];
242:         xrb = user->right[j-ys];
243:       } else {
244:         xr = x[j][i+1];
245:       }

247:       if (j+1==gys+gym) { /* top side */
248:         xt=user->top[i-xs+1];
249:         xlt = user->top[i-xs];
250:       }else {
251:         xt = x[j+1][i];
252:       }

254:       if (i>gxs && j+1<gys+gym) {
255:         xlt = x[j+1][i-1];
256:       }
257:       if (j>gys && i+1<gxs+gxm) {
258:         xrb = x[j-1][i+1];
259:       }

261:       d1 = (xc-xl);
262:       d2 = (xc-xr);
263:       d3 = (xc-xt);
264:       d4 = (xc-xb);
265:       d5 = (xr-xrb);
266:       d6 = (xrb-xb);
267:       d7 = (xlt-xl);
268:       d8 = (xt-xlt);

270:       df1dxc = d1*hydhx;
271:       df2dxc = (d1*hydhx + d4*hxdhy);
272:       df3dxc = d3*hxdhy;
273:       df4dxc = (d2*hydhx + d3*hxdhy);
274:       df5dxc = d2*hydhx;
275:       df6dxc = d4*hxdhy;

277:       d1 *= rhx;
278:       d2 *= rhx;
279:       d3 *= rhy;
280:       d4 *= rhy;
281:       d5 *= rhy;
282:       d6 *= rhx;
283:       d7 *= rhy;
284:       d8 *= rhx;

286:       f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
287:       f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
288:       f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
289:       f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
290:       f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
291:       f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);

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

295:       df1dxc /= f1;
296:       df2dxc /= f2;
297:       df3dxc /= f3;
298:       df4dxc /= f4;
299:       df5dxc /= f5;
300:       df6dxc /= f6;

302:       g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;

304:     }
305:   }

307:   /* Compute triangular areas along the border of the domain. */
308:   if (xs==0) { /* left side */
309:     for (j=ys; j<ys+ym; j++) {
310:       d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
311:       d2=(user->left[j-ys+1] - x[j][0]) *rhx;
312:       ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
313:     }
314:   }
315:   if (ys==0) { /* bottom side */
316:     for (i=xs; i<xs+xm; i++) {
317:       d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
318:       d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
319:       ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
320:     }
321:   }

323:   if (xs+xm==mx) { /* right side */
324:     for (j=ys; j< ys+ym; j++) {
325:       d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
326:       d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
327:       ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
328:     }
329:   }
330:   if (ys+ym==my) { /* top side */
331:     for (i=xs; i<xs+xm; i++) {
332:       d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
333:       d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
334:       ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
335:     }
336:   }

338:   if (ys==0 && xs==0) {
339:     d1=(user->left[0]-user->left[1])/hy;
340:     d2=(user->bottom[0]-user->bottom[1])*rhx;
341:     ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
342:   }
343:   if (ys+ym == my && xs+xm == mx) {
344:     d1=(user->right[ym+1] - user->right[ym])*rhy;
345:     d2=(user->top[xm+1] - user->top[xm])*rhx;
346:     ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
347:   }

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

352:   /* Restore vectors */
353:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
354:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

356:   /* Scatter values to global vector */
357:   DMRestoreLocalVector(user->dm,&localX);
358:   PetscLogFlops(67.0*xm*ym);
359:   return 0;
360: }

362: /* ------------------------------------------------------------------- */
363: /*
364:    FormHessian - Evaluates Hessian matrix.

366:    Input Parameters:
367: .  tao  - the Tao context
368: .  x    - input vector
369: .  ptr  - optional user-defined context, as set by TaoSetHessian()

371:    Output Parameters:
372: .  H    - Hessian matrix
373: .  Hpre - optionally different preconditioning matrix
374: .  flg  - flag indicating matrix structure

376: */
377: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
378: {
379:   AppCtx         *user = (AppCtx *) ptr;

381:   /* Evaluate the Hessian entries*/
382:   QuadraticH(user,X,H);
383:   return 0;
384: }

386: /* ------------------------------------------------------------------- */
387: /*
388:    QuadraticH - Evaluates Hessian matrix.

390:    Input Parameters:
391: .  user - user-defined context, as set by TaoSetHessian()
392: .  X    - input vector

394:    Output Parameter:
395: .  H    - Hessian matrix
396: */
397: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
398: {
399:   PetscInt i,j,k;
400:   PetscInt mx=user->mx, my=user->my;
401:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
402:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
403:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
404:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
405:   PetscReal **x, v[7];
406:   MatStencil col[7],row;
407:   Vec    localX;
408:   PetscBool assembled;

410:   /* Get local mesh boundaries */
411:   DMGetLocalVector(user->dm,&localX);

413:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
414:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

416:   /* Scatter ghost points to local vector */
417:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
418:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

420:   /* Get pointers to vector data */
421:   DMDAVecGetArray(user->dm,localX,(void**)&x);

423:   /* Initialize matrix entries to zero */
424:   MatAssembled(Hessian,&assembled);
425:   if (assembled) MatZeroEntries(Hessian);

427:   /* Set various matrix options */
428:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

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

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

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

436:       xc = x[j][i];
437:       xlt=xrb=xl=xr=xb=xt=xc;

439:       /* Left side */
440:       if (i==0) {
441:         xl  = user->left[j-ys+1];
442:         xlt = user->left[j-ys+2];
443:       } else {
444:         xl  = x[j][i-1];
445:       }

447:       if (j==0) {
448:         xb  = user->bottom[i-xs+1];
449:         xrb = user->bottom[i-xs+2];
450:       } else {
451:         xb  = x[j-1][i];
452:       }

454:       if (i+1 == mx) {
455:         xr  = user->right[j-ys+1];
456:         xrb = user->right[j-ys];
457:       } else {
458:         xr  = x[j][i+1];
459:       }

461:       if (j+1==my) {
462:         xt  = user->top[i-xs+1];
463:         xlt = user->top[i-xs];
464:       }else {
465:         xt  = x[j+1][i];
466:       }

468:       if (i>0 && j+1<my) {
469:         xlt = x[j+1][i-1];
470:       }
471:       if (j>0 && i+1<mx) {
472:         xrb = x[j-1][i+1];
473:       }

475:       d1 = (xc-xl)/hx;
476:       d2 = (xc-xr)/hx;
477:       d3 = (xc-xt)/hy;
478:       d4 = (xc-xb)/hy;
479:       d5 = (xrb-xr)/hy;
480:       d6 = (xrb-xb)/hx;
481:       d7 = (xlt-xl)/hy;
482:       d8 = (xlt-xt)/hx;

484:       f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
485:       f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
486:       f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
487:       f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
488:       f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
489:       f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);

491:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
492:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
493:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
494:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
495:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
496:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
497:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
498:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

503:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
504:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
505:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
506:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

508:       hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0;  hc/=2.0;

510:       row.j = j; row.i = i;
511:       k=0;
512:       if (j>0) {
513:         v[k]=hb;
514:         col[k].j = j - 1; col[k].i = i;
515:         k++;
516:       }

518:       if (j>0 && i < mx -1) {
519:         v[k]=hbr;
520:         col[k].j = j - 1; col[k].i = i+1;
521:         k++;
522:       }

524:       if (i>0) {
525:         v[k]= hl;
526:         col[k].j = j; col[k].i = i-1;
527:         k++;
528:       }

530:       v[k]= hc;
531:       col[k].j = j; col[k].i = i;
532:       k++;

534:       if (i < mx-1) {
535:         v[k]= hr;
536:         col[k].j = j; col[k].i = i+1;
537:         k++;
538:       }

540:       if (i>0 && j < my-1) {
541:         v[k]= htl;
542:         col[k].j = j+1; col[k].i = i-1;
543:         k++;
544:       }

546:       if (j < my-1) {
547:         v[k]= ht;
548:         col[k].j = j+1; col[k].i = i;
549:         k++;
550:       }

552:       /*
553:          Set matrix values using local numbering, which was defined
554:          earlier, in the main routine.
555:       */
556:       MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
557:     }
558:   }

560:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
561:   DMRestoreLocalVector(user->dm,&localX);

563:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
564:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

566:   PetscLogFlops(199.0*xm*ym);
567:   return 0;
568: }

570: /* ------------------------------------------------------------------- */
571: /*
572:    MSA_BoundaryConditions -  Calculates the boundary conditions for
573:    the region.

575:    Input Parameter:
576: .  user - user-defined application context

578:    Output Parameter:
579: .  user - user-defined application context
580: */
581: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
582: {
583:   PetscInt       i,j,k,limit=0,maxits=5;
584:   PetscInt       xs,ys,xm,ym,gxs,gys,gxm,gym;
585:   PetscInt       mx=user->mx,my=user->my;
586:   PetscInt       bsize=0, lsize=0, tsize=0, rsize=0;
587:   PetscReal      one=1.0, two=2.0, three=3.0, tol=1e-10;
588:   PetscReal      fnorm,det,hx,hy,xt=0,yt=0;
589:   PetscReal      u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
590:   PetscReal      b=-0.5, t=0.5, l=-0.5, r=0.5;
591:   PetscReal      *boundary;
592:   PetscBool      flg;

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

598:   bsize=xm+2;
599:   lsize=ym+2;
600:   rsize=ym+2;
601:   tsize=xm+2;

603:   PetscMalloc1(bsize,&user->bottom);
604:   PetscMalloc1(tsize,&user->top);
605:   PetscMalloc1(lsize,&user->left);
606:   PetscMalloc1(rsize,&user->right);

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

610:   for (j=0; j<4; j++) {
611:     if (j==0) {
612:       yt=b;
613:       xt=l+hx*xs;
614:       limit=bsize;
615:       boundary=user->bottom;
616:     } else if (j==1) {
617:       yt=t;
618:       xt=l+hx*xs;
619:       limit=tsize;
620:       boundary=user->top;
621:     } else if (j==2) {
622:       yt=b+hy*ys;
623:       xt=l;
624:       limit=lsize;
625:       boundary=user->left;
626:     } else { /* if (j==3) */
627:       yt=b+hy*ys;
628:       xt=r;
629:       limit=rsize;
630:       boundary=user->right;
631:     }

633:     for (i=0; i<limit; i++) {
634:       u1=xt;
635:       u2=-yt;
636:       for (k=0; k<maxits; k++) {
637:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
638:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
639:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
640:         if (fnorm <= tol) break;
641:         njac11=one+u2*u2-u1*u1;
642:         njac12=two*u1*u2;
643:         njac21=-two*u1*u2;
644:         njac22=-one - u1*u1 + u2*u2;
645:         det = njac11*njac22-njac21*njac12;
646:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
647:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
648:       }

650:       boundary[i]=u1*u1-u2*u2;
651:       if (j==0 || j==1) {
652:         xt=xt+hx;
653:       } else { /*  if (j==2 || j==3) */
654:         yt=yt+hy;
655:       }

657:     }

659:   }

661:   /* Scale the boundary if desired */
662:   if (1==1) {
663:     PetscReal scl = 1.0;

665:     PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
666:     if (flg) {
667:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
668:     }

670:     PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
671:     if (flg) {
672:       for (i=0;i<tsize;i++) user->top[i]*=scl;
673:     }

675:     PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
676:     if (flg) {
677:       for (i=0;i<rsize;i++) user->right[i]*=scl;
678:     }

680:     PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
681:     if (flg) {
682:       for (i=0;i<lsize;i++) user->left[i]*=scl;
683:     }
684:   }
685:   return 0;
686: }

688: /* ------------------------------------------------------------------- */
689: /*
690:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

692:    Input Parameters:
693: .  user - user-defined application context
694: .  X - vector for initial guess

696:    Output Parameters:
697: .  X - newly computed initial guess
698: */
699: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
700: {
701:   PetscInt       start2=-1,i,j;
702:   PetscReal      start1=0;
703:   PetscBool      flg1,flg2;

705:   PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
706:   PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);

708:   if (flg1) { /* The zero vector is reasonable */

710:     VecSet(X, start1);

712:   } else if (flg2 && start2>0) { /* Try a random start between -0.5 and 0.5 */
713:     PetscRandom rctx;  PetscReal np5=-0.5;

715:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
716:     VecSetRandom(X, rctx);
717:     PetscRandomDestroy(&rctx);
718:     VecShift(X, np5);

720:   } else { /* Take an average of the boundary conditions */
721:     PetscInt  xs,xm,ys,ym;
722:     PetscInt  mx=user->mx,my=user->my;
723:     PetscReal **x;

725:     /* Get local mesh boundaries */
726:     DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);

728:     /* Get pointers to vector data */
729:     DMDAVecGetArray(user->dm,X,(void**)&x);

731:     /* Perform local computations */
732:     for (j=ys; j<ys+ym; j++) {
733:       for (i=xs; i< xs+xm; i++) {
734:         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;
735:       }
736:     }
737:     DMDAVecRestoreArray(user->dm,X,(void**)&x);
738:     PetscLogFlops(9.0*xm*ym);
739:   }
740:   return 0;
741: }

743: /*-----------------------------------------------------------------------*/
744: PetscErrorCode My_Monitor(Tao tao, void *ctx)
745: {
746:   Vec            X;

748:   TaoGetSolution(tao,&X);
749:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
750:   return 0;
751: }

753: /*TEST

755:    build:
756:       requires: !complex

758:    test:
759:       args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
760:       requires: !single

762:    test:
763:       suffix: 2
764:       nsize: 2
765:       args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
766:       filter: grep -v "nls ksp"
767:       requires: !single

769:    test:
770:       suffix: 3
771:       nsize: 3
772:       args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
773:       requires: !single

775:    test:
776:       suffix: 5
777:       nsize: 2
778:       args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
779:       requires: !single

781: TEST*/