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