Actual source code: minsurf2.c
petsc-3.13.6 2020-09-29
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*/
36: /*
37: User-defined Section 1.5 Writing Application Codes with PETSc context - contains data needed by the
38: Section 1.5 Writing Application Codes with PETSc-provided call-back routines, FormFunctionGradient()
39: and FormHessian().
40: */
41: typedef struct {
42: PetscInt mx, my; /* discretization in x, y directions */
43: PetscReal *bottom, *top, *left, *right; /* boundary values */
44: DM dm; /* distributed array data structure */
45: Mat H; /* Hessian */
46: } AppCtx;
49: /* -------- User-defined Routines --------- */
51: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
52: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
53: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
54: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*);
55: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
56: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
57: 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 );if (ierr) return ierr;
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);
90: DMSetFromOptions(user.dm);
91: DMSetUp(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 -- Section 1.5 Writing Application Codes with PETSc 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 -- Section 1.5 Writing Application Codes with PETSc 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,NULL,"-tao_fddefault",&fddefault);
119: PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);
121: /*
122: Create a matrix data structure to store the Hessian and set
123: the Hessian evalution routine.
124: Set the matrix structure to be used for Hessian evalutions
125: */
126: DMCreateMatrix(user.dm,&user.H);
127: MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);
130: if (fdcoloring) {
131: DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
132: MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
133: ISColoringDestroy(&iscoloring);
134: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
135: MatFDColoringSetFromOptions(matfdcoloring);
136: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
137: } else if (fddefault){
138: TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
139: } else {
140: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
141: }
143: /*
144: If my_monitor option is in command line, then use the user-provided
145: monitoring function
146: */
147: PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
148: if (viewmat){
149: TaoSetMonitor(tao,My_Monitor,NULL,NULL);
150: }
152: /* Check for any tao command line options */
153: TaoSetFromOptions(tao);
155: /* SOLVE THE APPLICATION */
156: TaoSolve(tao);
158: TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
160: /* Free TAO data structures */
161: TaoDestroy(&tao);
163: /* Free PETSc data structures */
164: VecDestroy(&x);
165: MatDestroy(&user.H);
166: if (fdcoloring) {
167: MatFDColoringDestroy(&matfdcoloring);
168: }
169: PetscFree(user.bottom);
170: PetscFree(user.top);
171: PetscFree(user.left);
172: PetscFree(user.right);
173: DMDestroy(&user.dm);
174: PetscFinalize();
175: return ierr;
176: }
178: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
179: {
181: PetscReal fcn;
184: FormFunctionGradient(tao,X,&fcn,G,userCtx);
185: return(0);
186: }
188: /* -------------------------------------------------------------------- */
189: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
191: Input Parameters:
192: . tao - the Tao context
193: . XX - input vector
194: . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
196: Output Parameters:
197: . fcn - the newly evaluated function
198: . GG - vector containing the newly evaluated gradient
199: */
200: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){
202: AppCtx *user = (AppCtx *) userCtx;
204: PetscInt i,j;
205: PetscInt mx=user->mx, my=user->my;
206: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
207: PetscReal ft=0.0;
208: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
209: PetscReal rhx=mx+1, rhy=my+1;
210: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
211: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
212: PetscReal **g, **x;
213: Vec localX;
216: /* Get local mesh boundaries */
217: DMGetLocalVector(user->dm,&localX);
218: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
219: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
221: /* Scatter ghost points to local vector */
222: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
223: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
225: /* Get pointers to vector data */
226: DMDAVecGetArray(user->dm,localX,(void**)&x);
227: DMDAVecGetArray(user->dm,G,(void**)&g);
229: /* Compute function and gradient over the locally owned part of the mesh */
230: for (j=ys; j<ys+ym; j++){
231: for (i=xs; i< xs+xm; i++){
233: xc = x[j][i];
234: xlt=xrb=xl=xr=xb=xt=xc;
236: if (i==0){ /* left side */
237: xl= user->left[j-ys+1];
238: xlt = user->left[j-ys+2];
239: } else {
240: xl = x[j][i-1];
241: }
243: if (j==0){ /* bottom side */
244: xb=user->bottom[i-xs+1];
245: xrb =user->bottom[i-xs+2];
246: } else {
247: xb = x[j-1][i];
248: }
250: if (i+1 == gxs+gxm){ /* right side */
251: xr=user->right[j-ys+1];
252: xrb = user->right[j-ys];
253: } else {
254: xr = x[j][i+1];
255: }
257: if (j+1==gys+gym){ /* top side */
258: xt=user->top[i-xs+1];
259: xlt = user->top[i-xs];
260: }else {
261: xt = x[j+1][i];
262: }
264: if (i>gxs && j+1<gys+gym){
265: xlt = x[j+1][i-1];
266: }
267: if (j>gys && i+1<gxs+gxm){
268: xrb = x[j-1][i+1];
269: }
271: d1 = (xc-xl);
272: d2 = (xc-xr);
273: d3 = (xc-xt);
274: d4 = (xc-xb);
275: d5 = (xr-xrb);
276: d6 = (xrb-xb);
277: d7 = (xlt-xl);
278: d8 = (xt-xlt);
280: df1dxc = d1*hydhx;
281: df2dxc = ( d1*hydhx + d4*hxdhy );
282: df3dxc = d3*hxdhy;
283: df4dxc = ( d2*hydhx + d3*hxdhy );
284: df5dxc = d2*hydhx;
285: df6dxc = d4*hxdhy;
287: d1 *= rhx;
288: d2 *= rhx;
289: d3 *= rhy;
290: d4 *= rhy;
291: d5 *= rhy;
292: d6 *= rhx;
293: d7 *= rhy;
294: d8 *= rhx;
296: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
297: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
298: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
299: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
300: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
301: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
303: ft = ft + (f2 + f4);
305: df1dxc /= f1;
306: df2dxc /= f2;
307: df3dxc /= f3;
308: df4dxc /= f4;
309: df5dxc /= f5;
310: df6dxc /= f6;
312: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
314: }
315: }
317: /* Compute triangular areas along the border of the domain. */
318: if (xs==0){ /* left side */
319: for (j=ys; j<ys+ym; j++){
320: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
321: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
322: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
323: }
324: }
325: if (ys==0){ /* bottom side */
326: for (i=xs; i<xs+xm; i++){
327: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
328: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
329: ft = ft+PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
330: }
331: }
333: if (xs+xm==mx){ /* right side */
334: for (j=ys; j< ys+ym; j++){
335: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
336: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
337: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
338: }
339: }
340: if (ys+ym==my){ /* top side */
341: for (i=xs; i<xs+xm; i++){
342: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
343: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
344: ft = ft+PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
345: }
346: }
348: if (ys==0 && xs==0){
349: d1=(user->left[0]-user->left[1])/hy;
350: d2=(user->bottom[0]-user->bottom[1])*rhx;
351: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
352: }
353: if (ys+ym == my && xs+xm == mx){
354: d1=(user->right[ym+1] - user->right[ym])*rhy;
355: d2=(user->top[xm+1] - user->top[xm])*rhx;
356: ft +=PetscSqrtReal( 1.0 + d1*d1 + d2*d2);
357: }
359: ft=ft*area;
360: MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);
362: /* Restore vectors */
363: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
364: DMDAVecRestoreArray(user->dm,G,(void**)&g);
366: /* Scatter values to global vector */
367: DMRestoreLocalVector(user->dm,&localX);
368: PetscLogFlops(67.0*xm*ym);
369: return(0);
370: }
372: /* ------------------------------------------------------------------- */
373: /*
374: FormHessian - Evaluates Hessian matrix.
376: Input Parameters:
377: . tao - the Tao context
378: . x - input vector
379: . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
381: Output Parameters:
382: . H - Hessian matrix
383: . Hpre - optionally different preconditioning matrix
384: . flg - flag indicating matrix structure
386: */
387: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
388: {
390: AppCtx *user = (AppCtx *) ptr;
393: /* Evaluate the Hessian entries*/
394: QuadraticH(user,X,H);
395: return(0);
396: }
398: /* ------------------------------------------------------------------- */
399: /*
400: QuadraticH - Evaluates Hessian matrix.
402: Input Parameters:
403: . user - user-defined context, as set by TaoSetHessian()
404: . X - input vector
406: Output Parameter:
407: . H - Hessian matrix
408: */
409: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
410: {
411: PetscErrorCode ierr;
412: PetscInt i,j,k;
413: PetscInt mx=user->mx, my=user->my;
414: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
415: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
416: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
417: PetscReal hl,hr,ht,hb,hc,htl,hbr;
418: PetscReal **x, v[7];
419: MatStencil col[7],row;
420: Vec localX;
421: PetscBool assembled;
424: /* Get local mesh boundaries */
425: DMGetLocalVector(user->dm,&localX);
427: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
428: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
430: /* Scatter ghost points to local vector */
431: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
432: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
434: /* Get pointers to vector data */
435: DMDAVecGetArray(user->dm,localX,(void**)&x);
437: /* Initialize matrix entries to zero */
438: MatAssembled(Hessian,&assembled);
439: if (assembled){MatZeroEntries(Hessian);}
442: /* Set various matrix options */
443: MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
445: /* Compute Hessian over the locally owned part of the mesh */
447: for (j=ys; j<ys+ym; j++){
449: for (i=xs; i< xs+xm; i++){
451: xc = x[j][i];
452: xlt=xrb=xl=xr=xb=xt=xc;
454: /* Left side */
455: if (i==0){
456: xl = user->left[j-ys+1];
457: xlt = user->left[j-ys+2];
458: } else {
459: xl = x[j][i-1];
460: }
462: if (j==0){
463: xb = user->bottom[i-xs+1];
464: xrb = user->bottom[i-xs+2];
465: } else {
466: xb = x[j-1][i];
467: }
469: if (i+1 == mx){
470: xr = user->right[j-ys+1];
471: xrb = user->right[j-ys];
472: } else {
473: xr = x[j][i+1];
474: }
476: if (j+1==my){
477: xt = user->top[i-xs+1];
478: xlt = user->top[i-xs];
479: }else {
480: xt = x[j+1][i];
481: }
483: if (i>0 && j+1<my){
484: xlt = x[j+1][i-1];
485: }
486: if (j>0 && i+1<mx){
487: xrb = x[j-1][i+1];
488: }
491: d1 = (xc-xl)/hx;
492: d2 = (xc-xr)/hx;
493: d3 = (xc-xt)/hy;
494: d4 = (xc-xb)/hy;
495: d5 = (xrb-xr)/hy;
496: d6 = (xrb-xb)/hx;
497: d7 = (xlt-xl)/hy;
498: d8 = (xlt-xt)/hx;
500: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
501: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
502: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
503: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
504: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
505: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
508: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
509: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
510: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
511: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
512: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
513: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
514: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
515: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
517: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
518: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
520: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
521: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
522: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
523: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
525: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
527: row.j = j; row.i = i;
528: k=0;
529: if (j>0){
530: v[k]=hb;
531: col[k].j = j - 1; col[k].i = i;
532: k++;
533: }
535: if (j>0 && i < mx -1){
536: v[k]=hbr;
537: col[k].j = j - 1; col[k].i = i+1;
538: k++;
539: }
541: if (i>0){
542: v[k]= hl;
543: col[k].j = j; col[k].i = i-1;
544: k++;
545: }
547: v[k]= hc;
548: col[k].j = j; col[k].i = i;
549: k++;
551: if (i < mx-1 ){
552: v[k]= hr;
553: col[k].j = j; col[k].i = i+1;
554: k++;
555: }
557: if (i>0 && j < my-1 ){
558: v[k]= htl;
559: col[k].j = j+1; col[k].i = i-1;
560: k++;
561: }
563: if (j < my-1 ){
564: v[k]= ht;
565: col[k].j = j+1; col[k].i = i;
566: k++;
567: }
569: /*
570: Set matrix values using local numbering, which was defined
571: earlier, in the main routine.
572: */
573: MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
574: }
575: }
577: DMDAVecRestoreArray(user->dm,localX,(void**)&x);
578: DMRestoreLocalVector(user->dm,&localX);
580: MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
581: MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);
583: PetscLogFlops(199.0*xm*ym);
584: return(0);
585: }
587: /* ------------------------------------------------------------------- */
588: /*
589: MSA_BoundaryConditions - Calculates the boundary conditions for
590: the region.
592: Input Parameter:
593: . user - user-defined Section 1.5 Writing Application Codes with PETSc context
595: Output Parameter:
596: . user - user-defined Section 1.5 Writing Application Codes with PETSc context
597: */
598: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
599: {
601: PetscInt i,j,k,limit=0,maxits=5;
602: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
603: PetscInt mx=user->mx,my=user->my;
604: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
605: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
606: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
607: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
608: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
609: PetscReal *boundary;
610: PetscBool flg;
613: /* Get local mesh boundaries */
614: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
615: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
617: bsize=xm+2;
618: lsize=ym+2;
619: rsize=ym+2;
620: tsize=xm+2;
622: PetscMalloc1(bsize,&user->bottom);
623: PetscMalloc1(tsize,&user->top);
624: PetscMalloc1(lsize,&user->left);
625: PetscMalloc1(rsize,&user->right);
627: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
629: for (j=0; j<4; j++){
630: if (j==0){
631: yt=b;
632: xt=l+hx*xs;
633: limit=bsize;
634: boundary=user->bottom;
635: } else if (j==1){
636: yt=t;
637: xt=l+hx*xs;
638: limit=tsize;
639: boundary=user->top;
640: } else if (j==2){
641: yt=b+hy*ys;
642: xt=l;
643: limit=lsize;
644: boundary=user->left;
645: } else { /* if (j==3) */
646: yt=b+hy*ys;
647: xt=r;
648: limit=rsize;
649: boundary=user->right;
650: }
652: for (i=0; i<limit; i++){
653: u1=xt;
654: u2=-yt;
655: for (k=0; k<maxits; k++){
656: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
657: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
658: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
659: if (fnorm <= tol) break;
660: njac11=one+u2*u2-u1*u1;
661: njac12=two*u1*u2;
662: njac21=-two*u1*u2;
663: njac22=-one - u1*u1 + u2*u2;
664: det = njac11*njac22-njac21*njac12;
665: u1 = u1-(njac22*nf1-njac12*nf2)/det;
666: u2 = u2-(njac11*nf2-njac21*nf1)/det;
667: }
669: boundary[i]=u1*u1-u2*u2;
670: if (j==0 || j==1) {
671: xt=xt+hx;
672: } else { /* if (j==2 || j==3) */
673: yt=yt+hy;
674: }
676: }
678: }
680: /* Scale the boundary if desired */
681: if (1==1){
682: PetscReal scl = 1.0;
684: PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
685: if (flg){
686: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
687: }
689: PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
690: if (flg){
691: for (i=0;i<tsize;i++) user->top[i]*=scl;
692: }
694: PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
695: if (flg){
696: for (i=0;i<rsize;i++) user->right[i]*=scl;
697: }
699: PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
700: if (flg){
701: for (i=0;i<lsize;i++) user->left[i]*=scl;
702: }
703: }
704: return(0);
705: }
707: /* ------------------------------------------------------------------- */
708: /*
709: MSA_InitialPoint - Calculates the initial guess in one of three ways.
711: Input Parameters:
712: . user - user-defined Section 1.5 Writing Application Codes with PETSc context
713: . X - vector for initial guess
715: Output Parameters:
716: . X - newly computed initial guess
717: */
718: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
719: {
721: PetscInt start2=-1,i,j;
722: PetscReal start1=0;
723: PetscBool flg1,flg2;
726: PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
727: PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);
729: if (flg1){ /* The zero vector is reasonable */
731: VecSet(X, start1);
733: } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */
734: PetscRandom rctx; PetscReal np5=-0.5;
736: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
737: VecSetRandom(X, rctx);
738: PetscRandomDestroy(&rctx);
739: VecShift(X, np5);
741: } else { /* Take an average of the boundary conditions */
742: PetscInt xs,xm,ys,ym;
743: PetscInt mx=user->mx,my=user->my;
744: PetscReal **x;
746: /* Get local mesh boundaries */
747: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
749: /* Get pointers to vector data */
750: DMDAVecGetArray(user->dm,X,(void**)&x);
752: /* Perform local computations */
753: for (j=ys; j<ys+ym; j++){
754: for (i=xs; i< xs+xm; i++){
755: 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;
756: }
757: }
758: DMDAVecRestoreArray(user->dm,X,(void**)&x);
759: PetscLogFlops(9.0*xm*ym);
760: }
761: return(0);
762: }
764: /*-----------------------------------------------------------------------*/
765: PetscErrorCode My_Monitor(Tao tao, void *ctx)
766: {
768: Vec X;
771: TaoGetSolutionVector(tao,&X);
772: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
773: return(0);
774: }
777: /*TEST
779: build:
780: requires: !complex
782: test:
783: args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
784: requires: !single
786: test:
787: suffix: 2
788: nsize: 2
789: args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
790: filter: grep -v "nls ksp"
791: requires: !single
793: test:
794: suffix: 3
795: nsize: 3
796: args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
797: requires: !single
799: test:
800: suffix: 5
801: nsize: 2
802: args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
803: requires: !single
805: TEST*/