Actual source code: ex69.c


  2: static char help[] = "Tests recovery from domain errors in MatMult() and PCApply()\n\n";

  4: /*
  5:       See src/ksp/ksp/tutorials/ex19.c from which this was copied
  6: */

  8: #include <petscsnes.h>
  9: #include <petscdm.h>
 10: #include <petscdmda.h>

 12: /*
 13:    User-defined routines and data structures
 14: */
 15: typedef struct {
 16:   PetscScalar u,v,omega,temp;
 17: } Field;

 19: PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);

 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;

 29: typedef struct {
 30:   Mat Jmf;
 31: } MatShellCtx;

 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*);

 39: int main(int argc,char **argv)
 40: {
 41:   AppCtx         user;                /* user-defined work context */
 42:   PetscInt       mx,my;
 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;

 54:   PetscInitialize(&argc,&argv,(char*)0,help);
 55:   PetscOptionsGetBool(NULL,NULL,"-error_in_matmult",&errorinmatmult,NULL);
 56:   PetscOptionsGetBool(NULL,NULL,"-error_in_pcapply",&errorinpcapply,NULL);
 57:   PetscOptionsGetBool(NULL,NULL,"-error_in_pcsetup",&errorinpcsetup,NULL);
 58:   user.errorindomain = PETSC_FALSE;
 59:   PetscOptionsGetBool(NULL,NULL,"-error_in_domain",&user.errorindomain,NULL);
 60:   user.errorindomainmf = PETSC_FALSE;
 61:   PetscOptionsGetBool(NULL,NULL,"-error_in_domainmf",&user.errorindomainmf,NULL);

 63:   comm = PETSC_COMM_WORLD;
 64:   SNESCreate(comm,&user.snes);

 66:   /*
 67:       Create distributed array object to manage parallel grid and vectors
 68:       for principal unknowns (x) and governing residuals (f)
 69:   */
 70:   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:   DMSetFromOptions(da);
 72:   DMSetUp(da);
 73:   SNESSetDM(user.snes,da);

 75:   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;

 84:   PetscOptionsGetReal(NULL,NULL,"-lidvelocity",&user.lidvelocity,NULL);
 85:   PetscOptionsGetReal(NULL,NULL,"-prandtl",&user.prandtl,NULL);
 86:   PetscOptionsGetReal(NULL,NULL,"-grashof",&user.grashof,NULL);
 87:   PetscOptionsHasName(NULL,NULL,"-contours",&user.draw_contours);

 89:   DMDASetFieldName(da,0,"x_velocity");
 90:   DMDASetFieldName(da,1,"y_velocity");
 91:   DMDASetFieldName(da,2,"Omega");
 92:   DMDASetFieldName(da,3,"temperature");

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:      Create user context, set problem data, create vector data structures.
 96:      Also, compute the initial guess.
 97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100:      Create nonlinear solver context

102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   DMSetApplicationContext(da,&user);
104:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);

106:   if (errorinmatmult) {
107:     MatCreateSNESMF(user.snes,&Jmf);
108:     MatSetFromOptions(Jmf);
109:     MatGetLocalSize(Jmf,&mlocal,&nlocal);
110:     matshellctx.Jmf = Jmf;
111:     MatCreateShell(PetscObjectComm((PetscObject)Jmf),mlocal,nlocal,PETSC_DECIDE,PETSC_DECIDE,&matshellctx,&J);
112:     MatShellSetOperation(J,MATOP_MULT,(void (*)(void))MatMult_MyShell);
113:     MatShellSetOperation(J,MATOP_ASSEMBLY_END,(void (*)(void))MatAssemblyEnd_MyShell);
114:     SNESSetJacobian(user.snes,J,J,MatMFFDComputeJacobian,NULL);
115:   }

117:   SNESSetFromOptions(user.snes);
118:   PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %g\n",(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof);

120:   if (errorinpcapply) {
121:     SNESGetKSP(user.snes,&ksp);
122:     KSPGetPC(ksp,&pc);
123:     PCSetType(pc,PCSHELL);
124:     PCShellSetApply(pc,PCApply_MyShell);
125:   }

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Solve the nonlinear system
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   DMCreateGlobalVector(da,&x);
131:   FormInitialGuess(&user,da,x);

133:   if (errorinpcsetup) {
134:     SNESSetUp(user.snes);
135:     SNESSetJacobian(user.snes,NULL,NULL,SNESComputeJacobian_MyShell,NULL);
136:   }
137:   SNESSolve(user.snes,NULL,x);

139:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140:      Free work space.  All PETSc objects should be destroyed when they
141:      are no longer needed.
142:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143:   MatDestroy(&J);
144:   MatDestroy(&Jmf);
145:   VecDestroy(&x);
146:   DMDestroy(&da);
147:   SNESDestroy(&user.snes);
148:   PetscFinalize();
149:   return 0;
150: }

152: /* ------------------------------------------------------------------- */

154: /*
155:    FormInitialGuess - Forms initial approximation.

157:    Input Parameters:
158:    user - user-defined application context
159:    X - vector

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;

171:   grashof = user->grashof;

173:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
174:   dx   = 1.0/(mx-1);

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:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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:   DMDAVecGetArray(da,X,&x);

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:   }

205:   /*
206:      Restore vector
207:   */
208:   DMDAVecRestoreArray(da,X,&x);
209:   return 0;
210: }

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;

222:   if ((fail++ > 7 && user->errorindomainmf) || (fail++ > 36 && user->errorindomain)) {
223:     PetscMPIInt rank;
224:     MPI_Comm_rank(PetscObjectComm((PetscObject)user->snes),&rank);
225:     if (rank == 0) {
226:       SNESSetFunctionDomainError(user->snes);
227:     }
228:   }
229:   grashof = user->grashof;
230:   prandtl = user->prandtl;
231:   lid     = user->lidvelocity;

233:   /*
234:      Define mesh intervals ratios for uniform grid.

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.

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;

244:   xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;

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:   }

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:   }

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:   }

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:   }

298:   /* Compute over the interior points */
299:   for (j=yints; j<yinte; j++) {
300:     for (i=xints; i<xinte; i++) {

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);

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;

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;

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;

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:   }

339:   /*
340:      Flop count (multiply-adds are counted as 2 operations)
341:   */
342:   PetscLogFlops(84.0*info->ym*info->xm);
343:   return 0;
344: }

346: PetscErrorCode MatMult_MyShell(Mat A,Vec x,Vec y)
347: {
348:   MatShellCtx     *matshellctx;
349:   static PetscInt fail = 0;

351:   MatShellGetContext(A,&matshellctx);
352:   MatMult(matshellctx->Jmf,x,y);
353:   if (fail++ > 5) {
354:     PetscMPIInt rank;
355:     MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
356:     if (rank == 0) VecSetInf(y);
357:   }
358:   return 0;
359: }

361: PetscErrorCode MatAssemblyEnd_MyShell(Mat A,MatAssemblyType tp)
362: {
363:   MatShellCtx    *matshellctx;

365:   MatShellGetContext(A,&matshellctx);
366:   MatAssemblyEnd(matshellctx->Jmf,tp);
367:   return 0;
368: }

370: PetscErrorCode PCApply_MyShell(PC pc,Vec x,Vec y)
371: {
372:   static PetscInt fail = 0;

374:   VecCopy(x,y);
375:   if (fail++ > 3) {
376:     PetscMPIInt rank;
377:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
378:     if (rank == 0) VecSetInf(y);
379:   }
380:   return 0;
381: }

383: PETSC_EXTERN PetscErrorCode SNESComputeJacobian_DMDA(SNES,Vec,Mat,Mat,void*);

385: PetscErrorCode SNESComputeJacobian_MyShell(SNES snes,Vec X,Mat A,Mat B,void *ctx)
386: {
387:   static PetscInt fail = 0;

389:   SNESComputeJacobian_DMDA(snes,X,A,B,ctx);
390:   if (fail++ > 0) {
391:     MatZeroEntries(A);
392:   }
393:   return 0;
394: }

396: /*TEST

398:    test:
399:       args: -snes_converged_reason -ksp_converged_reason

401:    test:
402:       suffix: 2
403:       args: -snes_converged_reason -ksp_converged_reason -error_in_matmult

405:    test:
406:       suffix: 3
407:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply

409:    test:
410:       suffix: 4
411:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup

413:    test:
414:       suffix: 5
415:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi

417:    test:
418:       suffix: 5_fieldsplit
419:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type fieldsplit
420:       output_file: output/ex69_5.out

422:    test:
423:       suffix: 6
424:       args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf -pc_type none

426:    test:
427:       suffix: 7
428:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain

430:    test:
431:       suffix: 8
432:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf -pc_type none

434: TEST*/