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 in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n"; 4 5 /*F 6 This 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 For visualization, use 18 19 -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append 20 F*/ 21 22 #include <petscdmplex.h> 23 #include <petscsnes.h> 24 #include <petscts.h> 25 #include <petscds.h> 26 #include <petscbag.h> 27 28 typedef enum {SOL_QUADRATIC, SOL_CUBIC, SOL_CUBIC_TRIG, SOL_TAYLOR_GREEN, NUM_SOL_TYPES} SolType; 29 const char *solTypes[NUM_SOL_TYPES+1] = {"quadratic", "cubic", "cubic_trig", "taylor_green", "unknown"}; 30 31 typedef struct { 32 PetscReal nu; /* Kinematic viscosity */ 33 PetscReal alpha; /* Thermal diffusivity */ 34 PetscReal T_in; /* Inlet temperature*/ 35 } Parameter; 36 37 typedef struct { 38 /* Problem definition */ 39 PetscBag bag; /* Holds problem parameters */ 40 SolType solType; /* MMS solution type */ 41 /* Flow diagnostics */ 42 DM dmCell; /* A DM with piecewise constant discretization */ 43 } AppCtx; 44 45 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 46 { 47 PetscInt d; 48 for (d = 0; d < Nc; ++d) u[d] = 0.0; 49 return 0; 50 } 51 52 static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 53 { 54 PetscInt d; 55 for (d = 0; d < Nc; ++d) u[d] = 1.0; 56 return 0; 57 } 58 59 /* 60 CASE: quadratic 61 In 2D we use exact solution: 62 63 u = t + x^2 + y^2 64 v = t + 2x^2 - 2xy 65 p = x + y - 1 66 T = t + x + y 67 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> 68 Q = 1 + 2t + 3x^2 - 2xy + y^2 69 70 so that 71 72 \nabla \cdot u = 2x - 2x = 0 73 74 f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p 75 = <1, 1> + <t + x^2 + y^2, t + 2x^2 - 2xy> . <<2x, 4x - 2y>, <2y, -2x>> - \nu <4, 4> + <1, 1> 76 = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t (2x - 2y) + 2x^2y + 4xy^2 - 2y^3> + <-4 \nu + 2, -4\nu + 2> 77 = <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2 - 4\nu + 2, t (2x - 2y) + 4xy^2 + 2x^2y - 2y^3 - 4\nu + 2> 78 79 Q = dT/dt + u \cdot \nabla T - \alpha \Delta T 80 = 1 + <t + x^2 + y^2, t + 2x^2 - 2xy> . <1, 1> - \alpha 0 81 = 1 + 2t + 3x^2 - 2xy + y^2 82 */ 83 84 static PetscErrorCode quadratic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 85 { 86 u[0] = time + X[0]*X[0] + X[1]*X[1]; 87 u[1] = time + 2.0*X[0]*X[0] - 2.0*X[0]*X[1]; 88 return 0; 89 } 90 static PetscErrorCode quadratic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 91 { 92 u[0] = 1.0; 93 u[1] = 1.0; 94 return 0; 95 } 96 97 static PetscErrorCode quadratic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 98 { 99 p[0] = X[0] + X[1] - 1.0; 100 return 0; 101 } 102 103 static PetscErrorCode quadratic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 104 { 105 T[0] = time + X[0] + X[1]; 106 return 0; 107 } 108 static PetscErrorCode quadratic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 109 { 110 T[0] = 1.0; 111 return 0; 112 } 113 114 /* f0_v = du/dt - f */ 115 static void f0_quadratic_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 116 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 117 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 118 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 119 { 120 const PetscReal nu = PetscRealPart(constants[0]); 121 PetscInt Nc = dim; 122 PetscInt c, d; 123 124 for (d = 0; d<dim; ++d) f0[d] = u_t[uOff[0]+d]; 125 126 for (c = 0; c < Nc; ++c) { 127 for (d = 0; d < dim; ++d) f0[c] += u[d]*u_x[c*dim+d]; 128 } 129 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); 130 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); 131 } 132 133 /* f0_w = dT/dt + u.grad(T) - Q */ 134 static void f0_quadratic_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 135 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 136 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 137 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 138 { 139 PetscInt d; 140 f0[0] = 0; 141 for (d = 0; d < dim; ++d) f0[0] += u[uOff[0]+d]*u_x[uOff_x[2]+d]; 142 f0[0] += u_t[uOff[2]] - (2*t + 1 + 3*X[0]*X[0] - 2*X[0]*X[1] + X[1]*X[1]); 143 } 144 145 /* 146 CASE: cubic 147 In 2D we use exact solution: 148 149 u = t + x^3 + y^3 150 v = t + 2x^3 - 3x^2y 151 p = 3/2 x^2 + 3/2 y^2 - 1 152 T = t + 1/2 x^2 + 1/2 y^2 153 f = < t(3x^2 + 3y^2) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x + 1, 154 t(3x^2 - 6xy) + 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y + 1> 155 Q = x^4 + xy^3 + 2x^3y - 3x^2y^2 + xt + yt - 2\alpha + 1 156 157 so that 158 159 \nabla \cdot u = 3x^2 - 3x^2 = 0 160 161 du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p - f 162 = <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 163 164 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 165 */ 166 static PetscErrorCode cubic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 167 { 168 u[0] = time + X[0]*X[0]*X[0] + X[1]*X[1]*X[1]; 169 u[1] = time + 2.0*X[0]*X[0]*X[0] - 3.0*X[0]*X[0]*X[1]; 170 return 0; 171 } 172 static PetscErrorCode cubic_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 173 { 174 u[0] = 1.0; 175 u[1] = 1.0; 176 return 0; 177 } 178 179 static PetscErrorCode cubic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 180 { 181 p[0] = 3.0*X[0]*X[0]/2.0 + 3.0*X[1]*X[1]/2.0 - 1.0; 182 return 0; 183 } 184 185 static PetscErrorCode cubic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 186 { 187 T[0] = time + X[0]*X[0]/2.0 + X[1]*X[1]/2.0; 188 return 0; 189 } 190 static PetscErrorCode cubic_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 191 { 192 T[0] = 1.0; 193 return 0; 194 } 195 196 197 static void f0_cubic_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 198 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 199 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 200 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 201 { 202 PetscInt c, d; 203 PetscInt Nc = dim; 204 const PetscReal nu = PetscRealPart(constants[0]); 205 206 for (d=0; d<dim; ++d) f0[d] = u_t[uOff[0]+d]; 207 208 for (c=0; c<Nc; ++c) { 209 for (d=0; d<dim; ++d) f0[c] += u[d]*u_x[c*dim+d]; 210 } 211 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); 212 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); 213 } 214 215 static void f0_cubic_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 216 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 217 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 218 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 219 { 220 PetscInt d; 221 const PetscReal alpha = PetscRealPart(constants[1]); 222 223 for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0]+d]*u_x[uOff_x[2]+d]; 224 f0[0] += u_t[uOff[2]] - (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); 225 } 226 227 /* 228 CASE: cubic-trigonometric 229 In 2D we use exact solution: 230 231 u = beta cos t + x^3 + y^3 232 v = beta sin t + 2x^3 - 3x^2y 233 p = 3/2 x^2 + 3/2 y^2 - 1 234 T = 20 cos t + 1/2 x^2 + 1/2 y^2 235 f = < beta cos t 3x^2 + beta sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y) + 3x, 236 beta cos t (6x^2 - 6xy) - beta sin t (3x^2) + 3x^4y + 6x^2y^3 - 6xy^4 - \nu(12x - 6y) + 3y> 237 Q = beta cos t x + beta sin t (y - 1) + x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2\alpha 238 239 so that 240 241 \nabla \cdot u = 3x^2 - 3x^2 = 0 242 243 f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p 244 = <-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> 245 = <-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> 246 = <cos t (3x^2) + sin t (3y^2 - 1) + 3x^5 + 6x^3y^2 - 6x^2y^3 - \nu (6x + 6y) + 3x, 247 cos t (6x^2 - 6xy) - sin t (3x^2) + 3x^4y + 6x^2y^3 - 6xy^4 - \nu (12x - 6y) + 3y> 248 249 Q = dT/dt + u \cdot \nabla T - \alpha \Delta T 250 = -sin t + <cos t + x^3 + y^3, sin t + 2x^3 - 3x^2y> . <x, y> - 2 \alpha 251 = -sin t + cos t (x) + x^4 + xy^3 + sin t (y) + 2x^3y - 3x^2y^2 - 2 \alpha 252 = cos t x + sin t (y - 1) + (x^4 + 2x^3y - 3x^2y^2 + xy^3 - 2 \alpha) 253 */ 254 static PetscErrorCode cubic_trig_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 255 { 256 u[0] = 100.*PetscCosReal(time) + X[0]*X[0]*X[0] + X[1]*X[1]*X[1]; 257 u[1] = 100.*PetscSinReal(time) + 2.0*X[0]*X[0]*X[0] - 3.0*X[0]*X[0]*X[1]; 258 return 0; 259 } 260 static PetscErrorCode cubic_trig_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 261 { 262 u[0] = -100.*PetscSinReal(time); 263 u[1] = 100.*PetscCosReal(time); 264 return 0; 265 } 266 267 static PetscErrorCode cubic_trig_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 268 { 269 p[0] = 3.0*X[0]*X[0]/2.0 + 3.0*X[1]*X[1]/2.0 - 1.0; 270 return 0; 271 } 272 273 static PetscErrorCode cubic_trig_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 274 { 275 T[0] = 100.*PetscCosReal(time) + X[0]*X[0]/2.0 + X[1]*X[1]/2.0; 276 return 0; 277 } 278 static PetscErrorCode cubic_trig_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 279 { 280 T[0] = -100.*PetscSinReal(time); 281 return 0; 282 } 283 284 static void f0_cubic_trig_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 285 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 286 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 287 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 288 { 289 const PetscReal nu = PetscRealPart(constants[0]); 290 PetscInt Nc = dim; 291 PetscInt c, d; 292 293 for (d = 0; d < dim; ++d) f0[d] = u_t[uOff[0]+d]; 294 295 for (c = 0; c < Nc; ++c) { 296 for (d = 0; d < dim; ++d) f0[c] += u[d]*u_x[c*dim+d]; 297 } 298 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]; 299 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]; 300 } 301 302 static void f0_cubic_trig_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 303 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 304 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 305 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 306 { 307 const PetscReal alpha = PetscRealPart(constants[1]); 308 PetscInt d; 309 310 for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0]+d]*u_x[uOff_x[2]+d]; 311 f0[0] += u_t[uOff[2]] - (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); 312 } 313 314 /* 315 CASE: taylor-green vortex 316 In 2D we use exact solution: 317 318 u = 1 - cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t) 319 v = 1 + sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t) 320 p = -1/4 [cos(2 \pi(x - t)) + cos(2 \pi(y - t))] exp(-4 \pi^2 \nu t) 321 T = t + x + y 322 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)) > 323 Q = 3 + sin(\pi(x-y)) exp(-2\nu \pi^2 t) 324 325 so that 326 327 \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 328 329 f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p 330 = <-\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), 331 \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)> 332 + < \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), 333 \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)> 334 + <-\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), 335 -\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)> 336 + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t), 337 2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)> 338 + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t), 339 \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)> 340 = <-\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), 341 \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)> 342 + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t), 343 \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)> 344 + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t), 345 -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)> 346 + <-\pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t), 347 -\pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)> 348 + <-2\pi^2 cos(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t), 349 2\pi^2 sin(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)> 350 + < \pi/2 sin(2\pi(x - t)) exp(-4 \pi^2 \nu t), 351 \pi/2 sin(2\pi(y - t)) exp(-4 \pi^2 \nu t)> 352 = <-\pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t), 353 \pi (sin(\pi(x - t)) sin(\pi(y - t)) - cos(\pi(x - t)) cos(\pi(y - t))) exp(-2 \pi^2 \nu t)> 354 + < \pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t), 355 \pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t)> 356 + <-\pi cos(\pi(x - t)) cos(\pi(y - t)) exp(-2 \pi^2 \nu t), 357 -\pi sin(\pi(x - t)) sin(\pi(y - t)) exp(-2 \pi^2 \nu t)> 358 = < \pi cos(\pi(x - t)) cos(\pi(y - t)), 359 \pi sin(\pi(x - t)) sin(\pi(y - t))> 360 + <-\pi cos(\pi(x - t)) cos(\pi(y - t)), 361 -\pi sin(\pi(x - t)) sin(\pi(y - t))> = 0 362 Q = dT/dt + u \cdot \nabla T - \alpha \Delta T 363 = 1 + u \cdot <1, 1> - 0 364 = 1 + u + v 365 */ 366 367 static PetscErrorCode taylor_green_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 368 { 369 u[0] = 1 - PetscCosReal(PETSC_PI*(X[0]-time))*PetscSinReal(PETSC_PI*(X[1]-time))*PetscExpReal(-2*PETSC_PI*PETSC_PI*time); 370 u[1] = 1 + PetscSinReal(PETSC_PI*(X[0]-time))*PetscCosReal(PETSC_PI*(X[1]-time))*PetscExpReal(-2*PETSC_PI*PETSC_PI*time); 371 return 0; 372 } 373 static PetscErrorCode taylor_green_u_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 374 { 375 u[0] = -PETSC_PI*(PetscSinReal(PETSC_PI*(X[0]-time))*PetscSinReal(PETSC_PI*(X[1]-time)) 376 - PetscCosReal(PETSC_PI*(X[0]-time))*PetscCosReal(PETSC_PI*(X[1]-time)) 377 - 2*PETSC_PI*PetscCosReal(PETSC_PI*(X[0]-time))*PetscSinReal(PETSC_PI*(X[1]-time)))*PetscExpReal(-2*PETSC_PI*PETSC_PI*time); 378 u[1] = PETSC_PI*(PetscSinReal(PETSC_PI*(X[0]-time))*PetscSinReal(PETSC_PI*(X[1]-time)) 379 - PetscCosReal(PETSC_PI*(X[0]-time))*PetscCosReal(PETSC_PI*(X[1]-time)) 380 - 2*PETSC_PI*PetscSinReal(PETSC_PI*(X[0]-time))*PetscCosReal(PETSC_PI*(X[1]-time)))*PetscExpReal(-2*PETSC_PI*PETSC_PI*time); 381 return 0; 382 } 383 384 static PetscErrorCode taylor_green_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 385 { 386 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); 387 return 0; 388 } 389 390 static PetscErrorCode taylor_green_p_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 391 { 392 p[0] = PETSC_PI*(0.5*(PetscSinReal(2*PETSC_PI*(X[0]-time)) + PetscSinReal(2*PETSC_PI*(X[1]-time))) 393 + PETSC_PI*(PetscCosReal(2*PETSC_PI*(X[0]-time)) + PetscCosReal(2*PETSC_PI*(X[1]-time))))*PetscExpReal(-4*PETSC_PI*PETSC_PI*time); 394 return 0; 395 } 396 397 static PetscErrorCode taylor_green_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 398 { 399 T[0] = time + X[0] + X[1]; 400 return 0; 401 } 402 static PetscErrorCode taylor_green_T_t(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 403 { 404 T[0] = 1.0; 405 return 0; 406 } 407 408 static void f0_taylor_green_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 409 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 410 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 411 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 412 { 413 PetscInt Nc = dim; 414 PetscInt c, d; 415 416 for (d = 0; d < dim; ++d) f0[d] = u_t[uOff[0]+d]; 417 418 for (c = 0; c < Nc; ++c) { 419 for (d = 0; d < dim; ++d) f0[c] += u[d]*u_x[c*dim+d]; 420 } 421 } 422 423 static void f0_taylor_green_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 424 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 425 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 426 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 427 { 428 PetscScalar vel[2]; 429 PetscInt d; 430 431 taylor_green_u(dim, t, X, Nf, vel, NULL); 432 for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0]+d]*u_x[uOff_x[2]+d]; 433 f0[0] += u_t[uOff[2]] - (1.0 + vel[0] + vel[1]); 434 } 435 436 static void f0_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, 437 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 438 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 439 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 440 { 441 PetscInt d; 442 for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d*dim+d]; 443 } 444 445 /*f1_v = \nu[grad(u) + grad(u)^T] - pI */ 446 static void f1_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 447 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 448 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 449 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 450 { 451 const PetscReal nu = PetscRealPart(constants[0]); 452 const PetscInt Nc = dim; 453 PetscInt c, d; 454 455 for (c = 0; c < Nc; ++c) { 456 for (d = 0; d < dim; ++d) { 457 f1[c*dim+d] = nu*(u_x[c*dim+d] + u_x[d*dim+c]); 458 } 459 f1[c*dim+c] -= u[uOff[1]]; 460 } 461 } 462 463 static void f1_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 464 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 465 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 466 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 467 { 468 const PetscReal alpha = PetscRealPart(constants[1]); 469 PetscInt d; 470 for (d = 0; d < dim; ++d) f1[d] = alpha*u_x[uOff_x[2]+d]; 471 } 472 473 /*Jacobians*/ 474 static void g1_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 475 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 476 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 477 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 478 { 479 PetscInt d; 480 for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; 481 } 482 483 static void g0_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 484 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 485 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 486 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 487 { 488 PetscInt c, d; 489 const PetscInt Nc = dim; 490 491 for (d = 0; d < dim; ++d) g0[d*dim+d] = u_tShift; 492 493 for (c = 0; c < Nc; ++c) { 494 for (d = 0; d < dim; ++d) { 495 g0[c*Nc+d] += u_x[ c*Nc+d]; 496 } 497 } 498 } 499 500 static void g1_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 501 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 502 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 503 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 504 { 505 PetscInt NcI = dim; 506 PetscInt NcJ = dim; 507 PetscInt c, d, e; 508 509 for (c = 0; c < NcI; ++c) { 510 for (d = 0; d < NcJ; ++d) { 511 for (e = 0; e < dim; ++e) { 512 if (c == d) { 513 g1[(c*NcJ+d)*dim+e] += u[e]; 514 } 515 } 516 } 517 } 518 } 519 520 521 static void g2_vp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 522 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 523 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 524 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 525 { 526 PetscInt d; 527 for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; 528 } 529 530 static void g3_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 531 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 532 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 533 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 534 { 535 const PetscReal nu = PetscRealPart(constants[0]); 536 const PetscInt Nc = dim; 537 PetscInt c, d; 538 539 for (c = 0; c < Nc; ++c) { 540 for (d = 0; d < dim; ++d) { 541 g3[((c*Nc+c)*dim+d)*dim+d] += nu; 542 g3[((c*Nc+d)*dim+d)*dim+c] += nu; 543 } 544 } 545 } 546 547 static void g0_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 548 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 549 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 550 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 551 { 552 g0[0] = u_tShift; 553 } 554 555 static void g0_wu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 556 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 557 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 558 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 559 { 560 PetscInt d; 561 for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2]+d]; 562 } 563 564 static void g1_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 565 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 566 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 567 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 568 { 569 PetscInt d; 570 for (d = 0; d < dim; ++d) g1[d] = u[uOff[0]+d]; 571 } 572 573 static void g3_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 574 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 575 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 576 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 577 { 578 const PetscReal alpha = PetscRealPart(constants[1]); 579 PetscInt d; 580 581 for (d = 0; d < dim; ++d) g3[d*dim+d] = alpha; 582 } 583 584 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 585 { 586 PetscInt sol; 587 PetscErrorCode ierr; 588 589 590 PetscFunctionBeginUser; 591 options->solType = SOL_QUADRATIC; 592 593 ierr = PetscOptionsBegin(comm, "", "Low Mach flow Problem Options", "DMPLEX");CHKERRQ(ierr); 594 sol = options->solType; 595 ierr = PetscOptionsEList("-sol_type", "The solution type", "ex62.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL);CHKERRQ(ierr); 596 options->solType = (SolType) sol; 597 ierr = PetscOptionsEnd(); 598 PetscFunctionReturn(0); 599 } 600 601 static PetscErrorCode SetupParameters(AppCtx *user) 602 { 603 PetscBag bag; 604 Parameter *p; 605 PetscErrorCode ierr; 606 607 PetscFunctionBeginUser; 608 /* setup PETSc parameter bag */ 609 ierr = PetscBagGetData(user->bag, (void **) &p);CHKERRQ(ierr); 610 ierr = PetscBagSetName(user->bag, "par", "Low Mach flow parameters");CHKERRQ(ierr); 611 bag = user->bag; 612 ierr = PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity");CHKERRQ(ierr); 613 ierr = PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity");CHKERRQ(ierr); 614 ierr = PetscBagRegisterReal(bag, &p->T_in, 1.0, "T_in", "Inlet temperature");CHKERRQ(ierr); 615 PetscFunctionReturn(0); 616 } 617 618 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 619 { 620 PetscErrorCode ierr; 621 622 PetscFunctionBeginUser; 623 ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 624 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 625 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 626 PetscFunctionReturn(0); 627 } 628 629 static PetscErrorCode SetupProblem(DM dm, AppCtx *user) 630 { 631 PetscSimplePointFunc exactFuncs[3]; 632 PetscSimplePointFunc exactFuncs_t[3]; 633 PetscDS prob; 634 DMLabel label; 635 Parameter *ctx; 636 PetscInt id; 637 PetscErrorCode ierr; 638 639 PetscFunctionBeginUser; 640 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 641 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 642 switch(user->solType){ 643 case SOL_QUADRATIC: 644 ierr = PetscDSSetResidual(prob, 0, f0_quadratic_v, f1_v);CHKERRQ(ierr); 645 ierr = PetscDSSetResidual(prob, 2, f0_quadratic_w, f1_w);CHKERRQ(ierr); 646 647 exactFuncs[0] = quadratic_u; 648 exactFuncs[1] = quadratic_p; 649 exactFuncs[2] = quadratic_T; 650 exactFuncs_t[0] = quadratic_u_t; 651 exactFuncs_t[1] = NULL; 652 exactFuncs_t[2] = quadratic_T_t; 653 break; 654 case SOL_CUBIC: 655 ierr = PetscDSSetResidual(prob, 0, f0_cubic_v, f1_v);CHKERRQ(ierr); 656 ierr = PetscDSSetResidual(prob, 2, f0_cubic_w, f1_w);CHKERRQ(ierr); 657 658 exactFuncs[0] = cubic_u; 659 exactFuncs[1] = cubic_p; 660 exactFuncs[2] = cubic_T; 661 exactFuncs_t[0] = cubic_u_t; 662 exactFuncs_t[1] = NULL; 663 exactFuncs_t[2] = cubic_T_t; 664 break; 665 case SOL_CUBIC_TRIG: 666 ierr = PetscDSSetResidual(prob, 0, f0_cubic_trig_v, f1_v);CHKERRQ(ierr); 667 ierr = PetscDSSetResidual(prob, 2, f0_cubic_trig_w, f1_w);CHKERRQ(ierr); 668 669 exactFuncs[0] = cubic_trig_u; 670 exactFuncs[1] = cubic_trig_p; 671 exactFuncs[2] = cubic_trig_T; 672 exactFuncs_t[0] = cubic_trig_u_t; 673 exactFuncs_t[1] = NULL; 674 exactFuncs_t[2] = cubic_trig_T_t; 675 break; 676 case SOL_TAYLOR_GREEN: 677 ierr = PetscDSSetResidual(prob, 0, f0_taylor_green_v, f1_v);CHKERRQ(ierr); 678 ierr = PetscDSSetResidual(prob, 2, f0_taylor_green_w, f1_w);CHKERRQ(ierr); 679 680 exactFuncs[0] = taylor_green_u; 681 exactFuncs[1] = taylor_green_p; 682 exactFuncs[2] = taylor_green_T; 683 exactFuncs_t[0] = taylor_green_u_t; 684 exactFuncs_t[1] = taylor_green_p_t; 685 exactFuncs_t[2] = taylor_green_T_t; 686 break; 687 default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%D)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType); 688 } 689 690 ierr = PetscDSSetResidual(prob, 1, f0_q, NULL);CHKERRQ(ierr); 691 692 ierr = PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu);CHKERRQ(ierr); 693 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL);CHKERRQ(ierr); 694 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL);CHKERRQ(ierr); 695 ierr = PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL);CHKERRQ(ierr); 696 ierr = PetscDSSetJacobian(prob, 2, 2, g0_wT, g1_wT, NULL, g3_wT);CHKERRQ(ierr); 697 /* Setup constants */ 698 { 699 Parameter *param; 700 PetscScalar constants[3]; 701 702 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 703 704 constants[0] = param->nu; 705 constants[1] = param->alpha; 706 constants[2] = param->T_in; 707 ierr = PetscDSSetConstants(prob, 3, constants);CHKERRQ(ierr); 708 } 709 /* Setup Boundary Conditions */ 710 ierr = PetscBagGetData(user->bag, (void **) &ctx);CHKERRQ(ierr); 711 id = 3; 712 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], ctx, NULL);CHKERRQ(ierr); 713 id = 1; 714 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], ctx, NULL);CHKERRQ(ierr); 715 id = 2; 716 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], ctx, NULL);CHKERRQ(ierr); 717 id = 4; 718 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], ctx, NULL);CHKERRQ(ierr); 719 id = 3; 720 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], ctx, NULL);CHKERRQ(ierr); 721 id = 1; 722 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], ctx, NULL);CHKERRQ(ierr); 723 id = 2; 724 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], ctx, NULL);CHKERRQ(ierr); 725 id = 4; 726 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], ctx, NULL);CHKERRQ(ierr); 727 728 /*setup exact solution.*/ 729 ierr = PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx);CHKERRQ(ierr); 730 ierr = PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx);CHKERRQ(ierr); 731 ierr = PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx);CHKERRQ(ierr); 732 ierr = PetscDSSetExactSolutionTimeDerivative(prob, 0, exactFuncs_t[0], ctx);CHKERRQ(ierr); 733 ierr = PetscDSSetExactSolutionTimeDerivative(prob, 1, exactFuncs_t[1], ctx);CHKERRQ(ierr); 734 ierr = PetscDSSetExactSolutionTimeDerivative(prob, 2, exactFuncs_t[2], ctx);CHKERRQ(ierr); 735 PetscFunctionReturn(0); 736 } 737 738 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 739 { 740 DM cdm = dm; 741 PetscFE fe[3], fediv; 742 Parameter *param; 743 DMPolytopeType ct; 744 PetscInt dim, cStart; 745 PetscBool simplex; 746 PetscErrorCode ierr; 747 748 PetscFunctionBeginUser; 749 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 750 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 751 ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr); 752 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE; 753 /* Create finite element */ 754 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);CHKERRQ(ierr); 755 ierr = PetscObjectSetName((PetscObject) fe[0], "velocity");CHKERRQ(ierr); 756 757 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);CHKERRQ(ierr); 758 ierr = PetscFECopyQuadrature(fe[0], fe[1]);CHKERRQ(ierr); 759 ierr = PetscObjectSetName((PetscObject) fe[1], "pressure");CHKERRQ(ierr); 760 761 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]);CHKERRQ(ierr); 762 ierr = PetscFECopyQuadrature(fe[0], fe[2]);CHKERRQ(ierr); 763 ierr = PetscObjectSetName((PetscObject) fe[2], "temperature");CHKERRQ(ierr); 764 765 ierr = PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "div_", PETSC_DEFAULT, &fediv);CHKERRQ(ierr); 766 ierr = PetscFECopyQuadrature(fe[0], fediv);CHKERRQ(ierr); 767 ierr = PetscObjectSetName((PetscObject) fediv, "divergence");CHKERRQ(ierr); 768 769 /* Set discretization and boundary conditions for each mesh */ 770 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe[0]);CHKERRQ(ierr); 771 ierr = DMSetField(dm, 1, NULL, (PetscObject) fe[1]);CHKERRQ(ierr); 772 ierr = DMSetField(dm, 2, NULL, (PetscObject) fe[2]);CHKERRQ(ierr); 773 ierr = DMCreateDS(dm);CHKERRQ(ierr); 774 ierr = SetupProblem(dm, user);CHKERRQ(ierr); 775 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 776 while (cdm) { 777 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 778 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 779 } 780 ierr = PetscFEDestroy(&fe[0]);CHKERRQ(ierr); 781 ierr = PetscFEDestroy(&fe[1]);CHKERRQ(ierr); 782 ierr = PetscFEDestroy(&fe[2]);CHKERRQ(ierr); 783 784 ierr = DMClone(dm, &user->dmCell);CHKERRQ(ierr); 785 ierr = DMSetField(user->dmCell, 0, NULL, (PetscObject) fediv);CHKERRQ(ierr); 786 ierr = DMCreateDS(user->dmCell);CHKERRQ(ierr); 787 ierr = PetscFEDestroy(&fediv);CHKERRQ(ierr); 788 789 { 790 PetscObject pressure; 791 MatNullSpace nullspacePres; 792 793 ierr = DMGetField(dm, 1, NULL, &pressure);CHKERRQ(ierr); 794 ierr = MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres);CHKERRQ(ierr); 795 ierr = PetscObjectCompose(pressure, "nullspace", (PetscObject) nullspacePres);CHKERRQ(ierr); 796 ierr = MatNullSpaceDestroy(&nullspacePres);CHKERRQ(ierr); 797 } 798 799 PetscFunctionReturn(0); 800 } 801 802 static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace) 803 { 804 Vec vec; 805 PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero}; 806 PetscErrorCode ierr; 807 808 PetscFunctionBeginUser; 809 if (ofield != 1) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index 1, not %D", ofield); 810 funcs[nfield] = constant; 811 ierr = DMCreateGlobalVector(dm, &vec);CHKERRQ(ierr); 812 ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);CHKERRQ(ierr); 813 ierr = VecNormalize(vec, NULL);CHKERRQ(ierr); 814 ierr = PetscObjectSetName((PetscObject) vec, "Pressure Null Space");CHKERRQ(ierr); 815 ierr = VecViewFromOptions(vec, NULL, "-pressure_nullspace_view");CHKERRQ(ierr); 816 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_FALSE, 1, &vec, nullSpace);CHKERRQ(ierr); 817 ierr = VecDestroy(&vec);CHKERRQ(ierr); 818 PetscFunctionReturn(0); 819 } 820 821 static PetscErrorCode RemoveDiscretePressureNullspace_Private(TS ts, Vec u) 822 { 823 DM dm; 824 MatNullSpace nullsp; 825 PetscErrorCode ierr; 826 827 PetscFunctionBegin; 828 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 829 ierr = CreatePressureNullSpace(dm, 1, 1, &nullsp);CHKERRQ(ierr); 830 ierr = MatNullSpaceRemove(nullsp, u);CHKERRQ(ierr); 831 ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr); 832 PetscFunctionReturn(0); 833 } 834 835 /* Make the discrete pressure discretely divergence free */ 836 static PetscErrorCode RemoveDiscretePressureNullspace(TS ts) 837 { 838 Vec u; 839 PetscErrorCode ierr; 840 841 PetscFunctionBegin; 842 ierr = TSGetSolution(ts, &u);CHKERRQ(ierr); 843 ierr = RemoveDiscretePressureNullspace_Private(ts, u);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847 static PetscErrorCode SetInitialConditions(TS ts, Vec u) 848 { 849 DM dm; 850 PetscReal t; 851 PetscErrorCode ierr; 852 853 PetscFunctionBegin; 854 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 855 ierr = TSGetTime(ts, &t);CHKERRQ(ierr); 856 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr); 857 ierr = RemoveDiscretePressureNullspace_Private(ts, u);CHKERRQ(ierr); 858 PetscFunctionReturn(0); 859 } 860 861 static void divergence(PetscInt dim, PetscInt Nf, PetscInt NfAux, 862 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 863 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 864 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar divu[]) 865 { 866 PetscInt d; 867 868 divu[0] = 0.; 869 for (d = 0; d < dim; ++d) divu[0] += u_x[d*dim+d]; 870 } 871 872 static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) 873 { 874 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 875 void *ctxs[3]; 876 PetscPointFunc diagnostics[1] = {divergence}; 877 DM dm, dmCell = ((AppCtx *) ctx)->dmCell; 878 PetscDS ds; 879 Vec v, divu; 880 PetscReal ferrors[3], massFlux; 881 PetscInt f; 882 PetscErrorCode ierr; 883 884 PetscFunctionBeginUser; 885 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 886 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 887 888 for (f = 0; f < 3; ++f) {ierr = PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]);CHKERRQ(ierr);} 889 ierr = DMComputeL2FieldDiff(dm, crtime, exactFuncs, ctxs, u, ferrors);CHKERRQ(ierr); 890 ierr = DMGetGlobalVector(dmCell, &divu);CHKERRQ(ierr); 891 ierr = DMProjectField(dmCell, crtime, u, diagnostics, INSERT_VALUES, divu);CHKERRQ(ierr); 892 ierr = VecViewFromOptions(divu, NULL, "-divu_vec_view");CHKERRQ(ierr); 893 ierr = VecNorm(divu, NORM_2, &massFlux);CHKERRQ(ierr); 894 ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g, %2.3g] ||div u||: %2.3g\n", (int) step, (double) crtime, (double) ferrors[0], (double) ferrors[1], (double) ferrors[2], (double) massFlux);CHKERRQ(ierr); 895 896 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr); 897 898 ierr = DMGetGlobalVector(dm, &v);CHKERRQ(ierr); 899 ierr = DMProjectFunction(dm, crtime, exactFuncs, ctxs, INSERT_ALL_VALUES, v);CHKERRQ(ierr); 900 ierr = PetscObjectSetName((PetscObject) v, "Exact Solution");CHKERRQ(ierr); 901 ierr = VecViewFromOptions(v, NULL, "-exact_vec_view");CHKERRQ(ierr); 902 ierr = DMRestoreGlobalVector(dm, &v);CHKERRQ(ierr); 903 904 ierr = VecViewFromOptions(divu, NULL, "-div_vec_view");CHKERRQ(ierr); 905 ierr = DMRestoreGlobalVector(dmCell, &divu);CHKERRQ(ierr); 906 907 PetscFunctionReturn(0); 908 } 909 910 int main(int argc, char **argv) 911 { 912 DM dm; /* problem definition */ 913 TS ts; /* timestepper */ 914 Vec u; /* solution */ 915 AppCtx user; /* user-defined work context */ 916 PetscReal t; 917 PetscErrorCode ierr; 918 919 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 920 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 921 ierr = PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);CHKERRQ(ierr); 922 ierr = SetupParameters(&user);CHKERRQ(ierr); 923 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr); 924 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 925 ierr = TSSetDM(ts, dm);CHKERRQ(ierr); 926 ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr); 927 /* Setup problem */ 928 ierr = SetupDiscretization(dm, &user);CHKERRQ(ierr); 929 ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr); 930 931 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 932 ierr = PetscObjectSetName((PetscObject) u, "Numerical Solution");CHKERRQ(ierr); 933 ierr = DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace);CHKERRQ(ierr); 934 935 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user);CHKERRQ(ierr); 936 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user);CHKERRQ(ierr); 937 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user);CHKERRQ(ierr); 938 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 939 ierr = TSSetPreStep(ts, RemoveDiscretePressureNullspace);CHKERRQ(ierr); 940 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 941 942 ierr = TSSetComputeInitialCondition(ts, SetInitialConditions);CHKERRQ(ierr); /* Must come after SetFromOptions() */ 943 ierr = SetInitialConditions(ts, u);CHKERRQ(ierr); 944 ierr = TSGetTime(ts, &t);CHKERRQ(ierr); 945 ierr = DMSetOutputSequenceNumber(dm, 0, t);CHKERRQ(ierr); 946 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr); 947 ierr = TSMonitorSet(ts, MonitorError, &user, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 948 949 ierr = TSSolve(ts, u);CHKERRQ(ierr); 950 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr); 951 952 ierr = VecDestroy(&u);CHKERRQ(ierr); 953 ierr = DMDestroy(&user.dmCell);CHKERRQ(ierr); 954 ierr = DMDestroy(&dm);CHKERRQ(ierr); 955 ierr = TSDestroy(&ts);CHKERRQ(ierr); 956 ierr = PetscBagDestroy(&user.bag);CHKERRQ(ierr); 957 ierr = PetscFinalize(); 958 return ierr; 959 } 960 961 /*TEST 962 963 test: 964 suffix: 2d_tri_p2_p1_p1 965 requires: triangle !single 966 args: -dm_plex_separate_marker -sol_type quadratic -dm_refine 0 \ 967 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 968 -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \ 969 -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1 \ 970 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 971 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 972 -fieldsplit_0_pc_type lu \ 973 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 974 975 # TODO Need trig t for convergence in time, also need to refine in space 976 test: 977 # Using -dm_refine 5 -convest_num_refine 2 gives L_2 convergence rate: [0.89, 0.011, 1.0] 978 suffix: 2d_tri_p2_p1_p1_tconv 979 requires: triangle !single 980 args: -dm_plex_separate_marker -sol_type cubic_trig -dm_refine 0 \ 981 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 982 -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \ 983 -ts_max_steps 4 -ts_dt 0.1 -ts_convergence_estimate -convest_num_refine 1 \ 984 -snes_error_if_not_converged -snes_convergence_test correct_pressure \ 985 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 986 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 987 -fieldsplit_0_pc_type lu \ 988 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 989 990 test: 991 # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.5, 1.9] 992 suffix: 2d_tri_p2_p1_p1_sconv 993 requires: triangle !single 994 args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \ 995 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 996 -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \ 997 -ts_max_steps 1 -ts_dt 1e-4 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 998 -snes_error_if_not_converged -snes_convergence_test correct_pressure \ 999 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_atol 1e-16 -ksp_error_if_not_converged \ 1000 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 1001 -fieldsplit_0_pc_type lu \ 1002 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 1003 1004 test: 1005 suffix: 2d_tri_p3_p2_p2 1006 requires: triangle !single 1007 args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \ 1008 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -temp_petscspace_degree 2 \ 1009 -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \ 1010 -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1 \ 1011 -snes_convergence_test correct_pressure \ 1012 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 1013 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 1014 -fieldsplit_0_pc_type lu \ 1015 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 1016 1017 test: 1018 # Using -dm_refine 3 -convest_num_refine 3 gives L_2 convergence rate: [3.0, 2.1, 3.1] 1019 suffix: 2d_tri_p2_p1_p1_tg_sconv 1020 requires: triangle !single 1021 args: -dm_plex_separate_marker -sol_type taylor_green -dm_refine 0 \ 1022 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 1023 -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \ 1024 -ts_max_steps 1 -ts_dt 1e-8 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 1025 -snes_error_if_not_converged -snes_convergence_test correct_pressure \ 1026 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_atol 1e-16 -ksp_error_if_not_converged \ 1027 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 1028 -fieldsplit_0_pc_type lu \ 1029 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 1030 1031 test: 1032 # Using -dm_refine 3 -convest_num_refine 2 gives L_2 convergence rate: [1.2, 1.5, 1.2] 1033 suffix: 2d_tri_p2_p1_p1_tg_tconv 1034 requires: triangle !single 1035 args: -dm_plex_separate_marker -sol_type taylor_green -dm_refine 0 \ 1036 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 1037 -div_petscdualspace_lagrange_use_moments -div_petscdualspace_lagrange_moment_order 3 \ 1038 -ts_max_steps 4 -ts_dt 0.1 -ts_convergence_estimate -convest_num_refine 1 \ 1039 -snes_error_if_not_converged -snes_convergence_test correct_pressure \ 1040 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_atol 1e-16 -ksp_error_if_not_converged \ 1041 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 1042 -fieldsplit_0_pc_type lu \ 1043 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 1044 1045 TEST*/ 1046