Actual source code: ex18.c
2: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 2d.\n\
3: Uses 2-dimensional distributed arrays.\n\
4: A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
5: \n\
6: Solves the linear systems via multilevel methods \n\
7: \n\
8: The command line\n\
9: options are:\n\
10: -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
11: -tright <tr>, where <tr> indicates the right Diriclet BC \n\
12: -beta <beta>, where <beta> indicates the exponent in T \n\n";
14: /*T
15: Concepts: SNES^solving a system of nonlinear equations
16: Concepts: DA^using distributed arrays
17: Concepts: multigrid;
18: Processors: n
19: T*/
21: /*
22:
23: This example models the partial differential equation
24:
25: - Div(alpha* T^beta (GRAD T)) = 0.
26:
27: where beta = 2.5 and alpha = 1.0
28:
29: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
30:
31: in the unit square, which is uniformly discretized in each of x and
32: y in this simple encoding. The degrees of freedom are cell centered.
33:
34: A finite volume approximation with the usual 5-point stencil
35: is used to discretize the boundary value problem to obtain a
36: nonlinear system of equations.
38: This code was contributed by David Keyes
39:
40: */
42: #include petscsnes.h
43: #include petscda.h
44: #include petscmg.h
45: #include petscdmmg.h
47: /* User-defined application context */
49: typedef struct {
50: PetscReal tleft,tright; /* Dirichlet boundary conditions */
51: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
52: } AppCtx;
54: #define POWFLOP 5 /* assume a pow() takes five flops */
62: int main(int argc,char **argv)
63: {
64: DMMG *dmmg;
65: SNES snes;
66: AppCtx user;
68: PetscInt its,lits;
69: PetscReal litspit;
70: DA da;
72: PetscInitialize(&argc,&argv,PETSC_NULL,help);
74: /* set problem parameters */
75: user.tleft = 1.0;
76: user.tright = 0.1;
77: user.beta = 2.5;
78: PetscOptionsGetReal(PETSC_NULL,"-tleft",&user.tleft,PETSC_NULL);
79: PetscOptionsGetReal(PETSC_NULL,"-tright",&user.tright,PETSC_NULL);
80: PetscOptionsGetReal(PETSC_NULL,"-beta",&user.beta,PETSC_NULL);
81: user.bm1 = user.beta - 1.0;
82: user.coef = user.beta/2.0;
85: /*
86: Create the multilevel DA data structure
87: */
88: DMMGCreate(PETSC_COMM_WORLD,3,&user,&dmmg);
90: /*
91: Set the DA (grid structure) for the grids.
92: */
93: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
94: DMMGSetDM(dmmg,(DM)da);
95: DADestroy(da);
97: /*
98: Create the nonlinear solver, and tell the DMMG structure to use it
99: */
100: DMMGSetSNES(dmmg,FormFunction,FormJacobian);
101: DMMGSetFromOptions(dmmg);
103: /*
104: PreLoadBegin() means that the following section of code is run twice. The first time
105: through the flag PreLoading is on this the nonlinear solver is only run for a single step.
106: The second time through (the actually timed code) the maximum iterations is set to 10
107: Preload of the executable is done to eliminate from the timing the time spent bring the
108: executable into memory from disk (paging in).
109: */
110: PreLoadBegin(PETSC_TRUE,"Solve");
111: DMMGSetInitialGuess(dmmg,FormInitialGuess);
112: DMMGSolve(dmmg);
113: PreLoadEnd();
114: snes = DMMGGetSNES(dmmg);
115: SNESGetIterationNumber(snes,&its);
116: SNESGetLinearSolveIterations(snes,&lits);
117: litspit = ((PetscReal)lits)/((PetscReal)its);
118: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
119: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
120: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);
122: DMMGDestroy(dmmg);
123: PetscFinalize();
125: return 0;
126: }
127: /* -------------------- Form initial approximation ----------------- */
130: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
131: {
132: AppCtx *user = (AppCtx*)dmmg->user;
133: PetscInt i,j,xs,ys,xm,ym;
135: PetscReal tleft = user->tleft;
136: PetscScalar **x;
140: /* Get ghost points */
141: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
142: DAVecGetArray((DA)dmmg->dm,X,&x);
144: /* Compute initial guess */
145: for (j=ys; j<ys+ym; j++) {
146: for (i=xs; i<xs+xm; i++) {
147: x[j][i] = tleft;
148: }
149: }
150: DAVecRestoreArray((DA)dmmg->dm,X,&x);
151: return(0);
152: }
153: /* -------------------- Evaluate Function F(x) --------------------- */
156: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void* ptr)
157: {
158: DMMG dmmg = (DMMG)ptr;
159: AppCtx *user = (AppCtx*)dmmg->user;
161: PetscInt i,j,mx,my,xs,ys,xm,ym;
162: PetscScalar zero = 0.0,one = 1.0;
163: PetscScalar hx,hy,hxdhy,hydhx;
164: PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0;
165: PetscScalar tleft,tright,beta;
166: PetscScalar **x,**f;
167: Vec localX;
170: DAGetLocalVector((DA)dmmg->dm,&localX);
171: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
172: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
173: hxdhy = hx/hy; hydhx = hy/hx;
174: tleft = user->tleft; tright = user->tright;
175: beta = user->beta;
176:
177: /* Get ghost points */
178: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
179: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
180: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
181: DAVecGetArray((DA)dmmg->dm,localX,&x);
182: DAVecGetArray((DA)dmmg->dm,F,&f);
184: /* Evaluate function */
185: for (j=ys; j<ys+ym; j++) {
186: for (i=xs; i<xs+xm; i++) {
187: t0 = x[j][i];
189: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
191: /* general interior volume */
193: tw = x[j][i-1];
194: aw = 0.5*(t0 + tw);
195: dw = PetscPowScalar(aw,beta);
196: fw = dw*(t0 - tw);
198: te = x[j][i+1];
199: ae = 0.5*(t0 + te);
200: de = PetscPowScalar(ae,beta);
201: fe = de*(te - t0);
203: ts = x[j-1][i];
204: as = 0.5*(t0 + ts);
205: ds = PetscPowScalar(as,beta);
206: fs = ds*(t0 - ts);
207:
208: tn = x[j+1][i];
209: an = 0.5*(t0 + tn);
210: dn = PetscPowScalar(an,beta);
211: fn = dn*(tn - t0);
213: } else if (i == 0) {
215: /* left-hand boundary */
216: tw = tleft;
217: aw = 0.5*(t0 + tw);
218: dw = PetscPowScalar(aw,beta);
219: fw = dw*(t0 - tw);
221: te = x[j][i+1];
222: ae = 0.5*(t0 + te);
223: de = PetscPowScalar(ae,beta);
224: fe = de*(te - t0);
226: if (j > 0) {
227: ts = x[j-1][i];
228: as = 0.5*(t0 + ts);
229: ds = PetscPowScalar(as,beta);
230: fs = ds*(t0 - ts);
231: } else {
232: fs = zero;
233: }
235: if (j < my-1) {
236: tn = x[j+1][i];
237: an = 0.5*(t0 + tn);
238: dn = PetscPowScalar(an,beta);
239: fn = dn*(tn - t0);
240: } else {
241: fn = zero;
242: }
244: } else if (i == mx-1) {
246: /* right-hand boundary */
247: tw = x[j][i-1];
248: aw = 0.5*(t0 + tw);
249: dw = PetscPowScalar(aw,beta);
250: fw = dw*(t0 - tw);
251:
252: te = tright;
253: ae = 0.5*(t0 + te);
254: de = PetscPowScalar(ae,beta);
255: fe = de*(te - t0);
256:
257: if (j > 0) {
258: ts = x[j-1][i];
259: as = 0.5*(t0 + ts);
260: ds = PetscPowScalar(as,beta);
261: fs = ds*(t0 - ts);
262: } else {
263: fs = zero;
264: }
265:
266: if (j < my-1) {
267: tn = x[j+1][i];
268: an = 0.5*(t0 + tn);
269: dn = PetscPowScalar(an,beta);
270: fn = dn*(tn - t0);
271: } else {
272: fn = zero;
273: }
275: } else if (j == 0) {
277: /* bottom boundary,and i <> 0 or mx-1 */
278: tw = x[j][i-1];
279: aw = 0.5*(t0 + tw);
280: dw = PetscPowScalar(aw,beta);
281: fw = dw*(t0 - tw);
283: te = x[j][i+1];
284: ae = 0.5*(t0 + te);
285: de = PetscPowScalar(ae,beta);
286: fe = de*(te - t0);
288: fs = zero;
290: tn = x[j+1][i];
291: an = 0.5*(t0 + tn);
292: dn = PetscPowScalar(an,beta);
293: fn = dn*(tn - t0);
295: } else if (j == my-1) {
297: /* top boundary,and i <> 0 or mx-1 */
298: tw = x[j][i-1];
299: aw = 0.5*(t0 + tw);
300: dw = PetscPowScalar(aw,beta);
301: fw = dw*(t0 - tw);
303: te = x[j][i+1];
304: ae = 0.5*(t0 + te);
305: de = PetscPowScalar(ae,beta);
306: fe = de*(te - t0);
308: ts = x[j-1][i];
309: as = 0.5*(t0 + ts);
310: ds = PetscPowScalar(as,beta);
311: fs = ds*(t0 - ts);
313: fn = zero;
315: }
317: f[j][i] = - hydhx*(fe-fw) - hxdhy*(fn-fs);
319: }
320: }
321: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
322: DAVecRestoreArray((DA)dmmg->dm,F,&f);
323: DARestoreLocalVector((DA)dmmg->dm,&localX);
324: PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm);
325: return(0);
326: }
327: /* -------------------- Evaluate Jacobian F(x) --------------------- */
330: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
331: {
332: DMMG dmmg = (DMMG)ptr;
333: AppCtx *user = (AppCtx*)dmmg->user;
334: Mat jac = *J;
336: PetscInt i,j,mx,my,xs,ys,xm,ym;
337: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
338: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
339: PetscScalar tleft,tright,beta,bm1,coef;
340: PetscScalar v[5],**x;
341: Vec localX;
342: MatStencil col[5],row;
345: DAGetLocalVector((DA)dmmg->dm,&localX);
346: *flg = SAME_NONZERO_PATTERN;
347: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
348: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
349: hxdhy = hx/hy; hydhx = hy/hx;
350: tleft = user->tleft; tright = user->tright;
351: beta = user->beta; bm1 = user->bm1; coef = user->coef;
353: /* Get ghost points */
354: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
355: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
356: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
357: DAVecGetArray((DA)dmmg->dm,localX,&x);
359: /* Evaluate Jacobian of function */
360: for (j=ys; j<ys+ym; j++) {
361: for (i=xs; i<xs+xm; i++) {
362: t0 = x[j][i];
364: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
366: /* general interior volume */
368: tw = x[j][i-1];
369: aw = 0.5*(t0 + tw);
370: bw = PetscPowScalar(aw,bm1);
371: /* dw = bw * aw */
372: dw = PetscPowScalar(aw,beta);
373: gw = coef*bw*(t0 - tw);
375: te = x[j][i+1];
376: ae = 0.5*(t0 + te);
377: be = PetscPowScalar(ae,bm1);
378: /* de = be * ae; */
379: de = PetscPowScalar(ae,beta);
380: ge = coef*be*(te - t0);
382: ts = x[j-1][i];
383: as = 0.5*(t0 + ts);
384: bs = PetscPowScalar(as,bm1);
385: /* ds = bs * as; */
386: ds = PetscPowScalar(as,beta);
387: gs = coef*bs*(t0 - ts);
388:
389: tn = x[j+1][i];
390: an = 0.5*(t0 + tn);
391: bn = PetscPowScalar(an,bm1);
392: /* dn = bn * an; */
393: dn = PetscPowScalar(an,beta);
394: gn = coef*bn*(tn - t0);
396: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
397: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
398: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
399: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
400: v[4] = - hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
401: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
403: } else if (i == 0) {
405: /* left-hand boundary */
406: tw = tleft;
407: aw = 0.5*(t0 + tw);
408: bw = PetscPowScalar(aw,bm1);
409: /* dw = bw * aw */
410: dw = PetscPowScalar(aw,beta);
411: gw = coef*bw*(t0 - tw);
412:
413: te = x[j][i + 1];
414: ae = 0.5*(t0 + te);
415: be = PetscPowScalar(ae,bm1);
416: /* de = be * ae; */
417: de = PetscPowScalar(ae,beta);
418: ge = coef*be*(te - t0);
419:
420: /* left-hand bottom boundary */
421: if (j == 0) {
423: tn = x[j+1][i];
424: an = 0.5*(t0 + tn);
425: bn = PetscPowScalar(an,bm1);
426: /* dn = bn * an; */
427: dn = PetscPowScalar(an,beta);
428: gn = coef*bn*(tn - t0);
429:
430: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
431: v[1] = - hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
432: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
433: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
434:
435: /* left-hand interior boundary */
436: } else if (j < my-1) {
438: ts = x[j-1][i];
439: as = 0.5*(t0 + ts);
440: bs = PetscPowScalar(as,bm1);
441: /* ds = bs * as; */
442: ds = PetscPowScalar(as,beta);
443: gs = coef*bs*(t0 - ts);
444:
445: tn = x[j+1][i];
446: an = 0.5*(t0 + tn);
447: bn = PetscPowScalar(an,bm1);
448: /* dn = bn * an; */
449: dn = PetscPowScalar(an,beta);
450: gn = coef*bn*(tn - t0);
451:
452: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
453: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
454: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
455: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
456: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
457: /* left-hand top boundary */
458: } else {
460: ts = x[j-1][i];
461: as = 0.5*(t0 + ts);
462: bs = PetscPowScalar(as,bm1);
463: /* ds = bs * as; */
464: ds = PetscPowScalar(as,beta);
465: gs = coef*bs*(t0 - ts);
466:
467: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
468: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
469: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
470: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
471: }
473: } else if (i == mx-1) {
474:
475: /* right-hand boundary */
476: tw = x[j][i-1];
477: aw = 0.5*(t0 + tw);
478: bw = PetscPowScalar(aw,bm1);
479: /* dw = bw * aw */
480: dw = PetscPowScalar(aw,beta);
481: gw = coef*bw*(t0 - tw);
482:
483: te = tright;
484: ae = 0.5*(t0 + te);
485: be = PetscPowScalar(ae,bm1);
486: /* de = be * ae; */
487: de = PetscPowScalar(ae,beta);
488: ge = coef*be*(te - t0);
489:
490: /* right-hand bottom boundary */
491: if (j == 0) {
493: tn = x[j+1][i];
494: an = 0.5*(t0 + tn);
495: bn = PetscPowScalar(an,bm1);
496: /* dn = bn * an; */
497: dn = PetscPowScalar(an,beta);
498: gn = coef*bn*(tn - t0);
499:
500: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
501: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
502: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
503: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
504:
505: /* right-hand interior boundary */
506: } else if (j < my-1) {
508: ts = x[j-1][i];
509: as = 0.5*(t0 + ts);
510: bs = PetscPowScalar(as,bm1);
511: /* ds = bs * as; */
512: ds = PetscPowScalar(as,beta);
513: gs = coef*bs*(t0 - ts);
514:
515: tn = x[j+1][i];
516: an = 0.5*(t0 + tn);
517: bn = PetscPowScalar(an,bm1);
518: /* dn = bn * an; */
519: dn = PetscPowScalar(an,beta);
520: gn = coef*bn*(tn - t0);
521:
522: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
523: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
524: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
525: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
526: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
527: /* right-hand top boundary */
528: } else {
530: ts = x[j-1][i];
531: as = 0.5*(t0 + ts);
532: bs = PetscPowScalar(as,bm1);
533: /* ds = bs * as; */
534: ds = PetscPowScalar(as,beta);
535: gs = coef*bs*(t0 - ts);
536:
537: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
538: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
539: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
540: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
541: }
543: /* bottom boundary,and i <> 0 or mx-1 */
544: } else if (j == 0) {
546: tw = x[j][i-1];
547: aw = 0.5*(t0 + tw);
548: bw = PetscPowScalar(aw,bm1);
549: /* dw = bw * aw */
550: dw = PetscPowScalar(aw,beta);
551: gw = coef*bw*(t0 - tw);
553: te = x[j][i+1];
554: ae = 0.5*(t0 + te);
555: be = PetscPowScalar(ae,bm1);
556: /* de = be * ae; */
557: de = PetscPowScalar(ae,beta);
558: ge = coef*be*(te - t0);
560: tn = x[j+1][i];
561: an = 0.5*(t0 + tn);
562: bn = PetscPowScalar(an,bm1);
563: /* dn = bn * an; */
564: dn = PetscPowScalar(an,beta);
565: gn = coef*bn*(tn - t0);
566:
567: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
568: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
569: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
570: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
571: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
572:
573: /* top boundary,and i <> 0 or mx-1 */
574: } else if (j == my-1) {
575:
576: tw = x[j][i-1];
577: aw = 0.5*(t0 + tw);
578: bw = PetscPowScalar(aw,bm1);
579: /* dw = bw * aw */
580: dw = PetscPowScalar(aw,beta);
581: gw = coef*bw*(t0 - tw);
583: te = x[j][i+1];
584: ae = 0.5*(t0 + te);
585: be = PetscPowScalar(ae,bm1);
586: /* de = be * ae; */
587: de = PetscPowScalar(ae,beta);
588: ge = coef*be*(te - t0);
590: ts = x[j-1][i];
591: as = 0.5*(t0 + ts);
592: bs = PetscPowScalar(as,bm1);
593: /* ds = bs * as; */
594: ds = PetscPowScalar(as,beta);
595: gs = coef*bs*(t0 - ts);
597: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
598: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
599: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
600: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
601: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
602:
603: }
604: }
605: }
606: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
607: DAVecRestoreArray((DA)dmmg->dm,localX,&x);
608: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
609: DARestoreLocalVector((DA)dmmg->dm,&localX);
611: PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym);
612: return(0);
613: }