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