1 static char help[] = "Time-dependent reactive low Mach Flow in 2d and 3d channels with finite elements.\n\ 2 We solve the reactive low Mach flow problem in a rectangular domain\n\ 3 using a parallel unstructured mesh (DMPLEX) to discretize the flow\n\ 4 and particles (DWSWARM) to discretize the chemical species.\n\n\n"; 5 6 /*F 7 This low Mach flow is time-dependent isoviscous Navier-Stokes flow. We discretize using the 8 finite element method on an unstructured mesh. The weak form equations are 9 10 \begin{align*} 11 < q, \nabla\cdot u > = 0 12 <v, du/dt> + <v, u \cdot \nabla u> + < \nabla v, \nu (\nabla u + {\nabla u}^T) > - < \nabla\cdot v, p > - < v, f > = 0 13 < w, u \cdot \nabla T > + < \nabla w, \alpha \nabla T > - < w, Q > = 0 14 \end{align*} 15 16 where $\nu$ is the kinematic viscosity and $\alpha$ is thermal diffusivity. 17 18 For visualization, use 19 20 -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append 21 22 The particles can be visualized using 23 24 -part_dm_view draw -part_dm_view_swarm_radius 0.03 25 26 F*/ 27 28 #include <petscdmplex.h> 29 #include <petscdmswarm.h> 30 #include <petscts.h> 31 #include <petscds.h> 32 #include <petscbag.h> 33 34 typedef enum {SOL_TRIG_TRIG, NUM_SOL_TYPES} SolType; 35 const char *solTypes[NUM_SOL_TYPES+1] = {"trig_trig", "unknown"}; 36 37 typedef enum {PART_LAYOUT_CELL, PART_LAYOUT_BOX, NUM_PART_LAYOUT_TYPES} PartLayoutType; 38 const char *partLayoutTypes[NUM_PART_LAYOUT_TYPES+1] = {"cell", "box", "unknown"}; 39 40 typedef struct { 41 PetscReal nu; /* Kinematic viscosity */ 42 PetscReal alpha; /* Thermal diffusivity */ 43 PetscReal T_in; /* Inlet temperature*/ 44 PetscReal omega; /* Rotation speed in MMS benchmark */ 45 } Parameter; 46 47 typedef struct { 48 /* Problem definition */ 49 PetscBag bag; /* Holds problem parameters */ 50 SolType solType; /* MMS solution type */ 51 PartLayoutType partLayout; /* Type of particle distribution */ 52 PetscInt Npc; /* The initial number of particles per cell */ 53 PetscReal partLower[3]; /* Lower left corner of particle box */ 54 PetscReal partUpper[3]; /* Upper right corner of particle box */ 55 PetscInt Npb; /* The initial number of particles per box dimension */ 56 } AppCtx; 57 58 typedef struct { 59 PetscReal ti; /* The time for ui, at the beginning of the advection solve */ 60 PetscReal tf; /* The time for uf, at the end of the advection solve */ 61 Vec ui; /* The PDE solution field at ti */ 62 Vec uf; /* The PDE solution field at tf */ 63 Vec x0; /* The initial particle positions at t = 0 */ 64 PetscErrorCode (*exact)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 65 AppCtx *ctx; /* Context for exact solution */ 66 } AdvCtx; 67 68 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 69 { 70 PetscInt d; 71 for (d = 0; d < Nc; ++d) u[d] = 0.0; 72 return 0; 73 } 74 75 static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 76 { 77 PetscInt d; 78 for (d = 0; d < Nc; ++d) u[d] = 1.0; 79 return 0; 80 } 81 82 /* 83 CASE: trigonometric-trigonometric 84 In 2D we use exact solution: 85 86 x = r0 cos(w t + theta0) r0 = sqrt(x0^2 + y0^2) 87 y = r0 sin(w t + theta0) theta0 = arctan(y0/x0) 88 u = -w r0 sin(theta0) = -w y 89 v = w r0 cos(theta0) = w x 90 p = x + y - 1 91 T = t + x + y 92 f = <1, 1> 93 Q = 1 + w (x - y)/r 94 95 so that 96 97 \nabla \cdot u = 0 + 0 = 0 98 99 f = du/dt + u \cdot \nabla u - \nu \Delta u + \nabla p 100 = <0, 0> + u_i d_i u_j - \nu 0 + <1, 1> 101 = <1, 1> + w^2 <-y, x> . <<0, 1>, <-1, 0>> 102 = <1, 1> + w^2 <-x, -y> 103 = <1, 1> - w^2 <x, y> 104 105 Q = dT/dt + u \cdot \nabla T - \alpha \Delta T 106 = 1 + <u, v> . <1, 1> - \alpha 0 107 = 1 + u + v 108 */ 109 static PetscErrorCode trig_trig_x(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *x, void *ctx) 110 { 111 const PetscReal x0 = X[0]; 112 const PetscReal y0 = X[1]; 113 const PetscReal R0 = PetscSqrtReal(x0*x0 + y0*y0); 114 const PetscReal theta0 = PetscAtan2Real(y0, x0); 115 Parameter *p = (Parameter *) ctx; 116 117 x[0] = R0*PetscCosReal(p->omega*time + theta0); 118 x[1] = R0*PetscSinReal(p->omega*time + theta0); 119 return 0; 120 } 121 static PetscErrorCode trig_trig_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 122 { 123 Parameter *p = (Parameter *) ctx; 124 125 u[0] = -p->omega*X[1]; 126 u[1] = p->omega*X[0]; 127 return 0; 128 } 129 static PetscErrorCode trig_trig_u_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx) 130 { 131 u[0] = 0.0; 132 u[1] = 0.0; 133 return 0; 134 } 135 136 static PetscErrorCode trig_trig_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx) 137 { 138 p[0] = X[0] + X[1] - 1.0; 139 return 0; 140 } 141 142 static PetscErrorCode trig_trig_T(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 143 { 144 T[0] = time + X[0] + X[1]; 145 return 0; 146 } 147 static PetscErrorCode trig_trig_T_t(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx) 148 { 149 T[0] = 1.0; 150 return 0; 151 } 152 153 static void f0_trig_trig_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 154 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 155 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 156 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 157 { 158 const PetscReal omega = PetscRealPart(constants[3]); 159 PetscInt Nc = dim; 160 PetscInt c, d; 161 162 for (d = 0; d < dim; ++d) f0[d] = u_t[uOff[0]+d]; 163 164 for (c = 0; c < Nc; ++c) { 165 for (d = 0; d < dim; ++d) f0[c] += u[d]*u_x[c*dim+d]; 166 } 167 f0[0] -= 1.0 - omega*omega*X[0]; 168 f0[1] -= 1.0 - omega*omega*X[1]; 169 } 170 171 static void f0_trig_trig_w(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 const PetscReal omega = PetscRealPart(constants[3]); 177 PetscInt d; 178 179 for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0]+d]*u_x[uOff_x[2]+d]; 180 f0[0] += u_t[uOff[2]] - (1.0 + omega*(X[0] - X[1])); 181 } 182 183 static void f0_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, 184 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 185 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 186 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 187 { 188 PetscInt d; 189 for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d*dim+d]; 190 } 191 192 /*f1_v = \nu[grad(u) + grad(u)^T] - pI */ 193 static void f1_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, 194 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 195 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 196 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 197 { 198 const PetscReal nu = PetscRealPart(constants[0]); 199 const PetscInt Nc = dim; 200 PetscInt c, d; 201 202 for (c = 0; c < Nc; ++c) { 203 for (d = 0; d < dim; ++d) { 204 f1[c*dim+d] = nu*(u_x[c*dim+d] + u_x[d*dim+c]); 205 //f1[c*dim+d] = nu*u_x[c*dim+d]; 206 } 207 f1[c*dim+c] -= u[uOff[1]]; 208 } 209 } 210 211 static void f1_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, 212 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 213 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 214 PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 215 { 216 const PetscReal alpha = PetscRealPart(constants[1]); 217 PetscInt d; 218 for (d = 0; d < dim; ++d) f1[d] = alpha*u_x[uOff_x[2]+d]; 219 } 220 221 /*Jacobians*/ 222 static void g1_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 223 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 224 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 225 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 226 { 227 PetscInt d; 228 for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; 229 } 230 231 static void g0_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 232 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 233 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 234 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 235 { 236 PetscInt c, d; 237 const PetscInt Nc = dim; 238 239 for (d = 0; d < dim; ++d) g0[d*dim+d] = u_tShift; 240 241 for (c = 0; c < Nc; ++c) { 242 for (d = 0; d < dim; ++d) { 243 g0[c*Nc+d] += u_x[c*Nc+d]; 244 } 245 } 246 } 247 248 static void g1_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 249 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 250 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 251 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 252 { 253 PetscInt NcI = dim; 254 PetscInt NcJ = dim; 255 PetscInt c, d, e; 256 257 for (c = 0; c < NcI; ++c) { 258 for (d = 0; d < NcJ; ++d) { 259 for (e = 0; e < dim; ++e) { 260 if (c == d) { 261 g1[(c*NcJ+d)*dim+e] += u[e]; 262 } 263 } 264 } 265 } 266 } 267 268 269 static void g2_vp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 270 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 271 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 272 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 273 { 274 PetscInt d; 275 for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; 276 } 277 278 static void g3_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 279 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 280 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 281 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 282 { 283 const PetscReal nu = PetscRealPart(constants[0]); 284 const PetscInt Nc = dim; 285 PetscInt c, d; 286 287 for (c = 0; c < Nc; ++c) { 288 for (d = 0; d < dim; ++d) { 289 g3[((c*Nc+c)*dim+d)*dim+d] += nu; // gradU 290 g3[((c*Nc+d)*dim+d)*dim+c] += nu; // gradU transpose 291 } 292 } 293 } 294 295 static void g0_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 296 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 297 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 298 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 299 { 300 PetscInt d; 301 for (d = 0; d < dim; ++d) g0[d] = u_tShift; 302 } 303 304 static void g0_wu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 305 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 306 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 307 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 308 { 309 PetscInt d; 310 for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2]+d]; 311 } 312 313 static void g1_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 314 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 315 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 316 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 317 { 318 PetscInt d; 319 for (d = 0; d < dim; ++d) g1[d] = u[uOff[0]+d]; 320 } 321 322 static void g3_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, 323 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 324 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 325 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 326 { 327 const PetscReal alpha = PetscRealPart(constants[1]); 328 PetscInt d; 329 330 for (d = 0; d < dim; ++d) g3[d*dim+d] = alpha; 331 } 332 333 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 334 { 335 PetscInt sol, pl, n; 336 PetscErrorCode ierr; 337 338 PetscFunctionBeginUser; 339 options->solType = SOL_TRIG_TRIG; 340 options->partLayout = PART_LAYOUT_CELL; 341 options->Npc = 1; 342 options->Npb = 1; 343 344 options->partLower[0] = options->partLower[1] = options->partLower[2] = 0.; 345 options->partUpper[0] = options->partUpper[1] = options->partUpper[2] = 1.; 346 ierr = PetscOptionsBegin(comm, "", "Low Mach flow Problem Options", "DMPLEX");CHKERRQ(ierr); 347 sol = options->solType; 348 ierr = PetscOptionsEList("-sol_type", "The solution type", "ex77.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL);CHKERRQ(ierr); 349 options->solType = (SolType) sol; 350 pl = options->partLayout; 351 ierr = PetscOptionsEList("-part_layout_type", "The particle layout type", "ex77.c", partLayoutTypes, NUM_PART_LAYOUT_TYPES, partLayoutTypes[options->partLayout], &pl, NULL);CHKERRQ(ierr); 352 options->partLayout = (PartLayoutType) pl; 353 ierr = PetscOptionsInt("-Npc", "The initial number of particles per cell", "ex77.c", options->Npc, &options->Npc, NULL);CHKERRQ(ierr); 354 n = 3; 355 ierr = PetscOptionsRealArray("-part_lower", "The lower left corner of the particle box", "ex77.c", options->partLower, &n, NULL);CHKERRQ(ierr); 356 n = 3; 357 ierr = PetscOptionsRealArray("-part_upper", "The upper right corner of the particle box", "ex77.c", options->partUpper, &n, NULL);CHKERRQ(ierr); 358 ierr = PetscOptionsInt("-Npb", "The initial number of particles per box dimension", "ex77.c", options->Npb, &options->Npb, NULL);CHKERRQ(ierr); 359 ierr = PetscOptionsEnd(); 360 PetscFunctionReturn(0); 361 } 362 363 static PetscErrorCode SetupParameters(AppCtx *user) 364 { 365 PetscBag bag; 366 Parameter *p; 367 PetscErrorCode ierr; 368 369 PetscFunctionBeginUser; 370 /* setup PETSc parameter bag */ 371 ierr = PetscBagGetData(user->bag, (void **) &p);CHKERRQ(ierr); 372 ierr = PetscBagSetName(user->bag, "par", "Low Mach flow parameters");CHKERRQ(ierr); 373 bag = user->bag; 374 ierr = PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity");CHKERRQ(ierr); 375 ierr = PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity");CHKERRQ(ierr); 376 ierr = PetscBagRegisterReal(bag, &p->T_in, 1.0, "T_in", "Inlet temperature");CHKERRQ(ierr); 377 ierr = PetscBagRegisterReal(bag, &p->omega, 1.0, "omega", "Rotation speed in MMS benchmark");CHKERRQ(ierr); 378 PetscFunctionReturn(0); 379 } 380 381 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 382 { 383 PetscErrorCode ierr; 384 385 PetscFunctionBeginUser; 386 ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 387 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 388 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 389 PetscFunctionReturn(0); 390 } 391 392 static PetscErrorCode SetupProblem(DM dm, AppCtx *user) 393 { 394 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 395 PetscErrorCode (*exactFuncs_t[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 396 PetscDS prob; 397 Parameter *ctx; 398 PetscInt id; 399 PetscErrorCode ierr; 400 401 PetscFunctionBeginUser; 402 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 403 switch(user->solType){ 404 case SOL_TRIG_TRIG: 405 ierr = PetscDSSetResidual(prob, 0, f0_trig_trig_v, f1_v);CHKERRQ(ierr); 406 ierr = PetscDSSetResidual(prob, 2, f0_trig_trig_w, f1_w);CHKERRQ(ierr); 407 408 exactFuncs[0] = trig_trig_u; 409 exactFuncs[1] = trig_trig_p; 410 exactFuncs[2] = trig_trig_T; 411 exactFuncs_t[0] = trig_trig_u_t; 412 exactFuncs_t[1] = NULL; 413 exactFuncs_t[2] = trig_trig_T_t; 414 break; 415 default: SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%D)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType); 416 } 417 418 ierr = PetscDSSetResidual(prob, 1, f0_q, NULL);CHKERRQ(ierr); 419 420 ierr = PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu);CHKERRQ(ierr); 421 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL);CHKERRQ(ierr); 422 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL);CHKERRQ(ierr); 423 ierr = PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL);CHKERRQ(ierr); 424 ierr = PetscDSSetJacobian(prob, 2, 2, g0_wT, g1_wT, NULL, g3_wT);CHKERRQ(ierr); 425 /* Setup constants */ 426 { 427 Parameter *param; 428 PetscScalar constants[4]; 429 430 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 431 432 constants[0] = param->nu; 433 constants[1] = param->alpha; 434 constants[2] = param->T_in; 435 constants[3] = param->omega; 436 ierr = PetscDSSetConstants(prob, 4, constants);CHKERRQ(ierr); 437 } 438 /* Setup Boundary Conditions */ 439 ierr = PetscBagGetData(user->bag, (void **) &ctx);CHKERRQ(ierr); 440 id = 3; 441 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall velocity", "marker", 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], 1, &id, ctx);CHKERRQ(ierr); 442 id = 1; 443 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall velocity", "marker", 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], 1, &id, ctx);CHKERRQ(ierr); 444 id = 2; 445 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall velocity", "marker", 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], 1, &id, ctx);CHKERRQ(ierr); 446 id = 4; 447 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall velocity", "marker", 0, 0, NULL, (void (*)(void)) exactFuncs[0], (void (*)(void)) exactFuncs_t[0], 1, &id, ctx);CHKERRQ(ierr); 448 id = 3; 449 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall temp", "marker", 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], 1, &id, ctx);CHKERRQ(ierr); 450 id = 1; 451 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall temp", "marker", 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], 1, &id, ctx);CHKERRQ(ierr); 452 id = 2; 453 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall temp", "marker", 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], 1, &id, ctx);CHKERRQ(ierr); 454 id = 4; 455 ierr = PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall temp", "marker", 2, 0, NULL, (void (*)(void)) exactFuncs[2], (void (*)(void)) exactFuncs_t[2], 1, &id, ctx);CHKERRQ(ierr); 456 457 /*setup exact solution.*/ 458 ierr = PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx);CHKERRQ(ierr); 459 ierr = PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx);CHKERRQ(ierr); 460 ierr = PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx);CHKERRQ(ierr); 461 ierr = PetscDSSetExactSolutionTimeDerivative(prob, 0, exactFuncs_t[0], ctx);CHKERRQ(ierr); 462 ierr = PetscDSSetExactSolutionTimeDerivative(prob, 1, exactFuncs_t[1], ctx);CHKERRQ(ierr); 463 ierr = PetscDSSetExactSolutionTimeDerivative(prob, 2, exactFuncs_t[2], ctx);CHKERRQ(ierr); 464 PetscFunctionReturn(0); 465 } 466 467 /* x_t = v 468 469 Note that here we use the velocity field at t_{n+1} to advect the particles from 470 t_n to t_{n+1}. If we use both of these fields, we could use Crank-Nicholson or 471 the method of characteristics. 472 */ 473 static PetscErrorCode FreeStreaming(TS ts, PetscReal t, Vec X, Vec F, void *ctx) 474 { 475 AdvCtx *adv = (AdvCtx *) ctx; 476 Vec u = adv->ui; 477 DM sdm, dm, vdm; 478 Vec vel, locvel, pvel; 479 IS vis; 480 DMInterpolationInfo ictx; 481 const PetscScalar *coords, *v; 482 PetscScalar *f; 483 PetscInt vf[1] = {0}; 484 PetscInt dim, Np; 485 PetscErrorCode ierr; 486 487 PetscFunctionBeginUser; 488 ierr = TSGetDM(ts, &sdm);CHKERRQ(ierr); 489 ierr = DMSwarmGetCellDM(sdm, &dm);CHKERRQ(ierr); 490 ierr = DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor);CHKERRQ(ierr); 491 ierr = DMGetGlobalVector(sdm, &pvel);CHKERRQ(ierr); 492 ierr = DMSwarmGetLocalSize(sdm, &Np);CHKERRQ(ierr); 493 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 494 /* Get local velocity */ 495 ierr = DMCreateSubDM(dm, 1, vf, &vis, &vdm);CHKERRQ(ierr); 496 ierr = VecGetSubVector(u, vis, &vel);CHKERRQ(ierr); 497 ierr = DMGetLocalVector(vdm, &locvel);CHKERRQ(ierr); 498 ierr = DMPlexInsertBoundaryValues(vdm, PETSC_TRUE, locvel, adv->ti, NULL, NULL, NULL);CHKERRQ(ierr); 499 ierr = DMGlobalToLocalBegin(vdm, vel, INSERT_VALUES, locvel);CHKERRQ(ierr); 500 ierr = DMGlobalToLocalEnd(vdm, vel, INSERT_VALUES, locvel);CHKERRQ(ierr); 501 ierr = VecRestoreSubVector(u, vis, &vel);CHKERRQ(ierr); 502 ierr = ISDestroy(&vis);CHKERRQ(ierr); 503 /* Interpolate velocity */ 504 ierr = DMInterpolationCreate(PETSC_COMM_SELF, &ictx);CHKERRQ(ierr); 505 ierr = DMInterpolationSetDim(ictx, dim);CHKERRQ(ierr); 506 ierr = DMInterpolationSetDof(ictx, dim);CHKERRQ(ierr); 507 ierr = VecGetArrayRead(X, &coords);CHKERRQ(ierr); 508 ierr = DMInterpolationAddPoints(ictx, Np, (PetscReal *) coords);CHKERRQ(ierr); 509 ierr = VecRestoreArrayRead(X, &coords);CHKERRQ(ierr); 510 /* Particles that lie outside the domain should be dropped, 511 whereas particles that move to another partition should trigger a migration */ 512 ierr = DMInterpolationSetUp(ictx, vdm, PETSC_FALSE, PETSC_TRUE);CHKERRQ(ierr); 513 ierr = VecSet(pvel, 0.);CHKERRQ(ierr); 514 ierr = DMInterpolationEvaluate(ictx, vdm, locvel, pvel);CHKERRQ(ierr); 515 ierr = DMInterpolationDestroy(&ictx);CHKERRQ(ierr); 516 ierr = DMRestoreLocalVector(vdm, &locvel);CHKERRQ(ierr); 517 ierr = DMDestroy(&vdm);CHKERRQ(ierr); 518 519 ierr = VecGetArray(F, &f);CHKERRQ(ierr); 520 ierr = VecGetArrayRead(pvel, &v);CHKERRQ(ierr); 521 ierr = PetscArraycpy(f, v, Np*dim);CHKERRQ(ierr); 522 ierr = VecRestoreArrayRead(pvel, &v);CHKERRQ(ierr); 523 ierr = VecRestoreArray(F, &f);CHKERRQ(ierr); 524 ierr = DMRestoreGlobalVector(sdm, &pvel);CHKERRQ(ierr); 525 PetscFunctionReturn(0); 526 } 527 528 static PetscErrorCode SetInitialParticleConditions(TS ts, Vec u) 529 { 530 AppCtx *user; 531 void *ctx; 532 DM dm; 533 PetscScalar *coords; 534 PetscReal x[3], dx[3]; 535 PetscInt n[3]; 536 PetscInt Np, dim, d, i, j, k; 537 PetscErrorCode ierr; 538 539 PetscFunctionBegin; 540 ierr = TSGetApplicationContext(ts, &ctx);CHKERRQ(ierr); 541 user = ((AdvCtx *) ctx)->ctx; 542 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 543 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 544 switch (user->partLayout) { 545 case PART_LAYOUT_CELL: 546 ierr = DMSwarmSetPointCoordinatesRandom(dm, user->Npc);CHKERRQ(ierr); 547 break; 548 case PART_LAYOUT_BOX: 549 Np = 1; 550 for (d = 0; d < dim; ++d) { 551 n[d] = user->Npb; 552 dx[d] = (user->partUpper[d] - user->partLower[d])/PetscMax(1, n[d] - 1); 553 Np *= n[d]; 554 } 555 ierr = VecGetArray(u, &coords);CHKERRQ(ierr); 556 switch (dim) { 557 case 2: 558 x[0] = user->partLower[0]; 559 for (i = 0; i < n[0]; ++i, x[0] += dx[0]) { 560 x[1] = user->partLower[1]; 561 for (j = 0; j < n[1]; ++j, x[1] += dx[1]) { 562 const PetscInt p = j*n[0] + i; 563 for (d = 0; d < dim; ++d) coords[p*dim + d] = x[d]; 564 } 565 } 566 break; 567 case 3: 568 x[0] = user->partLower[0]; 569 for (i = 0; i < n[0]; ++i, x[0] += dx[0]) { 570 x[1] = user->partLower[1]; 571 for (j = 0; j < n[1]; ++j, x[1] += dx[1]) { 572 x[2] = user->partLower[2]; 573 for (k = 0; k < n[2]; ++k, x[2] += dx[2]) { 574 const PetscInt p = (k*n[1] + j)*n[0] + i; 575 for (d = 0; d < dim; ++d) coords[p*dim + d] = x[d]; 576 } 577 } 578 } 579 break; 580 default: SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_SUP, "Do not support particle layout in dimension %D", dim); 581 } 582 ierr = VecRestoreArray(u, &coords);CHKERRQ(ierr); 583 break; 584 default: SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_ARG_WRONG, "Invalid particle layout type %s", partLayoutTypes[PetscMin(user->partLayout, NUM_PART_LAYOUT_TYPES)]); 585 } 586 PetscFunctionReturn(0); 587 } 588 589 static PetscErrorCode SetupDiscretization(DM dm, DM sdm, AppCtx *user) 590 { 591 DM cdm = dm; 592 PetscFE fe[3]; 593 Parameter *param; 594 PetscInt *cellid, n[3]; 595 PetscReal x[3], dx[3]; 596 PetscScalar *coords; 597 DMPolytopeType ct; 598 PetscInt dim, d, cStart, cEnd, c, Np, p, i, j, k; 599 PetscBool simplex; 600 MPI_Comm comm; 601 PetscErrorCode ierr; 602 603 PetscFunctionBeginUser; 604 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 605 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 606 ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr); 607 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE; 608 /* Create finite element */ 609 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 610 ierr = PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);CHKERRQ(ierr); 611 ierr = PetscObjectSetName((PetscObject) fe[0], "velocity");CHKERRQ(ierr); 612 613 ierr = PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);CHKERRQ(ierr); 614 ierr = PetscFECopyQuadrature(fe[0], fe[1]);CHKERRQ(ierr); 615 ierr = PetscObjectSetName((PetscObject) fe[1], "pressure");CHKERRQ(ierr); 616 617 ierr = PetscFECreateDefault(comm, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]);CHKERRQ(ierr); 618 ierr = PetscFECopyQuadrature(fe[0], fe[2]);CHKERRQ(ierr); 619 ierr = PetscObjectSetName((PetscObject) fe[2], "temperature");CHKERRQ(ierr); 620 621 /* Set discretization and boundary conditions for each mesh */ 622 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe[0]);CHKERRQ(ierr); 623 ierr = DMSetField(dm, 1, NULL, (PetscObject) fe[1]);CHKERRQ(ierr); 624 ierr = DMSetField(dm, 2, NULL, (PetscObject) fe[2]);CHKERRQ(ierr); 625 ierr = DMCreateDS(dm);CHKERRQ(ierr); 626 ierr = SetupProblem(dm, user);CHKERRQ(ierr); 627 ierr = PetscBagGetData(user->bag, (void **) ¶m);CHKERRQ(ierr); 628 while (cdm) { 629 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 630 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 631 } 632 ierr = PetscFEDestroy(&fe[0]);CHKERRQ(ierr); 633 ierr = PetscFEDestroy(&fe[1]);CHKERRQ(ierr); 634 ierr = PetscFEDestroy(&fe[2]);CHKERRQ(ierr); 635 636 { 637 PetscObject pressure; 638 MatNullSpace nullspacePres; 639 640 ierr = DMGetField(dm, 1, NULL, &pressure);CHKERRQ(ierr); 641 ierr = MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullspacePres);CHKERRQ(ierr); 642 ierr = PetscObjectCompose(pressure, "nullspace", (PetscObject) nullspacePres);CHKERRQ(ierr); 643 ierr = MatNullSpaceDestroy(&nullspacePres);CHKERRQ(ierr); 644 } 645 646 /* Setup particle information */ 647 ierr = DMSwarmSetType(sdm, DMSWARM_PIC);CHKERRQ(ierr); 648 ierr = DMSwarmRegisterPetscDatatypeField(sdm, "mass", 1, PETSC_REAL);CHKERRQ(ierr); 649 ierr = DMSwarmFinalizeFieldRegister(sdm);CHKERRQ(ierr); 650 switch (user->partLayout) { 651 case PART_LAYOUT_CELL: 652 ierr = DMSwarmSetLocalSizes(sdm, (cEnd - cStart) * user->Npc, 0);CHKERRQ(ierr); 653 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr); 654 ierr = DMSwarmGetField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 655 for (c = cStart; c < cEnd; ++c) { 656 for (p = 0; p < user->Npc; ++p) { 657 const PetscInt n = c*user->Npc + p; 658 659 cellid[n] = c; 660 } 661 } 662 ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 663 ierr = DMSwarmSetPointCoordinatesRandom(sdm, user->Npc);CHKERRQ(ierr); 664 break; 665 case PART_LAYOUT_BOX: 666 Np = 1; 667 for (d = 0; d < dim; ++d) { 668 n[d] = user->Npb; 669 dx[d] = (user->partUpper[d] - user->partLower[d])/PetscMax(1, n[d] - 1); 670 Np *= n[d]; 671 } 672 ierr = DMSwarmSetLocalSizes(sdm, Np, 0);CHKERRQ(ierr); 673 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr); 674 ierr = DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 675 switch (dim) { 676 case 2: 677 x[0] = user->partLower[0]; 678 for (i = 0; i < n[0]; ++i, x[0] += dx[0]) { 679 x[1] = user->partLower[1]; 680 for (j = 0; j < n[1]; ++j, x[1] += dx[1]) { 681 const PetscInt p = j*n[0] + i; 682 for (d = 0; d < dim; ++d) coords[p*dim + d] = x[d]; 683 } 684 } 685 break; 686 case 3: 687 x[0] = user->partLower[0]; 688 for (i = 0; i < n[0]; ++i, x[0] += dx[0]) { 689 x[1] = user->partLower[1]; 690 for (j = 0; j < n[1]; ++j, x[1] += dx[1]) { 691 x[2] = user->partLower[2]; 692 for (k = 0; k < n[2]; ++k, x[2] += dx[2]) { 693 const PetscInt p = (k*n[1] + j)*n[0] + i; 694 for (d = 0; d < dim; ++d) coords[p*dim + d] = x[d]; 695 } 696 } 697 } 698 break; 699 default: SETERRQ1(comm, PETSC_ERR_SUP, "Do not support particle layout in dimension %D", dim); 700 } 701 ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr); 702 ierr = DMSwarmGetField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 703 for (p = 0; p < Np; ++p) cellid[p] = 0; 704 ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr); 705 ierr = DMSwarmMigrate(sdm, PETSC_TRUE);CHKERRQ(ierr); 706 break; 707 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Invalid particle layout type %s", partLayoutTypes[PetscMin(user->partLayout, NUM_PART_LAYOUT_TYPES)]); 708 } 709 ierr = PetscObjectSetName((PetscObject) sdm, "Particles");CHKERRQ(ierr); 710 ierr = DMViewFromOptions(sdm, NULL, "-dm_view");CHKERRQ(ierr); 711 PetscFunctionReturn(0); 712 } 713 714 static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace) 715 { 716 Vec vec; 717 PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero}; 718 PetscErrorCode ierr; 719 720 PetscFunctionBeginUser; 721 if (ofield != 1) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index 1, not %D", ofield); 722 funcs[nfield] = constant; 723 ierr = DMCreateGlobalVector(dm, &vec);CHKERRQ(ierr); 724 ierr = DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);CHKERRQ(ierr); 725 ierr = VecNormalize(vec, NULL);CHKERRQ(ierr); 726 ierr = PetscObjectSetName((PetscObject) vec, "Pressure Null Space");CHKERRQ(ierr); 727 ierr = VecViewFromOptions(vec, NULL, "-pressure_nullspace_view");CHKERRQ(ierr); 728 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_FALSE, 1, &vec, nullSpace);CHKERRQ(ierr); 729 ierr = VecDestroy(&vec);CHKERRQ(ierr); 730 PetscFunctionReturn(0); 731 } 732 733 static PetscErrorCode RemoveDiscretePressureNullspace_Private(TS ts, Vec u) 734 { 735 DM dm; 736 MatNullSpace nullsp; 737 PetscErrorCode ierr; 738 739 PetscFunctionBegin; 740 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 741 ierr = CreatePressureNullSpace(dm, 1, 1, &nullsp);CHKERRQ(ierr); 742 ierr = MatNullSpaceRemove(nullsp, u);CHKERRQ(ierr); 743 ierr = MatNullSpaceDestroy(&nullsp);CHKERRQ(ierr); 744 PetscFunctionReturn(0); 745 } 746 747 /* Make the discrete pressure discretely divergence free */ 748 static PetscErrorCode RemoveDiscretePressureNullspace(TS ts) 749 { 750 Vec u; 751 PetscErrorCode ierr; 752 753 PetscFunctionBegin; 754 ierr = TSGetSolution(ts, &u);CHKERRQ(ierr); 755 ierr = RemoveDiscretePressureNullspace_Private(ts, u);CHKERRQ(ierr); 756 PetscFunctionReturn(0); 757 } 758 759 static PetscErrorCode SetInitialConditions(TS ts, Vec u) 760 { 761 DM dm; 762 PetscReal t; 763 PetscErrorCode ierr; 764 765 PetscFunctionBegin; 766 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 767 ierr = TSGetTime(ts, &t);CHKERRQ(ierr); 768 ierr = DMComputeExactSolution(dm, t, u, NULL);CHKERRQ(ierr); 769 ierr = RemoveDiscretePressureNullspace_Private(ts, u);CHKERRQ(ierr); 770 PetscFunctionReturn(0); 771 } 772 773 static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) 774 { 775 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 776 void *ctxs[3]; 777 DM dm; 778 PetscDS ds; 779 Vec v; 780 PetscReal ferrors[3]; 781 PetscInt tl, l, f; 782 PetscErrorCode ierr; 783 784 PetscFunctionBeginUser; 785 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 786 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 787 788 for (f = 0; f < 3; ++f) {ierr = PetscDSGetExactSolution(ds, f, &exactFuncs[f], &ctxs[f]);CHKERRQ(ierr);} 789 ierr = DMComputeL2FieldDiff(dm, crtime, exactFuncs, ctxs, u, ferrors);CHKERRQ(ierr); 790 ierr = PetscObjectGetTabLevel((PetscObject) ts, &tl);CHKERRQ(ierr); 791 for (l = 0; l < tl; ++l) {ierr = PetscPrintf(PETSC_COMM_WORLD, "\t");CHKERRQ(ierr);} 792 ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g, %2.3g]\n", (int) step, (double) crtime, (double) ferrors[0], (double) ferrors[1], (double) ferrors[2]);CHKERRQ(ierr); 793 794 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 795 //ierr = TSGetSolution(ts, &u);CHKERRQ(ierr); 796 ierr = PetscObjectSetName((PetscObject) u, "Numerical Solution");CHKERRQ(ierr); 797 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr); 798 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 799 800 ierr = DMGetGlobalVector(dm, &v);CHKERRQ(ierr); 801 // ierr = VecSet(v, 0.0);CHKERRQ(ierr); 802 ierr = DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, v);CHKERRQ(ierr); 803 ierr = PetscObjectSetName((PetscObject) v, "Exact Solution");CHKERRQ(ierr); 804 ierr = VecViewFromOptions(v, NULL, "-exact_vec_view");CHKERRQ(ierr); 805 ierr = DMRestoreGlobalVector(dm, &v);CHKERRQ(ierr); 806 807 PetscFunctionReturn(0); 808 } 809 810 static PetscErrorCode ComputeParticleError(TS ts, Vec u, Vec e) 811 { 812 AdvCtx *adv; 813 DM sdm; 814 Parameter *param; 815 const PetscScalar *xp0, *xp; 816 PetscScalar *ep; 817 PetscReal time; 818 PetscInt dim, Np, p; 819 MPI_Comm comm; 820 PetscErrorCode ierr; 821 822 PetscFunctionBeginUser; 823 ierr = TSGetTime(ts, &time);CHKERRQ(ierr); 824 ierr = TSGetApplicationContext(ts, (void **) &adv);CHKERRQ(ierr); 825 ierr = PetscBagGetData(adv->ctx->bag, (void **) ¶m);CHKERRQ(ierr); 826 ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr); 827 ierr = TSGetDM(ts, &sdm);CHKERRQ(ierr); 828 ierr = DMGetDimension(sdm, &dim);CHKERRQ(ierr); 829 ierr = DMSwarmGetLocalSize(sdm, &Np);CHKERRQ(ierr); 830 ierr = VecGetArrayRead(adv->x0, &xp0);CHKERRQ(ierr); 831 ierr = VecGetArrayRead(u, &xp);CHKERRQ(ierr); 832 ierr = VecGetArrayWrite(e, &ep);CHKERRQ(ierr); 833 for (p = 0; p < Np; ++p) { 834 PetscScalar x[3]; 835 PetscReal x0[3]; 836 PetscInt d; 837 838 for (d = 0; d < dim; ++d) x0[d] = PetscRealPart(xp0[p*dim+d]); 839 ierr = adv->exact(dim, time, x0, 1, x, param);CHKERRQ(ierr); 840 for (d = 0; d < dim; ++d) ep[p*dim+d] += x[d] - xp[p*dim+d]; 841 } 842 ierr = VecRestoreArrayRead(adv->x0, &xp0);CHKERRQ(ierr); 843 ierr = VecRestoreArrayRead(u, &xp);CHKERRQ(ierr); 844 ierr = VecRestoreArrayWrite(e, &ep);CHKERRQ(ierr); 845 PetscFunctionReturn(0); 846 } 847 848 static PetscErrorCode MonitorParticleError(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx) 849 { 850 AdvCtx *adv = (AdvCtx *) ctx; 851 DM sdm; 852 Parameter *param; 853 const PetscScalar *xp0, *xp; 854 PetscReal error = 0.0; 855 PetscInt dim, tl, l, Np, p; 856 MPI_Comm comm; 857 PetscErrorCode ierr; 858 859 PetscFunctionBeginUser; 860 ierr = PetscBagGetData(adv->ctx->bag, (void **) ¶m);CHKERRQ(ierr); 861 ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr); 862 ierr = TSGetDM(ts, &sdm);CHKERRQ(ierr); 863 ierr = DMGetDimension(sdm, &dim);CHKERRQ(ierr); 864 ierr = DMSwarmGetLocalSize(sdm, &Np);CHKERRQ(ierr); 865 ierr = VecGetArrayRead(adv->x0, &xp0);CHKERRQ(ierr); 866 ierr = VecGetArrayRead(u, &xp);CHKERRQ(ierr); 867 for (p = 0; p < Np; ++p) { 868 PetscScalar x[3]; 869 PetscReal x0[3]; 870 PetscReal perror = 0.0; 871 PetscInt d; 872 873 for (d = 0; d < dim; ++d) x0[d] = PetscRealPart(xp0[p*dim+d]); 874 ierr = adv->exact(dim, time, x0, 1, x, param);CHKERRQ(ierr); 875 for (d = 0; d < dim; ++d) perror += PetscSqr(PetscRealPart(x[d] - xp[p*dim+d])); 876 error += perror; 877 } 878 ierr = VecRestoreArrayRead(adv->x0, &xp0);CHKERRQ(ierr); 879 ierr = VecRestoreArrayRead(u, &xp);CHKERRQ(ierr); 880 ierr = PetscObjectGetTabLevel((PetscObject) ts, &tl);CHKERRQ(ierr); 881 for (l = 0; l < tl; ++l) {ierr = PetscPrintf(PETSC_COMM_WORLD, "\t");CHKERRQ(ierr);} 882 ierr = PetscPrintf(comm, "Timestep: %04d time = %-8.4g \t L_2 Particle Error: [%2.3g]\n", (int) step, (double) time, (double) error);CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 static PetscErrorCode AdvectParticles(TS ts) 887 { 888 TS sts; 889 DM sdm; 890 Vec p; 891 AdvCtx *adv; 892 PetscScalar *coord, *a; 893 const PetscScalar *ca; 894 PetscReal time; 895 PetscInt dim, n, np; 896 PetscErrorCode ierr; 897 898 PetscFunctionBeginUser; 899 ierr = PetscObjectQuery((PetscObject) ts, "_SwarmTS", (PetscObject *) &sts);CHKERRQ(ierr); 900 ierr = PetscObjectQuery((PetscObject) ts, "_SwarmSol", (PetscObject *) &p);CHKERRQ(ierr); 901 ierr = TSGetDM(sts, &sdm);CHKERRQ(ierr); 902 ierr = TSGetRHSFunction(sts, NULL, NULL, (void **) &adv);CHKERRQ(ierr); 903 ierr = DMGetDimension(sdm, &dim);CHKERRQ(ierr); 904 ierr = DMSwarmGetLocalSize(sdm, &np);CHKERRQ(ierr); 905 ierr = VecGetLocalSize(p, &n);CHKERRQ(ierr); 906 if (n != np*dim) { 907 ierr = DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor);CHKERRQ(ierr); 908 ierr = DMCreateGlobalVector(sdm, &p);CHKERRQ(ierr); 909 ierr = PetscObjectCompose((PetscObject) ts, "_SwarmSol", (PetscObject) p);CHKERRQ(ierr); 910 ierr = VecDestroy(&p);CHKERRQ(ierr); 911 ierr = PetscObjectQuery((PetscObject) ts, "_SwarmSol", (PetscObject *) &p);CHKERRQ(ierr); 912 ierr = VecGetLocalSize(p, &n);CHKERRQ(ierr); 913 ierr = TSReset(sts);CHKERRQ(ierr); 914 } 915 916 ierr = DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coord);CHKERRQ(ierr); 917 ierr = VecGetArray(p, &a);CHKERRQ(ierr); 918 ierr = PetscArraycpy(a, coord, n);CHKERRQ(ierr); 919 ierr = VecRestoreArray(p, &a);CHKERRQ(ierr); 920 ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coord);CHKERRQ(ierr); 921 ierr = TSGetTime(ts, &time);CHKERRQ(ierr); 922 ierr = TSSetMaxTime(sts, time);CHKERRQ(ierr); 923 adv->tf = time; 924 ierr = TSSolve(sts, p);CHKERRQ(ierr); 925 ierr = DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coord);CHKERRQ(ierr); 926 ierr = DMSwarmGetLocalSize(sdm, &np);CHKERRQ(ierr); 927 ierr = VecGetLocalSize(p, &n);CHKERRQ(ierr); 928 if (n != np*dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Particle coords %D != %D Cached size", np*dim, n); 929 ierr = VecGetArrayRead(p, &ca);CHKERRQ(ierr); 930 ierr = PetscArraycpy(coord, ca, n);CHKERRQ(ierr); 931 ierr = VecRestoreArrayRead(p, &ca);CHKERRQ(ierr); 932 ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coord);CHKERRQ(ierr); 933 934 ierr = VecCopy(adv->uf, adv->ui);CHKERRQ(ierr); 935 adv->ti = adv->tf; 936 937 ierr = DMSwarmMigrate(sdm, PETSC_TRUE);CHKERRQ(ierr); 938 ierr = DMViewFromOptions(sdm, NULL, "-dm_view");CHKERRQ(ierr); 939 PetscFunctionReturn(0); 940 } 941 942 int main(int argc, char **argv) 943 { 944 DM dm, sdm; 945 TS ts, sts; 946 Vec u, p, xtmp; 947 AppCtx user; 948 AdvCtx adv; 949 PetscReal t; 950 PetscInt dim; 951 PetscErrorCode ierr; 952 953 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 954 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 955 ierr = PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);CHKERRQ(ierr); 956 ierr = SetupParameters(&user);CHKERRQ(ierr); 957 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr); 958 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 959 ierr = TSSetDM(ts, dm);CHKERRQ(ierr); 960 ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr); 961 /* Discretize chemical species */ 962 ierr = DMCreate(PETSC_COMM_WORLD, &sdm);CHKERRQ(ierr); 963 ierr = PetscObjectSetOptionsPrefix((PetscObject) sdm, "part_");CHKERRQ(ierr); 964 ierr = DMSetType(sdm, DMSWARM);CHKERRQ(ierr); 965 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 966 ierr = DMSetDimension(sdm, dim);CHKERRQ(ierr); 967 ierr = DMSwarmSetCellDM(sdm, dm);CHKERRQ(ierr); 968 /* Setup problem */ 969 ierr = SetupDiscretization(dm, sdm, &user);CHKERRQ(ierr); 970 ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr); 971 972 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 973 ierr = DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace);CHKERRQ(ierr); 974 975 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user);CHKERRQ(ierr); 976 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user);CHKERRQ(ierr); 977 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user);CHKERRQ(ierr); 978 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 979 ierr = TSSetPreStep(ts, RemoveDiscretePressureNullspace);CHKERRQ(ierr); 980 ierr = TSMonitorSet(ts, MonitorError, &user, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 981 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 982 983 ierr = TSSetComputeInitialCondition(ts, SetInitialConditions);CHKERRQ(ierr); /* Must come after SetFromOptions() */ 984 ierr = SetInitialConditions(ts, u);CHKERRQ(ierr); 985 ierr = TSGetTime(ts, &t);CHKERRQ(ierr); 986 ierr = DMSetOutputSequenceNumber(dm, 0, t);CHKERRQ(ierr); 987 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr); 988 989 /* Setup particle position integrator */ 990 ierr = TSCreate(PETSC_COMM_WORLD, &sts);CHKERRQ(ierr); 991 ierr = PetscObjectSetOptionsPrefix((PetscObject) sts, "part_");CHKERRQ(ierr); 992 ierr = PetscObjectIncrementTabLevel((PetscObject) sts, (PetscObject) ts, 1);CHKERRQ(ierr); 993 ierr = TSSetDM(sts, sdm);CHKERRQ(ierr); 994 ierr = TSSetProblemType(sts, TS_NONLINEAR);CHKERRQ(ierr); 995 ierr = TSSetExactFinalTime(sts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 996 ierr = TSMonitorSet(sts, MonitorParticleError, &adv, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 997 ierr = TSSetFromOptions(sts);CHKERRQ(ierr); 998 ierr = TSSetApplicationContext(sts, &adv);CHKERRQ(ierr); 999 ierr = TSSetComputeExactError(sts, ComputeParticleError);CHKERRQ(ierr); 1000 ierr = TSSetComputeInitialCondition(sts, SetInitialParticleConditions);CHKERRQ(ierr); 1001 adv.ti = t; 1002 adv.uf = u; 1003 ierr = VecDuplicate(adv.uf, &adv.ui); 1004 ierr = VecCopy(u, adv.ui);CHKERRQ(ierr); 1005 ierr = TSSetRHSFunction(sts, NULL, FreeStreaming, &adv);CHKERRQ(ierr); 1006 ierr = TSSetPostStep(ts, AdvectParticles);CHKERRQ(ierr); 1007 ierr = PetscObjectCompose((PetscObject) ts, "_SwarmTS", (PetscObject) sts);CHKERRQ(ierr); 1008 ierr = DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor);CHKERRQ(ierr); 1009 ierr = DMCreateGlobalVector(sdm, &p);CHKERRQ(ierr); 1010 ierr = PetscObjectCompose((PetscObject) ts, "_SwarmSol", (PetscObject) p);CHKERRQ(ierr); 1011 ierr = VecDestroy(&p);CHKERRQ(ierr); 1012 ierr = DMCreateGlobalVector(sdm, &adv.x0);CHKERRQ(ierr); 1013 ierr = DMSwarmCreateGlobalVectorFromField(sdm, DMSwarmPICField_coor, &xtmp);CHKERRQ(ierr); 1014 ierr = VecCopy(xtmp, adv.x0);CHKERRQ(ierr); 1015 ierr = DMSwarmDestroyGlobalVectorFromField(sdm, DMSwarmPICField_coor, &xtmp);CHKERRQ(ierr); 1016 switch(user.solType){ 1017 case SOL_TRIG_TRIG: adv.exact = trig_trig_x;break; 1018 default: SETERRQ2(PetscObjectComm((PetscObject) sdm), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%D)", solTypes[PetscMin(user.solType, NUM_SOL_TYPES)], user.solType); 1019 } 1020 adv.ctx = &user; 1021 1022 ierr = TSSolve(ts, u);CHKERRQ(ierr); 1023 ierr = DMTSCheckFromOptions(ts, u);CHKERRQ(ierr); 1024 ierr = PetscObjectSetName((PetscObject) u, "Numerical Solution");CHKERRQ(ierr); 1025 1026 ierr = VecDestroy(&u);CHKERRQ(ierr); 1027 ierr = VecDestroy(&adv.x0);CHKERRQ(ierr); 1028 ierr = VecDestroy(&adv.ui);CHKERRQ(ierr); 1029 ierr = DMDestroy(&dm);CHKERRQ(ierr); 1030 ierr = DMDestroy(&sdm);CHKERRQ(ierr); 1031 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1032 ierr = TSDestroy(&sts);CHKERRQ(ierr); 1033 ierr = PetscBagDestroy(&user.bag);CHKERRQ(ierr); 1034 ierr = PetscFinalize(); 1035 return ierr; 1036 } 1037 1038 /*TEST 1039 1040 # Swarm does not work with complex 1041 test: 1042 suffix: 2d_tri_p2_p1_p1_tconvp 1043 requires: triangle !single !complex 1044 args: -dm_plex_separate_marker -sol_type trig_trig -dm_refine 2 \ 1045 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 1046 -dmts_check .001 -ts_max_steps 4 -ts_dt 0.1 -ts_monitor_cancel \ 1047 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 1048 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 1049 -fieldsplit_0_pc_type lu \ 1050 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \ 1051 -omega 0.5 -part_layout_type box -part_lower 0.25,0.25 -part_upper 0.75,0.75 -Npb 5 \ 1052 -part_ts_max_steps 2 -part_ts_dt 0.05 -part_ts_convergence_estimate -convest_num_refine 1 -part_ts_monitor_cancel 1053 test: 1054 suffix: 2d_tri_p2_p1_p1_exit 1055 requires: triangle !single !complex 1056 args: -dm_plex_separate_marker -sol_type trig_trig -dm_refine 1 \ 1057 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \ 1058 -dmts_check .001 -ts_max_steps 10 -ts_dt 0.1 \ 1059 -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \ 1060 -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \ 1061 -fieldsplit_0_pc_type lu \ 1062 -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi \ 1063 -omega 0.5 -part_layout_type box -part_lower 0.25,0.25 -part_upper 0.75,0.75 -Npb 5 \ 1064 -part_ts_max_steps 20 -part_ts_dt 0.05 1065 1066 TEST*/ 1067