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