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
analytic_phi(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,PetscCtx ctx)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
velocity(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nf,PetscScalar * u,PetscCtx ctx)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> */
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[])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> */
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[])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> */
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[])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> */
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 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> */
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 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> */
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[])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
ProcessOptions(MPI_Comm comm,AppCtx * options)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
CreateMesh(MPI_Comm comm,DM * dm,AppCtx * ctx)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
SetupProblem(DM dm,AppCtx * ctx)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
SetupVelocity(DM dm,DM dmAux,AppCtx * user)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
SetupAuxDM(DM dm,PetscFE feAux,AppCtx * user)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
SetupDiscretization(DM dm,AppCtx * ctx)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
MonitorError(KSP ksp,PetscInt it,PetscReal rnorm,PetscCtx ctx)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
MyTSMonitorError(TS ts,PetscInt step,PetscReal crtime,Vec u,PetscCtx ctx)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
main(int argc,char ** argv)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