Actual source code: ex10.c
petsc-3.3-p6 2013-02-11
1: #include <petscsnes.h>
2: #include <../src/snes/impls/vi/viimpl.h>
4: static char help[] =
5: "This example is copied from the TAO package\n\
6: (src/complementarity/examples/tutorials/minsurf1.c). It solves a\n\
7: system of nonlinear equations in mixed complementarity form using\n\
8: semismooth newton algorithm.This example is based on a\n\
9: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and\n\
10: boundary values along the edges of the domain, the objective is to find the\n\
11: surface with the minimal area that satisfies the boundary conditions.\n\
12: This application solves this problem using complimentarity -- We are actually\n\
13: solving the system (grad f)_i >= 0, if x_i == l_i \n\
14: (grad f)_i = 0, if l_i < x_i < u_i \n\
15: (grad f)_i <= 0, if x_i == u_i \n\
16: where f is the function to be minimized. \n\
17: \n\
18: The command line options are:\n\
19: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
20: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
21: -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise\n\
22: -lb <value>, lower bound on the variables\n\
23: -ub <value>, upper bound on the variables\n\n";
25: /*
26: User-defined application context - contains data needed by the
27: application-provided call-back routines, FormJacobian() and
28: FormFunction().
29: */
31: typedef struct {
32: PetscInt mx, my;
33: PetscScalar *bottom, *top, *left, *right;
34: } AppCtx;
37: /* -------- User-defined Routines --------- */
39: extern PetscErrorCode MSA_BoundaryConditions(AppCtx *);
40: extern PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
41: extern PetscErrorCode FormGradient(SNES, Vec, Vec, void *);
42: extern PetscErrorCode FormJacobian(SNES, Vec, Mat *, Mat*, MatStructure*,void *);
46: int main(int argc, char **argv)
47: {
48: PetscErrorCode info; /* used to check for functions returning nonzeros */
49: Vec x,r; /* solution and residual vectors */
50: Vec xl,xu; /* Bounds on the variables */
51: PetscBool flg; /* A return variable when checking for user options */
52: SNES snes; /* nonlinear solver context */
53: Mat J; /* Jacobian matrix */
54: PetscInt N; /* Number of elements in vector */
55: PetscScalar lb = SNES_VI_NINF; /* lower bound constant */
56: PetscScalar ub = SNES_VI_INF; /* upper bound constant */
57: AppCtx user; /* user-defined work context */
59: /* Initialize PETSc */
60: PetscInitialize(&argc, &argv, (char *)0, help );
62: #if defined(PETSC_USE_COMPLEX)
63: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This example does not work for scalar type complex\n");
64: #endif
65: /* Specify default dimension of the problem */
66: user.mx = 4; user.my = 4;
68: /* Check for any command line arguments that override defaults */
69: info = PetscOptionsGetInt(PETSC_NULL, "-mx", &user.mx, &flg); CHKERRQ(info);
70: info = PetscOptionsGetInt(PETSC_NULL, "-my", &user.my, &flg); CHKERRQ(info);
71: info = PetscOptionsGetScalar(PETSC_NULL, "-lb", &lb, &flg);CHKERRQ(info);
72: info = PetscOptionsGetScalar(PETSC_NULL, "-ub", &ub, &flg);CHKERRQ(info);
75: /* Calculate any derived values from parameters */
76: N = user.mx*user.my;
77:
78: PetscPrintf(PETSC_COMM_SELF,"\n---- Minimum Surface Area Problem -----\n");
79: PetscPrintf(PETSC_COMM_SELF,"mx:%d, my:%d\n", user.mx,user.my);
81: /* Create appropriate vectors and matrices */
82: info = VecCreateSeq(MPI_COMM_SELF, N, &x);
83: info = VecDuplicate(x, &r); CHKERRQ(info);
84: info = MatCreateSeqAIJ(MPI_COMM_SELF, N, N, 7, PETSC_NULL, &J); CHKERRQ(info);
86: /* Create nonlinear solver context */
87: info = SNESCreate(PETSC_COMM_SELF,&snes); CHKERRQ(info);
90: /* Set function evaluation and Jacobian evaluation routines */
91: info = SNESSetFunction(snes,r,FormGradient,&user);CHKERRQ(info);
92: info = SNESSetJacobian(snes,J,J,FormJacobian,&user);CHKERRQ(info);
94: /* Set the variable bounds */
95: info = MSA_BoundaryConditions(&user); CHKERRQ(info);
97: /* Set initial solution guess */
98: info = MSA_InitialPoint(&user, x); CHKERRQ(info);
100: info = SNESSetFromOptions(snes);CHKERRQ(info);
102: /* Set Bounds on variables */
103: info = VecDuplicate(x, &xl); CHKERRQ(info);
104: info = VecDuplicate(x, &xu); CHKERRQ(info);
105: info = VecSet(xl, lb); CHKERRQ(info);
106: info = VecSet(xu, ub); CHKERRQ(info);
108: info = SNESVISetVariableBounds(snes,xl,xu);CHKERRQ(info);
110: /* Solve the application */
111: info = SNESSolve(snes,PETSC_NULL,x);CHKERRQ(info);
113: info = VecView(x,PETSC_VIEWER_STDOUT_SELF); CHKERRQ(info);
115: /* Free memory */
116: info = VecDestroy(&x); CHKERRQ(info);
117: info = VecDestroy(&xl); CHKERRQ(info);
118: info = VecDestroy(&xu); CHKERRQ(info);
119: info = VecDestroy(&r); CHKERRQ(info);
120: info = MatDestroy(&J); CHKERRQ(info);
121: info = SNESDestroy(&snes); CHKERRQ(info);
123: /* Free user-created data structures */
124: info = PetscFree(user.bottom); CHKERRQ(info);
125: info = PetscFree(user.top); CHKERRQ(info);
126: info = PetscFree(user.left); CHKERRQ(info);
127: info = PetscFree(user.right); CHKERRQ(info);
129: info = PetscFinalize();
131: return 0;
132: }
134: /* -------------------------------------------------------------------- */
138: /* FormGradient - Evaluates gradient of f.
140: Input Parameters:
141: . snes - the SNES context
142: . X - input vector
143: . ptr - optional user-defined context, as set by SNESSetFunction()
144:
145: Output Parameters:
146: . G - vector containing the newly evaluated gradient
147: */
148: int FormGradient(SNES snes, Vec X, Vec G, void *ptr){
149: AppCtx *user = (AppCtx *) ptr;
150: int info;
151: PetscInt i,j,row;
152: PetscInt mx=user->mx, my=user->my;
153: PetscScalar hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
154: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
155: PetscScalar df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
156: PetscScalar zero=0.0;
157: PetscScalar *g, *x;
159: /* Initialize vector to zero */
160: info = VecSet(G, zero); CHKERRQ(info);
162: /* Get pointers to vector data */
163: info = VecGetArray(X, &x); CHKERRQ(info);
164: info = VecGetArray(G, &g); CHKERRQ(info);
166: /* Compute function over the locally owned part of the mesh */
167: for (j=0; j<my; j++){
168: for (i=0; i< mx; i++){
169: row= j*mx + i;
170:
171: xc = x[row];
172: xlt=xrb=xl=xr=xb=xt=xc;
173:
174: if (i==0){ /* left side */
175: xl= user->left[j+1];
176: xlt = user->left[j+2];
177: } else {
178: xl = x[row-1];
179: }
181: if (j==0){ /* bottom side */
182: xb=user->bottom[i+1];
183: xrb = user->bottom[i+2];
184: } else {
185: xb = x[row-mx];
186: }
187:
188: if (i+1 == mx){ /* right side */
189: xr=user->right[j+1];
190: xrb = user->right[j];
191: } else {
192: xr = x[row+1];
193: }
195: if (j+1==0+my){ /* top side */
196: xt=user->top[i+1];
197: xlt = user->top[i];
198: }else {
199: xt = x[row+mx];
200: }
202: if (i>0 && j+1<my){
203: xlt = x[row-1+mx];
204: }
205: if (j>0 && i+1<mx){
206: xrb = x[row+1-mx];
207: }
209: d1 = (xc-xl);
210: d2 = (xc-xr);
211: d3 = (xc-xt);
212: d4 = (xc-xb);
213: d5 = (xr-xrb);
214: d6 = (xrb-xb);
215: d7 = (xlt-xl);
216: d8 = (xt-xlt);
217:
218: df1dxc = d1*hydhx;
219: df2dxc = ( d1*hydhx + d4*hxdhy );
220: df3dxc = d3*hxdhy;
221: df4dxc = ( d2*hydhx + d3*hxdhy );
222: df5dxc = d2*hydhx;
223: df6dxc = d4*hxdhy;
225: d1 /= hx;
226: d2 /= hx;
227: d3 /= hy;
228: d4 /= hy;
229: d5 /= hy;
230: d6 /= hx;
231: d7 /= hy;
232: d8 /= hx;
234: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
235: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
236: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
237: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
238: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
239: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
240:
241: df1dxc /= f1;
242: df2dxc /= f2;
243: df3dxc /= f3;
244: df4dxc /= f4;
245: df5dxc /= f5;
246: df6dxc /= f6;
248: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc )/2.0;
249:
250: }
251: }
252:
253: /* Restore vectors */
254: info = VecRestoreArray(X, &x); CHKERRQ(info);
255: info = VecRestoreArray(G, &g); CHKERRQ(info);
256: info = PetscLogFlops(67*mx*my); CHKERRQ(info);
257: return 0;
258: }
260: /* ------------------------------------------------------------------- */
263: /*
264: FormJacobian - Evaluates Jacobian matrix.
266: Input Parameters:
267: . snes - SNES context
268: . X - input vector
269: . ptr - optional user-defined context, as set by SNESSetJacobian()
271: Output Parameters:
272: . tH - Jacobian matrix
274: */
275: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat *tH, Mat* tHPre, MatStructure* flag, void *ptr)
276: {
277: AppCtx *user = (AppCtx *) ptr;
278: Mat H = *tH;
279: PetscErrorCode info;
280: PetscInt i,j,k,row;
281: PetscInt mx=user->mx, my=user->my;
282: PetscInt col[7];
283: PetscScalar hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
284: PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
285: PetscScalar hl,hr,ht,hb,hc,htl,hbr;
286: PetscScalar *x, v[7];
287: PetscBool assembled;
289: /* Set various matrix options */
290: info = MatSetOption(H,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE); CHKERRQ(info);
291: info = MatAssembled(H,&assembled); CHKERRQ(info);
292: if (assembled){info = MatZeroEntries(H); CHKERRQ(info);}
293: *flag=SAME_NONZERO_PATTERN;
295: /* Get pointers to vector data */
296: info = VecGetArray(X, &x); CHKERRQ(info);
298: /* Compute Jacobian over the locally owned part of the mesh */
299: for (i=0; i< mx; i++){
300: for (j=0; j<my; j++){
301: row= j*mx + i;
302:
303: xc = x[row];
304: xlt=xrb=xl=xr=xb=xt=xc;
306: /* Left side */
307: if (i==0){
308: xl= user->left[j+1];
309: xlt = user->left[j+2];
310: } else {
311: xl = x[row-1];
312: }
313:
314: if (j==0){
315: xb=user->bottom[i+1];
316: xrb = user->bottom[i+2];
317: } else {
318: xb = x[row-mx];
319: }
320:
321: if (i+1 == mx){
322: xr=user->right[j+1];
323: xrb = user->right[j];
324: } else {
325: xr = x[row+1];
326: }
328: if (j+1==my){
329: xt=user->top[i+1];
330: xlt = user->top[i];
331: }else {
332: xt = x[row+mx];
333: }
335: if (i>0 && j+1<my){
336: xlt = x[row-1+mx];
337: }
338: if (j>0 && i+1<mx){
339: xrb = x[row+1-mx];
340: }
343: d1 = (xc-xl)/hx;
344: d2 = (xc-xr)/hx;
345: d3 = (xc-xt)/hy;
346: d4 = (xc-xb)/hy;
347: d5 = (xrb-xr)/hy;
348: d6 = (xrb-xb)/hx;
349: d7 = (xlt-xl)/hy;
350: d8 = (xlt-xt)/hx;
351:
352: f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
353: f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
354: f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
355: f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
356: f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
357: f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);
360: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
361: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
362: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
363: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
364: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
365: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
366: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
367: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
369: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
370: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
372: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
373: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
374: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
375: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
377: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
379: k=0;
380: if (j>0){
381: v[k]=hb; col[k]=row - mx; k++;
382: }
383:
384: if (j>0 && i < mx -1){
385: v[k]=hbr; col[k]=row - mx+1; k++;
386: }
387:
388: if (i>0){
389: v[k]= hl; col[k]=row - 1; k++;
390: }
391:
392: v[k]= hc; col[k]=row; k++;
393:
394: if (i < mx-1 ){
395: v[k]= hr; col[k]=row+1; k++;
396: }
397:
398: if (i>0 && j < my-1 ){
399: v[k]= htl; col[k] = row+mx-1; k++;
400: }
401:
402: if (j < my-1 ){
403: v[k]= ht; col[k] = row+mx; k++;
404: }
405:
406: /*
407: Set matrix values using local numbering, which was defined
408: earlier, in the main routine.
409: */
410: info = MatSetValues(H,1,&row,k,col,v,INSERT_VALUES);
411: CHKERRQ(info);
412: }
413: }
414:
415: /* Restore vectors */
416: info = VecRestoreArray(X,&x); CHKERRQ(info);
418: /* Assemble the matrix */
419: info = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
420: info = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
421: info = PetscLogFlops(199*mx*my); CHKERRQ(info);
422: return 0;
423: }
425: /* ------------------------------------------------------------------- */
428: /*
429: MSA_BoundaryConditions - Calculates the boundary conditions for
430: the region.
432: Input Parameter:
433: . user - user-defined application context
435: Output Parameter:
436: . user - user-defined application context
437: */
438: PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
439: {
440: PetscErrorCode info;
441: PetscInt i,j,k,limit=0,maxits=5;
442: PetscInt mx=user->mx,my=user->my;
443: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
444: PetscScalar one=1.0, two=2.0, three=3.0, tol=1e-10;
445: PetscScalar fnorm,det,hx,hy,xt=0,yt=0;
446: PetscScalar u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
447: PetscScalar b=-0.5, t=0.5, l=-0.5, r=0.5;
448: PetscScalar *boundary;
450: bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2;
452: info = PetscMalloc(bsize*sizeof(PetscScalar), &user->bottom);CHKERRQ(info);
453: info = PetscMalloc(tsize*sizeof(PetscScalar), &user->top);CHKERRQ(info);
454: info = PetscMalloc(lsize*sizeof(PetscScalar), &user->left);CHKERRQ(info);
455: info = PetscMalloc(rsize*sizeof(PetscScalar), &user->right);CHKERRQ(info);
457: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
459: for (j=0; j<4; j++){
460: if (j==0){
461: yt=b;
462: xt=l;
463: limit=bsize;
464: boundary=user->bottom;
465: } else if (j==1){
466: yt=t;
467: xt=l;
468: limit=tsize;
469: boundary=user->top;
470: } else if (j==2){
471: yt=b;
472: xt=l;
473: limit=lsize;
474: boundary=user->left;
475: } else { // if (j==3)
476: yt=b;
477: xt=r;
478: limit=rsize;
479: boundary=user->right;
480: }
482: for (i=0; i<limit; i++){
483: u1=xt;
484: u2=-yt;
485: for (k=0; k<maxits; k++){
486: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
487: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
488: fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
489: if (fnorm <= tol) break;
490: njac11=one+u2*u2-u1*u1;
491: njac12=two*u1*u2;
492: njac21=-two*u1*u2;
493: njac22=-one - u1*u1 + u2*u2;
494: det = njac11*njac22-njac21*njac12;
495: u1 = u1-(njac22*nf1-njac12*nf2)/det;
496: u2 = u2-(njac11*nf2-njac21*nf1)/det;
497: }
499: boundary[i]=u1*u1-u2*u2;
500: if (j==0 || j==1) {
501: xt=xt+hx;
502: } else { // if (j==2 || j==3)
503: yt=yt+hy;
504: }
505: }
506: }
508: return 0;
509: }
511: /* ------------------------------------------------------------------- */
514: /*
515: MSA_InitialPoint - Calculates the initial guess in one of three ways.
517: Input Parameters:
518: . user - user-defined application context
519: . X - vector for initial guess
521: Output Parameters:
522: . X - newly computed initial guess
523: */
524: PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
525: {
526: PetscErrorCode info;
527: PetscInt start=-1,i,j;
528: PetscScalar zero=0.0;
529: PetscBool flg;
531: info = PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); CHKERRQ(info);
533: if (flg && start==0){ /* The zero vector is reasonable */
534:
535: info = VecSet(X, zero); CHKERRQ(info);
536: /* PLogInfo(user,"Min. Surface Area Problem: Start with 0 vector \n"); */
539: } else { /* Take an average of the boundary conditions */
541: PetscInt row;
542: PetscInt mx=user->mx,my=user->my;
543: PetscScalar *x;
544:
545: /* Get pointers to vector data */
546: info = VecGetArray(X,&x); CHKERRQ(info);
548: /* Perform local computations */
549: for (j=0; j<my; j++){
550: for (i=0; i< mx; i++){
551: row=(j)*mx + (i);
552: x[row] = ( ((j+1)*user->bottom[i+1]+(my-j+1)*user->top[i+1])/(my+2)+
553: ((i+1)*user->left[j+1]+(mx-i+1)*user->right[j+1])/(mx+2))/2.0;
554: }
555: }
556:
557: /* Restore vectors */
558: info = VecRestoreArray(X,&x); CHKERRQ(info);
559:
560: }
561: return 0;
562: }