Actual source code: minsurf2.c

petsc-master 2017-05-25
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 *);

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

 69:   /* Initialize TAO */
 70:   PetscInitialize( &argc, &argv,(char *)0,help );if (ierr) return ierr;

 72:   /* Specify dimension of the problem */
 73:   user.mx = 10; user.my = 10;

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

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


 83:   /* Let PETSc determine the vector distribution */
 84:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

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

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

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

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

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

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


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

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

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

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

156:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

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

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

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

182:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
183:   return(0);
184: }

186: /* -------------------------------------------------------------------- */
187: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

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

194:     Output Parameters:
195: .   fcn     - the newly evaluated function
196: .   GG       - vector containing the newly evaluated gradient
197: */
198: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){

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

214:   /* Get local mesh boundaries */
215:   DMGetLocalVector(user->dm,&localX);
216:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
217:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

219:   /* Scatter ghost points to local vector */
220:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
221:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

223:   /* Get pointers to vector data */
224:   DMDAVecGetArray(user->dm,localX,(void**)&x);
225:   DMDAVecGetArray(user->dm,G,(void**)&g);

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

231:       xc = x[j][i];
232:       xlt=xrb=xl=xr=xb=xt=xc;

234:       if (i==0){ /* left side */
235:         xl= user->left[j-ys+1];
236:         xlt = user->left[j-ys+2];
237:       } else {
238:         xl = x[j][i-1];
239:       }

241:       if (j==0){ /* bottom side */
242:         xb=user->bottom[i-xs+1];
243:         xrb =user->bottom[i-xs+2];
244:       } else {
245:         xb = x[j-1][i];
246:       }

248:       if (i+1 == gxs+gxm){ /* right side */
249:         xr=user->right[j-ys+1];
250:         xrb = user->right[j-ys];
251:       } else {
252:         xr = x[j][i+1];
253:       }

255:       if (j+1==gys+gym){ /* top side */
256:         xt=user->top[i-xs+1];
257:         xlt = user->top[i-xs];
258:       }else {
259:         xt = x[j+1][i];
260:       }

262:       if (i>gxs && j+1<gys+gym){
263:         xlt = x[j+1][i-1];
264:       }
265:       if (j>gys && i+1<gxs+gxm){
266:         xrb = x[j-1][i+1];
267:       }

269:       d1 = (xc-xl);
270:       d2 = (xc-xr);
271:       d3 = (xc-xt);
272:       d4 = (xc-xb);
273:       d5 = (xr-xrb);
274:       d6 = (xrb-xb);
275:       d7 = (xlt-xl);
276:       d8 = (xt-xlt);

278:       df1dxc = d1*hydhx;
279:       df2dxc = ( d1*hydhx + d4*hxdhy );
280:       df3dxc = d3*hxdhy;
281:       df4dxc = ( d2*hydhx + d3*hxdhy );
282:       df5dxc = d2*hydhx;
283:       df6dxc = d4*hxdhy;

285:       d1 *= rhx;
286:       d2 *= rhx;
287:       d3 *= rhy;
288:       d4 *= rhy;
289:       d5 *= rhy;
290:       d6 *= rhx;
291:       d7 *= rhy;
292:       d8 *= rhx;

294:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
295:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
296:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
297:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
298:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
299:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);

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

303:       df1dxc /= f1;
304:       df2dxc /= f2;
305:       df3dxc /= f3;
306:       df4dxc /= f4;
307:       df5dxc /= f5;
308:       df6dxc /= f6;

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

312:     }
313:   }

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

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

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

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

360:   /* Restore vectors */
361:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
362:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

364:   /* Scatter values to global vector */
365:   DMRestoreLocalVector(user->dm,&localX);
366:   PetscLogFlops(67*xm*ym);
367:   return(0);
368: }

370: /* ------------------------------------------------------------------- */
371: /*
372:    FormHessian - Evaluates Hessian matrix.

374:    Input Parameters:
375: .  tao  - the Tao context
376: .  x    - input vector
377: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

379:    Output Parameters:
380: .  H    - Hessian matrix
381: .  Hpre - optionally different preconditioning matrix
382: .  flg  - flag indicating matrix structure

384: */
385: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
386: {
388:   AppCtx         *user = (AppCtx *) ptr;

391:   /* Evaluate the Hessian entries*/
392:   QuadraticH(user,X,H);
393:   return(0);
394: }

396: /* ------------------------------------------------------------------- */
397: /*
398:    QuadraticH - Evaluates Hessian matrix.

400:    Input Parameters:
401: .  user - user-defined context, as set by TaoSetHessian()
402: .  X    - input vector

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

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

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

428:   /* Scatter ghost points to local vector */
429:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
430:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

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

435:   /* Initialize matrix entries to zero */
436:   MatAssembled(Hessian,&assembled);
437:   if (assembled){MatZeroEntries(Hessian);}


440:   /* Set various matrix options */
441:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

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

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

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

449:       xc = x[j][i];
450:       xlt=xrb=xl=xr=xb=xt=xc;

452:       /* Left side */
453:       if (i==0){
454:         xl  = user->left[j-ys+1];
455:         xlt = user->left[j-ys+2];
456:       } else {
457:         xl  = x[j][i-1];
458:       }

460:       if (j==0){
461:         xb  = user->bottom[i-xs+1];
462:         xrb = user->bottom[i-xs+2];
463:       } else {
464:         xb  = x[j-1][i];
465:       }

467:       if (i+1 == mx){
468:         xr  = user->right[j-ys+1];
469:         xrb = user->right[j-ys];
470:       } else {
471:         xr  = x[j][i+1];
472:       }

474:       if (j+1==my){
475:         xt  = user->top[i-xs+1];
476:         xlt = user->top[i-xs];
477:       }else {
478:         xt  = x[j+1][i];
479:       }

481:       if (i>0 && j+1<my){
482:         xlt = x[j+1][i-1];
483:       }
484:       if (j>0 && i+1<mx){
485:         xrb = x[j-1][i+1];
486:       }


489:       d1 = (xc-xl)/hx;
490:       d2 = (xc-xr)/hx;
491:       d3 = (xc-xt)/hy;
492:       d4 = (xc-xb)/hy;
493:       d5 = (xrb-xr)/hy;
494:       d6 = (xrb-xb)/hx;
495:       d7 = (xlt-xl)/hy;
496:       d8 = (xlt-xt)/hx;

498:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
499:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
500:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
501:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
502:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
503:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);


506:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
507:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
508:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
509:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
510:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
511:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
512:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
513:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

518:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
519:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
520:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
521:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

525:       row.j = j; row.i = i;
526:       k=0;
527:       if (j>0){
528:         v[k]=hb;
529:         col[k].j = j - 1; col[k].i = i;
530:         k++;
531:       }

533:       if (j>0 && i < mx -1){
534:         v[k]=hbr;
535:         col[k].j = j - 1; col[k].i = i+1;
536:         k++;
537:       }

539:       if (i>0){
540:         v[k]= hl;
541:         col[k].j = j; col[k].i = i-1;
542:         k++;
543:       }

545:       v[k]= hc;
546:       col[k].j = j; col[k].i = i;
547:       k++;

549:       if (i < mx-1 ){
550:         v[k]= hr;
551:         col[k].j = j; col[k].i = i+1;
552:         k++;
553:       }

555:       if (i>0 && j < my-1 ){
556:         v[k]= htl;
557:         col[k].j = j+1; col[k].i = i-1;
558:         k++;
559:       }

561:       if (j < my-1 ){
562:         v[k]= ht;
563:         col[k].j = j+1; col[k].i = i;
564:         k++;
565:       }

567:       /*
568:          Set matrix values using local numbering, which was defined
569:          earlier, in the main routine.
570:       */
571:       MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
572:     }
573:   }

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

578:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
579:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

581:   PetscLogFlops(199*xm*ym);
582:   return(0);
583: }

585: /* ------------------------------------------------------------------- */
586: /*
587:    MSA_BoundaryConditions -  Calculates the boundary conditions for
588:    the region.

590:    Input Parameter:
591: .  user - user-defined application context

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

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

615:   bsize=xm+2;
616:   lsize=ym+2;
617:   rsize=ym+2;
618:   tsize=xm+2;

620:   PetscMalloc1(bsize,&user->bottom);
621:   PetscMalloc1(tsize,&user->top);
622:   PetscMalloc1(lsize,&user->left);
623:   PetscMalloc1(rsize,&user->right);

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

627:   for (j=0; j<4; j++){
628:     if (j==0){
629:       yt=b;
630:       xt=l+hx*xs;
631:       limit=bsize;
632:       boundary=user->bottom;
633:     } else if (j==1){
634:       yt=t;
635:       xt=l+hx*xs;
636:       limit=tsize;
637:       boundary=user->top;
638:     } else if (j==2){
639:       yt=b+hy*ys;
640:       xt=l;
641:       limit=lsize;
642:       boundary=user->left;
643:     } else { /* if (j==3) */
644:       yt=b+hy*ys;
645:       xt=r;
646:       limit=rsize;
647:       boundary=user->right;
648:     }

650:     for (i=0; i<limit; i++){
651:       u1=xt;
652:       u2=-yt;
653:       for (k=0; k<maxits; k++){
654:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
655:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
656:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
657:         if (fnorm <= tol) break;
658:         njac11=one+u2*u2-u1*u1;
659:         njac12=two*u1*u2;
660:         njac21=-two*u1*u2;
661:         njac22=-one - u1*u1 + u2*u2;
662:         det = njac11*njac22-njac21*njac12;
663:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
664:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
665:       }

667:       boundary[i]=u1*u1-u2*u2;
668:       if (j==0 || j==1) {
669:         xt=xt+hx;
670:       } else { /*  if (j==2 || j==3) */
671:         yt=yt+hy;
672:       }

674:     }

676:   }

678:   /* Scale the boundary if desired */
679:   if (1==1){
680:     PetscReal scl = 1.0;

682:     PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
683:     if (flg){
684:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
685:     }

687:     PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
688:     if (flg){
689:       for (i=0;i<tsize;i++) user->top[i]*=scl;
690:     }

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

697:     PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
698:     if (flg){
699:       for (i=0;i<lsize;i++) user->left[i]*=scl;
700:     }
701:   }
702:   return(0);
703: }

705: /* ------------------------------------------------------------------- */
706: /*
707:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

709:    Input Parameters:
710: .  user - user-defined application context
711: .  X - vector for initial guess

713:    Output Parameters:
714: .  X - newly computed initial guess
715: */
716: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
717: {
719:   PetscInt       start2=-1,i,j;
720:   PetscReal      start1=0;
721:   PetscBool      flg1,flg2;

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

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

729:     VecSet(X, start1);

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

734:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
735:     VecSetRandom(X, rctx);
736:     PetscRandomDestroy(&rctx);
737:     VecShift(X, np5);

739:   } else { /* Take an average of the boundary conditions */
740:     PetscInt  xs,xm,ys,ym;
741:     PetscInt  mx=user->mx,my=user->my;
742:     PetscReal **x;

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

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

750:     /* Perform local computations */
751:     for (j=ys; j<ys+ym; j++){
752:       for (i=xs; i< xs+xm; i++){
753:         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;
754:       }
755:     }
756:     DMDAVecRestoreArray(user->dm,X,(void**)&x);
757:     PetscLogFlops(9*xm*ym);
758:   }
759:   return(0);
760: }

762: /*-----------------------------------------------------------------------*/
763: PetscErrorCode My_Monitor(Tao tao, void *ctx)
764: {
766:   Vec            X;

769:   TaoGetSolutionVector(tao,&X);
770:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
771:   return(0);
772: }