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:    KSP        ksp_coarse;
 58:    KSP        ksp_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:   KSP          ksp;
 88:   PC            pc;
 89:   KSP           ksp;

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

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

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

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

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

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

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

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

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

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

160:   /* set two level additive Schwarz preconditioner */
161:   SNESGetKSP(snes,&ksp);
162:   KSPGetPC(ksp,&pc);
163:   PCSetType(pc,PCMG);
164:   MGSetLevels(pc,2,PETSC_NULL);
165:   MGSetType(pc,MGADDITIVE);

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

171:   /* Create coarse level */
172:   MGGetCoarseSolve(pc,&user.ksp_coarse);
173:   KSPSetOptionsPrefix(user.ksp_coarse,"coarse_");
174:   KSPSetFromOptions(user.ksp_coarse);
175:   KSPSetOperators(user.ksp_coarse,user.coarse.J,user.coarse.J,DIFFERENT_NONZERO_PATTERN);
176:   MGSetX(pc,COARSE_LEVEL,user.coarse.x);
177:   MGSetRhs(pc,COARSE_LEVEL,user.coarse.b);
178:   if (user.redundant_build) {
179:     PC  rpc;
180:     KSP rksp;
181:     KSPGetPC(user.ksp_coarse,&rpc);
182:     PCRedundantSetScatter(rpc,user.tolocalall,user.fromlocalall);
183:   }

185:   /* Create fine level */
186:   MGGetSmoother(pc,FINE_LEVEL,&user.ksp_fine);
187:   KSPSetOptionsPrefix(user.ksp_fine,"fine_");
188:   KSPSetFromOptions(user.ksp_fine);
189:   KSPSetOperators(user.ksp_fine,user.fine.J,user.fine.J,DIFFERENT_NONZERO_PATTERN);
190:   MGSetR(pc,FINE_LEVEL,user.fine.r);
191:   MGSetResidual(pc,FINE_LEVEL,MGDefaultResidual,user.fine.J);

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

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

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

215:   MatDestroy(user.fine.J);
216:   VecDestroy(user.fine.x);
217:   VecDestroy(user.fine.r);
218:   VecDestroy(user.fine.b);
219:   DADestroy(user.fine.da);
220:   VecDestroy(user.fine.localX);
221:   VecDestroy(user.fine.localF);

223:   MatDestroy(user.coarse.J);
224:   VecDestroy(user.coarse.x);
225:   VecDestroy(user.coarse.b);
226:   DADestroy(user.coarse.da);

228:   SNESDestroy(snes);
229:   MatDestroy(user.R);
230:   VecDestroy(user.Rscale);
231:   PetscFinalize();

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

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

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

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

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

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

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

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

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

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

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

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

336:   mx = grid->mx;            my = grid->my;            lambda = user->param;
337:   hx = one/(double)(mx-1);  hy = one/(double)(my-1);
338:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

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

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

374:   return 0;
375: }

377: /*
378:         Computes the ENTIRE Jacobian associated with the ENTIRE grid sequentially
379:     This is for generating the coarse grid redundantly.

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

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

394:   mx = grid->mx;            my = grid->my;            lambda = user->param;
395:   hx = one/(double)(mx-1);  hy = one/(double)(my-1);
396:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

398:   VecGetArray(X,&x);

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

425:   return 0;
426: }

428: /* --------------------  Evaluate Jacobian F'(x) --------------------- */
431: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
432: {
433:   AppCtx     *user = (AppCtx *) ptr;
434:   int        ierr;
435:   KSP       ksp;
436:   PC         pc;
437:   PetscTruth ismg;
438:   KSP        ksp;

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

443:   /* create coarse grid jacobian for preconditioner */
444:   SNESGetKSP(snes,&ksp);
445:   KSPGetPC(ksp,&pc);
446: 
447:   PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
448:   if (ismg) {

450:     KSPSetOperators(user->ksp_fine,user->fine.J,user->fine.J,SAME_NONZERO_PATTERN);

452:     /* restrict X to coarse grid */
453:     MatMult(user->R,X,user->coarse.x);
454:     VecPointwiseMult(user->Rscale,user->coarse.x,user->coarse.x);

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

463:     } else {
464:       /* coarse grid Jacobian computed in parallel */
465:       FormJacobian_Grid(user,&user->coarse,user->coarse.x,&user->coarse.J,&user->coarse.J);
466:     }
467:     KSPSetOperators(user->ksp_coarse,user->coarse.J,user->coarse.J,SAME_NONZERO_PATTERN);
468:   }

470:   return 0;
471: }


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

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

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

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

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

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

556:   VecDuplicate(user->coarse.x,&Rscale);
557:   VecSet(&one,user->fine.x);
558:   MatMultTranspose(mat,user->fine.x,Rscale);
559:   VecReciprocal(Rscale);
560:   user->Rscale = Rscale;
561:   user->R = mat;
562:   return 0;
563: }