Actual source code: plate2.c

      1: /*$Id: plate2.c 1.74 05/05/11 08:44:55-05:00 sarich@zorak.(none) $*/

      3: #include "petscda.h"
     4:  #include tao.h

      6: static  char help[] =
      7: "This example demonstrates use of the TAO package to \n\
      8: solve an unconstrained minimization problem.  This example is based on a \n\
      9: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain, \n\
     10: boundary values along the edges of the domain, and a plate represented by \n\
     11: lower boundary conditions, the objective is to find the\n\
     12: surface with the minimal area that satisfies the boundary conditions.\n\
     13: The command line options are:\n\
     14:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
     15:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
     16:   -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
     17:   -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
     18:   -bheight <ht>, where <ht> = height of the plate\n\
     19:   -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
     20:                for an average of the boundary conditions\n\n";

     22: /*T
     23:    Concepts: TAO - Solving a bound constrained minimization problem
     24:    Routines: TaoInitialize(); TaoFinalize();
     25:    Routines: TaoApplicationCreate(); TaoAppDestroy();
     26:    Routines: TaoAppSetObjectiveAndGradientRoutine(); TaoAppSetHessianRoutine();
     27:    Routines: TaoAppSetInitialSolutionVec(); TaoAppSetHessianMat();
     28:    Routines: TaoAppSetVariableBounds();
     29:    Routines: TaoCreate();  TaoDestroy();
     30:    Routines: TaoSetOptions();
     31:    Processors: n
     32: T*/

     35: /*
     36:    User-defined application context - contains data needed by the
     37:    application-provided call-back routines, FormFunctionGradient(),
     38:    FormHessian().
     39: */
     40: typedef struct {
     41:   /* problem parameters */
     42:   PetscReal      bheight;                  /* Height of plate under the surface */
     43:   int         mx, my;                   /* discretization in x, y directions */
     44:   int         bmx,bmy;                  /* Size of plate under the surface */
     45:   Vec         Bottom, Top, Left, Right; /* boundary values */
     46:
     47:   /* Working space */
     48:   Vec         localX, localV;           /* ghosted local vector */
     49:   DA          da;                       /* distributed array data structure */
     50:   Mat         H;
     51: } AppCtx;

     53: /* -------- User-defined Routines --------- */

     55: static int MSA_BoundaryConditions(AppCtx*);
     56: static int MSA_InitialPoint(AppCtx*,Vec);
     57: static int MSA_Plate(TAO_APPLICATION,Vec,Vec,void*);
     58: int FormFunctionGradient(TAO_APPLICATION,Vec,double*,Vec,void*);
     59: int FormHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure*,void*);

     63: int main( int argc, char **argv )
     64: {
     65:   int        info;                 /* used to check for functions returning nonzeros */
     66:   int        Nx, Ny;               /* number of processors in x- and y- directions */
     67:   int        m, N;                 /* number of local and global elements in vectors */
     68:   Vec        x;                    /* solution vector */
     69:   PetscTruth   flg;                /* A return variable when checking for user options */
     70:   TAO_SOLVER tao;                  /* TAO_SOLVER solver context */
     71:   TAO_APPLICATION plateapp;        /* The PETSc application */
     72:   TaoMethod  method = "tao_blmvm"; /* default minimization method */
     73:   double     ff,gnorm,cnorm;       /* iteration information */
     74:   int        iter;
     75:   ISLocalToGlobalMapping isltog;   /* local-to-global mapping object */
     76:   TaoTerminateReason reason;
     77:   AppCtx     user;                 /* user-defined work context */

     79:   /* Initialize PETSc, TAO */
     80:   PetscInitialize( &argc, &argv,(char *)0,help );
     81:   TaoInitialize( &argc, &argv,(char *)0,help );

     83:   /* Specify default dimension of the problem */
     84:   user.mx = 10; user.my = 10; user.bheight=0.1;

     86:   /* Check for any command line arguments that override defaults */
     87:   info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg); CHKERRQ(info);
     88:   info = PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg); CHKERRQ(info);
     89:   info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user.bheight,&flg); CHKERRQ(info);

     91:   user.bmx = user.mx/2; user.bmy = user.my/2;
     92:   info = PetscOptionsGetInt(PETSC_NULL,"-bmx",&user.bmx,&flg); CHKERRQ(info);
     93:   info = PetscOptionsGetInt(PETSC_NULL,"-bmy",&user.bmy,&flg); CHKERRQ(info);

     95:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
     96:   PetscPrintf(PETSC_COMM_WORLD,"mx:%d, my:%d, bmx:%d, bmy:%d, height:%4.2f\n",
     97:               user.mx,user.my,user.bmx,user.bmy,user.bheight);

     99:   /* Calculate any derived values from parameters */
    100:   N    = user.mx*user.my;

    102:   /* Let Petsc determine the dimensions of the local vectors */
    103:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

    105:   /*
    106:      A two dimensional distributed array will help define this problem,
    107:      which derives from an elliptic PDE on two dimensional domain.  From
    108:      the distributed array, Create the vectors.
    109:   */
    110:   info = DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
    111:                     user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da); CHKERRQ(info);

    113:   /*
    114:      Extract global and local vectors from DA; The local vectors are
    115:      used solely as work space for the evaluation of the function,
    116:      gradient, and Hessian.  Duplicate for remaining vectors that are
    117:      the same types.
    118:   */
    119:   info = DACreateGlobalVector(user.da,&x); CHKERRQ(info); /* Solution */
    120:   info = DACreateLocalVector(user.da,&user.localX); CHKERRQ(info);
    121:   info = VecDuplicate(user.localX,&user.localV); CHKERRQ(info);

    123:   /*
    124:      Create a matrix data structure to store the Hessian.
    125:      Here we (optionally) also
    126:      associate the local numbering scheme with the matrix so that
    127:      later we can use local indices for matrix assembly.  We could
    128:      alternatively use global indices for matrix assembly.
    129:   */
    130:   info = VecGetLocalSize(x,&m); CHKERRQ(info);
    131:   info = MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL,
    132:                          3,PETSC_NULL,&(user.H)); CHKERRQ(info);
    133:   info = MatSetOption(user.H,MAT_SYMMETRIC); CHKERRQ(info);

    135:   info = DAGetISLocalToGlobalMapping(user.da,&isltog); CHKERRQ(info);
    136:   info = MatSetLocalToGlobalMapping(user.H,isltog); CHKERRQ(info);

    139:   /* The TAO code begins here */

    141:   /*
    142:      Create TAO solver and set desired solution method
    143:      The method must either be 'tao_tron' or 'tao_blmvm'
    144:      If blmvm is used, then hessian function is not called.
    145:   */
    146:   info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
    147:   info = TaoApplicationCreate(MPI_COMM_WORLD,&plateapp); CHKERRQ(info);

    149:   /* Set initial solution guess; */
    150:   info = MSA_BoundaryConditions(&user); CHKERRQ(info);
    151:   info = MSA_InitialPoint(&user,x); CHKERRQ(info);
    152:   info = TaoAppSetInitialSolutionVec(plateapp,x); CHKERRQ(info);

    154:   /* Set routines for function, gradient and hessian evaluation */
    155:   info = TaoAppSetObjectiveAndGradientRoutine(plateapp,FormFunctionGradient,(void*) &user);
    156:   CHKERRQ(info);
    157:
    158:   info = TaoAppSetHessianMat(plateapp,user.H,user.H); CHKERRQ(info);
    159:   info = TaoAppSetHessianRoutine(plateapp,FormHessian,(void*)&user);
    160:   CHKERRQ(info);

    162:   /* Set Variable bounds */
    163:   info = TaoAppSetVariableBoundsRoutine(plateapp,MSA_Plate,(void*)&user); CHKERRQ(info);

    165:   /* Check for any tao command line options */
    166:   info = TaoSetOptions(plateapp,tao); CHKERRQ(info);

    168:   /* SOLVE THE APPLICATION */
    169:   info = TaoSolveApplication(plateapp,tao);  CHKERRQ(info);

    171:   /* Get information on termination */
    172:   info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,&cnorm,0,&reason); CHKERRQ(info);
    173:   if (reason <= 0){
    174:     PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
    175:     PetscPrintf(MPI_COMM_WORLD,"Iteration: %d, f: %4.2e, Residual: %4.2e, Infeas: %4.2e\n",iter,ff,gnorm,cnorm);
    176:   }
    177:   /* Free TAO data structures */
    178:   info = TaoDestroy(tao); CHKERRQ(info);
    179:   info = TaoAppDestroy(plateapp); CHKERRQ(info);

    181:   /* Free PETSc data structures */
    182:   info = VecDestroy(x); CHKERRQ(info);
    183:   info = MatDestroy(user.H); CHKERRQ(info);
    184:   info = VecDestroy(user.localX); CHKERRQ(info);
    185:   info = VecDestroy(user.localV); CHKERRQ(info);
    186:   info = VecDestroy(user.Bottom); CHKERRQ(info);
    187:   info = VecDestroy(user.Top); CHKERRQ(info);
    188:   info = VecDestroy(user.Left); CHKERRQ(info);
    189:   info = VecDestroy(user.Right); CHKERRQ(info);
    190:   info = DADestroy(user.da); CHKERRQ(info);

    192:   /* Finalize TAO and PETSc */
    193:   PetscFinalize();
    194:   TaoFinalize();
    195:
    196:   return 0;
    197: }

    201: /*  FormFunctionGradient - Evaluates f(x) and gradient g(x).

    203:     Input Parameters:
    204: .   taoapp     - the TAO_APPLICATION context
    205: .   X      - input vector
    206: .   userCtx - optional user-defined context, as set by TaoSetPetscFunctionGradient()
    207:
    208:     Output Parameters:
    209: .   fcn     - the function value
    210: .   G      - vector containing the newly evaluated gradient

    212:    Notes:
    213:    In this case, we discretize the domain and Create triangles. The
    214:    surface of each triangle is planar, whose surface area can be easily
    215:    computed.  The total surface area is found by sweeping through the grid
    216:    and computing the surface area of the two triangles that have their
    217:    right angle at the grid point.  The diagonal line segments on the
    218:    grid that define the triangles run from top left to lower right.
    219:    The numbering of points starts at the lower left and runs left to
    220:    right, then bottom to top.
    221: */
    222: int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn, Vec G,void *userCtx)
    223: {
    224:   AppCtx * user = (AppCtx *) userCtx;
    225:   int    info,i,j,row;
    226:   int    mx=user->mx, my=user->my;
    227:   int    xs,xm,gxs,gxm,ys,ym,gys,gym;
    228:   double ft=0;
    229:   PetscReal zero=0.0;
    230:   PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
    231:   PetscReal rhx=mx+1, rhy=my+1;
    232:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
    233:   PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
    234:   PetscReal *g, *x,*left,*right,*bottom,*top;
    235:   Vec    localX = user->localX, localG = user->localV;

    237:   /* Get local mesh boundaries */
    238:   info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
    239:   info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);

    241:   /* Scatter ghost points to local vector */
    242:   info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
    243:   info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);

    245:   /* Initialize vector to zero */
    246:   info = VecSet(localG, zero); CHKERRQ(info);

    248:   /* Get pointers to vector data */
    249:   info = VecGetArray(localX,&x); CHKERRQ(info);
    250:   info = VecGetArray(localG,&g); CHKERRQ(info);
    251:   info = VecGetArray(user->Top,&top); CHKERRQ(info);
    252:   info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
    253:   info = VecGetArray(user->Left,&left); CHKERRQ(info);
    254:   info = VecGetArray(user->Right,&right); CHKERRQ(info);

    256:   /* Compute function over the locally owned part of the mesh */
    257:   for (j=ys; j<ys+ym; j++){
    258:     for (i=xs; i< xs+xm; i++){
    259:       row=(j-gys)*gxm + (i-gxs);
    260:
    261:       xc = x[row];
    262:       xlt=xrb=xl=xr=xb=xt=xc;
    263:
    264:       if (i==0){ /* left side */
    265:         xl= left[j-ys+1];
    266:         xlt = left[j-ys+2];
    267:       } else {
    268:         xl = x[row-1];
    269:       }

    271:       if (j==0){ /* bottom side */
    272:         xb=bottom[i-xs+1];
    273:         xrb = bottom[i-xs+2];
    274:       } else {
    275:         xb = x[row-gxm];
    276:       }
    277:
    278:       if (i+1 == gxs+gxm){ /* right side */
    279:         xr=right[j-ys+1];
    280:         xrb = right[j-ys];
    281:       } else {
    282:         xr = x[row+1];
    283:       }

    285:       if (j+1==gys+gym){ /* top side */
    286:         xt=top[i-xs+1];
    287:         xlt = top[i-xs];
    288:       }else {
    289:         xt = x[row+gxm];
    290:       }

    292:       if (i>gxs && j+1<gys+gym){
    293:         xlt = x[row-1+gxm];
    294:       }
    295:       if (j>gys && i+1<gxs+gxm){
    296:         xrb = x[row+1-gxm];
    297:       }

    299:       d1 = (xc-xl);
    300:       d2 = (xc-xr);
    301:       d3 = (xc-xt);
    302:       d4 = (xc-xb);
    303:       d5 = (xr-xrb);
    304:       d6 = (xrb-xb);
    305:       d7 = (xlt-xl);
    306:       d8 = (xt-xlt);
    307:
    308:       df1dxc = d1*hydhx;
    309:       df2dxc = ( d1*hydhx + d4*hxdhy );
    310:       df3dxc = d3*hxdhy;
    311:       df4dxc = ( d2*hydhx + d3*hxdhy );
    312:       df5dxc = d2*hydhx;
    313:       df6dxc = d4*hxdhy;

    315:       d1 *= rhx;
    316:       d2 *= rhx;
    317:       d3 *= rhy;
    318:       d4 *= rhy;
    319:       d5 *= rhy;
    320:       d6 *= rhx;
    321:       d7 *= rhy;
    322:       d8 *= rhx;

    324:       f1 = sqrt( 1.0 + d1*d1 + d7*d7);
    325:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
    326:       f3 = sqrt( 1.0 + d3*d3 + d8*d8);
    327:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);
    328:       f5 = sqrt( 1.0 + d2*d2 + d5*d5);
    329:       f6 = sqrt( 1.0 + d4*d4 + d6*d6);
    330:
    331:       ft = ft + (f2 + f4);

    333:       df1dxc /= f1;
    334:       df2dxc /= f2;
    335:       df3dxc /= f3;
    336:       df4dxc /= f4;
    337:       df5dxc /= f5;
    338:       df6dxc /= f6;

    340:       g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
    341:
    342:     }
    343:   }

    346:   /* Compute triangular areas along the border of the domain. */
    347:   if (xs==0){ /* left side */
    348:     for (j=ys; j<ys+ym; j++){
    349:       d3=(left[j-ys+1] - left[j-ys+2])*rhy;
    350:       d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
    351:       ft = ft+sqrt( 1.0 + d3*d3 + d2*d2);
    352:     }
    353:   }
    354:   if (ys==0){ /* bottom side */
    355:     for (i=xs; i<xs+xm; i++){
    356:       d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
    357:       d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
    358:       ft = ft+sqrt( 1.0 + d3*d3 + d2*d2);
    359:     }
    360:   }

    362:   if (xs+xm==mx){ /* right side */
    363:     for (j=ys; j< ys+ym; j++){
    364:       d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
    365:       d4=(right[j-ys]-right[j-ys+1])*rhy;
    366:       ft = ft+sqrt( 1.0 + d1*d1 + d4*d4);
    367:     }
    368:   }
    369:   if (ys+ym==my){ /* top side */
    370:     for (i=xs; i<xs+xm; i++){
    371:       d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
    372:       d4=(top[i-xs+1] - top[i-xs])*rhx;
    373:       ft = ft+sqrt( 1.0 + d1*d1 + d4*d4);
    374:     }
    375:   }

    377:   if (ys==0 && xs==0){
    378:     d1=(left[0]-left[1])*rhy;
    379:     d2=(bottom[0]-bottom[1])*rhx;
    380:     ft +=sqrt( 1.0 + d1*d1 + d2*d2);
    381:   }
    382:   if (ys+ym == my && xs+xm == mx){
    383:     d1=(right[ym+1] - right[ym])*rhy;
    384:     d2=(top[xm+1] - top[xm])*rhx;
    385:     ft +=sqrt( 1.0 + d1*d1 + d2*d2);
    386:   }

    388:   ft=ft*area;
    389:   info = MPI_Allreduce(&ft,fcn,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);CHKERRQ(info);

    391:
    392:   /* Restore vectors */
    393:   info = VecRestoreArray(localX,&x); CHKERRQ(info);
    394:   info = VecRestoreArray(localG,&g); CHKERRQ(info);
    395:   info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
    396:   info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
    397:   info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
    398:   info = VecRestoreArray(user->Right,&right); CHKERRQ(info);

    400:   /* Scatter values to global vector */
    401:   info = DALocalToGlobal(user->da,localG,INSERT_VALUES,G); CHKERRQ(info);

    403:   info = PetscLogFlops(70*xm*ym); CHKERRQ(info);

    405:   return 0;
    406: }

    408: /* ------------------------------------------------------------------- */
    411: /*
    412:    FormHessian - Evaluates Hessian matrix.

    414:    Input Parameters:
    415: .  taoapp  - the TAO_APPLICATION context
    416: .  x    - input vector
    417: .  ptr  - optional user-defined context, as set by TaoSetHessian()

    419:    Output Parameters:
    420: .  A    - Hessian matrix
    421: .  B    - optionally different preconditioning matrix
    422: .  flag - flag indicating matrix structure

    424:    Notes:
    425:    Due to mesh point reordering with DAs, we must always work
    426:    with the local mesh points, and then transform them to the new
    427:    global numbering with the local-to-global mapping.  We cannot work
    428:    directly with the global numbers for the original uniprocessor mesh!

    430:    Two methods are available for imposing this transformation
    431:    when setting matrix entries:
    432:      (A) MatSetValuesLocal(), using the local ordering (including
    433:          ghost points!)
    434:          - Do the following two steps once, before calling TaoSolve()
    435:            - Use DAGetISLocalToGlobalMapping() to extract the
    436:              local-to-global map from the DA
    437:            - Associate this map with the matrix by calling
    438:              MatSetLocalToGlobalMapping()
    439:          - Then set matrix entries using the local ordering
    440:            by calling MatSetValuesLocal()
    441:      (B) MatSetValues(), using the global ordering
    442:          - Use DAGetGlobalIndices() to extract the local-to-global map
    443:          - Then apply this map explicitly yourself
    444:          - Set matrix entries using the global ordering by calling
    445:            MatSetValues()
    446:    Option (A) seems cleaner/easier in many cases, and is the procedure
    447:    used in this example.
    448: */
    449: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *Hptr, Mat *Hpc, MatStructure *flag, void *ptr)
    450: {
    451:   int    info;
    452:   AppCtx *user = (AppCtx *) ptr;
    453:   Mat Hessian = *Hpc;
    454:   int    i,j,k,row;
    455:   int    mx=user->mx, my=user->my;
    456:   int    xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
    457:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
    458:   PetscReal rhx=mx+1, rhy=my+1;
    459:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
    460:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
    461:   PetscReal *x,*left,*right,*bottom,*top;
    462:   PetscReal v[7];
    463:   Vec    localX = user->localX;
    464:   PetscTruth assembled;

    467:   /* Set various matrix options */
    468:   info = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES); CHKERRQ(info);
    469:   info = MatSetOption(Hessian,MAT_COLUMNS_SORTED); CHKERRQ(info);
    470:   info = MatSetOption(Hessian,MAT_ROWS_SORTED); CHKERRQ(info);

    472:   /* Initialize matrix entries to zero */
    473:   info = MatAssembled(Hessian,&assembled); CHKERRQ(info);
    474:   if (assembled){info = MatZeroEntries(Hessian);  CHKERRQ(info);}

    476:   /* Get local mesh boundaries */
    477:   info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
    478:   info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);

    480:   /* Scatter ghost points to local vector */
    481:   info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
    482:   info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);

    484:   /* Get pointers to vector data */
    485:   info = VecGetArray(localX,&x); CHKERRQ(info);
    486:   info = VecGetArray(user->Top,&top); CHKERRQ(info);
    487:   info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
    488:   info = VecGetArray(user->Left,&left); CHKERRQ(info);
    489:   info = VecGetArray(user->Right,&right); CHKERRQ(info);

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

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

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

    497:       row=(j-gys)*gxm + (i-gxs);
    498:
    499:       xc = x[row];
    500:       xlt=xrb=xl=xr=xb=xt=xc;

    502:       /* Left side */
    503:       if (i==gxs){
    504:         xl= left[j-ys+1];
    505:         xlt = left[j-ys+2];
    506:       } else {
    507:         xl = x[row-1];
    508:       }
    509:
    510:       if (j==gys){
    511:         xb=bottom[i-xs+1];
    512:         xrb = bottom[i-xs+2];
    513:       } else {
    514:         xb = x[row-gxm];
    515:       }
    516:
    517:       if (i+1 == gxs+gxm){
    518:         xr=right[j-ys+1];
    519:         xrb = right[j-ys];
    520:       } else {
    521:         xr = x[row+1];
    522:       }

    524:       if (j+1==gys+gym){
    525:         xt=top[i-xs+1];
    526:         xlt = top[i-xs];
    527:       }else {
    528:         xt = x[row+gxm];
    529:       }

    531:       if (i>gxs && j+1<gys+gym){
    532:         xlt = x[row-1+gxm];
    533:       }
    534:       if (j>gys && i+1<gxs+gxm){
    535:         xrb = x[row+1-gxm];
    536:       }

    539:       d1 = (xc-xl)*rhx;
    540:       d2 = (xc-xr)*rhx;
    541:       d3 = (xc-xt)*rhy;
    542:       d4 = (xc-xb)*rhy;
    543:       d5 = (xrb-xr)*rhy;
    544:       d6 = (xrb-xb)*rhx;
    545:       d7 = (xlt-xl)*rhy;
    546:       d8 = (xlt-xt)*rhx;
    547:
    548:       f1 = sqrt( 1.0 + d1*d1 + d7*d7);
    549:       f2 = sqrt( 1.0 + d1*d1 + d4*d4);
    550:       f3 = sqrt( 1.0 + d3*d3 + d8*d8);
    551:       f4 = sqrt( 1.0 + d3*d3 + d2*d2);
    552:       f5 = sqrt( 1.0 + d2*d2 + d5*d5);
    553:       f6 = sqrt( 1.0 + d4*d4 + d6*d6);

    556:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
    557:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
    558:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
    559:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
    560:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
    561:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
    562:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
    563:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

    568:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
    569:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
    570:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
    571:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

    573:       hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5;  hc*=0.5;

    575:       k=0;
    576:       if (j>0){
    577:         v[k]=hb; col[k]=row - gxm; k++;
    578:       }
    579:
    580:       if (j>0 && i < mx -1){
    581:         v[k]=hbr; col[k]=row - gxm+1; k++;
    582:       }
    583:
    584:       if (i>0){
    585:         v[k]= hl; col[k]=row - 1; k++;
    586:       }
    587:
    588:       v[k]= hc; col[k]=row; k++;
    589:
    590:       if (i < mx-1 ){
    591:         v[k]= hr; col[k]=row+1; k++;
    592:       }
    593:
    594:       if (i>0 && j < my-1 ){
    595:         v[k]= htl; col[k] = row+gxm-1; k++;
    596:       }
    597:
    598:       if (j < my-1 ){
    599:         v[k]= ht; col[k] = row+gxm; k++;
    600:       }
    601:
    602:       /*
    603:          Set matrix values using local numbering, which was defined
    604:          earlier, in the main routine.
    605:       */
    606:       info = MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);
    607:       CHKERRQ(info);
    608:
    609:     }
    610:   }
    611:
    612:   /* Restore vectors */
    613:   info = VecRestoreArray(localX,&x); CHKERRQ(info);
    614:   info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
    615:   info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
    616:   info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
    617:   info = VecRestoreArray(user->Right,&right); CHKERRQ(info);

    619:   /* Assemble the matrix */
    620:   info = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    621:   info = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);

    623:   info = PetscLogFlops(199*xm*ym); CHKERRQ(info);
    624:   return 0;
    625: }

    627: /* ------------------------------------------------------------------- */
    630: /*
    631:    MSA_BoundaryConditions -  Calculates the boundary conditions for
    632:    the region.

    634:    Input Parameter:
    635: .  user - user-defined application context

    637:    Output Parameter:
    638: .  user - user-defined application context
    639: */
    640: static int MSA_BoundaryConditions(AppCtx * user)
    641: {
    642:   int        i,j,k,limit=0,info,maxits=5;
    643:   int        xs,ys,xm,ym,gxs,gys,gxm,gym;
    644:   int        mx=user->mx,my=user->my;
    645:   int        bsize=0, lsize=0, tsize=0, rsize=0;
    646:   PetscReal     one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
    647:   PetscReal     fnorm,det,hx,hy,xt=0,yt=0;
    648:   PetscReal     u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
    649:   PetscReal     b=-0.5, t=0.5, l=-0.5, r=0.5;
    650:   PetscReal     *boundary;
    651:   PetscTruth   flg;
    652:   Vec        Bottom,Top,Right,Left;

    654:   /* Get local mesh boundaries */
    655:   info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
    656:   info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);

    659:   bsize=xm+2;
    660:   lsize=ym+2;
    661:   rsize=ym+2;
    662:   tsize=xm+2;

    664:   info = VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom); CHKERRQ(info);
    665:   info = VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top); CHKERRQ(info);
    666:   info = VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left); CHKERRQ(info);
    667:   info = VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right); CHKERRQ(info);

    669:   user->Top=Top;
    670:   user->Left=Left;
    671:   user->Bottom=Bottom;
    672:   user->Right=Right;

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

    676:   for (j=0; j<4; j++){
    677:     if (j==0){
    678:       yt=b;
    679:       xt=l+hx*xs;
    680:       limit=bsize;
    681:       VecGetArray(Bottom,&boundary);
    682:     } else if (j==1){
    683:       yt=t;
    684:       xt=l+hx*xs;
    685:       limit=tsize;
    686:       VecGetArray(Top,&boundary);
    687:     } else if (j==2){
    688:       yt=b+hy*ys;
    689:       xt=l;
    690:       limit=lsize;
    691:       VecGetArray(Left,&boundary);
    692:     } else if (j==3){
    693:       yt=b+hy*ys;
    694:       xt=r;
    695:       limit=rsize;
    696:       VecGetArray(Right,&boundary);
    697:     }

    699:     for (i=0; i<limit; i++){
    700:       u1=xt;
    701:       u2=-yt;
    702:       for (k=0; k<maxits; k++){
    703:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
    704:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
    705:         fnorm=sqrt(nf1*nf1+nf2*nf2);
    706:         if (fnorm <= tol) break;
    707:         njac11=one+u2*u2-u1*u1;
    708:         njac12=two*u1*u2;
    709:         njac21=-two*u1*u2;
    710:         njac22=-one - u1*u1 + u2*u2;
    711:         det = njac11*njac22-njac21*njac12;
    712:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
    713:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
    714:       }

    716:       boundary[i]=u1*u1-u2*u2;
    717:       if (j==0 || j==1) {
    718:         xt=xt+hx;
    719:       } else if (j==2 || j==3){
    720:         yt=yt+hy;
    721:       }
    722:
    723:     }
    724:
    725:     if (j==0){
    726:       info = VecRestoreArray(Bottom,&boundary); CHKERRQ(info);
    727:     } else if (j==1){
    728:       info = VecRestoreArray(Top,&boundary); CHKERRQ(info);
    729:     } else if (j==2){
    730:       info = VecRestoreArray(Left,&boundary); CHKERRQ(info);
    731:     } else if (j==3){
    732:       info = VecRestoreArray(Right,&boundary); CHKERRQ(info);
    733:     }

    735:   }

    737:   /* Scale the boundary if desired */

    739:   info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
    740:   CHKERRQ(info);
    741:   if (flg){
    742:     info = VecScale(Bottom, scl); CHKERRQ(info);
    743:   }
    744:
    745:   info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
    746:   CHKERRQ(info);
    747:   if (flg){
    748:     info = VecScale(Top, scl); CHKERRQ(info);
    749:   }
    750:
    751:   info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
    752:   CHKERRQ(info);
    753:   if (flg){
    754:     info = VecScale(Right, scl); CHKERRQ(info);
    755:   }
    756:
    757:   info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
    758:   CHKERRQ(info);
    759:   if (flg){
    760:     info = VecScale(Left, scl); CHKERRQ(info);
    761:   }

    763:   return 0;
    764: }

    767: /* ------------------------------------------------------------------- */
    770: /*
    771:    MSA_Plate -  Calculates an obstacle for surface to stretch over.

    773:    Input Parameter:
    774: .  user - user-defined application context

    776:    Output Parameter:
    777: .  user - user-defined application context
    778: */
    779: static int MSA_Plate(TAO_APPLICATION taoapp, Vec XL,Vec XU,void *ctx){

    781:   AppCtx *user=(AppCtx *)ctx;
    782:   int i,j,row,info;
    783:   int xs,ys,xm,ym;
    784:   int mx=user->mx, my=user->my, bmy, bmx;
    785:   PetscReal t1,t2,t3;
    786:   PetscReal *xl, lb=TAO_NINFINITY, ub=TAO_INFINITY;
    787:   PetscTruth cylinder;

    789:   user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
    790:   user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
    791:   bmy=user->bmy, bmx=user->bmx;

    793:   info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);

    795:   info = VecSet(XL, lb); CHKERRQ(info);
    796:   info = VecSet(XU, ub); CHKERRQ(info);

    798:   info = VecGetArray(XL,&xl); CHKERRQ(info);

    800:   info = PetscOptionsHasName(PETSC_NULL,"-cylinder",&cylinder); CHKERRQ(info);
    801:   /* Compute the optional lower box */
    802:   if (cylinder){
    803:     for (i=xs; i< xs+xm; i++){
    804:       for (j=ys; j<ys+ym; j++){
    805:         row=(j-ys)*xm + (i-xs);
    806:         t1=(2.0*i-mx)*bmy;
    807:         t2=(2.0*j-my)*bmx;
    808:         t3=bmx*bmx*bmy*bmy;
    809:         if ( t1*t1 + t2*t2 <= t3 ){
    810:           xl[row] = user->bheight;
    811:         }
    812:       }
    813:     }
    814:   } else {
    815:     /* Compute the optional lower box */
    816:     for (i=xs; i< xs+xm; i++){
    817:       for (j=ys; j<ys+ym; j++){
    818:         row=(j-ys)*xm + (i-xs);
    819:         if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
    820:             j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
    821:           xl[row] = user->bheight;
    822:         }
    823:       }
    824:     }
    825:   }
    826:     info = VecRestoreArray(XL,&xl); CHKERRQ(info);

    828:   return 0;
    829: }

    832: /* ------------------------------------------------------------------- */
    835: /*
    836:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

    838:    Input Parameters:
    839: .  user - user-defined application context
    840: .  X - vector for initial guess

    842:    Output Parameters:
    843: .  X - newly computed initial guess
    844: */
    845: static int MSA_InitialPoint(AppCtx * user, Vec X)
    846: {
    847:   int      start=-1,i,j,info;
    848:   PetscReal   zero=0.0;
    849:   PetscTruth flg;

    851:   info = PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); CHKERRQ(info);

    853:   if (flg && start==0){ /* The zero vector is reasonable */
    854:
    855:     info = VecSet(X, zero); CHKERRQ(info);
    856:     /* PetscLogInfo((user,"Min. Surface Area Problem: Start with 0 vector \n")); */

    858:   } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */

    860:     PetscRandom rctx;  PetscReal np5=-0.5;

    862:     info = PetscRandomCreate(MPI_COMM_WORLD,RANDOM_DEFAULT,&rctx);
    863:     CHKERRQ(info);
    864:     for (i=0; i<start; i++){
    865:       info = VecSetRandom(X, rctx); CHKERRQ(info);
    866:     }
    867:     info = PetscRandomDestroy(rctx); CHKERRQ(info);
    868:     info = VecShift(X, np5);
    869:     /* PetscLogInfo((user,"Min. Surface Area Problem: Start with random vector %d\n",start)); */

    871:   } else { /* Take an average of the boundary conditions */

    873:     int    row,xs,xm,gxs,gxm,ys,ym,gys,gym;
    874:     int    mx=user->mx,my=user->my;
    875:     PetscReal *x,*left,*right,*bottom,*top;
    876:     Vec    localX = user->localX;
    877:
    878:     /* Get local mesh boundaries */
    879:     info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
    880:     info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);
    881:
    882:     /* Get pointers to vector data */
    883:     info = VecGetArray(user->Top,&top); CHKERRQ(info);
    884:     info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info);
    885:     info = VecGetArray(user->Left,&left); CHKERRQ(info);
    886:     info = VecGetArray(user->Right,&right); CHKERRQ(info);

    888:     info = VecGetArray(localX,&x); CHKERRQ(info);
    889:     /* Perform local computations */
    890:     for (j=ys; j<ys+ym; j++){
    891:       for (i=xs; i< xs+xm; i++){
    892:         row=(j-gys)*gxm + (i-gxs);
    893:         x[row] = ( (j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+
    894:                    (i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
    895:       }
    896:     }
    897:
    898:     /* Restore vectors */
    899:     info = VecRestoreArray(localX,&x); CHKERRQ(info);

    901:     info = VecRestoreArray(user->Left,&left); CHKERRQ(info);
    902:     info = VecRestoreArray(user->Top,&top); CHKERRQ(info);
    903:     info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info);
    904:     info = VecRestoreArray(user->Right,&right); CHKERRQ(info);
    905:
    906:     /* Scatter values into global vector */
    907:     info = DALocalToGlobal(user->da,localX,INSERT_VALUES,X); CHKERRQ(info);
    908:
    909:   }
    910:   return 0;
    911: }