1 2 static char help[] = "Tests recovery from domain errors in MatMult() and PCApply()\n\n"; 3 4 /* 5 See src/ksp/ksp/tutorials/ex19.c from which this was copied 6 */ 7 8 #include <petscsnes.h> 9 #include <petscdm.h> 10 #include <petscdmda.h> 11 12 /* 13 User-defined routines and data structures 14 */ 15 typedef struct { 16 PetscScalar u, v, omega, temp; 17 } Field; 18 19 PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field **, Field **, void *); 20 21 typedef struct { 22 PetscReal lidvelocity, prandtl, grashof; /* physical parameters */ 23 PetscBool draw_contours; /* flag - 1 indicates drawing contours */ 24 PetscBool errorindomain; 25 PetscBool errorindomainmf; 26 SNES snes; 27 } AppCtx; 28 29 typedef struct { 30 Mat Jmf; 31 } MatShellCtx; 32 33 extern PetscErrorCode FormInitialGuess(AppCtx *, DM, Vec); 34 extern PetscErrorCode MatMult_MyShell(Mat, Vec, Vec); 35 extern PetscErrorCode MatAssemblyEnd_MyShell(Mat, MatAssemblyType); 36 extern PetscErrorCode PCApply_MyShell(PC, Vec, Vec); 37 extern PetscErrorCode SNESComputeJacobian_MyShell(SNES, Vec, Mat, Mat, void *); 38 39 int main(int argc, char **argv) 40 { 41 AppCtx user; /* user-defined work context */ 42 PetscInt mx, my; 43 MPI_Comm comm; 44 DM da; 45 Vec x; 46 Mat J = NULL, Jmf = NULL; 47 MatShellCtx matshellctx; 48 PetscInt mlocal, nlocal; 49 PC pc; 50 KSP ksp; 51 PetscBool errorinmatmult = PETSC_FALSE, errorinpcapply = PETSC_FALSE, errorinpcsetup = PETSC_FALSE; 52 53 PetscFunctionBeginUser; 54 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 55 PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_matmult", &errorinmatmult, NULL)); 56 PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_pcapply", &errorinpcapply, NULL)); 57 PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_pcsetup", &errorinpcsetup, NULL)); 58 user.errorindomain = PETSC_FALSE; 59 PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_domain", &user.errorindomain, NULL)); 60 user.errorindomainmf = PETSC_FALSE; 61 PetscCall(PetscOptionsGetBool(NULL, NULL, "-error_in_domainmf", &user.errorindomainmf, NULL)); 62 63 comm = PETSC_COMM_WORLD; 64 PetscCall(SNESCreate(comm, &user.snes)); 65 66 /* 67 Create distributed array object to manage parallel grid and vectors 68 for principal unknowns (x) and governing residuals (f) 69 */ 70 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da)); 71 PetscCall(DMSetFromOptions(da)); 72 PetscCall(DMSetUp(da)); 73 PetscCall(SNESSetDM(user.snes, da)); 74 75 PetscCall(DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 76 /* 77 Problem parameters (velocity of lid, prandtl, and grashof numbers) 78 */ 79 user.lidvelocity = 1.0 / (mx * my); 80 user.prandtl = 1.0; 81 user.grashof = 1.0; 82 83 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lidvelocity", &user.lidvelocity, NULL)); 84 PetscCall(PetscOptionsGetReal(NULL, NULL, "-prandtl", &user.prandtl, NULL)); 85 PetscCall(PetscOptionsGetReal(NULL, NULL, "-grashof", &user.grashof, NULL)); 86 PetscCall(PetscOptionsHasName(NULL, NULL, "-contours", &user.draw_contours)); 87 88 PetscCall(DMDASetFieldName(da, 0, "x_velocity")); 89 PetscCall(DMDASetFieldName(da, 1, "y_velocity")); 90 PetscCall(DMDASetFieldName(da, 2, "Omega")); 91 PetscCall(DMDASetFieldName(da, 3, "temperature")); 92 93 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 94 Create user context, set problem data, create vector data structures. 95 Also, compute the initial guess. 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 98 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 Create nonlinear solver context 100 101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 102 PetscCall(DMSetApplicationContext(da, &user)); 103 PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, &user)); 104 105 if (errorinmatmult) { 106 PetscCall(MatCreateSNESMF(user.snes, &Jmf)); 107 PetscCall(MatSetFromOptions(Jmf)); 108 PetscCall(MatGetLocalSize(Jmf, &mlocal, &nlocal)); 109 matshellctx.Jmf = Jmf; 110 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)Jmf), mlocal, nlocal, PETSC_DECIDE, PETSC_DECIDE, &matshellctx, &J)); 111 PetscCall(MatShellSetOperation(J, MATOP_MULT, (void (*)(void))MatMult_MyShell)); 112 PetscCall(MatShellSetOperation(J, MATOP_ASSEMBLY_END, (void (*)(void))MatAssemblyEnd_MyShell)); 113 PetscCall(SNESSetJacobian(user.snes, J, J, MatMFFDComputeJacobian, NULL)); 114 } 115 116 PetscCall(SNESSetFromOptions(user.snes)); 117 PetscCall(PetscPrintf(comm, "lid velocity = %g, prandtl # = %g, grashof # = %g\n", (double)user.lidvelocity, (double)user.prandtl, (double)user.grashof)); 118 119 if (errorinpcapply) { 120 PetscCall(SNESGetKSP(user.snes, &ksp)); 121 PetscCall(KSPGetPC(ksp, &pc)); 122 PetscCall(PCSetType(pc, PCSHELL)); 123 PetscCall(PCShellSetApply(pc, PCApply_MyShell)); 124 } 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Solve the nonlinear system 128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 129 PetscCall(DMCreateGlobalVector(da, &x)); 130 PetscCall(FormInitialGuess(&user, da, x)); 131 132 if (errorinpcsetup) { 133 PetscCall(SNESSetUp(user.snes)); 134 PetscCall(SNESSetJacobian(user.snes, NULL, NULL, SNESComputeJacobian_MyShell, NULL)); 135 } 136 PetscCall(SNESSolve(user.snes, NULL, x)); 137 138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139 Free work space. All PETSc objects should be destroyed when they 140 are no longer needed. 141 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 142 PetscCall(MatDestroy(&J)); 143 PetscCall(MatDestroy(&Jmf)); 144 PetscCall(VecDestroy(&x)); 145 PetscCall(DMDestroy(&da)); 146 PetscCall(SNESDestroy(&user.snes)); 147 PetscCall(PetscFinalize()); 148 return 0; 149 } 150 151 /* 152 FormInitialGuess - Forms initial approximation. 153 154 Input Parameters: 155 user - user-defined application context 156 X - vector 157 158 Output Parameter: 159 X - vector 160 */ 161 PetscErrorCode FormInitialGuess(AppCtx *user, DM da, Vec X) 162 { 163 PetscInt i, j, mx, xs, ys, xm, ym; 164 PetscReal grashof, dx; 165 Field **x; 166 167 PetscFunctionBeginUser; 168 grashof = user->grashof; 169 170 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 171 dx = 1.0 / (mx - 1); 172 173 /* 174 Get local grid boundaries (for 2-dimensional DMDA): 175 xs, ys - starting grid indices (no ghost points) 176 xm, ym - widths of local grid (no ghost points) 177 */ 178 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 179 180 /* 181 Get a pointer to vector data. 182 - For default PETSc vectors, VecGetArray() returns a pointer to 183 the data array. Otherwise, the routine is implementation dependent. 184 - You MUST call VecRestoreArray() when you no longer need access to 185 the array. 186 */ 187 PetscCall(DMDAVecGetArray(da, X, &x)); 188 189 /* 190 Compute initial guess over the locally owned part of the grid 191 Initial condition is motionless fluid and equilibrium temperature 192 */ 193 for (j = ys; j < ys + ym; j++) { 194 for (i = xs; i < xs + xm; i++) { 195 x[j][i].u = 0.0; 196 x[j][i].v = 0.0; 197 x[j][i].omega = 0.0; 198 x[j][i].temp = (grashof > 0) * i * dx; 199 } 200 } 201 202 /* 203 Restore vector 204 */ 205 PetscCall(DMDAVecRestoreArray(da, X, &x)); 206 PetscFunctionReturn(PETSC_SUCCESS); 207 } 208 209 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, void *ptr) 210 { 211 AppCtx *user = (AppCtx *)ptr; 212 PetscInt xints, xinte, yints, yinte, i, j; 213 PetscReal hx, hy, dhx, dhy, hxdhy, hydhx; 214 PetscReal grashof, prandtl, lid; 215 PetscScalar u, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym; 216 static PetscInt fail = 0; 217 218 PetscFunctionBeginUser; 219 if ((fail++ > 7 && user->errorindomainmf) || (fail++ > 36 && user->errorindomain)) { 220 PetscMPIInt rank; 221 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)user->snes), &rank)); 222 if (rank == 0) PetscCall(SNESSetFunctionDomainError(user->snes)); 223 } 224 grashof = user->grashof; 225 prandtl = user->prandtl; 226 lid = user->lidvelocity; 227 228 /* 229 Define mesh intervals ratios for uniform grid. 230 231 Note: FD formulae below are normalized by multiplying through by 232 local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions. 233 234 */ 235 dhx = (PetscReal)(info->mx - 1); 236 dhy = (PetscReal)(info->my - 1); 237 hx = 1.0 / dhx; 238 hy = 1.0 / dhy; 239 hxdhy = hx * dhy; 240 hydhx = hy * dhx; 241 242 xints = info->xs; 243 xinte = info->xs + info->xm; 244 yints = info->ys; 245 yinte = info->ys + info->ym; 246 247 /* Test whether we are on the bottom edge of the global array */ 248 if (yints == 0) { 249 j = 0; 250 yints = yints + 1; 251 /* bottom edge */ 252 for (i = info->xs; i < info->xs + info->xm; i++) { 253 f[j][i].u = x[j][i].u; 254 f[j][i].v = x[j][i].v; 255 f[j][i].omega = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy; 256 f[j][i].temp = x[j][i].temp - x[j + 1][i].temp; 257 } 258 } 259 260 /* Test whether we are on the top edge of the global array */ 261 if (yinte == info->my) { 262 j = info->my - 1; 263 yinte = yinte - 1; 264 /* top edge */ 265 for (i = info->xs; i < info->xs + info->xm; i++) { 266 f[j][i].u = x[j][i].u - lid; 267 f[j][i].v = x[j][i].v; 268 f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy; 269 f[j][i].temp = x[j][i].temp - x[j - 1][i].temp; 270 } 271 } 272 273 /* Test whether we are on the left edge of the global array */ 274 if (xints == 0) { 275 i = 0; 276 xints = xints + 1; 277 /* left edge */ 278 for (j = info->ys; j < info->ys + info->ym; j++) { 279 f[j][i].u = x[j][i].u; 280 f[j][i].v = x[j][i].v; 281 f[j][i].omega = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx; 282 f[j][i].temp = x[j][i].temp; 283 } 284 } 285 286 /* Test whether we are on the right edge of the global array */ 287 if (xinte == info->mx) { 288 i = info->mx - 1; 289 xinte = xinte - 1; 290 /* right edge */ 291 for (j = info->ys; j < info->ys + info->ym; j++) { 292 f[j][i].u = x[j][i].u; 293 f[j][i].v = x[j][i].v; 294 f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx; 295 f[j][i].temp = x[j][i].temp - (PetscReal)(grashof > 0); 296 } 297 } 298 299 /* Compute over the interior points */ 300 for (j = yints; j < yinte; j++) { 301 for (i = xints; i < xinte; i++) { 302 /* 303 convective coefficients for upwinding 304 */ 305 vx = x[j][i].u; 306 avx = PetscAbsScalar(vx); 307 vxp = .5 * (vx + avx); 308 vxm = .5 * (vx - avx); 309 vy = x[j][i].v; 310 avy = PetscAbsScalar(vy); 311 vyp = .5 * (vy + avy); 312 vym = .5 * (vy - avy); 313 314 /* U velocity */ 315 u = x[j][i].u; 316 uxx = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx; 317 uyy = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy; 318 f[j][i].u = uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx; 319 320 /* V velocity */ 321 u = x[j][i].v; 322 uxx = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx; 323 uyy = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy; 324 f[j][i].v = uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy; 325 326 /* Omega */ 327 u = x[j][i].omega; 328 uxx = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx; 329 uyy = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy; 330 f[j][i].omega = uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (vyp * (u - x[j - 1][i].omega) + vym * (x[j + 1][i].omega - u)) * hx - .5 * grashof * (x[j][i + 1].temp - x[j][i - 1].temp) * hy; 331 332 /* Temperature */ 333 u = x[j][i].temp; 334 uxx = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx; 335 uyy = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy; 336 f[j][i].temp = uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) * hy + (vyp * (u - x[j - 1][i].temp) + vym * (x[j + 1][i].temp - u)) * hx); 337 } 338 } 339 340 /* 341 Flop count (multiply-adds are counted as 2 operations) 342 */ 343 PetscCall(PetscLogFlops(84.0 * info->ym * info->xm)); 344 PetscFunctionReturn(PETSC_SUCCESS); 345 } 346 347 PetscErrorCode MatMult_MyShell(Mat A, Vec x, Vec y) 348 { 349 MatShellCtx *matshellctx; 350 static PetscInt fail = 0; 351 352 PetscFunctionBegin; 353 PetscCall(MatShellGetContext(A, &matshellctx)); 354 PetscCall(MatMult(matshellctx->Jmf, x, y)); 355 if (fail++ > 5) { 356 PetscMPIInt rank; 357 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank)); 358 if (rank == 0) PetscCall(VecSetInf(y)); 359 } 360 PetscFunctionReturn(PETSC_SUCCESS); 361 } 362 363 PetscErrorCode MatAssemblyEnd_MyShell(Mat A, MatAssemblyType tp) 364 { 365 MatShellCtx *matshellctx; 366 367 PetscFunctionBegin; 368 PetscCall(MatShellGetContext(A, &matshellctx)); 369 PetscCall(MatAssemblyEnd(matshellctx->Jmf, tp)); 370 PetscFunctionReturn(PETSC_SUCCESS); 371 } 372 373 PetscErrorCode PCApply_MyShell(PC pc, Vec x, Vec y) 374 { 375 static PetscInt fail = 0; 376 377 PetscFunctionBegin; 378 PetscCall(VecCopy(x, y)); 379 if (fail++ > 3) { 380 PetscMPIInt rank; 381 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 382 if (rank == 0) PetscCall(VecSetInf(y)); 383 } 384 PetscFunctionReturn(PETSC_SUCCESS); 385 } 386 387 PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES, Vec, Mat, Mat, void *); 388 389 PetscErrorCode SNESComputeJacobian_MyShell(SNES snes, Vec X, Mat A, Mat B, void *ctx) 390 { 391 static PetscInt fail = 0; 392 393 PetscFunctionBegin; 394 PetscCall(SNESComputeJacobian_DMDA(snes, X, A, B, ctx)); 395 if (fail++ > 0) PetscCall(MatZeroEntries(A)); 396 PetscFunctionReturn(PETSC_SUCCESS); 397 } 398 399 /*TEST 400 401 test: 402 args: -snes_converged_reason -ksp_converged_reason 403 404 test: 405 suffix: 2 406 args: -snes_converged_reason -ksp_converged_reason -error_in_matmult 407 408 test: 409 suffix: 3 410 args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply 411 412 test: 413 suffix: 4 414 args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup 415 416 test: 417 suffix: 5 418 args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi 419 420 test: 421 suffix: 5_fieldsplit 422 args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type fieldsplit 423 output_file: output/ex69_5.out 424 425 test: 426 suffix: 6 427 args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf -pc_type none 428 429 test: 430 suffix: 7 431 args: -snes_converged_reason -ksp_converged_reason -error_in_domain 432 433 test: 434 suffix: 8 435 args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf -pc_type none 436 437 TEST*/ 438