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