Actual source code: ex71.c
1: static char help[] = "Poiseuille Flow in 2d and 3d channels with finite elements.\n\
2: We solve the Poiseuille flow problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
5: /*F
6: A Poiseuille flow is a steady-state isoviscous Stokes flow in a pipe of constant cross-section. We discretize using the
7: finite element method on an unstructured mesh. The weak form equations are
8: \begin{align*}
9: < \nabla v, \nu (\nabla u + {\nabla u}^T) > - < \nabla\cdot v, p > + < v, \Delta \hat n >_{\Gamma_o} = 0
10: < q, \nabla\cdot u > = 0
11: \end{align*}
12: where $\nu$ is the kinematic viscosity, $\Delta$ is the pressure drop per unit length, assuming that pressure is 0 on
13: the left edge, and $\Gamma_o$ is the outlet boundary at the right edge of the pipe. The normal velocity will be zero at
14: the wall, but we will allow a fixed tangential velocity $u_0$.
16: In order to test our global to local basis transformation, we will allow the pipe to be at an angle $\alpha$ to the
17: coordinate axes.
19: For visualization, use
21: -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append
22: F*/
24: #include <petscdmplex.h>
25: #include <petscsnes.h>
26: #include <petscds.h>
27: #include <petscbag.h>
29: typedef struct {
30: PetscReal Delta; /* Pressure drop per unit length */
31: PetscReal nu; /* Kinematic viscosity */
32: PetscReal u_0; /* Tangential velocity at the wall */
33: PetscReal alpha; /* Angle of pipe wall to x-axis */
34: } Parameter;
36: typedef struct {
37: PetscBag bag; /* Holds problem parameters */
38: } AppCtx;
40: /*
41: In 2D, plane Poiseuille flow has exact solution:
43: u = \Delta/(2 \nu) y (1 - y) + u_0
44: v = 0
45: p = -\Delta x
46: f = 0
48: so that
50: -\nu \Delta u + \nabla p + f = <\Delta, 0> + <-\Delta, 0> + <0, 0> = 0
51: \nabla \cdot u = 0 + 0 = 0
53: In 3D we use exact solution:
55: u = \Delta/(4 \nu) (y (1 - y) + z (1 - z)) + u_0
56: v = 0
57: w = 0
58: p = -\Delta x
59: f = 0
61: so that
63: -\nu \Delta u + \nabla p + f = <Delta, 0, 0> + <-Delta, 0, 0> + <0, 0, 0> = 0
64: \nabla \cdot u = 0 + 0 + 0 = 0
66: Note that these functions use coordinates X in the global (rotated) frame
67: */
68: PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
69: {
70: Parameter *param = (Parameter *)ctx;
71: PetscReal Delta = param->Delta;
72: PetscReal nu = param->nu;
73: PetscReal u_0 = param->u_0;
74: PetscReal fac = (PetscReal)(dim - 1);
75: PetscInt d;
77: u[0] = u_0;
78: for (d = 1; d < dim; ++d) u[0] += Delta / (fac * 2.0 * nu) * X[d] * (1.0 - X[d]);
79: for (d = 1; d < dim; ++d) u[d] = 0.0;
80: return 0;
81: }
83: PetscErrorCode linear_p(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
84: {
85: Parameter *param = (Parameter *)ctx;
86: PetscReal Delta = param->Delta;
88: p[0] = -Delta * X[0];
89: return 0;
90: }
92: PetscErrorCode wall_velocity(PetscInt dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
93: {
94: Parameter *param = (Parameter *)ctx;
95: PetscReal u_0 = param->u_0;
96: PetscInt d;
98: u[0] = u_0;
99: for (d = 1; d < dim; ++d) u[d] = 0.0;
100: return 0;
101: }
103: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z}
104: u[Ncomp] = {p} */
105: void f1_u(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[])
106: {
107: const PetscReal nu = PetscRealPart(constants[1]);
108: const PetscInt Nc = dim;
109: PetscInt c, d;
111: for (c = 0; c < Nc; ++c) {
112: for (d = 0; d < dim; ++d) {
113: /* f1[c*dim+d] = 0.5*nu*(u_x[c*dim+d] + u_x[d*dim+c]); */
114: f1[c * dim + d] = nu * u_x[c * dim + d];
115: }
116: f1[c * dim + c] -= u[uOff[1]];
117: }
118: }
120: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z} */
121: void f0_p(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[])
122: {
123: PetscInt d;
124: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
125: }
127: /* Residual functions are in reference coordinates */
128: static void f0_bd_u(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[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
129: {
130: const PetscReal Delta = PetscRealPart(constants[0]);
131: PetscReal alpha = PetscRealPart(constants[3]);
132: PetscReal X = PetscCosReal(alpha) * x[0] + PetscSinReal(alpha) * x[1];
133: PetscInt d;
135: for (d = 0; d < dim; ++d) f0[d] = -Delta * X * n[d];
136: }
138: /* < q, \nabla\cdot u >
139: NcompI = 1, NcompJ = dim */
140: void g1_pu(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[])
141: {
142: PetscInt d;
143: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
144: }
146: /* -< \nabla\cdot v, p >
147: NcompI = dim, NcompJ = 1 */
148: void g2_up(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[])
149: {
150: PetscInt d;
151: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
152: }
154: /* < \nabla v, \nabla u + {\nabla u}^T >
155: This just gives \nabla u, give the perdiagonal for the transpose */
156: void g3_uu(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 g3[])
157: {
158: const PetscReal nu = PetscRealPart(constants[1]);
159: const PetscInt Nc = dim;
160: PetscInt c, d;
162: for (c = 0; c < Nc; ++c) {
163: for (d = 0; d < dim; ++d) g3[((c * Nc + c) * dim + d) * dim + d] = nu;
164: }
165: }
167: static PetscErrorCode SetupParameters(AppCtx *user)
168: {
169: PetscBag bag;
170: Parameter *p;
173: /* setup PETSc parameter bag */
174: PetscBagGetData(user->bag, (void **)&p);
175: PetscBagSetName(user->bag, "par", "Poiseuille flow parameters");
176: bag = user->bag;
177: PetscBagRegisterReal(bag, &p->Delta, 1.0, "Delta", "Pressure drop per unit length");
178: PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity");
179: PetscBagRegisterReal(bag, &p->u_0, 0.0, "u_0", "Tangential velocity at the wall");
180: PetscBagRegisterReal(bag, &p->alpha, 0.0, "alpha", "Angle of pipe wall to x-axis");
181: return 0;
182: }
184: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
185: {
187: DMCreate(comm, dm);
188: DMSetType(*dm, DMPLEX);
189: DMSetFromOptions(*dm);
190: {
191: Parameter *param;
192: Vec coordinates;
193: PetscScalar *coords;
194: PetscReal alpha;
195: PetscInt cdim, N, bs, i;
197: DMGetCoordinateDim(*dm, &cdim);
198: DMGetCoordinates(*dm, &coordinates);
199: VecGetLocalSize(coordinates, &N);
200: VecGetBlockSize(coordinates, &bs);
202: VecGetArray(coordinates, &coords);
203: PetscBagGetData(user->bag, (void **)¶m);
204: alpha = param->alpha;
205: for (i = 0; i < N; i += cdim) {
206: PetscScalar x = coords[i + 0];
207: PetscScalar y = coords[i + 1];
209: coords[i + 0] = PetscCosReal(alpha) * x - PetscSinReal(alpha) * y;
210: coords[i + 1] = PetscSinReal(alpha) * x + PetscCosReal(alpha) * y;
211: }
212: VecRestoreArray(coordinates, &coords);
213: DMSetCoordinates(*dm, coordinates);
214: }
215: DMViewFromOptions(*dm, NULL, "-dm_view");
216: return 0;
217: }
219: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
220: {
221: PetscDS ds;
222: PetscWeakForm wf;
223: DMLabel label;
224: Parameter *ctx;
225: PetscInt id, bd;
228: PetscBagGetData(user->bag, (void **)&ctx);
229: DMGetDS(dm, &ds);
230: PetscDSSetResidual(ds, 0, NULL, f1_u);
231: PetscDSSetResidual(ds, 1, f0_p, NULL);
232: PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);
233: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL);
234: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL);
236: id = 2;
237: DMGetLabel(dm, "marker", &label);
238: DMAddBoundary(dm, DM_BC_NATURAL, "right wall", label, 1, &id, 0, 0, NULL, NULL, NULL, ctx, &bd);
239: PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
240: PetscWeakFormSetIndexBdResidual(wf, label, id, 0, 0, 0, f0_bd_u, 0, NULL);
241: /* Setup constants */
242: {
243: Parameter *param;
244: PetscScalar constants[4];
246: PetscBagGetData(user->bag, (void **)¶m);
248: constants[0] = param->Delta;
249: constants[1] = param->nu;
250: constants[2] = param->u_0;
251: constants[3] = param->alpha;
252: PetscDSSetConstants(ds, 4, constants);
253: }
254: /* Setup Boundary Conditions */
255: id = 3;
256: DMAddBoundary(dm, DM_BC_ESSENTIAL, "top wall", label, 1, &id, 0, 0, NULL, (void (*)(void))wall_velocity, NULL, ctx, NULL);
257: id = 1;
258: DMAddBoundary(dm, DM_BC_ESSENTIAL, "bottom wall", label, 1, &id, 0, 0, NULL, (void (*)(void))wall_velocity, NULL, ctx, NULL);
259: /* Setup exact solution */
260: PetscDSSetExactSolution(ds, 0, quadratic_u, ctx);
261: PetscDSSetExactSolution(ds, 1, linear_p, ctx);
262: return 0;
263: }
265: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
266: {
267: DM cdm = dm;
268: PetscFE fe[2];
269: Parameter *param;
270: PetscBool simplex;
271: PetscInt dim;
272: MPI_Comm comm;
275: DMGetDimension(dm, &dim);
276: DMPlexIsSimplex(dm, &simplex);
277: PetscObjectGetComm((PetscObject)dm, &comm);
278: PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);
279: PetscObjectSetName((PetscObject)fe[0], "velocity");
280: PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
281: PetscFECopyQuadrature(fe[0], fe[1]);
282: PetscObjectSetName((PetscObject)fe[1], "pressure");
283: /* Set discretization and boundary conditions for each mesh */
284: DMSetField(dm, 0, NULL, (PetscObject)fe[0]);
285: DMSetField(dm, 1, NULL, (PetscObject)fe[1]);
286: DMCreateDS(dm);
287: SetupProblem(dm, user);
288: PetscBagGetData(user->bag, (void **)¶m);
289: while (cdm) {
290: DMCopyDisc(dm, cdm);
291: DMPlexCreateBasisRotation(cdm, param->alpha, 0.0, 0.0);
292: DMGetCoarseDM(cdm, &cdm);
293: }
294: PetscFEDestroy(&fe[0]);
295: PetscFEDestroy(&fe[1]);
296: return 0;
297: }
299: int main(int argc, char **argv)
300: {
301: SNES snes; /* nonlinear solver */
302: DM dm; /* problem definition */
303: Vec u, r; /* solution and residual */
304: AppCtx user; /* user-defined work context */
307: PetscInitialize(&argc, &argv, NULL, help);
308: PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);
309: SetupParameters(&user);
310: PetscBagSetFromOptions(user.bag);
311: SNESCreate(PETSC_COMM_WORLD, &snes);
312: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
313: SNESSetDM(snes, dm);
314: DMSetApplicationContext(dm, &user);
315: /* Setup problem */
316: SetupDiscretization(dm, &user);
317: DMPlexCreateClosureIndex(dm, NULL);
319: DMCreateGlobalVector(dm, &u);
320: VecDuplicate(u, &r);
322: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
324: SNESSetFromOptions(snes);
326: {
327: PetscDS ds;
328: PetscSimplePointFunc exactFuncs[2];
329: void *ctxs[2];
331: DMGetDS(dm, &ds);
332: PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]);
333: PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]);
334: DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u);
335: PetscObjectSetName((PetscObject)u, "Exact Solution");
336: VecViewFromOptions(u, NULL, "-exact_vec_view");
337: }
338: DMSNESCheckFromOptions(snes, u);
339: VecSet(u, 0.0);
340: PetscObjectSetName((PetscObject)u, "Solution");
341: SNESSolve(snes, NULL, u);
342: VecViewFromOptions(u, NULL, "-sol_vec_view");
344: VecDestroy(&u);
345: VecDestroy(&r);
346: DMDestroy(&dm);
347: SNESDestroy(&snes);
348: PetscBagDestroy(&user.bag);
349: PetscFinalize();
350: return 0;
351: }
353: /*TEST
355: # Convergence
356: test:
357: suffix: 2d_quad_q1_p0_conv
358: requires: !single
359: args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 \
360: -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
361: -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
362: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
363: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
364: -fieldsplit_velocity_pc_type lu \
365: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
366: test:
367: suffix: 2d_quad_q1_p0_conv_u0
368: requires: !single
369: args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 -u_0 0.125 \
370: -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
371: -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
372: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
373: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
374: -fieldsplit_velocity_pc_type lu \
375: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
376: test:
377: suffix: 2d_quad_q1_p0_conv_u0_alpha
378: requires: !single
379: args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_refine 1 -u_0 0.125 -alpha 0.3927 \
380: -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
381: -snes_convergence_estimate -convest_num_refine 2 -snes_error_if_not_converged \
382: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
383: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
384: -fieldsplit_velocity_pc_type lu \
385: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
386: test:
387: suffix: 2d_quad_q1_p0_conv_gmg_vanka
388: requires: !single long_runtime
389: args: -dm_plex_simplex 0 -dm_plex_separate_marker -dm_plex_box_faces 2,2 -dm_refine_hierarchy 1 \
390: -vel_petscspace_degree 1 -pres_petscspace_degree 0 \
391: -snes_convergence_estimate -convest_num_refine 1 -snes_error_if_not_converged \
392: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
393: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
394: -fieldsplit_velocity_pc_type mg \
395: -fieldsplit_velocity_mg_levels_pc_type patch -fieldsplit_velocity_mg_levels_pc_patch_exclude_subspaces 1 \
396: -fieldsplit_velocity_mg_levels_pc_patch_construct_codim 0 -fieldsplit_velocity_mg_levels_pc_patch_construct_type vanka \
397: -fieldsplit_pressure_ksp_rtol 1e-5 -fieldsplit_pressure_pc_type jacobi
398: test:
399: suffix: 2d_tri_p2_p1_conv
400: requires: triangle !single
401: args: -dm_plex_separate_marker -dm_refine 1 \
402: -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
403: -dmsnes_check .001 -snes_error_if_not_converged \
404: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
405: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
406: -fieldsplit_velocity_pc_type lu \
407: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
408: test:
409: suffix: 2d_tri_p2_p1_conv_u0_alpha
410: requires: triangle !single
411: args: -dm_plex_separate_marker -dm_refine 0 -u_0 0.125 -alpha 0.3927 \
412: -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
413: -dmsnes_check .001 -snes_error_if_not_converged \
414: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
415: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
416: -fieldsplit_velocity_pc_type lu \
417: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
418: test:
419: suffix: 2d_tri_p2_p1_conv_gmg_vcycle
420: requires: triangle !single
421: args: -dm_plex_separate_marker -dm_plex_box_faces 2,2 -dm_refine_hierarchy 1 \
422: -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
423: -dmsnes_check .001 -snes_error_if_not_converged \
424: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
425: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
426: -fieldsplit_velocity_pc_type mg \
427: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
428: TEST*/