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