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++) { 121 x[k][j][i] = user->tleft; 122 } 123 } 124 } 125 PetscCall(DMDAVecRestoreArray(da,X,&x)); 126 PetscFunctionReturn(0); 127 } 128 /* -------------------- Evaluate Function F(x) --------------------- */ 129 PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr) 130 { 131 AppCtx *user = (AppCtx*)ptr; 132 PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm; 133 PetscScalar zero = 0.0,one = 1.0; 134 PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy; 135 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; 136 PetscScalar tleft,tright,beta,td,ad,dd,fd=0.0,tu,au,du=0.0,fu=0.0; 137 PetscScalar ***x,***f; 138 Vec localX; 139 DM da; 140 141 PetscFunctionBeginUser; 142 PetscCall(SNESGetDM(snes,&da)); 143 PetscCall(DMGetLocalVector(da,&localX)); 144 PetscCall(DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0)); 145 hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1); 146 hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy; 147 tleft = user->tleft; 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 165 /* general interior volume */ 166 167 tw = x[k][j][i-1]; 168 aw = 0.5*(t0 + tw); 169 dw = PetscPowScalar(aw,beta); 170 fw = dw*(t0 - tw); 171 172 te = x[k][j][i+1]; 173 ae = 0.5*(t0 + te); 174 de = PetscPowScalar(ae,beta); 175 fe = de*(te - t0); 176 177 ts = x[k][j-1][i]; 178 as = 0.5*(t0 + ts); 179 ds = PetscPowScalar(as,beta); 180 fs = ds*(t0 - ts); 181 182 tn = x[k][j+1][i]; 183 an = 0.5*(t0 + tn); 184 dn = PetscPowScalar(an,beta); 185 fn = dn*(tn - t0); 186 187 td = x[k-1][j][i]; 188 ad = 0.5*(t0 + td); 189 dd = PetscPowScalar(ad,beta); 190 fd = dd*(t0 - td); 191 192 tu = x[k+1][j][i]; 193 au = 0.5*(t0 + tu); 194 du = PetscPowScalar(au,beta); 195 fu = du*(tu - t0); 196 197 } else if (i == 0) { 198 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 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 297 /* bottom (south) boundary, and i <> 0 or mx-1 */ 298 tw = x[k][j][i-1]; 299 aw = 0.5*(t0 + tw); 300 dw = PetscPowScalar(aw,beta); 301 fw = dw*(t0 - tw); 302 303 te = x[k][j][i+1]; 304 ae = 0.5*(t0 + te); 305 de = PetscPowScalar(ae,beta); 306 fe = de*(te - t0); 307 308 fs = zero; 309 310 tn = x[k][j+1][i]; 311 an = 0.5*(t0 + tn); 312 dn = PetscPowScalar(an,beta); 313 fn = dn*(tn - t0); 314 315 if (k > 0) { 316 td = x[k-1][j][i]; 317 ad = 0.5*(t0 + td); 318 dd = PetscPowScalar(ad,beta); 319 fd = dd*(t0 - td); 320 } else { 321 fd = zero; 322 } 323 324 if (k < mz-1) { 325 tu = x[k+1][j][i]; 326 au = 0.5*(t0 + tu); 327 du = PetscPowScalar(au,beta); 328 fu = du*(tu - t0); 329 } else { 330 fu = zero; 331 } 332 333 } else if (j == my-1) { 334 335 /* top (north) boundary, and i <> 0 or mx-1 */ 336 tw = x[k][j][i-1]; 337 aw = 0.5*(t0 + tw); 338 dw = PetscPowScalar(aw,beta); 339 fw = dw*(t0 - tw); 340 341 te = x[k][j][i+1]; 342 ae = 0.5*(t0 + te); 343 de = PetscPowScalar(ae,beta); 344 fe = de*(te - t0); 345 346 ts = x[k][j-1][i]; 347 as = 0.5*(t0 + ts); 348 ds = PetscPowScalar(as,beta); 349 fs = ds*(t0 - ts); 350 351 fn = zero; 352 353 if (k > 0) { 354 td = x[k-1][j][i]; 355 ad = 0.5*(t0 + td); 356 dd = PetscPowScalar(ad,beta); 357 fd = dd*(t0 - td); 358 } else { 359 fd = zero; 360 } 361 362 if (k < mz-1) { 363 tu = x[k+1][j][i]; 364 au = 0.5*(t0 + tu); 365 du = PetscPowScalar(au,beta); 366 fu = du*(tu - t0); 367 } else { 368 fu = zero; 369 } 370 371 } else if (k == 0) { 372 373 /* down boundary (interior only) */ 374 tw = x[k][j][i-1]; 375 aw = 0.5*(t0 + tw); 376 dw = PetscPowScalar(aw,beta); 377 fw = dw*(t0 - tw); 378 379 te = x[k][j][i+1]; 380 ae = 0.5*(t0 + te); 381 de = PetscPowScalar(ae,beta); 382 fe = de*(te - t0); 383 384 ts = x[k][j-1][i]; 385 as = 0.5*(t0 + ts); 386 ds = PetscPowScalar(as,beta); 387 fs = ds*(t0 - ts); 388 389 tn = x[k][j+1][i]; 390 an = 0.5*(t0 + tn); 391 dn = PetscPowScalar(an,beta); 392 fn = dn*(tn - t0); 393 394 fd = zero; 395 396 tu = x[k+1][j][i]; 397 au = 0.5*(t0 + tu); 398 du = PetscPowScalar(au,beta); 399 fu = du*(tu - t0); 400 401 } else if (k == mz-1) { 402 403 /* up boundary (interior only) */ 404 tw = x[k][j][i-1]; 405 aw = 0.5*(t0 + tw); 406 dw = PetscPowScalar(aw,beta); 407 fw = dw*(t0 - tw); 408 409 te = x[k][j][i+1]; 410 ae = 0.5*(t0 + te); 411 de = PetscPowScalar(ae,beta); 412 fe = de*(te - t0); 413 414 ts = x[k][j-1][i]; 415 as = 0.5*(t0 + ts); 416 ds = PetscPowScalar(as,beta); 417 fs = ds*(t0 - ts); 418 419 tn = x[k][j+1][i]; 420 an = 0.5*(t0 + tn); 421 dn = PetscPowScalar(an,beta); 422 fn = dn*(tn - t0); 423 424 td = x[k-1][j][i]; 425 ad = 0.5*(t0 + td); 426 dd = PetscPowScalar(ad,beta); 427 fd = dd*(t0 - td); 428 429 fu = zero; 430 } 431 432 f[k][j][i] = -hyhzdhx*(fe-fw) - hzhxdhy*(fn-fs) - hxhydhz*(fu-fd); 433 } 434 } 435 } 436 PetscCall(DMDAVecRestoreArray(da,localX,&x)); 437 PetscCall(DMDAVecRestoreArray(da,F,&f)); 438 PetscCall(DMRestoreLocalVector(da,&localX)); 439 PetscCall(PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm)); 440 PetscFunctionReturn(0); 441 } 442 /* -------------------- Evaluate Jacobian F(x) --------------------- */ 443 PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr) 444 { 445 AppCtx *user = (AppCtx*)ptr; 446 PetscInt i,j,k,mx,my,mz,xs,ys,zs,xm,ym,zm; 447 PetscScalar one = 1.0; 448 PetscScalar hx,hy,hz,hxhydhz,hyhzdhx,hzhxdhy; 449 PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw; 450 PetscScalar tleft,tright,beta,td,ad,dd,tu,au,du,v[7],bm1,coef; 451 PetscScalar ***x,bn,bs,be,bw,bu,bd,gn,gs,ge,gw,gu,gd; 452 Vec localX; 453 MatStencil c[7],row; 454 DM da; 455 456 PetscFunctionBeginUser; 457 PetscCall(SNESGetDM(snes,&da)); 458 PetscCall(DMGetLocalVector(da,&localX)); 459 PetscCall(DMDAGetInfo(da,NULL,&mx,&my,&mz,0,0,0,0,0,0,0,0,0)); 460 hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); hz = one/(PetscReal)(mz-1); 461 hxhydhz = hx*hy/hz; hyhzdhx = hy*hz/hx; hzhxdhy = hz*hx/hy; 462 tleft = user->tleft; tright = user->tright; 463 beta = user->beta; bm1 = user->bm1; coef = user->coef; 464 465 /* Get ghost points */ 466 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 467 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 468 PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm)); 469 PetscCall(DMDAVecGetArray(da,localX,&x)); 470 471 /* Evaluate Jacobian of function */ 472 for (k=zs; k<zs+zm; k++) { 473 for (j=ys; j<ys+ym; j++) { 474 for (i=xs; i<xs+xm; i++) { 475 t0 = x[k][j][i]; 476 row.k = k; row.j = j; row.i = i; 477 if (i > 0 && i < mx-1 && j > 0 && j < my-1 && k > 0 && k < mz-1) { 478 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; c[0].j = j; c[0].i = i; v[0] = -hxhydhz*(dd - gd); 524 c[1].k = k; c[1].j = j-1; c[1].i = i; 525 v[1] = -hzhxdhy*(ds - gs); 526 c[2].k = k; c[2].j = j; c[2].i = i-1; 527 v[2] = -hyhzdhx*(dw - gw); 528 c[3].k = k; c[3].j = j; c[3].i = i; 529 v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 530 c[4].k = k; c[4].j = j; c[4].i = i+1; 531 v[4] = -hyhzdhx*(de + ge); 532 c[5].k = k; c[5].j = j+1; c[5].i = i; 533 v[5] = -hzhxdhy*(dn + gn); 534 c[6].k = k+1; c[6].j = j; c[6].i = i; 535 v[6] = -hxhydhz*(du + gu); 536 PetscCall(MatSetValuesStencil(jac,1,&row,7,c,v,INSERT_VALUES)); 537 538 } else if (i == 0) { 539 540 /* left-hand plane boundary */ 541 tw = tleft; 542 aw = 0.5*(t0 + tw); 543 bw = PetscPowScalar(aw,bm1); 544 /* dw = bw * aw */ 545 dw = PetscPowScalar(aw,beta); 546 gw = coef*bw*(t0 - tw); 547 548 te = x[k][j][i+1]; 549 ae = 0.5*(t0 + te); 550 be = PetscPowScalar(ae,bm1); 551 /* de = be * ae; */ 552 de = PetscPowScalar(ae,beta); 553 ge = coef*be*(te - t0); 554 555 /* left-hand bottom edge */ 556 if (j == 0) { 557 558 tn = x[k][j+1][i]; 559 an = 0.5*(t0 + tn); 560 bn = PetscPowScalar(an,bm1); 561 /* dn = bn * an; */ 562 dn = PetscPowScalar(an,beta); 563 gn = coef*bn*(tn - t0); 564 565 /* left-hand bottom down corner */ 566 if (k == 0) { 567 568 tu = x[k+1][j][i]; 569 au = 0.5*(t0 + tu); 570 bu = PetscPowScalar(au,bm1); 571 /* du = bu * au; */ 572 du = PetscPowScalar(au,beta); 573 gu = coef*bu*(tu - t0); 574 575 c[0].k = k; c[0].j = j; c[0].i = i; 576 v[0] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 577 c[1].k = k; c[1].j = j; c[1].i = i+1; 578 v[1] = -hyhzdhx*(de + ge); 579 c[2].k = k; c[2].j = j+1; c[2].i = i; 580 v[2] = -hzhxdhy*(dn + gn); 581 c[3].k = k+1; c[3].j = j; c[3].i = i; 582 v[3] = -hxhydhz*(du + gu); 583 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 584 585 /* left-hand bottom interior edge */ 586 } else if (k < mz-1) { 587 588 tu = x[k+1][j][i]; 589 au = 0.5*(t0 + tu); 590 bu = PetscPowScalar(au,bm1); 591 /* du = bu * au; */ 592 du = PetscPowScalar(au,beta); 593 gu = coef*bu*(tu - t0); 594 595 td = x[k-1][j][i]; 596 ad = 0.5*(t0 + td); 597 bd = PetscPowScalar(ad,bm1); 598 /* dd = bd * ad; */ 599 dd = PetscPowScalar(ad,beta); 600 gd = coef*bd*(td - t0); 601 602 c[0].k = k-1; c[0].j = j; c[0].i = i; 603 v[0] = -hxhydhz*(dd - gd); 604 c[1].k = k; c[1].j = j; c[1].i = i; 605 v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 606 c[2].k = k; c[2].j = j; c[2].i = i+1; 607 v[2] = -hyhzdhx*(de + ge); 608 c[3].k = k; c[3].j = j+1; c[3].i = i; 609 v[3] = -hzhxdhy*(dn + gn); 610 c[4].k = k+1; c[4].j = j; c[4].i = i; 611 v[4] = -hxhydhz*(du + gu); 612 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 613 614 /* left-hand bottom up corner */ 615 } else { 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; c[0].j = j; c[0].i = i; 625 v[0] = -hxhydhz*(dd - gd); 626 c[1].k = k; c[1].j = j; c[1].i = i; 627 v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 628 c[2].k = k; c[2].j = j; c[2].i = i+1; 629 v[2] = -hyhzdhx*(de + ge); 630 c[3].k = k; c[3].j = j+1; c[3].i = i; 631 v[3] = -hzhxdhy*(dn + gn); 632 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 633 } 634 635 /* left-hand top edge */ 636 } else if (j == my-1) { 637 638 ts = x[k][j-1][i]; 639 as = 0.5*(t0 + ts); 640 bs = PetscPowScalar(as,bm1); 641 /* ds = bs * as; */ 642 ds = PetscPowScalar(as,beta); 643 gs = coef*bs*(ts - t0); 644 645 /* left-hand top down corner */ 646 if (k == 0) { 647 648 tu = x[k+1][j][i]; 649 au = 0.5*(t0 + tu); 650 bu = PetscPowScalar(au,bm1); 651 /* du = bu * au; */ 652 du = PetscPowScalar(au,beta); 653 gu = coef*bu*(tu - t0); 654 655 c[0].k = k; c[0].j = j-1; c[0].i = i; 656 v[0] = -hzhxdhy*(ds - gs); 657 c[1].k = k; c[1].j = j; c[1].i = i; 658 v[1] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 659 c[2].k = k; c[2].j = j; c[2].i = i+1; 660 v[2] = -hyhzdhx*(de + ge); 661 c[3].k = k+1; c[3].j = j; c[3].i = i; 662 v[3] = -hxhydhz*(du + gu); 663 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 664 665 /* left-hand top interior edge */ 666 } else if (k < mz-1) { 667 668 tu = x[k+1][j][i]; 669 au = 0.5*(t0 + tu); 670 bu = PetscPowScalar(au,bm1); 671 /* du = bu * au; */ 672 du = PetscPowScalar(au,beta); 673 gu = coef*bu*(tu - t0); 674 675 td = x[k-1][j][i]; 676 ad = 0.5*(t0 + td); 677 bd = PetscPowScalar(ad,bm1); 678 /* dd = bd * ad; */ 679 dd = PetscPowScalar(ad,beta); 680 gd = coef*bd*(td - t0); 681 682 c[0].k = k-1; c[0].j = j; c[0].i = i; 683 v[0] = -hxhydhz*(dd - gd); 684 c[1].k = k; c[1].j = j-1; c[1].i = i; 685 v[1] = -hzhxdhy*(ds - gs); 686 c[2].k = k; c[2].j = j; c[2].i = i; 687 v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 688 c[3].k = k; c[3].j = j; c[3].i = i+1; 689 v[3] = -hyhzdhx*(de + ge); 690 c[4].k = k+1; c[4].j = j; c[4].i = i; 691 v[4] = -hxhydhz*(du + gu); 692 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 693 694 /* left-hand top up corner */ 695 } else { 696 697 td = x[k-1][j][i]; 698 ad = 0.5*(t0 + td); 699 bd = PetscPowScalar(ad,bm1); 700 /* dd = bd * ad; */ 701 dd = PetscPowScalar(ad,beta); 702 gd = coef*bd*(td - t0); 703 704 c[0].k = k-1; c[0].j = j; c[0].i = i; 705 v[0] = -hxhydhz*(dd - gd); 706 c[1].k = k; c[1].j = j-1; c[1].i = i; 707 v[1] = -hzhxdhy*(ds - gs); 708 c[2].k = k; c[2].j = j; c[2].i = i; 709 v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 710 c[3].k = k; c[3].j = j; c[3].i = i+1; 711 v[3] = -hyhzdhx*(de + ge); 712 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 713 } 714 715 } else { 716 717 ts = x[k][j-1][i]; 718 as = 0.5*(t0 + ts); 719 bs = PetscPowScalar(as,bm1); 720 /* ds = bs * as; */ 721 ds = PetscPowScalar(as,beta); 722 gs = coef*bs*(t0 - ts); 723 724 tn = x[k][j+1][i]; 725 an = 0.5*(t0 + tn); 726 bn = PetscPowScalar(an,bm1); 727 /* dn = bn * an; */ 728 dn = PetscPowScalar(an,beta); 729 gn = coef*bn*(tn - t0); 730 731 /* left-hand down interior edge */ 732 if (k == 0) { 733 734 tu = x[k+1][j][i]; 735 au = 0.5*(t0 + tu); 736 bu = PetscPowScalar(au,bm1); 737 /* du = bu * au; */ 738 du = PetscPowScalar(au,beta); 739 gu = coef*bu*(tu - t0); 740 741 c[0].k = k; c[0].j = j-1; c[0].i = i; 742 v[0] = -hzhxdhy*(ds - gs); 743 c[1].k = k; c[1].j = j; c[1].i = i; 744 v[1] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 745 c[2].k = k; c[2].j = j; c[2].i = i+1; 746 v[2] = -hyhzdhx*(de + ge); 747 c[3].k = k; c[3].j = j+1; c[3].i = i; 748 v[3] = -hzhxdhy*(dn + gn); 749 c[4].k = k+1; c[4].j = j; c[4].i = i; 750 v[4] = -hxhydhz*(du + gu); 751 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 752 753 } else if (k == mz-1) { /* left-hand up interior edge */ 754 755 td = x[k-1][j][i]; 756 ad = 0.5*(t0 + td); 757 bd = PetscPowScalar(ad,bm1); 758 /* dd = bd * ad; */ 759 dd = PetscPowScalar(ad,beta); 760 gd = coef*bd*(t0 - td); 761 762 c[0].k = k-1; c[0].j = j; c[0].i = i; 763 v[0] = -hxhydhz*(dd - gd); 764 c[1].k = k; c[1].j = j-1; c[1].i = i; 765 v[1] = -hzhxdhy*(ds - gs); 766 c[2].k = k; c[2].j = j; c[2].i = i; 767 v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 768 c[3].k = k; c[3].j = j; c[3].i = i+1; 769 v[3] = -hyhzdhx*(de + ge); 770 c[4].k = k; c[4].j = j+1; c[4].i = i; 771 v[4] = -hzhxdhy*(dn + gn); 772 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 773 } else { /* left-hand interior plane */ 774 775 td = x[k-1][j][i]; 776 ad = 0.5*(t0 + td); 777 bd = PetscPowScalar(ad,bm1); 778 /* dd = bd * ad; */ 779 dd = PetscPowScalar(ad,beta); 780 gd = coef*bd*(t0 - td); 781 782 tu = x[k+1][j][i]; 783 au = 0.5*(t0 + tu); 784 bu = PetscPowScalar(au,bm1); 785 /* du = bu * au; */ 786 du = PetscPowScalar(au,beta); 787 gu = coef*bu*(tu - t0); 788 789 c[0].k = k-1; c[0].j = j; c[0].i = i; 790 v[0] = -hxhydhz*(dd - gd); 791 c[1].k = k; c[1].j = j-1; c[1].i = i; 792 v[1] = -hzhxdhy*(ds - gs); 793 c[2].k = k; c[2].j = j; c[2].i = i; 794 v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 795 c[3].k = k; c[3].j = j; c[3].i = i+1; 796 v[3] = -hyhzdhx*(de + ge); 797 c[4].k = k; c[4].j = j+1; c[4].i = i; 798 v[4] = -hzhxdhy*(dn + gn); 799 c[5].k = k+1; c[5].j = j; c[5].i = i; 800 v[5] = -hxhydhz*(du + gu); 801 PetscCall(MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES)); 802 } 803 } 804 805 } else if (i == mx-1) { 806 807 /* right-hand plane boundary */ 808 tw = x[k][j][i-1]; 809 aw = 0.5*(t0 + tw); 810 bw = PetscPowScalar(aw,bm1); 811 /* dw = bw * aw */ 812 dw = PetscPowScalar(aw,beta); 813 gw = coef*bw*(t0 - tw); 814 815 te = tright; 816 ae = 0.5*(t0 + te); 817 be = PetscPowScalar(ae,bm1); 818 /* de = be * ae; */ 819 de = PetscPowScalar(ae,beta); 820 ge = coef*be*(te - t0); 821 822 /* right-hand bottom edge */ 823 if (j == 0) { 824 825 tn = x[k][j+1][i]; 826 an = 0.5*(t0 + tn); 827 bn = PetscPowScalar(an,bm1); 828 /* dn = bn * an; */ 829 dn = PetscPowScalar(an,beta); 830 gn = coef*bn*(tn - t0); 831 832 /* right-hand bottom down corner */ 833 if (k == 0) { 834 835 tu = x[k+1][j][i]; 836 au = 0.5*(t0 + tu); 837 bu = PetscPowScalar(au,bm1); 838 /* du = bu * au; */ 839 du = PetscPowScalar(au,beta); 840 gu = coef*bu*(tu - t0); 841 842 c[0].k = k; c[0].j = j; c[0].i = i-1; 843 v[0] = -hyhzdhx*(dw - gw); 844 c[1].k = k; c[1].j = j; c[1].i = i; 845 v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 846 c[2].k = k; c[2].j = j+1; c[2].i = i; 847 v[2] = -hzhxdhy*(dn + gn); 848 c[3].k = k+1; c[3].j = j; c[3].i = i; 849 v[3] = -hxhydhz*(du + gu); 850 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 851 852 /* right-hand bottom interior edge */ 853 } else if (k < mz-1) { 854 855 tu = x[k+1][j][i]; 856 au = 0.5*(t0 + tu); 857 bu = PetscPowScalar(au,bm1); 858 /* du = bu * au; */ 859 du = PetscPowScalar(au,beta); 860 gu = coef*bu*(tu - t0); 861 862 td = x[k-1][j][i]; 863 ad = 0.5*(t0 + td); 864 bd = PetscPowScalar(ad,bm1); 865 /* dd = bd * ad; */ 866 dd = PetscPowScalar(ad,beta); 867 gd = coef*bd*(td - t0); 868 869 c[0].k = k-1; c[0].j = j; c[0].i = i; 870 v[0] = -hxhydhz*(dd - gd); 871 c[1].k = k; c[1].j = j; c[1].i = i-1; 872 v[1] = -hyhzdhx*(dw - gw); 873 c[2].k = k; c[2].j = j; c[2].i = i; 874 v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 875 c[3].k = k; c[3].j = j+1; c[3].i = i; 876 v[3] = -hzhxdhy*(dn + gn); 877 c[4].k = k+1; c[4].j = j; c[4].i = i; 878 v[4] = -hxhydhz*(du + gu); 879 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 880 881 /* right-hand bottom up corner */ 882 } else { 883 884 td = x[k-1][j][i]; 885 ad = 0.5*(t0 + td); 886 bd = PetscPowScalar(ad,bm1); 887 /* dd = bd * ad; */ 888 dd = PetscPowScalar(ad,beta); 889 gd = coef*bd*(td - t0); 890 891 c[0].k = k-1; c[0].j = j; c[0].i = i; 892 v[0] = -hxhydhz*(dd - gd); 893 c[1].k = k; c[1].j = j; c[1].i = i-1; 894 v[1] = -hyhzdhx*(dw - gw); 895 c[2].k = k; c[2].j = j; c[2].i = i; 896 v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 897 c[3].k = k; c[3].j = j+1; c[3].i = i; 898 v[3] = -hzhxdhy*(dn + gn); 899 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 900 } 901 902 /* right-hand top edge */ 903 } else if (j == my-1) { 904 905 ts = x[k][j-1][i]; 906 as = 0.5*(t0 + ts); 907 bs = PetscPowScalar(as,bm1); 908 /* ds = bs * as; */ 909 ds = PetscPowScalar(as,beta); 910 gs = coef*bs*(ts - t0); 911 912 /* right-hand top down corner */ 913 if (k == 0) { 914 915 tu = x[k+1][j][i]; 916 au = 0.5*(t0 + tu); 917 bu = PetscPowScalar(au,bm1); 918 /* du = bu * au; */ 919 du = PetscPowScalar(au,beta); 920 gu = coef*bu*(tu - t0); 921 922 c[0].k = k; c[0].j = j-1; c[0].i = i; 923 v[0] = -hzhxdhy*(ds - gs); 924 c[1].k = k; c[1].j = j; c[1].i = i-1; 925 v[1] = -hyhzdhx*(dw - gw); 926 c[2].k = k; c[2].j = j; c[2].i = i; 927 v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 928 c[3].k = k+1; c[3].j = j; c[3].i = i; 929 v[3] = -hxhydhz*(du + gu); 930 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 931 932 /* right-hand top interior edge */ 933 } else if (k < mz-1) { 934 935 tu = x[k+1][j][i]; 936 au = 0.5*(t0 + tu); 937 bu = PetscPowScalar(au,bm1); 938 /* du = bu * au; */ 939 du = PetscPowScalar(au,beta); 940 gu = coef*bu*(tu - t0); 941 942 td = x[k-1][j][i]; 943 ad = 0.5*(t0 + td); 944 bd = PetscPowScalar(ad,bm1); 945 /* dd = bd * ad; */ 946 dd = PetscPowScalar(ad,beta); 947 gd = coef*bd*(td - t0); 948 949 c[0].k = k-1; c[0].j = j; c[0].i = i; 950 v[0] = -hxhydhz*(dd - gd); 951 c[1].k = k; c[1].j = j-1; c[1].i = i; 952 v[1] = -hzhxdhy*(ds - gs); 953 c[2].k = k; c[2].j = j; c[2].i = i-1; 954 v[2] = -hyhzdhx*(dw - gw); 955 c[3].k = k; c[3].j = j; c[3].i = i; 956 v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 957 c[4].k = k+1; c[4].j = j; c[4].i = i; 958 v[4] = -hxhydhz*(du + gu); 959 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 960 961 /* right-hand top up corner */ 962 } else { 963 964 td = x[k-1][j][i]; 965 ad = 0.5*(t0 + td); 966 bd = PetscPowScalar(ad,bm1); 967 /* dd = bd * ad; */ 968 dd = PetscPowScalar(ad,beta); 969 gd = coef*bd*(td - t0); 970 971 c[0].k = k-1; c[0].j = j; c[0].i = i; 972 v[0] = -hxhydhz*(dd - gd); 973 c[1].k = k; c[1].j = j-1; c[1].i = i; 974 v[1] = -hzhxdhy*(ds - gs); 975 c[2].k = k; c[2].j = j; c[2].i = i-1; 976 v[2] = -hyhzdhx*(dw - gw); 977 c[3].k = k; c[3].j = j; c[3].i = i; 978 v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 979 PetscCall(MatSetValuesStencil(jac,1,&row,4,c,v,INSERT_VALUES)); 980 } 981 982 } else { 983 984 ts = x[k][j-1][i]; 985 as = 0.5*(t0 + ts); 986 bs = PetscPowScalar(as,bm1); 987 /* ds = bs * as; */ 988 ds = PetscPowScalar(as,beta); 989 gs = coef*bs*(t0 - ts); 990 991 tn = x[k][j+1][i]; 992 an = 0.5*(t0 + tn); 993 bn = PetscPowScalar(an,bm1); 994 /* dn = bn * an; */ 995 dn = PetscPowScalar(an,beta); 996 gn = coef*bn*(tn - t0); 997 998 /* right-hand down interior edge */ 999 if (k == 0) { 1000 1001 tu = x[k+1][j][i]; 1002 au = 0.5*(t0 + tu); 1003 bu = PetscPowScalar(au,bm1); 1004 /* du = bu * au; */ 1005 du = PetscPowScalar(au,beta); 1006 gu = coef*bu*(tu - t0); 1007 1008 c[0].k = k; c[0].j = j-1; c[0].i = i; 1009 v[0] = -hzhxdhy*(ds - gs); 1010 c[1].k = k; c[1].j = j; c[1].i = i-1; 1011 v[1] = -hyhzdhx*(dw - gw); 1012 c[2].k = k; c[2].j = j; c[2].i = i; 1013 v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 1014 c[3].k = k; c[3].j = j+1; c[3].i = i; 1015 v[3] = -hzhxdhy*(dn + gn); 1016 c[4].k = k+1; c[4].j = j; c[4].i = i; 1017 v[4] = -hxhydhz*(du + gu); 1018 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 1019 1020 } else if (k == mz-1) { /* right-hand up interior edge */ 1021 1022 td = x[k-1][j][i]; 1023 ad = 0.5*(t0 + td); 1024 bd = PetscPowScalar(ad,bm1); 1025 /* dd = bd * ad; */ 1026 dd = PetscPowScalar(ad,beta); 1027 gd = coef*bd*(t0 - td); 1028 1029 c[0].k = k-1; c[0].j = j; c[0].i = i; 1030 v[0] = -hxhydhz*(dd - gd); 1031 c[1].k = k; c[1].j = j-1; c[1].i = i; 1032 v[1] = -hzhxdhy*(ds - gs); 1033 c[2].k = k; c[2].j = j; c[2].i = i-1; 1034 v[2] = -hyhzdhx*(dw - gw); 1035 c[3].k = k; c[3].j = j; c[3].i = i; 1036 v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 1037 c[4].k = k; c[4].j = j+1; c[4].i = i; 1038 v[4] = -hzhxdhy*(dn + gn); 1039 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 1040 1041 } else { /* right-hand interior plane */ 1042 1043 td = x[k-1][j][i]; 1044 ad = 0.5*(t0 + td); 1045 bd = PetscPowScalar(ad,bm1); 1046 /* dd = bd * ad; */ 1047 dd = PetscPowScalar(ad,beta); 1048 gd = coef*bd*(t0 - td); 1049 1050 tu = x[k+1][j][i]; 1051 au = 0.5*(t0 + tu); 1052 bu = PetscPowScalar(au,bm1); 1053 /* du = bu * au; */ 1054 du = PetscPowScalar(au,beta); 1055 gu = coef*bu*(tu - t0); 1056 1057 c[0].k = k-1; c[0].j = j; c[0].i = i; 1058 v[0] = -hxhydhz*(dd - gd); 1059 c[1].k = k; c[1].j = j-1; c[1].i = i; 1060 v[1] = -hzhxdhy*(ds - gs); 1061 c[2].k = k; c[2].j = j; c[2].i = i-1; 1062 v[2] = -hyhzdhx*(dw - gw); 1063 c[3].k = k; c[3].j = j; c[3].i = i; 1064 v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 1065 c[4].k = k; c[4].j = j+1; c[4].i = i; 1066 v[4] = -hzhxdhy*(dn + gn); 1067 c[5].k = k+1; c[5].j = j; c[5].i = i; 1068 v[5] = -hxhydhz*(du + gu); 1069 PetscCall(MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES)); 1070 } 1071 } 1072 1073 } else if (j == 0) { 1074 1075 tw = x[k][j][i-1]; 1076 aw = 0.5*(t0 + tw); 1077 bw = PetscPowScalar(aw,bm1); 1078 /* dw = bw * aw */ 1079 dw = PetscPowScalar(aw,beta); 1080 gw = coef*bw*(t0 - tw); 1081 1082 te = x[k][j][i+1]; 1083 ae = 0.5*(t0 + te); 1084 be = PetscPowScalar(ae,bm1); 1085 /* de = be * ae; */ 1086 de = PetscPowScalar(ae,beta); 1087 ge = coef*be*(te - t0); 1088 1089 tn = x[k][j+1][i]; 1090 an = 0.5*(t0 + tn); 1091 bn = PetscPowScalar(an,bm1); 1092 /* dn = bn * an; */ 1093 dn = PetscPowScalar(an,beta); 1094 gn = coef*bn*(tn - t0); 1095 1096 /* bottom down interior edge */ 1097 if (k == 0) { 1098 1099 tu = x[k+1][j][i]; 1100 au = 0.5*(t0 + tu); 1101 bu = PetscPowScalar(au,bm1); 1102 /* du = bu * au; */ 1103 du = PetscPowScalar(au,beta); 1104 gu = coef*bu*(tu - t0); 1105 1106 c[0].k = k; c[0].j = j; c[0].i = i-1; 1107 v[0] = -hyhzdhx*(dw - gw); 1108 c[1].k = k; c[1].j = j; c[1].i = i; 1109 v[1] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 1110 c[2].k = k; c[2].j = j; c[2].i = i+1; 1111 v[2] = -hyhzdhx*(de + ge); 1112 c[3].k = k; c[3].j = j+1; c[3].i = i; 1113 v[3] = -hzhxdhy*(dn + gn); 1114 c[4].k = k+1; c[4].j = j; c[4].i = i; 1115 v[4] = -hxhydhz*(du + gu); 1116 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 1117 1118 } else if (k == mz-1) { /* bottom up interior edge */ 1119 1120 td = x[k-1][j][i]; 1121 ad = 0.5*(t0 + td); 1122 bd = PetscPowScalar(ad,bm1); 1123 /* dd = bd * ad; */ 1124 dd = PetscPowScalar(ad,beta); 1125 gd = coef*bd*(td - t0); 1126 1127 c[0].k = k-1; c[0].j = j; c[0].i = i; 1128 v[0] = -hxhydhz*(dd - gd); 1129 c[1].k = k; c[1].j = j; c[1].i = i-1; 1130 v[1] = -hyhzdhx*(dw - gw); 1131 c[2].k = k; c[2].j = j; c[2].i = i; 1132 v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 1133 c[3].k = k; c[3].j = j; c[3].i = i+1; 1134 v[3] = -hyhzdhx*(de + ge); 1135 c[4].k = k; c[4].j = j+1; c[4].i = i; 1136 v[4] = -hzhxdhy*(dn + gn); 1137 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 1138 1139 } else { /* bottom interior plane */ 1140 1141 tu = x[k+1][j][i]; 1142 au = 0.5*(t0 + tu); 1143 bu = PetscPowScalar(au,bm1); 1144 /* du = bu * au; */ 1145 du = PetscPowScalar(au,beta); 1146 gu = coef*bu*(tu - t0); 1147 1148 td = x[k-1][j][i]; 1149 ad = 0.5*(t0 + td); 1150 bd = PetscPowScalar(ad,bm1); 1151 /* dd = bd * ad; */ 1152 dd = PetscPowScalar(ad,beta); 1153 gd = coef*bd*(td - t0); 1154 1155 c[0].k = k-1; c[0].j = j; c[0].i = i; 1156 v[0] = -hxhydhz*(dd - gd); 1157 c[1].k = k; c[1].j = j; c[1].i = i-1; 1158 v[1] = -hyhzdhx*(dw - gw); 1159 c[2].k = k; c[2].j = j; c[2].i = i; 1160 v[2] = hzhxdhy*(dn - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 1161 c[3].k = k; c[3].j = j; c[3].i = i+1; 1162 v[3] = -hyhzdhx*(de + ge); 1163 c[4].k = k; c[4].j = j+1; c[4].i = i; 1164 v[4] = -hzhxdhy*(dn + gn); 1165 c[5].k = k+1; c[5].j = j; c[5].i = i; 1166 v[5] = -hxhydhz*(du + gu); 1167 PetscCall(MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES)); 1168 } 1169 1170 } else if (j == my-1) { 1171 1172 tw = x[k][j][i-1]; 1173 aw = 0.5*(t0 + tw); 1174 bw = PetscPowScalar(aw,bm1); 1175 /* dw = bw * aw */ 1176 dw = PetscPowScalar(aw,beta); 1177 gw = coef*bw*(t0 - tw); 1178 1179 te = x[k][j][i+1]; 1180 ae = 0.5*(t0 + te); 1181 be = PetscPowScalar(ae,bm1); 1182 /* de = be * ae; */ 1183 de = PetscPowScalar(ae,beta); 1184 ge = coef*be*(te - t0); 1185 1186 ts = x[k][j-1][i]; 1187 as = 0.5*(t0 + ts); 1188 bs = PetscPowScalar(as,bm1); 1189 /* ds = bs * as; */ 1190 ds = PetscPowScalar(as,beta); 1191 gs = coef*bs*(t0 - ts); 1192 1193 /* top down interior edge */ 1194 if (k == 0) { 1195 1196 tu = x[k+1][j][i]; 1197 au = 0.5*(t0 + tu); 1198 bu = PetscPowScalar(au,bm1); 1199 /* du = bu * au; */ 1200 du = PetscPowScalar(au,beta); 1201 gu = coef*bu*(tu - t0); 1202 1203 c[0].k = k; c[0].j = j-1; c[0].i = i; 1204 v[0] = -hzhxdhy*(ds - gs); 1205 c[1].k = k; c[1].j = j; c[1].i = i-1; 1206 v[1] = -hyhzdhx*(dw - gw); 1207 c[2].k = k; c[2].j = j; c[2].i = i; 1208 v[2] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 1209 c[3].k = k; c[3].j = j; c[3].i = i+1; 1210 v[3] = -hyhzdhx*(de + ge); 1211 c[4].k = k+1; c[4].j = j; c[4].i = i; 1212 v[4] = -hxhydhz*(du + gu); 1213 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 1214 1215 } else if (k == mz-1) { /* top up interior edge */ 1216 1217 td = x[k-1][j][i]; 1218 ad = 0.5*(t0 + td); 1219 bd = PetscPowScalar(ad,bm1); 1220 /* dd = bd * ad; */ 1221 dd = PetscPowScalar(ad,beta); 1222 gd = coef*bd*(td - t0); 1223 1224 c[0].k = k-1; c[0].j = j; c[0].i = i; 1225 v[0] = -hxhydhz*(dd - gd); 1226 c[1].k = k; c[1].j = j-1; c[1].i = i; 1227 v[1] = -hzhxdhy*(ds - gs); 1228 c[2].k = k; c[2].j = j; c[2].i = i-1; 1229 v[2] = -hyhzdhx*(dw - gw); 1230 c[3].k = k; c[3].j = j; c[3].i = i; 1231 v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 1232 c[4].k = k; c[4].j = j; c[4].i = i+1; 1233 v[4] = -hyhzdhx*(de + ge); 1234 PetscCall(MatSetValuesStencil(jac,1,&row,5,c,v,INSERT_VALUES)); 1235 1236 } else { /* top interior plane */ 1237 1238 tu = x[k+1][j][i]; 1239 au = 0.5*(t0 + tu); 1240 bu = PetscPowScalar(au,bm1); 1241 /* du = bu * au; */ 1242 du = PetscPowScalar(au,beta); 1243 gu = coef*bu*(tu - t0); 1244 1245 td = x[k-1][j][i]; 1246 ad = 0.5*(t0 + td); 1247 bd = PetscPowScalar(ad,bm1); 1248 /* dd = bd * ad; */ 1249 dd = PetscPowScalar(ad,beta); 1250 gd = coef*bd*(td - t0); 1251 1252 c[0].k = k-1; c[0].j = j; c[0].i = i; 1253 v[0] = -hxhydhz*(dd - gd); 1254 c[1].k = k; c[1].j = j-1; c[1].i = i; 1255 v[1] = -hzhxdhy*(ds - gs); 1256 c[2].k = k; c[2].j = j; c[2].i = i-1; 1257 v[2] = -hyhzdhx*(dw - gw); 1258 c[3].k = k; c[3].j = j; c[3].i = i; 1259 v[3] = hzhxdhy*(ds + gs) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + du + gd - gu); 1260 c[4].k = k; c[4].j = j; c[4].i = i+1; 1261 v[4] = -hyhzdhx*(de + ge); 1262 c[5].k = k+1; c[5].j = j; c[5].i = i; 1263 v[5] = -hxhydhz*(du + gu); 1264 PetscCall(MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES)); 1265 } 1266 1267 } else if (k == 0) { 1268 1269 /* down interior plane */ 1270 1271 tw = x[k][j][i-1]; 1272 aw = 0.5*(t0 + tw); 1273 bw = PetscPowScalar(aw,bm1); 1274 /* dw = bw * aw */ 1275 dw = PetscPowScalar(aw,beta); 1276 gw = coef*bw*(t0 - tw); 1277 1278 te = x[k][j][i+1]; 1279 ae = 0.5*(t0 + te); 1280 be = PetscPowScalar(ae,bm1); 1281 /* de = be * ae; */ 1282 de = PetscPowScalar(ae,beta); 1283 ge = coef*be*(te - t0); 1284 1285 ts = x[k][j-1][i]; 1286 as = 0.5*(t0 + ts); 1287 bs = PetscPowScalar(as,bm1); 1288 /* ds = bs * as; */ 1289 ds = PetscPowScalar(as,beta); 1290 gs = coef*bs*(t0 - ts); 1291 1292 tn = x[k][j+1][i]; 1293 an = 0.5*(t0 + tn); 1294 bn = PetscPowScalar(an,bm1); 1295 /* dn = bn * an; */ 1296 dn = PetscPowScalar(an,beta); 1297 gn = coef*bn*(tn - t0); 1298 1299 tu = x[k+1][j][i]; 1300 au = 0.5*(t0 + tu); 1301 bu = PetscPowScalar(au,bm1); 1302 /* du = bu * au; */ 1303 du = PetscPowScalar(au,beta); 1304 gu = coef*bu*(tu - t0); 1305 1306 c[0].k = k; c[0].j = j-1; c[0].i = i; 1307 v[0] = -hzhxdhy*(ds - gs); 1308 c[1].k = k; c[1].j = j; c[1].i = i-1; 1309 v[1] = -hyhzdhx*(dw - gw); 1310 c[2].k = k; c[2].j = j; c[2].i = i; 1311 v[2] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(du - gu); 1312 c[3].k = k; c[3].j = j; c[3].i = i+1; 1313 v[3] = -hyhzdhx*(de + ge); 1314 c[4].k = k; c[4].j = j+1; c[4].i = i; 1315 v[4] = -hzhxdhy*(dn + gn); 1316 c[5].k = k+1; c[5].j = j; c[5].i = i; 1317 v[5] = -hxhydhz*(du + gu); 1318 PetscCall(MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES)); 1319 1320 } else if (k == mz-1) { 1321 1322 /* up interior plane */ 1323 1324 tw = x[k][j][i-1]; 1325 aw = 0.5*(t0 + tw); 1326 bw = PetscPowScalar(aw,bm1); 1327 /* dw = bw * aw */ 1328 dw = PetscPowScalar(aw,beta); 1329 gw = coef*bw*(t0 - tw); 1330 1331 te = x[k][j][i+1]; 1332 ae = 0.5*(t0 + te); 1333 be = PetscPowScalar(ae,bm1); 1334 /* de = be * ae; */ 1335 de = PetscPowScalar(ae,beta); 1336 ge = coef*be*(te - t0); 1337 1338 ts = x[k][j-1][i]; 1339 as = 0.5*(t0 + ts); 1340 bs = PetscPowScalar(as,bm1); 1341 /* ds = bs * as; */ 1342 ds = PetscPowScalar(as,beta); 1343 gs = coef*bs*(t0 - ts); 1344 1345 tn = x[k][j+1][i]; 1346 an = 0.5*(t0 + tn); 1347 bn = PetscPowScalar(an,bm1); 1348 /* dn = bn * an; */ 1349 dn = PetscPowScalar(an,beta); 1350 gn = coef*bn*(tn - t0); 1351 1352 td = x[k-1][j][i]; 1353 ad = 0.5*(t0 + td); 1354 bd = PetscPowScalar(ad,bm1); 1355 /* dd = bd * ad; */ 1356 dd = PetscPowScalar(ad,beta); 1357 gd = coef*bd*(t0 - td); 1358 1359 c[0].k = k-1; c[0].j = j; c[0].i = i; 1360 v[0] = -hxhydhz*(dd - gd); 1361 c[1].k = k; c[1].j = j-1; c[1].i = i; 1362 v[1] = -hzhxdhy*(ds - gs); 1363 c[2].k = k; c[2].j = j; c[2].i = i-1; 1364 v[2] = -hyhzdhx*(dw - gw); 1365 c[3].k = k; c[3].j = j; c[3].i = i; 1366 v[3] = hzhxdhy*(ds + dn + gs - gn) + hyhzdhx*(dw + de + gw - ge) + hxhydhz*(dd + gd); 1367 c[4].k = k; c[4].j = j; c[4].i = i+1; 1368 v[4] = -hyhzdhx*(de + ge); 1369 c[5].k = k; c[5].j = j+1; c[5].i = i; 1370 v[5] = -hzhxdhy*(dn + gn); 1371 PetscCall(MatSetValuesStencil(jac,1,&row,6,c,v,INSERT_VALUES)); 1372 } 1373 } 1374 } 1375 } 1376 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 1377 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 1378 if (jac != J) { 1379 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 1380 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 1381 } 1382 PetscCall(DMDAVecRestoreArray(da,localX,&x)); 1383 PetscCall(DMRestoreLocalVector(da,&localX)); 1384 1385 PetscCall(PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym)); 1386 PetscFunctionReturn(0); 1387 } 1388 1389 /*TEST 1390 1391 test: 1392 nsize: 4 1393 args: -snes_monitor_short -pc_mg_type full -ksp_type fgmres -pc_type mg -snes_view -pc_mg_levels 2 -pc_mg_galerkin pmat 1394 requires: !single 1395 1396 TEST*/ 1397