1 2 #include <petscsnes.h> 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 6 static const char help[] = "Parallel version of the minimum surface area problem in 2D using DMDA.\n\ 7 It solves a system of nonlinear equations in mixed\n\ 8 complementarity form.This example is based on a\n\ 9 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and\n\ 10 boundary values along the edges of the domain, the objective is to find the\n\ 11 surface with the minimal area that satisfies the boundary conditions.\n\ 12 This application solves this problem using complimentarity -- We are actually\n\ 13 solving the system (grad f)_i >= 0, if x_i == l_i \n\ 14 (grad f)_i = 0, if l_i < x_i < u_i \n\ 15 (grad f)_i <= 0, if x_i == u_i \n\ 16 where f is the function to be minimized. \n\ 17 \n\ 18 The command line options are:\n\ 19 -da_grid_x <nx>, where <nx> = number of grid points in the 1st coordinate direction\n\ 20 -da_grid_y <ny>, where <ny> = number of grid points in the 2nd coordinate direction\n\ 21 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise\n\ 22 -lb <value>, lower bound on the variables\n\ 23 -ub <value>, upper bound on the variables\n\n"; 24 25 /* 26 User-defined application context - contains data needed by the 27 application-provided call-back routines, FormJacobian() and 28 FormFunction(). 29 */ 30 31 /* 32 This is a new version of the ../tests/ex8.c code 33 34 Run, for example, with the options ./ex58 -snes_vi_monitor -ksp_monitor -mg_levels_ksp_monitor -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin pmat -ksp_type fgmres 35 36 Or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of 37 multigrid levels, it will be determined automatically based on the number of refinements done) 38 39 ./ex58 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3 40 -mg_levels_ksp_monitor -snes_vi_monitor -mg_levels_pc_type sor -pc_mg_type full 41 42 */ 43 44 typedef struct { 45 PetscScalar *bottom, *top, *left, *right; 46 PetscScalar lb, ub; 47 } AppCtx; 48 49 /* -------- User-defined Routines --------- */ 50 51 extern PetscErrorCode FormBoundaryConditions(SNES, AppCtx **); 52 extern PetscErrorCode DestroyBoundaryConditions(AppCtx **); 53 extern PetscErrorCode ComputeInitialGuess(SNES, Vec, void *); 54 extern PetscErrorCode FormGradient(SNES, Vec, Vec, void *); 55 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 56 extern PetscErrorCode FormBounds(SNES, Vec, Vec); 57 58 int main(int argc, char **argv) 59 { 60 Vec x, r; /* solution and residual vectors */ 61 SNES snes; /* nonlinear solver context */ 62 Mat J; /* Jacobian matrix */ 63 DM da; 64 65 PetscFunctionBeginUser; 66 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 67 68 /* Create distributed array to manage the 2d grid */ 69 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da)); 70 PetscCall(DMSetFromOptions(da)); 71 PetscCall(DMSetUp(da)); 72 73 /* Extract global vectors from DMDA; */ 74 PetscCall(DMCreateGlobalVector(da, &x)); 75 PetscCall(VecDuplicate(x, &r)); 76 77 PetscCall(DMSetMatType(da, MATAIJ)); 78 PetscCall(DMCreateMatrix(da, &J)); 79 80 /* Create nonlinear solver context */ 81 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 82 PetscCall(SNESSetDM(snes, da)); 83 84 /* Set function evaluation and Jacobian evaluation routines */ 85 PetscCall(SNESSetFunction(snes, r, FormGradient, NULL)); 86 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL)); 87 88 PetscCall(SNESSetComputeApplicationContext(snes, (PetscErrorCode(*)(SNES, void **))FormBoundaryConditions, (PetscErrorCode(*)(void **))DestroyBoundaryConditions)); 89 90 PetscCall(SNESSetComputeInitialGuess(snes, ComputeInitialGuess, NULL)); 91 92 PetscCall(SNESVISetComputeVariableBounds(snes, FormBounds)); 93 94 PetscCall(SNESSetFromOptions(snes)); 95 96 /* Solve the application */ 97 PetscCall(SNESSolve(snes, NULL, x)); 98 99 /* Free memory */ 100 PetscCall(VecDestroy(&x)); 101 PetscCall(VecDestroy(&r)); 102 PetscCall(MatDestroy(&J)); 103 PetscCall(SNESDestroy(&snes)); 104 105 /* Free user-created data structures */ 106 PetscCall(DMDestroy(&da)); 107 108 PetscCall(PetscFinalize()); 109 return 0; 110 } 111 112 /* -------------------------------------------------------------------- */ 113 114 /* FormBounds - sets the upper and lower bounds 115 116 Input Parameters: 117 . snes - the SNES context 118 119 Output Parameters: 120 . xl - lower bounds 121 . xu - upper bounds 122 */ 123 PetscErrorCode FormBounds(SNES snes, Vec xl, Vec xu) 124 { 125 AppCtx *ctx; 126 127 PetscFunctionBeginUser; 128 PetscCall(SNESGetApplicationContext(snes, &ctx)); 129 PetscCall(VecSet(xl, ctx->lb)); 130 PetscCall(VecSet(xu, ctx->ub)); 131 PetscFunctionReturn(0); 132 } 133 134 /* -------------------------------------------------------------------- */ 135 136 /* FormGradient - Evaluates gradient of f. 137 138 Input Parameters: 139 . snes - the SNES context 140 . X - input vector 141 . ptr - optional user-defined context, as set by SNESSetFunction() 142 143 Output Parameters: 144 . G - vector containing the newly evaluated gradient 145 */ 146 PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr) 147 { 148 AppCtx *user; 149 PetscInt i, j; 150 PetscInt mx, my; 151 PetscScalar hx, hy, hydhx, hxdhy; 152 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 153 PetscScalar df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc; 154 PetscScalar **g, **x; 155 PetscInt xs, xm, ys, ym; 156 Vec localX; 157 DM da; 158 159 PetscFunctionBeginUser; 160 PetscCall(SNESGetDM(snes, &da)); 161 PetscCall(SNESGetApplicationContext(snes, &user)); 162 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 163 hx = 1.0 / (mx + 1); 164 hy = 1.0 / (my + 1); 165 hydhx = hy / hx; 166 hxdhy = hx / hy; 167 168 PetscCall(VecSet(G, 0.0)); 169 170 /* Get local vector */ 171 PetscCall(DMGetLocalVector(da, &localX)); 172 /* Get ghost points */ 173 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 174 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 175 /* Get pointer to local vector data */ 176 PetscCall(DMDAVecGetArray(da, localX, &x)); 177 PetscCall(DMDAVecGetArray(da, G, &g)); 178 179 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 180 /* Compute function over the locally owned part of the mesh */ 181 for (j = ys; j < ys + ym; j++) { 182 for (i = xs; i < xs + xm; i++) { 183 xc = x[j][i]; 184 xlt = xrb = xl = xr = xb = xt = xc; 185 186 if (i == 0) { /* left side */ 187 xl = user->left[j + 1]; 188 xlt = user->left[j + 2]; 189 } else xl = x[j][i - 1]; 190 191 if (j == 0) { /* bottom side */ 192 xb = user->bottom[i + 1]; 193 xrb = user->bottom[i + 2]; 194 } else xb = x[j - 1][i]; 195 196 if (i + 1 == mx) { /* right side */ 197 xr = user->right[j + 1]; 198 xrb = user->right[j]; 199 } else xr = x[j][i + 1]; 200 201 if (j + 1 == 0 + my) { /* top side */ 202 xt = user->top[i + 1]; 203 xlt = user->top[i]; 204 } else xt = x[j + 1][i]; 205 206 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; /* left top side */ 207 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; /* right bottom */ 208 209 d1 = (xc - xl); 210 d2 = (xc - xr); 211 d3 = (xc - xt); 212 d4 = (xc - xb); 213 d5 = (xr - xrb); 214 d6 = (xrb - xb); 215 d7 = (xlt - xl); 216 d8 = (xt - xlt); 217 218 df1dxc = d1 * hydhx; 219 df2dxc = (d1 * hydhx + d4 * hxdhy); 220 df3dxc = d3 * hxdhy; 221 df4dxc = (d2 * hydhx + d3 * hxdhy); 222 df5dxc = d2 * hydhx; 223 df6dxc = d4 * hxdhy; 224 225 d1 /= hx; 226 d2 /= hx; 227 d3 /= hy; 228 d4 /= hy; 229 d5 /= hy; 230 d6 /= hx; 231 d7 /= hy; 232 d8 /= hx; 233 234 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 235 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 236 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 237 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 238 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 239 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 240 241 df1dxc /= f1; 242 df2dxc /= f2; 243 df3dxc /= f3; 244 df4dxc /= f4; 245 df5dxc /= f5; 246 df6dxc /= f6; 247 248 g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0; 249 } 250 } 251 252 /* Restore vectors */ 253 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 254 PetscCall(DMDAVecRestoreArray(da, G, &g)); 255 PetscCall(DMRestoreLocalVector(da, &localX)); 256 PetscCall(PetscLogFlops(67.0 * mx * my)); 257 PetscFunctionReturn(0); 258 } 259 260 /* ------------------------------------------------------------------- */ 261 /* 262 FormJacobian - Evaluates Jacobian matrix. 263 264 Input Parameters: 265 . snes - SNES context 266 . X - input vector 267 . ptr - optional user-defined context, as set by SNESSetJacobian() 268 269 Output Parameters: 270 . tH - Jacobian matrix 271 272 */ 273 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat H, Mat tHPre, void *ptr) 274 { 275 AppCtx *user; 276 PetscInt i, j, k; 277 PetscInt mx, my; 278 MatStencil row, col[7]; 279 PetscScalar hx, hy, hydhx, hxdhy; 280 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 281 PetscScalar hl, hr, ht, hb, hc, htl, hbr; 282 PetscScalar **x, v[7]; 283 PetscBool assembled; 284 PetscInt xs, xm, ys, ym; 285 Vec localX; 286 DM da; 287 288 PetscFunctionBeginUser; 289 PetscCall(SNESGetDM(snes, &da)); 290 PetscCall(SNESGetApplicationContext(snes, &user)); 291 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 292 hx = 1.0 / (mx + 1); 293 hy = 1.0 / (my + 1); 294 hydhx = hy / hx; 295 hxdhy = hx / hy; 296 297 /* Set various matrix options */ 298 PetscCall(MatAssembled(H, &assembled)); 299 if (assembled) PetscCall(MatZeroEntries(H)); 300 301 /* Get local vector */ 302 PetscCall(DMGetLocalVector(da, &localX)); 303 /* Get ghost points */ 304 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 305 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 306 307 /* Get pointers to vector data */ 308 PetscCall(DMDAVecGetArray(da, localX, &x)); 309 310 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 311 /* Compute Jacobian over the locally owned part of the mesh */ 312 for (j = ys; j < ys + ym; j++) { 313 for (i = xs; i < xs + xm; i++) { 314 xc = x[j][i]; 315 xlt = xrb = xl = xr = xb = xt = xc; 316 317 /* Left */ 318 if (i == 0) { 319 xl = user->left[j + 1]; 320 xlt = user->left[j + 2]; 321 } else xl = x[j][i - 1]; 322 323 /* Bottom */ 324 if (j == 0) { 325 xb = user->bottom[i + 1]; 326 xrb = user->bottom[i + 2]; 327 } else xb = x[j - 1][i]; 328 329 /* Right */ 330 if (i + 1 == mx) { 331 xr = user->right[j + 1]; 332 xrb = user->right[j]; 333 } else xr = x[j][i + 1]; 334 335 /* Top */ 336 if (j + 1 == my) { 337 xt = user->top[i + 1]; 338 xlt = user->top[i]; 339 } else xt = x[j + 1][i]; 340 341 /* Top left */ 342 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; 343 344 /* Bottom right */ 345 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; 346 347 d1 = (xc - xl) / hx; 348 d2 = (xc - xr) / hx; 349 d3 = (xc - xt) / hy; 350 d4 = (xc - xb) / hy; 351 d5 = (xrb - xr) / hy; 352 d6 = (xrb - xb) / hx; 353 d7 = (xlt - xl) / hy; 354 d8 = (xlt - xt) / hx; 355 356 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 357 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 358 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 359 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 360 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 361 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 362 363 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 364 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 365 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 366 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 367 368 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 369 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 370 371 hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2.0 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2.0 * d2 * d3) / (f4 * f4 * f4); 372 373 hl /= 2.0; 374 hr /= 2.0; 375 ht /= 2.0; 376 hb /= 2.0; 377 hbr /= 2.0; 378 htl /= 2.0; 379 hc /= 2.0; 380 381 k = 0; 382 row.i = i; 383 row.j = j; 384 /* Bottom */ 385 if (j > 0) { 386 v[k] = hb; 387 col[k].i = i; 388 col[k].j = j - 1; 389 k++; 390 } 391 392 /* Bottom right */ 393 if (j > 0 && i < mx - 1) { 394 v[k] = hbr; 395 col[k].i = i + 1; 396 col[k].j = j - 1; 397 k++; 398 } 399 400 /* left */ 401 if (i > 0) { 402 v[k] = hl; 403 col[k].i = i - 1; 404 col[k].j = j; 405 k++; 406 } 407 408 /* Centre */ 409 v[k] = hc; 410 col[k].i = row.i; 411 col[k].j = row.j; 412 k++; 413 414 /* Right */ 415 if (i < mx - 1) { 416 v[k] = hr; 417 col[k].i = i + 1; 418 col[k].j = j; 419 k++; 420 } 421 422 /* Top left */ 423 if (i > 0 && j < my - 1) { 424 v[k] = htl; 425 col[k].i = i - 1; 426 col[k].j = j + 1; 427 k++; 428 } 429 430 /* Top */ 431 if (j < my - 1) { 432 v[k] = ht; 433 col[k].i = i; 434 col[k].j = j + 1; 435 k++; 436 } 437 438 PetscCall(MatSetValuesStencil(H, 1, &row, k, col, v, INSERT_VALUES)); 439 } 440 } 441 442 /* Assemble the matrix */ 443 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 444 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 445 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 446 PetscCall(DMRestoreLocalVector(da, &localX)); 447 448 PetscCall(PetscLogFlops(199.0 * mx * my)); 449 PetscFunctionReturn(0); 450 } 451 452 /* ------------------------------------------------------------------- */ 453 /* 454 FormBoundaryConditions - Calculates the boundary conditions for 455 the region. 456 457 Input Parameter: 458 . user - user-defined application context 459 460 Output Parameter: 461 . user - user-defined application context 462 */ 463 PetscErrorCode FormBoundaryConditions(SNES snes, AppCtx **ouser) 464 { 465 PetscInt i, j, k, limit = 0, maxits = 5; 466 PetscInt mx, my; 467 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 468 PetscScalar one = 1.0, two = 2.0, three = 3.0; 469 PetscScalar det, hx, hy, xt = 0, yt = 0; 470 PetscReal fnorm, tol = 1e-10; 471 PetscScalar u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 472 PetscScalar b = -0.5, t = 0.5, l = -0.5, r = 0.5; 473 PetscScalar *boundary; 474 AppCtx *user; 475 DM da; 476 477 PetscFunctionBeginUser; 478 PetscCall(SNESGetDM(snes, &da)); 479 PetscCall(PetscNew(&user)); 480 *ouser = user; 481 user->lb = .05; 482 user->ub = PETSC_INFINITY; 483 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 484 485 /* Check if lower and upper bounds are set */ 486 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lb", &user->lb, 0)); 487 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ub", &user->ub, 0)); 488 bsize = mx + 2; 489 lsize = my + 2; 490 rsize = my + 2; 491 tsize = mx + 2; 492 493 PetscCall(PetscMalloc1(bsize, &user->bottom)); 494 PetscCall(PetscMalloc1(tsize, &user->top)); 495 PetscCall(PetscMalloc1(lsize, &user->left)); 496 PetscCall(PetscMalloc1(rsize, &user->right)); 497 498 hx = (r - l) / (mx + 1.0); 499 hy = (t - b) / (my + 1.0); 500 501 for (j = 0; j < 4; j++) { 502 if (j == 0) { 503 yt = b; 504 xt = l; 505 limit = bsize; 506 boundary = user->bottom; 507 } else if (j == 1) { 508 yt = t; 509 xt = l; 510 limit = tsize; 511 boundary = user->top; 512 } else if (j == 2) { 513 yt = b; 514 xt = l; 515 limit = lsize; 516 boundary = user->left; 517 } else { /* if (j==3) */ 518 yt = b; 519 xt = r; 520 limit = rsize; 521 boundary = user->right; 522 } 523 524 for (i = 0; i < limit; i++) { 525 u1 = xt; 526 u2 = -yt; 527 for (k = 0; k < maxits; k++) { 528 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 529 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 530 fnorm = PetscRealPart(PetscSqrtScalar(nf1 * nf1 + nf2 * nf2)); 531 if (fnorm <= tol) break; 532 njac11 = one + u2 * u2 - u1 * u1; 533 njac12 = two * u1 * u2; 534 njac21 = -two * u1 * u2; 535 njac22 = -one - u1 * u1 + u2 * u2; 536 det = njac11 * njac22 - njac21 * njac12; 537 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 538 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 539 } 540 541 boundary[i] = u1 * u1 - u2 * u2; 542 if (j == 0 || j == 1) xt = xt + hx; 543 else yt = yt + hy; /* if (j==2 || j==3) */ 544 } 545 } 546 PetscFunctionReturn(0); 547 } 548 549 PetscErrorCode DestroyBoundaryConditions(AppCtx **ouser) 550 { 551 AppCtx *user = *ouser; 552 553 PetscFunctionBeginUser; 554 PetscCall(PetscFree(user->bottom)); 555 PetscCall(PetscFree(user->top)); 556 PetscCall(PetscFree(user->left)); 557 PetscCall(PetscFree(user->right)); 558 PetscCall(PetscFree(*ouser)); 559 PetscFunctionReturn(0); 560 } 561 562 /* ------------------------------------------------------------------- */ 563 /* 564 ComputeInitialGuess - Calculates the initial guess 565 566 Input Parameters: 567 . user - user-defined application context 568 . X - vector for initial guess 569 570 Output Parameters: 571 . X - newly computed initial guess 572 */ 573 PetscErrorCode ComputeInitialGuess(SNES snes, Vec X, void *dummy) 574 { 575 PetscInt i, j, mx, my; 576 DM da; 577 AppCtx *user; 578 PetscScalar **x; 579 PetscInt xs, xm, ys, ym; 580 581 PetscFunctionBeginUser; 582 PetscCall(SNESGetDM(snes, &da)); 583 PetscCall(SNESGetApplicationContext(snes, &user)); 584 585 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 586 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 587 588 /* Get pointers to vector data */ 589 PetscCall(DMDAVecGetArray(da, X, &x)); 590 /* Perform local computations */ 591 for (j = ys; j < ys + ym; j++) { 592 for (i = xs; i < xs + xm; i++) x[j][i] = (((j + 1.0) * user->bottom[i + 1] + (my - j + 1.0) * user->top[i + 1]) / (my + 2.0) + ((i + 1.0) * user->left[j + 1] + (mx - i + 1.0) * user->right[j + 1]) / (mx + 2.0)) / 2.0; 593 } 594 /* Restore vectors */ 595 PetscCall(DMDAVecRestoreArray(da, X, &x)); 596 PetscFunctionReturn(0); 597 } 598 599 /*TEST 600 601 test: 602 args: -snes_type vinewtonrsls -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason 603 requires: !single 604 605 test: 606 suffix: 2 607 args: -snes_type vinewtonssls -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason 608 requires: !single 609 610 TEST*/ 611