Actual source code: ex15.c
petsc-3.7.3 2016-08-01
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,NULL,"-ecc",&user.ecc,&flg);
58: PetscOptionsGetReal(NULL,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);
156: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
158: /* Compute the linear term in the objective function */
159: for (i=xs; i<xs+xm; i++) {
160: temp=PetscSinReal((i+1)*hx);
161: for (j=ys; j<ys+ym; j++) b[j][i] = -ehxhy*temp;
162: }
163: /* Restore vectors */
164: DMDAVecRestoreArray(user->da,user->B,&b);
165: PetscLogFlops(5*xm*ym+3*xm);
166: return(0);
167: }
171: PetscErrorCode FormGradient(SNES snes, Vec X, Vec G,void *ctx)
172: {
173: AppCtx *user=(AppCtx*)ctx;
175: PetscInt i,j,k,kk;
176: PetscInt row[5],col[5];
177: PetscInt nx,ny,xs,xm,ys,ym;
178: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
179: PetscReal hx,hy,hxhy,hxhx,hyhy;
180: PetscReal xi,v[5];
181: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
182: PetscReal vmiddle, vup, vdown, vleft, vright;
183: PetscReal tt;
184: PetscReal **x,**g;
185: PetscReal zero=0.0;
186: Vec localX;
189: nx = user->nx;
190: ny = user->ny;
191: hx = two*pi/(nx+1.0);
192: hy = two*user->b/(ny+1.0);
193: hxhy = hx*hy;
194: hxhx = one/(hx*hx);
195: hyhy = one/(hy*hy);
197: VecSet(G, zero);
199: /* Get local vector */
200: DMGetLocalVector(user->da,&localX);
201: /* Get ghoist points */
202: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
203: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
204: /* Get pointer to vector data */
205: DMDAVecGetArrayRead(user->da,localX,&x);
206: DMDAVecGetArray(user->da,G,&g);
208: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
210: for (i=xs; i< xs+xm; i++) {
211: xi = (i+1)*hx;
212: trule1 = hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
213: trule2 = hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
214: trule3 = hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
215: trule4 = hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
216: trule5 = trule1; /* L(i,j-1) */
217: trule6 = trule2; /* U(i,j+1) */
219: vdown = -(trule5+trule2)*hyhy;
220: vleft = -hxhx*(trule2+trule4);
221: vright = -hxhx*(trule1+trule3);
222: vup = -hyhy*(trule1+trule6);
223: vmiddle = (hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
225: for (j=ys; j<ys+ym; j++) {
227: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
229: k=0;
230: if (j > 0) {
231: v[k]=vdown; row[k] = i; col[k] = j-1; k++;
232: }
234: if (i > 0) {
235: v[k]= vleft; row[k] = i-1; col[k] = j; k++;
236: }
238: v[k]= vmiddle; row[k] = i; col[k] = j; k++;
240: if (i+1 < nx) {
241: v[k]= vright; row[k] = i+1; col[k] = j; k++;
242: }
244: if (j+1 < ny) {
245: v[k]= vup; row[k] = i; col[k] = j+1; k++;
246: }
247: tt=0;
248: for (kk=0; kk<k; kk++) tt+=v[kk]*x[col[kk]][row[kk]];
249: g[j][i] = tt;
251: }
253: }
255: /* Restore vectors */
256: DMDAVecRestoreArrayRead(user->da,localX, &x);
257: DMDAVecRestoreArray(user->da,G, &g);
258: DMRestoreLocalVector(user->da,&localX);
260: VecAXPY(G, one, user->B);
262: PetscLogFlops((91 + 10*ym) * xm);
263: return(0);
264: }
270: /*
271: FormHessian computes the quadratic term in the quadratic objective function
272: Notice that the objective function in this problem is quadratic (therefore a constant
273: hessian). If using a nonquadratic solver, then you might want to reconsider this function
274: */
275: PetscErrorCode FormHessian(SNES snes,Vec X,Mat H, Mat Hpre, void *ptr)
276: {
277: AppCtx *user=(AppCtx*)ptr;
279: PetscInt i,j,k;
280: MatStencil row,col[5];
281: PetscInt nx,ny,xs,xm,ys,ym;
282: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
283: PetscReal hx,hy,hxhy,hxhx,hyhy;
284: PetscReal xi,v[5];
285: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
286: PetscReal vmiddle, vup, vdown, vleft, vright;
287: Mat hes=H;
288: PetscBool assembled;
289: PetscReal **x;
290: Vec localX;
293: nx = user->nx;
294: ny = user->ny;
295: hx = two*pi/(nx+1.0);
296: hy = two*user->b/(ny+1.0);
297: hxhy = hx*hy;
298: hxhx = one/(hx*hx);
299: hyhy = one/(hy*hy);
301: MatAssembled(hes,&assembled);
302: if (assembled) {MatZeroEntries(hes);}
304: /* Get local vector */
305: DMGetLocalVector(user->da,&localX);
306: /* Get ghost points */
307: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
308: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
310: /* Get pointers to vector data */
311: DMDAVecGetArrayRead(user->da,localX, &x);
313: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
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: MatSetValuesStencil(hes,1,&row,k,col,v,INSERT_VALUES);
353: }
354: }
356: MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY);
357: DMDAVecRestoreArrayRead(user->da,localX,&x);
358: MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY);
359: DMRestoreLocalVector(user->da,&localX);
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: MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
366: MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE);
368: PetscLogFlops(9*xm*ym+49*xm);
369: return(0);
370: }