1 static char help[] = "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 a steady-state 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, 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 <petscds.h> 25 #include <petscbag.h> 26 27 typedef enum {SOL_QUADRATIC, SOL_CUBIC, NUM_SOL_TYPES} SolType; 28 const char *solTypes[NUM_SOL_TYPES+1] = {"quadratic", "cubic", "unknown"}; 29 30 typedef struct { 31 PetscReal nu; /* Kinematic viscosity */ 32 PetscReal theta; /* Angle of pipe wall to x-axis */ 33 PetscReal alpha; /* Thermal diffusivity */ 34 PetscReal T_in; /* Inlet temperature*/ 35 } Parameter; 36 37 typedef struct { 38 PetscBool showError; 39 PetscBag bag; 40 SolType solType; 41 } AppCtx; 42 43 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 44 { 45 PetscInt d; 46 for (d = 0; d < Nc; ++d) u[d] = 0.0; 47 return 0; 48 } 49 50 static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 51 { 52 PetscInt d; 53 for (d = 0; d < Nc; ++d) u[d] = 1.0; 54 return 0; 55 } 56 57 /* 58 CASE: quadratic 59 In 2D we use exact solution: 60 61 u = x^2 + y^2 62 v = 2x^2 - 2xy 63 p = x + y - 1 64 T = x + y 65 f = <2x^3 + 4x^2y - 2xy^2 -4\nu + 1, 4xy^2 + 2x^2y - 2y^3 -4\nu + 1> 66 Q = 3x^2 + y^2 - 2xy 67 68 so that 69 70 (1) \nabla \cdot u = 2x - 2x = 0 71 72 (2) u \cdot \nabla u - \nu \Delta u + \nabla p - f 73 = <2x^3 + 4x^2y -2xy^2, 4xy^2 + 2x^2y - 2y^3> -\nu <4, 4> + <1, 1> - <2x^3 + 4x^2y - 2xy^2 -4\nu + 1, 4xy^2 + 2x^2y - 2y^3 - 4\nu + 1> = 0 74 75 (3) u \cdot \nabla T - \alpha \Delta T - Q = 3x^2 + y^2 - 2xy - \alpha*0 - 3x^2 - y^2 + 2xy = 0 76 */ 77 78 static PetscErrorCode quadratic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 79 { 80 u[0] = X[0]*X[0] + X[1]*X[1]; 81 u[1] = 2.0*X[0]*X[0] - 2.0*X[0]*X[1]; 82 return 0; 83 } 84 85 static PetscErrorCode linear_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 86 { 87 p[0] = X[0] + X[1] - 1.0; 88 return 0; 89 } 90 91 static PetscErrorCode linear_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 92 { 93 T[0] = X[0] + X[1]; 94 return 0; 95 } 96 97 static void f0_quadratic_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 98 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 99 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 100 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 101 { 102 PetscInt c, d; 103 PetscInt Nc = dim; 104 const PetscReal nu = PetscRealPart(constants[0]); 105 106 for (c=0; c<Nc; ++c) { 107 for (d=0; d<dim; ++d) f0[c] += u[d]*u_x[c*dim+d]; 108 } 109 f0[0] -= (2*X[0]*X[0]*X[0] + 4*X[0]*X[0]*X[1] - 2*X[0]*X[1]*X[1] - 4.0*nu + 1); 110 f0[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 + 1); 111 } 112 113 static void f0_quadratic_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 114 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 115 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 116 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 117 { 118 PetscInt d; 119 for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0]+d]*u_x[uOff_x[2]+d]; 120 f0[0] -= (3*X[0]*X[0] + X[1]*X[1] - 2*X[0]*X[1]); 121 } 122 123 /* 124 CASE: cubic 125 In 2D we use exact solution: 126 127 u = x^3 + y^3 128 v = 2x^3 - 3x^2y 129 p = 3/2 x^2 + 3/2 y^2 - 1 130 T = 1/2 x^2 + 1/2 y^2 131 f = <3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y), 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y> 132 Q = x^4 + xy^3 + 2x^3y - 3x^2y^2 - 2 133 134 so that 135 136 \nabla \cdot u = 3x^2 - 3x^2 = 0 137 138 u \cdot \nabla u - \nu \Delta u + \nabla p - f 139 = <3x^5 + 6x^3y^2 - 6x^2y^3, 6x^2y^3 + 3x^4y - 6xy^4> - \nu<6x + 6y, 12x - 6y> + <3x, 3y> - <3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y), 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y> = 0 140 141 u \cdot \nabla T - \alpha\Delta T - Q = (x^3 + y^3) x + (2x^3 - 3x^2y) y - 2*\alpha - (x^4 + xy^3 + 2x^3y - 3x^2y^2 - 2) = 0 142 */ 143 144 static PetscErrorCode cubic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 145 { 146 u[0] = X[0]*X[0]*X[0] + X[1]*X[1]*X[1]; 147 u[1] = 2.0*X[0]*X[0]*X[0] - 3.0*X[0]*X[0]*X[1]; 148 return 0; 149 } 150 151 static PetscErrorCode quadratic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 152 { 153 p[0] = 3.0*X[0]*X[0]/2.0 + 3.0*X[1]*X[1]/2.0 - 1.0; 154 return 0; 155 } 156 157 static PetscErrorCode quadratic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 158 { 159 T[0] = X[0]*X[0]/2.0 + X[1]*X[1]/2.0; 160 return 0; 161 } 162 163 static void f0_cubic_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 164 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 165 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 166 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 167 { 168 PetscInt c, d; 169 PetscInt Nc = dim; 170 const PetscReal nu = PetscRealPart(constants[0]); 171 172 for (c=0; c<Nc; ++c) { 173 for (d=0; d<dim; ++d) f0[c] += u[d]*u_x[c*dim+d]; 174 } 175 f0[0] -= (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]); 176 f0[1] -= (6*X[0]*X[0]*X[1]*X[1]*X[1] + 3*X[0]*X[0]*X[0]*X[0]*X[1] - 6*X[0]*X[1]*X[1]*X[1]*X[1] - (12*X[0] - 6*X[1])*nu + 3*X[1]); 177 } 178 179 static void f0_cubic_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 180 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 181 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 182 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 183 { 184 const PetscReal alpha = PetscRealPart(constants[1]); 185 PetscInt d; 186 187 for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0]+d]*u_x[uOff_x[2]+d]; 188 f0[0] -= (X[0]*X[0]*X[0]*X[0] + X[0]*X[1]*X[1]*X[1] + 2.0*X[0]*X[0]*X[0]*X[1] - 3.0*X[0]*X[0]*X[1]*X[1] - 2.0*alpha); 189 } 190 191 static void f0_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, 192 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 193 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 194 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 195 { 196 PetscInt d; 197 for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d*dim+d]; 198 } 199 200 static void f1_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 201 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 202 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 203 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 204 { 205 const PetscReal nu = PetscRealPart(constants[0]); 206 const PetscInt Nc = dim; 207 PetscInt c, d; 208 209 for (c = 0; c < Nc; ++c) { 210 for (d = 0; d < dim; ++d) { 211 f1[c*dim+d] = nu*(u_x[c*dim+d] + u_x[d*dim+c]); 212 //f1[c*dim+d] = nu*u_x[c*dim+d]; 213 } 214 f1[c*dim+c] -= u[uOff[1]]; 215 } 216 } 217 218 static void f1_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 219 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 220 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 221 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 222 { 223 const PetscReal alpha = PetscRealPart(constants[1]); 224 PetscInt d; 225 for (d = 0; d < dim; ++d) f1[d] = alpha*u_x[uOff_x[2]+d]; 226 } 227 228 /*Jacobians*/ 229 static void g1_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 230 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 231 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 232 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 233 { 234 PetscInt d; 235 for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; 236 } 237 238 static void g0_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 239 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 240 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 241 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 242 { 243 const PetscInt Nc = dim; 244 PetscInt c, d; 245 246 for (c = 0; c < Nc; ++c) { 247 for (d = 0; d < dim; ++d) { 248 g0[c*Nc+d] = u_x[ c*Nc+d]; 249 } 250 } 251 } 252 253 static void g1_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 254 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 255 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 256 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 257 { 258 PetscInt NcI = dim; 259 PetscInt NcJ = dim; 260 PetscInt c, d, e; 261 262 for (c = 0; c < NcI; ++c) { 263 for (d = 0; d < NcJ; ++d) { 264 for (e = 0; e < dim; ++e) { 265 if (c == d) { 266 g1[(c*NcJ+d)*dim+e] = u[e]; 267 } 268 } 269 } 270 } 271 } 272 273 static void g2_vp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 274 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 275 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 276 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 277 { 278 PetscInt d; 279 for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; 280 } 281 282 static void g3_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 283 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 284 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 285 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 286 { 287 const PetscReal nu = PetscRealPart(constants[0]); 288 const PetscInt Nc = dim; 289 PetscInt c, d; 290 291 for (c = 0; c < Nc; ++c) { 292 for (d = 0; d < dim; ++d) { 293 g3[((c*Nc+c)*dim+d)*dim+d] += nu; // gradU 294 g3[((c*Nc+d)*dim+d)*dim+c] += nu; // gradU transpose 295 } 296 } 297 } 298 299 static void g0_wu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 300 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 301 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 302 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 303 { 304 PetscInt d; 305 for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2]+d]; 306 } 307 308 static void g1_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 309 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 310 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 311 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 312 { 313 PetscInt d; 314 for (d = 0; d < dim; ++d) g1[d] = u[uOff[0]+d]; 315 } 316 317 static void g3_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 318 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 319 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 320 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 321 { 322 const PetscReal alpha = PetscRealPart(constants[1]); 323 PetscInt d; 324 325 for (d = 0; d < dim; ++d) g3[d*dim+d] = alpha; 326 } 327 328 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 329 { 330 PetscInt sol; 331 PetscErrorCode ierr; 332 333 PetscFunctionBeginUser; 334 options->solType = SOL_QUADRATIC; 335 options->showError = PETSC_FALSE; 336 337 ierr = PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");CHKERRQ(ierr); 338 sol = options->solType; 339 ierr = PetscOptionsEList("-sol_type", "The solution type", "ex62.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL);CHKERRQ(ierr); 340 options->solType = (SolType) sol; 341 ierr = PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL);CHKERRQ(ierr); 342 ierr = PetscOptionsEnd();CHKERRQ(ierr); 343 PetscFunctionReturn(0); 344 } 345 346 static PetscErrorCode SetupParameters(AppCtx *user) 347 { 348 PetscBag bag; 349 Parameter *p; 350 PetscErrorCode ierr; 351 352 PetscFunctionBeginUser; 353 /* setup PETSc parameter bag */ 354 ierr = PetscBagGetData(user->bag, (void **) &p);CHKERRQ(ierr); 355 ierr = PetscBagSetName(user->bag, "par", "Poiseuille flow parameters");CHKERRQ(ierr); 356 bag = user->bag; 357 ierr = PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity");CHKERRQ(ierr); 358 ierr = PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity");CHKERRQ(ierr); 359 ierr = PetscBagRegisterReal(bag, &p->theta, 0.0, "theta", "Angle of pipe wall to x-axis");CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 364 { 365 PetscErrorCode ierr; 366 367 PetscFunctionBeginUser; 368 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 369 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 370 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 371 { 372 Parameter *param; 373 Vec coordinates; 374 PetscScalar *coords; 375 PetscReal theta; 376 PetscInt cdim, N, bs, i; 377 378 ierr = DMGetCoordinateDim(*dm, &cdim);CHKERRQ(ierr); 379 ierr = DMGetCoordinates(*dm, &coordinates);CHKERRQ(ierr); 380 ierr = VecGetLocalSize(coordinates, &N);CHKERRQ(ierr); 381 ierr = VecGetBlockSize(coordinates, &bs);CHKERRQ(ierr); 382 PetscCheckFalse(bs != cdim,comm, PETSC_ERR_ARG_WRONG, "Invalid coordinate blocksize %D != embedding dimension %D", bs, cdim); 383 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 384 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 385 theta = param->theta; 386 for (i = 0; i < N; i += cdim) { 387 PetscScalar x = coords[i+0]; 388 PetscScalar y = coords[i+1]; 389 390 coords[i+0] = PetscCosReal(theta)*x - PetscSinReal(theta)*y; 391 coords[i+1] = PetscSinReal(theta)*x + PetscCosReal(theta)*y; 392 } 393 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 394 ierr = DMSetCoordinates(*dm, coordinates);CHKERRQ(ierr); 395 } 396 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 397 PetscFunctionReturn(0); 398 } 399 400 static PetscErrorCode SetupProblem(DM dm, AppCtx *user) 401 { 402 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 403 PetscDS prob; 404 DMLabel label; 405 Parameter *ctx; 406 PetscInt id; 407 PetscErrorCode ierr; 408 409 PetscFunctionBeginUser; 410 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 411 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 412 switch(user->solType) { 413 case SOL_QUADRATIC: 414 ierr = PetscDSSetResidual(prob, 0, f0_quadratic_v, f1_v);CHKERRQ(ierr); 415 ierr = PetscDSSetResidual(prob, 2, f0_quadratic_w, f1_w);CHKERRQ(ierr); 416 417 exactFuncs[0] = quadratic_u; 418 exactFuncs[1] = linear_p; 419 exactFuncs[2] = linear_T; 420 break; 421 case SOL_CUBIC: 422 ierr = PetscDSSetResidual(prob, 0, f0_cubic_v, f1_v);CHKERRQ(ierr); 423 ierr = PetscDSSetResidual(prob, 2, f0_cubic_w, f1_w);CHKERRQ(ierr); 424 425 exactFuncs[0] = cubic_u; 426 exactFuncs[1] = quadratic_p; 427 exactFuncs[2] = quadratic_T; 428 break; 429 default: SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%D)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType); 430 } 431 432 ierr = PetscDSSetResidual(prob, 1, f0_q, NULL);CHKERRQ(ierr); 433 434 ierr = PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu);CHKERRQ(ierr); 435 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL);CHKERRQ(ierr); 436 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL);CHKERRQ(ierr); 437 ierr = PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL);CHKERRQ(ierr); 438 ierr = PetscDSSetJacobian(prob, 2, 2, NULL, g1_wT, NULL, g3_wT);CHKERRQ(ierr); 439 /* Setup constants */ 440 { 441 Parameter *param; 442 PetscScalar constants[3]; 443 444 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 445 446 constants[0] = param->nu; 447 constants[1] = param->alpha; 448 constants[2] = param->theta; 449 ierr = PetscDSSetConstants(prob, 3, constants);CHKERRQ(ierr); 450 } 451 /* Setup Boundary Conditions */ 452 ierr = PetscBagGetData(user->bag, (void **) &ctx);CHKERRQ(ierr); 453 id = 3; 454 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], NULL, ctx, NULL);CHKERRQ(ierr); 455 id = 1; 456 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], NULL, ctx, NULL);CHKERRQ(ierr); 457 id = 2; 458 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], NULL, ctx, NULL);CHKERRQ(ierr); 459 id = 4; 460 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], NULL, ctx, NULL);CHKERRQ(ierr); 461 id = 3; 462 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], NULL, ctx, NULL);CHKERRQ(ierr); 463 id = 1; 464 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], NULL, ctx, NULL);CHKERRQ(ierr); 465 id = 2; 466 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], NULL, ctx, NULL);CHKERRQ(ierr); 467 id = 4; 468 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void)) exactFuncs[2], NULL, ctx, NULL);CHKERRQ(ierr); 469 470 /*setup exact solution.*/ 471 ierr = PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx);CHKERRQ(ierr); 472 ierr = PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx);CHKERRQ(ierr); 473 ierr = PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx);CHKERRQ(ierr); 474 PetscFunctionReturn(0); 475 } 476 477 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 478 { 479 DM cdm = dm; 480 PetscFE fe[3]; 481 Parameter *param; 482 MPI_Comm comm; 483 PetscInt dim; 484 PetscBool simplex; 485 PetscErrorCode ierr; 486 487 PetscFunctionBeginUser; 488 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 489 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 490 /* Create finite element */ 491 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 492 ierr = PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);CHKERRQ(ierr); 493 ierr = PetscObjectSetName((PetscObject) fe[0], "velocity");CHKERRQ(ierr); 494 495 ierr = PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);CHKERRQ(ierr); 496 ierr = PetscFECopyQuadrature(fe[0], fe[1]);CHKERRQ(ierr); 497 ierr = PetscObjectSetName((PetscObject) fe[1], "pressure");CHKERRQ(ierr); 498 499 ierr = PetscFECreateDefault(comm, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]);CHKERRQ(ierr); 500 ierr = PetscFECopyQuadrature(fe[0], fe[2]);CHKERRQ(ierr); 501 ierr = PetscObjectSetName((PetscObject) fe[2], "temperature");CHKERRQ(ierr); 502 503 /* Set discretization and boundary conditions for each mesh */ 504 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe[0]);CHKERRQ(ierr); 505 ierr = DMSetField(dm, 1, NULL, (PetscObject) fe[1]);CHKERRQ(ierr); 506 ierr = DMSetField(dm, 2, NULL, (PetscObject) fe[2]);CHKERRQ(ierr); 507 ierr = DMCreateDS(dm);CHKERRQ(ierr); 508 ierr = SetupProblem(dm, user);CHKERRQ(ierr); 509 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 510 while (cdm) { 511 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 512 ierr = DMPlexCreateBasisRotation(cdm, param->theta, 0.0, 0.0);CHKERRQ(ierr); 513 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 514 } 515 ierr = PetscFEDestroy(&fe[0]);CHKERRQ(ierr); 516 ierr = PetscFEDestroy(&fe[1]);CHKERRQ(ierr); 517 ierr = PetscFEDestroy(&fe[2]);CHKERRQ(ierr); 518 PetscFunctionReturn(0); 519 } 520 521 static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace) 522 { 523 Vec vec; 524 PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero}; 525 PetscErrorCode ierr; 526 527 PetscFunctionBeginUser; 528 PetscCheckFalse(ofield != 1,PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index 1, not %D", ofield); 529 funcs[nfield] = constant; 530 ierr = DMCreateGlobalVector(dm, &vec);CHKERRQ(ierr); 531 ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);CHKERRQ(ierr); 532 ierr = VecNormalize(vec, NULL);CHKERRQ(ierr); 533 ierr = PetscObjectSetName((PetscObject) vec, "Pressure Null Space");CHKERRQ(ierr); 534 ierr = VecViewFromOptions(vec, NULL, "-pressure_nullspace_view");CHKERRQ(ierr); 535 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_FALSE, 1, &vec, nullSpace);CHKERRQ(ierr); 536 ierr = VecDestroy(&vec);CHKERRQ(ierr); 537 PetscFunctionReturn(0); 538 } 539 540 int main(int argc, char **argv) 541 { 542 SNES snes; /* nonlinear solver */ 543 DM dm; /* problem definition */ 544 Vec u, r; /* solution, residual vectors */ 545 AppCtx user; /* user-defined work context */ 546 PetscErrorCode ierr; 547 548 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 549 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 550 ierr = PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);CHKERRQ(ierr); 551 ierr = SetupParameters(&user);CHKERRQ(ierr); 552 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 553 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 554 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 555 ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr); 556 /* Setup problem */ 557 ierr = SetupDiscretization(dm, &user);CHKERRQ(ierr); 558 ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr); 559 560 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 561 ierr = PetscObjectSetName((PetscObject) u, "Solution");CHKERRQ(ierr); 562 ierr = VecDuplicate(u, &r);CHKERRQ(ierr); 563 564 ierr = DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace);CHKERRQ(ierr); 565 ierr = DMPlexSetSNESLocalFEM(dm,&user,&user,&user);CHKERRQ(ierr); 566 567 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 568 { 569 PetscDS ds; 570 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 571 void *ctxs[3]; 572 573 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 574 ierr = PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]);CHKERRQ(ierr); 575 ierr = PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]);CHKERRQ(ierr); 576 ierr = PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]);CHKERRQ(ierr); 577 ierr = DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u);CHKERRQ(ierr); 578 ierr = PetscObjectSetName((PetscObject) u, "Exact Solution");CHKERRQ(ierr); 579 ierr = VecViewFromOptions(u, NULL, "-exact_vec_view");CHKERRQ(ierr); 580 } 581 ierr = DMSNESCheckFromOptions(snes, u);CHKERRQ(ierr); 582 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 583 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 584 585 if (user.showError) { 586 PetscDS ds; 587 Vec r; 588 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 589 void *ctxs[3]; 590 591 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 592 ierr = PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]);CHKERRQ(ierr); 593 ierr = PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]);CHKERRQ(ierr); 594 ierr = PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]);CHKERRQ(ierr); 595 ierr = DMGetGlobalVector(dm, &r);CHKERRQ(ierr); 596 ierr = DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, r);CHKERRQ(ierr); 597 ierr = VecAXPY(r, -1.0, u);CHKERRQ(ierr); 598 ierr = PetscObjectSetName((PetscObject) r, "Solution Error");CHKERRQ(ierr); 599 ierr = VecViewFromOptions(r, NULL, "-error_vec_view");CHKERRQ(ierr); 600 ierr = DMRestoreGlobalVector(dm, &r);CHKERRQ(ierr); 601 } 602 ierr = PetscObjectSetName((PetscObject) u, "Numerical Solution");CHKERRQ(ierr); 603 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr); 604 605 ierr = VecDestroy(&u);CHKERRQ(ierr); 606 ierr = VecDestroy(&r);CHKERRQ(ierr); 607 ierr = DMDestroy(&dm);CHKERRQ(ierr); 608 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 609 ierr = PetscBagDestroy(&user.bag);CHKERRQ(ierr); 610 ierr = PetscFinalize(); 611 return ierr; 612 } 613 614 /*TEST 615 616 test: 617 suffix: 2d_tri_p2_p1_p1 618 requires: triangle !single 619 args: -dm_plex_separate_marker -sol_type quadratic -dm_refine 0 \ 620 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 621 -dmsnes_check .001 -snes_error_if_not_converged \ 622 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 623 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 624 -fieldsplit_0_pc_type lu \ 625 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 626 627 test: 628 # Using -dm_refine 2 -convest_num_refine 3 gives L_2 convergence rate: [2.9, 2.3, 1.9] 629 suffix: 2d_tri_p2_p1_p1_conv 630 requires: triangle !single 631 args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \ 632 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 633 -snes_error_if_not_converged -snes_convergence_test correct_pressure -snes_convergence_estimate -convest_num_refine 1 \ 634 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 635 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 636 -fieldsplit_0_pc_type lu \ 637 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 638 639 test: 640 suffix: 2d_tri_p3_p2_p2 641 requires: triangle !single 642 args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \ 643 -vel_petscspace_degree 3 -pres_petscspace_degree 2 -temp_petscspace_degree 2 \ 644 -dmsnes_check .001 -snes_error_if_not_converged \ 645 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 646 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 647 -fieldsplit_0_pc_type lu \ 648 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi 649 650 TEST*/ 651