Actual source code: ex18.c
1: /* $Id: ex18.c,v 1.23 2001/08/07 21:31:17 bsmith Exp $ */
4: static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 2d.n
5: Uses 2-dimensional distributed arrays.n
6: A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. n
7: n
8: Solves the linear systems via multilevel methods n
9: n
10: The command linen
11: options are:n
12: -tleft <tl>, where <tl> indicates the left Diriclet BC n
13: -tright <tr>, where <tr> indicates the right Diriclet BC n
14: -beta <beta>, where <beta> indicates the exponent in T nn";
16: /*T
17: Concepts: SNES^solving a system of nonlinear equations
18: Concepts: DA^using distributed arrays
19: Concepts: multigrid;
20: Processors: n
21: T*/
23: /*
24:
25: This example models the partial differential equation
26:
27: - Div(alpha* T^beta (GRAD T)) = 0.
28:
29: where beta = 2.5 and alpha = 1.0
30:
31: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
32:
33: in the unit square, which is uniformly discretized in each of x and
34: y in this simple encoding. The degrees of freedom are cell centered.
35:
36: A finite volume approximation with the usual 5-point stencil
37: is used to discretize the boundary value problem to obtain a
38: nonlinear system of equations.
40: This code was contributed by David Keyes
41:
42: */
44: #include petscsnes.h
45: #include petscda.h
46: #include petscmg.h
48: /* User-defined application context */
50: typedef struct {
51: PetscReal tleft,tright; /* Dirichlet boundary conditions */
52: PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */
53: } AppCtx;
55: #define POWFLOP 5 /* assume a pow() takes five flops */
57: extern int FormInitialGuess(SNES,Vec,void*);
58: extern int FormFunction(SNES,Vec,Vec,void*);
59: extern int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
61: #undef __FUNCT__
63: int main(int argc,char **argv)
64: {
65: DMMG *dmmg;
66: SNES snes;
67: AppCtx user;
68: int ierr,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);
102: /*
103: PreLoadBegin() means that the following section of code is run twice. The first time
104: through the flag PreLoading is on this the nonlinear solver is only run for a single step.
105: The second time through (the actually timed code) the maximum iterations is set to 10
106: Preload of the executable is done to eliminate from the timing the time spent bring the
107: executable into memory from disk (paging in).
108: */
109: PreLoadBegin(PETSC_TRUE,"Solve");
110: DMMGSetInitialGuess(dmmg,FormInitialGuess);
111: DMMGSolve(dmmg);
112: PreLoadEnd();
113: snes = DMMGGetSNES(dmmg);
114: SNESGetIterationNumber(snes,&its);
115: SNESGetNumberLinearIterations(snes,&lits);
116: litspit = ((PetscReal)lits)/((PetscReal)its);
117: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %dn",its);
118: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %dn",lits);
119: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %en",litspit);
121: DMMGDestroy(dmmg);
122: PetscFinalize();
124: return 0;
125: }
126: /* -------------------- Form initial approximation ----------------- */
127: #undef __FUNCT__
129: int FormInitialGuess(SNES snes,Vec X,void *ptr)
130: {
131: DMMG dmmg = (DMMG)ptr;
132: AppCtx *user = (AppCtx*)dmmg->user;
133: int i,j,ierr,xs,ys,xm,ym;
134: PetscReal tleft = user->tleft;
135: PetscScalar **x;
139: /* Get ghost points */
140: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
141: DAVecGetArray((DA)dmmg->dm,X,(void**)&x);
143: /* Compute initial guess */
144: for (j=ys; j<ys+ym; j++) {
145: for (i=xs; i<xs+xm; i++) {
146: x[j][i] = tleft;
147: }
148: }
149: DAVecRestoreArray((DA)dmmg->dm,X,(void**)&x);
150: return(0);
151: }
152: /* -------------------- Evaluate Function F(x) --------------------- */
153: #undef __FUNCT__
155: int FormFunction(SNES snes,Vec X,Vec F,void* ptr)
156: {
157: DMMG dmmg = (DMMG)ptr;
158: AppCtx *user = (AppCtx*)dmmg->user;
159: int ierr,i,j,mx,my,xs,ys,xm,ym;
160: PetscScalar zero = 0.0,one = 1.0;
161: PetscScalar hx,hy,hxdhy,hydhx;
162: 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;
163: PetscScalar tleft,tright,beta;
164: PetscScalar **x,**f;
165: Vec localX;
168: DAGetLocalVector((DA)dmmg->dm,&localX);
169: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
170: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
171: hxdhy = hx/hy; hydhx = hy/hx;
172: tleft = user->tleft; tright = user->tright;
173: beta = user->beta;
174:
175: /* Get ghost points */
176: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
177: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
178: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
179: DAVecGetArray((DA)dmmg->dm,localX,(void**)&x);
180: DAVecGetArray((DA)dmmg->dm,F,(void**)&f);
182: /* Evaluate function */
183: for (j=ys; j<ys+ym; j++) {
184: for (i=xs; i<xs+xm; i++) {
185: t0 = x[j][i];
187: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
189: /* general interior volume */
191: tw = x[j][i-1];
192: aw = 0.5*(t0 + tw);
193: dw = PetscPowScalar(aw,beta);
194: fw = dw*(t0 - tw);
196: te = x[j][i+1];
197: ae = 0.5*(t0 + te);
198: de = PetscPowScalar(ae,beta);
199: fe = de*(te - t0);
201: ts = x[j-1][i];
202: as = 0.5*(t0 + ts);
203: ds = PetscPowScalar(as,beta);
204: fs = ds*(t0 - ts);
205:
206: tn = x[j+1][i];
207: an = 0.5*(t0 + tn);
208: dn = PetscPowScalar(an,beta);
209: fn = dn*(tn - t0);
211: } else if (i == 0) {
213: /* left-hand boundary */
214: tw = tleft;
215: aw = 0.5*(t0 + tw);
216: dw = PetscPowScalar(aw,beta);
217: fw = dw*(t0 - tw);
219: te = x[j][i+1];
220: ae = 0.5*(t0 + te);
221: de = PetscPowScalar(ae,beta);
222: fe = de*(te - t0);
224: if (j > 0) {
225: ts = x[j-1][i];
226: as = 0.5*(t0 + ts);
227: ds = PetscPowScalar(as,beta);
228: fs = ds*(t0 - ts);
229: } else {
230: fs = zero;
231: }
233: if (j < my-1) {
234: tn = x[j+1][i];
235: an = 0.5*(t0 + tn);
236: dn = PetscPowScalar(an,beta);
237: fn = dn*(tn - t0);
238: } else {
239: fn = zero;
240: }
242: } else if (i == mx-1) {
244: /* right-hand boundary */
245: tw = x[j][i-1];
246: aw = 0.5*(t0 + tw);
247: dw = PetscPowScalar(aw,beta);
248: fw = dw*(t0 - tw);
249:
250: te = tright;
251: ae = 0.5*(t0 + te);
252: de = PetscPowScalar(ae,beta);
253: fe = de*(te - t0);
254:
255: if (j > 0) {
256: ts = x[j-1][i];
257: as = 0.5*(t0 + ts);
258: ds = PetscPowScalar(as,beta);
259: fs = ds*(t0 - ts);
260: } else {
261: fs = zero;
262: }
263:
264: if (j < my-1) {
265: tn = x[j+1][i];
266: an = 0.5*(t0 + tn);
267: dn = PetscPowScalar(an,beta);
268: fn = dn*(tn - t0);
269: } else {
270: fn = zero;
271: }
273: } else if (j == 0) {
275: /* bottom boundary,and i <> 0 or mx-1 */
276: tw = x[j][i-1];
277: aw = 0.5*(t0 + tw);
278: dw = PetscPowScalar(aw,beta);
279: fw = dw*(t0 - tw);
281: te = x[j][i+1];
282: ae = 0.5*(t0 + te);
283: de = PetscPowScalar(ae,beta);
284: fe = de*(te - t0);
286: fs = zero;
288: tn = x[j+1][i];
289: an = 0.5*(t0 + tn);
290: dn = PetscPowScalar(an,beta);
291: fn = dn*(tn - t0);
293: } else if (j == my-1) {
295: /* top boundary,and i <> 0 or mx-1 */
296: tw = x[j][i-1];
297: aw = 0.5*(t0 + tw);
298: dw = PetscPowScalar(aw,beta);
299: fw = dw*(t0 - tw);
301: te = x[j][i+1];
302: ae = 0.5*(t0 + te);
303: de = PetscPowScalar(ae,beta);
304: fe = de*(te - t0);
306: ts = x[j-1][i];
307: as = 0.5*(t0 + ts);
308: ds = PetscPowScalar(as,beta);
309: fs = ds*(t0 - ts);
311: fn = zero;
313: }
315: f[j][i] = - hydhx*(fe-fw) - hxdhy*(fn-fs);
317: }
318: }
319: DAVecRestoreArray((DA)dmmg->dm,localX,(void**)&x);
320: DAVecRestoreArray((DA)dmmg->dm,F,(void**)&f);
321: DARestoreLocalVector((DA)dmmg->dm,&localX);
322: PetscLogFlops((22 + 4*POWFLOP)*ym*xm);
323: return(0);
324: }
325: /* -------------------- Evaluate Jacobian F(x) --------------------- */
326: #undef __FUNCT__
328: int FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flg,void *ptr)
329: {
330: DMMG dmmg = (DMMG)ptr;
331: AppCtx *user = (AppCtx*)dmmg->user;
332: Mat jac = *J;
333: int ierr,i,j,mx,my,xs,ys,xm,ym;
334: PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw;
335: PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw;
336: PetscScalar tleft,tright,beta,bm1,coef;
337: PetscScalar v[5],**x;
338: Vec localX;
339: MatStencil col[5],row;
342: DAGetLocalVector((DA)dmmg->dm,&localX);
343: *flg = SAME_NONZERO_PATTERN;
344: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
345: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
346: hxdhy = hx/hy; hydhx = hy/hx;
347: tleft = user->tleft; tright = user->tright;
348: beta = user->beta; bm1 = user->bm1; coef = user->coef;
350: /* Get ghost points */
351: DAGlobalToLocalBegin((DA)dmmg->dm,X,INSERT_VALUES,localX);
352: DAGlobalToLocalEnd((DA)dmmg->dm,X,INSERT_VALUES,localX);
353: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
354: DAVecGetArray((DA)dmmg->dm,localX,(void**)&x);
356: /* Evaluate Jacobian of function */
357: for (j=ys; j<ys+ym; j++) {
358: for (i=xs; i<xs+xm; i++) {
359: t0 = x[j][i];
361: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
363: /* general interior volume */
365: tw = x[j][i-1];
366: aw = 0.5*(t0 + tw);
367: bw = PetscPowScalar(aw,bm1);
368: /* dw = bw * aw */
369: dw = PetscPowScalar(aw,beta);
370: gw = coef*bw*(t0 - tw);
372: te = x[j][i+1];
373: ae = 0.5*(t0 + te);
374: be = PetscPowScalar(ae,bm1);
375: /* de = be * ae; */
376: de = PetscPowScalar(ae,beta);
377: ge = coef*be*(te - t0);
379: ts = x[j-1][i];
380: as = 0.5*(t0 + ts);
381: bs = PetscPowScalar(as,bm1);
382: /* ds = bs * as; */
383: ds = PetscPowScalar(as,beta);
384: gs = coef*bs*(t0 - ts);
385:
386: tn = x[j+1][i];
387: an = 0.5*(t0 + tn);
388: bn = PetscPowScalar(an,bm1);
389: /* dn = bn * an; */
390: dn = PetscPowScalar(an,beta);
391: gn = coef*bn*(tn - t0);
393: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
394: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
395: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
396: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
397: v[4] = - hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i;
398: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
400: } else if (i == 0) {
402: /* left-hand boundary */
403: tw = tleft;
404: aw = 0.5*(t0 + tw);
405: bw = PetscPowScalar(aw,bm1);
406: /* dw = bw * aw */
407: dw = PetscPowScalar(aw,beta);
408: gw = coef*bw*(t0 - tw);
409:
410: te = x[j][i + 1];
411: ae = 0.5*(t0 + te);
412: be = PetscPowScalar(ae,bm1);
413: /* de = be * ae; */
414: de = PetscPowScalar(ae,beta);
415: ge = coef*be*(te - t0);
416:
417: /* left-hand bottom boundary */
418: if (j == 0) {
420: tn = x[j+1][i];
421: an = 0.5*(t0 + tn);
422: bn = PetscPowScalar(an,bm1);
423: /* dn = bn * an; */
424: dn = PetscPowScalar(an,beta);
425: gn = coef*bn*(tn - t0);
426:
427: v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i;
428: v[1] = - hydhx*(de + ge); col[1].j = j; col[1].i = i+1;
429: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
430: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
431:
432: /* left-hand interior boundary */
433: } else if (j < my-1) {
435: ts = x[j-1][i];
436: as = 0.5*(t0 + ts);
437: bs = PetscPowScalar(as,bm1);
438: /* ds = bs * as; */
439: ds = PetscPowScalar(as,beta);
440: gs = coef*bs*(t0 - ts);
441:
442: tn = x[j+1][i];
443: an = 0.5*(t0 + tn);
444: bn = PetscPowScalar(an,bm1);
445: /* dn = bn * an; */
446: dn = PetscPowScalar(an,beta);
447: gn = coef*bn*(tn - t0);
448:
449: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
450: v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
451: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
452: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
453: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
454: /* left-hand top boundary */
455: } else {
457: ts = x[j-1][i];
458: as = 0.5*(t0 + ts);
459: bs = PetscPowScalar(as,bm1);
460: /* ds = bs * as; */
461: ds = PetscPowScalar(as,beta);
462: gs = coef*bs*(t0 - ts);
463:
464: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
465: v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
466: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
467: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
468: }
470: } else if (i == mx-1) {
471:
472: /* right-hand boundary */
473: tw = x[j][i-1];
474: aw = 0.5*(t0 + tw);
475: bw = PetscPowScalar(aw,bm1);
476: /* dw = bw * aw */
477: dw = PetscPowScalar(aw,beta);
478: gw = coef*bw*(t0 - tw);
479:
480: te = tright;
481: ae = 0.5*(t0 + te);
482: be = PetscPowScalar(ae,bm1);
483: /* de = be * ae; */
484: de = PetscPowScalar(ae,beta);
485: ge = coef*be*(te - t0);
486:
487: /* right-hand bottom boundary */
488: if (j == 0) {
490: tn = x[j+1][i];
491: an = 0.5*(t0 + tn);
492: bn = PetscPowScalar(an,bm1);
493: /* dn = bn * an; */
494: dn = PetscPowScalar(an,beta);
495: gn = coef*bn*(tn - t0);
496:
497: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
498: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
499: v[2] = - hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i;
500: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
501:
502: /* right-hand interior boundary */
503: } else if (j < my-1) {
505: ts = x[j-1][i];
506: as = 0.5*(t0 + ts);
507: bs = PetscPowScalar(as,bm1);
508: /* ds = bs * as; */
509: ds = PetscPowScalar(as,beta);
510: gs = coef*bs*(t0 - ts);
511:
512: tn = x[j+1][i];
513: an = 0.5*(t0 + tn);
514: bn = PetscPowScalar(an,bm1);
515: /* dn = bn * an; */
516: dn = PetscPowScalar(an,beta);
517: gn = coef*bn*(tn - t0);
518:
519: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
520: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
521: v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
522: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
523: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
524: /* right-hand top boundary */
525: } else {
527: ts = x[j-1][i];
528: as = 0.5*(t0 + ts);
529: bs = PetscPowScalar(as,bm1);
530: /* ds = bs * as; */
531: ds = PetscPowScalar(as,beta);
532: gs = coef*bs*(t0 - ts);
533:
534: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
535: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
536: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
537: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
538: }
540: /* bottom boundary,and i <> 0 or mx-1 */
541: } else if (j == 0) {
543: tw = x[j][i-1];
544: aw = 0.5*(t0 + tw);
545: bw = PetscPowScalar(aw,bm1);
546: /* dw = bw * aw */
547: dw = PetscPowScalar(aw,beta);
548: gw = coef*bw*(t0 - tw);
550: te = x[j][i+1];
551: ae = 0.5*(t0 + te);
552: be = PetscPowScalar(ae,bm1);
553: /* de = be * ae; */
554: de = PetscPowScalar(ae,beta);
555: ge = coef*be*(te - t0);
557: tn = x[j+1][i];
558: an = 0.5*(t0 + tn);
559: bn = PetscPowScalar(an,bm1);
560: /* dn = bn * an; */
561: dn = PetscPowScalar(an,beta);
562: gn = coef*bn*(tn - t0);
563:
564: v[0] = - hydhx*(dw - gw); col[0].j = j; col[0].i = i-1;
565: v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i;
566: v[2] = - hydhx*(de + ge); col[2].j = j; col[2].i = i+1;
567: v[3] = - hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i;
568: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
569:
570: /* top boundary,and i <> 0 or mx-1 */
571: } else if (j == my-1) {
572:
573: tw = x[j][i-1];
574: aw = 0.5*(t0 + tw);
575: bw = PetscPowScalar(aw,bm1);
576: /* dw = bw * aw */
577: dw = PetscPowScalar(aw,beta);
578: gw = coef*bw*(t0 - tw);
580: te = x[j][i+1];
581: ae = 0.5*(t0 + te);
582: be = PetscPowScalar(ae,bm1);
583: /* de = be * ae; */
584: de = PetscPowScalar(ae,beta);
585: ge = coef*be*(te - t0);
587: ts = x[j-1][i];
588: as = 0.5*(t0 + ts);
589: bs = PetscPowScalar(as,bm1);
590: /* ds = bs * as; */
591: ds = PetscPowScalar(as,beta);
592: gs = coef*bs*(t0 - ts);
594: v[0] = - hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i;
595: v[1] = - hydhx*(dw - gw); col[1].j = j; col[1].i = i-1;
596: v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i;
597: v[3] = - hydhx*(de + ge); col[3].j = j; col[3].i = i+1;
598: MatSetValuesStencil(jac,1,&row,4,col,v,INSERT_VALUES);
599:
600: }
601: }
602: }
603: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
604: DAVecRestoreArray((DA)dmmg->dm,localX,(void**)&x);
605: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
606: DARestoreLocalVector((DA)dmmg->dm,&localX);
608: PetscLogFlops((41 + 8*POWFLOP)*xm*ym);
609: return(0);
610: }