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