Actual source code: ex69.c

petsc-3.10.5 2019-03-28
Report Typos and Errors

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

  4: /*
  5:       See src/ksp/ksp/examples/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);


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

153: /* ------------------------------------------------------------------- */


156: /*
157:    FormInitialGuess - Forms initial approximation.

159:    Input Parameters:
160:    user - user-defined application context
161:    X - vector

163:    Output Parameter:
164:    X - vector
165: */
166: PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
167: {
168:   PetscInt       i,j,mx,xs,ys,xm,ym;
170:   PetscReal      grashof,dx;
171:   Field          **x;

174:   grashof = user->grashof;

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

179:   /*
180:      Get local grid boundaries (for 2-dimensional DMDA):
181:        xs, ys   - starting grid indices (no ghost points)
182:        xm, ym   - widths of local grid (no ghost points)
183:   */
184:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

186:   /*
187:      Get a pointer to vector data.
188:        - For default PETSc vectors, VecGetArray() returns a pointer to
189:          the data array.  Otherwise, the routine is implementation dependent.
190:        - You MUST call VecRestoreArray() when you no longer need access to
191:          the array.
192:   */
193:   DMDAVecGetArray(da,X,&x);

195:   /*
196:      Compute initial guess over the locally owned part of the grid
197:      Initial condition is motionless fluid and equilibrium temperature
198:   */
199:   for (j=ys; j<ys+ym; j++) {
200:     for (i=xs; i<xs+xm; i++) {
201:       x[j][i].u     = 0.0;
202:       x[j][i].v     = 0.0;
203:       x[j][i].omega = 0.0;
204:       x[j][i].temp  = (grashof>0)*i*dx;
205:     }
206:   }

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

215: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
216: {
217:   AppCtx          *user = (AppCtx*)ptr;
218:   PetscErrorCode  ierr;
219:   PetscInt        xints,xinte,yints,yinte,i,j;
220:   PetscReal       hx,hy,dhx,dhy,hxdhy,hydhx;
221:   PetscReal       grashof,prandtl,lid;
222:   PetscScalar     u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
223:   static PetscInt fail = 0;

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

237:   /*
238:      Define mesh intervals ratios for uniform grid.

240:      Note: FD formulae below are normalized by multiplying through by
241:      local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.


244:   */
245:   dhx   = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
246:   hx    = 1.0/dhx;                   hy = 1.0/dhy;
247:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

251:   /* Test whether we are on the bottom edge of the global array */
252:   if (yints == 0) {
253:     j     = 0;
254:     yints = yints + 1;
255:     /* bottom edge */
256:     for (i=info->xs; i<info->xs+info->xm; i++) {
257:       f[j][i].u     = x[j][i].u;
258:       f[j][i].v     = x[j][i].v;
259:       f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
260:       f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
261:     }
262:   }

264:   /* Test whether we are on the top edge of the global array */
265:   if (yinte == info->my) {
266:     j     = info->my - 1;
267:     yinte = yinte - 1;
268:     /* top edge */
269:     for (i=info->xs; i<info->xs+info->xm; i++) {
270:       f[j][i].u     = x[j][i].u - lid;
271:       f[j][i].v     = x[j][i].v;
272:       f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
273:       f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
274:     }
275:   }

277:   /* Test whether we are on the left edge of the global array */
278:   if (xints == 0) {
279:     i     = 0;
280:     xints = xints + 1;
281:     /* left edge */
282:     for (j=info->ys; j<info->ys+info->ym; j++) {
283:       f[j][i].u     = x[j][i].u;
284:       f[j][i].v     = x[j][i].v;
285:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
286:       f[j][i].temp  = x[j][i].temp;
287:     }
288:   }

290:   /* Test whether we are on the right edge of the global array */
291:   if (xinte == info->mx) {
292:     i     = info->mx - 1;
293:     xinte = xinte - 1;
294:     /* right edge */
295:     for (j=info->ys; j<info->ys+info->ym; j++) {
296:       f[j][i].u     = x[j][i].u;
297:       f[j][i].v     = x[j][i].v;
298:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
299:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
300:     }
301:   }

303:   /* Compute over the interior points */
304:   for (j=yints; j<yinte; j++) {
305:     for (i=xints; i<xinte; i++) {

307:       /*
308:        convective coefficients for upwinding
309:       */
310:       vx  = x[j][i].u; avx = PetscAbsScalar(vx);
311:       vxp = .5*(vx+avx); vxm = .5*(vx-avx);
312:       vy  = x[j][i].v; avy = PetscAbsScalar(vy);
313:       vyp = .5*(vy+avy); vym = .5*(vy-avy);

315:       /* U velocity */
316:       u         = x[j][i].u;
317:       uxx       = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
318:       uyy       = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
319:       f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

321:       /* V velocity */
322:       u         = x[j][i].v;
323:       uxx       = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
324:       uyy       = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
325:       f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

327:       /* Omega */
328:       u             = x[j][i].omega;
329:       uxx           = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
330:       uyy           = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
331:       f[j][i].omega = uxx + uyy + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u))*hy +
332:                       (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u))*hx -
333:                       .5*grashof*(x[j][i+1].temp - x[j][i-1].temp)*hy;

335:       /* Temperature */
336:       u            = x[j][i].temp;
337:       uxx          = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
338:       uyy          = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
339:       f[j][i].temp =  uxx + uyy  + prandtl*((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u))*hy +
340:                                             (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u))*hx);
341:     }
342:   }

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

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

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

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

374:   MatShellGetContext(A,&matshellctx);
375:   MatAssemblyEnd(matshellctx->Jmf,tp);
376:   return(0);
377: }

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

385:   VecCopy(x,y);
386:   if (fail++ > 3) {
387:     PetscMPIInt rank;
388:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
389:     if (!rank) {VecSetInf(y);}
390:   }
391:   return(0);
392: }

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

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

402:   SNESComputeJacobian_DMDA(snes,X,A,B,ctx);
403:   if (fail++ > 0) {
404:     MatZeroEntries(A);
405:   }
406:   return(0);
407: }


410: /*TEST

412:    test:
413:       args: -snes_converged_reason -ksp_converged_reason

415:    test:
416:       suffix: 2
417:       args: -snes_converged_reason -ksp_converged_reason -error_in_matmult

419:    test:
420:       suffix: 3
421:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcapply

423:    test:
424:       suffix: 4
425:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup

427:    test:
428:       suffix: 5
429:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type bjacobi

431:    test:
432:       suffix: 5_fieldsplit
433:       args: -snes_converged_reason -ksp_converged_reason -error_in_pcsetup -pc_type fieldsplit
434:       output_file: output/ex69_5.out

436:    test:
437:       suffix: 6
438:       args: -snes_converged_reason -ksp_converged_reason -error_in_domainmf -snes_mf

440:    test:
441:       suffix: 7
442:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain

444:    test:
445:       suffix: 8
446:       args: -snes_converged_reason -ksp_converged_reason -error_in_domain -snes_mf
447:       TODO: Need to determine if deprecated

449: TEST*/