Actual source code: minsurf3.c

  1: /*$Id: minsurf3.c,v 1.1 2002/05/13 19:32:57 benson Exp $*/

  3: /* Program usage: mpirun -np <proc> minsurf3 [-help] [all TAO options] */

  5: /* minsurf1: boundary values like in COPS example:
  6:       - parabola for top and bottom boundaries
  7:       - zero o.w.  */

  9: /*
 10:   Include "tao.h" so we can use TAO solvers.
 11:   petscda.h for distributed array
 12:   ad_deriv.h for AD gradient
 13: */

 15: #include "petscda.h"
 16: #include "tao.h"
 17: #include "taodaapplication.h"

 19: static char  help[] = 
 20: "Given a rectangular 2-D domain and boundary values along \n\
 21: the edges of the domain, the objective is to find the surface with \n\
 22: the minimal area that satisfies the boundary conditions.\n\
 23: The command line options are:\n\
 24:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 25:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 26:   -nlevels <nlevels>, where <nlevels> = number of levels in multigrid\n\
 27:   -byelement, if computation is made by functions on rectangular elements\n\
 28:   -adic, if AD is used (AD is not used by default)\n\
 29:   -bottom <b>, where <b> = bottom bound on the rectangular domain\n\
 30:   -top <t>, where <t> = bottom bound on the rectangular domain\n\
 31:   -left <l>, where <l> = left bound on the rectangular domain\n\
 32:   -right <r>, where <r> = right bound on the rectangular domain\n\
 33:   -cops <r>, if COPS boundaries should be use (default MINPACK bounds)\n\
 34:   -obst <obst>, where <obst> = 1 is obstacle is present, zero otherwise \n\n";

 36: /*T
 37:    Concepts: TAO - Solving a bounded minimization problem
 38:    Routines: TaoInitialize(); TaoFinalize();
 39:    Routines: TaoCreate(); TaoDestroy();
 40:    Routines: DAAppSetVariableBoundsRoutine();
 41:    Routines: DAAppSetElementObjectiveAndGradientRoutine();
 42:    Routines: DAAppSetElementHessianRoutine();
 43:    Routines: DAAppSetObjectiveAndGradientRoutine();
 44:    Routines: DAAppSetADElementFunctionGradient();
 45:    Routines: DAAppSetHessianRoutine();
 46:    Routines: TaoSetOptions();
 47:    Routines: TaoAppGetSolutionStatus(); TaoDAAppSolve();
 48:    Routines: DAAppSetBeforeMonitor(); TaoView();
 49:    Routines: DAAppGetSolution();
 50:    Routines: DAAppGetInterpolationMatrix();
 51:    Processors: n
 52: T*/

 54: /*
 55:    User-defined application context - contains data needed by the
 56:    application-provided call-back routines.
 57: */


 60: /*  
 61:     This structure is used only when an ADIC generated gradient is used.
 62:     An InactiveDouble type is a double 
 63: */
 64: typedef struct {
 65:   
 66:   InactiveDouble      hx, hy;        /* increment size in both directions */
 67:   InactiveDouble      area;          /* area of the triangles */

 69: } ADFGCtx;
 70: int ad_MSurfLocalFunction(PetscInt[2], DERIV_TYPE[4], DERIV_TYPE*, void*);

 72: typedef struct {
 73:   PetscReal      b, t, l, r;    /* domain boundaries */
 74:   PetscReal      bheight;       /* height of obstacle    */
 75:   PetscReal      fx, fy;        /* relative size of obstacle */
 76:   double      hx, hy;        /* increment size in both directions */
 77:   double      area;          /* area of the triangles */
 78:   PetscInt         mx, my;        /* discretization including boundaries */
 79:   PetscInt         bmx, bmy;      /* discretization of obstacle */
 80:   ADFGCtx     fgctx;         /* Used only when an ADIC generated gradient is used */
 81: } AppCtx;

 83: /* User-defined routines */
 84: static int AppCtxInitialize(void *);

 86: static int MSurfLocalFunctionGradient(PetscInt[2], double[4], double *, double[4], void *);
 87: static int WholeMSurfFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);

 89: static int MSurfLocalHessian(PetscInt[2], double[4], double[4][4], void *);
 90: static int WholeMSurfHessian(TAO_APPLICATION,DA,Vec,Mat,void*);

 92: static int COPS_Bounds(TAO_APPLICATION, DA, Vec, Vec, void*);
 93: static int MINPACK_Bounds(TAO_APPLICATION, DA, Vec, Vec, void *);

 95: static int MyGridMonitorBefore(TAO_APPLICATION, DA, PetscInt, void *);


100: int main( int argc, char **argv ) {

102:   int                  info;         /* used to check for functions returning nonzeros */
103:   PetscInt                  iter,nlevels;
104:   PetscInt                  Nx,Ny;
105:   double               ff,gnorm;
106:   DA                   DAarray[20];
107:   Vec                  X;
108:   PetscTruth           flg, PreLoad = PETSC_TRUE;                             /* flags */
109:   TaoMethod            method = "tao_bnls";                     /* minimization method */
110:   TaoTerminateReason   reason;
111:   TAO_SOLVER           tao;                               /* TAO_SOLVER solver context */
112:   TAO_APPLICATION   MSurfApp;                              /* The PETSc application */
113:   AppCtx               user;                              /* user-defined work context */
114:   KSP ksp;
115:   PC  pc;

117:   /* Initialize TAO */
118:   PetscInitialize(&argc, &argv, (char *)0, help);
119:   TaoInitialize(&argc, &argv, (char *)0, help);

121:   PreLoadBegin(PreLoad,"Solve");

123:   info = AppCtxInitialize((void*)&user); CHKERRQ(info);

125:   nlevels=4;
126:   info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
127:   if (PreLoadIt == 0) {
128:     nlevels = 1; user.mx = 11; user.my = 11;
129:   }

131:   PetscPrintf(MPI_COMM_WORLD,"\n---- Minimal Surface Problem (simple boundary) -----\n\n");

133:   /* Let PETSc determine the vector distribution */
134:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

136:   /* Create distributed array (DA) to manage parallel grid and vectors  */
137:   info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
138:                     user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
139:   for (iter=1;iter<nlevels;iter++){
140:     info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
141:   }

143:   /* Create TAO solver and set desired solution method */
144:   info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
145:   info = TaoApplicationCreate(PETSC_COMM_WORLD, &MSurfApp); CHKERRQ(info);
146:   info = TaoAppSetDAApp(MSurfApp,DAarray,nlevels); CHKERRQ(info);

148:   /* Sets routines for function, gradient and bounds evaluation */
149:   info = PetscOptionsHasName(TAO_NULL, "-cops", &flg); CHKERRQ(info);
150:   if (flg){
151:     info = DAAppSetVariableBoundsRoutine(MSurfApp,COPS_Bounds,(void *)&user); CHKERRQ(info);
152:   } else {
153:     info = DAAppSetVariableBoundsRoutine(MSurfApp,MINPACK_Bounds,(void *)&user); CHKERRQ(info);
154:   }

156:   info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
157:   if (flg) {

159:     /* Sets routines for function and gradient evaluation, element by element */
160:     info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
161:     if (flg) {
162:       info = DAAppSetADElementFunctionGradient(MSurfApp,ad_MSurfLocalFunction,150,(void *)&user.fgctx); CHKERRQ(info);
163:     } else {
164:       info = DAAppSetElementObjectiveAndGradientRoutine(MSurfApp,MSurfLocalFunctionGradient,36,(void *)&user); CHKERRQ(info);
165:     }

167:     /* Sets routines for Hessian evaluation, element by element */
168:     info = DAAppSetElementHessianRoutine(MSurfApp,MSurfLocalHessian,87,(void*)&user); CHKERRQ(info);

170:   } else {

172:     /* Sets routines for function and gradient evaluation, all in one routine */
173:     info = DAAppSetObjectiveAndGradientRoutine(MSurfApp,WholeMSurfFunctionGradient,(void *)&user); CHKERRQ(info);

175:     /* Sets routines for Hessian evaluation, all in one routine */
176:     info = DAAppSetHessianRoutine(MSurfApp,WholeMSurfHessian,(void*)&user); CHKERRQ(info);    
177:     
178:   }

180:   info = DAAppSetBeforeMonitor(MSurfApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
181:   info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
182:   if (flg){
183:     info = DAAppPrintInterpolationError(MSurfApp); CHKERRQ(info);
184:     info = DAAppPrintStageTimes(MSurfApp); CHKERRQ(info);
185:   }

187:   info = TaoAppSetRelativeTolerance(MSurfApp,1.0e-6); CHKERRQ(info);
188:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
189:   info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);

191:   info = TaoAppGetKSP(MSurfApp,&ksp);CHKERRQ(info);
192:   info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
193:   info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100);CHKERRQ(info);
194:   info = KSPGetPC(ksp,&pc);CHKERRQ(info);
195:   info = PCSetType(pc,PCBJACOBI);CHKERRQ(info);

197:   /* Check for any tao command line options */
198:   info = TaoSetOptions(MSurfApp,tao); CHKERRQ(info);

200:   /* SOLVE THE APPLICATION */
201:   info = TaoDAAppSolve(MSurfApp,tao);  CHKERRQ(info);

203:   /* Get information on termination */
204:   info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
205:   if (reason <= 0 ){
206:     PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
207:     PetscPrintf(MPI_COMM_WORLD," Iterations: %d,  Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
208:   }

210:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
211:   if (flg){
212:     info = DAAppGetSolution(MSurfApp,nlevels-1,&X); CHKERRQ(info);
213:     info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
214:   }

216:   /*  To View TAO solver information */
217:   // info = TaoView(tao); CHKERRQ(info);

219:   /* Free TAO data structures */
220:   info = TaoDestroy(tao); CHKERRQ(info);
221:   info = TaoAppDestroy(MSurfApp); CHKERRQ(info);

223:   /* Free PETSc data structures */
224:   for (iter=0;iter<nlevels;iter++){
225:     info = DADestroy(DAarray[iter]); CHKERRQ(info);
226:   }

228:   PreLoadEnd();

230:   /* Finalize TAO */
231:   TaoFinalize();
232:   PetscFinalize();

234:   return 0;
235: } /* main */



239: /*----- The following two routines
240:   MyGridMonitorBefore    MyGridMonitorAfter
241:   help diplay info of iterations at every grid level 
242: */

246: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx) {

248:   AppCtx *user = (AppCtx*)ctx;
249:   int info;
250:   PetscInt mx,my;

252:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
253:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
254:   user->mx = mx;
255:   user->my = my;
256:   user->hx = (user->r - user->l) / (user->mx - 1);
257:   user->hy = (user->t - user->b) / (user->my - 1);
258:   user->area = 0.5 * user->hx * user->hy;
259:   user->fgctx.hx   = user->hx;
260:   user->fgctx.hy   = user->hy;
261:   user->fgctx.area = user->area;

263:   user->bmx=(PetscInt)((mx+1)*user->fx); user->bmy=(PetscInt)((my+1)*user->fy); 

265:   PetscPrintf(MPI_COMM_WORLD,"Grid: %d,    mx: %d     my: %d   \n",level,mx,my);

267:   return 0;
268: }

272: /*
273:   AppCtxInitialize - Sets initial values for the application context parameters

275:   Input:
276:     ptr - void user-defined application context

278:   Output:
279:     ptr - user-defined application context with the default or user-provided
280:              parameters
281: */
282: static int AppCtxInitialize(void *ptr) {

284:   AppCtx *user = (AppCtx*)ptr;
285:   PetscTruth flg;
286:   int info;

288:   /* Specify default parameters */
289:   user->mx = user->my = 11;
290:   user->b = -0.5;
291:   user->t = 0.5;
292:   user->l = -0.5;
293:   user->r = 0.5;
294:   user->fx=0.5;
295:   user->fy=0.5;
296:   user->bheight=0.0;

298:   /* Check for command line arguments that override defaults */
299:   info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
300:   info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
301:   info = PetscOptionsGetReal(TAO_NULL, "-bottom", &user->b, &flg); CHKERRQ(info);
302:   info = PetscOptionsGetReal(TAO_NULL, "-top", &user->t, &flg); CHKERRQ(info);
303:   info = PetscOptionsGetReal(TAO_NULL, "-left", &user->l, &flg); CHKERRQ(info);
304:   info = PetscOptionsGetReal(TAO_NULL, "-right", &user->r, &flg); CHKERRQ(info);
305:   info = PetscOptionsGetReal(PETSC_NULL,"-bmx",&user->fx,&flg); CHKERRQ(info);
306:   info = PetscOptionsGetReal(PETSC_NULL,"-bmy",&user->fy,&flg); CHKERRQ(info);
307:   info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user->bheight,&flg); CHKERRQ(info);

309:   user->hx = (user->r - user->l) / (user->mx - 1);
310:   user->hy = (user->t - user->b) / (user->my - 1);
311:   user->area = 0.5 * user->hx * user->hy;
312:   info = PetscLogFlops(8); CHKERRQ(info);

314:   return 0;

316: } /* AppCtxInitialize */


319: /*------- USER-DEFINED: set the upper and lower bounds for the variables  -------*/

323: /*
324:   FormBounds - Forms bounds on the variables

326:   Input:
327:     user - user-defined application context

329:   Output:
330:     XL - vector of lower bounds
331:     XU - vector of upper bounds

333: */
334: static int COPS_Bounds(TAO_APPLICATION tao, DA da, Vec XL, Vec XU, void *ptr) {

336:   AppCtx *user = (AppCtx*)ptr;
337:   PetscTruth flg;
338:   int info;
339:   PetscInt i, j, mx, my, xs, xm, ys, ym;
340:   double lb = -TAO_INFINITY;
341:   double ub = TAO_INFINITY;
342:   double **xl, **xu;
343:   double xi, xi1, xi2;
344:   double cx, cy, radx, rady, height = 1.0;

346:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
347:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
348:   user->hx = (user->r - user->l) / (mx - 1);
349:   user->hy = (user->t - user->b) / (my - 1);
350:   user->area = 0.5 * user->hx * user->hy;

352:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
353:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
354:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

356:   /* Compute default bounds */
357:   for (j = ys; j < ys+ym; j++){
358:     for (i = xs; i < xs+xm; i++){

360:       if (j == 0 || j == my - 1) {
361:         xi = user->l + i * user->hx;
362:         xl[j][i] = xu[j][i] =  -4 * (xi - user->l) * (xi - user->r);
363:       } else if (i == 0 || i == mx - 1) {
364:         xl[j][i] = xu[j][i] = 0.0;
365:       } else {
366:         xl[j][i] = lb;
367:         xu[j][i] = ub;
368:       }

370:     }
371:   }

373:   /* Adjust lower bound if obstacle is present */
374:   info = PetscOptionsHasName(PETSC_NULL, "-obst", &flg); CHKERRQ(info);
375:   if (flg) {
376:     radx = (user->r - user->l) * 0.25;
377:     cx = user->l + 2.0 * radx;
378:     rady = (user->t - user->b) * 0.25;
379:     cy = user->b + 2.0 * rady;
380:     for (j = ys; j < ys+ym; j++){
381:       for (i = xs; i < xs+xm; i++){
382:         xi1 = user->l + i * user->hx;
383:         xi2 = user->b + j * user->hy;
384:         if ( fabs(xi1 - cx) <= radx && fabs(xi2 - cy) <= rady ) {
385:           xl[j][i] = height;
386:         }
387:       }
388:     }
389:     info = PetscLogFlops(8 + xm * ym * 6); CHKERRQ(info);
390:   }

392:   info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
393:   info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);

395:   info = PetscLogFlops(12 * ym + 6); CHKERRQ(info);

397:   return 0;

399: } /* DAGetBounds2d */


404: static int MINPACK_Bounds(TAO_APPLICATION tao, DA da, Vec XL,Vec XU, void *user1){

406:   AppCtx *user=(AppCtx*)user1;
407:   int info;
408:   PetscInt i,j,k,limit=0,maxits=5;
409:   PetscInt xs,ys,xm,ym;
410:   PetscInt mx, my, bmy, bmx;
411:   PetscInt row=0, bsize=0, lsize=0, tsize=0, rsize=0;
412:   double bheight=user->bheight;
413:   double one=1.0, two=2.0, three=3.0, tol=1e-10;
414:   double fnorm,det,hx,hy,xt=0,yt=0;
415:   double u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
416:   double b=-0.5, t=0.5, l=-0.5, r=0.5;
417:   PetscScalar scl,*xl,*xu, lb=TAO_NINFINITY, ub=TAO_INFINITY;
418:   PetscTruth flg;
419:   double **xll;

421:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
422:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);

424:   bheight=user->bheight,lb=-1000; ub=1000;
425:   bmx=user->bmx; bmy=user->bmy;

427:   info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);

429:   info = VecSet(XL, lb); CHKERRQ(info);
430:   info = VecSet(XU, ub); CHKERRQ(info);

432:   /*
433:     Get pointers to vector data
434:   */
435:   info = DAVecGetArray(da,XL,(void**)&xll); CHKERRQ(info);

437:   if (ys==0) bsize=xm;
438:   if (xs==0) lsize=ym;
439:   if (xs+xm==mx) rsize=ym;
440:   if (ys+ym==my) tsize=xm;
441:   
442:   hx= (r-l)/(mx-1); hy=(t-b)/(my-1);

444:   /* Compute the optional lower box */
445:   for (i=xs; i< xs+xm; i++){    
446:     for (j=ys; j<ys+ym; j++){
447:       
448:       if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 && j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
449:         xll[j][i] = bheight;
450:       }

452:     }
453:   }

455:   info = DAVecRestoreArray(da,XL,(void**)&xll); CHKERRQ(info);

457:   /* Boundary Values */
458:   info = VecGetArray(XL,&xl); CHKERRQ(info);
459:   info = VecGetArray(XU,&xu); CHKERRQ(info);

461:   for (j=0; j<4; j++){
462:     if (j==0){
463:       yt=b;
464:       xt=l+hx*xs;
465:       limit=bsize;
466:       scl=1.0;
467:       info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg); 
468:     } else if (j==1){
469:       yt=t;
470:       xt=l+hx*xs;
471:       limit=tsize;
472:       scl=1.0;
473:       info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg); 
474:     } else if (j==2){
475:       yt=b+hy*ys;
476:       xt=l;
477:       limit=lsize;
478:       scl=1.0;
479:       info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg); 
480:     } else if (j==3){
481:       yt=b+hy*ys;
482:       xt=r;
483:       limit=rsize;
484:       scl=1.0;
485:       info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg); 
486:     }

488:     for (i=0; i<limit; i++){
489:       u1=xt;
490:       u2=-yt;
491:       for (k=0; k<maxits; k++){
492:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
493:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
494:         fnorm=sqrt(nf1*nf1+nf2*nf2);
495:         if (fnorm <= tol) break;
496:         njac11=one+u2*u2-u1*u1;
497:         njac12=two*u1*u2;
498:         njac21=-two*u1*u2;
499:         njac22=-one - u1*u1 + u2*u2;
500:         det = njac11*njac22-njac21*njac12;
501:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
502:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
503:       }

505:       if (j==0){
506:         row = i;
507:       } else if (j==1){
508:         row = (ym-1)*xm+i;
509:       } else if (j==2){
510:         row = (xm*i);
511:       } else if (j==3){
512:         row = xm*(i+1)-1;
513:       }
514:       
515:       xl[row]=(u1*u1-u2*u2)*scl;
516:       xu[row]=(u1*u1-u2*u2)*scl;

518:       if (j==0 || j==1) {
519:         xt=xt+hx;
520:       } else if (j==2 || j==3){
521:         yt=yt+hy;
522:       }
523:       
524:     }
525:     
526:   }

528:   info = VecRestoreArray(XL,&xl); CHKERRQ(info);
529:   info = VecRestoreArray(XU,&xu); CHKERRQ(info);

531:   return 0;
532: }

534: /*------- USER-DEFINED: routine to evaluate the function and gradient
535:            at a local (rectangular element) level              -------*/

539: /*
540:   MSurfLocalFunctionGradient - Evaluates function and gradient over the 
541:       local rectangular element

543:   Input:
544:     coor - vector with the indices of the position of current element
545:              in the first, second and third directions
546:     x - current point (values over the current rectangular element)
547:     ptr - user-defined application context

549:   Output:
550:     f - value of the objective funtion at the local rectangular element
551:     g - gradient of the local function

553: */
554: static int MSurfLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {

556:   AppCtx *user = (AppCtx*)ptr;

558:   double hx, hy, area;
559:   double dvdx, dvdy, flow, fup;
560:   double d1,d2;

562:   hx = user->hx;
563:   hy = user->hy;
564:   area = user->area;

566:   /* lower triangle contribution */
567:   dvdx = (x[0] - x[1]);
568:   dvdy = (x[0] - x[2]);
569:   g[1] = (dvdx * (hy/hx))/(-2);
570:   g[2] = (dvdy * (hx/hy))/(-2);
571:   dvdx /= hx;
572:   dvdy /= hy;
573:   flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
574:   *f=flow;
575:   g[1] /= flow;
576:   g[2] /= flow;
577:   g[0] = -(g[1]+g[2]);

579:   /* upper triangle contribution */
580:   dvdx = (x[3] - x[2]);
581:   dvdy = (x[3] - x[1]);
582:   d1 = (dvdy*(hx/hy))/(-2);
583:   d2 = (dvdx*(hy/hx))/(-2);
584:   dvdx /= hx;
585:   dvdy /= hy;
586:   fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
587:   *f += fup;
588:   g[1] += d1/fup;
589:   g[2] += d2/fup; 
590:   g[3] = -(d1+d2)/fup;

592:   *f *= area;

594:   return 0;
595: } /* MSurfLocalFunctionGradient */



601: /*
602:   MSurfLocalHessian - Computes the Hessian of the local (partial) function
603:          defined over the current rectangle

605:   Input:
606:     coor - vector with the indices of the position of current element
607:              in the first, second and third directions
608:     x - current local solution (over the rectangle only)
609:     ptr - user-defined application context

611:   Output:
612:     H - Hessian matrix of the local function (wrt the four
613:            points of the rectangle only)

615: */
616: static int MSurfLocalHessian(PetscInt coor[2], double x[4], double H[4][4],void *ptr) {

618:   AppCtx *user = (AppCtx*)ptr;
619:   double hx, hy, area, byhxhx, byhyhy;
620:   double dvdx, dvdy, flow, fup;
621:   double areadivf, areadivf3;

623:   hx = user->hx;
624:   hy = user->hy;
625:   area = user->area;
626:   
627:   byhxhx = 1.0 / (hx * hx);
628:   byhyhy = 1.0 / (hy * hy);

630:   /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
631:   dvdx = (x[0] - x[1]) / hx;  /* lower triangle contribution */
632:   dvdy = (x[0] - x[2]) / hy;
633:   flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
634:   dvdx = dvdx / hx;
635:   dvdy = dvdy / hy;
636:   areadivf = area / flow;
637:   areadivf3 = areadivf / (flow * flow);
638:   H[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
639:   H[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
640:   H[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
641:   H[0][3] = 0.0;
642:   H[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
643:   H[1][2] = areadivf3 * (-dvdx) * dvdy;
644:   H[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;

646:   /* upper triangle contribution */
647:   dvdx = (x[3] - x[2]) / hx;
648:   dvdy = (x[3] - x[1]) / hy;
649:   fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
650:   dvdx = dvdx / hx;
651:   dvdy = dvdy / hy;
652:   areadivf = area / fup;
653:   areadivf3 = areadivf / (fup * fup);
654:   H[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
655:   H[1][2] += areadivf3 * (-dvdy) * dvdx;
656:   H[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
657:   H[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
658:   H[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
659:   H[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);

661:   H[1][0] = H[0][1];
662:   H[2][0] = H[0][2];
663:   H[3][0] = H[0][3];
664:   H[2][1] = H[1][2];
665:   H[3][1] = H[1][3];
666:   H[3][2] = H[2][3];

668:   return 0;

670: } /* MSurfLocalHessian */


673: /*------- USER-DEFINED: routine to evaluate the function 
674:           and gradient at the whole grid             -------*/

678: /*
679:   WholeMSurfFunctionGradient - Evaluates function and gradient over the 
680:       whole grid

682:   Input:
683:     daapplication - TAO application object
684:     da  - distributed array
685:     X   - the current point, at which the function and gradient are evaluated
686:     ptr - user-defined application context

688:   Output:
689:     f - value of the objective funtion at X
690:     G - gradient at X
691: */
692: static int WholeMSurfFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

694:   AppCtx *user = (AppCtx*)ptr;
695:   Vec localX, localG;
696:   PetscInt i, j;
697:   int info;
698:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
699:   double **x, **g;
700:   double floc = 0.0;
701:   PetscScalar zero = 0.0;

703:   double hx, hy, area;
704:   double dvdx, dvdy, flow, fup;
705:   double areadivf;

707:   hx = user->hx;
708:   hy = user->hy;
709:   area = user->area;

711:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
712:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
713:   info = VecSet(G, zero); CHKERRQ(info);
714:   info = VecSet(localG, zero); CHKERRQ(info);

716:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
717:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

719:   info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
720:   info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);

722:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
723:   info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);

725:   xe = gxs + gxm - 1;
726:   ye = gys + gym - 1;
727:   for (j = ys; j < ye; j++) {
728:     for (i = xs; i < xe; i++) {

730:       /* lower triangle contribution */
731:       dvdx = (x[j][i] - x[j][i+1]) / hx;  
732:       dvdy = (x[j][i] - x[j+1][i]) / hy;
733:       flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
734:       areadivf = area / flow;
735:       g[j][i] += (dvdx / hx + dvdy / hy) * areadivf;
736:       g[j][i+1] += (-dvdx / hx) * areadivf;
737:       g[j+1][i] += (-dvdy / hy) * areadivf;

739:       /* upper triangle contribution */
740:       dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
741:       dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
742:       fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
743:       areadivf = area / fup;
744:       g[j][i+1] += (-dvdy / hy) * areadivf;
745:       g[j+1][i] += (-dvdx / hx) * areadivf;
746:       g[j+1][i+1] += (dvdx / hx + dvdy / hy) * areadivf;

748:       floc += area * (flow + fup);

750:     }
751:   }

753:   info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);

755:   info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
756:   info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);

758:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
759:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

761:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
762:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

764:   info = PetscLogFlops((xe-xs) * (ye-ys) * 42); CHKERRQ(info);

766:   return 0;
767: } /* WholeMSurfFunctionGradient  */


770: /*------- USER-DEFINED: routine to evaluate the Hessian 
771:           at the whole grid             -------*/

775: /*
776:   WholeMSurfHessian - Evaluates Hessian over the whole grid

778:   Input:
779:     daapplication - TAO application object
780:     da  - distributed array
781:     X   - the current point, at which the function and gradient are evaluated
782:     ptr - user-defined application context

784:   Output:
785:     H - Hessian at X
786: */
787: static int WholeMSurfHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

789:   AppCtx *user = (AppCtx*)ptr;
790:   Vec localX;
791:   int info;
792:   PetscInt  i, j, ind[4];
793:   PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
794:   double smallH[4][4];
795:   double **x;

797:   double hx, hy, area, byhxhx, byhyhy;
798:   double dvdx, dvdy, flow, fup;
799:   double areadivf, areadivf3;
800:   PetscTruth assembled;

802:   hx = user->hx;
803:   hy = user->hy;
804:   area = user->area;
805:   
806:   byhxhx = 1.0 / (hx * hx);
807:   byhyhy = 1.0 / (hy * hy);

809:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
810:   info = MatAssembled(H,&assembled); CHKERRQ(info);
811:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}

813:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
814:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

816:   info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);

818:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
819:   info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);

821:   xe = gxs + gxm - 1;
822:   ye = gys + gym - 1;
823:   for (j = ys; j < ye; j++) {
824:     for (i = xs; i < xe; i++) {

826:       /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
827:       dvdx = (x[j][i] - x[j][i+1]) / hx;  /* lower triangle contribution */
828:       dvdy = (x[j][i] - x[j+1][i]) / hy;
829:       flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
830:       dvdx = dvdx / hx;
831:       dvdy = dvdy / hy;
832:       areadivf = area / flow;
833:       areadivf3 = areadivf / (flow * flow);
834:       smallH[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
835:       smallH[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
836:       smallH[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
837:       smallH[0][3] = 0.0;
838:       smallH[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
839:       smallH[1][2] = areadivf3 * (-dvdx) * dvdy;
840:       smallH[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;

842:       /* upper triangle contribution */
843:       dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
844:       dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
845:       fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
846:       dvdx = dvdx / hx;
847:       dvdy = dvdy / hy;
848:       areadivf = area / fup;
849:       areadivf3 = areadivf / (fup * fup);
850:       smallH[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
851:       smallH[1][2] += areadivf3 * (-dvdy) * dvdx;
852:       smallH[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
853:       smallH[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
854:       smallH[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
855:       smallH[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);

857:       smallH[1][0] = smallH[0][1];
858:       smallH[2][0] = smallH[0][2];
859:       smallH[3][0] = smallH[0][3];
860:       smallH[2][1] = smallH[1][2];
861:       smallH[3][1] = smallH[1][3];
862:       smallH[3][2] = smallH[2][3];

864:       ind[0] = (j-gys) * gxm + (i-gxs);
865:       ind[1] = ind[0] + 1;
866:       ind[2] = ind[0] + gxm;
867:       ind[3] = ind[2] + 1;
868:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

870:     }
871:   }

873:   info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);

875:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
876:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
877:   info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);

879:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);

881:   info = PetscLogFlops((xe-xs) * (ye-ys) * 83 + 4); CHKERRQ(info);
882:   return 0;

884: } /* WholeMSurfHessian */