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