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