Actual source code: minsurf2.c
1: /* Program usage: mpirun -np <proc> minsurf2 [-help] [all TAO options] */
3: /*
4: Include "taosolver.h" so we can use TAO solvers.
5: petscdm.h for distributed array
6: */
7: #include "taosolver.h"
8: #include "petscdm.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: TaoInitialize(); TaoFinalize();
25: Routines: TaoCreate(); TaoSetType();
26: Routines: TaoSetInitialVector();
27: Routines: TaoSetObjectiveAndGradientRoutine();
28: Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
29: Routines: TaoSetMonitor();
30: Routines: TaoSolve(); TaoView();
31: Routines: TaoGetTerminationReason(); TaoDestroy();
32: Processors: n
33: T*/
35: /*
36: User-defined application context - contains data needed by the
37: application-provided call-back routines, FormFunctionGradient()
38: and FormHessian().
39: */
40: typedef struct {
41: PetscInt mx, my; /* discretization in x, y directions */
42: PetscReal *bottom, *top, *left, *right; /* boundary values */
43: DM dm; /* distributed array data structure */
44: Mat H; /* Hessian */
45: } AppCtx;
48: /* -------- User-defined Routines --------- */
50: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
51: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
52: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
53: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal *,Vec,void*);
54: PetscErrorCode FormGradient(TaoSolver,Vec,Vec,void*);
55: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure *,void*);
56: PetscErrorCode My_Monitor(TaoSolver, void *);
60: int main( int argc, char **argv )
61: {
62: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
63: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
64: Vec x; /* solution, gradient vectors */
65: PetscBool flg, viewmat; /* flags */
66: PetscBool fddefault, fdcoloring; /* flags */
67: TaoSolverTerminationReason reason;
68: TaoSolver tao; /* TAO solver context */
69: AppCtx user; /* user-defined work context */
70: ISColoring iscoloring;
71: MatFDColoring matfdcoloring;
73: /* Initialize TAO */
74: PetscInitialize( &argc, &argv,(char *)0,help );
75: TaoInitialize( &argc, &argv,(char *)0,help );
77: /* Specify dimension of the problem */
78: user.mx = 10; user.my = 10;
80: /* Check for any command line arguments that override defaults */
81: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg);
82: PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg);
84: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
85: PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my);
88: /* Let PETSc determine the vector distribution */
89: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
91: /* Create distributed array (DM) to manage parallel grid and vectors */
92: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,
93: DMDA_STENCIL_BOX,
94: user.mx, user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,
95: &user.dm);
96:
98: /* Create TAO solver and set desired solution method.*/
99: TaoCreate(PETSC_COMM_WORLD,&tao);
100: TaoSetType(tao,"tao_cg");
102: /*
103: Extract global vector from DA for the vector of variables -- PETSC routine
104: Compute the initial solution -- application specific, see below
105: Set this vector for use by TAO -- TAO routine
106: */
107: DMCreateGlobalVector(user.dm,&x);
108: MSA_BoundaryConditions(&user);
109: MSA_InitialPoint(&user,x);
110: TaoSetInitialVector(tao,x);
112: /*
113: Initialize the Application context for use in function evaluations -- application specific, see below.
114: Set routines for function and gradient evaluation
115: */
116: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
118: /*
119: Given the command line arguments, calculate the hessian with either the user-
120: provided function FormHessian, or the default finite-difference driven Hessian
121: functions
122: */
123: PetscOptionsHasName(PETSC_NULL,"-tao_fddefault",&fddefault);
124: PetscOptionsHasName(PETSC_NULL,"-tao_fdcoloring",&fdcoloring);
127: /*
128: Create a matrix data structure to store the Hessian and set
129: the Hessian evalution routine.
130: Set the matrix structure to be used for Hessian evalutions
131: */
132: DMGetMatrix(user.dm,MATAIJ,&user.H);
133: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
136: if (fdcoloring) {
137: DMGetColoring(user.dm,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
138:
140: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
141:
143: ISColoringDestroy(&iscoloring);
144: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
145: MatFDColoringSetFromOptions(matfdcoloring);
146:
147: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
149: } else if (fddefault){
150: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)PETSC_NULL);
152: } else {
153: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
154: }
157: /*
158: If my_monitor option is in command line, then use the user-provided
159: monitoring function
160: */
161: PetscOptionsHasName(PETSC_NULL,"-my_monitor",&viewmat);
162: if (viewmat){
163: TaoSetMonitor(tao,My_Monitor,PETSC_NULL,PETSC_NULL);
164: }
166: /* Check for any tao command line options */
167: TaoSetFromOptions(tao);
169: /* SOLVE THE APPLICATION */
170: TaoSolve(tao);
172: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
174: /* Get information on termination */
175: TaoGetTerminationReason(tao,&reason);
176: if (reason <= 0 ){
177: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method \n");
178: }
181: /* Free TAO data structures */
182: TaoDestroy(&tao);
184: /* Free PETSc data structures */
185: VecDestroy(&x);
186: MatDestroy(&user.H);
187: if (fdcoloring) {
188: MatFDColoringDestroy(&matfdcoloring);
189: }
190: PetscFree(user.bottom);
191: PetscFree(user.top);
192: PetscFree(user.left);
193: PetscFree(user.right);
194: DMDestroy(&user.dm);
196: /* Finalize TAO */
197: TaoFinalize();
198: PetscFinalize();
199:
200: return 0;
201: }
205: PetscErrorCode FormGradient(TaoSolver tao, Vec X, Vec G,void *userCtx){
207: PetscReal fcn;
209: FormFunctionGradient(tao,X,&fcn,G,userCtx);
210: return(0);
211: }
213: /* -------------------------------------------------------------------- */
216: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
218: Input Parameters:
219: . tao - the TaoSolver context
220: . XX - input vector
221: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
222:
223: Output Parameters:
224: . fcn - the newly evaluated function
225: . GG - vector containing the newly evaluated gradient
226: */
227: PetscErrorCode FormFunctionGradient(TaoSolver tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){
229: AppCtx * user = (AppCtx *) userCtx;
230: PetscErrorCode ierr;
231: PetscInt i,j;
232: PetscInt mx=user->mx, my=user->my;
233: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
234: PetscReal ft=0.0;
235: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
236: PetscReal rhx=mx+1, rhy=my+1;
237: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
238: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
239: PetscReal **g, **x;
240: Vec localX;
243: /* Get local mesh boundaries */
244: DMGetLocalVector(user->dm,&localX);
246: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
247: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
249: /* Scatter ghost points to local vector */
250: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
251: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
253: /* Get pointers to vector data */
254: DMDAVecGetArray(user->dm,localX,(void**)&x);
255: DMDAVecGetArray(user->dm,G,(void**)&g);
257: /* Compute function and gradient over the locally owned part of the mesh */
258: for (j=ys; j<ys+ym; j++){
259: for (i=xs; i< xs+xm; i++){
260:
261: xc = x[j][i];
262: xlt=xrb=xl=xr=xb=xt=xc;
263:
264: if (i==0){ /* left side */
265: xl= user->left[j-ys+1];
266: xlt = user->left[j-ys+2];
267: } else {
268: xl = x[j][i-1];
269: }
271: if (j==0){ /* bottom side */
272: xb=user->bottom[i-xs+1];
273: xrb =user->bottom[i-xs+2];
274: } else {
275: xb = x[j-1][i];
276: }
277:
278: if (i+1 == gxs+gxm){ /* right side */
279: xr=user->right[j-ys+1];
280: xrb = user->right[j-ys];
281: } else {
282: xr = x[j][i+1];
283: }
285: if (j+1==gys+gym){ /* top side */
286: xt=user->top[i-xs+1];
287: xlt = user->top[i-xs];
288: }else {
289: xt = x[j+1][i];
290: }
292: if (i>gxs && j+1<gys+gym){
293: xlt = x[j+1][i-1];
294: }
295: if (j>gys && i+1<gxs+gxm){
296: xrb = x[j-1][i+1];
297: }
299: d1 = (xc-xl);
300: d2 = (xc-xr);
301: d3 = (xc-xt);
302: d4 = (xc-xb);
303: d5 = (xr-xrb);
304: d6 = (xrb-xb);
305: d7 = (xlt-xl);
306: d8 = (xt-xlt);
307:
308: df1dxc = d1*hydhx;
309: df2dxc = ( d1*hydhx + d4*hxdhy );
310: df3dxc = d3*hxdhy;
311: df4dxc = ( d2*hydhx + d3*hxdhy );
312: df5dxc = d2*hydhx;
313: df6dxc = d4*hxdhy;
315: d1 *= rhx;
316: d2 *= rhx;
317: d3 *= rhy;
318: d4 *= rhy;
319: d5 *= rhy;
320: d6 *= rhx;
321: d7 *= rhy;
322: d8 *= rhx;
324: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
325: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
326: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
327: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
328: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
329: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
330:
331: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
332: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
334: ft = ft + (f2 + f4);
336: df1dxc /= f1;
337: df2dxc /= f2;
338: df3dxc /= f3;
339: df4dxc /= f4;
340: df5dxc /= f5;
341: df6dxc /= f6;
343: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
344:
345: }
346: }
348: /* Compute triangular areas along the border of the domain. */
349: if (xs==0){ /* left side */
350: for (j=ys; j<ys+ym; j++){
351: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
352: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
353: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
354: }
355: }
356: if (ys==0){ /* bottom side */
357: for (i=xs; i<xs+xm; i++){
358: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
359: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
360: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
361: }
362: }
364: if (xs+xm==mx){ /* right side */
365: for (j=ys; j< ys+ym; j++){
366: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
367: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
368: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
369: }
370: }
371: if (ys+ym==my){ /* top side */
372: for (i=xs; i<xs+xm; i++){
373: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
374: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
375: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
376: }
377: }
379: if (ys==0 && xs==0){
380: d1=(user->left[0]-user->left[1])/hy;
381: d2=(user->bottom[0]-user->bottom[1])*rhx;
382: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
383: }
384: if (ys+ym == my && xs+xm == mx){
385: d1=(user->right[ym+1] - user->right[ym])*rhy;
386: d2=(user->top[xm+1] - user->top[xm])*rhx;
387: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
388: }
390: ft=ft*area;
391: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
393: /* Restore vectors */
394: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
395: DMDAVecRestoreArray(user->dm,G,(void**)&g);
397: /* Scatter values to global vector */
398: DMRestoreLocalVector(user->dm,&localX);
400: PetscLogFlops(67*xm*ym);
402: return(0);
403: }
405: /* ------------------------------------------------------------------- */
408: /*
409: FormHessian - Evaluates Hessian matrix.
411: Input Parameters:
412: . tao - the TaoSolver context
413: . x - input vector
414: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
416: Output Parameters:
417: . H - Hessian matrix
418: . Hpre - optionally different preconditioning matrix
419: . flg - flag indicating matrix structure
421: */
422: PetscErrorCode FormHessian(TaoSolver tao,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
423: {
424: PetscErrorCode ierr;
425: AppCtx *user = (AppCtx *) ptr;
428: /* Evaluate the Hessian entries*/
429: QuadraticH(user,X,*H);
432: /* Indicate that this matrix has the same sparsity pattern during
433: successive iterations; setting this flag can save significant work
434: in computing the preconditioner for some methods. */
435: *flg=SAME_NONZERO_PATTERN;
436: return(0);
437: }
439: /* ------------------------------------------------------------------- */
442: /*
443: QuadraticH - Evaluates Hessian matrix.
445: Input Parameters:
446: . user - user-defined context, as set by TaoSetHessian()
447: . X - input vector
449: Output Parameter:
450: . H - Hessian matrix
451: */
452: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
453: {
454: PetscErrorCode ierr;
455: PetscInt i,j,k;
456: PetscInt mx=user->mx, my=user->my;
457: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
458: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
459: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
460: PetscReal hl,hr,ht,hb,hc,htl,hbr;
461: PetscReal **x, v[7];
462: MatStencil col[7],row;
463: Vec localX;
464: PetscBool assembled;
467: /* Get local mesh boundaries */
468: DMGetLocalVector(user->dm,&localX);
470: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
471: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
473: /* Scatter ghost points to local vector */
474: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
475: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
477: /* Get pointers to vector data */
478: DMDAVecGetArray(user->dm,localX,(void**)&x);
480: /* Initialize matrix entries to zero */
481: MatAssembled(Hessian,&assembled);
482: if (assembled){MatZeroEntries(Hessian); }
485: /* Set various matrix options */
486: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
488: /* Compute Hessian over the locally owned part of the mesh */
490: for (j=ys; j<ys+ym; j++){
491:
492: for (i=xs; i< xs+xm; i++){
494: xc = x[j][i];
495: xlt=xrb=xl=xr=xb=xt=xc;
497: /* Left side */
498: if (i==0){
499: xl = user->left[j-ys+1];
500: xlt = user->left[j-ys+2];
501: } else {
502: xl = x[j][i-1];
503: }
504:
505: if (j==0){
506: xb = user->bottom[i-xs+1];
507: xrb = user->bottom[i-xs+2];
508: } else {
509: xb = x[j-1][i];
510: }
511:
512: if (i+1 == mx){
513: xr = user->right[j-ys+1];
514: xrb = user->right[j-ys];
515: } else {
516: xr = x[j][i+1];
517: }
519: if (j+1==my){
520: xt = user->top[i-xs+1];
521: xlt = user->top[i-xs];
522: }else {
523: xt = x[j+1][i];
524: }
526: if (i>0 && j+1<my){
527: xlt = x[j+1][i-1];
528: }
529: if (j>0 && i+1<mx){
530: xrb = x[j-1][i+1];
531: }
534: d1 = (xc-xl)/hx;
535: d2 = (xc-xr)/hx;
536: d3 = (xc-xt)/hy;
537: d4 = (xc-xb)/hy;
538: d5 = (xrb-xr)/hy;
539: d6 = (xrb-xb)/hx;
540: d7 = (xlt-xl)/hy;
541: d8 = (xlt-xt)/hx;
542:
543: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
544: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
545: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
546: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
547: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
548: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
551: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
552: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
553: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
554: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
555: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
556: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
557: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
558: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
560: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
561: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
563: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
564: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
565: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
566: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
568: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
570: row.j = j; row.i = i;
571: k=0;
572: if (j>0){
573: v[k]=hb;
574: col[k].j = j - 1; col[k].i = i;
575: k++;
576: }
577:
578: if (j>0 && i < mx -1){
579: v[k]=hbr;
580: col[k].j = j - 1; col[k].i = i+1;
581: k++;
582: }
583:
584: if (i>0){
585: v[k]= hl;
586: col[k].j = j; col[k].i = i-1;
587: k++;
588: }
589:
590: v[k]= hc;
591: col[k].j = j; col[k].i = i;
592: k++;
593:
594: if (i < mx-1 ){
595: v[k]= hr;
596: col[k].j = j; col[k].i = i+1;
597: k++;
598: }
599:
600: if (i>0 && j < my-1 ){
601: v[k]= htl;
602: col[k].j = j+1; col[k].i = i-1;
603: k++;
604: }
605:
606: if (j < my-1 ){
607: v[k]= ht;
608: col[k].j = j+1; col[k].i = i;
609: k++;
610: }
611:
612: /*
613: Set matrix values using local numbering, which was defined
614: earlier, in the main routine.
615: */
616: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
617:
618:
619: }
620: }
621:
622: /* Restore vectors */
623: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
625: DMRestoreLocalVector(user->dm,&localX);
627: /* Assemble the matrix */
628: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
629: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
631: PetscLogFlops(199*xm*ym);
632: return(0);
633: }
635: /* ------------------------------------------------------------------- */
638: /*
639: MSA_BoundaryConditions - Calculates the boundary conditions for
640: the region.
642: Input Parameter:
643: . user - user-defined application context
645: Output Parameter:
646: . user - user-defined application context
647: */
648: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
649: {
650: PetscErrorCode ierr;
651: PetscInt i,j,k,limit=0,maxits=5;
652: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
653: PetscInt mx=user->mx,my=user->my;
654: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
655: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
656: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
657: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
658: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
659: PetscReal *boundary;
660: PetscBool flg;
663: /* Get local mesh boundaries */
664: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
665: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
667: bsize=xm+2;
668: lsize=ym+2;
669: rsize=ym+2;
670: tsize=xm+2;
672: PetscMalloc(bsize*sizeof(PetscReal),&user->bottom);
673: PetscMalloc(tsize*sizeof(PetscReal),&user->top);
674: PetscMalloc(lsize*sizeof(PetscReal),&user->left);
675: PetscMalloc(rsize*sizeof(PetscReal),&user->right);
677: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
679: for (j=0; j<4; j++){
680: if (j==0){
681: yt=b;
682: xt=l+hx*xs;
683: limit=bsize;
684: boundary=user->bottom;
685: } else if (j==1){
686: yt=t;
687: xt=l+hx*xs;
688: limit=tsize;
689: boundary=user->top;
690: } else if (j==2){
691: yt=b+hy*ys;
692: xt=l;
693: limit=lsize;
694: boundary=user->left;
695: } else { /* if (j==3) */
696: yt=b+hy*ys;
697: xt=r;
698: limit=rsize;
699: boundary=user->right;
700: }
702: for (i=0; i<limit; i++){
703: u1=xt;
704: u2=-yt;
705: for (k=0; k<maxits; k++){
706: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
707: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
708: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
709: if (fnorm <= tol) break;
710: njac11=one+u2*u2-u1*u1;
711: njac12=two*u1*u2;
712: njac21=-two*u1*u2;
713: njac22=-one - u1*u1 + u2*u2;
714: det = njac11*njac22-njac21*njac12;
715: u1 = u1-(njac22*nf1-njac12*nf2)/det;
716: u2 = u2-(njac11*nf2-njac21*nf1)/det;
717: }
719: boundary[i]=u1*u1-u2*u2;
720: if (j==0 || j==1) {
721: xt=xt+hx;
722: } else { /* if (j==2 || j==3) */
723: yt=yt+hy;
724: }
725:
726: }
728: }
730: /* Scale the boundary if desired */
731: if (1==1){
732: PetscReal scl = 1.0;
734: PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
735:
736: if (flg){
737: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
738: }
740: PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
741:
742: if (flg){
743: for (i=0;i<tsize;i++) user->top[i]*=scl;
744: }
746: PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
747:
748: if (flg){
749: for (i=0;i<rsize;i++) user->right[i]*=scl;
750: }
752: PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
753:
754: if (flg){
755: for (i=0;i<lsize;i++) user->left[i]*=scl;
756: }
757: }
758:
759: return(0);
760: }
762: /* ------------------------------------------------------------------- */
765: /*
766: MSA_InitialPoint - Calculates the initial guess in one of three ways.
768: Input Parameters:
769: . user - user-defined application context
770: . X - vector for initial guess
772: Output Parameters:
773: . X - newly computed initial guess
774: */
775: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
776: {
777: PetscErrorCode ierr;
778: PetscInt start2=-1,i,j;
779: PetscReal start1=0;
780: PetscBool flg1,flg2;
783: PetscOptionsGetReal(PETSC_NULL,"-start",&start1,&flg1);
784: PetscOptionsGetInt(PETSC_NULL,"-random",&start2,&flg2);
786: if (flg1){ /* The zero vector is reasonable */
787:
788: VecSet(X, start1);
790: } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */
792: PetscRandom rctx; PetscReal np5=-0.5;
794: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
795: PetscRandomSetType(rctx,PETSCRAND);
796:
797: for (i=0; i<start2; i++){
798: VecSetRandom(X, rctx);
799: }
800: PetscRandomDestroy(&rctx);
801: VecShift(X, np5);
803: } else { /* Take an average of the boundary conditions */
805: PetscInt xs,xm,ys,ym;
806: PetscInt mx=user->mx,my=user->my;
807: PetscReal **x;
808:
809: /* Get local mesh boundaries */
810: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
811:
812: /* Get pointers to vector data */
813: DMDAVecGetArray(user->dm,X,(void**)&x);
815: /* Perform local computations */
816: for (j=ys; j<ys+ym; j++){
817: for (i=xs; i< xs+xm; i++){
818: x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+
819: ((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
820: }
821: }
822:
823: /* Restore vectors */
824: DMDAVecRestoreArray(user->dm,X,(void**)&x);
826: PetscLogFlops(9*xm*ym);
827:
828: }
829: return(0);
830: }
832: /*-----------------------------------------------------------------------*/
835: PetscErrorCode My_Monitor(TaoSolver tao, void *ctx){
837: Vec X;
839: TaoGetSolutionVector(tao,&X);
840: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
841: return(0);
842: }