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";
12
13 /*
14
15 This example models the partial differential equation
16
17 - Div(alpha* T^beta (GRAD T)) = 0.
18
19 where beta = 2.5 and alpha = 1.0
20
21 BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
22
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.
25
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.
29
30 This code was contributed by David Keyes
31
32 */
33
34 #include <petscsnes.h>
35 #include <petscdm.h>
36 #include <petscdmda.h>
37
38 /* User-defined application context */
39
40 typedef struct {
41 PetscReal tleft, tright; /* Dirichlet boundary conditions */
42 PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
43 } AppCtx;
44
45 #define POWFLOP 5 /* assume a pow() takes five flops */
46
47 extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
48 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
49 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
50
main(int argc,char ** argv)51 int main(int argc, char **argv)
52 {
53 SNES snes;
54 AppCtx user;
55 PetscInt its, lits;
56 PetscReal litspit;
57 DM da;
58
59 PetscFunctionBeginUser;
60 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
61
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;
71
72 /*
73 Create the multilevel DM data structure
74 */
75 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
76
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));
85
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));
93
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));
101
102 PetscCall(DMDestroy(&da));
103 PetscCall(SNESDestroy(&snes));
104 PetscCall(PetscFinalize());
105 return 0;
106 }
107 /* -------------------- Form initial approximation ----------------- */
FormInitialGuess(SNES snes,Vec X,PetscCtx ctx)108 PetscErrorCode FormInitialGuess(SNES snes, Vec X, PetscCtx ctx)
109 {
110 AppCtx *user;
111 PetscInt i, j, xs, ys, xm, ym;
112 PetscReal tleft;
113 PetscScalar **x;
114 DM da;
115
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));
123
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) --------------------- */
FormFunction(SNES snes,Vec X,Vec F,void * ptr)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;
143
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;
155
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));
162
163 /* Evaluate function */
164 for (j = ys; j < ys + ym; j++) {
165 for (i = xs; i < xs + xm; i++) {
166 t0 = x[j][i];
167
168 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
169 /* general interior volume */
170
171 tw = x[j][i - 1];
172 aw = 0.5 * (t0 + tw);
173 dw = PetscPowScalar(aw, beta);
174 fw = dw * (t0 - tw);
175
176 te = x[j][i + 1];
177 ae = 0.5 * (t0 + te);
178 de = PetscPowScalar(ae, beta);
179 fe = de * (te - t0);
180
181 ts = x[j - 1][i];
182 as = 0.5 * (t0 + ts);
183 ds = PetscPowScalar(as, beta);
184 fs = ds * (t0 - ts);
185
186 tn = x[j + 1][i];
187 an = 0.5 * (t0 + tn);
188 dn = PetscPowScalar(an, beta);
189 fn = dn * (tn - t0);
190
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);
197
198 te = x[j][i + 1];
199 ae = 0.5 * (t0 + te);
200 de = PetscPowScalar(ae, beta);
201 fe = de * (te - t0);
202
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;
209
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;
216
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);
223
224 te = tright;
225 ae = 0.5 * (t0 + te);
226 de = PetscPowScalar(ae, beta);
227 fe = de * (te - t0);
228
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;
235
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;
242
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);
249
250 te = x[j][i + 1];
251 ae = 0.5 * (t0 + te);
252 de = PetscPowScalar(ae, beta);
253 fe = de * (te - t0);
254
255 fs = zero;
256
257 tn = x[j + 1][i];
258 an = 0.5 * (t0 + tn);
259 dn = PetscPowScalar(an, beta);
260 fn = dn * (tn - t0);
261
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);
268
269 te = x[j][i + 1];
270 ae = 0.5 * (t0 + te);
271 de = PetscPowScalar(ae, beta);
272 fe = de * (te - t0);
273
274 ts = x[j - 1][i];
275 as = 0.5 * (t0 + ts);
276 ds = PetscPowScalar(as, beta);
277 fs = ds * (t0 - ts);
278
279 fn = zero;
280 }
281
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) --------------------- */
FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void * ptr)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;
303
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;
317
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));
323
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];
328
329 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
330 /* general interior volume */
331
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);
338
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);
345
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);
352
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);
359
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));
376
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);
385
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);
392
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);
401
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));
412
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);
421
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);
428
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);
450
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 }
462
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);
471
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);
478
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);
487
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));
498
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);
507
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);
514
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);
536
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 }
548
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);
557
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);
564
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);
571
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));
585
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);
594
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);
601
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);
608
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 }
633
634 PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym));
635 PetscFunctionReturn(PETSC_SUCCESS);
636 }
637
638 /*TEST
639
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
644
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
649
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
654
655 TEST*/
656