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