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