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