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