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