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(int[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:   int         mx, my;        /* discretization including boundaries */
 79:   int         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(int[2], double[4], double *, double[4], void *);
 87: static int WholeMSurfFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);

 89: static int MSurfLocalHessian(int[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, int, void *);


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

102:   int                  info;         /* used to check for functions returning nonzeros */
103:   int                  iter,nlevels;
104:   int                  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 */

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

119:   PreLoadBegin(PreLoad,"Solve");

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

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

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

131:   /* Let PETSc determine the vector distribution */
132:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

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

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

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

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

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

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

168:   } else {

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

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

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

185:   info = TaoAppSetRelativeTolerance(MSurfApp,1.0e-8); CHKERRQ(info);
186:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
187:   info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);

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

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

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

202:   info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
203:   if (flg){
204:     info = DAAppGetSolution(MSurfApp,nlevels-1,&X); CHKERRQ(info);
205:     info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
206:   }

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

211:   /* Free TAO data structures */
212:   info = TaoDestroy(tao); CHKERRQ(info);
213:   info = TaoAppDestroy(MSurfApp); CHKERRQ(info);

215:   /* Free PETSc data structures */
216:   for (iter=0;iter<nlevels;iter++){
217:     info = DADestroy(DAarray[iter]); CHKERRQ(info);
218:   }

220:   PreLoadEnd();

222:   /* Finalize TAO */
223:   TaoFinalize();
224:   PetscFinalize();

226:   return 0;
227: } /* main */



231: /*----- The following two routines
232:   MyGridMonitorBefore    MyGridMonitorAfter
233:   help diplay info of iterations at every grid level 
234: */

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

240:   AppCtx *user = (AppCtx*)ctx;
241:   int info,mx,my;

243:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
244:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
245:   user->mx = mx;
246:   user->my = my;
247:   user->hx = (user->r - user->l) / (user->mx - 1);
248:   user->hy = (user->t - user->b) / (user->my - 1);
249:   user->area = 0.5 * user->hx * user->hy;
250:   user->fgctx.hx   = user->hx;
251:   user->fgctx.hy   = user->hy;
252:   user->fgctx.area = user->area;

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

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

258:   return 0;
259: }

263: /*
264:   AppCtxInitialize - Sets initial values for the application context parameters

266:   Input:
267:     ptr - void user-defined application context

269:   Output:
270:     ptr - user-defined application context with the default or user-provided
271:              parameters
272: */
273: static int AppCtxInitialize(void *ptr) {

275:   AppCtx *user = (AppCtx*)ptr;
276:   PetscTruth flg;
277:   int info;

279:   /* Specify default parameters */
280:   user->mx = user->my = 11;
281:   user->b = -0.5;
282:   user->t = 0.5;
283:   user->l = -0.5;
284:   user->r = 0.5;
285:   user->fx=0.5;
286:   user->fy=0.5;
287:   user->bheight=0.0;

289:   /* Check for command line arguments that override defaults */
290:   info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
291:   info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
292:   info = PetscOptionsGetReal(TAO_NULL, "-bottom", &user->b, &flg); CHKERRQ(info);
293:   info = PetscOptionsGetReal(TAO_NULL, "-top", &user->t, &flg); CHKERRQ(info);
294:   info = PetscOptionsGetReal(TAO_NULL, "-left", &user->l, &flg); CHKERRQ(info);
295:   info = PetscOptionsGetReal(TAO_NULL, "-right", &user->r, &flg); CHKERRQ(info);
296:   info = PetscOptionsGetReal(PETSC_NULL,"-bmx",&user->fx,&flg); CHKERRQ(info);
297:   info = PetscOptionsGetReal(PETSC_NULL,"-bmy",&user->fy,&flg); CHKERRQ(info);
298:   info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user->bheight,&flg); CHKERRQ(info);

300:   user->hx = (user->r - user->l) / (user->mx - 1);
301:   user->hy = (user->t - user->b) / (user->my - 1);
302:   user->area = 0.5 * user->hx * user->hy;
303:   info = PetscLogFlops(8); CHKERRQ(info);

305:   return 0;

307: } /* AppCtxInitialize */


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

314: /*
315:   FormBounds - Forms bounds on the variables

317:   Input:
318:     user - user-defined application context

320:   Output:
321:     XL - vector of lower bounds
322:     XU - vector of upper bounds

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

327:   AppCtx *user = (AppCtx*)ptr;
328:   PetscTruth flg;
329:   int i, j, mx, my, info, xs, xm, ys, ym;
330:   double lb = -TAO_INFINITY;
331:   double ub = TAO_INFINITY;
332:   double **xl, **xu;
333:   double xi, xi1, xi2;
334:   double cx, cy, radx, rady, height = 1.0;

336:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
337:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
338:   user->hx = (user->r - user->l) / (mx - 1);
339:   user->hy = (user->t - user->b) / (my - 1);
340:   user->area = 0.5 * user->hx * user->hy;

342:   info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
343:   info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
344:   info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);

346:   /* Compute default bounds */
347:   for (j = ys; j < ys+ym; j++){
348:     for (i = xs; i < xs+xm; i++){

350:       if (j == 0 || j == my - 1) {
351:         xi = user->l + i * user->hx;
352:         xl[j][i] = xu[j][i] =  -4 * (xi - user->l) * (xi - user->r);
353:       } else if (i == 0 || i == mx - 1) {
354:         xl[j][i] = xu[j][i] = 0.0;
355:       } else {
356:         xl[j][i] = lb;
357:         xu[j][i] = ub;
358:       }

360:     }
361:   }

363:   /* Adjust lower bound if obstacle is present */
364:   info = PetscOptionsHasName(PETSC_NULL, "-obst", &flg); CHKERRQ(info);
365:   if (flg) {
366:     radx = (user->r - user->l) * 0.25;
367:     cx = user->l + 2.0 * radx;
368:     rady = (user->t - user->b) * 0.25;
369:     cy = user->b + 2.0 * rady;
370:     for (j = ys; j < ys+ym; j++){
371:       for (i = xs; i < xs+xm; i++){
372:         xi1 = user->l + i * user->hx;
373:         xi2 = user->b + j * user->hy;
374:         if ( fabs(xi1 - cx) <= radx && fabs(xi2 - cy) <= rady ) {
375:           xl[j][i] = height;
376:         }
377:       }
378:     }
379:     info = PetscLogFlops(8 + xm * ym * 6); CHKERRQ(info);
380:   }

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

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

387:   return 0;

389: } /* DAGetBounds2d */


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

396:   AppCtx *user=(AppCtx*)user1;
397:   int i,j,k,limit=0,info,maxits=5;
398:   int xs,ys,xm,ym;
399:   int mx, my, bmy, bmx;
400:   int row=0, bsize=0, lsize=0, tsize=0, rsize=0;
401:   double bheight=user->bheight;
402:   double one=1.0, two=2.0, three=3.0, tol=1e-10;
403:   double fnorm,det,hx,hy,xt=0,yt=0;
404:   double u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
405:   double b=-0.5, t=0.5, l=-0.5, r=0.5;
406:   PetscScalar scl,*xl,*xu, lb=TAO_NINFINITY, ub=TAO_INFINITY;
407:   PetscTruth flg;
408:   double **xll;

410:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
411:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);

413:   bheight=user->bheight,lb=-1000; ub=1000;
414:   bmx=user->bmx; bmy=user->bmy;

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

418:   info = VecSet(XL, lb); CHKERRQ(info);
419:   info = VecSet(XU, ub); CHKERRQ(info);

421:   /*
422:     Get pointers to vector data
423:   */
424:   info = DAVecGetArray(da,XL,(void**)&xll); CHKERRQ(info);

426:   if (ys==0) bsize=xm;
427:   if (xs==0) lsize=ym;
428:   if (xs+xm==mx) rsize=ym;
429:   if (ys+ym==my) tsize=xm;
430: 
431:   hx= (r-l)/(mx-1); hy=(t-b)/(my-1);

433:   /* Compute the optional lower box */
434:   for (i=xs; i< xs+xm; i++){
435:     for (j=ys; j<ys+ym; j++){
436: 
437:       if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 && j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
438:         xll[j][i] = bheight;
439:       }

441:     }
442:   }

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

446:   /* Boundary Values */
447:   info = VecGetArray(XL,&xl); CHKERRQ(info);
448:   info = VecGetArray(XU,&xu); CHKERRQ(info);

450:   for (j=0; j<4; j++){
451:     if (j==0){
452:       yt=b;
453:       xt=l+hx*xs;
454:       limit=bsize;
455:       scl=1.0;
456:       info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
457:     } else if (j==1){
458:       yt=t;
459:       xt=l+hx*xs;
460:       limit=tsize;
461:       scl=1.0;
462:       info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
463:     } else if (j==2){
464:       yt=b+hy*ys;
465:       xt=l;
466:       limit=lsize;
467:       scl=1.0;
468:       info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
469:     } else if (j==3){
470:       yt=b+hy*ys;
471:       xt=r;
472:       limit=rsize;
473:       scl=1.0;
474:       info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
475:     }

477:     for (i=0; i<limit; i++){
478:       u1=xt;
479:       u2=-yt;
480:       for (k=0; k<maxits; k++){
481:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
482:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
483:         fnorm=sqrt(nf1*nf1+nf2*nf2);
484:         if (fnorm <= tol) break;
485:         njac11=one+u2*u2-u1*u1;
486:         njac12=two*u1*u2;
487:         njac21=-two*u1*u2;
488:         njac22=-one - u1*u1 + u2*u2;
489:         det = njac11*njac22-njac21*njac12;
490:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
491:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
492:       }

494:       if (j==0){
495:         row = i;
496:       } else if (j==1){
497:         row = (ym-1)*xm+i;
498:       } else if (j==2){
499:         row = (xm*i);
500:       } else if (j==3){
501:         row = xm*(i+1)-1;
502:       }
503: 
504:       xl[row]=(u1*u1-u2*u2)*scl;
505:       xu[row]=(u1*u1-u2*u2)*scl;

507:       if (j==0 || j==1) {
508:         xt=xt+hx;
509:       } else if (j==2 || j==3){
510:         yt=yt+hy;
511:       }
512: 
513:     }
514: 
515:   }

517:   info = VecRestoreArray(XL,&xl); CHKERRQ(info);
518:   info = VecRestoreArray(XU,&xu); CHKERRQ(info);

520:   return 0;
521: }

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

528: /*
529:   MSurfLocalFunctionGradient - Evaluates function and gradient over the 
530:       local rectangular element

532:   Input:
533:     coor - vector with the indices of the position of current element
534:              in the first, second and third directions
535:     x - current point (values over the current rectangular element)
536:     ptr - user-defined application context

538:   Output:
539:     f - value of the objective funtion at the local rectangular element
540:     g - gradient of the local function

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

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

547:   double hx, hy, area;
548:   double dvdx, dvdy, flow, fup;
549:   double d1,d2;

551:   hx = user->hx;
552:   hy = user->hy;
553:   area = user->area;

555:   /* lower triangle contribution */
556:   dvdx = (x[0] - x[1]);
557:   dvdy = (x[0] - x[2]);
558:   g[1] = (dvdx * (hy/hx))/(-2);
559:   g[2] = (dvdy * (hx/hy))/(-2);
560:   dvdx /= hx;
561:   dvdy /= hy;
562:   flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
563:   *f=flow;
564:   g[1] /= flow;
565:   g[2] /= flow;
566:   g[0] = -(g[1]+g[2]);

568:   /* upper triangle contribution */
569:   dvdx = (x[3] - x[2]);
570:   dvdy = (x[3] - x[1]);
571:   d1 = (dvdy*(hx/hy))/(-2);
572:   d2 = (dvdx*(hy/hx))/(-2);
573:   dvdx /= hx;
574:   dvdy /= hy;
575:   fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
576:   *f += fup;
577:   g[1] += d1/fup;
578:   g[2] += d2/fup;
579:   g[3] = -(d1+d2)/fup;

581:   *f *= area;

583:   return 0;
584: } /* MSurfLocalFunctionGradient */



590: /*
591:   MSurfLocalHessian - Computes the Hessian of the local (partial) function
592:          defined over the current rectangle

594:   Input:
595:     coor - vector with the indices of the position of current element
596:              in the first, second and third directions
597:     x - current local solution (over the rectangle only)
598:     ptr - user-defined application context

600:   Output:
601:     H - Hessian matrix of the local function (wrt the four
602:            points of the rectangle only)

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

607:   AppCtx *user = (AppCtx*)ptr;
608:   double hx, hy, area, byhxhx, byhyhy;
609:   double dvdx, dvdy, flow, fup;
610:   double areadivf, areadivf3;

612:   hx = user->hx;
613:   hy = user->hy;
614:   area = user->area;
615: 
616:   byhxhx = 1.0 / (hx * hx);
617:   byhyhy = 1.0 / (hy * hy);

619:   /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
620:   dvdx = (x[0] - x[1]) / hx;  /* lower triangle contribution */
621:   dvdy = (x[0] - x[2]) / hy;
622:   flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
623:   dvdx = dvdx / hx;
624:   dvdy = dvdy / hy;
625:   areadivf = area / flow;
626:   areadivf3 = areadivf / (flow * flow);
627:   H[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
628:   H[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
629:   H[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
630:   H[0][3] = 0.0;
631:   H[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
632:   H[1][2] = areadivf3 * (-dvdx) * dvdy;
633:   H[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;

635:   /* upper triangle contribution */
636:   dvdx = (x[3] - x[2]) / hx;
637:   dvdy = (x[3] - x[1]) / hy;
638:   fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
639:   dvdx = dvdx / hx;
640:   dvdy = dvdy / hy;
641:   areadivf = area / fup;
642:   areadivf3 = areadivf / (fup * fup);
643:   H[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
644:   H[1][2] += areadivf3 * (-dvdy) * dvdx;
645:   H[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
646:   H[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
647:   H[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
648:   H[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);

650:   H[1][0] = H[0][1];
651:   H[2][0] = H[0][2];
652:   H[3][0] = H[0][3];
653:   H[2][1] = H[1][2];
654:   H[3][1] = H[1][3];
655:   H[3][2] = H[2][3];

657:   return 0;

659: } /* MSurfLocalHessian */


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

667: /*
668:   WholeMSurfFunctionGradient - Evaluates function and gradient over the 
669:       whole grid

671:   Input:
672:     daapplication - TAO application object
673:     da  - distributed array
674:     X   - the current point, at which the function and gradient are evaluated
675:     ptr - user-defined application context

677:   Output:
678:     f - value of the objective funtion at X
679:     G - gradient at X
680: */
681: static int WholeMSurfFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {

683:   AppCtx *user = (AppCtx*)ptr;
684:   Vec localX, localG;
685:   int info, i, j;
686:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
687:   double **x, **g;
688:   double floc = 0.0;
689:   PetscScalar zero = 0.0;

691:   double hx, hy, area;
692:   double dvdx, dvdy, flow, fup;
693:   double areadivf;

695:   hx = user->hx;
696:   hy = user->hy;
697:   area = user->area;

699:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
700:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
701:   info = VecSet(G, zero); CHKERRQ(info);
702:   info = VecSet(localG, zero); CHKERRQ(info);

704:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
705:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

713:   xe = gxs + gxm - 1;
714:   ye = gys + gym - 1;
715:   for (j = ys; j < ye; j++) {
716:     for (i = xs; i < xe; i++) {

718:       /* lower triangle contribution */
719:       dvdx = (x[j][i] - x[j][i+1]) / hx;
720:       dvdy = (x[j][i] - x[j+1][i]) / hy;
721:       flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
722:       areadivf = area / flow;
723:       g[j][i] += (dvdx / hx + dvdy / hy) * areadivf;
724:       g[j][i+1] += (-dvdx / hx) * areadivf;
725:       g[j+1][i] += (-dvdy / hy) * areadivf;

727:       /* upper triangle contribution */
728:       dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
729:       dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
730:       fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
731:       areadivf = area / fup;
732:       g[j][i+1] += (-dvdy / hy) * areadivf;
733:       g[j+1][i] += (-dvdx / hx) * areadivf;
734:       g[j+1][i+1] += (dvdx / hx + dvdy / hy) * areadivf;

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

738:     }
739:   }

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

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

746:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
747:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

749:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
750:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);

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

754:   return 0;
755: } /* WholeMSurfFunctionGradient  */


758: /*------- USER-DEFINED: routine to evaluate the Hessian 
759:           at the whole grid             -------*/

763: /*
764:   WholeMSurfHessian - Evaluates Hessian over the whole grid

766:   Input:
767:     daapplication - TAO application object
768:     da  - distributed array
769:     X   - the current point, at which the function and gradient are evaluated
770:     ptr - user-defined application context

772:   Output:
773:     H - Hessian at X
774: */
775: static int WholeMSurfHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {

777:   AppCtx *user = (AppCtx*)ptr;
778:   Vec localX;
779:   int info, i, j, ind[4];
780:   int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
781:   double smallH[4][4];
782:   double **x;

784:   double hx, hy, area, byhxhx, byhyhy;
785:   double dvdx, dvdy, flow, fup;
786:   double areadivf, areadivf3;
787:   PetscTruth assembled;

789:   hx = user->hx;
790:   hy = user->hy;
791:   area = user->area;
792: 
793:   byhxhx = 1.0 / (hx * hx);
794:   byhyhy = 1.0 / (hy * hy);

796:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
797:   info = MatAssembled(H,&assembled); CHKERRQ(info);
798:   if (assembled){info = MatZeroEntries(H);  CHKERRQ(info);}

800:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
801:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

808:   xe = gxs + gxm - 1;
809:   ye = gys + gym - 1;
810:   for (j = ys; j < ye; j++) {
811:     for (i = xs; i < xe; i++) {

813:       /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
814:       dvdx = (x[j][i] - x[j][i+1]) / hx;  /* lower triangle contribution */
815:       dvdy = (x[j][i] - x[j+1][i]) / hy;
816:       flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
817:       dvdx = dvdx / hx;
818:       dvdy = dvdy / hy;
819:       areadivf = area / flow;
820:       areadivf3 = areadivf / (flow * flow);
821:       smallH[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
822:       smallH[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
823:       smallH[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
824:       smallH[0][3] = 0.0;
825:       smallH[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
826:       smallH[1][2] = areadivf3 * (-dvdx) * dvdy;
827:       smallH[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;

829:       /* upper triangle contribution */
830:       dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
831:       dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
832:       fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
833:       dvdx = dvdx / hx;
834:       dvdy = dvdy / hy;
835:       areadivf = area / fup;
836:       areadivf3 = areadivf / (fup * fup);
837:       smallH[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
838:       smallH[1][2] += areadivf3 * (-dvdy) * dvdx;
839:       smallH[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
840:       smallH[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
841:       smallH[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
842:       smallH[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);

844:       smallH[1][0] = smallH[0][1];
845:       smallH[2][0] = smallH[0][2];
846:       smallH[3][0] = smallH[0][3];
847:       smallH[2][1] = smallH[1][2];
848:       smallH[3][1] = smallH[1][3];
849:       smallH[3][2] = smallH[2][3];

851:       ind[0] = (j-gys) * gxm + (i-gxs);
852:       ind[1] = ind[0] + 1;
853:       ind[2] = ind[0] + gxm;
854:       ind[3] = ind[2] + 1;
855:       info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);

857:     }
858:   }

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

862:   info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
863:   info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
864:   info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);

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

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

871: } /* WholeMSurfHessian */