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