Lines Matching +full:- +full:j
2 Uses 3-dimensional distributed arrays.\n\
3 A 3-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \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";
17 - Div(alpha* T^beta (GRAD T)) = 0.
26 A finite volume approximation with the usual 7-point stencil
38 /* User-defined application context */
68 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL)); in main()
69 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL)); in main()
70 PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL)); in main()
71 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL)); in main()
72 PetscCall(PetscOptionsGetBool(NULL, NULL, "-mark_converged", &user.converged, NULL)); in main()
73 user.bm1 = user.beta - 1.0; in main()
108 /* -------------------- Form initial approximation ----------------- */
112 PetscInt i, j, k, xs, ys, xm, ym, zs, zm; in FormInitialGuess() local
124 for (j = ys; j < ys + ym; j++) { in FormInitialGuess()
125 for (i = xs; i < xs + xm; i++) x[k][j][i] = user->tleft; in FormInitialGuess()
131 /* -------------------- Evaluate Function F(x) --------------------- */
135 PetscInt i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm; in FormFunction() local
148 hx = one / (PetscReal)(mx - 1); in FormFunction()
149 hy = one / (PetscReal)(my - 1); in FormFunction()
150 hz = one / (PetscReal)(mz - 1); in FormFunction()
154 tleft = user->tleft; in FormFunction()
155 tright = user->tright; in FormFunction()
156 beta = user->beta; in FormFunction()
167 for (j = ys; j < ys + ym; j++) { in FormFunction()
169 t0 = x[k][j][i]; in FormFunction()
171 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) { in FormFunction()
174 tw = x[k][j][i - 1]; in FormFunction()
177 fw = dw * (t0 - tw); in FormFunction()
179 te = x[k][j][i + 1]; in FormFunction()
182 fe = de * (te - t0); in FormFunction()
184 ts = x[k][j - 1][i]; in FormFunction()
187 fs = ds * (t0 - ts); in FormFunction()
189 tn = x[k][j + 1][i]; in FormFunction()
192 fn = dn * (tn - t0); in FormFunction()
194 td = x[k - 1][j][i]; in FormFunction()
197 fd = dd * (t0 - td); in FormFunction()
199 tu = x[k + 1][j][i]; in FormFunction()
202 fu = du * (tu - t0); in FormFunction()
205 /* left-hand (west) boundary */ in FormFunction()
209 fw = dw * (t0 - tw); in FormFunction()
211 te = x[k][j][i + 1]; in FormFunction()
214 fe = de * (te - t0); in FormFunction()
216 if (j > 0) { in FormFunction()
217 ts = x[k][j - 1][i]; in FormFunction()
220 fs = ds * (t0 - ts); in FormFunction()
225 if (j < my - 1) { in FormFunction()
226 tn = x[k][j + 1][i]; in FormFunction()
229 fn = dn * (tn - t0); in FormFunction()
235 td = x[k - 1][j][i]; in FormFunction()
238 fd = dd * (t0 - td); in FormFunction()
243 if (k < mz - 1) { in FormFunction()
244 tu = x[k + 1][j][i]; in FormFunction()
247 fu = du * (tu - t0); in FormFunction()
252 } else if (i == mx - 1) { in FormFunction()
253 /* right-hand (east) boundary */ in FormFunction()
254 tw = x[k][j][i - 1]; in FormFunction()
257 fw = dw * (t0 - tw); in FormFunction()
262 fe = de * (te - t0); in FormFunction()
264 if (j > 0) { in FormFunction()
265 ts = x[k][j - 1][i]; in FormFunction()
268 fs = ds * (t0 - ts); in FormFunction()
273 if (j < my - 1) { in FormFunction()
274 tn = x[k][j + 1][i]; in FormFunction()
277 fn = dn * (tn - t0); in FormFunction()
283 td = x[k - 1][j][i]; in FormFunction()
286 fd = dd * (t0 - td); in FormFunction()
291 if (k < mz - 1) { in FormFunction()
292 tu = x[k + 1][j][i]; in FormFunction()
295 fu = du * (tu - t0); in FormFunction()
300 } else if (j == 0) { in FormFunction()
301 /* bottom (south) boundary, and i <> 0 or mx-1 */ in FormFunction()
302 tw = x[k][j][i - 1]; in FormFunction()
305 fw = dw * (t0 - tw); in FormFunction()
307 te = x[k][j][i + 1]; in FormFunction()
310 fe = de * (te - t0); in FormFunction()
314 tn = x[k][j + 1][i]; in FormFunction()
317 fn = dn * (tn - t0); in FormFunction()
320 td = x[k - 1][j][i]; in FormFunction()
323 fd = dd * (t0 - td); in FormFunction()
328 if (k < mz - 1) { in FormFunction()
329 tu = x[k + 1][j][i]; in FormFunction()
332 fu = du * (tu - t0); in FormFunction()
337 } else if (j == my - 1) { in FormFunction()
338 /* top (north) boundary, and i <> 0 or mx-1 */ in FormFunction()
339 tw = x[k][j][i - 1]; in FormFunction()
342 fw = dw * (t0 - tw); in FormFunction()
344 te = x[k][j][i + 1]; in FormFunction()
347 fe = de * (te - t0); in FormFunction()
349 ts = x[k][j - 1][i]; in FormFunction()
352 fs = ds * (t0 - ts); in FormFunction()
357 td = x[k - 1][j][i]; in FormFunction()
360 fd = dd * (t0 - td); in FormFunction()
365 if (k < mz - 1) { in FormFunction()
366 tu = x[k + 1][j][i]; in FormFunction()
369 fu = du * (tu - t0); in FormFunction()
376 tw = x[k][j][i - 1]; in FormFunction()
379 fw = dw * (t0 - tw); in FormFunction()
381 te = x[k][j][i + 1]; in FormFunction()
384 fe = de * (te - t0); in FormFunction()
386 ts = x[k][j - 1][i]; in FormFunction()
389 fs = ds * (t0 - ts); in FormFunction()
391 tn = x[k][j + 1][i]; in FormFunction()
394 fn = dn * (tn - t0); in FormFunction()
398 tu = x[k + 1][j][i]; in FormFunction()
401 fu = du * (tu - t0); in FormFunction()
403 } else if (k == mz - 1) { in FormFunction()
405 tw = x[k][j][i - 1]; in FormFunction()
408 fw = dw * (t0 - tw); in FormFunction()
410 te = x[k][j][i + 1]; in FormFunction()
413 fe = de * (te - t0); in FormFunction()
415 ts = x[k][j - 1][i]; in FormFunction()
418 fs = ds * (t0 - ts); in FormFunction()
420 tn = x[k][j + 1][i]; in FormFunction()
423 fn = dn * (tn - t0); in FormFunction()
425 td = x[k - 1][j][i]; in FormFunction()
428 fd = dd * (t0 - td); in FormFunction()
433 f[k][j][i] = -hyhzdhx * (fe - fw) - hzhxdhy * (fn - fs) - hxhydhz * (fu - fd); in FormFunction()
443 /* -------------------- Evaluate Jacobian F(x) --------------------- */
444 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr) in FormJacobian() argument
447 PetscInt i, j, k, mx, my, mz, xs, ys, zs, xm, ym, zm; in FormJacobian() local
461 hx = one / (PetscReal)(mx - 1); in FormJacobian()
462 hy = one / (PetscReal)(my - 1); in FormJacobian()
463 hz = one / (PetscReal)(mz - 1); in FormJacobian()
467 tleft = user->tleft; in FormJacobian()
468 tright = user->tright; in FormJacobian()
469 beta = user->beta; in FormJacobian()
470 bm1 = user->bm1; in FormJacobian()
471 coef = user->coef; in FormJacobian()
481 for (j = ys; j < ys + ym; j++) { in FormJacobian()
483 t0 = x[k][j][i]; in FormJacobian()
485 row.j = j; in FormJacobian()
487 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1 && k > 0 && k < mz - 1) { in FormJacobian()
490 tw = x[k][j][i - 1]; in FormJacobian()
495 gw = coef * bw * (t0 - tw); in FormJacobian()
497 te = x[k][j][i + 1]; in FormJacobian()
502 ge = coef * be * (te - t0); in FormJacobian()
504 ts = x[k][j - 1][i]; in FormJacobian()
509 gs = coef * bs * (t0 - ts); in FormJacobian()
511 tn = x[k][j + 1][i]; in FormJacobian()
516 gn = coef * bn * (tn - t0); in FormJacobian()
518 td = x[k - 1][j][i]; in FormJacobian()
523 gd = coef * bd * (t0 - td); in FormJacobian()
525 tu = x[k + 1][j][i]; in FormJacobian()
530 gu = coef * bu * (tu - t0); in FormJacobian()
532 c[0].k = k - 1; in FormJacobian()
533 c[0].j = j; in FormJacobian()
535 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
537 c[1].j = j - 1; in FormJacobian()
539 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
541 c[2].j = j; in FormJacobian()
542 c[2].i = i - 1; in FormJacobian()
543 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
545 c[3].j = j; in FormJacobian()
547 …v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd -… in FormJacobian()
549 c[4].j = j; in FormJacobian()
551 v[4] = -hyhzdhx * (de + ge); in FormJacobian()
553 c[5].j = j + 1; in FormJacobian()
555 v[5] = -hzhxdhy * (dn + gn); in FormJacobian()
557 c[6].j = j; in FormJacobian()
559 v[6] = -hxhydhz * (du + gu); in FormJacobian()
563 /* left-hand plane boundary */ in FormJacobian()
569 gw = coef * bw * (t0 - tw); in FormJacobian()
571 te = x[k][j][i + 1]; in FormJacobian()
576 ge = coef * be * (te - t0); in FormJacobian()
578 /* left-hand bottom edge */ in FormJacobian()
579 if (j == 0) { in FormJacobian()
580 tn = x[k][j + 1][i]; in FormJacobian()
585 gn = coef * bn * (tn - t0); in FormJacobian()
587 /* left-hand bottom down corner */ in FormJacobian()
589 tu = x[k + 1][j][i]; in FormJacobian()
594 gu = coef * bu * (tu - t0); in FormJacobian()
597 c[0].j = j; in FormJacobian()
599 v[0] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
601 c[1].j = j; in FormJacobian()
603 v[1] = -hyhzdhx * (de + ge); in FormJacobian()
605 c[2].j = j + 1; in FormJacobian()
607 v[2] = -hzhxdhy * (dn + gn); in FormJacobian()
609 c[3].j = j; in FormJacobian()
611 v[3] = -hxhydhz * (du + gu); in FormJacobian()
614 /* left-hand bottom interior edge */ in FormJacobian()
615 } else if (k < mz - 1) { in FormJacobian()
616 tu = x[k + 1][j][i]; in FormJacobian()
621 gu = coef * bu * (tu - t0); in FormJacobian()
623 td = x[k - 1][j][i]; in FormJacobian()
628 gd = coef * bd * (td - t0); in FormJacobian()
630 c[0].k = k - 1; in FormJacobian()
631 c[0].j = j; in FormJacobian()
633 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
635 c[1].j = j; in FormJacobian()
637 … v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu); in FormJacobian()
639 c[2].j = j; in FormJacobian()
641 v[2] = -hyhzdhx * (de + ge); in FormJacobian()
643 c[3].j = j + 1; in FormJacobian()
645 v[3] = -hzhxdhy * (dn + gn); in FormJacobian()
647 c[4].j = j; in FormJacobian()
649 v[4] = -hxhydhz * (du + gu); in FormJacobian()
652 /* left-hand bottom up corner */ in FormJacobian()
654 td = x[k - 1][j][i]; in FormJacobian()
659 gd = coef * bd * (td - t0); in FormJacobian()
661 c[0].k = k - 1; in FormJacobian()
662 c[0].j = j; in FormJacobian()
664 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
666 c[1].j = j; in FormJacobian()
668 v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
670 c[2].j = j; in FormJacobian()
672 v[2] = -hyhzdhx * (de + ge); in FormJacobian()
674 c[3].j = j + 1; in FormJacobian()
676 v[3] = -hzhxdhy * (dn + gn); in FormJacobian()
680 /* left-hand top edge */ in FormJacobian()
681 } else if (j == my - 1) { in FormJacobian()
682 ts = x[k][j - 1][i]; in FormJacobian()
687 gs = coef * bs * (ts - t0); in FormJacobian()
689 /* left-hand top down corner */ in FormJacobian()
691 tu = x[k + 1][j][i]; in FormJacobian()
696 gu = coef * bu * (tu - t0); in FormJacobian()
699 c[0].j = j - 1; in FormJacobian()
701 v[0] = -hzhxdhy * (ds - gs); in FormJacobian()
703 c[1].j = j; in FormJacobian()
705 v[1] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
707 c[2].j = j; in FormJacobian()
709 v[2] = -hyhzdhx * (de + ge); in FormJacobian()
711 c[3].j = j; in FormJacobian()
713 v[3] = -hxhydhz * (du + gu); in FormJacobian()
716 /* left-hand top interior edge */ in FormJacobian()
717 } else if (k < mz - 1) { in FormJacobian()
718 tu = x[k + 1][j][i]; in FormJacobian()
723 gu = coef * bu * (tu - t0); in FormJacobian()
725 td = x[k - 1][j][i]; in FormJacobian()
730 gd = coef * bd * (td - t0); in FormJacobian()
732 c[0].k = k - 1; in FormJacobian()
733 c[0].j = j; in FormJacobian()
735 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
737 c[1].j = j - 1; in FormJacobian()
739 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
741 c[2].j = j; in FormJacobian()
743 … v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu); in FormJacobian()
745 c[3].j = j; in FormJacobian()
747 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
749 c[4].j = j; in FormJacobian()
751 v[4] = -hxhydhz * (du + gu); in FormJacobian()
754 /* left-hand top up corner */ in FormJacobian()
756 td = x[k - 1][j][i]; in FormJacobian()
761 gd = coef * bd * (td - t0); in FormJacobian()
763 c[0].k = k - 1; in FormJacobian()
764 c[0].j = j; in FormJacobian()
766 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
768 c[1].j = j - 1; in FormJacobian()
770 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
772 c[2].j = j; in FormJacobian()
774 v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
776 c[3].j = j; in FormJacobian()
778 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
783 ts = x[k][j - 1][i]; in FormJacobian()
788 gs = coef * bs * (t0 - ts); in FormJacobian()
790 tn = x[k][j + 1][i]; in FormJacobian()
795 gn = coef * bn * (tn - t0); in FormJacobian()
797 /* left-hand down interior edge */ in FormJacobian()
799 tu = x[k + 1][j][i]; in FormJacobian()
804 gu = coef * bu * (tu - t0); in FormJacobian()
807 c[0].j = j - 1; in FormJacobian()
809 v[0] = -hzhxdhy * (ds - gs); in FormJacobian()
811 c[1].j = j; in FormJacobian()
813 … v[1] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
815 c[2].j = j; in FormJacobian()
817 v[2] = -hyhzdhx * (de + ge); in FormJacobian()
819 c[3].j = j + 1; in FormJacobian()
821 v[3] = -hzhxdhy * (dn + gn); in FormJacobian()
823 c[4].j = j; in FormJacobian()
825 v[4] = -hxhydhz * (du + gu); in FormJacobian()
828 } else if (k == mz - 1) { /* left-hand up interior edge */ in FormJacobian()
830 td = x[k - 1][j][i]; in FormJacobian()
835 gd = coef * bd * (t0 - td); in FormJacobian()
837 c[0].k = k - 1; in FormJacobian()
838 c[0].j = j; in FormJacobian()
840 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
842 c[1].j = j - 1; in FormJacobian()
844 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
846 c[2].j = j; in FormJacobian()
848 … v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
850 c[3].j = j; in FormJacobian()
852 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
854 c[4].j = j + 1; in FormJacobian()
856 v[4] = -hzhxdhy * (dn + gn); in FormJacobian()
858 } else { /* left-hand interior plane */ in FormJacobian()
860 td = x[k - 1][j][i]; in FormJacobian()
865 gd = coef * bd * (t0 - td); in FormJacobian()
867 tu = x[k + 1][j][i]; in FormJacobian()
872 gu = coef * bu * (tu - t0); in FormJacobian()
874 c[0].k = k - 1; in FormJacobian()
875 c[0].j = j; in FormJacobian()
877 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
879 c[1].j = j - 1; in FormJacobian()
881 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
883 c[2].j = j; in FormJacobian()
885 …v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd -… in FormJacobian()
887 c[3].j = j; in FormJacobian()
889 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
891 c[4].j = j + 1; in FormJacobian()
893 v[4] = -hzhxdhy * (dn + gn); in FormJacobian()
895 c[5].j = j; in FormJacobian()
897 v[5] = -hxhydhz * (du + gu); in FormJacobian()
902 } else if (i == mx - 1) { in FormJacobian()
903 /* right-hand plane boundary */ in FormJacobian()
904 tw = x[k][j][i - 1]; in FormJacobian()
909 gw = coef * bw * (t0 - tw); in FormJacobian()
916 ge = coef * be * (te - t0); in FormJacobian()
918 /* right-hand bottom edge */ in FormJacobian()
919 if (j == 0) { in FormJacobian()
920 tn = x[k][j + 1][i]; in FormJacobian()
925 gn = coef * bn * (tn - t0); in FormJacobian()
927 /* right-hand bottom down corner */ in FormJacobian()
929 tu = x[k + 1][j][i]; in FormJacobian()
934 gu = coef * bu * (tu - t0); in FormJacobian()
937 c[0].j = j; in FormJacobian()
938 c[0].i = i - 1; in FormJacobian()
939 v[0] = -hyhzdhx * (dw - gw); in FormJacobian()
941 c[1].j = j; in FormJacobian()
943 v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
945 c[2].j = j + 1; in FormJacobian()
947 v[2] = -hzhxdhy * (dn + gn); in FormJacobian()
949 c[3].j = j; in FormJacobian()
951 v[3] = -hxhydhz * (du + gu); in FormJacobian()
954 /* right-hand bottom interior edge */ in FormJacobian()
955 } else if (k < mz - 1) { in FormJacobian()
956 tu = x[k + 1][j][i]; in FormJacobian()
961 gu = coef * bu * (tu - t0); in FormJacobian()
963 td = x[k - 1][j][i]; in FormJacobian()
968 gd = coef * bd * (td - t0); in FormJacobian()
970 c[0].k = k - 1; in FormJacobian()
971 c[0].j = j; in FormJacobian()
973 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
975 c[1].j = j; in FormJacobian()
976 c[1].i = i - 1; in FormJacobian()
977 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
979 c[2].j = j; in FormJacobian()
981 … v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu); in FormJacobian()
983 c[3].j = j + 1; in FormJacobian()
985 v[3] = -hzhxdhy * (dn + gn); in FormJacobian()
987 c[4].j = j; in FormJacobian()
989 v[4] = -hxhydhz * (du + gu); in FormJacobian()
992 /* right-hand bottom up corner */ in FormJacobian()
994 td = x[k - 1][j][i]; in FormJacobian()
999 gd = coef * bd * (td - t0); in FormJacobian()
1001 c[0].k = k - 1; in FormJacobian()
1002 c[0].j = j; in FormJacobian()
1004 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1006 c[1].j = j; in FormJacobian()
1007 c[1].i = i - 1; in FormJacobian()
1008 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
1010 c[2].j = j; in FormJacobian()
1012 v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
1014 c[3].j = j + 1; in FormJacobian()
1016 v[3] = -hzhxdhy * (dn + gn); in FormJacobian()
1020 /* right-hand top edge */ in FormJacobian()
1021 } else if (j == my - 1) { in FormJacobian()
1022 ts = x[k][j - 1][i]; in FormJacobian()
1027 gs = coef * bs * (ts - t0); in FormJacobian()
1029 /* right-hand top down corner */ in FormJacobian()
1031 tu = x[k + 1][j][i]; in FormJacobian()
1036 gu = coef * bu * (tu - t0); in FormJacobian()
1039 c[0].j = j - 1; in FormJacobian()
1041 v[0] = -hzhxdhy * (ds - gs); in FormJacobian()
1043 c[1].j = j; in FormJacobian()
1044 c[1].i = i - 1; in FormJacobian()
1045 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
1047 c[2].j = j; in FormJacobian()
1049 v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
1051 c[3].j = j; in FormJacobian()
1053 v[3] = -hxhydhz * (du + gu); in FormJacobian()
1056 /* right-hand top interior edge */ in FormJacobian()
1057 } else if (k < mz - 1) { in FormJacobian()
1058 tu = x[k + 1][j][i]; in FormJacobian()
1063 gu = coef * bu * (tu - t0); in FormJacobian()
1065 td = x[k - 1][j][i]; in FormJacobian()
1070 gd = coef * bd * (td - t0); in FormJacobian()
1072 c[0].k = k - 1; in FormJacobian()
1073 c[0].j = j; in FormJacobian()
1075 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1077 c[1].j = j - 1; in FormJacobian()
1079 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
1081 c[2].j = j; in FormJacobian()
1082 c[2].i = i - 1; in FormJacobian()
1083 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
1085 c[3].j = j; in FormJacobian()
1087 … v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu); in FormJacobian()
1089 c[4].j = j; in FormJacobian()
1091 v[4] = -hxhydhz * (du + gu); in FormJacobian()
1094 /* right-hand top up corner */ in FormJacobian()
1096 td = x[k - 1][j][i]; in FormJacobian()
1101 gd = coef * bd * (td - t0); in FormJacobian()
1103 c[0].k = k - 1; in FormJacobian()
1104 c[0].j = j; in FormJacobian()
1106 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1108 c[1].j = j - 1; in FormJacobian()
1110 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
1112 c[2].j = j; in FormJacobian()
1113 c[2].i = i - 1; in FormJacobian()
1114 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
1116 c[3].j = j; in FormJacobian()
1118 v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
1123 ts = x[k][j - 1][i]; in FormJacobian()
1128 gs = coef * bs * (t0 - ts); in FormJacobian()
1130 tn = x[k][j + 1][i]; in FormJacobian()
1135 gn = coef * bn * (tn - t0); in FormJacobian()
1137 /* right-hand down interior edge */ in FormJacobian()
1139 tu = x[k + 1][j][i]; in FormJacobian()
1144 gu = coef * bu * (tu - t0); in FormJacobian()
1147 c[0].j = j - 1; in FormJacobian()
1149 v[0] = -hzhxdhy * (ds - gs); in FormJacobian()
1151 c[1].j = j; in FormJacobian()
1152 c[1].i = i - 1; in FormJacobian()
1153 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
1155 c[2].j = j; in FormJacobian()
1157 … v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
1159 c[3].j = j + 1; in FormJacobian()
1161 v[3] = -hzhxdhy * (dn + gn); in FormJacobian()
1163 c[4].j = j; in FormJacobian()
1165 v[4] = -hxhydhz * (du + gu); in FormJacobian()
1168 } else if (k == mz - 1) { /* right-hand up interior edge */ in FormJacobian()
1170 td = x[k - 1][j][i]; in FormJacobian()
1175 gd = coef * bd * (t0 - td); in FormJacobian()
1177 c[0].k = k - 1; in FormJacobian()
1178 c[0].j = j; in FormJacobian()
1180 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1182 c[1].j = j - 1; in FormJacobian()
1184 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
1186 c[2].j = j; in FormJacobian()
1187 c[2].i = i - 1; in FormJacobian()
1188 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
1190 c[3].j = j; in FormJacobian()
1192 … v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
1194 c[4].j = j + 1; in FormJacobian()
1196 v[4] = -hzhxdhy * (dn + gn); in FormJacobian()
1199 } else { /* right-hand interior plane */ in FormJacobian()
1201 td = x[k - 1][j][i]; in FormJacobian()
1206 gd = coef * bd * (t0 - td); in FormJacobian()
1208 tu = x[k + 1][j][i]; in FormJacobian()
1213 gu = coef * bu * (tu - t0); in FormJacobian()
1215 c[0].k = k - 1; in FormJacobian()
1216 c[0].j = j; in FormJacobian()
1218 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1220 c[1].j = j - 1; in FormJacobian()
1222 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
1224 c[2].j = j; in FormJacobian()
1225 c[2].i = i - 1; in FormJacobian()
1226 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
1228 c[3].j = j; in FormJacobian()
1230 …v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd -… in FormJacobian()
1232 c[4].j = j + 1; in FormJacobian()
1234 v[4] = -hzhxdhy * (dn + gn); in FormJacobian()
1236 c[5].j = j; in FormJacobian()
1238 v[5] = -hxhydhz * (du + gu); in FormJacobian()
1243 } else if (j == 0) { in FormJacobian()
1244 tw = x[k][j][i - 1]; in FormJacobian()
1249 gw = coef * bw * (t0 - tw); in FormJacobian()
1251 te = x[k][j][i + 1]; in FormJacobian()
1256 ge = coef * be * (te - t0); in FormJacobian()
1258 tn = x[k][j + 1][i]; in FormJacobian()
1263 gn = coef * bn * (tn - t0); in FormJacobian()
1267 tu = x[k + 1][j][i]; in FormJacobian()
1272 gu = coef * bu * (tu - t0); in FormJacobian()
1275 c[0].j = j; in FormJacobian()
1276 c[0].i = i - 1; in FormJacobian()
1277 v[0] = -hyhzdhx * (dw - gw); in FormJacobian()
1279 c[1].j = j; in FormJacobian()
1281 v[1] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
1283 c[2].j = j; in FormJacobian()
1285 v[2] = -hyhzdhx * (de + ge); in FormJacobian()
1287 c[3].j = j + 1; in FormJacobian()
1289 v[3] = -hzhxdhy * (dn + gn); in FormJacobian()
1291 c[4].j = j; in FormJacobian()
1293 v[4] = -hxhydhz * (du + gu); in FormJacobian()
1296 } else if (k == mz - 1) { /* bottom up interior edge */ in FormJacobian()
1298 td = x[k - 1][j][i]; in FormJacobian()
1303 gd = coef * bd * (td - t0); in FormJacobian()
1305 c[0].k = k - 1; in FormJacobian()
1306 c[0].j = j; in FormJacobian()
1308 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1310 c[1].j = j; in FormJacobian()
1311 c[1].i = i - 1; in FormJacobian()
1312 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
1314 c[2].j = j; in FormJacobian()
1316 v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
1318 c[3].j = j; in FormJacobian()
1320 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
1322 c[4].j = j + 1; in FormJacobian()
1324 v[4] = -hzhxdhy * (dn + gn); in FormJacobian()
1329 tu = x[k + 1][j][i]; in FormJacobian()
1334 gu = coef * bu * (tu - t0); in FormJacobian()
1336 td = x[k - 1][j][i]; in FormJacobian()
1341 gd = coef * bd * (td - t0); in FormJacobian()
1343 c[0].k = k - 1; in FormJacobian()
1344 c[0].j = j; in FormJacobian()
1346 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1348 c[1].j = j; in FormJacobian()
1349 c[1].i = i - 1; in FormJacobian()
1350 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
1352 c[2].j = j; in FormJacobian()
1354 … v[2] = hzhxdhy * (dn - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu); in FormJacobian()
1356 c[3].j = j; in FormJacobian()
1358 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
1360 c[4].j = j + 1; in FormJacobian()
1362 v[4] = -hzhxdhy * (dn + gn); in FormJacobian()
1364 c[5].j = j; in FormJacobian()
1366 v[5] = -hxhydhz * (du + gu); in FormJacobian()
1370 } else if (j == my - 1) { in FormJacobian()
1371 tw = x[k][j][i - 1]; in FormJacobian()
1376 gw = coef * bw * (t0 - tw); in FormJacobian()
1378 te = x[k][j][i + 1]; in FormJacobian()
1383 ge = coef * be * (te - t0); in FormJacobian()
1385 ts = x[k][j - 1][i]; in FormJacobian()
1390 gs = coef * bs * (t0 - ts); in FormJacobian()
1394 tu = x[k + 1][j][i]; in FormJacobian()
1399 gu = coef * bu * (tu - t0); in FormJacobian()
1402 c[0].j = j - 1; in FormJacobian()
1404 v[0] = -hzhxdhy * (ds - gs); in FormJacobian()
1406 c[1].j = j; in FormJacobian()
1407 c[1].i = i - 1; in FormJacobian()
1408 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
1410 c[2].j = j; in FormJacobian()
1412 v[2] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
1414 c[3].j = j; in FormJacobian()
1416 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
1418 c[4].j = j; in FormJacobian()
1420 v[4] = -hxhydhz * (du + gu); in FormJacobian()
1423 } else if (k == mz - 1) { /* top up interior edge */ in FormJacobian()
1425 td = x[k - 1][j][i]; in FormJacobian()
1430 gd = coef * bd * (td - t0); in FormJacobian()
1432 c[0].k = k - 1; in FormJacobian()
1433 c[0].j = j; in FormJacobian()
1435 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1437 c[1].j = j - 1; in FormJacobian()
1439 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
1441 c[2].j = j; in FormJacobian()
1442 c[2].i = i - 1; in FormJacobian()
1443 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
1445 c[3].j = j; in FormJacobian()
1447 v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
1449 c[4].j = j; in FormJacobian()
1451 v[4] = -hyhzdhx * (de + ge); in FormJacobian()
1456 tu = x[k + 1][j][i]; in FormJacobian()
1461 gu = coef * bu * (tu - t0); in FormJacobian()
1463 td = x[k - 1][j][i]; in FormJacobian()
1468 gd = coef * bd * (td - t0); in FormJacobian()
1470 c[0].k = k - 1; in FormJacobian()
1471 c[0].j = j; in FormJacobian()
1473 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1475 c[1].j = j - 1; in FormJacobian()
1477 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
1479 c[2].j = j; in FormJacobian()
1480 c[2].i = i - 1; in FormJacobian()
1481 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
1483 c[3].j = j; in FormJacobian()
1485 … v[3] = hzhxdhy * (ds + gs) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + du + gd - gu); in FormJacobian()
1487 c[4].j = j; in FormJacobian()
1489 v[4] = -hyhzdhx * (de + ge); in FormJacobian()
1491 c[5].j = j; in FormJacobian()
1493 v[5] = -hxhydhz * (du + gu); in FormJacobian()
1500 tw = x[k][j][i - 1]; in FormJacobian()
1505 gw = coef * bw * (t0 - tw); in FormJacobian()
1507 te = x[k][j][i + 1]; in FormJacobian()
1512 ge = coef * be * (te - t0); in FormJacobian()
1514 ts = x[k][j - 1][i]; in FormJacobian()
1519 gs = coef * bs * (t0 - ts); in FormJacobian()
1521 tn = x[k][j + 1][i]; in FormJacobian()
1526 gn = coef * bn * (tn - t0); in FormJacobian()
1528 tu = x[k + 1][j][i]; in FormJacobian()
1533 gu = coef * bu * (tu - t0); in FormJacobian()
1536 c[0].j = j - 1; in FormJacobian()
1538 v[0] = -hzhxdhy * (ds - gs); in FormJacobian()
1540 c[1].j = j; in FormJacobian()
1541 c[1].i = i - 1; in FormJacobian()
1542 v[1] = -hyhzdhx * (dw - gw); in FormJacobian()
1544 c[2].j = j; in FormJacobian()
1546 … v[2] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (du - gu); in FormJacobian()
1548 c[3].j = j; in FormJacobian()
1550 v[3] = -hyhzdhx * (de + ge); in FormJacobian()
1552 c[4].j = j + 1; in FormJacobian()
1554 v[4] = -hzhxdhy * (dn + gn); in FormJacobian()
1556 c[5].j = j; in FormJacobian()
1558 v[5] = -hxhydhz * (du + gu); in FormJacobian()
1561 } else if (k == mz - 1) { in FormJacobian()
1564 tw = x[k][j][i - 1]; in FormJacobian()
1569 gw = coef * bw * (t0 - tw); in FormJacobian()
1571 te = x[k][j][i + 1]; in FormJacobian()
1576 ge = coef * be * (te - t0); in FormJacobian()
1578 ts = x[k][j - 1][i]; in FormJacobian()
1583 gs = coef * bs * (t0 - ts); in FormJacobian()
1585 tn = x[k][j + 1][i]; in FormJacobian()
1590 gn = coef * bn * (tn - t0); in FormJacobian()
1592 td = x[k - 1][j][i]; in FormJacobian()
1597 gd = coef * bd * (t0 - td); in FormJacobian()
1599 c[0].k = k - 1; in FormJacobian()
1600 c[0].j = j; in FormJacobian()
1602 v[0] = -hxhydhz * (dd - gd); in FormJacobian()
1604 c[1].j = j - 1; in FormJacobian()
1606 v[1] = -hzhxdhy * (ds - gs); in FormJacobian()
1608 c[2].j = j; in FormJacobian()
1609 c[2].i = i - 1; in FormJacobian()
1610 v[2] = -hyhzdhx * (dw - gw); in FormJacobian()
1612 c[3].j = j; in FormJacobian()
1614 … v[3] = hzhxdhy * (ds + dn + gs - gn) + hyhzdhx * (dw + de + gw - ge) + hxhydhz * (dd + gd); in FormJacobian()
1616 c[4].j = j; in FormJacobian()
1618 v[4] = -hyhzdhx * (de + ge); in FormJacobian()
1620 c[5].j = j + 1; in FormJacobian()
1622 v[5] = -hzhxdhy * (dn + gn); in FormJacobian()
1630 if (jac != J) { in FormJacobian()
1631 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); in FormJacobian()
1632 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); in FormJacobian()
1646 if (user->converged) *reason = SNES_CONVERGED_USER; in TestConvergence()
1655 …args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2…
1660 args: -snes_converged_reason -use_convergence_test -mark_converged
1664 args: -snes_converged_reason -use_convergence_test -mark_converged 0