Actual source code: ex15.c
petsc-3.4.5 2014-06-29
1: #include <petscsnes.h>
2: #include <petscdmda.h>
3: #include <../src/snes/impls/vi/viimpl.h>
5: static char help[]=
6: "This example is an implementation of the journal bearing problem from TAO package\n\
7: (src/bound/examples/tutorials/jbearing.c). This example is based on \n\
8: the problem DPJB from the MINPACK-2 test suite. This pressure journal \n\
9: bearing problem is an example of elliptic variational problem defined over \n\
10: a two dimensional rectangle. By discretizing the domain into triangular \n\
11: elements, the pressure surrounding the journal bearing is defined as the \n\
12: minimum of a quadratic function whose variables are bounded below by zero.\n";
14: typedef struct {
15: /* problem parameters */
16: PetscReal ecc; /* test problem parameter */
17: PetscReal b; /* A dimension of journal bearing */
18: PetscInt nx,ny; /* discretization in x, y directions */
19: DM da; /* distributed array data structure */
20: Mat A; /* Quadratic Objective term */
21: Vec B; /* Linear Objective term */
22: } AppCtx;
24: /* User-defined routines */
25: static PetscReal p(PetscReal xi,PetscReal ecc);
26: PetscErrorCode FormGradient(SNES,Vec,Vec,void*);
27: PetscErrorCode FormHessian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
28: PetscErrorCode ComputeB(AppCtx*);
32: int main(int argc, char **argv)
33: {
34: PetscErrorCode info; /* used to check for functions returning nonzeros */
35: Vec x; /* variables vector */
36: Vec xl,xu; /* lower and upper bound on variables */
37: PetscBool flg; /* A return variable when checking for user options */
38: SNESConvergedReason reason;
39: AppCtx user; /* user-defined work context */
40: SNES snes;
41: Vec r;
42: PetscReal zero=0.0,thnd=1000;
45: /* Initialize PETSC */
46: PetscInitialize(&argc, &argv,(char*)0,help);
48: #if defined(PETSC_USE_COMPLEX)
49: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example does not work for scalar type complex\n");
50: #endif
52: /* Set the default values for the problem parameters */
53: user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;
55: /* Check for any command line arguments that override defaults */
56: info = PetscOptionsGetReal(NULL,"-ecc",&user.ecc,&flg);CHKERRQ(info);
57: info = PetscOptionsGetReal(NULL,"-b",&user.b,&flg);CHKERRQ(info);
59: /*
60: A two dimensional distributed array will help define this problem,
61: which derives from an elliptic PDE on two dimensional domain. From
62: the distributed array, Create the vectors.
63: */
64: info = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-50,-50,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&user.da);CHKERRQ(info);
65: info = DMDAGetInfo(user.da,PETSC_IGNORE,&user.nx,&user.ny,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(info);
67: PetscPrintf(PETSC_COMM_WORLD,"\n---- Journal Bearing Problem -----\n");
68: PetscPrintf(PETSC_COMM_WORLD,"mx: %d, my: %d, ecc: %4.3f, b:%3.1f \n",
69: user.nx,user.ny,user.ecc,user.b);
70: /*
71: Extract global and local vectors from DA; the vector user.B is
72: used solely as work space for the evaluation of the function,
73: gradient, and Hessian. Duplicate for remaining vectors that are
74: the same types.
75: */
76: info = DMCreateGlobalVector(user.da,&x);CHKERRQ(info); /* Solution */
77: info = VecDuplicate(x,&user.B);CHKERRQ(info); /* Linear objective */
78: info = VecDuplicate(x,&r);CHKERRQ(info);
80: /* Create matrix user.A to store quadratic, Create a local ordering scheme. */
81: info = DMCreateMatrix(user.da,MATAIJ,&user.A);CHKERRQ(info);
83: /* User defined function -- compute linear term of quadratic */
84: info = ComputeB(&user);CHKERRQ(info);
86: /* Create nonlinear solver context */
87: info = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(info);
89: /* Set function evaluation and Jacobian evaluation routines */
90: info = SNESSetFunction(snes,r,FormGradient,&user);CHKERRQ(info);
91: info = SNESSetJacobian(snes,user.A,user.A,FormHessian,&user);CHKERRQ(info);
93: /* Set the initial solution guess */
94: info = VecSet(x, zero);CHKERRQ(info);
96: info = SNESSetFromOptions(snes);CHKERRQ(info);
98: /* Set variable bounds */
99: info = VecDuplicate(x,&xl);CHKERRQ(info);
100: info = VecDuplicate(x,&xu);CHKERRQ(info);
101: info = VecSet(xl,zero);CHKERRQ(info);
102: info = VecSet(xu,thnd);CHKERRQ(info);
103: info = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(info);
105: /* Solve the application */
106: info = SNESSolve(snes,NULL,x);CHKERRQ(info);
108: info = SNESGetConvergedReason(snes,&reason);CHKERRQ(info);
109: if (reason <= 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"The SNESVI solver did not converge, adjust some parameters, or check the function evaluation routines\n");
111: /* Free memory */
112: info = VecDestroy(&x);CHKERRQ(info);
113: info = VecDestroy(&xl);CHKERRQ(info);
114: info = VecDestroy(&xu);CHKERRQ(info);
115: info = VecDestroy(&r);CHKERRQ(info);
116: info = MatDestroy(&user.A);CHKERRQ(info);
117: info = VecDestroy(&user.B);CHKERRQ(info);
118: info = DMDestroy(&user.da);CHKERRQ(info);
119: info = SNESDestroy(&snes);CHKERRQ(info);
121: info = PetscFinalize();
123: return 0;
124: }
126: static PetscReal p(PetscReal xi, PetscReal ecc)
127: {
128: PetscReal t=1.0+ecc*cos(xi);
129: return(t*t*t);
130: }
134: PetscErrorCode ComputeB(AppCtx *user)
135: {
136: PetscErrorCode info;
137: PetscInt i,j;
138: PetscInt nx,ny,xs,xm,ys,ym;
139: PetscReal two=2.0, pi=4.0*atan(1.0);
140: PetscReal hx,hy,ehxhy;
141: PetscReal temp;
142: PetscReal ecc=user->ecc;
143: PetscReal **b;
146: nx = user->nx;
147: ny = user->ny;
148: hx = two*pi/(nx+1.0);
149: hy = two*user->b/(ny+1.0);
150: ehxhy = ecc*hx*hy;
152: /* Get pointer to local vector data */
153: info = DMDAVecGetArray(user->da,user->B, &b);CHKERRQ(info);
155: info = DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(info);
157: /* Compute the linear term in the objective function */
158: for (i=xs; i<xs+xm; i++) {
159: temp=sin((i+1)*hx);
160: for (j=ys; j<ys+ym; j++) b[j][i] = -ehxhy*temp;
161: }
162: /* Restore vectors */
163: info = DMDAVecRestoreArray(user->da,user->B,&b);CHKERRQ(info);
164: info = PetscLogFlops(5*xm*ym+3*xm);CHKERRQ(info);
165: return(0);
166: }
170: PetscErrorCode FormGradient(SNES snes, Vec X, Vec G,void *ctx)
171: {
172: AppCtx *user=(AppCtx*)ctx;
173: PetscErrorCode info;
174: PetscInt i,j,k,kk;
175: PetscInt row[5],col[5];
176: PetscInt nx,ny,xs,xm,ys,ym;
177: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
178: PetscReal hx,hy,hxhy,hxhx,hyhy;
179: PetscReal xi,v[5];
180: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
181: PetscReal vmiddle, vup, vdown, vleft, vright;
182: PetscReal tt;
183: PetscReal **x,**g;
184: PetscReal zero=0.0;
185: Vec localX;
188: nx = user->nx;
189: ny = user->ny;
190: hx = two*pi/(nx+1.0);
191: hy = two*user->b/(ny+1.0);
192: hxhy = hx*hy;
193: hxhx = one/(hx*hx);
194: hyhy = one/(hy*hy);
196: info = VecSet(G, zero);CHKERRQ(info);
198: /* Get local vector */
199: info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
200: /* Get ghoist points */
201: info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
202: info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
203: /* Get pointer to vector data */
204: info = DMDAVecGetArray(user->da,localX,&x);CHKERRQ(info);
205: info = DMDAVecGetArray(user->da,G,&g);CHKERRQ(info);
207: info = DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(info);
209: for (i=xs; i< xs+xm; i++) {
210: xi = (i+1)*hx;
211: trule1 = hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
212: trule2 = hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
213: trule3 = hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
214: trule4 = hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
215: trule5 = trule1; /* L(i,j-1) */
216: trule6 = trule2; /* U(i,j+1) */
218: vdown = -(trule5+trule2)*hyhy;
219: vleft = -hxhx*(trule2+trule4);
220: vright = -hxhx*(trule1+trule3);
221: vup = -hyhy*(trule1+trule6);
222: vmiddle = (hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
224: for (j=ys; j<ys+ym; j++) {
226: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
228: k=0;
229: if (j > 0) {
230: v[k]=vdown; row[k] = i; col[k] = j-1; k++;
231: }
233: if (i > 0) {
234: v[k]= vleft; row[k] = i-1; col[k] = j; k++;
235: }
237: v[k]= vmiddle; row[k] = i; col[k] = j; k++;
239: if (i+1 < nx) {
240: v[k]= vright; row[k] = i+1; col[k] = j; k++;
241: }
243: if (j+1 < ny) {
244: v[k]= vup; row[k] = i; col[k] = j+1; k++;
245: }
246: tt=0;
247: for (kk=0; kk<k; kk++) tt+=v[kk]*x[col[kk]][row[kk]];
248: g[j][i] = tt;
250: }
252: }
254: /* Restore vectors */
255: info = DMDAVecRestoreArray(user->da,localX, &x);CHKERRQ(info);
256: info = DMDAVecRestoreArray(user->da,G, &g);CHKERRQ(info);
257: info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);
259: info = VecAXPY(G, one, user->B);CHKERRQ(info);
261: info = PetscLogFlops((91 + 10*ym) * xm);CHKERRQ(info);
262: return(0);
263: }
269: /*
270: FormHessian computes the quadratic term in the quadratic objective function
271: Notice that the objective function in this problem is quadratic (therefore a constant
272: hessian). If using a nonquadratic solver, then you might want to reconsider this function
273: */
274: PetscErrorCode FormHessian(SNES snes,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
275: {
276: AppCtx *user=(AppCtx*)ptr;
277: PetscErrorCode info;
278: PetscInt i,j,k;
279: MatStencil row,col[5];
280: PetscInt nx,ny,xs,xm,ys,ym;
281: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
282: PetscReal hx,hy,hxhy,hxhx,hyhy;
283: PetscReal xi,v[5];
284: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
285: PetscReal vmiddle, vup, vdown, vleft, vright;
286: Mat hes=*H;
287: PetscBool assembled;
288: PetscReal **x;
289: Vec localX;
292: nx = user->nx;
293: ny = user->ny;
294: hx = two*pi/(nx+1.0);
295: hy = two*user->b/(ny+1.0);
296: hxhy = hx*hy;
297: hxhx = one/(hx*hx);
298: hyhy = one/(hy*hy);
300: info = MatAssembled(hes,&assembled);CHKERRQ(info);
301: if (assembled) {info = MatZeroEntries(hes);CHKERRQ(info);}
302: *flg=SAME_NONZERO_PATTERN;
304: /* Get local vector */
305: info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
306: /* Get ghost points */
307: info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
308: info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
310: /* Get pointers to vector data */
311: info = DMDAVecGetArray(user->da,localX, &x);CHKERRQ(info);
313: info = DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(info);
315: for (i=xs; i< xs+xm; i++) {
316: xi = (i+1)*hx;
317: trule1 = hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
318: trule2 = hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
319: trule3 = hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
320: trule4 = hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
321: trule5 = trule1; /* L(i,j-1) */
322: trule6 = trule2; /* U(i,j+1) */
324: vdown = -(trule5+trule2)*hyhy;
325: vleft = -hxhx*(trule2+trule4);
326: vright = -hxhx*(trule1+trule3);
327: vup = -hyhy*(trule1+trule6);
328: vmiddle = (hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
330: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
332: for (j=ys; j<ys+ym; j++) {
333: k =0;
334: row.i = i; row.j = j;
335: if (j > 0) {
336: v[k]=vdown; col[k].i=i;col[k].j = j-1; k++;
337: }
339: if (i > 0) {
340: v[k]= vleft; col[k].i= i-1; col[k].j = j;k++;
341: }
343: v[k]= vmiddle; col[k].i=i; col[k].j = j;k++;
345: if (i+1 < nx) {
346: v[k]= vright; col[k].i = i+1; col[k].j = j; k++;
347: }
349: if (j+1 < ny) {
350: v[k]= vup; col[k].i = i; col[k].j = j+1; k++;
351: }
352: info = MatSetValuesStencil(hes,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(info);
353: }
354: }
356: info = MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
357: info = DMDAVecRestoreArray(user->da,localX,&x);CHKERRQ(info);
358: info = MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
359: info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);
361: /*
362: Tell the matrix we will never add a new nonzero location to the
363: matrix. If we do it will generate an error.
364: */
365: info = MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(info);
366: info = MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(info);
368: info = PetscLogFlops(9*xm*ym+49*xm);CHKERRQ(info);
369: return(0);
370: }