Actual source code: ex18.c
1: static char help[] = "Nonlinear Radiative Transport PDE with multigrid in 2d.\n\
2: Uses 2-dimensional distributed arrays.\n\
3: A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
4: \n\
5: Solves the linear systems via multilevel methods \n\
6: \n\
7: The command line\n\
8: options are:\n\
9: -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
10: -tright <tr>, where <tr> indicates the right Diriclet BC \n\
11: -beta <beta>, where <beta> indicates the exponent in T \n\n";
13: /*
15: This example models the partial differential equation
17: - Div(alpha* T^beta (GRAD T)) = 0.
19: where beta = 2.5 and alpha = 1.0
21: BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
23: in the unit square, which is uniformly discretized in each of x and
24: y in this simple encoding. The degrees of freedom are cell centered.
26: A finite volume approximation with the usual 5-point stencil
27: is used to discretize the boundary value problem to obtain a
28: nonlinear system of equations.
30: This code was contributed by David Keyes
32: */
34: #include <petscsnes.h>
35: #include <petscdm.h>
36: #include <petscdmda.h>
38: /* User-defined application context */
40: typedef struct {
41: PetscReal tleft, tright; /* Dirichlet boundary conditions */
42: PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
43: } AppCtx;
45: #define POWFLOP 5 /* assume a pow() takes five flops */
47: extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
48: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
49: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
51: int main(int argc, char **argv)
52: {
53: SNES snes;
54: AppCtx user;
55: PetscInt its, lits;
56: PetscReal litspit;
57: DM da;
59: PetscFunctionBeginUser;
60: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
62: /* set problem parameters */
63: user.tleft = 1.0;
64: user.tright = 0.1;
65: user.beta = 2.5;
66: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL));
67: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL));
68: PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL));
69: user.bm1 = user.beta - 1.0;
70: user.coef = user.beta / 2.0;
72: /*
73: Create the multilevel DM data structure
74: */
75: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
77: /*
78: Set the DMDA (grid structure) for the grids.
79: */
80: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da));
81: PetscCall(DMSetFromOptions(da));
82: PetscCall(DMSetUp(da));
83: PetscCall(DMSetApplicationContext(da, &user));
84: PetscCall(SNESSetDM(snes, (DM)da));
86: /*
87: Create the nonlinear solver, and tell it the functions to use
88: */
89: PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user));
90: PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user));
91: PetscCall(SNESSetFromOptions(snes));
92: PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL));
94: PetscCall(SNESSolve(snes, NULL, NULL));
95: PetscCall(SNESGetIterationNumber(snes, &its));
96: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
97: litspit = ((PetscReal)lits) / ((PetscReal)its);
98: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
99: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits));
100: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit));
102: PetscCall(DMDestroy(&da));
103: PetscCall(SNESDestroy(&snes));
104: PetscCall(PetscFinalize());
105: return 0;
106: }
107: /* -------------------- Form initial approximation ----------------- */
108: PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx)
109: {
110: AppCtx *user;
111: PetscInt i, j, xs, ys, xm, ym;
112: PetscReal tleft;
113: PetscScalar **x;
114: DM da;
116: PetscFunctionBeginUser;
117: PetscCall(SNESGetDM(snes, &da));
118: PetscCall(DMGetApplicationContext(da, &user));
119: tleft = user->tleft;
120: /* Get ghost points */
121: PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
122: PetscCall(DMDAVecGetArray(da, X, &x));
124: /* Compute initial guess */
125: for (j = ys; j < ys + ym; j++) {
126: for (i = xs; i < xs + xm; i++) x[j][i] = tleft;
127: }
128: PetscCall(DMDAVecRestoreArray(da, X, &x));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
131: /* -------------------- Evaluate Function F(x) --------------------- */
132: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
133: {
134: AppCtx *user = (AppCtx *)ptr;
135: PetscInt i, j, mx, my, xs, ys, xm, ym;
136: PetscScalar zero = 0.0, one = 1.0;
137: PetscScalar hx, hy, hxdhy, hydhx;
138: 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;
139: PetscScalar tleft, tright, beta;
140: PetscScalar **x, **f;
141: Vec localX;
142: DM da;
144: PetscFunctionBeginUser;
145: PetscCall(SNESGetDM(snes, &da));
146: PetscCall(DMGetLocalVector(da, &localX));
147: PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
148: hx = one / (PetscReal)(mx - 1);
149: hy = one / (PetscReal)(my - 1);
150: hxdhy = hx / hy;
151: hydhx = hy / hx;
152: tleft = user->tleft;
153: tright = user->tright;
154: beta = user->beta;
156: /* Get ghost points */
157: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
158: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
159: PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
160: PetscCall(DMDAVecGetArray(da, localX, &x));
161: PetscCall(DMDAVecGetArray(da, F, &f));
163: /* Evaluate function */
164: for (j = ys; j < ys + ym; j++) {
165: for (i = xs; i < xs + xm; i++) {
166: t0 = x[j][i];
168: if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
169: /* general interior volume */
171: tw = x[j][i - 1];
172: aw = 0.5 * (t0 + tw);
173: dw = PetscPowScalar(aw, beta);
174: fw = dw * (t0 - tw);
176: te = x[j][i + 1];
177: ae = 0.5 * (t0 + te);
178: de = PetscPowScalar(ae, beta);
179: fe = de * (te - t0);
181: ts = x[j - 1][i];
182: as = 0.5 * (t0 + ts);
183: ds = PetscPowScalar(as, beta);
184: fs = ds * (t0 - ts);
186: tn = x[j + 1][i];
187: an = 0.5 * (t0 + tn);
188: dn = PetscPowScalar(an, beta);
189: fn = dn * (tn - t0);
191: } else if (i == 0) {
192: /* left-hand boundary */
193: tw = tleft;
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: if (j > 0) {
204: ts = x[j - 1][i];
205: as = 0.5 * (t0 + ts);
206: ds = PetscPowScalar(as, beta);
207: fs = ds * (t0 - ts);
208: } else fs = zero;
210: if (j < my - 1) {
211: tn = x[j + 1][i];
212: an = 0.5 * (t0 + tn);
213: dn = PetscPowScalar(an, beta);
214: fn = dn * (tn - t0);
215: } else fn = zero;
217: } else if (i == mx - 1) {
218: /* right-hand boundary */
219: tw = x[j][i - 1];
220: aw = 0.5 * (t0 + tw);
221: dw = PetscPowScalar(aw, beta);
222: fw = dw * (t0 - tw);
224: te = tright;
225: ae = 0.5 * (t0 + te);
226: de = PetscPowScalar(ae, beta);
227: fe = de * (te - t0);
229: if (j > 0) {
230: ts = x[j - 1][i];
231: as = 0.5 * (t0 + ts);
232: ds = PetscPowScalar(as, beta);
233: fs = ds * (t0 - ts);
234: } else fs = zero;
236: if (j < my - 1) {
237: tn = x[j + 1][i];
238: an = 0.5 * (t0 + tn);
239: dn = PetscPowScalar(an, beta);
240: fn = dn * (tn - t0);
241: } else fn = zero;
243: } else if (j == 0) {
244: /* bottom boundary,and i <> 0 or mx-1 */
245: tw = x[j][i - 1];
246: aw = 0.5 * (t0 + tw);
247: dw = PetscPowScalar(aw, beta);
248: fw = dw * (t0 - tw);
250: te = x[j][i + 1];
251: ae = 0.5 * (t0 + te);
252: de = PetscPowScalar(ae, beta);
253: fe = de * (te - t0);
255: fs = zero;
257: tn = x[j + 1][i];
258: an = 0.5 * (t0 + tn);
259: dn = PetscPowScalar(an, beta);
260: fn = dn * (tn - t0);
262: } else if (j == my - 1) {
263: /* top boundary,and i <> 0 or mx-1 */
264: tw = x[j][i - 1];
265: aw = 0.5 * (t0 + tw);
266: dw = PetscPowScalar(aw, beta);
267: fw = dw * (t0 - tw);
269: te = x[j][i + 1];
270: ae = 0.5 * (t0 + te);
271: de = PetscPowScalar(ae, beta);
272: fe = de * (te - t0);
274: ts = x[j - 1][i];
275: as = 0.5 * (t0 + ts);
276: ds = PetscPowScalar(as, beta);
277: fs = ds * (t0 - ts);
279: fn = zero;
280: }
282: f[j][i] = -hydhx * (fe - fw) - hxdhy * (fn - fs);
283: }
284: }
285: PetscCall(DMDAVecRestoreArray(da, localX, &x));
286: PetscCall(DMDAVecRestoreArray(da, F, &f));
287: PetscCall(DMRestoreLocalVector(da, &localX));
288: PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm));
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
291: /* -------------------- Evaluate Jacobian F(x) --------------------- */
292: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat jac, Mat B, void *ptr)
293: {
294: AppCtx *user = (AppCtx *)ptr;
295: PetscInt i, j, mx, my, xs, ys, xm, ym;
296: PetscScalar one = 1.0, hx, hy, hxdhy, hydhx, t0, tn, ts, te, tw;
297: PetscScalar dn, ds, de, dw, an, as, ae, aw, bn, bs, be, bw, gn, gs, ge, gw;
298: PetscScalar tleft, tright, beta, bm1, coef;
299: PetscScalar v[5], **x;
300: Vec localX;
301: MatStencil col[5], row;
302: DM da;
304: PetscFunctionBeginUser;
305: PetscCall(SNESGetDM(snes, &da));
306: PetscCall(DMGetLocalVector(da, &localX));
307: PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
308: hx = one / (PetscReal)(mx - 1);
309: hy = one / (PetscReal)(my - 1);
310: hxdhy = hx / hy;
311: hydhx = hy / hx;
312: tleft = user->tleft;
313: tright = user->tright;
314: beta = user->beta;
315: bm1 = user->bm1;
316: coef = user->coef;
318: /* Get ghost points */
319: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
320: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
321: PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
322: PetscCall(DMDAVecGetArray(da, localX, &x));
324: /* Evaluate Jacobian of function */
325: for (j = ys; j < ys + ym; j++) {
326: for (i = xs; i < xs + xm; i++) {
327: t0 = x[j][i];
329: if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
330: /* general interior volume */
332: tw = x[j][i - 1];
333: aw = 0.5 * (t0 + tw);
334: bw = PetscPowScalar(aw, bm1);
335: /* dw = bw * aw */
336: dw = PetscPowScalar(aw, beta);
337: gw = coef * bw * (t0 - tw);
339: te = x[j][i + 1];
340: ae = 0.5 * (t0 + te);
341: be = PetscPowScalar(ae, bm1);
342: /* de = be * ae; */
343: de = PetscPowScalar(ae, beta);
344: ge = coef * be * (te - t0);
346: ts = x[j - 1][i];
347: as = 0.5 * (t0 + ts);
348: bs = PetscPowScalar(as, bm1);
349: /* ds = bs * as; */
350: ds = PetscPowScalar(as, beta);
351: gs = coef * bs * (t0 - ts);
353: tn = x[j + 1][i];
354: an = 0.5 * (t0 + tn);
355: bn = PetscPowScalar(an, bm1);
356: /* dn = bn * an; */
357: dn = PetscPowScalar(an, beta);
358: gn = coef * bn * (tn - t0);
360: v[0] = -hxdhy * (ds - gs);
361: col[0].j = j - 1;
362: col[0].i = i;
363: v[1] = -hydhx * (dw - gw);
364: col[1].j = j;
365: col[1].i = i - 1;
366: v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
367: col[2].j = row.j = j;
368: col[2].i = row.i = i;
369: v[3] = -hydhx * (de + ge);
370: col[3].j = j;
371: col[3].i = i + 1;
372: v[4] = -hxdhy * (dn + gn);
373: col[4].j = j + 1;
374: col[4].i = i;
375: PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES));
377: } else if (i == 0) {
378: /* left-hand boundary */
379: tw = tleft;
380: aw = 0.5 * (t0 + tw);
381: bw = PetscPowScalar(aw, bm1);
382: /* dw = bw * aw */
383: dw = PetscPowScalar(aw, beta);
384: gw = coef * bw * (t0 - tw);
386: te = x[j][i + 1];
387: ae = 0.5 * (t0 + te);
388: be = PetscPowScalar(ae, bm1);
389: /* de = be * ae; */
390: de = PetscPowScalar(ae, beta);
391: ge = coef * be * (te - t0);
393: /* left-hand bottom boundary */
394: if (j == 0) {
395: tn = x[j + 1][i];
396: an = 0.5 * (t0 + tn);
397: bn = PetscPowScalar(an, bm1);
398: /* dn = bn * an; */
399: dn = PetscPowScalar(an, beta);
400: gn = coef * bn * (tn - t0);
402: v[0] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
403: col[0].j = row.j = j;
404: col[0].i = row.i = i;
405: v[1] = -hydhx * (de + ge);
406: col[1].j = j;
407: col[1].i = i + 1;
408: v[2] = -hxdhy * (dn + gn);
409: col[2].j = j + 1;
410: col[2].i = i;
411: PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
413: /* left-hand interior boundary */
414: } else if (j < my - 1) {
415: ts = x[j - 1][i];
416: as = 0.5 * (t0 + ts);
417: bs = PetscPowScalar(as, bm1);
418: /* ds = bs * as; */
419: ds = PetscPowScalar(as, beta);
420: gs = coef * bs * (t0 - ts);
422: tn = x[j + 1][i];
423: an = 0.5 * (t0 + tn);
424: bn = PetscPowScalar(an, bm1);
425: /* dn = bn * an; */
426: dn = PetscPowScalar(an, beta);
427: gn = coef * bn * (tn - t0);
429: v[0] = -hxdhy * (ds - gs);
430: col[0].j = j - 1;
431: col[0].i = i;
432: v[1] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
433: col[1].j = row.j = j;
434: col[1].i = row.i = i;
435: v[2] = -hydhx * (de + ge);
436: col[2].j = j;
437: col[2].i = i + 1;
438: v[3] = -hxdhy * (dn + gn);
439: col[3].j = j + 1;
440: col[3].i = i;
441: PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
442: /* left-hand top boundary */
443: } else {
444: ts = x[j - 1][i];
445: as = 0.5 * (t0 + ts);
446: bs = PetscPowScalar(as, bm1);
447: /* ds = bs * as; */
448: ds = PetscPowScalar(as, beta);
449: gs = coef * bs * (t0 - ts);
451: v[0] = -hxdhy * (ds - gs);
452: col[0].j = j - 1;
453: col[0].i = i;
454: v[1] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
455: col[1].j = row.j = j;
456: col[1].i = row.i = i;
457: v[2] = -hydhx * (de + ge);
458: col[2].j = j;
459: col[2].i = i + 1;
460: PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
461: }
463: } else if (i == mx - 1) {
464: /* right-hand boundary */
465: tw = x[j][i - 1];
466: aw = 0.5 * (t0 + tw);
467: bw = PetscPowScalar(aw, bm1);
468: /* dw = bw * aw */
469: dw = PetscPowScalar(aw, beta);
470: gw = coef * bw * (t0 - tw);
472: te = tright;
473: ae = 0.5 * (t0 + te);
474: be = PetscPowScalar(ae, bm1);
475: /* de = be * ae; */
476: de = PetscPowScalar(ae, beta);
477: ge = coef * be * (te - t0);
479: /* right-hand bottom boundary */
480: if (j == 0) {
481: tn = x[j + 1][i];
482: an = 0.5 * (t0 + tn);
483: bn = PetscPowScalar(an, bm1);
484: /* dn = bn * an; */
485: dn = PetscPowScalar(an, beta);
486: gn = coef * bn * (tn - t0);
488: v[0] = -hydhx * (dw - gw);
489: col[0].j = j;
490: col[0].i = i - 1;
491: v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
492: col[1].j = row.j = j;
493: col[1].i = row.i = i;
494: v[2] = -hxdhy * (dn + gn);
495: col[2].j = j + 1;
496: col[2].i = i;
497: PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
499: /* right-hand interior boundary */
500: } else if (j < my - 1) {
501: ts = x[j - 1][i];
502: as = 0.5 * (t0 + ts);
503: bs = PetscPowScalar(as, bm1);
504: /* ds = bs * as; */
505: ds = PetscPowScalar(as, beta);
506: gs = coef * bs * (t0 - ts);
508: tn = x[j + 1][i];
509: an = 0.5 * (t0 + tn);
510: bn = PetscPowScalar(an, bm1);
511: /* dn = bn * an; */
512: dn = PetscPowScalar(an, beta);
513: gn = coef * bn * (tn - t0);
515: v[0] = -hxdhy * (ds - gs);
516: col[0].j = j - 1;
517: col[0].i = i;
518: v[1] = -hydhx * (dw - gw);
519: col[1].j = j;
520: col[1].i = i - 1;
521: v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
522: col[2].j = row.j = j;
523: col[2].i = row.i = i;
524: v[3] = -hxdhy * (dn + gn);
525: col[3].j = j + 1;
526: col[3].i = i;
527: PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
528: /* right-hand top boundary */
529: } 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);
537: v[0] = -hxdhy * (ds - gs);
538: col[0].j = j - 1;
539: col[0].i = i;
540: v[1] = -hydhx * (dw - gw);
541: col[1].j = j;
542: col[1].i = i - 1;
543: v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
544: col[2].j = row.j = j;
545: col[2].i = row.i = i;
546: PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
547: }
549: /* bottom boundary,and i <> 0 or mx-1 */
550: } else if (j == 0) {
551: tw = x[j][i - 1];
552: aw = 0.5 * (t0 + tw);
553: bw = PetscPowScalar(aw, bm1);
554: /* dw = bw * aw */
555: dw = PetscPowScalar(aw, beta);
556: gw = coef * bw * (t0 - tw);
558: te = x[j][i + 1];
559: ae = 0.5 * (t0 + te);
560: be = PetscPowScalar(ae, bm1);
561: /* de = be * ae; */
562: de = PetscPowScalar(ae, beta);
563: ge = coef * be * (te - t0);
565: tn = x[j + 1][i];
566: an = 0.5 * (t0 + tn);
567: bn = PetscPowScalar(an, bm1);
568: /* dn = bn * an; */
569: dn = PetscPowScalar(an, beta);
570: gn = coef * bn * (tn - t0);
572: v[0] = -hydhx * (dw - gw);
573: col[0].j = j;
574: col[0].i = i - 1;
575: v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
576: col[1].j = row.j = j;
577: col[1].i = row.i = i;
578: v[2] = -hydhx * (de + ge);
579: col[2].j = j;
580: col[2].i = i + 1;
581: v[3] = -hxdhy * (dn + gn);
582: col[3].j = j + 1;
583: col[3].i = i;
584: PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
586: /* top boundary,and i <> 0 or mx-1 */
587: } else if (j == my - 1) {
588: tw = x[j][i - 1];
589: aw = 0.5 * (t0 + tw);
590: bw = PetscPowScalar(aw, bm1);
591: /* dw = bw * aw */
592: dw = PetscPowScalar(aw, beta);
593: gw = coef * bw * (t0 - tw);
595: te = x[j][i + 1];
596: ae = 0.5 * (t0 + te);
597: be = PetscPowScalar(ae, bm1);
598: /* de = be * ae; */
599: de = PetscPowScalar(ae, beta);
600: ge = coef * be * (te - t0);
602: ts = x[j - 1][i];
603: as = 0.5 * (t0 + ts);
604: bs = PetscPowScalar(as, bm1);
605: /* ds = bs * as; */
606: ds = PetscPowScalar(as, beta);
607: gs = coef * bs * (t0 - ts);
609: v[0] = -hxdhy * (ds - gs);
610: col[0].j = j - 1;
611: col[0].i = i;
612: v[1] = -hydhx * (dw - gw);
613: col[1].j = j;
614: col[1].i = i - 1;
615: v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
616: col[2].j = row.j = j;
617: col[2].i = row.i = i;
618: v[3] = -hydhx * (de + ge);
619: col[3].j = j;
620: col[3].i = i + 1;
621: PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
622: }
623: }
624: }
625: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
626: PetscCall(DMDAVecRestoreArray(da, localX, &x));
627: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
628: PetscCall(DMRestoreLocalVector(da, &localX));
629: if (jac != B) {
630: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
631: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
632: }
634: PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /*TEST
640: test:
641: suffix: 1
642: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
643: requires: !single
645: test:
646: suffix: 2
647: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc
648: requires: !single
650: test:
651: suffix: 3
652: args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontr -snes_tr_fallback_type dogleg
653: requires: !single
655: TEST*/