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