Actual source code: plate2.c
petsc-3.6.1 2015-08-06
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: TaoGetConvergedReason(); TaoDestroy();
30: Processors: n
31: T*/
34: /*
35: User-defined application context - contains data needed by the
36: application-provided call-back routines, FormFunctionGradient(),
37: FormHessian().
38: */
39: typedef struct {
40: /* problem parameters */
41: PetscReal bheight; /* Height of plate under the surface */
42: PetscInt mx, my; /* discretization in x, y directions */
43: PetscInt bmx,bmy; /* Size of plate under the surface */
44: Vec Bottom, Top, Left, Right; /* boundary values */
46: /* Working space */
47: Vec localX, localV; /* ghosted local vector */
48: DM dm; /* distributed array data structure */
49: Mat H;
50: } AppCtx;
52: /* -------- User-defined Routines --------- */
54: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
55: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
56: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
57: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
58: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
60: /* For testing matrix free submatrices */
61: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
62: PetscErrorCode MyMatMult(Mat,Vec,Vec);
66: int main( int argc, char **argv )
67: {
68: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
69: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
70: PetscInt m, N; /* number of local and global elements in vectors */
71: Vec x,xl,xu; /* solution vector and bounds*/
72: PetscBool flg; /* A return variable when checking for user options */
73: Tao tao; /* Tao solver context */
74: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
75: TaoConvergedReason reason;
76: Mat H_shell; /* to test matrix-free submatrices */
77: AppCtx user; /* user-defined work context */
79: /* Initialize PETSc, TAO */
80: PetscInitialize( &argc, &argv,(char *)0,help );
82: /* Specify default dimension of the problem */
83: user.mx = 10; user.my = 10; user.bheight=0.1;
85: /* Check for any command line arguments that override defaults */
86: PetscOptionsGetInt(NULL,"-mx",&user.mx,&flg);
87: PetscOptionsGetInt(NULL,"-my",&user.my,&flg);
88: PetscOptionsGetReal(NULL,"-bheight",&user.bheight,&flg);
90: user.bmx = user.mx/2; user.bmy = user.my/2;
91: PetscOptionsGetInt(NULL,"-bmx",&user.bmx,&flg);
92: PetscOptionsGetInt(NULL,"-bmy",&user.bmy,&flg);
94: PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
95: 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);
97: /* Calculate any derived values from parameters */
98: N = user.mx*user.my;
100: /* Let Petsc determine the dimensions of the local vectors */
101: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
103: /*
104: A two dimensional distributed array will help define this problem,
105: which derives from an elliptic PDE on two dimensional domain. From
106: the distributed array, Create the vectors.
107: */
108: DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
109: DMDA_STENCIL_BOX,user.mx,user.my,Nx,Ny,1,1,
110: NULL,NULL,&user.dm);
112: /*
113: Extract global and local vectors from DM; The local vectors are
114: used solely as work space for the evaluation of the function,
115: gradient, and Hessian. Duplicate for remaining vectors that are
116: the same types.
117: */
118: DMCreateGlobalVector(user.dm,&x); /* Solution */
119: DMCreateLocalVector(user.dm,&user.localX);
120: VecDuplicate(user.localX,&user.localV);
122: VecDuplicate(x,&xl);
123: VecDuplicate(x,&xu);
127: /* The TAO code begins here */
129: /*
130: Create TAO solver and set desired solution method
131: The method must either be TAOTRON or TAOBLMVM
132: If TAOBLMVM is used, then hessian function is not called.
133: */
134: TaoCreate(PETSC_COMM_WORLD,&tao);
135: TaoSetType(tao,TAOBLMVM);
137: /* Set initial solution guess; */
138: MSA_BoundaryConditions(&user);
139: MSA_InitialPoint(&user,x);
140: TaoSetInitialVector(tao,x);
142: /* Set routines for function, gradient and hessian evaluation */
143: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);
145: VecGetLocalSize(x,&m);
146: MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
147: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
149: DMGetLocalToGlobalMapping(user.dm,&isltog);
150: MatSetLocalToGlobalMapping(user.H,isltog,isltog);
151: PetscOptionsHasName(NULL,"-matrixfree",&flg);
152: if (flg) {
153: MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
154: MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
155: MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
156: TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
157: } else {
158: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
159: }
161: /* Set Variable bounds */
162: MSA_Plate(xl,xu,(void*)&user);
163: TaoSetVariableBounds(tao,xl,xu);
165: /* Check for any tao command line options */
166: TaoSetFromOptions(tao);
168: /* SOLVE THE APPLICATION */
169: TaoSolve(tao);
171: /* Get ierrrmation on converged */
172: TaoGetConvergedReason(tao,&reason);
173: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
174: if (reason <= 0){
175: PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
176: }
177: /* Free TAO data structures */
178: TaoDestroy(&tao);
180: /* Free PETSc data structures */
181: VecDestroy(&x);
182: VecDestroy(&xl);
183: VecDestroy(&xu);
184: MatDestroy(&user.H);
185: VecDestroy(&user.localX);
186: VecDestroy(&user.localV);
187: VecDestroy(&user.Bottom);
188: VecDestroy(&user.Top);
189: VecDestroy(&user.Left);
190: VecDestroy(&user.Right);
191: DMDestroy(&user.dm);
192: if (flg) {
193: MatDestroy(&H_shell);
194: }
195: PetscFinalize();
196: return 0;
197: }
201: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
203: Input Parameters:
204: . tao - the Tao context
205: . X - input vector
206: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
208: Output Parameters:
209: . fcn - the function value
210: . G - vector containing the newly evaluated gradient
212: Notes:
213: In this case, we discretize the domain and Create triangles. The
214: surface of each triangle is planar, whose surface area can be easily
215: computed. The total surface area is found by sweeping through the grid
216: and computing the surface area of the two triangles that have their
217: right angle at the grid point. The diagonal line segments on the
218: grid that define the triangles run from top left to lower right.
219: The numbering of points starts at the lower left and runs left to
220: right, then bottom to top.
221: */
222: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
223: {
224: AppCtx *user = (AppCtx *) userCtx;
226: PetscInt i,j,row;
227: PetscInt mx=user->mx, my=user->my;
228: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
229: PetscReal ft=0;
230: PetscReal zero=0.0;
231: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
232: PetscReal rhx=mx+1, rhy=my+1;
233: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
234: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
235: PetscReal *g, *x,*left,*right,*bottom,*top;
236: Vec localX = user->localX, localG = user->localV;
238: /* Get local mesh boundaries */
239: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
240: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
242: /* Scatter ghost points to local vector */
243: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
244: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
246: /* Initialize vector to zero */
247: VecSet(localG, zero);
249: /* Get pointers to vector data */
250: VecGetArray(localX,&x);
251: VecGetArray(localG,&g);
252: VecGetArray(user->Top,&top);
253: VecGetArray(user->Bottom,&bottom);
254: VecGetArray(user->Left,&left);
255: VecGetArray(user->Right,&right);
257: /* Compute function over the locally owned part of the mesh */
258: for (j=ys; j<ys+ym; j++){
259: for (i=xs; i< xs+xm; i++){
260: row=(j-gys)*gxm + (i-gxs);
262: xc = x[row];
263: xlt=xrb=xl=xr=xb=xt=xc;
265: if (i==0){ /* left side */
266: xl= left[j-ys+1];
267: xlt = left[j-ys+2];
268: } else {
269: xl = x[row-1];
270: }
272: if (j==0){ /* bottom side */
273: xb=bottom[i-xs+1];
274: xrb = bottom[i-xs+2];
275: } else {
276: xb = x[row-gxm];
277: }
279: if (i+1 == gxs+gxm){ /* right side */
280: xr=right[j-ys+1];
281: xrb = right[j-ys];
282: } else {
283: xr = x[row+1];
284: }
286: if (j+1==gys+gym){ /* top side */
287: xt=top[i-xs+1];
288: xlt = top[i-xs];
289: }else {
290: xt = x[row+gxm];
291: }
293: if (i>gxs && j+1<gys+gym){
294: xlt = x[row-1+gxm];
295: }
296: if (j>gys && i+1<gxs+gxm){
297: xrb = x[row+1-gxm];
298: }
300: d1 = (xc-xl);
301: d2 = (xc-xr);
302: d3 = (xc-xt);
303: d4 = (xc-xb);
304: d5 = (xr-xrb);
305: d6 = (xrb-xb);
306: d7 = (xlt-xl);
307: d8 = (xt-xlt);
309: df1dxc = d1*hydhx;
310: df2dxc = ( d1*hydhx + d4*hxdhy );
311: df3dxc = d3*hxdhy;
312: df4dxc = ( d2*hydhx + d3*hxdhy );
313: df5dxc = d2*hydhx;
314: df6dxc = d4*hxdhy;
316: d1 *= rhx;
317: d2 *= rhx;
318: d3 *= rhy;
319: d4 *= rhy;
320: d5 *= rhy;
321: d6 *= rhx;
322: d7 *= rhy;
323: d8 *= rhx;
325: f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
326: f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
327: f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
328: f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
329: f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
330: f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);
332: ft = ft + (f2 + f4);
334: df1dxc /= f1;
335: df2dxc /= f2;
336: df3dxc /= f3;
337: df4dxc /= f4;
338: df5dxc /= f5;
339: df6dxc /= f6;
341: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
343: }
344: }
347: /* Compute triangular areas along the border of the domain. */
348: if (xs==0){ /* left side */
349: for (j=ys; j<ys+ym; j++){
350: d3=(left[j-ys+1] - left[j-ys+2])*rhy;
351: d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
352: ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
353: }
354: }
355: if (ys==0){ /* bottom side */
356: for (i=xs; i<xs+xm; i++){
357: d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
358: d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
359: ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
360: }
361: }
363: if (xs+xm==mx){ /* right side */
364: for (j=ys; j< ys+ym; j++){
365: d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
366: d4=(right[j-ys]-right[j-ys+1])*rhy;
367: ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
368: }
369: }
370: if (ys+ym==my){ /* top side */
371: for (i=xs; i<xs+xm; i++){
372: d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
373: d4=(top[i-xs+1] - top[i-xs])*rhx;
374: ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
375: }
376: }
378: if (ys==0 && xs==0){
379: d1=(left[0]-left[1])*rhy;
380: d2=(bottom[0]-bottom[1])*rhx;
381: ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
382: }
383: if (ys+ym == my && xs+xm == mx){
384: d1=(right[ym+1] - right[ym])*rhy;
385: d2=(top[xm+1] - top[xm])*rhx;
386: ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
387: }
389: ft=ft*area;
390: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
393: /* Restore vectors */
394: VecRestoreArray(localX,&x);
395: VecRestoreArray(localG,&g);
396: VecRestoreArray(user->Left,&left);
397: VecRestoreArray(user->Top,&top);
398: VecRestoreArray(user->Bottom,&bottom);
399: VecRestoreArray(user->Right,&right);
401: /* Scatter values to global vector */
402: DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
403: DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);
405: PetscLogFlops(70*xm*ym);
407: return 0;
408: }
410: /* ------------------------------------------------------------------- */
413: /*
414: FormHessian - Evaluates Hessian matrix.
416: Input Parameters:
417: . tao - the Tao context
418: . x - input vector
419: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
421: Output Parameters:
422: . A - Hessian matrix
423: . B - optionally different preconditioning matrix
425: Notes:
426: Due to mesh point reordering with DMs, we must always work
427: with the local mesh points, and then transform them to the new
428: global numbering with the local-to-global mapping. We cannot work
429: directly with the global numbers for the original uniprocessor mesh!
431: Two methods are available for imposing this transformation
432: when setting matrix entries:
433: (A) MatSetValuesLocal(), using the local ordering (including
434: ghost points!)
435: - Do the following two steps once, before calling TaoSolve()
436: - Use DMGetISLocalToGlobalMapping() to extract the
437: local-to-global map from the DM
438: - Associate this map with the matrix by calling
439: MatSetLocalToGlobalMapping()
440: - Then set matrix entries using the local ordering
441: by calling MatSetValuesLocal()
442: (B) MatSetValues(), using the global ordering
443: - Use DMGetGlobalIndices() to extract the local-to-global map
444: - Then apply this map explicitly yourself
445: - Set matrix entries using the global ordering by calling
446: MatSetValues()
447: Option (A) seems cleaner/easier in many cases, and is the procedure
448: used in this example.
449: */
450: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
451: {
453: AppCtx *user = (AppCtx *) ptr;
454: PetscInt i,j,k,row;
455: PetscInt mx=user->mx, my=user->my;
456: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
457: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
458: PetscReal rhx=mx+1, rhy=my+1;
459: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
460: PetscReal hl,hr,ht,hb,hc,htl,hbr;
461: PetscReal *x,*left,*right,*bottom,*top;
462: PetscReal v[7];
463: Vec localX = user->localX;
464: PetscBool assembled;
467: /* Set various matrix options */
468: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
470: /* Initialize matrix entries to zero */
471: MatAssembled(Hessian,&assembled);
472: if (assembled){MatZeroEntries(Hessian);}
474: /* Get local mesh boundaries */
475: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
476: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
478: /* Scatter ghost points to local vector */
479: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
480: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
482: /* Get pointers to vector data */
483: VecGetArray(localX,&x);
484: VecGetArray(user->Top,&top);
485: VecGetArray(user->Bottom,&bottom);
486: VecGetArray(user->Left,&left);
487: VecGetArray(user->Right,&right);
489: /* Compute Hessian over the locally owned part of the mesh */
491: for (i=xs; i< xs+xm; i++){
493: for (j=ys; j<ys+ym; j++){
495: row=(j-gys)*gxm + (i-gxs);
497: xc = x[row];
498: xlt=xrb=xl=xr=xb=xt=xc;
500: /* Left side */
501: if (i==gxs){
502: xl= left[j-ys+1];
503: xlt = left[j-ys+2];
504: } else {
505: xl = x[row-1];
506: }
508: if (j==gys){
509: xb=bottom[i-xs+1];
510: xrb = bottom[i-xs+2];
511: } else {
512: xb = x[row-gxm];
513: }
515: if (i+1 == gxs+gxm){
516: xr=right[j-ys+1];
517: xrb = right[j-ys];
518: } else {
519: xr = x[row+1];
520: }
522: if (j+1==gys+gym){
523: xt=top[i-xs+1];
524: xlt = top[i-xs];
525: }else {
526: xt = x[row+gxm];
527: }
529: if (i>gxs && j+1<gys+gym){
530: xlt = x[row-1+gxm];
531: }
532: if (j>gys && i+1<gxs+gxm){
533: xrb = x[row+1-gxm];
534: }
537: d1 = (xc-xl)*rhx;
538: d2 = (xc-xr)*rhx;
539: d3 = (xc-xt)*rhy;
540: d4 = (xc-xb)*rhy;
541: d5 = (xrb-xr)*rhy;
542: d6 = (xrb-xb)*rhx;
543: d7 = (xlt-xl)*rhy;
544: d8 = (xlt-xt)*rhx;
546: f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
547: f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
548: f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
549: f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
550: f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
551: f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);
554: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
555: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
556: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
557: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
558: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
559: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
560: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
561: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
563: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
564: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
566: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
567: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
568: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
569: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
571: hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
573: k=0;
574: if (j>0){
575: v[k]=hb; col[k]=row - gxm; k++;
576: }
578: if (j>0 && i < mx -1){
579: v[k]=hbr; col[k]=row - gxm+1; k++;
580: }
582: if (i>0){
583: v[k]= hl; col[k]=row - 1; k++;
584: }
586: v[k]= hc; col[k]=row; k++;
588: if (i < mx-1 ){
589: v[k]= hr; col[k]=row+1; k++;
590: }
592: if (i>0 && j < my-1 ){
593: v[k]= htl; col[k] = row+gxm-1; k++;
594: }
596: if (j < my-1 ){
597: v[k]= ht; col[k] = row+gxm; k++;
598: }
600: /*
601: Set matrix values using local numbering, which was defined
602: earlier, in the main routine.
603: */
604: MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);
606: }
607: }
609: /* Restore vectors */
610: VecRestoreArray(localX,&x);
611: VecRestoreArray(user->Left,&left);
612: VecRestoreArray(user->Top,&top);
613: VecRestoreArray(user->Bottom,&bottom);
614: VecRestoreArray(user->Right,&right);
616: /* Assemble the matrix */
617: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
618: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
620: PetscLogFlops(199*xm*ym);
621: return 0;
622: }
624: /* ------------------------------------------------------------------- */
627: /*
628: MSA_BoundaryConditions - Calculates the boundary conditions for
629: the region.
631: Input Parameter:
632: . user - user-defined application context
634: Output Parameter:
635: . user - user-defined application context
636: */
637: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
638: {
639: int ierr;
640: PetscInt i,j,k,maxits=5,limit=0;
641: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
642: PetscInt mx=user->mx,my=user->my;
643: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
644: PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
645: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
646: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
647: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
648: PetscReal *boundary;
649: PetscBool flg;
650: Vec Bottom,Top,Right,Left;
652: /* Get local mesh boundaries */
653: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
654: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
657: bsize=xm+2;
658: lsize=ym+2;
659: rsize=ym+2;
660: tsize=xm+2;
662: VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
663: VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
664: VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
665: VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);
667: user->Top=Top;
668: user->Left=Left;
669: user->Bottom=Bottom;
670: user->Right=Right;
672: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
674: for (j=0; j<4; j++){
675: if (j==0){
676: yt=b;
677: xt=l+hx*xs;
678: limit=bsize;
679: VecGetArray(Bottom,&boundary);
680: } else if (j==1){
681: yt=t;
682: xt=l+hx*xs;
683: limit=tsize;
684: VecGetArray(Top,&boundary);
685: } else if (j==2){
686: yt=b+hy*ys;
687: xt=l;
688: limit=lsize;
689: VecGetArray(Left,&boundary);
690: } else if (j==3){
691: yt=b+hy*ys;
692: xt=r;
693: limit=rsize;
694: VecGetArray(Right,&boundary);
695: }
697: for (i=0; i<limit; i++){
698: u1=xt;
699: u2=-yt;
700: for (k=0; k<maxits; k++){
701: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
702: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
703: fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
704: if (fnorm <= tol) break;
705: njac11=one+u2*u2-u1*u1;
706: njac12=two*u1*u2;
707: njac21=-two*u1*u2;
708: njac22=-one - u1*u1 + u2*u2;
709: det = njac11*njac22-njac21*njac12;
710: u1 = u1-(njac22*nf1-njac12*nf2)/det;
711: u2 = u2-(njac11*nf2-njac21*nf1)/det;
712: }
714: boundary[i]=u1*u1-u2*u2;
715: if (j==0 || j==1) {
716: xt=xt+hx;
717: } else if (j==2 || j==3){
718: yt=yt+hy;
719: }
720: }
721: if (j==0){
722: VecRestoreArray(Bottom,&boundary);
723: } else if (j==1){
724: VecRestoreArray(Top,&boundary);
725: } else if (j==2){
726: VecRestoreArray(Left,&boundary);
727: } else if (j==3){
728: VecRestoreArray(Right,&boundary);
729: }
730: }
732: /* Scale the boundary if desired */
734: PetscOptionsGetReal(NULL,"-bottom",&scl,&flg);
735: if (flg){
736: VecScale(Bottom, scl);
737: }
738: PetscOptionsGetReal(NULL,"-top",&scl,&flg);
739: if (flg){
740: VecScale(Top, scl);
741: }
742: PetscOptionsGetReal(NULL,"-right",&scl,&flg);
743: if (flg){
744: VecScale(Right, scl);
745: }
747: PetscOptionsGetReal(NULL,"-left",&scl,&flg);
748: if (flg){
749: VecScale(Left, scl);
750: }
751: return 0;
752: }
755: /* ------------------------------------------------------------------- */
758: /*
759: MSA_Plate - Calculates an obstacle for surface to stretch over.
761: Input Parameter:
762: . user - user-defined application context
764: Output Parameter:
765: . user - user-defined application context
766: */
767: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx){
769: AppCtx *user=(AppCtx *)ctx;
771: PetscInt i,j,row;
772: PetscInt xs,ys,xm,ym;
773: PetscInt mx=user->mx, my=user->my, bmy, bmx;
774: PetscReal t1,t2,t3;
775: PetscReal *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
776: PetscBool cylinder;
778: user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
779: user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
780: bmy=user->bmy, bmx=user->bmx;
782: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
784: VecSet(XL, lb);
785: VecSet(XU, ub);
787: VecGetArray(XL,&xl);
789: PetscOptionsHasName(NULL,"-cylinder",&cylinder);
790: /* Compute the optional lower box */
791: if (cylinder){
792: for (i=xs; i< xs+xm; i++){
793: for (j=ys; j<ys+ym; j++){
794: row=(j-ys)*xm + (i-xs);
795: t1=(2.0*i-mx)*bmy;
796: t2=(2.0*j-my)*bmx;
797: t3=bmx*bmx*bmy*bmy;
798: if ( t1*t1 + t2*t2 <= t3 ){
799: xl[row] = user->bheight;
800: }
801: }
802: }
803: } else {
804: /* Compute the optional lower box */
805: for (i=xs; i< xs+xm; i++){
806: for (j=ys; j<ys+ym; j++){
807: row=(j-ys)*xm + (i-xs);
808: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
809: j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
810: xl[row] = user->bheight;
811: }
812: }
813: }
814: }
815: VecRestoreArray(XL,&xl);
817: return 0;
818: }
821: /* ------------------------------------------------------------------- */
824: /*
825: MSA_InitialPoint - Calculates the initial guess in one of three ways.
827: Input Parameters:
828: . user - user-defined application context
829: . X - vector for initial guess
831: Output Parameters:
832: . X - newly computed initial guess
833: */
834: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
835: {
837: PetscInt start=-1,i,j;
838: PetscReal zero=0.0;
839: PetscBool flg;
841: PetscOptionsGetInt(NULL,"-start",&start,&flg);
842: if (flg && start==0){ /* The zero vector is reasonable */
843: VecSet(X, zero);
844: } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */
845: PetscRandom rctx; PetscReal np5=-0.5;
847: PetscRandomCreate(MPI_COMM_WORLD,&rctx);
848: for (i=0; i<start; i++){
849: VecSetRandom(X, rctx);
850: }
851: PetscRandomDestroy(&rctx);
852: VecShift(X, np5);
854: } else { /* Take an average of the boundary conditions */
856: PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
857: PetscInt mx=user->mx,my=user->my;
858: PetscReal *x,*left,*right,*bottom,*top;
859: Vec localX = user->localX;
861: /* Get local mesh boundaries */
862: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
863: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
865: /* Get pointers to vector data */
866: VecGetArray(user->Top,&top);
867: VecGetArray(user->Bottom,&bottom);
868: VecGetArray(user->Left,&left);
869: VecGetArray(user->Right,&right);
871: VecGetArray(localX,&x);
872: /* Perform local computations */
873: for (j=ys; j<ys+ym; j++){
874: for (i=xs; i< xs+xm; i++){
875: row=(j-gys)*gxm + (i-gxs);
876: x[row] = ( (j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+
877: (i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
878: }
879: }
881: /* Restore vectors */
882: VecRestoreArray(localX,&x);
884: VecRestoreArray(user->Left,&left);
885: VecRestoreArray(user->Top,&top);
886: VecRestoreArray(user->Bottom,&bottom);
887: VecRestoreArray(user->Right,&right);
889: /* Scatter values into global vector */
890: DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
891: DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);
893: }
894: return 0;
895: }
897: /* For testing matrix free submatrices */
900: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
901: {
903: AppCtx *user = (AppCtx*)ptr;
905: FormHessian(tao,x,user->H,user->H,ptr);
906: return(0);
907: }
910: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
911: {
913: void *ptr;
914: AppCtx *user;
916: MatShellGetContext(H_shell,&ptr);
917: user = (AppCtx*)ptr;
918: MatMult(user->H,X,Y);
919: return(0);
920: }