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: FunctionalFn - 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 (*FunctionalFn)(DM, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
50: typedef struct _n_Functional *Functional;
51: struct _n_Functional {
52: char *name;
53: FunctionalFn 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, PetscCtx 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];
114: PetscFunctionBeginUser;
115: PetscOptionsBegin(comm, "", "Simulation Monitor Options", "DMPLEX");
116: options->numMonitorFuncs = PETSC_STATIC_ARRAY_LENGTH(names);
117: PetscCall(PetscOptionsStringArray("-monitor", "List of functionals to monitor", "", names, &options->numMonitorFuncs, NULL));
118: PetscCall(PetscMalloc1(options->numMonitorFuncs, &options->monitorFuncs));
119: for (PetscInt f = 0; f < options->numMonitorFuncs; ++f) {
120: for (func = options->functionalRegistry; func; func = func->next) {
121: PetscBool match;
123: PetscCall(PetscStrcasecmp(names[f], func->name, &match));
124: if (match) break;
125: }
126: PetscCheck(func, comm, PETSC_ERR_USER, "No known functional '%s'", names[f]);
127: options->monitorFuncs[f] = func;
128: /* Jed inserts a de-duplication of functionals here */
129: PetscCall(PetscFree(names[f]));
130: }
131: /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
132: options->maxMonitorFunc = -1;
133: for (func = options->functionalRegistry; func; func = func->next) {
134: for (PetscInt f = 0; f < options->numMonitorFuncs; ++f) {
135: Functional call = options->monitorFuncs[f];
137: if (func->func == call->func && func->ctx == call->ctx) options->maxMonitorFunc = PetscMax(options->maxMonitorFunc, func->offset);
138: }
139: }
140: PetscOptionsEnd();
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode FunctionalRegister(Functional *functionalRegistry, const char name[], PetscInt *offset, FunctionalFn func, PetscCtx ctx)
145: {
146: Functional *ptr, f;
147: PetscInt lastoffset = -1;
149: PetscFunctionBeginUser;
150: for (ptr = functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
151: PetscCall(PetscNew(&f));
152: PetscCall(PetscStrallocpy(name, &f->name));
153: f->offset = lastoffset + 1;
154: f->func = func;
155: f->ctx = ctx;
156: f->next = NULL;
157: *ptr = f;
158: *offset = f->offset;
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode FunctionalDestroy(Functional *link)
163: {
164: Functional next, l;
166: PetscFunctionBeginUser;
167: if (!link) PetscFunctionReturn(PETSC_SUCCESS);
168: l = *link;
169: *link = NULL;
170: for (; l; l = next) {
171: next = l->next;
172: PetscCall(PetscFree(l->name));
173: PetscCall(PetscFree(l));
174: }
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: 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[])
179: {
180: PetscInt comp;
181: for (comp = 0; comp < dim; ++comp) f0[comp] = u[comp];
182: }
184: 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[])
185: {
186: PetscScalar wind[3] = {0.0, 0.0, 0.0};
187: PetscInt comp;
189: PetscCallAbort(PETSC_COMM_SELF, constant_u_2d(dim, t, x, Nf, wind, NULL));
190: for (comp = 0; comp < dim && comp < 3; ++comp) f0[comp] = u[comp] - wind[comp];
191: }
193: 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[])
194: {
195: PetscInt comp;
196: for (comp = 0; comp < dim * dim; ++comp) f1[comp] = 0.0;
197: }
199: 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[])
200: {
201: PetscInt d;
202: for (d = 0; d < dim; ++d) g0[d * dim + d] = 1.0;
203: }
205: 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[])
206: {
207: g0[0] = 1.0;
208: }
210: 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[])
211: {
212: for (PetscInt comp = 0; comp < dim; ++comp) f0[comp] = 4.0;
213: }
215: 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[])
216: {
217: for (PetscInt comp = 0; comp < dim; ++comp)
218: for (PetscInt d = 0; d < dim; ++d) f1[comp * dim + d] = u_x[comp * dim + d];
219: }
221: 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[])
222: {
223: f0[0] = -PetscSinReal(2.0 * PETSC_PI * x[0]);
224: f0[1] = 2.0 * PETSC_PI * x[1] * PetscCosReal(2.0 * PETSC_PI * x[0]);
225: }
227: 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[])
228: {
229: f0[0] = -2.0 * PetscSinReal(2.0 * PETSC_PI * x[0]) * PetscCosReal(2.0 * PETSC_PI * x[1]);
230: f0[1] = 2.0 * PetscSinReal(2.0 * PETSC_PI * x[1]) * PetscCosReal(2.0 * PETSC_PI * x[0]);
231: }
233: 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[])
234: {
235: const PetscInt Ncomp = dim;
236: PetscInt d;
238: for (PetscInt compI = 0; compI < Ncomp; ++compI) {
239: for (d = 0; d < dim; ++d) g3[((compI * Ncomp + compI) * dim + d) * dim + d] = 1.0;
240: }
241: }
243: /* \frac{\partial\phi}{\partial t} + \nabla\phi \cdot \mathbf{u} + \phi \nabla \cdot \mathbf{u} = 0 */
244: 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[])
245: {
246: PetscInt d;
247: f0[0] = u_t[dim];
248: for (d = 0; d < dim; ++d) f0[0] += u[dim] * u_x[d * dim + d] + u_x[dim * dim + d] * u[d];
249: }
251: 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[])
252: {
253: PetscInt d;
254: for (d = 0; d < dim; ++d) f1[0] = 0.0;
255: }
257: 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[])
258: {
259: PetscInt d;
260: g0[0] = u_tShift;
261: for (d = 0; d < dim; ++d) g0[0] += u_x[d * dim + d];
262: }
264: 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[])
265: {
266: PetscInt d;
267: for (d = 0; d < dim; ++d) g1[d] = u[d];
268: }
270: 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[])
271: {
272: for (PetscInt d = 0; d < dim; ++d) g0[0] += u_x[dim * dim + d];
273: }
275: 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[])
276: {
277: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = u[dim];
278: }
280: 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, PetscCtx ctx)
281: {
282: PetscReal wind[3] = {0.0, 1.0, 0.0};
283: PetscReal wn = DMPlex_DotRealD_Internal(PetscMin(dim, 3), wind, n);
285: flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn;
286: }
288: 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, PetscCtx ctx)
289: {
290: PetscReal wn = DMPlex_DotD_Internal(dim, uL, n);
292: #if 1
293: flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn;
294: #else
295: /* 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]); */
296: /* Smear it out */
297: flux[0] = 0.5 * ((uL[dim] + uR[dim]) + (uL[dim] - uR[dim]) * tanh(1.0e5 * wn)) * wn;
298: #endif
299: }
301: static PetscErrorCode zero_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
302: {
303: u[0] = 0.0;
304: u[1] = 0.0;
305: return PETSC_SUCCESS;
306: }
308: static PetscErrorCode constant_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
309: {
310: u[0] = 0.0;
311: u[1] = 1.0;
312: return PETSC_SUCCESS;
313: }
315: /* Coordinates of the point which was at x at t = 0 */
316: static PetscErrorCode constant_x_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
317: {
318: const PetscReal t = *((PetscReal *)ctx);
319: u[0] = x[0];
320: u[1] = x[1] + t;
321: #if 0
322: PetscCall(DMLocalizeCoordinate(globalUser->dm, u, PETSC_FALSE, u));
323: #else
324: u[1] = u[1] - (int)PetscRealPart(u[1]);
325: #endif
326: return PETSC_SUCCESS;
327: }
329: /*
330: In 2D we use the exact solution:
332: u = x^2 + y^2
333: v = 2 x^2 - 2xy
334: phi = h(x + y + (u + v) t)
335: f_x = f_y = 4
337: so that
339: -\Delta u + f = <-4, -4> + <4, 4> = 0
340: {\partial\phi}{\partial t} - \nabla\cdot \phi u = 0
341: h_t(x + y + (u + v) t) - u . grad phi - phi div u
342: = u h' + v h' - u h_x - v h_y
343: = 0
345: We will conserve phi since
347: \nabla \cdot u = 2x - 2x = 0
349: Also try h((x + ut)^2 + (y + vt)^2), so that
351: h_t((x + ut)^2 + (y + vt)^2) - u . grad phi - phi div u
352: = 2 h' (u (x + ut) + v (y + vt)) - u h_x - v h_y
353: = 2 h' (u (x + ut) + v (y + vt)) - u h' 2 (x + u t) - v h' 2 (y + vt)
354: = 2 h' (u (x + ut) + v (y + vt) - u (x + u t) - v (y + vt))
355: = 0
357: */
358: static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
359: {
360: u[0] = x[0] * x[0] + x[1] * x[1];
361: u[1] = 2.0 * x[0] * x[0] - 2.0 * x[0] * x[1];
362: return PETSC_SUCCESS;
363: }
365: /*
366: In 2D we use the exact, periodic solution:
368: u = sin(2 pi x)/4 pi^2
369: v = -y cos(2 pi x)/2 pi
370: phi = h(x + y + (u + v) t)
371: f_x = -sin(2 pi x)
372: f_y = 2 pi y cos(2 pi x)
374: so that
376: -\Delta u + f = <sin(2pi x), -2pi y cos(2pi x)> + <-sin(2pi x), 2pi y cos(2pi x)> = 0
378: We will conserve phi since
380: \nabla \cdot u = cos(2pi x)/2pi - cos(2pi x)/2pi = 0
381: */
382: static PetscErrorCode periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
383: {
384: u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]) / PetscSqr(2.0 * PETSC_PI);
385: u[1] = -x[1] * PetscCosReal(2.0 * PETSC_PI * x[0]) / (2.0 * PETSC_PI);
386: return PETSC_SUCCESS;
387: }
389: /*
390: In 2D we use the exact, doubly periodic solution:
392: u = sin(2 pi x) cos(2 pi y)/4 pi^2
393: v = -sin(2 pi y) cos(2 pi x)/4 pi^2
394: phi = h(x + y + (u + v) t)
395: f_x = -2sin(2 pi x) cos(2 pi y)
396: f_y = 2sin(2 pi y) cos(2 pi x)
398: so that
400: -\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
402: We will conserve phi since
404: \nabla \cdot u = cos(2pi x) cos(2pi y)/2pi - cos(2pi y) cos(2pi x)/2pi = 0
405: */
406: static PetscErrorCode doubly_periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
407: {
408: u[0] = PetscSinReal(2.0 * PETSC_PI * x[0]) * PetscCosReal(2.0 * PETSC_PI * x[1]) / PetscSqr(2.0 * PETSC_PI);
409: u[1] = -PetscSinReal(2.0 * PETSC_PI * x[1]) * PetscCosReal(2.0 * PETSC_PI * x[0]) / PetscSqr(2.0 * PETSC_PI);
410: return PETSC_SUCCESS;
411: }
413: static PetscErrorCode shear_bc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
414: {
415: u[0] = x[1] - 0.5;
416: u[1] = 0.0;
417: return PETSC_SUCCESS;
418: }
420: static PetscErrorCode initialVelocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
421: {
422: for (PetscInt d = 0; d < dim; ++d) u[d] = 0.0;
423: return PETSC_SUCCESS;
424: }
426: static PetscErrorCode zero_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
427: {
428: u[0] = 0.0;
429: return PETSC_SUCCESS;
430: }
432: static PetscErrorCode constant_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
433: {
434: u[0] = 1.0;
435: return PETSC_SUCCESS;
436: }
438: static PetscErrorCode delta_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
439: {
440: PetscReal x0[2];
441: PetscScalar xn[2];
443: x0[0] = globalUser->source[0];
444: x0[1] = globalUser->source[1];
445: PetscCall(constant_x_2d(dim, time, x0, Nf, xn, ctx));
446: {
447: const PetscReal xi = x[0] - PetscRealPart(xn[0]);
448: const PetscReal eta = x[1] - PetscRealPart(xn[1]);
449: const PetscReal r2 = xi * xi + eta * eta;
451: u[0] = r2 < 1.0e-7 ? 1.0 : 0.0;
452: }
453: return PETSC_SUCCESS;
454: }
456: /*
457: Gaussian blob, initially centered on (0.5, 0.5)
459: xi = x(t) - x0, eta = y(t) - y0
461: where x(t), y(t) are the integral curves of v(t),
463: dx/dt . grad f = v . f
465: Check: constant v(t) = {v0, w0}, x(t) = {x0 + v0 t, y0 + w0 t}
467: v0 f_x + w0 f_y = v . f
468: */
469: static PetscErrorCode gaussian_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
470: {
471: const PetscReal x0[2] = {0.5, 0.5};
472: const PetscReal sigma = 1.0 / 6.0;
473: PetscScalar xn[2];
475: PetscCall(constant_x_2d(dim, time, x0, Nf, xn, ctx));
476: {
477: /* const PetscReal xi = x[0] + (sin(2.0*PETSC_PI*x[0])/(4.0*PETSC_PI*PETSC_PI))*t - x0[0]; */
478: /* const PetscReal eta = x[1] + (-x[1]*cos(2.0*PETSC_PI*x[0])/(2.0*PETSC_PI))*t - x0[1]; */
479: const PetscReal xi = x[0] - PetscRealPart(xn[0]);
480: const PetscReal eta = x[1] - PetscRealPart(xn[1]);
481: const PetscReal r2 = xi * xi + eta * eta;
483: u[0] = PetscExpReal(-r2 / (2.0 * sigma * sigma)) / (sigma * PetscSqrtReal(2.0 * PETSC_PI));
484: }
485: return PETSC_SUCCESS;
486: }
488: static PetscErrorCode tilted_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
489: {
490: PetscReal x0[3];
491: const PetscReal wind[3] = {0.0, 1.0, 0.0};
492: const PetscReal t = *((PetscReal *)ctx);
494: DMPlex_WaxpyD_Internal(2, -t, wind, x, x0);
495: if (x0[1] > 0) u[0] = 1.0 * x[0] + 3.0 * x[1];
496: else u[0] = -2.0; /* Inflow state */
497: return PETSC_SUCCESS;
498: }
500: static PetscErrorCode tilted_phi_coupled_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
501: {
502: PetscReal ur[3];
503: PetscReal x0[3];
504: const PetscReal t = *((PetscReal *)ctx);
506: ur[0] = PetscRealPart(u[0]);
507: ur[1] = PetscRealPart(u[1]);
508: ur[2] = PetscRealPart(u[2]);
509: DMPlex_WaxpyD_Internal(2, -t, ur, x, x0);
510: if (x0[1] > 0) u[0] = 1.0 * x[0] + 3.0 * x[1];
511: else u[0] = -2.0; /* Inflow state */
512: return PETSC_SUCCESS;
513: }
515: static PetscErrorCode advect_inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, PetscCtx ctx)
516: {
517: AppCtx *user = (AppCtx *)ctx;
519: PetscFunctionBeginUser;
520: xG[0] = user->inflowState;
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: static PetscErrorCode advect_outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, PetscCtx ctx)
525: {
526: PetscFunctionBeginUser;
527: //xG[0] = xI[dim];
528: xG[0] = xI[2];
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: static PetscErrorCode ExactSolution(DM dm, PetscReal time, const PetscReal *x, PetscScalar *u, PetscCtx ctx)
533: {
534: AppCtx *user = (AppCtx *)ctx;
535: PetscInt dim;
537: PetscFunctionBeginUser;
538: PetscCall(DMGetDimension(dm, &dim));
539: switch (user->porosityDist) {
540: case TILTED:
541: if (user->velocityDist == VEL_ZERO) PetscCall(tilted_phi_2d(dim, time, x, 2, u, (void *)&time));
542: else PetscCall(tilted_phi_coupled_2d(dim, time, x, 2, u, (void *)&time));
543: break;
544: case GAUSSIAN:
545: PetscCall(gaussian_phi_2d(dim, time, x, 2, u, (void *)&time));
546: break;
547: case DELTA:
548: PetscCall(delta_phi_2d(dim, time, x, 2, u, (void *)&time));
549: break;
550: default:
551: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
552: }
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static PetscErrorCode Functional_Error(DM dm, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, PetscCtx ctx)
557: {
558: AppCtx *user = (AppCtx *)ctx;
559: PetscScalar yexact[3] = {0, 0, 0};
561: PetscFunctionBeginUser;
562: PetscCall(ExactSolution(dm, time, x, yexact, ctx));
563: f[user->errorFunctional] = PetscAbsScalar(y[0] - yexact[0]);
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
568: {
569: PetscFunctionBeginUser;
570: PetscCall(DMCreate(comm, dm));
571: PetscCall(DMSetType(*dm, DMPLEX));
572: PetscCall(DMSetFromOptions(*dm));
573: PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view"));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: static PetscErrorCode SetupBC(DM dm, AppCtx *user)
578: {
579: PetscErrorCode (*exactFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
580: DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
581: PetscDS prob;
582: DMLabel label;
583: PetscBool check;
584: PetscInt dim, n = 3;
585: const char *prefix;
587: PetscFunctionBeginUser;
588: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
589: PetscCall(PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *)bdt, &n, NULL));
590: PetscCall(DMGetDimension(dm, &dim));
591: /* Set initial guesses and exact solutions */
592: switch (dim) {
593: case 2:
594: user->initialGuess[0] = initialVelocity;
595: switch (user->porosityDist) {
596: case ZERO:
597: user->initialGuess[1] = zero_phi;
598: break;
599: case CONSTANT:
600: user->initialGuess[1] = constant_phi;
601: break;
602: case GAUSSIAN:
603: user->initialGuess[1] = gaussian_phi_2d;
604: break;
605: case DELTA:
606: user->initialGuess[1] = delta_phi_2d;
607: break;
608: case TILTED:
609: if (user->velocityDist == VEL_ZERO) user->initialGuess[1] = tilted_phi_2d;
610: else user->initialGuess[1] = tilted_phi_coupled_2d;
611: break;
612: }
613: break;
614: default:
615: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", dim);
616: }
617: exactFuncs[0] = user->initialGuess[0];
618: exactFuncs[1] = user->initialGuess[1];
619: switch (dim) {
620: case 2:
621: switch (user->velocityDist) {
622: case VEL_ZERO:
623: exactFuncs[0] = zero_u_2d;
624: break;
625: case VEL_CONSTANT:
626: exactFuncs[0] = constant_u_2d;
627: break;
628: case VEL_HARMONIC:
629: switch (bdt[0]) {
630: case DM_BOUNDARY_PERIODIC:
631: switch (bdt[1]) {
632: case DM_BOUNDARY_PERIODIC:
633: exactFuncs[0] = doubly_periodic_u_2d;
634: break;
635: default:
636: exactFuncs[0] = periodic_u_2d;
637: break;
638: }
639: break;
640: default:
641: exactFuncs[0] = quadratic_u_2d;
642: break;
643: }
644: break;
645: case VEL_SHEAR:
646: exactFuncs[0] = shear_bc;
647: break;
648: default:
649: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %" PetscInt_FMT, dim);
650: }
651: break;
652: default:
653: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %" PetscInt_FMT " not supported", dim);
654: }
655: {
656: PetscBool isImplicit = PETSC_FALSE;
658: PetscCall(PetscOptionsHasName(NULL, "", "-use_implicit", &isImplicit));
659: if (user->velocityDist == VEL_CONSTANT && !isImplicit) user->initialGuess[0] = exactFuncs[0];
660: }
661: PetscCall(PetscOptionsHasName(NULL, NULL, "-dmts_check", &check));
662: if (check) {
663: user->initialGuess[0] = exactFuncs[0];
664: user->initialGuess[1] = exactFuncs[1];
665: }
666: /* Set BC */
667: PetscCall(DMGetDS(dm, &prob));
668: PetscCall(DMGetLabel(dm, "marker", &label));
669: PetscCall(PetscDSSetExactSolution(prob, 0, exactFuncs[0], user));
670: PetscCall(PetscDSSetExactSolution(prob, 1, exactFuncs[1], user));
671: if (label) {
672: const PetscInt id = 1;
674: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, user, NULL));
675: }
676: PetscCall(DMGetLabel(dm, "Face Sets", &label));
677: if (label && user->useFV) {
678: const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
680: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 1, 0, NULL, (PetscVoidFn *)advect_inflow, NULL, user, NULL));
681: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 1, 0, NULL, (PetscVoidFn *)advect_outflow, NULL, user, NULL));
682: }
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
687: {
688: DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
689: PetscDS prob;
690: PetscInt n = 3;
691: const char *prefix;
693: PetscFunctionBeginUser;
694: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
695: PetscCall(PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *)bdt, &n, NULL));
696: PetscCall(DMGetDS(dm, &prob));
697: switch (user->velocityDist) {
698: case VEL_ZERO:
699: PetscCall(PetscDSSetResidual(prob, 0, f0_zero_u, f1_constant_u));
700: break;
701: case VEL_CONSTANT:
702: PetscCall(PetscDSSetResidual(prob, 0, f0_constant_u, f1_constant_u));
703: PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_constant_uu, NULL, NULL, NULL));
704: PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_constant_pp, NULL, NULL, NULL));
705: break;
706: case VEL_HARMONIC:
707: switch (bdt[0]) {
708: case DM_BOUNDARY_PERIODIC:
709: switch (bdt[1]) {
710: case DM_BOUNDARY_PERIODIC:
711: PetscCall(PetscDSSetResidual(prob, 0, f0_lap_doubly_periodic_u, f1_lap_u));
712: break;
713: default:
714: PetscCall(PetscDSSetResidual(prob, 0, f0_lap_periodic_u, f1_lap_u));
715: break;
716: }
717: break;
718: default:
719: PetscCall(PetscDSSetResidual(prob, 0, f0_lap_u, f1_lap_u));
720: break;
721: }
722: PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu));
723: break;
724: case VEL_SHEAR:
725: PetscCall(PetscDSSetResidual(prob, 0, f0_zero_u, f1_lap_u));
726: PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu));
727: break;
728: }
729: PetscCall(PetscDSSetResidual(prob, 1, f0_advection, f1_advection));
730: PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_adv_pp, g1_adv_pp, NULL, NULL));
731: PetscCall(PetscDSSetJacobian(prob, 1, 0, g0_adv_pu, g1_adv_pu, NULL, NULL));
732: if (user->velocityDist == VEL_ZERO) PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_advection));
733: else PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_coupled_advection));
735: PetscCall(FunctionalRegister(&user->functionalRegistry, "Error", &user->errorFunctional, Functional_Error, user));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
740: {
741: DM cdm = dm;
742: PetscQuadrature q;
743: PetscFE fe[2];
744: PetscFV fv;
745: MPI_Comm comm;
746: PetscInt dim;
748: PetscFunctionBeginUser;
749: /* Create finite element */
750: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
751: PetscCall(DMGetDimension(dm, &dim));
752: PetscCall(PetscFECreateDefault(comm, dim, dim, PETSC_FALSE, "velocity_", PETSC_DEFAULT, &fe[0]));
753: PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
754: PetscCall(PetscFECreateDefault(comm, dim, 1, PETSC_FALSE, "porosity_", PETSC_DEFAULT, &fe[1]));
755: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
756: PetscCall(PetscObjectSetName((PetscObject)fe[1], "porosity"));
758: PetscCall(PetscFVCreate(PetscObjectComm((PetscObject)dm), &fv));
759: PetscCall(PetscObjectSetName((PetscObject)fv, "porosity"));
760: PetscCall(PetscFVSetFromOptions(fv));
761: PetscCall(PetscFVSetNumComponents(fv, 1));
762: PetscCall(PetscFVSetSpatialDimension(fv, dim));
763: PetscCall(PetscFEGetQuadrature(fe[0], &q));
764: PetscCall(PetscFVSetQuadrature(fv, q));
766: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
767: if (user->useFV) PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fv));
768: else PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
769: PetscCall(DMCreateDS(dm));
770: PetscCall(SetupProblem(dm, user));
772: /* Set discretization and boundary conditions for each mesh */
773: while (cdm) {
774: PetscCall(DMCopyDisc(dm, cdm));
775: PetscCall(DMGetCoarseDM(cdm, &cdm));
776: /* Coordinates were never localized for coarse meshes */
777: if (cdm) PetscCall(DMLocalizeCoordinates(cdm));
778: }
779: PetscCall(PetscFEDestroy(&fe[0]));
780: PetscCall(PetscFEDestroy(&fe[1]));
781: PetscCall(PetscFVDestroy(&fv));
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
785: static PetscErrorCode CreateDM(MPI_Comm comm, AppCtx *user, DM *dm)
786: {
787: PetscFunctionBeginUser;
788: PetscCall(CreateMesh(comm, user, dm));
789: /* Handle refinement, etc. */
790: PetscCall(DMSetFromOptions(*dm));
791: /* Construct ghost cells */
792: if (user->useFV) {
793: DM gdm;
795: PetscCall(DMPlexConstructGhostCells(*dm, NULL, NULL, &gdm));
796: PetscCall(DMDestroy(dm));
797: *dm = gdm;
798: }
799: /* Localize coordinates */
800: PetscCall(DMLocalizeCoordinates(*dm));
801: PetscCall(PetscObjectSetName((PetscObject)*dm, "Mesh"));
802: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
803: /* Setup problem */
804: PetscCall(SetupDiscretization(*dm, user));
805: /* Setup BC */
806: PetscCall(SetupBC(*dm, user));
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: static PetscErrorCode SetInitialConditionFVM(DM dm, Vec X, PetscInt field, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscCtx ctx)
811: {
812: PetscDS prob;
813: DM dmCell;
814: Vec cellgeom;
815: const PetscScalar *cgeom;
816: PetscScalar *x;
817: PetscInt dim, Nf, cStart, cEnd, c;
819: PetscFunctionBeginUser;
820: PetscCall(DMGetDS(dm, &prob));
821: PetscCall(DMGetDimension(dm, &dim));
822: PetscCall(PetscDSGetNumFields(prob, &Nf));
823: PetscCall(DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL));
824: PetscCall(VecGetDM(cellgeom, &dmCell));
825: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
826: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
827: PetscCall(VecGetArray(X, &x));
828: for (c = cStart; c < cEnd; ++c) {
829: PetscFVCellGeom *cg;
830: PetscScalar *xc;
832: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
833: PetscCall(DMPlexPointGlobalFieldRef(dm, c, field, x, &xc));
834: if (xc) PetscCall((*func)(dim, 0.0, cg->centroid, Nf, xc, ctx));
835: }
836: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
837: PetscCall(VecRestoreArray(X, &x));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: static PetscErrorCode MonitorFunctionals(TS ts, PetscInt stepnum, PetscReal time, Vec X, PetscCtx ctx)
842: {
843: AppCtx *user = (AppCtx *)ctx;
844: char *ftable = NULL;
845: DM dm;
846: PetscSection s;
847: Vec cellgeom;
848: const PetscScalar *x;
849: PetscScalar *a;
850: PetscReal *xnorms;
851: PetscInt pStart, pEnd, p, Nf, f;
853: PetscFunctionBeginUser;
854: PetscCall(VecViewFromOptions(X, (PetscObject)ts, "-view_solution"));
855: PetscCall(VecGetDM(X, &dm));
856: PetscCall(DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL));
857: PetscCall(DMGetLocalSection(dm, &s));
858: PetscCall(PetscSectionGetNumFields(s, &Nf));
859: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
860: PetscCall(PetscCalloc1(Nf * 2, &xnorms));
861: PetscCall(VecGetArrayRead(X, &x));
862: for (p = pStart; p < pEnd; ++p) {
863: for (f = 0; f < Nf; ++f) {
864: PetscInt dof, cdof, d;
866: PetscCall(PetscSectionGetFieldDof(s, p, f, &dof));
867: PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof));
868: PetscCall(DMPlexPointGlobalFieldRead(dm, p, f, x, &a));
869: /* TODO Use constrained indices here */
870: for (d = 0; d < dof - cdof; ++d) xnorms[f * 2 + 0] = PetscMax(xnorms[f * 2 + 0], PetscAbsScalar(a[d]));
871: for (d = 0; d < dof - cdof; ++d) xnorms[f * 2 + 1] += PetscAbsScalar(a[d]);
872: }
873: }
874: PetscCall(VecRestoreArrayRead(X, &x));
875: if (stepnum >= 0) { /* No summary for final time */
876: DM dmCell, *fdm;
877: Vec *fv;
878: const PetscScalar *cgeom;
879: PetscScalar **fx;
880: PetscReal *fmin, *fmax, *fint, *ftmp, t;
881: PetscInt cStart, cEnd, c, fcount, f, num;
883: size_t ftableused, ftablealloc;
885: /* Functionals have indices after registering, this is an upper bound */
886: fcount = user->numMonitorFuncs;
887: PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fint, fcount, &ftmp));
888: PetscCall(PetscMalloc3(fcount, &fdm, fcount, &fv, fcount, &fx));
889: for (f = 0; f < fcount; ++f) {
890: PetscSection fs;
891: const char *name = user->monitorFuncs[f]->name;
893: fmin[f] = PETSC_MAX_REAL;
894: fmax[f] = PETSC_MIN_REAL;
895: fint[f] = 0;
896: /* Make monitor vecs */
897: PetscCall(DMClone(dm, &fdm[f]));
898: PetscCall(DMGetOutputSequenceNumber(dm, &num, &t));
899: PetscCall(DMSetOutputSequenceNumber(fdm[f], num, t));
900: PetscCall(PetscSectionClone(s, &fs));
901: PetscCall(PetscSectionSetFieldName(fs, 0, NULL));
902: PetscCall(PetscSectionSetFieldName(fs, 1, name));
903: PetscCall(DMSetLocalSection(fdm[f], fs));
904: PetscCall(PetscSectionDestroy(&fs));
905: PetscCall(DMGetGlobalVector(fdm[f], &fv[f]));
906: PetscCall(PetscObjectSetName((PetscObject)fv[f], name));
907: PetscCall(VecGetArray(fv[f], &fx[f]));
908: }
909: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
910: PetscCall(VecGetDM(cellgeom, &dmCell));
911: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
912: PetscCall(VecGetArrayRead(X, &x));
913: for (c = cStart; c < cEnd; ++c) {
914: PetscFVCellGeom *cg;
915: PetscScalar *cx;
917: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
918: PetscCall(DMPlexPointGlobalFieldRead(dm, c, 1, x, &cx));
919: if (!cx) continue; /* not a global cell */
920: for (f = 0; f < user->numMonitorFuncs; ++f) {
921: Functional func = user->monitorFuncs[f];
922: PetscScalar *fxc;
924: PetscCall(DMPlexPointGlobalFieldRef(dm, c, 1, fx[f], &fxc));
925: /* I need to make it easier to get interpolated values here */
926: PetscCall((*func->func)(dm, time, cg->centroid, cx, ftmp, func->ctx));
927: fxc[0] = ftmp[user->monitorFuncs[f]->offset];
928: }
929: for (f = 0; f < fcount; ++f) {
930: fmin[f] = PetscMin(fmin[f], ftmp[f]);
931: fmax[f] = PetscMax(fmax[f], ftmp[f]);
932: fint[f] += cg->volume * ftmp[f];
933: }
934: }
935: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
936: PetscCall(VecRestoreArrayRead(X, &x));
937: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmin, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
938: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmax, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
939: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fint, (PetscMPIInt)fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
940: /* Output functional data */
941: ftablealloc = fcount * 100;
942: ftableused = 0;
943: PetscCall(PetscCalloc1(ftablealloc, &ftable));
944: for (f = 0; f < user->numMonitorFuncs; ++f) {
945: Functional func = user->monitorFuncs[f];
946: PetscInt id = func->offset;
947: char newline[] = "\n";
948: char buffer[256], *p, *prefix;
949: size_t countused, len;
951: /* Create string with functional outputs */
952: if (f % 3) {
953: PetscCall(PetscArraycpy(buffer, " ", 2));
954: p = buffer + 2;
955: } else if (f) {
956: PetscCall(PetscArraycpy(buffer, newline, sizeof(newline) - 1));
957: p = buffer + sizeof(newline) - 1;
958: } else {
959: p = buffer;
960: }
961: 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]));
962: countused += p - buffer;
963: /* reallocate */
964: if (countused > ftablealloc - ftableused - 1) {
965: char *ftablenew;
967: ftablealloc = 2 * ftablealloc + countused;
968: PetscCall(PetscMalloc1(ftablealloc, &ftablenew));
969: PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
970: PetscCall(PetscFree(ftable));
971: ftable = ftablenew;
972: }
973: PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
974: ftableused += countused;
975: ftable[ftableused] = 0;
976: /* Output vecs */
977: PetscCall(VecRestoreArray(fv[f], &fx[f]));
978: PetscCall(PetscStrlen(func->name, &len));
979: PetscCall(PetscMalloc1(len + 2, &prefix));
980: PetscCall(PetscStrncpy(prefix, func->name, len + 2));
981: PetscCall(PetscStrlcat(prefix, "_", len + 2));
982: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fv[f], prefix));
983: PetscCall(VecViewFromOptions(fv[f], NULL, "-vec_view"));
984: PetscCall(PetscFree(prefix));
985: PetscCall(DMRestoreGlobalVector(fdm[f], &fv[f]));
986: PetscCall(DMDestroy(&fdm[f]));
987: }
988: PetscCall(PetscFree4(fmin, fmax, fint, ftmp));
989: PetscCall(PetscFree3(fdm, fv, fx));
990: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT " time %8.4g |x| (", stepnum, (double)time));
991: for (f = 0; f < Nf; ++f) {
992: if (f > 0) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", "));
993: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%8.4g", (double)xnorms[f * 2 + 0]));
994: }
995: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ") |x|_1 ("));
996: for (f = 0; f < Nf; ++f) {
997: if (f > 0) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ", "));
998: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%8.4g", (double)xnorms[f * 2 + 1]));
999: }
1000: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), ") %s\n", ftable ? ftable : ""));
1001: PetscCall(PetscFree(ftable));
1002: }
1003: PetscCall(PetscFree(xnorms));
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: int main(int argc, char **argv)
1008: {
1009: MPI_Comm comm;
1010: TS ts;
1011: DM dm;
1012: Vec u;
1013: AppCtx user;
1014: PetscReal t0, t = 0.0;
1015: void *ctxs[2] = {&t, &t};
1017: PetscFunctionBeginUser;
1018: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1019: comm = PETSC_COMM_WORLD;
1020: user.functionalRegistry = NULL;
1021: globalUser = &user;
1022: PetscCall(ProcessOptions(comm, &user));
1023: PetscCall(TSCreate(comm, &ts));
1024: PetscCall(TSSetType(ts, TSBEULER));
1025: PetscCall(CreateDM(comm, &user, &dm));
1026: PetscCall(TSSetDM(ts, dm));
1027: PetscCall(ProcessMonitorOptions(comm, &user));
1029: PetscCall(DMCreateGlobalVector(dm, &u));
1030: PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
1031: if (user.useFV) {
1032: PetscBool isImplicit = PETSC_FALSE;
1034: PetscCall(PetscOptionsHasName(NULL, "", "-use_implicit", &isImplicit));
1035: if (isImplicit) {
1036: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
1037: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
1038: }
1039: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
1040: PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &user));
1041: } else {
1042: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
1043: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
1044: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
1045: }
1046: if (user.useFV) PetscCall(TSMonitorSet(ts, MonitorFunctionals, &user, NULL));
1047: PetscCall(TSSetMaxSteps(ts, 1));
1048: PetscCall(TSSetMaxTime(ts, 2.0));
1049: PetscCall(TSSetTimeStep(ts, 0.01));
1050: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1051: PetscCall(TSSetFromOptions(ts));
1053: PetscCall(DMProjectFunction(dm, 0.0, user.initialGuess, ctxs, INSERT_VALUES, u));
1054: if (user.useFV) PetscCall(SetInitialConditionFVM(dm, u, 1, user.initialGuess[1], ctxs[1]));
1055: PetscCall(VecViewFromOptions(u, NULL, "-init_vec_view"));
1056: PetscCall(TSGetTime(ts, &t));
1057: t0 = t;
1058: PetscCall(DMTSCheckFromOptions(ts, u));
1059: PetscCall(TSSolve(ts, u));
1060: PetscCall(TSGetTime(ts, &t));
1061: if (t > t0) PetscCall(DMTSCheckFromOptions(ts, u));
1062: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
1063: {
1064: PetscReal ftime;
1065: PetscInt nsteps;
1066: TSConvergedReason reason;
1068: PetscCall(TSGetSolveTime(ts, &ftime));
1069: PetscCall(TSGetStepNumber(ts, &nsteps));
1070: PetscCall(TSGetConvergedReason(ts, &reason));
1071: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
1072: }
1074: PetscCall(VecDestroy(&u));
1075: PetscCall(DMDestroy(&dm));
1076: PetscCall(TSDestroy(&ts));
1077: PetscCall(PetscFree(user.monitorFuncs));
1078: PetscCall(FunctionalDestroy(&user.functionalRegistry));
1079: PetscCall(PetscFinalize());
1080: return 0;
1081: }
1083: /*TEST
1085: testset:
1086: args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3
1088: # 2D harmonic velocity, no porosity
1089: test:
1090: suffix: p1p1
1091: requires: !complex !single
1092: 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
1093: test:
1094: suffix: p1p1_xper
1095: requires: !complex !single
1096: 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
1097: test:
1098: suffix: p1p1_xper_ref
1099: requires: !complex !single
1100: 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
1101: test:
1102: suffix: p1p1_xyper
1103: requires: !complex !single
1104: 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
1105: test:
1106: suffix: p1p1_xyper_ref
1107: requires: !complex !single
1108: 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
1109: test:
1110: suffix: p2p1
1111: requires: !complex !single
1112: 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
1113: test:
1114: suffix: p2p1_xyper
1115: requires: !complex !single
1116: 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
1118: test:
1119: suffix: adv_1
1120: requires: !complex !single
1121: args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_time_step 0.993392 -bc_inflow 1,2,4 -bc_outflow 3 -ts_view -dm_view
1123: test:
1124: suffix: adv_2
1125: requires: !complex
1126: TODO: broken memory corruption
1127: args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_time_step 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
1129: test:
1130: suffix: adv_3
1131: requires: !complex
1132: TODO: broken memory corruption
1133: 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_time_step 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
1135: test:
1136: suffix: adv_3_ex
1137: requires: !complex
1138: 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_time_step 0.1 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -ksp_max_it 100 -ts_view -dm_view
1140: test:
1141: suffix: adv_4
1142: requires: !complex
1143: TODO: broken memory corruption
1144: args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_time_step 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
1146: # 2D Advection, box, delta
1147: test:
1148: suffix: adv_delta_yper_0
1149: requires: !complex
1150: TODO: broken
1151: 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_time_step 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error
1153: test:
1154: suffix: adv_delta_yper_1
1155: requires: !complex
1156: TODO: broken
1157: 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_time_step 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 1 -source_loc 0.416666,0.416666
1159: test:
1160: suffix: adv_delta_yper_2
1161: requires: !complex
1162: TODO: broken
1163: 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_time_step 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 2 -source_loc 0.458333,0.458333
1165: test:
1166: suffix: adv_delta_yper_fim_0
1167: requires: !complex
1168: TODO: broken
1169: 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_time_step 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
1171: test:
1172: suffix: adv_delta_yper_fim_1
1173: requires: !complex
1174: TODO: broken
1175: 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_time_step 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
1177: test:
1178: suffix: adv_delta_yper_fim_2
1179: requires: !complex
1180: TODO: broken
1181: 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_time_step 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
1183: test:
1184: suffix: adv_delta_yper_im_0
1185: requires: !complex
1186: TODO: broken
1187: 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_time_step 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
1189: test:
1190: suffix: adv_delta_yper_im_1
1191: requires: !complex
1192: TODO: broken
1193: 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_time_step 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
1195: test:
1196: suffix: adv_delta_yper_im_2
1197: requires: !complex
1198: TODO: broken
1199: 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_time_step 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
1201: test:
1202: suffix: adv_delta_yper_im_3
1203: requires: !complex
1204: TODO: broken
1205: 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_time_step 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
1207: # I believe the nullspace is sin(pi y)
1208: test:
1209: suffix: adv_delta_yper_im_4
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 40 -ts_time_step 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
1214: test:
1215: suffix: adv_delta_yper_im_5
1216: requires: !complex
1217: TODO: broken
1218: 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_time_step 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
1220: test:
1221: suffix: adv_delta_yper_im_6
1222: requires: !complex
1223: TODO: broken
1224: 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_time_step 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
1225: # 2D Advection, magma benchmark 1
1227: test:
1228: suffix: adv_delta_shear_im_0
1229: requires: !complex
1230: TODO: broken
1231: 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_time_step 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
1232: # 2D Advection, box, gaussian
1234: test:
1235: suffix: adv_gauss
1236: requires: !complex
1237: TODO: broken
1238: 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_time_step 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view
1240: test:
1241: suffix: adv_gauss_im
1242: requires: !complex
1243: TODO: broken
1244: 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_time_step 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
1246: test:
1247: suffix: adv_gauss_im_1
1248: requires: !complex
1249: TODO: broken
1250: 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_time_step 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
1252: test:
1253: suffix: adv_gauss_im_2
1254: requires: !complex
1255: TODO: broken
1256: 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_time_step 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
1258: test:
1259: suffix: adv_gauss_corner
1260: requires: !complex
1261: TODO: broken
1262: 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_time_step 0.01 -bc_inflow 1 -bc_outflow 2 -ts_view -dm_view
1264: # 2D Advection+Harmonic 12-
1265: test:
1266: suffix: adv_harm_0
1267: requires: !complex
1268: TODO: broken memory corruption
1269: 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_time_step 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
1271: # Must check that FV BCs propagate to coarse meshes
1272: # Must check that FV BC ids propagate to coarse meshes
1273: # Must check that FE+FV BCs work at the same time
1274: # 2D Advection, matching wind in ex11 8-11
1275: # NOTE implicit solves are limited by accuracy of FD Jacobian
1276: test:
1277: suffix: adv_0
1278: requires: !complex !single exodusii
1279: 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_time_step 0.993392 -ts_view -dm_view
1281: test:
1282: suffix: adv_0_im
1283: requires: !complex exodusii
1284: TODO: broken memory corruption
1285: 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_time_step 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu
1287: test:
1288: suffix: adv_0_im_2
1289: requires: !complex exodusii
1290: TODO: broken
1291: 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_time_step 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
1293: test:
1294: suffix: adv_0_im_3
1295: requires: !complex exodusii
1296: TODO: broken
1297: 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_time_step 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
1299: test:
1300: suffix: adv_0_im_4
1301: requires: !complex exodusii
1302: TODO: broken
1303: 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_time_step 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
1304: # 2D Advection, misc
1306: TEST*/