1 static char help[] = "Nonlinear Radiative Transport PDE with multigrid in 3d.\n\
2 Uses 3-dimensional distributed arrays.\n\
3 A 3-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 = dT/dn_up = dT/dn_down = 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 7-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 Nickolas Jovanovic based on ex18.c
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 PetscBool converged; /* For user convergence test, whether to mark as converged */
44 } AppCtx;
45
46 #define POWFLOP 5 /* assume a pow() takes five flops */
47
48 extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
49 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
50 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
51 extern PetscErrorCode TestConvergence(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);
52
main(int argc,char ** argv)53 int main(int argc, char **argv)
54 {
55 SNES snes;
56 AppCtx user;
57 PetscInt its, lits;
58 PetscReal litspit = 0; /* avoid uninitialized warning */
59 DM da;
60 PetscBool use_convergence_test = PETSC_FALSE;
61
62 PetscFunctionBeginUser;
63 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
64 /* set problem parameters */
65 user.tleft = 1.0;
66 user.tright = 0.1;
67 user.beta = 2.5;
68 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL));
69 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL));
70 PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL));
71 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL));
72 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mark_converged", &user.converged, NULL));
73 user.bm1 = user.beta - 1.0;
74 user.coef = user.beta / 2.0;
75
76 /*
77 Set the DMDA (grid structure) for the grids.
78 */
79 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, 5, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, 0, &da));
80 PetscCall(DMSetFromOptions(da));
81 PetscCall(DMSetUp(da));
82 PetscCall(DMSetApplicationContext(da, &user));
83
84 /*
85 Create the nonlinear solver
86 */
87 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
88 PetscCall(SNESSetDM(snes, da));
89 PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user));
90 PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user));
91 if (use_convergence_test) PetscCall(SNESSetConvergenceTest(snes, TestConvergence, &user, NULL));
92 PetscCall(SNESSetFromOptions(snes));
93 PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL));
94
95 PetscCall(SNESSolve(snes, NULL, NULL));
96 PetscCall(SNESGetIterationNumber(snes, &its));
97 PetscCall(SNESGetLinearSolveIterations(snes, &lits));
98 if (its) litspit = ((PetscReal)lits) / ((PetscReal)its);
99 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
100 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits));
101 if (its) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit));
102
103 PetscCall(SNESDestroy(&snes));
104 PetscCall(DMDestroy(&da));
105 PetscCall(PetscFinalize());
106 return 0;
107 }
108 /* -------------------- Form initial approximation ----------------- */
FormInitialGuess(SNES snes,Vec X,PetscCtx ctx)109 PetscErrorCode FormInitialGuess(SNES snes, Vec X, PetscCtx ctx)
110 {
111 AppCtx *user;
112 PetscInt i, j, k, xs, ys, xm, ym, zs, zm;
113 PetscScalar ***x;
114 DM da;
115
116 PetscFunctionBeginUser;
117 PetscCall(SNESGetDM(snes, &da));
118 PetscCall(DMGetApplicationContext(da, &user));
119 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
120 PetscCall(DMDAVecGetArray(da, X, &x));
121
122 /* Compute initial guess */
123 for (k = zs; k < zs + zm; k++) {
124 for (j = ys; j < ys + ym; j++) {
125 for (i = xs; i < xs + xm; i++) x[k][j][i] = user->tleft;
126 }
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, k, mx, my, mz, xs, ys, zs, xm, ym, zm;
136 PetscScalar zero = 0.0, one = 1.0;
137 PetscScalar hx, hy, hz, hxhydhz, hyhzdhx, hzhxdhy;
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, td, ad, dd, fd = 0.0, tu, au, du = 0.0, fu = 0.0;
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, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
148 hx = one / (PetscReal)(mx - 1);
149 hy = one / (PetscReal)(my - 1);
150 hz = one / (PetscReal)(mz - 1);
151 hxhydhz = hx * hy / hz;
152 hyhzdhx = hy * hz / hx;
153 hzhxdhy = hz * hx / hy;
154 tleft = user->tleft;
155 tright = user->tright;
156 beta = user->beta;
157
158 /* Get ghost points */
159 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
160 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
161 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
162 PetscCall(DMDAVecGetArray(da, localX, &x));
163 PetscCall(DMDAVecGetArray(da, F, &f));
164
165 /* Evaluate function */
166 for (k = zs; k < zs + zm; k++) {
167 for (j = ys; j < ys + ym; j++) {
168 for (i = xs; i < xs + xm; i++) {
169 t0 = x[k][j][i];
170
171 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) {
172 /* general interior volume */
173
174 tw = x[k][j][i - 1];
175 aw = 0.5 * (t0 + tw);
176 dw = PetscPowScalar(aw, beta);
177 fw = dw * (t0 - tw);
178
179 te = x[k][j][i + 1];
180 ae = 0.5 * (t0 + te);
181 de = PetscPowScalar(ae, beta);
182 fe = de * (te - t0);
183
184 ts = x[k][j - 1][i];
185 as = 0.5 * (t0 + ts);
186 ds = PetscPowScalar(as, beta);
187 fs = ds * (t0 - ts);
188
189 tn = x[k][j + 1][i];
190 an = 0.5 * (t0 + tn);
191 dn = PetscPowScalar(an, beta);
192 fn = dn * (tn - t0);
193
194 td = x[k - 1][j][i];
195 ad = 0.5 * (t0 + td);
196 dd = PetscPowScalar(ad, beta);
197 fd = dd * (t0 - td);
198
199 tu = x[k + 1][j][i];
200 au = 0.5 * (t0 + tu);
201 du = PetscPowScalar(au, beta);
202 fu = du * (tu - t0);
203
204 } else if (i == 0) {
205 /* left-hand (west) boundary */
206 tw = tleft;
207 aw = 0.5 * (t0 + tw);
208 dw = PetscPowScalar(aw, beta);
209 fw = dw * (t0 - tw);
210
211 te = x[k][j][i + 1];
212 ae = 0.5 * (t0 + te);
213 de = PetscPowScalar(ae, beta);
214 fe = de * (te - t0);
215
216 if (j > 0) {
217 ts = x[k][j - 1][i];
218 as = 0.5 * (t0 + ts);
219 ds = PetscPowScalar(as, beta);
220 fs = ds * (t0 - ts);
221 } else {
222 fs = zero;
223 }
224
225 if (j < my - 1) {
226 tn = x[k][j + 1][i];
227 an = 0.5 * (t0 + tn);
228 dn = PetscPowScalar(an, beta);
229 fn = dn * (tn - t0);
230 } else {
231 fn = zero;
232 }
233
234 if (k > 0) {
235 td = x[k - 1][j][i];
236 ad = 0.5 * (t0 + td);
237 dd = PetscPowScalar(ad, beta);
238 fd = dd * (t0 - td);
239 } else {
240 fd = zero;
241 }
242
243 if (k < mz - 1) {
244 tu = x[k + 1][j][i];
245 au = 0.5 * (t0 + tu);
246 du = PetscPowScalar(au, beta);
247 fu = du * (tu - t0);
248 } else {
249 fu = zero;
250 }
251
252 } else if (i == mx - 1) {
253 /* right-hand (east) boundary */
254 tw = x[k][j][i - 1];
255 aw = 0.5 * (t0 + tw);
256 dw = PetscPowScalar(aw, beta);
257 fw = dw * (t0 - tw);
258
259 te = tright;
260 ae = 0.5 * (t0 + te);
261 de = PetscPowScalar(ae, beta);
262 fe = de * (te - t0);
263
264 if (j > 0) {
265 ts = x[k][j - 1][i];
266 as = 0.5 * (t0 + ts);
267 ds = PetscPowScalar(as, beta);
268 fs = ds * (t0 - ts);
269 } else {
270 fs = zero;
271 }
272
273 if (j < my - 1) {
274 tn = x[k][j + 1][i];
275 an = 0.5 * (t0 + tn);
276 dn = PetscPowScalar(an, beta);
277 fn = dn * (tn - t0);
278 } else {
279 fn = zero;
280 }
281
282 if (k > 0) {
283 td = x[k - 1][j][i];
284 ad = 0.5 * (t0 + td);
285 dd = PetscPowScalar(ad, beta);
286 fd = dd * (t0 - td);
287 } else {
288 fd = zero;
289 }
290
291 if (k < mz - 1) {
292 tu = x[k + 1][j][i];
293 au = 0.5 * (t0 + tu);
294 du = PetscPowScalar(au, beta);
295 fu = du * (tu - t0);
296 } else {
297 fu = zero;
298 }
299
300 } else if (j == 0) {
301 /* bottom (south) boundary, and i <> 0 or mx-1 */
302 tw = x[k][j][i - 1];
303 aw = 0.5 * (t0 + tw);
304 dw = PetscPowScalar(aw, beta);
305 fw = dw * (t0 - tw);
306
307 te = x[k][j][i + 1];
308 ae = 0.5 * (t0 + te);
309 de = PetscPowScalar(ae, beta);
310 fe = de * (te - t0);
311
312 fs = zero;
313
314 tn = x[k][j + 1][i];
315 an = 0.5 * (t0 + tn);
316 dn = PetscPowScalar(an, beta);
317 fn = dn * (tn - t0);
318
319 if (k > 0) {
320 td = x[k - 1][j][i];
321 ad = 0.5 * (t0 + td);
322 dd = PetscPowScalar(ad, beta);
323 fd = dd * (t0 - td);
324 } else {
325 fd = zero;
326 }
327
328 if (k < mz - 1) {
329 tu = x[k + 1][j][i];
330 au = 0.5 * (t0 + tu);
331 du = PetscPowScalar(au, beta);
332 fu = du * (tu - t0);
333 } else {
334 fu = zero;
335 }
336
337 } else if (j == my - 1) {
338 /* top (north) boundary, and i <> 0 or mx-1 */
339 tw = x[k][j][i - 1];
340 aw = 0.5 * (t0 + tw);
341 dw = PetscPowScalar(aw, beta);
342 fw = dw * (t0 - tw);
343
344 te = x[k][j][i + 1];
345 ae = 0.5 * (t0 + te);
346 de = PetscPowScalar(ae, beta);
347 fe = de * (te - t0);
348
349 ts = x[k][j - 1][i];
350 as = 0.5 * (t0 + ts);
351 ds = PetscPowScalar(as, beta);
352 fs = ds * (t0 - ts);
353
354 fn = zero;
355
356 if (k > 0) {
357 td = x[k - 1][j][i];
358 ad = 0.5 * (t0 + td);
359 dd = PetscPowScalar(ad, beta);
360 fd = dd * (t0 - td);
361 } else {
362 fd = zero;
363 }
364
365 if (k < mz - 1) {
366 tu = x[k + 1][j][i];
367 au = 0.5 * (t0 + tu);
368 du = PetscPowScalar(au, beta);
369 fu = du * (tu - t0);
370 } else {
371 fu = zero;
372 }
373
374 } else if (k == 0) {
375 /* down boundary (interior only) */
376 tw = x[k][j][i - 1];
377 aw = 0.5 * (t0 + tw);
378 dw = PetscPowScalar(aw, beta);
379 fw = dw * (t0 - tw);
380
381 te = x[k][j][i + 1];
382 ae = 0.5 * (t0 + te);
383 de = PetscPowScalar(ae, beta);
384 fe = de * (te - t0);
385
386 ts = x[k][j - 1][i];
387 as = 0.5 * (t0 + ts);
388 ds = PetscPowScalar(as, beta);
389 fs = ds * (t0 - ts);
390
391 tn = x[k][j + 1][i];
392 an = 0.5 * (t0 + tn);
393 dn = PetscPowScalar(an, beta);
394 fn = dn * (tn - t0);
395
396 fd = zero;
397
398 tu = x[k + 1][j][i];
399 au = 0.5 * (t0 + tu);
400 du = PetscPowScalar(au, beta);
401 fu = du * (tu - t0);
402
403 } else if (k == mz - 1) {
404 /* up boundary (interior only) */
405 tw = x[k][j][i - 1];
406 aw = 0.5 * (t0 + tw);
407 dw = PetscPowScalar(aw, beta);
408 fw = dw * (t0 - tw);
409
410 te = x[k][j][i + 1];
411 ae = 0.5 * (t0 + te);
412 de = PetscPowScalar(ae, beta);
413 fe = de * (te - t0);
414
415 ts = x[k][j - 1][i];
416 as = 0.5 * (t0 + ts);
417 ds = PetscPowScalar(as, beta);
418 fs = ds * (t0 - ts);
419
420 tn = x[k][j + 1][i];
421 an = 0.5 * (t0 + tn);
422 dn = PetscPowScalar(an, beta);
423 fn = dn * (tn - t0);
424
425 td = x[k - 1][j][i];
426 ad = 0.5 * (t0 + td);
427 dd = PetscPowScalar(ad, beta);
428 fd = dd * (t0 - td);
429
430 fu = zero;
431 }
432
433 f[k][j][i] = -hyhzdhx * (fe - fw) - hzhxdhy * (fn - fs) - hxhydhz * (fu - fd);
434 }
435 }
436 }
437 PetscCall(DMDAVecRestoreArray(da, localX, &x));
438 PetscCall(DMDAVecRestoreArray(da, F, &f));
439 PetscCall(DMRestoreLocalVector(da, &localX));
440 PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm));
441 PetscFunctionReturn(PETSC_SUCCESS);
442 }
443 /* -------------------- Evaluate Jacobian F(x) --------------------- */
FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void * ptr)444 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
445 {
446 AppCtx *user = (AppCtx *)ptr;
447 PetscInt i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm;
448 PetscScalar one = 1.0;
449 PetscScalar hx, hy, hz, hxhydhz, hyhzdhx, hzhxdhy;
450 PetscScalar t0, tn, ts, te, tw, an, as, ae, aw, dn, ds, de, dw;
451 PetscScalar tleft, tright, beta, td, ad, dd, tu, au, du, v[7], bm1, coef;
452 PetscScalar ***x, bn, bs, be, bw, bu, bd, gn, gs, ge, gw, gu, gd;
453 Vec localX;
454 MatStencil c[7], row;
455 DM da;
456
457 PetscFunctionBeginUser;
458 PetscCall(SNESGetDM(snes, &da));
459 PetscCall(DMGetLocalVector(da, &localX));
460 PetscCall(DMDAGetInfo(da, NULL, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
461 hx = one / (PetscReal)(mx - 1);
462 hy = one / (PetscReal)(my - 1);
463 hz = one / (PetscReal)(mz - 1);
464 hxhydhz = hx * hy / hz;
465 hyhzdhx = hy * hz / hx;
466 hzhxdhy = hz * hx / hy;
467 tleft = user->tleft;
468 tright = user->tright;
469 beta = user->beta;
470 bm1 = user->bm1;
471 coef = user->coef;
472
473 /* Get ghost points */
474 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
475 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
476 PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));
477 PetscCall(DMDAVecGetArray(da, localX, &x));
478
479 /* Evaluate Jacobian of function */
480 for (k = zs; k < zs + zm; k++) {
481 for (j = ys; j < ys + ym; j++) {
482 for (i = xs; i < xs + xm; i++) {
483 t0 = x[k][j][i];
484 row.k = k;
485 row.j = j;
486 row.i = i;
487 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) {
488 /* general interior volume */
489
490 tw = x[k][j][i - 1];
491 aw = 0.5 * (t0 + tw);
492 bw = PetscPowScalar(aw, bm1);
493 /* dw = bw * aw */
494 dw = PetscPowScalar(aw, beta);
495 gw = coef * bw * (t0 - tw);
496
497 te = x[k][j][i + 1];
498 ae = 0.5 * (t0 + te);
499 be = PetscPowScalar(ae, bm1);
500 /* de = be * ae; */
501 de = PetscPowScalar(ae, beta);
502 ge = coef * be * (te - t0);
503
504 ts = x[k][j - 1][i];
505 as = 0.5 * (t0 + ts);
506 bs = PetscPowScalar(as, bm1);
507 /* ds = bs * as; */
508 ds = PetscPowScalar(as, beta);
509 gs = coef * bs * (t0 - ts);
510
511 tn = x[k][j + 1][i];
512 an = 0.5 * (t0 + tn);
513 bn = PetscPowScalar(an, bm1);
514 /* dn = bn * an; */
515 dn = PetscPowScalar(an, beta);
516 gn = coef * bn * (tn - t0);
517
518 td = x[k - 1][j][i];
519 ad = 0.5 * (t0 + td);
520 bd = PetscPowScalar(ad, bm1);
521 /* dd = bd * ad; */
522 dd = PetscPowScalar(ad, beta);
523 gd = coef * bd * (t0 - td);
524
525 tu = x[k + 1][j][i];
526 au = 0.5 * (t0 + tu);
527 bu = PetscPowScalar(au, bm1);
528 /* du = bu * au; */
529 du = PetscPowScalar(au, beta);
530 gu = coef * bu * (tu - t0);
531
532 c[0].k = k - 1;
533 c[0].j = j;
534 c[0].i = i;
535 v[0] = -hxhydhz * (dd - gd);
536 c[1].k = k;
537 c[1].j = j - 1;
538 c[1].i = i;
539 v[1] = -hzhxdhy * (ds - gs);
540 c[2].k = k;
541 c[2].j = j;
542 c[2].i = i - 1;
543 v[2] = -hyhzdhx * (dw - gw);
544 c[3].k = k;
545 c[3].j = j;
546 c[3].i = i;
547 v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
548 c[4].k = k;
549 c[4].j = j;
550 c[4].i = i + 1;
551 v[4] = -hyhzdhx * (de + ge);
552 c[5].k = k;
553 c[5].j = j + 1;
554 c[5].i = i;
555 v[5] = -hzhxdhy * (dn + gn);
556 c[6].k = k + 1;
557 c[6].j = j;
558 c[6].i = i;
559 v[6] = -hxhydhz * (du + gu);
560 PetscCall(MatSetValuesStencil(jac, 1, &row, 7, c, v, INSERT_VALUES));
561
562 } else if (i == 0) {
563 /* left-hand plane boundary */
564 tw = tleft;
565 aw = 0.5 * (t0 + tw);
566 bw = PetscPowScalar(aw, bm1);
567 /* dw = bw * aw */
568 dw = PetscPowScalar(aw, beta);
569 gw = coef * bw * (t0 - tw);
570
571 te = x[k][j][i + 1];
572 ae = 0.5 * (t0 + te);
573 be = PetscPowScalar(ae, bm1);
574 /* de = be * ae; */
575 de = PetscPowScalar(ae, beta);
576 ge = coef * be * (te - t0);
577
578 /* left-hand bottom edge */
579 if (j == 0) {
580 tn = x[k][j + 1][i];
581 an = 0.5 * (t0 + tn);
582 bn = PetscPowScalar(an, bm1);
583 /* dn = bn * an; */
584 dn = PetscPowScalar(an, beta);
585 gn = coef * bn * (tn - t0);
586
587 /* left-hand bottom down corner */
588 if (k == 0) {
589 tu = x[k + 1][j][i];
590 au = 0.5 * (t0 + tu);
591 bu = PetscPowScalar(au, bm1);
592 /* du = bu * au; */
593 du = PetscPowScalar(au, beta);
594 gu = coef * bu * (tu - t0);
595
596 c[0].k = k;
597 c[0].j = j;
598 c[0].i = i;
599 v[0] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
600 c[1].k = k;
601 c[1].j = j;
602 c[1].i = i + 1;
603 v[1] = -hyhzdhx * (de + ge);
604 c[2].k = k;
605 c[2].j = j + 1;
606 c[2].i = i;
607 v[2] = -hzhxdhy * (dn + gn);
608 c[3].k = k + 1;
609 c[3].j = j;
610 c[3].i = i;
611 v[3] = -hxhydhz * (du + gu);
612 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
613
614 /* left-hand bottom interior edge */
615 } else if (k < mz - 1) {
616 tu = x[k + 1][j][i];
617 au = 0.5 * (t0 + tu);
618 bu = PetscPowScalar(au, bm1);
619 /* du = bu * au; */
620 du = PetscPowScalar(au, beta);
621 gu = coef * bu * (tu - t0);
622
623 td = x[k - 1][j][i];
624 ad = 0.5 * (t0 + td);
625 bd = PetscPowScalar(ad, bm1);
626 /* dd = bd * ad; */
627 dd = PetscPowScalar(ad, beta);
628 gd = coef * bd * (td - t0);
629
630 c[0].k = k - 1;
631 c[0].j = j;
632 c[0].i = i;
633 v[0] = -hxhydhz * (dd - gd);
634 c[1].k = k;
635 c[1].j = j;
636 c[1].i = i;
637 v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
638 c[2].k = k;
639 c[2].j = j;
640 c[2].i = i + 1;
641 v[2] = -hyhzdhx * (de + ge);
642 c[3].k = k;
643 c[3].j = j + 1;
644 c[3].i = i;
645 v[3] = -hzhxdhy * (dn + gn);
646 c[4].k = k + 1;
647 c[4].j = j;
648 c[4].i = i;
649 v[4] = -hxhydhz * (du + gu);
650 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
651
652 /* left-hand bottom up corner */
653 } else {
654 td = x[k - 1][j][i];
655 ad = 0.5 * (t0 + td);
656 bd = PetscPowScalar(ad, bm1);
657 /* dd = bd * ad; */
658 dd = PetscPowScalar(ad, beta);
659 gd = coef * bd * (td - t0);
660
661 c[0].k = k - 1;
662 c[0].j = j;
663 c[0].i = i;
664 v[0] = -hxhydhz * (dd - gd);
665 c[1].k = k;
666 c[1].j = j;
667 c[1].i = i;
668 v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
669 c[2].k = k;
670 c[2].j = j;
671 c[2].i = i + 1;
672 v[2] = -hyhzdhx * (de + ge);
673 c[3].k = k;
674 c[3].j = j + 1;
675 c[3].i = i;
676 v[3] = -hzhxdhy * (dn + gn);
677 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
678 }
679
680 /* left-hand top edge */
681 } else if (j == my - 1) {
682 ts = x[k][j - 1][i];
683 as = 0.5 * (t0 + ts);
684 bs = PetscPowScalar(as, bm1);
685 /* ds = bs * as; */
686 ds = PetscPowScalar(as, beta);
687 gs = coef * bs * (ts - t0);
688
689 /* left-hand top down corner */
690 if (k == 0) {
691 tu = x[k + 1][j][i];
692 au = 0.5 * (t0 + tu);
693 bu = PetscPowScalar(au, bm1);
694 /* du = bu * au; */
695 du = PetscPowScalar(au, beta);
696 gu = coef * bu * (tu - t0);
697
698 c[0].k = k;
699 c[0].j = j - 1;
700 c[0].i = i;
701 v[0] = -hzhxdhy * (ds - gs);
702 c[1].k = k;
703 c[1].j = j;
704 c[1].i = i;
705 v[1] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
706 c[2].k = k;
707 c[2].j = j;
708 c[2].i = i + 1;
709 v[2] = -hyhzdhx * (de + ge);
710 c[3].k = k + 1;
711 c[3].j = j;
712 c[3].i = i;
713 v[3] = -hxhydhz * (du + gu);
714 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
715
716 /* left-hand top interior edge */
717 } else if (k < mz - 1) {
718 tu = x[k + 1][j][i];
719 au = 0.5 * (t0 + tu);
720 bu = PetscPowScalar(au, bm1);
721 /* du = bu * au; */
722 du = PetscPowScalar(au, beta);
723 gu = coef * bu * (tu - t0);
724
725 td = x[k - 1][j][i];
726 ad = 0.5 * (t0 + td);
727 bd = PetscPowScalar(ad, bm1);
728 /* dd = bd * ad; */
729 dd = PetscPowScalar(ad, beta);
730 gd = coef * bd * (td - t0);
731
732 c[0].k = k - 1;
733 c[0].j = j;
734 c[0].i = i;
735 v[0] = -hxhydhz * (dd - gd);
736 c[1].k = k;
737 c[1].j = j - 1;
738 c[1].i = i;
739 v[1] = -hzhxdhy * (ds - gs);
740 c[2].k = k;
741 c[2].j = j;
742 c[2].i = i;
743 v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
744 c[3].k = k;
745 c[3].j = j;
746 c[3].i = i + 1;
747 v[3] = -hyhzdhx * (de + ge);
748 c[4].k = k + 1;
749 c[4].j = j;
750 c[4].i = i;
751 v[4] = -hxhydhz * (du + gu);
752 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
753
754 /* left-hand top up corner */
755 } else {
756 td = x[k - 1][j][i];
757 ad = 0.5 * (t0 + td);
758 bd = PetscPowScalar(ad, bm1);
759 /* dd = bd * ad; */
760 dd = PetscPowScalar(ad, beta);
761 gd = coef * bd * (td - t0);
762
763 c[0].k = k - 1;
764 c[0].j = j;
765 c[0].i = i;
766 v[0] = -hxhydhz * (dd - gd);
767 c[1].k = k;
768 c[1].j = j - 1;
769 c[1].i = i;
770 v[1] = -hzhxdhy * (ds - gs);
771 c[2].k = k;
772 c[2].j = j;
773 c[2].i = i;
774 v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
775 c[3].k = k;
776 c[3].j = j;
777 c[3].i = i + 1;
778 v[3] = -hyhzdhx * (de + ge);
779 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
780 }
781
782 } else {
783 ts = x[k][j - 1][i];
784 as = 0.5 * (t0 + ts);
785 bs = PetscPowScalar(as, bm1);
786 /* ds = bs * as; */
787 ds = PetscPowScalar(as, beta);
788 gs = coef * bs * (t0 - ts);
789
790 tn = x[k][j + 1][i];
791 an = 0.5 * (t0 + tn);
792 bn = PetscPowScalar(an, bm1);
793 /* dn = bn * an; */
794 dn = PetscPowScalar(an, beta);
795 gn = coef * bn * (tn - t0);
796
797 /* left-hand down interior edge */
798 if (k == 0) {
799 tu = x[k + 1][j][i];
800 au = 0.5 * (t0 + tu);
801 bu = PetscPowScalar(au, bm1);
802 /* du = bu * au; */
803 du = PetscPowScalar(au, beta);
804 gu = coef * bu * (tu - t0);
805
806 c[0].k = k;
807 c[0].j = j - 1;
808 c[0].i = i;
809 v[0] = -hzhxdhy * (ds - gs);
810 c[1].k = k;
811 c[1].j = j;
812 c[1].i = i;
813 v[1] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
814 c[2].k = k;
815 c[2].j = j;
816 c[2].i = i + 1;
817 v[2] = -hyhzdhx * (de + ge);
818 c[3].k = k;
819 c[3].j = j + 1;
820 c[3].i = i;
821 v[3] = -hzhxdhy * (dn + gn);
822 c[4].k = k + 1;
823 c[4].j = j;
824 c[4].i = i;
825 v[4] = -hxhydhz * (du + gu);
826 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
827
828 } else if (k == mz - 1) { /* left-hand up interior edge */
829
830 td = x[k - 1][j][i];
831 ad = 0.5 * (t0 + td);
832 bd = PetscPowScalar(ad, bm1);
833 /* dd = bd * ad; */
834 dd = PetscPowScalar(ad, beta);
835 gd = coef * bd * (t0 - td);
836
837 c[0].k = k - 1;
838 c[0].j = j;
839 c[0].i = i;
840 v[0] = -hxhydhz * (dd - gd);
841 c[1].k = k;
842 c[1].j = j - 1;
843 c[1].i = i;
844 v[1] = -hzhxdhy * (ds - gs);
845 c[2].k = k;
846 c[2].j = j;
847 c[2].i = i;
848 v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
849 c[3].k = k;
850 c[3].j = j;
851 c[3].i = i + 1;
852 v[3] = -hyhzdhx * (de + ge);
853 c[4].k = k;
854 c[4].j = j + 1;
855 c[4].i = i;
856 v[4] = -hzhxdhy * (dn + gn);
857 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
858 } else { /* left-hand interior plane */
859
860 td = x[k - 1][j][i];
861 ad = 0.5 * (t0 + td);
862 bd = PetscPowScalar(ad, bm1);
863 /* dd = bd * ad; */
864 dd = PetscPowScalar(ad, beta);
865 gd = coef * bd * (t0 - td);
866
867 tu = x[k + 1][j][i];
868 au = 0.5 * (t0 + tu);
869 bu = PetscPowScalar(au, bm1);
870 /* du = bu * au; */
871 du = PetscPowScalar(au, beta);
872 gu = coef * bu * (tu - t0);
873
874 c[0].k = k - 1;
875 c[0].j = j;
876 c[0].i = i;
877 v[0] = -hxhydhz * (dd - gd);
878 c[1].k = k;
879 c[1].j = j - 1;
880 c[1].i = i;
881 v[1] = -hzhxdhy * (ds - gs);
882 c[2].k = k;
883 c[2].j = j;
884 c[2].i = i;
885 v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
886 c[3].k = k;
887 c[3].j = j;
888 c[3].i = i + 1;
889 v[3] = -hyhzdhx * (de + ge);
890 c[4].k = k;
891 c[4].j = j + 1;
892 c[4].i = i;
893 v[4] = -hzhxdhy * (dn + gn);
894 c[5].k = k + 1;
895 c[5].j = j;
896 c[5].i = i;
897 v[5] = -hxhydhz * (du + gu);
898 PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
899 }
900 }
901
902 } else if (i == mx - 1) {
903 /* right-hand plane boundary */
904 tw = x[k][j][i - 1];
905 aw = 0.5 * (t0 + tw);
906 bw = PetscPowScalar(aw, bm1);
907 /* dw = bw * aw */
908 dw = PetscPowScalar(aw, beta);
909 gw = coef * bw * (t0 - tw);
910
911 te = tright;
912 ae = 0.5 * (t0 + te);
913 be = PetscPowScalar(ae, bm1);
914 /* de = be * ae; */
915 de = PetscPowScalar(ae, beta);
916 ge = coef * be * (te - t0);
917
918 /* right-hand bottom edge */
919 if (j == 0) {
920 tn = x[k][j + 1][i];
921 an = 0.5 * (t0 + tn);
922 bn = PetscPowScalar(an, bm1);
923 /* dn = bn * an; */
924 dn = PetscPowScalar(an, beta);
925 gn = coef * bn * (tn - t0);
926
927 /* right-hand bottom down corner */
928 if (k == 0) {
929 tu = x[k + 1][j][i];
930 au = 0.5 * (t0 + tu);
931 bu = PetscPowScalar(au, bm1);
932 /* du = bu * au; */
933 du = PetscPowScalar(au, beta);
934 gu = coef * bu * (tu - t0);
935
936 c[0].k = k;
937 c[0].j = j;
938 c[0].i = i - 1;
939 v[0] = -hyhzdhx * (dw - gw);
940 c[1].k = k;
941 c[1].j = j;
942 c[1].i = i;
943 v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
944 c[2].k = k;
945 c[2].j = j + 1;
946 c[2].i = i;
947 v[2] = -hzhxdhy * (dn + gn);
948 c[3].k = k + 1;
949 c[3].j = j;
950 c[3].i = i;
951 v[3] = -hxhydhz * (du + gu);
952 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
953
954 /* right-hand bottom interior edge */
955 } else if (k < mz - 1) {
956 tu = x[k + 1][j][i];
957 au = 0.5 * (t0 + tu);
958 bu = PetscPowScalar(au, bm1);
959 /* du = bu * au; */
960 du = PetscPowScalar(au, beta);
961 gu = coef * bu * (tu - t0);
962
963 td = x[k - 1][j][i];
964 ad = 0.5 * (t0 + td);
965 bd = PetscPowScalar(ad, bm1);
966 /* dd = bd * ad; */
967 dd = PetscPowScalar(ad, beta);
968 gd = coef * bd * (td - t0);
969
970 c[0].k = k - 1;
971 c[0].j = j;
972 c[0].i = i;
973 v[0] = -hxhydhz * (dd - gd);
974 c[1].k = k;
975 c[1].j = j;
976 c[1].i = i - 1;
977 v[1] = -hyhzdhx * (dw - gw);
978 c[2].k = k;
979 c[2].j = j;
980 c[2].i = i;
981 v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
982 c[3].k = k;
983 c[3].j = j + 1;
984 c[3].i = i;
985 v[3] = -hzhxdhy * (dn + gn);
986 c[4].k = k + 1;
987 c[4].j = j;
988 c[4].i = i;
989 v[4] = -hxhydhz * (du + gu);
990 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
991
992 /* right-hand bottom up corner */
993 } else {
994 td = x[k - 1][j][i];
995 ad = 0.5 * (t0 + td);
996 bd = PetscPowScalar(ad, bm1);
997 /* dd = bd * ad; */
998 dd = PetscPowScalar(ad, beta);
999 gd = coef * bd * (td - t0);
1000
1001 c[0].k = k - 1;
1002 c[0].j = j;
1003 c[0].i = i;
1004 v[0] = -hxhydhz * (dd - gd);
1005 c[1].k = k;
1006 c[1].j = j;
1007 c[1].i = i - 1;
1008 v[1] = -hyhzdhx * (dw - gw);
1009 c[2].k = k;
1010 c[2].j = j;
1011 c[2].i = i;
1012 v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1013 c[3].k = k;
1014 c[3].j = j + 1;
1015 c[3].i = i;
1016 v[3] = -hzhxdhy * (dn + gn);
1017 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1018 }
1019
1020 /* right-hand top edge */
1021 } else if (j == my - 1) {
1022 ts = x[k][j - 1][i];
1023 as = 0.5 * (t0 + ts);
1024 bs = PetscPowScalar(as, bm1);
1025 /* ds = bs * as; */
1026 ds = PetscPowScalar(as, beta);
1027 gs = coef * bs * (ts - t0);
1028
1029 /* right-hand top down corner */
1030 if (k == 0) {
1031 tu = x[k + 1][j][i];
1032 au = 0.5 * (t0 + tu);
1033 bu = PetscPowScalar(au, bm1);
1034 /* du = bu * au; */
1035 du = PetscPowScalar(au, beta);
1036 gu = coef * bu * (tu - t0);
1037
1038 c[0].k = k;
1039 c[0].j = j - 1;
1040 c[0].i = i;
1041 v[0] = -hzhxdhy * (ds - gs);
1042 c[1].k = k;
1043 c[1].j = j;
1044 c[1].i = i - 1;
1045 v[1] = -hyhzdhx * (dw - gw);
1046 c[2].k = k;
1047 c[2].j = j;
1048 c[2].i = i;
1049 v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1050 c[3].k = k + 1;
1051 c[3].j = j;
1052 c[3].i = i;
1053 v[3] = -hxhydhz * (du + gu);
1054 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1055
1056 /* right-hand top interior edge */
1057 } else if (k < mz - 1) {
1058 tu = x[k + 1][j][i];
1059 au = 0.5 * (t0 + tu);
1060 bu = PetscPowScalar(au, bm1);
1061 /* du = bu * au; */
1062 du = PetscPowScalar(au, beta);
1063 gu = coef * bu * (tu - t0);
1064
1065 td = x[k - 1][j][i];
1066 ad = 0.5 * (t0 + td);
1067 bd = PetscPowScalar(ad, bm1);
1068 /* dd = bd * ad; */
1069 dd = PetscPowScalar(ad, beta);
1070 gd = coef * bd * (td - t0);
1071
1072 c[0].k = k - 1;
1073 c[0].j = j;
1074 c[0].i = i;
1075 v[0] = -hxhydhz * (dd - gd);
1076 c[1].k = k;
1077 c[1].j = j - 1;
1078 c[1].i = i;
1079 v[1] = -hzhxdhy * (ds - gs);
1080 c[2].k = k;
1081 c[2].j = j;
1082 c[2].i = i - 1;
1083 v[2] = -hyhzdhx * (dw - gw);
1084 c[3].k = k;
1085 c[3].j = j;
1086 c[3].i = i;
1087 v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1088 c[4].k = k + 1;
1089 c[4].j = j;
1090 c[4].i = i;
1091 v[4] = -hxhydhz * (du + gu);
1092 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1093
1094 /* right-hand top up corner */
1095 } else {
1096 td = x[k - 1][j][i];
1097 ad = 0.5 * (t0 + td);
1098 bd = PetscPowScalar(ad, bm1);
1099 /* dd = bd * ad; */
1100 dd = PetscPowScalar(ad, beta);
1101 gd = coef * bd * (td - t0);
1102
1103 c[0].k = k - 1;
1104 c[0].j = j;
1105 c[0].i = i;
1106 v[0] = -hxhydhz * (dd - gd);
1107 c[1].k = k;
1108 c[1].j = j - 1;
1109 c[1].i = i;
1110 v[1] = -hzhxdhy * (ds - gs);
1111 c[2].k = k;
1112 c[2].j = j;
1113 c[2].i = i - 1;
1114 v[2] = -hyhzdhx * (dw - gw);
1115 c[3].k = k;
1116 c[3].j = j;
1117 c[3].i = i;
1118 v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1119 PetscCall(MatSetValuesStencil(jac, 1, &row, 4, c, v, INSERT_VALUES));
1120 }
1121
1122 } else {
1123 ts = x[k][j - 1][i];
1124 as = 0.5 * (t0 + ts);
1125 bs = PetscPowScalar(as, bm1);
1126 /* ds = bs * as; */
1127 ds = PetscPowScalar(as, beta);
1128 gs = coef * bs * (t0 - ts);
1129
1130 tn = x[k][j + 1][i];
1131 an = 0.5 * (t0 + tn);
1132 bn = PetscPowScalar(an, bm1);
1133 /* dn = bn * an; */
1134 dn = PetscPowScalar(an, beta);
1135 gn = coef * bn * (tn - t0);
1136
1137 /* right-hand down interior edge */
1138 if (k == 0) {
1139 tu = x[k + 1][j][i];
1140 au = 0.5 * (t0 + tu);
1141 bu = PetscPowScalar(au, bm1);
1142 /* du = bu * au; */
1143 du = PetscPowScalar(au, beta);
1144 gu = coef * bu * (tu - t0);
1145
1146 c[0].k = k;
1147 c[0].j = j - 1;
1148 c[0].i = i;
1149 v[0] = -hzhxdhy * (ds - gs);
1150 c[1].k = k;
1151 c[1].j = j;
1152 c[1].i = i - 1;
1153 v[1] = -hyhzdhx * (dw - gw);
1154 c[2].k = k;
1155 c[2].j = j;
1156 c[2].i = i;
1157 v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1158 c[3].k = k;
1159 c[3].j = j + 1;
1160 c[3].i = i;
1161 v[3] = -hzhxdhy * (dn + gn);
1162 c[4].k = k + 1;
1163 c[4].j = j;
1164 c[4].i = i;
1165 v[4] = -hxhydhz * (du + gu);
1166 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1167
1168 } else if (k == mz - 1) { /* right-hand up interior edge */
1169
1170 td = x[k - 1][j][i];
1171 ad = 0.5 * (t0 + td);
1172 bd = PetscPowScalar(ad, bm1);
1173 /* dd = bd * ad; */
1174 dd = PetscPowScalar(ad, beta);
1175 gd = coef * bd * (t0 - td);
1176
1177 c[0].k = k - 1;
1178 c[0].j = j;
1179 c[0].i = i;
1180 v[0] = -hxhydhz * (dd - gd);
1181 c[1].k = k;
1182 c[1].j = j - 1;
1183 c[1].i = i;
1184 v[1] = -hzhxdhy * (ds - gs);
1185 c[2].k = k;
1186 c[2].j = j;
1187 c[2].i = i - 1;
1188 v[2] = -hyhzdhx * (dw - gw);
1189 c[3].k = k;
1190 c[3].j = j;
1191 c[3].i = i;
1192 v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1193 c[4].k = k;
1194 c[4].j = j + 1;
1195 c[4].i = i;
1196 v[4] = -hzhxdhy * (dn + gn);
1197 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1198
1199 } else { /* right-hand interior plane */
1200
1201 td = x[k - 1][j][i];
1202 ad = 0.5 * (t0 + td);
1203 bd = PetscPowScalar(ad, bm1);
1204 /* dd = bd * ad; */
1205 dd = PetscPowScalar(ad, beta);
1206 gd = coef * bd * (t0 - td);
1207
1208 tu = x[k + 1][j][i];
1209 au = 0.5 * (t0 + tu);
1210 bu = PetscPowScalar(au, bm1);
1211 /* du = bu * au; */
1212 du = PetscPowScalar(au, beta);
1213 gu = coef * bu * (tu - t0);
1214
1215 c[0].k = k - 1;
1216 c[0].j = j;
1217 c[0].i = i;
1218 v[0] = -hxhydhz * (dd - gd);
1219 c[1].k = k;
1220 c[1].j = j - 1;
1221 c[1].i = i;
1222 v[1] = -hzhxdhy * (ds - gs);
1223 c[2].k = k;
1224 c[2].j = j;
1225 c[2].i = i - 1;
1226 v[2] = -hyhzdhx * (dw - gw);
1227 c[3].k = k;
1228 c[3].j = j;
1229 c[3].i = i;
1230 v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1231 c[4].k = k;
1232 c[4].j = j + 1;
1233 c[4].i = i;
1234 v[4] = -hzhxdhy * (dn + gn);
1235 c[5].k = k + 1;
1236 c[5].j = j;
1237 c[5].i = i;
1238 v[5] = -hxhydhz * (du + gu);
1239 PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1240 }
1241 }
1242
1243 } else if (j == 0) {
1244 tw = x[k][j][i - 1];
1245 aw = 0.5 * (t0 + tw);
1246 bw = PetscPowScalar(aw, bm1);
1247 /* dw = bw * aw */
1248 dw = PetscPowScalar(aw, beta);
1249 gw = coef * bw * (t0 - tw);
1250
1251 te = x[k][j][i + 1];
1252 ae = 0.5 * (t0 + te);
1253 be = PetscPowScalar(ae, bm1);
1254 /* de = be * ae; */
1255 de = PetscPowScalar(ae, beta);
1256 ge = coef * be * (te - t0);
1257
1258 tn = x[k][j + 1][i];
1259 an = 0.5 * (t0 + tn);
1260 bn = PetscPowScalar(an, bm1);
1261 /* dn = bn * an; */
1262 dn = PetscPowScalar(an, beta);
1263 gn = coef * bn * (tn - t0);
1264
1265 /* bottom down interior edge */
1266 if (k == 0) {
1267 tu = x[k + 1][j][i];
1268 au = 0.5 * (t0 + tu);
1269 bu = PetscPowScalar(au, bm1);
1270 /* du = bu * au; */
1271 du = PetscPowScalar(au, beta);
1272 gu = coef * bu * (tu - t0);
1273
1274 c[0].k = k;
1275 c[0].j = j;
1276 c[0].i = i - 1;
1277 v[0] = -hyhzdhx * (dw - gw);
1278 c[1].k = k;
1279 c[1].j = j;
1280 c[1].i = i;
1281 v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1282 c[2].k = k;
1283 c[2].j = j;
1284 c[2].i = i + 1;
1285 v[2] = -hyhzdhx * (de + ge);
1286 c[3].k = k;
1287 c[3].j = j + 1;
1288 c[3].i = i;
1289 v[3] = -hzhxdhy * (dn + gn);
1290 c[4].k = k + 1;
1291 c[4].j = j;
1292 c[4].i = i;
1293 v[4] = -hxhydhz * (du + gu);
1294 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1295
1296 } else if (k == mz - 1) { /* bottom up interior edge */
1297
1298 td = x[k - 1][j][i];
1299 ad = 0.5 * (t0 + td);
1300 bd = PetscPowScalar(ad, bm1);
1301 /* dd = bd * ad; */
1302 dd = PetscPowScalar(ad, beta);
1303 gd = coef * bd * (td - t0);
1304
1305 c[0].k = k - 1;
1306 c[0].j = j;
1307 c[0].i = i;
1308 v[0] = -hxhydhz * (dd - gd);
1309 c[1].k = k;
1310 c[1].j = j;
1311 c[1].i = i - 1;
1312 v[1] = -hyhzdhx * (dw - gw);
1313 c[2].k = k;
1314 c[2].j = j;
1315 c[2].i = i;
1316 v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1317 c[3].k = k;
1318 c[3].j = j;
1319 c[3].i = i + 1;
1320 v[3] = -hyhzdhx * (de + ge);
1321 c[4].k = k;
1322 c[4].j = j + 1;
1323 c[4].i = i;
1324 v[4] = -hzhxdhy * (dn + gn);
1325 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1326
1327 } else { /* bottom interior plane */
1328
1329 tu = x[k + 1][j][i];
1330 au = 0.5 * (t0 + tu);
1331 bu = PetscPowScalar(au, bm1);
1332 /* du = bu * au; */
1333 du = PetscPowScalar(au, beta);
1334 gu = coef * bu * (tu - t0);
1335
1336 td = x[k - 1][j][i];
1337 ad = 0.5 * (t0 + td);
1338 bd = PetscPowScalar(ad, bm1);
1339 /* dd = bd * ad; */
1340 dd = PetscPowScalar(ad, beta);
1341 gd = coef * bd * (td - t0);
1342
1343 c[0].k = k - 1;
1344 c[0].j = j;
1345 c[0].i = i;
1346 v[0] = -hxhydhz * (dd - gd);
1347 c[1].k = k;
1348 c[1].j = j;
1349 c[1].i = i - 1;
1350 v[1] = -hyhzdhx * (dw - gw);
1351 c[2].k = k;
1352 c[2].j = j;
1353 c[2].i = i;
1354 v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1355 c[3].k = k;
1356 c[3].j = j;
1357 c[3].i = i + 1;
1358 v[3] = -hyhzdhx * (de + ge);
1359 c[4].k = k;
1360 c[4].j = j + 1;
1361 c[4].i = i;
1362 v[4] = -hzhxdhy * (dn + gn);
1363 c[5].k = k + 1;
1364 c[5].j = j;
1365 c[5].i = i;
1366 v[5] = -hxhydhz * (du + gu);
1367 PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1368 }
1369
1370 } else if (j == my - 1) {
1371 tw = x[k][j][i - 1];
1372 aw = 0.5 * (t0 + tw);
1373 bw = PetscPowScalar(aw, bm1);
1374 /* dw = bw * aw */
1375 dw = PetscPowScalar(aw, beta);
1376 gw = coef * bw * (t0 - tw);
1377
1378 te = x[k][j][i + 1];
1379 ae = 0.5 * (t0 + te);
1380 be = PetscPowScalar(ae, bm1);
1381 /* de = be * ae; */
1382 de = PetscPowScalar(ae, beta);
1383 ge = coef * be * (te - t0);
1384
1385 ts = x[k][j - 1][i];
1386 as = 0.5 * (t0 + ts);
1387 bs = PetscPowScalar(as, bm1);
1388 /* ds = bs * as; */
1389 ds = PetscPowScalar(as, beta);
1390 gs = coef * bs * (t0 - ts);
1391
1392 /* top down interior edge */
1393 if (k == 0) {
1394 tu = x[k + 1][j][i];
1395 au = 0.5 * (t0 + tu);
1396 bu = PetscPowScalar(au, bm1);
1397 /* du = bu * au; */
1398 du = PetscPowScalar(au, beta);
1399 gu = coef * bu * (tu - t0);
1400
1401 c[0].k = k;
1402 c[0].j = j - 1;
1403 c[0].i = i;
1404 v[0] = -hzhxdhy * (ds - gs);
1405 c[1].k = k;
1406 c[1].j = j;
1407 c[1].i = i - 1;
1408 v[1] = -hyhzdhx * (dw - gw);
1409 c[2].k = k;
1410 c[2].j = j;
1411 c[2].i = i;
1412 v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1413 c[3].k = k;
1414 c[3].j = j;
1415 c[3].i = i + 1;
1416 v[3] = -hyhzdhx * (de + ge);
1417 c[4].k = k + 1;
1418 c[4].j = j;
1419 c[4].i = i;
1420 v[4] = -hxhydhz * (du + gu);
1421 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1422
1423 } else if (k == mz - 1) { /* top up interior edge */
1424
1425 td = x[k - 1][j][i];
1426 ad = 0.5 * (t0 + td);
1427 bd = PetscPowScalar(ad, bm1);
1428 /* dd = bd * ad; */
1429 dd = PetscPowScalar(ad, beta);
1430 gd = coef * bd * (td - t0);
1431
1432 c[0].k = k - 1;
1433 c[0].j = j;
1434 c[0].i = i;
1435 v[0] = -hxhydhz * (dd - gd);
1436 c[1].k = k;
1437 c[1].j = j - 1;
1438 c[1].i = i;
1439 v[1] = -hzhxdhy * (ds - gs);
1440 c[2].k = k;
1441 c[2].j = j;
1442 c[2].i = i - 1;
1443 v[2] = -hyhzdhx * (dw - gw);
1444 c[3].k = k;
1445 c[3].j = j;
1446 c[3].i = i;
1447 v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1448 c[4].k = k;
1449 c[4].j = j;
1450 c[4].i = i + 1;
1451 v[4] = -hyhzdhx * (de + ge);
1452 PetscCall(MatSetValuesStencil(jac, 1, &row, 5, c, v, INSERT_VALUES));
1453
1454 } else { /* top interior plane */
1455
1456 tu = x[k + 1][j][i];
1457 au = 0.5 * (t0 + tu);
1458 bu = PetscPowScalar(au, bm1);
1459 /* du = bu * au; */
1460 du = PetscPowScalar(au, beta);
1461 gu = coef * bu * (tu - t0);
1462
1463 td = x[k - 1][j][i];
1464 ad = 0.5 * (t0 + td);
1465 bd = PetscPowScalar(ad, bm1);
1466 /* dd = bd * ad; */
1467 dd = PetscPowScalar(ad, beta);
1468 gd = coef * bd * (td - t0);
1469
1470 c[0].k = k - 1;
1471 c[0].j = j;
1472 c[0].i = i;
1473 v[0] = -hxhydhz * (dd - gd);
1474 c[1].k = k;
1475 c[1].j = j - 1;
1476 c[1].i = i;
1477 v[1] = -hzhxdhy * (ds - gs);
1478 c[2].k = k;
1479 c[2].j = j;
1480 c[2].i = i - 1;
1481 v[2] = -hyhzdhx * (dw - gw);
1482 c[3].k = k;
1483 c[3].j = j;
1484 c[3].i = i;
1485 v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu);
1486 c[4].k = k;
1487 c[4].j = j;
1488 c[4].i = i + 1;
1489 v[4] = -hyhzdhx * (de + ge);
1490 c[5].k = k + 1;
1491 c[5].j = j;
1492 c[5].i = i;
1493 v[5] = -hxhydhz * (du + gu);
1494 PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1495 }
1496
1497 } else if (k == 0) {
1498 /* down interior plane */
1499
1500 tw = x[k][j][i - 1];
1501 aw = 0.5 * (t0 + tw);
1502 bw = PetscPowScalar(aw, bm1);
1503 /* dw = bw * aw */
1504 dw = PetscPowScalar(aw, beta);
1505 gw = coef * bw * (t0 - tw);
1506
1507 te = x[k][j][i + 1];
1508 ae = 0.5 * (t0 + te);
1509 be = PetscPowScalar(ae, bm1);
1510 /* de = be * ae; */
1511 de = PetscPowScalar(ae, beta);
1512 ge = coef * be * (te - t0);
1513
1514 ts = x[k][j - 1][i];
1515 as = 0.5 * (t0 + ts);
1516 bs = PetscPowScalar(as, bm1);
1517 /* ds = bs * as; */
1518 ds = PetscPowScalar(as, beta);
1519 gs = coef * bs * (t0 - ts);
1520
1521 tn = x[k][j + 1][i];
1522 an = 0.5 * (t0 + tn);
1523 bn = PetscPowScalar(an, bm1);
1524 /* dn = bn * an; */
1525 dn = PetscPowScalar(an, beta);
1526 gn = coef * bn * (tn - t0);
1527
1528 tu = x[k + 1][j][i];
1529 au = 0.5 * (t0 + tu);
1530 bu = PetscPowScalar(au, bm1);
1531 /* du = bu * au; */
1532 du = PetscPowScalar(au, beta);
1533 gu = coef * bu * (tu - t0);
1534
1535 c[0].k = k;
1536 c[0].j = j - 1;
1537 c[0].i = i;
1538 v[0] = -hzhxdhy * (ds - gs);
1539 c[1].k = k;
1540 c[1].j = j;
1541 c[1].i = i - 1;
1542 v[1] = -hyhzdhx * (dw - gw);
1543 c[2].k = k;
1544 c[2].j = j;
1545 c[2].i = i;
1546 v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu);
1547 c[3].k = k;
1548 c[3].j = j;
1549 c[3].i = i + 1;
1550 v[3] = -hyhzdhx * (de + ge);
1551 c[4].k = k;
1552 c[4].j = j + 1;
1553 c[4].i = i;
1554 v[4] = -hzhxdhy * (dn + gn);
1555 c[5].k = k + 1;
1556 c[5].j = j;
1557 c[5].i = i;
1558 v[5] = -hxhydhz * (du + gu);
1559 PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1560
1561 } else if (k == mz - 1) {
1562 /* up interior plane */
1563
1564 tw = x[k][j][i - 1];
1565 aw = 0.5 * (t0 + tw);
1566 bw = PetscPowScalar(aw, bm1);
1567 /* dw = bw * aw */
1568 dw = PetscPowScalar(aw, beta);
1569 gw = coef * bw * (t0 - tw);
1570
1571 te = x[k][j][i + 1];
1572 ae = 0.5 * (t0 + te);
1573 be = PetscPowScalar(ae, bm1);
1574 /* de = be * ae; */
1575 de = PetscPowScalar(ae, beta);
1576 ge = coef * be * (te - t0);
1577
1578 ts = x[k][j - 1][i];
1579 as = 0.5 * (t0 + ts);
1580 bs = PetscPowScalar(as, bm1);
1581 /* ds = bs * as; */
1582 ds = PetscPowScalar(as, beta);
1583 gs = coef * bs * (t0 - ts);
1584
1585 tn = x[k][j + 1][i];
1586 an = 0.5 * (t0 + tn);
1587 bn = PetscPowScalar(an, bm1);
1588 /* dn = bn * an; */
1589 dn = PetscPowScalar(an, beta);
1590 gn = coef * bn * (tn - t0);
1591
1592 td = x[k - 1][j][i];
1593 ad = 0.5 * (t0 + td);
1594 bd = PetscPowScalar(ad, bm1);
1595 /* dd = bd * ad; */
1596 dd = PetscPowScalar(ad, beta);
1597 gd = coef * bd * (t0 - td);
1598
1599 c[0].k = k - 1;
1600 c[0].j = j;
1601 c[0].i = i;
1602 v[0] = -hxhydhz * (dd - gd);
1603 c[1].k = k;
1604 c[1].j = j - 1;
1605 c[1].i = i;
1606 v[1] = -hzhxdhy * (ds - gs);
1607 c[2].k = k;
1608 c[2].j = j;
1609 c[2].i = i - 1;
1610 v[2] = -hyhzdhx * (dw - gw);
1611 c[3].k = k;
1612 c[3].j = j;
1613 c[3].i = i;
1614 v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd);
1615 c[4].k = k;
1616 c[4].j = j;
1617 c[4].i = i + 1;
1618 v[4] = -hyhzdhx * (de + ge);
1619 c[5].k = k;
1620 c[5].j = j + 1;
1621 c[5].i = i;
1622 v[5] = -hzhxdhy * (dn + gn);
1623 PetscCall(MatSetValuesStencil(jac, 1, &row, 6, c, v, INSERT_VALUES));
1624 }
1625 }
1626 }
1627 }
1628 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
1629 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
1630 if (jac != J) {
1631 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1632 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1633 }
1634 PetscCall(DMDAVecRestoreArray(da, localX, &x));
1635 PetscCall(DMRestoreLocalVector(da, &localX));
1636
1637 PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym));
1638 PetscFunctionReturn(PETSC_SUCCESS);
1639 }
1640
TestConvergence(SNES snes,PETSC_UNUSED PetscInt it,PETSC_UNUSED PetscReal xnorm,PETSC_UNUSED PetscReal snorm,PETSC_UNUSED PetscReal fnorm,SNESConvergedReason * reason,PetscCtx ctx)1641 PetscErrorCode TestConvergence(SNES snes, PETSC_UNUSED PetscInt it, PETSC_UNUSED PetscReal xnorm, PETSC_UNUSED PetscReal snorm, PETSC_UNUSED PetscReal fnorm, SNESConvergedReason *reason, PetscCtx ctx)
1642 {
1643 AppCtx *user = (AppCtx *)ctx;
1644
1645 PetscFunctionBeginUser;
1646 if (user->converged) *reason = SNES_CONVERGED_USER;
1647 else *reason = SNES_DIVERGED_USER;
1648 PetscFunctionReturn(PETSC_SUCCESS);
1649 }
1650
1651 /*TEST
1652
1653 test:
1654 nsize: 4
1655 args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat
1656 requires: !single
1657
1658 test:
1659 suffix: 2
1660 args: -snes_converged_reason -use_convergence_test -mark_converged
1661
1662 test:
1663 suffix: 3
1664 args: -snes_converged_reason -use_convergence_test -mark_converged 0
1665
1666 TEST*/
1667