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:     DMCompositeAddDA(app.pack,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:     DMCompositeAddDA(app.pack,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:     DMCompositeAddDA(app.pack,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);

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

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

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


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

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

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

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

219: /* ------------------------------------------------------------------- */


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

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

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

242:   return(0);
243: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

501: PetscErrorCode MyFormMatrix(DMMG fdmmg,Mat A,Mat B)
502: {
503:   AppCtx         *app = (AppCtx*)fdmmg->user;

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

513: /* 
514:    Setup for the custom preconditioner

516:  */
517: PetscErrorCode MyPCSetUp(void* ctx)
518: {
519:   AppCtx         *app = (AppCtx*)ctx;
521:   DA             da;

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

540: /* 
541:    Here is my custom preconditioner

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

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

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

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

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

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

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

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

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

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

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

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

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

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

623:   return(0);
624: }

628: PetscErrorCode MyPCDestroy(void* ctx)
629: {
630:   AppCtx         *app = (AppCtx*)ctx;

634:   DMMGDestroy(app->fdmmg);
635:   VecDestroy(app->c);
636:   return(0);
637: }