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: }
163: static PetscErrorCode FunctionalDestroy(Functional *link)
164: {
165: Functional next, l;
167: PetscFunctionBeginUser;
168: if (!link) PetscFunctionReturn(PETSC_SUCCESS);
169: l = *link;
170: *link = NULL;
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:
648: exactFuncs[0] = quadratic_u_2d;
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: }
736: PetscCall(PetscDSSetResidual(prob, 1, f0_advection, f1_advection));
737: PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_adv_pp, g1_adv_pp, NULL, NULL));
738: PetscCall(PetscDSSetJacobian(prob, 1, 0, g0_adv_pu, g1_adv_pu, NULL, NULL));
739: if (user->velocityDist == VEL_ZERO) PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_advection));
740: else PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_coupled_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;
749: PetscQuadrature q;
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]));
762: PetscCall(PetscFECopyQuadrature(fe[0], 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));
770: PetscCall(PetscFEGetQuadrature(fe[0], &q));
771: PetscCall(PetscFVSetQuadrature(fv, q));
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));
833: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
834: PetscCall(VecGetArray(X, &x));
835: for (c = cStart; c < cEnd; ++c) {
836: PetscFVCellGeom *cg;
837: PetscScalar *xc;
839: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
840: PetscCall(DMPlexPointGlobalFieldRef(dm, c, field, x, &xc));
841: if (xc) PetscCall((*func)(dim, 0.0, cg->centroid, Nf, xc, ctx));
842: }
843: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
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));
868: PetscCall(VecGetArrayRead(X, &x));
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: }
881: PetscCall(VecRestoreArrayRead(X, &x));
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));
918: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
919: PetscCall(VecGetArrayRead(X, &x));
920: for (c = cStart; c < cEnd; ++c) {
921: PetscFVCellGeom *cg;
922: PetscScalar *cx;
924: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
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: }
942: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
943: PetscCall(VecRestoreArrayRead(X, &x));
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:
1126: suffix: adv_1
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:
1131: suffix: adv_2
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:
1137: suffix: adv_3
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:
1143: suffix: adv_3_ex
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:
1148: suffix: adv_4
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:
1155: suffix: adv_delta_yper_0
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:
1161: suffix: adv_delta_yper_1
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:
1167: suffix: adv_delta_yper_2
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:
1173: suffix: adv_delta_yper_fim_0
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:
1179: suffix: adv_delta_yper_fim_1
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:
1185: suffix: adv_delta_yper_fim_2
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:
1191: suffix: adv_delta_yper_im_0
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:
1197: suffix: adv_delta_yper_im_1
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:
1203: suffix: adv_delta_yper_im_2
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:
1209: suffix: adv_delta_yper_im_3
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:
1216: suffix: adv_delta_yper_im_4
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:
1222: suffix: adv_delta_yper_im_5
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:
1228: suffix: adv_delta_yper_im_6
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:
1235: suffix: adv_delta_shear_im_0
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:
1242: suffix: adv_gauss
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:
1248: suffix: adv_gauss_im
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:
1254: suffix: adv_gauss_im_1
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:
1260: suffix: adv_gauss_im_2
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:
1266: suffix: adv_gauss_corner
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
1271: # 2D Advection+Harmonic 12-
1272: test:
1273: suffix: adv_harm_0
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:
1284: suffix: adv_0
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:
1295: suffix: adv_0_im_2
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:
1307: suffix: adv_0_im_4
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
1311: # 2D Advection, misc
1313: TEST*/