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);if (ierr) return ierr;
 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 ierr;
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;
168:   PetscReal      grashof,dx;
169:   Field          **x;

172:   grashof = user->grashof;

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

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

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

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

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

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;

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

235:   /*
236:      Define mesh intervals ratios for uniform grid.

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.

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;

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

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

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

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

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

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

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

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;

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;

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;

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

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

348: PetscErrorCode MatMult_MyShell(Mat A,Vec x,Vec y)
349: {
350:   PetscErrorCode  ierr;
351:   MatShellCtx     *matshellctx;
352:   static PetscInt fail = 0;

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

365: PetscErrorCode MatAssemblyEnd_MyShell(Mat A,MatAssemblyType tp)
366: {
368:   MatShellCtx    *matshellctx;

371:   MatShellGetContext(A,&matshellctx);
372:   MatAssemblyEnd(matshellctx->Jmf,tp);
373:   return(0);
374: }

376: PetscErrorCode PCApply_MyShell(PC pc,Vec x,Vec y)
377: {
379:   static PetscInt fail = 0;

382:   VecCopy(x,y);
383:   if (fail++ > 3) {
384:     PetscMPIInt rank;
385:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
386:     if (rank == 0) {VecSetInf(y);}
387:   }
388:   return(0);
389: }

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

393: PetscErrorCode SNESComputeJacobian_MyShell(SNES snes,Vec X,Mat A,Mat B,void *ctx)
394: {
395:   static PetscInt fail = 0;
396:   PetscErrorCode  ierr;

399:   SNESComputeJacobian_DMDA(snes,X,A,B,ctx);
400:   if (fail++ > 0) {
401:     MatZeroEntries(A);
402:   }
403:   return(0);
404: }

406: /*TEST

408:    test:
409:       args: -snes_converged_reason -ksp_converged_reason

411:    test:
412:       suffix: 2
413:       args: -snes_converged_reason -ksp_converged_reason -error_in_matmult

415:    test:
416:       suffix: 3
417:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply

419:    test:
420:       suffix: 4
421:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup

423:    test:
424:       suffix: 5
425:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi

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

432:    test:
433:       suffix: 6
434:       args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf -pc_type none

436:    test:
437:       suffix: 7
438:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain

440:    test:
441:       suffix: 8
442:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf -pc_type none

444: TEST*/