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, void *ctx)
107: {
108: PetscInt d;
109: for (d = 0; d < Nc; ++d) u[d] = 0.0;
110: return PETSC_SUCCESS;
111: }
113: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
114: {
115: PetscInt d;
116: for (d = 0; d < Nc; ++d) u[d] = 1.0;
117: return PETSC_SUCCESS;
118: }
120: /*
121: CASE: quadratic
122: In 2D we use exact solution:
124: u = t + x^2 + y^2
125: v = t + 2x^2 - 2xy
126: p = x + y - 1
127: T = t + x + y + 1
128: 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>
129: Q = 1 + 2t + 3x^2 - 2xy + y^2
131: so that
133: \nabla \cdot u = 2x - 2x = 0
135: f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
136: = <1, 1> + <t + x^2 + y^2, t + 2x^2 - 2xy> . <<2x, 4x - 2y>, <2y, -2x>> - \nu <4, 4> + <1, 1>
137: = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t (2x - 2y) + 2x^2y + 4xy^2 - 2y^3> + <-4 \nu + 2, -4\nu + 2>
138: = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2 - 4\nu + 2, t (2x - 2y) + 4xy^2 + 2x^2y - 2y^3 - 4\nu + 2>
140: Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
141: = 1 + <t + x^2 + y^2, t + 2x^2 - 2xy> . <1, 1> - \alpha 0
142: = 1 + 2t + 3x^2 - 2xy + y^2
143: */
145: static PetscErrorCode quadratic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
146: {
147: u[0] = time + X[0] * X[0] + X[1] * X[1];
148: u[1] = time + 2.0 * X[0] * X[0] - 2.0 * X[0] * X[1];
149: return PETSC_SUCCESS;
150: }
151: static PetscErrorCode quadratic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
152: {
153: u[0] = 1.0;
154: u[1] = 1.0;
155: return PETSC_SUCCESS;
156: }
158: static PetscErrorCode quadratic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
159: {
160: p[0] = X[0] + X[1] - 1.0;
161: return PETSC_SUCCESS;
162: }
164: static PetscErrorCode quadratic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
165: {
166: T[0] = time + X[0] + X[1] + 1.0;
167: return PETSC_SUCCESS;
168: }
169: static PetscErrorCode quadratic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
170: {
171: T[0] = 1.0;
172: return PETSC_SUCCESS;
173: }
175: 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[])
176: {
177: const PetscReal nu = PetscRealPart(constants[NU]);
179: 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;
180: 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;
181: }
183: 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[])
184: {
185: f0[0] -= 2 * t + 1 + 3 * X[0] * X[0] - 2 * X[0] * X[1] + X[1] * X[1];
186: }
188: /*
189: CASE: quadratic
190: In 2D we use exact solution:
192: u = t + x^2 + y^2
193: v = t + 2x^2 - 2xy
194: p = x + y - 1
195: T = t + x + y + 1
196: rho = p^{th} / T
198: so that
200: \nabla \cdot u = 2x - 2x = 0
201: grad u = <<2 x, 4x - 2y>, <2 y, -2x>>
202: epsilon(u) = 1/2 (grad u + grad u^T) = <<2x, 2x>, <2x, -2x>>
203: epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
204: div epsilon'(u) = <2, 2>
206: f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
207: = 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>
208: = 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>
210: g = S rho_t + div (rho u)
211: = -S pth T_t/T^2 + rho div (u) + u . grad rho
212: = -S pth 1/T^2 - pth u . grad T / T^2
213: = -pth / T^2 (S + 2t + 3 x^2 - 2xy + y^2)
215: Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
216: = c_p S pth / T + c_p pth (2t + 3 x^2 - 2xy + y^2) / T - k/Pe 0
217: = c_p pth / T (S + 2t + 3 x^2 - 2xy + y^2)
218: */
219: 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[])
220: {
221: const PetscReal S = PetscRealPart(constants[STROUHAL]);
222: const PetscReal F = PetscRealPart(constants[FROUDE]);
223: const PetscReal Re = PetscRealPart(constants[REYNOLDS]);
224: const PetscReal mu = PetscRealPart(constants[MU]);
225: const PetscReal p_th = PetscRealPart(constants[P_TH]);
226: const PetscReal rho = p_th / (t + X[0] + X[1] + 1.);
227: const PetscInt gd = (PetscInt)PetscRealPart(constants[G_DIR]);
229: 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.;
230: 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.;
231: f0[gd] -= rho / PetscSqr(F);
232: }
234: 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[])
235: {
236: const PetscReal S = PetscRealPart(constants[STROUHAL]);
237: const PetscReal p_th = PetscRealPart(constants[P_TH]);
239: 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.);
240: }
242: 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[])
243: {
244: const PetscReal S = PetscRealPart(constants[STROUHAL]);
245: const PetscReal c_p = PetscRealPart(constants[C_P]);
246: const PetscReal p_th = PetscRealPart(constants[P_TH]);
248: 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.);
249: }
251: /*
252: CASE: cubic
253: In 2D we use exact solution:
255: u = t + x^3 + y^3
256: v = t + 2x^3 - 3x^2y
257: p = 3/2 x^2 + 3/2 y^2 - 1
258: T = t + 1/2 x^2 + 1/2 y^2
259: f = < t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x + 1,
260: t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y + 1>
261: Q = x^4 + xy^3 + 2x^3y - 3x^2y^2 + xt + yt - 2\alpha + 1
263: so that
265: \nabla \cdot u = 3x^2 - 3x^2 = 0
267: du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p - f
268: = <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
270: 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
271: */
272: static PetscErrorCode cubic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
273: {
274: u[0] = time + X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
275: u[1] = time + 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
276: return PETSC_SUCCESS;
277: }
278: static PetscErrorCode cubic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
279: {
280: u[0] = 1.0;
281: u[1] = 1.0;
282: return PETSC_SUCCESS;
283: }
285: static PetscErrorCode cubic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
286: {
287: p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
288: return PETSC_SUCCESS;
289: }
291: static PetscErrorCode cubic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
292: {
293: T[0] = time + X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
294: return PETSC_SUCCESS;
295: }
296: static PetscErrorCode cubic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
297: {
298: T[0] = 1.0;
299: return PETSC_SUCCESS;
300: }
302: 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[])
303: {
304: const PetscReal nu = PetscRealPart(constants[NU]);
306: 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);
307: 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);
308: }
310: 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[])
311: {
312: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
314: 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;
315: }
317: /*
318: CASE: cubic-trigonometric
319: In 2D we use exact solution:
321: u = beta cos t + x^3 + y^3
322: v = beta sin t + 2x^3 - 3x^2y
323: p = 3/2 x^2 + 3/2 y^2 - 1
324: T = 20 cos t + 1/2 x^2 + 1/2 y^2
325: f = < beta cos t 3x^2 + beta sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x,
326: beta cos t (6x^2 - 6xy) - beta sin t (3x^2) + 3x^4y + 6x^2y^3 - 6xy^4 - \nu(12x - 6y) + 3y>
327: Q = beta cos t x + beta sin t (y - 1) + x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2\alpha
329: so that
331: \nabla \cdot u = 3x^2 - 3x^2 = 0
333: f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
334: = <-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>
335: = <-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>
336: = <cos t (3x^2) + sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu (6x + 6y) + 3x,
337: cos t (6x^2 - 6xy) - sin t (3x^2) + 3x^4y + 6x^2y^3 - 6xy^4 - \nu (12x - 6y) + 3y>
339: Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
340: = -sin t + <cos t + x^3 + y^3, sin t + 2x^3 - 3x^2y> . <x, y> - 2 \alpha
341: = -sin t + cos t (x) + x^4 + xy^3 + sin t (y) + 2x^3y - 3x^2y^2 - 2 \alpha
342: = cos t x + sin t (y - 1) + (x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2 \alpha)
343: */
344: static PetscErrorCode cubic_trig_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
345: {
346: u[0] = 100. * PetscCosReal(time) + X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
347: u[1] = 100. * PetscSinReal(time) + 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
348: return PETSC_SUCCESS;
349: }
350: static PetscErrorCode cubic_trig_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
351: {
352: u[0] = -100. * PetscSinReal(time);
353: u[1] = 100. * PetscCosReal(time);
354: return PETSC_SUCCESS;
355: }
357: static PetscErrorCode cubic_trig_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
358: {
359: p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
360: return PETSC_SUCCESS;
361: }
363: static PetscErrorCode cubic_trig_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
364: {
365: T[0] = 100. * PetscCosReal(time) + X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
366: return PETSC_SUCCESS;
367: }
368: static PetscErrorCode cubic_trig_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
369: {
370: T[0] = -100. * PetscSinReal(time);
371: return PETSC_SUCCESS;
372: }
374: 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[])
375: {
376: const PetscReal nu = PetscRealPart(constants[NU]);
378: 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];
379: 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];
380: }
382: 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[])
383: {
384: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
386: 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;
387: }
389: /*
390: CASE: Taylor-Green vortex
391: In 2D we use exact solution:
393: u = 1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)
394: v = 1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)
395: p = -1/4 [cos(2 \pi(x - t)) + cos(2 \pi(y - t))] exp(-4 \pi^2 \nu t)
396: T = t + x + y
397: 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)) >
398: Q = 3 + sin(\pi(x-y)) exp(-2\nu \pi^2 t)
400: so that
402: \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
404: f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p
405: = <-\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),
406: \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)>
407: + < \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),
408: \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)>
409: + <-\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),
410: -\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)>
411: + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
412: 2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
413: + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
414: \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
415: = <-\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),
416: \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)>
417: + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
418: \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
419: + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
420: -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
421: + <-\pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
422: -\pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
423: + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
424: 2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
425: + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t),
426: \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)>
427: = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t),
428: \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)>
429: + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t),
430: \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)>
431: + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t),
432: -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)>
433: = < \pi cos(\pi(x - t)) cos(\pi(y - t)),
434: \pi sin(\pi(x - t)) sin(\pi(y - t))>
435: + <-\pi cos(\pi(x - t)) cos(\pi(y - t)),
436: -\pi sin(\pi(x - t)) sin(\pi(y - t))> = 0
437: Q = dT/dt + u \cdot \nabla T - \alpha \Delta T
438: = 1 + u \cdot <1, 1> - 0
439: = 1 + u + v
440: */
442: static PetscErrorCode taylor_green_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
443: {
444: u[0] = 1 - PetscCosReal(PETSC_PI * (X[0] - time)) * PetscSinReal(PETSC_PI * (X[1] - time)) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
445: u[1] = 1 + PetscSinReal(PETSC_PI * (X[0] - time)) * PetscCosReal(PETSC_PI * (X[1] - time)) * PetscExpReal(-2 * PETSC_PI * PETSC_PI * time);
446: return PETSC_SUCCESS;
447: }
448: static PetscErrorCode taylor_green_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
449: {
450: 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);
451: 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);
452: return PETSC_SUCCESS;
453: }
455: static PetscErrorCode taylor_green_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
456: {
457: 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);
458: return PETSC_SUCCESS;
459: }
461: static PetscErrorCode taylor_green_p_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
462: {
463: 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);
464: return PETSC_SUCCESS;
465: }
467: static PetscErrorCode taylor_green_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
468: {
469: T[0] = time + X[0] + X[1];
470: return PETSC_SUCCESS;
471: }
472: static PetscErrorCode taylor_green_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
473: {
474: T[0] = 1.0;
475: return PETSC_SUCCESS;
476: }
478: 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[])
479: {
480: PetscScalar vel[2];
482: PetscCallAbort(PETSC_COMM_SELF, taylor_green_u(dim, t, X, Nf, vel, NULL));
483: f0[0] -= 1.0 + vel[0] + vel[1];
484: }
486: /*
487: CASE: Pipe flow
488: Poiseuille flow, with the incoming fluid having a parabolic temperature profile and the side walls being held at T_in
490: u = \Delta Re/(2 mu) y (1 - y)
491: v = 0
492: p = -\Delta x
493: T = y (1 - y) + T_in
494: rho = p^{th} / T
496: so that
498: \nabla \cdot u = 0 - 0 = 0
499: grad u = \Delta Re/(2 mu) <<0, 0>, <1 - 2y, 0>>
500: epsilon(u) = 1/2 (grad u + grad u^T) = \Delta Re/(4 mu) <<0, 1 - 2y>, <<1 - 2y, 0>>
501: epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
502: div epsilon'(u) = -\Delta Re/(2 mu) <1, 0>
504: f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
505: = 0 + 0 - div (2\mu/Re \epsilon'(u) - pI) + rho / F^2 \hat y
506: = -\Delta div <<x, (1 - 2y)/2>, <<(1 - 2y)/2, x>> + rho / F^2 \hat y
507: = \Delta <1, 0> - \Delta <1, 0> + rho/F^2 <0, 1>
508: = rho/F^2 <0, 1>
510: g = S rho_t + div (rho u)
511: = 0 + rho div (u) + u . grad rho
512: = 0 + 0 - pth u . grad T / T^2
513: = 0
515: Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
516: = 0 + c_p pth / T 0 + 2 k/Pe
517: = 2 k/Pe
519: The boundary conditions on the top and bottom are zero velocity and T_in temperature. The boundary term is
521: (2\mu/Re \epsilon'(u) - p I) . n = \Delta <<x, (1 - 2y)/2>, <<(1 - 2y)/2, x>> . n
523: so that
525: x = 0: \Delta <<0, (1 - 2y)/2>, <<(1 - 2y)/2, 0>> . <-1, 0> = <0, (2y - 1)/2>
526: x = 1: \Delta <<1, (1 - 2y)/2>, <<(1 - 2y)/2, 1>> . <1, 0> = <1, (1 - 2y)/2>
527: */
528: static PetscErrorCode pipe_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
529: {
530: Parameter *param = (Parameter *)ctx;
532: u[0] = (0.5 * param->Reynolds / param->mu) * X[1] * (1.0 - X[1]);
533: u[1] = 0.0;
534: return PETSC_SUCCESS;
535: }
536: static PetscErrorCode pipe_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
537: {
538: u[0] = 0.0;
539: u[1] = 0.0;
540: return PETSC_SUCCESS;
541: }
543: static PetscErrorCode pipe_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
544: {
545: p[0] = -X[0];
546: return PETSC_SUCCESS;
547: }
548: static PetscErrorCode pipe_p_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
549: {
550: p[0] = 0.0;
551: return PETSC_SUCCESS;
552: }
554: static PetscErrorCode pipe_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
555: {
556: Parameter *param = (Parameter *)ctx;
558: T[0] = X[1] * (1.0 - X[1]) + param->T_in;
559: return PETSC_SUCCESS;
560: }
561: static PetscErrorCode pipe_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
562: {
563: T[0] = 0.0;
564: return PETSC_SUCCESS;
565: }
567: 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[])
568: {
569: const PetscReal F = PetscRealPart(constants[FROUDE]);
570: const PetscReal p_th = PetscRealPart(constants[P_TH]);
571: const PetscReal T_in = PetscRealPart(constants[T_IN]);
572: const PetscReal rho = p_th / (X[1] * (1. - X[1]) + T_in);
573: const PetscInt gd = (PetscInt)PetscRealPart(constants[G_DIR]);
575: f0[gd] -= rho / PetscSqr(F);
576: }
578: 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[])
579: {
580: PetscReal sigma[4] = {X[0], (PetscReal)(0.5 * (1. - 2. * X[1])), (PetscReal)(0.5 * (1. - 2. * X[1])), X[0]};
581: PetscInt d, e;
583: for (d = 0; d < dim; ++d) {
584: for (e = 0; e < dim; ++e) f0[d] -= sigma[d * dim + e] * n[e];
585: }
586: }
588: 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[])
589: {
590: f0[0] += 0.0;
591: }
593: 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[])
594: {
595: const PetscReal k = PetscRealPart(constants[K]);
596: const PetscReal Pe = PetscRealPart(constants[PECLET]);
598: f0[0] -= 2 * k / Pe;
599: }
601: /*
602: CASE: Wiggly pipe flow
603: Perturbed Poiseuille flow, with the incoming fluid having a perturbed parabolic temperature profile and the side walls being held at T_in
605: u = \Delta Re/(2 mu) [y (1 - y) + a sin(pi y)]
606: v = 0
607: p = -\Delta x
608: T = y (1 - y) + a sin(pi y) + T_in
609: rho = p^{th} / T
611: so that
613: \nabla \cdot u = 0 - 0 = 0
614: grad u = \Delta Re/(2 mu) <<0, 0>, <1 - 2y + a pi cos(pi y), 0>>
615: 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>>
616: epsilon'(u) = epsilon(u) - 1/3 (div u) I = epsilon(u)
617: div epsilon'(u) = -\Delta Re/(2 mu) <1 + a pi^2/2 sin(pi y), 0>
619: f = rho S du/dt + rho u \cdot \nabla u - 2\mu/Re div \epsilon'(u) + \nabla p + rho / F^2 \hat y
620: = 0 + 0 - div (2\mu/Re \epsilon'(u) - pI) + rho / F^2 \hat y
621: = -\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
622: = -\Delta <1 - 1 - a pi^2/2 sin(pi y), 0> + rho/F^2 <0, 1>
623: = a \Delta pi^2/2 sin(pi y) <1, 0> + rho/F^2 <0, 1>
625: g = S rho_t + div (rho u)
626: = 0 + rho div (u) + u . grad rho
627: = 0 + 0 - pth u . grad T / T^2
628: = 0
630: Q = rho c_p S dT/dt + rho c_p u . grad T - 1/Pe div k grad T
631: = 0 + c_p pth / T 0 - k/Pe div <0, 1 - 2y + a pi cos(pi y)>
632: = - k/Pe (-2 - a pi^2 sin(pi y))
633: = 2 k/Pe (1 + a pi^2/2 sin(pi y))
635: The boundary conditions on the top and bottom are zero velocity and T_in temperature. The boundary term is
637: (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
639: so that
641: x = 0: \Delta <<0, (1 - 2y)/2>, <<(1 - 2y)/2, 0>> . <-1, 0> = <0, (2y - 1)/2 - a pi/2 cos(pi y)>
642: x = 1: \Delta <<1, (1 - 2y)/2>, <<(1 - 2y)/2, 1>> . < 1, 0> = <1, (1 - 2y)/2 + a pi/2 cos(pi y)>
643: */
644: static PetscErrorCode pipe_wiggly_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
645: {
646: Parameter *param = (Parameter *)ctx;
648: u[0] = (0.5 * param->Reynolds / param->mu) * (X[1] * (1.0 - X[1]) + param->epsilon * PetscSinReal(PETSC_PI * X[1]));
649: u[1] = 0.0;
650: return PETSC_SUCCESS;
651: }
652: static PetscErrorCode pipe_wiggly_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
653: {
654: u[0] = 0.0;
655: u[1] = 0.0;
656: return PETSC_SUCCESS;
657: }
659: static PetscErrorCode pipe_wiggly_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
660: {
661: p[0] = -X[0];
662: return PETSC_SUCCESS;
663: }
664: static PetscErrorCode pipe_wiggly_p_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
665: {
666: p[0] = 0.0;
667: return PETSC_SUCCESS;
668: }
670: static PetscErrorCode pipe_wiggly_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
671: {
672: Parameter *param = (Parameter *)ctx;
674: T[0] = X[1] * (1.0 - X[1]) + param->epsilon * PetscSinReal(PETSC_PI * X[1]) + param->T_in;
675: return PETSC_SUCCESS;
676: }
677: static PetscErrorCode pipe_wiggly_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
678: {
679: T[0] = 0.0;
680: return PETSC_SUCCESS;
681: }
683: 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[])
684: {
685: const PetscReal F = PetscRealPart(constants[FROUDE]);
686: const PetscReal p_th = PetscRealPart(constants[P_TH]);
687: const PetscReal T_in = PetscRealPart(constants[T_IN]);
688: const PetscReal eps = PetscRealPart(constants[EPSILON]);
689: const PetscReal rho = p_th / (X[1] * (1. - X[1]) + T_in);
690: const PetscInt gd = (PetscInt)PetscRealPart(constants[G_DIR]);
692: f0[0] -= eps * 0.5 * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * X[1]);
693: f0[gd] -= rho / PetscSqr(F);
694: }
696: 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[])
697: {
698: const PetscReal eps = PetscRealPart(constants[EPSILON]);
699: 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]};
700: PetscInt d, e;
702: for (d = 0; d < dim; ++d) {
703: for (e = 0; e < dim; ++e) f0[d] -= sigma[d * dim + e] * n[e];
704: }
705: }
707: 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[])
708: {
709: f0[0] += 0.0;
710: }
712: 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[])
713: {
714: const PetscReal k = PetscRealPart(constants[K]);
715: const PetscReal Pe = PetscRealPart(constants[PECLET]);
716: const PetscReal eps = PetscRealPart(constants[EPSILON]);
718: f0[0] -= 2 * k / Pe * (1.0 + eps * 0.5 * PetscSqr(PETSC_PI) * PetscSinReal(PETSC_PI * X[1]));
719: }
721: /* Physics Kernels */
723: 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[])
724: {
725: PetscInt d;
726: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
727: }
729: /* -\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 */
730: 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[])
731: {
732: const PetscReal S = PetscRealPart(constants[STROUHAL]);
733: const PetscReal p_th = PetscRealPart(constants[P_TH]);
734: PetscInt d;
736: // -\frac{S p^{th}}{T^2} \frac{\partial T}{\partial t}
737: f0[0] += -u_t[uOff[TEMP]] * S * p_th / PetscSqr(u[uOff[TEMP]]);
739: // \frac{p^{th}}{T} \nabla \cdot \vb{u}
740: for (d = 0; d < dim; ++d) f0[0] += p_th / u[uOff[TEMP]] * u_x[uOff_x[VEL] + d * dim + d];
742: // - \frac{p^{th}}{T^2} \vb{u} \cdot \nabla T
743: 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];
745: // Add in any fixed source term
746: if (NfAux > 0) f0[0] += a[aOff[MASS]];
747: }
749: /* \vb{u}_t + \vb{u} \cdot \nabla\vb{u} */
750: 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[])
751: {
752: const PetscInt Nc = dim;
753: PetscInt c, d;
755: for (c = 0; c < Nc; ++c) {
756: /* \vb{u}_t */
757: f0[c] += u_t[uOff[VEL] + c];
758: /* \vb{u} \cdot \nabla\vb{u} */
759: for (d = 0; d < dim; ++d) f0[c] += u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
760: }
761: }
763: /* \rho S \frac{\partial \vb{u}}{\partial t} + \rho \vb{u} \cdot \nabla \vb{u} + \rho \frac{\hat{\vb{z}}}{F^2} */
764: 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[])
765: {
766: const PetscReal S = PetscRealPart(constants[STROUHAL]);
767: const PetscReal F = PetscRealPart(constants[FROUDE]);
768: const PetscReal p_th = PetscRealPart(constants[P_TH]);
769: const PetscReal rho = p_th / PetscRealPart(u[uOff[TEMP]]);
770: const PetscInt gdir = (PetscInt)PetscRealPart(constants[G_DIR]);
771: PetscInt Nc = dim;
772: PetscInt c, d;
774: // \rho S \frac{\partial \vb{u}}{\partial t}
775: for (d = 0; d < dim; ++d) f0[d] = rho * S * u_t[uOff[VEL] + d];
777: // \rho \vb{u} \cdot \nabla \vb{u}
778: for (c = 0; c < Nc; ++c) {
779: for (d = 0; d < dim; ++d) f0[c] += rho * u[uOff[VEL] + d] * u_x[uOff_x[VEL] + c * dim + d];
780: }
782: // rho \hat{z}/F^2
783: f0[gdir] += rho / (F * F);
785: // Add in any fixed source term
786: if (NfAux > 0) {
787: for (d = 0; d < dim; ++d) f0[d] += a[aOff[MOMENTUM] + d];
788: }
789: }
791: /*f1_v = \nu[grad(u) + grad(u)^T] - pI */
792: 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[])
793: {
794: const PetscReal nu = PetscRealPart(constants[NU]);
795: const PetscInt Nc = dim;
796: PetscInt c, d;
798: for (c = 0; c < Nc; ++c) {
799: for (d = 0; d < dim; ++d) f1[c * dim + d] = nu * (u_x[c * dim + d] + u_x[d * dim + c]);
800: f1[c * dim + c] -= u[uOff[1]];
801: }
802: }
804: /* 2 \mu/Re (1/2 (\nabla \vb{u} + \nabla \vb{u}^T) - 1/3 (\nabla \cdot \vb{u}) I) - p I */
805: 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[])
806: {
807: const PetscReal Re = PetscRealPart(constants[REYNOLDS]);
808: const PetscReal mu = PetscRealPart(constants[MU]);
809: const PetscReal coef = mu / Re;
810: PetscReal u_div = 0.0;
811: const PetscInt Nc = dim;
812: PetscInt c, d;
814: for (c = 0; c < Nc; ++c) u_div += PetscRealPart(u_x[uOff_x[VEL] + c * dim + c]);
816: for (c = 0; c < Nc; ++c) {
817: // 2 \mu/Re 1/2 (\nabla \vb{u} + \nabla \vb{u}^T
818: 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]);
819: // -2/3 \mu/Re (\nabla \cdot \vb{u}) I
820: f1[c * dim + c] -= 2.0 * coef / 3.0 * u_div;
821: }
823: // -p I
824: for (c = 0; c < Nc; ++c) f1[c * dim + c] -= u[uOff[PRES]];
825: }
827: /* T_t + \vb{u} \cdot \nabla T */
828: 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[])
829: {
830: PetscInt d;
832: /* T_t */
833: f0[0] += u_t[uOff[TEMP]];
834: /* \vb{u} \cdot \nabla T */
835: for (d = 0; d < dim; ++d) f0[0] += u[uOff[VEL] + d] * u_x[uOff_x[TEMP] + d];
836: }
838: /* \frac{C_p S p^{th}}{T} \frac{\partial T}{\partial t} + \frac{C_p p^{th}}{T} \vb{u} \cdot \nabla T */
839: 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[])
840: {
841: const PetscReal S = PetscRealPart(constants[STROUHAL]);
842: const PetscReal c_p = PetscRealPart(constants[C_P]);
843: const PetscReal p_th = PetscRealPart(constants[P_TH]);
844: PetscInt d;
846: // \frac{C_p S p^{th}}{T} \frac{\partial T}{\partial t}
847: f0[0] = c_p * S * p_th / u[uOff[TEMP]] * u_t[uOff[TEMP]];
849: // \frac{C_p p^{th}}{T} \vb{u} \cdot \nabla T
850: 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];
852: // Add in any fixed source term
853: if (NfAux > 0) f0[0] += a[aOff[ENERGY]];
854: }
856: 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[])
857: {
858: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
859: PetscInt d;
861: for (d = 0; d < dim; ++d) f1[d] = alpha * u_x[uOff_x[2] + d];
862: }
864: /* \frac{k}{Pe} \nabla T */
865: 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[])
866: {
867: const PetscReal Pe = PetscRealPart(constants[PECLET]);
868: const PetscReal k = PetscRealPart(constants[K]);
869: PetscInt d;
871: // \frac{k}{Pe} \nabla T
872: for (d = 0; d < dim; ++d) f1[d] = k / Pe * u_x[uOff_x[TEMP] + d];
873: }
875: 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[])
876: {
877: PetscInt d;
878: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
879: }
881: 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[])
882: {
883: PetscInt c, d;
884: const PetscInt Nc = dim;
886: for (d = 0; d < dim; ++d) g0[d * dim + d] = u_tShift;
888: for (c = 0; c < Nc; ++c) {
889: for (d = 0; d < dim; ++d) g0[c * Nc + d] += u_x[c * Nc + d];
890: }
891: }
893: 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[])
894: {
895: PetscInt NcI = dim;
896: PetscInt NcJ = dim;
897: PetscInt c, d, e;
899: for (c = 0; c < NcI; ++c) {
900: for (d = 0; d < NcJ; ++d) {
901: for (e = 0; e < dim; ++e) {
902: if (c == d) g1[(c * NcJ + d) * dim + e] += u[e];
903: }
904: }
905: }
906: }
908: 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[])
909: {
910: const PetscReal p_th = PetscRealPart(constants[P_TH]);
911: PetscInt d;
913: // - \phi_i \frac{p^{th}}{T^2} \frac{\partial T}{\partial x_c} \psi_{j, u_c}
914: for (d = 0; d < dim; ++d) g0[d] = -p_th / PetscSqr(u[uOff[TEMP]]) * u_x[uOff_x[TEMP] + d];
915: }
917: 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[])
918: {
919: const PetscReal p_th = PetscRealPart(constants[P_TH]);
920: PetscInt d;
922: // \phi_i \frac{p^{th}}{T} \frac{\partial \psi_{u_c,j}}{\partial x_c}
923: for (d = 0; d < dim; ++d) g1[d * dim + d] = p_th / u[uOff[TEMP]];
924: }
926: 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[])
927: {
928: const PetscReal S = PetscRealPart(constants[STROUHAL]);
929: const PetscReal p_th = PetscRealPart(constants[P_TH]);
930: PetscInt d;
932: // - \phi_i \frac{S p^{th}}{T^2} \psi_j
933: g0[0] -= S * p_th / PetscSqr(u[uOff[TEMP]]) * u_tShift;
934: // \phi_i 2 \frac{S p^{th}}{T^3} T_t \psi_j
935: g0[0] += 2.0 * S * p_th / PetscPowScalarInt(u[uOff[TEMP]], 3) * u_t[uOff[TEMP]];
936: // \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)
937: 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]);
938: }
940: 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[])
941: {
942: const PetscReal p_th = PetscRealPart(constants[P_TH]);
943: PetscInt d;
945: // - \phi_i \frac{p^{th}}{T^2} \vb{u} \cdot \nabla \psi_j
946: for (d = 0; d < dim; ++d) g1[d] = -p_th / PetscSqr(u[uOff[TEMP]]) * u[uOff[VEL] + d];
947: }
949: 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[])
950: {
951: PetscInt d;
952: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
953: }
955: 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[])
956: {
957: const PetscReal nu = PetscRealPart(constants[NU]);
958: const PetscInt Nc = dim;
959: PetscInt c, d;
961: for (c = 0; c < Nc; ++c) {
962: for (d = 0; d < dim; ++d) {
963: g3[((c * Nc + c) * dim + d) * dim + d] += nu;
964: g3[((c * Nc + d) * dim + d) * dim + c] += nu;
965: }
966: }
967: }
969: 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[])
970: {
971: const PetscReal S = PetscRealPart(constants[STROUHAL]);
972: const PetscReal F = PetscRealPart(constants[FROUDE]);
973: const PetscReal p_th = PetscRealPart(constants[P_TH]);
974: const PetscInt gdir = (PetscInt)PetscRealPart(constants[G_DIR]);
975: const PetscInt Nc = dim;
976: PetscInt c, d;
978: // - \vb{\phi}_i \cdot \vb{u}_t \frac{p^{th} S}{T^2} \psi_j
979: for (d = 0; d < dim; ++d) g0[d] -= p_th * S / PetscSqr(u[uOff[TEMP]]) * u_t[uOff[VEL] + d];
981: // - \vb{\phi}_i \cdot \vb{u} \cdot \nabla \vb{u} \frac{p^{th}}{T^2} \psi_j
982: for (c = 0; c < Nc; ++c) {
983: 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];
984: }
986: // - \vb{\phi}_i \cdot \vu{z} \frac{p^{th}}{T^2 F^2} \psi_j
987: g0[gdir] -= p_th / PetscSqr(u[uOff[TEMP]] * F);
988: }
990: 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[])
991: {
992: const PetscReal S = PetscRealPart(constants[STROUHAL]);
993: const PetscReal p_th = PetscRealPart(constants[P_TH]);
994: const PetscInt Nc = dim;
995: PetscInt c, d;
997: // \vb{\phi}_i \cdot S \rho \psi_j
998: for (d = 0; d < dim; ++d) g0[d * dim + d] = S * p_th / u[uOff[TEMP]] * u_tShift;
1000: // \phi^c_i \cdot \rho \frac{\partial u^c}{\partial x^d} \psi^d_j
1001: for (c = 0; c < Nc; ++c) {
1002: for (d = 0; d < dim; ++d) g0[c * Nc + d] += p_th / u[uOff[TEMP]] * u_x[uOff_x[VEL] + c * Nc + d];
1003: }
1004: }
1006: 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[])
1007: {
1008: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1009: const PetscInt NcI = dim;
1010: const PetscInt NcJ = dim;
1011: PetscInt c, d, e;
1013: // \phi^c_i \rho u^e \frac{\partial \psi^d_j}{\partial x^e}
1014: for (c = 0; c < NcI; ++c) {
1015: for (d = 0; d < NcJ; ++d) {
1016: for (e = 0; e < dim; ++e) {
1017: if (c == d) g1[(c * NcJ + d) * dim + e] += p_th / u[uOff[TEMP]] * u[uOff[VEL] + e];
1018: }
1019: }
1020: }
1021: }
1023: 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[])
1024: {
1025: const PetscReal Re = PetscRealPart(constants[REYNOLDS]);
1026: const PetscReal mu = PetscRealPart(constants[MU]);
1027: const PetscInt Nc = dim;
1028: PetscInt c, d;
1030: for (c = 0; c < Nc; ++c) {
1031: for (d = 0; d < dim; ++d) {
1032: // \frac{\partial \phi^c_i}{\partial x^d} \mu/Re \frac{\partial \psi^c_i}{\partial x^d}
1033: g3[((c * Nc + c) * dim + d) * dim + d] += mu / Re; // gradU
1034: // \frac{\partial \phi^c_i}{\partial x^d} \mu/Re \frac{\partial \psi^d_i}{\partial x^c}
1035: g3[((c * Nc + d) * dim + d) * dim + c] += mu / Re; // gradU transpose
1036: // \frac{\partial \phi^c_i}{\partial x^d} -2/3 \mu/Re \frac{\partial \psi^d_i}{\partial x^c}
1037: g3[((c * Nc + d) * dim + c) * dim + d] -= 2.0 / 3.0 * mu / Re;
1038: }
1039: }
1040: }
1042: 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[])
1043: {
1044: PetscInt d;
1045: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
1046: }
1048: 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[])
1049: {
1050: g0[0] = u_tShift;
1051: }
1053: 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[])
1054: {
1055: PetscInt d;
1056: for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2] + d];
1057: }
1059: 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[])
1060: {
1061: PetscInt d;
1062: for (d = 0; d < dim; ++d) g1[d] = u[uOff[0] + d];
1063: }
1065: 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[])
1066: {
1067: const PetscReal alpha = PetscRealPart(constants[ALPHA]);
1068: PetscInt d;
1070: for (d = 0; d < dim; ++d) g3[d * dim + d] = alpha;
1071: }
1073: 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[])
1074: {
1075: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1076: const PetscReal c_p = PetscRealPart(constants[C_P]);
1077: PetscInt d;
1079: // \phi_i \frac{C_p p^{th}}{T} \nabla T \cdot \psi_j
1080: for (d = 0; d < dim; ++d) g0[d] = c_p * p_th / u[uOff[TEMP]] * u_x[uOff_x[TEMP] + d];
1081: }
1083: 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[])
1084: {
1085: const PetscReal S = PetscRealPart(constants[STROUHAL]);
1086: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1087: const PetscReal c_p = PetscRealPart(constants[C_P]);
1088: PetscInt d;
1090: // \psi_i C_p S p^{th}\T \psi_{j}
1091: g0[0] += c_p * S * p_th / u[uOff[TEMP]] * u_tShift;
1092: // - \phi_i C_p S p^{th}/T^2 T_t \psi_j
1093: g0[0] -= c_p * S * p_th / PetscSqr(u[uOff[TEMP]]) * u_t[uOff[TEMP]];
1094: // - \phi_i C_p p^{th}/T^2 \vb{u} \cdot \nabla T \psi_j
1095: for (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];
1096: }
1098: 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[])
1099: {
1100: const PetscReal p_th = PetscRealPart(constants[P_TH]);
1101: const PetscReal c_p = PetscRealPart(constants[C_P]);
1102: PetscInt d;
1104: // \phi_i C_p p^{th}/T \vb{u} \cdot \nabla \psi_j
1105: for (d = 0; d < dim; ++d) g1[d] += c_p * p_th / u[uOff[TEMP]] * u[uOff[VEL] + d];
1106: }
1108: 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[])
1109: {
1110: const PetscReal Pe = PetscRealPart(constants[PECLET]);
1111: const PetscReal k = PetscRealPart(constants[K]);
1112: PetscInt d;
1114: // \nabla \phi_i \frac{k}{Pe} \nabla \phi_j
1115: for (d = 0; d < dim; ++d) g3[d * dim + d] = k / Pe;
1116: }
1118: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
1119: {
1120: PetscInt mod, sol;
1122: PetscFunctionBeginUser;
1123: options->modType = MOD_INCOMPRESSIBLE;
1124: options->solType = SOL_QUADRATIC;
1125: options->hasNullSpace = PETSC_TRUE;
1126: options->dmCell = NULL;
1128: PetscOptionsBegin(comm, "", "Low Mach flow Problem Options", "DMPLEX");
1129: mod = options->modType;
1130: PetscCall(PetscOptionsEList("-mod_type", "The model type", "ex76.c", modTypes, NUM_MOD_TYPES, modTypes[options->modType], &mod, NULL));
1131: options->modType = (ModType)mod;
1132: sol = options->solType;
1133: PetscCall(PetscOptionsEList("-sol_type", "The solution type", "ex76.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL));
1134: options->solType = (SolType)sol;
1135: PetscOptionsEnd();
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }
1139: static PetscErrorCode SetupParameters(DM dm, AppCtx *user)
1140: {
1141: PetscBag bag;
1142: Parameter *p;
1143: PetscReal dir;
1144: PetscInt dim;
1146: PetscFunctionBeginUser;
1147: PetscCall(DMGetDimension(dm, &dim));
1148: dir = (PetscReal)(dim - 1);
1149: /* setup PETSc parameter bag */
1150: PetscCall(PetscBagGetData(user->bag, (void **)&p));
1151: PetscCall(PetscBagSetName(user->bag, "par", "Low Mach flow parameters"));
1152: bag = user->bag;
1153: PetscCall(PetscBagRegisterReal(bag, &p->Strouhal, 1.0, "S", "Strouhal number"));
1154: PetscCall(PetscBagRegisterReal(bag, &p->Froude, 1.0, "Fr", "Froude number"));
1155: PetscCall(PetscBagRegisterReal(bag, &p->Reynolds, 1.0, "Re", "Reynolds number"));
1156: PetscCall(PetscBagRegisterReal(bag, &p->Peclet, 1.0, "Pe", "Peclet number"));
1157: PetscCall(PetscBagRegisterReal(bag, &p->p_th, 1.0, "p_th", "Thermodynamic pressure"));
1158: PetscCall(PetscBagRegisterReal(bag, &p->mu, 1.0, "mu", "Dynamic viscosity"));
1159: PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
1160: PetscCall(PetscBagRegisterReal(bag, &p->c_p, 1.0, "c_p", "Specific heat at constant pressure"));
1161: PetscCall(PetscBagRegisterReal(bag, &p->k, 1.0, "k", "Thermal conductivity"));
1162: PetscCall(PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity"));
1163: PetscCall(PetscBagRegisterReal(bag, &p->T_in, 1.0, "T_in", "Inlet temperature"));
1164: PetscCall(PetscBagRegisterReal(bag, &p->g_dir, dir, "g_dir", "Gravity direction"));
1165: PetscCall(PetscBagRegisterReal(bag, &p->epsilon, 1.0, "epsilon", "Perturbation strength"));
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
1170: {
1171: PetscFunctionBeginUser;
1172: PetscCall(DMCreate(comm, dm));
1173: PetscCall(DMSetType(*dm, DMPLEX));
1174: PetscCall(DMSetFromOptions(*dm));
1175: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: static PetscErrorCode UniformBoundaryConditions(DM dm, DMLabel label, PetscSimplePointFn *exactFuncs[], PetscSimplePointFn *exactFuncs_t[], AppCtx *user)
1180: {
1181: PetscDS ds;
1182: PetscInt id;
1183: void *ctx;
1185: PetscFunctionBeginUser;
1186: PetscCall(DMGetDS(dm, &ds));
1187: PetscCall(PetscBagGetData(user->bag, &ctx));
1188: id = 3;
1189: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1190: id = 1;
1191: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1192: id = 2;
1193: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1194: id = 4;
1195: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1196: id = 3;
1197: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1198: id = 1;
1199: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1200: id = 2;
1201: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1202: id = 4;
1203: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1204: PetscFunctionReturn(PETSC_SUCCESS);
1205: }
1207: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
1208: {
1209: PetscSimplePointFn *exactFuncs[3];
1210: PetscSimplePointFn *exactFuncs_t[3];
1211: PetscDS ds;
1212: PetscWeakForm wf;
1213: DMLabel label;
1214: Parameter *ctx;
1215: PetscInt id, bd;
1217: PetscFunctionBeginUser;
1218: PetscCall(DMGetLabel(dm, "marker", &label));
1219: PetscCall(DMGetDS(dm, &ds));
1220: PetscCall(PetscDSGetWeakForm(ds, &wf));
1222: switch (user->modType) {
1223: case MOD_INCOMPRESSIBLE:
1224: PetscCall(PetscDSSetResidual(ds, VEL, f0_v, f1_v));
1225: PetscCall(PetscDSSetResidual(ds, PRES, f0_q, NULL));
1226: PetscCall(PetscDSSetResidual(ds, TEMP, f0_w, f1_w));
1228: PetscCall(PetscDSSetJacobian(ds, VEL, VEL, g0_vu, g1_vu, NULL, g3_vu));
1229: PetscCall(PetscDSSetJacobian(ds, VEL, PRES, NULL, NULL, g2_vp, NULL));
1230: PetscCall(PetscDSSetJacobian(ds, PRES, VEL, NULL, g1_qu, NULL, NULL));
1231: PetscCall(PetscDSSetJacobian(ds, TEMP, VEL, g0_wu, NULL, NULL, NULL));
1232: PetscCall(PetscDSSetJacobian(ds, TEMP, TEMP, g0_wT, g1_wT, NULL, g3_wT));
1234: switch (user->solType) {
1235: case SOL_QUADRATIC:
1236: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_quadratic_v, 0, NULL));
1237: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_quadratic_w, 0, NULL));
1239: exactFuncs[VEL] = quadratic_u;
1240: exactFuncs[PRES] = quadratic_p;
1241: exactFuncs[TEMP] = quadratic_T;
1242: exactFuncs_t[VEL] = quadratic_u_t;
1243: exactFuncs_t[PRES] = NULL;
1244: exactFuncs_t[TEMP] = quadratic_T_t;
1246: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1247: break;
1248: case SOL_CUBIC:
1249: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_cubic_v, 0, NULL));
1250: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_cubic_w, 0, NULL));
1252: exactFuncs[VEL] = cubic_u;
1253: exactFuncs[PRES] = cubic_p;
1254: exactFuncs[TEMP] = cubic_T;
1255: exactFuncs_t[VEL] = cubic_u_t;
1256: exactFuncs_t[PRES] = NULL;
1257: exactFuncs_t[TEMP] = cubic_T_t;
1259: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1260: break;
1261: case SOL_CUBIC_TRIG:
1262: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_cubic_trig_v, 0, NULL));
1263: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_cubic_trig_w, 0, NULL));
1265: exactFuncs[VEL] = cubic_trig_u;
1266: exactFuncs[PRES] = cubic_trig_p;
1267: exactFuncs[TEMP] = cubic_trig_T;
1268: exactFuncs_t[VEL] = cubic_trig_u_t;
1269: exactFuncs_t[PRES] = NULL;
1270: exactFuncs_t[TEMP] = cubic_trig_T_t;
1272: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1273: break;
1274: case SOL_TAYLOR_GREEN:
1275: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_taylor_green_w, 0, NULL));
1277: exactFuncs[VEL] = taylor_green_u;
1278: exactFuncs[PRES] = taylor_green_p;
1279: exactFuncs[TEMP] = taylor_green_T;
1280: exactFuncs_t[VEL] = taylor_green_u_t;
1281: exactFuncs_t[PRES] = taylor_green_p_t;
1282: exactFuncs_t[TEMP] = taylor_green_T_t;
1284: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1285: break;
1286: default:
1287: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
1288: }
1289: break;
1290: case MOD_CONDUCTING:
1291: PetscCall(PetscDSSetResidual(ds, VEL, f0_conduct_v, f1_conduct_v));
1292: PetscCall(PetscDSSetResidual(ds, PRES, f0_conduct_q, NULL));
1293: PetscCall(PetscDSSetResidual(ds, TEMP, f0_conduct_w, f1_conduct_w));
1295: PetscCall(PetscDSSetJacobian(ds, VEL, VEL, g0_conduct_vu, g1_conduct_vu, NULL, g3_conduct_vu));
1296: PetscCall(PetscDSSetJacobian(ds, VEL, PRES, NULL, NULL, g2_conduct_vp, NULL));
1297: PetscCall(PetscDSSetJacobian(ds, VEL, TEMP, g0_conduct_vT, NULL, NULL, NULL));
1298: PetscCall(PetscDSSetJacobian(ds, PRES, VEL, g0_conduct_qu, g1_conduct_qu, NULL, NULL));
1299: PetscCall(PetscDSSetJacobian(ds, PRES, TEMP, g0_conduct_qT, g1_conduct_qT, NULL, NULL));
1300: PetscCall(PetscDSSetJacobian(ds, TEMP, VEL, g0_conduct_wu, NULL, NULL, NULL));
1301: PetscCall(PetscDSSetJacobian(ds, TEMP, TEMP, g0_conduct_wT, g1_conduct_wT, NULL, g3_conduct_wT));
1303: switch (user->solType) {
1304: case SOL_QUADRATIC:
1305: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_quadratic_v, 0, NULL));
1306: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_quadratic_q, 0, NULL));
1307: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_quadratic_w, 0, NULL));
1309: exactFuncs[VEL] = quadratic_u;
1310: exactFuncs[PRES] = quadratic_p;
1311: exactFuncs[TEMP] = quadratic_T;
1312: exactFuncs_t[VEL] = quadratic_u_t;
1313: exactFuncs_t[PRES] = NULL;
1314: exactFuncs_t[TEMP] = quadratic_T_t;
1316: PetscCall(UniformBoundaryConditions(dm, label, exactFuncs, exactFuncs_t, user));
1317: break;
1318: case SOL_PIPE:
1319: user->hasNullSpace = PETSC_FALSE;
1320: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_pipe_v, 0, NULL));
1321: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_pipe_q, 0, NULL));
1322: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_pipe_w, 0, NULL));
1324: exactFuncs[VEL] = pipe_u;
1325: exactFuncs[PRES] = pipe_p;
1326: exactFuncs[TEMP] = pipe_T;
1327: exactFuncs_t[VEL] = pipe_u_t;
1328: exactFuncs_t[PRES] = pipe_p_t;
1329: exactFuncs_t[TEMP] = pipe_T_t;
1331: PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
1332: id = 2;
1333: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1334: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1335: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_v, 0, NULL));
1336: id = 4;
1337: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1338: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1339: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_v, 0, NULL));
1340: id = 4;
1341: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temperature", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1342: id = 3;
1343: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1344: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temperature", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1345: id = 1;
1346: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1347: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temperature", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1348: break;
1349: case SOL_PIPE_WIGGLY:
1350: user->hasNullSpace = PETSC_FALSE;
1351: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, VEL, 0, 1, f0_conduct_pipe_wiggly_v, 0, NULL));
1352: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, PRES, 0, 1, f0_conduct_pipe_wiggly_q, 0, NULL));
1353: PetscCall(PetscWeakFormSetIndexResidual(wf, NULL, 0, TEMP, 0, 1, f0_conduct_pipe_wiggly_w, 0, NULL));
1355: exactFuncs[VEL] = pipe_wiggly_u;
1356: exactFuncs[PRES] = pipe_wiggly_p;
1357: exactFuncs[TEMP] = pipe_wiggly_T;
1358: exactFuncs_t[VEL] = pipe_wiggly_u_t;
1359: exactFuncs_t[PRES] = pipe_wiggly_p_t;
1360: exactFuncs_t[TEMP] = pipe_wiggly_T_t;
1362: PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
1363: id = 2;
1364: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1365: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1366: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_wiggly_v, 0, NULL));
1367: id = 4;
1368: PetscCall(DMAddBoundary(dm, DM_BC_NATURAL, "left wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd));
1369: PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1370: PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, id, VEL, 0, 0, f0_conduct_bd_pipe_wiggly_v, 0, NULL));
1371: id = 4;
1372: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "left wall temperature", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1373: id = 3;
1374: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1375: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "top wall temperature", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1376: id = 1;
1377: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, VEL, 0, NULL, (void (*)(void))exactFuncs[VEL], (void (*)(void))exactFuncs_t[VEL], ctx, NULL));
1378: PetscCall(PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "bottom wall temperature", label, 1, &id, TEMP, 0, NULL, (void (*)(void))exactFuncs[TEMP], (void (*)(void))exactFuncs_t[TEMP], ctx, NULL));
1379: break;
1380: default:
1381: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
1382: }
1383: break;
1384: default:
1385: SETERRQ(PetscObjectComm((PetscObject)ds), PETSC_ERR_ARG_WRONG, "Unsupported model type: %s (%d)", modTypes[PetscMin(user->modType, NUM_MOD_TYPES)], user->modType);
1386: }
1387: /* Setup constants */
1388: {
1389: Parameter *param;
1390: PetscScalar constants[13];
1392: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1394: constants[STROUHAL] = param->Strouhal;
1395: constants[FROUDE] = param->Froude;
1396: constants[REYNOLDS] = param->Reynolds;
1397: constants[PECLET] = param->Peclet;
1398: constants[P_TH] = param->p_th;
1399: constants[MU] = param->mu;
1400: constants[NU] = param->nu;
1401: constants[C_P] = param->c_p;
1402: constants[K] = param->k;
1403: constants[ALPHA] = param->alpha;
1404: constants[T_IN] = param->T_in;
1405: constants[G_DIR] = param->g_dir;
1406: constants[EPSILON] = param->epsilon;
1407: PetscCall(PetscDSSetConstants(ds, 13, constants));
1408: }
1410: PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
1411: PetscCall(PetscDSSetExactSolution(ds, VEL, exactFuncs[VEL], ctx));
1412: PetscCall(PetscDSSetExactSolution(ds, PRES, exactFuncs[PRES], ctx));
1413: PetscCall(PetscDSSetExactSolution(ds, TEMP, exactFuncs[TEMP], ctx));
1414: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, VEL, exactFuncs_t[VEL], ctx));
1415: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, PRES, exactFuncs_t[PRES], ctx));
1416: PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, TEMP, exactFuncs_t[TEMP], ctx));
1417: PetscFunctionReturn(PETSC_SUCCESS);
1418: }
1420: static PetscErrorCode CreateCellDM(DM dm, AppCtx *user)
1421: {
1422: PetscFE fe, fediv;
1423: DMPolytopeType ct;
1424: PetscInt dim, cStart;
1425: PetscBool simplex;
1427: PetscFunctionBeginUser;
1428: PetscCall(DMGetDimension(dm, &dim));
1429: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
1430: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1431: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1433: PetscCall(DMGetField(dm, VEL, NULL, (PetscObject *)&fe));
1434: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "div_", PETSC_DEFAULT, &fediv));
1435: PetscCall(PetscFECopyQuadrature(fe, fediv));
1436: PetscCall(PetscObjectSetName((PetscObject)fediv, "divergence"));
1438: PetscCall(DMDestroy(&user->dmCell));
1439: PetscCall(DMClone(dm, &user->dmCell));
1440: PetscCall(DMSetField(user->dmCell, 0, NULL, (PetscObject)fediv));
1441: PetscCall(DMCreateDS(user->dmCell));
1442: PetscCall(PetscFEDestroy(&fediv));
1443: PetscFunctionReturn(PETSC_SUCCESS);
1444: }
1446: static PetscErrorCode GetCellDM(DM dm, AppCtx *user, DM *dmCell)
1447: {
1448: PetscInt cStart, cEnd, cellStart = -1, cellEnd = -1;
1450: PetscFunctionBeginUser;
1451: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1452: if (user->dmCell) PetscCall(DMPlexGetSimplexOrBoxCells(user->dmCell, 0, &cellStart, &cellEnd));
1453: if (cStart != cellStart || cEnd != cellEnd) PetscCall(CreateCellDM(dm, user));
1454: *dmCell = user->dmCell;
1455: PetscFunctionReturn(PETSC_SUCCESS);
1456: }
1458: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
1459: {
1460: DM cdm = dm;
1461: PetscFE fe[3];
1462: Parameter *param;
1463: DMPolytopeType ct;
1464: PetscInt dim, cStart;
1465: PetscBool simplex;
1467: PetscFunctionBeginUser;
1468: PetscCall(DMGetDimension(dm, &dim));
1469: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
1470: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
1471: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
1472: /* Create finite element */
1473: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]));
1474: PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
1476: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
1477: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
1478: PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
1480: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]));
1481: PetscCall(PetscFECopyQuadrature(fe[0], fe[2]));
1482: PetscCall(PetscObjectSetName((PetscObject)fe[2], "temperature"));
1484: /* Set discretization and boundary conditions for each mesh */
1485: PetscCall(DMSetField(dm, VEL, NULL, (PetscObject)fe[VEL]));
1486: PetscCall(DMSetField(dm, PRES, NULL, (PetscObject)fe[PRES]));
1487: PetscCall(DMSetField(dm, TEMP, NULL, (PetscObject)fe[TEMP]));
1488: PetscCall(DMCreateDS(dm));
1489: PetscCall(SetupProblem(dm, user));
1490: PetscCall(PetscBagGetData(user->bag, (void **)¶m));
1491: while (cdm) {
1492: PetscCall(DMCopyDisc(dm, cdm));
1493: PetscCall(DMGetCoarseDM(cdm, &cdm));
1494: }
1495: PetscCall(PetscFEDestroy(&fe[VEL]));
1496: PetscCall(PetscFEDestroy(&fe[PRES]));
1497: PetscCall(PetscFEDestroy(&fe[TEMP]));
1499: if (user->hasNullSpace) {
1500: PetscObject pressure;
1501: MatNullSpace nullspacePres;
1503: PetscCall(DMGetField(dm, PRES, NULL, &pressure));
1504: PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres));
1505: PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullspacePres));
1506: PetscCall(MatNullSpaceDestroy(&nullspacePres));
1507: }
1508: PetscFunctionReturn(PETSC_SUCCESS);
1509: }
1511: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
1512: {
1513: Vec vec;
1514: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero};
1516: PetscFunctionBeginUser;
1517: 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);
1518: funcs[nfield] = constant;
1519: PetscCall(DMCreateGlobalVector(dm, &vec));
1520: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
1521: PetscCall(VecNormalize(vec, NULL));
1522: PetscCall(PetscObjectSetName((PetscObject)vec, "Pressure Null Space"));
1523: PetscCall(VecViewFromOptions(vec, NULL, "-pressure_nullspace_view"));
1524: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
1525: PetscCall(VecDestroy(&vec));
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: static PetscErrorCode RemoveDiscretePressureNullspace_Private(TS ts, Vec u)
1530: {
1531: DM dm;
1532: AppCtx *user;
1533: MatNullSpace nullsp;
1535: PetscFunctionBeginUser;
1536: PetscCall(TSGetDM(ts, &dm));
1537: PetscCall(DMGetApplicationContext(dm, &user));
1538: if (!user->hasNullSpace) PetscFunctionReturn(PETSC_SUCCESS);
1539: PetscCall(CreatePressureNullSpace(dm, 1, 1, &nullsp));
1540: PetscCall(MatNullSpaceRemove(nullsp, u));
1541: PetscCall(MatNullSpaceDestroy(&nullsp));
1542: PetscFunctionReturn(PETSC_SUCCESS);
1543: }
1545: /* Make the discrete pressure discretely divergence free */
1546: static PetscErrorCode RemoveDiscretePressureNullspace(TS ts)
1547: {
1548: Vec u;
1550: PetscFunctionBeginUser;
1551: PetscCall(TSGetSolution(ts, &u));
1552: PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: 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[])
1557: {
1558: PetscInt d;
1560: divu[0] = 0.;
1561: for (d = 0; d < dim; ++d) divu[0] += u_x[d * dim + d];
1562: }
1564: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
1565: {
1566: AppCtx *user;
1567: DM dm;
1568: PetscReal t;
1570: PetscFunctionBeginUser;
1571: PetscCall(TSGetDM(ts, &dm));
1572: PetscCall(TSGetTime(ts, &t));
1573: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1574: PetscCall(DMGetApplicationContext(dm, &user));
1575: PetscCall(RemoveDiscretePressureNullspace_Private(ts, u));
1576: PetscFunctionReturn(PETSC_SUCCESS);
1577: }
1579: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
1580: {
1581: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
1582: void *ctxs[3];
1583: PetscPointFunc diagnostics[1] = {divergence};
1584: DM dm, dmCell = NULL;
1585: PetscDS ds;
1586: Vec v, divu;
1587: PetscReal ferrors[3], massFlux;
1588: PetscInt f;
1590: PetscFunctionBeginUser;
1591: PetscCall(TSGetDM(ts, &dm));
1592: PetscCall(DMGetDS(dm, &ds));
1594: for (f = 0; f < 3; ++f) PetscCall(PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]));
1595: PetscCall(DMComputeL2FieldDiff(dm, crtime, exactFuncs, ctxs, u, ferrors));
1596: PetscCall(GetCellDM(dm, (AppCtx *)ctx, &dmCell));
1597: PetscCall(DMGetGlobalVector(dmCell, &divu));
1598: PetscCall(DMProjectField(dmCell, crtime, u, diagnostics, INSERT_VALUES, divu));
1599: PetscCall(VecViewFromOptions(divu, NULL, "-divu_vec_view"));
1600: PetscCall(VecNorm(divu, NORM_2, &massFlux));
1601: 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));
1603: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
1605: PetscCall(DMGetGlobalVector(dm, &v));
1606: PetscCall(DMProjectFunction(dm, crtime, exactFuncs, ctxs, INSERT_ALL_VALUES, v));
1607: PetscCall(PetscObjectSetName((PetscObject)v, "Exact Solution"));
1608: PetscCall(VecViewFromOptions(v, NULL, "-exact_vec_view"));
1609: PetscCall(DMRestoreGlobalVector(dm, &v));
1611: PetscCall(VecViewFromOptions(divu, NULL, "-div_vec_view"));
1612: PetscCall(DMRestoreGlobalVector(dmCell, &divu));
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: int main(int argc, char **argv)
1617: {
1618: DM dm; /* problem definition */
1619: TS ts; /* timestepper */
1620: Vec u; /* solution */
1621: AppCtx user; /* user-defined work context */
1622: PetscReal t;
1624: PetscFunctionBeginUser;
1625: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1626: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
1627: PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
1628: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
1629: PetscCall(SetupParameters(dm, &user));
1630: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1631: PetscCall(TSSetDM(ts, dm));
1632: PetscCall(DMSetApplicationContext(dm, &user));
1633: /* Setup problem */
1634: PetscCall(SetupDiscretization(dm, &user));
1635: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
1637: PetscCall(DMCreateGlobalVector(dm, &u));
1638: PetscCall(PetscObjectSetName((PetscObject)u, "Numerical Solution"));
1639: if (user.hasNullSpace) PetscCall(DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace));
1641: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user));
1642: PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user));
1643: PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user));
1644: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
1645: PetscCall(TSSetPreStep(ts, RemoveDiscretePressureNullspace));
1646: PetscCall(TSSetFromOptions(ts));
1648: PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions)); /* Must come after SetFromOptions() */
1649: PetscCall(SetInitialConditions(ts, u));
1650: PetscCall(TSGetTime(ts, &t));
1651: PetscCall(DMSetOutputSequenceNumber(dm, 0, t));
1652: PetscCall(DMTSCheckFromOptions(ts, u));
1653: PetscCall(TSMonitorSet(ts, MonitorError, &user, NULL));
1655: PetscCall(TSSolve(ts, u));
1656: PetscCall(DMTSCheckFromOptions(ts, u));
1658: PetscCall(VecDestroy(&u));
1659: PetscCall(DMDestroy(&user.dmCell));
1660: PetscCall(DMDestroy(&dm));
1661: PetscCall(TSDestroy(&ts));
1662: PetscCall(PetscBagDestroy(&user.bag));
1663: PetscCall(PetscFinalize());
1664: return 0;
1665: }
1667: /*TEST
1669: testset:
1670: requires: triangle !single
1671: args: -dm_plex_separate_marker \
1672: -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \
1673: -snes_error_if_not_converged -snes_convergence_test correct_pressure \
1674: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
1675: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 \
1676: -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
1677: -fieldsplit_0_pc_type lu \
1678: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1680: test:
1681: suffix: 2d_tri_p2_p1_p1
1682: args: -sol_type quadratic \
1683: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1684: -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1
1686: test:
1687: # Using -dm_refine 5 -convest_num_refine 2 gives L_2 convergence rate: [0.89, 0.011, 1.0]
1688: suffix: 2d_tri_p2_p1_p1_tconv
1689: args: -sol_type cubic_trig \
1690: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1691: -ts_max_steps 4 -ts_dt 0.1 -ts_convergence_estimate -convest_num_refine 1
1693: test:
1694: # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.5, 1.9]
1695: suffix: 2d_tri_p2_p1_p1_sconv
1696: args: -sol_type cubic \
1697: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1698: -ts_max_steps 1 -ts_dt 1e-4 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
1700: test:
1701: suffix: 2d_tri_p3_p2_p2
1702: args: -sol_type cubic \
1703: -vel_petscspace_degree 3 -pres_petscspace_degree 2 -temp_petscspace_degree 2 \
1704: -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1
1706: test:
1707: # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.1, 3.1]
1708: suffix: 2d_tri_p2_p1_p1_tg_sconv
1709: args: -sol_type taylor_green \
1710: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1711: -ts_max_steps 1 -ts_dt 1e-8 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1
1713: test:
1714: # Using -dm_refine 3 -convest_num_refine 2 gives L_2 convergence rate: [1.2, 1.5, 1.2]
1715: suffix: 2d_tri_p2_p1_p1_tg_tconv
1716: args: -sol_type taylor_green \
1717: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1718: -ts_max_steps 4 -ts_dt 0.1 -ts_convergence_estimate -convest_num_refine 1
1720: testset:
1721: requires: triangle !single
1722: args: -dm_plex_separate_marker -mod_type conducting \
1723: -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \
1724: -snes_error_if_not_converged -snes_max_linear_solve_fail 5 \
1725: -ksp_type fgmres -ksp_max_it 2 -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 \
1726: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 \
1727: -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
1728: -fieldsplit_0_pc_type lu \
1729: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
1731: test:
1732: # At this resolution, the rhs is inconsistent on some Newton steps
1733: suffix: 2d_tri_p2_p1_p1_conduct
1734: args: -sol_type quadratic \
1735: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
1736: -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1 \
1737: -pc_fieldsplit_schur_precondition full \
1738: -fieldsplit_pressure_ksp_max_it 2 -fieldsplit_pressure_pc_type svd
1740: test:
1741: suffix: 2d_tri_p2_p1_p2_conduct_pipe
1742: args: -sol_type pipe \
1743: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 2 \
1744: -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1
1746: test:
1747: suffix: 2d_tri_p2_p1_p2_conduct_pipe_wiggly_sconv
1748: args: -sol_type pipe_wiggly -Fr 1e10 \
1749: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 2 \
1750: -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \
1751: -ts_max_steps 1 -ts_dt 1e10 \
1752: -ksp_atol 1e-12 -ksp_max_it 300 \
1753: -fieldsplit_pressure_ksp_atol 1e-14
1755: TEST*/