Actual source code: ex11.c

petsc-dev 2014-04-19
Report Typos and Errors
  1: static char help[] = "Second Order TVD Finite Volume Example.\n";

We use a second order TVD finite volume method to evolve a system of PDEs. Our simple upwinded residual evaluation loops
over all mesh faces and uses a Riemann solver to produce the flux given the face geometry and cell values,
\begin{equation}
f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R)
\end{equation}
and then update the cell values given the cell volume.
\begin{eqnarray}
f^L_i &-=& \frac{f_i}{vol^L} \\
f^R_i &+=& \frac{f_i}{vol^R}
\end{eqnarray}

As an example, we can consider the shallow water wave equation,
\begin{eqnarray}
h_t + \nabla\cdot \left( uh \right) &=& 0 \\
(uh)_t + \nabla\cdot \left( u\otimes uh + \frac{g h^2}{2} I \right) &=& 0
\end{eqnarray}
where $h$ is wave height, $u$ is wave velocity, and $g$ is the acceleration due to gravity.

A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
\begin{eqnarray}
f^{L,R}_h &=& uh^{L,R} \cdot \hat n \\
f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
c^{L,R} &=& \sqrt{g h^{L,R}} \\
s &=& \max\left( \left|\frac{uh^L \cdot \hat n}{h^L}\right| + c^L, \left|\frac{uh^R \cdot \hat n}{h^R}\right| + c^R \right) \\
f_i &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
\end{eqnarray}
where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.

The more sophisticated residual evaluation in RHSFunctionLocal_LS() uses a least-squares fit to a quadratic polynomial
over a neighborhood of the given element.

The mesh is read in from an ExodusII file, usually generated by Cubit.
 37: #include <petscts.h>
 38: #include <petscdmplex.h>
 39: #include <petscsf.h>
 40: #include <petscblaslapack.h>

 42: #define DIM 2                   /* Geometric dimension */
 43: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))

 45: static PetscFunctionList PhysicsList;
 46: static PetscFunctionList LimitList;

 48: /* Represents continuum physical equations. */
 49: typedef struct _n_Physics *Physics;

 51: /* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
 52:  * discretization-independent, but its members depend on the scenario being solved. */
 53: typedef struct _n_Model *Model;

 55: /* 'User' implements a discretization of a continuous model. */
 56: typedef struct _n_User *User;

 58: typedef PetscErrorCode (*RiemannFunction)(Physics,const PetscReal*,const PetscReal*,const PetscScalar*,const PetscScalar*,PetscScalar*);
 59: typedef PetscErrorCode (*SolutionFunction)(Model,PetscReal,const PetscReal*,PetscScalar*,void*);
 60: typedef PetscErrorCode (*BoundaryFunction)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*);
 61: typedef PetscErrorCode (*FunctionalFunction)(Model,PetscReal,const PetscReal*,const PetscScalar*,PetscReal*,void*);
 62: typedef PetscErrorCode (*SetupFields)(Physics,PetscSection);
 63: static PetscErrorCode ModelBoundaryRegister(Model,const char*,BoundaryFunction,void*,PetscInt,const PetscInt*);
 64: static PetscErrorCode ModelSolutionSetDefault(Model,SolutionFunction,void*);
 65: static PetscErrorCode ModelFunctionalRegister(Model,const char*,PetscInt*,FunctionalFunction,void*);
 66: static PetscErrorCode OutputVTK(DM,const char*,PetscViewer*);

 68: struct FieldDescription {
 69:   const char *name;
 70:   PetscInt dof;
 71: };

 73: typedef struct _n_BoundaryLink *BoundaryLink;
 74: struct _n_BoundaryLink {
 75:   char             *name;
 76:   BoundaryFunction func;
 77:   void             *ctx;
 78:   PetscInt         numids;
 79:   PetscInt         *ids;
 80:   BoundaryLink     next;
 81: };

 83: typedef struct _n_FunctionalLink *FunctionalLink;
 84: struct _n_FunctionalLink {
 85:   char               *name;
 86:   FunctionalFunction func;
 87:   void               *ctx;
 88:   PetscInt           offset;
 89:   FunctionalLink     next;
 90: };

 92: struct _n_Physics {
 93:   RiemannFunction riemann;
 94:   PetscInt        dof;          /* number of degrees of freedom per cell */
 95:   PetscReal       maxspeed;     /* kludge to pick initial time step, need to add monitoring and step control */
 96:   void            *data;
 97:   PetscInt        nfields;
 98:   const struct FieldDescription *field_desc;
 99: };

101: struct _n_Model {
102:   MPI_Comm         comm;        /* Does not do collective communicaton, but some error conditions can be collective */
103:   Physics          physics;
104:   BoundaryLink     boundary;
105:   FunctionalLink   functionalRegistry;
106:   PetscInt         maxComputed;
107:   PetscInt         numMonitored;
108:   FunctionalLink   *functionalMonitored;
109:   PetscInt         numCall;
110:   FunctionalLink   *functionalCall;
111:   SolutionFunction solution;
112:   void             *solutionctx;
113:   PetscReal        maxspeed;    /* estimate of global maximum speed (for CFL calculation) */
114: };

116: struct _n_User {
117:   PetscReal      (*Limit)(PetscReal);
118:   PetscBool      reconstruct;
119:   PetscInt       numGhostCells, numSplitFaces;
120:   PetscInt       cEndInterior; /* First boundary ghost cell */
121:   Vec            cellgeom, facegeom;
122:   DM             dmGrad;
123:   PetscReal      minradius;
124:   PetscInt       vtkInterval;   /* For monitor */
125:   Model          model;
126:   struct {
127:     PetscScalar *flux;
128:     PetscScalar *state0;
129:     PetscScalar *state1;
130:   } work;
131: };

133: typedef struct {
134:   PetscScalar normal[DIM];              /* Area-scaled normals */
135:   PetscScalar centroid[DIM];            /* Location of centroid (quadrature point) */
136:   PetscScalar grad[2][DIM];             /* Face contribution to gradient in left and right cell */
137: } FaceGeom;

139: typedef struct {
140:   PetscScalar centroid[DIM];
141:   PetscScalar volume;
142: } CellGeom;


145: PETSC_STATIC_INLINE PetscScalar DotDIM(const PetscScalar *x,const PetscScalar *y)
146: {
147:   PetscInt    i;
148:   PetscScalar prod=0.0;

150:   for (i=0; i<DIM; i++) prod += x[i]*y[i];
151:   return prod;
152: }
153: PETSC_STATIC_INLINE PetscReal NormDIM(const PetscScalar *x) { return PetscSqrtReal(PetscAbsScalar(DotDIM(x,x))); }
154: PETSC_STATIC_INLINE void axDIM(const PetscScalar a,PetscScalar *x)
155: {
156:   PetscInt i;
157:   for (i=0; i<DIM; i++) x[i] *= a;
158: }
159: PETSC_STATIC_INLINE void waxDIM(const PetscScalar a,const PetscScalar *x, PetscScalar *w)
160: {
161:   PetscInt i;
162:   for (i=0; i<DIM; i++) w[i] = x[i]*a;
163: }
164: PETSC_STATIC_INLINE void NormalSplitDIM(const PetscReal *n,const PetscScalar *x,PetscScalar *xn,PetscScalar *xt)
165: {                               /* Split x into normal and tangential components */
166:   PetscInt    i;
167:   PetscScalar c;
168:   c = DotDIM(x,n)/DotDIM(n,n);
169:   for (i=0; i<DIM; i++) {
170:     xn[i] = c*n[i];
171:     xt[i] = x[i]-xn[i];
172:   }
173: }

175: PETSC_STATIC_INLINE PetscScalar Dot2(const PetscScalar *x,const PetscScalar *y) { return x[0]*y[0] + x[1]*y[1];}
176: PETSC_STATIC_INLINE PetscReal Norm2(const PetscScalar *x) { return PetscSqrtReal(PetscAbsScalar(Dot2(x,x)));}
177: PETSC_STATIC_INLINE void Normalize2(PetscScalar *x) { PetscReal a = 1./Norm2(x); x[0] *= a; x[1] *= a; }
178: PETSC_STATIC_INLINE void Waxpy2(PetscScalar a,const PetscScalar *x,const PetscScalar *y,PetscScalar *w) { w[0] = a*x[0] + y[0]; w[1] = a*x[1] + y[1]; }
179: PETSC_STATIC_INLINE void Scale2(PetscScalar a,const PetscScalar *x,PetscScalar *y) { y[0] = a*x[0]; y[1] = a*x[1]; }

181: PETSC_STATIC_INLINE void WaxpyD(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];}
182: PETSC_STATIC_INLINE 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;}
183: PETSC_STATIC_INLINE PetscReal NormD(PetscInt dim, const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(DotD(dim,x,x)));}

185: PETSC_STATIC_INLINE void NormalSplit(const PetscReal *n,const PetscScalar *x,PetscScalar *xn,PetscScalar *xt)
186: {                               /* Split x into normal and tangential components */
187:   Scale2(Dot2(x,n)/Dot2(n,n),n,xn);
188:   Waxpy2(-1,xn,x,xt);
189: }

191: /* Limiters given in symmetric form following Berger, Aftosmis, and Murman 2005
192:  *
193:  * The classical flux-limited formulation is psi(r) where
194:  *
195:  * r = (u[0] - u[-1]) / (u[1] - u[0])
196:  *
197:  * The second order TVD region is bounded by
198:  *
199:  * psi_minmod(r) = min(r,1)      and        psi_superbee(r) = min(2, 2r, max(1,r))
200:  *
201:  * where all limiters are implicitly clipped to be non-negative. A more convenient slope-limited form is psi(r) =
202:  * phi(r)(r+1)/2 in which the reconstructed interface values are
203:  *
204:  * u(v) = u[0] + phi(r) (grad u)[0] v
205:  *
206:  * where v is the vector from centroid to quadrature point. In these variables, the usual limiters become
207:  *
208:  * 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))
209:  *
210:  * For a nicer symmetric formulation, rewrite in terms of
211:  *
212:  * f = (u[0] - u[-1]) / (u[1] - u[-1])
213:  *
214:  * where r(f) = f/(1-f). Not that r(1-f) = (1-f)/f = 1/r(f) so the symmetry condition
215:  *
216:  * phi(r) = phi(1/r)
217:  *
218:  * becomes
219:  *
220:  * w(f) = w(1-f).
221:  *
222:  * The limiters below implement this final form w(f). The reference methods are
223:  *
224:  * w_minmod(f) = 2 min(f,(1-f))             w_superbee(r) = 4 min((1-f), f)
225:  * */
226: static PetscReal Limit_Zero(PetscReal f) { return 0; }
227: static PetscReal Limit_None(PetscReal f) { return 1; }
228: static PetscReal Limit_Minmod(PetscReal f) { return PetscMax(0,PetscMin(f,(1-f))*2); }
229: static PetscReal Limit_VanLeer(PetscReal f) { return PetscMax(0,4*f*(1-f)); }
230: static PetscReal Limit_VanAlbada(PetscReal f) { return PetscMax(0,2*f*(1-f) / (PetscSqr(f) + PetscSqr(1-f))); }
231: static PetscReal Limit_Sin(PetscReal f)
232: {
233:   PetscReal fclip = PetscMax(0,PetscMin(f,1));
234:   return PetscSinReal(PETSC_PI*fclip);
235: }
236: static PetscReal Limit_Superbee(PetscReal f) { return 2*Limit_Minmod(f); }
237: static PetscReal Limit_MC(PetscReal f) { return PetscMin(1,Limit_Superbee(f)); } /* aka Barth-Jespersen */

239: /******************* Advect ********************/
240: typedef enum {ADVECT_SOL_TILTED,ADVECT_SOL_BUMP} AdvectSolType;
241: static const char *const AdvectSolTypes[] = {"TILTED","BUMP","AdvectSolType","ADVECT_SOL_",0};
242: typedef enum {ADVECT_SOL_BUMP_CONE,ADVECT_SOL_BUMP_COS} AdvectSolBumpType;
243: static const char *const AdvectSolBumpTypes[] = {"CONE","COS","AdvectSolBumpType","ADVECT_SOL_BUMP_",0};

245: typedef struct {
246:   PetscReal wind[DIM];
247: } Physics_Advect_Tilted;
248: typedef struct {
249:   PetscReal         center[DIM];
250:   PetscReal         radius;
251:   AdvectSolBumpType type;
252: } Physics_Advect_Bump;

254: typedef struct {
255:   PetscReal     inflowState;
256:   AdvectSolType soltype;
257:   union {
258:     Physics_Advect_Tilted tilted;
259:     Physics_Advect_Bump   bump;
260:   } sol;
261:   struct {
262:     PetscInt Error;
263:   } functional;
264: } Physics_Advect;

266: static const struct FieldDescription PhysicsFields_Advect[] = {{"U",1},{NULL,0}};

270: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
271: {
272:   Physics        phys    = (Physics)ctx;
273:   Physics_Advect *advect = (Physics_Advect*)phys->data;

276:   xG[0] = advect->inflowState;
277:   return(0);
278: }

282: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
283: {
285:   xG[0] = xI[0];
286:   return(0);
287: }

291: static PetscErrorCode PhysicsRiemann_Advect(Physics phys, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux)
292: {
293:   Physics_Advect *advect = (Physics_Advect*)phys->data;
294:   PetscReal      wind[DIM],wn;

297:   switch (advect->soltype) {
298:   case ADVECT_SOL_TILTED: {
299:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
300:     wind[0] = tilted->wind[0];
301:     wind[1] = tilted->wind[1];
302:   } break;
303:   case ADVECT_SOL_BUMP:
304:     wind[0] = -qp[1];
305:     wind[1] = qp[0];
306:     break;
307:   default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]);
308:   }
309:   wn      = Dot2(wind, n);
310:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
311:   return(0);
312: }

316: static PetscErrorCode PhysicsSolution_Advect(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
317: {
318:   Physics        phys    = (Physics)ctx;
319:   Physics_Advect *advect = (Physics_Advect*)phys->data;

322:   switch (advect->soltype) {
323:   case ADVECT_SOL_TILTED: {
324:     PetscReal             x0[DIM];
325:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
326:     Waxpy2(-time,tilted->wind,x,x0);
327:     if (x0[1] > 0) u[0] = 1.*x[0] + 3.*x[1];
328:     else u[0] = advect->inflowState;
329:   } break;
330:   case ADVECT_SOL_BUMP: {
331:     Physics_Advect_Bump *bump = &advect->sol.bump;
332:     PetscReal           x0[DIM],v[DIM],r,cost,sint;
333:     cost  = PetscCosReal(time);
334:     sint  = PetscSinReal(time);
335:     x0[0] = cost*x[0] + sint*x[1];
336:     x0[1] = -sint*x[0] + cost*x[1];
337:     Waxpy2(-1,bump->center,x0,v);
338:     r = Norm2(v);
339:     switch (bump->type) {
340:     case ADVECT_SOL_BUMP_CONE:
341:       u[0] = PetscMax(1 - r/bump->radius,0);
342:       break;
343:     case ADVECT_SOL_BUMP_COS:
344:       u[0] = 0.5 + 0.5*PetscCosReal(PetscMin(r/bump->radius,1)*PETSC_PI);
345:       break;
346:     }
347:   } break;
348:   default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown solution type");
349:   }
350:   return(0);
351: }

355: static PetscErrorCode PhysicsFunctional_Advect(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *y,PetscReal *f,void *ctx)
356: {
357:   Physics        phys    = (Physics)ctx;
358:   Physics_Advect *advect = (Physics_Advect*)phys->data;
359:   PetscScalar    yexact[1];

363:   PhysicsSolution_Advect(mod,time,x,yexact,phys);
364:   f[advect->functional.Error] = PetscAbsScalar(y[0]-yexact[0]);
365:   return(0);
366: }

370: static PetscErrorCode PhysicsCreate_Advect(DM dm, Model mod,Physics phys)
371: {
372:   Physics_Advect *advect;

376:   phys->field_desc = PhysicsFields_Advect;
377:   phys->riemann = PhysicsRiemann_Advect;
378:   PetscNew(&advect);
379:   phys->data = advect;
380:   PetscOptionsHead("Advect options");
381:   {
382:     PetscInt two = 2,dof = 1;
383:     advect->soltype = ADVECT_SOL_TILTED;
384:     PetscOptionsEnum("-advect_sol_type","solution type","",AdvectSolTypes,(PetscEnum)advect->soltype,(PetscEnum*)&advect->soltype,NULL);
385:     switch (advect->soltype) {
386:     case ADVECT_SOL_TILTED: {
387:       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
388:       two = 2;
389:       tilted->wind[0] = 0.0;
390:       tilted->wind[1] = 1.0;
391:       PetscOptionsRealArray("-advect_tilted_wind","background wind vx,vy","",tilted->wind,&two,NULL);
392:       advect->inflowState = -2.0;
393:       PetscOptionsRealArray("-advect_tilted_inflow","Inflow state","",&advect->inflowState,&dof,NULL);
394:       phys->maxspeed = Norm2(tilted->wind);
395:     } break;
396:     case ADVECT_SOL_BUMP: {
397:       Physics_Advect_Bump *bump = &advect->sol.bump;
398:       two = 2;
399:       bump->center[0] = 2.;
400:       bump->center[1] = 0.;
401:       PetscOptionsRealArray("-advect_bump_center","location of center of bump x,y","",bump->center,&two,NULL);
402:       bump->radius = 0.9;
403:       PetscOptionsReal("-advect_bump_radius","radius of bump","",bump->radius,&bump->radius,NULL);
404:       bump->type = ADVECT_SOL_BUMP_CONE;
405:       PetscOptionsEnum("-advect_bump_type","type of bump","",AdvectSolBumpTypes,(PetscEnum)bump->type,(PetscEnum*)&bump->type,NULL);
406:       phys->maxspeed = 3.;       /* radius of mesh, kludge */
407:     } break;
408:     }
409:   }
410:   PetscOptionsTail();
411:   {
412:     const PetscInt inflowids[] = {100,200,300},outflowids[] = {101};
413:     /* Register "canned" boundary conditions and defaults for where to apply. */
414:     ModelBoundaryRegister(mod,"inflow",PhysicsBoundary_Advect_Inflow,phys,ALEN(inflowids),inflowids);
415:     ModelBoundaryRegister(mod,"outflow",PhysicsBoundary_Advect_Outflow,phys,ALEN(outflowids),outflowids);
416:     DMPlexAddBoundary(dm, PETSC_TRUE, "inflow",  0, (void (*)()) PhysicsBoundary_Advect_Inflow,  ALEN(inflowids),  inflowids,  phys);
417:     DMPlexAddBoundary(dm, PETSC_TRUE, "outflow", 0, (void (*)()) PhysicsBoundary_Advect_Outflow, ALEN(outflowids), outflowids, phys);
418:     /* Initial/transient solution with default boundary conditions */
419:     ModelSolutionSetDefault(mod,PhysicsSolution_Advect,phys);
420:     /* Register "canned" functionals */
421:     ModelFunctionalRegister(mod,"Error",&advect->functional.Error,PhysicsFunctional_Advect,phys);
422:   }
423:   return(0);
424: }

426: /******************* Shallow Water ********************/
427: typedef struct {
428:   PetscReal gravity;
429:   PetscReal boundaryHeight;
430:   struct {
431:     PetscInt Height;
432:     PetscInt Speed;
433:     PetscInt Energy;
434:   } functional;
435: } Physics_SW;
436: typedef struct {
437:   PetscScalar vals[0];
438:   PetscScalar h;
439:   PetscScalar uh[DIM];
440: } SWNode;

442: static const struct FieldDescription PhysicsFields_SW[] = {{"Height",1},{"Momentum",DIM},{NULL,0}};

446: /*
447:  * h_t + div(uh) = 0
448:  * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
449:  *
450:  * */
451: static PetscErrorCode SWFlux(Physics phys,const PetscReal *n,const SWNode *x,SWNode *f)
452: {
453:   Physics_SW  *sw = (Physics_SW*)phys->data;
454:   PetscScalar uhn,u[DIM];
455:   PetscInt    i;

458:   Scale2(1./x->h,x->uh,u);
459:   uhn  = Dot2(x->uh,n);
460:   f->h = uhn;
461:   for (i=0; i<DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
462:   return(0);
463: }

467: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
468: {
470:   xG[0] = xI[0];
471:   xG[1] = -xI[1];
472:   xG[2] = -xI[2];
473:   return(0);
474: }

478: static PetscErrorCode PhysicsRiemann_SW(Physics phys, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux)
479: {
480:   Physics_SW   *sw = (Physics_SW*)phys->data;
481:   PetscReal    cL,cR,speed,nn[DIM];
482:   const SWNode *uL = (const SWNode*)xL,*uR = (const SWNode*)xR;
483:   SWNode       fL,fR;
484:   PetscInt     i;

487:   if (uL->h < 0 || uR->h < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative");
488:   nn[0] = n[0];
489:   nn[1] = n[1];
490:   Normalize2(nn);
491:   SWFlux(phys,nn,uL,&fL);
492:   SWFlux(phys,nn,uR,&fR);
493:   cL    = PetscSqrtReal(sw->gravity*PetscRealPart(uL->h));
494:   cR    = PetscSqrtReal(sw->gravity*PetscRealPart(uR->h)); /* gravity wave speed */
495:   speed = PetscMax(PetscAbsScalar(Dot2(uL->uh,nn)/uL->h) + cL,PetscAbsScalar(Dot2(uR->uh,nn)/uR->h) + cR);
496:   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);
497:   return(0);
498: }

502: static PetscErrorCode PhysicsSolution_SW(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
503: {
504:   PetscReal dx[2],r,sigma;

507:   if (time != 0.0) SETERRQ1(mod->comm,PETSC_ERR_SUP,"No solution known for time %g",(double)time);
508:   dx[0] = x[0] - 1.5;
509:   dx[1] = x[1] - 1.0;
510:   r     = Norm2(dx);
511:   sigma = 0.5;
512:   u[0]  = 1 + 2*PetscExpScalar(-PetscSqr(r)/(2*PetscSqr(sigma)));
513:   u[1]  = 0.0;
514:   u[2]  = 0.0;
515:   return(0);
516: }

520: static PetscErrorCode PhysicsFunctional_SW(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)
521: {
522:   Physics      phys = (Physics)ctx;
523:   Physics_SW   *sw  = (Physics_SW*)phys->data;
524:   const SWNode *x   = (const SWNode*)xx;
525:   PetscScalar  u[2];
526:   PetscReal    h;

529:   h = PetscRealPart(x->h);
530:   Scale2(1./x->h,x->uh,u);
531:   f[sw->functional.Height] = h;
532:   f[sw->functional.Speed]  = Norm2(u) + PetscSqrtReal(sw->gravity*h);
533:   f[sw->functional.Energy] = 0.5*(Dot2(x->uh,u) + sw->gravity*PetscSqr(h));
534:   return(0);
535: }

539: static PetscErrorCode PhysicsCreate_SW(DM dm, Model mod,Physics phys)
540: {
541:   Physics_SW     *sw;

545:   phys->field_desc = PhysicsFields_SW;
546:   phys->riemann = PhysicsRiemann_SW;
547:   PetscNew(&sw);
548:   phys->data    = sw;
549:   PetscOptionsHead("SW options");
550:   {
551:     sw->gravity = 1.0;
552:     PetscOptionsReal("-sw_gravity","Gravitational constant","",sw->gravity,&sw->gravity,NULL);
553:   }
554:   PetscOptionsTail();
555:   phys->maxspeed = PetscSqrtReal(2.0*sw->gravity); /* Mach 1 for depth of 2 */

557:   {
558:     const PetscInt wallids[] = {100,101,200,300};
559:     ModelBoundaryRegister(mod,"wall",PhysicsBoundary_SW_Wall,phys,ALEN(wallids),wallids);
560:     DMPlexAddBoundary(dm, PETSC_TRUE, "wall", 0, (void (*)()) PhysicsBoundary_SW_Wall, ALEN(wallids), wallids, phys);
561:     ModelSolutionSetDefault(mod,PhysicsSolution_SW,phys);
562:     ModelFunctionalRegister(mod,"Height",&sw->functional.Height,PhysicsFunctional_SW,phys);
563:     ModelFunctionalRegister(mod,"Speed",&sw->functional.Speed,PhysicsFunctional_SW,phys);
564:     ModelFunctionalRegister(mod,"Energy",&sw->functional.Energy,PhysicsFunctional_SW,phys);
565:   }
566:   return(0);
567: }

569: /******************* Euler ********************/
570: typedef struct {
571:   PetscScalar vals[0];
572:   PetscScalar r;
573:   PetscScalar ru[DIM];
574:   PetscScalar e;
575: } EulerNode;
576: typedef PetscErrorCode (*EquationOfState)(const PetscReal*, const EulerNode*, PetscScalar*);
577: typedef struct {
578:   PetscInt        npars;
579:   PetscReal       pars[DIM];
580:   EquationOfState pressure;
581:   EquationOfState sound;
582:   struct {
583:     PetscInt Density;
584:     PetscInt Momentum;
585:     PetscInt Energy;
586:     PetscInt Pressure;
587:     PetscInt Speed;
588:   } monitor;
589: } Physics_Euler;

591: static const struct FieldDescription PhysicsFields_Euler[] = {{"Density",1},{"Momentum",DIM},{"Energy",1},{NULL,0}};

595: static PetscErrorCode Pressure_PG(const PetscReal *pars,const EulerNode *x,PetscScalar *p)
596: {
597:   PetscScalar ru2;

600:   ru2  = DotDIM(x->ru,x->ru);
601:   ru2 /= x->r;
602:   /* kinematic dof = params[0] */
603:   (*p)=2.0*(x->e-0.5*ru2)/pars[0];
604:   return(0);
605: }

609: static PetscErrorCode SpeedOfSound_PG(const PetscReal *pars,const EulerNode *x,PetscScalar *c)
610: {
611:   PetscScalar p;

614:   /* TODO remove direct usage of Pressure_PG */
615:   Pressure_PG(pars,x,&p);
616:   /* TODO check the sign of p */
617:   /* pars[1] = heat capacity ratio */
618:   (*c)=PetscSqrtScalar(pars[1]*p/x->r);
619:   return(0);
620: }

624: /*
625:  * x = (rho,rho*(u_1),...,rho*e)^T
626:  * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
627:  *
628:  * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
629:  *
630:  * */
631: static PetscErrorCode EulerFlux(Physics phys,const PetscReal *n,const EulerNode *x,EulerNode *f)
632: {
633:   Physics_Euler *eu = (Physics_Euler*)phys->data;
634:   PetscScalar   u,nu,p;
635:   PetscInt      i;

638:   u  = DotDIM(x->ru,x->ru);
639:   u /= (x->r * x->r);
640:   nu = DotDIM(x->ru,n);
641:   /* TODO check the sign of p */
642:   eu->pressure(eu->pars,x,&p);
643:   f->r = nu * x->r;
644:   for (i=0; i<DIM; i++) f->ru[i] = nu * x->ru[i] + n[i]*p;
645:   f->e = nu*(x->e+p);
646:   return(0);
647: }

649: /* PetscReal* => EulerNode* conversion */
652: static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
653: {
654:   PetscInt    i;
655:   PetscScalar xn[DIM],xt[DIM];

658:   xG[0] = xI[0];
659:   NormalSplitDIM(n,xI+1,xn,xt);
660:   for (i=0; i<DIM; i++) xG[i+1] = -xn[i]+xt[i];
661:   xG[DIM+1] = xI[DIM+1];
662:   return(0);
663: }

665: /* PetscReal* => EulerNode* conversion */
668: static PetscErrorCode PhysicsRiemann_Euler_Rusanov(Physics phys, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux)
669: {
670:   Physics_Euler   *eu = (Physics_Euler*)phys->data;
671:   PetscScalar     cL,cR,speed;
672:   const EulerNode *uL = (const EulerNode*)xL,*uR = (const EulerNode*)xR;
673:   EulerNode       fL,fR;
674:   PetscInt        i;

677:   if (uL->r < 0 || uR->r < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed density is negative");
678:   EulerFlux(phys,n,uL,&fL);
679:   EulerFlux(phys,n,uR,&fR);
680:   eu->sound(eu->pars,uL,&cL);
681:   eu->sound(eu->pars,uR,&cR);
682:   speed = PetscMax(cL,cR)+PetscMax(PetscAbsScalar(DotDIM(uL->ru,n)/NormDIM(n)),PetscAbsScalar(DotDIM(uR->ru,n)/NormDIM(n)));
683:   for (i=0; i<2+DIM; i++) flux[i] = 0.5*(fL.vals[i]+fR.vals[i])+0.5*speed*(xL[i]-xR[i]);
684:   return(0);
685: }

689: static PetscErrorCode PhysicsSolution_Euler(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
690: {
691:   PetscInt i;

694:   if (time != 0.0) SETERRQ1(mod->comm,PETSC_ERR_SUP,"No solution known for time %g",(double)time);
695:   u[0]     = 1.0;
696:   u[DIM+1] = 1.0+PetscAbsReal(x[0]);
697:   for (i=1; i<DIM+1; i++) u[i] = 0.0;
698:   return(0);
699: }

703: static PetscErrorCode PhysicsFunctional_Euler(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)
704: {
705:   Physics         phys = (Physics)ctx;
706:   Physics_Euler   *eu  = (Physics_Euler*)phys->data;
707:   const EulerNode *x   = (const EulerNode*)xx;
708:   PetscScalar     p;

711:   f[eu->monitor.Density]  = x->r;
712:   f[eu->monitor.Momentum] = NormDIM(x->ru);
713:   f[eu->monitor.Energy]   = x->e;
714:   f[eu->monitor.Speed]    = NormDIM(x->ru)/x->r;
715:   eu->pressure(eu->pars, x, &p);
716:   f[eu->monitor.Pressure] = p;
717:   return(0);
718: }

722: static PetscErrorCode PhysicsCreate_Euler(DM dm, Model mod,Physics phys)
723: {
724:   Physics_Euler   *eu;
725:   PetscErrorCode  ierr;

728:   phys->field_desc = PhysicsFields_Euler;
729:   phys->riemann = PhysicsRiemann_Euler_Rusanov;
730:   PetscNew(&eu);
731:   phys->data    = eu;
732:   PetscOptionsHead("Euler options");
733:   {
734:     eu->pars[0] = 3.0;
735:     eu->pars[1] = 1.67;
736:     PetscOptionsReal("-eu_f","Degrees of freedom","",eu->pars[0],&eu->pars[0],NULL);
737:     PetscOptionsReal("-eu_gamma","Heat capacity ratio","",eu->pars[1],&eu->pars[1],NULL);
738:   }
739:   PetscOptionsTail();
740:   eu->pressure = Pressure_PG;
741:   eu->sound    = SpeedOfSound_PG;
742:   phys->maxspeed = 1.0;
743:   {
744:     const PetscInt wallids[] = {100,101,200,300};
745:     ModelBoundaryRegister(mod,"wall",PhysicsBoundary_Euler_Wall,phys,ALEN(wallids),wallids);
746:     DMPlexAddBoundary(dm, PETSC_TRUE, "wall", 0, (void (*)()) PhysicsBoundary_Euler_Wall, ALEN(wallids), wallids, phys);
747:     ModelSolutionSetDefault(mod,PhysicsSolution_Euler,phys);
748:     ModelFunctionalRegister(mod,"Speed",&eu->monitor.Speed,PhysicsFunctional_Euler,phys);
749:     ModelFunctionalRegister(mod,"Energy",&eu->monitor.Energy,PhysicsFunctional_Euler,phys);
750:     ModelFunctionalRegister(mod,"Density",&eu->monitor.Density,PhysicsFunctional_Euler,phys);
751:     ModelFunctionalRegister(mod,"Momentum",&eu->monitor.Momentum,PhysicsFunctional_Euler,phys);
752:     ModelFunctionalRegister(mod,"Pressure",&eu->monitor.Pressure,PhysicsFunctional_Euler,phys);
753:   }
754:   return(0);
755: }

759: PetscErrorCode ConstructCellBoundary(DM dm, User user)
760: {
761:   const char     *name   = "Cell Sets";
762:   const char     *bdname = "split faces";
763:   IS             regionIS, innerIS;
764:   const PetscInt *regions, *cells;
765:   PetscInt       numRegions, innerRegion, numCells, c;

767:   PetscInt cStart, cEnd, fStart, fEnd;

769:   PetscBool      hasLabel;

773:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
774:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);

776:   DMPlexHasLabel(dm, name, &hasLabel);
777:   if (!hasLabel) return(0);
778:   DMPlexGetLabelSize(dm, name, &numRegions);
779:   if (numRegions != 2) return(0);
780:   /* Get the inner id */
781:   DMPlexGetLabelIdIS(dm, name, &regionIS);
782:   ISGetIndices(regionIS, &regions);
783:   innerRegion = regions[0];
784:   ISRestoreIndices(regionIS, &regions);
785:   ISDestroy(&regionIS);
786:   /* Find the faces between cells in different regions, could call DMPlexCreateNeighborCSR() */
787:   DMPlexGetStratumIS(dm, name, innerRegion, &innerIS);
788:   ISGetLocalSize(innerIS, &numCells);
789:   ISGetIndices(innerIS, &cells);
790:   DMPlexCreateLabel(dm, bdname);
791:   for (c = 0; c < numCells; ++c) {
792:     const PetscInt cell = cells[c];
793:     const PetscInt *faces;
794:     PetscInt       numFaces, f;

796:     if ((cell < cStart) || (cell >= cEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", cell);
797:     DMPlexGetConeSize(dm, cell, &numFaces);
798:     DMPlexGetCone(dm, cell, &faces);
799:     for (f = 0; f < numFaces; ++f) {
800:       const PetscInt face = faces[f];
801:       const PetscInt *neighbors;
802:       PetscInt       nC, regionA, regionB;

804:       if ((face < fStart) || (face >= fEnd)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a face", face);
805:       DMPlexGetSupportSize(dm, face, &nC);
806:       if (nC != 2) continue;
807:       DMPlexGetSupport(dm, face, &neighbors);
808:       if ((neighbors[0] >= user->cEndInterior) || (neighbors[1] >= user->cEndInterior)) continue;
809:       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]);
810:       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]);
811:       DMPlexGetLabelValue(dm, name, neighbors[0], &regionA);
812:       DMPlexGetLabelValue(dm, name, neighbors[1], &regionB);
813:       if (regionA < 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[0]);
814:       if (regionB < 0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[1]);
815:       if (regionA != regionB) {
816:         DMPlexSetLabelValue(dm, bdname, faces[f], 1);
817:       }
818:     }
819:   }
820:   ISRestoreIndices(innerIS, &cells);
821:   ISDestroy(&innerIS);
822:   {
823:     DMLabel label;

825:     PetscViewerASCIISynchronizedAllow(PETSC_VIEWER_STDOUT_WORLD, PETSC_TRUE);
826:     DMPlexGetLabel(dm, bdname, &label);
827:     DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD);
828:   }
829:   return(0);
830: }

834: /* Right now, I have just added duplicate faces, which see both cells. We can
835: - Add duplicate vertices and decouple the face cones
836: - Disconnect faces from cells across the rotation gap
837: */
838: PetscErrorCode SplitFaces(DM *dmSplit, const char labelName[], User user)
839: {
840:   DM             dm = *dmSplit, sdm;
841:   PetscSF        sfPoint, gsfPoint;
842:   PetscSection   coordSection, newCoordSection;
843:   Vec            coordinates;
844:   IS             idIS;
845:   const PetscInt *ids;
846:   PetscInt       *newpoints;
847:   PetscInt       dim, depth, maxConeSize, maxSupportSize, numLabels;
848:   PetscInt       numFS, fs, pStart, pEnd, p, vStart, vEnd, v, fStart, fEnd, newf, d, l;
849:   PetscBool      hasLabel;

853:   DMPlexHasLabel(dm, labelName, &hasLabel);
854:   if (!hasLabel) return(0);
855:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
856:   DMSetType(sdm, DMPLEX);
857:   DMPlexGetDimension(dm, &dim);
858:   DMPlexSetDimension(sdm, dim);

860:   DMPlexGetLabelIdIS(dm, labelName, &idIS);
861:   ISGetLocalSize(idIS, &numFS);
862:   ISGetIndices(idIS, &ids);

864:   user->numSplitFaces = 0;
865:   for (fs = 0; fs < numFS; ++fs) {
866:     PetscInt numBdFaces;

868:     DMPlexGetStratumSize(dm, labelName, ids[fs], &numBdFaces);
869:     user->numSplitFaces += numBdFaces;
870:   }
871:   DMPlexGetChart(dm, &pStart, &pEnd);
872:   pEnd += user->numSplitFaces;
873:   DMPlexSetChart(sdm, pStart, pEnd);
874:   /* Set cone and support sizes */
875:   DMPlexGetDepth(dm, &depth);
876:   for (d = 0; d <= depth; ++d) {
877:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
878:     for (p = pStart; p < pEnd; ++p) {
879:       PetscInt newp = p;
880:       PetscInt size;

882:       DMPlexGetConeSize(dm, p, &size);
883:       DMPlexSetConeSize(sdm, newp, size);
884:       DMPlexGetSupportSize(dm, p, &size);
885:       DMPlexSetSupportSize(sdm, newp, size);
886:     }
887:   }
888:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
889:   for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
890:     IS             faceIS;
891:     const PetscInt *faces;
892:     PetscInt       numFaces, f;

894:     DMPlexGetStratumIS(dm, labelName, ids[fs], &faceIS);
895:     ISGetLocalSize(faceIS, &numFaces);
896:     ISGetIndices(faceIS, &faces);
897:     for (f = 0; f < numFaces; ++f, ++newf) {
898:       PetscInt size;

900:       /* Right now I think that both faces should see both cells */
901:       DMPlexGetConeSize(dm, faces[f], &size);
902:       DMPlexSetConeSize(sdm, newf, size);
903:       DMPlexGetSupportSize(dm, faces[f], &size);
904:       DMPlexSetSupportSize(sdm, newf, size);
905:     }
906:     ISRestoreIndices(faceIS, &faces);
907:     ISDestroy(&faceIS);
908:   }
909:   DMSetUp(sdm);
910:   /* Set cones and supports */
911:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
912:   PetscMalloc(PetscMax(maxConeSize, maxSupportSize) * sizeof(PetscInt), &newpoints);
913:   DMPlexGetChart(dm, &pStart, &pEnd);
914:   for (p = pStart; p < pEnd; ++p) {
915:     const PetscInt *points, *orientations;
916:     PetscInt       size, i, newp = p;

918:     DMPlexGetConeSize(dm, p, &size);
919:     DMPlexGetCone(dm, p, &points);
920:     DMPlexGetConeOrientation(dm, p, &orientations);
921:     for (i = 0; i < size; ++i) newpoints[i] = points[i];
922:     DMPlexSetCone(sdm, newp, newpoints);
923:     DMPlexSetConeOrientation(sdm, newp, orientations);
924:     DMPlexGetSupportSize(dm, p, &size);
925:     DMPlexGetSupport(dm, p, &points);
926:     for (i = 0; i < size; ++i) newpoints[i] = points[i];
927:     DMPlexSetSupport(sdm, newp, newpoints);
928:   }
929:   PetscFree(newpoints);
930:   for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
931:     IS             faceIS;
932:     const PetscInt *faces;
933:     PetscInt       numFaces, f;

935:     DMPlexGetStratumIS(dm, labelName, ids[fs], &faceIS);
936:     ISGetLocalSize(faceIS, &numFaces);
937:     ISGetIndices(faceIS, &faces);
938:     for (f = 0; f < numFaces; ++f, ++newf) {
939:       const PetscInt *points;

941:       DMPlexGetCone(dm, faces[f], &points);
942:       DMPlexSetCone(sdm, newf, points);
943:       DMPlexGetSupport(dm, faces[f], &points);
944:       DMPlexSetSupport(sdm, newf, points);
945:     }
946:     ISRestoreIndices(faceIS, &faces);
947:     ISDestroy(&faceIS);
948:   }
949:   ISRestoreIndices(idIS, &ids);
950:   ISDestroy(&idIS);
951:   DMPlexStratify(sdm);
952:   /* Convert coordinates */
953:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
954:   DMGetCoordinateSection(dm, &coordSection);
955:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
956:   PetscSectionSetNumFields(newCoordSection, 1);
957:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
958:   PetscSectionSetChart(newCoordSection, vStart, vEnd);
959:   for (v = vStart; v < vEnd; ++v) {
960:     PetscSectionSetDof(newCoordSection, v, dim);
961:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
962:   }
963:   PetscSectionSetUp(newCoordSection);
964:   DMSetCoordinateSection(sdm, newCoordSection);
965:   PetscSectionDestroy(&newCoordSection); /* relinquish our reference */
966:   DMGetCoordinatesLocal(dm, &coordinates);
967:   DMSetCoordinatesLocal(sdm, coordinates);
968:   /* Convert labels */
969:   DMPlexGetNumLabels(dm, &numLabels);
970:   for (l = 0; l < numLabels; ++l) {
971:     const char *lname;
972:     PetscBool  isDepth;

974:     DMPlexGetLabelName(dm, l, &lname);
975:     PetscStrcmp(lname, "depth", &isDepth);
976:     if (isDepth) continue;
977:     DMPlexCreateLabel(sdm, lname);
978:     DMPlexGetLabelIdIS(dm, lname, &idIS);
979:     ISGetLocalSize(idIS, &numFS);
980:     ISGetIndices(idIS, &ids);
981:     for (fs = 0; fs < numFS; ++fs) {
982:       IS             pointIS;
983:       const PetscInt *points;
984:       PetscInt       numPoints;

986:       DMPlexGetStratumIS(dm, lname, ids[fs], &pointIS);
987:       ISGetLocalSize(pointIS, &numPoints);
988:       ISGetIndices(pointIS, &points);
989:       for (p = 0; p < numPoints; ++p) {
990:         PetscInt newpoint = points[p];

992:         DMPlexSetLabelValue(sdm, lname, newpoint, ids[fs]);
993:       }
994:       ISRestoreIndices(pointIS, &points);
995:       ISDestroy(&pointIS);
996:     }
997:     ISRestoreIndices(idIS, &ids);
998:     ISDestroy(&idIS);
999:   }
1000:   /* Convert pointSF */
1001:   const PetscSFNode *remotePoints;
1002:   PetscSFNode       *gremotePoints;
1003:   const PetscInt    *localPoints;
1004:   PetscInt          *glocalPoints,*newLocation,*newRemoteLocation;
1005:   PetscInt          numRoots, numLeaves;
1006:   PetscMPIInt       numProcs;

1008:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
1009:   DMGetPointSF(dm, &sfPoint);
1010:   DMGetPointSF(sdm, &gsfPoint);
1011:   DMPlexGetChart(dm,&pStart,&pEnd);
1012:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
1013:   if (numRoots >= 0) {
1014:     PetscMalloc2(numRoots,&newLocation,pEnd-pStart,&newRemoteLocation);
1015:     for (l=0; l<numRoots; l++) newLocation[l] = l; /* + (l >= cEnd ? user->numGhostCells : 0); */
1016:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
1017:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation);
1018:     PetscMalloc1(numLeaves,    &glocalPoints);
1019:     PetscMalloc1(numLeaves, &gremotePoints);
1020:     for (l = 0; l < numLeaves; ++l) {
1021:       glocalPoints[l]        = localPoints[l]; /* localPoints[l] >= cEnd ? localPoints[l] + user->numGhostCells : localPoints[l]; */
1022:       gremotePoints[l].rank  = remotePoints[l].rank;
1023:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
1024:     }
1025:     PetscFree2(newLocation,newRemoteLocation);
1026:     PetscSFSetGraph(gsfPoint, numRoots+user->numGhostCells, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
1027:   }
1028:   DMDestroy(dmSplit);
1029:   *dmSplit = sdm;
1030:   return(0);
1031: }

1035: static PetscErrorCode IsExteriorOrGhostFace(DM dm,PetscInt face,PetscBool *isghost)
1036: {
1038:   PetscInt       ghost,boundary;

1041:   *isghost = PETSC_FALSE;
1042:   DMPlexGetLabelValue(dm, "ghost", face, &ghost);
1043:   DMPlexGetLabelValue(dm, "Face Sets", face, &boundary);
1044:   if (ghost >= 0 || boundary >= 0) *isghost = PETSC_TRUE;
1045:   return(0);
1046: }

1050: /* Overwrites A. Can only handle full-rank problems with m>=n */
1051: static PetscErrorCode PseudoInverse(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1052: {
1053:   PetscBool      debug = PETSC_FALSE;
1055:   PetscBLASInt   M,N,K,lda,ldb,ldwork,info;
1056:   PetscScalar    *R,*Q,*Aback,Alpha;

1059:   if (debug) {
1060:     PetscMalloc1(m*n,&Aback);
1061:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1062:   }

1064:   PetscBLASIntCast(m,&M);
1065:   PetscBLASIntCast(n,&N);
1066:   PetscBLASIntCast(mstride,&lda);
1067:   PetscBLASIntCast(worksize,&ldwork);
1068:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1069:   LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
1070:   PetscFPTrapPop();
1071:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
1072:   R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */

1074:   /* Extract an explicit representation of Q */
1075:   Q    = Ainv;
1076:   PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
1077:   K    = N;                     /* full rank */
1078:   LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
1079:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");

1081:   /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
1082:   Alpha = 1.0;
1083:   ldb   = lda;
1084:   BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
1085:   /* Ainv is Q, overwritten with inverse */

1087:   if (debug) {                      /* Check that pseudo-inverse worked */
1088:     PetscScalar Beta = 0.0;
1089:     PetscInt    ldc;
1090:     K   = N;
1091:     ldc = N;
1092:     BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc);
1093:     PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);
1094:     PetscFree(Aback);
1095:   }
1096:   return(0);
1097: }

1101: static PetscErrorCode PseudoInverseGetWorkRequired(PetscInt maxFaces,PetscInt *work)
1102: {
1103:   PetscInt m,n,nrhs,minwork;

1106:   m       = maxFaces;
1107:   n       = DIM;
1108:   nrhs    = maxFaces;
1109:   minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
1110:   *work   = 5*minwork;          /* We can afford to be extra generous */
1111:   return(0);
1112: }

1116: /* Overwrites A. Can handle degenerate problems and m<n. */
1117: static PetscErrorCode PseudoInverseSVD(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
1118: {
1119:   PetscBool      debug = PETSC_FALSE;
1121:   PetscInt       i,j,maxmn;
1122:   PetscBLASInt   M,N,nrhs,lda,ldb,irank,ldwork,info;
1123:   PetscScalar    rcond,*tmpwork,*Brhs,*Aback;

1126:   if (debug) {
1127:     PetscMalloc1(m*n,&Aback);
1128:     PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));
1129:   }

1131:   /* initialize to identity */
1132:   tmpwork = Ainv;
1133:   Brhs = work;
1134:   maxmn = PetscMax(m,n);
1135:   for (j=0; j<maxmn; j++) {
1136:     for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
1137:   }

1139:   PetscBLASIntCast(m,&M);
1140:   PetscBLASIntCast(n,&N);
1141:   nrhs  = M;
1142:   PetscBLASIntCast(mstride,&lda);
1143:   PetscBLASIntCast(maxmn,&ldb);
1144:   PetscBLASIntCast(worksize,&ldwork);
1145:   rcond = -1;
1146:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1147:   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb,tau,&rcond,&irank,tmpwork,&ldwork,&info);
1148:   PetscFPTrapPop();
1149:   if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error");
1150:   /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */
1151:   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");

1153:   /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn.
1154:    * Here we transpose to (N,nrhs) row-major rowdim=mstride. */
1155:   for (i=0; i<n; i++) {
1156:     for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn];
1157:   }
1158:   return(0);
1159: }

1163: /* Build least squares gradient reconstruction operators */
1164: static PetscErrorCode BuildLeastSquares(DM dm,PetscInt cEndInterior,DM dmFace,PetscScalar *fgeom,DM dmCell,PetscScalar *cgeom)
1165: {
1167:   PetscInt       c,cStart,cEnd,maxNumFaces,worksize;
1168:   PetscScalar    *B,*Binv,*work,*tau,**gref;

1171:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
1172:   DMPlexGetMaxSizes(dm,&maxNumFaces,NULL);
1173:   PseudoInverseGetWorkRequired(maxNumFaces,&worksize);
1174:   PetscMalloc5(maxNumFaces*DIM,&B,worksize,&Binv,worksize,&work,maxNumFaces,&tau,maxNumFaces,&gref);
1175:   for (c=cStart; c<cEndInterior; c++) {
1176:     const PetscInt *faces;
1177:     PetscInt       numFaces,usedFaces,f,i,j;
1178:     const CellGeom *cg;
1179:     PetscBool      ghost;
1180:     DMPlexGetConeSize(dm,c,&numFaces);
1181:     if (numFaces < DIM) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction",c,numFaces);
1182:     DMPlexGetCone(dm,c,&faces);
1183:     DMPlexPointLocalRead(dmCell,c,cgeom,&cg);
1184:     for (f=0,usedFaces=0; f<numFaces; f++) {
1185:       const PetscInt *fcells;
1186:       PetscInt       ncell,side;
1187:       FaceGeom       *fg;
1188:       const CellGeom *cg1;
1189:       IsExteriorOrGhostFace(dm,faces[f],&ghost);
1190:       if (ghost) continue;
1191:       DMPlexGetSupport(dm,faces[f],&fcells);
1192:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
1193:       ncell = fcells[!side];   /* the neighbor */
1194:       DMPlexPointLocalRef(dmFace,faces[f],fgeom,&fg);
1195:       DMPlexPointLocalRead(dmCell,ncell,cgeom,&cg1);
1196:       for (j=0; j<DIM; j++) B[j*numFaces+usedFaces] = cg1->centroid[j] - cg->centroid[j];
1197:       gref[usedFaces++] = fg->grad[side];  /* Gradient reconstruction term will go here */
1198:     }
1199:     if (!usedFaces) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Mesh contains isolated cell (no neighbors). Is it intentional?");
1200:     /* Overwrites B with garbage, returns Binv in row-major format */
1201:     if (0) {
1202:       PseudoInverse(usedFaces,numFaces,DIM,B,Binv,tau,worksize,work);
1203:     } else {
1204:       PseudoInverseSVD(usedFaces,numFaces,DIM,B,Binv,tau,worksize,work);
1205:     }
1206:     for (f=0,i=0; f<numFaces; f++) {
1207:       IsExteriorOrGhostFace(dm,faces[f],&ghost);
1208:       if (ghost) continue;
1209:       for (j=0; j<DIM; j++) gref[i][j] = Binv[j*numFaces+i];
1210:       i++;
1211:     }

1213: #if 1
1214:     if (0) {
1215:       PetscReal grad[2] = {0,0};
1216:       for (f=0; f<numFaces; f++) {
1217:         const PetscInt *fcells;
1218:         const CellGeom *cg1;
1219:         const FaceGeom *fg;
1220:         DMPlexGetSupport(dm,faces[f],&fcells);
1221:         DMPlexPointLocalRead(dmFace,faces[f],fgeom,&fg);
1222:         for (i=0; i<2; i++) {
1223:           if (fcells[i] == c) continue;
1224:           DMPlexPointLocalRead(dmCell,fcells[i],cgeom,&cg1);
1225:           PetscScalar du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]);
1226:           grad[0] += fg->grad[!i][0] * du;
1227:           grad[1] += fg->grad[!i][1] * du;
1228:         }
1229:       }
1230:       printf("cell[%d] grad (%g,%g)\n",c,grad[0],grad[1]);
1231:     }
1232: #endif
1233:   }
1234:   PetscFree5(B,Binv,work,tau,gref);
1235:   return(0);
1236: }

1240: /* Set up face data and cell data */
1241: PetscErrorCode ConstructGeometry(DM dm, Vec *facegeom, Vec *cellgeom, User user)
1242: {
1243:   DM             dmFace, dmCell;
1244:   DMLabel        ghostLabel;
1245:   PetscSection   sectionFace, sectionCell;
1246:   PetscSection   coordSection;
1247:   Vec            coordinates;
1248:   PetscReal      minradius;
1249:   PetscScalar    *fgeom, *cgeom;
1250:   PetscInt       dim, cStart, cEnd, c, fStart, fEnd, f;

1254:   DMPlexGetDimension(dm, &dim);
1255:   if (dim != DIM) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No support for dim %D != DIM %D",dim,DIM);
1256:   DMGetCoordinateSection(dm, &coordSection);
1257:   DMGetCoordinatesLocal(dm, &coordinates);

1259:   /* Make cell centroids and volumes */
1260:   DMClone(dm, &dmCell);
1261:   DMSetCoordinateSection(dmCell, coordSection);
1262:   DMSetCoordinatesLocal(dmCell, coordinates);
1263:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);
1264:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1265:   PetscSectionSetChart(sectionCell, cStart, cEnd);
1266:   for (c = cStart; c < cEnd; ++c) {
1267:     PetscSectionSetDof(sectionCell, c, sizeof(CellGeom)/sizeof(PetscScalar));
1268:   }
1269:   PetscSectionSetUp(sectionCell);
1270:   DMSetDefaultSection(dmCell, sectionCell);
1271:   PetscSectionDestroy(&sectionCell); /* relinquish our reference */

1273:   DMCreateLocalVector(dmCell, cellgeom);
1274:   VecGetArray(*cellgeom, &cgeom);
1275:   for (c = cStart; c < user->cEndInterior; ++c) {
1276:     CellGeom *cg;

1278:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
1279:     PetscMemzero(cg,sizeof(*cg));
1280:     DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
1281:   }
1282:   /* Compute face normals and minimum cell radius */
1283:   DMClone(dm, &dmFace);
1284:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionFace);
1285:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1286:   PetscSectionSetChart(sectionFace, fStart, fEnd);
1287:   for (f = fStart; f < fEnd; ++f) {
1288:     PetscSectionSetDof(sectionFace, f, sizeof(FaceGeom)/sizeof(PetscScalar));
1289:   }
1290:   PetscSectionSetUp(sectionFace);
1291:   DMSetDefaultSection(dmFace, sectionFace);
1292:   PetscSectionDestroy(&sectionFace);
1293:   DMCreateLocalVector(dmFace, facegeom);
1294:   VecGetArray(*facegeom, &fgeom);
1295:   DMPlexGetLabel(dm, "ghost", &ghostLabel);
1296:   minradius = PETSC_MAX_REAL;
1297:   for (f = fStart; f < fEnd; ++f) {
1298:     FaceGeom *fg;
1299:     PetscReal area;
1300:     PetscInt  ghost, d;

1302:     DMLabelGetValue(ghostLabel, f, &ghost);
1303:     if (ghost >= 0) continue;
1304:     DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
1305:     DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
1306:     for (d = 0; d < dim; ++d) fg->normal[d] *= area;
1307:     /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
1308:     {
1309:       CellGeom       *cL, *cR;
1310:       const PetscInt *cells;
1311:       PetscReal      *lcentroid, *rcentroid;
1312:       PetscScalar     v[3];

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:   MPI_Allreduce(&minradius, &user->minradius, 1, MPIU_SCALAR, MPI_MIN, PetscObjectComm((PetscObject)dm));
1343:   /* Compute centroids of ghost cells */
1344:   for (c = user->cEndInterior; c < cEnd; ++c) {
1345:     FaceGeom       *fg;
1346:     const PetscInt *cone,    *support;
1347:     PetscInt        coneSize, supportSize, s;

1349:     DMPlexGetConeSize(dmCell, c, &coneSize);
1350:     if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
1351:     DMPlexGetCone(dmCell, c, &cone);
1352:     DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
1353:     if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize);
1354:     DMPlexGetSupport(dmCell, cone[0], &support);
1355:     DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
1356:     for (s = 0; s < 2; ++s) {
1357:       /* Reflect ghost centroid across plane of face */
1358:       if (support[s] == c) {
1359:         const CellGeom *ci;
1360:         CellGeom       *cg;
1361:         PetscScalar     c2f[3], a;

1363:         DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
1364:         WaxpyD(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
1365:         a    = DotD(dim, c2f, fg->normal)/DotD(dim, fg->normal, fg->normal);
1366:         DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
1367:         WaxpyD(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
1368:         cg->volume = ci->volume;
1369:       }
1370:     }
1371:   }
1372:   VecRestoreArray(*facegeom, &fgeom);
1373:   VecRestoreArray(*cellgeom, &cgeom);
1374:   DMDestroy(&dmCell);
1375:   DMDestroy(&dmFace);
1376:   return(0);
1377: }

1381: PetscErrorCode ReconstructionSetUp(DM dm, User user)
1382: {
1383:   DM             dmFace, dmCell;
1384:   PetscSection   sectionGrad;
1385:   PetscScalar   *fgeom, *cgeom;
1386:   PetscInt       cStart, cEnd, c;

1390:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1391:   VecGetDM(user->facegeom, &dmFace);
1392:   VecGetDM(user->cellgeom, &dmCell);
1393:   VecGetArray(user->facegeom, &fgeom);
1394:   VecGetArray(user->cellgeom, &cgeom);
1395:   BuildLeastSquares(dm, user->cEndInterior, dmFace, fgeom, dmCell, cgeom);
1396:   DMClone(dm, &user->dmGrad);
1397:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionGrad);
1398:   PetscSectionSetChart(sectionGrad, cStart, cEnd);
1399:   for (c = cStart; c < cEnd; ++c) {
1400:     PetscSectionSetDof(sectionGrad, c, user->model->physics->dof*DIM);
1401:   }
1402:   PetscSectionSetUp(sectionGrad);
1403:   DMSetDefaultSection(user->dmGrad, sectionGrad);
1404:   PetscSectionDestroy(&sectionGrad);
1405:   VecRestoreArray(user->facegeom, &fgeom);
1406:   VecRestoreArray(user->cellgeom, &cgeom);
1407:   return(0);
1408: }

1412: PetscErrorCode CreateMesh(MPI_Comm comm, const char *filename, PetscInt overlap, PetscBool splitFaces, DM *newdm, User user)
1413: {
1414:   DM             dm, dmDist;

1418:   DMPlexCreateExodusFromFile(comm, filename, PETSC_TRUE, &dm);
1419:   /* Distribute mesh */
1420:   DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE);
1421:   DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE);
1422:   DMPlexDistribute(dm, "chaco", overlap, NULL, &dmDist);
1423:   if (dmDist) {
1424:     DMDestroy(&dm);
1425:     dm   = dmDist;
1426:   }
1427:   DMSetFromOptions(dm);
1428:   {
1429:     DM gdm;

1431:     DMPlexGetHeightStratum(dm, 0, NULL, &user->cEndInterior);
1432:     DMPlexConstructGhostCells(dm, NULL, &user->numGhostCells, &gdm);
1433:     DMDestroy(&dm);
1434:     dm   = gdm;
1435:     DMViewFromOptions(dm, NULL, "-dm_view");
1436:   }
1437:   if (splitFaces) {ConstructCellBoundary(dm, user);}
1438:   SplitFaces(&dm, "split faces", user);
1439:   ConstructGeometry(dm, &user->facegeom, &user->cellgeom, user);
1440:   PetscObjectCompose((PetscObject) dm, "FaceGeometry", (PetscObject) user->facegeom);
1441:   if (0) {VecView(user->cellgeom, PETSC_VIEWER_STDOUT_WORLD);}
1442:   *newdm = dm;
1443:   return(0);
1444: }

1448: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
1449: {
1450:   PetscSF        sfPoint;
1451:   PetscSection   coordSection;
1452:   Vec            coordinates;
1453:   PetscSection   sectionCell;
1454:   PetscScalar    *part;
1455:   PetscInt       cStart, cEnd, c;
1456:   PetscMPIInt    rank;

1460:   DMGetCoordinateSection(dm, &coordSection);
1461:   DMGetCoordinatesLocal(dm, &coordinates);
1462:   DMClone(dm, dmCell);
1463:   DMGetPointSF(dm, &sfPoint);
1464:   DMSetPointSF(*dmCell, sfPoint);
1465:   DMSetCoordinateSection(*dmCell, coordSection);
1466:   DMSetCoordinatesLocal(*dmCell, coordinates);
1467:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1468:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);
1469:   DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd);
1470:   PetscSectionSetChart(sectionCell, cStart, cEnd);
1471:   for (c = cStart; c < cEnd; ++c) {
1472:     PetscSectionSetDof(sectionCell, c, 1);
1473:   }
1474:   PetscSectionSetUp(sectionCell);
1475:   DMSetDefaultSection(*dmCell, sectionCell);
1476:   PetscSectionDestroy(&sectionCell);
1477:   DMCreateLocalVector(*dmCell, partition);
1478:   PetscObjectSetName((PetscObject)*partition, "partition");
1479:   VecGetArray(*partition, &part);
1480:   for (c = cStart; c < cEnd; ++c) {
1481:     PetscScalar *p;

1483:     DMPlexPointLocalRef(*dmCell, c, part, &p);
1484:     p[0] = rank;
1485:   }
1486:   VecRestoreArray(*partition, &part);
1487:   return(0);
1488: }

1492: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
1493: {
1494:   DM                dmMass, dmFace, dmCell, dmCoord;
1495:   PetscSection      coordSection;
1496:   Vec               coordinates;
1497:   PetscSection      sectionMass;
1498:   PetscScalar       *m;
1499:   const PetscScalar *facegeom, *cellgeom, *coords;
1500:   PetscInt          vStart, vEnd, v;
1501:   PetscErrorCode    ierr;

1504:   DMGetCoordinateSection(dm, &coordSection);
1505:   DMGetCoordinatesLocal(dm, &coordinates);
1506:   DMClone(dm, &dmMass);
1507:   DMSetCoordinateSection(dmMass, coordSection);
1508:   DMSetCoordinatesLocal(dmMass, coordinates);
1509:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass);
1510:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1511:   PetscSectionSetChart(sectionMass, vStart, vEnd);
1512:   for (v = vStart; v < vEnd; ++v) {
1513:     PetscInt numFaces;

1515:     DMPlexGetSupportSize(dmMass, v, &numFaces);
1516:     PetscSectionSetDof(sectionMass, v, numFaces*numFaces);
1517:   }
1518:   PetscSectionSetUp(sectionMass);
1519:   DMSetDefaultSection(dmMass, sectionMass);
1520:   PetscSectionDestroy(&sectionMass);
1521:   DMGetLocalVector(dmMass, massMatrix);
1522:   VecGetArray(*massMatrix, &m);
1523:   VecGetDM(user->facegeom, &dmFace);
1524:   VecGetArrayRead(user->facegeom, &facegeom);
1525:   VecGetDM(user->cellgeom, &dmCell);
1526:   VecGetArrayRead(user->cellgeom, &cellgeom);
1527:   DMGetCoordinateDM(dm, &dmCoord);
1528:   VecGetArrayRead(coordinates, &coords);
1529:   for (v = vStart; v < vEnd; ++v) {
1530:     const PetscInt    *faces;
1531:     const FaceGeom    *fgA, *fgB, *cg;
1532:     const PetscScalar *vertex;
1533:     PetscInt          numFaces, sides[2], f, g;

1535:     DMPlexPointLocalRead(dmCoord, v, coords, &vertex);
1536:     DMPlexGetSupportSize(dmMass, v, &numFaces);
1537:     DMPlexGetSupport(dmMass, v, &faces);
1538:     for (f = 0; f < numFaces; ++f) {
1539:       sides[0] = faces[f];
1540:       DMPlexPointLocalRead(dmFace, faces[f], facegeom, &fgA);
1541:       for (g = 0; g < numFaces; ++g) {
1542:         const PetscInt *cells = NULL;;
1543:         PetscReal      area   = 0.0;
1544:         PetscInt       numCells;

1546:         sides[1] = faces[g];
1547:         DMPlexPointLocalRead(dmFace, faces[g], facegeom, &fgB);
1548:         DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells);
1549:         if (numCells != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
1550:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
1551:         area += PetscAbsScalar((vertex[0] - cg->centroid[0])*(fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1])*(fgA->centroid[0] - cg->centroid[0]));
1552:         area += PetscAbsScalar((vertex[0] - cg->centroid[0])*(fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1])*(fgB->centroid[0] - cg->centroid[0]));
1553:         m[f*numFaces+g] = Dot2(fgA->normal, fgB->normal)*area*0.5;
1554:         DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells);
1555:       }
1556:     }
1557:   }
1558:   VecRestoreArrayRead(user->facegeom, &facegeom);
1559:   VecRestoreArrayRead(user->facegeom, &facegeom);
1560:   VecRestoreArrayRead(coordinates, &coords);
1561:   VecRestoreArray(*massMatrix, &m);
1562:   DMDestroy(&dmMass);
1563:   return(0);
1564: }

1568: PetscErrorCode SetUpLocalSpace(DM dm, User user)
1569: {
1570:   PetscSection   stateSection;
1571:   Physics        phys;
1572:   PetscInt       dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, c, i;

1576:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1577:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection);
1578:   phys = user->model->physics;
1579:   PetscSectionSetNumFields(stateSection,phys->nfields);
1580:   for (i=0; i<phys->nfields; i++) {
1581:     PetscSectionSetFieldName(stateSection,i,phys->field_desc[i].name);
1582:     PetscSectionSetFieldComponents(stateSection,i,phys->field_desc[i].dof);
1583:   }
1584:   PetscSectionSetChart(stateSection, cStart, cEnd);
1585:   for (c = cStart; c < cEnd; ++c) {
1586:     for (i=0; i<phys->nfields; i++) {
1587:       PetscSectionSetFieldDof(stateSection,c,i,phys->field_desc[i].dof);
1588:     }
1589:     PetscSectionSetDof(stateSection, c, dof);
1590:   }
1591:   for (c = user->cEndInterior; c < cEnd; ++c) {
1592:     PetscSectionSetConstraintDof(stateSection, c, dof);
1593:   }
1594:   PetscSectionSetUp(stateSection);
1595:   PetscMalloc1(dof, &cind);
1596:   for (d = 0; d < dof; ++d) cind[d] = d;
1597: #if 0
1598:   for (c = cStart; c < cEnd; ++c) {
1599:     PetscInt val;

1601:     DMPlexGetLabelValue(dm, "vtk", c, &val);
1602:     if (val < 0) {PetscSectionSetConstraintIndices(stateSection, c, cind);}
1603:   }
1604: #endif
1605:   for (c = user->cEndInterior; c < cEnd; ++c) {
1606:     PetscSectionSetConstraintIndices(stateSection, c, cind);
1607:   }
1608:   PetscFree(cind);
1609:   PetscSectionGetStorageSize(stateSection, &stateSize);
1610:   DMSetDefaultSection(dm,stateSection);
1611:   PetscSectionDestroy(&stateSection);
1612:   return(0);
1613: }

1617: PetscErrorCode SetUpBoundaries(DM dm, User user)
1618: {
1619:   Model          mod = user->model;
1621:   BoundaryLink   b;

1624:   PetscOptionsBegin(PetscObjectComm((PetscObject)dm),NULL,"Boundary condition options","");
1625:   for (b = mod->boundary; b; b=b->next) {
1626:     char      optname[512];
1627:     PetscInt  ids[512],len = 512;
1628:     PetscBool flg;
1629:     PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name);
1630:     PetscMemzero(ids,sizeof(ids));
1631:     PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg);
1632:     if (flg) {
1633:       /* TODO: check all IDs to make sure they exist in the mesh */
1634:       PetscFree(b->ids);
1635:       b->numids = len;
1636:       PetscMalloc1(len,&b->ids);
1637:       PetscMemcpy(b->ids,ids,len*sizeof(PetscInt));
1638:     }
1639:   }
1640:   PetscOptionsEnd();
1641:   return(0);
1642: }

1646: /* The ids are just defaults, can be overridden at command line. I expect to set defaults based on names in the future. */
1647: static PetscErrorCode ModelBoundaryRegister(Model mod,const char *name,BoundaryFunction bcFunc,void *ctx,PetscInt numids,const PetscInt *ids)
1648: {
1650:   BoundaryLink   link;

1653:   PetscNew(&link);
1654:   PetscStrallocpy(name,&link->name);
1655:   link->numids  = numids;
1656:   PetscMalloc1(numids,&link->ids);
1657:   PetscMemcpy(link->ids,ids,numids*sizeof(PetscInt));
1658:   link->func    = bcFunc;
1659:   link->ctx     = ctx;
1660:   link->next    = mod->boundary;
1661:   mod->boundary = link;
1662:   return(0);
1663: }

1667: static PetscErrorCode BoundaryLinkDestroy(BoundaryLink *link)
1668: {
1670:   BoundaryLink   l,next;

1673:   if (!link) return(0);
1674:   l     = *link;
1675:   *link = NULL;
1676:   for (; l; l=next) {
1677:     next = l->next;
1678:     PetscFree(l->ids);
1679:     PetscFree(l->name);
1680:     PetscFree(l);
1681:   }
1682:   return(0);
1683: }

1687: static PetscErrorCode ModelBoundaryFind(Model mod,PetscInt id,BoundaryFunction *bcFunc,void **ctx)
1688: {
1689:   BoundaryLink link;
1690:   PetscInt     i;

1693:   *bcFunc = NULL;
1694:   for (link=mod->boundary; link; link=link->next) {
1695:     for (i=0; i<link->numids; i++) {
1696:       if (link->ids[i] == id) {
1697:         *bcFunc = link->func;
1698:         *ctx    = link->ctx;
1699:         return(0);
1700:       }
1701:     }
1702:   }
1703:   SETERRQ1(mod->comm,PETSC_ERR_USER,"Boundary ID %D not associated with any registered boundary condition",id);
1704:   return(0);
1705: }
1708: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
1709: static PetscErrorCode ModelSolutionSetDefault(Model mod,SolutionFunction func,void *ctx)
1710: {
1712:   mod->solution    = func;
1713:   mod->solutionctx = ctx;
1714:   return(0);
1715: }

1719: static PetscErrorCode ModelFunctionalRegister(Model mod,const char *name,PetscInt *offset,FunctionalFunction func,void *ctx)
1720: {
1722:   FunctionalLink link,*ptr;
1723:   PetscInt       lastoffset = -1;

1726:   for (ptr=&mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
1727:   PetscNew(&link);
1728:   PetscStrallocpy(name,&link->name);
1729:   link->offset = lastoffset + 1;
1730:   link->func   = func;
1731:   link->ctx    = ctx;
1732:   link->next   = NULL;
1733:   *ptr         = link;
1734:   *offset      = link->offset;
1735:   return(0);
1736: }

1740: static PetscErrorCode ModelFunctionalSetFromOptions(Model mod)
1741: {
1743:   PetscInt       i,j;
1744:   FunctionalLink link;
1745:   char           *names[256];

1748:   mod->numMonitored = ALEN(names);
1749:   PetscOptionsStringArray("-monitor","list of functionals to monitor","",names,&mod->numMonitored,NULL);
1750:   /* Create list of functionals that will be computed somehow */
1751:   PetscMalloc1(mod->numMonitored,&mod->functionalMonitored);
1752:   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
1753:   PetscMalloc1(mod->numMonitored,&mod->functionalCall);
1754:   mod->numCall = 0;
1755:   for (i=0; i<mod->numMonitored; i++) {
1756:     for (link=mod->functionalRegistry; link; link=link->next) {
1757:       PetscBool match;
1758:       PetscStrcasecmp(names[i],link->name,&match);
1759:       if (match) break;
1760:     }
1761:     if (!link) SETERRQ1(mod->comm,PETSC_ERR_USER,"No known functional '%s'",names[i]);
1762:     mod->functionalMonitored[i] = link;
1763:     for (j=0; j<i; j++) {
1764:       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
1765:     }
1766:     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
1767: next_name:
1768:     PetscFree(names[i]);
1769:   }

1771:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
1772:   mod->maxComputed = -1;
1773:   for (link=mod->functionalRegistry; link; link=link->next) {
1774:     for (i=0; i<mod->numCall; i++) {
1775:       FunctionalLink call = mod->functionalCall[i];
1776:       if (link->func == call->func && link->ctx == call->ctx) {
1777:         mod->maxComputed = PetscMax(mod->maxComputed,link->offset);
1778:       }
1779:     }
1780:   }
1781:   return(0);
1782: }

1786: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
1787: {
1789:   FunctionalLink l,next;

1792:   if (!link) return(0);
1793:   l     = *link;
1794:   *link = NULL;
1795:   for (; l; l=next) {
1796:     next = l->next;
1797:     PetscFree(l->name);
1798:     PetscFree(l);
1799:   }
1800:   return(0);
1801: }

1805: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
1806: {
1807:   DM                dmCell;
1808:   Model             mod = user->model;
1809:   const PetscScalar *cellgeom;
1810:   PetscScalar       *x;
1811:   PetscInt          cStart, cEnd, cEndInterior = user->cEndInterior, c;
1812:   PetscErrorCode    ierr;

1815:   VecGetDM(user->cellgeom, &dmCell);
1816:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1817:   VecGetArrayRead(user->cellgeom, &cellgeom);
1818:   VecGetArray(X, &x);
1819:   for (c = cStart; c < cEndInterior; ++c) {
1820:     const CellGeom *cg;
1821:     PetscScalar    *xc;

1823:     DMPlexPointLocalRead(dmCell,c,cellgeom,&cg);
1824:     DMPlexPointGlobalRef(dm,c,x,&xc);
1825:     if (xc) {(*mod->solution)(mod,0.0,cg->centroid,xc,mod->solutionctx);}
1826:   }
1827:   VecRestoreArrayRead(user->cellgeom, &cellgeom);
1828:   VecRestoreArray(X, &x);
1829:   return(0);
1830: }

1834: static PetscErrorCode RHSFunctionLocal_Upwind(DM dm,PetscReal time,Vec locX,Vec F,User user)
1835: {
1836:   Physics           phys = user->model->physics;
1837:   DM                dmFace, dmCell;
1838:   DMLabel           ghostLabel;
1839:   PetscErrorCode    ierr;
1840:   const PetscScalar *facegeom, *cellgeom, *x;
1841:   PetscScalar       *f;
1842:   PetscInt          fStart, fEnd, face;

1845:   VecGetDM(user->facegeom,&dmFace);
1846:   VecGetDM(user->cellgeom,&dmCell);
1847:   DMPlexInsertBoundaryValuesFVM(dm, time, locX);
1848:   VecGetArrayRead(user->facegeom,&facegeom);
1849:   VecGetArrayRead(user->cellgeom,&cellgeom);
1850:   VecGetArrayRead(locX,&x);
1851:   VecGetArray(F,&f);
1852:   DMPlexGetLabel(dm, "ghost", &ghostLabel);
1853:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1854:   for (face = fStart; face < fEnd; ++face) {
1855:     const PetscInt    *cells;
1856:     PetscInt          i,ghost;
1857:     PetscScalar       *flux = user->work.flux,*fL,*fR;
1858:     const FaceGeom    *fg;
1859:     const CellGeom    *cgL,*cgR;
1860:     const PetscScalar *xL,*xR;

1862:     DMLabelGetValue(ghostLabel, face, &ghost);
1863:     if (ghost >= 0) continue;
1864:     DMPlexGetSupport(dm, face, &cells);
1865:     DMPlexPointLocalRead(dmFace,face,facegeom,&fg);
1866:     DMPlexPointLocalRead(dmCell,cells[0],cellgeom,&cgL);
1867:     DMPlexPointLocalRead(dmCell,cells[1],cellgeom,&cgR);
1868:     DMPlexPointLocalRead(dm,cells[0],x,&xL);
1869:     DMPlexPointLocalRead(dm,cells[1],x,&xR);
1870:     DMPlexPointGlobalRef(dm,cells[0],f,&fL);
1871:     DMPlexPointGlobalRef(dm,cells[1],f,&fR);
1872:     (*phys->riemann)(phys, fg->centroid, fg->normal, xL, xR, flux);
1873:     for (i=0; i<phys->dof; i++) {
1874:       if (fL) fL[i] -= flux[i] / cgL->volume;
1875:       if (fR) fR[i] += flux[i] / cgR->volume;
1876:     }
1877:   }
1878:   VecRestoreArrayRead(user->facegeom,&facegeom);
1879:   VecRestoreArrayRead(user->cellgeom,&cellgeom);
1880:   VecRestoreArrayRead(locX,&x);
1881:   VecRestoreArray(F,&f);
1882:   return(0);
1883: }

1887: static PetscErrorCode RHSFunctionLocal_LS(DM dm,PetscReal time,Vec locX,Vec F,User user)
1888: {
1889:   DM                dmFace, dmCell;
1890:   DM                dmGrad = user->dmGrad;
1891:   Model             mod    = user->model;
1892:   Physics           phys   = mod->physics;
1893:   const PetscInt    dof    = phys->dof;
1894:   PetscErrorCode    ierr;
1895:   const PetscScalar *facegeom, *cellgeom, *x;
1896:   PetscScalar       *f;
1897:   PetscInt          fStart, fEnd, face, cStart, cell;
1898:   Vec               locGrad,Grad;

1901:   VecGetDM(user->facegeom,&dmFace);
1902:   VecGetDM(user->cellgeom,&dmCell);
1903:   DMPlexInsertBoundaryValuesFVM(dm, time, locX);
1904:   DMGetGlobalVector(dmGrad,&Grad);
1905:   VecZeroEntries(Grad);
1906:   VecGetArrayRead(user->facegeom,&facegeom);
1907:   VecGetArrayRead(user->cellgeom,&cellgeom);
1908:   VecGetArrayRead(locX,&x);
1909:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1910:   DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
1911:   {
1912:     PetscScalar *grad;
1913:     VecGetArray(Grad,&grad);
1914:     /* Reconstruct gradients */
1915:     for (face=fStart; face<fEnd; ++face) {
1916:       const PetscInt    *cells;
1917:       const PetscScalar *cx[2];
1918:       const FaceGeom    *fg;
1919:       PetscScalar       *cgrad[2];
1920:       PetscInt          i,j;
1921:       PetscBool         ghost;

1923:       IsExteriorOrGhostFace(dm,face,&ghost);
1924:       if (ghost) continue;
1925:       DMPlexGetSupport(dm,face,&cells);
1926:       DMPlexPointLocalRead(dmFace,face,facegeom,&fg);
1927:       for (i=0; i<2; i++) {
1928:         DMPlexPointLocalRead(dm,cells[i],x,&cx[i]);
1929:         DMPlexPointGlobalRef(dmGrad,cells[i],grad,&cgrad[i]);
1930:       }
1931:       for (i=0; i<dof; i++) {
1932:         PetscScalar delta = cx[1][i] - cx[0][i];
1933:         for (j=0; j<DIM; j++) {
1934:           if (cgrad[0]) cgrad[0][i*DIM+j] += fg->grad[0][j] * delta;
1935:           if (cgrad[1]) cgrad[1][i*DIM+j] -= fg->grad[1][j] * delta;
1936:         }
1937:       }
1938:     }
1939:     /* Limit interior gradients. Using cell-based loop because it generalizes better to vector limiters. */
1940:     for (cell=cStart; cell<user->cEndInterior; cell++) {
1941:       const PetscInt    *faces;
1942:       PetscInt          numFaces,f;
1943:       PetscReal         *cellPhi = user->work.state0; /* Scalar limiter applied to each component separately */
1944:       const PetscScalar *cx;
1945:       const CellGeom    *cg;
1946:       PetscScalar       *cgrad;
1947:       PetscInt          i;
1948:       DMPlexGetConeSize(dm,cell,&numFaces);
1949:       DMPlexGetCone(dm,cell,&faces);
1950:       DMPlexPointLocalRead(dm,cell,x,&cx);
1951:       DMPlexPointLocalRead(dmCell,cell,cellgeom,&cg);
1952:       DMPlexPointGlobalRef(dmGrad,cell,grad,&cgrad);
1953:       if (!cgrad) continue;     /* ghost cell, we don't compute */
1954:       /* Limiter will be minimum value over all neighbors */
1955:       for (i=0; i<dof; i++) cellPhi[i] = PETSC_MAX_REAL;
1956:       for (f=0; f<numFaces; f++) {
1957:         const PetscScalar *ncx;
1958:         const CellGeom    *ncg;
1959:         const PetscInt    *fcells;
1960:         PetscInt          face = faces[f],ncell;
1961:         PetscScalar       v[DIM];
1962:         PetscBool         ghost;
1963:         IsExteriorOrGhostFace(dm,face,&ghost);
1964:         if (ghost) continue;
1965:         DMPlexGetSupport(dm,face,&fcells);
1966:         ncell = cell == fcells[0] ? fcells[1] : fcells[0];
1967:         DMPlexPointLocalRead(dm,ncell,x,&ncx);
1968:         DMPlexPointLocalRead(dmCell,ncell,cellgeom,&ncg);
1969:         Waxpy2(-1,cg->centroid,ncg->centroid,v);
1970:         for (i=0; i<dof; i++) {
1971:           /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
1972:           PetscScalar phi,flim = 0.5 * (ncx[i] - cx[i]) / Dot2(&cgrad[i*DIM],v);
1973:           phi        = (*user->Limit)(flim);
1974:           cellPhi[i] = PetscMin(cellPhi[i],phi);
1975:         }
1976:       }
1977:       /* Apply limiter to gradient */
1978:       for (i=0; i<dof; i++) Scale2(cellPhi[i],&cgrad[i*DIM],&cgrad[i*DIM]);
1979:     }
1980:     VecRestoreArray(Grad,&grad);
1981:   }
1982:   DMGetLocalVector(dmGrad,&locGrad);
1983:   DMGlobalToLocalBegin(dmGrad,Grad,INSERT_VALUES,locGrad);
1984:   DMGlobalToLocalEnd(dmGrad,Grad,INSERT_VALUES,locGrad);
1985:   DMRestoreGlobalVector(dmGrad,&Grad);

1987:   {
1988:     DMLabel            ghostLabel, faceLabel;
1989:     const PetscScalar *grad;
1990:     DMPlexGetLabel(dm, "ghost", &ghostLabel);
1991:     DMPlexGetLabel(dm, "Face Sets", &faceLabel);
1992:     VecGetArrayRead(locGrad,&grad);
1993:     VecGetArray(F,&f);
1994:     for (face=fStart; face<fEnd; ++face) {
1995:       const PetscInt    *cells;
1996:       PetscInt          ghost,i,j,bset;
1997:       PetscScalar       *flux = user->work.flux,*fx[2] = {user->work.state0,user->work.state1},*cf[2];
1998:       const FaceGeom    *fg;
1999:       const CellGeom    *cg[2];
2000:       const PetscScalar *cx[2],*cgrad[2];

2002:       DMLabelGetValue(ghostLabel, face, &ghost);
2003:       if (ghost >= 0) continue;
2004:       DMPlexGetSupport(dm, face, &cells);
2005:       DMPlexPointLocalRead(dmFace,face,facegeom,&fg);
2006:       for (i=0; i<2; i++) {
2007:         PetscScalar dx[DIM];
2008:         DMPlexPointLocalRead(dmCell,cells[i],cellgeom,&cg[i]);
2009:         DMPlexPointLocalRead(dm,cells[i],x,&cx[i]);
2010:         DMPlexPointLocalRead(dmGrad,cells[i],grad,&cgrad[i]);
2011:         DMPlexPointGlobalRef(dm,cells[i],f,&cf[i]);
2012:         Waxpy2(-1,cg[i]->centroid,fg->centroid,dx);
2013:         for (j=0; j<dof; j++) fx[i][j] = cx[i][j] + Dot2(cgrad[i],dx);
2014:       }
2015:       DMLabelGetValue(faceLabel, face, &bset);
2016:       if (bset != -1) {
2017:         BoundaryFunction bcFunc;
2018:         void             *bcCtx;
2019:         ModelBoundaryFind(mod,bset,&bcFunc,&bcCtx);
2020:         (*bcFunc)(time,fg->centroid,fg->normal,fx[0],fx[1],bcCtx);
2021:       }
2022:       (*phys->riemann)(phys, fg->centroid, fg->normal, fx[0], fx[1], flux);
2023:       for (i=0; i<phys->dof; i++) {
2024:         if (cf[0]) cf[0][i] -= flux[i] / cg[0]->volume;
2025:         if (cf[1]) cf[1][i] += flux[i] / cg[1]->volume;
2026:       }
2027:     }
2028:     VecRestoreArrayRead(locGrad,&grad);
2029:     VecRestoreArray(F,&f);
2030:   }
2031:   VecRestoreArrayRead(user->facegeom,&facegeom);
2032:   VecRestoreArrayRead(user->cellgeom,&cellgeom);
2033:   VecRestoreArrayRead(locX,&x);
2034:   DMRestoreLocalVector(dmGrad,&locGrad);
2035:   return(0);
2036: }

2040: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
2041: {

2045:   PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);
2046:   PetscViewerSetType(*viewer, PETSCVIEWERVTK);
2047:   PetscViewerFileSetName(*viewer, filename);
2048:   return(0);
2049: }

2053: static PetscErrorCode MonitorVTK(TS ts,PetscInt stepnum,PetscReal time,Vec X,void *ctx)
2054: {
2055:   User           user = (User)ctx;
2056:   DM             dm;
2057:   PetscViewer    viewer;
2058:   char           filename[PETSC_MAX_PATH_LEN],*ftable = NULL;
2059:   PetscReal      xnorm;

2063:   PetscObjectSetName((PetscObject) X, "solution");
2064:   VecGetDM(X,&dm);
2065:   VecNorm(X,NORM_INFINITY,&xnorm);
2066:   if (stepnum >= 0) {           /* No summary for final time */
2067:     Model             mod = user->model;
2068:     PetscInt          c,cStart,cEnd,fcount,i;
2069:     size_t            ftableused,ftablealloc;
2070:     const PetscScalar *cellgeom,*x;
2071:     DM                dmCell;
2072:     PetscReal         *fmin,*fmax,*fintegral,*ftmp;
2073:     fcount = mod->maxComputed+1;
2074:     PetscMalloc4(fcount,&fmin,fcount,&fmax,fcount,&fintegral,fcount,&ftmp);
2075:     for (i=0; i<fcount; i++) {
2076:       fmin[i]      = PETSC_MAX_REAL;
2077:       fmax[i]      = PETSC_MIN_REAL;
2078:       fintegral[i] = 0;
2079:     }
2080:     DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
2081:     VecGetDM(user->cellgeom,&dmCell);
2082:     VecGetArrayRead(user->cellgeom,&cellgeom);
2083:     VecGetArrayRead(X,&x);
2084:     for (c=cStart; c<user->cEndInterior; c++) {
2085:       const CellGeom    *cg;
2086:       const PetscScalar *cx;
2087:       DMPlexPointLocalRead(dmCell,c,cellgeom,&cg);
2088:       DMPlexPointGlobalRead(dm,c,x,&cx);
2089:       if (!cx) continue;        /* not a global cell */
2090:       for (i=0; i<mod->numCall; i++) {
2091:         FunctionalLink flink = mod->functionalCall[i];
2092:         (*flink->func)(mod,time,cg->centroid,cx,ftmp,flink->ctx);
2093:       }
2094:       for (i=0; i<fcount; i++) {
2095:         fmin[i]       = PetscMin(fmin[i],ftmp[i]);
2096:         fmax[i]       = PetscMax(fmax[i],ftmp[i]);
2097:         fintegral[i] += cg->volume * ftmp[i];
2098:       }
2099:     }
2100:     VecRestoreArrayRead(user->cellgeom,&cellgeom);
2101:     VecRestoreArrayRead(X,&x);
2102:     MPI_Allreduce(MPI_IN_PLACE,fmin,fcount,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)ts));
2103:     MPI_Allreduce(MPI_IN_PLACE,fmax,fcount,MPIU_REAL,MPI_MAX,PetscObjectComm((PetscObject)ts));
2104:     MPI_Allreduce(MPI_IN_PLACE,fintegral,fcount,MPIU_REAL,MPI_SUM,PetscObjectComm((PetscObject)ts));

2106:     ftablealloc = fcount * 100;
2107:     ftableused  = 0;
2108:     PetscMalloc1(ftablealloc,&ftable);
2109:     for (i=0; i<mod->numMonitored; i++) {
2110:       size_t         countused;
2111:       char           buffer[256],*p;
2112:       FunctionalLink flink = mod->functionalMonitored[i];
2113:       PetscInt       id    = flink->offset;
2114:       if (i % 3) {
2115:         PetscMemcpy(buffer,"  ",2);
2116:         p    = buffer + 2;
2117:       } else if (i) {
2118:         char newline[] = "\n";
2119:         PetscMemcpy(buffer,newline,sizeof newline-1);
2120:         p    = buffer + sizeof newline - 1;
2121:       } else {
2122:         p = buffer;
2123:       }
2124:       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]);
2125:       countused += p - buffer;
2126:       if (countused > ftablealloc-ftableused-1) { /* reallocate */
2127:         char *ftablenew;
2128:         ftablealloc = 2*ftablealloc + countused;
2129:         PetscMalloc(ftablealloc,&ftablenew);
2130:         PetscMemcpy(ftablenew,ftable,ftableused);
2131:         PetscFree(ftable);
2132:         ftable = ftablenew;
2133:       }
2134:       PetscMemcpy(ftable+ftableused,buffer,countused);
2135:       ftableused += countused;
2136:       ftable[ftableused] = 0;
2137:     }
2138:     PetscFree4(fmin,fmax,fintegral,ftmp);

2140:     PetscPrintf(PetscObjectComm((PetscObject)ts),"% 3D  time %8.4g  |x| %8.4g  %s\n",stepnum,(double)time,(double)xnorm,ftable ? ftable : "");
2141:     PetscFree(ftable);
2142:   }
2143:   if (user->vtkInterval < 1) return(0);
2144:   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
2145:     if (stepnum == -1) {        /* Final time is not multiple of normal time interval, write it anyway */
2146:       TSGetTimeStepNumber(ts,&stepnum);
2147:     }
2148:     PetscSNPrintf(filename,sizeof filename,"ex11-%03D.vtu",stepnum);
2149:     OutputVTK(dm,filename,&viewer);
2150:     VecView(X,viewer);
2151:     PetscViewerDestroy(&viewer);
2152:   }
2153:   return(0);
2154: }

2158: int main(int argc, char **argv)
2159: {
2160:   MPI_Comm          comm;
2161:   User              user;
2162:   Model             mod;
2163:   Physics           phys;
2164:   DM                dm;
2165:   PetscReal         ftime,cfl,dt;
2166:   PetscInt          dim, nsteps;
2167:   TS                ts;
2168:   TSConvergedReason reason;
2169:   Vec               X;
2170:   PetscViewer       viewer;
2171:   PetscMPIInt       rank;
2172:   PetscBool         vtkCellGeom, splitFaces;
2173:   PetscInt          overlap;
2174:   char              filename[PETSC_MAX_PATH_LEN] = "sevenside.exo";
2175:   PetscErrorCode    ierr;

2177:   PetscInitialize(&argc, &argv, (char*) 0, help);
2178:   comm = PETSC_COMM_WORLD;
2179:   MPI_Comm_rank(comm, &rank);

2181:   PetscNew(&user);
2182:   PetscNew(&user->model);
2183:   PetscNew(&user->model->physics);
2184:   mod  = user->model;
2185:   phys = mod->physics;
2186:   mod->comm = comm;

2188:   /* Register physical models to be available on the command line */
2189:   PetscFunctionListAdd(&PhysicsList,"advect"          ,PhysicsCreate_Advect);
2190:   PetscFunctionListAdd(&PhysicsList,"sw"              ,PhysicsCreate_SW);
2191:   PetscFunctionListAdd(&PhysicsList,"euler"           ,PhysicsCreate_Euler);

2193:   PetscFunctionListAdd(&LimitList,"zero"              ,Limit_Zero);
2194:   PetscFunctionListAdd(&LimitList,"none"              ,Limit_None);
2195:   PetscFunctionListAdd(&LimitList,"minmod"            ,Limit_Minmod);
2196:   PetscFunctionListAdd(&LimitList,"vanleer"           ,Limit_VanLeer);
2197:   PetscFunctionListAdd(&LimitList,"vanalbada"         ,Limit_VanAlbada);
2198:   PetscFunctionListAdd(&LimitList,"sin"               ,Limit_Sin);
2199:   PetscFunctionListAdd(&LimitList,"superbee"          ,Limit_Superbee);
2200:   PetscFunctionListAdd(&LimitList,"mc"                ,Limit_MC);

2202:   PetscOptionsBegin(comm,NULL,"Unstructured Finite Volume Mesh Options","");
2203:   {
2204:     cfl  = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
2205:     PetscOptionsReal("-ufv_cfl","CFL number per step","",cfl,&cfl,NULL);
2206:     PetscOptionsString("-f","Exodus.II filename to read","",filename,filename,sizeof(filename),NULL);
2207:     splitFaces = PETSC_FALSE;
2208:     PetscOptionsBool("-ufv_split_faces","Split faces between cell sets","",splitFaces,&splitFaces,NULL);
2209:     overlap = 1;
2210:     PetscOptionsInt("-ufv_mesh_overlap","Number of cells to overlap partitions","",overlap,&overlap,NULL);
2211:     user->vtkInterval = 1;
2212:     PetscOptionsInt("-ufv_vtk_interval","VTK output interval (0 to disable)","",user->vtkInterval,&user->vtkInterval,NULL);
2213:     vtkCellGeom = PETSC_FALSE;
2214:     PetscOptionsBool("-ufv_vtk_cellgeom","Write cell geometry (for debugging)","",vtkCellGeom,&vtkCellGeom,NULL);
2215:   }
2216:   PetscOptionsEnd();
2217:   CreateMesh(comm, filename, overlap, splitFaces, &dm, user);
2218:   DMPlexGetDimension(dm, &dim);
2219:   PetscOptionsBegin(comm,NULL,"Unstructured Finite Volume Physics Options","");
2220:   {
2221:     PetscErrorCode (*physcreate)(DM,Model,Physics);
2222:     char             physname[256]  = "advect";
2223:     char             limitname[256] = "minmod";

2225:     PetscOptionsFList("-physics","Physics module to solve","",PhysicsList,physname,physname,sizeof physname,NULL);
2226:     PetscFunctionListFind(PhysicsList,physname,&physcreate);
2227:     PetscMemzero(phys,sizeof(struct _n_Physics));
2228:     (*physcreate)(dm,mod,phys);
2229:     mod->maxspeed = phys->maxspeed;
2230:     /* Count number of fields and dofs */
2231:     for (phys->nfields=0,phys->dof=0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;

2233:     if (mod->maxspeed <= 0) SETERRQ1(comm,PETSC_ERR_ARG_WRONGSTATE,"Physics '%s' did not set maxspeed",physname);
2234:     if (phys->dof <= 0) SETERRQ1(comm,PETSC_ERR_ARG_WRONGSTATE,"Physics '%s' did not set dof",physname);
2235:     PetscMalloc3(phys->dof,&user->work.flux,phys->dof,&user->work.state0,phys->dof,&user->work.state1);
2236:     user->reconstruct = PETSC_FALSE;
2237:     PetscOptionsBool("-ufv_reconstruct","Reconstruct gradients for a second order method (grows stencil)","",user->reconstruct,&user->reconstruct,NULL);
2238:     if (user->reconstruct) {
2239:       PetscOptionsFList("-ufv_limit","Limiter to apply to reconstructed solution","",LimitList,limitname,limitname,sizeof limitname,NULL);
2240:       PetscFunctionListFind(LimitList,limitname,&user->Limit);
2241:       ReconstructionSetUp(dm, user);
2242:     }
2243:     ModelFunctionalSetFromOptions(mod);
2244:   }
2245:   PetscOptionsEnd();

2247:   /* Set up DM with section describing local vector and configure local vector. */
2248:   SetUpLocalSpace(dm, user);
2249:   SetUpBoundaries(dm, user);

2251:   DMCreateGlobalVector(dm, &X);
2252:   PetscObjectSetName((PetscObject) X, "solution");
2253:   SetInitialCondition(dm, X, user);
2254:   if (vtkCellGeom) {
2255:     DM  dmCell;
2256:     Vec partition;

2258:     OutputVTK(dm, "ex11-cellgeom.vtk", &viewer);
2259:     VecView(user->cellgeom, viewer);
2260:     PetscViewerDestroy(&viewer);
2261:     CreatePartitionVec(dm, &dmCell, &partition);
2262:     OutputVTK(dmCell, "ex11-partition.vtk", &viewer);
2263:     VecView(partition, viewer);
2264:     PetscViewerDestroy(&viewer);
2265:     VecDestroy(&partition);
2266:     DMDestroy(&dmCell);
2267:   }

2269:   TSCreate(comm, &ts);
2270:   TSSetType(ts, TSSSP);
2271:   TSSetDM(ts, dm);
2272:   TSMonitorSet(ts,MonitorVTK,user,NULL);
2273:   DMTSSetRHSFunctionLocal(dm, (PetscErrorCode (*)(DM,PetscReal,Vec,Vec,void*)) (user->reconstruct ? RHSFunctionLocal_LS : RHSFunctionLocal_Upwind), user);
2274:   TSSetDuration(ts,1000,2.0);
2275:   dt   = cfl * user->minradius / user->model->maxspeed;
2276:   TSSetInitialTimeStep(ts,0.0,dt);
2277:   TSSetFromOptions(ts);
2278:   TSSolve(ts,X);
2279:   TSGetSolveTime(ts,&ftime);
2280:   TSGetTimeStepNumber(ts,&nsteps);
2281:   TSGetConvergedReason(ts,&reason);
2282:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,nsteps);
2283:   TSDestroy(&ts);

2285:   VecDestroy(&user->cellgeom);
2286:   VecDestroy(&user->facegeom);
2287:   DMDestroy(&user->dmGrad);
2288:   PetscFunctionListDestroy(&PhysicsList);
2289:   PetscFunctionListDestroy(&LimitList);
2290:   BoundaryLinkDestroy(&user->model->boundary);
2291:   FunctionalLinkDestroy(&user->model->functionalRegistry);
2292:   PetscFree(user->model->functionalMonitored);
2293:   PetscFree(user->model->functionalCall);
2294:   PetscFree(user->model->physics->data);
2295:   PetscFree(user->model->physics);
2296:   PetscFree(user->model);
2297:   PetscFree3(user->work.flux,user->work.state0,user->work.state1);
2298:   PetscFree(user);
2299:   VecDestroy(&X);
2300:   DMDestroy(&dm);
2301:   PetscFinalize();
2302:   return(0);
2303: }