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(int[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: int 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(int[2], double x[4], double *f, double g[4], void *ptr);
83: static int JBearLocalHessian(int[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, int, void *);
91: static int MyGridMonitorAfter1(TAO_APPLICATION, DA, int, void *);
95: int main( int argc, char **argv ) {
97: int info,iter; /* used to check for functions returning nonzeros */
98: int nlevels; /* multigrid levels */
99: int 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;
110: /* Initialize TAO */
111: PetscInitialize(&argc, &argv, (char *)0, help);
112: TaoInitialize(&argc, &argv, (char *)0, help);
114: PreLoadBegin(PreLoad,"Solve");
115:
116: info = AppCtxInitialize((void*)&user); CHKERRQ(info);
117:
118: nlevels=5;
119: info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
120: if (PreLoadIt == 0) {
121: nlevels = 1; user.mx = 11; user.my = 11; }
123: PetscPrintf(MPI_COMM_WORLD,"\n---- Journal Bearing Problem -----\n\n");
125: /* Let PETSc determine the vector distribution */
126: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
128: /* Create distributed array (DA) to manage parallel grid and vectors */
129: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
130: user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
131: for (iter=1;iter<nlevels;iter++){
132: info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
133: }
135: /* Create TAO solver and set desired solution method */
136: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
137: info = TaoApplicationCreate(PETSC_COMM_WORLD, &JBearApp); CHKERRQ(info);
138: info = TaoAppSetDAApp(JBearApp,DAarray,nlevels); CHKERRQ(info);
140: /* Sets routine bounds evaluation */
141: info = DAAppSetVariableBoundsRoutine(JBearApp,DASetBounds,(void *)&user); CHKERRQ(info);
143: info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
144: if (flg) {
146: /* Sets routines for function and gradient evaluation, element by element */
147: info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
148: if (flg) {
149: info = DAAppSetADElementFunctionGradient(JBearApp,ad_JBearLocalFunction,248,(void *)&user.fgctx); CHKERRQ(info);
150: } else {
151: info = DAAppSetElementObjectiveAndGradientRoutine(JBearApp,JBearLocalFunctionGradient,63,(void *)&user); CHKERRQ(info);
152: }
153: /* Sets routines for Hessian evaluation, element by element */
154: info = DAAppSetElementHessianRoutine(JBearApp,JBearLocalHessian,16,(void*)&user); CHKERRQ(info);
156: } else {
158: /* Sets routines for function and gradient evaluation, all in one routine */
159: info = DAAppSetObjectiveAndGradientRoutine(JBearApp,WholeJBearFunctionGradient,(void *)&user); CHKERRQ(info);
161: /* Sets routines for Hessian evaluation, all in one routine */
162: info = DAAppSetHessianRoutine(JBearApp,WholeJBearHessian,(void*)&user); CHKERRQ(info);
163:
164: }
166: info = DAAppSetBeforeMonitor(JBearApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
167: info = DAAppSetAfterMonitor(JBearApp,MyGridMonitorAfter1,(void*)&user); CHKERRQ(info);
168: info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
169: if (flg){
170: info = DAAppPrintStageTimes(JBearApp); CHKERRQ(info);
171: info = DAAppPrintInterpolationError(JBearApp); CHKERRQ(info);
172: }
173: /* Check for any tao command line options */
174: info = TaoAppSetRelativeTolerance(JBearApp,1.0e-8); CHKERRQ(info);
175: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
176: info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);
178: info = TaoSetOptions(JBearApp,tao); CHKERRQ(info);
180: info = DAAppGetSolution(JBearApp,0,&X); CHKERRQ(info);
181: info = FormInitialGuess(DAarray[0],X); CHKERRQ(info);
182: info = DAAppSetInitialSolution(JBearApp,X); CHKERRQ(info);
183: /* SOLVE THE APPLICATION */
184: info = TaoDAAppSolve(JBearApp,tao); CHKERRQ(info);
186: /* Get information on termination */
187: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
188: if (reason <= 0 ){
189: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
190: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
191: }
193: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
194: if (flg){
195: info = DAAppGetSolution(JBearApp,nlevels-1,&X); CHKERRQ(info);
196: info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
197: }
199: /* To View TAO solver information */
200: // info = TaoView(tao); CHKERRQ(info);
202: /* Free TAO data structures */
203: info = TaoDestroy(tao); CHKERRQ(info);
204: info = TaoAppDestroy(JBearApp); CHKERRQ(info);
206: /* Free PETSc data structures */
207: for (iter=0;iter<nlevels;iter++){
208: info = DADestroy(DAarray[iter]); CHKERRQ(info);
209: }
211: PreLoadEnd();
213: /* Finalize TAO */
214: TaoFinalize();
215: PetscFinalize();
217: return 0;
218: } /* main */
222: /*----- The following two routines
223:
224: MyGridMonitorBefore MyGridMonitorAfter
226: help diplay info of iterations at every grid level -------*/
230: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, int level, void *ctx) {
232: AppCtx *user = (AppCtx*)ctx;
233: int info,mx,my;
234: double t;
236: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
237: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
238: user->mx = mx;
239: user->my = my;
241: user->hx = (2.0 * 3.14159265358979) / (user->mx - 1);
242: user->hy = (2.0 * user->b) / (user->my - 1);
243: user->area = 0.5 * user->hx * user->hy;
245: user->wq = new double[user->mx];
246: user->wl = new double[user->mx];
247: for (int i=0; i<user->mx; i++) {
248: t = 1.0 + user->ecc * cos(i*user->hx);
249: user->wq[i] = t*t*t;
250: user->wl[i] = user->ecc * sin(i*user->hx);
251: }
252: info = PetscLogFlops(8 + 7 * user->mx); CHKERRQ(info);
254: user->fgctx.hx = user->hx;
255: user->fgctx.hy = user->hy;
256: user->fgctx.area = user->area;
257: user->fgctx.wq = user->wq;
258: user->fgctx.wl = user->wl;
260: PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my);
262: return 0;
263: }
267: static int MyGridMonitorAfter1(TAO_APPLICATION myapp, DA da, int level, void *ctx){
268:
269: AppCtx *user = (AppCtx*)ctx;
270: delete [] user->wq; delete [] user->wl;
271: return 0;
272: }
276: /*------- USER-DEFINED: initialize the application context information -------*/
280: /*
281: AppCtxInitialize - Sets initial values for the application context parameters
283: Input:
284: ptr - void user-defined application context
286: Output:
287: ptr - user-defined application context with the default or user-provided
288: parameters
289: */
290: static int AppCtxInitialize(void *ptr) {
292: AppCtx *user = (AppCtx*)ptr;
293: PetscTruth flg; /* flag for PETSc calls */
294: int info;
296: /* Specify default parameters */
297: user->ecc = 0.1;
298: user->b = 10.0;
299: user->mx = user->my = 11;
301: /* Check for command line arguments that override defaults */
302: info = PetscOptionsGetReal(TAO_NULL, "-ecc", &user->ecc, &flg); CHKERRQ(info);
303: info = PetscOptionsGetReal(TAO_NULL, "-b", &user->b, &flg); CHKERRQ(info);
304: info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
305: info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
307: return 0;
308: } /* AppCtxInitialize */
313: static int FormInitialGuess(DA da, Vec X)
314: {
315: int info, i, j, mx;
316: int xs, ys, xm, ym, xe, ye;
317: PetscReal hx, val;
318: double **x;
320: /* Get local mesh boundaries */
321: info = DAGetInfo(da,PETSC_NULL,&mx,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
322: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
323: hx = 2.0*4.0*atan(1.0)/(mx-1);
325: info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
326: xe = xs+xm; ye = ys+ym;
328: info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
329: /* Compute initial guess over locally owned part of mesh */
330: for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */
331: for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */
332: val = PetscMax(sin((i+1.0)*hx),0.0);
333: x[j][i] = val;
334: x[j][i] = 0;
335: }
336: }
337: info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);
339: return 0;
340: }
343: /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/
346: /*
347: FormBounds - Forms bounds on the variables
349: Input:
350: user - user-defined application context
352: Output:
353: XL - vector of lower bounds
354: XU - vector of upper bounds
356: */
357: static int DASetBounds(TAO_APPLICATION daapplication, DA da, Vec XL, Vec XU, void *ptr) {
359: AppCtx *user = (AppCtx*)ptr;
360: int i, j, info, mx, my;
361: int xs, xm, ys, ym;
362: double **xl, **xu;
364: mx = user->mx;
365: my = user->my;
367: info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
368: info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
369: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
371: for (j = ys; j < ys+ym; j++){
372: for (i = xs; i < xs+xm; i++){
373: xl[j][i] = 0.0;
374: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
375: xu[j][i] = 0.0;
376: } else {
377: xu[j][i] = TAO_INFINITY;
378: }
379: }
380: }
382: info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
383: info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);
385: return 0;
387: } /* DASetBounds */
393: /*
394: JBearLocalFunctionGradient - Evaluates function and gradient over the
395: local rectangular element
397: Input:
398: coor - vector with the indices of the position of current element
399: in the first, second and third directions
400: x - current point (values over the current rectangular element)
401: ptr - user-defined application context
403: Output:
404: f - value of the objective funtion at the local rectangular element
405: g - gradient of the local function
407: */
408: static int JBearLocalFunctionGradient(int coor[2], double x[4], double *f, double g[4], void *ptr) {
410: AppCtx *user = (AppCtx*)ptr;
412: double avgWq, sqGrad, avgWV, fl, fu;
413: double hx, hy, area, aread3, *wq, *wl;
414: double dvdx, dvdy;
415: int i;
417: hx = user->hx;
418: hy = user->hy;
419: area = user->area;
420: aread3 = area / 3.0;
421: wq = user->wq;
422: wl = user->wl;
423: i = coor[0];
425: /* lower triangle contribution */
426: dvdx = (x[0] - x[1]) / hx;
427: dvdy = (x[0] - x[2]) / hy;
428: sqGrad = dvdx * dvdx + dvdy * dvdy;
429: avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
430: avgWV = (wl[i]*x[0] + wl[i+1]*x[1] + wl[i]*x[2]) / 3.0;
431: fl = avgWq * sqGrad - avgWV;
433: dvdx = dvdx * hy * avgWq;
434: dvdy = dvdy * hx * avgWq;
435: g[0] = ( dvdx + dvdy ) - wl[i] * aread3;
436: g[1] = ( -dvdx ) - wl[i+1] * aread3;
437: g[2] = ( -dvdy ) - wl[i] * aread3;
439: /* upper triangle contribution */
440: dvdx = (x[3] - x[2]) / hx;
441: dvdy = (x[3] - x[1]) / hy;
442: sqGrad = dvdx * dvdx + dvdy * dvdy;
443: avgWq = (2.0 * wq[i+1] + wq[i]) / 6.0;
444: avgWV = (wl[i+1]*x[1] + wl[i]*x[2] + wl[i+1]*x[3]) / 3.0;
445: fu = avgWq * sqGrad - avgWV;
447: dvdx = dvdx * hy * avgWq;
448: dvdy = dvdy * hx * avgWq;
449: g[1] += (-dvdy) - wl[i+1] * aread3;
450: g[2] += (-dvdx) - wl[i] * aread3;
451: g[3] = ( dvdx + dvdy ) - wl[i+1] * aread3;
453: *f = area * (fl + fu);
455: return 0;
456: } /* JBearLocalFunctionGradient */
459: /*------- USER-DEFINED: routine to evaluate the Hessian
460: at a local (rectangular element) level -------*/
463: /*
464: JBearLocalHessian - Computes the Hessian of the local (partial) function
465: defined over the current rectangle
467: Input:
468: coor - vector with the indices of the position of current element
469: in the first, second and third directions
470: x - current local solution (over the rectangle only)
471: ptr - user-defined application context
473: Output:
474: H - Hessian matrix of the local function (wrt the four
475: points of the rectangle only)
477: */
478: static int JBearLocalHessian(int coor[2], double x[4], double H[4][4],void *ptr) {
480: AppCtx *user = (AppCtx*)ptr;
481: double wql, wqu;
482: double hx, hy, dydx, dxdy;
483: double *wq;
484: double wqldydx,wqldxdy,wqudydx,wqudxdy;
485: int i;
487: hx = user->hx;
488: hy = user->hy;
489: wq = user->wq;
490: i = coor[0];
492: dxdy = hx / hy;
493: dydx = hy / hx;
494: wql = (2.0 * wq[i] + wq[i+1]) / 6.0;
495: wqu = (wq[i] + 2.0 * wq[i+1]) / 6.0;
496: wqldydx = wql * dydx;
497: wqldxdy = wql * dxdy;
498: wqudydx = wqu * dydx;
499: wqudxdy = wqu * dxdy;
501: /* Hessian contribution at 0,0 */
502: H[0][0] = wqldxdy + wqldydx;
503: H[0][1] = H[1][0] = -wqldydx;
504: H[0][2] = H[2][0] = -wqldxdy;
505: H[0][3] = H[3][0] = 0.0;
507: /* Hessian contribution at 1,0 */
508: H[1][1] = wqldydx + wqudxdy;
509: H[1][2] = H[2][1] = 0.0;
510: H[1][3] = H[3][1] = -wqudxdy;
512: /* Hessian contribution at 0,1 */
513: H[2][2] = wqldxdy + wqudydx;
514: H[2][3] = H[3][2] = -wqudydx;
516: /* Hessian contribution at 1,1 */
517: H[3][3] = wqudydx + wqudxdy;
519: return 0;
521: } /* JBearLocalHessian */
524: /*------- USER-DEFINED: routine to evaluate the function
525: and gradient at the whole grid -------*/
529: /*
530: WholeJBearFunctionGradient - Evaluates function and gradient over the
531: whole grid
533: Input:
534: daapplication - TAO application object
535: da - distributed array
536: X - the current point, at which the function and gradient are evaluated
537: ptr - user-defined application context
539: Output:
540: f - value of the objective funtion at X
541: G - gradient at X
542: */
543: static int WholeJBearFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {
545: AppCtx *user = (AppCtx*)ptr;
546: Vec localX, localG;
547: int info, i, j;
548: int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
549: double **x, **g;
550: double floc = 0.0;
551: PetscScalar zero = 0.0;
553: double avgWq, sqGrad, avgWV, fl, fu;
554: double hx, hy, area, aread3, *wq, *wl;
555: double dvdx, dvdy;
557: hx = user->hx;
558: hy = user->hy;
559: area = user->area;
560: aread3= area/3.0;
561: wq = user->wq;
562: wl = user->wl;
564: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
565: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
566: info = VecSet(G, zero); CHKERRQ(info);
567: info = VecSet(localG, zero); CHKERRQ(info);
569: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
570: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
572: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
573: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
575: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
576: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
578: xe = gxs + gxm - 1;
579: ye = gys + gym - 1;
580: for (j = ys; j < ye; j++) {
581: for (i = xs; i < xe; i++) {
583: /* lower triangle contribution */
584: dvdx = (x[j][i] - x[j][i+1]) / hx;
585: dvdy = (x[j][i] - x[j+1][i]) / hy;
586: sqGrad = dvdx * dvdx + dvdy * dvdy;
587: avgWq = (2.0 * wq[i] + wq[i+1]) / 6.0;
588: avgWV = (wl[i]*x[j][i] + wl[i+1]*x[j][i+1] + wl[i]*x[j+1][i]) / 3.0;
589: fl = avgWq * sqGrad - avgWV;
591: dvdx = dvdx * hy * avgWq;
592: dvdy = dvdy * hx * avgWq;
593: g[j][i] += ( dvdx + dvdy ) - wl[i] * aread3;
594: g[j][i+1] += ( -dvdx ) - wl[i+1] * aread3;
595: g[j+1][i] += ( -dvdy ) - wl[i] * aread3;
597: /* upper triangle contribution */
598: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
599: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
600: sqGrad = dvdx * dvdx + dvdy * dvdy;
601: avgWq = (2.0 * wq[i+1] + wq[i]) / 6.0;
602: avgWV = (wl[i+1]*x[j][i+1] + wl[i]*x[j+1][i] + wl[i+1]*x[j+1][i+1]) / 3.0;
603: fu = avgWq * sqGrad - avgWV;
605: dvdx = dvdx * hy * avgWq;
606: dvdy = dvdy * hx * avgWq;
607: g[j][i+1] += (-dvdy) - wl[i+1] * aread3;
608: g[j+1][i] += (-dvdx) - wl[i] * aread3;
609: g[j+1][i+1] += ( dvdx + dvdy ) - wl[i+1] * aread3;
611: floc += area * (fl + fu);
613: }
614: }
616: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);
618: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
619: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
621: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
622: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
624: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
625: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
627: info = PetscLogFlops((xe-xs) * (ye-ys) * 67 + 1); CHKERRQ(info);
628: return 0;
629: } /* WholeJBearFunctionGradient */
634: /*
635: WholeJBearHessian - Evaluates Hessian over the whole grid
637: Input:
638: daapplication - TAO application object
639: da - distributed array
640: X - the current point, at which the function and gradient are evaluated
641: ptr - user-defined application context
643: Output:
644: H - Hessian at X
645: */
646: static int WholeJBearHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {
648: AppCtx *user = (AppCtx*)ptr;
649: int info, i, j, ind[4];
650: int xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
651: double smallH[4][4];
652: double wql, wqu;
653: double dydx, dxdy;
654: double *wq;
655: double wqldydx,wqldxdy,wqudydx,wqudxdy;
656: PetscTruth assembled;
658: wq = user->wq;
659: dydx = user->hy / user->hx;
660: dxdy = user->hx / user->hy;
662: info = MatAssembled(H,&assembled); CHKERRQ(info);
663: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
666: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
667: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
669: xe = gxs + gxm - 1;
670: ye = gys + gym - 1;
671: for (j = ys; j < ye; j++) {
672: for (i = xs; i < xe; i++) {
674: wql = (2.0 * wq[i] + wq[i+1]) / 6.0;
675: wqu = (wq[i] + 2.0 * wq[i+1]) / 6.0;
677: wqldydx = wql * dydx;
678: wqldxdy = wql * dxdy;
679: wqudydx = wqu * dydx;
680: wqudxdy = wqu * dxdy;
682: /* Hessian contribution at 0,0 */
683: smallH[0][0] = wqldxdy + wqldydx;
684: smallH[0][1] = smallH[1][0] = -wqldydx;
685: smallH[0][2] = smallH[2][0] = -wqldxdy;
686: smallH[0][3] = smallH[3][0] = 0.0;
688: /* Hessian contribution at 1,0 */
689: smallH[1][1] = (wqldydx + wqudxdy);
690: smallH[1][2] = smallH[2][1] = 0.0;
691: smallH[1][3] = smallH[3][1] = -wqudxdy;
693: /* Hessian contribution at 0,1 */
694: smallH[2][2] = (wqldxdy + wqudydx);
695: smallH[2][3] = smallH[3][2] = -wqudydx;
697: /* Hessian contribution at 1,1 */
698: smallH[3][3] = wqudydx + wqudxdy;
700: ind[0] = (j-gys) * gxm + (i-gxs);
701: ind[1] = ind[0] + 1;
702: ind[2] = ind[0] + gxm;
703: ind[3] = ind[2] + 1;
704: info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);
706: }
707: }
710: info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
711: info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
712: info = MatSetOption(H, MAT_SYMMETRIC); CHKERRQ(info);
715: info = PetscLogFlops((xe-xs) * (ye-ys) * 14 + 2); CHKERRQ(info);
716: return 0;
718: } /* WholeJBearHessian */