Actual source code: ex76.c
1: static char help[] = "Time-dependent Low Mach Flow in 2d and 3d channels with finite elements.\n\
2: We solve the Low Mach flow problem for both conducting and non-conducting fluids,\n\
3: using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
5: /*F
6: The non-conducting Low Mach flow is time-dependent isoviscous Navier-Stokes flow. We discretize using the
7: finite element method on an unstructured mesh. The weak form equations are
9: \begin{align*}
10: < q, \nabla\cdot u > = 0
11: <v, du/dt> + <v, u \cdot \nabla u> + < \nabla v, \nu (\nabla u + {\nabla u}^T) > - < \nabla\cdot v, p > - < v, f > = 0
12: < w, u \cdot \nabla T > + < \nabla w, \alpha \nabla T > - < w, Q > = 0
13: \end{align*}
15: where $\nu$ is the kinematic viscosity and $\alpha$ is thermal diffusivity.
17: The conducting form is given in the ABLATE documentation [1,2] and derived in Principe and Codina [2].
19: For visualization, use
21: -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append
23: To look at nonlinear solver convergence, use
25: -dm_refine <k> -ts_max_steps 1 \
26: -ts_view -ts_monitor -snes_monitor -snes_converged_reason -ksp_converged_reason -fieldsplit_pressure_ksp_converged_reason
28: [1] https://ubchrest.github.io/ablate/content/formulations/lowMachFlow/
29: [2] https://github.com/UBCHREST/ablate/blob/main/ablateCore/flow/lowMachFlow.c
30: [3] J. Principe and R. Codina, "Mathematical models for thermally coupled low speed flows", Adv. in Theo. and App. Mech., 2(1), pp.93--112, 2009.
31: F*/
33: #include <petscdmplex.h>
34: #include <petscsnes.h>
35: #include <petscts.h>
36: #include <petscds.h>
37: #include <petscbag.h>
39: typedef enum {
40: MOD_INCOMPRESSIBLE,
41: MOD_CONDUCTING,
42: NUM_MOD_TYPES
43: } ModType;
44: const char *modTypes[NUM_MOD_TYPES + 1] = {"incompressible", "conducting", "unknown"};
46: typedef enum {
47: SOL_QUADRATIC,
48: SOL_CUBIC,
49: SOL_CUBIC_TRIG,
50: SOL_TAYLOR_GREEN,
51: SOL_PIPE,
52: SOL_PIPE_WIGGLY,
53: NUM_SOL_TYPES
54: } SolType;
55: const char *solTypes[NUM_SOL_TYPES + 1] = {"quadratic", "cubic", "cubic_trig", "taylor_green", "pipe", "pipe_wiggly", "unknown"};
57: /* Fields */
58: const PetscInt VEL = 0;
59: const PetscInt PRES = 1;
60: const PetscInt TEMP = 2;
61: /* Sources */
62: const PetscInt MOMENTUM = 0;
63: const PetscInt MASS = 1;
64: const PetscInt ENERGY = 2;
65: /* Constants */
66: const PetscInt STROUHAL = 0;
67: const PetscInt FROUDE = 1;
68: const PetscInt REYNOLDS = 2;
69: const PetscInt PECLET = 3;
70: const PetscInt P_TH = 4;
71: const PetscInt MU = 5;
72: const PetscInt NU = 6;
73: const PetscInt C_P = 7;
74: const PetscInt K = 8;
75: const PetscInt ALPHA = 9;
76: const PetscInt T_IN = 10;
77: const PetscInt G_DIR = 11;
78: const PetscInt EPSILON = 12;
80: typedef struct {
81: PetscReal Strouhal; /* Strouhal number */
82: PetscReal Froude; /* Froude number */
83: PetscReal Reynolds; /* Reynolds number */
84: PetscReal Peclet; /* Peclet number */
85: PetscReal p_th; /* Thermodynamic pressure */
86: PetscReal mu; /* Dynamic viscosity */
87: PetscReal nu; /* Kinematic viscosity */
88: PetscReal c_p; /* Specific heat at constant pressure */
89: PetscReal k; /* Thermal conductivity */
90: PetscReal alpha; /* Thermal diffusivity */
91: PetscReal T_in; /* Inlet temperature */
92: PetscReal g_dir; /* Gravity direction */
93: PetscReal epsilon; /* Strength of perturbation */
94: } Parameter;
96: typedef struct {
97: /* Problem definition */
98: PetscBag bag; /* Holds problem parameters */
99: ModType modType; /* Model type */
100: SolType solType; /* MMS solution type */
101: PetscBool hasNullSpace; /* Problem has the constant null space for pressure */
102: /* Flow diagnostics */
103: DM dmCell; /* A DM with piecewise constant discretization */
104: } AppCtx;
106: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
107: {
108: for (PetscInt d = 0; d < Nc; ++d) u[d] = 0.0;
109: return PETSC_SUCCESS;
110: }
112: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
113: {
114: for (PetscInt d = 0; d < Nc; ++d) u[d] = 1.0;
115: return PETSC_SUCCESS;
116: }
118: /*
119: CASE: quadratic
120: In 2D we use exact solution:
122: u = t + x^2 + y^2
123: v = t + 2x^2 - 2xy
124: p = x + y - 1
125: T = t + x + y + 1
126: f = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2 -4\nu + 2, t (2x - 2y) + 4xy^2 + 2x^2y - 2y^3 -4\nu + 2>
127: Q = 1 + 2t + 3x^2 - 2xy + y^2
129: so that
131: \nabla \cdot u = 2x - 2x = 0
133: f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
134: = <1, 1> + <t + x^2 + y^2, t + 2x^2 - 2xy> . <<2x, 4x - 2y>, <2y, -2x>> - \nu <4, 4> + <1, 1>
135: = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t (2x - 2y) + 2x^2y + 4xy^2 - 2y^3> + <-4 \nu + 2, -4\nu + 2>
136: = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2 - 4\nu + 2, t (2x - 2y) + 4xy^2 + 2x^2y - 2y^3 - 4\nu + 2>
138: Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
139: = 1 + <t + x^2 + y^2, t + 2x^2 - 2xy> . <1, 1> - \alpha 0
140: = 1 + 2t + 3x^2 - 2xy + y^2
141: */
143: static PetscErrorCode quadratic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
144: {
145: u[0] = time + X[0] * X[0] + X[1] * X[1];
146: u[1] = time + 2.0 * X[0] * X[0] - 2.0 * X[0] * X[1];
147: return PETSC_SUCCESS;
148: }
149: static PetscErrorCode quadratic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
150: {
151: u[0] = 1.0;
152: u[1] = 1.0;
153: return PETSC_SUCCESS;
154: }
156: static PetscErrorCode quadratic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
157: {
158: p[0] = X[0] + X[1] - 1.0;
159: return PETSC_SUCCESS;
160: }
162: static PetscErrorCode quadratic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
163: {
164: T[0] = time + X[0] + X[1] + 1.0;
165: return PETSC_SUCCESS;
166: }
167: static PetscErrorCode quadratic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
168: {
169: T[0] = 1.0;
170: return PETSC_SUCCESS;
171: }
173: static void f0_quadratic_v(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[])
174: {
175: const PetscReal nu = PetscRealPart(constants[NU]);
177: f0[0] -= t * (2 * X[0] + 2 * X[1]) + 2 * X[0] * X[0] * X[0] + 4 * X[0] * X[0] * X[1] - 2 * X[0] * X[1] * X[1] - 4.0 * nu + 2;
178: f0[1] -= t * (2 * X[0] - 2 * X[1]) + 4 * X[0] * X[1] * X[1] + 2 * X[0] * X[0] * X[1] - 2 * X[1] * X[1] * X[1] - 4.0 * nu + 2;
179: }
181: static void f0_quadratic_w(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[])
182: {
183: f0[0] -= 2 * t + 1 + 3 * X[0] * X[0] - 2 * X[0] * X[1] + X[1] * X[1];
184: }
186: /*
187: CASE: quadratic
188: In 2D we use exact solution:
190: u = t + x^2 + y^2
191: v = t + 2x^2 - 2xy
192: p = x + y - 1
193: T = t + x + y + 1
194: rho = p^{th} / T
196: so that
198: \nabla \cdot u = 2x - 2x = 0
199: grad u = <<2 x, 4x - 2y>, <2 y, -2x>>
200: epsilon(u) = 1/2 (grad u + grad u^T) = <<2x, 2x>, <2x, -2x>>
201: epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
202: div epsilon'(u) = <2, 2>
204: f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
205: = rho S <1, 1> + rho <t + x^2 + y^2, t + 2x^2 - 2xy> . <<2x, 4x - 2y>, <2y, -2x>> - 2\mu/Re <2, 2> + <1, 1> + rho/F^2 <0, 1>
206: = rho S <1, 1> + rho <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t (2x - 2y) + 2x^2y + 4xy^2 - 2y^3> - mu/Re <4, 4> + <1, 1> + rho/F^2 <0, 1>
208: g = S rho_t + div (rho u)
209: = -S pth T_t/T^2 + rho div (u) + u . grad rho
210: = -S pth 1/T^2 - pth u . grad T / T^2
211: = -pth / T^2 (S + 2t + 3 x^2 - 2xy + y^2)
213: Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
214: = c_p S pth / T + c_p pth (2t + 3 x^2 - 2xy + y^2) / T - k/Pe 0
215: = c_p pth / T (S + 2t + 3 x^2 - 2xy + y^2)
216: */
217: static void f0_conduct_quadratic_v(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[])
218: {
219: const PetscReal S = PetscRealPart(constants[STROUHAL]);
220: const PetscReal F = PetscRealPart(constants[FROUDE]);
221: const PetscReal Re = PetscRealPart(constants[REYNOLDS]);
222: const PetscReal mu = PetscRealPart(constants[MU]);
223: const PetscReal p_th = PetscRealPart(constants[P_TH]);
224: const PetscReal rho = p_th / (t + X[0] + X[1] + 1.);
225: const PetscInt gd = (PetscInt)PetscRealPart(constants[G_DIR]);
227: f0[0] -= rho * S + rho * (2. * t * (X[0] + X[1]) + 2. * X[0] * X[0] * X[0] + 4. * X[0] * X[0] * X[1] - 2. * X[0] * X[1] * X[1]) - 4. * mu / Re + 1.;
228: f0[1] -= rho * S + rho * (2. * t * (X[0] - X[1]) + 2. * X[0] * X[0] * X[1] + 4. * X[0] * X[1] * X[1] - 2. * X[1] * X[1] * X[1]) - 4. * mu / Re + 1.;
229: f0[gd] -= rho / PetscSqr(F);
230: }
232: static void f0_conduct_quadratic_q(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: const PetscReal S = PetscRealPart(constants[STROUHAL]);
235: const PetscReal p_th = PetscRealPart(constants[P_TH]);
237: f0[0] += p_th * (S + 2. * t + 3. * X[0] * X[0] - 2. * X[0] * X[1] + X[1] * X[1]) / PetscSqr(t + X[0] + X[1] + 1.);
238: }
240: static void f0_conduct_quadratic_w(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[])
241: {
242: const PetscReal S = PetscRealPart(constants[STROUHAL]);
243: const PetscReal c_p = PetscRealPart(constants[C_P]);
244: const PetscReal p_th = PetscRealPart(constants[P_TH]);
246: f0[0] -= c_p * p_th * (S + 2. * t + 3. * X[0] * X[0] - 2. * X[0] * X[1] + X[1] * X[1]) / (t + X[0] + X[1] + 1.);
247: }
249: /*
250: CASE: cubic
251: In 2D we use exact solution:
253: u = t + x^3 + y^3
254: v = t + 2x^3 - 3x^2y
255: p = 3/2 x^2 + 3/2 y^2 - 1
256: T = t + 1/2 x^2 + 1/2 y^2
257: f = < t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x + 1,
258: t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y + 1>
259: Q = x^4 + xy^3 + 2x^3y - 3x^2y^2 + xt + yt - 2\alpha + 1
261: so that
263: \nabla \cdot u = 3x^2 - 3x^2 = 0
265: du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p - f
266: = <1,1> + <t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3, t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4> - \nu<6x + 6y, 12x - 6y> + <3x, 3y> - <t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x + 1, t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y + 1> = 0
268: dT/dt + u \cdot \nabla T - \alpha \Delta T - Q = 1 + (x^3 + y^3) x + (2x^3 - 3x^2y) y - 2*\alpha - (x^4 + xy^3 + 2x^3y - 3x^2y^2 - 2*\alpha +1) = 0
269: */
270: static PetscErrorCode cubic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
271: {
272: u[0] = time + X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
273: u[1] = time + 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
274: return PETSC_SUCCESS;
275: }
276: static PetscErrorCode cubic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
277: {
278: u[0] = 1.0;
279: u[1] = 1.0;
280: return PETSC_SUCCESS;
281: }
283: static PetscErrorCode cubic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
284: {
285: p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
286: return PETSC_SUCCESS;
287: }
289: static PetscErrorCode cubic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
290: {
291: T[0] = time + X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
292: return PETSC_SUCCESS;
293: }
294: static PetscErrorCode cubic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
295: {
296: T[0] = 1.0;
297: return PETSC_SUCCESS;
298: }
300: static void f0_cubic_v(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[])
301: {
302: const PetscReal nu = PetscRealPart(constants[NU]);
304: f0[0] -= (t * (3 * X[0] * X[0] + 3 * X[1] * X[1]) + 3 * X[0] * X[0] * X[0] * X[0] * X[0] + 6 * X[0] * X[0] * X[0] * X[1] * X[1] - 6 * X[0] * X[0] * X[1] * X[1] * X[1] - (6 * X[0] + 6 * X[1]) * nu + 3 * X[0] + 1);
305: f0[1] -= (t * (3 * X[0] * X[0] - 6 * X[0] * X[1]) + 3 * X[0] * X[0] * X[0] * X[0] * X[1] + 6 * X[0] * X[0] * X[1] * X[1] * X[1] - 6 * X[0] * X[1] * X[1] * X[1] * X[1] - (12 * X[0] - 6 * X[1]) * nu + 3 * X[1] + 1);
306: }
308: static void f0_cubic_w(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[])
309: {
310: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
312: f0[0] -= X[0] * X[0] * X[0] * X[0] + 2.0 * X[0] * X[0] * X[0] * X[1] - 3.0 * X[0] * X[0] * X[1] * X[1] + X[0] * X[1] * X[1] * X[1] + X[0] * t + X[1] * t - 2.0 * alpha + 1;
313: }
315: /*
316: CASE: cubic-trigonometric
317: In 2D we use exact solution:
319: u = beta cos t + x^3 + y^3
320: v = beta sin t + 2x^3 - 3x^2y
321: p = 3/2 x^2 + 3/2 y^2 - 1
322: T = 20 cos t + 1/2 x^2 + 1/2 y^2
323: f = < beta cos t 3x^2 + beta sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x,
324: beta cos t (6x^2 - 6xy) - beta sin t (3x^2) + 3x^4y + 6x^2y^3 - 6xy^4 - \nu(12x - 6y) + 3y>
325: Q = beta cos t x + beta sin t (y - 1) + x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2\alpha
327: so that
329: \nabla \cdot u = 3x^2 - 3x^2 = 0
331: f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
332: = <-sin t, cos t> + <cos t + x^3 + y^3, sin t + 2x^3 - 3x^2y> <<3x^2, 6x^2 - 6xy>, <3y^2, -3x^2>> - \nu <6x + 6y, 12x - 6y> + <3x, 3y>
333: = <-sin t, cos t> + <cos t 3x^2 + 3x^5 + 3x^2y^3 + sin t 3y^2 + 6x^3y^2 - 9x^2y^3, cos t (6x^2 - 6xy) + 6x^5 - 6x^4y + 6x^2y^3 - 6xy^4 + sin t (-3x^2) - 6x^5 + 9x^4y> - \nu <6x + 6y, 12x - 6y> + <3x, 3y>
334: = <cos t (3x^2) + sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu (6x + 6y) + 3x,
335: cos t (6x^2 - 6xy) - sin t (3x^2) + 3x^4y + 6x^2y^3 - 6xy^4 - \nu (12x - 6y) + 3y>
337: Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
338: = -sin t + <cos t + x^3 + y^3, sin t + 2x^3 - 3x^2y> . <x, y> - 2 \alpha
339: = -sin t + cos t (x) + x^4 + xy^3 + sin t (y) + 2x^3y - 3x^2y^2 - 2 \alpha
340: = cos t x + sin t (y - 1) + (x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2 \alpha)
341: */
342: static PetscErrorCode cubic_trig_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
343: {
344: u[0] = 100. * PetscCosReal(time) + X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
345: u[1] = 100. * PetscSinReal(time) + 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
346: return PETSC_SUCCESS;
347: }
348: static PetscErrorCode cubic_trig_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
349: {
350: u[0] = -100. * PetscSinReal(time);
351: u[1] = 100. * PetscCosReal(time);
352: return PETSC_SUCCESS;
353: }
355: static PetscErrorCode cubic_trig_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
356: {
357: p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
358: return PETSC_SUCCESS;
359: }
361: static PetscErrorCode cubic_trig_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
362: {
363: T[0] = 100. * PetscCosReal(time) + X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
364: return PETSC_SUCCESS;
365: }
366: static PetscErrorCode cubic_trig_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
367: {
368: T[0] = -100. * PetscSinReal(time);
369: return PETSC_SUCCESS;
370: }
372: static void f0_cubic_trig_v(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[])
373: {
374: const PetscReal nu = PetscRealPart(constants[NU]);
376: f0[0] -= 100. * PetscCosReal(t) * (3 * X[0] * X[0]) + 100. * PetscSinReal(t) * (3 * X[1] * X[1] - 1.) + 3 * X[0] * X[0] * X[0] * X[0] * X[0] + 6 * X[0] * X[0] * X[0] * X[1] * X[1] - 6 * X[0] * X[0] * X[1] * X[1] * X[1] - (6 * X[0] + 6 * X[1]) * nu + 3 * X[0];
377: f0[1] -= 100. * PetscCosReal(t) * (6 * X[0] * X[0] - 6 * X[0] * X[1]) - 100. * PetscSinReal(t) * (3 * X[0] * X[0]) + 3 * X[0] * X[0] * X[0] * X[0] * X[1] + 6 * X[0] * X[0] * X[1] * X[1] * X[1] - 6 * X[0] * X[1] * X[1] * X[1] * X[1] - (12 * X[0] - 6 * X[1]) * nu + 3 * X[1];
378: }
380: static void f0_cubic_trig_w(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[])
381: {
382: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
384: f0[0] -= 100. * PetscCosReal(t) * X[0] + 100. * PetscSinReal(t) * (X[1] - 1.) + X[0] * X[0] * X[0] * X[0] + 2.0 * X[0] * X[0] * X[0] * X[1] - 3.0 * X[0] * X[0] * X[1] * X[1] + X[0] * X[1] * X[1] * X[1] - 2.0 * alpha;
385: }
387: /*
388: CASE: Taylor-Green vortex
389: In 2D we use exact solution:
391: u = 1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)
392: v = 1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)
393: p = -1/4 [cos(2 \pi(x - t)) + cos(2 \pi(y - t))] exp(-4 \pi^2 \nu t)
394: T = t + x + y
395: f = <\nu \pi^2 exp(-2\nu \pi^2 t) cos(\pi(x-t)) sin(\pi(y-t)), -\nu \pi^2 exp(-2\nu \pi^2 t) sin(\pi(x-t)) cos(\pi(y-t)) >
396: Q = 3 + sin(\pi(x-y)) exp(-2\nu \pi^2 t)
398: so that
400: \nabla \cdot u = \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t) - \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t) = 0
402: f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
403: = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi cos(\pi(x - t)) sin(\pi(y - t))) exp(-2 \pi^2 \nu t),
404: \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi sin(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)>
405: + < \pi (1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)) sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
406: \pi (1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)) cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
407: + <-\pi (1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)) cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
408: -\pi (1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)) sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
409: + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
410: 2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
411: + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
412: \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
413: = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi cos(\pi(x - t)) sin(\pi(y - t))) exp(-2 \pi^2 \nu t),
414: \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t)) - 2\pi sin(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)>
415: + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
416: \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
417: + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
418: -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
419: + <-\pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
420: -\pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
421: + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
422: 2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
423: + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
424: \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
425: = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t),
426: \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)>
427: + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
428: \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
429: + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
430: -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
431: = < \pi cos(\pi(x - t)) cos(\pi(y - t)),
432: \pi sin(\pi(x - t)) sin(\pi(y - t))>
433: + <-\pi cos(\pi(x - t)) cos(\pi(y - t)),
434: -\pi sin(\pi(x - t)) sin(\pi(y - t))> = 0
435: Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
436: = 1 + u \cdot <1, 1> - 0
437: = 1 + u + v
438: */
440: static PetscErrorCode taylor_green_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
441: {
442: u[0] = 1 - PetscCosReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time)) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
443: u[1] = 1 + PetscSinReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time)) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
444: return PETSC_SUCCESS;
445: }
446: static PetscErrorCode taylor_green_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
447: {
448: u[0] = -PETSC_PI * (PetscSinReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time)) - PetscCosReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time)) - 2 * PETSC_PI * PetscCosReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time))) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
449: u[1] = PETSC_PI * (PetscSinReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time)) - PetscCosReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time)) - 2 * PETSC_PI * PetscSinReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time))) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
450: return PETSC_SUCCESS;
451: }
453: static PetscErrorCode taylor_green_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
454: {
455: p[0] = -0.25 * (PetscCosReal(2 * PETSC_PI * (X[0] - time)) + PetscCosReal(2 * PETSC_PI * (X[1] - time))) * PetscExpReal(-4 * PETSC_PI * PETSC_PI * time);
456: return PETSC_SUCCESS;
457: }
459: static PetscErrorCode taylor_green_p_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
460: {
461: p[0] = PETSC_PI * (0.5 * (PetscSinReal(2 * PETSC_PI * (X[0] - time)) + PetscSinReal(2 * PETSC_PI * (X[1] - time))) + PETSC_PI * (PetscCosReal(2 * PETSC_PI * (X[0] - time)) + PetscCosReal(2 * PETSC_PI * (X[1] - time)))) * PetscExpReal(-4 * PETSC_PI * PETSC_PI * time);
462: return PETSC_SUCCESS;
463: }
465: static PetscErrorCode taylor_green_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
466: {
467: T[0] = time + X[0] + X[1];
468: return PETSC_SUCCESS;
469: }
470: static PetscErrorCode taylor_green_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
471: {
472: T[0] = 1.0;
473: return PETSC_SUCCESS;
474: }
476: static void f0_taylor_green_w(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[])
477: {
478: PetscScalar vel[2];
480: PetscCallAbort(PETSC_COMM_SELF, taylor_green_u(dim, t, X, Nf, vel, NULL));
481: f0[0] -= 1.0 + vel[0] + vel[1];
482: }
484: /*
485: CASE: Pipe flow
486: Poiseuille flow, with the incoming fluid having a parabolic temperature profile and the side walls being held at T_in
488: u = \Delta Re/(2 mu) y (1 - y)
489: v = 0
490: p = -\Delta x
491: T = y (1 - y) + T_in
492: rho = p^{th} / T
494: so that
496: \nabla \cdot u = 0 - 0 = 0
497: grad u = \Delta Re/(2 mu) <<0, 0>, <1 - 2y, 0>>
498: epsilon(u) = 1/2 (grad u + grad u^T) = \Delta Re/(4 mu) <<0, 1 - 2y>, <<1 - 2y, 0>>
499: epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
500: div epsilon'(u) = -\Delta Re/(2 mu) <1, 0>
502: f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
503: = 0 + 0 - div (2\mu/Re \epsilon'(u) - pI) + rho / F^2 \hat y
504: = -\Delta div <<x, (1 - 2y)/2>, <<(1 - 2y)/2, x>> + rho / F^2 \hat y
505: = \Delta <1, 0> - \Delta <1, 0> + rho/F^2 <0, 1>
506: = rho/F^2 <0, 1>
508: g = S rho_t + div (rho u)
509: = 0 + rho div (u) + u . grad rho
510: = 0 + 0 - pth u . grad T / T^2
511: = 0
513: Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
514: = 0 + c_p pth / T 0 + 2 k/Pe
515: = 2 k/Pe
517: The boundary conditions on the top and bottom are zero velocity and T_in temperature. The boundary term is
519: (2\mu/Re \epsilon'(u) - p I) . n = \Delta <<x, (1 - 2y)/2>, <<(1 - 2y)/2, x>> . n
521: so that
523: x = 0: \Delta <<0, (1 - 2y)/2>, <<(1 - 2y)/2, 0>> . <-1, 0> = <0, (2y - 1)/2>
524: x = 1: \Delta <<1, (1 - 2y)/2>, <<(1 - 2y)/2, 1>> . <1, 0> = <1, (1 - 2y)/2>
525: */
526: static PetscErrorCode pipe_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
527: {
528: Parameter *param = (Parameter *)ctx;
530: u[0] = (0.5 * param->Reynolds / param->mu) * X[1] * (1.0 - X[1]);
531: u[1] = 0.0;
532: return PETSC_SUCCESS;
533: }
534: static PetscErrorCode pipe_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
535: {
536: u[0] = 0.0;
537: u[1] = 0.0;
538: return PETSC_SUCCESS;
539: }
541: static PetscErrorCode pipe_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
542: {
543: p[0] = -X[0];
544: return PETSC_SUCCESS;
545: }
546: static PetscErrorCode pipe_p_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
547: {
548: p[0] = 0.0;
549: return PETSC_SUCCESS;
550: }
552: static PetscErrorCode pipe_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
553: {
554: Parameter *param = (Parameter *)ctx;
556: T[0] = X[1] * (1.0 - X[1]) + param->T_in;
557: return PETSC_SUCCESS;
558: }
559: static PetscErrorCode pipe_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
560: {
561: T[0] = 0.0;
562: return PETSC_SUCCESS;
563: }
565: static void f0_conduct_pipe_v(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[])
566: {
567: const PetscReal F = PetscRealPart(constants[FROUDE]);
568: const PetscReal p_th = PetscRealPart(constants[P_TH]);
569: const PetscReal T_in = PetscRealPart(constants[T_IN]);
570: const PetscReal rho = p_th / (X[1] * (1. - X[1]) + T_in);
571: const PetscInt gd = (PetscInt)PetscRealPart(constants[G_DIR]);
573: f0[gd] -= rho / PetscSqr(F);
574: }
576: static void f0_conduct_bd_pipe_v(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
577: {
578: PetscReal sigma[4] = {X[0], (PetscReal)(0.5 * (1. - 2. * X[1])), (PetscReal)(0.5 * (1. - 2. * X[1])), X[0]};
579: PetscInt d, e;
581: for (d = 0; d < dim; ++d) {
582: for (e = 0; e < dim; ++e) f0[d] -= sigma[d * dim + e] * n[e];
583: }
584: }
586: static void f0_conduct_pipe_q(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[])
587: {
588: f0[0] += 0.0;
589: }
591: static void f0_conduct_pipe_w(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[])
592: {
593: const PetscReal k = PetscRealPart(constants[K]);
594: const PetscReal Pe = PetscRealPart(constants[PECLET]);
596: f0[0] -= 2 * k / Pe;
597: }
599: /*
600: CASE: Wiggly pipe flow
601: Perturbed Poiseuille flow, with the incoming fluid having a perturbed parabolic temperature profile and the side walls being held at T_in
603: u = \Delta Re/(2 mu) [y (1 - y) + a sin(pi y)]
604: v = 0
605: p = -\Delta x
606: T = y (1 - y) + a sin(pi y) + T_in
607: rho = p^{th} / T
609: so that
611: \nabla \cdot u = 0 - 0 = 0
612: grad u = \Delta Re/(2 mu) <<0, 0>, <1 - 2y + a pi cos(pi y), 0>>
613: epsilon(u) = 1/2 (grad u + grad u^T) = \Delta Re/(4 mu) <<0, 1 - 2y + a pi cos(pi y)>, <<1 - 2y + a pi cos(pi y), 0>>
614: epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
615: div epsilon'(u) = -\Delta Re/(2 mu) <1 + a pi^2/2 sin(pi y), 0>
617: f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
618: = 0 + 0 - div (2\mu/Re \epsilon'(u) - pI) + rho / F^2 \hat y
619: = -\Delta div <<x, (1 - 2y)/2 + a pi/2 cos(pi y)>, <<(1 - 2y)/2 + a pi/2 cos(pi y), x>> + rho / F^2 \hat y
620: = -\Delta <1 - 1 - a pi^2/2 sin(pi y), 0> + rho/F^2 <0, 1>
621: = a \Delta pi^2/2 sin(pi y) <1, 0> + rho/F^2 <0, 1>
623: g = S rho_t + div (rho u)
624: = 0 + rho div (u) + u . grad rho
625: = 0 + 0 - pth u . grad T / T^2
626: = 0
628: Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
629: = 0 + c_p pth / T 0 - k/Pe div <0, 1 - 2y + a pi cos(pi y)>
630: = - k/Pe (-2 - a pi^2 sin(pi y))
631: = 2 k/Pe (1 + a pi^2/2 sin(pi y))
633: The boundary conditions on the top and bottom are zero velocity and T_in temperature. The boundary term is
635: (2\mu/Re \epsilon'(u) - p I) . n = \Delta <<x, (1 - 2y)/2 + a pi/2 cos(pi y)>, <<(1 - 2y)/2 + a pi/2 cos(pi y), x>> . n
637: so that
639: x = 0: \Delta <<0, (1 - 2y)/2>, <<(1 - 2y)/2, 0>> . <-1, 0> = <0, (2y - 1)/2 - a pi/2 cos(pi y)>
640: x = 1: \Delta <<1, (1 - 2y)/2>, <<(1 - 2y)/2, 1>> . < 1, 0> = <1, (1 - 2y)/2 + a pi/2 cos(pi y)>
641: */
642: static PetscErrorCode pipe_wiggly_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
643: {
644: Parameter *param = (Parameter *)ctx;
646: u[0] = (0.5 * param->Reynolds / param->mu) * (X[1] * (1.0 - X[1]) + param->epsilon * PetscSinReal(PETSC_PI * X[1]));
647: u[1] = 0.0;
648: return PETSC_SUCCESS;
649: }
650: static PetscErrorCode pipe_wiggly_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
651: {
652: u[0] = 0.0;
653: u[1] = 0.0;
654: return PETSC_SUCCESS;
655: }
657: static PetscErrorCode pipe_wiggly_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
658: {
659: p[0] = -X[0];
660: return PETSC_SUCCESS;
661: }
662: static PetscErrorCode pipe_wiggly_p_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
663: {
664: p[0] = 0.0;
665: return PETSC_SUCCESS;
666: }
668: static PetscErrorCode pipe_wiggly_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
669: {
670: Parameter *param = (Parameter *)ctx;
672: T[0] = X[1] * (1.0 - X[1]) + param->epsilon * PetscSinReal(PETSC_PI * X[1]) + param->T_in;
673: return PETSC_SUCCESS;
674: }
675: static PetscErrorCode pipe_wiggly_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
676: {
677: T[0] = 0.0;
678: return PETSC_SUCCESS;
679: }
681: static void f0_conduct_pipe_wiggly_v(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[])
682: {
683: const PetscReal F = PetscRealPart(constants[FROUDE]);
684: const PetscReal p_th = PetscRealPart(constants[P_TH]);
685: const PetscReal T_in = PetscRealPart(constants[T_IN]);
686: const PetscReal eps = PetscRealPart(constants[EPSILON]);
687: const PetscReal rho = p_th / (X[1] * (1. - X[1]) + T_in);
688: const PetscInt gd = (PetscInt)PetscRealPart(constants[G_DIR]);
690: f0[0] -= eps * 0.5 * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * X[1]);
691: f0[gd] -= rho / PetscSqr(F);
692: }
694: static void f0_conduct_bd_pipe_wiggly_v(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
695: {
696: const PetscReal eps = PetscRealPart(constants[EPSILON]);
697: PetscReal sigma[4] = {X[0], (PetscReal)(0.5 * (1. - 2. * X[1]) + eps * 0.5 * PETSC_PI * PetscCosReal(PETSC_PI * X[1])), (PetscReal)(0.5 * (1. - 2. * X[1]) + eps * 0.5 * PETSC_PI * PetscCosReal(PETSC_PI * X[1])), X[0]};
698: PetscInt d, e;
700: for (d = 0; d < dim; ++d) {
701: for (e = 0; e < dim; ++e) f0[d] -= sigma[d * dim + e] * n[e];
702: }
703: }
705: static void f0_conduct_pipe_wiggly_q(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[])
706: {
707: f0[0] += 0.0;
708: }
710: static void f0_conduct_pipe_wiggly_w(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[])
711: {
712: const PetscReal k = PetscRealPart(constants[K]);
713: const PetscReal Pe = PetscRealPart(constants[PECLET]);
714: const PetscReal eps = PetscRealPart(constants[EPSILON]);
716: f0[0] -= 2 * k / Pe * (1.0 + eps * 0.5 * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * X[1]));
717: }
719: /* Physics Kernels */
721: static void f0_q(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[])
722: {
723: PetscInt d;
724: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
725: }
727: /* -\frac{Sp^{th}}{T^2} \frac{\partial T}{\partial t} + \frac{p^{th}}{T} \nabla \cdot \vb{u} - \frac{p^{th}}{T^2} \vb{u} \cdot \nabla T */
728: static void f0_conduct_q(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[])
729: {
730: const PetscReal S = PetscRealPart(constants[STROUHAL]);
731: const PetscReal p_th = PetscRealPart(constants[P_TH]);
732: PetscInt d;
734: // -\frac{S p^{th}}{T^2} \frac{\partial T}{\partial t}
735: f0[0] += -u_t[uOff[TEMP]] * S * p_th / PetscSqr(u[uOff[TEMP]]);
737: // \frac{p^{th}}{T} \nabla \cdot \vb{u}
738: for (d = 0; d < dim; ++d) f0[0] += p_th / u[uOff[TEMP]] * u_x[uOff_x[VEL] + d * dim + d];
740: // - \frac{p^{th}}{T^2} \vb{u} \cdot \nabla T
741: for (d = 0; d < dim; ++d) f0[0] -= p_th / (u[uOff[TEMP]] * u[uOff[TEMP]]) * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];
743: // Add in any fixed source term
744: if (NfAux > 0) f0[0] += a[aOff[MASS]];
745: }
747: /* \vb{u}_t + \vb{u} \cdot \nabla\vb{u} */
748: static void f0_v(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[])
749: {
750: const PetscInt Nc = dim;
751: PetscInt c, d;
753: for (c = 0; c < Nc; ++c) {
754: /* \vb{u}_t */
755: f0[c] += u_t[uOff[VEL] + c];
756: /* \vb{u} \cdot \nabla\vb{u} */
757: for (d = 0; d < dim; ++d) f0[c] += u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
758: }
759: }
761: /* \rho S \frac{\partial \vb{u}}{\partial t} + \rho \vb{u} \cdot \nabla \vb{u} + \rho \frac{\hat{\vb{z}}}{F^2} */
762: static void f0_conduct_v(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[])
763: {
764: const PetscReal S = PetscRealPart(constants[STROUHAL]);
765: const PetscReal F = PetscRealPart(constants[FROUDE]);
766: const PetscReal p_th = PetscRealPart(constants[P_TH]);
767: const PetscReal rho = p_th / PetscRealPart(u[uOff[TEMP]]);
768: const PetscInt gdir = (PetscInt)PetscRealPart(constants[G_DIR]);
769: PetscInt Nc = dim;
770: PetscInt c, d;
772: // \rho S \frac{\partial \vb{u}}{\partial t}
773: for (d = 0; d < dim; ++d) f0[d] = rho * S * u_t[uOff[VEL] + d];
775: // \rho \vb{u} \cdot \nabla \vb{u}
776: for (c = 0; c < Nc; ++c) {
777: for (d = 0; d < dim; ++d) f0[c] += rho * u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
778: }
780: // rho \hat{z}/F^2
781: f0[gdir] += rho / (F * F);
783: // Add in any fixed source term
784: if (NfAux > 0) {
785: for (d = 0; d < dim; ++d) f0[d] += a[aOff[MOMENTUM] + d];
786: }
787: }
789: /*f1_v = \nu[grad(u) + grad(u)^T] - pI */
790: static void f1_v(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[])
791: {
792: const PetscReal nu = PetscRealPart(constants[NU]);
793: const PetscInt Nc = dim;
794: PetscInt c, d;
796: for (c = 0; c < Nc; ++c) {
797: for (d = 0; d < dim; ++d) f1[c * dim + d] = nu * (u_x[c * dim + d] + u_x[d * dim + c]);
798: f1[c * dim + c] -= u[uOff[1]];
799: }
800: }
802: /* 2 \mu/Re (1/2 (\nabla \vb{u} + \nabla \vb{u}^T) - 1/3 (\nabla \cdot \vb{u}) I) - p I */
803: static void f1_conduct_v(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[])
804: {
805: const PetscReal Re = PetscRealPart(constants[REYNOLDS]);
806: const PetscReal mu = PetscRealPart(constants[MU]);
807: const PetscReal coef = mu / Re;
808: PetscReal u_div = 0.0;
809: const PetscInt Nc = dim;
810: PetscInt c, d;
812: for (c = 0; c < Nc; ++c) u_div += PetscRealPart(u_x[uOff_x[VEL] + c * dim + c]);
814: for (c = 0; c < Nc; ++c) {
815: // 2 \mu/Re 1/2 (\nabla \vb{u} + \nabla \vb{u}^T
816: for (d = 0; d < dim; ++d) f1[c * dim + d] += coef * (u_x[uOff_x[VEL] + c * dim + d] + u_x[uOff_x[VEL] + d * dim + c]);
817: // -2/3 \mu/Re (\nabla \cdot \vb{u}) I
818: f1[c * dim + c] -= 2.0 * coef / 3.0 * u_div;
819: }
821: // -p I
822: for (c = 0; c < Nc; ++c) f1[c * dim + c] -= u[uOff[PRES]];
823: }
825: /* T_t + \vb{u} \cdot \nabla T */
826: static void f0_w(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[])
827: {
828: PetscInt d;
830: /* T_t */
831: f0[0] += u_t[uOff[TEMP]];
832: /* \vb{u} \cdot \nabla T */
833: for (d = 0; d < dim; ++d) f0[0] += u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];
834: }
836: /* \frac{C_p S p^{th}}{T} \frac{\partial T}{\partial t} + \frac{C_p p^{th}}{T} \vb{u} \cdot \nabla T */
837: static void f0_conduct_w(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[])
838: {
839: const PetscReal S = PetscRealPart(constants[STROUHAL]);
840: const PetscReal c_p = PetscRealPart(constants[C_P]);
841: const PetscReal p_th = PetscRealPart(constants[P_TH]);
842: PetscInt d;
844: // \frac{C_p S p^{th}}{T} \frac{\partial T}{\partial t}
845: f0[0] = c_p * S * p_th / u[uOff[TEMP]] * u_t[uOff[TEMP]];
847: // \frac{C_p p^{th}}{T} \vb{u} \cdot \nabla T
848: for (d = 0; d < dim; ++d) f0[0] += c_p * p_th / u[uOff[TEMP]] * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];
850: // Add in any fixed source term
851: if (NfAux > 0) f0[0] += a[aOff[ENERGY]];
852: }
854: static void f1_w(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[])
855: {
856: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
857: PetscInt d;
859: for (d = 0; d < dim; ++d) f1[d] = alpha * u_x[uOff_x[2] + d];
860: }
862: /* \frac{k}{Pe} \nabla T */
863: static void f1_conduct_w(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[])
864: {
865: const PetscReal Pe = PetscRealPart(constants[PECLET]);
866: const PetscReal k = PetscRealPart(constants[K]);
867: PetscInt d;
869: // \frac{k}{Pe} \nabla T
870: for (d = 0; d < dim; ++d) f1[d] = k / Pe * u_x[uOff_x[TEMP] + d];
871: }
873: static void g1_qu(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[])
874: {
875: PetscInt d;
876: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
877: }
879: static void g0_vu(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[])
880: {
881: PetscInt c, d;
882: const PetscInt Nc = dim;
884: for (d = 0; d < dim; ++d) g0[d * dim + d] = u_tShift;
886: for (c = 0; c < Nc; ++c) {
887: for (d = 0; d < dim; ++d) g0[c * Nc + d] += u_x[c * Nc + d];
888: }
889: }
891: static void g1_vu(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[])
892: {
893: PetscInt NcI = dim;
894: PetscInt NcJ = dim;
895: PetscInt c, d, e;
897: for (c = 0; c < NcI; ++c) {
898: for (d = 0; d < NcJ; ++d) {
899: for (e = 0; e < dim; ++e) {
900: if (c == d) g1[(c * NcJ + d) * dim + e] += u[e];
901: }
902: }
903: }
904: }
906: static void g0_conduct_qu(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[])
907: {
908: const PetscReal p_th = PetscRealPart(constants[P_TH]);
909: PetscInt d;
911: // - \phi_i \frac{p^{th}}{T^2} \frac{\partial T}{\partial x_c} \psi_{j, u_c}
912: for (d = 0; d < dim; ++d) g0[d] = -p_th / PetscSqr(u[uOff[TEMP]]) * u_x[uOff_x[TEMP] + d];
913: }
915: static void g1_conduct_qu(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[])
916: {
917: const PetscReal p_th = PetscRealPart(constants[P_TH]);
918: PetscInt d;
920: // \phi_i \frac{p^{th}}{T} \frac{\partial \psi_{u_c,j}}{\partial x_c}
921: for (d = 0; d < dim; ++d) g1[d * dim + d] = p_th / u[uOff[TEMP]];
922: }
924: static void g0_conduct_qT(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[])
925: {
926: const PetscReal S = PetscRealPart(constants[STROUHAL]);
927: const PetscReal p_th = PetscRealPart(constants[P_TH]);
928: PetscInt d;
930: // - \phi_i \frac{S p^{th}}{T^2} \psi_j
931: g0[0] -= S * p_th / PetscSqr(u[uOff[TEMP]]) * u_tShift;
932: // \phi_i 2 \frac{S p^{th}}{T^3} T_t \psi_j
933: g0[0] += 2.0 * S * p_th / PetscPowScalarInt(u[uOff[TEMP]], 3) * u_t[uOff[TEMP]];
934: // \phi_i \frac{p^{th}}{T^2} \left( - \nabla \cdot \vb{u} \psi_j + \frac{2}{T} \vb{u} \cdot \nabla T \psi_j \right)
935: for (d = 0; d < dim; ++d) g0[0] += p_th / PetscSqr(u[uOff[TEMP]]) * (-u_x[uOff_x[VEL] + d * dim + d] + 2.0 / u[uOff[TEMP]] * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d]);
936: }
938: static void g1_conduct_qT(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[])
939: {
940: const PetscReal p_th = PetscRealPart(constants[P_TH]);
941: PetscInt d;
943: // - \phi_i \frac{p^{th}}{T^2} \vb{u} \cdot \nabla \psi_j
944: for (d = 0; d < dim; ++d) g1[d] = -p_th / PetscSqr(u[uOff[TEMP]]) * u[uOff[VEL] + d];
945: }
947: static void g2_vp(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 g2[])
948: {
949: PetscInt d;
950: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
951: }
953: static void g3_vu(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[])
954: {
955: const PetscReal nu = PetscRealPart(constants[NU]);
956: const PetscInt Nc = dim;
957: PetscInt c, d;
959: for (c = 0; c < Nc; ++c) {
960: for (d = 0; d < dim; ++d) {
961: g3[((c * Nc + c) * dim + d) * dim + d] += nu;
962: g3[((c * Nc + d) * dim + d) * dim + c] += nu;
963: }
964: }
965: }
967: static void g0_conduct_vT(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[])
968: {
969: const PetscReal S = PetscRealPart(constants[STROUHAL]);
970: const PetscReal F = PetscRealPart(constants[FROUDE]);
971: const PetscReal p_th = PetscRealPart(constants[P_TH]);
972: const PetscInt gdir = (PetscInt)PetscRealPart(constants[G_DIR]);
973: const PetscInt Nc = dim;
974: PetscInt c, d;
976: // - \vb{\phi}_i \cdot \vb{u}_t \frac{p^{th} S}{T^2} \psi_j
977: for (d = 0; d < dim; ++d) g0[d] -= p_th * S / PetscSqr(u[uOff[TEMP]]) * u_t[uOff[VEL] + d];
979: // - \vb{\phi}_i \cdot \vb{u} \cdot \nabla \vb{u} \frac{p^{th}}{T^2} \psi_j
980: for (c = 0; c < Nc; ++c) {
981: for (d = 0; d < dim; ++d) g0[c] -= p_th / PetscSqr(u[uOff[TEMP]]) * u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
982: }
984: // - \vb{\phi}_i \cdot \vu{z} \frac{p^{th}}{T^2 F^2} \psi_j
985: g0[gdir] -= p_th / PetscSqr(u[uOff[TEMP]] * F);
986: }
988: static void g0_conduct_vu(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[])
989: {
990: const PetscReal S = PetscRealPart(constants[STROUHAL]);
991: const PetscReal p_th = PetscRealPart(constants[P_TH]);
992: const PetscInt Nc = dim;
993: PetscInt c, d;
995: // \vb{\phi}_i \cdot S \rho \psi_j
996: for (d = 0; d < dim; ++d) g0[d * dim + d] = S * p_th / u[uOff[TEMP]] * u_tShift;
998: // \phi^c_i \cdot \rho \frac{\partial u^c}{\partial x^d} \psi^d_j
999: for (c = 0; c < Nc; ++c) {
1000: for (d = 0; d < dim; ++d) g0[c * Nc + d] += p_th / u[uOff[TEMP]] * u_x[uOff_x[VEL] + c * Nc + d];
1001: }
1002: }
1004: static void g1_conduct_vu(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[])
1005: {
1006: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1007: const PetscInt NcI = dim;
1008: const PetscInt NcJ = dim;
1010: // \phi^c_i \rho u^e \frac{\partial \psi^d_j}{\partial x^e}
1011: for (PetscInt c = 0; c < NcI; ++c) {
1012: for (PetscInt d = 0; d < NcJ; ++d) {
1013: for (PetscInt e = 0; e < dim; ++e) {
1014: if (c == d) g1[(c * NcJ + d) * dim + e] += p_th / u[uOff[TEMP]] * u[uOff[VEL] + e];
1015: }
1016: }
1017: }
1018: }
1020: static void g3_conduct_vu(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[])
1021: {
1022: const PetscReal Re = PetscRealPart(constants[REYNOLDS]);
1023: const PetscReal mu = PetscRealPart(constants[MU]);
1024: const PetscInt Nc = dim;
1026: for (PetscInt c = 0; c < Nc; ++c) {
1027: for (PetscInt d = 0; d < dim; ++d) {
1028: // \frac{\partial \phi^c_i}{\partial x^d} \mu/Re \frac{\partial \psi^c_i}{\partial x^d}
1029: g3[((c * Nc + c) * dim + d) * dim + d] += mu / Re; // gradU
1030: // \frac{\partial \phi^c_i}{\partial x^d} \mu/Re \frac{\partial \psi^d_i}{\partial x^c}
1031: g3[((c * Nc + d) * dim + d) * dim + c] += mu / Re; // gradU transpose
1032: // \frac{\partial \phi^c_i}{\partial x^d} -2/3 \mu/Re \frac{\partial \psi^d_i}{\partial x^c}
1033: g3[((c * Nc + d) * dim + c) * dim + d] -= 2.0 / 3.0 * mu / Re;
1034: }
1035: }
1036: }
1038: static void g2_conduct_vp(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 g2[])
1039: {
1040: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
1041: }
1043: static void g0_wT(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[])
1044: {
1045: g0[0] = u_tShift;
1046: }
1048: static void g0_wu(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[])
1049: {
1050: PetscInt d;
1051: for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2] + d];
1052: }
1054: static void g1_wT(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[])
1055: {
1056: PetscInt d;
1057: for (d = 0; d < dim; ++d) g1[d] = u[uOff[0] + d];
1058: }
1060: static void g3_wT(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[])
1061: {
1062: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
1063: PetscInt d;
1065: for (d = 0; d < dim; ++d) g3[d * dim + d] = alpha;
1066: }
1068: static void g0_conduct_wu(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[])
1069: {
1070: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1071: const PetscReal c_p = PetscRealPart(constants[C_P]);
1072: PetscInt d;
1074: // \phi_i \frac{C_p p^{th}}{T} \nabla T \cdot \psi_j
1075: for (d = 0; d < dim; ++d) g0[d] = c_p * p_th / u[uOff[TEMP]] * u_x[uOff_x[TEMP] + d];
1076: }
1078: static void g0_conduct_wT(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[])
1079: {
1080: const PetscReal S = PetscRealPart(constants[STROUHAL]);
1081: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1082: const PetscReal c_p = PetscRealPart(constants[C_P]);
1084: // \psi_i C_p S p^{th}\T \psi_{j}
1085: g0[0] += c_p * S * p_th / u[uOff[TEMP]] * u_tShift;
1086: // - \phi_i C_p S p^{th}/T^2 T_t \psi_j
1087: g0[0] -= c_p * S * p_th / PetscSqr(u[uOff[TEMP]]) * u_t[uOff[TEMP]];
1088: // - \phi_i C_p p^{th}/T^2 \vb{u} \cdot \nabla T \psi_j
1089: for (PetscInt d = 0; d < dim; ++d) g0[0] -= c_p * p_th / PetscSqr(u[uOff[TEMP]]) * u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];
1090: }
1092: static void g1_conduct_wT(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[])
1093: {
1094: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1095: const PetscReal c_p = PetscRealPart(constants[C_P]);
1096: // \phi_i C_p p^{th}/T \vb{u} \cdot \nabla \psi_j
1097: for (PetscInt d = 0; d < dim; ++d) g1[d] += c_p * p_th / u[uOff[TEMP]] * u[uOff[VEL] + d];
1098: }
1100: static void g3_conduct_wT(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[])
1101: {
1102: const PetscReal Pe = PetscRealPart(constants[PECLET]);
1103: const PetscReal k = PetscRealPart(constants[K]);
1104: // \nabla \phi_i \frac{k}{Pe} \nabla \phi_j
1105: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = k / Pe;
1106: }
1108: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
1109: {
1110: PetscInt mod, sol;
1112: PetscFunctionBeginUser;
1113: options->modType = MOD_INCOMPRESSIBLE;
1114: options->solType = SOL_QUADRATIC;
1115: options->hasNullSpace = PETSC_TRUE;
1116: options->dmCell = NULL;
1118: PetscOptionsBegin(comm, "", "Low Mach flow Problem Options", "DMPLEX");
1119: mod = options->modType;
1120: PetscCall(PetscOptionsEList("-mod_type", "The model type", "ex76.c", modTypes, NUM_MOD_TYPES, modTypes[options->modType], &mod, NULL));
1121: options->modType = (ModType)mod;
1122: sol = options->solType;
1123: PetscCall(PetscOptionsEList("-sol_type", "The solution type", "ex76.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL));
1124: options->solType = (SolType)sol;
1125: PetscOptionsEnd();
1126: PetscFunctionReturn(PETSC_SUCCESS);
1127: }
1129: static PetscErrorCode SetupParameters(DM dm, AppCtx *user)
1130: {
1131: PetscBag bag;
1132: Parameter *p;
1133: PetscReal dir;
1134: PetscInt dim;
1136: PetscFunctionBeginUser;
1137: PetscCall(DMGetDimension(dm, &dim));
1138: dir = (PetscReal)(dim - 1);
1139: /* setup PETSc parameter bag */
1140: PetscCall(PetscBagGetData(user->bag, &p));
1141: PetscCall(PetscBagSetName(user->bag, "par", "Low Mach flow parameters"));
1142: bag = user->bag;
1143: PetscCall(PetscBagRegisterReal(bag, &p->Strouhal, 1.0, "S", "Strouhal number"));
1144: PetscCall(PetscBagRegisterReal(bag, &p->Froude, 1.0, "Fr", "Froude number"));
1145: PetscCall(PetscBagRegisterReal(bag, &p->Reynolds, 1.0, "Re", "Reynolds number"));
1146: PetscCall(PetscBagRegisterReal(bag, &p->Peclet, 1.0, "Pe", "Peclet number"));
1147: PetscCall(PetscBagRegisterReal(bag, &p->p_th, 1.0, "p_th", "Thermodynamic pressure"));
1148: PetscCall(PetscBagRegisterReal(bag, &p->mu, 1.0, "mu", "Dynamic viscosity"));
1149: PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
1150: PetscCall(PetscBagRegisterReal(bag, &p->c_p, 1.0, "c_p", "Specific heat at constant pressure"));
1151: PetscCall(PetscBagRegisterReal(bag, &p->k, 1.0, "k", "Thermal conductivity"));
1152: PetscCall(PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity"));
1153: PetscCall(PetscBagRegisterReal(bag, &p->T_in, 1.0, "T_in", "Inlet temperature"));
1154: PetscCall(PetscBagRegisterReal(bag, &p->g_dir, dir, "g_dir", "Gravity direction"));
1155: PetscCall(PetscBagRegisterReal(bag, &p->epsilon, 1.0, "epsilon", "Perturbation strength"));
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
1160: {
1161: PetscFunctionBeginUser;
1162: PetscCall(DMCreate(comm, dm));
1163: PetscCall(DMSetType(*dm, DMPLEX));
1164: PetscCall(DMSetFromOptions(*dm));
1165: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: static PetscErrorCode UniformBoundaryConditions(DM dm, DMLabel label, PetscSimplePointFn *exactFuncs[], PetscSimplePointFn *exactFuncs_t[], AppCtx *user)
1170: {
1171: PetscDS ds;
1172: PetscInt id;
1173: void *ctx;
1175: PetscFunctionBeginUser;
1176: PetscCall(DMGetDS(dm, &ds));
1177: PetscCall(PetscBagGetData(user->bag, &ctx));
1178: id = 3;
1179: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1180: id = 1;
1181: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1182: id = 2;
1183: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1184: id = 4;
1185: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1186: id = 3;
1187: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1188: id = 1;
1189: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1190: id = 2;
1191: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1192: id = 4;
1193: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
1198: {
1199: PetscSimplePointFn *exactFuncs[3];
1200: PetscSimplePointFn *exactFuncs_t[3];
1201: PetscDS ds;
1202: PetscWeakForm wf;
1203: DMLabel label;
1204: Parameter *ctx;
1205: PetscInt id, bd;
1207: PetscFunctionBeginUser;
1208: PetscCall(DMGetLabel(dm, "marker", &label));
1209: PetscCall(DMGetDS(dm, &ds));
1210: PetscCall(PetscDSGetWeakForm(ds, &wf));
1212: switch (user->modType) {
1213: case MOD_INCOMPRESSIBLE:
1214: PetscCall(PetscDSSetResidual(ds, VEL, f0_v, f1_v));
1215: PetscCall(PetscDSSetResidual(ds, PRES, f0_q, NULL));
1216: PetscCall(PetscDSSetResidual(ds, TEMP, f0_w, f1_w));
1218: PetscCall(PetscDSSetJacobian(ds, VEL, VEL, g0_vu, g1_vu, NULL, g3_vu));
1219: PetscCall(PetscDSSetJacobian(ds, VEL, PRES, NULL, NULL, g2_vp, NULL));
1220: PetscCall(PetscDSSetJacobian(ds, PRES, VEL, NULL, g1_qu, NULL, NULL));
1221: PetscCall(PetscDSSetJacobian(ds, TEMP, VEL, g0_wu, NULL, NULL, NULL));
1222: PetscCall(PetscDSSetJacobian(ds, TEMP, TEMP, g0_wT, g1_wT, NULL, g3_wT));
1224: switch (user->solType) {
1225: case SOL_QUADRATIC:
1226: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_quadratic_v, 0, NULL));
1227: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_quadratic_w, 0, NULL));
1229: exactFuncs[VEL] = quadratic_u;
1230: exactFuncs[PRES] = quadratic_p;
1231: exactFuncs[TEMP] = quadratic_T;
1232: exactFuncs_t[VEL] = quadratic_u_t;
1233: exactFuncs_t[PRES] = NULL;
1234: exactFuncs_t[TEMP] = quadratic_T_t;
1236: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1237: break;
1238: case SOL_CUBIC:
1239: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_cubic_v, 0, NULL));
1240: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_cubic_w, 0, NULL));
1242: exactFuncs[VEL] = cubic_u;
1243: exactFuncs[PRES] = cubic_p;
1244: exactFuncs[TEMP] = cubic_T;
1245: exactFuncs_t[VEL] = cubic_u_t;
1246: exactFuncs_t[PRES] = NULL;
1247: exactFuncs_t[TEMP] = cubic_T_t;
1249: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1250: break;
1251: case SOL_CUBIC_TRIG:
1252: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_cubic_trig_v, 0, NULL));
1253: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_cubic_trig_w, 0, NULL));
1255: exactFuncs[VEL] = cubic_trig_u;
1256: exactFuncs[PRES] = cubic_trig_p;
1257: exactFuncs[TEMP] = cubic_trig_T;
1258: exactFuncs_t[VEL] = cubic_trig_u_t;
1259: exactFuncs_t[PRES] = NULL;
1260: exactFuncs_t[TEMP] = cubic_trig_T_t;
1262: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1263: break;
1264: case SOL_TAYLOR_GREEN:
1265: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_taylor_green_w, 0, NULL));
1267: exactFuncs[VEL] = taylor_green_u;
1268: exactFuncs[PRES] = taylor_green_p;
1269: exactFuncs[TEMP] = taylor_green_T;
1270: exactFuncs_t[VEL] = taylor_green_u_t;
1271: exactFuncs_t[PRES] = taylor_green_p_t;
1272: exactFuncs_t[TEMP] = taylor_green_T_t;
1274: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1275: break;
1276: default:
1277: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
1278: }
1279: break;
1280: case MOD_CONDUCTING:
1281: PetscCall(PetscDSSetResidual(ds, VEL, f0_conduct_v, f1_conduct_v));
1282: PetscCall(PetscDSSetResidual(ds, PRES, f0_conduct_q, NULL));
1283: PetscCall(PetscDSSetResidual(ds, TEMP, f0_conduct_w, f1_conduct_w));
1285: PetscCall(PetscDSSetJacobian(ds, VEL, VEL, g0_conduct_vu, g1_conduct_vu, NULL, g3_conduct_vu));
1286: PetscCall(PetscDSSetJacobian(ds, VEL, PRES, NULL, NULL, g2_conduct_vp, NULL));
1287: PetscCall(PetscDSSetJacobian(ds, VEL, TEMP, g0_conduct_vT, NULL, NULL, NULL));
1288: PetscCall(PetscDSSetJacobian(ds, PRES, VEL, g0_conduct_qu, g1_conduct_qu, NULL, NULL));
1289: PetscCall(PetscDSSetJacobian(ds, PRES, TEMP, g0_conduct_qT, g1_conduct_qT, NULL, NULL));
1290: PetscCall(PetscDSSetJacobian(ds, TEMP, VEL, g0_conduct_wu, NULL, NULL, NULL));
1291: PetscCall(PetscDSSetJacobian(ds, TEMP, TEMP, g0_conduct_wT, g1_conduct_wT, NULL, g3_conduct_wT));
1293: switch (user->solType) {
1294: case SOL_QUADRATIC:
1295: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_quadratic_v, 0, NULL));
1296: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_quadratic_q, 0, NULL));
1297: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_quadratic_w, 0, NULL));
1299: exactFuncs[VEL] = quadratic_u;
1300: exactFuncs[PRES] = quadratic_p;
1301: exactFuncs[TEMP] = quadratic_T;
1302: exactFuncs_t[VEL] = quadratic_u_t;
1303: exactFuncs_t[PRES] = NULL;
1304: exactFuncs_t[TEMP] = quadratic_T_t;
1306: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1307: break;
1308: case SOL_PIPE:
1309: user->hasNullSpace = PETSC_FALSE;
1310: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_pipe_v, 0, NULL));
1311: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_pipe_q, 0, NULL));
1312: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_pipe_w, 0, NULL));
1314: exactFuncs[VEL] = pipe_u;
1315: exactFuncs[PRES] = pipe_p;
1316: exactFuncs[TEMP] = pipe_T;
1317: exactFuncs_t[VEL] = pipe_u_t;
1318: exactFuncs_t[PRES] = pipe_p_t;
1319: exactFuncs_t[TEMP] = pipe_T_t;
1321: PetscCall(PetscBagGetData(user->bag, &ctx));
1322: id = 2;
1323: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1324: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1325: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_v, 0, NULL));
1326: id = 4;
1327: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1328: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1329: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_v, 0, NULL));
1330: id = 4;
1331: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1332: id = 3;
1333: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1334: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1335: id = 1;
1336: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1337: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1338: break;
1339: case SOL_PIPE_WIGGLY:
1340: user->hasNullSpace = PETSC_FALSE;
1341: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_pipe_wiggly_v, 0, NULL));
1342: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_pipe_wiggly_q, 0, NULL));
1343: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_pipe_wiggly_w, 0, NULL));
1345: exactFuncs[VEL] = pipe_wiggly_u;
1346: exactFuncs[PRES] = pipe_wiggly_p;
1347: exactFuncs[TEMP] = pipe_wiggly_T;
1348: exactFuncs_t[VEL] = pipe_wiggly_u_t;
1349: exactFuncs_t[PRES] = pipe_wiggly_p_t;
1350: exactFuncs_t[TEMP] = pipe_wiggly_T_t;
1352: PetscCall(PetscBagGetData(user->bag, &ctx));
1353: id = 2;
1354: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1355: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1356: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_wiggly_v, 0, NULL));
1357: id = 4;
1358: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1359: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1360: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_wiggly_v, 0, NULL));
1361: id = 4;
1362: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1363: id = 3;
1364: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1365: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1366: id = 1;
1367: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (PetscVoidFn *)exactFuncs[VEL], (PetscVoidFn *)exactFuncs_t[VEL], ctx, NULL));
1368: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temperature", label, 1, &id, TEMP, 0, NULL, (PetscVoidFn *)exactFuncs[TEMP], (PetscVoidFn *)exactFuncs_t[TEMP], ctx, NULL));
1369: break;
1370: default:
1371: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
1372: }
1373: break;
1374: default:
1375: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported model type: %s (%d)", modTypes[PetscMin(user->modType, NUM_MOD_TYPES)], user->modType);
1376: }
1377: /* Setup constants */
1378: {
1379: Parameter *param;
1380: PetscScalar constants[13];
1382: PetscCall(PetscBagGetData(user->bag, ¶m));
1384: constants[STROUHAL] = param->Strouhal;
1385: constants[FROUDE] = param->Froude;
1386: constants[REYNOLDS] = param->Reynolds;
1387: constants[PECLET] = param->Peclet;
1388: constants[P_TH] = param->p_th;
1389: constants[MU] = param->mu;
1390: constants[NU] = param->nu;
1391: constants[C_P] = param->c_p;
1392: constants[K] = param->k;
1393: constants[ALPHA] = param->alpha;
1394: constants[T_IN] = param->T_in;
1395: constants[G_DIR] = param->g_dir;
1396: constants[EPSILON] = param->epsilon;
1397: PetscCall(PetscDSSetConstants(ds, 13, constants));
1398: }
1400: PetscCall(PetscBagGetData(user->bag, &ctx));
1401: PetscCall(PetscDSSetExactSolution(ds, VEL, exactFuncs[VEL], ctx));
1402: PetscCall(PetscDSSetExactSolution(ds, PRES, exactFuncs[PRES], ctx));
1403: PetscCall(PetscDSSetExactSolution(ds, TEMP, exactFuncs[TEMP], ctx));
1404: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, VEL, exactFuncs_t[VEL], ctx));
1405: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, PRES, exactFuncs_t[PRES], ctx));
1406: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, TEMP, exactFuncs_t[TEMP], ctx));
1407: PetscFunctionReturn(PETSC_SUCCESS);
1408: }
1410: static PetscErrorCode CreateCellDM(DM dm, AppCtx *user)
1411: {
1412: PetscFE fe, fediv;
1413: DMPolytopeType ct;
1414: PetscInt dim, cStart;
1415: PetscBool simplex;
1417: PetscFunctionBeginUser;
1418: PetscCall(DMGetDimension(dm, &dim));
1419: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
1420: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1421: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1423: PetscCall(DMGetField(dm, VEL, NULL, (PetscObject *)&fe));
1424: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "div_", PETSC_DEFAULT, &fediv));
1425: PetscCall(PetscFECopyQuadrature(fe, fediv));
1426: PetscCall(PetscObjectSetName((PetscObject)fediv, "divergence"));
1428: PetscCall(DMDestroy(&user->dmCell));
1429: PetscCall(DMClone(dm, &user->dmCell));
1430: PetscCall(DMSetField(user->dmCell, 0, NULL, (PetscObject)fediv));
1431: PetscCall(DMCreateDS(user->dmCell));
1432: PetscCall(PetscFEDestroy(&fediv));
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: static PetscErrorCode GetCellDM(DM dm, AppCtx *user, DM *dmCell)
1437: {
1438: PetscInt cStart, cEnd, cellStart = -1, cellEnd = -1;
1440: PetscFunctionBeginUser;
1441: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1442: if (user->dmCell) PetscCall(DMPlexGetSimplexOrBoxCells(user->dmCell, 0, &cellStart, &cellEnd));
1443: if (cStart != cellStart || cEnd != cellEnd) PetscCall(CreateCellDM(dm, user));
1444: *dmCell = user->dmCell;
1445: PetscFunctionReturn(PETSC_SUCCESS);
1446: }
1448: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
1449: {
1450: DM cdm = dm;
1451: PetscFE fe[3];
1452: Parameter *param;
1453: DMPolytopeType ct;
1454: PetscInt dim, cStart;
1455: PetscBool simplex;
1457: PetscFunctionBeginUser;
1458: PetscCall(DMGetDimension(dm, &dim));
1459: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
1460: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1461: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1462: /* Create finite element */
1463: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]));
1464: PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
1466: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
1467: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
1468: PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
1470: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]));
1471: PetscCall(PetscFECopyQuadrature(fe[0], fe[2]));
1472: PetscCall(PetscObjectSetName((PetscObject)fe[2], "temperature"));
1474: /* Set discretization and boundary conditions for each mesh */
1475: PetscCall(DMSetField(dm, VEL, NULL, (PetscObject)fe[VEL]));
1476: PetscCall(DMSetField(dm, PRES, NULL, (PetscObject)fe[PRES]));
1477: PetscCall(DMSetField(dm, TEMP, NULL, (PetscObject)fe[TEMP]));
1478: PetscCall(DMCreateDS(dm));
1479: PetscCall(SetupProblem(dm, user));
1480: PetscCall(PetscBagGetData(user->bag, ¶m));
1481: while (cdm) {
1482: PetscCall(DMCopyDisc(dm, cdm));
1483: PetscCall(DMGetCoarseDM(cdm, &cdm));
1484: }
1485: PetscCall(PetscFEDestroy(&fe[VEL]));
1486: PetscCall(PetscFEDestroy(&fe[PRES]));
1487: PetscCall(PetscFEDestroy(&fe[TEMP]));
1489: if (user->hasNullSpace) {
1490: PetscObject pressure;
1491: MatNullSpace nullspacePres;
1493: PetscCall(DMGetField(dm, PRES, NULL, &pressure));
1494: PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres));
1495: PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres));
1496: PetscCall(MatNullSpaceDestroy(&nullspacePres));
1497: }
1498: PetscFunctionReturn(PETSC_SUCCESS);
1499: }
1501: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
1502: {
1503: Vec vec;
1504: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero};
1506: PetscFunctionBeginUser;
1507: PetscCheck(ofield == PRES, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index %" PetscInt_FMT ", not %" PetscInt_FMT, PRES, ofield);
1508: funcs[nfield] = constant;
1509: PetscCall(DMCreateGlobalVector(dm, &vec));
1510: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
1511: PetscCall(VecNormalize(vec, NULL));
1512: PetscCall(PetscObjectSetName((PetscObject)vec, "Pressure Null Space"));
1513: PetscCall(VecViewFromOptions(vec, NULL, "-pressure_nullspace_view"));
1514: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
1515: PetscCall(VecDestroy(&vec));
1516: PetscFunctionReturn(PETSC_SUCCESS);
1517: }
1519: static PetscErrorCode RemoveDiscretePressureNullspace_Private(TS ts, Vec u)
1520: {
1521: DM dm;
1522: AppCtx *user;
1523: MatNullSpace nullsp;
1525: PetscFunctionBeginUser;
1526: PetscCall(TSGetDM(ts, &dm));
1527: PetscCall(DMGetApplicationContext(dm, &user));
1528: if (!user->hasNullSpace) PetscFunctionReturn(PETSC_SUCCESS);
1529: PetscCall(CreatePressureNullSpace(dm, 1, 1, &nullsp));
1530: PetscCall(MatNullSpaceRemove(nullsp, u));
1531: PetscCall(MatNullSpaceDestroy(&nullsp));
1532: PetscFunctionReturn(PETSC_SUCCESS);
1533: }
1535: /* Make the discrete pressure discretely divergence free */
1536: static PetscErrorCode RemoveDiscretePressureNullspace(TS ts)
1537: {
1538: Vec u;
1540: PetscFunctionBeginUser;
1541: PetscCall(TSGetSolution(ts, &u));
1542: PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
1543: PetscFunctionReturn(PETSC_SUCCESS);
1544: }
1546: static void divergence(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 divu[])
1547: {
1548: divu[0] = 0.;
1549: for (PetscInt d = 0; d < dim; ++d) divu[0] += u_x[d * dim + d];
1550: }
1552: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
1553: {
1554: AppCtx *user;
1555: DM dm;
1556: PetscReal t;
1558: PetscFunctionBeginUser;
1559: PetscCall(TSGetDM(ts, &dm));
1560: PetscCall(TSGetTime(ts, &t));
1561: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1562: PetscCall(DMGetApplicationContext(dm, &user));
1563: PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
1564: PetscFunctionReturn(PETSC_SUCCESS);
1565: }
1567: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, PetscCtx ctx)
1568: {
1569: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
1570: void *ctxs[3];
1571: PetscPointFn *diagnostics[1] = {divergence};
1572: DM dm, dmCell = NULL;
1573: PetscDS ds;
1574: Vec v, divu;
1575: PetscReal ferrors[3], massFlux;
1577: PetscFunctionBeginUser;
1578: PetscCall(TSGetDM(ts, &dm));
1579: PetscCall(DMGetDS(dm, &ds));
1581: for (PetscInt f = 0; f < 3; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]));
1582: PetscCall(DMComputeL2FieldDiff(dm, crtime, exactFuncs, ctxs, u, ferrors));
1583: PetscCall(GetCellDM(dm, (AppCtx *)ctx, &dmCell));
1584: PetscCall(DMGetGlobalVector(dmCell, &divu));
1585: PetscCall(DMProjectField(dmCell, crtime, u, diagnostics, INSERT_VALUES, divu));
1586: PetscCall(VecViewFromOptions(divu, NULL, "-divu_vec_view"));
1587: PetscCall(VecNorm(divu, NORM_2, &massFlux));
1588: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g, %2.3g] ||div u||: %2.3g\n", (int)step, (double)crtime, (double)ferrors[0], (double)ferrors[1], (double)ferrors[2], (double)massFlux));
1590: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
1592: PetscCall(DMGetGlobalVector(dm, &v));
1593: PetscCall(DMProjectFunction(dm, crtime, exactFuncs, ctxs, INSERT_ALL_VALUES, v));
1594: PetscCall(PetscObjectSetName((PetscObject)v, "Exact Solution"));
1595: PetscCall(VecViewFromOptions(v, NULL, "-exact_vec_view"));
1596: PetscCall(DMRestoreGlobalVector(dm, &v));
1598: PetscCall(VecViewFromOptions(divu, NULL, "-div_vec_view"));
1599: PetscCall(DMRestoreGlobalVector(dmCell, &divu));
1600: PetscFunctionReturn(PETSC_SUCCESS);
1601: }
1603: int main(int argc, char **argv)
1604: {
1605: DM dm; /* problem definition */
1606: TS ts; /* timestepper */
1607: Vec u; /* solution */
1608: AppCtx user; /* user-defined work context */
1609: PetscReal t;
1611: PetscFunctionBeginUser;
1612: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1613: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1614: PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
1615: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1616: PetscCall(SetupParameters(dm, &user));
1617: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1618: PetscCall(TSSetDM(ts, dm));
1619: PetscCall(DMSetApplicationContext(dm, &user));
1620: /* Setup problem */
1621: PetscCall(SetupDiscretization(dm, &user));
1622: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
1624: PetscCall(DMCreateGlobalVector(dm, &u));
1625: PetscCall(PetscObjectSetName((PetscObject)u, "Numerical Solution"));
1626: if (user.hasNullSpace) PetscCall(DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace));
1628: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
1629: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
1630: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
1631: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1632: PetscCall(TSSetPreStep(ts, RemoveDiscretePressureNullspace));
1633: PetscCall(TSSetFromOptions(ts));
1635: PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions)); /* Must come after SetFromOptions() */
1636: PetscCall(SetInitialConditions(ts, u));
1637: PetscCall(TSGetTime(ts, &t));
1638: PetscCall(DMSetOutputSequenceNumber(dm, 0, t));
1639: PetscCall(DMTSCheckFromOptions(ts, u));
1640: PetscCall(TSMonitorSet(ts, MonitorError, &user, NULL));
1642: PetscCall(TSSolve(ts, u));
1643: PetscCall(DMTSCheckFromOptions(ts, u));
1645: PetscCall(VecDestroy(&u));
1646: PetscCall(DMDestroy(&user.dmCell));
1647: PetscCall(DMDestroy(&dm));
1648: PetscCall(TSDestroy(&ts));
1649: PetscCall(PetscBagDestroy(&user.bag));
1650: PetscCall(PetscFinalize());
1651: return 0;
1652: }
1654: /*TEST
1656: testset:
1657: requires: triangle !single
1658: args: -dm_plex_separate_marker \
1659: -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \
1660: -snes_error_if_not_converged -snes_convergence_test correct_pressure \
1661: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
1662: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 \
1663: -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
1664: -fieldsplit_0_pc_type lu \
1665: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1667: test:
1668: suffix: 2d_tri_p2_p1_p1
1669: args: -sol_type quadratic \
1670: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1671: -dmts_check .001 -ts_max_steps 4 -ts_time_step 0.1
1673: test:
1674: # Using -dm_refine 5 -convest_num_refine 2 gives L_2 convergence rate: [0.89, 0.011, 1.0]
1675: suffix: 2d_tri_p2_p1_p1_tconv
1676: args: -sol_type cubic_trig \
1677: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1678: -ts_max_steps 4 -ts_time_step 0.1 -ts_convergence_estimate -convest_num_refine 1
1680: test:
1681: # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.5, 1.9]
1682: suffix: 2d_tri_p2_p1_p1_sconv
1683: args: -sol_type cubic \
1684: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1685: -ts_max_steps 1 -ts_time_step 1e-4 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
1687: test:
1688: suffix: 2d_tri_p3_p2_p2
1689: args: -sol_type cubic \
1690: -vel_petscspace_degree 3 -pres_petscspace_degree 2 -temp_petscspace_degree 2 \
1691: -dmts_check .001 -ts_max_steps 4 -ts_time_step 0.1
1693: test:
1694: # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.1, 3.1]
1695: suffix: 2d_tri_p2_p1_p1_tg_sconv
1696: args: -sol_type taylor_green \
1697: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1698: -ts_max_steps 1 -ts_time_step 1e-8 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
1700: test:
1701: # Using -dm_refine 3 -convest_num_refine 2 gives L_2 convergence rate: [1.2, 1.5, 1.2]
1702: suffix: 2d_tri_p2_p1_p1_tg_tconv
1703: args: -sol_type taylor_green \
1704: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1705: -ts_max_steps 4 -ts_time_step 0.1 -ts_convergence_estimate -convest_num_refine 1
1707: testset:
1708: requires: triangle !single
1709: args: -dm_plex_separate_marker -mod_type conducting \
1710: -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \
1711: -snes_error_if_not_converged -snes_max_linear_solve_fail 5 \
1712: -ksp_type fgmres -ksp_max_it 2 -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 \
1713: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 \
1714: -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
1715: -fieldsplit_0_pc_type lu \
1716: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1718: test:
1719: # At this resolution, the rhs is inconsistent on some Newton steps
1720: suffix: 2d_tri_p2_p1_p1_conduct
1721: args: -sol_type quadratic \
1722: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1723: -dmts_check .001 -ts_max_steps 4 -ts_time_step 0.1 \
1724: -pc_fieldsplit_schur_precondition full \
1725: -fieldsplit_pressure_ksp_max_it 2 -fieldsplit_pressure_pc_type svd
1727: test:
1728: suffix: 2d_tri_p2_p1_p2_conduct_pipe
1729: args: -sol_type pipe \
1730: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 2 \
1731: -dmts_check .001 -ts_max_steps 4 -ts_time_step 0.1
1733: test:
1734: suffix: 2d_tri_p2_p1_p2_conduct_pipe_wiggly_sconv
1735: args: -sol_type pipe_wiggly -Fr 1e10 \
1736: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 2 \
1737: -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
1738: -ts_max_steps 1 -ts_time_step 1e10 \
1739: -ksp_atol 1e-12 -ksp_max_it 300 \
1740: -fieldsplit_pressure_ksp_atol 1e-14
1742: TEST*/