Actual source code: minsurf2.c
petsc-3.7.7 2017-09-25
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 *);
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: Tao tao; /* TAO solver context */
67: AppCtx user; /* user-defined work context */
68: ISColoring iscoloring;
69: MatFDColoring matfdcoloring;
71: /* Initialize TAO */
72: PetscInitialize( &argc, &argv,(char *)0,help );
74: /* Specify dimension of the problem */
75: user.mx = 10; user.my = 10;
77: /* Check for any command line arguments that override defaults */
78: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
79: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
81: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
82: PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
85: /* Let PETSc determine the vector distribution */
86: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
88: /* Create distributed array (DM) to manage parallel grid and vectors */
89: 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);
92: /* Create TAO solver and set desired solution method.*/
93: TaoCreate(PETSC_COMM_WORLD,&tao);
94: TaoSetType(tao,TAOCG);
96: /*
97: Extract global vector from DA for the vector of variables -- PETSC routine
98: Compute the initial solution -- application specific, see below
99: Set this vector for use by TAO -- TAO routine
100: */
101: DMCreateGlobalVector(user.dm,&x);
102: MSA_BoundaryConditions(&user);
103: MSA_InitialPoint(&user,x);
104: TaoSetInitialVector(tao,x);
106: /*
107: Initialize the Application context for use in function evaluations -- application specific, see below.
108: Set routines for function and gradient evaluation
109: */
110: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
112: /*
113: Given the command line arguments, calculate the hessian with either the user-
114: provided function FormHessian, or the default finite-difference driven Hessian
115: functions
116: */
117: PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault);
118: PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);
120: /*
121: Create a matrix data structure to store the Hessian and set
122: the Hessian evalution routine.
123: Set the matrix structure to be used for Hessian evalutions
124: */
125: DMCreateMatrix(user.dm,&user.H);
126: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
129: if (fdcoloring) {
130: DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
131: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
132: ISColoringDestroy(&iscoloring);
133: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
134: MatFDColoringSetFromOptions(matfdcoloring);
135: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
136: } else if (fddefault){
137: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
138: } else {
139: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
140: }
142: /*
143: If my_monitor option is in command line, then use the user-provided
144: monitoring function
145: */
146: PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
147: if (viewmat){
148: TaoSetMonitor(tao,My_Monitor,NULL,NULL);
149: }
151: /* Check for any tao command line options */
152: TaoSetFromOptions(tao);
154: /* SOLVE THE APPLICATION */
155: TaoSolve(tao);
157: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
159: /* Free TAO data structures */
160: TaoDestroy(&tao);
162: /* Free PETSc data structures */
163: VecDestroy(&x);
164: MatDestroy(&user.H);
165: if (fdcoloring) {
166: MatFDColoringDestroy(&matfdcoloring);
167: }
168: PetscFree(user.bottom);
169: PetscFree(user.top);
170: PetscFree(user.left);
171: PetscFree(user.right);
172: DMDestroy(&user.dm);
173: PetscFinalize();
174: return 0;
175: }
179: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
180: {
182: PetscReal fcn;
185: FormFunctionGradient(tao,X,&fcn,G,userCtx);
186: return(0);
187: }
189: /* -------------------------------------------------------------------- */
192: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
194: Input Parameters:
195: . tao - the Tao context
196: . XX - input vector
197: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
199: Output Parameters:
200: . fcn - the newly evaluated function
201: . GG - vector containing the newly evaluated gradient
202: */
203: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){
205: AppCtx *user = (AppCtx *) userCtx;
207: PetscInt i,j;
208: PetscInt mx=user->mx, my=user->my;
209: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
210: PetscReal ft=0.0;
211: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
212: PetscReal rhx=mx+1, rhy=my+1;
213: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
214: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
215: PetscReal **g, **x;
216: Vec localX;
219: /* Get local mesh boundaries */
220: DMGetLocalVector(user->dm,&localX);
221: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
222: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
224: /* Scatter ghost points to local vector */
225: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
226: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
228: /* Get pointers to vector data */
229: DMDAVecGetArray(user->dm,localX,(void**)&x);
230: DMDAVecGetArray(user->dm,G,(void**)&g);
232: /* Compute function and gradient over the locally owned part of the mesh */
233: for (j=ys; j<ys+ym; j++){
234: for (i=xs; i< xs+xm; i++){
236: xc = x[j][i];
237: xlt=xrb=xl=xr=xb=xt=xc;
239: if (i==0){ /* left side */
240: xl= user->left[j-ys+1];
241: xlt = user->left[j-ys+2];
242: } else {
243: xl = x[j][i-1];
244: }
246: if (j==0){ /* bottom side */
247: xb=user->bottom[i-xs+1];
248: xrb =user->bottom[i-xs+2];
249: } else {
250: xb = x[j-1][i];
251: }
253: if (i+1 == gxs+gxm){ /* right side */
254: xr=user->right[j-ys+1];
255: xrb = user->right[j-ys];
256: } else {
257: xr = x[j][i+1];
258: }
260: if (j+1==gys+gym){ /* top side */
261: xt=user->top[i-xs+1];
262: xlt = user->top[i-xs];
263: }else {
264: xt = x[j+1][i];
265: }
267: if (i>gxs && j+1<gys+gym){
268: xlt = x[j+1][i-1];
269: }
270: if (j>gys && i+1<gxs+gxm){
271: xrb = x[j-1][i+1];
272: }
274: d1 = (xc-xl);
275: d2 = (xc-xr);
276: d3 = (xc-xt);
277: d4 = (xc-xb);
278: d5 = (xr-xrb);
279: d6 = (xrb-xb);
280: d7 = (xlt-xl);
281: d8 = (xt-xlt);
283: df1dxc = d1*hydhx;
284: df2dxc = ( d1*hydhx + d4*hxdhy );
285: df3dxc = d3*hxdhy;
286: df4dxc = ( d2*hydhx + d3*hxdhy );
287: df5dxc = d2*hydhx;
288: df6dxc = d4*hxdhy;
290: d1 *= rhx;
291: d2 *= rhx;
292: d3 *= rhy;
293: d4 *= rhy;
294: d5 *= rhy;
295: d6 *= rhx;
296: d7 *= rhy;
297: d8 *= rhx;
299: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
300: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
301: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
302: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
303: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
304: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
306: ft = ft + (f2 + f4);
308: df1dxc /= f1;
309: df2dxc /= f2;
310: df3dxc /= f3;
311: df4dxc /= f4;
312: df5dxc /= f5;
313: df6dxc /= f6;
315: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
317: }
318: }
320: /* Compute triangular areas along the border of the domain. */
321: if (xs==0){ /* left side */
322: for (j=ys; j<ys+ym; j++){
323: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
324: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
325: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
326: }
327: }
328: if (ys==0){ /* bottom side */
329: for (i=xs; i<xs+xm; i++){
330: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
331: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
332: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
333: }
334: }
336: if (xs+xm==mx){ /* right side */
337: for (j=ys; j< ys+ym; j++){
338: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
339: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
340: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
341: }
342: }
343: if (ys+ym==my){ /* top side */
344: for (i=xs; i<xs+xm; i++){
345: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
346: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
347: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
348: }
349: }
351: if (ys==0 && xs==0){
352: d1=(user->left[0]-user->left[1])/hy;
353: d2=(user->bottom[0]-user->bottom[1])*rhx;
354: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
355: }
356: if (ys+ym == my && xs+xm == mx){
357: d1=(user->right[ym+1] - user->right[ym])*rhy;
358: d2=(user->top[xm+1] - user->top[xm])*rhx;
359: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
360: }
362: ft=ft*area;
363: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
365: /* Restore vectors */
366: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
367: DMDAVecRestoreArray(user->dm,G,(void**)&g);
369: /* Scatter values to global vector */
370: DMRestoreLocalVector(user->dm,&localX);
371: PetscLogFlops(67*xm*ym);
372: return(0);
373: }
375: /* ------------------------------------------------------------------- */
378: /*
379: FormHessian - Evaluates Hessian matrix.
381: Input Parameters:
382: . tao - the Tao context
383: . x - input vector
384: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
386: Output Parameters:
387: . H - Hessian matrix
388: . Hpre - optionally different preconditioning matrix
389: . flg - flag indicating matrix structure
391: */
392: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
393: {
395: AppCtx *user = (AppCtx *) ptr;
398: /* Evaluate the Hessian entries*/
399: QuadraticH(user,X,H);
400: return(0);
401: }
403: /* ------------------------------------------------------------------- */
406: /*
407: QuadraticH - Evaluates Hessian matrix.
409: Input Parameters:
410: . user - user-defined context, as set by TaoSetHessian()
411: . X - input vector
413: Output Parameter:
414: . H - Hessian matrix
415: */
416: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
417: {
418: PetscErrorCode ierr;
419: PetscInt i,j,k;
420: PetscInt mx=user->mx, my=user->my;
421: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
422: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
423: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
424: PetscReal hl,hr,ht,hb,hc,htl,hbr;
425: PetscReal **x, v[7];
426: MatStencil col[7],row;
427: Vec localX;
428: PetscBool assembled;
431: /* Get local mesh boundaries */
432: DMGetLocalVector(user->dm,&localX);
434: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
435: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
437: /* Scatter ghost points to local vector */
438: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
439: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
441: /* Get pointers to vector data */
442: DMDAVecGetArray(user->dm,localX,(void**)&x);
444: /* Initialize matrix entries to zero */
445: MatAssembled(Hessian,&assembled);
446: if (assembled){MatZeroEntries(Hessian);}
449: /* Set various matrix options */
450: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
452: /* Compute Hessian over the locally owned part of the mesh */
454: for (j=ys; j<ys+ym; j++){
456: for (i=xs; i< xs+xm; i++){
458: xc = x[j][i];
459: xlt=xrb=xl=xr=xb=xt=xc;
461: /* Left side */
462: if (i==0){
463: xl = user->left[j-ys+1];
464: xlt = user->left[j-ys+2];
465: } else {
466: xl = x[j][i-1];
467: }
469: if (j==0){
470: xb = user->bottom[i-xs+1];
471: xrb = user->bottom[i-xs+2];
472: } else {
473: xb = x[j-1][i];
474: }
476: if (i+1 == mx){
477: xr = user->right[j-ys+1];
478: xrb = user->right[j-ys];
479: } else {
480: xr = x[j][i+1];
481: }
483: if (j+1==my){
484: xt = user->top[i-xs+1];
485: xlt = user->top[i-xs];
486: }else {
487: xt = x[j+1][i];
488: }
490: if (i>0 && j+1<my){
491: xlt = x[j+1][i-1];
492: }
493: if (j>0 && i+1<mx){
494: xrb = x[j-1][i+1];
495: }
498: d1 = (xc-xl)/hx;
499: d2 = (xc-xr)/hx;
500: d3 = (xc-xt)/hy;
501: d4 = (xc-xb)/hy;
502: d5 = (xrb-xr)/hy;
503: d6 = (xrb-xb)/hx;
504: d7 = (xlt-xl)/hy;
505: d8 = (xlt-xt)/hx;
507: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
508: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
509: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
510: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
511: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
512: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
515: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
516: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
517: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
518: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
519: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
520: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
521: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
522: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
524: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
525: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
527: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
528: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
529: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
530: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
532: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
534: row.j = j; row.i = i;
535: k=0;
536: if (j>0){
537: v[k]=hb;
538: col[k].j = j - 1; col[k].i = i;
539: k++;
540: }
542: if (j>0 && i < mx -1){
543: v[k]=hbr;
544: col[k].j = j - 1; col[k].i = i+1;
545: k++;
546: }
548: if (i>0){
549: v[k]= hl;
550: col[k].j = j; col[k].i = i-1;
551: k++;
552: }
554: v[k]= hc;
555: col[k].j = j; col[k].i = i;
556: k++;
558: if (i < mx-1 ){
559: v[k]= hr;
560: col[k].j = j; col[k].i = i+1;
561: k++;
562: }
564: if (i>0 && j < my-1 ){
565: v[k]= htl;
566: col[k].j = j+1; col[k].i = i-1;
567: k++;
568: }
570: if (j < my-1 ){
571: v[k]= ht;
572: col[k].j = j+1; col[k].i = i;
573: k++;
574: }
576: /*
577: Set matrix values using local numbering, which was defined
578: earlier, in the main routine.
579: */
580: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
581: }
582: }
584: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
585: DMRestoreLocalVector(user->dm,&localX);
587: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
588: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
590: PetscLogFlops(199*xm*ym);
591: return(0);
592: }
594: /* ------------------------------------------------------------------- */
597: /*
598: MSA_BoundaryConditions - Calculates the boundary conditions for
599: the region.
601: Input Parameter:
602: . user - user-defined application context
604: Output Parameter:
605: . user - user-defined application context
606: */
607: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
608: {
610: PetscInt i,j,k,limit=0,maxits=5;
611: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
612: PetscInt mx=user->mx,my=user->my;
613: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
614: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
615: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
616: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
617: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
618: PetscReal *boundary;
619: PetscBool flg;
622: /* Get local mesh boundaries */
623: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
624: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
626: bsize=xm+2;
627: lsize=ym+2;
628: rsize=ym+2;
629: tsize=xm+2;
631: PetscMalloc1(bsize,&user->bottom);
632: PetscMalloc1(tsize,&user->top);
633: PetscMalloc1(lsize,&user->left);
634: PetscMalloc1(rsize,&user->right);
636: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
638: for (j=0; j<4; j++){
639: if (j==0){
640: yt=b;
641: xt=l+hx*xs;
642: limit=bsize;
643: boundary=user->bottom;
644: } else if (j==1){
645: yt=t;
646: xt=l+hx*xs;
647: limit=tsize;
648: boundary=user->top;
649: } else if (j==2){
650: yt=b+hy*ys;
651: xt=l;
652: limit=lsize;
653: boundary=user->left;
654: } else { /* if (j==3) */
655: yt=b+hy*ys;
656: xt=r;
657: limit=rsize;
658: boundary=user->right;
659: }
661: for (i=0; i<limit; i++){
662: u1=xt;
663: u2=-yt;
664: for (k=0; k<maxits; k++){
665: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
666: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
667: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
668: if (fnorm <= tol) break;
669: njac11=one+u2*u2-u1*u1;
670: njac12=two*u1*u2;
671: njac21=-two*u1*u2;
672: njac22=-one - u1*u1 + u2*u2;
673: det = njac11*njac22-njac21*njac12;
674: u1 = u1-(njac22*nf1-njac12*nf2)/det;
675: u2 = u2-(njac11*nf2-njac21*nf1)/det;
676: }
678: boundary[i]=u1*u1-u2*u2;
679: if (j==0 || j==1) {
680: xt=xt+hx;
681: } else { /* if (j==2 || j==3) */
682: yt=yt+hy;
683: }
685: }
687: }
689: /* Scale the boundary if desired */
690: if (1==1){
691: PetscReal scl = 1.0;
693: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
694: if (flg){
695: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
696: }
698: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
699: if (flg){
700: for (i=0;i<tsize;i++) user->top[i]*=scl;
701: }
703: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
704: if (flg){
705: for (i=0;i<rsize;i++) user->right[i]*=scl;
706: }
708: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
709: if (flg){
710: for (i=0;i<lsize;i++) user->left[i]*=scl;
711: }
712: }
713: return(0);
714: }
716: /* ------------------------------------------------------------------- */
719: /*
720: MSA_InitialPoint - Calculates the initial guess in one of three ways.
722: Input Parameters:
723: . user - user-defined application context
724: . X - vector for initial guess
726: Output Parameters:
727: . X - newly computed initial guess
728: */
729: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
730: {
732: PetscInt start2=-1,i,j;
733: PetscReal start1=0;
734: PetscBool flg1,flg2;
737: PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
738: PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);
740: if (flg1){ /* The zero vector is reasonable */
742: VecSet(X, start1);
744: } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */
745: PetscRandom rctx; PetscReal np5=-0.5;
747: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
748: VecSetRandom(X, rctx);
749: PetscRandomDestroy(&rctx);
750: VecShift(X, np5);
752: } else { /* Take an average of the boundary conditions */
753: PetscInt xs,xm,ys,ym;
754: PetscInt mx=user->mx,my=user->my;
755: PetscReal **x;
757: /* Get local mesh boundaries */
758: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
760: /* Get pointers to vector data */
761: DMDAVecGetArray(user->dm,X,(void**)&x);
763: /* Perform local computations */
764: for (j=ys; j<ys+ym; j++){
765: for (i=xs; i< xs+xm; i++){
766: 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;
767: }
768: }
769: DMDAVecRestoreArray(user->dm,X,(void**)&x);
770: PetscLogFlops(9*xm*ym);
771: }
772: return(0);
773: }
775: /*-----------------------------------------------------------------------*/
778: PetscErrorCode My_Monitor(Tao tao, void *ctx)
779: {
781: Vec X;
784: TaoGetSolutionVector(tao,&X);
785: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
786: return(0);
787: }