Actual source code: ex11.c

  1: /*$Id: ex11.c,v 1.38 1999/11/05 14:47:16 bsmith Exp bsmith $*/

  3: static char help[] =
  4: "This program demonstrates use of the SNES package to solve systems of\n\
  5: nonlinear equations in parallel, using 2-dimensional distributed arrays.\n\
  6: The 2-dim Bratu (SFI - solid fuel ignition) test problem is used, where\n\
  7: analytic formation of the Jacobian is the default.  \n\
  8: \n\
  9:   Solves the linear systems via 2 level additive Schwarz \n\
 10: \n\
 11: The command line\n\
 12: options are:\n\
 13:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
 14:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
 15:   -Mx <xg>, where <xg> = number of grid points in the x-direction on coarse grid\n\
 16:   -My <yg>, where <yg> = number of grid points in the y-direction on coarse grid\n\n";

 18: /*  
 19:     1) Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 20:     the partial differential equation
 21:   
 22:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1 ,
 23:   
 24:     with boundary conditions
 25:    
 26:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 27:   
 28:     A finite difference approximation with the usual 5-point stencil
 29:     is used to discretize the boundary value problem to obtain a nonlinear 
 30:     system of equations.

 32:    The code has two cases for multilevel solver
 33:      I. the coarse grid Jacobian is computed in parallel 
 34:      II. the coarse grid Jacobian is computed sequentially on each processor
 35:    in both cases the coarse problem is SOLVED redundantly.

 37: */

 39:  #include petscsnes.h
 40:  #include petscda.h
 41:  #include petscmg.h

 43: /* User-defined application contexts */

 45: typedef struct {
 46:    int        mx,my;            /* number grid points in x and y direction */
 47:    Vec        localX,localF;    /* local vectors with ghost region */
 48:    DA         da;
 49:    Vec        x,b,r;            /* global vectors */
 50:    Mat        J;                /* Jacobian on grid */
 51: } GridCtx;

 53: typedef struct {
 54:    double      param;           /* test problem parameter */
 55:    GridCtx     fine;
 56:    GridCtx     coarse;
 57:    SLES        sles_coarse;
 58:    SLES        sles_fine;
 59:    int         ratio;
 60:    Mat         R;               /* restriction fine to coarse */
 61:    Vec         Rscale;
 62:    PetscTruth  redundant_build; /* build coarse matrix redundantly */
 63:    Vec         localall;        /* contains entire coarse vector on each processor in NATURAL order*/
 64:    VecScatter  tolocalall;      /* maps from parallel "global" coarse vector to localall */
 65:    VecScatter  fromlocalall;    /* maps from localall vector back to global coarse vector */
 66: } AppCtx;

 68: #define COARSE_LEVEL 0
 69: #define FINE_LEVEL   1

 71: extern int FormFunction(SNES,Vec,Vec,void*), FormInitialGuess1(AppCtx*,Vec);
 72: extern int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 73: extern int FormInterpolation(AppCtx *);

 75: /*
 76:       Mm_ratio - ration of grid lines between fine and coarse grids.
 77: */
 80: int main( int argc, char **argv )
 81: {
 82:   SNES          snes;
 83:   AppCtx        user;
 84:   int           ierr, its, N, n, Nx = PETSC_DECIDE, Ny = PETSC_DECIDE;
 85:   int           size, nlocal,Nlocal;
 86:   double        bratu_lambda_max = 6.81, bratu_lambda_min = 0.;
 87:   SLES          sles;
 88:   PC            pc;

 90:   /*
 91:       Initialize PETSc, note that default options in ex11options can be 
 92:       overridden at the command line
 93:   */
 94:   PetscInitialize( &argc, &argv,"ex11options",help );

 96:   user.ratio = 2;
 97:   user.coarse.mx = 5; user.coarse.my = 5; user.param = 6.0;
 98:   PetscOptionsGetInt(PETSC_NULL,"-Mx",&user.coarse.mx,PETSC_NULL);
 99:   PetscOptionsGetInt(PETSC_NULL,"-My",&user.coarse.my,PETSC_NULL);
100:   PetscOptionsGetInt(PETSC_NULL,"-ratio",&user.ratio,PETSC_NULL);
101:   user.fine.mx = user.ratio*(user.coarse.mx-1)+1; user.fine.my = user.ratio*(user.coarse.my-1)+1;

103:   PetscOptionsHasName(PETSC_NULL,"-redundant_build",&user.redundant_build);
104:   if (user.redundant_build) {
105:     PetscPrintf(PETSC_COMM_WORLD,"Building coarse Jacobian redundantly\n");
106:   }

108:   PetscPrintf(PETSC_COMM_WORLD,"Coarse grid size %d by %d\n",user.coarse.mx,user.coarse.my);
109:   PetscPrintf(PETSC_COMM_WORLD,"Fine grid size %d by %d\n",user.fine.mx,user.fine.my);

111:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
112:   if (user.param >= bratu_lambda_max || user.param < bratu_lambda_min) {
113:     SETERRQ(1,"Lambda is out of range");
114:   }
115:   n = user.fine.mx*user.fine.my; N = user.coarse.mx*user.coarse.my;

117:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
118:   PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
119:   PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);

121:   /* Set up distributed array for fine grid */
122:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.fine.mx,
123:                     user.fine.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.fine.da);
124:   DACreateGlobalVector(user.fine.da,&user.fine.x);
125:   VecDuplicate(user.fine.x,&user.fine.r);
126:   VecDuplicate(user.fine.x,&user.fine.b);
127:   VecGetLocalSize(user.fine.x,&nlocal);
128:   DACreateLocalVector(user.fine.da,&user.fine.localX);
129:   VecDuplicate(user.fine.localX,&user.fine.localF);
130:   MatCreateMPIAIJ(PETSC_COMM_WORLD,nlocal,nlocal,n,n,5,PETSC_NULL,3,PETSC_NULL,&user.fine.J);

132:   /* Set up distributed array for coarse grid */
133:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.coarse.mx,
134:                     user.coarse.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.coarse.da);
135:   DACreateGlobalVector(user.coarse.da,&user.coarse.x);
136:   VecDuplicate(user.coarse.x,&user.coarse.b);
137:   if (user.redundant_build) {
138:     /* Create scatter from parallel global numbering to redundant with natural ordering */
139:     DAGlobalToNaturalAllCreate(user.coarse.da,&user.tolocalall);
140:     DANaturalAllToGlobalCreate(user.coarse.da,&user.fromlocalall);
141:     VecCreateSeq(PETSC_COMM_SELF,N,&user.localall);
142:     /* Create sequential matrix to hold entire coarse grid Jacobian on each processor */
143:     MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,5,PETSC_NULL,&user.coarse.J);
144:   } else {
145:     VecGetLocalSize(user.coarse.x,&Nlocal);
146:     DACreateLocalVector(user.coarse.da,&user.coarse.localX);
147:     VecDuplicate(user.coarse.localX,&user.coarse.localF);
148:     /* We will compute the coarse Jacobian in parallel */
149:     MatCreateMPIAIJ(PETSC_COMM_WORLD,Nlocal,Nlocal,N,N,5,PETSC_NULL,3,PETSC_NULL,&user.coarse.J);
150:   }

152:   /* Create nonlinear solver */
153:   SNESCreate(PETSC_COMM_WORLD,&snes);

155:   /* provide user function and Jacobian */
156:   SNESSetFunction(snes,user.fine.b,FormFunction,&user);
157:   SNESSetJacobian(snes,user.fine.J,user.fine.J,FormJacobian,&user);

159:   /* set two level additive Schwarz preconditioner */
160:   SNESGetSLES(snes,&sles);
161:   SLESGetPC(sles,&pc);
162:   PCSetType(pc,PCMG);
163:   MGSetLevels(pc,2,PETSC_NULL);
164:   MGSetType(pc,MGADDITIVE);

166:   /* always solve the coarse problem redundantly with direct LU solver */
167:   PetscOptionsSetValue("-coarse_pc_type","redundant");
168:   PetscOptionsSetValue("-coarse_redundant_pc_type","lu");

170:   /* Create coarse level */
171:   MGGetCoarseSolve(pc,&user.sles_coarse);
172:   SLESSetOptionsPrefix(user.sles_coarse,"coarse_");
173:   SLESSetFromOptions(user.sles_coarse);
174:   SLESSetOperators(user.sles_coarse,user.coarse.J,user.coarse.J,DIFFERENT_NONZERO_PATTERN);
175:   MGSetX(pc,COARSE_LEVEL,user.coarse.x);
176:   MGSetRhs(pc,COARSE_LEVEL,user.coarse.b);
177:   if (user.redundant_build) {
178:     PC rpc;
179:     SLESGetPC(user.sles_coarse,&rpc);
180:     PCRedundantSetScatter(rpc,user.tolocalall,user.fromlocalall);
181:   }

183:   /* Create fine level */
184:   MGGetSmoother(pc,FINE_LEVEL,&user.sles_fine);
185:   SLESSetOptionsPrefix(user.sles_fine,"fine_");
186:   SLESSetFromOptions(user.sles_fine);
187:   SLESSetOperators(user.sles_fine,user.fine.J,user.fine.J,DIFFERENT_NONZERO_PATTERN);
188:   MGSetR(pc,FINE_LEVEL,user.fine.r);
189:   MGSetResidual(pc,FINE_LEVEL,MGDefaultResidual,user.fine.J);

191:   /* Create interpolation between the levels */
192:   FormInterpolation(&user);
193:   MGSetInterpolate(pc,FINE_LEVEL,user.R);
194:   MGSetRestriction(pc,FINE_LEVEL,user.R);

196:   /* Set options, then solve nonlinear system */
197:   SNESSetFromOptions(snes);
198:   FormInitialGuess1(&user,user.fine.x);
199:   SNESSolve(snes,user.fine.x,&its);
200:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n", its );

202:   /* Free data structures */
203:   if (user.redundant_build) {
204:     VecScatterDestroy(user.tolocalall);
205:     VecScatterDestroy(user.fromlocalall);
206:     VecDestroy(user.localall);
207:   } else {
208:     VecDestroy(user.coarse.localX);
209:     VecDestroy(user.coarse.localF);
210:   }

212:   MatDestroy(user.fine.J);
213:   VecDestroy(user.fine.x);
214:   VecDestroy(user.fine.r);
215:   VecDestroy(user.fine.b);
216:   DADestroy(user.fine.da);
217:   VecDestroy(user.fine.localX);
218:   VecDestroy(user.fine.localF);

220:   MatDestroy(user.coarse.J);
221:   VecDestroy(user.coarse.x);
222:   VecDestroy(user.coarse.b);
223:   DADestroy(user.coarse.da);

225:   SNESDestroy(snes);
226:   MatDestroy(user.R);
227:   VecDestroy(user.Rscale);
228:   PetscFinalize();

230:   return 0;
231: }/* --------------------  Form initial approximation ----------------- */
234: int FormInitialGuess1(AppCtx *user,Vec X)
235: {
236:   int     i, j, row, mx, my, ierr, xs, ys, xm, ym, Xm, Ym, Xs, Ys;
237:   double  one = 1.0, lambda, temp1, temp, hx, hy, hxdhy, hydhx,sc;
238:   PetscScalar  *x;
239:   Vec     localX = user->fine.localX;

241:   mx = user->fine.mx;       my = user->fine.my;            lambda = user->param;
242:   hx = one/(double)(mx-1);  hy = one/(double)(my-1);
243:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

245:   temp1 = lambda/(lambda + one);

247:   /* Get ghost points */
248:   DAGetCorners(user->fine.da,&xs,&ys,0,&xm,&ym,0);
249:   DAGetGhostCorners(user->fine.da,&Xs,&Ys,0,&Xm,&Ym,0);
250:   VecGetArray(localX,&x);

252:   /* Compute initial guess */
253:   for (j=ys; j<ys+ym; j++) {
254:     temp = (double)(PetscMin(j,my-j-1))*hy;
255:     for (i=xs; i<xs+xm; i++) {
256:       row = i - Xs + (j - Ys)*Xm;
257:       if (i == 0 || j == 0 || i == mx-1 || j == my-1 ) {
258:         x[row] = 0.0;
259:         continue;
260:       }
261:       x[row] = temp1*sqrt( PetscMin( (double)(PetscMin(i,mx-i-1))*hx,temp) );
262:     }
263:   }
264:   VecRestoreArray(localX,&x);

266:   /* Insert values into global vector */
267:   DALocalToGlobal(user->fine.da,localX,INSERT_VALUES,X);
268:   return 0;
269: }

271:  /* --------------------  Evaluate Function F(x) --------------------- */
274: int FormFunction(SNES snes,Vec X,Vec F,void *ptr)
275: {
276:   AppCtx  *user = (AppCtx *) ptr;
277:   int     ierr, i, j, row, mx, my, xs, ys, xm, ym, Xs, Ys, Xm, Ym;
278:   double  two = 2.0, one = 1.0, lambda,hx, hy, hxdhy, hydhx,sc;
279:   PetscScalar  u, uxx, uyy, *x,*f;
280:   Vec     localX = user->fine.localX, localF = user->fine.localF;

282:   mx = user->fine.mx;       my = user->fine.my;       lambda = user->param;
283:   hx = one/(double)(mx-1);  hy = one/(double)(my-1);
284:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

286:   /* Get ghost points */
287:   DAGlobalToLocalBegin(user->fine.da,X,INSERT_VALUES,localX);
288:   DAGlobalToLocalEnd(user->fine.da,X,INSERT_VALUES,localX);
289:   DAGetCorners(user->fine.da,&xs,&ys,0,&xm,&ym,0);
290:   DAGetGhostCorners(user->fine.da,&Xs,&Ys,0,&Xm,&Ym,0);
291:   VecGetArray(localX,&x);
292:   VecGetArray(localF,&f);

294:   /* Evaluate function */
295:   for (j=ys; j<ys+ym; j++) {
296:     row = (j - Ys)*Xm + xs - Xs - 1;
297:     for (i=xs; i<xs+xm; i++) {
298:       row++;
299:       if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
300:         u = x[row];
301:         uxx = (two*u - x[row-1] - x[row+1])*hydhx;
302:         uyy = (two*u - x[row-Xm] - x[row+Xm])*hxdhy;
303:         f[row] = uxx + uyy - sc*exp(u);
304:       } else if ((i > 0 && i < mx-1) || (j > 0 && j < my-1)){
305:         f[row] = .5*two*(hydhx + hxdhy)*x[row];
306:       } else {
307:         f[row] = .25*two*(hydhx + hxdhy)*x[row];
308:       }
309:     }
310:   }
311:   VecRestoreArray(localX,&x);
312:   VecRestoreArray(localF,&f);

314:   /* Insert values into global vector */
315:   DALocalToGlobal(user->fine.da,localF,INSERT_VALUES,F);
316:   PetscLogFlops(11*ym*xm);
317:   return 0;
318: }

320: /*
321:         Computes the part of the Jacobian associated with this processor 
322: */
325: int FormJacobian_Grid(AppCtx *user,GridCtx *grid,Vec X, Mat *J,Mat *B)
326: {
327:   Mat     jac = *J;
328:   int     ierr, i, j, row, mx, my, xs, ys, xm, ym, Xs, Ys, Xm, Ym, col[5];
329:   int     nloc, *ltog, grow;
330:   PetscScalar  two = 2.0, one = 1.0, lambda, v[5], hx, hy, hxdhy, hydhx, sc, *x, value;
331:   Vec     localX = grid->localX;

333:   mx = grid->mx;            my = grid->my;            lambda = user->param;
334:   hx = one/(double)(mx-1);  hy = one/(double)(my-1);
335:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

337:   /* Get ghost points */
338:   DAGlobalToLocalBegin(grid->da,X,INSERT_VALUES,localX);
339:   DAGlobalToLocalEnd(grid->da,X,INSERT_VALUES,localX);
340:   DAGetCorners(grid->da,&xs,&ys,0,&xm,&ym,0);
341:   DAGetGhostCorners(grid->da,&Xs,&Ys,0,&Xm,&Ym,0);
342:   DAGetGlobalIndices(grid->da,&nloc,&ltog);
343:   VecGetArray(localX,&x);

345:   /* Evaluate Jacobian of function */
346:   for (j=ys; j<ys+ym; j++) {
347:     row = (j - Ys)*Xm + xs - Xs - 1;
348:     for (i=xs; i<xs+xm; i++) {
349:       row++;
350:       grow = ltog[row];
351:       if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
352:         v[0] = -hxdhy; col[0] = ltog[row - Xm];
353:         v[1] = -hydhx; col[1] = ltog[row - 1];
354:         v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = grow;
355:         v[3] = -hydhx; col[3] = ltog[row + 1];
356:         v[4] = -hxdhy; col[4] = ltog[row + Xm];
357:         MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
358:       } else if ((i > 0 && i < mx-1) || (j > 0 && j < my-1)){
359:         value = .5*two*(hydhx + hxdhy);
360:         MatSetValues(jac,1,&grow,1,&grow,&value,INSERT_VALUES);
361:       } else {
362:         value = .25*two*(hydhx + hxdhy);
363:         MatSetValues(jac,1,&grow,1,&grow,&value,INSERT_VALUES);
364:       }
365:     }
366:   }
367:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
368:   VecRestoreArray(localX,&x);
369:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

371:   return 0;
372: }

374: /*
375:         Computes the ENTIRE Jacobian associated with the ENTIRE grid sequentially
376:     This is for generating the coarse grid redundantly.

378:           This is BAD code duplication, since the bulk of this routine is the
379:     same as the routine above

381:        Note the numbering of the rows/columns is the NATURAL numbering
382: */
385: int FormJacobian_Coarse(AppCtx *user,GridCtx *grid,Vec X, Mat *J,Mat *B)
386: {
387:   Mat     jac = *J;
388:   int     ierr, i, j, row, mx, my, col[5];
389:   PetscScalar  two = 2.0, one = 1.0, lambda, v[5], hx, hy, hxdhy, hydhx, sc, *x, value;

391:   mx = grid->mx;            my = grid->my;            lambda = user->param;
392:   hx = one/(double)(mx-1);  hy = one/(double)(my-1);
393:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

395:   VecGetArray(X,&x);

397:   /* Evaluate Jacobian of function */
398:   for (j=0; j<my; j++) {
399:     row = j*mx - 1;
400:     for (i=0; i<mx; i++) {
401:       row++;
402:       if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
403:         v[0] = -hxdhy; col[0] = row - mx;
404:         v[1] = -hydhx; col[1] = row - 1;
405:         v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = row;
406:         v[3] = -hydhx; col[3] = row + 1;
407:         v[4] = -hxdhy; col[4] = row + mx;
408:         MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
409:       } else if ((i > 0 && i < mx-1) || (j > 0 && j < my-1)){
410:         value = .5*two*(hydhx + hxdhy);
411:         MatSetValues(jac,1,&row,1,&row,&value,INSERT_VALUES);
412:       } else {
413:         value = .25*two*(hydhx + hxdhy);
414:         MatSetValues(jac,1,&row,1,&row,&value,INSERT_VALUES);
415:       }
416:     }
417:   }
418:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
419:   VecRestoreArray(X,&x);
420:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

422:   return 0;
423: }

425: /* --------------------  Evaluate Jacobian F'(x) --------------------- */
428: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
429: {
430:   AppCtx     *user = (AppCtx *) ptr;
431:   int        ierr;
432:   SLES       sles;
433:   PC         pc;
434:   PetscTruth ismg;

436:   *flag = SAME_NONZERO_PATTERN;
437:   FormJacobian_Grid(user,&user->fine,X,J,B);

439:   /* create coarse grid jacobian for preconditioner */
440:   SNESGetSLES(snes,&sles);
441:   SLESGetPC(sles,&pc);
442: 
443:   PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
444:   if (ismg) {

446:     SLESSetOperators(user->sles_fine,user->fine.J,user->fine.J,SAME_NONZERO_PATTERN);

448:     /* restrict X to coarse grid */
449:     MatMult(user->R,X,user->coarse.x);
450:     VecPointwiseMult(user->Rscale,user->coarse.x,user->coarse.x);

452:     /* form Jacobian on coarse grid */
453:     if (user->redundant_build) {
454:       /* get copy of coarse X onto each processor */
455:       VecScatterBegin(user->coarse.x,user->localall,INSERT_VALUES,SCATTER_FORWARD,user->tolocalall);
456:       VecScatterEnd(user->coarse.x,user->localall,INSERT_VALUES,SCATTER_FORWARD,user->tolocalall);
457:       FormJacobian_Coarse(user,&user->coarse,user->localall,&user->coarse.J,&user->coarse.J);

459:     } else {
460:       /* coarse grid Jacobian computed in parallel */
461:       FormJacobian_Grid(user,&user->coarse,user->coarse.x,&user->coarse.J,&user->coarse.J);
462:     }
463:     SLESSetOperators(user->sles_coarse,user->coarse.J,user->coarse.J,SAME_NONZERO_PATTERN);
464:   }

466:   return 0;
467: }


472: /*
473:       Forms the interpolation (and restriction) operator from 
474: coarse grid to fine.
475: */
476: int FormInterpolation(AppCtx *user)
477: {
478:   int      ierr,i,j,i_start,m_fine,j_start,m,n,M,Mx = user->coarse.mx,My = user->coarse.my,*idx;
479:   int      m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,m_coarse;
480:   int      row,i_start_ghost,j_start_ghost,cols[4],mx = user->fine.mx, m_c,my = user->fine.my;
481:   int      c0,c1,c2,c3,nc,ratio = user->ratio,i_end,i_end_ghost,m_c_local,m_fine_local;
482:   int      i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col;
483:   PetscScalar   v[4],x,y, one = 1.0;
484:   Mat      mat;
485:   Vec      Rscale;
486: 
487:   DAGetCorners(user->fine.da,&i_start,&j_start,0,&m,&n,0);
488:   DAGetGhostCorners(user->fine.da,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
489:   DAGetGlobalIndices(user->fine.da,PETSC_NULL,&idx);

491:   DAGetCorners(user->coarse.da,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
492:   DAGetGhostCorners(user->coarse.da,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
493:   DAGetGlobalIndices(user->coarse.da,PETSC_NULL,&idx_c);

495:   /* create interpolation matrix */
496:   VecGetLocalSize(user->fine.x,&m_fine_local);
497:   VecGetLocalSize(user->coarse.x,&m_c_local);
498:   VecGetSize(user->fine.x,&m_fine);
499:   VecGetSize(user->coarse.x,&m_coarse);
500:   MatCreateMPIAIJ(PETSC_COMM_WORLD,m_fine_local,m_c_local,m_fine,m_coarse,
501:                          5,0,3,0,&mat);

503:   /* loop over local fine grid nodes setting interpolation for those*/
504:   for ( j=j_start; j<j_start+n; j++ ) {
505:     for ( i=i_start; i<i_start+m; i++ ) {
506:       /* convert to local "natural" numbering and then to PETSc global numbering */
507:       row    = idx[m_ghost*(j-j_start_ghost) + (i-i_start_ghost)];

509:       i_c = (i/ratio);    /* coarse grid node to left of fine grid node */
510:       j_c = (j/ratio);    /* coarse grid node below fine grid node */

512:       /* 
513:          Only include those interpolation points that are truly 
514:          nonzero. Note this is very important for final grid lines
515:          in x and y directions; since they have no right/top neighbors
516:       */
517:       x  = ((double)(i - i_c*ratio))/((double)ratio);
518:       y  = ((double)(j - j_c*ratio))/((double)ratio);
519:       /* printf("i j %d %d %g %g\n",i,j,x,y); */
520:       nc = 0;
521:       /* one left and below; or we are right on it */
522:       if (j_c < j_start_ghost_c || j_c > j_start_ghost_c+n_ghost_c) {
523:         SETERRQ3(1,"Sorry j %d %d %d",j_c,j_start_ghost_c,j_start_ghost_c+n_ghost_c);
524:       }
525:       if (i_c < i_start_ghost_c || i_c > i_start_ghost_c+m_ghost_c) {
526:         SETERRQ3(1,"Sorry i %d %d %d",i_c,i_start_ghost_c,i_start_ghost_c+m_ghost_c);
527:       }
528:       col      = m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c);
529:       cols[nc] = idx_c[col];
530:       v[nc++]  = x*y - x - y + 1.0;
531:       /* one right and below */
532:       if (i_c*ratio != i) {
533:         cols[nc] = idx_c[col+1];
534:         v[nc++]  = -x*y + x;
535:       }
536:       /* one left and above */
537:       if (j_c*ratio != j) {
538:         cols[nc] = idx_c[col+m_ghost_c];
539:         v[nc++]  = -x*y + y;
540:       }
541:       /* one right and above */
542:       if (j_c*ratio != j && i_c*ratio != i) {
543:         cols[nc] = idx_c[col+m_ghost_c+1];
544:         v[nc++]  = x*y;
545:       }
546:       MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);
547:     }
548:   }
549:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
550:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

552:   VecDuplicate(user->coarse.x,&Rscale);
553:   VecSet(&one,user->fine.x);
554:   MatMultTranspose(mat,user->fine.x,Rscale);
555:   VecReciprocal(Rscale);
556:   user->Rscale = Rscale;
557:   user->R = mat;
558:   return 0;
559: }