Actual source code: ex31.c

  2: static char help[] = "Model multi-physics solver\n\n";

  4: /*
  5:      A model "multi-physics" solver based on the Vincent Mousseau's reactor core pilot code.

  7:      There are three grids:

  9:             --------------------- DMDA1        

 11:                     nyv  -->       --------------------- DMDA2
 12:                                    |                    | 
 13:                                    |                    | 
 14:                                    |                    |   
 15:                                    |                    | 
 16:                     nyvf-1 -->     |                    |         --------------------- DMDA3
 17:                                    |                    |        |                    |
 18:                                    |                    |        |                    |
 19:                                    |                    |        |                    |
 20:                                    |                    |        |                    |
 21:                          0 -->     ---------------------          ---------------------

 23:                    nxv                     nxv                          nxv

 25:             Fluid                     Thermal conduction              Fission (core)
 26:                                     (cladding and core)

 28:     Notes:
 29:      * The discretization approach used is to have ghost nodes OUTSIDE the physical domain
 30:       that are used to apply the stencil near the boundary; in order to implement this with
 31:       PETSc DMDAs we simply define the DMDAs to have periodic boundary conditions and use those
 32:       periodic ghost points to store the needed extra variables (which do not equations associated
 33:       with them). Note that these periodic ghost nodes have NOTHING to do with the ghost nodes
 34:       used for parallel computing.

 36:    This can be run in two modes:

 38:         -snes_mf -pc_type shell * for matrix free with "physics based" preconditioning or
 39:         -snes_mf *                for matrix free with an explicit matrix based preconditioner where the explicit
 40:                                   matrix is computed via finite differences igoring the coupling between the models or
 41:                    * for "approximate Newton" where the Jacobian is not used but rather the approximate Jacobian as
 42:                    computed in the alternative above.

 44: */

 46: #include <petscdmmg.h>
 47: #include <petscdmcomposite.h>

 49: typedef struct {
 50:   PetscScalar pri,ugi,ufi,agi,vgi,vfi;              /* initial conditions for fluid variables */
 51:   PetscScalar prin,ugin,ufin,agin,vgin,vfin;        /* inflow boundary conditions for fluid */
 52:   PetscScalar prout,ugout,ufout,agout,vgout;        /* outflow boundary conditions for fluid */

 54:   PetscScalar twi;                                  /* initial condition for tempature */

 56:   PetscScalar phii;                                 /* initial conditions for fuel */
 57:   PetscScalar prei;

 59:   PetscInt    nxv,nyv,nyvf;

 61:   MPI_Comm    comm;

 63:   DM          pack;

 65:   DMMG        *fdmmg;                              /* used by PCShell to solve diffusion problem */
 66:   Vec         dx,dy;
 67:   Vec         c;
 68: } AppCtx;

 70: typedef struct {                 /* Fluid unknowns */
 71:   PetscScalar prss;
 72:   PetscScalar ergg;
 73:   PetscScalar ergf;
 74:   PetscScalar alfg;
 75:   PetscScalar velg;
 76:   PetscScalar velf;
 77: } FluidField;

 79: typedef struct {                 /* Fuel unknowns */
 80:   PetscScalar phii;
 81:   PetscScalar prei;
 82: } FuelField;


 92: int main(int argc,char **argv)
 93: {
 94:   DMMG           *dmmg;               /* multilevel grid structure */
 96:   DM             da;
 97:   AppCtx         app;
 98:   PC             pc;
 99:   KSP            ksp;
100:   PetscBool      isshell;
101:   PetscViewer    v1;

103:   PetscInitialize(&argc,&argv,(char *)0,help);

105:   PetscPreLoadBegin(PETSC_TRUE,"SetUp");

107:     app.comm = PETSC_COMM_WORLD;
108:     app.nxv  = 6;
109:     app.nyvf = 3;
110:     app.nyv  = app.nyvf + 2;
111:     PetscOptionsBegin(app.comm,PETSC_NULL,"Options for Grid Sizes",PETSC_NULL);
112:       PetscOptionsInt("-nxv","Grid spacing in X direction",PETSC_NULL,app.nxv,&app.nxv,PETSC_NULL);
113:       PetscOptionsInt("-nyvf","Grid spacing in Y direction of Fuel",PETSC_NULL,app.nyvf,&app.nyvf,PETSC_NULL);
114:       PetscOptionsInt("-nyv","Total Grid spacing in Y direction of",PETSC_NULL,app.nyv,&app.nyv,PETSC_NULL);
115:     PetscOptionsEnd();

117:     PetscViewerDrawOpen(app.comm,PETSC_NULL,"",-1,-1,-1,-1,&v1);

119:     /*
120:        Create the DMComposite object to manage the three grids/physics. 
121:        We use a 1d decomposition along the y direction (since one of the grids is 1d).

123:     */
124:     DMCompositeCreate(app.comm,&app.pack);

126:     /* 6 fluid unknowns, 3 ghost points on each end for either periodicity or simply boundary conditions */
127:     DMDACreate1d(app.comm,DMDA_BOUNDARY_PERIODIC,app.nxv,6,3,0,&da);
128:     DMDASetFieldName(da,0,"prss");
129:     DMDASetFieldName(da,1,"ergg");
130:     DMDASetFieldName(da,2,"ergf");
131:     DMDASetFieldName(da,3,"alfg");
132:     DMDASetFieldName(da,4,"velg");
133:     DMDASetFieldName(da,5,"velf");
134:     DMCompositeAddDM(app.pack,(DM)da);
135:     DMDestroy(&da);

137:     DMDACreate2d(app.comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,app.nxv,app.nyv,PETSC_DETERMINE,1,1,1,0,0,&da);
138:     DMDASetFieldName(da,0,"Tempature");
139:     DMCompositeAddDM(app.pack,(DM)da);
140:     DMDestroy(&da);

142:     DMDACreate2d(app.comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,app.nxv,app.nyvf,PETSC_DETERMINE,1,2,1,0,0,&da);
143:     DMDASetFieldName(da,0,"Phi");
144:     DMDASetFieldName(da,1,"Pre");
145:     DMCompositeAddDM(app.pack,(DM)da);
146:     DMDestroy(&da);
147: 
148:     app.pri = 1.0135e+5;
149:     app.ugi = 2.5065e+6;
150:     app.ufi = 4.1894e+5;
151:     app.agi = 1.00e-1;
152:     app.vgi = 1.0e-1 ;
153:     app.vfi = 1.0e-1;

155:     app.prin = 1.0135e+5;
156:     app.ugin = 2.5065e+6;
157:     app.ufin = 4.1894e+5;
158:     app.agin = 1.00e-1;
159:     app.vgin = 1.0e-1 ;
160:     app.vfin = 1.0e-1;

162:     app.prout = 1.0135e+5;
163:     app.ugout = 2.5065e+6;
164:     app.ufout = 4.1894e+5;
165:     app.agout = 3.0e-1;

167:     app.twi = 373.15e+0;

169:     app.phii = 1.0e+0;
170:     app.prei = 1.0e-5;

172:     /*
173:        Create the solver object and attach the grid/physics info 
174:     */
175:     DMMGCreate(app.comm,1,0,&dmmg);
176:     DMMGSetDM(dmmg,(DM)app.pack);
177:     DMMGSetUser(dmmg,0,&app);
178:     DMMGSetISColoringType(dmmg,IS_COLORING_GLOBAL);
179:     CHKMEMQ;


182:     DMMGSetInitialGuess(dmmg,FormInitialGuess);
183:     DMMGSetSNES(dmmg,FormFunction,0);
184:     DMMGSetFromOptions(dmmg);

186:     /* Supply custom shell preconditioner if requested */
187:     SNESGetKSP(DMMGGetSNES(dmmg),&ksp);
188:     KSPGetPC(ksp,&pc);
189:     PetscTypeCompare((PetscObject)pc,PCSHELL,&isshell);
190:     if (isshell) {
191:       PCShellSetContext(pc,&app);
192:       PCShellSetSetUp(pc,MyPCSetUp);
193:       PCShellSetApply(pc,MyPCApply);
194:       PCShellSetDestroy(pc,MyPCDestroy);
195:     }

197:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198:        Solve the nonlinear system
199:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

201:   PetscPreLoadStage("Solve");
202:     DMMGSolve(dmmg);


205:     VecView(DMMGGetx(dmmg),v1);

207:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208:        Free work space.  All PETSc objects should be destroyed when they
209:        are no longer needed.
210:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

212:     PetscViewerDestroy(&v1);
213:     DMDestroy(&app.pack);
214:     DMMGDestroy(dmmg);
215:   PetscPreLoadEnd();

217:   PetscFinalize();
218:   return 0;
219: }

221: /* ------------------------------------------------------------------- */


224: /* 
225:    FormInitialGuessLocal* Forms the initial SOLUTION for the fluid, cladding and fuel

227:  */
230: PetscErrorCode FormInitialGuessLocalFluid(AppCtx *app,DMDALocalInfo *info,FluidField *f)
231: {
232:   PetscInt       i;

235:   for (i=info->xs; i<info->xs+info->xm; i++) {
236:     f[i].prss = app->pri;
237:     f[i].ergg = app->ugi;
238:     f[i].ergf = app->ufi;
239:     f[i].alfg = app->agi;
240:     f[i].velg = app->vgi;
241:     f[i].velf = app->vfi;
242:   }

244:   return(0);
245: }

249: PetscErrorCode FormInitialGuessLocalThermal(AppCtx *app,DMDALocalInfo *info2,PetscScalar **T)
250: {
251:   PetscInt i,j;

254:   for (i=info2->xs; i<info2->xs+info2->xm; i++) {
255:     for (j=info2->ys;j<info2->ys+info2->ym; j++) {
256:       T[j][i] = app->twi;
257:     }
258:   }
259:   return(0);
260: }

264: PetscErrorCode FormInitialGuessLocalFuel(AppCtx *app,DMDALocalInfo *info2,FuelField **F)
265: {
266:   PetscInt i,j;

269:   for (i=info2->xs; i<info2->xs+info2->xm; i++) {
270:     for (j=info2->ys;j<info2->ys+info2->ym; j++) {
271:       F[j][i].phii = app->phii;
272:       F[j][i].prei = app->prei;
273:     }
274:   }
275:   return(0);
276: }

278: /* 
279:    FormFunctionLocal* - Forms user provided function

281: */
284: PetscErrorCode FormFunctionLocalFluid(DMDALocalInfo *info,FluidField *u,FluidField *f)
285: {
286:   PetscInt       i;

289:   for (i=info->xs; i<info->xs+info->xm; i++) {
290:     f[i].prss = u[i].prss;
291:     f[i].ergg = u[i].ergg;
292:     f[i].ergf = u[i].ergf;
293:     f[i].alfg = u[i].alfg;
294:     f[i].velg = u[i].velg;
295:     f[i].velf = u[i].velf;
296:   }
297:   return(0);
298: }

302: PetscErrorCode FormFunctionLocalThermal(DMDALocalInfo *info,PetscScalar **T,PetscScalar **f)
303: {
304:   PetscInt i,j;

307:   for (i=info->xs; i<info->xs+info->xm; i++) {
308:     for (j=info->ys;j<info->ys+info->ym; j++) {
309:       f[j][i] = T[j][i];
310:     }
311:   }
312:   return(0);
313: }

317: PetscErrorCode FormFunctionLocalFuel(DMDALocalInfo *info,FuelField **U,FuelField **F)
318: {
319:   PetscInt i,j;

322:   for (i=info->xs; i<info->xs+info->xm; i++) {
323:     for (j=info->ys;j<info->ys+info->ym; j++) {
324:       F[j][i].phii = U[j][i].phii;
325:       F[j][i].prei = U[j][i].prei;
326:     }
327:   }
328:   return(0);
329: }

331: 
334: /* 
335:    FormInitialGuess  - Unwraps the global solution vector and passes its local pieces into the user function

337:  */
338: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
339: {
340:   DM             dm = (DM)dmmg->dm;
341:   DMDALocalInfo    info1,info2,info3;
342:   DM             da1,da2,da3;
343:   FluidField     *x1;
344:   PetscScalar    **x2;
345:   FuelField      **x3;
346:   Vec            X1,X2,X3;
348:   AppCtx         *app = (AppCtx*)dmmg->user;

351:   DMCompositeGetEntries(dm,&da1,&da2,&da3);
352:   DMDAGetLocalInfo(da1,&info1);
353:   DMDAGetLocalInfo(da2,&info2);
354:   DMDAGetLocalInfo(da3,&info3);

356:   /* Access the three subvectors in X */
357:   DMCompositeGetAccess(dm,X,&X1,&X2,&X3);

359:   /* Access the arrays inside the subvectors of X */
360:   DMDAVecGetArray(da1,X1,(void**)&x1);
361:   DMDAVecGetArray(da2,X2,(void**)&x2);
362:   DMDAVecGetArray(da3,X3,(void**)&x3);

364:   /* Evaluate local user provided function */
365:   FormInitialGuessLocalFluid(app,&info1,x1);
366:   FormInitialGuessLocalThermal(app,&info2,x2);
367:   FormInitialGuessLocalFuel(app,&info3,x3);

369:   DMDAVecRestoreArray(da1,X1,(void**)&x1);
370:   DMDAVecRestoreArray(da2,X2,(void**)&x2);
371:   DMDAVecRestoreArray(da3,X3,(void**)&x3);
372:   DMCompositeRestoreAccess(dm,X,&X1,&X2,&X3);
373:   return(0);
374: }

378: /* 
379:    FormFunction  - Unwraps the input vector and passes its local ghosted pieces into the user function

381:  */
382: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
383: {
384:   DMMG           dmmg = (DMMG)ctx;
385:   DM             dm = (DM)dmmg->dm;
386:   DMDALocalInfo    info1,info2,info3;
387:   DM             da1,da2,da3;
388:   FluidField     *x1,*f1;
389:   PetscScalar    **x2,**f2;
390:   FuelField      **x3,**f3;
391:   Vec            X1,X2,X3,F1,F2,F3;
393:   PetscInt       i,j;
394:   AppCtx         *app = (AppCtx*)dmmg->user;

397:   DMCompositeGetEntries(dm,&da1,&da2,&da3);
398:   DMDAGetLocalInfo(da1,&info1);
399:   DMDAGetLocalInfo(da2,&info2);
400:   DMDAGetLocalInfo(da3,&info3);

402:   /* Get local vectors to hold ghosted parts of X; then fill in the ghosted vectors from the unghosted global vector X */
403:   DMCompositeGetLocalVectors(dm,&X1,&X2,&X3);
404:   DMCompositeScatter(dm,X,X1,X2,X3);

406:   /* Access the arrays inside the subvectors of X */
407:   DMDAVecGetArray(da1,X1,(void**)&x1);
408:   DMDAVecGetArray(da2,X2,(void**)&x2);
409:   DMDAVecGetArray(da3,X3,(void**)&x3);

411:    /*
412:     Ghost points for periodicity are used to "force" inflow/outflow fluid boundary conditions 
413:     Note that there is no periodicity; we define periodicity to "cheat" and have ghost spaces to store "exterior to boundary" values

415:   */
416:   /* FLUID */
417:   if (info1.gxs == -3) {                   /* 3 points at left end of line */
418:     for (i=-3; i<0; i++) {
419:       x1[i].prss = app->prin;
420:       x1[i].ergg = app->ugin;
421:       x1[i].ergf = app->ufin;
422:       x1[i].alfg = app->agin;
423:       x1[i].velg = app->vgin;
424:       x1[i].velf = app->vfin;
425:     }
426:   }
427:   if (info1.gxs+info1.gxm == info1.mx+3) { /* 3 points at right end of line */
428:     for (i=info1.mx; i<info1.mx+3; i++) {
429:       x1[i].prss = app->prout;
430:       x1[i].ergg = app->ugout;
431:       x1[i].ergf = app->ufout;
432:       x1[i].alfg = app->agout;
433:       x1[i].velg = app->vgi;
434:       x1[i].velf = app->vfi;
435:     }
436:   }

438:   /* Thermal */
439:   if (info2.gxs == -1) {                                      /* left side of domain */
440:     for (j=info2.gys;j<info2.gys+info2.gym; j++) {
441:       x2[j][-1] = app->twi;
442:     }
443:   }
444:   if (info2.gxs+info2.gxm == info2.mx+1) {                   /* right side of domain */
445:     for (j=info2.gys;j<info2.gys+info2.gym; j++) {
446:       x2[j][info2.mx] = app->twi;
447:     }
448:   }

450:   /* FUEL */
451:   if (info3.gxs == -1) {                                      /* left side of domain */
452:     for (j=info3.gys;j<info3.gys+info3.gym; j++) {
453:       x3[j][-1].phii = app->phii;
454:       x3[j][-1].prei = app->prei;
455:     }
456:   }
457:   if (info3.gxs+info3.gxm == info3.mx+1) {                   /* right side of domain */
458:     for (j=info3.gys;j<info3.gys+info3.gym; j++) {
459:       x3[j][info3.mx].phii = app->phii;
460:       x3[j][info3.mx].prei = app->prei;
461:     }
462:   }
463:   if (info3.gys == -1) {                                      /* bottom of domain */
464:     for (i=info3.gxs;i<info3.gxs+info3.gxm; i++) {
465:       x3[-1][i].phii = app->phii;
466:       x3[-1][i].prei = app->prei;
467:     }
468:   }
469:   if (info3.gys+info3.gym == info3.my+1) {                   /* top of domain */
470:     for (i=info3.gxs;i<info3.gxs+info3.gxm; i++) {
471:       x3[info3.my][i].phii = app->phii;
472:       x3[info3.my][i].prei = app->prei;
473:     }
474:   }

476:   /* Access the three subvectors in F: these are not ghosted and directly access the memory locations in F */
477:   DMCompositeGetAccess(dm,F,&F1,&F2,&F3);

479:   /* Access the arrays inside the subvectors of F */
480:   DMDAVecGetArray(da1,F1,(void**)&f1);
481:   DMDAVecGetArray(da2,F2,(void**)&f2);
482:   DMDAVecGetArray(da3,F3,(void**)&f3);

484:   /* Evaluate local user provided function */
485:   FormFunctionLocalFluid(&info1,x1,f1);
486:   FormFunctionLocalThermal(&info2,x2,f2);
487:   FormFunctionLocalFuel(&info3,x3,f3);

489:   DMDAVecRestoreArray(da1,X1,(void**)&x1);
490:   DMDAVecRestoreArray(da2,X2,(void**)&x2);
491:   DMDAVecRestoreArray(da3,X3,(void**)&x3);
492:   DMCompositeRestoreLocalVectors(dm,&X1,&X2,&X3);

494:   DMDAVecRestoreArray(da1,F1,(void**)&f1);
495:   DMDAVecRestoreArray(da2,F2,(void**)&f2);
496:   DMDAVecRestoreArray(da3,F3,(void**)&f3);
497:   DMCompositeRestoreAccess(dm,F,&F1,&F2,&F3);
498:   return(0);
499: }

503: PetscErrorCode MyFormMatrix(DMMG fdmmg,Mat A,Mat B)
504: {

508:   MatShift(A,1.0);
509:   return(0);
510: }

514: /* 
515:    Setup for the custom preconditioner

517:  */
518: PetscErrorCode MyPCSetUp(PC pc)
519: {
520:   AppCtx         *app;
522:   DM             da;

525:   PCShellGetContext(pc,(void**)&app);
526:   /* create the linear solver for the Neutron diffusion */
527:   DMMGCreate(app->comm,1,0,&app->fdmmg);
528:   DMMGSetOptionsPrefix(app->fdmmg,"phi_");
529:   DMMGSetUser(app->fdmmg,0,app);
530:   DMDACreate2d(app->comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,app->nxv,app->nyvf,PETSC_DETERMINE,1,1,1,0,0,&da);
531:   DMMGSetDM(app->fdmmg,(DM)da);
532:   DMMGSetKSP(app->fdmmg,PETSC_NULL,MyFormMatrix);
533:   app->dx = DMMGGetRHS(app->fdmmg);
534:   app->dy = DMMGGetx(app->fdmmg);
535:   VecDuplicate(app->dy,&app->c);
536:   DMDestroy(&da);
537:   return(0);
538: }

542: /* 
543:    Here is my custom preconditioner

545:     Capital vectors: X, X1 are global vectors
546:     Small vectors: x, x1 are local ghosted vectors
547:     Prefixed a: ax1, aY1 are arrays that access the vector values (either local (ax1) or global aY1)

549:  */
550: PetscErrorCode MyPCApply(PC pc,Vec X,Vec Y)
551: {
552:   AppCtx         *app;
554:   Vec            X1,X2,X3,x1,x2,Y1,Y2,Y3;
555:   DMDALocalInfo    info1,info2,info3;
556:   DM             da1,da2,da3;
557:   PetscInt       i,j;
558:   FluidField     *ax1,*aY1;
559:   PetscScalar    **ax2,**aY2;

562:   PCShellGetContext(pc,(void**)&app);
563:   /* obtain information about the three meshes */
564:   DMCompositeGetEntries(app->pack,&da1,&da2,&da3);
565:   DMDAGetLocalInfo(da1,&info1);
566:   DMDAGetLocalInfo(da2,&info2);
567:   DMDAGetLocalInfo(da3,&info3);

569:   /* get ghosted version of fluid and thermal conduction, global for phi and C */
570:   DMCompositeGetAccess(app->pack,X,&X1,&X2,&X3);
571:   DMCompositeGetLocalVectors(app->pack,&x1,&x2,PETSC_NULL);
572:   DMGlobalToLocalBegin(da1,X1,INSERT_VALUES,x1);
573:   DMGlobalToLocalEnd(da1,X1,INSERT_VALUES,x1);
574:   DMGlobalToLocalBegin(da2,X2,INSERT_VALUES,x2);
575:   DMGlobalToLocalEnd(da2,X2,INSERT_VALUES,x2);

577:   /* get global version of result vector */
578:   DMCompositeGetAccess(app->pack,Y,&Y1,&Y2,&Y3);

580:   /* pull out the phi and C values */
581:   VecStrideGather(X3,0,app->dx,INSERT_VALUES);
582:   VecStrideGather(X3,1,app->c,INSERT_VALUES);

584:   /* update C via formula 38; put back into return vector */
585:   VecAXPY(app->c,0.0,app->dx);
586:   VecScale(app->c,1.0);
587:   VecStrideScatter(app->c,1,Y3,INSERT_VALUES);

589:   /* form the right hand side of the phi equation; solve system; put back into return vector */
590:   VecAXPBY(app->dx,0.0,1.0,app->c);
591:   DMMGSolve(app->fdmmg);
592:   VecStrideScatter(app->dy,0,Y3,INSERT_VALUES);

594:   /* access the ghosted x1 and x2 as arrays */
595:   DMDAVecGetArray(da1,x1,&ax1);
596:   DMDAVecGetArray(da2,x2,&ax2);

598:   /* access global y1 and y2 as arrays */
599:   DMDAVecGetArray(da1,Y1,&aY1);
600:   DMDAVecGetArray(da2,Y2,&aY2);

602:   for (i=info1.xs; i<info1.xs+info1.xm; i++) {
603:     aY1[i].prss = ax1[i].prss;
604:     aY1[i].ergg = ax1[i].ergg;
605:     aY1[i].ergf = ax1[i].ergf;
606:     aY1[i].alfg = ax1[i].alfg;
607:     aY1[i].velg = ax1[i].velg;
608:     aY1[i].velf = ax1[i].velf;
609:   }

611:   for (j=info2.ys; j<info2.ys+info2.ym; j++) {
612:     for (i=info2.xs; i<info2.xs+info2.xm; i++) {
613:       aY2[j][i] = ax2[j][i];
614:     }
615:   }

617:   DMDAVecRestoreArray(da1,x1,&ax1);
618:   DMDAVecRestoreArray(da2,x2,&ax2);
619:   DMDAVecRestoreArray(da1,Y1,&aY1);
620:   DMDAVecRestoreArray(da2,Y2,&aY2);

622:   DMCompositeRestoreLocalVectors(app->pack,&x1,&x2,PETSC_NULL);
623:   DMCompositeRestoreAccess(app->pack,X,&X1,&X2,&X3);
624:   DMCompositeRestoreAccess(app->pack,Y,&Y1,&Y2,&Y3);

626:   return(0);
627: }

631: PetscErrorCode MyPCDestroy(PC pc)
632: {
633:   AppCtx         *app;

637:   PCShellGetContext(pc,(void**)&app);
638:   DMMGDestroy(app->fdmmg);
639:   VecDestroy(&app->c);
640:   return(0);
641: }