1 static char help[] = "Pure advection with finite elements.\n\ 2 We solve the hyperbolic problem in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n"; 4 5 /* 6 The continuity equation (https://en.wikipedia.org/wiki/Continuity_equation) for advection 7 (https://en.wikipedia.org/wiki/Advection) of a conserved scalar quantity phi, with source q, 8 9 phi_t + div (phi u) = q 10 11 if used with a solenoidal velocity field u (div u = 0) is given by 12 13 phi_t + u . grad phi = q 14 15 For a vector quantity a, we likewise have 16 17 a_t + u . grad a = q 18 */ 19 20 /* 21 r1: 8 SOR 22 r2: 1128 SOR 23 r3: > 10000 SOR 24 25 SOR is completely unreliable as a smoother, use Jacobi 26 r1: 8 MG 27 r2: 28 */ 29 30 #include <petscdmplex.h> 31 #include <petscts.h> 32 #include <petscds.h> 33 34 typedef enum {PRIMITIVE, INT_BY_PARTS} WeakFormType; 35 36 typedef struct { 37 /* Domain and mesh definition */ 38 PetscInt dim; /* The topological mesh dimension */ 39 PetscBool simplex; /* Use simplices or tensor product cells */ 40 /* Problem definition */ 41 WeakFormType formType; 42 PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 43 } AppCtx; 44 45 /* MMS1: 46 47 2D: 48 u = <1, 1> 49 phi = x + y - 2t 50 51 phi_t + u . grad phi = -2 + <1, 1> . <1, 1> = 0 52 53 3D: 54 u = <1, 1, 1> 55 phi = x + y + z - 3t 56 57 phi_t + u . grad phi = -3 + <1, 1, 1> . <1, 1, 1> = 0 58 */ 59 60 static PetscErrorCode analytic_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 61 { 62 PetscInt d; 63 64 *u = -dim*time; 65 for (d = 0; d < dim; ++d) *u += x[d]; 66 return 0; 67 } 68 69 static PetscErrorCode velocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 70 { 71 PetscInt d; 72 for (d = 0; d < dim; ++d) u[d] = 1.0; 73 return 0; 74 } 75 76 /* <psi, phi_t> + <psi, u . grad phi> */ 77 static void f0_prim_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 78 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 79 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 80 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 81 { 82 PetscInt d; 83 84 f0[0] = u_t[0]; 85 for (d = 0; d < dim; ++d) f0[0] += a[d] * u_x[d]; 86 } 87 88 /* <psi, phi_t> */ 89 static void f0_ibp_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 90 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 91 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 92 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 93 { 94 f0[0] = u_t[0]; 95 } 96 97 /* <grad psi, u phi> */ 98 static void f1_ibp_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 99 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 100 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 101 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 102 { 103 PetscInt d; 104 for (d = 0; d < dim; ++d) f1[d] = a[d]*u[0]; 105 } 106 107 /* <psi, phi_t> */ 108 static void g0_prim_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 109 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 110 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 111 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 112 { 113 g0[0] = u_tShift*1.0; 114 } 115 116 /* <psi, u . grad phi> */ 117 static void g1_prim_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 118 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 119 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 120 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 121 { 122 PetscInt d; 123 for (d = 0; d < dim; ++d) g1[d] = a[d]; 124 } 125 126 /* <grad psi, u phi> */ 127 static void g2_ibp_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 128 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 129 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 130 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 131 { 132 PetscInt d; 133 for (d = 0; d < dim; ++d) g2[d] = a[d]; 134 } 135 136 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 137 { 138 const char *formTypes[2] = {"primitive", "int_by_parts"}; 139 PetscInt form; 140 PetscErrorCode ierr; 141 142 PetscFunctionBeginUser; 143 options->dim = 2; 144 options->simplex = PETSC_TRUE; 145 options->formType = PRIMITIVE; 146 147 ierr = PetscOptionsBegin(comm, "", "Advection Equation Options", "DMPLEX");CHKERRQ(ierr); 148 ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex47.c", options->dim, &options->dim, NULL);CHKERRQ(ierr); 149 ierr = PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex47.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 150 form = options->formType; 151 ierr = PetscOptionsEList("-form_type", "The weak form type", "ex47.c", formTypes, 2, formTypes[options->formType], &form, NULL);CHKERRQ(ierr); 152 options->formType = (WeakFormType) form; 153 ierr = PetscOptionsEnd();CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) 158 { 159 DM plex; 160 DMLabel label; 161 PetscErrorCode ierr; 162 163 PetscFunctionBeginUser; 164 ierr = DMCreateLabel(dm, name);CHKERRQ(ierr); 165 ierr = DMGetLabel(dm, name, &label);CHKERRQ(ierr); 166 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 167 ierr = DMPlexMarkBoundaryFaces(dm, 1, label);CHKERRQ(ierr); 168 ierr = DMDestroy(&plex);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) 173 { 174 DM pdm = NULL; 175 const PetscInt dim = ctx->dim; 176 PetscBool hasLabel; 177 PetscErrorCode ierr; 178 179 PetscFunctionBeginUser; 180 ierr = DMPlexCreateBoxMesh(comm, dim, ctx->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 181 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 182 /* If no boundary marker exists, mark the whole boundary */ 183 ierr = DMHasLabel(*dm, "marker", &hasLabel);CHKERRQ(ierr); 184 if (!hasLabel) {ierr = CreateBCLabel(*dm, "marker");CHKERRQ(ierr);} 185 /* Distribute mesh over processes */ 186 ierr = DMPlexDistribute(*dm, 0, NULL, &pdm);CHKERRQ(ierr); 187 if (pdm) { 188 ierr = DMDestroy(dm);CHKERRQ(ierr); 189 *dm = pdm; 190 } 191 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 192 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) 197 { 198 PetscDS prob; 199 const PetscInt id = 1; 200 PetscErrorCode ierr; 201 202 PetscFunctionBeginUser; 203 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 204 switch (ctx->formType) { 205 case PRIMITIVE: 206 ierr = PetscDSSetResidual(prob, 0, f0_prim_phi, NULL);CHKERRQ(ierr); 207 ierr = PetscDSSetJacobian(prob, 0, 0, g0_prim_phi, g1_prim_phi, NULL, NULL);CHKERRQ(ierr); 208 break; 209 case INT_BY_PARTS: 210 ierr = PetscDSSetResidual(prob, 0, f0_ibp_phi, f1_ibp_phi);CHKERRQ(ierr); 211 ierr = PetscDSSetJacobian(prob, 0, 0, g0_prim_phi, NULL, g2_ibp_phi, NULL);CHKERRQ(ierr); 212 break; 213 } 214 ctx->exactFuncs[0] = analytic_phi; 215 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) ctx->exactFuncs[0], NULL, 1, &id, ctx);CHKERRQ(ierr); 216 PetscFunctionReturn(0); 217 } 218 219 static PetscErrorCode SetupVelocity(DM dm, DM dmAux, AppCtx *user) 220 { 221 PetscErrorCode (*funcs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {velocity}; 222 Vec v; 223 PetscErrorCode ierr; 224 225 PetscFunctionBeginUser; 226 ierr = DMCreateLocalVector(dmAux, &v);CHKERRQ(ierr); 227 ierr = DMProjectFunctionLocal(dmAux, 0.0, funcs, NULL, INSERT_ALL_VALUES, v);CHKERRQ(ierr); 228 ierr = PetscObjectCompose((PetscObject) dm, "A", (PetscObject) v);CHKERRQ(ierr); 229 ierr = VecDestroy(&v);CHKERRQ(ierr); 230 PetscFunctionReturn(0); 231 } 232 233 static PetscErrorCode SetupAuxDM(DM dm, PetscFE feAux, AppCtx *user) 234 { 235 DM dmAux, coordDM; 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */ 240 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 241 if (!feAux) PetscFunctionReturn(0); 242 ierr = DMClone(dm, &dmAux);CHKERRQ(ierr); 243 ierr = PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);CHKERRQ(ierr); 244 ierr = DMSetCoordinateDM(dmAux, coordDM);CHKERRQ(ierr); 245 ierr = DMSetField(dmAux, 0, NULL, (PetscObject) feAux);CHKERRQ(ierr); 246 ierr = DMCreateDS(dmAux);CHKERRQ(ierr); 247 ierr = SetupVelocity(dm, dmAux, user);CHKERRQ(ierr); 248 ierr = DMDestroy(&dmAux);CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 252 static PetscErrorCode SetupDiscretization(DM dm, AppCtx* ctx) 253 { 254 DM cdm = dm; 255 const PetscInt dim = ctx->dim; 256 PetscFE fe, feAux; 257 MPI_Comm comm; 258 PetscErrorCode ierr; 259 260 PetscFunctionBeginUser; 261 /* Create finite element */ 262 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 263 ierr = PetscFECreateDefault(comm, dim, 1, ctx->simplex, "phi_", -1, &fe);CHKERRQ(ierr); 264 ierr = PetscObjectSetName((PetscObject) fe, "phi");CHKERRQ(ierr); 265 /* Create velocity */ 266 ierr = PetscFECreateDefault(comm, dim, dim, ctx->simplex, "vel_", -1, &feAux);CHKERRQ(ierr); 267 ierr = PetscFECopyQuadrature(fe, feAux);CHKERRQ(ierr); 268 /* Set discretization and boundary conditions for each mesh */ 269 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 270 ierr = DMCreateDS(dm);CHKERRQ(ierr); 271 ierr = SetupProblem(dm, ctx);CHKERRQ(ierr); 272 while (cdm) { 273 PetscBool hasLabel; 274 275 ierr = SetupAuxDM(cdm, feAux, ctx);CHKERRQ(ierr); 276 ierr = DMHasLabel(cdm, "marker", &hasLabel);CHKERRQ(ierr); 277 if (!hasLabel) {ierr = CreateBCLabel(cdm, "marker");CHKERRQ(ierr);} 278 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 279 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 280 } 281 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 282 ierr = PetscFEDestroy(&feAux);CHKERRQ(ierr); 283 PetscFunctionReturn(0); 284 } 285 286 static PetscErrorCode MonitorError(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx) 287 { 288 AppCtx *user = (AppCtx *) ctx; 289 DM dm; 290 Vec u, r, error; 291 PetscReal time = 0.5, res; 292 PetscErrorCode ierr; 293 294 PetscFunctionBeginUser; 295 ierr = KSPGetDM(ksp, &dm);CHKERRQ(ierr); 296 ierr = DMSetOutputSequenceNumber(dm, it, time);CHKERRQ(ierr); 297 /* Calculate residual */ 298 ierr = KSPBuildResidual(ksp, NULL, NULL, &r);CHKERRQ(ierr); 299 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 300 ierr = DMSetOutputSequenceNumber(dm, it, res);CHKERRQ(ierr); 301 ierr = PetscObjectSetName((PetscObject) r, "residual");CHKERRQ(ierr); 302 ierr = VecViewFromOptions(r, NULL, "-res_vec_view");CHKERRQ(ierr); 303 ierr = VecDestroy(&r);CHKERRQ(ierr); 304 /* Calculate error */ 305 ierr = KSPBuildSolution(ksp, NULL, &u);CHKERRQ(ierr); 306 ierr = DMGetGlobalVector(dm, &error);CHKERRQ(ierr); 307 ierr = DMProjectFunction(dm, time, user->exactFuncs, NULL, INSERT_ALL_VALUES, error);CHKERRQ(ierr); 308 ierr = VecAXPY(error, -1.0, u);CHKERRQ(ierr); 309 ierr = PetscObjectSetName((PetscObject) error, "error");CHKERRQ(ierr); 310 ierr = VecViewFromOptions(error, NULL, "-err_vec_view");CHKERRQ(ierr); 311 ierr = DMRestoreGlobalVector(dm, &error);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 static PetscErrorCode MyTSMonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) 316 { 317 AppCtx *user = (AppCtx *) ctx; 318 DM dm; 319 PetscReal error; 320 PetscErrorCode ierr; 321 322 PetscFunctionBeginUser; 323 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 324 ierr = DMComputeL2Diff(dm, crtime, user->exactFuncs, NULL, u, &error);CHKERRQ(ierr); 325 ierr = PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: %2.5g\n", (int) step, (double) crtime, (double) error);CHKERRQ(ierr); 326 PetscFunctionReturn(0); 327 } 328 329 int main(int argc, char **argv) 330 { 331 AppCtx ctx; 332 DM dm; 333 TS ts; 334 Vec u, r; 335 PetscReal t = 0.0; 336 PetscErrorCode ierr; 337 338 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 339 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 340 ierr = CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);CHKERRQ(ierr); 341 ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr); 342 ierr = PetscMalloc1(1, &ctx.exactFuncs);CHKERRQ(ierr); 343 ierr = SetupDiscretization(dm, &ctx);CHKERRQ(ierr); 344 345 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 346 ierr = PetscObjectSetName((PetscObject) u, "phi");CHKERRQ(ierr); 347 ierr = VecDuplicate(u, &r);CHKERRQ(ierr); 348 349 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr); 350 ierr = TSMonitorSet(ts, MyTSMonitorError, &ctx, NULL);CHKERRQ(ierr); 351 ierr = TSSetDM(ts, dm);CHKERRQ(ierr); 352 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr); 353 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr); 354 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr); 355 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 356 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 357 358 ierr = DMProjectFunction(dm, t, ctx.exactFuncs, NULL, INSERT_ALL_VALUES, u);CHKERRQ(ierr); 359 { 360 SNES snes; 361 KSP ksp; 362 363 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 364 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 365 ierr = KSPMonitorSet(ksp, MonitorError, &ctx, NULL);CHKERRQ(ierr); 366 } 367 ierr = TSSolve(ts, u);CHKERRQ(ierr); 368 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr); 369 370 ierr = VecDestroy(&u);CHKERRQ(ierr); 371 ierr = VecDestroy(&r);CHKERRQ(ierr); 372 ierr = TSDestroy(&ts);CHKERRQ(ierr); 373 ierr = DMDestroy(&dm);CHKERRQ(ierr); 374 ierr = PetscFree(ctx.exactFuncs);CHKERRQ(ierr); 375 ierr = PetscFinalize(); 376 return ierr; 377 } 378 379 /*TEST 380 381 # Full solves 382 test: 383 suffix: 2d_p1p1_r1 384 requires: triangle 385 args: -dm_refine 1 -phi_petscspace_degree 1 -vel_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -snes_monitor_short -snes_converged_reason -ts_monitor 386 387 test: 388 suffix: 2d_p1p1_sor_r1 389 requires: triangle !single 390 args: -dm_refine 1 -phi_petscspace_degree 1 -vel_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ksp_rtol 1.0e-9 -pc_type sor -snes_monitor_short -snes_converged_reason -ksp_monitor_short -ts_monitor 391 392 test: 393 suffix: 2d_p1p1_mg_r1 394 requires: triangle !single 395 args: -dm_refine_hierarchy 1 -phi_petscspace_degree 1 -vel_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ksp_type fgmres -ksp_rtol 1.0e-9 -pc_type mg -pc_mg_levels 2 -snes_monitor_short -snes_converged_reason -snes_view -ksp_monitor_true_residual -ts_monitor 396 397 TEST*/ 398