Actual source code: plate2.c
1: #include "petscdm.h" /*I "petscdm.h" I*/
2: #include "taosolver.h" /*I "taosolver.h" I*/
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: TaoInitialize(); TaoFinalize();
23: Routines: TaoCreate();
24: Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
25: Routines: TaoSetHessianRoutine();
26: Routines: TaoSetInitialVector();
27: Routines: TaoSetVariableBounds();
28: Routines: TaoSetFromOptions();
29: Routines: TaoSolve(); TaoView();
30: Routines: TaoGetTerminationReason(); TaoDestroy();
31: Processors: n
32: T*/
35: /*
36: User-defined application context - contains data needed by the
37: application-provided call-back routines, FormFunctionGradient(),
38: FormHessian().
39: */
40: typedef struct {
41: /* problem parameters */
42: PetscReal bheight; /* Height of plate under the surface */
43: PetscInt mx, my; /* discretization in x, y directions */
44: PetscInt bmx,bmy; /* Size of plate under the surface */
45: Vec Bottom, Top, Left, Right; /* boundary values */
46:
47: /* Working space */
48: Vec localX, localV; /* ghosted local vector */
49: DM dm; /* distributed array data structure */
50: Mat H;
51: } AppCtx;
53: /* -------- User-defined Routines --------- */
55: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
56: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
57: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
58: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal*,Vec,void*);
59: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure*,void*);
61: /* For testing matrix free submatrices */
62: PetscErrorCode MatrixFreeHessian(TaoSolver,Vec,Mat*, Mat*,MatStructure*,void*);
63: PetscErrorCode MyMatMult(Mat,Vec,Vec);
67: int main( int argc, char **argv )
68: {
70: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
71: PetscInt m, N; /* number of local and global elements in vectors */
72: Vec x,xl,xu; /* solution vector and bounds*/
73: PetscBool flg; /* A return variable when checking for user options */
74: TaoSolver tao; /* TaoSolver solver context */
75: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
76: TaoSolverTerminationReason reason;
77: Mat H_shell; /* to test matrix-free submatrices */
78: AppCtx user; /* user-defined work context */
80: /* Initialize PETSc, TAO */
81: PetscInitialize( &argc, &argv,(char *)0,help );
82: TaoInitialize( &argc, &argv,(char *)0,help );
84: /* Specify default dimension of the problem */
85: user.mx = 10; user.my = 10; user.bheight=0.1;
87: /* Check for any command line arguments that override defaults */
88: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg);
89: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg);
90: PetscOptionsGetReal(PETSC_NULL,"-bheight",&user.bheight,&flg);
92: user.bmx = user.mx/2; user.bmy = user.my/2;
93: PetscOptionsGetInt(PETSC_NULL,"-bmx",&user.bmx,&flg);
94: PetscOptionsGetInt(PETSC_NULL,"-bmy",&user.bmy,&flg);
96: PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
97: PetscPrintf(PETSC_COMM_WORLD,"mx:%D, my:%D, bmx:%D, bmy:%D, height:%G\n",
98: user.mx,user.my,user.bmx,user.bmy,user.bheight);
100: /* Calculate any derived values from parameters */
101: N = user.mx*user.my;
103: /* Let Petsc determine the dimensions of the local vectors */
104: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
106: /*
107: A two dimensional distributed array will help define this problem,
108: which derives from an elliptic PDE on two dimensional domain. From
109: the distributed array, Create the vectors.
110: */
111: DMDACreate2d(MPI_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,
112: DMDA_STENCIL_BOX,user.mx,user.my,Nx,Ny,1,1,
113: PETSC_NULL,PETSC_NULL,&user.dm);
115: /*
116: Extract global and local vectors from DM; The local vectors are
117: used solely as work space for the evaluation of the function,
118: gradient, and Hessian. Duplicate for remaining vectors that are
119: the same types.
120: */
121: DMCreateGlobalVector(user.dm,&x); /* Solution */
122: DMCreateLocalVector(user.dm,&user.localX);
123: VecDuplicate(user.localX,&user.localV);
125: VecDuplicate(x,&xl);
126: VecDuplicate(x,&xu);
130: /* The TAO code begins here */
132: /*
133: Create TAO solver and set desired solution method
134: The method must either be 'tao_tron' or 'tao_blmvm'
135: If blmvm is used, then hessian function is not called.
136: */
137: TaoCreate(PETSC_COMM_WORLD,&tao);
138: TaoSetType(tao,"tao_blmvm");
140: /* Set initial solution guess; */
141: MSA_BoundaryConditions(&user);
142: MSA_InitialPoint(&user,x);
143: TaoSetInitialVector(tao,x);
144:
145: /* Set routines for function, gradient and hessian evaluation */
146: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);
147:
149: VecGetLocalSize(x,&m);
150: MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL,
151: 3,PETSC_NULL,&(user.H));
152: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
153:
154: DMGetLocalToGlobalMapping(user.dm,&isltog);
155: MatSetLocalToGlobalMapping(user.H,isltog,isltog);
156: PetscOptionsHasName(PETSC_NULL,"-matrixfree",&flg);
157: if (flg) {
158: MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
159: MatShellSetOperation(H_shell,MATOP_MULT,(void(*)())MyMatMult);
160: MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
161: TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
162: } else {
163: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
164: }
165:
167: /* Set Variable bounds */
168: MSA_Plate(xl,xu,(void*)&user);
169: TaoSetVariableBounds(tao,xl,xu);
171: /* Check for any tao command line options */
172: TaoSetFromOptions(tao);
174: /* SOLVE THE APPLICATION */
175: TaoSolve(tao);
177: /* Get ierrrmation on termination */
178: TaoGetTerminationReason(tao,&reason);
179: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
180: if (reason <= 0){
181: PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
182: }
183: /* Free TAO data structures */
184: TaoDestroy(&tao);
186: /* Free PETSc data structures */
187: VecDestroy(&x);
188: VecDestroy(&xl);
189: VecDestroy(&xu);
190: MatDestroy(&user.H);
191: VecDestroy(&user.localX);
192: VecDestroy(&user.localV);
193: VecDestroy(&user.Bottom);
194: VecDestroy(&user.Top);
195: VecDestroy(&user.Left);
196: VecDestroy(&user.Right);
197: DMDestroy(&user.dm);
198: if (flg) {
199: MatDestroy(&H_shell);
200: }
201: /* Finalize TAO and PETSc */
202: TaoFinalize();
203: PetscFinalize();
204:
205: return 0;
206: }
210: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
212: Input Parameters:
213: . tao - the TaoSolver context
214: . X - input vector
215: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
216:
217: Output Parameters:
218: . fcn - the function value
219: . G - vector containing the newly evaluated gradient
221: Notes:
222: In this case, we discretize the domain and Create triangles. The
223: surface of each triangle is planar, whose surface area can be easily
224: computed. The total surface area is found by sweeping through the grid
225: and computing the surface area of the two triangles that have their
226: right angle at the grid point. The diagonal line segments on the
227: grid that define the triangles run from top left to lower right.
228: The numbering of points starts at the lower left and runs left to
229: right, then bottom to top.
230: */
231: PetscErrorCode FormFunctionGradient(TaoSolver tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
232: {
233: AppCtx * user = (AppCtx *) userCtx;
234: PetscErrorCode ierr;
235: PetscInt i,j,row;
236: PetscInt mx=user->mx, my=user->my;
237: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
238: PetscReal ft=0;
239: PetscReal zero=0.0;
240: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
241: PetscReal rhx=mx+1, rhy=my+1;
242: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
243: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
244: PetscReal *g, *x,*left,*right,*bottom,*top;
245: Vec localX = user->localX, localG = user->localV;
247: /* Get local mesh boundaries */
248: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
249: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
251: /* Scatter ghost points to local vector */
252: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
253: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
255: /* Initialize vector to zero */
256: VecSet(localG, zero);
258: /* Get pointers to vector data */
259: VecGetArray(localX,&x);
260: VecGetArray(localG,&g);
261: VecGetArray(user->Top,&top);
262: VecGetArray(user->Bottom,&bottom);
263: VecGetArray(user->Left,&left);
264: VecGetArray(user->Right,&right);
266: /* Compute function over the locally owned part of the mesh */
267: for (j=ys; j<ys+ym; j++){
268: for (i=xs; i< xs+xm; i++){
269: row=(j-gys)*gxm + (i-gxs);
270:
271: xc = x[row];
272: xlt=xrb=xl=xr=xb=xt=xc;
273:
274: if (i==0){ /* left side */
275: xl= left[j-ys+1];
276: xlt = left[j-ys+2];
277: } else {
278: xl = x[row-1];
279: }
281: if (j==0){ /* bottom side */
282: xb=bottom[i-xs+1];
283: xrb = bottom[i-xs+2];
284: } else {
285: xb = x[row-gxm];
286: }
287:
288: if (i+1 == gxs+gxm){ /* right side */
289: xr=right[j-ys+1];
290: xrb = right[j-ys];
291: } else {
292: xr = x[row+1];
293: }
295: if (j+1==gys+gym){ /* top side */
296: xt=top[i-xs+1];
297: xlt = top[i-xs];
298: }else {
299: xt = x[row+gxm];
300: }
302: if (i>gxs && j+1<gys+gym){
303: xlt = x[row-1+gxm];
304: }
305: if (j>gys && i+1<gxs+gxm){
306: xrb = x[row+1-gxm];
307: }
309: d1 = (xc-xl);
310: d2 = (xc-xr);
311: d3 = (xc-xt);
312: d4 = (xc-xb);
313: d5 = (xr-xrb);
314: d6 = (xrb-xb);
315: d7 = (xlt-xl);
316: d8 = (xt-xlt);
317:
318: df1dxc = d1*hydhx;
319: df2dxc = ( d1*hydhx + d4*hxdhy );
320: df3dxc = d3*hxdhy;
321: df4dxc = ( d2*hydhx + d3*hxdhy );
322: df5dxc = d2*hydhx;
323: df6dxc = d4*hxdhy;
325: d1 *= rhx;
326: d2 *= rhx;
327: d3 *= rhy;
328: d4 *= rhy;
329: d5 *= rhy;
330: d6 *= rhx;
331: d7 *= rhy;
332: d8 *= rhx;
334: f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
335: f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
336: f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
337: f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
338: f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
339: f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);
340:
341: ft = ft + (f2 + f4);
343: df1dxc /= f1;
344: df2dxc /= f2;
345: df3dxc /= f3;
346: df4dxc /= f4;
347: df5dxc /= f5;
348: df6dxc /= f6;
350: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
351:
352: }
353: }
356: /* Compute triangular areas along the border of the domain. */
357: if (xs==0){ /* left side */
358: for (j=ys; j<ys+ym; j++){
359: d3=(left[j-ys+1] - left[j-ys+2])*rhy;
360: d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
361: ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
362: }
363: }
364: if (ys==0){ /* bottom side */
365: for (i=xs; i<xs+xm; i++){
366: d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
367: d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
368: ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
369: }
370: }
372: if (xs+xm==mx){ /* right side */
373: for (j=ys; j< ys+ym; j++){
374: d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
375: d4=(right[j-ys]-right[j-ys+1])*rhy;
376: ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
377: }
378: }
379: if (ys+ym==my){ /* top side */
380: for (i=xs; i<xs+xm; i++){
381: d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
382: d4=(top[i-xs+1] - top[i-xs])*rhx;
383: ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
384: }
385: }
387: if (ys==0 && xs==0){
388: d1=(left[0]-left[1])*rhy;
389: d2=(bottom[0]-bottom[1])*rhx;
390: ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
391: }
392: if (ys+ym == my && xs+xm == mx){
393: d1=(right[ym+1] - right[ym])*rhy;
394: d2=(top[xm+1] - top[xm])*rhx;
395: ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
396: }
398: ft=ft*area;
399: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
401:
402: /* Restore vectors */
403: VecRestoreArray(localX,&x);
404: VecRestoreArray(localG,&g);
405: VecRestoreArray(user->Left,&left);
406: VecRestoreArray(user->Top,&top);
407: VecRestoreArray(user->Bottom,&bottom);
408: VecRestoreArray(user->Right,&right);
410: /* Scatter values to global vector */
411: DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
412: DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);
414: PetscLogFlops(70*xm*ym);
416: return 0;
417: }
419: /* ------------------------------------------------------------------- */
422: /*
423: FormHessian - Evaluates Hessian matrix.
425: Input Parameters:
426: . tao - the TaoSolver context
427: . x - input vector
428: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
430: Output Parameters:
431: . A - Hessian matrix
432: . B - optionally different preconditioning matrix
433: . flag - flag indicating matrix structure
435: Notes:
436: Due to mesh point reordering with DMs, we must always work
437: with the local mesh points, and then transform them to the new
438: global numbering with the local-to-global mapping. We cannot work
439: directly with the global numbers for the original uniprocessor mesh!
441: Two methods are available for imposing this transformation
442: when setting matrix entries:
443: (A) MatSetValuesLocal(), using the local ordering (including
444: ghost points!)
445: - Do the following two steps once, before calling TaoSolve()
446: - Use DMGetISLocalToGlobalMapping() to extract the
447: local-to-global map from the DM
448: - Associate this map with the matrix by calling
449: MatSetLocalToGlobalMapping()
450: - Then set matrix entries using the local ordering
451: by calling MatSetValuesLocal()
452: (B) MatSetValues(), using the global ordering
453: - Use DMGetGlobalIndices() to extract the local-to-global map
454: - Then apply this map explicitly yourself
455: - Set matrix entries using the global ordering by calling
456: MatSetValues()
457: Option (A) seems cleaner/easier in many cases, and is the procedure
458: used in this example.
459: */
460: PetscErrorCode FormHessian(TaoSolver tao,Vec X,Mat *Hptr, Mat *Hpc, MatStructure *flag, void *ptr)
461: {
462: PetscErrorCode ierr;
463: AppCtx *user = (AppCtx *) ptr;
464: Mat Hessian = *Hpc;
465: PetscInt i,j,k,row;
466: PetscInt mx=user->mx, my=user->my;
467: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
468: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
469: PetscReal rhx=mx+1, rhy=my+1;
470: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
471: PetscReal hl,hr,ht,hb,hc,htl,hbr;
472: PetscReal *x,*left,*right,*bottom,*top;
473: PetscReal v[7];
474: Vec localX = user->localX;
475: PetscBool assembled;
478: /* Set various matrix options */
479: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
481: /* Initialize matrix entries to zero */
482: MatAssembled(Hessian,&assembled);
483: if (assembled){MatZeroEntries(Hessian); }
485: /* Get local mesh boundaries */
486: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
487: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
489: /* Scatter ghost points to local vector */
490: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
491: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
493: /* Get pointers to vector data */
494: VecGetArray(localX,&x);
495: VecGetArray(user->Top,&top);
496: VecGetArray(user->Bottom,&bottom);
497: VecGetArray(user->Left,&left);
498: VecGetArray(user->Right,&right);
500: /* Compute Hessian over the locally owned part of the mesh */
502: for (i=xs; i< xs+xm; i++){
504: for (j=ys; j<ys+ym; j++){
506: row=(j-gys)*gxm + (i-gxs);
507:
508: xc = x[row];
509: xlt=xrb=xl=xr=xb=xt=xc;
511: /* Left side */
512: if (i==gxs){
513: xl= left[j-ys+1];
514: xlt = left[j-ys+2];
515: } else {
516: xl = x[row-1];
517: }
518:
519: if (j==gys){
520: xb=bottom[i-xs+1];
521: xrb = bottom[i-xs+2];
522: } else {
523: xb = x[row-gxm];
524: }
525:
526: if (i+1 == gxs+gxm){
527: xr=right[j-ys+1];
528: xrb = right[j-ys];
529: } else {
530: xr = x[row+1];
531: }
533: if (j+1==gys+gym){
534: xt=top[i-xs+1];
535: xlt = top[i-xs];
536: }else {
537: xt = x[row+gxm];
538: }
540: if (i>gxs && j+1<gys+gym){
541: xlt = x[row-1+gxm];
542: }
543: if (j>gys && i+1<gxs+gxm){
544: xrb = x[row+1-gxm];
545: }
548: d1 = (xc-xl)*rhx;
549: d2 = (xc-xr)*rhx;
550: d3 = (xc-xt)*rhy;
551: d4 = (xc-xb)*rhy;
552: d5 = (xrb-xr)*rhy;
553: d6 = (xrb-xb)*rhx;
554: d7 = (xlt-xl)*rhy;
555: d8 = (xlt-xt)*rhx;
556:
557: f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
558: f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
559: f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
560: f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
561: f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
562: f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);
565: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
566: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
567: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
568: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
569: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
570: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
571: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
572: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
574: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
575: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
577: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
578: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
579: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
580: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
582: hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
584: k=0;
585: if (j>0){
586: v[k]=hb; col[k]=row - gxm; k++;
587: }
588:
589: if (j>0 && i < mx -1){
590: v[k]=hbr; col[k]=row - gxm+1; k++;
591: }
592:
593: if (i>0){
594: v[k]= hl; col[k]=row - 1; k++;
595: }
596:
597: v[k]= hc; col[k]=row; k++;
598:
599: if (i < mx-1 ){
600: v[k]= hr; col[k]=row+1; k++;
601: }
602:
603: if (i>0 && j < my-1 ){
604: v[k]= htl; col[k] = row+gxm-1; k++;
605: }
606:
607: if (j < my-1 ){
608: v[k]= ht; col[k] = row+gxm; k++;
609: }
610:
611: /*
612: Set matrix values using local numbering, which was defined
613: earlier, in the main routine.
614: */
615: MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);
616:
617:
618: }
619: }
620:
621: /* Restore vectors */
622: VecRestoreArray(localX,&x);
623: VecRestoreArray(user->Left,&left);
624: VecRestoreArray(user->Top,&top);
625: VecRestoreArray(user->Bottom,&bottom);
626: VecRestoreArray(user->Right,&right);
628: /* Assemble the matrix */
629: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
630: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
632: PetscLogFlops(199*xm*ym);
633: return 0;
634: }
636: /* ------------------------------------------------------------------- */
639: /*
640: MSA_BoundaryConditions - Calculates the boundary conditions for
641: the region.
643: Input Parameter:
644: . user - user-defined application context
646: Output Parameter:
647: . user - user-defined application context
648: */
649: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
650: {
651: int ierr;
652: PetscInt i,j,k,maxits=5,limit=0;
653: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
654: PetscInt mx=user->mx,my=user->my;
655: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
656: PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
657: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
658: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
659: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
660: PetscReal *boundary;
661: PetscBool flg;
662: Vec Bottom,Top,Right,Left;
664: /* Get local mesh boundaries */
665: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
666: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
669: bsize=xm+2;
670: lsize=ym+2;
671: rsize=ym+2;
672: tsize=xm+2;
674: VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
675: VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
676: VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
677: VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);
679: user->Top=Top;
680: user->Left=Left;
681: user->Bottom=Bottom;
682: user->Right=Right;
684: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
686: for (j=0; j<4; j++){
687: if (j==0){
688: yt=b;
689: xt=l+hx*xs;
690: limit=bsize;
691: VecGetArray(Bottom,&boundary);
692: } else if (j==1){
693: yt=t;
694: xt=l+hx*xs;
695: limit=tsize;
696: VecGetArray(Top,&boundary);
697: } else if (j==2){
698: yt=b+hy*ys;
699: xt=l;
700: limit=lsize;
701: VecGetArray(Left,&boundary);
702: } else if (j==3){
703: yt=b+hy*ys;
704: xt=r;
705: limit=rsize;
706: VecGetArray(Right,&boundary);
707: }
709: for (i=0; i<limit; i++){
710: u1=xt;
711: u2=-yt;
712: for (k=0; k<maxits; k++){
713: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
714: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
715: fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
716: if (fnorm <= tol) break;
717: njac11=one+u2*u2-u1*u1;
718: njac12=two*u1*u2;
719: njac21=-two*u1*u2;
720: njac22=-one - u1*u1 + u2*u2;
721: det = njac11*njac22-njac21*njac12;
722: u1 = u1-(njac22*nf1-njac12*nf2)/det;
723: u2 = u2-(njac11*nf2-njac21*nf1)/det;
724: }
726: boundary[i]=u1*u1-u2*u2;
727: if (j==0 || j==1) {
728: xt=xt+hx;
729: } else if (j==2 || j==3){
730: yt=yt+hy;
731: }
732:
733: }
734:
735: if (j==0){
736: VecRestoreArray(Bottom,&boundary);
737: } else if (j==1){
738: VecRestoreArray(Top,&boundary);
739: } else if (j==2){
740: VecRestoreArray(Left,&boundary);
741: } else if (j==3){
742: VecRestoreArray(Right,&boundary);
743: }
745: }
747: /* Scale the boundary if desired */
749: PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
750:
751: if (flg){
752: VecScale(Bottom, scl);
753: }
754:
755: PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
756:
757: if (flg){
758: VecScale(Top, scl);
759: }
760:
761: PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
762:
763: if (flg){
764: VecScale(Right, scl);
765: }
766:
767: PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
768:
769: if (flg){
770: VecScale(Left, scl);
771: }
773: return 0;
774: }
777: /* ------------------------------------------------------------------- */
780: /*
781: MSA_Plate - Calculates an obstacle for surface to stretch over.
783: Input Parameter:
784: . user - user-defined application context
786: Output Parameter:
787: . user - user-defined application context
788: */
789: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx){
791: AppCtx *user=(AppCtx *)ctx;
793: PetscInt i,j,row;
794: PetscInt xs,ys,xm,ym;
795: PetscInt mx=user->mx, my=user->my, bmy, bmx;
796: PetscReal t1,t2,t3;
797: PetscReal *xl, lb=TAO_NINFINITY, ub=TAO_INFINITY;
798: PetscBool cylinder;
800: user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
801: user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
802: bmy=user->bmy, bmx=user->bmx;
804: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
806: VecSet(XL, lb);
807: VecSet(XU, ub);
809: VecGetArray(XL,&xl);
811: PetscOptionsHasName(PETSC_NULL,"-cylinder",&cylinder);
812: /* Compute the optional lower box */
813: if (cylinder){
814: for (i=xs; i< xs+xm; i++){
815: for (j=ys; j<ys+ym; j++){
816: row=(j-ys)*xm + (i-xs);
817: t1=(2.0*i-mx)*bmy;
818: t2=(2.0*j-my)*bmx;
819: t3=bmx*bmx*bmy*bmy;
820: if ( t1*t1 + t2*t2 <= t3 ){
821: xl[row] = user->bheight;
822: }
823: }
824: }
825: } else {
826: /* Compute the optional lower box */
827: for (i=xs; i< xs+xm; i++){
828: for (j=ys; j<ys+ym; j++){
829: row=(j-ys)*xm + (i-xs);
830: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
831: j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
832: xl[row] = user->bheight;
833: }
834: }
835: }
836: }
837: VecRestoreArray(XL,&xl);
839: return 0;
840: }
843: /* ------------------------------------------------------------------- */
846: /*
847: MSA_InitialPoint - Calculates the initial guess in one of three ways.
849: Input Parameters:
850: . user - user-defined application context
851: . X - vector for initial guess
853: Output Parameters:
854: . X - newly computed initial guess
855: */
856: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
857: {
858: PetscErrorCode ierr;
859: PetscInt start=-1,i,j;
860: PetscReal zero=0.0;
861: PetscBool flg;
863: PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg);
865: if (flg && start==0){ /* The zero vector is reasonable */
866:
867: VecSet(X, zero);
869: } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */
871: PetscRandom rctx; PetscReal np5=-0.5;
873: PetscRandomCreate(MPI_COMM_WORLD,&rctx);
874:
875: for (i=0; i<start; i++){
876: VecSetRandom(X, rctx);
877: }
878: PetscRandomDestroy(&rctx);
879: VecShift(X, np5);
881: } else { /* Take an average of the boundary conditions */
883: PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
884: PetscInt mx=user->mx,my=user->my;
885: PetscReal *x,*left,*right,*bottom,*top;
886: Vec localX = user->localX;
887:
888: /* Get local mesh boundaries */
889: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
890: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
891:
892: /* Get pointers to vector data */
893: VecGetArray(user->Top,&top);
894: VecGetArray(user->Bottom,&bottom);
895: VecGetArray(user->Left,&left);
896: VecGetArray(user->Right,&right);
898: VecGetArray(localX,&x);
899: /* Perform local computations */
900: for (j=ys; j<ys+ym; j++){
901: for (i=xs; i< xs+xm; i++){
902: row=(j-gys)*gxm + (i-gxs);
903: x[row] = ( (j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+
904: (i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
905: }
906: }
907:
908: /* Restore vectors */
909: VecRestoreArray(localX,&x);
911: VecRestoreArray(user->Left,&left);
912: VecRestoreArray(user->Top,&top);
913: VecRestoreArray(user->Bottom,&bottom);
914: VecRestoreArray(user->Right,&right);
915:
916: /* Scatter values into global vector */
917: DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
918: DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);
919:
920: }
921: return 0;
922: }
924: /* For testing matrix free submatrices */
927: PetscErrorCode MatrixFreeHessian(TaoSolver tao, Vec x, Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
928: {
930: AppCtx *user = (AppCtx*)ptr;
932: FormHessian(tao,x,&user->H,&user->H,flg,ptr);
933: return(0);
934: }
937: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
938: {
940: void *ptr;
941: AppCtx *user;
943: MatShellGetContext(H_shell,&ptr);
944: user = (AppCtx*)ptr;
945: MatMult(user->H,X,Y);
946: return(0);
947: }