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