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 */