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 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 ----------------- */ 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) --------------------- */ 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) --------------------- */ 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 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