Actual source code: ex11.c petsc-3.4.5 2014-06-29
1: static char help[] = "Second Order TVD Finite Volume Example.\n" ;
37: #include <petscts.h>
38: #include <petscdmplex.h>
39: #include <petscsf.h>
40: #include <petscblaslapack.h>
41: #if defined(PETSC_HAVE_EXODUSII)
42: #include <exodusII.h>
43: #else
44: #error This example requires ExodusII support. Reconfigure using --download-exodusii
45: #endif
47: #define DIM 2 /* Geometric dimension */ 48: #define ALEN(a) (sizeof(a)/sizeof((a)[0])) 50: static PetscFunctionList PhysicsList;
51: static PetscFunctionList LimitList;
53: /* Represents continuum physical equations. */
54: typedef struct _n_Physics *Physics;
56: /* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
57: * discretization-independent, but its members depend on the scenario being solved. */
58: typedef struct _n_Model *Model;
60: /* 'User' implements a discretization of a continuous model. */
61: typedef struct _n_User *User;
63: typedef PetscErrorCode (*RiemannFunction)(Physics,const PetscReal *,const PetscReal *,const PetscScalar *,const PetscScalar *,PetscScalar *) 64: typedef PetscErrorCode (*SolutionFunction)(Model,PetscReal ,const PetscReal *,PetscScalar *,void*) 65: typedef PetscErrorCode (*BoundaryFunction)(Model,PetscReal ,const PetscReal *,const PetscReal *,const PetscScalar *,PetscScalar *,void*) 66: typedef PetscErrorCode (*FunctionalFunction)(Model,PetscReal ,const PetscReal *,const PetscScalar *,PetscReal *,void*) 67: typedef PetscErrorCode (*SetupFields)(Physics,PetscSection ) 68: static PetscErrorCode ModelBoundaryRegister(Model,const char*,BoundaryFunction,void*,PetscInt ,const PetscInt *) 69: static PetscErrorCode ModelSolutionSetDefault(Model,SolutionFunction,void*) 70: static PetscErrorCode ModelFunctionalRegister(Model,const char*,PetscInt *,FunctionalFunction,void*) 71: static PetscErrorCode OutputVTK(DM,const char*,PetscViewer *) 73: 74: const char *name;
75: PetscInt dof;
76: };
78: typedef struct _n_BoundaryLink *BoundaryLink;
79: 80: char *name;
81: BoundaryFunction func;
82: void *ctx;
83: PetscInt numids;
84: PetscInt *ids;
85: BoundaryLink next;
86: };
88: typedef struct _n_FunctionalLink *FunctionalLink;
89: 90: char *name;
91: FunctionalFunction func;
92: void *ctx;
93: PetscInt offset;
94: FunctionalLink next;
95: };
97: 98: RiemannFunction riemann;
99: PetscInt dof; /* number of degrees of freedom per cell */
100: PetscReal maxspeed; /* kludge to pick initial time step, need to add monitoring and step control */
101: void *data;
102: PetscInt nfields;
103: const struct FieldDescription *field_desc;
104: };
106: 107: MPI_Comm comm; /* Does not do collective communicaton, but some error conditions can be collective */
108: Physics physics;
109: BoundaryLink boundary;
110: FunctionalLink functionalRegistry;
111: PetscInt maxComputed;
112: PetscInt numMonitored;
113: FunctionalLink *functionalMonitored;
114: PetscInt numCall;
115: FunctionalLink *functionalCall;
116: SolutionFunction solution;
117: void *solutionctx;
118: PetscReal maxspeed; /* estimate of global maximum speed (for CFL calculation) */
119: };
121: 122: PetscErrorCode (*RHSFunctionLocal)(DM,DM,DM,PetscReal ,Vec ,Vec ,User);
123: PetscReal (*Limit)(PetscReal );
124: PetscBool reconstruct;
125: PetscInt numGhostCells, numSplitFaces;
126: PetscInt cEndInterior; /* First boundary ghost cell */
127: Vec cellgeom, facegeom;
128: DM dmGrad;
129: PetscReal minradius;
130: PetscInt vtkInterval; /* For monitor */
131: Model model;
132: struct {
133: PetscScalar *flux;
134: PetscScalar *state0;
135: PetscScalar *state1;
136: } work;
137: };
139: typedef struct {
140: PetscScalar normal[DIM]; /* Area-scaled normals */
141: PetscScalar centroid[DIM]; /* Location of centroid (quadrature point) */
142: PetscScalar grad[2][DIM]; /* Face contribution to gradient in left and right cell */
143: } FaceGeom;
145: typedef struct {
146: PetscScalar centroid[DIM];
147: PetscScalar volume;
148: } CellGeom;
151: PetscScalar DotDIM(const PetscScalar *x,const PetscScalar *y)152: {
153: PetscInt i;
154: PetscScalar prod=0.0;
156: for (i=0; i<DIM; i++) prod += x[i]*y[i];
157: return prod;
158: }
159: PetscReal NormDIM(const PetscScalar *x)return PetscSqrtReal(PetscAbsScalar(DotDIM(x,x))); }
160: PetscScalar a,PetscScalar *x)161: {
162: PetscInt i;
163: for (i=0; i<DIM; i++) x[i] *= a;
164: }
165: PetscScalar a,const PetscScalar *x, PetscScalar *w)166: {
167: PetscInt i;
168: for (i=0; i<DIM; i++) w[i] = x[i]*a;
169: }
170: PetscReal *n,const PetscScalar *x,PetscScalar *xn,PetscScalar *xt)171: { /* Split x into normal and tangential components */
172: PetscInt i;
173: PetscScalar c;
174: c = DotDIM(x,n)/DotDIM(n,n);
175: for (i=0; i<DIM; i++) {
176: xn[i] = c*n[i];
177: xt[i] = x[i]-xn[i];
178: }
179: }
181: PetscScalar Dot2(const PetscScalar *x,const PetscScalar *y)return x[0]*y[0] + x[1]*y[1];}
182: PetscReal Norm2(const PetscScalar *x)return PetscSqrtReal(PetscAbsScalar(Dot2(x,x)));}
183: PetscScalar *x)PetscReal a = 1./Norm2(x); x[0] *= a; x[1] *= a; }
184: PetscScalar a,const PetscScalar *x,const PetscScalar *y,PetscScalar *w)185: PetscScalar a,const PetscScalar *x,PetscScalar *y)187: PetscInt dim, PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];}
188: PetscScalar DotD(PetscInt dim, const PetscScalar *x, const PetscScalar *y)PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;}
189: PetscReal NormD(PetscInt dim, const PetscScalar *x)return PetscSqrtReal(PetscAbsScalar(DotD(dim,x,x)));}
191: PetscReal *n,const PetscScalar *x,PetscScalar *xn,PetscScalar *xt)192: { /* Split x into normal and tangential components */
193: Scale2(Dot2(x,n)/Dot2(n,n),n,xn);
194: Waxpy2(-1,xn,x,xt);
195: }
197: /* Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
198: *
199: * The classical flux-limited formulation is psi(r) where
200: *
201: * r = (u[0] - u[-1]) / (u[1] - u[0])
202: *
203: * The second order TVD region is bounded by
204: *
205: * psi_minmod(r) = min(r,1) and psi_superbee(r) = min(2, 2r, max(1,r))
206: *
207: * where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
208: * phi(r)(r+1)/2 in which the reconstructed interface values are
209: *
210: * u(v) = u[0] + phi(r) (grad u)[0] v
211: *
212: * where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
213: *
214: * phi_minmod(r) = 2 min(1/(1+r),r/(1+r)) phi_superbee(r) = 2 min(2/(1+r), 2r/(1+r), max(1,r)/(1+r))
215: *
216: * For a nicer symmetric formulation, rewrite in terms of
217: *
218: * f = (u[0] - u[-1]) / (u[1] - u[-1])
219: *
220: * where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
221: *
222: * phi(r) = phi(1/r)
223: *
224: * becomes
225: *
226: * w(f) = w(1-f).
227: *
228: * The limiters below implement this final form w(f). The reference methods are
229: *
230: * w_minmod(f) = 2 min(f,(1-f)) w_superbee(r) = 4 min((1-f), f)
231: * */
232: PetscReal Limit_Zero(PetscReal f)return 0; }
233: PetscReal Limit_None(PetscReal f)return 1; }
234: PetscReal Limit_Minmod(PetscReal f)return PetscMax(0,PetscMin(f,(1-f))*2); }
235: PetscReal Limit_VanLeer(PetscReal f)return PetscMax(0,4*f*(1-f)); }
236: PetscReal Limit_VanAlbada(PetscReal f)return PetscMax(0,2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f))); }
237: PetscReal Limit_Sin(PetscReal f)238: {
239: PetscReal fclip = PetscMax(0,PetscMin(f,1));
240: return PetscSinReal(PETSC_PI*fclip);
241: }
242: PetscReal Limit_Superbee(PetscReal f)return 2*Limit_Minmod(f); }
243: PetscReal Limit_MC(PetscReal f)return PetscMin(1,Limit_Superbee(f)); } /* aka Barth-Jespersen */
245: /******************* Advect ********************/
246: typedef enum {ADVECT_SOL_TILTED,ADVECT_SOL_BUMP} AdvectSolType;
247: static const char *const AdvectSolTypes[] = {"TILTED" ,"BUMP" ,"AdvectSolType" ,"ADVECT_SOL_" ,0};
248: typedef enum {ADVECT_SOL_BUMP_CONE,ADVECT_SOL_BUMP_COS} AdvectSolBumpType;
249: static const char *const AdvectSolBumpTypes[] = {"CONE" ,"COS" ,"AdvectSolBumpType" ,"ADVECT_SOL_BUMP_" ,0};
251: typedef struct {
252: PetscReal wind[DIM];
253: } Physics_Advect_Tilted;
254: typedef struct {
255: PetscReal center[DIM];
256: PetscReal radius;
257: AdvectSolBumpType type;
258: } Physics_Advect_Bump;
260: typedef struct {
261: PetscReal inflowState;
262: AdvectSolType soltype;
263: union {
264: Physics_Advect_Tilted tilted;
265: Physics_Advect_Bump bump;
266: } sol;
267: struct {
268: PetscInt Error;
269: } functional;
270: } Physics_Advect;
272: static const struct FieldDescription PhysicsFields_Advect[] = {{"U" ,1},{NULL,0}};
276: PetscErrorCode PhysicsBoundary_Advect_Inflow(Model mod, PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)277: {
278: Physics phys = (Physics)ctx;
279: Physics_Advect *advect = (Physics_Advect*)phys->data;
282: xG[0] = advect->inflowState;
283: return (0);
284: }
288: PetscErrorCode PhysicsBoundary_Advect_Outflow(Model mod, PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)289: {
291: xG[0] = xI[0];
292: return (0);
293: }
297: PetscErrorCode PhysicsRiemann_Advect(Physics phys, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux)298: {
299: Physics_Advect *advect = (Physics_Advect*)phys->data;
300: PetscReal wind[DIM],wn;
303: switch (advect->soltype) {
304: case ADVECT_SOL_TILTED: {
305: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
306: wind[0] = tilted->wind[0];
307: wind[1] = tilted->wind[1];
308: } break ;
309: case ADVECT_SOL_BUMP:
310: wind[0] = -qp[1];
311: wind[1] = qp[0];
312: break ;
313: default: SETERRQ1 (PETSC_COMM_SELF ,PETSC_ERR_SUP,"No support for solution type %s" ,AdvectSolBumpTypes[advect->soltype]);
314: }
315: wn = Dot2(wind, n);
316: flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
317: return (0);
318: }
322: PetscErrorCode PhysicsSolution_Advect(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)323: {
324: Physics phys = (Physics)ctx;
325: Physics_Advect *advect = (Physics_Advect*)phys->data;
328: switch (advect->soltype) {
329: case ADVECT_SOL_TILTED: {
330: PetscReal x0[DIM];
331: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
332: Waxpy2(-time,tilted->wind,x,x0);
333: if (x0[1] > 0) u[0] = 1.*x[0] + 3.*x[1];
334: else u[0] = advect->inflowState;
335: } break ;
336: case ADVECT_SOL_BUMP: {
337: Physics_Advect_Bump *bump = &advect->sol.bump;
338: PetscReal x0[DIM],v[DIM],r,cost,sint;
339: cost = PetscCosReal(time);
340: sint = PetscSinReal(time);
341: x0[0] = cost*x[0] + sint*x[1];
342: x0[1] = -sint*x[0] + cost*x[1];
343: Waxpy2(-1,bump->center,x0,v);
344: r = Norm2(v);
345: switch (bump->type) {
346: case ADVECT_SOL_BUMP_CONE:
347: u[0] = PetscMax(1 - r/bump->radius,0);
348: break ;
349: case ADVECT_SOL_BUMP_COS:
350: u[0] = 0.5 + 0.5*PetscCosReal(PetscMin(r/bump->radius,1)*PETSC_PI);
351: break ;
352: }
353: } break ;
354: default: SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_SUP,"Unknown solution type" );
355: }
356: return (0);
357: }
361: PetscErrorCode PhysicsFunctional_Advect(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *y,PetscReal *f,void *ctx)362: {
363: Physics phys = (Physics)ctx;
364: Physics_Advect *advect = (Physics_Advect*)phys->data;
365: PetscScalar yexact[1];
369: PhysicsSolution_Advect(mod,time,x,yexact,phys);
370: f[advect->functional.Error] = PetscAbsScalar(y[0]-yexact[0]);
371: return (0);
372: }
376: PetscErrorCode PhysicsCreate_Advect(Model mod,Physics phys)377: {
378: Physics_Advect *advect = (Physics_Advect*)phys->data;
382: phys->field_desc = PhysicsFields_Advect;
383: phys->riemann = PhysicsRiemann_Advect;
384: PetscNew (Physics_Advect,&phys->data);
385: advect = phys->data;
386: PetscOptionsHead ("Advect options" );
387: {
388: PetscInt two = 2,dof = 1;
389: advect->soltype = ADVECT_SOL_TILTED;
390: PetscOptionsEnum ("-advect_sol_type" ,"solution type" ,"" ,AdvectSolTypes,(PetscEnum )advect->soltype,(PetscEnum *)&advect->soltype,NULL);
391: switch (advect->soltype) {
392: case ADVECT_SOL_TILTED: {
393: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
394: two = 2;
395: tilted->wind[0] = 0.0;
396: tilted->wind[1] = 1.0;
397: PetscOptionsRealArray ("-advect_tilted_wind" ,"background wind vx,vy" ,"" ,tilted->wind,&two,NULL);
398: advect->inflowState = -2.0;
399: PetscOptionsRealArray ("-advect_tilted_inflow" ,"Inflow state" ,"" ,&advect->inflowState,&dof,NULL);
400: phys->maxspeed = Norm2(tilted->wind);
401: } break ;
402: case ADVECT_SOL_BUMP: {
403: Physics_Advect_Bump *bump = &advect->sol.bump;
404: two = 2;
405: bump->center[0] = 2.;
406: bump->center[1] = 0.;
407: PetscOptionsRealArray ("-advect_bump_center" ,"location of center of bump x,y" ,"" ,bump->center,&two,NULL);
408: bump->radius = 0.9;
409: PetscOptionsReal ("-advect_bump_radius" ,"radius of bump" ,"" ,bump->radius,&bump->radius,NULL);
410: bump->type = ADVECT_SOL_BUMP_CONE;
411: PetscOptionsEnum ("-advect_bump_type" ,"type of bump" ,"" ,AdvectSolBumpTypes,(PetscEnum )bump->type,(PetscEnum *)&bump->type,NULL);
412: phys->maxspeed = 3.; /* radius of mesh, kludge */
413: } break ;
414: }
415: }
416: PetscOptionsTail ();
418: {
419: const PetscInt inflowids[] = {100,200,300},outflowids[] = {101};
420: /* Register "canned" boundary conditions and defaults for where to apply. */
421: ModelBoundaryRegister(mod,"inflow" ,PhysicsBoundary_Advect_Inflow,phys,ALEN(inflowids),inflowids);
422: ModelBoundaryRegister(mod,"outflow" ,PhysicsBoundary_Advect_Outflow,phys,ALEN(outflowids),outflowids);
423: /* Initial/transient solution with default boundary conditions */
424: ModelSolutionSetDefault(mod,PhysicsSolution_Advect,phys);
425: /* Register "canned" functionals */
426: ModelFunctionalRegister(mod,"Error" ,&advect->functional.Error,PhysicsFunctional_Advect,phys);
427: }
428: return (0);
429: }
431: /******************* Shallow Water ********************/
432: typedef struct {
433: PetscReal gravity;
434: PetscReal boundaryHeight;
435: struct {
436: PetscInt Height;
437: PetscInt Speed;
438: PetscInt Energy;
439: } functional;
440: } Physics_SW;
441: typedef struct {
442: PetscScalar vals[0];
443: PetscScalar h;
444: PetscScalar uh[DIM];
445: } SWNode;
447: static const struct FieldDescription PhysicsFields_SW[] = {{"Height" ,1},{"Momentum" ,DIM},{NULL,0}};
451: /*
452: * h_t + div(uh) = 0
453: * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
454: *
455: * */
456: PetscErrorCode SWFlux(Physics phys,const PetscReal *n,const SWNode *x,SWNode *f)457: {
458: Physics_SW *sw = (Physics_SW*)phys->data;
459: PetscScalar uhn,u[DIM];
460: PetscInt i;
463: Scale2(1./x->h,x->uh,u);
464: uhn = Dot2(x->uh,n);
465: f->h = uhn;
466: for (i=0; i<DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
467: return (0);
468: }
472: PetscErrorCode PhysicsBoundary_SW_Wall(Model mod, PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)473: {
475: xG[0] = xI[0];
476: xG[1] = -xI[1];
477: xG[2] = -xI[2];
478: return (0);
479: }
483: PetscErrorCode PhysicsRiemann_SW(Physics phys, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux)484: {
485: Physics_SW *sw = (Physics_SW*)phys->data;
486: PetscReal cL,cR,speed,nn[DIM];
487: const SWNode *uL = (const SWNode*)xL,*uR = (const SWNode*)xR;
488: SWNode fL,fR;
489: PetscInt i;
492: if (uL->h < 0 || uR->h < 0) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative" );
493: nn[0] = n[0];
494: nn[1] = n[1];
495: Normalize2(nn);
496: SWFlux(phys,nn,uL,&fL);
497: SWFlux(phys,nn,uR,&fR);
498: cL = PetscSqrtReal(sw->gravity*PetscRealPart(uL->h));
499: cR = PetscSqrtReal(sw->gravity*PetscRealPart(uR->h)); /* gravity wave speed */
500: speed = PetscMax(PetscAbsScalar(Dot2(uL->uh,nn)/uL->h) + cL,PetscAbsScalar(Dot2(uR->uh,nn)/uR->h) + cR);
501: for (i=0; i<1+DIM; i++) flux[i] = (0.5*(fL.vals[i] + fR.vals[i]) + 0.5*speed*(xL[i] - xR[i])) * Norm2(n);
502: return (0);
503: }
507: PetscErrorCode PhysicsSolution_SW(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)508: {
509: PetscReal dx[2],r,sigma;
512: if (time != 0.0) SETERRQ1 (mod->comm,PETSC_ERR_SUP,"No solution known for time %G" ,time);
513: dx[0] = x[0] - 1.5;
514: dx[1] = x[1] - 1.0;
515: r = Norm2(dx);
516: sigma = 0.5;
517: u[0] = 1 + 2*PetscExpScalar(-PetscSqr(r)/(2*PetscSqr(sigma)));
518: u[1] = 0.0;
519: u[2] = 0.0;
520: return (0);
521: }
525: PetscErrorCode PhysicsFunctional_SW(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)526: {
527: Physics phys = (Physics)ctx;
528: Physics_SW *sw = (Physics_SW*)phys->data;
529: const SWNode *x = (const SWNode*)xx;
530: PetscScalar u[2];
531: PetscReal h;
534: h = PetscRealPart(x->h);
535: Scale2(1./x->h,x->uh,u);
536: f[sw->functional.Height] = h;
537: f[sw->functional.Speed] = Norm2(u) + PetscSqrtReal(sw->gravity*h);
538: f[sw->functional.Energy] = 0.5*(Dot2(x->uh,u) + sw->gravity*PetscSqr(h));
539: return (0);
540: }
544: PetscErrorCode PhysicsCreate_SW(Model mod,Physics phys)545: {
546: Physics_SW *sw;
550: phys->field_desc = PhysicsFields_SW;
551: phys->riemann = PhysicsRiemann_SW;
552: PetscNew (Physics_SW,&phys->data);
553: sw = phys->data;
554: PetscOptionsHead ("SW options" );
555: {
556: sw->gravity = 1.0;
557: PetscOptionsReal ("-sw_gravity" ,"Gravitational constant" ,"" ,sw->gravity,&sw->gravity,NULL);
558: }
559: PetscOptionsTail ();
560: phys->maxspeed = PetscSqrtReal(2.0*sw->gravity); /* Mach 1 for depth of 2 */
562: {
563: const PetscInt wallids[] = {100,101,200,300};
564: ModelBoundaryRegister(mod,"wall" ,PhysicsBoundary_SW_Wall,phys,ALEN(wallids),wallids);
565: ModelSolutionSetDefault(mod,PhysicsSolution_SW,phys);
566: ModelFunctionalRegister(mod,"Height" ,&sw->functional.Height,PhysicsFunctional_SW,phys);
567: ModelFunctionalRegister(mod,"Speed" ,&sw->functional.Speed,PhysicsFunctional_SW,phys);
568: ModelFunctionalRegister(mod,"Energy" ,&sw->functional.Energy,PhysicsFunctional_SW,phys);
569: }
570: return (0);
571: }
573: /******************* Euler ********************/
574: typedef struct {
575: PetscScalar vals[0];
576: PetscScalar r;
577: PetscScalar ru[DIM];
578: PetscScalar e;
579: } EulerNode;
580: typedef PetscErrorCode (*EquationOfState)(const PetscReal *, const EulerNode*, PetscScalar *) 581: typedef struct {
582: PetscInt npars;
583: PetscReal pars[DIM];
584: EquationOfState pressure;
585: EquationOfState sound;
586: struct {
587: PetscInt Density;
588: PetscInt Momentum;
589: PetscInt Energy;
590: PetscInt Pressure;
591: PetscInt Speed;
592: } monitor;
593: } PhysicsEuler;
595: static const struct FieldDescription PhysicsFields_Euler[] = {{"Density" ,1},{"Momentum" ,DIM},{"Energy" ,1},{NULL,0}};
599: PetscErrorCode Pressure_PG(const PetscReal *pars,const EulerNode *x,PetscScalar *p)600: {
601: PetscScalar ru2;
604: ru2 = DotDIM(x->ru,x->ru);
605: ru2 /= x->r;
606: /* kinematic dof = params[0] */
607: (*p)=2.0*(x->e-0.5*ru2)/pars[0];
608: return (0);
609: }
613: PetscErrorCode SpeedOfSound_PG(const PetscReal *pars,const EulerNode *x,PetscScalar *c)614: {
615: PetscScalar p;
618: /* TODO remove direct usage of Pressure_PG */
619: Pressure_PG(pars,x,&p);
620: /* TODO check the sign of p */
621: /* pars[1] = heat capacity ratio */
622: (*c)=PetscSqrtScalar(pars[1]*p/x->r);
623: return (0);
624: }
628: /*
629: * x = (rho,rho*(u_1),...,rho*e)^T
630: * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
631: *
632: * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
633: *
634: * */
635: PetscErrorCode EulerFlux(Physics phys,const PetscReal *n,const EulerNode *x,EulerNode *f)636: {
637: PhysicsEuler *eu = (PhysicsEuler*)phys->data;
638: PetscScalar u,nu,p;
639: PetscInt i;
642: u = DotDIM(x->ru,x->ru);
643: u /= (x->r * x->r);
644: nu = DotDIM(x->ru,n);
645: /* TODO check the sign of p */
646: eu->pressure(eu->pars,x,&p);
647: f->r = nu * x->r;
648: for (i=0; i<DIM; i++) f->ru[i] = nu * x->ru[i] + n[i]*p;
649: f->e = nu*(x->e+p);
650: return (0);
651: }
653: /* PetscReal * => EulerNode* conversion */
656: PetscErrorCode PhysicsBoundary_Euler_Wall(Model mod, PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)657: {
658: PetscInt i;
659: PetscScalar xn[DIM],xt[DIM];
662: xG[0] = xI[0];
663: NormalSplitDIM(n,xI+1,xn,xt);
664: for (i=0; i<DIM; i++) xG[i+1] = -xn[i]+xt[i];
665: xG[DIM+1] = xI[DIM+1];
666: return (0);
667: }
669: /* PetscReal * => EulerNode* conversion */
672: PetscErrorCode PhysicsRiemann_Euler_Rusanov(Physics phys, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux)673: {
674: PhysicsEuler *eu = (PhysicsEuler*)phys->data;
675: PetscScalar cL,cR,speed;
676: const EulerNode *uL = (const EulerNode*)xL,*uR = (const EulerNode*)xR;
677: EulerNode fL,fR;
678: PetscInt i;
681: if (uL->r < 0 || uR->r < 0) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed density is negative" );
682: EulerFlux(phys,n,uL,&fL);
683: EulerFlux(phys,n,uR,&fR);
684: eu->sound(eu->pars,uL,&cL);
685: eu->sound(eu->pars,uR,&cR);
686: speed = PetscMax(cL,cR)+PetscMax(PetscAbsScalar(DotDIM(uL->ru,n)/NormDIM(n)),PetscAbsScalar(DotDIM(uR->ru,n)/NormDIM(n)));
687: for (i=0; i<2+DIM; i++) flux[i] = 0.5*(fL.vals[i]+fR.vals[i])+0.5*speed*(xL[i]-xR[i]);
688: return (0);
689: }
693: PetscErrorCode PhysicsSolution_Euler(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)694: {
695: PetscInt i;
698: if (time != 0.0) SETERRQ1 (mod->comm,PETSC_ERR_SUP,"No solution known for time %G" ,time);
699: u[0] = 1.0;
700: u[DIM+1] = 1.0+PetscAbsReal(x[0]);
701: for (i=1; i<DIM+1; i++) u[i] = 0.0;
702: return (0);
703: }
707: PetscErrorCode PhysicsFunctional_Euler(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)708: {
709: Physics phys = (Physics)ctx;
710: PhysicsEuler *eu = (PhysicsEuler*)phys->data;
711: const EulerNode *x = (const EulerNode*)xx;
712: PetscScalar p;
715: f[eu->monitor.Density] = x->r;
716: f[eu->monitor.Momentum] = NormDIM(x->ru);
717: f[eu->monitor.Energy] = x->e;
718: f[eu->monitor.Speed] = NormDIM(x->ru)/x->r;
719: eu->pressure(eu->pars, x, &p);
720: f[eu->monitor.Pressure] = p;
721: return (0);
722: }
726: PetscErrorCode PhysicsCreate_Euler(Model mod,Physics phys)727: {
728: PhysicsEuler *eu;
732: phys->field_desc = PhysicsFields_Euler;
733: phys->riemann = PhysicsRiemann_Euler_Rusanov;
734: PetscNew (PhysicsEuler,&phys->data);
735: eu = phys->data;
736: PetscOptionsHead ("Euler options" );
737: {
738: eu->pars[0] = 3.0;
739: eu->pars[1] = 1.67;
740: PetscOptionsReal ("-eu_f" ,"Degrees of freedom" ,"" ,eu->pars[0],&eu->pars[0],NULL);
741: PetscOptionsReal ("-eu_gamma" ,"Heat capacity ratio" ,"" ,eu->pars[1],&eu->pars[1],NULL);
742: }
743: PetscOptionsTail ();
744: eu->pressure = Pressure_PG;
745: eu->sound = SpeedOfSound_PG;
746: phys->maxspeed = 1.0;
747: {
748: const PetscInt wallids[] = {100,101,200,300};
749: ModelBoundaryRegister(mod,"wall" ,PhysicsBoundary_Euler_Wall,phys,ALEN(wallids),wallids);
750: ModelSolutionSetDefault(mod,PhysicsSolution_Euler,phys);
751: ModelFunctionalRegister(mod,"Speed" ,&eu->monitor.Speed,PhysicsFunctional_Euler,phys);
752: ModelFunctionalRegister(mod,"Energy" ,&eu->monitor.Energy,PhysicsFunctional_Euler,phys);
753: ModelFunctionalRegister(mod,"Density" ,&eu->monitor.Density,PhysicsFunctional_Euler,phys);
754: ModelFunctionalRegister(mod,"Momentum" ,&eu->monitor.Momentum,PhysicsFunctional_Euler,phys);
755: ModelFunctionalRegister(mod,"Pressure" ,&eu->monitor.Pressure,PhysicsFunctional_Euler,phys);
756: }
757: return (0);
758: }
762: PetscErrorCode ConstructCellBoundary(DM dm, User user)763: {
764: const char *name = "Cell Sets" ;
765: const char *bdname = "split faces" ;
766: IS regionIS, innerIS;
767: const PetscInt *regions, *cells;
768: PetscInt numRegions, innerRegion, numCells, c;
770: PetscInt cStart, cEnd, fStart, fEnd;
772: PetscBool hasLabel;
776: DMPlexGetHeightStratum (dm, 0, &cStart, &cEnd);
777: DMPlexGetHeightStratum (dm, 1, &fStart, &fEnd);
779: DMPlexHasLabel (dm, name, &hasLabel);
780: if (!hasLabel) return (0);
781: DMPlexGetLabelSize (dm, name, &numRegions);
782: if (numRegions != 2) return (0);
783: /* Get the inner id */
784: DMPlexGetLabelIdIS (dm, name, ®ionIS);
785: ISGetIndices (regionIS, ®ions);
786: innerRegion = regions[0];
787: ISRestoreIndices (regionIS, ®ions);
788: ISDestroy (®ionIS);
789: /* Find the faces between cells in different regions, could call DMPlexCreateNeighborCSR() */
790: DMPlexGetStratumIS (dm, name, innerRegion, &innerIS);
791: ISGetLocalSize (innerIS, &numCells);
792: ISGetIndices (innerIS, &cells);
793: DMPlexCreateLabel (dm, bdname);
794: for (c = 0; c < numCells; ++c) {
795: const PetscInt cell = cells[c];
796: const PetscInt *faces;
797: PetscInt numFaces, f;
799: if ((cell < cStart) || (cell >= cEnd)) SETERRQ1 (PETSC_COMM_SELF , PETSC_ERR_LIB, "Got invalid point %d which is not a cell" , cell);
800: DMPlexGetConeSize (dm, cell, &numFaces);
801: DMPlexGetCone (dm, cell, &faces);
802: for (f = 0; f < numFaces; ++f) {
803: const PetscInt face = faces[f];
804: const PetscInt *neighbors;
805: PetscInt nC, regionA, regionB;
807: if ((face < fStart) || (face >= fEnd)) SETERRQ1 (PETSC_COMM_SELF , PETSC_ERR_LIB, "Got invalid point %d which is not a face" , face);
808: DMPlexGetSupportSize (dm, face, &nC);
809: if (nC != 2) continue ;
810: DMPlexGetSupport (dm, face, &neighbors);
811: if ((neighbors[0] >= user->cEndInterior) || (neighbors[1] >= user->cEndInterior)) continue ;
812: if ((neighbors[0] < cStart) || (neighbors[0] >= cEnd)) SETERRQ1 (PETSC_COMM_SELF , PETSC_ERR_LIB, "Got invalid point %d which is not a cell" , neighbors[0]);
813: if ((neighbors[1] < cStart) || (neighbors[1] >= cEnd)) SETERRQ1 (PETSC_COMM_SELF , PETSC_ERR_LIB, "Got invalid point %d which is not a cell" , neighbors[1]);
814: DMPlexGetLabelValue (dm, name, neighbors[0], ®ionA);
815: DMPlexGetLabelValue (dm, name, neighbors[1], ®ionB);
816: if (regionA < 0) SETERRQ2 (PetscObjectComm ((PetscObject )dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value" , name, neighbors[0]);
817: if (regionB < 0) SETERRQ2 (PetscObjectComm ((PetscObject )dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value" , name, neighbors[1]);
818: if (regionA != regionB) {
819: DMPlexSetLabelValue (dm, bdname, faces[f], 1);
820: }
821: }
822: }
823: ISRestoreIndices (innerIS, &cells);
824: ISDestroy (&innerIS);
825: {
826: DMLabel label;
828: PetscViewerASCIISynchronizedAllow (PETSC_VIEWER_STDOUT_WORLD , PETSC_TRUE );
829: DMPlexGetLabel (dm, bdname, &label);
830: DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD );
831: }
832: return (0);
833: }
837: /* Right now, I have just added duplicate faces, which see both cells. We can
838: - Add duplicate vertices and decouple the face cones
839: - Disconnect faces from cells across the rotation gap
840: */
841: PetscErrorCode SplitFaces(DM *dmSplit, const char labelName[], User user)842: {
843: DM dm = *dmSplit, sdm;
844: PetscSF sfPoint, gsfPoint;
845: PetscSection coordSection, newCoordSection;
846: Vec coordinates;
847: IS idIS;
848: const PetscInt *ids;
849: PetscInt *newpoints;
850: PetscInt dim, depth, maxConeSize, maxSupportSize, numLabels;
851: PetscInt numFS, fs, pStart, pEnd, p, vStart, vEnd, v, fStart, fEnd, newf, d, l;
852: PetscBool hasLabel;
856: DMPlexHasLabel (dm, labelName, &hasLabel);
857: if (!hasLabel) return (0);
858: DMCreate (PetscObjectComm ((PetscObject )dm), &sdm);
859: DMSetType (sdm, DMPLEX );
860: DMPlexGetDimension (dm, &dim);
861: DMPlexSetDimension (sdm, dim);
863: DMPlexGetLabelIdIS (dm, labelName, &idIS);
864: ISGetLocalSize (idIS, &numFS);
865: ISGetIndices (idIS, &ids);
867: user->numSplitFaces = 0;
868: for (fs = 0; fs < numFS; ++fs) {
869: PetscInt numBdFaces;
871: DMPlexGetStratumSize (dm, labelName, ids[fs], &numBdFaces);
872: user->numSplitFaces += numBdFaces;
873: }
874: DMPlexGetChart (dm, &pStart, &pEnd);
875: pEnd += user->numSplitFaces;
876: DMPlexSetChart (sdm, pStart, pEnd);
877: /* Set cone and support sizes */
878: DMPlexGetDepth (dm, &depth);
879: for (d = 0; d <= depth; ++d) {
880: DMPlexGetDepthStratum (dm, d, &pStart, &pEnd);
881: for (p = pStart; p < pEnd; ++p) {
882: PetscInt newp = p;
883: PetscInt size;
885: DMPlexGetConeSize (dm, p, &size);
886: DMPlexSetConeSize (sdm, newp, size);
887: DMPlexGetSupportSize (dm, p, &size);
888: DMPlexSetSupportSize (sdm, newp, size);
889: }
890: }
891: DMPlexGetHeightStratum (dm, 1, &fStart, &fEnd);
892: for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
893: IS faceIS;
894: const PetscInt *faces;
895: PetscInt numFaces, f;
897: DMPlexGetStratumIS (dm, labelName, ids[fs], &faceIS);
898: ISGetLocalSize (faceIS, &numFaces);
899: ISGetIndices (faceIS, &faces);
900: for (f = 0; f < numFaces; ++f, ++newf) {
901: PetscInt size;
903: /* Right now I think that both faces should see both cells */
904: DMPlexGetConeSize (dm, faces[f], &size);
905: DMPlexSetConeSize (sdm, newf, size);
906: DMPlexGetSupportSize (dm, faces[f], &size);
907: DMPlexSetSupportSize (sdm, newf, size);
908: }
909: ISRestoreIndices (faceIS, &faces);
910: ISDestroy (&faceIS);
911: }
912: DMSetUp (sdm);
913: /* Set cones and supports */
914: DMPlexGetMaxSizes (dm, &maxConeSize, &maxSupportSize);
915: PetscMalloc (PetscMax(maxConeSize, maxSupportSize) * sizeof (PetscInt ), &newpoints);
916: DMPlexGetChart (dm, &pStart, &pEnd);
917: for (p = pStart; p < pEnd; ++p) {
918: const PetscInt *points, *orientations;
919: PetscInt size, i, newp = p;
921: DMPlexGetConeSize (dm, p, &size);
922: DMPlexGetCone (dm, p, &points);
923: DMPlexGetConeOrientation (dm, p, &orientations);
924: for (i = 0; i < size; ++i) newpoints[i] = points[i];
925: DMPlexSetCone (sdm, newp, newpoints);
926: DMPlexSetConeOrientation (sdm, newp, orientations);
927: DMPlexGetSupportSize (dm, p, &size);
928: DMPlexGetSupport (dm, p, &points);
929: for (i = 0; i < size; ++i) newpoints[i] = points[i];
930: DMPlexSetSupport (sdm, newp, newpoints);
931: }
932: PetscFree (newpoints);
933: for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
934: IS faceIS;
935: const PetscInt *faces;
936: PetscInt numFaces, f;
938: DMPlexGetStratumIS (dm, labelName, ids[fs], &faceIS);
939: ISGetLocalSize (faceIS, &numFaces);
940: ISGetIndices (faceIS, &faces);
941: for (f = 0; f < numFaces; ++f, ++newf) {
942: const PetscInt *points;
944: DMPlexGetCone (dm, faces[f], &points);
945: DMPlexSetCone (sdm, newf, points);
946: DMPlexGetSupport (dm, faces[f], &points);
947: DMPlexSetSupport (sdm, newf, points);
948: }
949: ISRestoreIndices (faceIS, &faces);
950: ISDestroy (&faceIS);
951: }
952: ISRestoreIndices (idIS, &ids);
953: ISDestroy (&idIS);
954: DMPlexStratify (sdm);
955: /* Convert coordinates */
956: DMPlexGetDepthStratum (dm, 0, &vStart, &vEnd);
957: DMPlexGetCoordinateSection (dm, &coordSection);
958: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), &newCoordSection);
959: PetscSectionSetNumFields (newCoordSection, 1);
960: PetscSectionSetFieldComponents (newCoordSection, 0, dim);
961: PetscSectionSetChart (newCoordSection, vStart, vEnd);
962: for (v = vStart; v < vEnd; ++v) {
963: PetscSectionSetDof (newCoordSection, v, dim);
964: PetscSectionSetFieldDof (newCoordSection, v, 0, dim);
965: }
966: PetscSectionSetUp (newCoordSection);
967: DMPlexSetCoordinateSection (sdm, newCoordSection);
968: PetscSectionDestroy (&newCoordSection); /* relinquish our reference */
969: DMGetCoordinatesLocal (dm, &coordinates);
970: DMSetCoordinatesLocal (sdm, coordinates);
971: /* Convert labels */
972: DMPlexGetNumLabels (dm, &numLabels);
973: for (l = 0; l < numLabels; ++l) {
974: const char *lname;
975: PetscBool isDepth;
977: DMPlexGetLabelName (dm, l, &lname);
978: PetscStrcmp (lname, "depth" , &isDepth);
979: if (isDepth) continue ;
980: DMPlexCreateLabel (sdm, lname);
981: DMPlexGetLabelIdIS (dm, lname, &idIS);
982: ISGetLocalSize (idIS, &numFS);
983: ISGetIndices (idIS, &ids);
984: for (fs = 0; fs < numFS; ++fs) {
985: IS pointIS;
986: const PetscInt *points;
987: PetscInt numPoints;
989: DMPlexGetStratumIS (dm, lname, ids[fs], &pointIS);
990: ISGetLocalSize (pointIS, &numPoints);
991: ISGetIndices (pointIS, &points);
992: for (p = 0; p < numPoints; ++p) {
993: PetscInt newpoint = points[p];
995: DMPlexSetLabelValue (sdm, lname, newpoint, ids[fs]);
996: }
997: ISRestoreIndices (pointIS, &points);
998: ISDestroy (&pointIS);
999: }
1000: ISRestoreIndices (idIS, &ids);
1001: ISDestroy (&idIS);
1002: }
1003: /* Convert pointSF */
1004: const PetscSFNode *remotePoints;
1005: PetscSFNode *gremotePoints;
1006: const PetscInt *localPoints;
1007: PetscInt *glocalPoints,*newLocation,*newRemoteLocation;
1008: PetscInt numRoots, numLeaves;
1009: PetscMPIInt numProcs;
1011: MPI_Comm_size (PetscObjectComm ((PetscObject )dm), &numProcs);
1012: DMGetPointSF (dm, &sfPoint);
1013: DMGetPointSF (sdm, &gsfPoint);
1014: DMPlexGetChart (dm,&pStart,&pEnd);
1015: PetscSFGetGraph (sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1016: if (numRoots >= 0) {
1017: PetscMalloc2 (numRoots,PetscInt ,&newLocation,pEnd-pStart,PetscInt ,&newRemoteLocation);
1018: for (l=0; l<numRoots; l++) newLocation[l] = l; /* + (l >= cEnd ? user->numGhostCells : 0); */
1019: PetscSFBcastBegin (sfPoint, MPIU_INT, newLocation, newRemoteLocation);
1020: PetscSFBcastEnd (sfPoint, MPIU_INT, newLocation, newRemoteLocation);
1021: PetscMalloc (numLeaves * sizeof (PetscInt ), &glocalPoints);
1022: PetscMalloc (numLeaves * sizeof (PetscSFNode ), &gremotePoints);
1023: for (l = 0; l < numLeaves; ++l) {
1024: glocalPoints[l] = localPoints[l]; /* localPoints[l] >= cEnd ? localPoints[l] + user->numGhostCells : localPoints[l]; */
1025: gremotePoints[l].rank = remotePoints[l].rank;
1026: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
1027: }
1028: PetscFree2 (newLocation,newRemoteLocation);
1029: PetscSFSetGraph (gsfPoint, numRoots+user->numGhostCells, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
1030: }
1031: DMDestroy (dmSplit);
1032: *dmSplit = sdm;
1033: return (0);
1034: }
1038: PetscErrorCode IsExteriorGhostFace(DM dm,PetscInt face,PetscBool *isghost)1039: {
1041: PetscInt ghost,boundary;
1044: *isghost = PETSC_FALSE ;
1045: DMPlexGetLabelValue (dm, "ghost" , face, &ghost);
1046: DMPlexGetLabelValue (dm, "Face Sets" , face, &boundary);
1047: if (ghost >= 0 || boundary >= 0) *isghost = PETSC_TRUE ;
1048: return (0);
1049: }
1053: /* Overwrites A. Can only handle full-rank problems with m>=n */
1054: PetscErrorCode PseudoInverse(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)1055: {
1056: PetscBool debug = PETSC_FALSE ;
1058: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
1059: PetscScalar *R,*Q,*Aback,Alpha;
1062: if (debug) {
1063: PetscMalloc (m*n*sizeof (PetscScalar ),&Aback);
1064: PetscMemcpy (Aback,A,m*n*sizeof (PetscScalar ));
1065: }
1067: PetscBLASIntCast(m,&M);
1068: PetscBLASIntCast(n,&N);
1069: PetscBLASIntCast(mstride,&lda);
1070: PetscBLASIntCast(worksize,&ldwork);
1071: PetscFPTrapPush (PETSC_FP_TRAP_OFF);
1072: LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
1073: PetscFPTrapPop ();
1074: if (info) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_LIB,"xGEQRF error" );
1075: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
1077: /* Extract an explicit representation of Q */
1078: Q = Ainv;
1079: PetscMemcpy (Q,A,mstride*n*sizeof (PetscScalar ));
1080: K = N; /* full rank */
1081: LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
1082: if (info) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_LIB,"xORGQR/xUNGQR error" );
1084: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1085: Alpha = 1.0;
1086: ldb = lda;
1087: BLAStrsm_("Right" ,"Upper" ,"ConjugateTranspose" ,"NotUnitTriangular" ,&M,&N,&Alpha,R,&lda,Q,&ldb);
1088: /* Ainv is Q, overwritten with inverse */
1090: if (debug) { /* Check that pseudo-inverse worked */
1091: PetscScalar Beta = 0.0;
1092: PetscInt ldc;
1093: K = N;
1094: ldc = N;
1095: BLASgemm_("ConjugateTranspose" ,"Normal" ,&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1096: PetscScalarView (n*n,work,PETSC_VIEWER_STDOUT_SELF );
1097: PetscFree (Aback);
1098: }
1099: return (0);
1100: }
1104: PetscErrorCode PseudoInverseGetWorkRequired(PetscInt maxFaces,PetscInt *work)1105: {
1106: PetscInt m,n,nrhs,minwork;
1109: m = maxFaces;
1110: n = DIM;
1111: nrhs = maxFaces;
1112: minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
1113: *work = 5*minwork; /* We can afford to be extra generous */
1114: return (0);
1115: }
1119: /* Overwrites A. Can handle degenerate problems and m<n. */
1120: PetscErrorCode PseudoInverseSVD(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)1121: {
1122: PetscBool debug = PETSC_FALSE ;
1124: PetscInt i,j,maxmn;
1125: PetscBLASInt M,N,nrhs,lda,ldb,irank,ldwork,info;
1126: PetscScalar rcond,*tmpwork,*Brhs,*Aback;
1129: if (debug) {
1130: PetscMalloc (m*n*sizeof (PetscScalar ),&Aback);
1131: PetscMemcpy (Aback,A,m*n*sizeof (PetscScalar ));
1132: }
1134: /* initialize to identity */
1135: tmpwork = Ainv;
1136: Brhs = work;
1137: maxmn = PetscMax(m,n);
1138: for (j=0; j<maxmn; j++) {
1139: for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1140: }
1142: PetscBLASIntCast(m,&M);
1143: PetscBLASIntCast(n,&N);
1144: nrhs = M;
1145: PetscBLASIntCast(mstride,&lda);
1146: PetscBLASIntCast(maxmn,&ldb);
1147: PetscBLASIntCast(worksize,&ldwork);
1148: rcond = -1;
1149: PetscFPTrapPush (PETSC_FP_TRAP_OFF);
1150: LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb,tau,&rcond,&irank,tmpwork,&ldwork,&info);
1151: PetscFPTrapPop ();
1152: if (info) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_LIB,"xGELSS error" );
1153: /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1154: if (irank < PetscMin(M,N)) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_USER,"Rank deficient least squares fit, indicates an isolated cell with two colinear points" );
1156: /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
1157: * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
1158: for (i=0; i<n; i++) {
1159: for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
1160: }
1161: return (0);
1162: }
1166: /* Build least squares gradient reconstruction operators */
1167: PetscErrorCode BuildLeastSquares(DM dm,PetscInt cEndInterior,DM dmFace,PetscScalar *fgeom,DM dmCell,PetscScalar *cgeom)1168: {
1170: PetscInt c,cStart,cEnd,maxNumFaces,worksize;
1171: PetscScalar *B,*Binv,*work,*tau,**gref;
1174: DMPlexGetHeightStratum (dm,0,&cStart,&cEnd);
1175: DMPlexGetMaxSizes (dm,&maxNumFaces,NULL);
1176: PseudoInverseGetWorkRequired(maxNumFaces,&worksize);
1177: PetscMalloc5 (maxNumFaces*DIM,PetscScalar ,&B,worksize,PetscScalar ,&Binv,worksize,PetscScalar ,&work,maxNumFaces,PetscScalar ,&tau,maxNumFaces,PetscScalar *,&gref);
1178: for (c=cStart; c<cEndInterior; c++) {
1179: const PetscInt *faces;
1180: PetscInt numFaces,usedFaces,f,i,j;
1181: const CellGeom *cg;
1182: PetscBool ghost;
1183: DMPlexGetConeSize (dm,c,&numFaces);
1184: if (numFaces < DIM) SETERRQ2 (PETSC_COMM_SELF ,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction" ,c,numFaces);
1185: DMPlexGetCone (dm,c,&faces);
1186: DMPlexPointLocalRead (dmCell,c,cgeom,&cg);
1187: for (f=0,usedFaces=0; f<numFaces; f++) {
1188: const PetscInt *fcells;
1189: PetscInt ncell,side;
1190: FaceGeom *fg;
1191: const CellGeom *cg1;
1192: IsExteriorGhostFace(dm,faces[f],&ghost);
1193: if (ghost) continue ;
1194: DMPlexGetSupport (dm,faces[f],&fcells);
1195: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
1196: ncell = fcells[!side]; /* the neighbor */
1197: DMPlexPointLocalRef (dmFace,faces[f],fgeom,&fg);
1198: DMPlexPointLocalRead (dmCell,ncell,cgeom,&cg1);
1199: for (j=0; j<DIM; j++) B[j*numFaces+usedFaces] = cg1->centroid[j] - cg->centroid[j];
1200: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
1201: }
1202: if (!usedFaces) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_USER,"Mesh contains isolated cell (no neighbors). Is it intentional?" );
1203: /* Overwrites B with garbage, returns Binv in row-major format */
1204: if (0) {
1205: PseudoInverse(usedFaces,numFaces,DIM,B,Binv,tau,worksize,work);
1206: } else {
1207: PseudoInverseSVD(usedFaces,numFaces,DIM,B,Binv,tau,worksize,work);
1208: }
1209: for (f=0,i=0; f<numFaces; f++) {
1210: IsExteriorGhostFace(dm,faces[f],&ghost);
1211: if (ghost) continue ;
1212: for (j=0; j<DIM; j++) gref[i][j] = Binv[j*numFaces+i];
1213: i++;
1214: }
1216: #if 1
1217: if (0) {
1218: PetscReal grad[2] = {0,0};
1219: for (f=0; f<numFaces; f++) {
1220: const PetscInt *fcells;
1221: const CellGeom *cg1;
1222: const FaceGeom *fg;
1223: DMPlexGetSupport (dm,faces[f],&fcells);
1224: DMPlexPointLocalRead (dmFace,faces[f],fgeom,&fg);
1225: for (i=0; i<2; i++) {
1226: if (fcells[i] == c) continue ;
1227: DMPlexPointLocalRead (dmCell,fcells[i],cgeom,&cg1);
1228: PetscScalar du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1229: grad[0] += fg->grad[!i][0] * du;
1230: grad[1] += fg->grad[!i][1] * du;
1231: }
1232: }
1233: printf("cell[%d] grad (%g,%g)\n" ,c,grad[0],grad[1]);
1234: }
1235: #endif
1236: }
1237: PetscFree5 (B,Binv,work,tau,gref);
1238: return (0);
1239: }
1243: /* Set up face data and cell data */
1244: PetscErrorCode ConstructGeometry(DM dm, Vec *facegeom, Vec *cellgeom, User user)1245: {
1246: DM dmFace, dmCell;
1247: PetscSection sectionFace, sectionCell;
1248: PetscSection coordSection;
1249: Vec coordinates;
1250: PetscReal minradius;
1251: PetscScalar *fgeom, *cgeom;
1252: PetscInt dim, cStart, cEnd, c, fStart, fEnd, f;
1256: DMPlexGetDimension (dm, &dim);
1257: if (dim != DIM) SETERRQ2 (PetscObjectComm ((PetscObject )dm),PETSC_ERR_SUP,"No support for dim %D != DIM %D" ,dim,DIM);
1258: DMPlexGetCoordinateSection (dm, &coordSection);
1259: DMGetCoordinatesLocal (dm, &coordinates);
1261: /* Make cell centroids and volumes */
1262: DMPlexClone (dm, &dmCell);
1263: DMPlexSetCoordinateSection (dmCell, coordSection);
1264: DMSetCoordinatesLocal (dmCell, coordinates);
1265: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), §ionCell);
1266: DMPlexGetHeightStratum (dm, 0, &cStart, &cEnd);
1267: PetscSectionSetChart (sectionCell, cStart, cEnd);
1268: for (c = cStart; c < cEnd; ++c) {
1269: PetscSectionSetDof (sectionCell, c, sizeof (CellGeom)/sizeof (PetscScalar ));
1270: }
1271: PetscSectionSetUp (sectionCell);
1272: DMSetDefaultSection (dmCell, sectionCell);
1273: PetscSectionDestroy (§ionCell); /* relinquish our reference */
1275: DMCreateLocalVector (dmCell, cellgeom);
1276: VecGetArray (*cellgeom, &cgeom);
1277: for (c = cStart; c < user->cEndInterior; ++c) {
1278: CellGeom *cg;
1280: DMPlexPointLocalRef (dmCell, c, cgeom, &cg);
1281: PetscMemzero (cg,sizeof (*cg));
1282: DMPlexComputeCellGeometryFVM (dmCell, c, &cg->volume, cg->centroid, NULL);
1283: }
1284: /* Compute face normals and minimum cell radius */
1285: DMPlexClone (dm, &dmFace);
1286: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), §ionFace);
1287: DMPlexGetHeightStratum (dm, 1, &fStart, &fEnd);
1288: PetscSectionSetChart (sectionFace, fStart, fEnd);
1289: for (f = fStart; f < fEnd; ++f) {
1290: PetscSectionSetDof (sectionFace, f, sizeof (FaceGeom)/sizeof (PetscScalar ));
1291: }
1292: PetscSectionSetUp (sectionFace);
1293: DMSetDefaultSection (dmFace, sectionFace);
1294: PetscSectionDestroy (§ionFace);
1295: DMCreateLocalVector (dmFace, facegeom);
1296: VecGetArray (*facegeom, &fgeom);
1297: minradius = PETSC_MAX_REAL;
1298: for (f = fStart; f < fEnd; ++f) {
1299: FaceGeom *fg;
1300: PetscInt ghost;
1302: DMPlexGetLabelValue (dm, "ghost" , f, &ghost);
1303: if (ghost >= 0) continue ;
1304: DMPlexPointLocalRef (dmFace, f, fgeom, &fg);
1305: DMPlexComputeCellGeometryFVM (dm, f, NULL, fg->centroid, fg->normal);
1306: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
1307: {
1308: CellGeom *cL, *cR;
1309: const PetscInt *cells;
1310: PetscReal *lcentroid, *rcentroid;
1311: PetscScalar v[3];
1312: PetscInt d;
1314: DMPlexGetSupport (dm, f, &cells);
1315: DMPlexPointLocalRead (dmCell, cells[0], cgeom, &cL);
1316: DMPlexPointLocalRead (dmCell, cells[1], cgeom, &cR);
1317: lcentroid = cells[0] >= user->cEndInterior ? fg->centroid : cL->centroid;
1318: rcentroid = cells[1] >= user->cEndInterior ? fg->centroid : cR->centroid;
1319: WaxpyD(dim, -1, lcentroid, rcentroid, v);
1320: if (DotD(dim, fg->normal, v) < 0) {
1321: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
1322: }
1323: if (DotD(dim, fg->normal, v) <= 0) {
1324: #if DIM == 2
1325: SETERRQ5(PETSC_COMM_SELF ,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)" , f, fg->normal[0], fg->normal[1], v[0], v[1]);
1326: #elif DIM == 3
1327: SETERRQ7(PETSC_COMM_SELF ,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)" , f, fg->normal[0], fg->normal[1], fg->normal[2], v[0], v[1], v[2]);
1328: #else
1329: # error DIM not supported
1330: #endif
1331: }
1332: if (cells[0] < user->cEndInterior) {
1333: WaxpyD(dim, -1, fg->centroid, cL->centroid, v);
1334: minradius = PetscMin(minradius, NormD(dim, v));
1335: }
1336: if (cells[1] < user->cEndInterior) {
1337: WaxpyD(dim, -1, fg->centroid, cR->centroid, v);
1338: minradius = PetscMin(minradius, NormD(dim, v));
1339: }
1340: }
1341: }
1342: /* Compute centroids of ghost cells */
1343: for (c = user->cEndInterior; c < cEnd; ++c) {
1344: FaceGeom *fg;
1345: const PetscInt *cone, *support;
1346: PetscInt coneSize, supportSize, s;
1348: DMPlexGetConeSize (dmCell, c, &coneSize);
1349: if (coneSize != 1) SETERRQ2 (PETSC_COMM_SELF , PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1" , c, coneSize);
1350: DMPlexGetCone (dmCell, c, &cone);
1351: DMPlexGetSupportSize (dmCell, cone[0], &supportSize);
1352: if (supportSize != 2) SETERRQ2 (PETSC_COMM_SELF , PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1" , cone[0], supportSize);
1353: DMPlexGetSupport (dmCell, cone[0], &support);
1354: DMPlexPointLocalRef (dmFace, cone[0], fgeom, &fg);
1355: for (s = 0; s < 2; ++s) {
1356: /* Reflect ghost centroid across plane of face */
1357: if (support[s] == c) {
1358: const CellGeom *ci;
1359: CellGeom *cg;
1360: PetscScalar c2f[3], a;
1362: DMPlexPointLocalRead (dmCell, support[(s+1)%2], cgeom, &ci);
1363: WaxpyD(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
1364: a = DotD(dim, c2f, fg->normal)/DotD(dim, fg->normal, fg->normal);
1365: DMPlexPointLocalRef (dmCell, support[s], cgeom, &cg);
1366: WaxpyD(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
1367: cg->volume = ci->volume;
1368: }
1369: }
1370: }
1371: if (user->reconstruct) {
1372: PetscSection sectionGrad;
1373: BuildLeastSquares(dm,user->cEndInterior,dmFace,fgeom,dmCell,cgeom);
1374: DMPlexClone (dm,&user->dmGrad);
1375: PetscSectionCreate (PetscObjectComm ((PetscObject )dm),§ionGrad);
1376: PetscSectionSetChart (sectionGrad,cStart,cEnd);
1377: for (c=cStart; c<cEnd; c++) {
1378: PetscSectionSetDof (sectionGrad,c,user->model->physics->dof*DIM);
1379: }
1380: PetscSectionSetUp (sectionGrad);
1381: DMSetDefaultSection (user->dmGrad,sectionGrad);
1382: PetscSectionDestroy (§ionGrad);
1383: }
1384: VecRestoreArray (*facegeom, &fgeom);
1385: VecRestoreArray (*cellgeom, &cgeom);
1386: MPI_Allreduce (&minradius, &user->minradius, 1, MPIU_SCALAR , MPI_MIN, PetscObjectComm ((PetscObject )dm));
1387: DMDestroy (&dmCell);
1388: DMDestroy (&dmFace);
1389: return (0);
1390: }
1394: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)1395: {
1396: PetscSF sfPoint;
1397: PetscSection coordSection;
1398: Vec coordinates;
1399: PetscSection sectionCell;
1400: PetscScalar *part;
1401: PetscInt cStart, cEnd, c;
1402: PetscMPIInt rank;
1406: DMPlexGetCoordinateSection (dm, &coordSection);
1407: DMGetCoordinatesLocal (dm, &coordinates);
1408: DMPlexClone (dm, dmCell);
1409: DMGetPointSF (dm, &sfPoint);
1410: DMSetPointSF (*dmCell, sfPoint);
1411: DMPlexSetCoordinateSection (*dmCell, coordSection);
1412: DMSetCoordinatesLocal (*dmCell, coordinates);
1413: MPI_Comm_rank (PetscObjectComm ((PetscObject )dm), &rank);
1414: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), §ionCell);
1415: DMPlexGetHeightStratum (*dmCell, 0, &cStart, &cEnd);
1416: PetscSectionSetChart (sectionCell, cStart, cEnd);
1417: for (c = cStart; c < cEnd; ++c) {
1418: PetscSectionSetDof (sectionCell, c, 1);
1419: }
1420: PetscSectionSetUp (sectionCell);
1421: DMSetDefaultSection (*dmCell, sectionCell);
1422: PetscSectionDestroy (§ionCell);
1423: DMCreateLocalVector (*dmCell, partition);
1424: PetscObjectSetName ((PetscObject )*partition, "partition" );
1425: VecGetArray (*partition, &part);
1426: for (c = cStart; c < cEnd; ++c) {
1427: PetscScalar *p;
1429: DMPlexPointLocalRef (*dmCell, c, part, &p);
1430: p[0] = rank;
1431: }
1432: VecRestoreArray (*partition, &part);
1433: return (0);
1434: }
1438: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)1439: {
1440: DM dmMass, dmFace, dmCell, dmCoord;
1441: PetscSection coordSection;
1442: Vec coordinates;
1443: PetscSection sectionMass;
1444: PetscScalar *m;
1445: const PetscScalar *facegeom, *cellgeom, *coords;
1446: PetscInt vStart, vEnd, v;
1447: PetscErrorCode ierr;
1450: DMPlexGetCoordinateSection (dm, &coordSection);
1451: DMGetCoordinatesLocal (dm, &coordinates);
1452: DMPlexClone (dm, &dmMass);
1453: DMPlexSetCoordinateSection (dmMass, coordSection);
1454: DMSetCoordinatesLocal (dmMass, coordinates);
1455: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), §ionMass);
1456: DMPlexGetDepthStratum (dm, 0, &vStart, &vEnd);
1457: PetscSectionSetChart (sectionMass, vStart, vEnd);
1458: for (v = vStart; v < vEnd; ++v) {
1459: PetscInt numFaces;
1461: DMPlexGetSupportSize (dmMass, v, &numFaces);
1462: PetscSectionSetDof (sectionMass, v, numFaces*numFaces);
1463: }
1464: PetscSectionSetUp (sectionMass);
1465: DMSetDefaultSection (dmMass, sectionMass);
1466: PetscSectionDestroy (§ionMass);
1467: DMGetLocalVector (dmMass, massMatrix);
1468: VecGetArray (*massMatrix, &m);
1469: VecGetDM (user->facegeom, &dmFace);
1470: VecGetArrayRead (user->facegeom, &facegeom);
1471: VecGetDM (user->cellgeom, &dmCell);
1472: VecGetArrayRead (user->cellgeom, &cellgeom);
1473: DMGetCoordinateDM (dm, &dmCoord);
1474: VecGetArrayRead (coordinates, &coords);
1475: for (v = vStart; v < vEnd; ++v) {
1476: const PetscInt *faces;
1477: const FaceGeom *fgA, *fgB, *cg;
1478: const PetscScalar *vertex;
1479: PetscInt numFaces, sides[2], f, g;
1481: DMPlexPointLocalRead (dmCoord, v, coords, &vertex);
1482: DMPlexGetSupportSize (dmMass, v, &numFaces);
1483: DMPlexGetSupport (dmMass, v, &faces);
1484: for (f = 0; f < numFaces; ++f) {
1485: sides[0] = faces[f];
1486: DMPlexPointLocalRead (dmFace, faces[f], facegeom, &fgA);
1487: for (g = 0; g < numFaces; ++g) {
1488: const PetscInt *cells = NULL;;
1489: PetscReal area = 0.0;
1490: PetscInt numCells;
1492: sides[1] = faces[g];
1493: DMPlexPointLocalRead (dmFace, faces[g], facegeom, &fgB);
1494: DMPlexGetJoin (dmMass, 2, sides, &numCells, &cells);
1495: if (numCells != 1) SETERRQ (PETSC_COMM_SELF , PETSC_ERR_LIB, "Invalid join for faces" );
1496: DMPlexPointLocalRead (dmCell, cells[0], cellgeom, &cg);
1497: area += PetscAbsScalar((vertex[0] - cg->centroid[0])*(fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1])*(fgA->centroid[0] - cg->centroid[0]));
1498: area += PetscAbsScalar((vertex[0] - cg->centroid[0])*(fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1])*(fgB->centroid[0] - cg->centroid[0]));
1499: m[f*numFaces+g] = Dot2(fgA->normal, fgB->normal)*area*0.5;
1500: DMPlexRestoreJoin (dmMass, 2, sides, &numCells, &cells);
1501: }
1502: }
1503: }
1504: VecRestoreArrayRead (user->facegeom, &facegeom);
1505: VecRestoreArrayRead (user->facegeom, &facegeom);
1506: VecRestoreArrayRead (coordinates, &coords);
1507: VecRestoreArray (*massMatrix, &m);
1508: DMDestroy (&dmMass);
1509: return (0);
1510: }
1514: PetscErrorCode SetUpLocalSpace(DM dm, User user)1515: {
1516: PetscSection stateSection;
1517: Physics phys;
1518: PetscInt dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, c, i;
1522: DMPlexGetHeightStratum (dm, 0, &cStart, &cEnd);
1523: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), &stateSection);
1524: phys = user->model->physics;
1525: PetscSectionSetNumFields (stateSection,phys->nfields);
1526: for (i=0; i<phys->nfields; i++) {
1527: PetscSectionSetFieldName (stateSection,i,phys->field_desc[i].name);
1528: PetscSectionSetFieldComponents (stateSection,i,phys->field_desc[i].dof);
1529: }
1530: PetscSectionSetChart (stateSection, cStart, cEnd);
1531: for (c = cStart; c < cEnd; ++c) {
1532: for (i=0; i<phys->nfields; i++) {
1533: PetscSectionSetFieldDof (stateSection,c,i,phys->field_desc[i].dof);
1534: }
1535: PetscSectionSetDof (stateSection, c, dof);
1536: }
1537: for (c = user->cEndInterior; c < cEnd; ++c) {
1538: PetscSectionSetConstraintDof (stateSection, c, dof);
1539: }
1540: PetscSectionSetUp (stateSection);
1541: PetscMalloc (dof * sizeof (PetscInt ), &cind);
1542: for (d = 0; d < dof; ++d) cind[d] = d;
1543: #if 0
1544: for (c = cStart; c < cEnd; ++c) {
1545: PetscInt val;
1547: DMPlexGetLabelValue (dm, "vtk" , c, &val);
1548: if (val < 0) {PetscSectionSetConstraintIndices(stateSection, c, cind);}
1549: }
1550: #endif
1551: for (c = user->cEndInterior; c < cEnd; ++c) {
1552: PetscSectionSetConstraintIndices(stateSection, c, cind);
1553: }
1554: PetscFree (cind);
1555: PetscSectionGetStorageSize (stateSection, &stateSize);
1556: DMSetDefaultSection (dm,stateSection);
1557: PetscSectionDestroy (&stateSection);
1558: return (0);
1559: }
1563: PetscErrorCode SetUpBoundaries(DM dm, User user)1564: {
1565: Model mod = user->model;
1567: BoundaryLink b;
1570: PetscOptionsBegin (PetscObjectComm ((PetscObject )dm),NULL,"Boundary condition options" ,"" );
1571: for (b = mod->boundary; b; b=b->next) {
1572: char optname[512];
1573: PetscInt ids[512],len = 512;
1574: PetscBool flg;
1575: PetscSNPrintf (optname,sizeof optname,"-bc_%s" ,b->name);
1576: PetscMemzero (ids,sizeof (ids));
1577: PetscOptionsIntArray (optname,"List of boundary IDs" ,"" ,ids,&len,&flg);
1578: if (flg) {
1579: /* TODO: check all IDs to make sure they exist in the mesh */
1580: PetscFree (b->ids);
1581: b->numids = len;
1582: PetscMalloc (len*sizeof (PetscInt ),&b->ids);
1583: PetscMemcpy (b->ids,ids,len*sizeof (PetscInt ));
1584: }
1585: }
1586: PetscOptionsEnd ();
1587: return (0);
1588: }
1592: /* The ids are just defaults, can be overridden at command line. I expect to set defaults based on names in the future. */
1593: PetscErrorCode ModelBoundaryRegister(Model mod,const char *name,BoundaryFunction bcFunc,void *ctx,PetscInt numids,const PetscInt *ids)1594: {
1596: BoundaryLink link;
1599: PetscNew (struct _n_BoundaryLink ,&link);
1600: PetscStrallocpy (name,&link->name);
1601: link->numids = numids;
1602: PetscMalloc (numids*sizeof (PetscInt ),&link->ids);
1603: PetscMemcpy (link->ids,ids,numids*sizeof (PetscInt ));
1604: link->func = bcFunc;
1605: link->ctx = ctx;
1606: link->next = mod->boundary;
1607: mod->boundary = link;
1608: return (0);
1609: }
1613: PetscErrorCode BoundaryLinkDestroy(BoundaryLink *link)1614: {
1616: BoundaryLink l,next;
1619: if (!link) return (0);
1620: l = *link;
1621: *link = NULL;
1622: for (; l; l=next) {
1623: next = l->next;
1624: PetscFree (l->ids);
1625: PetscFree (l->name);
1626: PetscFree (l);
1627: }
1628: return (0);
1629: }
1633: PetscErrorCode ModelBoundaryFind(Model mod,PetscInt id,BoundaryFunction *bcFunc,void **ctx)1634: {
1635: BoundaryLink link;
1636: PetscInt i;
1639: *bcFunc = NULL;
1640: for (link=mod->boundary; link; link=link->next) {
1641: for (i=0; i<link->numids; i++) {
1642: if (link->ids[i] == id) {
1643: *bcFunc = link->func;
1644: *ctx = link->ctx;
1645: return (0);
1646: }
1647: }
1648: }
1649: SETERRQ1 (mod->comm,PETSC_ERR_USER,"Boundary ID %D not associated with any registered boundary condition" ,id);
1650: return (0);
1651: }
1654: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
1655: PetscErrorCode ModelSolutionSetDefault(Model mod,SolutionFunction func,void *ctx)1656: {
1658: mod->solution = func;
1659: mod->solutionctx = ctx;
1660: return (0);
1661: }
1665: PetscErrorCode ModelFunctionalRegister(Model mod,const char *name,PetscInt *offset,FunctionalFunction func,void *ctx)1666: {
1668: FunctionalLink link,*ptr;
1669: PetscInt lastoffset = -1;
1672: for (ptr=&mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
1673: PetscNew (struct _n_FunctionalLink ,&link);
1674: PetscStrallocpy (name,&link->name);
1675: link->offset = lastoffset + 1;
1676: link->func = func;
1677: link->ctx = ctx;
1678: link->next = NULL;
1679: *ptr = link;
1680: *offset = link->offset;
1681: return (0);
1682: }
1686: PetscErrorCode ModelFunctionalSetFromOptions(Model mod)1687: {
1689: PetscInt i,j;
1690: FunctionalLink link;
1691: char *names[256];
1694: mod->numMonitored = ALEN(names);
1695: PetscOptionsStringArray ("-monitor" ,"list of functionals to monitor" ,"" ,names,&mod->numMonitored,NULL);
1696: /* Create list of functionals that will be computed somehow */
1697: PetscMalloc (mod->numMonitored*sizeof (FunctionalLink),&mod->functionalMonitored);
1698: /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
1699: PetscMalloc (mod->numMonitored*sizeof (FunctionalLink),&mod->functionalCall);
1700: mod->numCall = 0;
1701: for (i=0; i<mod->numMonitored; i++) {
1702: for (link=mod->functionalRegistry; link; link=link->next) {
1703: PetscBool match;
1704: PetscStrcasecmp (names[i],link->name,&match);
1705: if (match) break ;
1706: }
1707: if (!link) SETERRQ1 (mod->comm,PETSC_ERR_USER,"No known functional '%s'" ,names[i]);
1708: mod->functionalMonitored[i] = link;
1709: for (j=0; j<i; j++) {
1710: if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
1711: }
1712: mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
1713: next_name: 1714: PetscFree (names[i]);
1715: }
1717: /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
1718: mod->maxComputed = -1;
1719: for (link=mod->functionalRegistry; link; link=link->next) {
1720: for (i=0; i<mod->numCall; i++) {
1721: FunctionalLink call = mod->functionalCall[i];
1722: if (link->func == call->func && link->ctx == call->ctx) {
1723: mod->maxComputed = PetscMax(mod->maxComputed,link->offset);
1724: }
1725: }
1726: }
1727: return (0);
1728: }
1732: PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)1733: {
1735: FunctionalLink l,next;
1738: if (!link) return (0);
1739: l = *link;
1740: *link = NULL;
1741: for (; l; l=next) {
1742: next = l->next;
1743: PetscFree (l->name);
1744: PetscFree (l);
1745: }
1746: return (0);
1747: }
1751: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)1752: {
1753: DM dmCell;
1754: Model mod = user->model;
1755: const PetscScalar *cellgeom;
1756: PetscScalar *x;
1757: PetscInt cStart, cEnd, cEndInterior = user->cEndInterior, c;
1758: PetscErrorCode ierr;
1761: VecGetDM (user->cellgeom, &dmCell);
1762: DMPlexGetHeightStratum (dm, 0, &cStart, &cEnd);
1763: VecGetArrayRead (user->cellgeom, &cellgeom);
1764: VecGetArray (X, &x);
1765: for (c = cStart; c < cEndInterior; ++c) {
1766: const CellGeom *cg;
1767: PetscScalar *xc;
1769: DMPlexPointLocalRead (dmCell,c,cellgeom,&cg);
1770: DMPlexPointGlobalRef (dm,c,x,&xc);
1771: if (xc) {(*mod->solution)(mod,0.0,cg->centroid,xc,mod->solutionctx);}
1772: }
1773: VecRestoreArrayRead (user->cellgeom, &cellgeom);
1774: VecRestoreArray (X, &x);
1775: return (0);
1776: }
1780: PetscErrorCode ApplyBC(DM dm, PetscReal time, Vec locX, User user)1781: {
1782: Model mod = user->model;
1783: const char *name = "Face Sets" ;
1784: DM dmFace;
1785: IS idIS;
1786: const PetscInt *ids;
1787: PetscScalar *x;
1788: const PetscScalar *facegeom;
1789: PetscInt numFS, fs;
1790: PetscErrorCode ierr;
1793: VecGetDM (user->facegeom,&dmFace);
1794: DMPlexGetLabelIdIS (dm, name, &idIS);
1795: if (!idIS) return (0);
1796: ISGetLocalSize (idIS, &numFS);
1797: ISGetIndices (idIS, &ids);
1798: VecGetArrayRead (user->facegeom, &facegeom);
1799: VecGetArray (locX, &x);
1800: for (fs = 0; fs < numFS; ++fs) {
1801: BoundaryFunction bcFunc;
1802: void *bcCtx;
1803: IS faceIS;
1804: const PetscInt *faces;
1805: PetscInt numFaces, f;
1807: ModelBoundaryFind(mod,ids[fs],&bcFunc,&bcCtx);
1808: DMPlexGetStratumIS (dm, name, ids[fs], &faceIS);
1809: ISGetLocalSize (faceIS, &numFaces);
1810: ISGetIndices (faceIS, &faces);
1811: for (f = 0; f < numFaces; ++f) {
1812: const PetscInt face = faces[f], *cells;
1813: const PetscScalar *xI;
1814: PetscScalar *xG;
1815: const FaceGeom *fg;
1817: DMPlexPointLocalRead (dmFace, face, facegeom, &fg);
1818: DMPlexGetSupport (dm, face, &cells);
1819: DMPlexPointLocalRead (dm, cells[0], x, &xI);
1820: DMPlexPointLocalRef (dm, cells[1], x, &xG);
1821: (*bcFunc)(mod, time, fg->centroid, fg->normal, xI, xG, bcCtx);
1822: }
1823: ISRestoreIndices (faceIS, &faces);
1824: ISDestroy (&faceIS);
1825: }
1826: VecRestoreArray (locX, &x);
1827: VecRestoreArrayRead (user->facegeom,&facegeom);
1828: ISRestoreIndices (idIS, &ids);
1829: ISDestroy (&idIS);
1830: return (0);
1831: }
1835: PetscErrorCode RHSFunctionLocal_Upwind(DM dm,DM dmFace,DM dmCell,PetscReal time,Vec locX,Vec F,User user)1836: {
1837: Physics phys = user->model->physics;
1838: PetscErrorCode ierr;
1839: const PetscScalar *facegeom, *cellgeom, *x;
1840: PetscScalar *f;
1841: PetscInt fStart, fEnd, face;
1844: VecGetArrayRead (user->facegeom,&facegeom);
1845: VecGetArrayRead (user->cellgeom,&cellgeom);
1846: VecGetArrayRead (locX,&x);
1847: VecGetArray (F,&f);
1848: DMPlexGetHeightStratum (dm, 1, &fStart, &fEnd);
1849: for (face = fStart; face < fEnd; ++face) {
1850: const PetscInt *cells;
1851: PetscInt i,ghost;
1852: PetscScalar *flux = user->work.flux,*fL,*fR;
1853: const FaceGeom *fg;
1854: const CellGeom *cgL,*cgR;
1855: const PetscScalar *xL,*xR;
1857: DMPlexGetLabelValue (dm, "ghost" , face, &ghost);
1858: if (ghost >= 0) continue ;
1859: DMPlexGetSupport (dm, face, &cells);
1860: DMPlexPointLocalRead (dmFace,face,facegeom,&fg);
1861: DMPlexPointLocalRead (dmCell,cells[0],cellgeom,&cgL);
1862: DMPlexPointLocalRead (dmCell,cells[1],cellgeom,&cgR);
1863: DMPlexPointLocalRead (dm,cells[0],x,&xL);
1864: DMPlexPointLocalRead (dm,cells[1],x,&xR);
1865: DMPlexPointGlobalRef (dm,cells[0],f,&fL);
1866: DMPlexPointGlobalRef (dm,cells[1],f,&fR);
1867: (*phys->riemann)(phys, fg->centroid, fg->normal, xL, xR, flux);
1868: for (i=0; i<phys->dof; i++) {
1869: if (fL) fL[i] -= flux[i] / cgL->volume;
1870: if (fR) fR[i] += flux[i] / cgR->volume;
1871: }
1872: }
1873: VecRestoreArrayRead (user->facegeom,&facegeom);
1874: VecRestoreArrayRead (user->cellgeom,&cellgeom);
1875: VecRestoreArrayRead (locX,&x);
1876: VecRestoreArray (F,&f);
1877: return (0);
1878: }
1882: PetscErrorCode RHSFunctionLocal_LS(DM dm,DM dmFace,DM dmCell,PetscReal time,Vec locX,Vec F,User user)1883: {
1884: DM dmGrad = user->dmGrad;
1885: Model mod = user->model;
1886: Physics phys = mod->physics;
1887: const PetscInt dof = phys->dof;
1888: PetscErrorCode ierr;
1889: const PetscScalar *facegeom, *cellgeom, *x;
1890: PetscScalar *f;
1891: PetscInt fStart, fEnd, face, cStart, cell;
1892: Vec locGrad,Grad;
1895: DMGetGlobalVector (dmGrad,&Grad);
1896: VecZeroEntries (Grad);
1897: VecGetArrayRead (user->facegeom,&facegeom);
1898: VecGetArrayRead (user->cellgeom,&cellgeom);
1899: VecGetArrayRead (locX,&x);
1900: DMPlexGetHeightStratum (dm, 1, &fStart, &fEnd);
1901: DMPlexGetHeightStratum (dm, 0, &cStart, NULL);
1902: {
1903: PetscScalar *grad;
1904: VecGetArray (Grad,&grad);
1905: /* Reconstruct gradients */
1906: for (face=fStart; face<fEnd; ++face) {
1907: const PetscInt *cells;
1908: const PetscScalar *cx[2];
1909: const FaceGeom *fg;
1910: PetscScalar *cgrad[2];
1911: PetscInt i,j;
1912: PetscBool ghost;
1914: IsExteriorGhostFace(dm,face,&ghost);
1915: if (ghost) continue ;
1916: DMPlexGetSupport (dm,face,&cells);
1917: DMPlexPointLocalRead (dmFace,face,facegeom,&fg);
1918: for (i=0; i<2; i++) {
1919: DMPlexPointLocalRead (dm,cells[i],x,&cx[i]);
1920: DMPlexPointGlobalRef (dmGrad,cells[i],grad,&cgrad[i]);
1921: }
1922: for (i=0; i<dof; i++) {
1923: PetscScalar delta = cx[1][i] - cx[0][i];
1924: for (j=0; j<DIM; j++) {
1925: if (cgrad[0]) cgrad[0][i*DIM+j] += fg->grad[0][j] * delta;
1926: if (cgrad[1]) cgrad[1][i*DIM+j] -= fg->grad[1][j] * delta;
1927: }
1928: }
1929: }
1930: /* Limit interior gradients. Using cell-based loop because it generalizes better to vector limiters. */
1931: for (cell=cStart; cell<user->cEndInterior; cell++) {
1932: const PetscInt *faces;
1933: PetscInt numFaces,f;
1934: PetscReal *cellPhi = user->work.state0; /* Scalar limiter applied to each component separately */
1935: const PetscScalar *cx;
1936: const CellGeom *cg;
1937: PetscScalar *cgrad;
1938: PetscInt i;
1939: DMPlexGetConeSize (dm,cell,&numFaces);
1940: DMPlexGetCone (dm,cell,&faces);
1941: DMPlexPointLocalRead (dm,cell,x,&cx);
1942: DMPlexPointLocalRead (dmCell,cell,cellgeom,&cg);
1943: DMPlexPointGlobalRef (dmGrad,cell,grad,&cgrad);
1944: if (!cgrad) continue ; /* ghost cell, we don't compute */
1945: /* Limiter will be minimum value over all neighbors */
1946: for (i=0; i<dof; i++) cellPhi[i] = PETSC_MAX_REAL;
1947: for (f=0; f<numFaces; f++) {
1948: const PetscScalar *ncx;
1949: const CellGeom *ncg;
1950: const PetscInt *fcells;
1951: PetscInt face = faces[f],ncell;
1952: PetscScalar v[DIM];
1953: PetscBool ghost;
1954: IsExteriorGhostFace(dm,face,&ghost);
1955: if (ghost) continue ;
1956: DMPlexGetSupport (dm,face,&fcells);
1957: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
1958: DMPlexPointLocalRead (dm,ncell,x,&ncx);
1959: DMPlexPointLocalRead (dmCell,ncell,cellgeom,&ncg);
1960: Waxpy2(-1,cg->centroid,ncg->centroid,v);
1961: for (i=0; i<dof; i++) {
1962: /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
1963: PetscScalar phi,flim = 0.5 * (ncx[i] - cx[i]) / Dot2(&cgrad[i*DIM],v);
1964: phi = (*user->Limit)(flim);
1965: cellPhi[i] = PetscMin(cellPhi[i],phi);
1966: }
1967: }
1968: /* Apply limiter to gradient */
1969: for (i=0; i<dof; i++) Scale2(cellPhi[i],&cgrad[i*DIM],&cgrad[i*DIM]);
1970: }
1971: VecRestoreArray (Grad,&grad);
1972: }
1973: DMGetLocalVector (dmGrad,&locGrad);
1974: DMGlobalToLocalBegin (dmGrad,Grad,INSERT_VALUES ,locGrad);
1975: DMGlobalToLocalEnd (dmGrad,Grad,INSERT_VALUES ,locGrad);
1976: DMRestoreGlobalVector (dmGrad,&Grad);
1978: {
1979: const PetscScalar *grad;
1980: VecGetArrayRead (locGrad,&grad);
1981: VecGetArray (F,&f);
1982: for (face=fStart; face<fEnd; ++face) {
1983: const PetscInt *cells;
1984: PetscInt ghost,i,j,bset;
1985: PetscScalar *flux = user->work.flux,*fx[2] = {user->work.state0,user->work.state1},*cf[2];
1986: const FaceGeom *fg;
1987: const CellGeom *cg[2];
1988: const PetscScalar *cx[2],*cgrad[2];
1990: DMPlexGetLabelValue (dm, "ghost" , face, &ghost);
1991: if (ghost >= 0) continue ;
1992: DMPlexGetSupport (dm, face, &cells);
1993: DMPlexPointLocalRead (dmFace,face,facegeom,&fg);
1994: for (i=0; i<2; i++) {
1995: PetscScalar dx[DIM];
1996: DMPlexPointLocalRead (dmCell,cells[i],cellgeom,&cg[i]);
1997: DMPlexPointLocalRead (dm,cells[i],x,&cx[i]);
1998: DMPlexPointLocalRead (dmGrad,cells[i],grad,&cgrad[i]);
1999: DMPlexPointGlobalRef (dm,cells[i],f,&cf[i]);
2000: Waxpy2(-1,cg[i]->centroid,fg->centroid,dx);
2001: for (j=0; j<dof; j++) fx[i][j] = cx[i][j] + Dot2(cgrad[i],dx);
2002: }
2003: DMPlexGetLabelValue (dm, "Face Sets" , face, &bset);
2004: if (bset != -1) {
2005: BoundaryFunction bcFunc;
2006: void *bcCtx;
2007: ModelBoundaryFind(mod,bset,&bcFunc,&bcCtx);
2008: (*bcFunc)(mod,time,fg->centroid,fg->normal,fx[0],fx[1],bcCtx);
2009: }
2010: (*phys->riemann)(phys, fg->centroid, fg->normal, fx[0], fx[1], flux);
2011: for (i=0; i<phys->dof; i++) {
2012: if (cf[0]) cf[0][i] -= flux[i] / cg[0]->volume;
2013: if (cf[1]) cf[1][i] += flux[i] / cg[1]->volume;
2014: }
2015: }
2016: VecRestoreArrayRead (locGrad,&grad);
2017: VecRestoreArray (F,&f);
2018: }
2019: VecRestoreArrayRead (user->facegeom,&facegeom);
2020: VecRestoreArrayRead (user->cellgeom,&cellgeom);
2021: VecRestoreArrayRead (locX,&x);
2022: DMRestoreLocalVector (dmGrad,&locGrad);
2023: return (0);
2024: }
2028: PetscErrorCode RHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *ctx)2029: {
2030: User user = (User)ctx;
2031: DM dm, dmFace, dmCell;
2032: PetscSection section;
2033: Vec locX;
2037: TSGetDM (ts,&dm);
2038: VecGetDM (user->facegeom,&dmFace);
2039: VecGetDM (user->cellgeom,&dmCell);
2040: DMGetLocalVector (dm,&locX);
2041: DMGlobalToLocalBegin (dm, X, INSERT_VALUES , locX);
2042: DMGlobalToLocalEnd (dm, X, INSERT_VALUES , locX);
2043: DMGetDefaultSection (dm, §ion);
2045: ApplyBC(dm, time, locX, user);
2047: VecZeroEntries (F);
2048: (*user->RHSFunctionLocal)(dm,dmFace,dmCell,time,locX,F,user);
2049: DMRestoreLocalVector (dm,&locX);
2050: return (0);
2051: }
2055: PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)2056: {
2060: PetscViewerCreate (PetscObjectComm ((PetscObject )dm), viewer);
2061: PetscViewerSetType (*viewer, PETSCVIEWERVTK);
2062: PetscViewerFileSetName (*viewer, filename);
2063: return (0);
2064: }
2068: PetscErrorCode MonitorVTK(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx)2069: {
2070: User user = (User)ctx;
2071: DM dm;
2072: PetscViewer viewer;
2073: char filename[PETSC_MAX_PATH_LEN],*ftable = NULL;
2074: PetscReal xnorm;
2078: PetscObjectSetName ((PetscObject ) X, "solution" );
2079: VecGetDM (X,&dm);
2080: VecNorm (X,NORM_INFINITY ,&xnorm);
2081: if (stepnum >= 0) { /* No summary for final time */
2082: Model mod = user->model;
2083: PetscInt c,cStart,cEnd,fcount,i;
2084: size_t ftableused,ftablealloc;
2085: const PetscScalar *cellgeom,*x;
2086: DM dmCell;
2087: PetscReal *fmin,*fmax,*fintegral,*ftmp;
2088: fcount = mod->maxComputed+1;
2089: PetscMalloc4 (fcount,PetscReal ,&fmin,fcount,PetscReal ,&fmax,fcount,PetscReal ,&fintegral,fcount,PetscReal ,&ftmp);
2090: for (i=0; i<fcount; i++) {
2091: fmin[i] = PETSC_MAX_REAL;
2092: fmax[i] = PETSC_MIN_REAL;
2093: fintegral[i] = 0;
2094: }
2095: DMPlexGetHeightStratum (dm,0,&cStart,&cEnd);
2096: VecGetDM (user->cellgeom,&dmCell);
2097: VecGetArrayRead (user->cellgeom,&cellgeom);
2098: VecGetArrayRead (X,&x);
2099: for (c=cStart; c<user->cEndInterior; c++) {
2100: const CellGeom *cg;
2101: const PetscScalar *cx;
2102: DMPlexPointLocalRead (dmCell,c,cellgeom,&cg);
2103: DMPlexPointGlobalRead (dm,c,x,&cx);
2104: if (!cx) continue ; /* not a global cell */
2105: for (i=0; i<mod->numCall; i++) {
2106: FunctionalLink flink = mod->functionalCall[i];
2107: (*flink->func)(mod,time,cg->centroid,cx,ftmp,flink->ctx);
2108: }
2109: for (i=0; i<fcount; i++) {
2110: fmin[i] = PetscMin(fmin[i],ftmp[i]);
2111: fmax[i] = PetscMax(fmax[i],ftmp[i]);
2112: fintegral[i] += cg->volume * ftmp[i];
2113: }
2114: }
2115: VecRestoreArrayRead (user->cellgeom,&cellgeom);
2116: VecRestoreArrayRead (X,&x);
2117: MPI_Allreduce (MPI_IN_PLACE,fmin,fcount,MPIU_REAL,MPI_MIN,PetscObjectComm ((PetscObject )ts));
2118: MPI_Allreduce (MPI_IN_PLACE,fmax,fcount,MPIU_REAL,MPI_MAX,PetscObjectComm ((PetscObject )ts));
2119: MPI_Allreduce (MPI_IN_PLACE,fintegral,fcount,MPIU_REAL,MPI_SUM,PetscObjectComm ((PetscObject )ts));
2121: ftablealloc = fcount * 100;
2122: ftableused = 0;
2123: PetscMalloc (ftablealloc,&ftable);
2124: for (i=0; i<mod->numMonitored; i++) {
2125: size_t countused;
2126: char buffer[256],*p;
2127: FunctionalLink flink = mod->functionalMonitored[i];
2128: PetscInt id = flink->offset;
2129: if (i % 3) {
2130: PetscMemcpy (buffer," " ,2);
2131: p = buffer + 2;
2132: } else if (i) {
2133: char newline[] = "\n" ;
2134: PetscMemcpy (buffer,newline,sizeof newline-1);
2135: p = buffer + sizeof newline - 1;
2136: } else {
2137: p = buffer;
2138: }
2139: PetscSNPrintfCount (p,sizeof buffer-(p-buffer),"%12s [%10.7G,%10.7G] int %10.7G" ,&countused,flink->name,fmin[id],fmax[id],fintegral[id]);
2140: countused += p - buffer;
2141: if (countused > ftablealloc-ftableused-1) { /* reallocate */
2142: char *ftablenew;
2143: ftablealloc = 2*ftablealloc + countused;
2144: PetscMalloc (ftablealloc,&ftablenew);
2145: PetscMemcpy (ftablenew,ftable,ftableused);
2146: PetscFree (ftable);
2147: ftable = ftablenew;
2148: }
2149: PetscMemcpy (ftable+ftableused,buffer,countused);
2150: ftableused += countused;
2151: ftable[ftableused] = 0;
2152: }
2153: PetscFree4 (fmin,fmax,fintegral,ftmp);
2155: PetscPrintf (PetscObjectComm ((PetscObject )ts),"% 3D time %8.4G |x| %8.4G %s\n" ,stepnum,time,xnorm,ftable ? ftable : "" );
2156: PetscFree (ftable);
2157: }
2158: if (user->vtkInterval < 1) return (0);
2159: if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
2160: if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
2161: TSGetTimeStepNumber (ts,&stepnum);
2162: }
2163: PetscSNPrintf (filename,sizeof filename,"ex11-%03D.vtu" ,stepnum);
2164: OutputVTK(dm,filename,&viewer);
2165: VecView (X,viewer);
2166: PetscViewerDestroy (&viewer);
2167: }
2168: return (0);
2169: }
2173: 2174: {
2175: MPI_Comm comm;
2176: User user;
2177: Model mod;
2178: Physics phys;
2179: DM dm, dmDist;
2180: PetscReal ftime,cfl,dt;
2181: PetscInt dim, overlap, nsteps;
2182: int CPU_word_size = 0, IO_word_size = 0, exoid;
2183: float version;
2184: TS ts;
2185: TSConvergedReason reason;
2186: Vec X;
2187: PetscViewer viewer;
2188: PetscMPIInt rank;
2189: char filename[PETSC_MAX_PATH_LEN] = "sevenside.exo" ;
2190: PetscBool vtkCellGeom, splitFaces;
2191: PetscErrorCode ierr;
2193: PetscInitialize (&argc, &argv, (char*) 0, help);
2194: comm = PETSC_COMM_WORLD ;
2195: MPI_Comm_rank (comm, &rank);
2197: PetscNew (struct _n_User ,&user);
2198: PetscNew (struct _n_Model ,&user->model);
2199: PetscNew (struct _n_Physics ,&user->model->physics);
2200: mod = user->model;
2201: phys = mod->physics;
2202: mod->comm = comm;
2204: /* Register physical models to be available on the command line */
2205: PetscFunctionListAdd (&PhysicsList,"advect" ,PhysicsCreate_Advect);
2206: PetscFunctionListAdd (&PhysicsList,"sw" ,PhysicsCreate_SW);
2207: PetscFunctionListAdd (&PhysicsList,"euler" ,PhysicsCreate_Euler);
2209: PetscFunctionListAdd (&LimitList,"zero" ,Limit_Zero);
2210: PetscFunctionListAdd (&LimitList,"none" ,Limit_None);
2211: PetscFunctionListAdd (&LimitList,"minmod" ,Limit_Minmod);
2212: PetscFunctionListAdd (&LimitList,"vanleer" ,Limit_VanLeer);
2213: PetscFunctionListAdd (&LimitList,"vanalbada" ,Limit_VanAlbada);
2214: PetscFunctionListAdd (&LimitList,"sin" ,Limit_Sin);
2215: PetscFunctionListAdd (&LimitList,"superbee" ,Limit_Superbee);
2216: PetscFunctionListAdd (&LimitList,"mc" ,Limit_MC);
2218: PetscOptionsBegin (comm,NULL,"Unstructured Finite Volume Options" ,"" );
2219: {
2220: char physname[256] = "advect" ,limitname[256] = "minmod" ;
2221: PetscErrorCode (*physcreate)(Model,Physics);
2222: cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
2223: PetscOptionsReal ("-ufv_cfl" ,"CFL number per step" ,"" ,cfl,&cfl,NULL);
2224: PetscOptionsString ("-f" ,"Exodus.II filename to read" ,"" ,filename,filename,sizeof (filename),NULL);
2225: user->vtkInterval = 1;
2226: PetscOptionsInt ("-ufv_vtk_interval" ,"VTK output interval (0 to disable)" ,"" ,user->vtkInterval,&user->vtkInterval,NULL);
2227: overlap = 1;
2228: PetscOptionsInt ("-ufv_mesh_overlap" ,"Number of cells to overlap partitions" ,"" ,overlap,&overlap,NULL);
2229: vtkCellGeom = PETSC_FALSE ;
2231: PetscOptionsBool ("-ufv_vtk_cellgeom" ,"Write cell geometry (for debugging)" ,"" ,vtkCellGeom,&vtkCellGeom,NULL);
2232: PetscOptionsList ("-physics" ,"Physics module to solve" ,"" ,PhysicsList,physname,physname,sizeof physname,NULL);
2233: PetscFunctionListFind (PhysicsList,physname,&physcreate);
2234: PetscMemzero (phys,sizeof (struct _n_Physics ));
2235: (*physcreate)(mod,phys);
2236: mod->maxspeed = phys->maxspeed;
2237: /* Count number of fields and dofs */
2238: for (phys->nfields=0,phys->dof=0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
2240: if (mod->maxspeed <= 0) SETERRQ1 (comm,PETSC_ERR_ARG_WRONGSTATE,"Physics '%s' did not set maxspeed" ,physname);
2241: if (phys->dof <= 0) SETERRQ1 (comm,PETSC_ERR_ARG_WRONGSTATE,"Physics '%s' did not set dof" ,physname);
2242: PetscMalloc3 (phys->dof,PetscScalar ,&user->work.flux,phys->dof,PetscScalar ,&user->work.state0,phys->dof,PetscScalar ,&user->work.state1);
2243: user->reconstruct = PETSC_FALSE ;
2244: PetscOptionsBool ("-ufv_reconstruct" ,"Reconstruct gradients for a second order method (grows stencil)" ,"" ,user->reconstruct,&user->reconstruct,NULL);
2245: user->RHSFunctionLocal = user->reconstruct ? RHSFunctionLocal_LS : RHSFunctionLocal_Upwind;
2246: splitFaces = PETSC_FALSE ;
2247: PetscOptionsBool ("-ufv_split_faces" ,"Split faces between cell sets" ,"" ,splitFaces,&splitFaces,NULL);
2248: if (user->reconstruct) {
2249: PetscOptionsList ("-ufv_limit" ,"Limiter to apply to reconstructed solution" ,"" ,LimitList,limitname,limitname,sizeof limitname,NULL);
2250: PetscFunctionListFind (LimitList,limitname,&user->Limit);
2251: }
2252: ModelFunctionalSetFromOptions(mod);
2253: }
2254: PetscOptionsEnd ();
2256: if (!rank) {
2257: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
2258: if (exoid <= 0) SETERRQ1 (PETSC_COMM_SELF ,PETSC_ERR_LIB,"ex_open(\"%s\",...) did not return a valid file ID" ,filename);
2259: } else exoid = -1; /* Not used */
2260: DMPlexCreateExodus (comm, exoid, PETSC_TRUE , &dm);
2261: if (!rank) {ex_close(exoid);}
2262: /* Distribute mesh */
2263: DMPlexDistribute (dm, "chaco" , overlap, &dmDist);
2264: if (dmDist) {
2265: DMDestroy (&dm);
2266: dm = dmDist;
2267: }
2268: DMSetFromOptions (dm);
2269: {
2270: DM gdm;
2272: DMPlexGetHeightStratum (dm, 0, NULL, &user->cEndInterior);
2273: DMPlexConstructGhostCells (dm, NULL, &user->numGhostCells, &gdm);
2274: DMDestroy (&dm);
2275: dm = gdm;
2276: }
2277: if (splitFaces) {ConstructCellBoundary(dm, user);}
2278: SplitFaces(&dm, "split faces" , user);
2279: ConstructGeometry(dm, &user->facegeom, &user->cellgeom, user);
2280: if (0) {VecView (user->cellgeom, PETSC_VIEWER_STDOUT_WORLD );}
2281: DMPlexGetDimension (dm, &dim);
2282: DMPlexSetPreallocationCenterDimension(dm, 0);
2284: /* Set up DM with section describing local vector and configure local vector. */
2285: SetUpLocalSpace(dm, user);
2286: SetUpBoundaries(dm, user);
2288: DMCreateGlobalVector (dm, &X);
2289: PetscObjectSetName ((PetscObject ) X, "solution" );
2290: SetInitialCondition(dm, X, user);
2291: if (vtkCellGeom) {
2292: DM dmCell;
2293: Vec partition;
2295: OutputVTK(dm, "ex11-cellgeom.vtk" , &viewer);
2296: VecView (user->cellgeom, viewer);
2297: PetscViewerDestroy (&viewer);
2298: CreatePartitionVec(dm, &dmCell, &partition);
2299: OutputVTK(dmCell, "ex11-partition.vtk" , &viewer);
2300: VecView (partition, viewer);
2301: PetscViewerDestroy (&viewer);
2302: VecDestroy (&partition);
2303: DMDestroy (&dmCell);
2304: }
2306: TSCreate (comm, &ts);
2307: TSSetType (ts, TSSSP );
2308: TSSetDM (ts, dm);
2309: TSMonitorSet (ts,MonitorVTK,user,NULL);
2310: TSSetRHSFunction (ts,NULL,RHSFunction,user);
2311: TSSetDuration (ts,1000,2.0);
2312: dt = cfl * user->minradius / user->model->maxspeed;
2313: TSSetInitialTimeStep (ts,0.0,dt);
2314: TSSetFromOptions (ts);
2315: TSSolve (ts,X);
2316: TSGetSolveTime (ts,&ftime);
2317: TSGetTimeStepNumber (ts,&nsteps);
2318: TSGetConvergedReason (ts,&reason);
2319: PetscPrintf (PETSC_COMM_WORLD ,"%s at time %G after %D steps\n" ,TSConvergedReasons[reason],ftime,nsteps);
2320: TSDestroy (&ts);
2322: VecDestroy (&user->cellgeom);
2323: VecDestroy (&user->facegeom);
2324: DMDestroy (&user->dmGrad);
2325: PetscFunctionListDestroy (&PhysicsList);
2326: PetscFunctionListDestroy (&LimitList);
2327: BoundaryLinkDestroy(&user->model->boundary);
2328: FunctionalLinkDestroy(&user->model->functionalRegistry);
2329: PetscFree (user->model->functionalMonitored);
2330: PetscFree (user->model->functionalCall);
2331: PetscFree (user->model->physics->data);
2332: PetscFree (user->model->physics);
2333: PetscFree (user->model);
2334: PetscFree3 (user->work.flux,user->work.state0,user->work.state1);
2335: PetscFree (user);
2336: VecDestroy (&X);
2337: DMDestroy (&dm);
2338: PetscFinalize ();
2339: return (0);
2340: }