Actual source code: jbearing3.c
1: /*$Id: jbearing.c, v 1.1 2002/08/08 11:35 lopezca@mauddib.mcs.anl.gov $*/
3: /* Program usage: mpirun -np <proc> jbearing [-help] [all TAO options] */
5: /*
6: Include "tao.h" so we can use TAO solvers.
7: petscda.h for distributed array
8: ad_deriv.h for AD gradient
9: */
11: #include "petscda.h"
12: #include "tao.h"
13: #include "taodaapplication.h"
15: static char help[] ="Pressure distribution in a Journal Bearing. \n\
16: This example is based on the problem DPJB from the MINPACK-2 test suite.\n\
17: This pressure journal bearing problem is an example of elliptic variational\n\
18: problem defined over a two dimensional rectangle. By discretizing the domain \n\
19: into triangular elements, the pressure surrounding the journal bearing is\n\
20: defined as the minimum of a quadratic function whose variables are bounded\n\
21: below by zero. The command line options are:\n\
22: -ecc <ecc>, where <ecc> = epsilon parameter\n\
23: -b <b>, where <b> = half the upper limit in the 2nd coordinate direction\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\n";
30: /*T
31: Concepts: TAO - Solving a bounded minimization problem
32: Routines: TaoInitialize(); TaoFinalize();
33: Routines: TaoCreate(); TaoDestroy();
34: Routines: DAApplicationCreate(); DAApplicationDestroy();
35: Routines: DAAppSetVariableBoundsRoutine;
36: Routines: DAAppSetElementObjectiveAndGradientRoutine();
37: Routines: DAAppSetElementHessianRoutine();
38: Routines: DAAppSetObjectiveAndGradientRoutine();
39: Routines: DAAppSetADElementFunctionGradient();
40: Routines: DAAppSetHessianRoutine();
41: Routines: TaoSetOptions();
42: Routines: TaoGetSolutionStatus(); TaoDAAppSolve();
43: Routines: DAAppSetBeforeMonitor(); DAAppSetAfterMonitor
44: Routines: DAAppGetSolution(); TaoView();
45: Routines: DAAppGetInterpolationMatrix();
46: Processors: n
47: T*/
48:
49: /*
50: User-defined application context - contains data needed by the
51: application-provided call-back routines.
52: */
54: int ad_JBearLocalFunction(PetscInt[2] ,DERIV_TYPE[4], DERIV_TYPE *, void*);
55: typedef struct {
57: InactiveDouble *wq, *wl; /* vectors with the parameters w_q(x) and w_l(x) */
58: InactiveDouble hx, hy; /* increment size in both directions */
59: InactiveDouble area; /* area of the triangles */
61: } ADFGCtx;
64: typedef struct {
66: PetscReal ecc; /* epsilon value */
67: PetscReal b; /* 0.5 * upper limit for 2nd variable */
68: double *wq, *wl; /* vectors with the parameters w_q(x) and w_l(x) */
69: double hx, hy; /* increment size in both directions */
70: double area; /* area of the triangles */
72: PetscInt mx, my; /* discretization including boundaries */
74: ADFGCtx fgctx; /* Used only when an ADIC generated gradient is used */
76: } AppCtx;
78: /* User-defined routines found in this file */
79: static int AppCtxInitialize(void *ptr);
80: static int FormInitialGuess(DA, Vec);
82: static int JBearLocalFunctionGradient(PetscInt[2], double x[4], double *f, double g[4], void *ptr);
83: static int JBearLocalHessian(PetscInt[2], double x[4], double H[4][4], void *ptr);
85: static int WholeJBearFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);
86: static int WholeJBearHessian(TAO_APPLICATION,DA,Vec,Mat,void*);
88: static int DASetBounds(TAO_APPLICATION, DA, Vec, Vec, void*);
90: static int MyGridMonitorBefore(TAO_APPLICATION, DA, PetscInt, void *);
91: static int MyGridMonitorAfter1(TAO_APPLICATION, DA, PetscInt, void *);
95: int main( int argc, char **argv ) {
97: PetscInt info,iter; /* used to check for functions returning nonzeros */
98: PetscInt nlevels; /* multigrid levels */
99: PetscInt Nx,Ny;
100: double ff,gnorm;
101: DA DAarray[20];
102: Vec X;
103: PetscTruth flg, PreLoad = PETSC_TRUE; /* flags */
104: AppCtx user; /* user-defined work context */
105: TaoMethod method = "tao_gpcg"; /* minimization method */
106: TAO_SOLVER tao; /* TAO_SOLVER solver context */
107: TAO_APPLICATION JBearApp; /* The PETSc application */
108: TaoTerminateReason reason;
109: KSP ksp;
110: PC pc;
112: /* Initialize TAO */
113: PetscInitialize(&argc, &argv, (char *)0, help);
114: TaoInitialize(&argc, &argv, (char *)0, help);
116: PreLoadBegin(PreLoad,"Solve");
117:
118: info = AppCtxInitialize((void*)&user); CHKERRQ(info);
119:
120: nlevels=5;
121: info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
122: if (PreLoadIt == 0) {
123: nlevels = 1; user.mx = 11; user.my = 11; }
125: PetscPrintf(MPI_COMM_WORLD,"\n---- Journal Bearing Problem -----\n\n");
127: /* Let PETSc determine the vector distribution */
128: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
130: /* Create distributed array (DA) to manage parallel grid and vectors */
131: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
132: user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
133: for (iter=1;iter<nlevels;iter++){
134: info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
135: }
137: /* Create TAO solver and set desired solution method */
138: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
139: info = TaoApplicationCreate(PETSC_COMM_WORLD, &JBearApp); CHKERRQ(info);
140: info = TaoAppSetDAApp(JBearApp,DAarray,nlevels); CHKERRQ(info);
142: /* Sets routine bounds evaluation */
143: info = DAAppSetVariableBoundsRoutine(JBearApp,DASetBounds,(void *)&user); CHKERRQ(info);
145: info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
146: if (flg) {
148: /* Sets routines for function and gradient evaluation, element by element */
149: info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
150: if (flg) {
151: info = DAAppSetADElementFunctionGradient(JBearApp,ad_JBearLocalFunction,248,(void *)&user.fgctx); CHKERRQ(info);
152: } else {
153: info = DAAppSetElementObjectiveAndGradientRoutine(JBearApp,JBearLocalFunctionGradient,63,(void *)&user); CHKERRQ(info);
154: }
155: /* Sets routines for Hessian evaluation, element by element */
156: info = DAAppSetElementHessianRoutine(JBearApp,JBearLocalHessian,16,(void*)&user); CHKERRQ(info);
158: } else {
160: /* Sets routines for function and gradient evaluation, all in one routine */
161: info = DAAppSetObjectiveAndGradientRoutine(JBearApp,WholeJBearFunctionGradient,(void *)&user); CHKERRQ(info);
163: /* Sets routines for Hessian evaluation, all in one routine */
164: info = DAAppSetHessianRoutine(JBearApp,WholeJBearHessian,(void*)&user); CHKERRQ(info);
165:
166: }
168: info = DAAppSetBeforeMonitor(JBearApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
169: info = DAAppSetAfterMonitor(JBearApp,MyGridMonitorAfter1,(void*)&user); CHKERRQ(info);
170: info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
171: if (flg){
172: info = DAAppPrintInterpolationError(JBearApp); CHKERRQ(info);
173: info = DAAppPrintStageTimes(JBearApp); CHKERRQ(info);
174: }
175: /* Check for any tao command line options */
176: info = TaoAppSetRelativeTolerance(JBearApp,1.0e-6); CHKERRQ(info);
177: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
178: info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);
180: info = TaoAppGetKSP(JBearApp,&ksp);CHKERRQ(info);
181: info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
182: info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);CHKERRQ(info);
183: info = KSPGetPC(ksp,&pc);CHKERRQ(info);
184: info = PCSetType(pc,PCBJACOBI);CHKERRQ(info);
186: info = TaoSetOptions(JBearApp,tao); CHKERRQ(info);
188: info = DAAppGetSolution(JBearApp,0,&X); CHKERRQ(info);
189: info = FormInitialGuess(DAarray[0],X); CHKERRQ(info);
190: info = DAAppSetInitialSolution(JBearApp,X); CHKERRQ(info);
191: /* SOLVE THE APPLICATION */
192: info = TaoDAAppSolve(JBearApp,tao); CHKERRQ(info);
194: /* Get information on termination */
195: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
196: if (reason <= 0 ){
197: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
198: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
199: }
201: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
202: if (flg){
203: info = DAAppGetSolution(JBearApp,nlevels-1,&X); CHKERRQ(info);
204: info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
205: }
207: /* To View TAO solver information */
208: // info = TaoView(tao); CHKERRQ(info);
210: /* Free TAO data structures */
211: info = TaoDestroy(tao); CHKERRQ(info);
212: info = TaoAppDestroy(JBearApp); CHKERRQ(info);
214: /* Free PETSc data structures */
215: for (iter=0;iter<nlevels;iter++){
216: info = DADestroy(DAarray[iter]); CHKERRQ(info);
217: }
219: PreLoadEnd();
221: /* Finalize TAO */
222: TaoFinalize();
223: PetscFinalize();
225: return 0;
226: } /* main */
230: /*----- The following two routines
231:
232: MyGridMonitorBefore MyGridMonitorAfter
234: help diplay info of iterations at every grid level -------*/
238: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx) {
240: AppCtx *user = (AppCtx*)ctx;
241: int info;
242: PetscInt mx,my;
243: double t;
245: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
246: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
247: user->mx = mx;
248: user->my = my;
250: user->hx = (2.0 * 3.14159265358979) / (user->mx - 1);
251: user->hy = (2.0 * user->b) / (user->my - 1);
252: user->area = 0.5 * user->hx * user->hy;
254: user->wq = new double[user->mx];
255: user->wl = new double[user->mx];
256: for (PetscInt i=0; i<user->mx; i++) {
257: t = 1.0 + user->ecc * cos(i*user->hx);
258: user->wq[i] = t*t*t;
259: user->wl[i] = user->ecc * sin(i*user->hx);
260: }
261: info = PetscLogFlops(8 + 7 * user->mx); CHKERRQ(info);
263: user->fgctx.hx = user->hx;
264: user->fgctx.hy = user->hy;
265: user->fgctx.area = user->area;
266: user->fgctx.wq = user->wq;
267: user->fgctx.wl = user->wl;
269: PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my);
271: return 0;
272: }
276: static int MyGridMonitorAfter1(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx){
277:
278: AppCtx *user = (AppCtx*)ctx;
279: delete [] user->wq; delete [] user->wl;
280: return 0;
281: }
285: /*------- USER-DEFINED: initialize the application context information -------*/
289: /*
290: AppCtxInitialize - Sets initial values for the application context parameters
292: Input:
293: ptr - void user-defined application context
295: Output:
296: ptr - user-defined application context with the default or user-provided
297: parameters
298: */
299: static int AppCtxInitialize(void *ptr) {
301: AppCtx *user = (AppCtx*)ptr;
302: PetscTruth flg; /* flag for PETSc calls */
303: int info;
305: /* Specify default parameters */
306: user->ecc = 0.1;
307: user->b = 10.0;
308: user->mx = user->my = 11;
310: /* Check for command line arguments that override defaults */
311: info = PetscOptionsGetReal(TAO_NULL, "-ecc", &user->ecc, &flg); CHKERRQ(info);
312: info = PetscOptionsGetReal(TAO_NULL, "-b", &user->b, &flg); CHKERRQ(info);
313: info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
314: info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
316: return 0;
317: } /* AppCtxInitialize */
322: static int FormInitialGuess(DA da, Vec X)
323: {
324: int info;
325: PetscInt i, j, mx;
326: PetscInt xs, ys, xm, ym, xe, ye;
327: PetscReal hx, val;
328: double **x;
330: /* Get local mesh boundaries */
331: info = DAGetInfo(da,PETSC_NULL,&mx,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
332: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
333: hx = 2.0*4.0*atan(1.0)/(mx-1);
335: info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
336: xe = xs+xm; ye = ys+ym;
338: info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
339: /* Compute initial guess over locally owned part of mesh */
340: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
341: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
342: val = PetscMax(sin((i+1.0)*hx),0.0);
343: x[j][i] = val;
344: x[j][i] = 0;
345: }
346: }
347: info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);
349: return 0;
350: }
353: /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/
356: /*
357: FormBounds - Forms bounds on the variables
359: Input:
360: user - user-defined application context
362: Output:
363: XL - vector of lower bounds
364: XU - vector of upper bounds
366: */
367: static int DASetBounds(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr) {
369: AppCtx *user = (AppCtx*)ptr;
370: int info;
371: PetscInt i, j, mx, my;
372: PetscInt xs, xm, ys, ym;
373: double **xl, **xu;
375: mx = user->mx;
376: my = user->my;
378: info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
379: info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
380: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
382: for (j = ys; j < ys+ym; j++){
383: for (i = xs; i < xs+xm; i++){
384: xl[j][i] = 0.0;
385: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
386: xu[j][i] = 0.0;
387: } else {
388: xu[j][i] = TAO_INFINITY;
389: }
390: }
391: }
393: info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
394: info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);
396: return 0;
398: } /* DASetBounds */
404: /*
405: JBearLocalFunctionGradient - Evaluates function and gradient over the
406: local rectangular element
408: Input:
409: coor - vector with the indices of the position of current element
410: in the first, second and third directions
411: x - current point (values over the current rectangular element)
412: ptr - user-defined application context
414: Output:
415: f - value of the objective funtion at the local rectangular element
416: g - gradient of the local function
418: */
419: static int JBearLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {
421: AppCtx *user = (AppCtx*)ptr;
423: double avgWq, sqGrad, avgWV, fl, fu;
424: double hx, hy, area, aread3, *wq, *wl;
425: double dvdx, dvdy;
426: PetscInt i;
428: hx = user->hx;
429: hy = user->hy;
430: area = user->area;
431: aread3 = area / 3.0;
432: wq = user->wq;
433: wl = user->wl;
434: i = coor[0];
436: /* lower triangle contribution */
437: dvdx = (x[0] - x[1]) / hx;
438: dvdy = (x[0] - x[2]) / hy;
439: sqGrad = dvdx * dvdx + dvdy * dvdy;
440: avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
441: avgWV = (wl[i]*x[0] + wl[i+1]*x[1] + wl[i]*x[2]) / 3.0;
442: fl = avgWq * sqGrad - avgWV;
444: dvdx = dvdx * hy * avgWq;
445: dvdy = dvdy * hx * avgWq;
446: g[0] = ( dvdx + dvdy ) - wl[i] * aread3;
447: g[1] = ( -dvdx ) - wl[i+1] * aread3;
448: g[2] = ( -dvdy ) - wl[i] * aread3;
450: /* upper triangle contribution */
451: dvdx = (x[3] - x[2]) / hx;
452: dvdy = (x[3] - x[1]) / hy;
453: sqGrad = dvdx * dvdx + dvdy * dvdy;
454: avgWq = (2.0 * wq[i+1] + wq[i]) / 6.0;
455: avgWV = (wl[i+1]*x[1] + wl[i]*x[2] + wl[i+1]*x[3]) / 3.0;
456: fu = avgWq * sqGrad - avgWV;
458: dvdx = dvdx * hy * avgWq;
459: dvdy = dvdy * hx * avgWq;
460: g[1] += (-dvdy) - wl[i+1] * aread3;
461: g[2] += (-dvdx) - wl[i] * aread3;
462: g[3] = ( dvdx + dvdy ) - wl[i+1] * aread3;
464: *f = area * (fl + fu);
466: return 0;
467: } /* JBearLocalFunctionGradient */
470: /*------- USER-DEFINED: routine to evaluate the Hessian
471: at a local (rectangular element) level -------*/
474: /*
475: JBearLocalHessian - Computes the Hessian of the local (partial) function
476: defined over the current rectangle
478: Input:
479: coor - vector with the indices of the position of current element
480: in the first, second and third directions
481: x - current local solution (over the rectangle only)
482: ptr - user-defined application context
484: Output:
485: H - Hessian matrix of the local function (wrt the four
486: points of the rectangle only)
488: */
489: static int JBearLocalHessian(PetscInt coor[2], double x[4], double H[4][4],void *ptr) {
491: AppCtx *user = (AppCtx*)ptr;
492: double wql, wqu;
493: double hx, hy, dydx, dxdy;
494: double *wq;
495: double wqldydx,wqldxdy,wqudydx,wqudxdy;
496: PetscInt i;
498: hx = user->hx;
499: hy = user->hy;
500: wq = user->wq;
501: i = coor[0];
503: dxdy = hx / hy;
504: dydx = hy / hx;
505: wql = (2.0 * wq[i] + wq[i+1]) / 6.0;
506: wqu = (wq[i] + 2.0 * wq[i+1]) / 6.0;
507: wqldydx = wql * dydx;
508: wqldxdy = wql * dxdy;
509: wqudydx = wqu * dydx;
510: wqudxdy = wqu * dxdy;
512: /* Hessian contribution at 0,0 */
513: H[0][0] = wqldxdy + wqldydx;
514: H[0][1] = H[1][0] = -wqldydx;
515: H[0][2] = H[2][0] = -wqldxdy;
516: H[0][3] = H[3][0] = 0.0;
518: /* Hessian contribution at 1,0 */
519: H[1][1] = wqldydx + wqudxdy;
520: H[1][2] = H[2][1] = 0.0;
521: H[1][3] = H[3][1] = -wqudxdy;
523: /* Hessian contribution at 0,1 */
524: H[2][2] = wqldxdy + wqudydx;
525: H[2][3] = H[3][2] = -wqudydx;
527: /* Hessian contribution at 1,1 */
528: H[3][3] = wqudydx + wqudxdy;
530: return 0;
532: } /* JBearLocalHessian */
535: /*------- USER-DEFINED: routine to evaluate the function
536: and gradient at the whole grid -------*/
540: /*
541: WholeJBearFunctionGradient - Evaluates function and gradient over the
542: whole grid
544: Input:
545: daapplication - TAO application object
546: da - distributed array
547: X - the current point, at which the function and gradient are evaluated
548: ptr - user-defined application context
550: Output:
551: f - value of the objective funtion at X
552: G - gradient at X
553: */
554: static int WholeJBearFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {
556: AppCtx *user = (AppCtx*)ptr;
557: Vec localX, localG;
558: int info;
559: PetscInt i, j;
560: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
561: double **x, **g;
562: double floc = 0.0;
563: PetscScalar zero = 0.0;
565: double avgWq, sqGrad, avgWV, fl, fu;
566: double hx, hy, area, aread3, *wq, *wl;
567: double dvdx, dvdy;
569: hx = user->hx;
570: hy = user->hy;
571: area = user->area;
572: aread3= area/3.0;
573: wq = user->wq;
574: wl = user->wl;
576: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
577: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
578: info = VecSet(G, zero); CHKERRQ(info);
579: info = VecSet(localG, zero); CHKERRQ(info);
581: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
582: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
584: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
585: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
587: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
588: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
590: xe = gxs + gxm - 1;
591: ye = gys + gym - 1;
592: for (j = ys; j < ye; j++) {
593: for (i = xs; i < xe; i++) {
595: /* lower triangle contribution */
596: dvdx = (x[j][i] - x[j][i+1]) / hx;
597: dvdy = (x[j][i] - x[j+1][i]) / hy;
598: sqGrad = dvdx * dvdx + dvdy * dvdy;
599: avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
600: avgWV = (wl[i]*x[j][i] + wl[i+1]*x[j][i+1] + wl[i]*x[j+1][i]) / 3.0;
601: fl = avgWq * sqGrad - avgWV;
603: dvdx = dvdx * hy * avgWq;
604: dvdy = dvdy * hx * avgWq;
605: g[j][i] += ( dvdx + dvdy ) - wl[i] * aread3;
606: g[j][i+1] += ( -dvdx ) - wl[i+1] * aread3;
607: g[j+1][i] += ( -dvdy ) - wl[i] * aread3;
609: /* upper triangle contribution */
610: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
611: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
612: sqGrad = dvdx * dvdx + dvdy * dvdy;
613: avgWq = (2.0 * wq[i+1] + wq[i]) / 6.0;
614: avgWV = (wl[i+1]*x[j][i+1] + wl[i]*x[j+1][i] + wl[i+1]*x[j+1][i+1]) / 3.0;
615: fu = avgWq * sqGrad - avgWV;
617: dvdx = dvdx * hy * avgWq;
618: dvdy = dvdy * hx * avgWq;
619: g[j][i+1] += (-dvdy) - wl[i+1] * aread3;
620: g[j+1][i] += (-dvdx) - wl[i] * aread3;
621: g[j+1][i+1] += ( dvdx + dvdy ) - wl[i+1] * aread3;
623: floc += area * (fl + fu);
625: }
626: }
628: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);
630: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
631: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
633: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
634: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
636: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
637: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
639: info = PetscLogFlops((xe-xs) * (ye-ys) * 67 + 1); CHKERRQ(info);
640: return 0;
641: } /* WholeJBearFunctionGradient */
646: /*
647: WholeJBearHessian - Evaluates Hessian over the whole grid
649: Input:
650: daapplication - TAO application object
651: da - distributed array
652: X - the current point, at which the function and gradient are evaluated
653: ptr - user-defined application context
655: Output:
656: H - Hessian at X
657: */
658: static int WholeJBearHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {
660: AppCtx *user = (AppCtx*)ptr;
661: int info;
662: PetscInt i, j, ind[4];
663: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
664: double smallH[4][4];
665: double wql, wqu;
666: double dydx, dxdy;
667: double *wq;
668: double wqldydx,wqldxdy,wqudydx,wqudxdy;
669: PetscTruth assembled;
671: wq = user->wq;
672: dydx = user->hy / user->hx;
673: dxdy = user->hx / user->hy;
675: info = MatAssembled(H,&assembled); CHKERRQ(info);
676: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
679: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
680: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
682: xe = gxs + gxm - 1;
683: ye = gys + gym - 1;
684: for (j = ys; j < ye; j++) {
685: for (i = xs; i < xe; i++) {
687: wql = (2.0 * wq[i] + wq[i+1]) / 6.0;
688: wqu = (wq[i] + 2.0 * wq[i+1]) / 6.0;
690: wqldydx = wql * dydx;
691: wqldxdy = wql * dxdy;
692: wqudydx = wqu * dydx;
693: wqudxdy = wqu * dxdy;
695: /* Hessian contribution at 0,0 */
696: smallH[0][0] = wqldxdy + wqldydx;
697: smallH[0][1] = smallH[1][0] = -wqldydx;
698: smallH[0][2] = smallH[2][0] = -wqldxdy;
699: smallH[0][3] = smallH[3][0] = 0.0;
701: /* Hessian contribution at 1,0 */
702: smallH[1][1] = (wqldydx + wqudxdy);
703: smallH[1][2] = smallH[2][1] = 0.0;
704: smallH[1][3] = smallH[3][1] = -wqudxdy;
706: /* Hessian contribution at 0,1 */
707: smallH[2][2] = (wqldxdy + wqudydx);
708: smallH[2][3] = smallH[3][2] = -wqudydx;
710: /* Hessian contribution at 1,1 */
711: smallH[3][3] = wqudydx + wqudxdy;
713: ind[0] = (j-gys) * gxm + (i-gxs);
714: ind[1] = ind[0] + 1;
715: ind[2] = ind[0] + gxm;
716: ind[3] = ind[2] + 1;
717: info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);
719: }
720: }
723: info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
724: info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
725: info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);
728: info = PetscLogFlops((xe-xs) * (ye-ys) * 14 + 2); CHKERRQ(info);
729: return 0;
731: } /* WholeJBearHessian */