Actual source code: ex18.c

  1: static char help[] = "Hybrid Finite Element-Finite Volume Example.\n";
  2: /*F
  3:   Here we are advecting a passive tracer in a harmonic velocity field, defined by
  4: a forcing function $f$:
  5: \begin{align}
  6:   -\Delta \mathbf{u} + f &= 0 \\
  7:   \frac{\partial\phi}{\partial t} + \nabla\cdot \phi \mathbf{u} &= 0
  8: \end{align}
  9: F*/

 11: #include <petscdmplex.h>
 12: #include <petscds.h>
 13: #include <petscts.h>

 15: #include <petsc/private/dmpleximpl.h>

 17: typedef enum {
 18:   VEL_ZERO,
 19:   VEL_CONSTANT,
 20:   VEL_HARMONIC,
 21:   VEL_SHEAR
 22: } VelocityDistribution;

 24: typedef enum {
 25:   ZERO,
 26:   CONSTANT,
 27:   GAUSSIAN,
 28:   TILTED,
 29:   DELTA
 30: } PorosityDistribution;

 32: static PetscErrorCode constant_u_2d(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *);

 34: /*
 35:   FunctionalFunc - Calculates the value of a functional of the solution at a point

 37:   Input Parameters:
 38: + dm   - The DM
 39: . time - The TS time
 40: . x    - The coordinates of the evaluation point
 41: . u    - The field values at point x
 42: - ctx  - A user context, or NULL

 44:   Output Parameter:
 45: . f    - The value of the functional at point x

 47: */
 48: typedef PetscErrorCode (*FunctionalFunc)(DM, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);

 50: typedef struct _n_Functional *Functional;
 51: struct _n_Functional {
 52:   char          *name;
 53:   FunctionalFunc func;
 54:   void          *ctx;
 55:   PetscInt       offset;
 56:   Functional     next;
 57: };

 59: typedef struct {
 60:   /* Problem definition */
 61:   PetscBool useFV; /* Use a finite volume scheme for advection */
 62:   PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
 63:   VelocityDistribution velocityDist;
 64:   PorosityDistribution porosityDist;
 65:   PetscReal            inflowState;
 66:   PetscReal            source[3];
 67:   /* Monitoring */
 68:   PetscInt    numMonitorFuncs, maxMonitorFunc;
 69:   Functional *monitorFuncs;
 70:   PetscInt    errorFunctional;
 71:   Functional  functionalRegistry;
 72: } AppCtx;

 74: static AppCtx *globalUser;

 76: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 77: {
 78:   const char *velocityDist[4] = {"zero", "constant", "harmonic", "shear"};
 79:   const char *porosityDist[5] = {"zero", "constant", "gaussian", "tilted", "delta"};
 80:   PetscInt    vd, pd, d;
 81:   PetscBool   flg;

 84:   options->useFV           = PETSC_FALSE;
 85:   options->velocityDist    = VEL_HARMONIC;
 86:   options->porosityDist    = ZERO;
 87:   options->inflowState     = -2.0;
 88:   options->numMonitorFuncs = 0;
 89:   options->source[0]       = 0.5;
 90:   options->source[1]       = 0.5;
 91:   options->source[2]       = 0.5;

 93:   PetscOptionsBegin(comm, "", "Magma Dynamics Options", "DMPLEX");
 94:   PetscOptionsBool("-use_fv", "Use the finite volume method for advection", "ex18.c", options->useFV, &options->useFV, NULL);
 95:   vd = options->velocityDist;
 96:   PetscOptionsEList("-velocity_dist", "Velocity distribution type", "ex18.c", velocityDist, 4, velocityDist[options->velocityDist], &vd, NULL);
 97:   options->velocityDist = (VelocityDistribution)vd;
 98:   pd                    = options->porosityDist;
 99:   PetscOptionsEList("-porosity_dist", "Initial porosity distribution type", "ex18.c", porosityDist, 5, porosityDist[options->porosityDist], &pd, NULL);
100:   options->porosityDist = (PorosityDistribution)pd;
101:   PetscOptionsReal("-inflow_state", "The inflow state", "ex18.c", options->inflowState, &options->inflowState, NULL);
102:   d = 2;
103:   PetscOptionsRealArray("-source_loc", "The source location", "ex18.c", options->source, &d, &flg);
105:   PetscOptionsEnd();

107:   return 0;
108: }

110: static PetscErrorCode ProcessMonitorOptions(MPI_Comm comm, AppCtx *options)
111: {
112:   Functional func;
113:   char      *names[256];
114:   PetscInt   f;

117:   PetscOptionsBegin(comm, "", "Simulation Monitor Options", "DMPLEX");
118:   options->numMonitorFuncs = PETSC_STATIC_ARRAY_LENGTH(names);
119:   PetscOptionsStringArray("-monitor", "List of functionals to monitor", "", names, &options->numMonitorFuncs, NULL);
120:   PetscMalloc1(options->numMonitorFuncs, &options->monitorFuncs);
121:   for (f = 0; f < options->numMonitorFuncs; ++f) {
122:     for (func = options->functionalRegistry; func; func = func->next) {
123:       PetscBool match;

125:       PetscStrcasecmp(names[f], func->name, &match);
126:       if (match) break;
127:     }
129:     options->monitorFuncs[f] = func;
130:     /* Jed inserts a de-duplication of functionals here */
131:     PetscFree(names[f]);
132:   }
133:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
134:   options->maxMonitorFunc = -1;
135:   for (func = options->functionalRegistry; func; func = func->next) {
136:     for (f = 0; f < options->numMonitorFuncs; ++f) {
137:       Functional call = options->monitorFuncs[f];

139:       if (func->func == call->func && func->ctx == call->ctx) options->maxMonitorFunc = PetscMax(options->maxMonitorFunc, func->offset);
140:     }
141:   }
142:   PetscOptionsEnd();
143:   return 0;
144: }

146: static PetscErrorCode FunctionalRegister(Functional *functionalRegistry, const char name[], PetscInt *offset, FunctionalFunc func, void *ctx)
147: {
148:   Functional *ptr, f;
149:   PetscInt    lastoffset = -1;

152:   for (ptr = functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
153:   PetscNew(&f);
154:   PetscStrallocpy(name, &f->name);
155:   f->offset = lastoffset + 1;
156:   f->func   = func;
157:   f->ctx    = ctx;
158:   f->next   = NULL;
159:   *ptr      = f;
160:   *offset   = f->offset;
161:   return 0;
162: }

164: static PetscErrorCode FunctionalDestroy(Functional *link)
165: {
166:   Functional next, l;

169:   if (!link) return 0;
170:   l     = *link;
171:   *link = NULL;
172:   for (; l; l = next) {
173:     next = l->next;
174:     PetscFree(l->name);
175:     PetscFree(l);
176:   }
177:   return 0;
178: }

180: static void f0_zero_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
181: {
182:   PetscInt comp;
183:   for (comp = 0; comp < dim; ++comp) f0[comp] = u[comp];
184: }

186: static void f0_constant_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
187: {
188:   PetscScalar wind[3] = {0.0, 0.0, 0.0};
189:   PetscInt    comp;

191:   constant_u_2d(dim, t, x, Nf, wind, NULL);
192:   for (comp = 0; comp < dim && comp < 3; ++comp) f0[comp] = u[comp] - wind[comp];
193: }

195: static void f1_constant_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
196: {
197:   PetscInt comp;
198:   for (comp = 0; comp < dim * dim; ++comp) f1[comp] = 0.0;
199: }

201: static void g0_constant_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
202: {
203:   PetscInt d;
204:   for (d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
205: }

207: static void g0_constant_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
208: {
209:   g0[0] = 1.0;
210: }

212: static void f0_lap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
213: {
214:   PetscInt comp;
215:   for (comp = 0; comp < dim; ++comp) f0[comp] = 4.0;
216: }

218: static void f1_lap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
219: {
220:   PetscInt comp, d;
221:   for (comp = 0; comp < dim; ++comp) {
222:     for (d = 0; d < dim; ++d) f1[comp * dim + d] = u_x[comp * dim + d];
223:   }
224: }

226: static void f0_lap_periodic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
227: {
228:   f0[0] = -PetscSinReal(2.0 * PETSC_PI * x[0]);
229:   f0[1] = 2.0 * PETSC_PI * x[1] * PetscCosReal(2.0 * PETSC_PI * x[0]);
230: }

232: static void f0_lap_doubly_periodic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
233: {
234:   f0[0] = -2.0 * PetscSinReal(2.0 * PETSC_PI * x[0]) * PetscCosReal(2.0 * PETSC_PI * x[1]);
235:   f0[1] = 2.0 * PetscSinReal(2.0 * PETSC_PI * x[1]) * PetscCosReal(2.0 * PETSC_PI * x[0]);
236: }

238: void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
239: {
240:   const PetscInt Ncomp = dim;
241:   PetscInt       compI, d;

243:   for (compI = 0; compI < Ncomp; ++compI) {
244:     for (d = 0; d < dim; ++d) g3[((compI * Ncomp + compI) * dim + d) * dim + d] = 1.0;
245:   }
246: }

248: /* \frac{\partial\phi}{\partial t} + \nabla\phi \cdot \mathbf{u} + \phi \nabla \cdot \mathbf{u} = 0 */
249: static void f0_advection(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
250: {
251:   PetscInt d;
252:   f0[0] = u_t[dim];
253:   for (d = 0; d < dim; ++d) f0[0] += u[dim] * u_x[d * dim + d] + u_x[dim * dim + d] * u[d];
254: }

256: static void f1_advection(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
257: {
258:   PetscInt d;
259:   for (d = 0; d < dim; ++d) f1[0] = 0.0;
260: }

262: void g0_adv_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
263: {
264:   PetscInt d;
265:   g0[0] = u_tShift;
266:   for (d = 0; d < dim; ++d) g0[0] += u_x[d * dim + d];
267: }

269: void g1_adv_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
270: {
271:   PetscInt d;
272:   for (d = 0; d < dim; ++d) g1[d] = u[d];
273: }

275: void g0_adv_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
276: {
277:   PetscInt d;
278:   for (d = 0; d < dim; ++d) g0[0] += u_x[dim * dim + d];
279: }

281: void g1_adv_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
282: {
283:   PetscInt d;
284:   for (d = 0; d < dim; ++d) g1[d * dim + d] = u[dim];
285: }

287: static void riemann_advection(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, void *ctx)
288: {
289:   PetscReal wind[3] = {0.0, 1.0, 0.0};
290:   PetscReal wn      = DMPlex_DotRealD_Internal(PetscMin(dim, 3), wind, n);

292:   flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn;
293: }

295: static void riemann_coupled_advection(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, void *ctx)
296: {
297:   PetscReal wn = DMPlex_DotD_Internal(dim, uL, n);

299: #if 1
300:   flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn;
301: #else
302:   /* if (fabs(uL[0] - wind[0]) > 1.0e-7 || fabs(uL[1] - wind[1]) > 1.0e-7) PetscPrintf(PETSC_COMM_SELF, "wind (%g, %g) uL (%g, %g) uR (%g, %g)\n", wind[0], wind[1], uL[0], uL[1], uR[0], uR[1]); */
303:   /* Smear it out */
304:   flux[0] = 0.5 * ((uL[dim] + uR[dim]) + (uL[dim] - uR[dim]) * tanh(1.0e5 * wn)) * wn;
305: #endif
306: }

308: static PetscErrorCode zero_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
309: {
310:   u[0] = 0.0;
311:   u[1] = 0.0;
312:   return 0;
313: }

315: static PetscErrorCode constant_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
316: {
317:   u[0] = 0.0;
318:   u[1] = 1.0;
319:   return 0;
320: }

322: /* Coordinates of the point which was at x at t = 0 */
323: static PetscErrorCode constant_x_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
324: {
325:   const PetscReal t = *((PetscReal *)ctx);
326:   u[0]              = x[0];
327:   u[1]              = x[1] + t;
328: #if 0
329:   DMLocalizeCoordinate(globalUser->dm, u, PETSC_FALSE, u);
330: #else
331:   u[1]    = u[1] - (int)PetscRealPart(u[1]);
332: #endif
333:   return 0;
334: }

336: /*
337:   In 2D we use the exact solution:

339:     u   = x^2 + y^2
340:     v   = 2 x^2 - 2xy
341:     phi = h(x + y + (u + v) t)
342:     f_x = f_y = 4

344:   so that

346:     -\Delta u + f = <-4, -4> + <4, 4> = 0
347:     {\partial\phi}{\partial t} - \nabla\cdot \phi u = 0
348:     h_t(x + y + (u + v) t) - u . grad phi - phi div u
349:   = u h' + v h'              - u h_x - v h_y
350:   = 0

352: We will conserve phi since

354:     \nabla \cdot u = 2x - 2x = 0

356: Also try h((x + ut)^2 + (y + vt)^2), so that

358:     h_t((x + ut)^2 + (y + vt)^2) - u . grad phi - phi div u
359:   = 2 h' (u (x + ut) + v (y + vt)) - u h_x - v h_y
360:   = 2 h' (u (x + ut) + v (y + vt)) - u h' 2 (x + u t) - v h' 2 (y + vt)
361:   = 2 h' (u (x + ut) + v (y + vt)  - u (x + u t) - v (y + vt))
362:   = 0

364: */
365: static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
366: {
367:   u[0] = x[0] * x[0] + x[1] * x[1];
368:   u[1] = 2.0 * x[0] * x[0] - 2.0 * x[0] * x[1];
369:   return 0;
370: }

372: /*
373:   In 2D we use the exact, periodic solution:

375:     u   =  sin(2 pi x)/4 pi^2
376:     v   = -y cos(2 pi x)/2 pi
377:     phi = h(x + y + (u + v) t)
378:     f_x = -sin(2 pi x)
379:     f_y = 2 pi y cos(2 pi x)

381:   so that

383:     -\Delta u + f = <sin(2pi x),  -2pi y cos(2pi x)> + <-sin(2pi x), 2pi y cos(2pi x)> = 0

385: We will conserve phi since

387:     \nabla \cdot u = cos(2pi x)/2pi - cos(2pi x)/2pi = 0
388: */
389: static PetscErrorCode periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
390: {
391:   u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]) / PetscSqr(2.0 * PETSC_PI);
392:   u[1] = -x[1] * PetscCosReal(2.0 * PETSC_PI * x[0]) / (2.0 * PETSC_PI);
393:   return 0;
394: }

396: /*
397:   In 2D we use the exact, doubly periodic solution:

399:     u   =  sin(2 pi x) cos(2 pi y)/4 pi^2
400:     v   = -sin(2 pi y) cos(2 pi x)/4 pi^2
401:     phi = h(x + y + (u + v) t)
402:     f_x = -2sin(2 pi x) cos(2 pi y)
403:     f_y =  2sin(2 pi y) cos(2 pi x)

405:   so that

407:     -\Delta u + f = <2 sin(2pi x) cos(2pi y),  -2 sin(2pi y) cos(2pi x)> + <-2 sin(2pi x) cos(2pi y), 2 sin(2pi y) cos(2pi x)> = 0

409: We will conserve phi since

411:     \nabla \cdot u = cos(2pi x) cos(2pi y)/2pi - cos(2pi y) cos(2pi x)/2pi = 0
412: */
413: static PetscErrorCode doubly_periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
414: {
415:   u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]) * PetscCosReal(2.0 * PETSC_PI * x[1]) / PetscSqr(2.0 * PETSC_PI);
416:   u[1] = -PetscSinReal(2.0 * PETSC_PI * x[1]) * PetscCosReal(2.0 * PETSC_PI * x[0]) / PetscSqr(2.0 * PETSC_PI);
417:   return 0;
418: }

420: static PetscErrorCode shear_bc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
421: {
422:   u[0] = x[1] - 0.5;
423:   u[1] = 0.0;
424:   return 0;
425: }

427: static PetscErrorCode initialVelocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
428: {
429:   PetscInt d;
430:   for (d = 0; d < dim; ++d) u[d] = 0.0;
431:   return 0;
432: }

434: static PetscErrorCode zero_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
435: {
436:   u[0] = 0.0;
437:   return 0;
438: }

440: static PetscErrorCode constant_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
441: {
442:   u[0] = 1.0;
443:   return 0;
444: }

446: static PetscErrorCode delta_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
447: {
448:   PetscReal   x0[2];
449:   PetscScalar xn[2];

451:   x0[0] = globalUser->source[0];
452:   x0[1] = globalUser->source[1];
453:   constant_x_2d(dim, time, x0, Nf, xn, ctx);
454:   {
455:     const PetscReal xi  = x[0] - PetscRealPart(xn[0]);
456:     const PetscReal eta = x[1] - PetscRealPart(xn[1]);
457:     const PetscReal r2  = xi * xi + eta * eta;

459:     u[0] = r2 < 1.0e-7 ? 1.0 : 0.0;
460:   }
461:   return 0;
462: }

464: /*
465:   Gaussian blob, initially centered on (0.5, 0.5)

467:   xi = x(t) - x0, eta = y(t) - y0

469: where x(t), y(t) are the integral curves of v(t),

471:   dx/dt . grad f = v . f

473: Check: constant v(t) = {v0, w0}, x(t) = {x0 + v0 t, y0 + w0 t}

475:   v0 f_x + w0 f_y = v . f
476: */
477: static PetscErrorCode gaussian_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
478: {
479:   const PetscReal x0[2] = {0.5, 0.5};
480:   const PetscReal sigma = 1.0 / 6.0;
481:   PetscScalar     xn[2];

483:   constant_x_2d(dim, time, x0, Nf, xn, ctx);
484:   {
485:     /* const PetscReal xi  = x[0] + (sin(2.0*PETSC_PI*x[0])/(4.0*PETSC_PI*PETSC_PI))*t - x0[0]; */
486:     /* const PetscReal eta = x[1] + (-x[1]*cos(2.0*PETSC_PI*x[0])/(2.0*PETSC_PI))*t - x0[1]; */
487:     const PetscReal xi  = x[0] - PetscRealPart(xn[0]);
488:     const PetscReal eta = x[1] - PetscRealPart(xn[1]);
489:     const PetscReal r2  = xi * xi + eta * eta;

491:     u[0] = PetscExpReal(-r2 / (2.0 * sigma * sigma)) / (sigma * PetscSqrtReal(2.0 * PETSC_PI));
492:   }
493:   return 0;
494: }

496: static PetscErrorCode tilted_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
497: {
498:   PetscReal       x0[3];
499:   const PetscReal wind[3] = {0.0, 1.0, 0.0};
500:   const PetscReal t       = *((PetscReal *)ctx);

502:   DMPlex_WaxpyD_Internal(2, -t, wind, x, x0);
503:   if (x0[1] > 0) u[0] = 1.0 * x[0] + 3.0 * x[1];
504:   else u[0] = -2.0; /* Inflow state */
505:   return 0;
506: }

508: static PetscErrorCode tilted_phi_coupled_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
509: {
510:   PetscReal       ur[3];
511:   PetscReal       x0[3];
512:   const PetscReal t = *((PetscReal *)ctx);

514:   ur[0] = PetscRealPart(u[0]);
515:   ur[1] = PetscRealPart(u[1]);
516:   ur[2] = PetscRealPart(u[2]);
517:   DMPlex_WaxpyD_Internal(2, -t, ur, x, x0);
518:   if (x0[1] > 0) u[0] = 1.0 * x[0] + 3.0 * x[1];
519:   else u[0] = -2.0; /* Inflow state */
520:   return 0;
521: }

523: static PetscErrorCode advect_inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
524: {
525:   AppCtx *user = (AppCtx *)ctx;

528:   xG[0] = user->inflowState;
529:   return 0;
530: }

532: static PetscErrorCode advect_outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
533: {
535:   //xG[0] = xI[dim];
536:   xG[0] = xI[2];
537:   return 0;
538: }

540: static PetscErrorCode ExactSolution(DM dm, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
541: {
542:   AppCtx  *user = (AppCtx *)ctx;
543:   PetscInt dim;

546:   DMGetDimension(dm, &dim);
547:   switch (user->porosityDist) {
548:   case TILTED:
549:     if (user->velocityDist == VEL_ZERO) tilted_phi_2d(dim, time, x, 2, u, (void *)&time);
550:     else tilted_phi_coupled_2d(dim, time, x, 2, u, (void *)&time);
551:     break;
552:   case GAUSSIAN:
553:     gaussian_phi_2d(dim, time, x, 2, u, (void *)&time);
554:     break;
555:   case DELTA:
556:     delta_phi_2d(dim, time, x, 2, u, (void *)&time);
557:     break;
558:   default:
559:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
560:   }
561:   return 0;
562: }

564: static PetscErrorCode Functional_Error(DM dm, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx)
565: {
566:   AppCtx     *user      = (AppCtx *)ctx;
567:   PetscScalar yexact[3] = {0, 0, 0};

570:   ExactSolution(dm, time, x, yexact, ctx);
571:   f[user->errorFunctional] = PetscAbsScalar(y[0] - yexact[0]);
572:   return 0;
573: }

575: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
576: {
578:   DMCreate(comm, dm);
579:   DMSetType(*dm, DMPLEX);
580:   DMSetFromOptions(*dm);
581:   DMViewFromOptions(*dm, NULL, "-orig_dm_view");
582:   return 0;
583: }

585: static PetscErrorCode SetupBC(DM dm, AppCtx *user)
586: {
587:   PetscErrorCode (*exactFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
588:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
589:   PetscDS        prob;
590:   DMLabel        label;
591:   PetscBool      check;
592:   PetscInt       dim, n = 3;
593:   const char    *prefix;

596:   PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
597:   PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *)bdt, &n, NULL);
598:   DMGetDimension(dm, &dim);
599:   /* Set initial guesses and exact solutions */
600:   switch (dim) {
601:   case 2:
602:     user->initialGuess[0] = initialVelocity;
603:     switch (user->porosityDist) {
604:     case ZERO:
605:       user->initialGuess[1] = zero_phi;
606:       break;
607:     case CONSTANT:
608:       user->initialGuess[1] = constant_phi;
609:       break;
610:     case GAUSSIAN:
611:       user->initialGuess[1] = gaussian_phi_2d;
612:       break;
613:     case DELTA:
614:       user->initialGuess[1] = delta_phi_2d;
615:       break;
616:     case TILTED:
617:       if (user->velocityDist == VEL_ZERO) user->initialGuess[1] = tilted_phi_2d;
618:       else user->initialGuess[1] = tilted_phi_coupled_2d;
619:       break;
620:     }
621:     break;
622:   default:
623:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", dim);
624:   }
625:   exactFuncs[0] = user->initialGuess[0];
626:   exactFuncs[1] = user->initialGuess[1];
627:   switch (dim) {
628:   case 2:
629:     switch (user->velocityDist) {
630:     case VEL_ZERO:
631:       exactFuncs[0] = zero_u_2d;
632:       break;
633:     case VEL_CONSTANT:
634:       exactFuncs[0] = constant_u_2d;
635:       break;
636:     case VEL_HARMONIC:
637:       switch (bdt[0]) {
638:       case DM_BOUNDARY_PERIODIC:
639:         switch (bdt[1]) {
640:         case DM_BOUNDARY_PERIODIC:
641:           exactFuncs[0] = doubly_periodic_u_2d;
642:           break;
643:         default:
644:           exactFuncs[0] = periodic_u_2d;
645:           break;
646:         }
647:         break;
648:       default:
649:         exactFuncs[0] = quadratic_u_2d;
650:         break;
651:       }
652:       break;
653:     case VEL_SHEAR:
654:       exactFuncs[0] = shear_bc;
655:       break;
656:     default:
657:       SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %" PetscInt_FMT, dim);
658:     }
659:     break;
660:   default:
661:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", dim);
662:   }
663:   {
664:     PetscBool isImplicit = PETSC_FALSE;

666:     PetscOptionsHasName(NULL, "", "-use_implicit", &isImplicit);
667:     if (user->velocityDist == VEL_CONSTANT && !isImplicit) user->initialGuess[0] = exactFuncs[0];
668:   }
669:   PetscOptionsHasName(NULL, NULL, "-dmts_check", &check);
670:   if (check) {
671:     user->initialGuess[0] = exactFuncs[0];
672:     user->initialGuess[1] = exactFuncs[1];
673:   }
674:   /* Set BC */
675:   DMGetDS(dm, &prob);
676:   DMGetLabel(dm, "marker", &label);
677:   PetscDSSetExactSolution(prob, 0, exactFuncs[0], user);
678:   PetscDSSetExactSolution(prob, 1, exactFuncs[1], user);
679:   if (label) {
680:     const PetscInt id = 1;

682:     DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], NULL, user, NULL);
683:   }
684:   DMGetLabel(dm, "Face Sets", &label);
685:   if (label && user->useFV) {
686:     const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};

688:     DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 1, 0, NULL, (void (*)(void))advect_inflow, NULL, user, NULL);
689:     DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 1, 0, NULL, (void (*)(void))advect_outflow, NULL, user, NULL);
690:   }
691:   return 0;
692: }

694: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
695: {
696:   DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
697:   PetscDS        prob;
698:   PetscInt       n = 3;
699:   const char    *prefix;

702:   PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
703:   PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *)bdt, &n, NULL);
704:   DMGetDS(dm, &prob);
705:   switch (user->velocityDist) {
706:   case VEL_ZERO:
707:     PetscDSSetResidual(prob, 0, f0_zero_u, f1_constant_u);
708:     break;
709:   case VEL_CONSTANT:
710:     PetscDSSetResidual(prob, 0, f0_constant_u, f1_constant_u);
711:     PetscDSSetJacobian(prob, 0, 0, g0_constant_uu, NULL, NULL, NULL);
712:     PetscDSSetJacobian(prob, 1, 1, g0_constant_pp, NULL, NULL, NULL);
713:     break;
714:   case VEL_HARMONIC:
715:     switch (bdt[0]) {
716:     case DM_BOUNDARY_PERIODIC:
717:       switch (bdt[1]) {
718:       case DM_BOUNDARY_PERIODIC:
719:         PetscDSSetResidual(prob, 0, f0_lap_doubly_periodic_u, f1_lap_u);
720:         break;
721:       default:
722:         PetscDSSetResidual(prob, 0, f0_lap_periodic_u, f1_lap_u);
723:         break;
724:       }
725:       break;
726:     default:
727:       PetscDSSetResidual(prob, 0, f0_lap_u, f1_lap_u);
728:       break;
729:     }
730:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
731:     break;
732:   case VEL_SHEAR:
733:     PetscDSSetResidual(prob, 0, f0_zero_u, f1_lap_u);
734:     PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
735:     break;
736:   }
737:   PetscDSSetResidual(prob, 1, f0_advection, f1_advection);
738:   PetscDSSetJacobian(prob, 1, 1, g0_adv_pp, g1_adv_pp, NULL, NULL);
739:   PetscDSSetJacobian(prob, 1, 0, g0_adv_pu, g1_adv_pu, NULL, NULL);
740:   if (user->velocityDist == VEL_ZERO) PetscDSSetRiemannSolver(prob, 1, riemann_advection);
741:   else PetscDSSetRiemannSolver(prob, 1, riemann_coupled_advection);

743:   FunctionalRegister(&user->functionalRegistry, "Error", &user->errorFunctional, Functional_Error, user);
744:   return 0;
745: }

747: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
748: {
749:   DM              cdm = dm;
750:   PetscQuadrature q;
751:   PetscFE         fe[2];
752:   PetscFV         fv;
753:   MPI_Comm        comm;
754:   PetscInt        dim;

757:   /* Create finite element */
758:   PetscObjectGetComm((PetscObject)dm, &comm);
759:   DMGetDimension(dm, &dim);
760:   PetscFECreateDefault(comm, dim, dim, PETSC_FALSE, "velocity_", PETSC_DEFAULT, &fe[0]);
761:   PetscObjectSetName((PetscObject)fe[0], "velocity");
762:   PetscFECreateDefault(comm, dim, 1, PETSC_FALSE, "porosity_", PETSC_DEFAULT, &fe[1]);
763:   PetscFECopyQuadrature(fe[0], fe[1]);
764:   PetscObjectSetName((PetscObject)fe[1], "porosity");

766:   PetscFVCreate(PetscObjectComm((PetscObject)dm), &fv);
767:   PetscObjectSetName((PetscObject)fv, "porosity");
768:   PetscFVSetFromOptions(fv);
769:   PetscFVSetNumComponents(fv, 1);
770:   PetscFVSetSpatialDimension(fv, dim);
771:   PetscFEGetQuadrature(fe[0], &q);
772:   PetscFVSetQuadrature(fv, q);

774:   DMSetField(dm, 0, NULL, (PetscObject)fe[0]);
775:   if (user->useFV) DMSetField(dm, 1, NULL, (PetscObject)fv);
776:   else DMSetField(dm, 1, NULL, (PetscObject)fe[1]);
777:   DMCreateDS(dm);
778:   SetupProblem(dm, user);

780:   /* Set discretization and boundary conditions for each mesh */
781:   while (cdm) {
782:     DMCopyDisc(dm, cdm);
783:     DMGetCoarseDM(cdm, &cdm);
784:     /* Coordinates were never localized for coarse meshes */
785:     if (cdm) DMLocalizeCoordinates(cdm);
786:   }
787:   PetscFEDestroy(&fe[0]);
788:   PetscFEDestroy(&fe[1]);
789:   PetscFVDestroy(&fv);
790:   return 0;
791: }

793: static PetscErrorCode CreateDM(MPI_Comm comm, AppCtx *user, DM *dm)
794: {
796:   CreateMesh(comm, user, dm);
797:   /* Handle refinement, etc. */
798:   DMSetFromOptions(*dm);
799:   /* Construct ghost cells */
800:   if (user->useFV) {
801:     DM gdm;

803:     DMPlexConstructGhostCells(*dm, NULL, NULL, &gdm);
804:     DMDestroy(dm);
805:     *dm = gdm;
806:   }
807:   /* Localize coordinates */
808:   DMLocalizeCoordinates(*dm);
809:   PetscObjectSetName((PetscObject)(*dm), "Mesh");
810:   DMViewFromOptions(*dm, NULL, "-dm_view");
811:   /* Setup problem */
812:   SetupDiscretization(*dm, user);
813:   /* Setup BC */
814:   SetupBC(*dm, user);
815:   return 0;
816: }

818: static PetscErrorCode SetInitialConditionFVM(DM dm, Vec X, PetscInt field, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx)
819: {
820:   PetscDS            prob;
821:   DM                 dmCell;
822:   Vec                cellgeom;
823:   const PetscScalar *cgeom;
824:   PetscScalar       *x;
825:   PetscInt           dim, Nf, cStart, cEnd, c;

828:   DMGetDS(dm, &prob);
829:   DMGetDimension(dm, &dim);
830:   PetscDSGetNumFields(prob, &Nf);
831:   DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL);
832:   VecGetDM(cellgeom, &dmCell);
833:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
834:   VecGetArrayRead(cellgeom, &cgeom);
835:   VecGetArray(X, &x);
836:   for (c = cStart; c < cEnd; ++c) {
837:     PetscFVCellGeom *cg;
838:     PetscScalar     *xc;

840:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
841:     DMPlexPointGlobalFieldRef(dm, c, field, x, &xc);
842:     if (xc) (*func)(dim, 0.0, cg->centroid, Nf, xc, ctx);
843:   }
844:   VecRestoreArrayRead(cellgeom, &cgeom);
845:   VecRestoreArray(X, &x);
846:   return 0;
847: }

849: static PetscErrorCode MonitorFunctionals(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
850: {
851:   AppCtx            *user   = (AppCtx *)ctx;
852:   char              *ftable = NULL;
853:   DM                 dm;
854:   PetscSection       s;
855:   Vec                cellgeom;
856:   const PetscScalar *x;
857:   PetscScalar       *a;
858:   PetscReal         *xnorms;
859:   PetscInt           pStart, pEnd, p, Nf, f;

862:   VecViewFromOptions(X, (PetscObject)ts, "-view_solution");
863:   VecGetDM(X, &dm);
864:   DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL);
865:   DMGetLocalSection(dm, &s);
866:   PetscSectionGetNumFields(s, &Nf);
867:   PetscSectionGetChart(s, &pStart, &pEnd);
868:   PetscCalloc1(Nf * 2, &xnorms);
869:   VecGetArrayRead(X, &x);
870:   for (p = pStart; p < pEnd; ++p) {
871:     for (f = 0; f < Nf; ++f) {
872:       PetscInt dof, cdof, d;

874:       PetscSectionGetFieldDof(s, p, f, &dof);
875:       PetscSectionGetFieldConstraintDof(s, p, f, &cdof);
876:       DMPlexPointGlobalFieldRead(dm, p, f, x, &a);
877:       /* TODO Use constrained indices here */
878:       for (d = 0; d < dof - cdof; ++d) xnorms[f * 2 + 0] = PetscMax(xnorms[f * 2 + 0], PetscAbsScalar(a[d]));
879:       for (d = 0; d < dof - cdof; ++d) xnorms[f * 2 + 1] += PetscAbsScalar(a[d]);
880:     }
881:   }
882:   VecRestoreArrayRead(X, &x);
883:   if (stepnum >= 0) { /* No summary for final time */
884:     DM                 dmCell, *fdm;
885:     Vec               *fv;
886:     const PetscScalar *cgeom;
887:     PetscScalar      **fx;
888:     PetscReal         *fmin, *fmax, *fint, *ftmp, t;
889:     PetscInt           cStart, cEnd, c, fcount, f, num;

891:     size_t ftableused, ftablealloc;

893:     /* Functionals have indices after registering, this is an upper bound */
894:     fcount = user->numMonitorFuncs;
895:     PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fint, fcount, &ftmp);
896:     PetscMalloc3(fcount, &fdm, fcount, &fv, fcount, &fx);
897:     for (f = 0; f < fcount; ++f) {
898:       PetscSection fs;
899:       const char  *name = user->monitorFuncs[f]->name;

901:       fmin[f] = PETSC_MAX_REAL;
902:       fmax[f] = PETSC_MIN_REAL;
903:       fint[f] = 0;
904:       /* Make monitor vecs */
905:       DMClone(dm, &fdm[f]);
906:       DMGetOutputSequenceNumber(dm, &num, &t);
907:       DMSetOutputSequenceNumber(fdm[f], num, t);
908:       PetscSectionClone(s, &fs);
909:       PetscSectionSetFieldName(fs, 0, NULL);
910:       PetscSectionSetFieldName(fs, 1, name);
911:       DMSetLocalSection(fdm[f], fs);
912:       PetscSectionDestroy(&fs);
913:       DMGetGlobalVector(fdm[f], &fv[f]);
914:       PetscObjectSetName((PetscObject)fv[f], name);
915:       VecGetArray(fv[f], &fx[f]);
916:     }
917:     DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
918:     VecGetDM(cellgeom, &dmCell);
919:     VecGetArrayRead(cellgeom, &cgeom);
920:     VecGetArrayRead(X, &x);
921:     for (c = cStart; c < cEnd; ++c) {
922:       PetscFVCellGeom *cg;
923:       PetscScalar     *cx;

925:       DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
926:       DMPlexPointGlobalFieldRead(dm, c, 1, x, &cx);
927:       if (!cx) continue; /* not a global cell */
928:       for (f = 0; f < user->numMonitorFuncs; ++f) {
929:         Functional   func = user->monitorFuncs[f];
930:         PetscScalar *fxc;

932:         DMPlexPointGlobalFieldRef(dm, c, 1, fx[f], &fxc);
933:         /* I need to make it easier to get interpolated values here */
934:         (*func->func)(dm, time, cg->centroid, cx, ftmp, func->ctx);
935:         fxc[0] = ftmp[user->monitorFuncs[f]->offset];
936:       }
937:       for (f = 0; f < fcount; ++f) {
938:         fmin[f] = PetscMin(fmin[f], ftmp[f]);
939:         fmax[f] = PetscMax(fmax[f], ftmp[f]);
940:         fint[f] += cg->volume * ftmp[f];
941:       }
942:     }
943:     VecRestoreArrayRead(cellgeom, &cgeom);
944:     VecRestoreArrayRead(X, &x);
945:     MPI_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts));
946:     MPI_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
947:     MPI_Allreduce(MPI_IN_PLACE, fint, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));
948:     /* Output functional data */
949:     ftablealloc = fcount * 100;
950:     ftableused  = 0;
951:     PetscCalloc1(ftablealloc, &ftable);
952:     for (f = 0; f < user->numMonitorFuncs; ++f) {
953:       Functional func      = user->monitorFuncs[f];
954:       PetscInt   id        = func->offset;
955:       char       newline[] = "\n";
956:       char       buffer[256], *p, *prefix;
957:       size_t     countused, len;

959:       /* Create string with functional outputs */
960:       if (f % 3) {
961:         PetscArraycpy(buffer, "  ", 2);
962:         p = buffer + 2;
963:       } else if (f) {
964:         PetscArraycpy(buffer, newline, sizeof(newline) - 1);
965:         p = buffer + sizeof(newline) - 1;
966:       } else {
967:         p = buffer;
968:       }
969:       PetscSNPrintfCount(p, sizeof buffer - (p - buffer), "%12s [%12.6g,%12.6g] int %12.6g", &countused, func->name, (double)fmin[id], (double)fmax[id], (double)fint[id]);
970:       countused += p - buffer;
971:       /* reallocate */
972:       if (countused > ftablealloc - ftableused - 1) {
973:         char *ftablenew;

975:         ftablealloc = 2 * ftablealloc + countused;
976:         PetscMalloc1(ftablealloc, &ftablenew);
977:         PetscArraycpy(ftablenew, ftable, ftableused);
978:         PetscFree(ftable);
979:         ftable = ftablenew;
980:       }
981:       PetscArraycpy(ftable + ftableused, buffer, countused);
982:       ftableused += countused;
983:       ftable[ftableused] = 0;
984:       /* Output vecs */
985:       VecRestoreArray(fv[f], &fx[f]);
986:       PetscStrlen(func->name, &len);
987:       PetscMalloc1(len + 2, &prefix);
988:       PetscStrcpy(prefix, func->name);
989:       PetscStrcat(prefix, "_");
990:       PetscObjectSetOptionsPrefix((PetscObject)fv[f], prefix);
991:       VecViewFromOptions(fv[f], NULL, "-vec_view");
992:       PetscFree(prefix);
993:       DMRestoreGlobalVector(fdm[f], &fv[f]);
994:       DMDestroy(&fdm[f]);
995:     }
996:     PetscFree4(fmin, fmax, fint, ftmp);
997:     PetscFree3(fdm, fv, fx);
998:     PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| (", stepnum, (double)time);
999:     for (f = 0; f < Nf; ++f) {
1000:       if (f > 0) PetscPrintf(PetscObjectComm((PetscObject)ts), ", ");
1001:       PetscPrintf(PetscObjectComm((PetscObject)ts), "%8.4g", (double)xnorms[f * 2 + 0]);
1002:     }
1003:     PetscPrintf(PetscObjectComm((PetscObject)ts), ") |x|_1 (");
1004:     for (f = 0; f < Nf; ++f) {
1005:       if (f > 0) PetscPrintf(PetscObjectComm((PetscObject)ts), ", ");
1006:       PetscPrintf(PetscObjectComm((PetscObject)ts), "%8.4g", (double)xnorms[f * 2 + 1]);
1007:     }
1008:     PetscPrintf(PetscObjectComm((PetscObject)ts), ")  %s\n", ftable ? ftable : "");
1009:     PetscFree(ftable);
1010:   }
1011:   PetscFree(xnorms);
1012:   return 0;
1013: }

1015: int main(int argc, char **argv)
1016: {
1017:   MPI_Comm  comm;
1018:   TS        ts;
1019:   DM        dm;
1020:   Vec       u;
1021:   AppCtx    user;
1022:   PetscReal t0, t = 0.0;
1023:   void     *ctxs[2];

1025:   ctxs[0] = &t;
1026:   ctxs[1] = &t;
1028:   PetscInitialize(&argc, &argv, (char *)0, help);
1029:   comm                    = PETSC_COMM_WORLD;
1030:   user.functionalRegistry = NULL;
1031:   globalUser              = &user;
1032:   ProcessOptions(comm, &user);
1033:   TSCreate(comm, &ts);
1034:   TSSetType(ts, TSBEULER);
1035:   CreateDM(comm, &user, &dm);
1036:   TSSetDM(ts, dm);
1037:   ProcessMonitorOptions(comm, &user);

1039:   DMCreateGlobalVector(dm, &u);
1040:   PetscObjectSetName((PetscObject)u, "solution");
1041:   if (user.useFV) {
1042:     PetscBool isImplicit = PETSC_FALSE;

1044:     PetscOptionsHasName(NULL, "", "-use_implicit", &isImplicit);
1045:     if (isImplicit) {
1046:       DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user);
1047:       DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user);
1048:     }
1049:     DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user);
1050:     DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &user);
1051:   } else {
1052:     DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user);
1053:     DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user);
1054:     DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user);
1055:   }
1056:   if (user.useFV) TSMonitorSet(ts, MonitorFunctionals, &user, NULL);
1057:   TSSetMaxSteps(ts, 1);
1058:   TSSetMaxTime(ts, 2.0);
1059:   TSSetTimeStep(ts, 0.01);
1060:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1061:   TSSetFromOptions(ts);

1063:   DMProjectFunction(dm, 0.0, user.initialGuess, ctxs, INSERT_VALUES, u);
1064:   if (user.useFV) SetInitialConditionFVM(dm, u, 1, user.initialGuess[1], ctxs[1]);
1065:   VecViewFromOptions(u, NULL, "-init_vec_view");
1066:   TSGetTime(ts, &t);
1067:   t0 = t;
1068:   DMTSCheckFromOptions(ts, u);
1069:   TSSolve(ts, u);
1070:   TSGetTime(ts, &t);
1071:   if (t > t0) DMTSCheckFromOptions(ts, u);
1072:   VecViewFromOptions(u, NULL, "-sol_vec_view");
1073:   {
1074:     PetscReal         ftime;
1075:     PetscInt          nsteps;
1076:     TSConvergedReason reason;

1078:     TSGetSolveTime(ts, &ftime);
1079:     TSGetStepNumber(ts, &nsteps);
1080:     TSGetConvergedReason(ts, &reason);
1081:     PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps);
1082:   }

1084:   VecDestroy(&u);
1085:   DMDestroy(&dm);
1086:   TSDestroy(&ts);
1087:   PetscFree(user.monitorFuncs);
1088:   FunctionalDestroy(&user.functionalRegistry);
1089:   PetscFinalize();
1090:   return 0;
1091: }

1093: /*TEST

1095:   testset:
1096:     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3

1098:     # 2D harmonic velocity, no porosity
1099:     test:
1100:       suffix: p1p1
1101:       requires: !complex !single
1102:       args: -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_factor_shift_type nonzero -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
1103:     test:
1104:       suffix: p1p1_xper
1105:       requires: !complex !single
1106:       args: -dm_refine 1 -dm_plex_box_bd periodic,none -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
1107:     test:
1108:       suffix: p1p1_xper_ref
1109:       requires: !complex !single
1110:       args: -dm_refine 2 -dm_plex_box_bd periodic,none -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
1111:     test:
1112:       suffix: p1p1_xyper
1113:       requires: !complex !single
1114:       args: -dm_refine 1 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
1115:     test:
1116:       suffix: p1p1_xyper_ref
1117:       requires: !complex !single
1118:       args: -dm_refine 2 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
1119:     test:
1120:       suffix: p2p1
1121:       requires: !complex !single
1122:       args: -velocity_petscspace_degree 2 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ts_monitor   -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
1123:     test:
1124:       suffix: p2p1_xyper
1125:       requires: !complex !single
1126:       args: -dm_refine 1 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 2 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check

1128:     test:
1129:       suffix: adv_1
1130:       requires: !complex !single
1131:       args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 1,2,4 -bc_outflow 3 -ts_view -dm_view

1133:     test:
1134:       suffix: adv_2
1135:       requires: !complex
1136:       TODO: broken memory corruption
1137:       args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3,4 -bc_outflow 1,2 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason -ksp_converged_reason

1139:     test:
1140:       suffix: adv_3
1141:       requires: !complex
1142:       TODO: broken memory corruption
1143:       args: -dm_plex_box_bd periodic,none -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason

1145:     test:
1146:       suffix: adv_3_ex
1147:       requires: !complex
1148:       args: -dm_plex_box_bd periodic,none -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.1 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -ksp_max_it 100 -ts_view -dm_view

1150:     test:
1151:       suffix: adv_4
1152:       requires: !complex
1153:       TODO: broken memory corruption
1154:       args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason

1156:     # 2D Advection, box, delta
1157:     test:
1158:       suffix: adv_delta_yper_0
1159:       requires: !complex
1160:       TODO: broken
1161:       args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error

1163:     test:
1164:       suffix: adv_delta_yper_1
1165:       requires: !complex
1166:       TODO: broken
1167:       args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 1 -source_loc 0.416666,0.416666

1169:     test:
1170:       suffix: adv_delta_yper_2
1171:       requires: !complex
1172:       TODO: broken
1173:       args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 2 -source_loc 0.458333,0.458333

1175:     test:
1176:       suffix: adv_delta_yper_fim_0
1177:       requires: !complex
1178:       TODO: broken
1179:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view   -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason

1181:     test:
1182:       suffix: adv_delta_yper_fim_1
1183:       requires: !complex
1184:       TODO: broken
1185:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view   -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -snes_linesearch_type basic

1187:     test:
1188:       suffix: adv_delta_yper_fim_2
1189:       requires: !complex
1190:       TODO: broken
1191:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view   -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -snes_linesearch_type basic

1193:     test:
1194:       suffix: adv_delta_yper_im_0
1195:       requires: !complex
1196:       TODO: broken
1197:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason

1199:     test:
1200:       suffix: adv_delta_yper_im_1
1201:       requires: !complex
1202:       TODO: broken
1203:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 1 -source_loc 0.416666,0.416666

1205:     test:
1206:       suffix: adv_delta_yper_im_2
1207:       requires: !complex
1208:       TODO: broken
1209:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 2 -source_loc 0.458333,0.458333

1211:     test:
1212:       suffix: adv_delta_yper_im_3
1213:       requires: !complex
1214:       TODO: broken
1215:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason

1217:     #    I believe the nullspace is sin(pi y)
1218:     test:
1219:       suffix: adv_delta_yper_im_4
1220:       requires: !complex
1221:       TODO: broken
1222:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 1 -source_loc 0.416666,0.416666

1224:     test:
1225:       suffix: adv_delta_yper_im_5
1226:       requires: !complex
1227:       TODO: broken
1228:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 2 -source_loc 0.458333,0.458333

1230:     test:
1231:       suffix: adv_delta_yper_im_6
1232:       requires: !complex
1233:       TODO: broken
1234:       args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view   -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type svd -snes_converged_reason
1235:     # 2D Advection, magma benchmark 1

1237:     test:
1238:       suffix: adv_delta_shear_im_0
1239:       requires: !complex
1240:       TODO: broken
1241:       args: -dm_plex_box_bd periodic,none -dm_refine 2 -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist shear -porosity_dist   delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 1,3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -source_loc 0.458333,0.708333
1242:     # 2D Advection, box, gaussian

1244:     test:
1245:       suffix: adv_gauss
1246:       requires: !complex
1247:       TODO: broken
1248:       args: -use_fv -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type ssp -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view

1250:     test:
1251:       suffix: adv_gauss_im
1252:       requires: !complex
1253:       TODO: broken
1254:       args: -use_fv -use_implicit -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps   100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7

1256:     test:
1257:       suffix: adv_gauss_im_1
1258:       requires: !complex
1259:       TODO: broken
1260:       args: -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7

1262:     test:
1263:       suffix: adv_gauss_im_2
1264:       requires: !complex
1265:       TODO: broken
1266:       args: -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7

1268:     test:
1269:       suffix: adv_gauss_corner
1270:       requires: !complex
1271:       TODO: broken
1272:       args: -use_fv -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type ssp -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 2 -ts_view -dm_view

1274:     # 2D Advection+Harmonic 12-
1275:     test:
1276:       suffix: adv_harm_0
1277:       requires: !complex
1278:       TODO: broken memory corruption
1279:       args: -velocity_petscspace_degree 2 -use_fv -velocity_dist harmonic -porosity_dist gaussian -ts_type beuler -ts_max_time 2.0 -ts_max_steps   1000 -ts_dt 0.993392 -bc_inflow 1,2,4 -bc_outflow 3 -use_implicit -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it   100 -ts_view -dm_view -snes_converged_reason -ksp_converged_reason -snes_monitor -dmts_check

1281:   #   Must check that FV BCs propagate to coarse meshes
1282:   #   Must check that FV BC ids propagate to coarse meshes
1283:   #   Must check that FE+FV BCs work at the same time
1284:   # 2D Advection, matching wind in ex11 8-11
1285:   #   NOTE implicit solves are limited by accuracy of FD Jacobian
1286:   test:
1287:     suffix: adv_0
1288:     requires: !complex !single exodusii
1289:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view

1291:   test:
1292:     suffix: adv_0_im
1293:     requires: !complex exodusii
1294:     TODO: broken  memory corruption
1295:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu

1297:   test:
1298:     suffix: adv_0_im_2
1299:     requires: !complex exodusii
1300:     TODO: broken
1301:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -snes_rtol 1.0e-7

1303:   test:
1304:     suffix: adv_0_im_3
1305:     requires: !complex exodusii
1306:     TODO: broken
1307:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type svd -snes_rtol 1.0e-7

1309:   test:
1310:     suffix: adv_0_im_4
1311:     requires: !complex exodusii
1312:     TODO: broken
1313:     args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type svd -snes_rtol 1.0e-7
1314:   # 2D Advection, misc

1316: TEST*/