Actual source code: plate2.c
petsc-3.8.4 2018-03-24
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*/
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);
64: int main( int argc, char **argv )
65: {
66: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
67: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
68: PetscInt m, N; /* number of local and global elements in vectors */
69: Vec x,xl,xu; /* solution vector and bounds*/
70: PetscBool flg; /* A return variable when checking for user options */
71: Tao tao; /* Tao solver context */
72: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
73: Mat H_shell; /* to test matrix-free submatrices */
74: AppCtx user; /* user-defined work context */
76: /* Initialize PETSc, TAO */
77: PetscInitialize( &argc, &argv,(char *)0,help );if (ierr) return ierr;
79: /* Specify default dimension of the problem */
80: user.mx = 10; user.my = 10; user.bheight=0.1;
82: /* Check for any command line arguments that override defaults */
83: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
84: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
85: PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);
87: user.bmx = user.mx/2; user.bmy = user.my/2;
88: PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);
89: PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);
91: PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
92: 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);
94: /* Calculate any derived values from parameters */
95: N = user.mx*user.my;
97: /* Let Petsc determine the dimensions of the local vectors */
98: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
100: /*
101: A two dimensional distributed array will help define this problem,
102: which derives from an elliptic PDE on two dimensional domain. From
103: the distributed array, Create the vectors.
104: */
105: 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);
106: DMSetFromOptions(user.dm);
107: DMSetUp(user.dm);
108: /*
109: Extract global and local vectors from DM; The local vectors are
110: used solely as work space for the evaluation of the function,
111: gradient, and Hessian. Duplicate for remaining vectors that are
112: the same types.
113: */
114: DMCreateGlobalVector(user.dm,&x); /* Solution */
115: DMCreateLocalVector(user.dm,&user.localX);
116: VecDuplicate(user.localX,&user.localV);
118: VecDuplicate(x,&xl);
119: VecDuplicate(x,&xu);
121: /* The TAO code begins here */
123: /*
124: Create TAO solver and set desired solution method
125: The method must either be TAOTRON or TAOBLMVM
126: If TAOBLMVM is used, then hessian function is not called.
127: */
128: TaoCreate(PETSC_COMM_WORLD,&tao);
129: TaoSetType(tao,TAOBLMVM);
131: /* Set initial solution guess; */
132: MSA_BoundaryConditions(&user);
133: MSA_InitialPoint(&user,x);
134: TaoSetInitialVector(tao,x);
136: /* Set routines for function, gradient and hessian evaluation */
137: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);
139: VecGetLocalSize(x,&m);
140: MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
141: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
143: DMGetLocalToGlobalMapping(user.dm,&isltog);
144: MatSetLocalToGlobalMapping(user.H,isltog,isltog);
145: PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);
146: if (flg) {
147: MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
148: MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
149: MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
150: TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
151: } else {
152: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
153: }
155: /* Set Variable bounds */
156: MSA_Plate(xl,xu,(void*)&user);
157: TaoSetVariableBounds(tao,xl,xu);
159: /* Check for any tao command line options */
160: TaoSetFromOptions(tao);
162: /* SOLVE THE APPLICATION */
163: TaoSolve(tao);
165: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
167: /* Free TAO data structures */
168: TaoDestroy(&tao);
170: /* Free PETSc data structures */
171: VecDestroy(&x);
172: VecDestroy(&xl);
173: VecDestroy(&xu);
174: MatDestroy(&user.H);
175: VecDestroy(&user.localX);
176: VecDestroy(&user.localV);
177: VecDestroy(&user.Bottom);
178: VecDestroy(&user.Top);
179: VecDestroy(&user.Left);
180: VecDestroy(&user.Right);
181: DMDestroy(&user.dm);
182: if (flg) {
183: MatDestroy(&H_shell);
184: }
185: PetscFinalize();
186: return ierr;
187: }
189: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
191: Input Parameters:
192: . tao - the Tao context
193: . X - input vector
194: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
196: Output Parameters:
197: . fcn - the function value
198: . G - vector containing the newly evaluated gradient
200: Notes:
201: In this case, we discretize the domain and Create triangles. The
202: surface of each triangle is planar, whose surface area can be easily
203: computed. The total surface area is found by sweeping through the grid
204: and computing the surface area of the two triangles that have their
205: right angle at the grid point. The diagonal line segments on the
206: grid that define the triangles run from top left to lower right.
207: The numbering of points starts at the lower left and runs left to
208: right, then bottom to top.
209: */
210: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
211: {
212: AppCtx *user = (AppCtx *) userCtx;
214: PetscInt i,j,row;
215: PetscInt mx=user->mx, my=user->my;
216: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
217: PetscReal ft=0;
218: PetscReal zero=0.0;
219: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
220: PetscReal rhx=mx+1, rhy=my+1;
221: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
222: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
223: PetscReal *g, *x,*left,*right,*bottom,*top;
224: Vec localX = user->localX, localG = user->localV;
226: /* Get local mesh boundaries */
227: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
228: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
230: /* Scatter ghost points to local vector */
231: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
232: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
234: /* Initialize vector to zero */
235: VecSet(localG, zero);
237: /* Get pointers to vector data */
238: VecGetArray(localX,&x);
239: VecGetArray(localG,&g);
240: VecGetArray(user->Top,&top);
241: VecGetArray(user->Bottom,&bottom);
242: VecGetArray(user->Left,&left);
243: VecGetArray(user->Right,&right);
245: /* Compute function over the locally owned part of the mesh */
246: for (j=ys; j<ys+ym; j++){
247: for (i=xs; i< xs+xm; i++){
248: row=(j-gys)*gxm + (i-gxs);
250: xc = x[row];
251: xlt=xrb=xl=xr=xb=xt=xc;
253: if (i==0){ /* left side */
254: xl= left[j-ys+1];
255: xlt = left[j-ys+2];
256: } else {
257: xl = x[row-1];
258: }
260: if (j==0){ /* bottom side */
261: xb=bottom[i-xs+1];
262: xrb = bottom[i-xs+2];
263: } else {
264: xb = x[row-gxm];
265: }
267: if (i+1 == gxs+gxm){ /* right side */
268: xr=right[j-ys+1];
269: xrb = right[j-ys];
270: } else {
271: xr = x[row+1];
272: }
274: if (j+1==gys+gym){ /* top side */
275: xt=top[i-xs+1];
276: xlt = top[i-xs];
277: }else {
278: xt = x[row+gxm];
279: }
281: if (i>gxs && j+1<gys+gym){
282: xlt = x[row-1+gxm];
283: }
284: if (j>gys && i+1<gxs+gxm){
285: xrb = x[row+1-gxm];
286: }
288: d1 = (xc-xl);
289: d2 = (xc-xr);
290: d3 = (xc-xt);
291: d4 = (xc-xb);
292: d5 = (xr-xrb);
293: d6 = (xrb-xb);
294: d7 = (xlt-xl);
295: d8 = (xt-xlt);
297: df1dxc = d1*hydhx;
298: df2dxc = ( d1*hydhx + d4*hxdhy );
299: df3dxc = d3*hxdhy;
300: df4dxc = ( d2*hydhx + d3*hxdhy );
301: df5dxc = d2*hydhx;
302: df6dxc = d4*hxdhy;
304: d1 *= rhx;
305: d2 *= rhx;
306: d3 *= rhy;
307: d4 *= rhy;
308: d5 *= rhy;
309: d6 *= rhx;
310: d7 *= rhy;
311: d8 *= rhx;
313: f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
314: f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
315: f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
316: f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
317: f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
318: f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);
320: ft = ft + (f2 + f4);
322: df1dxc /= f1;
323: df2dxc /= f2;
324: df3dxc /= f3;
325: df4dxc /= f4;
326: df5dxc /= f5;
327: df6dxc /= f6;
329: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
331: }
332: }
335: /* Compute triangular areas along the border of the domain. */
336: if (xs==0){ /* left side */
337: for (j=ys; j<ys+ym; j++){
338: d3=(left[j-ys+1] - left[j-ys+2])*rhy;
339: d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
340: ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
341: }
342: }
343: if (ys==0){ /* bottom side */
344: for (i=xs; i<xs+xm; i++){
345: d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
346: d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
347: ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
348: }
349: }
351: if (xs+xm==mx){ /* right side */
352: for (j=ys; j< ys+ym; j++){
353: d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
354: d4=(right[j-ys]-right[j-ys+1])*rhy;
355: ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
356: }
357: }
358: if (ys+ym==my){ /* top side */
359: for (i=xs; i<xs+xm; i++){
360: d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
361: d4=(top[i-xs+1] - top[i-xs])*rhx;
362: ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
363: }
364: }
366: if (ys==0 && xs==0){
367: d1=(left[0]-left[1])*rhy;
368: d2=(bottom[0]-bottom[1])*rhx;
369: ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
370: }
371: if (ys+ym == my && xs+xm == mx){
372: d1=(right[ym+1] - right[ym])*rhy;
373: d2=(top[xm+1] - top[xm])*rhx;
374: ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
375: }
377: ft=ft*area;
378: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
381: /* Restore vectors */
382: VecRestoreArray(localX,&x);
383: VecRestoreArray(localG,&g);
384: VecRestoreArray(user->Left,&left);
385: VecRestoreArray(user->Top,&top);
386: VecRestoreArray(user->Bottom,&bottom);
387: VecRestoreArray(user->Right,&right);
389: /* Scatter values to global vector */
390: DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
391: DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);
393: PetscLogFlops(70*xm*ym);
395: return 0;
396: }
398: /* ------------------------------------------------------------------- */
399: /*
400: FormHessian - Evaluates Hessian matrix.
402: Input Parameters:
403: . tao - the Tao context
404: . x - input vector
405: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
407: Output Parameters:
408: . A - Hessian matrix
409: . B - optionally different preconditioning matrix
411: Notes:
412: Due to mesh point reordering with DMs, we must always work
413: with the local mesh points, and then transform them to the new
414: global numbering with the local-to-global mapping. We cannot work
415: directly with the global numbers for the original uniprocessor mesh!
417: Two methods are available for imposing this transformation
418: when setting matrix entries:
419: (A) MatSetValuesLocal(), using the local ordering (including
420: ghost points!)
421: - Do the following two steps once, before calling TaoSolve()
422: - Use DMGetISLocalToGlobalMapping() to extract the
423: local-to-global map from the DM
424: - Associate this map with the matrix by calling
425: MatSetLocalToGlobalMapping()
426: - Then set matrix entries using the local ordering
427: by calling MatSetValuesLocal()
428: (B) MatSetValues(), using the global ordering
429: - Use DMGetGlobalIndices() to extract the local-to-global map
430: - Then apply this map explicitly yourself
431: - Set matrix entries using the global ordering by calling
432: MatSetValues()
433: Option (A) seems cleaner/easier in many cases, and is the procedure
434: used in this example.
435: */
436: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
437: {
439: AppCtx *user = (AppCtx *) ptr;
440: PetscInt i,j,k,row;
441: PetscInt mx=user->mx, my=user->my;
442: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
443: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
444: PetscReal rhx=mx+1, rhy=my+1;
445: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
446: PetscReal hl,hr,ht,hb,hc,htl,hbr;
447: PetscReal *x,*left,*right,*bottom,*top;
448: PetscReal v[7];
449: Vec localX = user->localX;
450: PetscBool assembled;
453: /* Set various matrix options */
454: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
456: /* Initialize matrix entries to zero */
457: MatAssembled(Hessian,&assembled);
458: if (assembled){MatZeroEntries(Hessian);}
460: /* Get local mesh boundaries */
461: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
462: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
464: /* Scatter ghost points to local vector */
465: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
466: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
468: /* Get pointers to vector data */
469: VecGetArray(localX,&x);
470: VecGetArray(user->Top,&top);
471: VecGetArray(user->Bottom,&bottom);
472: VecGetArray(user->Left,&left);
473: VecGetArray(user->Right,&right);
475: /* Compute Hessian over the locally owned part of the mesh */
477: for (i=xs; i< xs+xm; i++){
479: for (j=ys; j<ys+ym; j++){
481: row=(j-gys)*gxm + (i-gxs);
483: xc = x[row];
484: xlt=xrb=xl=xr=xb=xt=xc;
486: /* Left side */
487: if (i==gxs){
488: xl= left[j-ys+1];
489: xlt = left[j-ys+2];
490: } else {
491: xl = x[row-1];
492: }
494: if (j==gys){
495: xb=bottom[i-xs+1];
496: xrb = bottom[i-xs+2];
497: } else {
498: xb = x[row-gxm];
499: }
501: if (i+1 == gxs+gxm){
502: xr=right[j-ys+1];
503: xrb = right[j-ys];
504: } else {
505: xr = x[row+1];
506: }
508: if (j+1==gys+gym){
509: xt=top[i-xs+1];
510: xlt = top[i-xs];
511: }else {
512: xt = x[row+gxm];
513: }
515: if (i>gxs && j+1<gys+gym){
516: xlt = x[row-1+gxm];
517: }
518: if (j>gys && i+1<gxs+gxm){
519: xrb = x[row+1-gxm];
520: }
523: d1 = (xc-xl)*rhx;
524: d2 = (xc-xr)*rhx;
525: d3 = (xc-xt)*rhy;
526: d4 = (xc-xb)*rhy;
527: d5 = (xrb-xr)*rhy;
528: d6 = (xrb-xb)*rhx;
529: d7 = (xlt-xl)*rhy;
530: d8 = (xlt-xt)*rhx;
532: f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
533: f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
534: f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
535: f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
536: f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
537: f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);
540: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
541: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
542: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
543: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
544: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
545: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
546: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
547: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
549: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
550: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
552: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
553: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
554: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
555: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
557: hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
559: k=0;
560: if (j>0){
561: v[k]=hb; col[k]=row - gxm; k++;
562: }
564: if (j>0 && i < mx -1){
565: v[k]=hbr; col[k]=row - gxm+1; k++;
566: }
568: if (i>0){
569: v[k]= hl; col[k]=row - 1; k++;
570: }
572: v[k]= hc; col[k]=row; k++;
574: if (i < mx-1 ){
575: v[k]= hr; col[k]=row+1; k++;
576: }
578: if (i>0 && j < my-1 ){
579: v[k]= htl; col[k] = row+gxm-1; k++;
580: }
582: if (j < my-1 ){
583: v[k]= ht; col[k] = row+gxm; k++;
584: }
586: /*
587: Set matrix values using local numbering, which was defined
588: earlier, in the main routine.
589: */
590: MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);
592: }
593: }
595: /* Restore vectors */
596: VecRestoreArray(localX,&x);
597: VecRestoreArray(user->Left,&left);
598: VecRestoreArray(user->Top,&top);
599: VecRestoreArray(user->Bottom,&bottom);
600: VecRestoreArray(user->Right,&right);
602: /* Assemble the matrix */
603: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
604: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
606: PetscLogFlops(199*xm*ym);
607: return 0;
608: }
610: /* ------------------------------------------------------------------- */
611: /*
612: MSA_BoundaryConditions - Calculates the boundary conditions for
613: the region.
615: Input Parameter:
616: . user - user-defined application context
618: Output Parameter:
619: . user - user-defined application context
620: */
621: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
622: {
623: int ierr;
624: PetscInt i,j,k,maxits=5,limit=0;
625: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
626: PetscInt mx=user->mx,my=user->my;
627: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
628: PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
629: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
630: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
631: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
632: PetscReal *boundary;
633: PetscBool flg;
634: Vec Bottom,Top,Right,Left;
636: /* Get local mesh boundaries */
637: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
638: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
641: bsize=xm+2;
642: lsize=ym+2;
643: rsize=ym+2;
644: tsize=xm+2;
646: VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
647: VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
648: VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
649: VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);
651: user->Top=Top;
652: user->Left=Left;
653: user->Bottom=Bottom;
654: user->Right=Right;
656: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
658: for (j=0; j<4; j++){
659: if (j==0){
660: yt=b;
661: xt=l+hx*xs;
662: limit=bsize;
663: VecGetArray(Bottom,&boundary);
664: } else if (j==1){
665: yt=t;
666: xt=l+hx*xs;
667: limit=tsize;
668: VecGetArray(Top,&boundary);
669: } else if (j==2){
670: yt=b+hy*ys;
671: xt=l;
672: limit=lsize;
673: VecGetArray(Left,&boundary);
674: } else if (j==3){
675: yt=b+hy*ys;
676: xt=r;
677: limit=rsize;
678: VecGetArray(Right,&boundary);
679: }
681: for (i=0; i<limit; i++){
682: u1=xt;
683: u2=-yt;
684: for (k=0; k<maxits; k++){
685: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
686: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
687: fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
688: if (fnorm <= tol) break;
689: njac11=one+u2*u2-u1*u1;
690: njac12=two*u1*u2;
691: njac21=-two*u1*u2;
692: njac22=-one - u1*u1 + u2*u2;
693: det = njac11*njac22-njac21*njac12;
694: u1 = u1-(njac22*nf1-njac12*nf2)/det;
695: u2 = u2-(njac11*nf2-njac21*nf1)/det;
696: }
698: boundary[i]=u1*u1-u2*u2;
699: if (j==0 || j==1) {
700: xt=xt+hx;
701: } else if (j==2 || j==3){
702: yt=yt+hy;
703: }
704: }
705: if (j==0){
706: VecRestoreArray(Bottom,&boundary);
707: } else if (j==1){
708: VecRestoreArray(Top,&boundary);
709: } else if (j==2){
710: VecRestoreArray(Left,&boundary);
711: } else if (j==3){
712: VecRestoreArray(Right,&boundary);
713: }
714: }
716: /* Scale the boundary if desired */
718: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
719: if (flg){
720: VecScale(Bottom, scl);
721: }
722: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
723: if (flg){
724: VecScale(Top, scl);
725: }
726: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
727: if (flg){
728: VecScale(Right, scl);
729: }
731: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
732: if (flg){
733: VecScale(Left, scl);
734: }
735: return 0;
736: }
739: /* ------------------------------------------------------------------- */
740: /*
741: MSA_Plate - Calculates an obstacle for surface to stretch over.
743: Input Parameter:
744: . user - user-defined application context
746: Output Parameter:
747: . user - user-defined application context
748: */
749: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx){
751: AppCtx *user=(AppCtx *)ctx;
753: PetscInt i,j,row;
754: PetscInt xs,ys,xm,ym;
755: PetscInt mx=user->mx, my=user->my, bmy, bmx;
756: PetscReal t1,t2,t3;
757: PetscReal *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
758: PetscBool cylinder;
760: user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
761: user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
762: bmy=user->bmy; bmx=user->bmx;
764: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
766: VecSet(XL, lb);
767: VecSet(XU, ub);
769: VecGetArray(XL,&xl);
771: PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);
772: /* Compute the optional lower box */
773: if (cylinder){
774: for (i=xs; i< xs+xm; i++){
775: for (j=ys; j<ys+ym; j++){
776: row=(j-ys)*xm + (i-xs);
777: t1=(2.0*i-mx)*bmy;
778: t2=(2.0*j-my)*bmx;
779: t3=bmx*bmx*bmy*bmy;
780: if ( t1*t1 + t2*t2 <= t3 ){
781: xl[row] = user->bheight;
782: }
783: }
784: }
785: } else {
786: /* Compute the optional lower box */
787: for (i=xs; i< xs+xm; i++){
788: for (j=ys; j<ys+ym; j++){
789: row=(j-ys)*xm + (i-xs);
790: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
791: j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
792: xl[row] = user->bheight;
793: }
794: }
795: }
796: }
797: VecRestoreArray(XL,&xl);
799: return 0;
800: }
803: /* ------------------------------------------------------------------- */
804: /*
805: MSA_InitialPoint - Calculates the initial guess in one of three ways.
807: Input Parameters:
808: . user - user-defined application context
809: . X - vector for initial guess
811: Output Parameters:
812: . X - newly computed initial guess
813: */
814: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
815: {
817: PetscInt start=-1,i,j;
818: PetscReal zero=0.0;
819: PetscBool flg;
821: PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
822: if (flg && start==0){ /* The zero vector is reasonable */
823: VecSet(X, zero);
824: } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */
825: PetscRandom rctx; PetscReal np5=-0.5;
827: PetscRandomCreate(MPI_COMM_WORLD,&rctx);
828: for (i=0; i<start; i++){
829: VecSetRandom(X, rctx);
830: }
831: PetscRandomDestroy(&rctx);
832: VecShift(X, np5);
834: } else { /* Take an average of the boundary conditions */
836: PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
837: PetscInt mx=user->mx,my=user->my;
838: PetscReal *x,*left,*right,*bottom,*top;
839: Vec localX = user->localX;
841: /* Get local mesh boundaries */
842: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
843: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
845: /* Get pointers to vector data */
846: VecGetArray(user->Top,&top);
847: VecGetArray(user->Bottom,&bottom);
848: VecGetArray(user->Left,&left);
849: VecGetArray(user->Right,&right);
851: VecGetArray(localX,&x);
852: /* Perform local computations */
853: for (j=ys; j<ys+ym; j++){
854: for (i=xs; i< xs+xm; i++){
855: row=(j-gys)*gxm + (i-gxs);
856: x[row] = ( (j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+
857: (i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
858: }
859: }
861: /* Restore vectors */
862: VecRestoreArray(localX,&x);
864: VecRestoreArray(user->Left,&left);
865: VecRestoreArray(user->Top,&top);
866: VecRestoreArray(user->Bottom,&bottom);
867: VecRestoreArray(user->Right,&right);
869: /* Scatter values into global vector */
870: DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
871: DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);
873: }
874: return 0;
875: }
877: /* For testing matrix free submatrices */
878: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
879: {
881: AppCtx *user = (AppCtx*)ptr;
883: FormHessian(tao,x,user->H,user->H,ptr);
884: return(0);
885: }
886: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
887: {
889: void *ptr;
890: AppCtx *user;
892: MatShellGetContext(H_shell,&ptr);
893: user = (AppCtx*)ptr;
894: MatMult(user->H,X,Y);
895: return(0);
896: }