Actual source code: ex15.c
petsc-3.3-p7 2013-05-11
1: #include <petscsnes.h>
2: #include <petscdmda.h>
3: #include <../src/snes/impls/vi/viimpl.h>
4: #include <math.h> /* for cos() sin(0), and atan() */
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 *, MatStructure *, void *);
29: PetscErrorCode ComputeB(AppCtx*);
33: int main( int argc, char **argv )
34: {
35: PetscErrorCode info; /* 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;
45:
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: info = PetscOptionsGetReal(PETSC_NULL,"-ecc",&user.ecc,&flg); CHKERRQ(info);
58: info = PetscOptionsGetReal(PETSC_NULL,"-b",&user.b,&flg); CHKERRQ(info);
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: info = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-50,-50,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&user.da);CHKERRQ(info);
66: 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);
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: info = DMCreateGlobalVector(user.da,&x); CHKERRQ(info); /* Solution */
78: info = VecDuplicate(x,&user.B); CHKERRQ(info); /* Linear objective */
79: info = VecDuplicate(x,&r);CHKERRQ(info);
81: /* Create matrix user.A to store quadratic, Create a local ordering scheme. */
82: info = DMCreateMatrix(user.da,MATAIJ,&user.A);CHKERRQ(info);
83:
84: /* User defined function -- compute linear term of quadratic */
85: info = ComputeB(&user); CHKERRQ(info);
87: /* Create nonlinear solver context */
88: info = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(info);
90: /* Set function evaluation and Jacobian evaluation routines */
91: info = SNESSetFunction(snes,r,FormGradient,&user);CHKERRQ(info);
92: info = SNESSetJacobian(snes,user.A,user.A,FormHessian,&user);CHKERRQ(info);
94: /* Set the initial solution guess */
95: info = VecSet(x, zero); CHKERRQ(info);
97: info = SNESSetFromOptions(snes);CHKERRQ(info);
99: /* Set variable bounds */
100: info = VecDuplicate(x,&xl);CHKERRQ(info);
101: info = VecDuplicate(x,&xu);CHKERRQ(info);
102: info = VecSet(xl,zero);CHKERRQ(info);
103: info = VecSet(xu,thnd);CHKERRQ(info);
104: info = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(info);
106: /* Solve the application */
107: info = SNESSolve(snes,PETSC_NULL,x);CHKERRQ(info);
109: info = SNESGetConvergedReason(snes,&reason); CHKERRQ(info);
110: if (reason <= 0) {
111: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"The SNESVI solver did not converge, adjust some parameters, or check the function evaluation routines\n");
112: }
113:
114: /* Free memory */
115: info = VecDestroy(&x); CHKERRQ(info);
116: info = VecDestroy(&xl);CHKERRQ(info);
117: info = VecDestroy(&xu);CHKERRQ(info);
118: info = VecDestroy(&r);CHKERRQ(info);
119: info = MatDestroy(&user.A); CHKERRQ(info);
120: info = VecDestroy(&user.B); CHKERRQ(info);
121: info = DMDestroy(&user.da); CHKERRQ(info);
122: info = SNESDestroy(&snes);CHKERRQ(info);
124: info = PetscFinalize();
126: return 0;
127: }
129: static PetscReal p(PetscReal xi, PetscReal ecc)
130: {
131: PetscReal t=1.0+ecc*cos(xi);
132: return(t*t*t);
133: }
137: PetscErrorCode ComputeB(AppCtx* user)
138: {
139: PetscErrorCode info;
140: PetscInt i,j;
141: PetscInt nx,ny,xs,xm,ys,ym;
142: PetscReal two=2.0, pi=4.0*atan(1.0);
143: PetscReal hx,hy,ehxhy;
144: PetscReal temp;
145: PetscReal ecc=user->ecc;
146: PetscReal **b;
149: nx=user->nx;
150: ny=user->ny;
151: hx=two*pi/(nx+1.0);
152: hy=two*user->b/(ny+1.0);
153: ehxhy = ecc*hx*hy;
155: /* Get pointer to local vector data */
156: info = DMDAVecGetArray(user->da,user->B, &b); CHKERRQ(info);
158: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
160: /* Compute the linear term in the objective function */
161: for (i=xs; i<xs+xm; i++){
162: temp=sin((i+1)*hx);
163: for (j=ys; j<ys+ym; j++){
164: b[j][i] = - ehxhy*temp;
165: }
166: }
167: /* Restore vectors */
168: info = DMDAVecRestoreArray(user->da,user->B,&b);CHKERRQ(info);
169: info = PetscLogFlops(5*xm*ym+3*xm); CHKERRQ(info);
171: return(0);
172: }
176: PetscErrorCode FormGradient(SNES snes, Vec X, Vec G,void *ctx)
177: {
178: AppCtx* user=(AppCtx*)ctx;
179: PetscErrorCode info;
180: PetscInt i,j,k,kk;
181: PetscInt row[5],col[5];
182: PetscInt nx,ny,xs,xm,ys,ym;
183: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
184: PetscReal hx,hy,hxhy,hxhx,hyhy;
185: PetscReal xi,v[5];
186: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
187: PetscReal vmiddle, vup, vdown, vleft, vright;
188: PetscReal tt;
189: PetscReal **x,**g;
190: PetscReal zero=0.0;
191: Vec localX;
194: nx=user->nx;
195: ny=user->ny;
196: hx=two*pi/(nx+1.0);
197: hy=two*user->b/(ny+1.0);
198: hxhy=hx*hy;
199: hxhx=one/(hx*hx);
200: hyhy=one/(hy*hy);
202: info = VecSet(G, zero); CHKERRQ(info);
203:
204: /* Get local vector */
205: info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
206: /* Get ghoist points */
207: info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
208: info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
209: /* Get pointer to vector data */
210: info = DMDAVecGetArray(user->da,localX,&x);CHKERRQ(info);
211: info = DMDAVecGetArray(user->da,G,&g);CHKERRQ(info);
213: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
215: for (i=xs; i< xs+xm; i++){
216: xi=(i+1)*hx;
217: trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
218: trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
219: trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
220: trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
221: trule5=trule1; /* L(i,j-1) */
222: trule6=trule2; /* U(i,j+1) */
224: vdown=-(trule5+trule2)*hyhy;
225: vleft=-hxhx*(trule2+trule4);
226: vright= -hxhx*(trule1+trule3);
227: vup=-hyhy*(trule1+trule6);
228: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
230: for (j=ys; j<ys+ym; j++){
231:
232: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
233:
234: k=0;
235: if (j > 0){
236: v[k]=vdown; row[k] = i; col[k] = j-1; k++;
237: }
238:
239: if (i > 0){
240: v[k]= vleft; row[k] = i-1; col[k] = j; k++;
241: }
243: v[k]= vmiddle; row[k] = i; col[k] = j; k++;
244:
245: if (i+1 < nx){
246: v[k]= vright; row[k] = i+1; col[k] = j; k++;
247: }
248:
249: if (j+1 < ny){
250: v[k]= vup; row[k] = i; col[k] = j+1; k++;
251: }
252: tt=0;
253: for (kk=0;kk<k;kk++){
254: tt+=v[kk]*x[col[kk]][row[kk]];
255: }
256: g[j][i] = tt;
258: }
260: }
262: /* Restore vectors */
263: info = DMDAVecRestoreArray(user->da,localX, &x); CHKERRQ(info);
264: info = DMDAVecRestoreArray(user->da,G, &g); CHKERRQ(info);
265: info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);
267: info = VecAXPY(G, one, user->B); CHKERRQ(info);
269: info = PetscLogFlops((91 + 10*ym) * xm); CHKERRQ(info);
270: return(0);
272: }
278: /*
279: FormHessian computes the quadratic term in the quadratic objective function
280: Notice that the objective function in this problem is quadratic (therefore a constant
281: hessian). If using a nonquadratic solver, then you might want to reconsider this function
282: */
283: PetscErrorCode FormHessian(SNES snes,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
284: {
285: AppCtx* user=(AppCtx*)ptr;
286: PetscErrorCode info;
287: PetscInt i,j,k;
288: MatStencil row,col[5];
289: PetscInt nx,ny,xs,xm,ys,ym;
290: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
291: PetscReal hx,hy,hxhy,hxhx,hyhy;
292: PetscReal xi,v[5];
293: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
294: PetscReal vmiddle, vup, vdown, vleft, vright;
295: Mat hes=*H;
296: PetscBool assembled;
297: PetscReal **x;
298: Vec localX;
301: nx=user->nx;
302: ny=user->ny;
303: hx=two*pi/(nx+1.0);
304: hy=two*user->b/(ny+1.0);
305: hxhy=hx*hy;
306: hxhx=one/(hx*hx);
307: hyhy=one/(hy*hy);
309: info = MatAssembled(hes,&assembled); CHKERRQ(info);
310: if (assembled){info = MatZeroEntries(hes); CHKERRQ(info);}
311: *flg=SAME_NONZERO_PATTERN;
313: /* Get local vector */
314: info = DMGetLocalVector(user->da,&localX);CHKERRQ(info);
315: /* Get ghost points */
316: info = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
317: info = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(info);
318:
319: /* Get pointers to vector data */
320: info = DMDAVecGetArray(user->da,localX, &x); CHKERRQ(info);
322: info = DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(info);
324: for (i=xs; i< xs+xm; i++){
325: xi=(i+1)*hx;
326: trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
327: trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
328: trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
329: trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
330: trule5=trule1; /* L(i,j-1) */
331: trule6=trule2; /* U(i,j+1) */
333: vdown=-(trule5+trule2)*hyhy;
334: vleft=-hxhx*(trule2+trule4);
335: vright= -hxhx*(trule1+trule3);
336: vup=-hyhy*(trule1+trule6);
337: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
338: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
340: for (j=ys; j<ys+ym; j++){
341: k=0;
342: row.i = i; row.j = j;
343: if (j > 0){
344: v[k]=vdown; col[k].i=i;col[k].j = j-1; k++;
345: }
346:
347: if (i > 0){
348: v[k]= vleft; col[k].i= i-1; col[k].j = j;k++;
349: }
351: v[k]= vmiddle; col[k].i=i; col[k].j = j;k++;
352:
353: if (i+1 < nx){
354: v[k]= vright; col[k].i = i+1; col[k].j = j; k++;
355: }
356:
357: if (j+1 < ny){
358: v[k]= vup; col[k].i = i; col[k].j = j+1; k++;
359: }
360: info = MatSetValuesStencil(hes,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info);
361: }
362: }
364: info = MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
365: info = DMDAVecRestoreArray(user->da,localX,&x);CHKERRQ(info);
366: info = MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
367: info = DMRestoreLocalVector(user->da,&localX);CHKERRQ(info);
369: /*
370: Tell the matrix we will never add a new nonzero location to the
371: matrix. If we do it will generate an error.
372: */
373: info = MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); CHKERRQ(info);
374: info = MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);
376: info = PetscLogFlops(9*xm*ym+49*xm); CHKERRQ(info);
378: return(0);
379: }