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: /*
16: This example models the partial differential equation
18: - Div(alpha* T^beta (GRAD T)) = 0.
20: where beta = 2.5 and alpha = 1.0
22: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
24: in the unit square, which is uniformly discretized in each of x and
25: y in this simple encoding. The degrees of freedom are cell centered.
27: A finite volume approximation with the usual 5-point stencil
28: is used to discretize the boundary value problem to obtain a
29: nonlinear system of equations.
31: This code was contributed by David Keyes
33: */
35: #include <petscsnes.h>
36: #include <petscdm.h>
37: #include <petscdmda.h>
39: /* User-defined application context */
41: typedef struct {
42: PetscReal tleft, tright; /* Dirichlet boundary conditions */
43: PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
44: } AppCtx;
46: #define POWFLOP 5 /* assume a pow() takes five flops */
48: extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
49: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
50: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
52: int main(int argc, char **argv)
53: {
54: SNES snes;
55: AppCtx user;
56: PetscInt its, lits;
57: PetscReal litspit;
58: DM da;
61: PetscInitialize(&argc, &argv, NULL, help);
63: /* set problem parameters */
64: user.tleft = 1.0;
65: user.tright = 0.1;
66: user.beta = 2.5;
67: PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL);
68: PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL);
69: PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL);
70: user.bm1 = user.beta - 1.0;
71: user.coef = user.beta / 2.0;
73: /*
74: Create the multilevel DM data structure
75: */
76: SNESCreate(PETSC_COMM_WORLD, &snes);
78: /*
79: Set the DMDA (grid structure) for the grids.
80: */
81: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da);
82: DMSetFromOptions(da);
83: DMSetUp(da);
84: DMSetApplicationContext(da, &user);
85: SNESSetDM(snes, (DM)da);
87: /*
88: Create the nonlinear solver, and tell it the functions to use
89: */
90: SNESSetFunction(snes, NULL, FormFunction, &user);
91: SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user);
92: SNESSetFromOptions(snes);
93: SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL);
95: SNESSolve(snes, NULL, NULL);
96: SNESGetIterationNumber(snes, &its);
97: SNESGetLinearSolveIterations(snes, &lits);
98: litspit = ((PetscReal)lits) / ((PetscReal)its);
99: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its);
100: PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits);
101: PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit);
103: DMDestroy(&da);
104: SNESDestroy(&snes);
105: PetscFinalize();
106: return 0;
107: }
108: /* -------------------- Form initial approximation ----------------- */
109: PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx)
110: {
111: AppCtx *user;
112: PetscInt i, j, xs, ys, xm, ym;
113: PetscReal tleft;
114: PetscScalar **x;
115: DM da;
118: SNESGetDM(snes, &da);
119: DMGetApplicationContext(da, &user);
120: tleft = user->tleft;
121: /* Get ghost points */
122: DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
123: DMDAVecGetArray(da, X, &x);
125: /* Compute initial guess */
126: for (j = ys; j < ys + ym; j++) {
127: for (i = xs; i < xs + xm; i++) x[j][i] = tleft;
128: }
129: DMDAVecRestoreArray(da, X, &x);
130: return 0;
131: }
132: /* -------------------- Evaluate Function F(x) --------------------- */
133: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
134: {
135: AppCtx *user = (AppCtx *)ptr;
136: PetscInt i, j, mx, my, xs, ys, xm, ym;
137: PetscScalar zero = 0.0, one = 1.0;
138: PetscScalar hx, hy, hxdhy, hydhx;
139: 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;
140: PetscScalar tleft, tright, beta;
141: PetscScalar **x, **f;
142: Vec localX;
143: DM da;
146: SNESGetDM(snes, &da);
147: DMGetLocalVector(da, &localX);
148: DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
149: hx = one / (PetscReal)(mx - 1);
150: hy = one / (PetscReal)(my - 1);
151: hxdhy = hx / hy;
152: hydhx = hy / hx;
153: tleft = user->tleft;
154: tright = user->tright;
155: beta = user->beta;
157: /* Get ghost points */
158: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
159: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
160: DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
161: DMDAVecGetArray(da, localX, &x);
162: DMDAVecGetArray(da, F, &f);
164: /* Evaluate function */
165: for (j = ys; j < ys + ym; j++) {
166: for (i = xs; i < xs + xm; i++) {
167: t0 = x[j][i];
169: if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
170: /* general interior volume */
172: tw = x[j][i - 1];
173: aw = 0.5 * (t0 + tw);
174: dw = PetscPowScalar(aw, beta);
175: fw = dw * (t0 - tw);
177: te = x[j][i + 1];
178: ae = 0.5 * (t0 + te);
179: de = PetscPowScalar(ae, beta);
180: fe = de * (te - t0);
182: ts = x[j - 1][i];
183: as = 0.5 * (t0 + ts);
184: ds = PetscPowScalar(as, beta);
185: fs = ds * (t0 - ts);
187: tn = x[j + 1][i];
188: an = 0.5 * (t0 + tn);
189: dn = PetscPowScalar(an, beta);
190: fn = dn * (tn - t0);
192: } else if (i == 0) {
193: /* left-hand boundary */
194: tw = tleft;
195: aw = 0.5 * (t0 + tw);
196: dw = PetscPowScalar(aw, beta);
197: fw = dw * (t0 - tw);
199: te = x[j][i + 1];
200: ae = 0.5 * (t0 + te);
201: de = PetscPowScalar(ae, beta);
202: fe = de * (te - t0);
204: if (j > 0) {
205: ts = x[j - 1][i];
206: as = 0.5 * (t0 + ts);
207: ds = PetscPowScalar(as, beta);
208: fs = ds * (t0 - ts);
209: } else fs = zero;
211: if (j < my - 1) {
212: tn = x[j + 1][i];
213: an = 0.5 * (t0 + tn);
214: dn = PetscPowScalar(an, beta);
215: fn = dn * (tn - t0);
216: } else fn = zero;
218: } else if (i == mx - 1) {
219: /* right-hand boundary */
220: tw = x[j][i - 1];
221: aw = 0.5 * (t0 + tw);
222: dw = PetscPowScalar(aw, beta);
223: fw = dw * (t0 - tw);
225: te = tright;
226: ae = 0.5 * (t0 + te);
227: de = PetscPowScalar(ae, beta);
228: fe = de * (te - t0);
230: if (j > 0) {
231: ts = x[j - 1][i];
232: as = 0.5 * (t0 + ts);
233: ds = PetscPowScalar(as, beta);
234: fs = ds * (t0 - ts);
235: } else fs = zero;
237: if (j < my - 1) {
238: tn = x[j + 1][i];
239: an = 0.5 * (t0 + tn);
240: dn = PetscPowScalar(an, beta);
241: fn = dn * (tn - t0);
242: } else fn = zero;
244: } else if (j == 0) {
245: /* bottom boundary,and i <> 0 or mx-1 */
246: tw = x[j][i - 1];
247: aw = 0.5 * (t0 + tw);
248: dw = PetscPowScalar(aw, beta);
249: fw = dw * (t0 - tw);
251: te = x[j][i + 1];
252: ae = 0.5 * (t0 + te);
253: de = PetscPowScalar(ae, beta);
254: fe = de * (te - t0);
256: fs = zero;
258: tn = x[j + 1][i];
259: an = 0.5 * (t0 + tn);
260: dn = PetscPowScalar(an, beta);
261: fn = dn * (tn - t0);
263: } else if (j == my - 1) {
264: /* top boundary,and i <> 0 or mx-1 */
265: tw = x[j][i - 1];
266: aw = 0.5 * (t0 + tw);
267: dw = PetscPowScalar(aw, beta);
268: fw = dw * (t0 - tw);
270: te = x[j][i + 1];
271: ae = 0.5 * (t0 + te);
272: de = PetscPowScalar(ae, beta);
273: fe = de * (te - t0);
275: ts = x[j - 1][i];
276: as = 0.5 * (t0 + ts);
277: ds = PetscPowScalar(as, beta);
278: fs = ds * (t0 - ts);
280: fn = zero;
281: }
283: f[j][i] = -hydhx * (fe - fw) - hxdhy * (fn - fs);
284: }
285: }
286: DMDAVecRestoreArray(da, localX, &x);
287: DMDAVecRestoreArray(da, F, &f);
288: DMRestoreLocalVector(da, &localX);
289: PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm);
290: return 0;
291: }
292: /* -------------------- Evaluate Jacobian F(x) --------------------- */
293: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat jac, Mat B, void *ptr)
294: {
295: AppCtx *user = (AppCtx *)ptr;
296: PetscInt i, j, mx, my, xs, ys, xm, ym;
297: PetscScalar one = 1.0, hx, hy, hxdhy, hydhx, t0, tn, ts, te, tw;
298: PetscScalar dn, ds, de, dw, an, as, ae, aw, bn, bs, be, bw, gn, gs, ge, gw;
299: PetscScalar tleft, tright, beta, bm1, coef;
300: PetscScalar v[5], **x;
301: Vec localX;
302: MatStencil col[5], row;
303: DM da;
306: SNESGetDM(snes, &da);
307: DMGetLocalVector(da, &localX);
308: DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
309: hx = one / (PetscReal)(mx - 1);
310: hy = one / (PetscReal)(my - 1);
311: hxdhy = hx / hy;
312: hydhx = hy / hx;
313: tleft = user->tleft;
314: tright = user->tright;
315: beta = user->beta;
316: bm1 = user->bm1;
317: coef = user->coef;
319: /* Get ghost points */
320: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
321: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
322: DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
323: DMDAVecGetArray(da, localX, &x);
325: /* Evaluate Jacobian of function */
326: for (j = ys; j < ys + ym; j++) {
327: for (i = xs; i < xs + xm; i++) {
328: t0 = x[j][i];
330: if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
331: /* general interior volume */
333: tw = x[j][i - 1];
334: aw = 0.5 * (t0 + tw);
335: bw = PetscPowScalar(aw, bm1);
336: /* dw = bw * aw */
337: dw = PetscPowScalar(aw, beta);
338: gw = coef * bw * (t0 - tw);
340: te = x[j][i + 1];
341: ae = 0.5 * (t0 + te);
342: be = PetscPowScalar(ae, bm1);
343: /* de = be * ae; */
344: de = PetscPowScalar(ae, beta);
345: ge = coef * be * (te - t0);
347: ts = x[j - 1][i];
348: as = 0.5 * (t0 + ts);
349: bs = PetscPowScalar(as, bm1);
350: /* ds = bs * as; */
351: ds = PetscPowScalar(as, beta);
352: gs = coef * bs * (t0 - ts);
354: tn = x[j + 1][i];
355: an = 0.5 * (t0 + tn);
356: bn = PetscPowScalar(an, bm1);
357: /* dn = bn * an; */
358: dn = PetscPowScalar(an, beta);
359: gn = coef * bn * (tn - t0);
361: v[0] = -hxdhy * (ds - gs);
362: col[0].j = j - 1;
363: col[0].i = i;
364: v[1] = -hydhx * (dw - gw);
365: col[1].j = j;
366: col[1].i = i - 1;
367: v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
368: col[2].j = row.j = j;
369: col[2].i = row.i = i;
370: v[3] = -hydhx * (de + ge);
371: col[3].j = j;
372: col[3].i = i + 1;
373: v[4] = -hxdhy * (dn + gn);
374: col[4].j = j + 1;
375: col[4].i = i;
376: MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES);
378: } else if (i == 0) {
379: /* left-hand boundary */
380: tw = tleft;
381: aw = 0.5 * (t0 + tw);
382: bw = PetscPowScalar(aw, bm1);
383: /* dw = bw * aw */
384: dw = PetscPowScalar(aw, beta);
385: gw = coef * bw * (t0 - tw);
387: te = x[j][i + 1];
388: ae = 0.5 * (t0 + te);
389: be = PetscPowScalar(ae, bm1);
390: /* de = be * ae; */
391: de = PetscPowScalar(ae, beta);
392: ge = coef * be * (te - t0);
394: /* left-hand bottom boundary */
395: if (j == 0) {
396: tn = x[j + 1][i];
397: an = 0.5 * (t0 + tn);
398: bn = PetscPowScalar(an, bm1);
399: /* dn = bn * an; */
400: dn = PetscPowScalar(an, beta);
401: gn = coef * bn * (tn - t0);
403: v[0] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
404: col[0].j = row.j = j;
405: col[0].i = row.i = i;
406: v[1] = -hydhx * (de + ge);
407: col[1].j = j;
408: col[1].i = i + 1;
409: v[2] = -hxdhy * (dn + gn);
410: col[2].j = j + 1;
411: col[2].i = i;
412: MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES);
414: /* left-hand interior boundary */
415: } else if (j < my - 1) {
416: ts = x[j - 1][i];
417: as = 0.5 * (t0 + ts);
418: bs = PetscPowScalar(as, bm1);
419: /* ds = bs * as; */
420: ds = PetscPowScalar(as, beta);
421: gs = coef * bs * (t0 - ts);
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);
430: v[0] = -hxdhy * (ds - gs);
431: col[0].j = j - 1;
432: col[0].i = i;
433: v[1] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
434: col[1].j = row.j = j;
435: col[1].i = row.i = i;
436: v[2] = -hydhx * (de + ge);
437: col[2].j = j;
438: col[2].i = i + 1;
439: v[3] = -hxdhy * (dn + gn);
440: col[3].j = j + 1;
441: col[3].i = i;
442: MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES);
443: /* left-hand top boundary */
444: } else {
445: ts = x[j - 1][i];
446: as = 0.5 * (t0 + ts);
447: bs = PetscPowScalar(as, bm1);
448: /* ds = bs * as; */
449: ds = PetscPowScalar(as, beta);
450: gs = coef * bs * (t0 - ts);
452: v[0] = -hxdhy * (ds - gs);
453: col[0].j = j - 1;
454: col[0].i = i;
455: v[1] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
456: col[1].j = row.j = j;
457: col[1].i = row.i = i;
458: v[2] = -hydhx * (de + ge);
459: col[2].j = j;
460: col[2].i = i + 1;
461: MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES);
462: }
464: } else if (i == mx - 1) {
465: /* right-hand boundary */
466: tw = x[j][i - 1];
467: aw = 0.5 * (t0 + tw);
468: bw = PetscPowScalar(aw, bm1);
469: /* dw = bw * aw */
470: dw = PetscPowScalar(aw, beta);
471: gw = coef * bw * (t0 - tw);
473: te = tright;
474: ae = 0.5 * (t0 + te);
475: be = PetscPowScalar(ae, bm1);
476: /* de = be * ae; */
477: de = PetscPowScalar(ae, beta);
478: ge = coef * be * (te - t0);
480: /* right-hand bottom boundary */
481: if (j == 0) {
482: tn = x[j + 1][i];
483: an = 0.5 * (t0 + tn);
484: bn = PetscPowScalar(an, bm1);
485: /* dn = bn * an; */
486: dn = PetscPowScalar(an, beta);
487: gn = coef * bn * (tn - t0);
489: v[0] = -hydhx * (dw - gw);
490: col[0].j = j;
491: col[0].i = i - 1;
492: v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
493: col[1].j = row.j = j;
494: col[1].i = row.i = i;
495: v[2] = -hxdhy * (dn + gn);
496: col[2].j = j + 1;
497: col[2].i = i;
498: MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES);
500: /* right-hand interior boundary */
501: } else if (j < my - 1) {
502: ts = x[j - 1][i];
503: as = 0.5 * (t0 + ts);
504: bs = PetscPowScalar(as, bm1);
505: /* ds = bs * as; */
506: ds = PetscPowScalar(as, beta);
507: gs = coef * bs * (t0 - ts);
509: tn = x[j + 1][i];
510: an = 0.5 * (t0 + tn);
511: bn = PetscPowScalar(an, bm1);
512: /* dn = bn * an; */
513: dn = PetscPowScalar(an, beta);
514: gn = coef * bn * (tn - t0);
516: v[0] = -hxdhy * (ds - gs);
517: col[0].j = j - 1;
518: col[0].i = i;
519: v[1] = -hydhx * (dw - gw);
520: col[1].j = j;
521: col[1].i = i - 1;
522: v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
523: col[2].j = row.j = j;
524: col[2].i = row.i = i;
525: v[3] = -hxdhy * (dn + gn);
526: col[3].j = j + 1;
527: col[3].i = i;
528: MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES);
529: /* right-hand top boundary */
530: } else {
531: ts = x[j - 1][i];
532: as = 0.5 * (t0 + ts);
533: bs = PetscPowScalar(as, bm1);
534: /* ds = bs * as; */
535: ds = PetscPowScalar(as, beta);
536: gs = coef * bs * (t0 - ts);
538: v[0] = -hxdhy * (ds - gs);
539: col[0].j = j - 1;
540: col[0].i = i;
541: v[1] = -hydhx * (dw - gw);
542: col[1].j = j;
543: col[1].i = i - 1;
544: v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
545: col[2].j = row.j = j;
546: col[2].i = row.i = i;
547: MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES);
548: }
550: /* bottom boundary,and i <> 0 or mx-1 */
551: } else if (j == 0) {
552: tw = x[j][i - 1];
553: aw = 0.5 * (t0 + tw);
554: bw = PetscPowScalar(aw, bm1);
555: /* dw = bw * aw */
556: dw = PetscPowScalar(aw, beta);
557: gw = coef * bw * (t0 - tw);
559: te = x[j][i + 1];
560: ae = 0.5 * (t0 + te);
561: be = PetscPowScalar(ae, bm1);
562: /* de = be * ae; */
563: de = PetscPowScalar(ae, beta);
564: ge = coef * be * (te - t0);
566: tn = x[j + 1][i];
567: an = 0.5 * (t0 + tn);
568: bn = PetscPowScalar(an, bm1);
569: /* dn = bn * an; */
570: dn = PetscPowScalar(an, beta);
571: gn = coef * bn * (tn - t0);
573: v[0] = -hydhx * (dw - gw);
574: col[0].j = j;
575: col[0].i = i - 1;
576: v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
577: col[1].j = row.j = j;
578: col[1].i = row.i = i;
579: v[2] = -hydhx * (de + ge);
580: col[2].j = j;
581: col[2].i = i + 1;
582: v[3] = -hxdhy * (dn + gn);
583: col[3].j = j + 1;
584: col[3].i = i;
585: MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES);
587: /* top boundary,and i <> 0 or mx-1 */
588: } else if (j == my - 1) {
589: tw = x[j][i - 1];
590: aw = 0.5 * (t0 + tw);
591: bw = PetscPowScalar(aw, bm1);
592: /* dw = bw * aw */
593: dw = PetscPowScalar(aw, beta);
594: gw = coef * bw * (t0 - tw);
596: te = x[j][i + 1];
597: ae = 0.5 * (t0 + te);
598: be = PetscPowScalar(ae, bm1);
599: /* de = be * ae; */
600: de = PetscPowScalar(ae, beta);
601: ge = coef * be * (te - t0);
603: ts = x[j - 1][i];
604: as = 0.5 * (t0 + ts);
605: bs = PetscPowScalar(as, bm1);
606: /* ds = bs * as; */
607: ds = PetscPowScalar(as, beta);
608: gs = coef * bs * (t0 - ts);
610: v[0] = -hxdhy * (ds - gs);
611: col[0].j = j - 1;
612: col[0].i = i;
613: v[1] = -hydhx * (dw - gw);
614: col[1].j = j;
615: col[1].i = i - 1;
616: v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
617: col[2].j = row.j = j;
618: col[2].i = row.i = i;
619: v[3] = -hydhx * (de + ge);
620: col[3].j = j;
621: col[3].i = i + 1;
622: MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES);
623: }
624: }
625: }
626: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
627: DMDAVecRestoreArray(da, localX, &x);
628: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
629: DMRestoreLocalVector(da, &localX);
630: if (jac != B) {
631: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
632: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
633: }
635: PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym);
636: return 0;
637: }
639: /*TEST
641: test:
642: suffix: 1
643: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
644: requires: !single
646: test:
647: suffix: 2
648: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc
649: requires: !single
651: TEST*/