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