Actual source code: ex10.c

  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 --------- */


 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: }