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