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;

83:   PetscFunctionBeginUser;
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:   PetscCall(PetscOptionsBool("-use_fv", "Use the finite volume method for advection", "ex18.c", options->useFV, &options->useFV, NULL));
95:   vd = options->velocityDist;
96:   PetscCall(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:   PetscCall(PetscOptionsEList("-porosity_dist", "Initial porosity distribution type", "ex18.c", porosityDist, 5, porosityDist[options->porosityDist], &pd, NULL));
100:   options->porosityDist = (PorosityDistribution)pd;
101:   PetscCall(PetscOptionsReal("-inflow_state", "The inflow state", "ex18.c", options->inflowState, &options->inflowState, NULL));
102:   d = 2;
103:   PetscCall(PetscOptionsRealArray("-source_loc", "The source location", "ex18.c", options->source, &d, &flg));
104:   PetscCheck(!flg || d == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must give dim coordinates for the source location, not %" PetscInt_FMT, d);
105:   PetscOptionsEnd();
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

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

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

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

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

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

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

164: {
165:   Functional next, l;

167:   PetscFunctionBeginUser;
171:   for (; l; l = next) {
172:     next = l->next;
173:     PetscCall(PetscFree(l->name));
174:     PetscCall(PetscFree(l));
175:   }
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: 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[])
180: {
181:   PetscInt comp;
182:   for (comp = 0; comp < dim; ++comp) f0[comp] = u[comp];
183: }

185: 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[])
186: {
187:   PetscScalar wind[3] = {0.0, 0.0, 0.0};
188:   PetscInt    comp;

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

194: 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[])
195: {
196:   PetscInt comp;
197:   for (comp = 0; comp < dim * dim; ++comp) f1[comp] = 0.0;
198: }

200: 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[])
201: {
202:   PetscInt d;
203:   for (d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
204: }

206: 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[])
207: {
208:   g0[0] = 1.0;
209: }

211: 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[])
212: {
213:   PetscInt comp;
214:   for (comp = 0; comp < dim; ++comp) f0[comp] = 4.0;
215: }

217: 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[])
218: {
219:   PetscInt comp, d;
220:   for (comp = 0; comp < dim; ++comp) {
221:     for (d = 0; d < dim; ++d) f1[comp * dim + d] = u_x[comp * dim + d];
222:   }
223: }

225: 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[])
226: {
227:   f0[0] = -PetscSinReal(2.0 * PETSC_PI * x[0]);
228:   f0[1] = 2.0 * PETSC_PI * x[1] * PetscCosReal(2.0 * PETSC_PI * x[0]);
229: }

231: 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[])
232: {
233:   f0[0] = -2.0 * PetscSinReal(2.0 * PETSC_PI * x[0]) * PetscCosReal(2.0 * PETSC_PI * x[1]);
234:   f0[1] = 2.0 * PetscSinReal(2.0 * PETSC_PI * x[1]) * PetscCosReal(2.0 * PETSC_PI * x[0]);
235: }

237: 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[])
238: {
239:   const PetscInt Ncomp = dim;
240:   PetscInt       compI, d;

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

247: /* \frac{\partial\phi}{\partial t} + \nabla\phi \cdot \mathbf{u} + \phi \nabla \cdot \mathbf{u} = 0 */
248: 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[])
249: {
250:   PetscInt d;
251:   f0[0] = u_t[dim];
252:   for (d = 0; d < dim; ++d) f0[0] += u[dim] * u_x[d * dim + d] + u_x[dim * dim + d] * u[d];
253: }

255: 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[])
256: {
257:   PetscInt d;
258:   for (d = 0; d < dim; ++d) f1[0] = 0.0;
259: }

261: 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[])
262: {
263:   PetscInt d;
264:   g0[0] = u_tShift;
265:   for (d = 0; d < dim; ++d) g0[0] += u_x[d * dim + d];
266: }

268: 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[])
269: {
270:   PetscInt d;
271:   for (d = 0; d < dim; ++d) g1[d] = u[d];
272: }

274: 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[])
275: {
276:   PetscInt d;
277:   for (d = 0; d < dim; ++d) g0[0] += u_x[dim * dim + d];
278: }

280: 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[])
281: {
282:   PetscInt d;
283:   for (d = 0; d < dim; ++d) g1[d * dim + d] = u[dim];
284: }

286: 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)
287: {
288:   PetscReal wind[3] = {0.0, 1.0, 0.0};
289:   PetscReal wn      = DMPlex_DotRealD_Internal(PetscMin(dim, 3), wind, n);

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

294: 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)
295: {
296:   PetscReal wn = DMPlex_DotD_Internal(dim, uL, n);

298: #if 1
299:   flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn;
300: #else
301:   /* 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]); */
302:   /* Smear it out */
303:   flux[0] = 0.5 * ((uL[dim] + uR[dim]) + (uL[dim] - uR[dim]) * tanh(1.0e5 * wn)) * wn;
304: #endif
305: }

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

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

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

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

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

343:   so that

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

351: We will conserve phi since

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

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

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

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

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

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

380:   so that

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

384: We will conserve phi since

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

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

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

404:   so that

406:     -\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

408: We will conserve phi since

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

526:   PetscFunctionBeginUser;
527:   xG[0] = user->inflowState;
528:   PetscFunctionReturn(PETSC_SUCCESS);
529: }

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

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

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

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

568:   PetscFunctionBeginUser;
569:   PetscCall(ExactSolution(dm, time, x, yexact, ctx));
570:   f[user->errorFunctional] = PetscAbsScalar(y[0] - yexact[0]);
571:   PetscFunctionReturn(PETSC_SUCCESS);
572: }

574: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
575: {
576:   PetscFunctionBeginUser;
577:   PetscCall(DMCreate(comm, dm));
578:   PetscCall(DMSetType(*dm, DMPLEX));
579:   PetscCall(DMSetFromOptions(*dm));
580:   PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view"));
581:   PetscFunctionReturn(PETSC_SUCCESS);
582: }

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

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

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

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

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

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

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

742:   PetscCall(FunctionalRegister(&user->functionalRegistry, "Error", &user->errorFunctional, Functional_Error, user));
743:   PetscFunctionReturn(PETSC_SUCCESS);
744: }

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

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

765:   PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)dm), &fv));
766:   PetscCall(PetscObjectSetName((PetscObject)fv, "porosity"));
767:   PetscCall(PetscFVSetFromOptions(fv));
768:   PetscCall(PetscFVSetNumComponents(fv, 1));
769:   PetscCall(PetscFVSetSpatialDimension(fv, dim));

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

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

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

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

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

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

840:     PetscCall(DMPlexPointGlobalFieldRef(dm, c, field, x, &xc));
841:     if (xc) PetscCall((*func)(dim, 0.0, cg->centroid, Nf, xc, ctx));
842:   }
844:   PetscCall(VecRestoreArray(X, &x));
845:   PetscFunctionReturn(PETSC_SUCCESS);
846: }

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

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

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

890:     size_t ftableused, ftablealloc;

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

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

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

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

958:       /* Create string with functional outputs */
959:       if (f % 3) {
960:         PetscCall(PetscArraycpy(buffer, "  ", 2));
961:         p = buffer + 2;
962:       } else if (f) {
963:         PetscCall(PetscArraycpy(buffer, newline, sizeof(newline) - 1));
964:         p = buffer + sizeof(newline) - 1;
965:       } else {
966:         p = buffer;
967:       }
968:       PetscCall(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]));
969:       countused += p - buffer;
970:       /* reallocate */
971:       if (countused > ftablealloc - ftableused - 1) {
972:         char *ftablenew;

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

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

1024:   PetscFunctionBeginUser;
1025:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1026:   comm                    = PETSC_COMM_WORLD;
1027:   user.functionalRegistry = NULL;
1028:   globalUser              = &user;
1029:   PetscCall(ProcessOptions(comm, &user));
1030:   PetscCall(TSCreate(comm, &ts));
1031:   PetscCall(TSSetType(ts, TSBEULER));
1032:   PetscCall(CreateDM(comm, &user, &dm));
1033:   PetscCall(TSSetDM(ts, dm));
1034:   PetscCall(ProcessMonitorOptions(comm, &user));

1036:   PetscCall(DMCreateGlobalVector(dm, &u));
1037:   PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
1038:   if (user.useFV) {
1039:     PetscBool isImplicit = PETSC_FALSE;

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

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

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

1081:   PetscCall(VecDestroy(&u));
1082:   PetscCall(DMDestroy(&dm));
1083:   PetscCall(TSDestroy(&ts));
1084:   PetscCall(PetscFree(user.monitorFuncs));
1085:   PetscCall(FunctionalDestroy(&user.functionalRegistry));
1086:   PetscCall(PetscFinalize());
1087:   return 0;
1088: }

1090: /*TEST

1092:   testset:
1093:     args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3

1095:     # 2D harmonic velocity, no porosity
1096:     test:
1097:       suffix: p1p1
1098:       requires: !complex !single
1099:       args: -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type ilu -pc_factor_shift_type nonzero -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check
1100:     test:
1101:       suffix: p1p1_xper
1102:       requires: !complex !single
1103:       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
1104:     test:
1105:       suffix: p1p1_xper_ref
1106:       requires: !complex !single
1107:       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
1108:     test:
1109:       suffix: p1p1_xyper
1110:       requires: !complex !single
1111:       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
1112:     test:
1113:       suffix: p1p1_xyper_ref
1114:       requires: !complex !single
1115:       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
1116:     test:
1117:       suffix: p2p1
1118:       requires: !complex !single
1119:       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
1120:     test:
1121:       suffix: p2p1_xyper
1122:       requires: !complex !single
1123:       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

1125:     test:
1127:       requires: !complex !single
1128:       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

1130:     test:
1132:       requires: !complex
1133:       TODO: broken memory corruption
1134:       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

1136:     test:
1138:       requires: !complex
1139:       TODO: broken memory corruption
1140:       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

1142:     test:
1144:       requires: !complex
1145:       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

1147:     test:
1149:       requires: !complex
1150:       TODO: broken memory corruption
1151:       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

1153:     # 2D Advection, box, delta
1154:     test:
1156:       requires: !complex
1157:       TODO: broken
1158:       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

1160:     test:
1162:       requires: !complex
1163:       TODO: broken
1164:       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

1166:     test:
1168:       requires: !complex
1169:       TODO: broken
1170:       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

1172:     test:
1174:       requires: !complex
1175:       TODO: broken
1176:       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

1178:     test:
1180:       requires: !complex
1181:       TODO: broken
1182:       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

1184:     test:
1186:       requires: !complex
1187:       TODO: broken
1188:       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

1190:     test:
1192:       requires: !complex
1193:       TODO: broken
1194:       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

1196:     test:
1198:       requires: !complex
1199:       TODO: broken
1200:       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

1202:     test:
1204:       requires: !complex
1205:       TODO: broken
1206:       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

1208:     test:
1210:       requires: !complex
1211:       TODO: broken
1212:       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

1214:     #    I believe the nullspace is sin(pi y)
1215:     test:
1217:       requires: !complex
1218:       TODO: broken
1219:       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

1221:     test:
1223:       requires: !complex
1224:       TODO: broken
1225:       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

1227:     test:
1229:       requires: !complex
1230:       TODO: broken
1231:       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
1232:     # 2D Advection, magma benchmark 1

1234:     test:
1236:       requires: !complex
1237:       TODO: broken
1238:       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
1239:     # 2D Advection, box, gaussian

1241:     test:
1243:       requires: !complex
1244:       TODO: broken
1245:       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

1247:     test:
1249:       requires: !complex
1250:       TODO: broken
1251:       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

1253:     test:
1255:       requires: !complex
1256:       TODO: broken
1257:       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

1259:     test:
1261:       requires: !complex
1262:       TODO: broken
1263:       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

1265:     test:
1267:       requires: !complex
1268:       TODO: broken
1269:       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

1272:     test:
1274:       requires: !complex
1275:       TODO: broken memory corruption
1276:       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

1278:   #   Must check that FV BCs propagate to coarse meshes
1279:   #   Must check that FV BC ids propagate to coarse meshes
1280:   #   Must check that FE+FV BCs work at the same time
1281:   # 2D Advection, matching wind in ex11 8-11
1282:   #   NOTE implicit solves are limited by accuracy of FD Jacobian
1283:   test:
1285:     requires: !complex !single exodusii
1286:     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 1288: test: 1289: suffix: adv_0_im 1290: requires: !complex exodusii 1291: TODO: broken memory corruption 1292: 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

1294:   test:
1296:     requires: !complex exodusii
1297:     TODO: broken
1298:     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 1300: test: 1301: suffix: adv_0_im_3 1302: requires: !complex exodusii 1303: TODO: broken 1304: 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

1306:   test:
1308:     requires: !complex exodusii
1309:     TODO: broken
1310:     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