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