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