1 static char help[] = "Time dependent Navier-Stokes problem in 2d and 3d with finite elements.\n\ 2 We solve the Navier-Stokes in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example supports discretized auxiliary fields (Re) as well as\n\ 5 multilevel nonlinear solvers.\n\ 6 Contributed by: Julian Andrej <juan@tf.uni-kiel.de>\n\n\n"; 7 8 #include <petscdmplex.h> 9 #include <petscsnes.h> 10 #include <petscts.h> 11 #include <petscds.h> 12 13 /* 14 Navier-Stokes equation: 15 16 du/dt + u . grad u - \Delta u - grad p = f 17 div u = 0 18 */ 19 20 typedef struct { 21 PetscInt dim; 22 PetscBool simplex; 23 PetscInt mms; 24 PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 25 } AppCtx; 26 27 #define REYN 400.0 28 29 /* MMS1 30 31 u = t + x^2 + y^2; 32 v = t + 2*x^2 - 2*x*y; 33 p = x + y - 1; 34 35 f_x = -2*t*(x + y) + 2*x*y^2 - 4*x^2*y - 2*x^3 + 4.0/Re - 1.0 36 f_y = -2*t*x + 2*y^3 - 4*x*y^2 - 2*x^2*y + 4.0/Re - 1.0 37 38 so that 39 40 u_t + u \cdot \nabla u - 1/Re \Delta u + \nabla p + f = <1, 1> + <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t 2x + 2x^2y + 4xy^2 - 2y^3> - 1/Re <4, 4> + <1, 1> 41 + <-t (2x + 2y) + 2xy^2 - 4x^2y - 2x^3 + 4/Re - 1, -2xt + 2y^3 - 4xy^2 - 2x^2y + 4/Re - 1> = 0 42 \nabla \cdot u = 2x - 2x = 0 43 44 where 45 46 <u, v> . <<u_x, v_x>, <u_y, v_y>> = <u u_x + v u_y, u v_x + v v_y> 47 */ 48 PetscErrorCode mms1_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 49 { 50 u[0] = time + x[0]*x[0] + x[1]*x[1]; 51 u[1] = time + 2.0*x[0]*x[0] - 2.0*x[0]*x[1]; 52 return 0; 53 } 54 55 PetscErrorCode mms1_p_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx) 56 { 57 *p = x[0] + x[1] - 1.0; 58 return 0; 59 } 60 61 /* MMS 2*/ 62 63 static PetscErrorCode mms2_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 64 { 65 u[0] = PetscSinReal(time + x[0])*PetscSinReal(time + x[1]); 66 u[1] = PetscCosReal(time + x[0])*PetscCosReal(time + x[1]); 67 return 0; 68 } 69 70 static PetscErrorCode mms2_p_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx) 71 { 72 *p = PetscSinReal(time + x[0] - x[1]); 73 return 0; 74 } 75 76 static void f0_mms1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 77 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 78 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 79 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 80 { 81 const PetscReal Re = REYN; 82 const PetscInt Ncomp = dim; 83 PetscInt c, d; 84 85 for (c = 0; c < Ncomp; ++c) { 86 for (d = 0; d < dim; ++d) { 87 f0[c] += u[d] * u_x[c*dim+d]; 88 } 89 } 90 f0[0] += u_t[0]; 91 f0[1] += u_t[1]; 92 93 f0[0] += -2.0*t*(x[0] + x[1]) + 2.0*x[0]*x[1]*x[1] - 4.0*x[0]*x[0]*x[1] - 2.0*x[0]*x[0]*x[0] + 4.0/Re - 1.0; 94 f0[1] += -2.0*t*x[0] + 2.0*x[1]*x[1]*x[1] - 4.0*x[0]*x[1]*x[1] - 2.0*x[0]*x[0]*x[1] + 4.0/Re - 1.0; 95 } 96 97 static void f0_mms2_u(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 const PetscReal Re = REYN; 103 const PetscInt Ncomp = dim; 104 PetscInt c, d; 105 106 for (c = 0; c < Ncomp; ++c) { 107 for (d = 0; d < dim; ++d) { 108 f0[c] += u[d] * u_x[c*dim+d]; 109 } 110 } 111 f0[0] += u_t[0]; 112 f0[1] += u_t[1]; 113 114 f0[0] -= ( Re*((1.0L/2.0L)*PetscSinReal(2*t + 2*x[0]) + PetscSinReal(2*t + x[0] + x[1]) + PetscCosReal(t + x[0] - x[1])) + 2.0*PetscSinReal(t + x[0])*PetscSinReal(t + x[1]))/Re; 115 f0[1] -= (-Re*((1.0L/2.0L)*PetscSinReal(2*t + 2*x[1]) + PetscSinReal(2*t + x[0] + x[1]) + PetscCosReal(t + x[0] - x[1])) + 2.0*PetscCosReal(t + x[0])*PetscCosReal(t + x[1]))/Re; 116 } 117 118 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 119 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 120 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 121 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 122 { 123 const PetscReal Re = REYN; 124 const PetscInt Ncomp = dim; 125 PetscInt comp, d; 126 127 for (comp = 0; comp < Ncomp; ++comp) { 128 for (d = 0; d < dim; ++d) { 129 f1[comp*dim+d] = 1.0/Re * u_x[comp*dim+d]; 130 } 131 f1[comp*dim+comp] -= u[Ncomp]; 132 } 133 } 134 135 static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, 136 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 137 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 138 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 139 { 140 PetscInt d; 141 for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d*dim+d]; 142 } 143 144 static void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, 145 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 146 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 147 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 148 { 149 PetscInt d; 150 for (d = 0; d < dim; ++d) f1[d] = 0.0; 151 } 152 153 /* 154 (psi_i, u_j grad_j u_i) ==> (\psi_i, \phi_j grad_j u_i) 155 */ 156 static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 157 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 158 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 159 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 160 { 161 PetscInt NcI = dim, NcJ = dim; 162 PetscInt fc, gc; 163 PetscInt d; 164 165 for (d = 0; d < dim; ++d) { 166 g0[d*dim+d] = u_tShift; 167 } 168 169 for (fc = 0; fc < NcI; ++fc) { 170 for (gc = 0; gc < NcJ; ++gc) { 171 g0[fc*NcJ+gc] += u_x[fc*NcJ+gc]; 172 } 173 } 174 } 175 176 /* 177 (psi_i, u_j grad_j u_i) ==> (\psi_i, \u_j grad_j \phi_i) 178 */ 179 static void g1_uu(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, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 183 { 184 PetscInt NcI = dim; 185 PetscInt NcJ = dim; 186 PetscInt fc, gc, dg; 187 for (fc = 0; fc < NcI; ++fc) { 188 for (gc = 0; gc < NcJ; ++gc) { 189 for (dg = 0; dg < dim; ++dg) { 190 /* kronecker delta */ 191 if (fc == gc) { 192 g1[(fc*NcJ+gc)*dim+dg] += u[dg]; 193 } 194 } 195 } 196 } 197 } 198 199 /* < q, \nabla\cdot u > 200 NcompI = 1, NcompJ = dim */ 201 static void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 202 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 203 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 204 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 205 { 206 PetscInt d; 207 for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */ 208 } 209 210 /* -< \nabla\cdot v, p > 211 NcompI = dim, NcompJ = 1 */ 212 static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, 213 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 214 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 215 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 216 { 217 PetscInt d; 218 for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */ 219 } 220 221 /* < \nabla v, \nabla u + {\nabla u}^T > 222 This just gives \nabla u, give the perdiagonal for the transpose */ 223 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 224 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 225 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 226 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 227 { 228 const PetscReal Re = REYN; 229 const PetscInt Ncomp = dim; 230 PetscInt compI, d; 231 232 for (compI = 0; compI < Ncomp; ++compI) { 233 for (d = 0; d < dim; ++d) { 234 g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0/Re; 235 } 236 } 237 } 238 239 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 240 { 241 PetscErrorCode ierr; 242 243 PetscFunctionBeginUser; 244 options->dim = 2; 245 options->simplex = PETSC_TRUE; 246 options->mms = 1; 247 248 ierr = PetscOptionsBegin(comm, "", "Navier-Stokes Equation Options", "DMPLEX");CHKERRQ(ierr); 249 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex46.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 250 ierr = PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex46.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 251 ierr = PetscOptionsInt("-mms", "The manufactured solution to use", "ex46.c", options->mms, &options->mms, NULL);CHKERRQ(ierr); 252 ierr = PetscOptionsEnd();CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) 257 { 258 DM plex; 259 DMLabel label; 260 PetscErrorCode ierr; 261 262 PetscFunctionBeginUser; 263 ierr = DMCreateLabel(dm, name);CHKERRQ(ierr); 264 ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr); 265 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 266 ierr = DMPlexMarkBoundaryFaces(plex, 1, label);CHKERRQ(ierr); 267 ierr = DMDestroy(&plex);CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 271 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) 272 { 273 DM pdm = NULL; 274 const PetscInt dim = ctx->dim; 275 PetscBool hasLabel; 276 PetscErrorCode ierr; 277 278 PetscFunctionBeginUser; 279 ierr = DMPlexCreateBoxMesh(comm, dim, ctx->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 280 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 281 /* If no boundary marker exists, mark the whole boundary */ 282 ierr = DMHasLabel(*dm, "marker", &hasLabel);CHKERRQ(ierr); 283 if (!hasLabel) {ierr = CreateBCLabel(*dm, "marker");CHKERRQ(ierr);} 284 /* Distribute mesh over processes */ 285 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 286 if (pdm) { 287 ierr = DMDestroy(dm);CHKERRQ(ierr); 288 *dm = pdm; 289 } 290 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 291 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) 296 { 297 PetscDS prob; 298 const PetscInt id = 1; 299 PetscErrorCode ierr; 300 301 PetscFunctionBeginUser; 302 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 303 switch (ctx->mms) { 304 case 1: 305 ierr = PetscDSSetResidual(prob, 0, f0_mms1_u, f1_u);CHKERRQ(ierr);break; 306 case 2: 307 ierr = PetscDSSetResidual(prob, 0, f0_mms2_u, f1_u);CHKERRQ(ierr);break; 308 } 309 ierr = PetscDSSetResidual(prob, 1, f0_p, f1_p);CHKERRQ(ierr); 310 ierr = PetscDSSetJacobian(prob, 0, 0, g0_uu, g1_uu, NULL, g3_uu);CHKERRQ(ierr); 311 ierr = PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_up, NULL);CHKERRQ(ierr); 312 ierr = PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu, NULL, NULL);CHKERRQ(ierr); 313 switch (ctx->dim) { 314 case 2: 315 switch (ctx->mms) { 316 case 1: 317 ctx->exactFuncs[0] = mms1_u_2d; 318 ctx->exactFuncs[1] = mms1_p_2d; 319 break; 320 case 2: 321 ctx->exactFuncs[0] = mms2_u_2d; 322 ctx->exactFuncs[1] = mms2_p_2d; 323 break; 324 default: 325 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid MMS %D", ctx->mms); 326 } 327 break; 328 default: 329 SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %D", ctx->dim); 330 } 331 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) ctx->exactFuncs[0], NULL, 1, &id, ctx);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) 336 { 337 DM cdm = dm; 338 const PetscInt dim = ctx->dim; 339 PetscFE fe[2]; 340 MPI_Comm comm; 341 PetscErrorCode ierr; 342 343 PetscFunctionBeginUser; 344 /* Create finite element */ 345 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 346 ierr = PetscFECreateDefault(comm, dim, dim, ctx->simplex, "vel_", PETSC_DEFAULT, &fe[0]);CHKERRQ(ierr); 347 ierr = PetscObjectSetName((PetscObject) fe[0], "velocity");CHKERRQ(ierr); 348 ierr = PetscFECreateDefault(comm, dim, 1, ctx->simplex, "pres_", PETSC_DEFAULT, &fe[1]);CHKERRQ(ierr); 349 ierr = PetscFECopyQuadrature(fe[0], fe[1]);CHKERRQ(ierr); 350 ierr = PetscObjectSetName((PetscObject) fe[1], "pressure");CHKERRQ(ierr); 351 /* Set discretization and boundary conditions for each mesh */ 352 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe[0]);CHKERRQ(ierr); 353 ierr = DMSetField(dm, 1, NULL, (PetscObject) fe[1]);CHKERRQ(ierr); 354 ierr = DMCreateDS(dm);CHKERRQ(ierr); 355 ierr = SetupProblem(dm, ctx);CHKERRQ(ierr); 356 while (cdm) { 357 PetscObject pressure; 358 MatNullSpace nsp; 359 PetscBool hasLabel; 360 361 ierr = DMGetField(cdm, 1, NULL, &pressure);CHKERRQ(ierr); 362 ierr = MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nsp);CHKERRQ(ierr); 363 ierr = PetscObjectCompose(pressure, "nullspace", (PetscObject) nsp);CHKERRQ(ierr); 364 ierr = MatNullSpaceDestroy(&nsp);CHKERRQ(ierr); 365 366 ierr = DMHasLabel(cdm, "marker", &hasLabel);CHKERRQ(ierr); 367 if (!hasLabel) {ierr = CreateBCLabel(cdm, "marker");CHKERRQ(ierr);} 368 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 369 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 370 } 371 ierr = PetscFEDestroy(&fe[0]);CHKERRQ(ierr); 372 ierr = PetscFEDestroy(&fe[1]);CHKERRQ(ierr); 373 PetscFunctionReturn(0); 374 } 375 376 static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) 377 { 378 AppCtx *user = (AppCtx *) ctx; 379 DM dm; 380 PetscReal ferrors[2]; 381 PetscErrorCode ierr; 382 383 PetscFunctionBeginUser; 384 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 385 ierr = DMComputeL2FieldDiff(dm, crtime, user->exactFuncs, NULL, u, ferrors);CHKERRQ(ierr); 386 ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g]\n", (int) step, (double) crtime, (double) ferrors[0], (double) ferrors[1]);CHKERRQ(ierr); 387 PetscFunctionReturn(0); 388 } 389 390 int main(int argc, char **argv) 391 { 392 AppCtx ctx; 393 DM dm; 394 TS ts; 395 Vec u, r; 396 PetscErrorCode ierr; 397 398 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 399 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 400 ierr = CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);CHKERRQ(ierr); 401 ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr); 402 ierr = PetscMalloc1(2, &ctx.exactFuncs);CHKERRQ(ierr); 403 ierr = SetupDiscretization(dm, &ctx);CHKERRQ(ierr); 404 ierr = DMPlexCreateClosureIndex(dm, NULL);CHKERRQ(ierr); 405 406 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 407 ierr = VecDuplicate(u, &r);CHKERRQ(ierr); 408 409 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr); 410 ierr = TSMonitorSet(ts, MonitorError, &ctx, NULL);CHKERRQ(ierr); 411 ierr = TSSetDM(ts, dm);CHKERRQ(ierr); 412 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr); 413 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr); 414 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr); 415 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 416 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 417 418 ierr = DMProjectFunction(dm, 0.0, ctx.exactFuncs, NULL, INSERT_ALL_VALUES, u);CHKERRQ(ierr); 419 ierr = TSSolve(ts, u);CHKERRQ(ierr); 420 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr); 421 422 ierr = VecDestroy(&u);CHKERRQ(ierr); 423 ierr = VecDestroy(&r);CHKERRQ(ierr); 424 ierr = TSDestroy(&ts);CHKERRQ(ierr); 425 ierr = DMDestroy(&dm);CHKERRQ(ierr); 426 ierr = PetscFree(ctx.exactFuncs);CHKERRQ(ierr); 427 ierr = PetscFinalize(); 428 return ierr; 429 } 430 431 /*TEST 432 433 # Full solves 434 test: 435 suffix: 2d_p2p1_r1 436 requires: !single triangle 437 filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" 438 args: -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor 439 test: 440 suffix: 2d_p2p1_r2 441 requires: !single triangle 442 filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" 443 args: -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor 444 test: 445 suffix: 2d_q2q1_r1 446 requires: !single 447 filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" -e "s~ 0\]~ 0.0\]~g" 448 args: -simplex 0 -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor 449 test: 450 suffix: 2d_q2q1_r2 451 requires: !single 452 filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" 453 args: -simplex 0 -dm_refine 2 -vel_petscspace_degree 2 -pres_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -fieldsplit_velocity_pc_type lu -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor 454 455 TEST*/ 456