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 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 AppCtx *ctx; 124 125 PetscFunctionBeginUser; 126 PetscCall(SNESGetApplicationContext(snes, &ctx)); 127 PetscCall(VecSet(xl, ctx->lb)); 128 PetscCall(VecSet(xu, ctx->ub)); 129 PetscFunctionReturn(0); 130 } 131 132 /* -------------------------------------------------------------------- */ 133 134 /* FormGradient - Evaluates gradient of f. 135 136 Input Parameters: 137 . snes - the SNES context 138 . X - input vector 139 . ptr - optional user-defined context, as set by SNESSetFunction() 140 141 Output Parameters: 142 . G - vector containing the newly evaluated gradient 143 */ 144 PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr) { 145 AppCtx *user; 146 PetscInt i, j; 147 PetscInt mx, my; 148 PetscScalar hx, hy, hydhx, hxdhy; 149 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 150 PetscScalar df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc; 151 PetscScalar **g, **x; 152 PetscInt xs, xm, ys, ym; 153 Vec localX; 154 DM da; 155 156 PetscFunctionBeginUser; 157 PetscCall(SNESGetDM(snes, &da)); 158 PetscCall(SNESGetApplicationContext(snes, &user)); 159 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)); 160 hx = 1.0 / (mx + 1); 161 hy = 1.0 / (my + 1); 162 hydhx = hy / hx; 163 hxdhy = hx / hy; 164 165 PetscCall(VecSet(G, 0.0)); 166 167 /* Get local vector */ 168 PetscCall(DMGetLocalVector(da, &localX)); 169 /* Get ghost points */ 170 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 171 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 172 /* Get pointer to local vector data */ 173 PetscCall(DMDAVecGetArray(da, localX, &x)); 174 PetscCall(DMDAVecGetArray(da, G, &g)); 175 176 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 177 /* Compute function over the locally owned part of the mesh */ 178 for (j = ys; j < ys + ym; j++) { 179 for (i = xs; i < xs + xm; i++) { 180 xc = x[j][i]; 181 xlt = xrb = xl = xr = xb = xt = xc; 182 183 if (i == 0) { /* left side */ 184 xl = user->left[j + 1]; 185 xlt = user->left[j + 2]; 186 } else xl = x[j][i - 1]; 187 188 if (j == 0) { /* bottom side */ 189 xb = user->bottom[i + 1]; 190 xrb = user->bottom[i + 2]; 191 } else xb = x[j - 1][i]; 192 193 if (i + 1 == mx) { /* right side */ 194 xr = user->right[j + 1]; 195 xrb = user->right[j]; 196 } else xr = x[j][i + 1]; 197 198 if (j + 1 == 0 + my) { /* top side */ 199 xt = user->top[i + 1]; 200 xlt = user->top[i]; 201 } else xt = x[j + 1][i]; 202 203 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; /* left top side */ 204 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; /* right bottom */ 205 206 d1 = (xc - xl); 207 d2 = (xc - xr); 208 d3 = (xc - xt); 209 d4 = (xc - xb); 210 d5 = (xr - xrb); 211 d6 = (xrb - xb); 212 d7 = (xlt - xl); 213 d8 = (xt - xlt); 214 215 df1dxc = d1 * hydhx; 216 df2dxc = (d1 * hydhx + d4 * hxdhy); 217 df3dxc = d3 * hxdhy; 218 df4dxc = (d2 * hydhx + d3 * hxdhy); 219 df5dxc = d2 * hydhx; 220 df6dxc = d4 * hxdhy; 221 222 d1 /= hx; 223 d2 /= hx; 224 d3 /= hy; 225 d4 /= hy; 226 d5 /= hy; 227 d6 /= hx; 228 d7 /= hy; 229 d8 /= hx; 230 231 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 232 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 233 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 234 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 235 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 236 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 237 238 df1dxc /= f1; 239 df2dxc /= f2; 240 df3dxc /= f3; 241 df4dxc /= f4; 242 df5dxc /= f5; 243 df6dxc /= f6; 244 245 g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0; 246 } 247 } 248 249 /* Restore vectors */ 250 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 251 PetscCall(DMDAVecRestoreArray(da, G, &g)); 252 PetscCall(DMRestoreLocalVector(da, &localX)); 253 PetscCall(PetscLogFlops(67.0 * mx * my)); 254 PetscFunctionReturn(0); 255 } 256 257 /* ------------------------------------------------------------------- */ 258 /* 259 FormJacobian - Evaluates Jacobian matrix. 260 261 Input Parameters: 262 . snes - SNES context 263 . X - input vector 264 . ptr - optional user-defined context, as set by SNESSetJacobian() 265 266 Output Parameters: 267 . tH - Jacobian matrix 268 269 */ 270 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat H, Mat tHPre, void *ptr) { 271 AppCtx *user; 272 PetscInt i, j, k; 273 PetscInt mx, my; 274 MatStencil row, col[7]; 275 PetscScalar hx, hy, hydhx, hxdhy; 276 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 277 PetscScalar hl, hr, ht, hb, hc, htl, hbr; 278 PetscScalar **x, v[7]; 279 PetscBool assembled; 280 PetscInt xs, xm, ys, ym; 281 Vec localX; 282 DM da; 283 284 PetscFunctionBeginUser; 285 PetscCall(SNESGetDM(snes, &da)); 286 PetscCall(SNESGetApplicationContext(snes, &user)); 287 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)); 288 hx = 1.0 / (mx + 1); 289 hy = 1.0 / (my + 1); 290 hydhx = hy / hx; 291 hxdhy = hx / hy; 292 293 /* Set various matrix options */ 294 PetscCall(MatAssembled(H, &assembled)); 295 if (assembled) PetscCall(MatZeroEntries(H)); 296 297 /* Get local vector */ 298 PetscCall(DMGetLocalVector(da, &localX)); 299 /* Get ghost points */ 300 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 301 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 302 303 /* Get pointers to vector data */ 304 PetscCall(DMDAVecGetArray(da, localX, &x)); 305 306 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 307 /* Compute Jacobian over the locally owned part of the mesh */ 308 for (j = ys; j < ys + ym; j++) { 309 for (i = xs; i < xs + xm; i++) { 310 xc = x[j][i]; 311 xlt = xrb = xl = xr = xb = xt = xc; 312 313 /* Left */ 314 if (i == 0) { 315 xl = user->left[j + 1]; 316 xlt = user->left[j + 2]; 317 } else xl = x[j][i - 1]; 318 319 /* Bottom */ 320 if (j == 0) { 321 xb = user->bottom[i + 1]; 322 xrb = user->bottom[i + 2]; 323 } else xb = x[j - 1][i]; 324 325 /* Right */ 326 if (i + 1 == mx) { 327 xr = user->right[j + 1]; 328 xrb = user->right[j]; 329 } else xr = x[j][i + 1]; 330 331 /* Top */ 332 if (j + 1 == my) { 333 xt = user->top[i + 1]; 334 xlt = user->top[i]; 335 } else xt = x[j + 1][i]; 336 337 /* Top left */ 338 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; 339 340 /* Bottom right */ 341 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; 342 343 d1 = (xc - xl) / hx; 344 d2 = (xc - xr) / hx; 345 d3 = (xc - xt) / hy; 346 d4 = (xc - xb) / hy; 347 d5 = (xrb - xr) / hy; 348 d6 = (xrb - xb) / hx; 349 d7 = (xlt - xl) / hy; 350 d8 = (xlt - xt) / hx; 351 352 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 353 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 354 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 355 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 356 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 357 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 358 359 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 360 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 361 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 362 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 363 364 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 365 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 366 367 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); 368 369 hl /= 2.0; 370 hr /= 2.0; 371 ht /= 2.0; 372 hb /= 2.0; 373 hbr /= 2.0; 374 htl /= 2.0; 375 hc /= 2.0; 376 377 k = 0; 378 row.i = i; 379 row.j = j; 380 /* Bottom */ 381 if (j > 0) { 382 v[k] = hb; 383 col[k].i = i; 384 col[k].j = j - 1; 385 k++; 386 } 387 388 /* Bottom right */ 389 if (j > 0 && i < mx - 1) { 390 v[k] = hbr; 391 col[k].i = i + 1; 392 col[k].j = j - 1; 393 k++; 394 } 395 396 /* left */ 397 if (i > 0) { 398 v[k] = hl; 399 col[k].i = i - 1; 400 col[k].j = j; 401 k++; 402 } 403 404 /* Centre */ 405 v[k] = hc; 406 col[k].i = row.i; 407 col[k].j = row.j; 408 k++; 409 410 /* Right */ 411 if (i < mx - 1) { 412 v[k] = hr; 413 col[k].i = i + 1; 414 col[k].j = j; 415 k++; 416 } 417 418 /* Top left */ 419 if (i > 0 && j < my - 1) { 420 v[k] = htl; 421 col[k].i = i - 1; 422 col[k].j = j + 1; 423 k++; 424 } 425 426 /* Top */ 427 if (j < my - 1) { 428 v[k] = ht; 429 col[k].i = i; 430 col[k].j = j + 1; 431 k++; 432 } 433 434 PetscCall(MatSetValuesStencil(H, 1, &row, k, col, v, INSERT_VALUES)); 435 } 436 } 437 438 /* Assemble the matrix */ 439 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 440 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 441 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 442 PetscCall(DMRestoreLocalVector(da, &localX)); 443 444 PetscCall(PetscLogFlops(199.0 * mx * my)); 445 PetscFunctionReturn(0); 446 } 447 448 /* ------------------------------------------------------------------- */ 449 /* 450 FormBoundaryConditions - Calculates the boundary conditions for 451 the region. 452 453 Input Parameter: 454 . user - user-defined application context 455 456 Output Parameter: 457 . user - user-defined application context 458 */ 459 PetscErrorCode FormBoundaryConditions(SNES snes, AppCtx **ouser) { 460 PetscInt i, j, k, limit = 0, maxits = 5; 461 PetscInt mx, my; 462 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 463 PetscScalar one = 1.0, two = 2.0, three = 3.0; 464 PetscScalar det, hx, hy, xt = 0, yt = 0; 465 PetscReal fnorm, tol = 1e-10; 466 PetscScalar u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 467 PetscScalar b = -0.5, t = 0.5, l = -0.5, r = 0.5; 468 PetscScalar *boundary; 469 AppCtx *user; 470 DM da; 471 472 PetscFunctionBeginUser; 473 PetscCall(SNESGetDM(snes, &da)); 474 PetscCall(PetscNew(&user)); 475 *ouser = user; 476 user->lb = .05; 477 user->ub = PETSC_INFINITY; 478 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)); 479 480 /* Check if lower and upper bounds are set */ 481 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lb", &user->lb, 0)); 482 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ub", &user->ub, 0)); 483 bsize = mx + 2; 484 lsize = my + 2; 485 rsize = my + 2; 486 tsize = mx + 2; 487 488 PetscCall(PetscMalloc1(bsize, &user->bottom)); 489 PetscCall(PetscMalloc1(tsize, &user->top)); 490 PetscCall(PetscMalloc1(lsize, &user->left)); 491 PetscCall(PetscMalloc1(rsize, &user->right)); 492 493 hx = (r - l) / (mx + 1.0); 494 hy = (t - b) / (my + 1.0); 495 496 for (j = 0; j < 4; j++) { 497 if (j == 0) { 498 yt = b; 499 xt = l; 500 limit = bsize; 501 boundary = user->bottom; 502 } else if (j == 1) { 503 yt = t; 504 xt = l; 505 limit = tsize; 506 boundary = user->top; 507 } else if (j == 2) { 508 yt = b; 509 xt = l; 510 limit = lsize; 511 boundary = user->left; 512 } else { /* if (j==3) */ 513 yt = b; 514 xt = r; 515 limit = rsize; 516 boundary = user->right; 517 } 518 519 for (i = 0; i < limit; i++) { 520 u1 = xt; 521 u2 = -yt; 522 for (k = 0; k < maxits; k++) { 523 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 524 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 525 fnorm = PetscRealPart(PetscSqrtScalar(nf1 * nf1 + nf2 * nf2)); 526 if (fnorm <= tol) break; 527 njac11 = one + u2 * u2 - u1 * u1; 528 njac12 = two * u1 * u2; 529 njac21 = -two * u1 * u2; 530 njac22 = -one - u1 * u1 + u2 * u2; 531 det = njac11 * njac22 - njac21 * njac12; 532 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 533 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 534 } 535 536 boundary[i] = u1 * u1 - u2 * u2; 537 if (j == 0 || j == 1) xt = xt + hx; 538 else yt = yt + hy; /* if (j==2 || j==3) */ 539 } 540 } 541 PetscFunctionReturn(0); 542 } 543 544 PetscErrorCode DestroyBoundaryConditions(AppCtx **ouser) { 545 AppCtx *user = *ouser; 546 547 PetscFunctionBeginUser; 548 PetscCall(PetscFree(user->bottom)); 549 PetscCall(PetscFree(user->top)); 550 PetscCall(PetscFree(user->left)); 551 PetscCall(PetscFree(user->right)); 552 PetscCall(PetscFree(*ouser)); 553 PetscFunctionReturn(0); 554 } 555 556 /* ------------------------------------------------------------------- */ 557 /* 558 ComputeInitialGuess - Calculates the initial guess 559 560 Input Parameters: 561 . user - user-defined application context 562 . X - vector for initial guess 563 564 Output Parameters: 565 . X - newly computed initial guess 566 */ 567 PetscErrorCode ComputeInitialGuess(SNES snes, Vec X, void *dummy) { 568 PetscInt i, j, mx, my; 569 DM da; 570 AppCtx *user; 571 PetscScalar **x; 572 PetscInt xs, xm, ys, ym; 573 574 PetscFunctionBeginUser; 575 PetscCall(SNESGetDM(snes, &da)); 576 PetscCall(SNESGetApplicationContext(snes, &user)); 577 578 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 579 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)); 580 581 /* Get pointers to vector data */ 582 PetscCall(DMDAVecGetArray(da, X, &x)); 583 /* Perform local computations */ 584 for (j = ys; j < ys + ym; j++) { 585 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; } 586 } 587 /* Restore vectors */ 588 PetscCall(DMDAVecRestoreArray(da, X, &x)); 589 PetscFunctionReturn(0); 590 } 591 592 /*TEST 593 594 test: 595 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 596 requires: !single 597 598 test: 599 suffix: 2 600 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 601 requires: !single 602 603 TEST*/ 604