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

 11:                     nyv  -->       --------------------- DA2
 12:                                    |                    | 
 13:                                    |                    | 
 14:                                    |                    |   
 15:                                    |                    | 
 16:                     nyvf-1 -->     |                    |         --------------------- DA3
 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 DAs we simply define the DAs 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

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

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

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

 58:   PetscInt    nxv,nyv,nyvf;

 60:   MPI_Comm    comm;

 62:   DMComposite pack;

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

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

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


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

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

104:   PreLoadBegin(PETSC_TRUE,"SetUp");

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

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

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

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

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

136:     DACreate2d(app.comm,DA_YPERIODIC,DA_STENCIL_STAR,app.nxv,app.nyv,PETSC_DETERMINE,1,1,1,0,0,&da);
137:     DASetFieldName(da,0,"Tempature");
138:     DMCompositeAddDM(app.pack,(DM)da);
139:     DADestroy(da);

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

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

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

166:     app.twi = 373.15e+0;

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

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


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

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

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

200:   PreLoadStage("Solve");
201:     DMMGSolve(dmmg);


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

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

211:     PetscViewerDestroy(v1);
212:     DMCompositeDestroy(app.pack);
213:     DMMGDestroy(dmmg);
214:   PreLoadEnd();

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

220: /* ------------------------------------------------------------------- */


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

226:  */
229: PetscErrorCode FormInitialGuessLocalFluid(AppCtx *app,DALocalInfo *info,FluidField *f)
230: {
231:   PetscInt       i;

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

243:   return(0);
244: }

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

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

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

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

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

280: */
283: PetscErrorCode FormFunctionLocalFluid(DALocalInfo *info,FluidField *u,FluidField *f)
284: {
285:   PetscInt       i;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

513: /* 
514:    Setup for the custom preconditioner

516:  */
517: PetscErrorCode MyPCSetUp(PC pc)
518: {
519:   AppCtx         *app;
521:   DA             da;

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

541: /* 
542:    Here is my custom preconditioner

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

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

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

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

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

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

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

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

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

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

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

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

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

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

625:   return(0);
626: }

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

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