Actual source code: ex47.c
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";
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,
9: phi_t + div (phi u) = q
11: if used with a solenoidal velocity field u (div u = 0) is given by
13: phi_t + u . grad phi = q
15: For a vector quantity a, we likewise have
17: a_t + u . grad a = q
18: */
20: /*
21: r1: 8 SOR
22: r2: 1128 SOR
23: r3: > 10000 SOR
25: SOR is completely unreliable as a smoother, use Jacobi
26: r1: 8 MG
27: r2:
28: */
30: #include <petscdmplex.h>
31: #include <petscts.h>
32: #include <petscds.h>
34: typedef enum {
35: PRIMITIVE,
36: INT_BY_PARTS
37: } WeakFormType;
39: typedef struct {
40: WeakFormType formType;
41: } AppCtx;
43: /* MMS1:
45: 2D:
46: u = <1, 1>
47: phi = x + y - 2t
49: phi_t + u . grad phi = -2 + <1, 1> . <1, 1> = 0
51: 3D:
52: u = <1, 1, 1>
53: phi = x + y + z - 3t
55: phi_t + u . grad phi = -3 + <1, 1, 1> . <1, 1, 1> = 0
56: */
58: static PetscErrorCode analytic_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
59: {
60: PetscInt d;
62: *u = -dim * time;
63: for (d = 0; d < dim; ++d) *u += x[d];
64: return 0;
65: }
67: static PetscErrorCode velocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
68: {
69: PetscInt d;
70: for (d = 0; d < dim; ++d) u[d] = 1.0;
71: return 0;
72: }
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;
79: f0[0] = u_t[0];
80: for (d = 0; d < dim; ++d) f0[0] += a[d] * u_x[d];
81: }
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: }
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: }
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: }
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: }
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: }
116: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
117: {
118: const char *formTypes[2] = {"primitive", "int_by_parts"};
119: PetscInt form;
122: options->formType = PRIMITIVE;
123: PetscOptionsBegin(comm, "", "Advection Equation Options", "DMPLEX");
124: form = options->formType;
125: PetscOptionsEList("-form_type", "The weak form type", "ex47.c", formTypes, 2, formTypes[options->formType], &form, NULL);
126: options->formType = (WeakFormType)form;
127: PetscOptionsEnd();
128: return 0;
129: }
131: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
132: {
134: DMCreate(comm, dm);
135: DMSetType(*dm, DMPLEX);
136: DMSetFromOptions(*dm);
137: DMViewFromOptions(*dm, NULL, "-dm_view");
138: return 0;
139: }
141: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
142: {
143: PetscDS ds;
144: DMLabel label;
145: const PetscInt id = 1;
148: DMGetDS(dm, &ds);
149: switch (ctx->formType) {
150: case PRIMITIVE:
151: PetscDSSetResidual(ds, 0, f0_prim_phi, NULL);
152: PetscDSSetJacobian(ds, 0, 0, g0_prim_phi, g1_prim_phi, NULL, NULL);
153: break;
154: case INT_BY_PARTS:
155: PetscDSSetResidual(ds, 0, f0_ibp_phi, f1_ibp_phi);
156: PetscDSSetJacobian(ds, 0, 0, g0_prim_phi, NULL, g2_ibp_phi, NULL);
157: break;
158: }
159: PetscDSSetExactSolution(ds, 0, analytic_phi, ctx);
160: DMGetLabel(dm, "marker", &label);
161: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))analytic_phi, NULL, ctx, NULL);
162: return 0;
163: }
165: static PetscErrorCode SetupVelocity(DM dm, DM dmAux, AppCtx *user)
166: {
167: PetscSimplePointFunc funcs[1] = {velocity};
168: Vec v;
171: DMCreateLocalVector(dmAux, &v);
172: DMProjectFunctionLocal(dmAux, 0.0, funcs, NULL, INSERT_ALL_VALUES, v);
173: DMSetAuxiliaryVec(dm, NULL, 0, 0, v);
174: VecDestroy(&v);
175: return 0;
176: }
178: static PetscErrorCode SetupAuxDM(DM dm, PetscFE feAux, AppCtx *user)
179: {
180: DM dmAux, coordDM;
183: /* MUST call DMGetCoordinateDM() in order to get p4est setup if present */
184: DMGetCoordinateDM(dm, &coordDM);
185: if (!feAux) return 0;
186: DMClone(dm, &dmAux);
187: DMSetCoordinateDM(dmAux, coordDM);
188: DMSetField(dmAux, 0, NULL, (PetscObject)feAux);
189: DMCreateDS(dmAux);
190: SetupVelocity(dm, dmAux, user);
191: DMDestroy(&dmAux);
192: return 0;
193: }
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;
204: DMGetDimension(dm, &dim);
205: DMPlexIsSimplex(dm, &simplex);
206: PetscObjectGetComm((PetscObject)dm, &comm);
207: PetscFECreateDefault(comm, dim, 1, simplex, "phi_", -1, &fe);
208: PetscObjectSetName((PetscObject)fe, "phi");
209: PetscFECreateDefault(comm, dim, dim, simplex, "vel_", -1, &feAux);
210: PetscFECopyQuadrature(fe, feAux);
211: DMSetField(dm, 0, NULL, (PetscObject)fe);
212: DMCreateDS(dm);
213: SetupProblem(dm, ctx);
214: while (cdm) {
215: SetupAuxDM(cdm, feAux, ctx);
216: DMCopyDisc(dm, cdm);
217: DMGetCoarseDM(cdm, &cdm);
218: }
219: PetscFEDestroy(&fe);
220: PetscFEDestroy(&feAux);
221: return 0;
222: }
224: static PetscErrorCode MonitorError(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
225: {
226: DM dm;
227: PetscDS ds;
228: PetscSimplePointFunc func[1];
229: void *ctxs[1];
230: Vec u, r, error;
231: PetscReal time = 0.5, res;
234: KSPGetDM(ksp, &dm);
235: DMSetOutputSequenceNumber(dm, it, time);
236: /* Calculate residual */
237: KSPBuildResidual(ksp, NULL, NULL, &r);
238: VecNorm(r, NORM_2, &res);
239: DMSetOutputSequenceNumber(dm, it, res);
240: PetscObjectSetName((PetscObject)r, "residual");
241: VecViewFromOptions(r, NULL, "-res_vec_view");
242: VecDestroy(&r);
243: /* Calculate error */
244: DMGetDS(dm, &ds);
245: PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]);
246: KSPBuildSolution(ksp, NULL, &u);
247: DMGetGlobalVector(dm, &error);
248: DMProjectFunction(dm, time, func, ctxs, INSERT_ALL_VALUES, error);
249: VecAXPY(error, -1.0, u);
250: PetscObjectSetName((PetscObject)error, "error");
251: VecViewFromOptions(error, NULL, "-err_vec_view");
252: DMRestoreGlobalVector(dm, &error);
253: return 0;
254: }
256: static PetscErrorCode MyTSMonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
257: {
258: DM dm;
259: PetscDS ds;
260: PetscSimplePointFunc func[1];
261: void *ctxs[1];
262: PetscReal error;
265: TSGetDM(ts, &dm);
266: DMGetDS(dm, &ds);
267: PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]);
268: DMComputeL2Diff(dm, crtime, func, ctxs, u, &error);
269: PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: %2.5g\n", (int)step, (double)crtime, (double)error);
270: return 0;
271: }
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;
282: PetscInitialize(&argc, &argv, NULL, help);
283: ProcessOptions(PETSC_COMM_WORLD, &ctx);
284: CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);
285: DMSetApplicationContext(dm, &ctx);
286: SetupDiscretization(dm, &ctx);
288: DMCreateGlobalVector(dm, &u);
289: PetscObjectSetName((PetscObject)u, "phi");
290: VecDuplicate(u, &r);
292: TSCreate(PETSC_COMM_WORLD, &ts);
293: TSMonitorSet(ts, MyTSMonitorError, &ctx, NULL);
294: TSSetDM(ts, dm);
295: DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
296: DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
297: DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
298: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
299: TSSetFromOptions(ts);
301: {
302: PetscDS ds;
303: PetscSimplePointFunc func[1];
304: void *ctxs[1];
306: DMGetDS(dm, &ds);
307: PetscDSGetExactSolution(ds, 0, &func[0], &ctxs[0]);
308: DMProjectFunction(dm, t, func, ctxs, INSERT_ALL_VALUES, u);
309: }
310: {
311: SNES snes;
312: KSP ksp;
314: TSGetSNES(ts, &snes);
315: SNESGetKSP(snes, &ksp);
316: KSPMonitorSet(ksp, MonitorError, &ctx, NULL);
317: }
318: TSSolve(ts, u);
319: VecViewFromOptions(u, NULL, "-sol_vec_view");
321: VecDestroy(&u);
322: VecDestroy(&r);
323: TSDestroy(&ts);
324: DMDestroy(&dm);
325: PetscFinalize();
326: return 0;
327: }
329: /*TEST
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_dt 0.1 -pc_type lu -snes_monitor_short -snes_converged_reason -ts_monitor
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_dt 0.1 -ksp_rtol 1.0e-9 -pc_type sor -snes_monitor_short -snes_converged_reason -ksp_monitor_short -ts_monitor
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_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
347: TEST*/