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