Actual source code: plate2.c
petsc-3.14.6 2021-03-30
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*/
36: /*
37: User-defined application context - contains data needed by the
38: application-provided call-back routines, FormFunctionGradient(),
39: FormHessian().
40: */
41: typedef struct {
42: /* problem parameters */
43: PetscReal bheight; /* Height of plate under the surface */
44: PetscInt mx, my; /* discretization in x, y directions */
45: PetscInt bmx,bmy; /* Size of plate under the surface */
46: Vec Bottom, Top, Left, Right; /* boundary values */
48: /* Working space */
49: Vec localX, localV; /* ghosted local vector */
50: DM dm; /* distributed array data structure */
51: Mat H;
52: } AppCtx;
54: /* -------- User-defined Routines --------- */
56: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
57: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
58: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
59: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
60: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
62: /* For testing matrix free submatrices */
63: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
64: 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: Mat H_shell; /* to test matrix-free submatrices */
76: AppCtx user; /* user-defined work context */
78: /* Initialize PETSc, TAO */
79: PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;
81: /* Specify default dimension of the problem */
82: user.mx = 10; user.my = 10; user.bheight=0.1;
84: /* Check for any command line arguments that override defaults */
85: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
86: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
87: PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);
89: user.bmx = user.mx/2; user.bmy = user.my/2;
90: PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);
91: PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);
93: PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
94: 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);
96: /* Calculate any derived values from parameters */
97: N = user.mx*user.my;
99: /* Let Petsc determine the dimensions of the local vectors */
100: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
102: /*
103: A two dimensional distributed array will help define this problem,
104: which derives from an elliptic PDE on two dimensional domain. From
105: the distributed array, Create the vectors.
106: */
107: 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);
108: DMSetFromOptions(user.dm);
109: DMSetUp(user.dm);
110: /*
111: Extract global and local vectors from DM; The local vectors are
112: used solely as work space for the evaluation of the function,
113: gradient, and Hessian. Duplicate for remaining vectors that are
114: the same types.
115: */
116: DMCreateGlobalVector(user.dm,&x); /* Solution */
117: DMCreateLocalVector(user.dm,&user.localX);
118: VecDuplicate(user.localX,&user.localV);
120: VecDuplicate(x,&xl);
121: VecDuplicate(x,&xu);
123: /* The TAO code begins here */
125: /*
126: Create TAO solver and set desired solution method
127: The method must either be TAOTRON or TAOBLMVM
128: If TAOBLMVM is used, then hessian function is not called.
129: */
130: TaoCreate(PETSC_COMM_WORLD,&tao);
131: TaoSetType(tao,TAOBLMVM);
133: /* Set initial solution guess; */
134: MSA_BoundaryConditions(&user);
135: MSA_InitialPoint(&user,x);
136: TaoSetInitialVector(tao,x);
138: /* Set routines for function, gradient and hessian evaluation */
139: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);
141: VecGetLocalSize(x,&m);
142: MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
143: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
145: DMGetLocalToGlobalMapping(user.dm,&isltog);
146: MatSetLocalToGlobalMapping(user.H,isltog,isltog);
147: PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);
148: if (flg) {
149: MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
150: MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
151: MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
152: TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
153: } else {
154: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
155: }
157: /* Set Variable bounds */
158: MSA_Plate(xl,xu,(void*)&user);
159: TaoSetVariableBounds(tao,xl,xu);
161: /* Check for any tao command line options */
162: TaoSetFromOptions(tao);
164: /* SOLVE THE APPLICATION */
165: TaoSolve(tao);
167: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
169: /* Free TAO data structures */
170: TaoDestroy(&tao);
172: /* Free PETSc data structures */
173: VecDestroy(&x);
174: VecDestroy(&xl);
175: VecDestroy(&xu);
176: MatDestroy(&user.H);
177: VecDestroy(&user.localX);
178: VecDestroy(&user.localV);
179: VecDestroy(&user.Bottom);
180: VecDestroy(&user.Top);
181: VecDestroy(&user.Left);
182: VecDestroy(&user.Right);
183: DMDestroy(&user.dm);
184: if (flg) {
185: MatDestroy(&H_shell);
186: }
187: PetscFinalize();
188: return ierr;
189: }
191: /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
193: Input Parameters:
194: . tao - the Tao context
195: . X - input vector
196: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
198: Output Parameters:
199: . fcn - the function value
200: . G - vector containing the newly evaluated gradient
202: Notes:
203: In this case, we discretize the domain and Create triangles. The
204: surface of each triangle is planar, whose surface area can be easily
205: computed. The total surface area is found by sweeping through the grid
206: and computing the surface area of the two triangles that have their
207: right angle at the grid point. The diagonal line segments on the
208: grid that define the triangles run from top left to lower right.
209: The numbering of points starts at the lower left and runs left to
210: right, then bottom to top.
211: */
212: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
213: {
214: AppCtx *user = (AppCtx *) userCtx;
216: PetscInt i,j,row;
217: PetscInt mx=user->mx, my=user->my;
218: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
219: PetscReal ft=0;
220: PetscReal zero=0.0;
221: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
222: PetscReal rhx=mx+1, rhy=my+1;
223: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
224: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
225: PetscReal *g, *x,*left,*right,*bottom,*top;
226: Vec localX = user->localX, localG = user->localV;
228: /* Get local mesh boundaries */
229: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
230: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
232: /* Scatter ghost points to local vector */
233: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
234: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
236: /* Initialize vector to zero */
237: VecSet(localG, zero);
239: /* Get pointers to vector data */
240: VecGetArray(localX,&x);
241: VecGetArray(localG,&g);
242: VecGetArray(user->Top,&top);
243: VecGetArray(user->Bottom,&bottom);
244: VecGetArray(user->Left,&left);
245: VecGetArray(user->Right,&right);
247: /* Compute function over the locally owned part of the mesh */
248: for (j=ys; j<ys+ym; j++){
249: for (i=xs; i< xs+xm; i++){
250: row=(j-gys)*gxm + (i-gxs);
252: xc = x[row];
253: xlt=xrb=xl=xr=xb=xt=xc;
255: if (i==0){ /* left side */
256: xl= left[j-ys+1];
257: xlt = left[j-ys+2];
258: } else {
259: xl = x[row-1];
260: }
262: if (j==0){ /* bottom side */
263: xb=bottom[i-xs+1];
264: xrb = bottom[i-xs+2];
265: } else {
266: xb = x[row-gxm];
267: }
269: if (i+1 == gxs+gxm){ /* right side */
270: xr=right[j-ys+1];
271: xrb = right[j-ys];
272: } else {
273: xr = x[row+1];
274: }
276: if (j+1==gys+gym){ /* top side */
277: xt=top[i-xs+1];
278: xlt = top[i-xs];
279: }else {
280: xt = x[row+gxm];
281: }
283: if (i>gxs && j+1<gys+gym){
284: xlt = x[row-1+gxm];
285: }
286: if (j>gys && i+1<gxs+gxm){
287: xrb = x[row+1-gxm];
288: }
290: d1 = (xc-xl);
291: d2 = (xc-xr);
292: d3 = (xc-xt);
293: d4 = (xc-xb);
294: d5 = (xr-xrb);
295: d6 = (xrb-xb);
296: d7 = (xlt-xl);
297: d8 = (xt-xlt);
299: df1dxc = d1*hydhx;
300: df2dxc = (d1*hydhx + d4*hxdhy);
301: df3dxc = d3*hxdhy;
302: df4dxc = (d2*hydhx + d3*hxdhy);
303: df5dxc = d2*hydhx;
304: df6dxc = d4*hxdhy;
306: d1 *= rhx;
307: d2 *= rhx;
308: d3 *= rhy;
309: d4 *= rhy;
310: d5 *= rhy;
311: d6 *= rhx;
312: d7 *= rhy;
313: d8 *= rhx;
315: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
316: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
317: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
318: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
319: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
320: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
322: ft = ft + (f2 + f4);
324: df1dxc /= f1;
325: df2dxc /= f2;
326: df3dxc /= f3;
327: df4dxc /= f4;
328: df5dxc /= f5;
329: df6dxc /= f6;
331: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;
333: }
334: }
337: /* Compute triangular areas along the border of the domain. */
338: if (xs==0){ /* left side */
339: for (j=ys; j<ys+ym; j++){
340: d3=(left[j-ys+1] - left[j-ys+2])*rhy;
341: d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
342: ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
343: }
344: }
345: if (ys==0){ /* bottom side */
346: for (i=xs; i<xs+xm; i++){
347: d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
348: d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
349: ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
350: }
351: }
353: if (xs+xm==mx){ /* right side */
354: for (j=ys; j< ys+ym; j++){
355: d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
356: d4=(right[j-ys]-right[j-ys+1])*rhy;
357: ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
358: }
359: }
360: if (ys+ym==my){ /* top side */
361: for (i=xs; i<xs+xm; i++){
362: d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
363: d4=(top[i-xs+1] - top[i-xs])*rhx;
364: ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
365: }
366: }
368: if (ys==0 && xs==0){
369: d1=(left[0]-left[1])*rhy;
370: d2=(bottom[0]-bottom[1])*rhx;
371: ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
372: }
373: if (ys+ym == my && xs+xm == mx){
374: d1=(right[ym+1] - right[ym])*rhy;
375: d2=(top[xm+1] - top[xm])*rhx;
376: ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
377: }
379: ft=ft*area;
380: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
383: /* Restore vectors */
384: VecRestoreArray(localX,&x);
385: VecRestoreArray(localG,&g);
386: VecRestoreArray(user->Left,&left);
387: VecRestoreArray(user->Top,&top);
388: VecRestoreArray(user->Bottom,&bottom);
389: VecRestoreArray(user->Right,&right);
391: /* Scatter values to global vector */
392: DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
393: DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);
395: PetscLogFlops(70.0*xm*ym);
397: return 0;
398: }
400: /* ------------------------------------------------------------------- */
401: /*
402: FormHessian - Evaluates Hessian matrix.
404: Input Parameters:
405: . tao - the Tao context
406: . x - input vector
407: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
409: Output Parameters:
410: . A - Hessian matrix
411: . B - optionally different preconditioning matrix
413: Notes:
414: Due to mesh point reordering with DMs, we must always work
415: with the local mesh points, and then transform them to the new
416: global numbering with the local-to-global mapping. We cannot work
417: directly with the global numbers for the original uniprocessor mesh!
419: Two methods are available for imposing this transformation
420: when setting matrix entries:
421: (A) MatSetValuesLocal(), using the local ordering (including
422: ghost points!)
423: - Do the following two steps once, before calling TaoSolve()
424: - Use DMGetISLocalToGlobalMapping() to extract the
425: local-to-global map from the DM
426: - Associate this map with the matrix by calling
427: MatSetLocalToGlobalMapping()
428: - Then set matrix entries using the local ordering
429: by calling MatSetValuesLocal()
430: (B) MatSetValues(), using the global ordering
431: - Use DMGetGlobalIndices() to extract the local-to-global map
432: - Then apply this map explicitly yourself
433: - Set matrix entries using the global ordering by calling
434: MatSetValues()
435: Option (A) seems cleaner/easier in many cases, and is the procedure
436: used in this example.
437: */
438: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
439: {
441: AppCtx *user = (AppCtx *) ptr;
442: PetscInt i,j,k,row;
443: PetscInt mx=user->mx, my=user->my;
444: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
445: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
446: PetscReal rhx=mx+1, rhy=my+1;
447: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
448: PetscReal hl,hr,ht,hb,hc,htl,hbr;
449: PetscReal *x,*left,*right,*bottom,*top;
450: PetscReal v[7];
451: Vec localX = user->localX;
452: PetscBool assembled;
455: /* Set various matrix options */
456: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
458: /* Initialize matrix entries to zero */
459: MatAssembled(Hessian,&assembled);
460: if (assembled){MatZeroEntries(Hessian);}
462: /* Get local mesh boundaries */
463: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
464: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
466: /* Scatter ghost points to local vector */
467: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
468: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
470: /* Get pointers to vector data */
471: VecGetArray(localX,&x);
472: VecGetArray(user->Top,&top);
473: VecGetArray(user->Bottom,&bottom);
474: VecGetArray(user->Left,&left);
475: VecGetArray(user->Right,&right);
477: /* Compute Hessian over the locally owned part of the mesh */
479: for (i=xs; i< xs+xm; i++){
481: for (j=ys; j<ys+ym; j++){
483: row=(j-gys)*gxm + (i-gxs);
485: xc = x[row];
486: xlt=xrb=xl=xr=xb=xt=xc;
488: /* Left side */
489: if (i==gxs){
490: xl= left[j-ys+1];
491: xlt = left[j-ys+2];
492: } else {
493: xl = x[row-1];
494: }
496: if (j==gys){
497: xb=bottom[i-xs+1];
498: xrb = bottom[i-xs+2];
499: } else {
500: xb = x[row-gxm];
501: }
503: if (i+1 == gxs+gxm){
504: xr=right[j-ys+1];
505: xrb = right[j-ys];
506: } else {
507: xr = x[row+1];
508: }
510: if (j+1==gys+gym){
511: xt=top[i-xs+1];
512: xlt = top[i-xs];
513: }else {
514: xt = x[row+gxm];
515: }
517: if (i>gxs && j+1<gys+gym){
518: xlt = x[row-1+gxm];
519: }
520: if (j>gys && i+1<gxs+gxm){
521: xrb = x[row+1-gxm];
522: }
525: d1 = (xc-xl)*rhx;
526: d2 = (xc-xr)*rhx;
527: d3 = (xc-xt)*rhy;
528: d4 = (xc-xb)*rhy;
529: d5 = (xrb-xr)*rhy;
530: d6 = (xrb-xb)*rhx;
531: d7 = (xlt-xl)*rhy;
532: d8 = (xlt-xt)*rhx;
534: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
535: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
536: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
537: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
538: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
539: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
542: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
543: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
544: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
545: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
546: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
547: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
548: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
549: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
551: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
552: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
554: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
555: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
556: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
557: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
559: hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
561: k=0;
562: if (j>0){
563: v[k]=hb; col[k]=row - gxm; k++;
564: }
566: if (j>0 && i < mx -1){
567: v[k]=hbr; col[k]=row - gxm+1; k++;
568: }
570: if (i>0){
571: v[k]= hl; col[k]=row - 1; k++;
572: }
574: v[k]= hc; col[k]=row; k++;
576: if (i < mx-1){
577: v[k]= hr; col[k]=row+1; k++;
578: }
580: if (i>0 && j < my-1){
581: v[k]= htl; col[k] = row+gxm-1; k++;
582: }
584: if (j < my-1){
585: v[k]= ht; col[k] = row+gxm; k++;
586: }
588: /*
589: Set matrix values using local numbering, which was defined
590: earlier, in the main routine.
591: */
592: MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);
594: }
595: }
597: /* Restore vectors */
598: VecRestoreArray(localX,&x);
599: VecRestoreArray(user->Left,&left);
600: VecRestoreArray(user->Top,&top);
601: VecRestoreArray(user->Bottom,&bottom);
602: VecRestoreArray(user->Right,&right);
604: /* Assemble the matrix */
605: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
606: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
608: PetscLogFlops(199.0*xm*ym);
609: return 0;
610: }
612: /* ------------------------------------------------------------------- */
613: /*
614: MSA_BoundaryConditions - Calculates the boundary conditions for
615: the region.
617: Input Parameter:
618: . user - user-defined application context
620: Output Parameter:
621: . user - user-defined application context
622: */
623: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
624: {
625: int ierr;
626: PetscInt i,j,k,maxits=5,limit=0;
627: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
628: PetscInt mx=user->mx,my=user->my;
629: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
630: PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
631: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
632: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
633: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
634: PetscReal *boundary;
635: PetscBool flg;
636: Vec Bottom,Top,Right,Left;
638: /* Get local mesh boundaries */
639: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
640: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
643: bsize=xm+2;
644: lsize=ym+2;
645: rsize=ym+2;
646: tsize=xm+2;
648: VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
649: VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
650: VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
651: VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);
653: user->Top=Top;
654: user->Left=Left;
655: user->Bottom=Bottom;
656: user->Right=Right;
658: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
660: for (j=0; j<4; j++){
661: if (j==0){
662: yt=b;
663: xt=l+hx*xs;
664: limit=bsize;
665: VecGetArray(Bottom,&boundary);
666: } else if (j==1){
667: yt=t;
668: xt=l+hx*xs;
669: limit=tsize;
670: VecGetArray(Top,&boundary);
671: } else if (j==2){
672: yt=b+hy*ys;
673: xt=l;
674: limit=lsize;
675: VecGetArray(Left,&boundary);
676: } else if (j==3){
677: yt=b+hy*ys;
678: xt=r;
679: limit=rsize;
680: VecGetArray(Right,&boundary);
681: }
683: for (i=0; i<limit; i++){
684: u1=xt;
685: u2=-yt;
686: for (k=0; k<maxits; k++){
687: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
688: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
689: fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
690: if (fnorm <= tol) break;
691: njac11=one+u2*u2-u1*u1;
692: njac12=two*u1*u2;
693: njac21=-two*u1*u2;
694: njac22=-one - u1*u1 + u2*u2;
695: det = njac11*njac22-njac21*njac12;
696: u1 = u1-(njac22*nf1-njac12*nf2)/det;
697: u2 = u2-(njac11*nf2-njac21*nf1)/det;
698: }
700: boundary[i]=u1*u1-u2*u2;
701: if (j==0 || j==1) {
702: xt=xt+hx;
703: } else if (j==2 || j==3){
704: yt=yt+hy;
705: }
706: }
707: if (j==0){
708: VecRestoreArray(Bottom,&boundary);
709: } else if (j==1){
710: VecRestoreArray(Top,&boundary);
711: } else if (j==2){
712: VecRestoreArray(Left,&boundary);
713: } else if (j==3){
714: VecRestoreArray(Right,&boundary);
715: }
716: }
718: /* Scale the boundary if desired */
720: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
721: if (flg){
722: VecScale(Bottom, scl);
723: }
724: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
725: if (flg){
726: VecScale(Top, scl);
727: }
728: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
729: if (flg){
730: VecScale(Right, scl);
731: }
733: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
734: if (flg){
735: VecScale(Left, scl);
736: }
737: return 0;
738: }
741: /* ------------------------------------------------------------------- */
742: /*
743: MSA_Plate - Calculates an obstacle for surface to stretch over.
745: Input Parameter:
746: . user - user-defined application context
748: Output Parameter:
749: . user - user-defined application context
750: */
751: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx){
753: AppCtx *user=(AppCtx *)ctx;
755: PetscInt i,j,row;
756: PetscInt xs,ys,xm,ym;
757: PetscInt mx=user->mx, my=user->my, bmy, bmx;
758: PetscReal t1,t2,t3;
759: PetscReal *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
760: PetscBool cylinder;
762: user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
763: user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
764: bmy=user->bmy; bmx=user->bmx;
766: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
768: VecSet(XL, lb);
769: VecSet(XU, ub);
771: VecGetArray(XL,&xl);
773: PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);
774: /* Compute the optional lower box */
775: if (cylinder){
776: for (i=xs; i< xs+xm; i++){
777: for (j=ys; j<ys+ym; j++){
778: row=(j-ys)*xm + (i-xs);
779: t1=(2.0*i-mx)*bmy;
780: t2=(2.0*j-my)*bmx;
781: t3=bmx*bmx*bmy*bmy;
782: if (t1*t1 + t2*t2 <= t3){
783: xl[row] = user->bheight;
784: }
785: }
786: }
787: } else {
788: /* Compute the optional lower box */
789: for (i=xs; i< xs+xm; i++){
790: for (j=ys; j<ys+ym; j++){
791: row=(j-ys)*xm + (i-xs);
792: if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
793: j>=(my-bmy)/2 && j<my-(my-bmy)/2){
794: xl[row] = user->bheight;
795: }
796: }
797: }
798: }
799: VecRestoreArray(XL,&xl);
801: return 0;
802: }
805: /* ------------------------------------------------------------------- */
806: /*
807: MSA_InitialPoint - Calculates the initial guess in one of three ways.
809: Input Parameters:
810: . user - user-defined application context
811: . X - vector for initial guess
813: Output Parameters:
814: . X - newly computed initial guess
815: */
816: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
817: {
819: PetscInt start=-1,i,j;
820: PetscReal zero=0.0;
821: PetscBool flg;
823: PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
824: if (flg && start==0){ /* The zero vector is reasonable */
825: VecSet(X, zero);
826: } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */
827: PetscRandom rctx; PetscReal np5=-0.5;
829: PetscRandomCreate(MPI_COMM_WORLD,&rctx);
830: for (i=0; i<start; i++){
831: VecSetRandom(X, rctx);
832: }
833: PetscRandomDestroy(&rctx);
834: VecShift(X, np5);
836: } else { /* Take an average of the boundary conditions */
838: PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
839: PetscInt mx=user->mx,my=user->my;
840: PetscReal *x,*left,*right,*bottom,*top;
841: Vec localX = user->localX;
843: /* Get local mesh boundaries */
844: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
845: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
847: /* Get pointers to vector data */
848: VecGetArray(user->Top,&top);
849: VecGetArray(user->Bottom,&bottom);
850: VecGetArray(user->Left,&left);
851: VecGetArray(user->Right,&right);
853: VecGetArray(localX,&x);
854: /* Perform local computations */
855: for (j=ys; j<ys+ym; j++){
856: for (i=xs; i< xs+xm; i++){
857: row=(j-gys)*gxm + (i-gxs);
858: 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;
859: }
860: }
862: /* Restore vectors */
863: VecRestoreArray(localX,&x);
865: VecRestoreArray(user->Left,&left);
866: VecRestoreArray(user->Top,&top);
867: VecRestoreArray(user->Bottom,&bottom);
868: VecRestoreArray(user->Right,&right);
870: /* Scatter values into global vector */
871: DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
872: DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);
874: }
875: return 0;
876: }
878: /* For testing matrix free submatrices */
879: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
880: {
882: AppCtx *user = (AppCtx*)ptr;
884: FormHessian(tao,x,user->H,user->H,ptr);
885: return(0);
886: }
887: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
888: {
890: void *ptr;
891: AppCtx *user;
893: MatShellGetContext(H_shell,&ptr);
894: user = (AppCtx*)ptr;
895: MatMult(user->H,X,Y);
896: return(0);
897: }
900: /*TEST
902: build:
903: requires: !complex
905: test:
906: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
907: requires: !single
909: test:
910: suffix: 2
911: nsize: 2
912: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
913: requires: !single
915: test:
916: suffix: 3
917: nsize: 3
918: args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
919: requires: !single
921: test:
922: suffix: 4
923: nsize: 3
924: 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
925: requires: !single
927: test:
928: suffix: 5
929: nsize: 3
930: 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
931: requires: !single
933: test:
934: suffix: 6
935: nsize: 3
936: 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
937: requires: !single
939: test:
940: suffix: 7
941: nsize: 3
942: 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
943: requires: !single
945: test:
946: suffix: 8
947: 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
948: requires: !single
950: test:
951: suffix: 9
952: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
953: requires: !single
955: test:
956: suffix: 10
957: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
958: requires: !single
960: test:
961: suffix: 11
962: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
963: requires: !single
965: test:
966: suffix: 12
967: args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
968: requires: !single
970: test:
971: suffix: 13
972: 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
973: requires: !single
975: test:
976: suffix: 14
977: 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
978: requires: !single
980: test:
981: suffix: 15
982: 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
983: requires: !single
985: test:
986: suffix: 16
987: args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
988: requires: !single
990: test:
991: suffix: 17
992: 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
993: requires: !single
995: TEST*/