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