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