Actual source code: minsurf3.c
1: /*$Id: minsurf3.c,v 1.1 2002/05/13 19:32:57 benson Exp $*/
3: /* Program usage: mpirun -np <proc> minsurf3 [-help] [all TAO options] */
5: /* minsurf1: boundary values like in COPS example:
6: - parabola for top and bottom boundaries
7: - zero o.w. */
9: /*
10: Include "tao.h" so we can use TAO solvers.
11: petscda.h for distributed array
12: ad_deriv.h for AD gradient
13: */
15: #include "petscda.h"
16: #include "tao.h"
17: #include "taodaapplication.h"
19: static char help[] =
20: "Given a rectangular 2-D domain and boundary values along \n\
21: the edges of the domain, the objective is to find the surface with \n\
22: the minimal area that satisfies the boundary conditions.\n\
23: The command line options are:\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\
29: -bottom <b>, where <b> = bottom bound on the rectangular domain\n\
30: -top <t>, where <t> = bottom bound on the rectangular domain\n\
31: -left <l>, where <l> = left bound on the rectangular domain\n\
32: -right <r>, where <r> = right bound on the rectangular domain\n\
33: -cops <r>, if COPS boundaries should be use (default MINPACK bounds)\n\
34: -obst <obst>, where <obst> = 1 is obstacle is present, zero otherwise \n\n";
36: /*T
37: Concepts: TAO - Solving a bounded minimization problem
38: Routines: TaoInitialize(); TaoFinalize();
39: Routines: TaoCreate(); TaoDestroy();
40: Routines: DAAppSetVariableBoundsRoutine();
41: Routines: DAAppSetElementObjectiveAndGradientRoutine();
42: Routines: DAAppSetElementHessianRoutine();
43: Routines: DAAppSetObjectiveAndGradientRoutine();
44: Routines: DAAppSetADElementFunctionGradient();
45: Routines: DAAppSetHessianRoutine();
46: Routines: TaoSetOptions();
47: Routines: TaoAppGetSolutionStatus(); TaoDAAppSolve();
48: Routines: DAAppSetBeforeMonitor(); TaoView();
49: Routines: DAAppGetSolution();
50: Routines: DAAppGetInterpolationMatrix();
51: Processors: n
52: T*/
54: /*
55: User-defined application context - contains data needed by the
56: application-provided call-back routines.
57: */
60: /*
61: This structure is used only when an ADIC generated gradient is used.
62: An InactiveDouble type is a double
63: */
64: typedef struct {
65:
66: InactiveDouble hx, hy; /* increment size in both directions */
67: InactiveDouble area; /* area of the triangles */
69: } ADFGCtx;
70: int ad_MSurfLocalFunction(PetscInt[2], DERIV_TYPE[4], DERIV_TYPE*, void*);
72: typedef struct {
73: PetscReal b, t, l, r; /* domain boundaries */
74: PetscReal bheight; /* height of obstacle */
75: PetscReal fx, fy; /* relative size of obstacle */
76: double hx, hy; /* increment size in both directions */
77: double area; /* area of the triangles */
78: PetscInt mx, my; /* discretization including boundaries */
79: PetscInt bmx, bmy; /* discretization of obstacle */
80: ADFGCtx fgctx; /* Used only when an ADIC generated gradient is used */
81: } AppCtx;
83: /* User-defined routines */
84: static int AppCtxInitialize(void *);
86: static int MSurfLocalFunctionGradient(PetscInt[2], double[4], double *, double[4], void *);
87: static int WholeMSurfFunctionGradient(TAO_APPLICATION,DA,Vec,double *,Vec,void*);
89: static int MSurfLocalHessian(PetscInt[2], double[4], double[4][4], void *);
90: static int WholeMSurfHessian(TAO_APPLICATION,DA,Vec,Mat,void*);
92: static int COPS_Bounds(TAO_APPLICATION, DA, Vec, Vec, void*);
93: static int MINPACK_Bounds(TAO_APPLICATION, DA, Vec, Vec, void *);
95: static int MyGridMonitorBefore(TAO_APPLICATION, DA, PetscInt, void *);
100: int main( int argc, char **argv ) {
102: int info; /* used to check for functions returning nonzeros */
103: PetscInt iter,nlevels;
104: PetscInt Nx,Ny;
105: double ff,gnorm;
106: DA DAarray[20];
107: Vec X;
108: PetscTruth flg, PreLoad = PETSC_TRUE; /* flags */
109: TaoMethod method = "tao_bnls"; /* minimization method */
110: TaoTerminateReason reason;
111: TAO_SOLVER tao; /* TAO_SOLVER solver context */
112: TAO_APPLICATION MSurfApp; /* The PETSc application */
113: AppCtx user; /* user-defined work context */
114: KSP ksp;
115: PC pc;
117: /* Initialize TAO */
118: PetscInitialize(&argc, &argv, (char *)0, help);
119: TaoInitialize(&argc, &argv, (char *)0, help);
121: PreLoadBegin(PreLoad,"Solve");
123: info = AppCtxInitialize((void*)&user); CHKERRQ(info);
125: nlevels=4;
126: info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);
127: if (PreLoadIt == 0) {
128: nlevels = 1; user.mx = 11; user.my = 11;
129: }
131: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimal Surface Problem (simple boundary) -----\n\n");
133: /* Let PETSc determine the vector distribution */
134: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
136: /* Create distributed array (DA) to manage parallel grid and vectors */
137: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
138: user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
139: for (iter=1;iter<nlevels;iter++){
140: info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
141: }
143: /* Create TAO solver and set desired solution method */
144: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
145: info = TaoApplicationCreate(PETSC_COMM_WORLD, &MSurfApp); CHKERRQ(info);
146: info = TaoAppSetDAApp(MSurfApp,DAarray,nlevels); CHKERRQ(info);
148: /* Sets routines for function, gradient and bounds evaluation */
149: info = PetscOptionsHasName(TAO_NULL, "-cops", &flg); CHKERRQ(info);
150: if (flg){
151: info = DAAppSetVariableBoundsRoutine(MSurfApp,COPS_Bounds,(void *)&user); CHKERRQ(info);
152: } else {
153: info = DAAppSetVariableBoundsRoutine(MSurfApp,MINPACK_Bounds,(void *)&user); CHKERRQ(info);
154: }
156: info = PetscOptionsHasName(TAO_NULL, "-byelement", &flg); CHKERRQ(info);
157: if (flg) {
159: /* Sets routines for function and gradient evaluation, element by element */
160: info = PetscOptionsHasName(TAO_NULL, "-adic", &flg); CHKERRQ(info);
161: if (flg) {
162: info = DAAppSetADElementFunctionGradient(MSurfApp,ad_MSurfLocalFunction,150,(void *)&user.fgctx); CHKERRQ(info);
163: } else {
164: info = DAAppSetElementObjectiveAndGradientRoutine(MSurfApp,MSurfLocalFunctionGradient,36,(void *)&user); CHKERRQ(info);
165: }
167: /* Sets routines for Hessian evaluation, element by element */
168: info = DAAppSetElementHessianRoutine(MSurfApp,MSurfLocalHessian,87,(void*)&user); CHKERRQ(info);
170: } else {
172: /* Sets routines for function and gradient evaluation, all in one routine */
173: info = DAAppSetObjectiveAndGradientRoutine(MSurfApp,WholeMSurfFunctionGradient,(void *)&user); CHKERRQ(info);
175: /* Sets routines for Hessian evaluation, all in one routine */
176: info = DAAppSetHessianRoutine(MSurfApp,WholeMSurfHessian,(void*)&user); CHKERRQ(info);
177:
178: }
180: info = DAAppSetBeforeMonitor(MSurfApp,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);
181: info = PetscOptionsHasName(TAO_NULL,"-tao_monitor", &flg); CHKERRQ(info);
182: if (flg){
183: info = DAAppPrintInterpolationError(MSurfApp); CHKERRQ(info);
184: info = DAAppPrintStageTimes(MSurfApp); CHKERRQ(info);
185: }
187: info = TaoAppSetRelativeTolerance(MSurfApp,1.0e-6); CHKERRQ(info);
188: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
189: info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);
191: info = TaoAppGetKSP(MSurfApp,&ksp);CHKERRQ(info);
192: info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
193: info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,100);CHKERRQ(info);
194: info = KSPGetPC(ksp,&pc);CHKERRQ(info);
195: info = PCSetType(pc,PCBJACOBI);CHKERRQ(info);
197: /* Check for any tao command line options */
198: info = TaoSetOptions(MSurfApp,tao); CHKERRQ(info);
200: /* SOLVE THE APPLICATION */
201: info = TaoDAAppSolve(MSurfApp,tao); CHKERRQ(info);
203: /* Get information on termination */
204: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
205: if (reason <= 0 ){
206: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
207: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
208: }
210: info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
211: if (flg){
212: info = DAAppGetSolution(MSurfApp,nlevels-1,&X); CHKERRQ(info);
213: info=VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
214: }
216: /* To View TAO solver information */
217: // info = TaoView(tao); CHKERRQ(info);
219: /* Free TAO data structures */
220: info = TaoDestroy(tao); CHKERRQ(info);
221: info = TaoAppDestroy(MSurfApp); CHKERRQ(info);
223: /* Free PETSc data structures */
224: for (iter=0;iter<nlevels;iter++){
225: info = DADestroy(DAarray[iter]); CHKERRQ(info);
226: }
228: PreLoadEnd();
230: /* Finalize TAO */
231: TaoFinalize();
232: PetscFinalize();
234: return 0;
235: } /* main */
239: /*----- The following two routines
240: MyGridMonitorBefore MyGridMonitorAfter
241: help diplay info of iterations at every grid level
242: */
246: static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, PetscInt level, void *ctx) {
248: AppCtx *user = (AppCtx*)ctx;
249: int info;
250: PetscInt mx,my;
252: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
253: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
254: user->mx = mx;
255: user->my = my;
256: user->hx = (user->r - user->l) / (user->mx - 1);
257: user->hy = (user->t - user->b) / (user->my - 1);
258: user->area = 0.5 * user->hx * user->hy;
259: user->fgctx.hx = user->hx;
260: user->fgctx.hy = user->hy;
261: user->fgctx.area = user->area;
263: user->bmx=(PetscInt)((mx+1)*user->fx); user->bmy=(PetscInt)((my+1)*user->fy);
265: PetscPrintf(MPI_COMM_WORLD,"Grid: %d, mx: %d my: %d \n",level,mx,my);
267: return 0;
268: }
272: /*
273: AppCtxInitialize - Sets initial values for the application context parameters
275: Input:
276: ptr - void user-defined application context
278: Output:
279: ptr - user-defined application context with the default or user-provided
280: parameters
281: */
282: static int AppCtxInitialize(void *ptr) {
284: AppCtx *user = (AppCtx*)ptr;
285: PetscTruth flg;
286: int info;
288: /* Specify default parameters */
289: user->mx = user->my = 11;
290: user->b = -0.5;
291: user->t = 0.5;
292: user->l = -0.5;
293: user->r = 0.5;
294: user->fx=0.5;
295: user->fy=0.5;
296: user->bheight=0.0;
298: /* Check for command line arguments that override defaults */
299: info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
300: info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
301: info = PetscOptionsGetReal(TAO_NULL, "-bottom", &user->b, &flg); CHKERRQ(info);
302: info = PetscOptionsGetReal(TAO_NULL, "-top", &user->t, &flg); CHKERRQ(info);
303: info = PetscOptionsGetReal(TAO_NULL, "-left", &user->l, &flg); CHKERRQ(info);
304: info = PetscOptionsGetReal(TAO_NULL, "-right", &user->r, &flg); CHKERRQ(info);
305: info = PetscOptionsGetReal(PETSC_NULL,"-bmx",&user->fx,&flg); CHKERRQ(info);
306: info = PetscOptionsGetReal(PETSC_NULL,"-bmy",&user->fy,&flg); CHKERRQ(info);
307: info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user->bheight,&flg); CHKERRQ(info);
309: user->hx = (user->r - user->l) / (user->mx - 1);
310: user->hy = (user->t - user->b) / (user->my - 1);
311: user->area = 0.5 * user->hx * user->hy;
312: info = PetscLogFlops(8); CHKERRQ(info);
314: return 0;
316: } /* AppCtxInitialize */
319: /*------- USER-DEFINED: set the upper and lower bounds for the variables -------*/
323: /*
324: FormBounds - Forms bounds on the variables
326: Input:
327: user - user-defined application context
329: Output:
330: XL - vector of lower bounds
331: XU - vector of upper bounds
333: */
334: static int COPS_Bounds(TAO_APPLICATION tao, DA da, Vec XL, Vec XU, void *ptr) {
336: AppCtx *user = (AppCtx*)ptr;
337: PetscTruth flg;
338: int info;
339: PetscInt i, j, mx, my, xs, xm, ys, ym;
340: double lb = -TAO_INFINITY;
341: double ub = TAO_INFINITY;
342: double **xl, **xu;
343: double xi, xi1, xi2;
344: double cx, cy, radx, rady, height = 1.0;
346: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
347: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
348: user->hx = (user->r - user->l) / (mx - 1);
349: user->hy = (user->t - user->b) / (my - 1);
350: user->area = 0.5 * user->hx * user->hy;
352: info = DAVecGetArray(da, XL, (void**)&xl); CHKERRQ(info);
353: info = DAVecGetArray(da, XU, (void**)&xu); CHKERRQ(info);
354: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
356: /* Compute default bounds */
357: for (j = ys; j < ys+ym; j++){
358: for (i = xs; i < xs+xm; i++){
360: if (j == 0 || j == my - 1) {
361: xi = user->l + i * user->hx;
362: xl[j][i] = xu[j][i] = -4 * (xi - user->l) * (xi - user->r);
363: } else if (i == 0 || i == mx - 1) {
364: xl[j][i] = xu[j][i] = 0.0;
365: } else {
366: xl[j][i] = lb;
367: xu[j][i] = ub;
368: }
370: }
371: }
373: /* Adjust lower bound if obstacle is present */
374: info = PetscOptionsHasName(PETSC_NULL, "-obst", &flg); CHKERRQ(info);
375: if (flg) {
376: radx = (user->r - user->l) * 0.25;
377: cx = user->l + 2.0 * radx;
378: rady = (user->t - user->b) * 0.25;
379: cy = user->b + 2.0 * rady;
380: for (j = ys; j < ys+ym; j++){
381: for (i = xs; i < xs+xm; i++){
382: xi1 = user->l + i * user->hx;
383: xi2 = user->b + j * user->hy;
384: if ( fabs(xi1 - cx) <= radx && fabs(xi2 - cy) <= rady ) {
385: xl[j][i] = height;
386: }
387: }
388: }
389: info = PetscLogFlops(8 + xm * ym * 6); CHKERRQ(info);
390: }
392: info = DAVecRestoreArray(da, XL, (void**)&xl); CHKERRQ(info);
393: info = DAVecRestoreArray(da, XU, (void**)&xu); CHKERRQ(info);
395: info = PetscLogFlops(12 * ym + 6); CHKERRQ(info);
397: return 0;
399: } /* DAGetBounds2d */
404: static int MINPACK_Bounds(TAO_APPLICATION tao, DA da, Vec XL,Vec XU, void *user1){
406: AppCtx *user=(AppCtx*)user1;
407: int info;
408: PetscInt i,j,k,limit=0,maxits=5;
409: PetscInt xs,ys,xm,ym;
410: PetscInt mx, my, bmy, bmx;
411: PetscInt row=0, bsize=0, lsize=0, tsize=0, rsize=0;
412: double bheight=user->bheight;
413: double one=1.0, two=2.0, three=3.0, tol=1e-10;
414: double fnorm,det,hx,hy,xt=0,yt=0;
415: double u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
416: double b=-0.5, t=0.5, l=-0.5, r=0.5;
417: PetscScalar scl,*xl,*xu, lb=TAO_NINFINITY, ub=TAO_INFINITY;
418: PetscTruth flg;
419: double **xll;
421: info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
422: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
424: bheight=user->bheight,lb=-1000; ub=1000;
425: bmx=user->bmx; bmy=user->bmy;
427: info = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
429: info = VecSet(XL, lb); CHKERRQ(info);
430: info = VecSet(XU, ub); CHKERRQ(info);
432: /*
433: Get pointers to vector data
434: */
435: info = DAVecGetArray(da,XL,(void**)&xll); CHKERRQ(info);
437: if (ys==0) bsize=xm;
438: if (xs==0) lsize=ym;
439: if (xs+xm==mx) rsize=ym;
440: if (ys+ym==my) tsize=xm;
441:
442: hx= (r-l)/(mx-1); hy=(t-b)/(my-1);
444: /* Compute the optional lower box */
445: for (i=xs; i< xs+xm; i++){
446: for (j=ys; j<ys+ym; j++){
447:
448: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 && j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
449: xll[j][i] = bheight;
450: }
452: }
453: }
455: info = DAVecRestoreArray(da,XL,(void**)&xll); CHKERRQ(info);
457: /* Boundary Values */
458: info = VecGetArray(XL,&xl); CHKERRQ(info);
459: info = VecGetArray(XU,&xu); CHKERRQ(info);
461: for (j=0; j<4; j++){
462: if (j==0){
463: yt=b;
464: xt=l+hx*xs;
465: limit=bsize;
466: scl=1.0;
467: info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
468: } else if (j==1){
469: yt=t;
470: xt=l+hx*xs;
471: limit=tsize;
472: scl=1.0;
473: info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
474: } else if (j==2){
475: yt=b+hy*ys;
476: xt=l;
477: limit=lsize;
478: scl=1.0;
479: info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
480: } else if (j==3){
481: yt=b+hy*ys;
482: xt=r;
483: limit=rsize;
484: scl=1.0;
485: info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
486: }
488: for (i=0; i<limit; i++){
489: u1=xt;
490: u2=-yt;
491: for (k=0; k<maxits; k++){
492: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
493: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
494: fnorm=sqrt(nf1*nf1+nf2*nf2);
495: if (fnorm <= tol) break;
496: njac11=one+u2*u2-u1*u1;
497: njac12=two*u1*u2;
498: njac21=-two*u1*u2;
499: njac22=-one - u1*u1 + u2*u2;
500: det = njac11*njac22-njac21*njac12;
501: u1 = u1-(njac22*nf1-njac12*nf2)/det;
502: u2 = u2-(njac11*nf2-njac21*nf1)/det;
503: }
505: if (j==0){
506: row = i;
507: } else if (j==1){
508: row = (ym-1)*xm+i;
509: } else if (j==2){
510: row = (xm*i);
511: } else if (j==3){
512: row = xm*(i+1)-1;
513: }
514:
515: xl[row]=(u1*u1-u2*u2)*scl;
516: xu[row]=(u1*u1-u2*u2)*scl;
518: if (j==0 || j==1) {
519: xt=xt+hx;
520: } else if (j==2 || j==3){
521: yt=yt+hy;
522: }
523:
524: }
525:
526: }
528: info = VecRestoreArray(XL,&xl); CHKERRQ(info);
529: info = VecRestoreArray(XU,&xu); CHKERRQ(info);
531: return 0;
532: }
534: /*------- USER-DEFINED: routine to evaluate the function and gradient
535: at a local (rectangular element) level -------*/
539: /*
540: MSurfLocalFunctionGradient - Evaluates function and gradient over the
541: local rectangular element
543: Input:
544: coor - vector with the indices of the position of current element
545: in the first, second and third directions
546: x - current point (values over the current rectangular element)
547: ptr - user-defined application context
549: Output:
550: f - value of the objective funtion at the local rectangular element
551: g - gradient of the local function
553: */
554: static int MSurfLocalFunctionGradient(PetscInt coor[2], double x[4], double *f, double g[4], void *ptr) {
556: AppCtx *user = (AppCtx*)ptr;
558: double hx, hy, area;
559: double dvdx, dvdy, flow, fup;
560: double d1,d2;
562: hx = user->hx;
563: hy = user->hy;
564: area = user->area;
566: /* lower triangle contribution */
567: dvdx = (x[0] - x[1]);
568: dvdy = (x[0] - x[2]);
569: g[1] = (dvdx * (hy/hx))/(-2);
570: g[2] = (dvdy * (hx/hy))/(-2);
571: dvdx /= hx;
572: dvdy /= hy;
573: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
574: *f=flow;
575: g[1] /= flow;
576: g[2] /= flow;
577: g[0] = -(g[1]+g[2]);
579: /* upper triangle contribution */
580: dvdx = (x[3] - x[2]);
581: dvdy = (x[3] - x[1]);
582: d1 = (dvdy*(hx/hy))/(-2);
583: d2 = (dvdx*(hy/hx))/(-2);
584: dvdx /= hx;
585: dvdy /= hy;
586: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
587: *f += fup;
588: g[1] += d1/fup;
589: g[2] += d2/fup;
590: g[3] = -(d1+d2)/fup;
592: *f *= area;
594: return 0;
595: } /* MSurfLocalFunctionGradient */
601: /*
602: MSurfLocalHessian - Computes the Hessian of the local (partial) function
603: defined over the current rectangle
605: Input:
606: coor - vector with the indices of the position of current element
607: in the first, second and third directions
608: x - current local solution (over the rectangle only)
609: ptr - user-defined application context
611: Output:
612: H - Hessian matrix of the local function (wrt the four
613: points of the rectangle only)
615: */
616: static int MSurfLocalHessian(PetscInt coor[2], double x[4], double H[4][4],void *ptr) {
618: AppCtx *user = (AppCtx*)ptr;
619: double hx, hy, area, byhxhx, byhyhy;
620: double dvdx, dvdy, flow, fup;
621: double areadivf, areadivf3;
623: hx = user->hx;
624: hy = user->hy;
625: area = user->area;
626:
627: byhxhx = 1.0 / (hx * hx);
628: byhyhy = 1.0 / (hy * hy);
630: /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
631: dvdx = (x[0] - x[1]) / hx; /* lower triangle contribution */
632: dvdy = (x[0] - x[2]) / hy;
633: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
634: dvdx = dvdx / hx;
635: dvdy = dvdy / hy;
636: areadivf = area / flow;
637: areadivf3 = areadivf / (flow * flow);
638: H[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
639: H[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
640: H[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
641: H[0][3] = 0.0;
642: H[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
643: H[1][2] = areadivf3 * (-dvdx) * dvdy;
644: H[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;
646: /* upper triangle contribution */
647: dvdx = (x[3] - x[2]) / hx;
648: dvdy = (x[3] - x[1]) / hy;
649: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
650: dvdx = dvdx / hx;
651: dvdy = dvdy / hy;
652: areadivf = area / fup;
653: areadivf3 = areadivf / (fup * fup);
654: H[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
655: H[1][2] += areadivf3 * (-dvdy) * dvdx;
656: H[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
657: H[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
658: H[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
659: H[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
661: H[1][0] = H[0][1];
662: H[2][0] = H[0][2];
663: H[3][0] = H[0][3];
664: H[2][1] = H[1][2];
665: H[3][1] = H[1][3];
666: H[3][2] = H[2][3];
668: return 0;
670: } /* MSurfLocalHessian */
673: /*------- USER-DEFINED: routine to evaluate the function
674: and gradient at the whole grid -------*/
678: /*
679: WholeMSurfFunctionGradient - Evaluates function and gradient over the
680: whole grid
682: Input:
683: daapplication - TAO application object
684: da - distributed array
685: X - the current point, at which the function and gradient are evaluated
686: ptr - user-defined application context
688: Output:
689: f - value of the objective funtion at X
690: G - gradient at X
691: */
692: static int WholeMSurfFunctionGradient(TAO_APPLICATION daapplication, DA da, Vec X, double *f, Vec G, void *ptr) {
694: AppCtx *user = (AppCtx*)ptr;
695: Vec localX, localG;
696: PetscInt i, j;
697: int info;
698: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
699: double **x, **g;
700: double floc = 0.0;
701: PetscScalar zero = 0.0;
703: double hx, hy, area;
704: double dvdx, dvdy, flow, fup;
705: double areadivf;
707: hx = user->hx;
708: hy = user->hy;
709: area = user->area;
711: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
712: info = DAGetLocalVector(da, &localG); CHKERRQ(info);
713: info = VecSet(G, zero); CHKERRQ(info);
714: info = VecSet(localG, zero); CHKERRQ(info);
716: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
717: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
719: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
720: info = DAVecGetArray(da, localG, (void**)&g); CHKERRQ(info);
722: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
723: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
725: xe = gxs + gxm - 1;
726: ye = gys + gym - 1;
727: for (j = ys; j < ye; j++) {
728: for (i = xs; i < xe; i++) {
730: /* lower triangle contribution */
731: dvdx = (x[j][i] - x[j][i+1]) / hx;
732: dvdy = (x[j][i] - x[j+1][i]) / hy;
733: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
734: areadivf = area / flow;
735: g[j][i] += (dvdx / hx + dvdy / hy) * areadivf;
736: g[j][i+1] += (-dvdx / hx) * areadivf;
737: g[j+1][i] += (-dvdy / hy) * areadivf;
739: /* upper triangle contribution */
740: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
741: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
742: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
743: areadivf = area / fup;
744: g[j][i+1] += (-dvdy / hy) * areadivf;
745: g[j+1][i] += (-dvdx / hx) * areadivf;
746: g[j+1][i+1] += (dvdx / hx + dvdy / hy) * areadivf;
748: floc += area * (flow + fup);
750: }
751: }
753: info = MPI_Allreduce(&floc, f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); CHKERRQ(info);
755: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
756: info = DAVecRestoreArray(da, localG, (void**)&g); CHKERRQ(info);
758: info = DALocalToGlobalBegin(da, localG, G); CHKERRQ(info);
759: info = DALocalToGlobalEnd(da, localG, G); CHKERRQ(info);
761: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
762: info = DARestoreLocalVector(da, &localG); CHKERRQ(info);
764: info = PetscLogFlops((xe-xs) * (ye-ys) * 42); CHKERRQ(info);
766: return 0;
767: } /* WholeMSurfFunctionGradient */
770: /*------- USER-DEFINED: routine to evaluate the Hessian
771: at the whole grid -------*/
775: /*
776: WholeMSurfHessian - Evaluates Hessian over the whole grid
778: Input:
779: daapplication - TAO application object
780: da - distributed array
781: X - the current point, at which the function and gradient are evaluated
782: ptr - user-defined application context
784: Output:
785: H - Hessian at X
786: */
787: static int WholeMSurfHessian(TAO_APPLICATION daapplication, DA da, Vec X, Mat H, void *ptr) {
789: AppCtx *user = (AppCtx*)ptr;
790: Vec localX;
791: int info;
792: PetscInt i, j, ind[4];
793: PetscInt xs, xm, gxs, gxm, xe, ys, ym, gys, gym, ye;
794: double smallH[4][4];
795: double **x;
797: double hx, hy, area, byhxhx, byhyhy;
798: double dvdx, dvdy, flow, fup;
799: double areadivf, areadivf3;
800: PetscTruth assembled;
802: hx = user->hx;
803: hy = user->hy;
804: area = user->area;
805:
806: byhxhx = 1.0 / (hx * hx);
807: byhyhy = 1.0 / (hy * hy);
809: info = DAGetLocalVector(da, &localX); CHKERRQ(info);
810: info = MatAssembled(H,&assembled); CHKERRQ(info);
811: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
813: info = DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX); CHKERRQ(info);
814: info = DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX); CHKERRQ(info);
816: info = DAVecGetArray(da, localX, (void**)&x); CHKERRQ(info);
818: info = DAGetCorners(da, &xs, &ys, TAO_NULL, &xm, &ym, TAO_NULL); CHKERRQ(info);
819: info = DAGetGhostCorners(da, &gxs, &gys, TAO_NULL, &gxm, &gym, TAO_NULL); CHKERRQ(info);
821: xe = gxs + gxm - 1;
822: ye = gys + gym - 1;
823: for (j = ys; j < ye; j++) {
824: for (i = xs; i < xe; i++) {
826: /* 0 is 0,0; 1 is 1,0; 2 is 0,1; 3 is 1,1 */
827: dvdx = (x[j][i] - x[j][i+1]) / hx; /* lower triangle contribution */
828: dvdy = (x[j][i] - x[j+1][i]) / hy;
829: flow = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
830: dvdx = dvdx / hx;
831: dvdy = dvdy / hy;
832: areadivf = area / flow;
833: areadivf3 = areadivf / (flow * flow);
834: smallH[0][0] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
835: smallH[0][1] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * (dvdx);
836: smallH[0][2] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * (dvdy);
837: smallH[0][3] = 0.0;
838: smallH[1][1] = areadivf * byhxhx - areadivf3 * dvdx * dvdx;
839: smallH[1][2] = areadivf3 * (-dvdx) * dvdy;
840: smallH[2][2] = areadivf * byhyhy - areadivf3 * dvdy * dvdy;
842: /* upper triangle contribution */
843: dvdx = (x[j+1][i+1] - x[j+1][i]) / hx;
844: dvdy = (x[j+1][i+1] - x[j][i+1]) / hy;
845: fup = sqrt( 1 + dvdx * dvdx + dvdy * dvdy );
846: dvdx = dvdx / hx;
847: dvdy = dvdy / hy;
848: areadivf = area / fup;
849: areadivf3 = areadivf / (fup * fup);
850: smallH[1][1] += areadivf * byhyhy - areadivf3 * dvdy * dvdy;
851: smallH[1][2] += areadivf3 * (-dvdy) * dvdx;
852: smallH[2][2] += areadivf * byhxhx - areadivf3 * (dvdx * dvdx);
853: smallH[1][3] = areadivf * (-byhyhy) + areadivf3 * (dvdx + dvdy) * dvdy;
854: smallH[2][3] = areadivf * (-byhxhx) + areadivf3 * (dvdx + dvdy) * dvdx;
855: smallH[3][3] = areadivf * (byhxhx + byhyhy) - areadivf3 * (dvdx + dvdy) * (dvdx + dvdy);
857: smallH[1][0] = smallH[0][1];
858: smallH[2][0] = smallH[0][2];
859: smallH[3][0] = smallH[0][3];
860: smallH[2][1] = smallH[1][2];
861: smallH[3][1] = smallH[1][3];
862: smallH[3][2] = smallH[2][3];
864: ind[0] = (j-gys) * gxm + (i-gxs);
865: ind[1] = ind[0] + 1;
866: ind[2] = ind[0] + gxm;
867: ind[3] = ind[2] + 1;
868: info = MatSetValuesLocal(H,4,ind,4,ind,(PetscScalar*)smallH,ADD_VALUES); CHKERRQ(info);
870: }
871: }
873: info = DAVecRestoreArray(da, localX, (void**)&x); CHKERRQ(info);
875: info = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
876: info = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); CHKERRQ(info);
877: info = MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(info);
879: info = DARestoreLocalVector(da, &localX); CHKERRQ(info);
881: info = PetscLogFlops((xe-xs) * (ye-ys) * 83 + 4); CHKERRQ(info);
882: return 0;
884: } /* WholeMSurfHessian */