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