Actual source code: adelement2.c

      1: #ifdef ad_GRAD_MAX
      2: #undef ad_GRAD_MAX
      3: #endif
      4: #define ad_GRAD_MAX 16

     6:  #include taodaapplication.h
     7:  #include taoapp.h

      9: //#include "ad_deriv.h"

     11: typedef struct {
     12:   PetscScalar x,y,vpotx,vpoty;
     13: } Field;

     16: /* Function Gradient Evaluation over single element  */
     17: typedef struct {
     18:   int (*computeadicfunctiongradient)(int[2],DERIV_TYPE[16],DERIV_TYPE*,void*);
     19:   void *elementfgctx;
     20:   int elementfgflops;
     21:   DERIV_TYPE adX[16];
     22: } TaoDA2D4DOFADICFGCtx;

     27: static int TaoDA2dLoopADFunctionGradient(TAO_APPLICATION tao, DA da, Vec X, double *f, Vec G, void * ctx) {

     29:   TaoDA2D4DOFADICFGCtx *myapp = (TaoDA2D4DOFADICFGCtx*) ctx;
     30:   MPI_Comm comm;
     31:   Vec localX, localG;
     32:   int info, i, j, coor[2];
     33:   int xs, xm, mx, xe, ys, ym, ye, my;
     34:   PetscScalar zero=0.0, floc = 0.0;
     35:   DERIV_TYPE adF,*adX=myapp->adX;
     36:   Field     **x,**g;

     38:   info = DAGetLocalVector(da, &localX); CHKERRQ(info);
     39:   info = DAGetLocalVector(da, &localG); CHKERRQ(info);
     40:   info = VecSet(G, zero); CHKERRQ(info);
     41:   info = VecSet(localG, zero); CHKERRQ(info);

     43:   info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
     44:                    PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);

     46:   info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
     47:   info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);

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

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

     54:   xe = xs + xm;
     55:   ye = ys + ym;

     57:   for (j = ys; j < ye; j++) {
     58:     for (i = xs; i < xe; i++) {

     60:         DERIV_val(adX[0]) = x[j][i].x;
     61:         DERIV_val(adX[1]) = x[j][i].y;
     62:         DERIV_val(adX[2]) = x[j][i].vpotx;
     63:         DERIV_val(adX[3]) = x[j][i].vpoty;

     65:         DERIV_val(adX[4]) = x[j][i+1].x;
     66:         DERIV_val(adX[5]) = x[j][i+1].y;
     67:         DERIV_val(adX[6]) = x[j][i+1].vpotx;
     68:         DERIV_val(adX[7]) = x[j][i+1].vpoty;

     70:         DERIV_val(adX[8]) = x[j+1][i].x;
     71:         DERIV_val(adX[9]) = x[j+1][i].y;
     72:         DERIV_val(adX[10]) = x[j+1][i].vpotx;

     74:         DERIV_val(adX[11]) = x[j+1][i].vpoty;

     76:         DERIV_val(adX[12]) = x[j+1][i+1].x;
     77:         DERIV_val(adX[13]) = x[j+1][i+1].y;
     78:         DERIV_val(adX[14]) = x[j+1][i+1].vpotx;
     79:         DERIV_val(adX[15]) = x[j+1][i+1].vpoty;

     81:         coor[0] = i; coor[1] = j;

     83:         info = myapp->computeadicfunctiongradient(coor,adX,&adF,myapp->elementfgctx);
     84:         CHKERRQ(info);

     86:         floc += DERIV_val(adF);

     88:         g[j][i].x += DERIV_grad(adF)[0];
     89:         g[j][i].y += DERIV_grad(adF)[1];
     90:         g[j][i].vpotx += DERIV_grad(adF)[2];
     91:         g[j][i].vpoty += DERIV_grad(adF)[3];

     93:         g[j][i+1].x += DERIV_grad(adF)[4];
     94:         g[j][i+1].y += DERIV_grad(adF)[5];
     95:         g[j][i+1].vpotx += DERIV_grad(adF)[6];
     96:         g[j][i+1].vpoty += DERIV_grad(adF)[7];

     98:         g[j+1][i].x += DERIV_grad(adF)[8];
     99:         g[j+1][i].y += DERIV_grad(adF)[9];
    100:         g[j+1][i].vpotx += DERIV_grad(adF)[10];

    102:         g[j+1][i].vpoty += DERIV_grad(adF)[11];

    104:         g[j+1][i+1].x += DERIV_grad(adF)[12];
    105:         g[j+1][i+1].y += DERIV_grad(adF)[13];
    106:         g[j+1][i+1].vpotx += DERIV_grad(adF)[14];
    107:         g[j+1][i+1].vpoty += DERIV_grad(adF)[15];

    109:     }
    110:   }

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

    116:   info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
    117:   info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
    118:   info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
    119:   info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);

    121:   PetscLogFlops((ye-ys)*(xe-xs)*(myapp->elementfgflops + 33));
    122:   PetscObjectGetComm((PetscObject)da,&comm); CHKERRQ(info);
    123:   info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, comm); CHKERRQ(info);

    125:   return(0);
    126: } /* TaoDA2dLoopADFunctionGradient */

    130: static int TaoShutDownADICQQQQ(void* ctx){
    132:   ad_AD_Final();
    133:   return(0);
    134: }

    138: /*@C
    139:    DAAppSetADElementFunctionGradient2 - Set routine that evaluates the
    140:    local part of a function on a 2-dimensional periodic DA with 4 degrees of freedom.

    142:    Collective on TAO_APPLICATION

    144:    Input Parameters:
    145: +  daapp - the TAO_DA_APPLICATION solver context
    146: .  funcgrad - local function gradient routine
    147: .  flops - the number of flops done performed in the funcgrad routine
    148: -  fgctx - [optional] user-defined context for private data for the evaluation.

    150:    Calling sequence of funcgrad:
    151: $     int funcgrad(int coordinates[2], DERIV_TYPE x[16], DERIV_TYPE *f, void* ctx)

    153: +    coord - the global coordinates [i j] in each direction of the DA
    154: .    x - the variables on the DA ( da[j][i], da[j][j+1], da[j+1][i], da[j+1][i+1] ) (bottom left, bottom right, top left, top right)
    155: .    g - the ADIC differentiated objective function with respect to each variable
    156: -    ctx - user defined context

    158:    Calling sequence of func before calling ADIC
    159: $     int func(int coordinates[2], double x[16], double *f, void* ctx)

    161: +    coord - the global coordinates [i j] in each direction of the DA
    162: .    x - the variables on the DA ( da[j][i], da[j][j+1], da[j+1][i], da[j+1][i+1] ) (bottom left, bottom right, top left, top right)
    163: .    f - the objective function
    164: -    ctx - user defined context

    167:    Level: intermediate

    169: Concepts: DA, gradient, ADIC

    171: .seealso: DAAppSetObjectiveAndGradientRoutine(), DAAppSetADElementFunctionGradient()
    172: @*/
    173: int DAAppSetADElementFunctionGradient2(TAO_APPLICATION daapplication,
    174:                                          int (*funcgrad)(int[2],DERIV_TYPE[16],DERIV_TYPE*,void*),
    175:                                          int flops, void *ctx){
    176:   int i,n,info;
    177:   int dim,dof,s;
    178:   DAStencilType st;
    179:   TaoDA2D4DOFADICFGCtx *fgctx;
    180:   DA da;

    183:   info=DAAppGetNumberOfDAGrids(daapplication,&n); CHKERRQ(info);
    184:   for (i=0;i<n;i++){
    185:     info = DAAppGetDA(daapplication, i, &da); CHKERRQ(info);
    186:     info = DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,&s,0,&st); CHKERRQ(info);
    187:     if (dim!=2){
    188:       SETERRQ(1,"TAO DA ERROR: DA must have dimension of 2");}
    189:     if (dof!=4){
    190:       SETERRQ(1,"TAO DA ERROR: DA must have exactly 4 degrees of freedom per node");}
    191:     if (s!=1){
    192:       SETERRQ(1,"TAO DA ERROR: DA stencil width must equal 1"); }
    193:     //    if (st!=DA_STENCIL_BOX){
    194:     //      SETERRQ(1,"TAO DA ERROR: DA stencil must be DA_STENCIL_BOX");}
    195:   }
    196:   PetscNew(TaoDA2D4DOFADICFGCtx,&fgctx);
    197:   ad_AD_Init(16);
    198:   //  ad_AD_Init(ad_GRAD_MAX);
    199:   ad_AD_SetIndepArray(fgctx->adX,16);
    200:   ad_AD_SetIndepDone();
    201:   fgctx->computeadicfunctiongradient=funcgrad;
    202:   fgctx->elementfgctx=ctx;
    203:   fgctx->elementfgflops=flops;
    204:   info = DAAppSetObjectiveAndGradientRoutine(daapplication, TaoDA2dLoopADFunctionGradient, (void*)fgctx);
    205:   CHKERRQ(info);
    206:   info = TaoAppSetDestroyRoutine(daapplication,TaoApplicationFreeMemory, (void*)fgctx); CHKERRQ(info);
    207:   info = TaoAppSetDestroyRoutine(daapplication, TaoShutDownADICQQQQ, 0); CHKERRQ(info);
    208:   info = PetscLogInfo((daapplication,"Set objective function pointer for TAO_DA_APPLICATION object.\n")); CHKERRQ(info);
    209:   return(0);
    210: }