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