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: /*T
21: Concepts: TAO^Solving a bound constrained minimization problem
22: Routines: TaoCreate();
23: Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
24: Routines: TaoSetHessianRoutine();
25: Routines: TaoSetInitialVector();
26: Routines: TaoSetVariableBounds();
27: Routines: TaoSetFromOptions();
28: Routines: TaoSolve(); TaoView();
29: Routines: TaoDestroy();
30: Processors: n
31: T*/
33: /*
34: User-defined application context - contains data needed by the
35: application-provided call-back routines, FormFunctionGradient(),
36: FormHessian().
37: */
38: typedef struct {
39: /* problem parameters */
40: PetscReal bheight; /* Height of plate under the surface */
41: PetscInt mx, my; /* discretization in x, y directions */
42: PetscInt bmx,bmy; /* Size of plate under the surface */
43: Vec Bottom, Top, Left, Right; /* boundary values */
45: /* Working space */
46: Vec localX, localV; /* ghosted local vector */
47: DM dm; /* distributed array data structure */
48: Mat H;
49: } AppCtx;
51: /* -------- User-defined Routines --------- */
53: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
54: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
55: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
56: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
57: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
59: /* For testing matrix free submatrices */
60: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
61: PetscErrorCode MyMatMult(Mat,Vec,Vec);
63: int main(int argc, char **argv)
64: {
65: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
66: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
67: PetscInt m, N; /* number of local and global elements in vectors */
68: Vec x,xl,xu; /* solution vector and bounds*/
69: PetscBool flg; /* A return variable when checking for user options */
70: Tao tao; /* Tao solver context */
71: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
72: Mat H_shell; /* to test matrix-free submatrices */
73: AppCtx user; /* user-defined work context */
75: /* Initialize PETSc, TAO */
76: PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;
78: /* Specify default dimension of the problem */
79: user.mx = 10; user.my = 10; user.bheight=0.1;
81: /* Check for any command line arguments that override defaults */
82: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
83: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
84: PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);
86: user.bmx = user.mx/2; user.bmy = user.my/2;
87: PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);
88: PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);
90: PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
91: 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);
93: /* Calculate any derived values from parameters */
94: N = user.mx*user.my;
96: /* Let Petsc determine the dimensions of the local vectors */
97: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
99: /*
100: A two dimensional distributed array will help define this problem,
101: which derives from an elliptic PDE on two dimensional domain. From
102: the distributed array, Create the vectors.
103: */
104: 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);
105: DMSetFromOptions(user.dm);
106: DMSetUp(user.dm);
107: /*
108: Extract global and local vectors from DM; The local vectors are
109: used solely as work space for the evaluation of the function,
110: gradient, and Hessian. Duplicate for remaining vectors that are
111: the same types.
112: */
113: DMCreateGlobalVector(user.dm,&x); /* Solution */
114: DMCreateLocalVector(user.dm,&user.localX);
115: VecDuplicate(user.localX,&user.localV);
117: VecDuplicate(x,&xl);
118: VecDuplicate(x,&xu);
120: /* The TAO code begins here */
122: /*
123: Create TAO solver and set desired solution method
124: The method must either be TAOTRON or TAOBLMVM
125: If TAOBLMVM is used, then hessian function is not called.
126: */
127: TaoCreate(PETSC_COMM_WORLD,&tao);
128: TaoSetType(tao,TAOBLMVM);
130: /* Set initial solution guess; */
131: MSA_BoundaryConditions(&user);
132: MSA_InitialPoint(&user,x);
133: TaoSetInitialVector(tao,x);
135: /* Set routines for function, gradient and hessian evaluation */
136: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);
138: VecGetLocalSize(x,&m);
139: MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
140: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
142: DMGetLocalToGlobalMapping(user.dm,&isltog);
143: MatSetLocalToGlobalMapping(user.H,isltog,isltog);
144: PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);
145: if (flg) {
146: MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
147: MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
148: MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
149: TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
150: } else {
151: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
152: }
154: /* Set Variable bounds */
155: MSA_Plate(xl,xu,(void*)&user);
156: TaoSetVariableBounds(tao,xl,xu);
158: /* Check for any tao command line options */
159: TaoSetFromOptions(tao);
161: /* SOLVE THE APPLICATION */
162: TaoSolve(tao);
164: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
166: /* Free TAO data structures */
167: TaoDestroy(&tao);
169: /* Free PETSc data structures */
170: VecDestroy(&x);
171: VecDestroy(&xl);
172: VecDestroy(&xu);
173: MatDestroy(&user.H);
174: VecDestroy(&user.localX);
175: VecDestroy(&user.localV);
176: VecDestroy(&user.Bottom);
177: VecDestroy(&user.Top);
178: VecDestroy(&user.Left);
179: VecDestroy(&user.Right);
180: DMDestroy(&user.dm);
181: if (flg) {
182: MatDestroy(&H_shell);
183: }
184: PetscFinalize();
185: return ierr;
186: }
188: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
190: Input Parameters:
191: . tao - the Tao context
192: . X - input vector
193: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
195: Output Parameters:
196: . fcn - the function value
197: . G - vector containing the newly evaluated gradient
199: Notes:
200: In this case, we discretize the domain and Create triangles. The
201: surface of each triangle is planar, whose surface area can be easily
202: computed. The total surface area is found by sweeping through the grid
203: and computing the surface area of the two triangles that have their
204: right angle at the grid point. The diagonal line segments on the
205: grid that define the triangles run from top left to lower right.
206: The numbering of points starts at the lower left and runs left to
207: right, then bottom to top.
208: */
209: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
210: {
211: AppCtx *user = (AppCtx *) userCtx;
213: PetscInt i,j,row;
214: PetscInt mx=user->mx, my=user->my;
215: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
216: PetscReal ft=0;
217: PetscReal zero=0.0;
218: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
219: PetscReal rhx=mx+1, rhy=my+1;
220: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
221: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
222: PetscReal *g, *x,*left,*right,*bottom,*top;
223: Vec localX = user->localX, localG = user->localV;
225: /* Get local mesh boundaries */
226: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
227: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
229: /* Scatter ghost points to local vector */
230: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
231: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
233: /* Initialize vector to zero */
234: VecSet(localG, zero);
236: /* Get pointers to vector data */
237: VecGetArray(localX,&x);
238: VecGetArray(localG,&g);
239: VecGetArray(user->Top,&top);
240: VecGetArray(user->Bottom,&bottom);
241: VecGetArray(user->Left,&left);
242: VecGetArray(user->Right,&right);
244: /* Compute function over the locally owned part of the mesh */
245: for (j=ys; j<ys+ym; j++) {
246: for (i=xs; i< xs+xm; i++) {
247: row=(j-gys)*gxm + (i-gxs);
249: xc = x[row];
250: xlt=xrb=xl=xr=xb=xt=xc;
252: if (i==0) { /* left side */
253: xl= left[j-ys+1];
254: xlt = left[j-ys+2];
255: } else {
256: xl = x[row-1];
257: }
259: if (j==0) { /* bottom side */
260: xb=bottom[i-xs+1];
261: xrb = bottom[i-xs+2];
262: } else {
263: xb = x[row-gxm];
264: }
266: if (i+1 == gxs+gxm) { /* right side */
267: xr=right[j-ys+1];
268: xrb = right[j-ys];
269: } else {
270: xr = x[row+1];
271: }
273: if (j+1==gys+gym) { /* top side */
274: xt=top[i-xs+1];
275: xlt = top[i-xs];
276: }else {
277: xt = x[row+gxm];
278: }
280: if (i>gxs && j+1<gys+gym) {
281: xlt = x[row-1+gxm];
282: }
283: if (j>gys && i+1<gxs+gxm) {
284: xrb = x[row+1-gxm];
285: }
287: d1 = (xc-xl);
288: d2 = (xc-xr);
289: d3 = (xc-xt);
290: d4 = (xc-xb);
291: d5 = (xr-xrb);
292: d6 = (xrb-xb);
293: d7 = (xlt-xl);
294: d8 = (xt-xlt);
296: df1dxc = d1*hydhx;
297: df2dxc = (d1*hydhx + d4*hxdhy);
298: df3dxc = d3*hxdhy;
299: df4dxc = (d2*hydhx + d3*hxdhy);
300: df5dxc = d2*hydhx;
301: df6dxc = d4*hxdhy;
303: d1 *= rhx;
304: d2 *= rhx;
305: d3 *= rhy;
306: d4 *= rhy;
307: d5 *= rhy;
308: d6 *= rhx;
309: d7 *= rhy;
310: d8 *= rhx;
312: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
313: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
314: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
315: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
316: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
317: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
319: ft = ft + (f2 + f4);
321: df1dxc /= f1;
322: df2dxc /= f2;
323: df3dxc /= f3;
324: df4dxc /= f4;
325: df5dxc /= f5;
326: df6dxc /= f6;
328: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;
330: }
331: }
333: /* Compute triangular areas along the border of the domain. */
334: if (xs==0) { /* left side */
335: for (j=ys; j<ys+ym; j++) {
336: d3=(left[j-ys+1] - left[j-ys+2])*rhy;
337: d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
338: ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
339: }
340: }
341: if (ys==0) { /* bottom side */
342: for (i=xs; i<xs+xm; i++) {
343: d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
344: d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
345: ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
346: }
347: }
349: if (xs+xm==mx) { /* right side */
350: for (j=ys; j< ys+ym; j++) {
351: d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
352: d4=(right[j-ys]-right[j-ys+1])*rhy;
353: ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
354: }
355: }
356: if (ys+ym==my) { /* top side */
357: for (i=xs; i<xs+xm; i++) {
358: d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
359: d4=(top[i-xs+1] - top[i-xs])*rhx;
360: ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
361: }
362: }
364: if (ys==0 && xs==0) {
365: d1=(left[0]-left[1])*rhy;
366: d2=(bottom[0]-bottom[1])*rhx;
367: ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
368: }
369: if (ys+ym == my && xs+xm == mx) {
370: d1=(right[ym+1] - right[ym])*rhy;
371: d2=(top[xm+1] - top[xm])*rhx;
372: ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
373: }
375: ft=ft*area;
376: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
378: /* Restore vectors */
379: VecRestoreArray(localX,&x);
380: VecRestoreArray(localG,&g);
381: VecRestoreArray(user->Left,&left);
382: VecRestoreArray(user->Top,&top);
383: VecRestoreArray(user->Bottom,&bottom);
384: VecRestoreArray(user->Right,&right);
386: /* Scatter values to global vector */
387: DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
388: DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);
390: PetscLogFlops(70.0*xm*ym);
392: return 0;
393: }
395: /* ------------------------------------------------------------------- */
396: /*
397: FormHessian - Evaluates Hessian matrix.
399: Input Parameters:
400: . tao - the Tao context
401: . x - input vector
402: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
404: Output Parameters:
405: . A - Hessian matrix
406: . B - optionally different preconditioning matrix
408: Notes:
409: Due to mesh point reordering with DMs, we must always work
410: with the local mesh points, and then transform them to the new
411: global numbering with the local-to-global mapping. We cannot work
412: directly with the global numbers for the original uniprocessor mesh!
414: Two methods are available for imposing this transformation
415: when setting matrix entries:
416: (A) MatSetValuesLocal(), using the local ordering (including
417: ghost points!)
418: - Do the following two steps once, before calling TaoSolve()
419: - Use DMGetISLocalToGlobalMapping() to extract the
420: local-to-global map from the DM
421: - Associate this map with the matrix by calling
422: MatSetLocalToGlobalMapping()
423: - Then set matrix entries using the local ordering
424: by calling MatSetValuesLocal()
425: (B) MatSetValues(), using the global ordering
426: - Use DMGetGlobalIndices() to extract the local-to-global map
427: - Then apply this map explicitly yourself
428: - Set matrix entries using the global ordering by calling
429: MatSetValues()
430: Option (A) seems cleaner/easier in many cases, and is the procedure
431: used in this example.
432: */
433: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
434: {
436: AppCtx *user = (AppCtx *) ptr;
437: PetscInt i,j,k,row;
438: PetscInt mx=user->mx, my=user->my;
439: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
440: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
441: PetscReal rhx=mx+1, rhy=my+1;
442: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
443: PetscReal hl,hr,ht,hb,hc,htl,hbr;
444: PetscReal *x,*left,*right,*bottom,*top;
445: PetscReal v[7];
446: Vec localX = user->localX;
447: PetscBool assembled;
449: /* Set various matrix options */
450: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
452: /* Initialize matrix entries to zero */
453: MatAssembled(Hessian,&assembled);
454: if (assembled) {MatZeroEntries(Hessian);}
456: /* Get local mesh boundaries */
457: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
458: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
460: /* Scatter ghost points to local vector */
461: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
462: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
464: /* Get pointers to vector data */
465: VecGetArray(localX,&x);
466: VecGetArray(user->Top,&top);
467: VecGetArray(user->Bottom,&bottom);
468: VecGetArray(user->Left,&left);
469: VecGetArray(user->Right,&right);
471: /* Compute Hessian over the locally owned part of the mesh */
473: for (i=xs; i< xs+xm; i++) {
475: for (j=ys; j<ys+ym; j++) {
477: row=(j-gys)*gxm + (i-gxs);
479: xc = x[row];
480: xlt=xrb=xl=xr=xb=xt=xc;
482: /* Left side */
483: if (i==gxs) {
484: xl= left[j-ys+1];
485: xlt = left[j-ys+2];
486: } else {
487: xl = x[row-1];
488: }
490: if (j==gys) {
491: xb=bottom[i-xs+1];
492: xrb = bottom[i-xs+2];
493: } else {
494: xb = x[row-gxm];
495: }
497: if (i+1 == gxs+gxm) {
498: xr=right[j-ys+1];
499: xrb = right[j-ys];
500: } else {
501: xr = x[row+1];
502: }
504: if (j+1==gys+gym) {
505: xt=top[i-xs+1];
506: xlt = top[i-xs];
507: }else {
508: xt = x[row+gxm];
509: }
511: if (i>gxs && j+1<gys+gym) {
512: xlt = x[row-1+gxm];
513: }
514: if (j>gys && i+1<gxs+gxm) {
515: xrb = x[row+1-gxm];
516: }
518: d1 = (xc-xl)*rhx;
519: d2 = (xc-xr)*rhx;
520: d3 = (xc-xt)*rhy;
521: d4 = (xc-xb)*rhy;
522: d5 = (xrb-xr)*rhy;
523: d6 = (xrb-xb)*rhx;
524: d7 = (xlt-xl)*rhy;
525: d8 = (xlt-xt)*rhx;
527: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
528: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
529: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
530: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
531: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
532: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
534: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
535: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
536: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
537: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
538: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
539: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
540: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
541: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
543: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
544: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
546: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
547: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
548: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
549: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
551: hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
553: k=0;
554: if (j>0) {
555: v[k]=hb; col[k]=row - gxm; k++;
556: }
558: if (j>0 && i < mx -1) {
559: v[k]=hbr; col[k]=row - gxm+1; k++;
560: }
562: if (i>0) {
563: v[k]= hl; col[k]=row - 1; k++;
564: }
566: v[k]= hc; col[k]=row; k++;
568: if (i < mx-1) {
569: v[k]= hr; col[k]=row+1; k++;
570: }
572: if (i>0 && j < my-1) {
573: v[k]= htl; col[k] = row+gxm-1; k++;
574: }
576: if (j < my-1) {
577: v[k]= ht; col[k] = row+gxm; k++;
578: }
580: /*
581: Set matrix values using local numbering, which was defined
582: earlier, in the main routine.
583: */
584: MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);
586: }
587: }
589: /* Restore vectors */
590: VecRestoreArray(localX,&x);
591: VecRestoreArray(user->Left,&left);
592: VecRestoreArray(user->Top,&top);
593: VecRestoreArray(user->Bottom,&bottom);
594: VecRestoreArray(user->Right,&right);
596: /* Assemble the matrix */
597: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
598: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
600: PetscLogFlops(199.0*xm*ym);
601: return 0;
602: }
604: /* ------------------------------------------------------------------- */
605: /*
606: MSA_BoundaryConditions - Calculates the boundary conditions for
607: the region.
609: Input Parameter:
610: . user - user-defined application context
612: Output Parameter:
613: . user - user-defined application context
614: */
615: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
616: {
617: int ierr;
618: PetscInt i,j,k,maxits=5,limit=0;
619: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
620: PetscInt mx=user->mx,my=user->my;
621: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
622: PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
623: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
624: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
625: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
626: PetscReal *boundary;
627: PetscBool flg;
628: Vec Bottom,Top,Right,Left;
630: /* Get local mesh boundaries */
631: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
632: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
634: bsize=xm+2;
635: lsize=ym+2;
636: rsize=ym+2;
637: tsize=xm+2;
639: VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
640: VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
641: VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
642: VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);
644: user->Top=Top;
645: user->Left=Left;
646: user->Bottom=Bottom;
647: user->Right=Right;
649: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
651: for (j=0; j<4; j++) {
652: if (j==0) {
653: yt=b;
654: xt=l+hx*xs;
655: limit=bsize;
656: VecGetArray(Bottom,&boundary);
657: } else if (j==1) {
658: yt=t;
659: xt=l+hx*xs;
660: limit=tsize;
661: VecGetArray(Top,&boundary);
662: } else if (j==2) {
663: yt=b+hy*ys;
664: xt=l;
665: limit=lsize;
666: VecGetArray(Left,&boundary);
667: } else if (j==3) {
668: yt=b+hy*ys;
669: xt=r;
670: limit=rsize;
671: VecGetArray(Right,&boundary);
672: }
674: for (i=0; i<limit; i++) {
675: u1=xt;
676: u2=-yt;
677: for (k=0; k<maxits; k++) {
678: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
679: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
680: fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
681: if (fnorm <= tol) break;
682: njac11=one+u2*u2-u1*u1;
683: njac12=two*u1*u2;
684: njac21=-two*u1*u2;
685: njac22=-one - u1*u1 + u2*u2;
686: det = njac11*njac22-njac21*njac12;
687: u1 = u1-(njac22*nf1-njac12*nf2)/det;
688: u2 = u2-(njac11*nf2-njac21*nf1)/det;
689: }
691: boundary[i]=u1*u1-u2*u2;
692: if (j==0 || j==1) {
693: xt=xt+hx;
694: } else if (j==2 || j==3) {
695: yt=yt+hy;
696: }
697: }
698: if (j==0) {
699: VecRestoreArray(Bottom,&boundary);
700: } else if (j==1) {
701: VecRestoreArray(Top,&boundary);
702: } else if (j==2) {
703: VecRestoreArray(Left,&boundary);
704: } else if (j==3) {
705: VecRestoreArray(Right,&boundary);
706: }
707: }
709: /* Scale the boundary if desired */
711: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
712: if (flg) {
713: VecScale(Bottom, scl);
714: }
715: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
716: if (flg) {
717: VecScale(Top, scl);
718: }
719: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
720: if (flg) {
721: VecScale(Right, scl);
722: }
724: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
725: if (flg) {
726: VecScale(Left, scl);
727: }
728: return 0;
729: }
731: /* ------------------------------------------------------------------- */
732: /*
733: MSA_Plate - Calculates an obstacle for surface to stretch over.
735: Input Parameter:
736: . user - user-defined application context
738: Output Parameter:
739: . user - user-defined application context
740: */
741: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx)
742: {
743: AppCtx *user=(AppCtx *)ctx;
745: PetscInt i,j,row;
746: PetscInt xs,ys,xm,ym;
747: PetscInt mx=user->mx, my=user->my, bmy, bmx;
748: PetscReal t1,t2,t3;
749: PetscReal *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
750: PetscBool cylinder;
752: user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
753: user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
754: bmy=user->bmy; bmx=user->bmx;
756: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
758: VecSet(XL, lb);
759: VecSet(XU, ub);
761: VecGetArray(XL,&xl);
763: PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);
764: /* Compute the optional lower box */
765: if (cylinder) {
766: for (i=xs; i< xs+xm; i++) {
767: for (j=ys; j<ys+ym; j++) {
768: row=(j-ys)*xm + (i-xs);
769: t1=(2.0*i-mx)*bmy;
770: t2=(2.0*j-my)*bmx;
771: t3=bmx*bmx*bmy*bmy;
772: if (t1*t1 + t2*t2 <= t3) {
773: xl[row] = user->bheight;
774: }
775: }
776: }
777: } else {
778: /* Compute the optional lower box */
779: for (i=xs; i< xs+xm; i++) {
780: for (j=ys; j<ys+ym; j++) {
781: row=(j-ys)*xm + (i-xs);
782: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
783: j>=(my-bmy)/2 && j<my-(my-bmy)/2) {
784: xl[row] = user->bheight;
785: }
786: }
787: }
788: }
789: VecRestoreArray(XL,&xl);
791: return 0;
792: }
794: /* ------------------------------------------------------------------- */
795: /*
796: MSA_InitialPoint - Calculates the initial guess in one of three ways.
798: Input Parameters:
799: . user - user-defined application context
800: . X - vector for initial guess
802: Output Parameters:
803: . X - newly computed initial guess
804: */
805: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
806: {
808: PetscInt start=-1,i,j;
809: PetscReal zero=0.0;
810: PetscBool flg;
812: PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
813: if (flg && start==0) { /* The zero vector is reasonable */
814: VecSet(X, zero);
815: } else if (flg && start>0) { /* Try a random start between -0.5 and 0.5 */
816: PetscRandom rctx; PetscReal np5=-0.5;
818: PetscRandomCreate(MPI_COMM_WORLD,&rctx);
819: for (i=0; i<start; i++) {
820: VecSetRandom(X, rctx);
821: }
822: PetscRandomDestroy(&rctx);
823: VecShift(X, np5);
825: } else { /* Take an average of the boundary conditions */
827: PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
828: PetscInt mx=user->mx,my=user->my;
829: PetscReal *x,*left,*right,*bottom,*top;
830: Vec localX = user->localX;
832: /* Get local mesh boundaries */
833: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
834: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
836: /* Get pointers to vector data */
837: VecGetArray(user->Top,&top);
838: VecGetArray(user->Bottom,&bottom);
839: VecGetArray(user->Left,&left);
840: VecGetArray(user->Right,&right);
842: VecGetArray(localX,&x);
843: /* Perform local computations */
844: for (j=ys; j<ys+ym; j++) {
845: for (i=xs; i< xs+xm; i++) {
846: row=(j-gys)*gxm + (i-gxs);
847: 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;
848: }
849: }
851: /* Restore vectors */
852: VecRestoreArray(localX,&x);
854: VecRestoreArray(user->Left,&left);
855: VecRestoreArray(user->Top,&top);
856: VecRestoreArray(user->Bottom,&bottom);
857: VecRestoreArray(user->Right,&right);
859: /* Scatter values into global vector */
860: DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
861: DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);
863: }
864: return 0;
865: }
867: /* For testing matrix free submatrices */
868: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
869: {
871: AppCtx *user = (AppCtx*)ptr;
873: FormHessian(tao,x,user->H,user->H,ptr);
874: return(0);
875: }
876: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
877: {
879: void *ptr;
880: AppCtx *user;
882: MatShellGetContext(H_shell,&ptr);
883: user = (AppCtx*)ptr;
884: MatMult(user->H,X,Y);
885: return(0);
886: }
888: /*TEST
890: build:
891: requires: !complex
893: test:
894: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
895: requires: !single
897: test:
898: suffix: 2
899: nsize: 2
900: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
901: requires: !single
903: test:
904: suffix: 3
905: nsize: 3
906: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
907: requires: !single
909: test:
910: suffix: 4
911: nsize: 3
912: 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
913: requires: !single
915: test:
916: suffix: 5
917: nsize: 3
918: 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
919: requires: !single
921: test:
922: suffix: 6
923: nsize: 3
924: 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
925: requires: !single
927: test:
928: suffix: 7
929: nsize: 3
930: 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
931: requires: !single
933: test:
934: suffix: 8
935: 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
936: requires: !single
938: test:
939: suffix: 9
940: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
941: requires: !single
943: test:
944: suffix: 10
945: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
946: requires: !single
948: test:
949: suffix: 11
950: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
951: requires: !single
953: test:
954: suffix: 12
955: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
956: requires: !single
958: test:
959: suffix: 13
960: 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
961: requires: !single
963: test:
964: suffix: 14
965: 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
966: requires: !single
968: test:
969: suffix: 15
970: 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
971: requires: !single
973: test:
974: suffix: 16
975: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
976: requires: !single
978: test:
979: suffix: 17
980: 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
981: requires: !single
983: test:
984: suffix: 18
985: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
986: requires: !single
988: test:
989: suffix: 19
990: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
991: requires: !single
993: test:
994: suffix: 20
995: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
997: TEST*/