Actual source code: minsurf2.c
petsc-3.6.1 2015-08-06
1: /* Program usage: mpirun -np <proc> minsurf2 [-help] [all TAO options] */
3: /*
4: Include "petsctao.h" so we can use TAO solvers.
5: petscdmda.h for distributed array
6: */
7: #include <petsctao.h>
8: #include <petscdmda.h>
10: static char help[] =
11: "This example demonstrates use of the TAO package to \n\
12: solve an unconstrained minimization problem. This example is based on a \n\
13: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
14: boundary values along the edges of the domain, the objective is to find the\n\
15: surface with the minimal area that satisfies the boundary conditions.\n\
16: The command line options are:\n\
17: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
18: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
19: -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
20: for an average of the boundary conditions\n\n";
22: /*T
23: Concepts: TAO^Solving an unconstrained minimization problem
24: Routines: TaoCreate(); TaoSetType();
25: Routines: TaoSetInitialVector();
26: Routines: TaoSetObjectiveAndGradientRoutine();
27: Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
28: Routines: TaoSetMonitor();
29: Routines: TaoSolve(); TaoView();
30: Routines: TaoGetConvergedReason(); TaoDestroy();
31: Processors: n
32: T*/
34: /*
35: User-defined application context - contains data needed by the
36: application-provided call-back routines, FormFunctionGradient()
37: and FormHessian().
38: */
39: typedef struct {
40: PetscInt mx, my; /* discretization in x, y directions */
41: PetscReal *bottom, *top, *left, *right; /* boundary values */
42: DM dm; /* distributed array data structure */
43: Mat H; /* Hessian */
44: } AppCtx;
47: /* -------- User-defined Routines --------- */
49: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
50: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
51: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
52: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*);
53: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
54: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
55: PetscErrorCode My_Monitor(Tao, void *);
59: int main( int argc, char **argv )
60: {
61: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
62: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
63: Vec x; /* solution, gradient vectors */
64: PetscBool flg, viewmat; /* flags */
65: PetscBool fddefault, fdcoloring; /* flags */
66: TaoConvergedReason reason;
67: Tao tao; /* TAO solver context */
68: AppCtx user; /* user-defined work context */
69: ISColoring iscoloring;
70: MatFDColoring matfdcoloring;
72: /* Initialize TAO */
73: PetscInitialize( &argc, &argv,(char *)0,help );
75: /* Specify dimension of the problem */
76: user.mx = 10; user.my = 10;
78: /* Check for any command line arguments that override defaults */
79: PetscOptionsGetInt(NULL,"-mx",&user.mx,&flg);
80: PetscOptionsGetInt(NULL,"-my",&user.my,&flg);
82: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
83: PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
86: /* Let PETSc determine the vector distribution */
87: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
89: /* Create distributed array (DM) to manage parallel grid and vectors */
90: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx, user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
93: /* Create TAO solver and set desired solution method.*/
94: TaoCreate(PETSC_COMM_WORLD,&tao);
95: TaoSetType(tao,TAOCG);
97: /*
98: Extract global vector from DA for the vector of variables -- PETSC routine
99: Compute the initial solution -- application specific, see below
100: Set this vector for use by TAO -- TAO routine
101: */
102: DMCreateGlobalVector(user.dm,&x);
103: MSA_BoundaryConditions(&user);
104: MSA_InitialPoint(&user,x);
105: TaoSetInitialVector(tao,x);
107: /*
108: Initialize the Application context for use in function evaluations -- application specific, see below.
109: Set routines for function and gradient evaluation
110: */
111: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
113: /*
114: Given the command line arguments, calculate the hessian with either the user-
115: provided function FormHessian, or the default finite-difference driven Hessian
116: functions
117: */
118: PetscOptionsHasName(NULL,"-tao_fddefault",&fddefault);
119: PetscOptionsHasName(NULL,"-tao_fdcoloring",&fdcoloring);
122: /*
123: Create a matrix data structure to store the Hessian and set
124: the Hessian evalution routine.
125: Set the matrix structure to be used for Hessian evalutions
126: */
127: DMCreateMatrix(user.dm,&user.H);
128: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
131: if (fdcoloring) {
132: DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
133:
135: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
136:
138: ISColoringDestroy(&iscoloring);
139: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
140: MatFDColoringSetFromOptions(matfdcoloring);
142: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
144: } else if (fddefault){
145: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
147: } else {
148: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
149: }
152: /*
153: If my_monitor option is in command line, then use the user-provided
154: monitoring function
155: */
156: PetscOptionsHasName(NULL,"-my_monitor",&viewmat);
157: if (viewmat){
158: TaoSetMonitor(tao,My_Monitor,NULL,NULL);
159: }
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: /* Get information on termination */
170: TaoGetConvergedReason(tao,&reason);
171: if (reason <= 0 ){
172: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method \n");
173: }
176: /* Free TAO data structures */
177: TaoDestroy(&tao);
179: /* Free PETSc data structures */
180: VecDestroy(&x);
181: MatDestroy(&user.H);
182: if (fdcoloring) {
183: MatFDColoringDestroy(&matfdcoloring);
184: }
185: PetscFree(user.bottom);
186: PetscFree(user.top);
187: PetscFree(user.left);
188: PetscFree(user.right);
189: DMDestroy(&user.dm);
191: PetscFinalize();
192: return 0;
193: }
197: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx){
199: PetscReal fcn;
201: FormFunctionGradient(tao,X,&fcn,G,userCtx);
202: return(0);
203: }
205: /* -------------------------------------------------------------------- */
208: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
210: Input Parameters:
211: . tao - the Tao context
212: . XX - input vector
213: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
215: Output Parameters:
216: . fcn - the newly evaluated function
217: . GG - vector containing the newly evaluated gradient
218: */
219: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){
221: AppCtx *user = (AppCtx *) userCtx;
223: PetscInt i,j;
224: PetscInt mx=user->mx, my=user->my;
225: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
226: PetscReal ft=0.0;
227: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
228: PetscReal rhx=mx+1, rhy=my+1;
229: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
230: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
231: PetscReal **g, **x;
232: Vec localX;
235: /* Get local mesh boundaries */
236: DMGetLocalVector(user->dm,&localX);
238: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
239: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
241: /* Scatter ghost points to local vector */
242: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
243: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
245: /* Get pointers to vector data */
246: DMDAVecGetArray(user->dm,localX,(void**)&x);
247: DMDAVecGetArray(user->dm,G,(void**)&g);
249: /* Compute function and gradient over the locally owned part of the mesh */
250: for (j=ys; j<ys+ym; j++){
251: for (i=xs; i< xs+xm; i++){
253: xc = x[j][i];
254: xlt=xrb=xl=xr=xb=xt=xc;
256: if (i==0){ /* left side */
257: xl= user->left[j-ys+1];
258: xlt = user->left[j-ys+2];
259: } else {
260: xl = x[j][i-1];
261: }
263: if (j==0){ /* bottom side */
264: xb=user->bottom[i-xs+1];
265: xrb =user->bottom[i-xs+2];
266: } else {
267: xb = x[j-1][i];
268: }
270: if (i+1 == gxs+gxm){ /* right side */
271: xr=user->right[j-ys+1];
272: xrb = user->right[j-ys];
273: } else {
274: xr = x[j][i+1];
275: }
277: if (j+1==gys+gym){ /* top side */
278: xt=user->top[i-xs+1];
279: xlt = user->top[i-xs];
280: }else {
281: xt = x[j+1][i];
282: }
284: if (i>gxs && j+1<gys+gym){
285: xlt = x[j+1][i-1];
286: }
287: if (j>gys && i+1<gxs+gxm){
288: xrb = x[j-1][i+1];
289: }
291: d1 = (xc-xl);
292: d2 = (xc-xr);
293: d3 = (xc-xt);
294: d4 = (xc-xb);
295: d5 = (xr-xrb);
296: d6 = (xrb-xb);
297: d7 = (xlt-xl);
298: d8 = (xt-xlt);
300: df1dxc = d1*hydhx;
301: df2dxc = ( d1*hydhx + d4*hxdhy );
302: df3dxc = d3*hxdhy;
303: df4dxc = ( d2*hydhx + d3*hxdhy );
304: df5dxc = d2*hydhx;
305: df6dxc = d4*hxdhy;
307: d1 *= rhx;
308: d2 *= rhx;
309: d3 *= rhy;
310: d4 *= rhy;
311: d5 *= rhy;
312: d6 *= rhx;
313: d7 *= rhy;
314: d8 *= rhx;
316: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
317: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
318: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
319: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
320: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
321: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
323: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
324: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
326: ft = ft + (f2 + f4);
328: df1dxc /= f1;
329: df2dxc /= f2;
330: df3dxc /= f3;
331: df4dxc /= f4;
332: df5dxc /= f5;
333: df6dxc /= f6;
335: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
337: }
338: }
340: /* Compute triangular areas along the border of the domain. */
341: if (xs==0){ /* left side */
342: for (j=ys; j<ys+ym; j++){
343: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
344: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
345: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
346: }
347: }
348: if (ys==0){ /* bottom side */
349: for (i=xs; i<xs+xm; i++){
350: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
351: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
352: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
353: }
354: }
356: if (xs+xm==mx){ /* right side */
357: for (j=ys; j< ys+ym; j++){
358: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
359: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
360: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
361: }
362: }
363: if (ys+ym==my){ /* top side */
364: for (i=xs; i<xs+xm; i++){
365: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
366: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
367: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
368: }
369: }
371: if (ys==0 && xs==0){
372: d1=(user->left[0]-user->left[1])/hy;
373: d2=(user->bottom[0]-user->bottom[1])*rhx;
374: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
375: }
376: if (ys+ym == my && xs+xm == mx){
377: d1=(user->right[ym+1] - user->right[ym])*rhy;
378: d2=(user->top[xm+1] - user->top[xm])*rhx;
379: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
380: }
382: ft=ft*area;
383: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
385: /* Restore vectors */
386: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
387: DMDAVecRestoreArray(user->dm,G,(void**)&g);
389: /* Scatter values to global vector */
390: DMRestoreLocalVector(user->dm,&localX);
392: PetscLogFlops(67*xm*ym);
394: return(0);
395: }
397: /* ------------------------------------------------------------------- */
400: /*
401: FormHessian - Evaluates Hessian matrix.
403: Input Parameters:
404: . tao - the Tao context
405: . x - input vector
406: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
408: Output Parameters:
409: . H - Hessian matrix
410: . Hpre - optionally different preconditioning matrix
411: . flg - flag indicating matrix structure
413: */
414: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
415: {
417: AppCtx *user = (AppCtx *) ptr;
420: /* Evaluate the Hessian entries*/
421: QuadraticH(user,X,H);
422: return(0);
423: }
425: /* ------------------------------------------------------------------- */
428: /*
429: QuadraticH - Evaluates Hessian matrix.
431: Input Parameters:
432: . user - user-defined context, as set by TaoSetHessian()
433: . X - input vector
435: Output Parameter:
436: . H - Hessian matrix
437: */
438: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
439: {
440: PetscErrorCode ierr;
441: PetscInt i,j,k;
442: PetscInt mx=user->mx, my=user->my;
443: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
444: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
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, v[7];
448: MatStencil col[7],row;
449: Vec localX;
450: PetscBool assembled;
453: /* Get local mesh boundaries */
454: DMGetLocalVector(user->dm,&localX);
456: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
457: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
459: /* Scatter ghost points to local vector */
460: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
461: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
463: /* Get pointers to vector data */
464: DMDAVecGetArray(user->dm,localX,(void**)&x);
466: /* Initialize matrix entries to zero */
467: MatAssembled(Hessian,&assembled);
468: if (assembled){MatZeroEntries(Hessian); }
471: /* Set various matrix options */
472: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
474: /* Compute Hessian over the locally owned part of the mesh */
476: for (j=ys; j<ys+ym; j++){
478: for (i=xs; i< xs+xm; i++){
480: xc = x[j][i];
481: xlt=xrb=xl=xr=xb=xt=xc;
483: /* Left side */
484: if (i==0){
485: xl = user->left[j-ys+1];
486: xlt = user->left[j-ys+2];
487: } else {
488: xl = x[j][i-1];
489: }
491: if (j==0){
492: xb = user->bottom[i-xs+1];
493: xrb = user->bottom[i-xs+2];
494: } else {
495: xb = x[j-1][i];
496: }
498: if (i+1 == mx){
499: xr = user->right[j-ys+1];
500: xrb = user->right[j-ys];
501: } else {
502: xr = x[j][i+1];
503: }
505: if (j+1==my){
506: xt = user->top[i-xs+1];
507: xlt = user->top[i-xs];
508: }else {
509: xt = x[j+1][i];
510: }
512: if (i>0 && j+1<my){
513: xlt = x[j+1][i-1];
514: }
515: if (j>0 && i+1<mx){
516: xrb = x[j-1][i+1];
517: }
520: d1 = (xc-xl)/hx;
521: d2 = (xc-xr)/hx;
522: d3 = (xc-xt)/hy;
523: d4 = (xc-xb)/hy;
524: d5 = (xrb-xr)/hy;
525: d6 = (xrb-xb)/hx;
526: d7 = (xlt-xl)/hy;
527: d8 = (xlt-xt)/hx;
529: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
530: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
531: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
532: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
533: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
534: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
537: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
538: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
539: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
540: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
541: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
542: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
543: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
544: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
546: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
547: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
549: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
550: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
551: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
552: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
554: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
556: row.j = j; row.i = i;
557: k=0;
558: if (j>0){
559: v[k]=hb;
560: col[k].j = j - 1; col[k].i = i;
561: k++;
562: }
564: if (j>0 && i < mx -1){
565: v[k]=hbr;
566: col[k].j = j - 1; col[k].i = i+1;
567: k++;
568: }
570: if (i>0){
571: v[k]= hl;
572: col[k].j = j; col[k].i = i-1;
573: k++;
574: }
576: v[k]= hc;
577: col[k].j = j; col[k].i = i;
578: k++;
580: if (i < mx-1 ){
581: v[k]= hr;
582: col[k].j = j; col[k].i = i+1;
583: k++;
584: }
586: if (i>0 && j < my-1 ){
587: v[k]= htl;
588: col[k].j = j+1; col[k].i = i-1;
589: k++;
590: }
592: if (j < my-1 ){
593: v[k]= ht;
594: col[k].j = j+1; col[k].i = i;
595: k++;
596: }
598: /*
599: Set matrix values using local numbering, which was defined
600: earlier, in the main routine.
601: */
602: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
603:
605: }
606: }
608: /* Restore vectors */
609: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
611: DMRestoreLocalVector(user->dm,&localX);
613: /* Assemble the matrix */
614: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
615: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
617: PetscLogFlops(199*xm*ym);
618: return(0);
619: }
621: /* ------------------------------------------------------------------- */
624: /*
625: MSA_BoundaryConditions - Calculates the boundary conditions for
626: the region.
628: Input Parameter:
629: . user - user-defined application context
631: Output Parameter:
632: . user - user-defined application context
633: */
634: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
635: {
636: PetscErrorCode ierr;
637: PetscInt i,j,k,limit=0,maxits=5;
638: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
639: PetscInt mx=user->mx,my=user->my;
640: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
641: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
642: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
643: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
644: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
645: PetscReal *boundary;
646: PetscBool flg;
649: /* Get local mesh boundaries */
650: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
651: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
653: bsize=xm+2;
654: lsize=ym+2;
655: rsize=ym+2;
656: tsize=xm+2;
658: PetscMalloc1(bsize,&user->bottom);
659: PetscMalloc1(tsize,&user->top);
660: PetscMalloc1(lsize,&user->left);
661: PetscMalloc1(rsize,&user->right);
663: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
665: for (j=0; j<4; j++){
666: if (j==0){
667: yt=b;
668: xt=l+hx*xs;
669: limit=bsize;
670: boundary=user->bottom;
671: } else if (j==1){
672: yt=t;
673: xt=l+hx*xs;
674: limit=tsize;
675: boundary=user->top;
676: } else if (j==2){
677: yt=b+hy*ys;
678: xt=l;
679: limit=lsize;
680: boundary=user->left;
681: } else { /* if (j==3) */
682: yt=b+hy*ys;
683: xt=r;
684: limit=rsize;
685: boundary=user->right;
686: }
688: for (i=0; i<limit; i++){
689: u1=xt;
690: u2=-yt;
691: for (k=0; k<maxits; k++){
692: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
693: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
694: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
695: if (fnorm <= tol) break;
696: njac11=one+u2*u2-u1*u1;
697: njac12=two*u1*u2;
698: njac21=-two*u1*u2;
699: njac22=-one - u1*u1 + u2*u2;
700: det = njac11*njac22-njac21*njac12;
701: u1 = u1-(njac22*nf1-njac12*nf2)/det;
702: u2 = u2-(njac11*nf2-njac21*nf1)/det;
703: }
705: boundary[i]=u1*u1-u2*u2;
706: if (j==0 || j==1) {
707: xt=xt+hx;
708: } else { /* if (j==2 || j==3) */
709: yt=yt+hy;
710: }
712: }
714: }
716: /* Scale the boundary if desired */
717: if (1==1){
718: PetscReal scl = 1.0;
720: PetscOptionsGetReal(NULL,"-bottom",&scl,&flg);
721:
722: if (flg){
723: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
724: }
726: PetscOptionsGetReal(NULL,"-top",&scl,&flg);
727:
728: if (flg){
729: for (i=0;i<tsize;i++) user->top[i]*=scl;
730: }
732: PetscOptionsGetReal(NULL,"-right",&scl,&flg);
733:
734: if (flg){
735: for (i=0;i<rsize;i++) user->right[i]*=scl;
736: }
738: PetscOptionsGetReal(NULL,"-left",&scl,&flg);
739:
740: if (flg){
741: for (i=0;i<lsize;i++) user->left[i]*=scl;
742: }
743: }
745: return(0);
746: }
748: /* ------------------------------------------------------------------- */
751: /*
752: MSA_InitialPoint - Calculates the initial guess in one of three ways.
754: Input Parameters:
755: . user - user-defined application context
756: . X - vector for initial guess
758: Output Parameters:
759: . X - newly computed initial guess
760: */
761: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
762: {
763: PetscErrorCode ierr;
764: PetscInt start2=-1,i,j;
765: PetscReal start1=0;
766: PetscBool flg1,flg2;
769: PetscOptionsGetReal(NULL,"-start",&start1,&flg1);
770: PetscOptionsGetInt(NULL,"-random",&start2,&flg2);
772: if (flg1){ /* The zero vector is reasonable */
774: VecSet(X, start1);
776: } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */
778: PetscRandom rctx; PetscReal np5=-0.5;
780: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
781: PetscRandomSetType(rctx,PETSCRAND);
782:
783: for (i=0; i<start2; i++){
784: VecSetRandom(X, rctx);
785: }
786: PetscRandomDestroy(&rctx);
787: VecShift(X, np5);
789: } else { /* Take an average of the boundary conditions */
791: PetscInt xs,xm,ys,ym;
792: PetscInt mx=user->mx,my=user->my;
793: PetscReal **x;
795: /* Get local mesh boundaries */
796: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
798: /* Get pointers to vector data */
799: DMDAVecGetArray(user->dm,X,(void**)&x);
801: /* Perform local computations */
802: for (j=ys; j<ys+ym; j++){
803: for (i=xs; i< xs+xm; i++){
804: x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+
805: ((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
806: }
807: }
809: /* Restore vectors */
810: DMDAVecRestoreArray(user->dm,X,(void**)&x);
812: PetscLogFlops(9*xm*ym);
814: }
815: return(0);
816: }
818: /*-----------------------------------------------------------------------*/
821: PetscErrorCode My_Monitor(Tao tao, void *ctx){
823: Vec X;
825: TaoGetSolutionVector(tao,&X);
826: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
827: return(0);
828: }