1 2 static char help[] = "Nonlinear Radiative Transport PDE with multigrid in 2d.\n\ 3 Uses 2-dimensional distributed arrays.\n\ 4 A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\ 5 \n\ 6 Solves the linear systems via multilevel methods \n\ 7 \n\ 8 The command line\n\ 9 options are:\n\ 10 -tleft <tl>, where <tl> indicates the left Diriclet BC \n\ 11 -tright <tr>, where <tr> indicates the right Diriclet BC \n\ 12 -beta <beta>, where <beta> indicates the exponent in T \n\n"; 13 14 /* 15 16 This example models the partial differential equation 17 18 - Div(alpha* T^beta (GRAD T)) = 0. 19 20 where beta = 2.5 and alpha = 1.0 21 22 BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0. 23 24 in the unit square, which is uniformly discretized in each of x and 25 y in this simple encoding. The degrees of freedom are cell centered. 26 27 A finite volume approximation with the usual 5-point stencil 28 is used to discretize the boundary value problem to obtain a 29 nonlinear system of equations. 30 31 This code was contributed by David Keyes 32 33 */ 34 35 #include <petscsnes.h> 36 #include <petscdm.h> 37 #include <petscdmda.h> 38 39 /* User-defined application context */ 40 41 typedef struct { 42 PetscReal tleft, tright; /* Dirichlet boundary conditions */ 43 PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */ 44 } AppCtx; 45 46 #define POWFLOP 5 /* assume a pow() takes five flops */ 47 48 extern PetscErrorCode FormInitialGuess(SNES, Vec, void *); 49 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 50 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 51 52 int main(int argc, char **argv) { 53 SNES snes; 54 AppCtx user; 55 PetscInt its, lits; 56 PetscReal litspit; 57 DM da; 58 59 PetscFunctionBeginUser; 60 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 61 62 /* set problem parameters */ 63 user.tleft = 1.0; 64 user.tright = 0.1; 65 user.beta = 2.5; 66 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL)); 67 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL)); 68 PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL)); 69 user.bm1 = user.beta - 1.0; 70 user.coef = user.beta / 2.0; 71 72 /* 73 Create the multilevel DM data structure 74 */ 75 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 76 77 /* 78 Set the DMDA (grid structure) for the grids. 79 */ 80 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da)); 81 PetscCall(DMSetFromOptions(da)); 82 PetscCall(DMSetUp(da)); 83 PetscCall(DMSetApplicationContext(da, &user)); 84 PetscCall(SNESSetDM(snes, (DM)da)); 85 86 /* 87 Create the nonlinear solver, and tell it the functions to use 88 */ 89 PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user)); 90 PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user)); 91 PetscCall(SNESSetFromOptions(snes)); 92 PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL)); 93 94 PetscCall(SNESSolve(snes, NULL, NULL)); 95 PetscCall(SNESGetIterationNumber(snes, &its)); 96 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 97 litspit = ((PetscReal)lits) / ((PetscReal)its); 98 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its)); 99 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits)); 100 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit)); 101 102 PetscCall(DMDestroy(&da)); 103 PetscCall(SNESDestroy(&snes)); 104 PetscCall(PetscFinalize()); 105 return 0; 106 } 107 /* -------------------- Form initial approximation ----------------- */ 108 PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx) { 109 AppCtx *user; 110 PetscInt i, j, xs, ys, xm, ym; 111 PetscReal tleft; 112 PetscScalar **x; 113 DM da; 114 115 PetscFunctionBeginUser; 116 PetscCall(SNESGetDM(snes, &da)); 117 PetscCall(DMGetApplicationContext(da, &user)); 118 tleft = user->tleft; 119 /* Get ghost points */ 120 PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 121 PetscCall(DMDAVecGetArray(da, X, &x)); 122 123 /* Compute initial guess */ 124 for (j = ys; j < ys + ym; j++) { 125 for (i = xs; i < xs + xm; i++) { x[j][i] = tleft; } 126 } 127 PetscCall(DMDAVecRestoreArray(da, X, &x)); 128 PetscFunctionReturn(0); 129 } 130 /* -------------------- Evaluate Function F(x) --------------------- */ 131 PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr) { 132 AppCtx *user = (AppCtx *)ptr; 133 PetscInt i, j, mx, my, xs, ys, xm, ym; 134 PetscScalar zero = 0.0, one = 1.0; 135 PetscScalar hx, hy, hxdhy, hydhx; 136 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; 137 PetscScalar tleft, tright, beta; 138 PetscScalar **x, **f; 139 Vec localX; 140 DM da; 141 142 PetscFunctionBeginUser; 143 PetscCall(SNESGetDM(snes, &da)); 144 PetscCall(DMGetLocalVector(da, &localX)); 145 PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 146 hx = one / (PetscReal)(mx - 1); 147 hy = one / (PetscReal)(my - 1); 148 hxdhy = hx / hy; 149 hydhx = hy / hx; 150 tleft = user->tleft; 151 tright = user->tright; 152 beta = user->beta; 153 154 /* Get ghost points */ 155 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 156 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 157 PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 158 PetscCall(DMDAVecGetArray(da, localX, &x)); 159 PetscCall(DMDAVecGetArray(da, F, &f)); 160 161 /* Evaluate function */ 162 for (j = ys; j < ys + ym; j++) { 163 for (i = xs; i < xs + xm; i++) { 164 t0 = x[j][i]; 165 166 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) { 167 /* general interior volume */ 168 169 tw = x[j][i - 1]; 170 aw = 0.5 * (t0 + tw); 171 dw = PetscPowScalar(aw, beta); 172 fw = dw * (t0 - tw); 173 174 te = x[j][i + 1]; 175 ae = 0.5 * (t0 + te); 176 de = PetscPowScalar(ae, beta); 177 fe = de * (te - t0); 178 179 ts = x[j - 1][i]; 180 as = 0.5 * (t0 + ts); 181 ds = PetscPowScalar(as, beta); 182 fs = ds * (t0 - ts); 183 184 tn = x[j + 1][i]; 185 an = 0.5 * (t0 + tn); 186 dn = PetscPowScalar(an, beta); 187 fn = dn * (tn - t0); 188 189 } else if (i == 0) { 190 /* left-hand boundary */ 191 tw = tleft; 192 aw = 0.5 * (t0 + tw); 193 dw = PetscPowScalar(aw, beta); 194 fw = dw * (t0 - tw); 195 196 te = x[j][i + 1]; 197 ae = 0.5 * (t0 + te); 198 de = PetscPowScalar(ae, beta); 199 fe = de * (te - t0); 200 201 if (j > 0) { 202 ts = x[j - 1][i]; 203 as = 0.5 * (t0 + ts); 204 ds = PetscPowScalar(as, beta); 205 fs = ds * (t0 - ts); 206 } else fs = zero; 207 208 if (j < my - 1) { 209 tn = x[j + 1][i]; 210 an = 0.5 * (t0 + tn); 211 dn = PetscPowScalar(an, beta); 212 fn = dn * (tn - t0); 213 } else fn = zero; 214 215 } else if (i == mx - 1) { 216 /* right-hand boundary */ 217 tw = x[j][i - 1]; 218 aw = 0.5 * (t0 + tw); 219 dw = PetscPowScalar(aw, beta); 220 fw = dw * (t0 - tw); 221 222 te = tright; 223 ae = 0.5 * (t0 + te); 224 de = PetscPowScalar(ae, beta); 225 fe = de * (te - t0); 226 227 if (j > 0) { 228 ts = x[j - 1][i]; 229 as = 0.5 * (t0 + ts); 230 ds = PetscPowScalar(as, beta); 231 fs = ds * (t0 - ts); 232 } else fs = zero; 233 234 if (j < my - 1) { 235 tn = x[j + 1][i]; 236 an = 0.5 * (t0 + tn); 237 dn = PetscPowScalar(an, beta); 238 fn = dn * (tn - t0); 239 } else fn = zero; 240 241 } else if (j == 0) { 242 /* bottom boundary,and i <> 0 or mx-1 */ 243 tw = x[j][i - 1]; 244 aw = 0.5 * (t0 + tw); 245 dw = PetscPowScalar(aw, beta); 246 fw = dw * (t0 - tw); 247 248 te = x[j][i + 1]; 249 ae = 0.5 * (t0 + te); 250 de = PetscPowScalar(ae, beta); 251 fe = de * (te - t0); 252 253 fs = zero; 254 255 tn = x[j + 1][i]; 256 an = 0.5 * (t0 + tn); 257 dn = PetscPowScalar(an, beta); 258 fn = dn * (tn - t0); 259 260 } else if (j == my - 1) { 261 /* top boundary,and i <> 0 or mx-1 */ 262 tw = x[j][i - 1]; 263 aw = 0.5 * (t0 + tw); 264 dw = PetscPowScalar(aw, beta); 265 fw = dw * (t0 - tw); 266 267 te = x[j][i + 1]; 268 ae = 0.5 * (t0 + te); 269 de = PetscPowScalar(ae, beta); 270 fe = de * (te - t0); 271 272 ts = x[j - 1][i]; 273 as = 0.5 * (t0 + ts); 274 ds = PetscPowScalar(as, beta); 275 fs = ds * (t0 - ts); 276 277 fn = zero; 278 } 279 280 f[j][i] = -hydhx * (fe - fw) - hxdhy * (fn - fs); 281 } 282 } 283 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 284 PetscCall(DMDAVecRestoreArray(da, F, &f)); 285 PetscCall(DMRestoreLocalVector(da, &localX)); 286 PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm)); 287 PetscFunctionReturn(0); 288 } 289 /* -------------------- Evaluate Jacobian F(x) --------------------- */ 290 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat jac, Mat B, void *ptr) { 291 AppCtx *user = (AppCtx *)ptr; 292 PetscInt i, j, mx, my, xs, ys, xm, ym; 293 PetscScalar one = 1.0, hx, hy, hxdhy, hydhx, t0, tn, ts, te, tw; 294 PetscScalar dn, ds, de, dw, an, as, ae, aw, bn, bs, be, bw, gn, gs, ge, gw; 295 PetscScalar tleft, tright, beta, bm1, coef; 296 PetscScalar v[5], **x; 297 Vec localX; 298 MatStencil col[5], row; 299 DM da; 300 301 PetscFunctionBeginUser; 302 PetscCall(SNESGetDM(snes, &da)); 303 PetscCall(DMGetLocalVector(da, &localX)); 304 PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 305 hx = one / (PetscReal)(mx - 1); 306 hy = one / (PetscReal)(my - 1); 307 hxdhy = hx / hy; 308 hydhx = hy / hx; 309 tleft = user->tleft; 310 tright = user->tright; 311 beta = user->beta; 312 bm1 = user->bm1; 313 coef = user->coef; 314 315 /* Get ghost points */ 316 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 317 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 318 PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 319 PetscCall(DMDAVecGetArray(da, localX, &x)); 320 321 /* Evaluate Jacobian of function */ 322 for (j = ys; j < ys + ym; j++) { 323 for (i = xs; i < xs + xm; i++) { 324 t0 = x[j][i]; 325 326 if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) { 327 /* general interior volume */ 328 329 tw = x[j][i - 1]; 330 aw = 0.5 * (t0 + tw); 331 bw = PetscPowScalar(aw, bm1); 332 /* dw = bw * aw */ 333 dw = PetscPowScalar(aw, beta); 334 gw = coef * bw * (t0 - tw); 335 336 te = x[j][i + 1]; 337 ae = 0.5 * (t0 + te); 338 be = PetscPowScalar(ae, bm1); 339 /* de = be * ae; */ 340 de = PetscPowScalar(ae, beta); 341 ge = coef * be * (te - t0); 342 343 ts = x[j - 1][i]; 344 as = 0.5 * (t0 + ts); 345 bs = PetscPowScalar(as, bm1); 346 /* ds = bs * as; */ 347 ds = PetscPowScalar(as, beta); 348 gs = coef * bs * (t0 - ts); 349 350 tn = x[j + 1][i]; 351 an = 0.5 * (t0 + tn); 352 bn = PetscPowScalar(an, bm1); 353 /* dn = bn * an; */ 354 dn = PetscPowScalar(an, beta); 355 gn = coef * bn * (tn - t0); 356 357 v[0] = -hxdhy * (ds - gs); 358 col[0].j = j - 1; 359 col[0].i = i; 360 v[1] = -hydhx * (dw - gw); 361 col[1].j = j; 362 col[1].i = i - 1; 363 v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 364 col[2].j = row.j = j; 365 col[2].i = row.i = i; 366 v[3] = -hydhx * (de + ge); 367 col[3].j = j; 368 col[3].i = i + 1; 369 v[4] = -hxdhy * (dn + gn); 370 col[4].j = j + 1; 371 col[4].i = i; 372 PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES)); 373 374 } else if (i == 0) { 375 /* left-hand boundary */ 376 tw = tleft; 377 aw = 0.5 * (t0 + tw); 378 bw = PetscPowScalar(aw, bm1); 379 /* dw = bw * aw */ 380 dw = PetscPowScalar(aw, beta); 381 gw = coef * bw * (t0 - tw); 382 383 te = x[j][i + 1]; 384 ae = 0.5 * (t0 + te); 385 be = PetscPowScalar(ae, bm1); 386 /* de = be * ae; */ 387 de = PetscPowScalar(ae, beta); 388 ge = coef * be * (te - t0); 389 390 /* left-hand bottom boundary */ 391 if (j == 0) { 392 tn = x[j + 1][i]; 393 an = 0.5 * (t0 + tn); 394 bn = PetscPowScalar(an, bm1); 395 /* dn = bn * an; */ 396 dn = PetscPowScalar(an, beta); 397 gn = coef * bn * (tn - t0); 398 399 v[0] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 400 col[0].j = row.j = j; 401 col[0].i = row.i = i; 402 v[1] = -hydhx * (de + ge); 403 col[1].j = j; 404 col[1].i = i + 1; 405 v[2] = -hxdhy * (dn + gn); 406 col[2].j = j + 1; 407 col[2].i = i; 408 PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 409 410 /* left-hand interior boundary */ 411 } else if (j < my - 1) { 412 ts = x[j - 1][i]; 413 as = 0.5 * (t0 + ts); 414 bs = PetscPowScalar(as, bm1); 415 /* ds = bs * as; */ 416 ds = PetscPowScalar(as, beta); 417 gs = coef * bs * (t0 - ts); 418 419 tn = x[j + 1][i]; 420 an = 0.5 * (t0 + tn); 421 bn = PetscPowScalar(an, bm1); 422 /* dn = bn * an; */ 423 dn = PetscPowScalar(an, beta); 424 gn = coef * bn * (tn - t0); 425 426 v[0] = -hxdhy * (ds - gs); 427 col[0].j = j - 1; 428 col[0].i = i; 429 v[1] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 430 col[1].j = row.j = j; 431 col[1].i = row.i = i; 432 v[2] = -hydhx * (de + ge); 433 col[2].j = j; 434 col[2].i = i + 1; 435 v[3] = -hxdhy * (dn + gn); 436 col[3].j = j + 1; 437 col[3].i = i; 438 PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 439 /* left-hand top boundary */ 440 } else { 441 ts = x[j - 1][i]; 442 as = 0.5 * (t0 + ts); 443 bs = PetscPowScalar(as, bm1); 444 /* ds = bs * as; */ 445 ds = PetscPowScalar(as, beta); 446 gs = coef * bs * (t0 - ts); 447 448 v[0] = -hxdhy * (ds - gs); 449 col[0].j = j - 1; 450 col[0].i = i; 451 v[1] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 452 col[1].j = row.j = j; 453 col[1].i = row.i = i; 454 v[2] = -hydhx * (de + ge); 455 col[2].j = j; 456 col[2].i = i + 1; 457 PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 458 } 459 460 } else if (i == mx - 1) { 461 /* right-hand boundary */ 462 tw = x[j][i - 1]; 463 aw = 0.5 * (t0 + tw); 464 bw = PetscPowScalar(aw, bm1); 465 /* dw = bw * aw */ 466 dw = PetscPowScalar(aw, beta); 467 gw = coef * bw * (t0 - tw); 468 469 te = tright; 470 ae = 0.5 * (t0 + te); 471 be = PetscPowScalar(ae, bm1); 472 /* de = be * ae; */ 473 de = PetscPowScalar(ae, beta); 474 ge = coef * be * (te - t0); 475 476 /* right-hand bottom boundary */ 477 if (j == 0) { 478 tn = x[j + 1][i]; 479 an = 0.5 * (t0 + tn); 480 bn = PetscPowScalar(an, bm1); 481 /* dn = bn * an; */ 482 dn = PetscPowScalar(an, beta); 483 gn = coef * bn * (tn - t0); 484 485 v[0] = -hydhx * (dw - gw); 486 col[0].j = j; 487 col[0].i = i - 1; 488 v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 489 col[1].j = row.j = j; 490 col[1].i = row.i = i; 491 v[2] = -hxdhy * (dn + gn); 492 col[2].j = j + 1; 493 col[2].i = i; 494 PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 495 496 /* right-hand interior boundary */ 497 } else if (j < my - 1) { 498 ts = x[j - 1][i]; 499 as = 0.5 * (t0 + ts); 500 bs = PetscPowScalar(as, bm1); 501 /* ds = bs * as; */ 502 ds = PetscPowScalar(as, beta); 503 gs = coef * bs * (t0 - ts); 504 505 tn = x[j + 1][i]; 506 an = 0.5 * (t0 + tn); 507 bn = PetscPowScalar(an, bm1); 508 /* dn = bn * an; */ 509 dn = PetscPowScalar(an, beta); 510 gn = coef * bn * (tn - t0); 511 512 v[0] = -hxdhy * (ds - gs); 513 col[0].j = j - 1; 514 col[0].i = i; 515 v[1] = -hydhx * (dw - gw); 516 col[1].j = j; 517 col[1].i = i - 1; 518 v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 519 col[2].j = row.j = j; 520 col[2].i = row.i = i; 521 v[3] = -hxdhy * (dn + gn); 522 col[3].j = j + 1; 523 col[3].i = i; 524 PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 525 /* right-hand top boundary */ 526 } else { 527 ts = x[j - 1][i]; 528 as = 0.5 * (t0 + ts); 529 bs = PetscPowScalar(as, bm1); 530 /* ds = bs * as; */ 531 ds = PetscPowScalar(as, beta); 532 gs = coef * bs * (t0 - ts); 533 534 v[0] = -hxdhy * (ds - gs); 535 col[0].j = j - 1; 536 col[0].i = i; 537 v[1] = -hydhx * (dw - gw); 538 col[1].j = j; 539 col[1].i = i - 1; 540 v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 541 col[2].j = row.j = j; 542 col[2].i = row.i = i; 543 PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 544 } 545 546 /* bottom boundary,and i <> 0 or mx-1 */ 547 } else if (j == 0) { 548 tw = x[j][i - 1]; 549 aw = 0.5 * (t0 + tw); 550 bw = PetscPowScalar(aw, bm1); 551 /* dw = bw * aw */ 552 dw = PetscPowScalar(aw, beta); 553 gw = coef * bw * (t0 - tw); 554 555 te = x[j][i + 1]; 556 ae = 0.5 * (t0 + te); 557 be = PetscPowScalar(ae, bm1); 558 /* de = be * ae; */ 559 de = PetscPowScalar(ae, beta); 560 ge = coef * be * (te - t0); 561 562 tn = x[j + 1][i]; 563 an = 0.5 * (t0 + tn); 564 bn = PetscPowScalar(an, bm1); 565 /* dn = bn * an; */ 566 dn = PetscPowScalar(an, beta); 567 gn = coef * bn * (tn - t0); 568 569 v[0] = -hydhx * (dw - gw); 570 col[0].j = j; 571 col[0].i = i - 1; 572 v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 573 col[1].j = row.j = j; 574 col[1].i = row.i = i; 575 v[2] = -hydhx * (de + ge); 576 col[2].j = j; 577 col[2].i = i + 1; 578 v[3] = -hxdhy * (dn + gn); 579 col[3].j = j + 1; 580 col[3].i = i; 581 PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 582 583 /* top boundary,and i <> 0 or mx-1 */ 584 } else if (j == my - 1) { 585 tw = x[j][i - 1]; 586 aw = 0.5 * (t0 + tw); 587 bw = PetscPowScalar(aw, bm1); 588 /* dw = bw * aw */ 589 dw = PetscPowScalar(aw, beta); 590 gw = coef * bw * (t0 - tw); 591 592 te = x[j][i + 1]; 593 ae = 0.5 * (t0 + te); 594 be = PetscPowScalar(ae, bm1); 595 /* de = be * ae; */ 596 de = PetscPowScalar(ae, beta); 597 ge = coef * be * (te - t0); 598 599 ts = x[j - 1][i]; 600 as = 0.5 * (t0 + ts); 601 bs = PetscPowScalar(as, bm1); 602 /* ds = bs * as; */ 603 ds = PetscPowScalar(as, beta); 604 gs = coef * bs * (t0 - ts); 605 606 v[0] = -hxdhy * (ds - gs); 607 col[0].j = j - 1; 608 col[0].i = i; 609 v[1] = -hydhx * (dw - gw); 610 col[1].j = j; 611 col[1].i = i - 1; 612 v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 613 col[2].j = row.j = j; 614 col[2].i = row.i = i; 615 v[3] = -hydhx * (de + ge); 616 col[3].j = j; 617 col[3].i = i + 1; 618 PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 619 } 620 } 621 } 622 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 623 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 624 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 625 PetscCall(DMRestoreLocalVector(da, &localX)); 626 if (jac != B) { 627 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 628 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 629 } 630 631 PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym)); 632 PetscFunctionReturn(0); 633 } 634 635 /*TEST 636 637 test: 638 suffix: 1 639 args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view 640 requires: !single 641 642 test: 643 suffix: 2 644 args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc 645 requires: !single 646 647 TEST*/ 648