Actual source code: ex11.c

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

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

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

 36: */

 38: #include <petscsnes.h>
 39: #include <petscdmda.h>
 40: #include <petscpcmg.h>

 42: /* User-defined application contexts */

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

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

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


 74: /*
 75:       Mm_ratio - ration of grid lines between fine and coarse grids.
 76: */
 79: int main( int argc, char **argv )
 80: {
 81:   SNES           snes;
 82:   AppCtx         user;
 84:   PetscInt       its, N, n, Nx = PETSC_DECIDE, Ny = PETSC_DECIDE, nlocal,Nlocal;
 85:   PetscMPIInt    size;
 86:   double         bratu_lambda_max = 6.81, bratu_lambda_min = 0.;
 87:   KSP            ksp;
 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) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
113:   n = user.fine.mx*user.fine.my; N = user.coarse.mx*user.coarse.my;

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

119:   /* Set up distributed array for fine grid */
120:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.fine.mx,
121:                     user.fine.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.fine.da);
122:   DMCreateGlobalVector(user.fine.da,&user.fine.x);
123:   VecDuplicate(user.fine.x,&user.fine.r);
124:   VecDuplicate(user.fine.x,&user.fine.b);
125:   VecGetLocalSize(user.fine.x,&nlocal);
126:   DMCreateLocalVector(user.fine.da,&user.fine.localX);
127:   VecDuplicate(user.fine.localX,&user.fine.localF);
128:   MatCreateMPIAIJ(PETSC_COMM_WORLD,nlocal,nlocal,n,n,5,PETSC_NULL,3,PETSC_NULL,&user.fine.J);

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

150:   /* Create nonlinear solver */
151:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

157:   /* set two level additive Schwarz preconditioner */
158:   SNESGetKSP(snes,&ksp);
159:   KSPGetPC(ksp,&pc);
160:   PCSetType(pc,PCMG);
161:   PCMGSetLevels(pc,2,PETSC_NULL);
162:   PCMGSetType(pc,PC_MG_ADDITIVE);

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

169:   /* Create coarse level */
170:   PCMGGetCoarseSolve(pc,&user.ksp_coarse);
171:   KSPSetOptionsPrefix(user.ksp_coarse,"coarse_");
172:   KSPSetFromOptions(user.ksp_coarse);
173:   KSPSetOperators(user.ksp_coarse,user.coarse.J,user.coarse.J,DIFFERENT_NONZERO_PATTERN);
174:   PCMGSetX(pc,COARSE_LEVEL,user.coarse.x);
175:   PCMGSetRhs(pc,COARSE_LEVEL,user.coarse.b);
176:   if (user.redundant_build) {
177:     PC  rpc;
178:     KSPGetPC(user.ksp_coarse,&rpc);
179:     PCRedundantSetScatter(rpc,user.tolocalall,user.fromlocalall);
180:   }

182:   /* Create fine level */
183:   PCMGGetSmoother(pc,FINE_LEVEL,&user.ksp_fine);
184:   KSPSetOptionsPrefix(user.ksp_fine,"fine_");
185:   KSPSetFromOptions(user.ksp_fine);
186:   KSPSetOperators(user.ksp_fine,user.fine.J,user.fine.J,DIFFERENT_NONZERO_PATTERN);
187:   PCMGSetR(pc,FINE_LEVEL,user.fine.r);
188:   PCMGSetResidual(pc,FINE_LEVEL,PCMGDefaultResidual,user.fine.J);

190:   /* Create interpolation between the levels */
191:   FormInterpolation(&user);
192:   PCMGSetInterpolation(pc,FINE_LEVEL,user.R);
193:   PCMGSetRestriction(pc,FINE_LEVEL,user.R);

195:   /* Set options, then solve nonlinear system */
196:   SNESSetFromOptions(snes);
197:   FormInitialGuess1(&user,user.fine.x);
198:   SNESSolve(snes,PETSC_NULL,user.fine.x);
199:   SNESGetIterationNumber(snes,&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:   DMDestroy(&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:   DMDestroy(&user.coarse.da);

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

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

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

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

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

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

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

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

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

290:   /* Get ghost points */
291:   DMGlobalToLocalBegin(user->fine.da,X,INSERT_VALUES,localX);
292:   DMGlobalToLocalEnd(user->fine.da,X,INSERT_VALUES,localX);
293:   DMDAGetCorners(user->fine.da,&xs,&ys,0,&xm,&ym,0);
294:   DMDAGetGhostCorners(user->fine.da,&Xs,&Ys,0,&Xm,&Ym,0);
295:   VecGetArray(localX,&x);
296:   VecGetArray(localF,&f);

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

318:   /* Insert values into global vector */
319:   DMLocalToGlobalBegin(user->fine.da,localF,INSERT_VALUES,F);
320:   DMLocalToGlobalEnd(user->fine.da,localF,INSERT_VALUES,F);
321:   PetscLogFlops(11.0*ym*xm);
322:   return 0;
323: }

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

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

342:   /* Get ghost points */
343:   DMGlobalToLocalBegin(grid->da,X,INSERT_VALUES,localX);
344:   DMGlobalToLocalEnd(grid->da,X,INSERT_VALUES,localX);
345:   DMDAGetCorners(grid->da,&xs,&ys,0,&xm,&ym,0);
346:   DMDAGetGhostCorners(grid->da,&Xs,&Ys,0,&Xm,&Ym,0);
347:   DMDAGetGlobalIndices(grid->da,&nloc,&ltog);
348:   VecGetArray(localX,&x);

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

376:   return 0;
377: }

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

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

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

397:   mx = grid->mx;            my = grid->my;            lambda = user->param;
398:   hx = one/(double)(mx-1);  hy = one/(double)(my-1);
399:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

401:   VecGetArray(X,&x);

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

428:   return 0;
429: }

431: /* --------------------  Evaluate Jacobian F'(x) --------------------- */
434: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
435: {
436:   AppCtx         *user = (AppCtx *) ptr;
438:   KSP            ksp;
439:   PC             pc;
440:   PetscBool      ismg;

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

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

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

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

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

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

472:   return 0;
473: }


478: /*
479:       Forms the interpolation (and restriction) operator from 
480: coarse grid to fine.
481: */
482: PetscErrorCode FormInterpolation(AppCtx *user)
483: {
485:   PetscInt       i,j,i_start,m_fine,j_start,m,n,*idx;
486:   PetscInt       m_ghost,n_ghost,*idx_c,m_ghost_c,n_ghost_c,m_coarse;
487:   PetscInt       row,i_start_ghost,j_start_ghost,cols[4], m_c;
488:   PetscInt       nc,ratio = user->ratio,m_c_local,m_fine_local;
489:   PetscInt       i_c,j_c,i_start_c,j_start_c,n_c,i_start_ghost_c,j_start_ghost_c,col;
490:   PetscScalar    v[4],x,y, one = 1.0;
491:   Mat            mat;
492:   Vec            Rscale;
493: 
494:   DMDAGetCorners(user->fine.da,&i_start,&j_start,0,&m,&n,0);
495:   DMDAGetGhostCorners(user->fine.da,&i_start_ghost,&j_start_ghost,0,&m_ghost,&n_ghost,0);
496:   DMDAGetGlobalIndices(user->fine.da,PETSC_NULL,&idx);

498:   DMDAGetCorners(user->coarse.da,&i_start_c,&j_start_c,0,&m_c,&n_c,0);
499:   DMDAGetGhostCorners(user->coarse.da,&i_start_ghost_c,&j_start_ghost_c,0,&m_ghost_c,&n_ghost_c,0);
500:   DMDAGetGlobalIndices(user->coarse.da,PETSC_NULL,&idx_c);

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

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

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

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

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