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: }