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 WeakFormType formType; 38 } AppCtx; 39 40 /* MMS1: 41 42 2D: 43 u = <1, 1> 44 phi = x + y - 2t 45 46 phi_t + u . grad phi = -2 + <1, 1> . <1, 1> = 0 47 48 3D: 49 u = <1, 1, 1> 50 phi = x + y + z - 3t 51 52 phi_t + u . grad phi = -3 + <1, 1, 1> . <1, 1, 1> = 0 53 */ 54 55 static PetscErrorCode analytic_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 56 { 57 PetscInt d; 58 59 *u = -dim*time; 60 for (d = 0; d < dim; ++d) *u += x[d]; 61 return 0; 62 } 63 64 static PetscErrorCode velocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 65 { 66 PetscInt d; 67 for (d = 0; d < dim; ++d) u[d] = 1.0; 68 return 0; 69 } 70 71 /* <psi, phi_t> + <psi, u . grad phi> */ 72 static void f0_prim_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 73 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 74 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 75 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 76 { 77 PetscInt d; 78 79 f0[0] = u_t[0]; 80 for (d = 0; d < dim; ++d) f0[0] += a[d] * u_x[d]; 81 } 82 83 /* <psi, phi_t> */ 84 static void f0_ibp_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 85 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 86 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 87 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 88 { 89 f0[0] = u_t[0]; 90 } 91 92 /* <grad psi, u phi> */ 93 static void f1_ibp_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 94 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 95 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 96 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 97 { 98 PetscInt d; 99 for (d = 0; d < dim; ++d) f1[d] = a[d]*u[0]; 100 } 101 102 /* <psi, phi_t> */ 103 static void g0_prim_phi(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, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 107 { 108 g0[0] = u_tShift*1.0; 109 } 110 111 /* <psi, u . grad phi> */ 112 static void g1_prim_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 113 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 114 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 115 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 116 { 117 PetscInt d; 118 for (d = 0; d < dim; ++d) g1[d] = a[d]; 119 } 120 121 /* <grad psi, u phi> */ 122 static void g2_ibp_phi(PetscInt dim, PetscInt Nf, PetscInt NfAux, 123 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 124 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 125 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]) 126 { 127 PetscInt d; 128 for (d = 0; d < dim; ++d) g2[d] = a[d]; 129 } 130 131 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 132 { 133 const char *formTypes[2] = {"primitive", "int_by_parts"}; 134 PetscInt form; 135 PetscErrorCode ierr; 136 137 PetscFunctionBeginUser; 138 options->formType = PRIMITIVE; 139 140 ierr = PetscOptionsBegin(comm, "", "Advection Equation Options", "DMPLEX");CHKERRQ(ierr); 141 form = options->formType; 142 ierr = PetscOptionsEList("-form_type", "The weak form type", "ex47.c", formTypes, 2, formTypes[options->formType], &form, NULL);CHKERRQ(ierr); 143 options->formType = (WeakFormType) form; 144 ierr = PetscOptionsEnd();CHKERRQ(ierr); 145 PetscFunctionReturn(0); 146 } 147 148 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) 149 { 150 PetscErrorCode ierr; 151 152 PetscFunctionBeginUser; 153 ierr = DMCreate(comm, dm);CHKERRQ(ierr); 154 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr); 155 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 156 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) 161 { 162 PetscDS ds; 163 DMLabel label; 164 const PetscInt id = 1; 165 PetscErrorCode ierr; 166 167 PetscFunctionBeginUser; 168 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 169 switch (ctx->formType) { 170 case PRIMITIVE: 171 ierr = PetscDSSetResidual(ds, 0, f0_prim_phi, NULL);CHKERRQ(ierr); 172 ierr = PetscDSSetJacobian(ds, 0, 0, g0_prim_phi, g1_prim_phi, NULL, NULL);CHKERRQ(ierr); 173 break; 174 case INT_BY_PARTS: 175 ierr = PetscDSSetResidual(ds, 0, f0_ibp_phi, f1_ibp_phi);CHKERRQ(ierr); 176 ierr = PetscDSSetJacobian(ds, 0, 0, g0_prim_phi, NULL, g2_ibp_phi, NULL);CHKERRQ(ierr); 177 break; 178 } 179 ierr = PetscDSSetExactSolution(ds, 0, analytic_phi, ctx);CHKERRQ(ierr); 180 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 181 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) analytic_phi, NULL, ctx, NULL);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 static PetscErrorCode SetupVelocity(DM dm, DM dmAux, AppCtx *user) 186 { 187 PetscSimplePointFunc funcs[1] = {velocity}; 188 Vec v; 189 PetscErrorCode ierr; 190 191 PetscFunctionBeginUser; 192 ierr = DMCreateLocalVector(dmAux, &v);CHKERRQ(ierr); 193 ierr = DMProjectFunctionLocal(dmAux, 0.0, funcs, NULL, INSERT_ALL_VALUES, v);CHKERRQ(ierr); 194 ierr = DMSetAuxiliaryVec(dm, NULL, 0, v);CHKERRQ(ierr); 195 ierr = VecDestroy(&v);CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 199 static PetscErrorCode SetupAuxDM(DM dm, PetscFE feAux, AppCtx *user) 200 { 201 DM dmAux, coordDM; 202 PetscErrorCode ierr; 203 204 PetscFunctionBegin; 205 /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */ 206 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 207 if (!feAux) PetscFunctionReturn(0); 208 ierr = DMClone(dm, &dmAux);CHKERRQ(ierr); 209 ierr = DMSetCoordinateDM(dmAux, coordDM);CHKERRQ(ierr); 210 ierr = DMSetField(dmAux, 0, NULL, (PetscObject) feAux);CHKERRQ(ierr); 211 ierr = DMCreateDS(dmAux);CHKERRQ(ierr); 212 ierr = SetupVelocity(dm, dmAux, user);CHKERRQ(ierr); 213 ierr = DMDestroy(&dmAux);CHKERRQ(ierr); 214 PetscFunctionReturn(0); 215 } 216 217 static PetscErrorCode SetupDiscretization(DM dm, AppCtx* ctx) 218 { 219 DM cdm = dm; 220 PetscFE fe, feAux; 221 MPI_Comm comm; 222 PetscInt dim; 223 PetscBool simplex; 224 PetscErrorCode ierr; 225 226 PetscFunctionBeginUser; 227 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 228 ierr = DMPlexIsSimplex(dm, &simplex);CHKERRQ(ierr); 229 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 230 ierr = PetscFECreateDefault(comm, dim, 1, simplex, "phi_", -1, &fe);CHKERRQ(ierr); 231 ierr = PetscObjectSetName((PetscObject) fe, "phi");CHKERRQ(ierr); 232 ierr = PetscFECreateDefault(comm, dim, dim, simplex, "vel_", -1, &feAux);CHKERRQ(ierr); 233 ierr = PetscFECopyQuadrature(fe, feAux);CHKERRQ(ierr); 234 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 235 ierr = DMCreateDS(dm);CHKERRQ(ierr); 236 ierr = SetupProblem(dm, ctx);CHKERRQ(ierr); 237 while (cdm) { 238 ierr = SetupAuxDM(cdm, feAux, ctx);CHKERRQ(ierr); 239 ierr = DMCopyDisc(dm, cdm);CHKERRQ(ierr); 240 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 241 } 242 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 243 ierr = PetscFEDestroy(&feAux);CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246 247 static PetscErrorCode MonitorError(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx) 248 { 249 DM dm; 250 PetscDS ds; 251 PetscSimplePointFunc func[1]; 252 void *ctxs[1]; 253 Vec u, r, error; 254 PetscReal time = 0.5, res; 255 PetscErrorCode ierr; 256 257 PetscFunctionBeginUser; 258 ierr = KSPGetDM(ksp, &dm);CHKERRQ(ierr); 259 ierr = DMSetOutputSequenceNumber(dm, it, time);CHKERRQ(ierr); 260 /* Calculate residual */ 261 ierr = KSPBuildResidual(ksp, NULL, NULL, &r);CHKERRQ(ierr); 262 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 263 ierr = DMSetOutputSequenceNumber(dm, it, res);CHKERRQ(ierr); 264 ierr = PetscObjectSetName((PetscObject) r, "residual");CHKERRQ(ierr); 265 ierr = VecViewFromOptions(r, NULL, "-res_vec_view");CHKERRQ(ierr); 266 ierr = VecDestroy(&r);CHKERRQ(ierr); 267 /* Calculate error */ 268 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 269 ierr = PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]);CHKERRQ(ierr); 270 ierr = KSPBuildSolution(ksp, NULL, &u);CHKERRQ(ierr); 271 ierr = DMGetGlobalVector(dm, &error);CHKERRQ(ierr); 272 ierr = DMProjectFunction(dm, time, func, ctxs, INSERT_ALL_VALUES, error);CHKERRQ(ierr); 273 ierr = VecAXPY(error, -1.0, u);CHKERRQ(ierr); 274 ierr = PetscObjectSetName((PetscObject) error, "error");CHKERRQ(ierr); 275 ierr = VecViewFromOptions(error, NULL, "-err_vec_view");CHKERRQ(ierr); 276 ierr = DMRestoreGlobalVector(dm, &error);CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 280 static PetscErrorCode MyTSMonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) 281 { 282 DM dm; 283 PetscDS ds; 284 PetscSimplePointFunc func[1]; 285 void *ctxs[1]; 286 PetscReal error; 287 PetscErrorCode ierr; 288 289 PetscFunctionBeginUser; 290 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 291 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 292 ierr = PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]);CHKERRQ(ierr); 293 ierr = DMComputeL2Diff(dm, crtime, func, ctxs, u, &error);CHKERRQ(ierr); 294 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); 295 PetscFunctionReturn(0); 296 } 297 298 int main(int argc, char **argv) 299 { 300 AppCtx ctx; 301 DM dm; 302 TS ts; 303 Vec u, r; 304 PetscReal t = 0.0; 305 PetscErrorCode ierr; 306 307 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 308 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 309 ierr = CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);CHKERRQ(ierr); 310 ierr = DMSetApplicationContext(dm, &ctx);CHKERRQ(ierr); 311 ierr = SetupDiscretization(dm, &ctx);CHKERRQ(ierr); 312 313 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 314 ierr = PetscObjectSetName((PetscObject) u, "phi");CHKERRQ(ierr); 315 ierr = VecDuplicate(u, &r);CHKERRQ(ierr); 316 317 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr); 318 ierr = TSMonitorSet(ts, MyTSMonitorError, &ctx, NULL);CHKERRQ(ierr); 319 ierr = TSSetDM(ts, dm);CHKERRQ(ierr); 320 ierr = DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);CHKERRQ(ierr); 321 ierr = DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);CHKERRQ(ierr); 322 ierr = DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);CHKERRQ(ierr); 323 ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 324 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 325 326 { 327 PetscDS ds; 328 PetscSimplePointFunc func[1]; 329 void *ctxs[1]; 330 331 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 332 ierr = PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]);CHKERRQ(ierr); 333 ierr = DMProjectFunction(dm, t, func, ctxs, INSERT_ALL_VALUES, u);CHKERRQ(ierr); 334 } 335 { 336 SNES snes; 337 KSP ksp; 338 339 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 340 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 341 ierr = KSPMonitorSet(ksp, MonitorError, &ctx, NULL);CHKERRQ(ierr); 342 } 343 ierr = TSSolve(ts, u);CHKERRQ(ierr); 344 ierr = VecViewFromOptions(u, NULL, "-sol_vec_view");CHKERRQ(ierr); 345 346 ierr = VecDestroy(&u);CHKERRQ(ierr); 347 ierr = VecDestroy(&r);CHKERRQ(ierr); 348 ierr = TSDestroy(&ts);CHKERRQ(ierr); 349 ierr = DMDestroy(&dm);CHKERRQ(ierr); 350 ierr = PetscFinalize(); 351 return ierr; 352 } 353 354 /*TEST 355 356 # Full solves 357 test: 358 suffix: 2d_p1p1_r1 359 requires: triangle 360 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 361 362 test: 363 suffix: 2d_p1p1_sor_r1 364 requires: triangle !single 365 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 366 367 test: 368 suffix: 2d_p1p1_mg_r1 369 requires: triangle !single 370 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 371 372 TEST*/ 373