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 complementarity -- 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 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 32 33 Or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of 34 multigrid levels, it will be determined automatically based on the number of refinements done) 35 36 ./ex58 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3 37 -mg_levels_ksp_monitor -snes_vi_monitor -mg_levels_pc_type sor -pc_mg_type full 38 39 */ 40 41 typedef struct { 42 PetscScalar *bottom, *top, *left, *right; 43 PetscScalar lb, ub; 44 } AppCtx; 45 46 /* -------- User-defined Routines --------- */ 47 48 extern PetscErrorCode FormBoundaryConditions(SNES, void **); 49 extern PetscErrorCode DestroyBoundaryConditions(void **); 50 extern PetscErrorCode ComputeInitialGuess(SNES, Vec, void *); 51 extern PetscErrorCode FormGradient(SNES, Vec, Vec, void *); 52 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 53 extern PetscErrorCode FormBounds(SNES, Vec, Vec); 54 55 int main(int argc, char **argv) 56 { 57 Vec x, r; /* solution and residual vectors */ 58 SNES snes; /* nonlinear solver context */ 59 Mat J; /* Jacobian matrix */ 60 DM da; 61 62 PetscFunctionBeginUser; 63 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 64 65 /* Create distributed array to manage the 2d grid */ 66 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)); 67 PetscCall(DMSetFromOptions(da)); 68 PetscCall(DMSetUp(da)); 69 70 /* Extract global vectors from DMDA; */ 71 PetscCall(DMCreateGlobalVector(da, &x)); 72 PetscCall(VecDuplicate(x, &r)); 73 74 PetscCall(DMSetMatType(da, MATAIJ)); 75 PetscCall(DMCreateMatrix(da, &J)); 76 77 /* Create nonlinear solver context */ 78 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 79 PetscCall(SNESSetDM(snes, da)); 80 81 /* Set function evaluation and Jacobian evaluation routines */ 82 PetscCall(SNESSetFunction(snes, r, FormGradient, NULL)); 83 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL)); 84 85 PetscCall(SNESSetComputeApplicationContext(snes, (PetscErrorCode (*)(SNES, void **))FormBoundaryConditions, (PetscCtxDestroyFn *)DestroyBoundaryConditions)); 86 87 PetscCall(SNESSetComputeInitialGuess(snes, ComputeInitialGuess, NULL)); 88 89 PetscCall(SNESVISetComputeVariableBounds(snes, FormBounds)); 90 91 PetscCall(SNESSetFromOptions(snes)); 92 93 /* Solve the application */ 94 PetscCall(SNESSolve(snes, NULL, x)); 95 96 /* Free memory */ 97 PetscCall(VecDestroy(&x)); 98 PetscCall(VecDestroy(&r)); 99 PetscCall(MatDestroy(&J)); 100 PetscCall(SNESDestroy(&snes)); 101 102 /* Free user-created data structures */ 103 PetscCall(DMDestroy(&da)); 104 105 PetscCall(PetscFinalize()); 106 return 0; 107 } 108 109 /* -------------------------------------------------------------------- */ 110 111 /* FormBounds - sets the upper and lower bounds 112 113 Input Parameters: 114 . snes - the SNES context 115 116 Output Parameters: 117 . xl - lower bounds 118 . xu - upper bounds 119 */ 120 PetscErrorCode FormBounds(SNES snes, Vec xl, Vec xu) 121 { 122 AppCtx *ctx; 123 124 PetscFunctionBeginUser; 125 PetscCall(SNESGetApplicationContext(snes, &ctx)); 126 PetscCall(VecSet(xl, ctx->lb)); 127 PetscCall(VecSet(xu, ctx->ub)); 128 PetscFunctionReturn(PETSC_SUCCESS); 129 } 130 131 /* -------------------------------------------------------------------- */ 132 133 /* FormGradient - Evaluates gradient of f. 134 135 Input Parameters: 136 . snes - the SNES context 137 . X - input vector 138 . ptr - optional user-defined context, as set by SNESSetFunction() 139 140 Output Parameters: 141 . G - vector containing the newly evaluated gradient 142 */ 143 PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr) 144 { 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(PETSC_SUCCESS); 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 { 272 AppCtx *user; 273 PetscInt i, j, k; 274 PetscInt mx, my; 275 MatStencil row, col[7]; 276 PetscScalar hx, hy, hydhx, hxdhy; 277 PetscScalar f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb; 278 PetscScalar hl, hr, ht, hb, hc, htl, hbr; 279 PetscScalar **x, v[7]; 280 PetscBool assembled; 281 PetscInt xs, xm, ys, ym; 282 Vec localX; 283 DM da; 284 285 PetscFunctionBeginUser; 286 PetscCall(SNESGetDM(snes, &da)); 287 PetscCall(SNESGetApplicationContext(snes, &user)); 288 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)); 289 hx = 1.0 / (mx + 1); 290 hy = 1.0 / (my + 1); 291 hydhx = hy / hx; 292 hxdhy = hx / hy; 293 294 /* Set various matrix options */ 295 PetscCall(MatAssembled(H, &assembled)); 296 if (assembled) PetscCall(MatZeroEntries(H)); 297 298 /* Get local vector */ 299 PetscCall(DMGetLocalVector(da, &localX)); 300 /* Get ghost points */ 301 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 302 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 303 304 /* Get pointers to vector data */ 305 PetscCall(DMDAVecGetArray(da, localX, &x)); 306 307 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 308 /* Compute Jacobian over the locally owned part of the mesh */ 309 for (j = ys; j < ys + ym; j++) { 310 for (i = xs; i < xs + xm; i++) { 311 xc = x[j][i]; 312 xlt = xrb = xl = xr = xb = xt = xc; 313 314 /* Left */ 315 if (i == 0) { 316 xl = user->left[j + 1]; 317 xlt = user->left[j + 2]; 318 } else xl = x[j][i - 1]; 319 320 /* Bottom */ 321 if (j == 0) { 322 xb = user->bottom[i + 1]; 323 xrb = user->bottom[i + 2]; 324 } else xb = x[j - 1][i]; 325 326 /* Right */ 327 if (i + 1 == mx) { 328 xr = user->right[j + 1]; 329 xrb = user->right[j]; 330 } else xr = x[j][i + 1]; 331 332 /* Top */ 333 if (j + 1 == my) { 334 xt = user->top[i + 1]; 335 xlt = user->top[i]; 336 } else xt = x[j + 1][i]; 337 338 /* Top left */ 339 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; 340 341 /* Bottom right */ 342 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; 343 344 d1 = (xc - xl) / hx; 345 d2 = (xc - xr) / hx; 346 d3 = (xc - xt) / hy; 347 d4 = (xc - xb) / hy; 348 d5 = (xrb - xr) / hy; 349 d6 = (xrb - xb) / hx; 350 d7 = (xlt - xl) / hy; 351 d8 = (xlt - xt) / hx; 352 353 f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7); 354 f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4); 355 f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8); 356 f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2); 357 f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5); 358 f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6); 359 360 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2); 361 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4); 362 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4); 363 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2); 364 365 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6); 366 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3); 367 368 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); 369 370 hl /= 2.0; 371 hr /= 2.0; 372 ht /= 2.0; 373 hb /= 2.0; 374 hbr /= 2.0; 375 htl /= 2.0; 376 hc /= 2.0; 377 378 k = 0; 379 row.i = i; 380 row.j = j; 381 /* Bottom */ 382 if (j > 0) { 383 v[k] = hb; 384 col[k].i = i; 385 col[k].j = j - 1; 386 k++; 387 } 388 389 /* Bottom right */ 390 if (j > 0 && i < mx - 1) { 391 v[k] = hbr; 392 col[k].i = i + 1; 393 col[k].j = j - 1; 394 k++; 395 } 396 397 /* left */ 398 if (i > 0) { 399 v[k] = hl; 400 col[k].i = i - 1; 401 col[k].j = j; 402 k++; 403 } 404 405 /* Centre */ 406 v[k] = hc; 407 col[k].i = row.i; 408 col[k].j = row.j; 409 k++; 410 411 /* Right */ 412 if (i < mx - 1) { 413 v[k] = hr; 414 col[k].i = i + 1; 415 col[k].j = j; 416 k++; 417 } 418 419 /* Top left */ 420 if (i > 0 && j < my - 1) { 421 v[k] = htl; 422 col[k].i = i - 1; 423 col[k].j = j + 1; 424 k++; 425 } 426 427 /* Top */ 428 if (j < my - 1) { 429 v[k] = ht; 430 col[k].i = i; 431 col[k].j = j + 1; 432 k++; 433 } 434 435 PetscCall(MatSetValuesStencil(H, 1, &row, k, col, v, INSERT_VALUES)); 436 } 437 } 438 439 /* Assemble the matrix */ 440 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 441 PetscCall(DMDAVecRestoreArray(da, localX, &x)); 442 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 443 PetscCall(DMRestoreLocalVector(da, &localX)); 444 445 PetscCall(PetscLogFlops(199.0 * mx * my)); 446 PetscFunctionReturn(PETSC_SUCCESS); 447 } 448 449 /* ------------------------------------------------------------------- */ 450 /* 451 FormBoundaryConditions - Calculates the boundary conditions for 452 the region. 453 454 Input Parameter: 455 . user - user-defined application context 456 457 Output Parameter: 458 . user - user-defined application context 459 */ 460 PetscErrorCode FormBoundaryConditions(SNES snes, void **inctx) 461 { 462 PetscInt i, j, k, limit = 0, maxits = 5; 463 PetscInt mx, my; 464 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0; 465 PetscScalar one = 1.0, two = 2.0, three = 3.0; 466 PetscScalar det, hx, hy, xt = 0, yt = 0; 467 PetscReal fnorm, tol = 1e-10; 468 PetscScalar u1, u2, nf1, nf2, njac11, njac12, njac21, njac22; 469 PetscScalar b = -0.5, t = 0.5, l = -0.5, r = 0.5; 470 PetscScalar *boundary; 471 AppCtx *user, **ouser = (AppCtx **)inctx; 472 DM da; 473 474 PetscFunctionBeginUser; 475 PetscCall(SNESGetDM(snes, &da)); 476 PetscCall(PetscNew(&user)); 477 *ouser = user; 478 user->lb = .05; 479 user->ub = PETSC_INFINITY; 480 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)); 481 482 /* Check if lower and upper bounds are set */ 483 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lb", &user->lb, 0)); 484 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ub", &user->ub, 0)); 485 bsize = mx + 2; 486 lsize = my + 2; 487 rsize = my + 2; 488 tsize = mx + 2; 489 490 PetscCall(PetscMalloc1(bsize, &user->bottom)); 491 PetscCall(PetscMalloc1(tsize, &user->top)); 492 PetscCall(PetscMalloc1(lsize, &user->left)); 493 PetscCall(PetscMalloc1(rsize, &user->right)); 494 495 hx = (r - l) / (mx + 1.0); 496 hy = (t - b) / (my + 1.0); 497 498 for (j = 0; j < 4; j++) { 499 if (j == 0) { 500 yt = b; 501 xt = l; 502 limit = bsize; 503 boundary = user->bottom; 504 } else if (j == 1) { 505 yt = t; 506 xt = l; 507 limit = tsize; 508 boundary = user->top; 509 } else if (j == 2) { 510 yt = b; 511 xt = l; 512 limit = lsize; 513 boundary = user->left; 514 } else { /* if (j==3) */ 515 yt = b; 516 xt = r; 517 limit = rsize; 518 boundary = user->right; 519 } 520 521 for (i = 0; i < limit; i++) { 522 u1 = xt; 523 u2 = -yt; 524 for (k = 0; k < maxits; k++) { 525 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt; 526 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt; 527 fnorm = PetscRealPart(PetscSqrtScalar(nf1 * nf1 + nf2 * nf2)); 528 if (fnorm <= tol) break; 529 njac11 = one + u2 * u2 - u1 * u1; 530 njac12 = two * u1 * u2; 531 njac21 = -two * u1 * u2; 532 njac22 = -one - u1 * u1 + u2 * u2; 533 det = njac11 * njac22 - njac21 * njac12; 534 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det; 535 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det; 536 } 537 538 boundary[i] = u1 * u1 - u2 * u2; 539 if (j == 0 || j == 1) xt = xt + hx; 540 else yt = yt + hy; /* if (j==2 || j==3) */ 541 } 542 } 543 PetscFunctionReturn(PETSC_SUCCESS); 544 } 545 546 PetscErrorCode DestroyBoundaryConditions(void **ouser) 547 { 548 AppCtx *user = (AppCtx *)*ouser; 549 550 PetscFunctionBeginUser; 551 PetscCall(PetscFree(user->bottom)); 552 PetscCall(PetscFree(user->top)); 553 PetscCall(PetscFree(user->left)); 554 PetscCall(PetscFree(user->right)); 555 PetscCall(PetscFree(*ouser)); 556 PetscFunctionReturn(PETSC_SUCCESS); 557 } 558 559 /* ------------------------------------------------------------------- */ 560 /* 561 ComputeInitialGuess - Calculates the initial guess 562 563 Input Parameters: 564 . user - user-defined application context 565 . X - vector for initial guess 566 567 Output Parameters: 568 . X - newly computed initial guess 569 */ 570 PetscErrorCode ComputeInitialGuess(SNES snes, Vec X, void *dummy) 571 { 572 PetscInt i, j, mx, my; 573 DM da; 574 AppCtx *user; 575 PetscScalar **x; 576 PetscInt xs, xm, ys, ym; 577 578 PetscFunctionBeginUser; 579 PetscCall(SNESGetDM(snes, &da)); 580 PetscCall(SNESGetApplicationContext(snes, &user)); 581 582 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 583 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)); 584 585 /* Get pointers to vector data */ 586 PetscCall(DMDAVecGetArray(da, X, &x)); 587 /* Perform local computations */ 588 for (j = ys; j < ys + ym; j++) { 589 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; 590 } 591 /* Restore vectors */ 592 PetscCall(DMDAVecRestoreArray(da, X, &x)); 593 PetscFunctionReturn(PETSC_SUCCESS); 594 } 595 596 /*TEST 597 598 test: 599 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 600 requires: !single 601 602 test: 603 suffix: 2 604 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 605 requires: !single 606 607 test: 608 suffix: 3 609 args: -snes_type vinewtonrsls -snes_linesearch_type bisection -snes_linesearch_monitor -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 610 requires: !single 611 612 test: 613 suffix: 4 614 args: -snes_type vinewtonssls -snes_linesearch_type bisection -snes_linesearch_monitor -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 615 requires: !single 616 617 test: 618 suffix: 5 619 args: -snes_type vinewtonssls -snes_linesearch_type cp -snes_linesearch_order 3 -snes_linesearch_monitor -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 620 requires: !single 621 622 TEST*/ 623