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 PETSC_SUCCESS;
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 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: 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;
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: }
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: }
141: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
142: {
143: PetscDS ds;
144: DMLabel label;
145: const PetscInt id = 1;
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, (void (*)(void))analytic_phi, NULL, ctx, NULL));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode SetupVelocity(DM dm, DM dmAux, AppCtx *user)
166: {
167: PetscSimplePointFunc funcs[1] = {velocity};
168: Vec v;
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: }
178: static PetscErrorCode SetupAuxDM(DM dm, PetscFE feAux, AppCtx *user)
179: {
180: DM dmAux, coordDM;
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: }
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;
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: }
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;
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: }
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;
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: }
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;
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));
288: PetscCall(DMCreateGlobalVector(dm, &u));
289: PetscCall(PetscObjectSetName((PetscObject)u, "phi"));
290: PetscCall(VecDuplicate(u, &r));
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));
301: {
302: PetscDS ds;
303: PetscSimplePointFunc func[1];
304: void *ctxs[1];
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;
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"));
321: PetscCall(VecDestroy(&u));
322: PetscCall(VecDestroy(&r));
323: PetscCall(TSDestroy(&ts));
324: PetscCall(DMDestroy(&dm));
325: PetscCall(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*/