Actual source code: plate2.c
1: #include <petscdmda.h>
2: #include <petsctao.h>
4: static char help[] =
5: "This example demonstrates use of the TAO package to \n\
6: solve an unconstrained minimization problem. This example is based on a \n\
7: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain, \n\
8: boundary values along the edges of the domain, and a plate represented by \n\
9: lower boundary conditions, the objective is to find the\n\
10: surface with the minimal area that satisfies the boundary conditions.\n\
11: The command line options are:\n\
12: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
13: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
14: -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
15: -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
16: -bheight <ht>, where <ht> = height of the plate\n\
17: -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
18: for an average of the boundary conditions\n\n";
20: /*
21: User-defined application context - contains data needed by the
22: application-provided call-back routines, FormFunctionGradient(),
23: FormHessian().
24: */
25: typedef struct {
26: /* problem parameters */
27: PetscReal bheight; /* Height of plate under the surface */
28: PetscInt mx, my; /* discretization in x, y directions */
29: PetscInt bmx,bmy; /* Size of plate under the surface */
30: Vec Bottom, Top, Left, Right; /* boundary values */
32: /* Working space */
33: Vec localX, localV; /* ghosted local vector */
34: DM dm; /* distributed array data structure */
35: Mat H;
36: } AppCtx;
38: /* -------- User-defined Routines --------- */
40: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
41: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
42: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
43: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
44: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
46: /* For testing matrix free submatrices */
47: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
48: PetscErrorCode MyMatMult(Mat,Vec,Vec);
50: int main(int argc, char **argv)
51: {
52: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
53: PetscInt m, N; /* number of local and global elements in vectors */
54: Vec x,xl,xu; /* solution vector and bounds*/
55: PetscBool flg; /* A return variable when checking for user options */
56: Tao tao; /* Tao solver context */
57: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
58: Mat H_shell; /* to test matrix-free submatrices */
59: AppCtx user; /* user-defined work context */
61: /* Initialize PETSc, TAO */
62: PetscInitialize(&argc, &argv,(char *)0,help);
64: /* Specify default dimension of the problem */
65: user.mx = 10; user.my = 10; user.bheight=0.1;
67: /* Check for any command line arguments that override defaults */
68: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
69: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
70: PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);
72: user.bmx = user.mx/2; user.bmy = user.my/2;
73: PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);
74: PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);
76: PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
77: PetscPrintf(PETSC_COMM_WORLD,"mx:%D, my:%D, bmx:%D, bmy:%D, height:%g\n",user.mx,user.my,user.bmx,user.bmy,(double)user.bheight);
79: /* Calculate any derived values from parameters */
80: N = user.mx*user.my;
82: /* Let Petsc determine the dimensions of the local vectors */
83: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
85: /*
86: A two dimensional distributed array will help define this problem,
87: which derives from an elliptic PDE on two dimensional domain. From
88: the distributed array, Create the vectors.
89: */
90: DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
91: DMSetFromOptions(user.dm);
92: DMSetUp(user.dm);
93: /*
94: Extract global and local vectors from DM; The local vectors are
95: used solely as work space for the evaluation of the function,
96: gradient, and Hessian. Duplicate for remaining vectors that are
97: the same types.
98: */
99: DMCreateGlobalVector(user.dm,&x); /* Solution */
100: DMCreateLocalVector(user.dm,&user.localX);
101: VecDuplicate(user.localX,&user.localV);
103: VecDuplicate(x,&xl);
104: VecDuplicate(x,&xu);
106: /* The TAO code begins here */
108: /*
109: Create TAO solver and set desired solution method
110: The method must either be TAOTRON or TAOBLMVM
111: If TAOBLMVM is used, then hessian function is not called.
112: */
113: TaoCreate(PETSC_COMM_WORLD,&tao);
114: TaoSetType(tao,TAOBLMVM);
116: /* Set initial solution guess; */
117: MSA_BoundaryConditions(&user);
118: MSA_InitialPoint(&user,x);
119: TaoSetSolution(tao,x);
121: /* Set routines for function, gradient and hessian evaluation */
122: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*) &user);
124: VecGetLocalSize(x,&m);
125: MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
126: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
128: DMGetLocalToGlobalMapping(user.dm,&isltog);
129: MatSetLocalToGlobalMapping(user.H,isltog,isltog);
130: PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);
131: if (flg) {
132: MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
133: MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
134: MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
135: TaoSetHessian(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
136: } else {
137: TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user);
138: }
140: /* Set Variable bounds */
141: MSA_Plate(xl,xu,(void*)&user);
142: TaoSetVariableBounds(tao,xl,xu);
144: /* Check for any tao command line options */
145: TaoSetFromOptions(tao);
147: /* SOLVE THE APPLICATION */
148: TaoSolve(tao);
150: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
152: /* Free TAO data structures */
153: TaoDestroy(&tao);
155: /* Free PETSc data structures */
156: VecDestroy(&x);
157: VecDestroy(&xl);
158: VecDestroy(&xu);
159: MatDestroy(&user.H);
160: VecDestroy(&user.localX);
161: VecDestroy(&user.localV);
162: VecDestroy(&user.Bottom);
163: VecDestroy(&user.Top);
164: VecDestroy(&user.Left);
165: VecDestroy(&user.Right);
166: DMDestroy(&user.dm);
167: if (flg) {
168: MatDestroy(&H_shell);
169: }
170: PetscFinalize();
171: return 0;
172: }
174: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
176: Input Parameters:
177: . tao - the Tao context
178: . X - input vector
179: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
181: Output Parameters:
182: . fcn - the function value
183: . G - vector containing the newly evaluated gradient
185: Notes:
186: In this case, we discretize the domain and Create triangles. The
187: surface of each triangle is planar, whose surface area can be easily
188: computed. The total surface area is found by sweeping through the grid
189: and computing the surface area of the two triangles that have their
190: right angle at the grid point. The diagonal line segments on the
191: grid that define the triangles run from top left to lower right.
192: The numbering of points starts at the lower left and runs left to
193: right, then bottom to top.
194: */
195: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
196: {
197: AppCtx *user = (AppCtx *) userCtx;
198: PetscInt i,j,row;
199: PetscInt mx=user->mx, my=user->my;
200: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
201: PetscReal ft=0;
202: PetscReal zero=0.0;
203: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
204: PetscReal rhx=mx+1, rhy=my+1;
205: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
206: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
207: PetscReal *g, *x,*left,*right,*bottom,*top;
208: Vec localX = user->localX, localG = user->localV;
210: /* Get local mesh boundaries */
211: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
212: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
214: /* Scatter ghost points to local vector */
215: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
216: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
218: /* Initialize vector to zero */
219: VecSet(localG, zero);
221: /* Get pointers to vector data */
222: VecGetArray(localX,&x);
223: VecGetArray(localG,&g);
224: VecGetArray(user->Top,&top);
225: VecGetArray(user->Bottom,&bottom);
226: VecGetArray(user->Left,&left);
227: VecGetArray(user->Right,&right);
229: /* Compute function over the locally owned part of the mesh */
230: for (j=ys; j<ys+ym; j++) {
231: for (i=xs; i< xs+xm; i++) {
232: row=(j-gys)*gxm + (i-gxs);
234: xc = x[row];
235: xlt=xrb=xl=xr=xb=xt=xc;
237: if (i==0) { /* left side */
238: xl= left[j-ys+1];
239: xlt = left[j-ys+2];
240: } else {
241: xl = x[row-1];
242: }
244: if (j==0) { /* bottom side */
245: xb=bottom[i-xs+1];
246: xrb = bottom[i-xs+2];
247: } else {
248: xb = x[row-gxm];
249: }
251: if (i+1 == gxs+gxm) { /* right side */
252: xr=right[j-ys+1];
253: xrb = right[j-ys];
254: } else {
255: xr = x[row+1];
256: }
258: if (j+1==gys+gym) { /* top side */
259: xt=top[i-xs+1];
260: xlt = top[i-xs];
261: }else {
262: xt = x[row+gxm];
263: }
265: if (i>gxs && j+1<gys+gym) {
266: xlt = x[row-1+gxm];
267: }
268: if (j>gys && i+1<gxs+gxm) {
269: xrb = x[row+1-gxm];
270: }
272: d1 = (xc-xl);
273: d2 = (xc-xr);
274: d3 = (xc-xt);
275: d4 = (xc-xb);
276: d5 = (xr-xrb);
277: d6 = (xrb-xb);
278: d7 = (xlt-xl);
279: d8 = (xt-xlt);
281: df1dxc = d1*hydhx;
282: df2dxc = (d1*hydhx + d4*hxdhy);
283: df3dxc = d3*hxdhy;
284: df4dxc = (d2*hydhx + d3*hxdhy);
285: df5dxc = d2*hydhx;
286: df6dxc = d4*hxdhy;
288: d1 *= rhx;
289: d2 *= rhx;
290: d3 *= rhy;
291: d4 *= rhy;
292: d5 *= rhy;
293: d6 *= rhx;
294: d7 *= rhy;
295: d8 *= rhx;
297: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
298: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
299: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
300: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
301: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
302: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
304: ft = ft + (f2 + f4);
306: df1dxc /= f1;
307: df2dxc /= f2;
308: df3dxc /= f3;
309: df4dxc /= f4;
310: df5dxc /= f5;
311: df6dxc /= f6;
313: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;
315: }
316: }
318: /* Compute triangular areas along the border of the domain. */
319: if (xs==0) { /* left side */
320: for (j=ys; j<ys+ym; j++) {
321: d3=(left[j-ys+1] - left[j-ys+2])*rhy;
322: d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
323: ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
324: }
325: }
326: if (ys==0) { /* bottom side */
327: for (i=xs; i<xs+xm; i++) {
328: d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
329: d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
330: ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
331: }
332: }
334: if (xs+xm==mx) { /* right side */
335: for (j=ys; j< ys+ym; j++) {
336: d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
337: d4=(right[j-ys]-right[j-ys+1])*rhy;
338: ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
339: }
340: }
341: if (ys+ym==my) { /* top side */
342: for (i=xs; i<xs+xm; i++) {
343: d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
344: d4=(top[i-xs+1] - top[i-xs])*rhx;
345: ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
346: }
347: }
349: if (ys==0 && xs==0) {
350: d1=(left[0]-left[1])*rhy;
351: d2=(bottom[0]-bottom[1])*rhx;
352: ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
353: }
354: if (ys+ym == my && xs+xm == mx) {
355: d1=(right[ym+1] - right[ym])*rhy;
356: d2=(top[xm+1] - top[xm])*rhx;
357: ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
358: }
360: ft=ft*area;
361: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
363: /* Restore vectors */
364: VecRestoreArray(localX,&x);
365: VecRestoreArray(localG,&g);
366: VecRestoreArray(user->Left,&left);
367: VecRestoreArray(user->Top,&top);
368: VecRestoreArray(user->Bottom,&bottom);
369: VecRestoreArray(user->Right,&right);
371: /* Scatter values to global vector */
372: DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
373: DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);
375: PetscLogFlops(70.0*xm*ym);
377: return 0;
378: }
380: /* ------------------------------------------------------------------- */
381: /*
382: FormHessian - Evaluates Hessian matrix.
384: Input Parameters:
385: . tao - the Tao context
386: . x - input vector
387: . ptr - optional user-defined context, as set by TaoSetHessian()
389: Output Parameters:
390: . A - Hessian matrix
391: . B - optionally different preconditioning matrix
393: Notes:
394: Due to mesh point reordering with DMs, we must always work
395: with the local mesh points, and then transform them to the new
396: global numbering with the local-to-global mapping. We cannot work
397: directly with the global numbers for the original uniprocessor mesh!
399: Two methods are available for imposing this transformation
400: when setting matrix entries:
401: (A) MatSetValuesLocal(), using the local ordering (including
402: ghost points!)
403: - Do the following two steps once, before calling TaoSolve()
404: - Use DMGetISLocalToGlobalMapping() to extract the
405: local-to-global map from the DM
406: - Associate this map with the matrix by calling
407: MatSetLocalToGlobalMapping()
408: - Then set matrix entries using the local ordering
409: by calling MatSetValuesLocal()
410: (B) MatSetValues(), using the global ordering
411: - Use DMGetGlobalIndices() to extract the local-to-global map
412: - Then apply this map explicitly yourself
413: - Set matrix entries using the global ordering by calling
414: MatSetValues()
415: Option (A) seems cleaner/easier in many cases, and is the procedure
416: used in this example.
417: */
418: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
419: {
420: AppCtx *user = (AppCtx *) ptr;
421: PetscInt i,j,k,row;
422: PetscInt mx=user->mx, my=user->my;
423: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
424: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
425: PetscReal rhx=mx+1, rhy=my+1;
426: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
427: PetscReal hl,hr,ht,hb,hc,htl,hbr;
428: PetscReal *x,*left,*right,*bottom,*top;
429: PetscReal v[7];
430: Vec localX = user->localX;
431: PetscBool assembled;
433: /* Set various matrix options */
434: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
436: /* Initialize matrix entries to zero */
437: MatAssembled(Hessian,&assembled);
438: if (assembled) MatZeroEntries(Hessian);
440: /* Get local mesh boundaries */
441: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
442: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
444: /* Scatter ghost points to local vector */
445: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
446: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
448: /* Get pointers to vector data */
449: VecGetArray(localX,&x);
450: VecGetArray(user->Top,&top);
451: VecGetArray(user->Bottom,&bottom);
452: VecGetArray(user->Left,&left);
453: VecGetArray(user->Right,&right);
455: /* Compute Hessian over the locally owned part of the mesh */
457: for (i=xs; i< xs+xm; i++) {
459: for (j=ys; j<ys+ym; j++) {
461: row=(j-gys)*gxm + (i-gxs);
463: xc = x[row];
464: xlt=xrb=xl=xr=xb=xt=xc;
466: /* Left side */
467: if (i==gxs) {
468: xl= left[j-ys+1];
469: xlt = left[j-ys+2];
470: } else {
471: xl = x[row-1];
472: }
474: if (j==gys) {
475: xb=bottom[i-xs+1];
476: xrb = bottom[i-xs+2];
477: } else {
478: xb = x[row-gxm];
479: }
481: if (i+1 == gxs+gxm) {
482: xr=right[j-ys+1];
483: xrb = right[j-ys];
484: } else {
485: xr = x[row+1];
486: }
488: if (j+1==gys+gym) {
489: xt=top[i-xs+1];
490: xlt = top[i-xs];
491: }else {
492: xt = x[row+gxm];
493: }
495: if (i>gxs && j+1<gys+gym) {
496: xlt = x[row-1+gxm];
497: }
498: if (j>gys && i+1<gxs+gxm) {
499: xrb = x[row+1-gxm];
500: }
502: d1 = (xc-xl)*rhx;
503: d2 = (xc-xr)*rhx;
504: d3 = (xc-xt)*rhy;
505: d4 = (xc-xb)*rhy;
506: d5 = (xrb-xr)*rhy;
507: d6 = (xrb-xb)*rhx;
508: d7 = (xlt-xl)*rhy;
509: d8 = (xlt-xt)*rhx;
511: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
512: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
513: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
514: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
515: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
516: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
518: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
519: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
520: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
521: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
522: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
523: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
524: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
525: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
527: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
528: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
530: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
531: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
532: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
533: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
535: hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
537: k=0;
538: if (j>0) {
539: v[k]=hb; col[k]=row - gxm; k++;
540: }
542: if (j>0 && i < mx -1) {
543: v[k]=hbr; col[k]=row - gxm+1; k++;
544: }
546: if (i>0) {
547: v[k]= hl; col[k]=row - 1; k++;
548: }
550: v[k]= hc; col[k]=row; k++;
552: if (i < mx-1) {
553: v[k]= hr; col[k]=row+1; k++;
554: }
556: if (i>0 && j < my-1) {
557: v[k]= htl; col[k] = row+gxm-1; k++;
558: }
560: if (j < my-1) {
561: v[k]= ht; col[k] = row+gxm; k++;
562: }
564: /*
565: Set matrix values using local numbering, which was defined
566: earlier, in the main routine.
567: */
568: MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);
570: }
571: }
573: /* Restore vectors */
574: VecRestoreArray(localX,&x);
575: VecRestoreArray(user->Left,&left);
576: VecRestoreArray(user->Top,&top);
577: VecRestoreArray(user->Bottom,&bottom);
578: VecRestoreArray(user->Right,&right);
580: /* Assemble the matrix */
581: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
582: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
584: PetscLogFlops(199.0*xm*ym);
585: return 0;
586: }
588: /* ------------------------------------------------------------------- */
589: /*
590: MSA_BoundaryConditions - Calculates the boundary conditions for
591: the region.
593: Input Parameter:
594: . user - user-defined application context
596: Output Parameter:
597: . user - user-defined application context
598: */
599: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
600: {
601: PetscInt i,j,k,maxits=5,limit=0;
602: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
603: PetscInt mx=user->mx,my=user->my;
604: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
605: PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
606: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
607: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
608: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
609: PetscReal *boundary;
610: PetscBool flg;
611: Vec Bottom,Top,Right,Left;
613: /* Get local mesh boundaries */
614: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
615: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
617: bsize=xm+2;
618: lsize=ym+2;
619: rsize=ym+2;
620: tsize=xm+2;
622: VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
623: VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
624: VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
625: VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);
627: user->Top=Top;
628: user->Left=Left;
629: user->Bottom=Bottom;
630: user->Right=Right;
632: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
634: for (j=0; j<4; j++) {
635: if (j==0) {
636: yt=b;
637: xt=l+hx*xs;
638: limit=bsize;
639: VecGetArray(Bottom,&boundary);
640: } else if (j==1) {
641: yt=t;
642: xt=l+hx*xs;
643: limit=tsize;
644: VecGetArray(Top,&boundary);
645: } else if (j==2) {
646: yt=b+hy*ys;
647: xt=l;
648: limit=lsize;
649: VecGetArray(Left,&boundary);
650: } else if (j==3) {
651: yt=b+hy*ys;
652: xt=r;
653: limit=rsize;
654: VecGetArray(Right,&boundary);
655: }
657: for (i=0; i<limit; i++) {
658: u1=xt;
659: u2=-yt;
660: for (k=0; k<maxits; k++) {
661: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
662: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
663: fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
664: if (fnorm <= tol) break;
665: njac11=one+u2*u2-u1*u1;
666: njac12=two*u1*u2;
667: njac21=-two*u1*u2;
668: njac22=-one - u1*u1 + u2*u2;
669: det = njac11*njac22-njac21*njac12;
670: u1 = u1-(njac22*nf1-njac12*nf2)/det;
671: u2 = u2-(njac11*nf2-njac21*nf1)/det;
672: }
674: boundary[i]=u1*u1-u2*u2;
675: if (j==0 || j==1) {
676: xt=xt+hx;
677: } else if (j==2 || j==3) {
678: yt=yt+hy;
679: }
680: }
681: if (j==0) {
682: VecRestoreArray(Bottom,&boundary);
683: } else if (j==1) {
684: VecRestoreArray(Top,&boundary);
685: } else if (j==2) {
686: VecRestoreArray(Left,&boundary);
687: } else if (j==3) {
688: VecRestoreArray(Right,&boundary);
689: }
690: }
692: /* Scale the boundary if desired */
694: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
695: if (flg) {
696: VecScale(Bottom, scl);
697: }
698: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
699: if (flg) {
700: VecScale(Top, scl);
701: }
702: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
703: if (flg) {
704: VecScale(Right, scl);
705: }
707: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
708: if (flg) {
709: VecScale(Left, scl);
710: }
711: return 0;
712: }
714: /* ------------------------------------------------------------------- */
715: /*
716: MSA_Plate - Calculates an obstacle for surface to stretch over.
718: Input Parameter:
719: . user - user-defined application context
721: Output Parameter:
722: . user - user-defined application context
723: */
724: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx)
725: {
726: AppCtx *user=(AppCtx *)ctx;
727: PetscInt i,j,row;
728: PetscInt xs,ys,xm,ym;
729: PetscInt mx=user->mx, my=user->my, bmy, bmx;
730: PetscReal t1,t2,t3;
731: PetscReal *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
732: PetscBool cylinder;
734: user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
735: user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
736: bmy=user->bmy; bmx=user->bmx;
738: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
740: VecSet(XL, lb);
741: VecSet(XU, ub);
743: VecGetArray(XL,&xl);
745: PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);
746: /* Compute the optional lower box */
747: if (cylinder) {
748: for (i=xs; i< xs+xm; i++) {
749: for (j=ys; j<ys+ym; j++) {
750: row=(j-ys)*xm + (i-xs);
751: t1=(2.0*i-mx)*bmy;
752: t2=(2.0*j-my)*bmx;
753: t3=bmx*bmx*bmy*bmy;
754: if (t1*t1 + t2*t2 <= t3) {
755: xl[row] = user->bheight;
756: }
757: }
758: }
759: } else {
760: /* Compute the optional lower box */
761: for (i=xs; i< xs+xm; i++) {
762: for (j=ys; j<ys+ym; j++) {
763: row=(j-ys)*xm + (i-xs);
764: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
765: j>=(my-bmy)/2 && j<my-(my-bmy)/2) {
766: xl[row] = user->bheight;
767: }
768: }
769: }
770: }
771: VecRestoreArray(XL,&xl);
773: return 0;
774: }
776: /* ------------------------------------------------------------------- */
777: /*
778: MSA_InitialPoint - Calculates the initial guess in one of three ways.
780: Input Parameters:
781: . user - user-defined application context
782: . X - vector for initial guess
784: Output Parameters:
785: . X - newly computed initial guess
786: */
787: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
788: {
789: PetscInt start=-1,i,j;
790: PetscReal zero=0.0;
791: PetscBool flg;
793: PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
794: if (flg && start==0) { /* The zero vector is reasonable */
795: VecSet(X, zero);
796: } else if (flg && start>0) { /* Try a random start between -0.5 and 0.5 */
797: PetscRandom rctx; PetscReal np5=-0.5;
799: PetscRandomCreate(MPI_COMM_WORLD,&rctx);
800: for (i=0; i<start; i++) {
801: VecSetRandom(X, rctx);
802: }
803: PetscRandomDestroy(&rctx);
804: VecShift(X, np5);
806: } else { /* Take an average of the boundary conditions */
808: PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
809: PetscInt mx=user->mx,my=user->my;
810: PetscReal *x,*left,*right,*bottom,*top;
811: Vec localX = user->localX;
813: /* Get local mesh boundaries */
814: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
815: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
817: /* Get pointers to vector data */
818: VecGetArray(user->Top,&top);
819: VecGetArray(user->Bottom,&bottom);
820: VecGetArray(user->Left,&left);
821: VecGetArray(user->Right,&right);
823: VecGetArray(localX,&x);
824: /* Perform local computations */
825: for (j=ys; j<ys+ym; j++) {
826: for (i=xs; i< xs+xm; i++) {
827: row=(j-gys)*gxm + (i-gxs);
828: x[row] = ((j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+(i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
829: }
830: }
832: /* Restore vectors */
833: VecRestoreArray(localX,&x);
835: VecRestoreArray(user->Left,&left);
836: VecRestoreArray(user->Top,&top);
837: VecRestoreArray(user->Bottom,&bottom);
838: VecRestoreArray(user->Right,&right);
840: /* Scatter values into global vector */
841: DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
842: DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);
844: }
845: return 0;
846: }
848: /* For testing matrix free submatrices */
849: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
850: {
851: AppCtx *user = (AppCtx*)ptr;
852: FormHessian(tao,x,user->H,user->H,ptr);
853: return 0;
854: }
855: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
856: {
857: void *ptr;
858: AppCtx *user;
859: MatShellGetContext(H_shell,&ptr);
860: user = (AppCtx*)ptr;
861: MatMult(user->H,X,Y);
862: return 0;
863: }
865: /*TEST
867: build:
868: requires: !complex
870: test:
871: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
872: requires: !single
874: test:
875: suffix: 2
876: nsize: 2
877: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
878: requires: !single
880: test:
881: suffix: 3
882: nsize: 3
883: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
884: requires: !single
886: test:
887: suffix: 4
888: nsize: 3
889: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
890: requires: !single
892: test:
893: suffix: 5
894: nsize: 3
895: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -pc_type none -tao_type tron -tao_gatol 1.e-5
896: requires: !single
898: test:
899: suffix: 6
900: nsize: 3
901: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
902: requires: !single
904: test:
905: suffix: 7
906: nsize: 3
907: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -pc_type none -tao_type gpcg -tao_gatol 1.e-5
908: requires: !single
910: test:
911: suffix: 8
912: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
913: requires: !single
915: test:
916: suffix: 9
917: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
918: requires: !single
920: test:
921: suffix: 10
922: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
923: requires: !single
925: test:
926: suffix: 11
927: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
928: requires: !single
930: test:
931: suffix: 12
932: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
933: requires: !single
935: test:
936: suffix: 13
937: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
938: requires: !single
940: test:
941: suffix: 14
942: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
943: requires: !single
945: test:
946: suffix: 15
947: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
948: requires: !single
950: test:
951: suffix: 16
952: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
953: requires: !single
955: test:
956: suffix: 17
957: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
958: requires: !single
960: test:
961: suffix: 18
962: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
963: requires: !single
965: test:
966: suffix: 19
967: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
968: requires: !single
970: test:
971: suffix: 20
972: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
974: TEST*/