Actual source code: ex11.c petsc-3.7.5 2017-01-01
1: static char help[] = "Second Order TVD Finite Volume Example.\n" ;
37: #include <petscdmplex.h>
38: #include <petscds.h>
39: #include <petscts.h>
40: #include <petscsf.h> /* For SplitFaces() */ 42: #define DIM 2 /* Geometric dimension */ 43: #define ALEN(a) (sizeof(a)/sizeof((a)[0])) 45: static PetscFunctionList PhysicsList;
47: /* Represents continuum physical equations. */
48: typedef struct _n_Physics *Physics;
50: /* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
51: * discretization-independent, but its members depend on the scenario being solved. */
52: typedef struct _n_Model *Model;
54: /* 'User' implements a discretization of a continuous model. */
55: typedef struct _n_User *User;
56: typedef PetscErrorCode (*SolutionFunction)(Model,PetscReal ,const PetscReal *,PetscScalar *,void*) 57: typedef PetscErrorCode (*SetUpBCFunction)(DM ,Physics) 58: typedef PetscErrorCode (*FunctionalFunction)(Model,PetscReal ,const PetscReal *,const PetscScalar *,PetscReal *,void*) 59: typedef PetscErrorCode (*SetupFields)(Physics,PetscSection ) 60: static PetscErrorCode ModelSolutionSetDefault(Model,SolutionFunction,void*) 61: static PetscErrorCode ModelFunctionalRegister(Model,const char*,PetscInt *,FunctionalFunction,void*) 62: static PetscErrorCode OutputVTK(DM ,const char*,PetscViewer *) 64: 65: const char *name;
66: PetscInt dof;
67: };
69: typedef struct _n_FunctionalLink *FunctionalLink;
70: 71: char *name;
72: FunctionalFunction func;
73: void *ctx;
74: PetscInt offset;
75: FunctionalLink next;
76: };
78: 79: PetscRiemannFunc riemann;
80: PetscInt dof; /* number of degrees of freedom per cell */
81: PetscReal maxspeed; /* kludge to pick initial time step, need to add monitoring and step control */
82: void *data;
83: PetscInt nfields;
84: const struct FieldDescription *field_desc;
85: };
87: 88: MPI_Comm comm; /* Does not do collective communicaton, but some error conditions can be collective */
89: Physics physics;
90: FunctionalLink functionalRegistry;
91: PetscInt maxComputed;
92: PetscInt numMonitored;
93: FunctionalLink *functionalMonitored;
94: PetscInt numCall;
95: FunctionalLink *functionalCall;
96: SolutionFunction solution;
97: SetUpBCFunction setupbc;
98: void *solutionctx;
99: PetscReal maxspeed; /* estimate of global maximum speed (for CFL calculation) */
100: PetscReal bounds[2*DIM];
101: DMBoundaryType bcs[3];
102: };
104: 105: PetscInt numSplitFaces;
106: PetscInt vtkInterval; /* For monitor */
107: Model model;
108: };
110: PetscScalar DotDIM(const PetscScalar *x,const PetscScalar *y)111: {
112: PetscInt i;
113: PetscScalar prod=0.0;
115: for (i=0; i<DIM; i++) prod += x[i]*y[i];
116: return prod;
117: }
118: PetscReal NormDIM(const PetscScalar *x)return PetscSqrtReal(PetscAbsScalar(DotDIM(x,x))); }
120: PetscScalar Dot2(const PetscScalar *x,const PetscScalar *y)return x[0]*y[0] + x[1]*y[1];}
121: PetscReal Norm2(const PetscScalar *x)return PetscSqrtReal(PetscAbsScalar(Dot2(x,x)));}
122: PetscScalar *x)PetscReal a = 1./Norm2(x); x[0] *= a; x[1] *= a; }
123: PetscScalar a,const PetscScalar *x,const PetscScalar *y,PetscScalar *w)124: PetscScalar a,const PetscScalar *x,PetscScalar *y)126: /******************* Advect ********************/
127: typedef enum {ADVECT_SOL_TILTED,ADVECT_SOL_BUMP} AdvectSolType;
128: static const char *const AdvectSolTypes[] = {"TILTED" ,"BUMP" ,"AdvectSolType" ,"ADVECT_SOL_" ,0};
129: typedef enum {ADVECT_SOL_BUMP_CONE,ADVECT_SOL_BUMP_COS} AdvectSolBumpType;
130: static const char *const AdvectSolBumpTypes[] = {"CONE" ,"COS" ,"AdvectSolBumpType" ,"ADVECT_SOL_BUMP_" ,0};
132: typedef struct {
133: PetscReal wind[DIM];
134: } Physics_Advect_Tilted;
135: typedef struct {
136: PetscReal center[DIM];
137: PetscReal radius;
138: AdvectSolBumpType type;
139: } Physics_Advect_Bump;
141: typedef struct {
142: PetscReal inflowState;
143: AdvectSolType soltype;
144: union {
145: Physics_Advect_Tilted tilted;
146: Physics_Advect_Bump bump;
147: } sol;
148: struct {
149: PetscInt Error;
150: } functional;
151: } Physics_Advect;
153: static const struct FieldDescription PhysicsFields_Advect[] = {{"U" ,1},{NULL,0}};
157: PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)158: {
159: Physics phys = (Physics)ctx;
160: Physics_Advect *advect = (Physics_Advect*)phys->data;
163: xG[0] = advect->inflowState;
164: return (0);
165: }
169: PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)170: {
172: xG[0] = xI[0];
173: return (0);
174: }
178: PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)179: {
180: Physics_Advect *advect = (Physics_Advect*)phys->data;
181: PetscReal wind[DIM],wn;
183: switch (advect->soltype) {
184: case ADVECT_SOL_TILTED: {
185: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
186: wind[0] = tilted->wind[0];
187: wind[1] = tilted->wind[1];
188: } break ;
189: case ADVECT_SOL_BUMP:
190: wind[0] = -qp[1];
191: wind[1] = qp[0];
192: break ;
193: /* default: SETERRQ1 (PETSC_COMM_SELF ,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
194: }
195: wn = Dot2(wind, n);
196: flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
197: }
201: PetscErrorCode PhysicsSolution_Advect(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)202: {
203: Physics phys = (Physics)ctx;
204: Physics_Advect *advect = (Physics_Advect*)phys->data;
207: switch (advect->soltype) {
208: case ADVECT_SOL_TILTED: {
209: PetscReal x0[DIM];
210: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
211: Waxpy2(-time,tilted->wind,x,x0);
212: if (x0[1] > 0) u[0] = 1.*x[0] + 3.*x[1];
213: else u[0] = advect->inflowState;
214: } break ;
215: case ADVECT_SOL_BUMP: {
216: Physics_Advect_Bump *bump = &advect->sol.bump;
217: PetscReal x0[DIM],v[DIM],r,cost,sint;
218: cost = PetscCosReal(time);
219: sint = PetscSinReal(time);
220: x0[0] = cost*x[0] + sint*x[1];
221: x0[1] = -sint*x[0] + cost*x[1];
222: Waxpy2(-1,bump->center,x0,v);
223: r = Norm2(v);
224: switch (bump->type) {
225: case ADVECT_SOL_BUMP_CONE:
226: u[0] = PetscMax (1 - r/bump->radius,0);
227: break ;
228: case ADVECT_SOL_BUMP_COS:
229: u[0] = 0.5 + 0.5*PetscCosReal(PetscMin (r/bump->radius,1)*PETSC_PI);
230: break ;
231: }
232: } break ;
233: default: SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_SUP,"Unknown solution type" );
234: }
235: return (0);
236: }
240: PetscErrorCode PhysicsFunctional_Advect(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *y,PetscReal *f,void *ctx)241: {
242: Physics phys = (Physics)ctx;
243: Physics_Advect *advect = (Physics_Advect*)phys->data;
244: PetscScalar yexact[1];
248: PhysicsSolution_Advect(mod,time,x,yexact,phys);
249: f[advect->functional.Error] = PetscAbsScalar(y[0]-yexact[0]);
250: return (0);
251: }
255: PetscErrorCode SetUpBC_Advect(DM dm, Physics phys)256: {
258: const PetscInt inflowids[] = {100,200,300},outflowids[] = {101};
260: /* Register "canned" boundary conditions and defaults for where to apply. */
261: DMAddBoundary (dm, PETSC_TRUE , "inflow" , "Face Sets" , 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Inflow, ALEN(inflowids), inflowids, phys);
262: DMAddBoundary (dm, PETSC_FALSE , "outflow" , "Face Sets" , 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Outflow, ALEN(outflowids), outflowids, phys);
263: return (0);
264: }
268: PetscErrorCode PhysicsCreate_Advect(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)269: {
270: Physics_Advect *advect;
274: phys->field_desc = PhysicsFields_Advect;
275: phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Advect;
276: PetscNew (&advect);
277: phys->data = advect;
278: mod->setupbc = SetUpBC_Advect;
280: PetscOptionsHead (PetscOptionsObject,"Advect options" );
281: {
282: PetscInt two = 2,dof = 1;
283: advect->soltype = ADVECT_SOL_TILTED;
284: PetscOptionsEnum ("-advect_sol_type" ,"solution type" ,"" ,AdvectSolTypes,(PetscEnum )advect->soltype,(PetscEnum *)&advect->soltype,NULL);
285: switch (advect->soltype) {
286: case ADVECT_SOL_TILTED: {
287: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
288: two = 2;
289: tilted->wind[0] = 0.0;
290: tilted->wind[1] = 1.0;
291: PetscOptionsRealArray ("-advect_tilted_wind" ,"background wind vx,vy" ,"" ,tilted->wind,&two,NULL);
292: advect->inflowState = -2.0;
293: PetscOptionsRealArray ("-advect_tilted_inflow" ,"Inflow state" ,"" ,&advect->inflowState,&dof,NULL);
294: phys->maxspeed = Norm2(tilted->wind);
295: } break ;
296: case ADVECT_SOL_BUMP: {
297: Physics_Advect_Bump *bump = &advect->sol.bump;
298: two = 2;
299: bump->center[0] = 2.;
300: bump->center[1] = 0.;
301: PetscOptionsRealArray ("-advect_bump_center" ,"location of center of bump x,y" ,"" ,bump->center,&two,NULL);
302: bump->radius = 0.9;
303: PetscOptionsReal ("-advect_bump_radius" ,"radius of bump" ,"" ,bump->radius,&bump->radius,NULL);
304: bump->type = ADVECT_SOL_BUMP_CONE;
305: PetscOptionsEnum ("-advect_bump_type" ,"type of bump" ,"" ,AdvectSolBumpTypes,(PetscEnum )bump->type,(PetscEnum *)&bump->type,NULL);
306: phys->maxspeed = 3.; /* radius of mesh, kludge */
307: } break ;
308: }
309: }
310: PetscOptionsTail ();
311: /* Initial/transient solution with default boundary conditions */
312: ModelSolutionSetDefault(mod,PhysicsSolution_Advect,phys);
313: /* Register "canned" functionals */
314: ModelFunctionalRegister(mod,"Error" ,&advect->functional.Error,PhysicsFunctional_Advect,phys);
315: mod->bcs[0] = mod->bcs[1] = mod->bcs[2] = DM_BOUNDARY_GHOSTED;
316: return (0);
317: }
319: /******************* Shallow Water ********************/
320: typedef struct {
321: PetscReal gravity;
322: PetscReal boundaryHeight;
323: struct {
324: PetscInt Height;
325: PetscInt Speed;
326: PetscInt Energy;
327: } functional;
328: } Physics_SW;
329: typedef struct {
330: PetscScalar vals[0];
331: PetscScalar h;
332: PetscScalar uh[DIM];
333: } SWNode;
335: static const struct FieldDescription PhysicsFields_SW[] = {{"Height" ,1},{"Momentum" ,DIM},{NULL,0}};
339: /*
340: * h_t + div(uh) = 0
341: * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
342: *
343: * */
344: PetscErrorCode SWFlux(Physics phys,const PetscReal *n,const SWNode *x,SWNode *f)345: {
346: Physics_SW *sw = (Physics_SW*)phys->data;
347: PetscScalar uhn,u[DIM];
348: PetscInt i;
351: Scale2(1./x->h,x->uh,u);
352: uhn = Dot2(x->uh,n);
353: f->h = uhn;
354: for (i=0; i<DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr (x->h) * n[i];
355: return (0);
356: }
360: PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)361: {
363: xG[0] = xI[0];
364: xG[1] = -xI[1];
365: xG[2] = -xI[2];
366: return (0);
367: }
371: PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)372: {
373: Physics_SW *sw = (Physics_SW*)phys->data;
374: PetscReal cL,cR,speed,nn[DIM];
375: const SWNode *uL = (const SWNode*)xL,*uR = (const SWNode*)xR;
376: SWNode fL,fR;
377: PetscInt i;
379: if (uL->h < 0 || uR->h < 0) {for (i=0; i<1+dim; i++) flux[i] = NAN; return ;} /* SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); */
380: nn[0] = n[0];
381: nn[1] = n[1];
382: Normalize2(nn);
383: SWFlux(phys,nn,uL,&fL);
384: SWFlux(phys,nn,uR,&fR);
385: cL = PetscSqrtReal(sw->gravity*PetscRealPart(uL->h));
386: cR = PetscSqrtReal(sw->gravity*PetscRealPart(uR->h)); /* gravity wave speed */
387: speed = PetscMax (PetscAbsScalar(Dot2(uL->uh,nn)/uL->h) + cL,PetscAbsScalar(Dot2(uR->uh,nn)/uR->h) + cR);
388: 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);
389: }
393: PetscErrorCode PhysicsSolution_SW(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)394: {
395: PetscReal dx[2],r,sigma;
398: if (time != 0.0) SETERRQ1 (mod->comm,PETSC_ERR_SUP,"No solution known for time %g" ,(double)time);
399: dx[0] = x[0] - 1.5;
400: dx[1] = x[1] - 1.0;
401: r = Norm2(dx);
402: sigma = 0.5;
403: u[0] = 1 + 2*PetscExpScalar(-PetscSqr (r)/(2*PetscSqr (sigma)));
404: u[1] = 0.0;
405: u[2] = 0.0;
406: return (0);
407: }
411: PetscErrorCode PhysicsFunctional_SW(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)412: {
413: Physics phys = (Physics)ctx;
414: Physics_SW *sw = (Physics_SW*)phys->data;
415: const SWNode *x = (const SWNode*)xx;
416: PetscScalar u[2];
417: PetscReal h;
420: h = PetscRealPart(x->h);
421: Scale2(1./x->h,x->uh,u);
422: f[sw->functional.Height] = h;
423: f[sw->functional.Speed] = Norm2(u) + PetscSqrtReal(sw->gravity*h);
424: f[sw->functional.Energy] = 0.5*(Dot2(x->uh,u) + sw->gravity*PetscSqr (h));
425: return (0);
426: }
430: PetscErrorCode SetUpBC_SW(DM dm,Physics phys)431: {
433: const PetscInt wallids[] = {100,101,200,300};
435: DMAddBoundary (dm, PETSC_FALSE , "wall" , "Face Sets" , 0, 0, NULL, (void (*)()) PhysicsBoundary_SW_Wall, ALEN(wallids), wallids, phys);
436: return (0);
437: }
441: PetscErrorCode PhysicsCreate_SW(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)442: {
443: Physics_SW *sw;
447: phys->field_desc = PhysicsFields_SW;
448: phys->riemann = (PetscRiemannFunc) PhysicsRiemann_SW;
449: PetscNew (&sw);
450: phys->data = sw;
451: mod->setupbc = SetUpBC_SW;
453: PetscOptionsHead (PetscOptionsObject,"SW options" );
454: {
455: sw->gravity = 1.0;
456: PetscOptionsReal ("-sw_gravity" ,"Gravitational constant" ,"" ,sw->gravity,&sw->gravity,NULL);
457: }
458: PetscOptionsTail ();
459: phys->maxspeed = PetscSqrtReal(2.0*sw->gravity); /* Mach 1 for depth of 2 */
461: ModelSolutionSetDefault(mod,PhysicsSolution_SW,phys);
462: ModelFunctionalRegister(mod,"Height" ,&sw->functional.Height,PhysicsFunctional_SW,phys);
463: ModelFunctionalRegister(mod,"Speed" ,&sw->functional.Speed,PhysicsFunctional_SW,phys);
464: ModelFunctionalRegister(mod,"Energy" ,&sw->functional.Energy,PhysicsFunctional_SW,phys);
466: mod->bcs[0] = mod->bcs[1] = mod->bcs[2] = DM_BOUNDARY_GHOSTED;
468: return (0);
469: }
471: /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
472: /* An initialâvalue and selfâsimilar solutions of the compressible Euler equations */
473: /* Ravi Samtaney and D. I. Pullin */
474: /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
475: typedef enum {EULER_PAR_GAMMA,EULER_PAR_RHOR,EULER_PAR_AMACH,EULER_PAR_ITANA,EULER_PAR_SIZE} EulerParamIdx;
476: typedef enum {EULER_IV_SHOCK,EULER_SS_SHOCK,EULER_SHOCK_TUBE,EULER_LINEAR_WAVE} EulerType;
477: typedef struct {
478: PetscScalar vals[0];
479: PetscScalar r;
480: PetscScalar ru[DIM];
481: PetscScalar E;
482: } EulerNode;
483: typedef PetscErrorCode (*EquationOfState)(const PetscReal *, const EulerNode*, PetscScalar *) 484: typedef struct {
485: EulerType type;
486: PetscReal pars[EULER_PAR_SIZE];
487: EquationOfState sound;
488: struct {
489: PetscInt Density;
490: PetscInt Momentum;
491: PetscInt Energy;
492: PetscInt Pressure;
493: PetscInt Speed;
494: } monitor;
495: } Physics_Euler;
497: static const struct FieldDescription PhysicsFields_Euler[] = {{"Density" ,1},{"Momentum" ,DIM},{"Energy" ,1},{NULL,0}};
499: /* initial condition */
500: int initLinearWave(EulerNode *ux, const PetscScalar gamma, const PetscReal coord[], const PetscReal Lx) 503: PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)504: {
505: PetscInt i;
506: Physics phys = (Physics)ctx;
507: Physics_Euler *eu = (Physics_Euler*)phys->data;
508: EulerNode *uu = (EulerNode*)u;
509: PetscScalar p0,gamma,c;
511: if (time != 0.0) SETERRQ1 (mod->comm,PETSC_ERR_SUP,"No solution known for time %g" ,(double)time);
513: for (i=0; i<DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */
514: /* set E and rho */
515: gamma = eu->pars[EULER_PAR_GAMMA];
517: if (eu->type==EULER_IV_SHOCK || eu->type==EULER_SS_SHOCK) {
518: /******************* Euler Density Shock ********************/
519: /* On initialâvalue and selfâsimilar solutions of the compressible Euler equations */
520: /* Ravi Samtaney and D. I. Pullin */
521: /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
522: /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity, */
523: p0 = 1.;
524: if (x[0] < 0.0 + x[1]*eu->pars[EULER_PAR_ITANA]) {
525: if (x[0] < mod->bounds[0]*0.5) { /* left of shock (1) */
526: PetscScalar amach,rho,press,gas1,p1;
527: amach = eu->pars[EULER_PAR_AMACH];
528: rho = 1.;
529: press = p0;
530: p1 = press*(1.0+2.0*gamma/(gamma+1.0)*(amach*amach-1.0));
531: gas1 = (gamma-1.0)/(gamma+1.0);
532: uu->r = rho*(p1/press+gas1)/(gas1*p1/press+1.0);
533: uu->ru[0] = ((uu->r - rho)*sqrt(gamma*press/rho)*amach);
534: uu->E = p1/(gamma-1.0) + .5/uu->r*uu->ru[0]*uu->ru[0];
535: }
536: else { /* left of discontinuity (0) */
537: uu->r = 1.; /* rho = 1 */
538: uu->E = p0/(gamma-1.0);
539: }
540: }
541: else { /* right of discontinuity (2) */
542: uu->r = eu->pars[EULER_PAR_RHOR];
543: uu->E = p0/(gamma-1.0);
544: }
545: }
546: else if (eu->type==EULER_SHOCK_TUBE) {
547: /* For (x<x0) set (rho,u,p)=(8,0,10) and for (x>x0) set (rho,u,p)=(1,0,1). Choose x0 to the midpoint of the domain in the x-direction. */
548: if (x[0] < 0.0 ) {
549: uu->r = 8.;
550: uu->E = 10./(gamma-1.);
551: }
552: else {
553: uu->r = 1.;
554: uu->E = 1./(gamma-1.);
555: }
556: }
557: else if (eu->type==EULER_LINEAR_WAVE) {
558: initLinearWave( uu, gamma, x, mod->bounds[1] - mod->bounds[0]);
559: }
560: else SETERRQ1 (mod->comm,PETSC_ERR_SUP,"Unknown type %d" ,eu->type);
562: // set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main;
563: eu->sound(&gamma,uu,&c);
564: c = PetscAbsScalar(uu->ru[0]/uu->r) + c;
565: if (c > phys->maxspeed) phys->maxspeed = c;
567: return (0);
568: }
572: PetscErrorCode Pressure_PG(const PetscReal gamma,const EulerNode *x,PetscScalar *p)573: {
574: PetscScalar ru2;
576: ru2 = DotDIM(x->ru,x->ru);
577: (*p)=(x->E - 0.5*ru2/x->r)*(gamma - 1.0); /* (E - rho V^2/2)(gamma-1) = e rho (gamma-1) */
578: return (0);
579: }
583: PetscErrorCode SpeedOfSound_PG(const PetscReal *gamma, const EulerNode *x, PetscScalar *c)584: {
585: PetscScalar p;
588: Pressure_PG(*gamma,x,&p);
589: if (p<0.) SETERRQ1 (PETSC_COMM_WORLD ,PETSC_ERR_SUP,"negative pressure time %g -- NEED TO FIX!!!!!!" ,p);
590: /* pars[EULER_PAR_GAMMA] = heat capacity ratio */
591: (*c)=PetscSqrtScalar(*gamma * p / x->r);
592: return (0);
593: }
597: /*
598: * x = (rho,rho*(u_1),...,rho*e)^T
599: * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
600: *
601: * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
602: *
603: */
604: PetscErrorCode EulerFlux(Physics phys,const PetscReal *n,const EulerNode *x,EulerNode *f)605: {
606: Physics_Euler *eu = (Physics_Euler*)phys->data;
607: PetscScalar nu,p;
608: PetscInt i;
611: Pressure_PG(eu->pars[EULER_PAR_GAMMA],x,&p);
612: nu = DotDIM(x->ru,n);
613: f->r = nu; /* A rho u */
614: nu /= x->r; /* A u */
615: for (i=0; i<DIM; i++) f->ru[i] = nu * x->ru[i] + n[i]*p; /* r u^2 + p */
616: f->E = nu * (x->E + p); /* u(e+p) */
617: return (0);
618: }
620: /* PetscReal * => EulerNode* conversion */
623: PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)624: {
625: PetscInt i;
626: const EulerNode *xI = (const EulerNode*)a_xI;
627: EulerNode *xG = (EulerNode*)a_xG;
628: Physics phys = (Physics)ctx;
629: Physics_Euler *eu = (Physics_Euler*)phys->data;
631: xG->r = xI->r; /* ghost cell density - same */
632: xG->E = xI->E; /* ghost cell energy - same */
633: if (n[1] != 0.) { /* top and bottom */
634: xG->ru[0] = xI->ru[0]; /* copy tang to wall */
635: xG->ru[1] = -xI->ru[1]; /* reflect perp to t/b wall */
636: }
637: else { /* sides */
638: for (i=0; i<DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */
639: }
640: if (eu->type == EULER_LINEAR_WAVE) { // debug
641: // PetscPrintf (PETSC_COMM_WORLD ,"%s coord=%g,%g\n" ,__FUNCT__,c[0],c[1]);
642: }
643: return (0);
644: }
645: int godunovflux( const PetscScalar *ul, const PetscScalar *ur, PetscScalar *flux, const PetscReal *nn, const int *ndim, const PetscReal *gamma) 646: /* PetscReal * => EulerNode* conversion */
649: PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n,650: const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys) 651: {
652: Physics_Euler *eu = (Physics_Euler*)phys->data;
653: PetscReal cL,cR,speed,velL,velR,nn[DIM],s2;
654: PetscInt i;
655: PetscErrorCode ierr;
658: for (i=0,s2=0.; i<DIM; i++) {
659: nn[i] = n[i];
660: s2 += nn[i]*nn[i];
661: }
662: s2 = PetscSqrtReal(s2); /* |n|_2 = sum(n^2)^1/2 */
663: for (i=0.; i<DIM; i++) nn[i] /= s2;
664: if (0) { /* Rusanov */
665: const EulerNode *uL = (const EulerNode*)xL,*uR = (const EulerNode*)xR;
666: EulerNode fL,fR;
667: EulerFlux(phys,nn,uL,&fL);
668: EulerFlux(phys,nn,uR,&fR);
669: eu->sound(&eu->pars[EULER_PAR_GAMMA],uL,&cL);if (ierr) exit(13);
670: eu->sound(&eu->pars[EULER_PAR_GAMMA],uR,&cR);if (ierr) exit(14);
671: velL = DotDIM(uL->ru,nn)/uL->r;
672: velR = DotDIM(uR->ru,nn)/uR->r;
673: speed = PetscMax (PetscAbsScalar(velR) + cR,PetscAbsScalar(velL) + cL);
674: for (i=0; i<2+dim; i++) flux[i] = 0.5*((fL.vals[i]+fR.vals[i]) + speed*(xL[i] - xR[i]))*s2;
675: }
676: else {
677: int dim = DIM;
678: /* int iwave = */
679: godunovflux(xL, xR, flux, nn, &dim, &eu->pars[EULER_PAR_GAMMA]);
680: for (i=0; i<2+dim; i++) flux[i] *= s2;
681: }
682: PetscFunctionReturnVoid();
683: }
687: PetscErrorCode PhysicsFunctional_Euler(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)688: {
689: Physics phys = (Physics)ctx;
690: Physics_Euler *eu = (Physics_Euler*)phys->data;
691: const EulerNode *x = (const EulerNode*)xx;
692: PetscScalar p;
695: f[eu->monitor.Density] = x->r;
696: f[eu->monitor.Momentum] = NormDIM(x->ru);
697: f[eu->monitor.Energy] = x->E;
698: f[eu->monitor.Speed] = NormDIM(x->ru)/x->r;
699: Pressure_PG(eu->pars[EULER_PAR_GAMMA], x, &p);
700: f[eu->monitor.Pressure] = p;
701: return (0);
702: }
706: PetscErrorCode SetUpBC_Euler(DM dm,Physics phys)707: {
708: PetscErrorCode ierr;
709: Physics_Euler *eu = phys->data;
710: if (eu->type == EULER_LINEAR_WAVE) {
711: const PetscInt wallids[] = {100,101};
712: DMAddBoundary (dm, PETSC_FALSE , "wall" , "Face Sets" , 0, 0, NULL, (void (*)()) PhysicsBoundary_Euler_Wall, ALEN(wallids), wallids, phys);
713: }
714: else {
715: const PetscInt wallids[] = {100,101,200,300};
716: DMAddBoundary (dm, PETSC_FALSE , "wall" , "Face Sets" , 0, 0, NULL, (void (*)()) PhysicsBoundary_Euler_Wall, ALEN(wallids), wallids, phys);
717: }
718: return (0);
719: }
723: PetscErrorCode PhysicsCreate_Euler(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)724: {
725: Physics_Euler *eu;
726: PetscErrorCode ierr;
729: phys->field_desc = PhysicsFields_Euler;
730: phys->riemann = (PetscRiemannFunc) PhysicsRiemann_Euler_Godunov;
731: PetscNew (&eu);
732: phys->data = eu;
733: mod->setupbc = SetUpBC_Euler;
734: PetscOptionsHead (PetscOptionsObject,"Euler options" );
735: {
736: PetscReal alpha;
737: char type[64] = "linear_wave" ;
738: PetscBool is;
739: mod->bcs[0] = mod->bcs[1] = mod->bcs[2] = DM_BOUNDARY_GHOSTED;
740: eu->pars[EULER_PAR_GAMMA] = 1.4;
741: eu->pars[EULER_PAR_AMACH] = 2.02;
742: eu->pars[EULER_PAR_RHOR] = 3.0;
743: eu->pars[EULER_PAR_ITANA] = 0.57735026918963; /* angle of Euler self similar (SS) shock */
744: PetscOptionsReal ("-eu_gamma" ,"Heat capacity ratio" ,"" ,eu->pars[EULER_PAR_GAMMA],&eu->pars[EULER_PAR_GAMMA],NULL);
745: PetscOptionsReal ("-eu_amach" ,"Shock speed (Mach)" ,"" ,eu->pars[EULER_PAR_AMACH],&eu->pars[EULER_PAR_AMACH],NULL);
746: PetscOptionsReal ("-eu_rho2" ,"Density right of discontinuity" ,"" ,eu->pars[EULER_PAR_RHOR],&eu->pars[EULER_PAR_RHOR],NULL);
747: alpha = 60.;
748: PetscOptionsReal ("-eu_alpha" ,"Angle of discontinuity" ,"" ,alpha,&alpha,NULL);
749: if (alpha<=0. || alpha>90.) SETERRQ1 (PETSC_COMM_WORLD ,PETSC_ERR_SUP,"Alpha bust be > 0 and <= 90 (%g)" ,alpha);
750: eu->pars[EULER_PAR_ITANA] = 1./tan ( alpha * M_PI / 180.0 );
751: PetscOptionsString ("-eu_type" ,"Type of Euler test" ,"" ,type,type,sizeof (type),NULL);
752: PetscStrcmp (type,"linear_wave" , &is);
753: if (is) {
754: eu->type = EULER_LINEAR_WAVE;
755: mod->bcs[0] = mod->bcs[1] = mod->bcs[2] = DM_BOUNDARY_PERIODIC;
756: mod->bcs[1] = DM_BOUNDARY_GHOSTED; /* debug */
757: PetscPrintf (PETSC_COMM_WORLD ,"%s set Euler type: %s\n" ,__FUNCT__,"linear_wave" );
758: }
759: else {
760: if (DIM != 2) SETERRQ1 (PETSC_COMM_WORLD ,PETSC_ERR_SUP,"DIM must be 2 unless linear wave test %s" ,type);
761: PetscStrcmp (type,"iv_shock" , &is);
762: if (is) {
763: eu->type = EULER_IV_SHOCK;
764: PetscPrintf (PETSC_COMM_WORLD ,"%s set Euler type: %s\n" ,__FUNCT__,"iv_shock" );
765: }
766: else {
767: PetscStrcmp (type,"ss_shock" , &is);
768: if (is) {
769: eu->type = EULER_SS_SHOCK;
770: PetscPrintf (PETSC_COMM_WORLD ,"%s set Euler type: %s\n" ,__FUNCT__,"ss_shock" );
771: }
772: else {
773: PetscStrcmp (type,"shock_tube" , &is);
774: if (is) eu->type = EULER_SHOCK_TUBE;
775: else SETERRQ1 (PETSC_COMM_WORLD ,PETSC_ERR_SUP,"Unknown Euler type %s" ,type);
776: PetscPrintf (PETSC_COMM_WORLD ,"%s set Euler type: %s\n" ,__FUNCT__,"shock_tube" );
777: }
778: }
779: }
780: }
781: PetscOptionsTail ();
782: eu->sound = SpeedOfSound_PG;
783: phys->maxspeed = 0.; /* will get set in solution */
784: ModelSolutionSetDefault(mod,PhysicsSolution_Euler,phys);
785: ModelFunctionalRegister(mod,"Speed" ,&eu->monitor.Speed,PhysicsFunctional_Euler,phys);
786: ModelFunctionalRegister(mod,"Energy" ,&eu->monitor.Energy,PhysicsFunctional_Euler,phys);
787: ModelFunctionalRegister(mod,"Density" ,&eu->monitor.Density,PhysicsFunctional_Euler,phys);
788: ModelFunctionalRegister(mod,"Momentum" ,&eu->monitor.Momentum,PhysicsFunctional_Euler,phys);
789: ModelFunctionalRegister(mod,"Pressure" ,&eu->monitor.Pressure,PhysicsFunctional_Euler,phys);
791: return (0);
792: }
796: PetscErrorCode ConstructCellBoundary(DM dm, User user)797: {
798: const char *name = "Cell Sets" ;
799: const char *bdname = "split faces" ;
800: IS regionIS, innerIS;
801: const PetscInt *regions, *cells;
802: PetscInt numRegions, innerRegion, numCells, c;
803: PetscInt cStart, cEnd, cEndInterior, fStart, fEnd;
804: PetscBool hasLabel;
808: DMPlexGetHeightStratum (dm, 0, &cStart, &cEnd);
809: DMPlexGetHeightStratum (dm, 1, &fStart, &fEnd);
810: DMPlexGetHybridBounds (dm, &cEndInterior, NULL, NULL, NULL);
812: DMHasLabel (dm, name, &hasLabel);
813: if (!hasLabel) return (0);
814: DMGetLabelSize (dm, name, &numRegions);
815: if (numRegions != 2) return (0);
816: /* Get the inner id */
817: DMGetLabelIdIS (dm, name, ®ionIS);
818: ISGetIndices (regionIS, ®ions);
819: innerRegion = regions[0];
820: ISRestoreIndices (regionIS, ®ions);
821: ISDestroy (®ionIS);
822: /* Find the faces between cells in different regions, could call DMPlexCreateNeighborCSR() */
823: DMGetStratumIS (dm, name, innerRegion, &innerIS);
824: ISGetLocalSize (innerIS, &numCells);
825: ISGetIndices (innerIS, &cells);
826: DMCreateLabel (dm, bdname);
827: for (c = 0; c < numCells; ++c) {
828: const PetscInt cell = cells[c];
829: const PetscInt *faces;
830: PetscInt numFaces, f;
832: if ((cell < cStart) || (cell >= cEnd)) SETERRQ1 (PETSC_COMM_SELF , PETSC_ERR_LIB, "Got invalid point %d which is not a cell" , cell);
833: DMPlexGetConeSize (dm, cell, &numFaces);
834: DMPlexGetCone (dm, cell, &faces);
835: for (f = 0; f < numFaces; ++f) {
836: const PetscInt face = faces[f];
837: const PetscInt *neighbors;
838: PetscInt nC, regionA, regionB;
840: if ((face < fStart) || (face >= fEnd)) SETERRQ1 (PETSC_COMM_SELF , PETSC_ERR_LIB, "Got invalid point %d which is not a face" , face);
841: DMPlexGetSupportSize (dm, face, &nC);
842: if (nC != 2) continue ;
843: DMPlexGetSupport (dm, face, &neighbors);
844: if ((neighbors[0] >= cEndInterior) || (neighbors[1] >= cEndInterior)) continue ;
845: 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]);
846: 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]);
847: DMGetLabelValue (dm, name, neighbors[0], ®ionA);
848: DMGetLabelValue (dm, name, neighbors[1], ®ionB);
849: if (regionA < 0) SETERRQ2 (PetscObjectComm ((PetscObject )dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value" , name, neighbors[0]);
850: if (regionB < 0) SETERRQ2 (PetscObjectComm ((PetscObject )dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value" , name, neighbors[1]);
851: if (regionA != regionB) {
852: DMSetLabelValue (dm, bdname, faces[f], 1);
853: }
854: }
855: }
856: ISRestoreIndices (innerIS, &cells);
857: ISDestroy (&innerIS);
858: {
859: DMLabel label;
861: DMGetLabel (dm, bdname, &label);
862: DMLabelView (label, PETSC_VIEWER_STDOUT_WORLD );
863: }
864: return (0);
865: }
869: /* Right now, I have just added duplicate faces, which see both cells. We can
870: - Add duplicate vertices and decouple the face cones
871: - Disconnect faces from cells across the rotation gap
872: */
873: PetscErrorCode SplitFaces(DM *dmSplit, const char labelName[], User user)874: {
875: DM dm = *dmSplit, sdm;
876: PetscSF sfPoint, gsfPoint;
877: PetscSection coordSection, newCoordSection;
878: Vec coordinates;
879: IS idIS;
880: const PetscInt *ids;
881: PetscInt *newpoints;
882: PetscInt dim, depth, maxConeSize, maxSupportSize, numLabels, numGhostCells;
883: PetscInt numFS, fs, pStart, pEnd, p, cEnd, cEndInterior, vStart, vEnd, v, fStart, fEnd, newf, d, l;
884: PetscBool hasLabel;
888: DMHasLabel (dm, labelName, &hasLabel);
889: if (!hasLabel) return (0);
890: DMCreate (PetscObjectComm ((PetscObject )dm), &sdm);
891: DMSetType (sdm, DMPLEX );
892: DMGetDimension (dm, &dim);
893: DMSetDimension (sdm, dim);
895: DMGetLabelIdIS (dm, labelName, &idIS);
896: ISGetLocalSize (idIS, &numFS);
897: ISGetIndices (idIS, &ids);
899: user->numSplitFaces = 0;
900: for (fs = 0; fs < numFS; ++fs) {
901: PetscInt numBdFaces;
903: DMGetStratumSize (dm, labelName, ids[fs], &numBdFaces);
904: user->numSplitFaces += numBdFaces;
905: }
906: DMPlexGetChart (dm, &pStart, &pEnd);
907: pEnd += user->numSplitFaces;
908: DMPlexSetChart (sdm, pStart, pEnd);
909: DMPlexGetHybridBounds (dm, &cEndInterior, NULL, NULL, NULL);
910: DMPlexSetHybridBounds (sdm, cEndInterior, PETSC_DETERMINE , PETSC_DETERMINE , PETSC_DETERMINE );
911: DMPlexGetHeightStratum (dm, 0, NULL, &cEnd);
912: numGhostCells = cEnd - cEndInterior;
913: /* Set cone and support sizes */
914: DMPlexGetDepth (dm, &depth);
915: for (d = 0; d <= depth; ++d) {
916: DMPlexGetDepthStratum (dm, d, &pStart, &pEnd);
917: for (p = pStart; p < pEnd; ++p) {
918: PetscInt newp = p;
919: PetscInt size;
921: DMPlexGetConeSize (dm, p, &size);
922: DMPlexSetConeSize (sdm, newp, size);
923: DMPlexGetSupportSize (dm, p, &size);
924: DMPlexSetSupportSize (sdm, newp, size);
925: }
926: }
927: DMPlexGetHeightStratum (dm, 1, &fStart, &fEnd);
928: for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
929: IS faceIS;
930: const PetscInt *faces;
931: PetscInt numFaces, f;
933: DMGetStratumIS (dm, labelName, ids[fs], &faceIS);
934: ISGetLocalSize (faceIS, &numFaces);
935: ISGetIndices (faceIS, &faces);
936: for (f = 0; f < numFaces; ++f, ++newf) {
937: PetscInt size;
939: /* Right now I think that both faces should see both cells */
940: DMPlexGetConeSize (dm, faces[f], &size);
941: DMPlexSetConeSize (sdm, newf, size);
942: DMPlexGetSupportSize (dm, faces[f], &size);
943: DMPlexSetSupportSize (sdm, newf, size);
944: }
945: ISRestoreIndices (faceIS, &faces);
946: ISDestroy (&faceIS);
947: }
948: DMSetUp (sdm);
949: /* Set cones and supports */
950: DMPlexGetMaxSizes (dm, &maxConeSize, &maxSupportSize);
951: PetscMalloc (PetscMax (maxConeSize, maxSupportSize) * sizeof (PetscInt ), &newpoints);
952: DMPlexGetChart (dm, &pStart, &pEnd);
953: for (p = pStart; p < pEnd; ++p) {
954: const PetscInt *points, *orientations;
955: PetscInt size, i, newp = p;
957: DMPlexGetConeSize (dm, p, &size);
958: DMPlexGetCone (dm, p, &points);
959: DMPlexGetConeOrientation (dm, p, &orientations);
960: for (i = 0; i < size; ++i) newpoints[i] = points[i];
961: DMPlexSetCone (sdm, newp, newpoints);
962: DMPlexSetConeOrientation (sdm, newp, orientations);
963: DMPlexGetSupportSize (dm, p, &size);
964: DMPlexGetSupport (dm, p, &points);
965: for (i = 0; i < size; ++i) newpoints[i] = points[i];
966: DMPlexSetSupport (sdm, newp, newpoints);
967: }
968: PetscFree (newpoints);
969: for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
970: IS faceIS;
971: const PetscInt *faces;
972: PetscInt numFaces, f;
974: DMGetStratumIS (dm, labelName, ids[fs], &faceIS);
975: ISGetLocalSize (faceIS, &numFaces);
976: ISGetIndices (faceIS, &faces);
977: for (f = 0; f < numFaces; ++f, ++newf) {
978: const PetscInt *points;
980: DMPlexGetCone (dm, faces[f], &points);
981: DMPlexSetCone (sdm, newf, points);
982: DMPlexGetSupport (dm, faces[f], &points);
983: DMPlexSetSupport (sdm, newf, points);
984: }
985: ISRestoreIndices (faceIS, &faces);
986: ISDestroy (&faceIS);
987: }
988: ISRestoreIndices (idIS, &ids);
989: ISDestroy (&idIS);
990: DMPlexStratify (sdm);
991: /* Convert coordinates */
992: DMPlexGetDepthStratum (dm, 0, &vStart, &vEnd);
993: DMGetCoordinateSection (dm, &coordSection);
994: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), &newCoordSection);
995: PetscSectionSetNumFields (newCoordSection, 1);
996: PetscSectionSetFieldComponents (newCoordSection, 0, dim);
997: PetscSectionSetChart (newCoordSection, vStart, vEnd);
998: for (v = vStart; v < vEnd; ++v) {
999: PetscSectionSetDof (newCoordSection, v, dim);
1000: PetscSectionSetFieldDof (newCoordSection, v, 0, dim);
1001: }
1002: PetscSectionSetUp (newCoordSection);
1003: DMSetCoordinateSection (sdm, PETSC_DETERMINE , newCoordSection);
1004: PetscSectionDestroy (&newCoordSection); /* relinquish our reference */
1005: DMGetCoordinatesLocal (dm, &coordinates);
1006: DMSetCoordinatesLocal (sdm, coordinates);
1007: /* Convert labels */
1008: DMGetNumLabels (dm, &numLabels);
1009: for (l = 0; l < numLabels; ++l) {
1010: const char *lname;
1011: PetscBool isDepth;
1013: DMGetLabelName (dm, l, &lname);
1014: PetscStrcmp (lname, "depth" , &isDepth);
1015: if (isDepth) continue ;
1016: DMCreateLabel (sdm, lname);
1017: DMGetLabelIdIS (dm, lname, &idIS);
1018: ISGetLocalSize (idIS, &numFS);
1019: ISGetIndices (idIS, &ids);
1020: for (fs = 0; fs < numFS; ++fs) {
1021: IS pointIS;
1022: const PetscInt *points;
1023: PetscInt numPoints;
1025: DMGetStratumIS (dm, lname, ids[fs], &pointIS);
1026: ISGetLocalSize (pointIS, &numPoints);
1027: ISGetIndices (pointIS, &points);
1028: for (p = 0; p < numPoints; ++p) {
1029: PetscInt newpoint = points[p];
1031: DMSetLabelValue (sdm, lname, newpoint, ids[fs]);
1032: }
1033: ISRestoreIndices (pointIS, &points);
1034: ISDestroy (&pointIS);
1035: }
1036: ISRestoreIndices (idIS, &ids);
1037: ISDestroy (&idIS);
1038: }
1039: /* Convert pointSF */
1040: const PetscSFNode *remotePoints;
1041: PetscSFNode *gremotePoints;
1042: const PetscInt *localPoints;
1043: PetscInt *glocalPoints,*newLocation,*newRemoteLocation;
1044: PetscInt numRoots, numLeaves;
1045: PetscMPIInt numProcs;
1047: MPI_Comm_size (PetscObjectComm ((PetscObject )dm), &numProcs);
1048: DMGetPointSF (dm, &sfPoint);
1049: DMGetPointSF (sdm, &gsfPoint);
1050: DMPlexGetChart (dm,&pStart,&pEnd);
1051: PetscSFGetGraph (sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1052: if (numRoots >= 0) {
1053: PetscMalloc2 (numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
1054: for (l=0; l<numRoots; l++) newLocation[l] = l; /* + (l >= cEnd ? numGhostCells : 0); */
1055: PetscSFBcastBegin (sfPoint, MPIU_INT, newLocation, newRemoteLocation);
1056: PetscSFBcastEnd (sfPoint, MPIU_INT, newLocation, newRemoteLocation);
1057: PetscMalloc1 (numLeaves, &glocalPoints);
1058: PetscMalloc1 (numLeaves, &gremotePoints);
1059: for (l = 0; l < numLeaves; ++l) {
1060: glocalPoints[l] = localPoints[l]; /* localPoints[l] >= cEnd ? localPoints[l] + numGhostCells : localPoints[l]; */
1061: gremotePoints[l].rank = remotePoints[l].rank;
1062: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
1063: }
1064: PetscFree2 (newLocation,newRemoteLocation);
1065: PetscSFSetGraph (gsfPoint, numRoots+numGhostCells, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
1066: }
1067: DMDestroy (dmSplit);
1068: *dmSplit = sdm;
1069: return (0);
1070: }
1074: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)1075: {
1076: PetscSF sfPoint;
1077: PetscSection coordSection;
1078: Vec coordinates;
1079: PetscSection sectionCell;
1080: PetscScalar *part;
1081: PetscInt cStart, cEnd, c;
1082: PetscMPIInt rank;
1086: DMGetCoordinateSection (dm, &coordSection);
1087: DMGetCoordinatesLocal (dm, &coordinates);
1088: DMClone (dm, dmCell);
1089: DMGetPointSF (dm, &sfPoint);
1090: DMSetPointSF (*dmCell, sfPoint);
1091: DMSetCoordinateSection (*dmCell, PETSC_DETERMINE , coordSection);
1092: DMSetCoordinatesLocal (*dmCell, coordinates);
1093: MPI_Comm_rank (PetscObjectComm ((PetscObject )dm), &rank);
1094: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), §ionCell);
1095: DMPlexGetHeightStratum (*dmCell, 0, &cStart, &cEnd);
1096: PetscSectionSetChart (sectionCell, cStart, cEnd);
1097: for (c = cStart; c < cEnd; ++c) {
1098: PetscSectionSetDof (sectionCell, c, 1);
1099: }
1100: PetscSectionSetUp (sectionCell);
1101: DMSetDefaultSection (*dmCell, sectionCell);
1102: PetscSectionDestroy (§ionCell);
1103: DMCreateLocalVector (*dmCell, partition);
1104: PetscObjectSetName ((PetscObject )*partition, "partition" );
1105: VecGetArray (*partition, &part);
1106: for (c = cStart; c < cEnd; ++c) {
1107: PetscScalar *p;
1109: DMPlexPointLocalRef (*dmCell, c, part, &p);
1110: p[0] = rank;
1111: }
1112: VecRestoreArray (*partition, &part);
1113: return (0);
1114: }
1118: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)1119: {
1120: DM dmMass, dmFace, dmCell, dmCoord;
1121: PetscSection coordSection;
1122: Vec coordinates, facegeom, cellgeom;
1123: PetscSection sectionMass;
1124: PetscScalar *m;
1125: const PetscScalar *fgeom, *cgeom, *coords;
1126: PetscInt vStart, vEnd, v;
1127: PetscErrorCode ierr;
1130: DMGetCoordinateSection (dm, &coordSection);
1131: DMGetCoordinatesLocal (dm, &coordinates);
1132: DMClone (dm, &dmMass);
1133: DMSetCoordinateSection (dmMass, PETSC_DETERMINE , coordSection);
1134: DMSetCoordinatesLocal (dmMass, coordinates);
1135: PetscSectionCreate (PetscObjectComm ((PetscObject )dm), §ionMass);
1136: DMPlexGetDepthStratum (dm, 0, &vStart, &vEnd);
1137: PetscSectionSetChart (sectionMass, vStart, vEnd);
1138: for (v = vStart; v < vEnd; ++v) {
1139: PetscInt numFaces;
1141: DMPlexGetSupportSize (dmMass, v, &numFaces);
1142: PetscSectionSetDof (sectionMass, v, numFaces*numFaces);
1143: }
1144: PetscSectionSetUp (sectionMass);
1145: DMSetDefaultSection (dmMass, sectionMass);
1146: PetscSectionDestroy (§ionMass);
1147: DMGetLocalVector (dmMass, massMatrix);
1148: VecGetArray (*massMatrix, &m);
1149: DMPlexTSGetGeometryFVM (dm, &facegeom, &cellgeom, NULL);
1150: VecGetDM (facegeom, &dmFace);
1151: VecGetArrayRead (facegeom, &fgeom);
1152: VecGetDM (cellgeom, &dmCell);
1153: VecGetArrayRead (cellgeom, &cgeom);
1154: DMGetCoordinateDM (dm, &dmCoord);
1155: VecGetArrayRead (coordinates, &coords);
1156: for (v = vStart; v < vEnd; ++v) {
1157: const PetscInt *faces;
1158: const PetscFVFaceGeom *fgA, *fgB, *cg;
1159: const PetscScalar *vertex;
1160: PetscInt numFaces, sides[2], f, g;
1162: DMPlexPointLocalRead (dmCoord, v, coords, &vertex);
1163: DMPlexGetSupportSize (dmMass, v, &numFaces);
1164: DMPlexGetSupport (dmMass, v, &faces);
1165: for (f = 0; f < numFaces; ++f) {
1166: sides[0] = faces[f];
1167: DMPlexPointLocalRead (dmFace, faces[f], fgeom, &fgA);
1168: for (g = 0; g < numFaces; ++g) {
1169: const PetscInt *cells = NULL;;
1170: PetscReal area = 0.0;
1171: PetscInt numCells;
1173: sides[1] = faces[g];
1174: DMPlexPointLocalRead (dmFace, faces[g], fgeom, &fgB);
1175: DMPlexGetJoin (dmMass, 2, sides, &numCells, &cells);
1176: if (numCells != 1) SETERRQ (PETSC_COMM_SELF , PETSC_ERR_LIB, "Invalid join for faces" );
1177: DMPlexPointLocalRead (dmCell, cells[0], cgeom, &cg);
1178: area += PetscAbsScalar((vertex[0] - cg->centroid[0])*(fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1])*(fgA->centroid[0] - cg->centroid[0]));
1179: area += PetscAbsScalar((vertex[0] - cg->centroid[0])*(fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1])*(fgB->centroid[0] - cg->centroid[0]));
1180: m[f*numFaces+g] = Dot2(fgA->normal, fgB->normal)*area*0.5;
1181: DMPlexRestoreJoin (dmMass, 2, sides, &numCells, &cells);
1182: }
1183: }
1184: }
1185: VecRestoreArrayRead (facegeom, &fgeom);
1186: VecRestoreArrayRead (cellgeom, &cgeom);
1187: VecRestoreArrayRead (coordinates, &coords);
1188: VecRestoreArray (*massMatrix, &m);
1189: DMDestroy (&dmMass);
1190: return (0);
1191: }
1195: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
1196: PetscErrorCode ModelSolutionSetDefault(Model mod,SolutionFunction func,void *ctx)1197: {
1199: mod->solution = func;
1200: mod->solutionctx = ctx;
1201: return (0);
1202: }
1206: PetscErrorCode ModelFunctionalRegister(Model mod,const char *name,PetscInt *offset,FunctionalFunction func,void *ctx)1207: {
1209: FunctionalLink link,*ptr;
1210: PetscInt lastoffset = -1;
1213: for (ptr=&mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
1214: PetscNew (&link);
1215: PetscStrallocpy (name,&link->name);
1216: link->offset = lastoffset + 1;
1217: link->func = func;
1218: link->ctx = ctx;
1219: link->next = NULL;
1220: *ptr = link;
1221: *offset = link->offset;
1222: return (0);
1223: }
1227: PetscErrorCode ModelFunctionalSetFromOptions(Model mod,PetscOptionItems *PetscOptionsObject)1228: {
1230: PetscInt i,j;
1231: FunctionalLink link;
1232: char *names[256];
1235: mod->numMonitored = ALEN(names);
1236: PetscOptionsStringArray ("-monitor" ,"list of functionals to monitor" ,"" ,names,&mod->numMonitored,NULL);
1237: /* Create list of functionals that will be computed somehow */
1238: PetscMalloc1 (mod->numMonitored,&mod->functionalMonitored);
1239: /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
1240: PetscMalloc1 (mod->numMonitored,&mod->functionalCall);
1241: mod->numCall = 0;
1242: for (i=0; i<mod->numMonitored; i++) {
1243: for (link=mod->functionalRegistry; link; link=link->next) {
1244: PetscBool match;
1245: PetscStrcasecmp (names[i],link->name,&match);
1246: if (match) break ;
1247: }
1248: if (!link) SETERRQ1 (mod->comm,PETSC_ERR_USER,"No known functional '%s'" ,names[i]);
1249: mod->functionalMonitored[i] = link;
1250: for (j=0; j<i; j++) {
1251: if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
1252: }
1253: mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
1254: next_name: 1255: PetscFree (names[i]);
1256: }
1258: /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
1259: mod->maxComputed = -1;
1260: for (link=mod->functionalRegistry; link; link=link->next) {
1261: for (i=0; i<mod->numCall; i++) {
1262: FunctionalLink call = mod->functionalCall[i];
1263: if (link->func == call->func && link->ctx == call->ctx) {
1264: mod->maxComputed = PetscMax (mod->maxComputed,link->offset);
1265: }
1266: }
1267: }
1268: return (0);
1269: }
1273: PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)1274: {
1276: FunctionalLink l,next;
1279: if (!link) return (0);
1280: l = *link;
1281: *link = NULL;
1282: for (; l; l=next) {
1283: next = l->next;
1284: PetscFree (l->name);
1285: PetscFree (l);
1286: }
1287: return (0);
1288: }
1292: /* put the solution callback into a functional callback */
1293: PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)1294: {
1295: Model mod;
1298: mod = (Model) modctx;
1299: (*mod->solution)(mod, time, x, u, mod->solutionctx);
1300: return (0);
1301: }
1305: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)1306: {
1307: PetscErrorCode (*func[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
1308: void *ctx[1];
1309: Model mod = user->model;
1310: PetscErrorCode ierr;
1313: func[0] = SolutionFunctional;
1314: ctx[0] = (void *) mod;
1315: DMProjectFunction (dm,0.0,func,ctx,INSERT_ALL_VALUES,X);
1316: return (0);
1317: }
1321: PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)1322: {
1326: PetscViewerCreate (PetscObjectComm ((PetscObject )dm), viewer);
1327: PetscViewerSetType (*viewer, PETSCVIEWERVTK);
1328: PetscViewerFileSetName (*viewer, filename);
1329: return (0);
1330: }
1334: PetscErrorCode MonitorVTK(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx)1335: {
1336: User user = (User)ctx;
1337: DM dm;
1338: Vec cellgeom;
1339: PetscViewer viewer;
1340: char filename[PETSC_MAX_PATH_LEN],*ftable = NULL;
1341: PetscReal xnorm;
1342: PetscInt cEndInterior;
1346: PetscObjectSetName ((PetscObject ) X, "u" );
1347: VecGetDM (X,&dm);
1348: DMPlexTSGetGeometryFVM (dm, NULL, &cellgeom, NULL);
1349: VecNorm (X,NORM_INFINITY ,&xnorm);
1351: if (stepnum >= 0) { /* No summary for final time */
1352: Model mod = user->model;
1353: PetscInt c,cStart,cEnd,fcount,i;
1354: size_t ftableused,ftablealloc;
1355: const PetscScalar *cgeom,*x;
1356: DM dmCell;
1357: PetscReal *fmin,*fmax,*fintegral,*ftmp;
1358: fcount = mod->maxComputed+1;
1359: PetscMalloc4 (fcount,&fmin,fcount,&fmax,fcount,&fintegral,fcount,&ftmp);
1360: for (i=0; i<fcount; i++) {
1361: fmin[i] = PETSC_MAX_REAL;
1362: fmax[i] = PETSC_MIN_REAL;
1363: fintegral[i] = 0;
1364: }
1365: VecGetDM (cellgeom,&dmCell);
1366: DMPlexGetHybridBounds (dmCell, &cEndInterior, NULL, NULL, NULL);
1367: DMPlexGetHeightStratum (dmCell,0,&cStart,&cEnd);
1368: VecGetArrayRead (cellgeom,&cgeom);
1369: VecGetArrayRead (X,&x);
1370: for (c = cStart; c < cEndInterior; ++c) {
1371: const PetscFVCellGeom *cg;
1372: const PetscScalar *cx;
1373: /* not that these two routines as currently implemented work for any dm with a
1374: * defaultSection/defaultGlobalSection */
1375: DMPlexPointLocalRead (dmCell,c,cgeom,&cg);
1376: DMPlexPointGlobalRead (dm,c,x,&cx);
1377: if (!cx) continue ; /* not a global cell */
1378: for (i=0; i<mod->numCall; i++) {
1379: FunctionalLink flink = mod->functionalCall[i];
1380: (*flink->func)(mod,time,cg->centroid,cx,ftmp,flink->ctx);
1381: }
1382: for (i=0; i<fcount; i++) {
1383: fmin[i] = PetscMin (fmin[i],ftmp[i]);
1384: fmax[i] = PetscMax (fmax[i],ftmp[i]);
1385: fintegral[i] += cg->volume * ftmp[i];
1386: }
1387: }
1388: VecRestoreArrayRead (cellgeom,&cgeom);
1389: VecRestoreArrayRead (X,&x);
1390: MPI_Allreduce (MPI_IN_PLACE,fmin,fcount,MPIU_REAL,MPIU_MIN,PetscObjectComm ((PetscObject )ts));
1391: MPI_Allreduce (MPI_IN_PLACE,fmax,fcount,MPIU_REAL,MPIU_MAX,PetscObjectComm ((PetscObject )ts));
1392: MPI_Allreduce (MPI_IN_PLACE,fintegral,fcount,MPIU_REAL,MPIU_SUM,PetscObjectComm ((PetscObject )ts));
1394: ftablealloc = fcount * 100;
1395: ftableused = 0;
1396: PetscMalloc1 (ftablealloc,&ftable);
1397: for (i=0; i<mod->numMonitored; i++) {
1398: size_t countused;
1399: char buffer[256],*p;
1400: FunctionalLink flink = mod->functionalMonitored[i];
1401: PetscInt id = flink->offset;
1402: if (i % 3) {
1403: PetscMemcpy (buffer," " ,2);
1404: p = buffer + 2;
1405: } else if (i) {
1406: char newline[] = "\n" ;
1407: PetscMemcpy (buffer,newline,sizeof newline-1);
1408: p = buffer + sizeof newline - 1;
1409: } else {
1410: p = buffer;
1411: }
1412: PetscSNPrintfCount (p,sizeof buffer-(p-buffer),"%12s [%10.7g,%10.7g] int %10.7g" ,&countused,flink->name,(double)fmin[id],(double)fmax[id],(double)fintegral[id]);
1413: countused += p - buffer;
1414: if (countused > ftablealloc-ftableused-1) { /* reallocate */
1415: char *ftablenew;
1416: ftablealloc = 2*ftablealloc + countused;
1417: PetscMalloc (ftablealloc,&ftablenew);
1418: PetscMemcpy (ftablenew,ftable,ftableused);
1419: PetscFree (ftable);
1420: ftable = ftablenew;
1421: }
1422: PetscMemcpy (ftable+ftableused,buffer,countused);
1423: ftableused += countused;
1424: ftable[ftableused] = 0;
1425: }
1426: PetscFree4 (fmin,fmax,fintegral,ftmp);
1428: PetscPrintf (PetscObjectComm ((PetscObject )ts),"% 3D time %8.4g |x| %8.4g %s\n" ,stepnum,(double)time,(double)xnorm,ftable ? ftable : "" );
1429: PetscFree (ftable);
1430: }
1431: if (user->vtkInterval < 1) return (0);
1432: if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1433: if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
1434: TSGetTimeStepNumber (ts,&stepnum);
1435: }
1436: PetscSNPrintf (filename,sizeof filename,"ex11-%03D.vtu" ,stepnum);
1437: OutputVTK(dm,filename,&viewer);
1438: VecView (X,viewer);
1439: PetscViewerDestroy (&viewer);
1440: }
1441: return (0);
1442: }
1446: 1447: {
1448: MPI_Comm comm;
1449: PetscDS prob;
1450: PetscFV fvm;
1451: User user;
1452: Model mod;
1453: Physics phys;
1454: DM dm;
1455: PetscReal ftime, cfl, dt, minRadius;
1456: PetscInt dim, nsteps;
1457: TS ts;
1458: TSConvergedReason reason;
1459: Vec X;
1460: PetscViewer viewer;
1461: PetscBool vtkCellGeom, splitFaces;
1462: PetscInt overlap;
1463: char filename[PETSC_MAX_PATH_LEN] = "sevenside.exo" ;
1464: char physname[256] = "advect" ;
1465: PetscErrorCode ierr;
1467: PetscInitialize (&argc, &argv, (char*) 0, help);
1468: comm = PETSC_COMM_WORLD ;
1470: PetscNew (&user);
1471: PetscNew (&user->model);
1472: PetscNew (&user->model->physics);
1473: mod = user->model;
1474: phys = mod->physics;
1475: mod->comm = comm;
1477: /* Register physical models to be available on the command line */
1478: PetscFunctionListAdd (&PhysicsList,"advect" ,PhysicsCreate_Advect);
1479: PetscFunctionListAdd (&PhysicsList,"sw" ,PhysicsCreate_SW);
1480: PetscFunctionListAdd (&PhysicsList,"euler" ,PhysicsCreate_Euler);
1482: PetscOptionsBegin (comm,NULL,"Unstructured Finite Volume Mesh Options" ,"" );
1483: {
1484: cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1485: PetscOptionsReal ("-ufv_cfl" ,"CFL number per step" ,"" ,cfl,&cfl,NULL);
1486: PetscOptionsString ("-f" ,"Exodus.II filename to read" ,"" ,filename,filename,sizeof (filename),NULL);
1487: splitFaces = PETSC_FALSE ;
1488: PetscOptionsBool ("-ufv_split_faces" ,"Split faces between cell sets" ,"" ,splitFaces,&splitFaces,NULL);
1489: overlap = 1;
1490: PetscOptionsInt ("-ufv_mesh_overlap" ,"Number of cells to overlap partitions" ,"" ,overlap,&overlap,NULL);
1491: user->vtkInterval = 1;
1492: PetscOptionsInt ("-ufv_vtk_interval" ,"VTK output interval (0 to disable)" ,"" ,user->vtkInterval,&user->vtkInterval,NULL);
1493: vtkCellGeom = PETSC_FALSE ;
1494: PetscOptionsBool ("-ufv_vtk_cellgeom" ,"Write cell geometry (for debugging)" ,"" ,vtkCellGeom,&vtkCellGeom,NULL);
1495: }
1496: PetscOptionsEnd ();
1498: PetscOptionsBegin (comm,NULL,"Unstructured Finite Volume Physics Options" ,"" );
1499: {
1500: PetscErrorCode (*physcreate)(Model,Physics,PetscOptionItems*);
1501: PetscOptionsFList ("-physics" ,"Physics module to solve" ,"" ,PhysicsList,physname,physname,sizeof physname,NULL);
1502: PetscFunctionListFind (PhysicsList,physname,&physcreate);
1503: PetscMemzero (phys,sizeof (struct _n_Physics ));
1504: (*physcreate)(mod,phys,PetscOptionsObject);
1505: /* Count number of fields and dofs */
1506: for (phys->nfields=0,phys->dof=0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1507: if (phys->dof <= 0) SETERRQ1 (comm,PETSC_ERR_ARG_WRONGSTATE,"Physics '%s' did not set dof" ,physname);
1508: ModelFunctionalSetFromOptions(mod,PetscOptionsObject);
1509: }
1510: PetscOptionsEnd ();
1512: {
1513: size_t len,i;
1514: for (i = 0; i < DIM; i++) { mod->bounds[2*i] = 0.; mod->bounds[2*i+1] = 1.;};
1515: PetscStrlen (filename,&len);
1516: dim = DIM;
1517: if (!len) { /* a null name means just do a hex box */
1518: PetscInt cells[3] = {1, 1, 1}; /* coarse mesh is one cell; refine from there */
1519: PetscBool flg1, flg2;
1520: PetscInt nret1 = DIM;
1521: PetscInt nret2 = 2*DIM;
1522: PetscOptionsBegin (comm,NULL,"Rectangular mesh options" ,"" );
1523: PetscOptionsIntArray ("-grid_size" ,"number of cells in each direction" ,"" ,cells,&nret1,&flg1);
1524: PetscOptionsRealArray ("-grid_bounds" ,"bounds of the mesh in each direction (e.g., x_min,x_max,y_min,y_max" ,"" ,mod->bounds,&nret2,&flg2);
1525: PetscOptionsEnd ();
1526: if (flg1) {
1527: dim = nret1;
1528: if (dim != DIM) SETERRQ1 (comm,PETSC_ERR_ARG_SIZ,"Dim wrong size %D in -grid_size" ,dim);
1529: }
1530: DMPlexCreateHexBoxMesh (comm, dim, cells, mod->bcs[0], mod->bcs[1], mod->bcs[2], &dm);
1531: if (flg2) {
1532: PetscInt dimEmbed, i;
1533: PetscInt nCoords;
1534: PetscScalar *coords;
1535: Vec coordinates;
1537: DMGetCoordinatesLocal (dm,&coordinates);
1538: DMGetCoordinateDim (dm,&dimEmbed);
1539: VecGetLocalSize (coordinates,&nCoords);
1540: if (nCoords % dimEmbed) SETERRQ (PETSC_COMM_SELF ,PETSC_ERR_ARG_SIZ,"Coordinate vector the wrong size" );
1541: VecGetArray (coordinates,&coords);
1542: for (i = 0; i < nCoords; i += dimEmbed) {
1543: PetscInt j;
1545: PetscScalar *coord = &coords[i];
1546: for (j = 0; j < dimEmbed; j++) {
1547: coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1548: if (dim==2 && cells[1]==1 && j==0 && 0) {
1549: if (cells[0]==2 && coord[j]==mod->bounds[3] && i==8) {
1550: coord[j] *= (1.57735026918963); /* hack to get 60 deg skewed mesh */
1551: }
1552: else if (cells[0]==3) {
1553: if (i==2 || i==10) coord[j] = mod->bounds[1]/4.;
1554: else if (i==4) coord[j] = mod->bounds[1]/2.;
1555: else if (i==12) coord[j] = 1.57735026918963*mod->bounds[1]/2.;
1556: }
1557: }
1558: }
1559: }
1560: VecRestoreArray (coordinates,&coords);
1561: DMSetCoordinatesLocal (dm,coordinates);
1562: }
1563: }
1564: else {
1565: DMPlexCreateFromFile (comm, filename, PETSC_TRUE , &dm);
1566: }
1567: }
1568: DMViewFromOptions(dm, NULL, "-dm_view" );
1569: DMGetDimension (dm, &dim);
1571: /* set up BCs, functions, tags */
1572: DMCreateLabel (dm, "Face Sets" );
1573: (*mod->setupbc)(dm,phys);
1575: {
1576: DM dmDist;
1578: DMPlexSetAdjacencyUseCone (dm, PETSC_TRUE );
1579: DMPlexSetAdjacencyUseClosure (dm, PETSC_FALSE );
1580: DMPlexDistribute (dm, overlap, NULL, &dmDist);
1581: if (dmDist) {
1582: DMDestroy (&dm);
1583: dm = dmDist;
1584: }
1585: }
1587: DMSetFromOptions (dm);
1589: {
1590: DM gdm;
1592: DMPlexConstructGhostCells (dm, NULL, NULL, &gdm);
1593: DMDestroy (&dm);
1594: dm = gdm;
1595: DMViewFromOptions(dm, NULL, "-dm_view" );
1596: }
1597: if (splitFaces) {ConstructCellBoundary(dm, user);}
1598: SplitFaces(&dm, "split faces" , user);
1600: {
1601: char convType[256];
1602: PetscBool flg;
1604: PetscOptionsBegin (comm, "" , "Mesh conversion options" , "DMPLEX " );
1605: PetscOptionsFList ("-dm_plex_convert_type" ,"Convert DMPlex to another format" ,"ex12" ,DMList,DMPLEX ,convType,256,&flg);
1606: PetscOptionsEnd ();
1607: if (flg) {
1608: DM dmConv;
1610: DMConvert (dm,convType,&dmConv);
1611: if (dmConv) {
1612: DMViewFromOptions(dmConv, NULL, "-dm_conv_view" );
1613: DMDestroy (&dm);
1614: dm = dmConv;
1615: DMSetFromOptions (dm);
1616: }
1617: }
1618: }
1620: PetscFVCreate (comm, &fvm);
1621: PetscFVSetFromOptions (fvm);
1622: PetscFVSetNumComponents (fvm, phys->dof);
1623: PetscFVSetSpatialDimension (fvm, dim);
1624: PetscObjectSetName ((PetscObject ) fvm,"" );
1625: {
1626: PetscInt f, dof;
1627: for (f=0,dof=0; f < phys->nfields; f++) {
1628: PetscInt newDof = phys->field_desc[f].dof;
1630: if (newDof == 1) {
1631: PetscFVSetComponentName (fvm,dof,phys->field_desc[f].name);
1632: }
1633: else {
1634: PetscInt j;
1636: for (j = 0; j < newDof; j++) {
1637: char compName[256] = "Unknown" ;
1639: PetscSNPrintf (compName,sizeof (compName),"%s_%d" ,phys->field_desc[f].name,j);
1640: PetscFVSetComponentName (fvm,dof+j,compName);
1641: }
1642: }
1643: dof += newDof;
1644: }
1645: }
1646: DMGetDS (dm, &prob);
1647: /* FV is now structured with one field having all physics as components */
1648: PetscDSAddDiscretization (prob, (PetscObject ) fvm);
1649: PetscDSSetRiemannSolver (prob, 0, user->model->physics->riemann);
1650: PetscDSSetContext(prob, 0, user->model->physics);
1652: TSCreate (comm, &ts);
1653: TSSetType (ts, TSSSP );
1654: TSSetDM (ts, dm);
1655: TSMonitorSet (ts,MonitorVTK,user,NULL);
1656: DMTSSetBoundaryLocal (dm, DMPlexTSComputeBoundary , user);
1657: DMTSSetRHSFunctionLocal (dm, DMPlexTSComputeRHSFunctionFVM , user);
1659: DMCreateGlobalVector (dm, &X);
1660: PetscObjectSetName ((PetscObject ) X, "solution" );
1661: SetInitialCondition(dm, X, user);
1662: if (vtkCellGeom) {
1663: DM dmCell;
1664: Vec cellgeom, partition;
1666: DMPlexTSGetGeometryFVM (dm, NULL, &cellgeom, NULL);
1667: OutputVTK(dm, "ex11-cellgeom.vtk" , &viewer);
1668: VecView (cellgeom, viewer);
1669: PetscViewerDestroy (&viewer);
1670: CreatePartitionVec(dm, &dmCell, &partition);
1671: OutputVTK(dmCell, "ex11-partition.vtk" , &viewer);
1672: VecView (partition, viewer);
1673: PetscViewerDestroy (&viewer);
1674: VecDestroy (&partition);
1675: DMDestroy (&dmCell);
1676: }
1678: DMPlexTSGetGeometryFVM (dm, NULL, NULL, &minRadius);
1679: TSSetDuration (ts,1000,2.0);
1680: /* collect max maxspeed from all processes -- todo */
1681: mod->maxspeed = phys->maxspeed;
1682: if (mod->maxspeed <= 0) SETERRQ1 (comm,PETSC_ERR_ARG_WRONGSTATE,"Physics '%s' did not set maxspeed" ,physname);
1683: TSSetExactFinalTime (ts,TS_EXACTFINALTIME_STEPOVER);
1684: dt = cfl * minRadius / user->model->maxspeed;
1685: TSSetInitialTimeStep (ts,0.0,dt);
1686: TSSetFromOptions (ts);
1687: TSSolve (ts,X);
1688: TSGetSolveTime (ts,&ftime);
1689: TSGetTimeStepNumber (ts,&nsteps);
1690: TSGetConvergedReason (ts,&reason);
1691: PetscPrintf (PETSC_COMM_WORLD ,"%s at time %g after %D steps\n" ,TSConvergedReasons[reason],(double)ftime,nsteps);
1692: TSDestroy (&ts);
1694: PetscFunctionListDestroy (&PhysicsList);
1695: FunctionalLinkDestroy(&user->model->functionalRegistry);
1696: PetscFree (user->model->functionalMonitored);
1697: PetscFree (user->model->functionalCall);
1698: PetscFree (user->model->physics->data);
1699: PetscFree (user->model->physics);
1700: PetscFree (user->model);
1701: PetscFree (user);
1702: VecDestroy (&X);
1703: PetscFVDestroy (&fvm);
1704: DMDestroy (&dm);
1705: PetscFinalize ();
1706: return (0);
1707: }
1709: /* Godunov fluxs */
1710: PetscScalar cvmgp_(PetscScalar *a, PetscScalar *b, PetscScalar *test)1711: {
1712: /* System generated locals */
1713: PetscScalar ret_val;
1715: if (*test > 0.) {
1716: goto L10;
1717: }
1718: ret_val = *b;
1719: return ret_val;
1720: L10: 1721: ret_val = *a;
1722: return ret_val;
1723: } /* cvmgp_ */
1725: PetscScalar cvmgm_(PetscScalar *a, PetscScalar *b, PetscScalar *test)1726: {
1727: /* System generated locals */
1728: PetscScalar ret_val;
1730: if (*test < 0.) {
1731: goto L10;
1732: }
1733: ret_val = *b;
1734: return ret_val;
1735: L10: 1736: ret_val = *a;
1737: return ret_val;
1738: } /* cvmgm_ */
1742: PetscScalar *gaml, PetscScalar *gamr, PetscScalar *rl, PetscScalar *pl,1743: PetscScalar *uxl, PetscScalar *rr, PetscScalar *pr, 1744: PetscScalar *uxr, PetscScalar *rstarl, PetscScalar *rstarr, PetscScalar * 1745: pstar, PetscScalar *ustar) 1746: {
1747: /* Initialized data */
1749: static PetscScalar smallp = 1e-8;
1751: /* System generated locals */
1752: int i__1;
1753: PetscScalar d__1, d__2;
1755: /* Local variables */
1756: static int i0;
1757: static PetscScalar cl, cr, wl, zl, wr, zr, pst, durl, skpr1, skpr2;
1758: static int iwave;
1759: static PetscScalar csqrl, csqrr, gascl1, gascl2, gascl3, gascl4, gascr1, gascr2, gascr3, gascr4, cstarl, dpstar, cstarr;
1760: static int iterno;
1761: static PetscScalar ustarl, ustarr, rarepr1, rarepr2;
1765: gascl1 = *gaml - 1.;
1766: gascl2 = (*gaml + 1.) * .5;
1767: gascl3 = gascl2 / *gaml;
1768: gascl4 = 1.f / (*gaml - 1.);
1770: gascr1 = *gamr - 1.;
1771: gascr2 = (*gamr + 1.) * .5;
1772: gascr3 = gascr2 / *gamr;
1773: gascr4 = 1. / (*gamr - 1.);
1774: iterno = 10;
1775: /* find pstar: */
1776: cl = sqrt(*gaml * *pl / *rl);
1777: cr = sqrt(*gamr * *pr / *rr);
1778: wl = *rl * cl;
1779: wr = *rr * cr;
1780: csqrl = wl * wl;
1781: csqrr = wr * wr;
1782: *pstar = (wl * *pr + wr * *pl) / (wl + wr);
1783: *pstar = PetscMax (*pstar,smallp);
1784: pst = *pl / *pr;
1785: skpr1 = cr * (pst - 1.) * sqrt(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1786: d__1 = (*gamr - 1.) / (*gamr * 2.);
1787: rarepr2 = gascr4 * 2. * cr * (1. - PetscPowScalar(pst, d__1));
1788: pst = *pr / *pl;
1789: skpr2 = cl * (pst - 1.) * sqrt(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1790: d__1 = (*gaml - 1.) / (*gaml * 2.);
1791: rarepr1 = gascl4 * 2. * cl * (1. - PetscPowScalar(pst, d__1));
1792: durl = *uxr - *uxl;
1793: if (*pr < *pl) {
1794: if (durl >= rarepr1) {
1795: iwave = 100;
1796: } else if (durl <= -skpr1) {
1797: iwave = 300;
1798: } else {
1799: iwave = 400;
1800: }
1801: } else {
1802: if (durl >= rarepr2) {
1803: iwave = 100;
1804: } else if (durl <= -skpr2) {
1805: iwave = 300;
1806: } else {
1807: iwave = 200;
1808: }
1809: }
1810: if (iwave == 100) {
1811: /* 1-wave: rarefaction wave, 3-wave: rarefaction wave */
1812: /* case (100) */
1813: i__1 = iterno;
1814: for (i0 = 1; i0 <= i__1; ++i0) {
1815: d__1 = *pstar / *pl;
1816: d__2 = 1. / *gaml;
1817: *rstarl = *rl * PetscPowScalar(d__1, d__2);
1818: cstarl = sqrt(*gaml * *pstar / *rstarl);
1819: ustarl = *uxl - gascl4 * 2. * (cstarl - cl);
1820: zl = *rstarl * cstarl;
1821: d__1 = *pstar / *pr;
1822: d__2 = 1. / *gamr;
1823: *rstarr = *rr * PetscPowScalar(d__1, d__2);
1824: cstarr = sqrt(*gamr * *pstar / *rstarr);
1825: ustarr = *uxr + gascr4 * 2. * (cstarr - cr);
1826: zr = *rstarr * cstarr;
1827: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
1828: *pstar -= dpstar;
1829: *pstar = PetscMax (*pstar,smallp);
1830: if (PetscAbsScalar(dpstar) / *pstar <= 1e-8) {
1831: //break ;
1832: }
1833: }
1834: /* 1-wave: shock wave, 3-wave: rarefaction wave */
1835: } else if (iwave == 200) {
1836: /* case (200) */
1837: i__1 = iterno;
1838: for (i0 = 1; i0 <= i__1; ++i0) {
1839: pst = *pstar / *pl;
1840: ustarl = *uxl - (pst - 1.) * cl * sqrt(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1841: zl = *pl / cl * sqrt(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
1842: d__1 = *pstar / *pr;
1843: d__2 = 1. / *gamr;
1844: *rstarr = *rr * PetscPowScalar(d__1, d__2);
1845: cstarr = sqrt(*gamr * *pstar / *rstarr);
1846: zr = *rstarr * cstarr;
1847: ustarr = *uxr + gascr4 * 2. * (cstarr - cr);
1848: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
1849: *pstar -= dpstar;
1850: *pstar = PetscMax (*pstar,smallp);
1851: if (PetscAbsScalar(dpstar) / *pstar <= 1e-8) {
1852: //break ;
1853: }
1854: }
1855: /* 1-wave: shock wave, 3-wave: shock */
1856: } else if (iwave == 300) {
1857: /* case (300) */
1858: i__1 = iterno;
1859: for (i0 = 1; i0 <= i__1; ++i0) {
1860: pst = *pstar / *pl;
1861: ustarl = *uxl - (pst - 1.) * cl * sqrt(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1862: zl = *pl / cl * sqrt(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
1863: pst = *pstar / *pr;
1864: ustarr = *uxr + (pst - 1.) * cr * sqrt(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1865: zr = *pr / cr * sqrt(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
1866: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
1867: *pstar -= dpstar;
1868: *pstar = PetscMax (*pstar,smallp);
1869: if (PetscAbsScalar(dpstar) / *pstar <= 1e-8) {
1870: //break ;
1871: }
1872: }
1873: /* 1-wave: rarefaction wave, 3-wave: shock */
1874: } else if (iwave == 400) {
1875: /* case (400) */
1876: i__1 = iterno;
1877: for (i0 = 1; i0 <= i__1; ++i0) {
1878: d__1 = *pstar / *pl;
1879: d__2 = 1. / *gaml;
1880: *rstarl = *rl * PetscPowScalar(d__1, d__2);
1881: cstarl = sqrt(*gaml * *pstar / *rstarl);
1882: ustarl = *uxl - gascl4 * 2. * (cstarl - cl);
1883: zl = *rstarl * cstarl;
1884: pst = *pstar / *pr;
1885: ustarr = *uxr + (pst - 1.) * cr * sqrt(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1886: zr = *pr / cr * sqrt(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
1887: dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
1888: *pstar -= dpstar;
1889: *pstar = PetscMax (*pstar,smallp);
1890: if (PetscAbsScalar(dpstar) / *pstar <= 1e-8) {
1891: //break ;
1892: }
1893: }
1894: }
1896: *ustar = (zl * ustarr + zr * ustarl) / (zl + zr);
1897: if (*pstar > *pl) {
1898: pst = *pstar / *pl;
1899: *rstarl = ((*gaml + 1.) * pst + *gaml - 1.) / ((*gaml - 1.) * pst + *
1900: gaml + 1.) * *rl;
1901: }
1902: if (*pstar > *pr) {
1903: pst = *pstar / *pr;
1904: *rstarr = ((*gamr + 1.) * pst + *gamr - 1.) / ((*gamr - 1.) * pst + *
1905: gamr + 1.) * *rr;
1906: }
1907: return iwave;
1908: }
1910: PetscScalar sign(PetscScalar x)1911: {
1912: if (x > 0) return 1.0;
1913: if (x < 0) return -1.0;
1914: return 0.0;
1915: }
1916: /* Riemann Solver */
1919: /* -------------------------------------------------------------------- */
1920: PetscScalar *xcen, PetscScalar *xp,1921: PetscScalar *dtt, PetscScalar *rl, PetscScalar *uxl, PetscScalar *pl, 1922: PetscScalar *utl, PetscScalar *ubl, PetscScalar *gaml, PetscScalar *rho1l, 1923: PetscScalar *rr, PetscScalar *uxr, PetscScalar *pr, PetscScalar *utr, 1924: PetscScalar *ubr, PetscScalar *gamr, PetscScalar *rho1r, PetscScalar *rx, 1925: PetscScalar *uxm, PetscScalar *px, PetscScalar *utx, PetscScalar *ubx, 1926: PetscScalar *gam, PetscScalar *rho1) 1927: {
1928: /* System generated locals */
1929: PetscScalar d__1, d__2;
1931: /* Local variables */
1932: static PetscScalar s, c0, p0, r0, u0, w0, x0, x2, ri, cx, sgn0, wsp0, gasc1, gasc2, gasc3, gasc4;
1933: static PetscScalar cstar, pstar, rstar, ustar, xstar, wspst, ushock, streng, rstarl, rstarr, rstars;
1935: if (*rl == *rr && *pr == *pl && *uxl == *uxr && *gaml == *gamr) {
1936: *rx = *rl;
1937: *px = *pl;
1938: *uxm = *uxl;
1939: *gam = *gaml;
1940: x2 = *xcen + *uxm * *dtt;
1942: if (*xp >= x2) {
1943: *utx = *utr;
1944: *ubx = *ubr;
1945: *rho1 = *rho1r;
1946: } else {
1947: *utx = *utl;
1948: *ubx = *ubl;
1949: *rho1 = *rho1l;
1950: }
1951: return 0;
1952: }
1953: int iwave = riem1mdt(gaml, gamr, rl, pl, uxl, rr, pr, uxr, &rstarl, &rstarr, &pstar, &ustar);
1955: x2 = *xcen + ustar * *dtt;
1956: d__1 = *xp - x2;
1957: sgn0 = sign(d__1);
1958: /* x is in 3-wave if sgn0 = 1 */
1959: /* x is in 1-wave if sgn0 = -1 */
1960: r0 = cvmgm_(rl, rr, &sgn0);
1961: p0 = cvmgm_(pl, pr, &sgn0);
1962: u0 = cvmgm_(uxl, uxr, &sgn0);
1963: *gam = cvmgm_(gaml, gamr, &sgn0);
1964: gasc1 = *gam - 1.;
1965: gasc2 = (*gam + 1.) * .5;
1966: gasc3 = gasc2 / *gam;
1967: gasc4 = 1. / (*gam - 1.);
1968: c0 = sqrt(*gam * p0 / r0);
1969: streng = pstar - p0;
1970: w0 = *gam * r0 * p0 * (gasc3 * streng / p0 + 1.);
1971: rstars = r0 / (1. - r0 * streng / w0);
1972: d__1 = p0 / pstar;
1973: d__2 = -1. / *gam;
1974: rstarr = r0 * PetscPowScalar(d__1, d__2);
1975: rstar = cvmgm_(&rstarr, &rstars, &streng);
1976: w0 = sqrt(w0);
1977: cstar = sqrt(*gam * pstar / rstar);
1978: wsp0 = u0 + sgn0 * c0;
1979: wspst = ustar + sgn0 * cstar;
1980: ushock = ustar + sgn0 * w0 / rstar;
1981: wspst = cvmgp_(&ushock, &wspst, &streng);
1982: wsp0 = cvmgp_(&ushock, &wsp0, &streng);
1983: x0 = *xcen + wsp0 * *dtt;
1984: xstar = *xcen + wspst * *dtt;
1985: /* using gas formula to evaluate rarefaction wave */
1986: /* ri : reiman invariant */
1987: ri = u0 - sgn0 * 2. * gasc4 * c0;
1988: cx = sgn0 * .5 * gasc1 / gasc2 * ((*xp - *xcen) / *dtt - ri);
1989: *uxm = ri + sgn0 * 2. * gasc4 * cx;
1990: s = p0 / PetscPowScalar(r0, *gam);
1991: d__1 = cx * cx / (*gam * s);
1992: *rx = PetscPowScalar(d__1, gasc4);
1993: *px = cx * cx * *rx / *gam;
1994: d__1 = sgn0 * (x0 - *xp);
1995: *rx = cvmgp_(rx, &r0, &d__1);
1996: d__1 = sgn0 * (x0 - *xp);
1997: *px = cvmgp_(px, &p0, &d__1);
1998: d__1 = sgn0 * (x0 - *xp);
1999: *uxm = cvmgp_(uxm, &u0, &d__1);
2000: d__1 = sgn0 * (xstar - *xp);
2001: *rx = cvmgm_(rx, &rstar, &d__1);
2002: d__1 = sgn0 * (xstar - *xp);
2003: *px = cvmgm_(px, &pstar, &d__1);
2004: d__1 = sgn0 * (xstar - *xp);
2005: *uxm = cvmgm_(uxm, &ustar, &d__1);
2006: if (*xp >= x2) {
2007: *utx = *utr;
2008: *ubx = *ubr;
2009: *rho1 = *rho1r;
2010: } else {
2011: *utx = *utl;
2012: *ubx = *ubl;
2013: *rho1 = *rho1l;
2014: }
2015: return iwave;
2016: }
2019: PetscScalar *ul, const PetscScalar *ur,2020: PetscScalar *flux, const PetscScalar *nn, const int *ndim, 2021: const PetscScalar *gamma) 2022: {
2023: /* System generated locals */
2024: int i__1,iwave;
2025: PetscScalar d__1, d__2, d__3;
2027: /* Local variables */
2028: static int k;
2029: static PetscScalar bn[3], fn, ft, tg[3], pl, rl, pm, pr, rr, xp, ubl, ubm,
2030: ubr, dtt, unm, tmp, utl, utm, uxl, utr, uxr, gaml, gamm, gamr,
2031: xcen, rhom, rho1l, rho1m, rho1r;
2032: /* Parameter adjustments */
2033: --nn;
2034: --flux;
2035: --ur;
2036: --ul;
2038: /* Function Body */
2039: xcen = 0.;
2040: xp = 0.;
2041: i__1 = *ndim;
2042: for (k = 1; k <= i__1; ++k) {
2043: tg[k - 1] = 0.;
2044: bn[k - 1] = 0.;
2045: }
2046: dtt = 1.;
2047: if (*ndim == 3) {
2048: if (nn[1] == 0. && nn[2] == 0.) {
2049: tg[0] = 1.;
2050: } else {
2051: tg[0] = -nn[2];
2052: tg[1] = nn[1];
2053: }
2054: /* tmp=dsqrt(tg(1)**2+tg(2)**2) */
2055: /* tg=tg/tmp */
2056: bn[0] = -nn[3] * tg[1];
2057: bn[1] = nn[3] * tg[0];
2058: bn[2] = nn[1] * tg[1] - nn[2] * tg[0];
2059: /* Computing 2nd power */
2060: d__1 = bn[0];
2061: /* Computing 2nd power */
2062: d__2 = bn[1];
2063: /* Computing 2nd power */
2064: d__3 = bn[2];
2065: tmp = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
2066: i__1 = *ndim;
2067: for (k = 1; k <= i__1; ++k) {
2068: bn[k - 1] /= tmp;
2069: }
2070: } else if (*ndim == 2) {
2071: tg[0] = -nn[2];
2072: tg[1] = nn[1];
2073: /* tmp=dsqrt(tg(1)**2+tg(2)**2) */
2074: /* tg=tg/tmp */
2075: bn[0] = 0.;
2076: bn[1] = 0.;
2077: bn[2] = 1.;
2078: }
2079: rl = ul[1];
2080: rr = ur[1];
2081: uxl = 0.;
2082: uxr = 0.;
2083: utl = 0.;
2084: utr = 0.;
2085: ubl = 0.;
2086: ubr = 0.;
2087: i__1 = *ndim;
2088: for (k = 1; k <= i__1; ++k) {
2089: uxl += ul[k + 1] * nn[k];
2090: uxr += ur[k + 1] * nn[k];
2091: utl += ul[k + 1] * tg[k - 1];
2092: utr += ur[k + 1] * tg[k - 1];
2093: ubl += ul[k + 1] * bn[k - 1];
2094: ubr += ur[k + 1] * bn[k - 1];
2095: }
2096: uxl /= rl;
2097: uxr /= rr;
2098: utl /= rl;
2099: utr /= rr;
2100: ubl /= rl;
2101: ubr /= rr;
2103: gaml = *gamma;
2104: gamr = *gamma;
2105: /* Computing 2nd power */
2106: d__1 = uxl;
2107: /* Computing 2nd power */
2108: d__2 = utl;
2109: /* Computing 2nd power */
2110: d__3 = ubl;
2111: pl = (*gamma - 1.) * (ul[*ndim + 2] - rl * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
2112: /* Computing 2nd power */
2113: d__1 = uxr;
2114: /* Computing 2nd power */
2115: d__2 = utr;
2116: /* Computing 2nd power */
2117: d__3 = ubr;
2118: pr = (*gamma - 1.) * (ur[*ndim + 2] - rr * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
2119: rho1l = rl;
2120: rho1r = rr;
2122: iwave = riemannsolver(&xcen, &xp, &dtt, &rl, &uxl, &pl, &utl, &ubl, &gaml, &
2123: rho1l, &rr, &uxr, &pr, &utr, &ubr, &gamr, &rho1r, &rhom, &unm, &
2124: pm, &utm, &ubm, &gamm, &rho1m);
2126: flux[1] = rhom * unm;
2127: fn = rhom * unm * unm + pm;
2128: ft = rhom * unm * utm;
2129: /* flux(2)=fn*nn(1)+ft*nn(2) */
2130: /* flux(3)=fn*tg(1)+ft*tg(2) */
2131: flux[2] = fn * nn[1] + ft * tg[0];
2132: flux[3] = fn * nn[2] + ft * tg[1];
2133: /* flux(2)=rhom*unm*(unm)+pm */
2134: /* flux(3)=rhom*(unm)*utm */
2135: if (*ndim == 3) {
2136: flux[4] = rhom * unm * ubm;
2137: }
2138: flux[*ndim + 2] = (rhom * .5 * (unm * unm + utm * utm + ubm * ubm) + gamm / (gamm - 1.) * pm) * unm;
2139: return iwave;
2140: } /* godunovflux_ */
2142: /* Subroutine to set up the initial conditions for the */
2143: /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
2144: /* ----------------------------------------------------------------------- */
2145: int projecteqstate(PetscScalar wc[], const PetscScalar ueq[], const PetscScalar lv[][3])
2146: {
2147: int j,k;
2148: /* Wc=matmul(lv,Ueq) 3 vars */
2149: for (k = 0; k < 3; ++k) {
2150: wc[k] = 0.;
2151: for (j = 0; j < 3; ++j) {
2152: wc[k] += lv[k][j]*ueq[j];
2153: }
2154: }
2155: return 0;
2156: }
2157: /* ----------------------------------------------------------------------- */
2158: int projecttoprim(PetscScalar v[], const PetscScalar wc[], PetscScalar rv[][3])
2159: {
2160: int k,j;
2161: /* V=matmul(rv,WC) 3 vars */
2162: for (k = 0; k < 3; ++k) {
2163: v[k] = 0.;
2164: for (j = 0; j < 3; ++j) {
2165: v[k] += rv[k][j]*wc[j];
2166: }
2167: }
2168: return 0;
2169: }
2170: /* ---------------------------------------------------------------------- */
2171: int eigenvectors(PetscScalar rv[][3], PetscScalar lv[][3], const PetscScalar ueq[], PetscScalar gamma)
2172: {
2173: int j,k;
2174: PetscScalar rho,csnd,p0,u;
2176: for (k = 0; k < 3; ++k) for (j = 0; j < 3; ++j) { lv[k][j] = 0.; rv[k][j] = 0.; }
2177: rho = ueq[0];
2178: u = ueq[1];
2179: p0 = ueq[2];
2180: csnd = PetscSqrtScalar(gamma * p0 / rho);
2181: lv[0][1] = rho * .5;
2182: lv[0][2] = -.5 / csnd;
2183: lv[1][0] = csnd;
2184: lv[1][2] = -1. / csnd;
2185: lv[2][1] = rho * .5;
2186: lv[2][2] = .5 / csnd;
2187: rv[0][0] = -1. / csnd;
2188: rv[1][0] = 1. / rho;
2189: rv[2][0] = -csnd;
2190: rv[0][1] = 1. / csnd;
2191: rv[0][2] = 1. / csnd;
2192: rv[1][2] = 1. / rho;
2193: rv[2][2] = csnd;
2194: return 0;
2195: }
2198: PetscScalar gamma, const PetscReal coord[], const PetscReal Lx)2199: {
2200: PetscScalar p0,u0,wcp[3],wc[3];
2201: PetscScalar lv[3][3];
2202: PetscScalar vp[3];
2203: PetscScalar rv[3][3];
2204: PetscScalar eps, ueq[3], rho0, twopi;
2206: /* Function Body */
2207: twopi = 2.*M_PI;
2208: eps = 1e-4; /* perturbation */
2209: rho0 = 1e3; /* density of water */
2210: p0 = 101325.; /* init pressure of 1 atm (?) */
2211: u0 = 0.;
2212: ueq[0] = rho0;
2213: ueq[1] = u0;
2214: ueq[2] = p0;
2215: /* Project initial state to characteristic variables */
2216: eigenvectors(rv, lv, ueq, gamma);
2217: projecteqstate(wc, ueq, lv);
2218: wcp[0] = wc[0];
2219: wcp[1] = wc[1];
2220: wcp[2] = wc[2] + eps * cos(coord[0] * 2. * twopi / Lx);
2221: projecttoprim(vp, wcp, rv);
2222: ux->r = vp[0]; /* density */
2223: ux->ru[0] = vp[0] * vp[1]; /* x momentum */
2224: ux->ru[1] = 0.;
2225: #if defined DIM > 2
2226: if (dim>2) ux->ru[2] = 0.;
2227: #endif
2228: /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
2229: ux->E = vp[2]/(gamma - 1.) + 0.5*vp[0]*vp[1]*vp[1];
2230: return 0;
2231: }