Actual source code: minsurf2.c
petsc-3.8.4 2018-03-24
1: /* Program usage: mpiexec -n <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: 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 *);
57: int main( int argc, char **argv )
58: {
59: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
60: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
61: Vec x; /* solution, gradient vectors */
62: PetscBool flg, viewmat; /* flags */
63: PetscBool fddefault, fdcoloring; /* flags */
64: Tao tao; /* TAO solver context */
65: AppCtx user; /* user-defined work context */
66: ISColoring iscoloring;
67: MatFDColoring matfdcoloring;
69: /* Initialize TAO */
70: PetscInitialize( &argc, &argv,(char *)0,help );if (ierr) return ierr;
72: /* Specify dimension of the problem */
73: user.mx = 10; user.my = 10;
75: /* Check for any command line arguments that override defaults */
76: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
77: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
79: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
80: PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
83: /* Let PETSc determine the vector distribution */
84: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
86: /* Create distributed array (DM) to manage parallel grid and vectors */
87: 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);
88: DMSetFromOptions(user.dm);
89: DMSetUp(user.dm);
91: /* Create TAO solver and set desired solution method.*/
92: TaoCreate(PETSC_COMM_WORLD,&tao);
93: TaoSetType(tao,TAOCG);
95: /*
96: Extract global vector from DA for the vector of variables -- PETSC routine
97: Compute the initial solution -- application specific, see below
98: Set this vector for use by TAO -- TAO routine
99: */
100: DMCreateGlobalVector(user.dm,&x);
101: MSA_BoundaryConditions(&user);
102: MSA_InitialPoint(&user,x);
103: TaoSetInitialVector(tao,x);
105: /*
106: Initialize the Application context for use in function evaluations -- application specific, see below.
107: Set routines for function and gradient evaluation
108: */
109: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
111: /*
112: Given the command line arguments, calculate the hessian with either the user-
113: provided function FormHessian, or the default finite-difference driven Hessian
114: functions
115: */
116: PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault);
117: PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);
119: /*
120: Create a matrix data structure to store the Hessian and set
121: the Hessian evalution routine.
122: Set the matrix structure to be used for Hessian evalutions
123: */
124: DMCreateMatrix(user.dm,&user.H);
125: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
128: if (fdcoloring) {
129: DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
130: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
131: ISColoringDestroy(&iscoloring);
132: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
133: MatFDColoringSetFromOptions(matfdcoloring);
134: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
135: } else if (fddefault){
136: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
137: } else {
138: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
139: }
141: /*
142: If my_monitor option is in command line, then use the user-provided
143: monitoring function
144: */
145: PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
146: if (viewmat){
147: TaoSetMonitor(tao,My_Monitor,NULL,NULL);
148: }
150: /* Check for any tao command line options */
151: TaoSetFromOptions(tao);
153: /* SOLVE THE APPLICATION */
154: TaoSolve(tao);
156: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
158: /* Free TAO data structures */
159: TaoDestroy(&tao);
161: /* Free PETSc data structures */
162: VecDestroy(&x);
163: MatDestroy(&user.H);
164: if (fdcoloring) {
165: MatFDColoringDestroy(&matfdcoloring);
166: }
167: PetscFree(user.bottom);
168: PetscFree(user.top);
169: PetscFree(user.left);
170: PetscFree(user.right);
171: DMDestroy(&user.dm);
172: PetscFinalize();
173: return ierr;
174: }
176: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
177: {
179: PetscReal fcn;
182: FormFunctionGradient(tao,X,&fcn,G,userCtx);
183: return(0);
184: }
186: /* -------------------------------------------------------------------- */
187: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
189: Input Parameters:
190: . tao - the Tao context
191: . XX - input vector
192: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
194: Output Parameters:
195: . fcn - the newly evaluated function
196: . GG - vector containing the newly evaluated gradient
197: */
198: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){
200: AppCtx *user = (AppCtx *) userCtx;
202: PetscInt i,j;
203: PetscInt mx=user->mx, my=user->my;
204: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
205: PetscReal ft=0.0;
206: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
207: PetscReal rhx=mx+1, rhy=my+1;
208: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
209: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
210: PetscReal **g, **x;
211: Vec localX;
214: /* Get local mesh boundaries */
215: DMGetLocalVector(user->dm,&localX);
216: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
217: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
219: /* Scatter ghost points to local vector */
220: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
221: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
223: /* Get pointers to vector data */
224: DMDAVecGetArray(user->dm,localX,(void**)&x);
225: DMDAVecGetArray(user->dm,G,(void**)&g);
227: /* Compute function and gradient over the locally owned part of the mesh */
228: for (j=ys; j<ys+ym; j++){
229: for (i=xs; i< xs+xm; i++){
231: xc = x[j][i];
232: xlt=xrb=xl=xr=xb=xt=xc;
234: if (i==0){ /* left side */
235: xl= user->left[j-ys+1];
236: xlt = user->left[j-ys+2];
237: } else {
238: xl = x[j][i-1];
239: }
241: if (j==0){ /* bottom side */
242: xb=user->bottom[i-xs+1];
243: xrb =user->bottom[i-xs+2];
244: } else {
245: xb = x[j-1][i];
246: }
248: if (i+1 == gxs+gxm){ /* right side */
249: xr=user->right[j-ys+1];
250: xrb = user->right[j-ys];
251: } else {
252: xr = x[j][i+1];
253: }
255: if (j+1==gys+gym){ /* top side */
256: xt=user->top[i-xs+1];
257: xlt = user->top[i-xs];
258: }else {
259: xt = x[j+1][i];
260: }
262: if (i>gxs && j+1<gys+gym){
263: xlt = x[j+1][i-1];
264: }
265: if (j>gys && i+1<gxs+gxm){
266: xrb = x[j-1][i+1];
267: }
269: d1 = (xc-xl);
270: d2 = (xc-xr);
271: d3 = (xc-xt);
272: d4 = (xc-xb);
273: d5 = (xr-xrb);
274: d6 = (xrb-xb);
275: d7 = (xlt-xl);
276: d8 = (xt-xlt);
278: df1dxc = d1*hydhx;
279: df2dxc = ( d1*hydhx + d4*hxdhy );
280: df3dxc = d3*hxdhy;
281: df4dxc = ( d2*hydhx + d3*hxdhy );
282: df5dxc = d2*hydhx;
283: df6dxc = d4*hxdhy;
285: d1 *= rhx;
286: d2 *= rhx;
287: d3 *= rhy;
288: d4 *= rhy;
289: d5 *= rhy;
290: d6 *= rhx;
291: d7 *= rhy;
292: d8 *= rhx;
294: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
295: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
296: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
297: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
298: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
299: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
301: ft = ft + (f2 + f4);
303: df1dxc /= f1;
304: df2dxc /= f2;
305: df3dxc /= f3;
306: df4dxc /= f4;
307: df5dxc /= f5;
308: df6dxc /= f6;
310: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
312: }
313: }
315: /* Compute triangular areas along the border of the domain. */
316: if (xs==0){ /* left side */
317: for (j=ys; j<ys+ym; j++){
318: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
319: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
320: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
321: }
322: }
323: if (ys==0){ /* bottom side */
324: for (i=xs; i<xs+xm; i++){
325: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
326: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
327: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
328: }
329: }
331: if (xs+xm==mx){ /* right side */
332: for (j=ys; j< ys+ym; j++){
333: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
334: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
335: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
336: }
337: }
338: if (ys+ym==my){ /* top side */
339: for (i=xs; i<xs+xm; i++){
340: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
341: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
342: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
343: }
344: }
346: if (ys==0 && xs==0){
347: d1=(user->left[0]-user->left[1])/hy;
348: d2=(user->bottom[0]-user->bottom[1])*rhx;
349: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
350: }
351: if (ys+ym == my && xs+xm == mx){
352: d1=(user->right[ym+1] - user->right[ym])*rhy;
353: d2=(user->top[xm+1] - user->top[xm])*rhx;
354: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
355: }
357: ft=ft*area;
358: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
360: /* Restore vectors */
361: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
362: DMDAVecRestoreArray(user->dm,G,(void**)&g);
364: /* Scatter values to global vector */
365: DMRestoreLocalVector(user->dm,&localX);
366: PetscLogFlops(67*xm*ym);
367: return(0);
368: }
370: /* ------------------------------------------------------------------- */
371: /*
372: FormHessian - Evaluates Hessian matrix.
374: Input Parameters:
375: . tao - the Tao context
376: . x - input vector
377: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
379: Output Parameters:
380: . H - Hessian matrix
381: . Hpre - optionally different preconditioning matrix
382: . flg - flag indicating matrix structure
384: */
385: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
386: {
388: AppCtx *user = (AppCtx *) ptr;
391: /* Evaluate the Hessian entries*/
392: QuadraticH(user,X,H);
393: return(0);
394: }
396: /* ------------------------------------------------------------------- */
397: /*
398: QuadraticH - Evaluates Hessian matrix.
400: Input Parameters:
401: . user - user-defined context, as set by TaoSetHessian()
402: . X - input vector
404: Output Parameter:
405: . H - Hessian matrix
406: */
407: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
408: {
409: PetscErrorCode ierr;
410: PetscInt i,j,k;
411: PetscInt mx=user->mx, my=user->my;
412: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
413: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
414: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
415: PetscReal hl,hr,ht,hb,hc,htl,hbr;
416: PetscReal **x, v[7];
417: MatStencil col[7],row;
418: Vec localX;
419: PetscBool assembled;
422: /* Get local mesh boundaries */
423: DMGetLocalVector(user->dm,&localX);
425: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
426: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
428: /* Scatter ghost points to local vector */
429: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
430: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
432: /* Get pointers to vector data */
433: DMDAVecGetArray(user->dm,localX,(void**)&x);
435: /* Initialize matrix entries to zero */
436: MatAssembled(Hessian,&assembled);
437: if (assembled){MatZeroEntries(Hessian);}
440: /* Set various matrix options */
441: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
443: /* Compute Hessian over the locally owned part of the mesh */
445: for (j=ys; j<ys+ym; j++){
447: for (i=xs; i< xs+xm; i++){
449: xc = x[j][i];
450: xlt=xrb=xl=xr=xb=xt=xc;
452: /* Left side */
453: if (i==0){
454: xl = user->left[j-ys+1];
455: xlt = user->left[j-ys+2];
456: } else {
457: xl = x[j][i-1];
458: }
460: if (j==0){
461: xb = user->bottom[i-xs+1];
462: xrb = user->bottom[i-xs+2];
463: } else {
464: xb = x[j-1][i];
465: }
467: if (i+1 == mx){
468: xr = user->right[j-ys+1];
469: xrb = user->right[j-ys];
470: } else {
471: xr = x[j][i+1];
472: }
474: if (j+1==my){
475: xt = user->top[i-xs+1];
476: xlt = user->top[i-xs];
477: }else {
478: xt = x[j+1][i];
479: }
481: if (i>0 && j+1<my){
482: xlt = x[j+1][i-1];
483: }
484: if (j>0 && i+1<mx){
485: xrb = x[j-1][i+1];
486: }
489: d1 = (xc-xl)/hx;
490: d2 = (xc-xr)/hx;
491: d3 = (xc-xt)/hy;
492: d4 = (xc-xb)/hy;
493: d5 = (xrb-xr)/hy;
494: d6 = (xrb-xb)/hx;
495: d7 = (xlt-xl)/hy;
496: d8 = (xlt-xt)/hx;
498: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
499: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
500: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
501: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
502: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
503: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
506: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
507: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
508: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
509: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
510: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
511: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
512: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
513: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
515: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
516: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
518: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
519: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
520: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
521: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
523: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
525: row.j = j; row.i = i;
526: k=0;
527: if (j>0){
528: v[k]=hb;
529: col[k].j = j - 1; col[k].i = i;
530: k++;
531: }
533: if (j>0 && i < mx -1){
534: v[k]=hbr;
535: col[k].j = j - 1; col[k].i = i+1;
536: k++;
537: }
539: if (i>0){
540: v[k]= hl;
541: col[k].j = j; col[k].i = i-1;
542: k++;
543: }
545: v[k]= hc;
546: col[k].j = j; col[k].i = i;
547: k++;
549: if (i < mx-1 ){
550: v[k]= hr;
551: col[k].j = j; col[k].i = i+1;
552: k++;
553: }
555: if (i>0 && j < my-1 ){
556: v[k]= htl;
557: col[k].j = j+1; col[k].i = i-1;
558: k++;
559: }
561: if (j < my-1 ){
562: v[k]= ht;
563: col[k].j = j+1; col[k].i = i;
564: k++;
565: }
567: /*
568: Set matrix values using local numbering, which was defined
569: earlier, in the main routine.
570: */
571: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
572: }
573: }
575: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
576: DMRestoreLocalVector(user->dm,&localX);
578: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
579: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
581: PetscLogFlops(199*xm*ym);
582: return(0);
583: }
585: /* ------------------------------------------------------------------- */
586: /*
587: MSA_BoundaryConditions - Calculates the boundary conditions for
588: the region.
590: Input Parameter:
591: . user - user-defined application context
593: Output Parameter:
594: . user - user-defined application context
595: */
596: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
597: {
599: PetscInt i,j,k,limit=0,maxits=5;
600: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
601: PetscInt mx=user->mx,my=user->my;
602: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
603: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
604: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
605: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
606: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
607: PetscReal *boundary;
608: PetscBool flg;
611: /* Get local mesh boundaries */
612: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
613: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
615: bsize=xm+2;
616: lsize=ym+2;
617: rsize=ym+2;
618: tsize=xm+2;
620: PetscMalloc1(bsize,&user->bottom);
621: PetscMalloc1(tsize,&user->top);
622: PetscMalloc1(lsize,&user->left);
623: PetscMalloc1(rsize,&user->right);
625: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
627: for (j=0; j<4; j++){
628: if (j==0){
629: yt=b;
630: xt=l+hx*xs;
631: limit=bsize;
632: boundary=user->bottom;
633: } else if (j==1){
634: yt=t;
635: xt=l+hx*xs;
636: limit=tsize;
637: boundary=user->top;
638: } else if (j==2){
639: yt=b+hy*ys;
640: xt=l;
641: limit=lsize;
642: boundary=user->left;
643: } else { /* if (j==3) */
644: yt=b+hy*ys;
645: xt=r;
646: limit=rsize;
647: boundary=user->right;
648: }
650: for (i=0; i<limit; i++){
651: u1=xt;
652: u2=-yt;
653: for (k=0; k<maxits; k++){
654: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
655: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
656: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
657: if (fnorm <= tol) break;
658: njac11=one+u2*u2-u1*u1;
659: njac12=two*u1*u2;
660: njac21=-two*u1*u2;
661: njac22=-one - u1*u1 + u2*u2;
662: det = njac11*njac22-njac21*njac12;
663: u1 = u1-(njac22*nf1-njac12*nf2)/det;
664: u2 = u2-(njac11*nf2-njac21*nf1)/det;
665: }
667: boundary[i]=u1*u1-u2*u2;
668: if (j==0 || j==1) {
669: xt=xt+hx;
670: } else { /* if (j==2 || j==3) */
671: yt=yt+hy;
672: }
674: }
676: }
678: /* Scale the boundary if desired */
679: if (1==1){
680: PetscReal scl = 1.0;
682: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
683: if (flg){
684: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
685: }
687: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
688: if (flg){
689: for (i=0;i<tsize;i++) user->top[i]*=scl;
690: }
692: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
693: if (flg){
694: for (i=0;i<rsize;i++) user->right[i]*=scl;
695: }
697: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
698: if (flg){
699: for (i=0;i<lsize;i++) user->left[i]*=scl;
700: }
701: }
702: return(0);
703: }
705: /* ------------------------------------------------------------------- */
706: /*
707: MSA_InitialPoint - Calculates the initial guess in one of three ways.
709: Input Parameters:
710: . user - user-defined application context
711: . X - vector for initial guess
713: Output Parameters:
714: . X - newly computed initial guess
715: */
716: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
717: {
719: PetscInt start2=-1,i,j;
720: PetscReal start1=0;
721: PetscBool flg1,flg2;
724: PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
725: PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);
727: if (flg1){ /* The zero vector is reasonable */
729: VecSet(X, start1);
731: } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */
732: PetscRandom rctx; PetscReal np5=-0.5;
734: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
735: VecSetRandom(X, rctx);
736: PetscRandomDestroy(&rctx);
737: VecShift(X, np5);
739: } else { /* Take an average of the boundary conditions */
740: PetscInt xs,xm,ys,ym;
741: PetscInt mx=user->mx,my=user->my;
742: PetscReal **x;
744: /* Get local mesh boundaries */
745: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
747: /* Get pointers to vector data */
748: DMDAVecGetArray(user->dm,X,(void**)&x);
750: /* Perform local computations */
751: for (j=ys; j<ys+ym; j++){
752: for (i=xs; i< xs+xm; i++){
753: x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
754: }
755: }
756: DMDAVecRestoreArray(user->dm,X,(void**)&x);
757: PetscLogFlops(9*xm*ym);
758: }
759: return(0);
760: }
762: /*-----------------------------------------------------------------------*/
763: PetscErrorCode My_Monitor(Tao tao, void *ctx)
764: {
766: Vec X;
769: TaoGetSolutionVector(tao,&X);
770: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
771: return(0);
772: }