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