Actual source code: ex11.c
petsc-3.3-p6 2013-02-11
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
70: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*), FormInitialGuess1(AppCtx*,Vec);
71: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
72: extern PetscErrorCode FormInterpolation(AppCtx *);
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: MatCreateAIJ(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: MatCreateAIJ(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 SNES 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,<og);
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: PetscObjectTypeCompare((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: MatCreateAIJ(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: }