Actual source code: minsurf2.c

petsc-3.7.1 2016-05-15
Report Typos and Errors
  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: TaoSetInitialVector();
 26:    Routines: TaoSetObjectiveAndGradientRoutine();
 27:    Routines: TaoSetHessianRoutine(); 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;


 47: /* -------- User-defined Routines --------- */

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

 59: int main( int argc, char **argv )
 60: {
 61:   PetscErrorCode     ierr;                /* used to check for functions returning nonzeros */
 62:   PetscInt           Nx, Ny;              /* number of processors in x- and y- directions */
 63:   Vec                x;                   /* solution, gradient vectors */
 64:   PetscBool          flg, viewmat;        /* flags */
 65:   PetscBool          fddefault, fdcoloring;   /* flags */
 66:   Tao                tao;                 /* TAO solver context */
 67:   AppCtx             user;                /* user-defined work context */
 68:   ISColoring         iscoloring;
 69:   MatFDColoring      matfdcoloring;

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

 74:   /* Specify dimension of the problem */
 75:   user.mx = 10; user.my = 10;

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

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


 85:   /* Let PETSc determine the vector distribution */
 86:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 88:   /* Create distributed array (DM) to manage parallel grid and vectors  */
 89:   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);


 92:   /* Create TAO solver and set desired solution method.*/
 93:   TaoCreate(PETSC_COMM_WORLD,&tao);
 94:   TaoSetType(tao,TAOCG);

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

106:   /*
107:      Initialize the Application context for use in function evaluations  --  application specific, see below.
108:      Set routines for function and gradient evaluation
109:   */
110:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);

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

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


129:   if (fdcoloring) {
130:     DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
131:     MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
132:     ISColoringDestroy(&iscoloring);
133:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
134:     MatFDColoringSetFromOptions(matfdcoloring);
135:     TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
136:   } else if (fddefault){
137:     TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
138:   } else {
139:     TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
140:   }

142:   /*
143:      If my_monitor option is in command line, then use the user-provided
144:      monitoring function
145:   */
146:   PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
147:   if (viewmat){
148:     TaoSetMonitor(tao,My_Monitor,NULL,NULL);
149:   }

151:   /* Check for any tao command line options */
152:   TaoSetFromOptions(tao);

154:   /* SOLVE THE APPLICATION */
155:   TaoSolve(tao);

157:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

159:   /* Free TAO data structures */
160:   TaoDestroy(&tao);

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

179: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
180: {
182:   PetscReal      fcn;

185:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
186:   return(0);
187: }

189: /* -------------------------------------------------------------------- */
192: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

194:     Input Parameters:
195: .   tao     - the Tao context
196: .   XX      - input vector
197: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

199:     Output Parameters:
200: .   fcn     - the newly evaluated function
201: .   GG       - vector containing the newly evaluated gradient
202: */
203: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){

205:   AppCtx         *user = (AppCtx *) userCtx;
207:   PetscInt       i,j;
208:   PetscInt       mx=user->mx, my=user->my;
209:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
210:   PetscReal      ft=0.0;
211:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
212:   PetscReal      rhx=mx+1, rhy=my+1;
213:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
214:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
215:   PetscReal      **g, **x;
216:   Vec            localX;

219:   /* Get local mesh boundaries */
220:   DMGetLocalVector(user->dm,&localX);
221:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
222:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

224:   /* Scatter ghost points to local vector */
225:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
226:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

228:   /* Get pointers to vector data */
229:   DMDAVecGetArray(user->dm,localX,(void**)&x);
230:   DMDAVecGetArray(user->dm,G,(void**)&g);

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

236:       xc = x[j][i];
237:       xlt=xrb=xl=xr=xb=xt=xc;

239:       if (i==0){ /* left side */
240:         xl= user->left[j-ys+1];
241:         xlt = user->left[j-ys+2];
242:       } else {
243:         xl = x[j][i-1];
244:       }

246:       if (j==0){ /* bottom side */
247:         xb=user->bottom[i-xs+1];
248:         xrb =user->bottom[i-xs+2];
249:       } else {
250:         xb = x[j-1][i];
251:       }

253:       if (i+1 == gxs+gxm){ /* right side */
254:         xr=user->right[j-ys+1];
255:         xrb = user->right[j-ys];
256:       } else {
257:         xr = x[j][i+1];
258:       }

260:       if (j+1==gys+gym){ /* top side */
261:         xt=user->top[i-xs+1];
262:         xlt = user->top[i-xs];
263:       }else {
264:         xt = x[j+1][i];
265:       }

267:       if (i>gxs && j+1<gys+gym){
268:         xlt = x[j+1][i-1];
269:       }
270:       if (j>gys && i+1<gxs+gxm){
271:         xrb = x[j-1][i+1];
272:       }

274:       d1 = (xc-xl);
275:       d2 = (xc-xr);
276:       d3 = (xc-xt);
277:       d4 = (xc-xb);
278:       d5 = (xr-xrb);
279:       d6 = (xrb-xb);
280:       d7 = (xlt-xl);
281:       d8 = (xt-xlt);

283:       df1dxc = d1*hydhx;
284:       df2dxc = ( d1*hydhx + d4*hxdhy );
285:       df3dxc = d3*hxdhy;
286:       df4dxc = ( d2*hydhx + d3*hxdhy );
287:       df5dxc = d2*hydhx;
288:       df6dxc = d4*hxdhy;

290:       d1 *= rhx;
291:       d2 *= rhx;
292:       d3 *= rhy;
293:       d4 *= rhy;
294:       d5 *= rhy;
295:       d6 *= rhx;
296:       d7 *= rhy;
297:       d8 *= rhx;

299:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
300:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
301:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
302:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
303:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
304:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);

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

308:       df1dxc /= f1;
309:       df2dxc /= f2;
310:       df3dxc /= f3;
311:       df4dxc /= f4;
312:       df5dxc /= f5;
313:       df6dxc /= f6;

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

317:     }
318:   }

320:   /* Compute triangular areas along the border of the domain. */
321:   if (xs==0){ /* left side */
322:     for (j=ys; j<ys+ym; j++){
323:       d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
324:       d2=(user->left[j-ys+1] - x[j][0]) *rhx;
325:       ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
326:     }
327:   }
328:   if (ys==0){ /* bottom side */
329:     for (i=xs; i<xs+xm; i++){
330:       d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
331:       d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
332:       ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
333:     }
334:   }

336:   if (xs+xm==mx){ /* right side */
337:     for (j=ys; j< ys+ym; j++){
338:       d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
339:       d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
340:       ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
341:     }
342:   }
343:   if (ys+ym==my){ /* top side */
344:     for (i=xs; i<xs+xm; i++){
345:       d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
346:       d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
347:       ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
348:     }
349:   }

351:   if (ys==0 && xs==0){
352:     d1=(user->left[0]-user->left[1])/hy;
353:     d2=(user->bottom[0]-user->bottom[1])*rhx;
354:     ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
355:   }
356:   if (ys+ym == my && xs+xm == mx){
357:     d1=(user->right[ym+1] - user->right[ym])*rhy;
358:     d2=(user->top[xm+1] - user->top[xm])*rhx;
359:     ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
360:   }

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

365:   /* Restore vectors */
366:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
367:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

369:   /* Scatter values to global vector */
370:   DMRestoreLocalVector(user->dm,&localX);
371:   PetscLogFlops(67*xm*ym);
372:   return(0);
373: }

375: /* ------------------------------------------------------------------- */
378: /*
379:    FormHessian - Evaluates Hessian matrix.

381:    Input Parameters:
382: .  tao  - the Tao context
383: .  x    - input vector
384: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

386:    Output Parameters:
387: .  H    - Hessian matrix
388: .  Hpre - optionally different preconditioning matrix
389: .  flg  - flag indicating matrix structure

391: */
392: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
393: {
395:   AppCtx         *user = (AppCtx *) ptr;

398:   /* Evaluate the Hessian entries*/
399:   QuadraticH(user,X,H);
400:   return(0);
401: }

403: /* ------------------------------------------------------------------- */
406: /*
407:    QuadraticH - Evaluates Hessian matrix.

409:    Input Parameters:
410: .  user - user-defined context, as set by TaoSetHessian()
411: .  X    - input vector

413:    Output Parameter:
414: .  H    - Hessian matrix
415: */
416: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
417: {
418:   PetscErrorCode    ierr;
419:   PetscInt i,j,k;
420:   PetscInt mx=user->mx, my=user->my;
421:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
422:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
423:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
424:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
425:   PetscReal **x, v[7];
426:   MatStencil col[7],row;
427:   Vec    localX;
428:   PetscBool assembled;

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

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

437:   /* Scatter ghost points to local vector */
438:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
439:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

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

444:   /* Initialize matrix entries to zero */
445:   MatAssembled(Hessian,&assembled);
446:   if (assembled){MatZeroEntries(Hessian);}


449:   /* Set various matrix options */
450:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

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

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

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

458:       xc = x[j][i];
459:       xlt=xrb=xl=xr=xb=xt=xc;

461:       /* Left side */
462:       if (i==0){
463:         xl  = user->left[j-ys+1];
464:         xlt = user->left[j-ys+2];
465:       } else {
466:         xl  = x[j][i-1];
467:       }

469:       if (j==0){
470:         xb  = user->bottom[i-xs+1];
471:         xrb = user->bottom[i-xs+2];
472:       } else {
473:         xb  = x[j-1][i];
474:       }

476:       if (i+1 == mx){
477:         xr  = user->right[j-ys+1];
478:         xrb = user->right[j-ys];
479:       } else {
480:         xr  = x[j][i+1];
481:       }

483:       if (j+1==my){
484:         xt  = user->top[i-xs+1];
485:         xlt = user->top[i-xs];
486:       }else {
487:         xt  = x[j+1][i];
488:       }

490:       if (i>0 && j+1<my){
491:         xlt = x[j+1][i-1];
492:       }
493:       if (j>0 && i+1<mx){
494:         xrb = x[j-1][i+1];
495:       }


498:       d1 = (xc-xl)/hx;
499:       d2 = (xc-xr)/hx;
500:       d3 = (xc-xt)/hy;
501:       d4 = (xc-xb)/hy;
502:       d5 = (xrb-xr)/hy;
503:       d6 = (xrb-xb)/hx;
504:       d7 = (xlt-xl)/hy;
505:       d8 = (xlt-xt)/hx;

507:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
508:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
509:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
510:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
511:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
512:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);


515:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
516:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
517:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
518:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
519:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
520:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
521:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
522:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

527:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
528:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
529:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
530:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

534:       row.j = j; row.i = i;
535:       k=0;
536:       if (j>0){
537:         v[k]=hb;
538:         col[k].j = j - 1; col[k].i = i;
539:         k++;
540:       }

542:       if (j>0 && i < mx -1){
543:         v[k]=hbr;
544:         col[k].j = j - 1; col[k].i = i+1;
545:         k++;
546:       }

548:       if (i>0){
549:         v[k]= hl;
550:         col[k].j = j; col[k].i = i-1;
551:         k++;
552:       }

554:       v[k]= hc;
555:       col[k].j = j; col[k].i = i;
556:       k++;

558:       if (i < mx-1 ){
559:         v[k]= hr;
560:         col[k].j = j; col[k].i = i+1;
561:         k++;
562:       }

564:       if (i>0 && j < my-1 ){
565:         v[k]= htl;
566:         col[k].j = j+1; col[k].i = i-1;
567:         k++;
568:       }

570:       if (j < my-1 ){
571:         v[k]= ht;
572:         col[k].j = j+1; col[k].i = i;
573:         k++;
574:       }

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

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

587:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
588:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

590:   PetscLogFlops(199*xm*ym);
591:   return(0);
592: }

594: /* ------------------------------------------------------------------- */
597: /*
598:    MSA_BoundaryConditions -  Calculates the boundary conditions for
599:    the region.

601:    Input Parameter:
602: .  user - user-defined application context

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

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

626:   bsize=xm+2;
627:   lsize=ym+2;
628:   rsize=ym+2;
629:   tsize=xm+2;

631:   PetscMalloc1(bsize,&user->bottom);
632:   PetscMalloc1(tsize,&user->top);
633:   PetscMalloc1(lsize,&user->left);
634:   PetscMalloc1(rsize,&user->right);

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

638:   for (j=0; j<4; j++){
639:     if (j==0){
640:       yt=b;
641:       xt=l+hx*xs;
642:       limit=bsize;
643:       boundary=user->bottom;
644:     } else if (j==1){
645:       yt=t;
646:       xt=l+hx*xs;
647:       limit=tsize;
648:       boundary=user->top;
649:     } else if (j==2){
650:       yt=b+hy*ys;
651:       xt=l;
652:       limit=lsize;
653:       boundary=user->left;
654:     } else { /* if (j==3) */
655:       yt=b+hy*ys;
656:       xt=r;
657:       limit=rsize;
658:       boundary=user->right;
659:     }

661:     for (i=0; i<limit; i++){
662:       u1=xt;
663:       u2=-yt;
664:       for (k=0; k<maxits; k++){
665:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
666:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
667:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
668:         if (fnorm <= tol) break;
669:         njac11=one+u2*u2-u1*u1;
670:         njac12=two*u1*u2;
671:         njac21=-two*u1*u2;
672:         njac22=-one - u1*u1 + u2*u2;
673:         det = njac11*njac22-njac21*njac12;
674:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
675:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
676:       }

678:       boundary[i]=u1*u1-u2*u2;
679:       if (j==0 || j==1) {
680:         xt=xt+hx;
681:       } else { /*  if (j==2 || j==3) */
682:         yt=yt+hy;
683:       }

685:     }

687:   }

689:   /* Scale the boundary if desired */
690:   if (1==1){
691:     PetscReal scl = 1.0;

693:     PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
694:     if (flg){
695:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
696:     }

698:     PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
699:     if (flg){
700:       for (i=0;i<tsize;i++) user->top[i]*=scl;
701:     }

703:     PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
704:     if (flg){
705:       for (i=0;i<rsize;i++) user->right[i]*=scl;
706:     }

708:     PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
709:     if (flg){
710:       for (i=0;i<lsize;i++) user->left[i]*=scl;
711:     }
712:   }
713:   return(0);
714: }

716: /* ------------------------------------------------------------------- */
719: /*
720:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

722:    Input Parameters:
723: .  user - user-defined application context
724: .  X - vector for initial guess

726:    Output Parameters:
727: .  X - newly computed initial guess
728: */
729: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
730: {
732:   PetscInt       start2=-1,i,j;
733:   PetscReal      start1=0;
734:   PetscBool      flg1,flg2;

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

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

742:     VecSet(X, start1);

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

747:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
748:     VecSetRandom(X, rctx);
749:     PetscRandomDestroy(&rctx);
750:     VecShift(X, np5);

752:   } else { /* Take an average of the boundary conditions */
753:     PetscInt  xs,xm,ys,ym;
754:     PetscInt  mx=user->mx,my=user->my;
755:     PetscReal **x;

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

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

763:     /* Perform local computations */
764:     for (j=ys; j<ys+ym; j++){
765:       for (i=xs; i< xs+xm; i++){
766:         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;
767:       }
768:     }
769:     DMDAVecRestoreArray(user->dm,X,(void**)&x);
770:     PetscLogFlops(9*xm*ym);
771:   }
772:   return(0);
773: }

775: /*-----------------------------------------------------------------------*/
778: PetscErrorCode My_Monitor(Tao tao, void *ctx)
779: {
781:   Vec            X;

784:   TaoGetSolutionVector(tao,&X);
785:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
786:   return(0);
787: }