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