Actual source code: ex46.c
1: static char help[] = "Time dependent Navier-Stokes problem in 2d and 3d with finite elements.\n\
2: We solve the Navier-Stokes in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
4: This example supports discretized auxiliary fields (Re) as well as\n\
5: multilevel nonlinear solvers.\n\
6: Contributed by: Julian Andrej <juan@tf.uni-kiel.de>\n\n\n";
8: #include <petscdmplex.h>
9: #include <petscsnes.h>
10: #include <petscts.h>
11: #include <petscds.h>
13: /*
14: Navier-Stokes equation:
16: du/dt + u . grad u - \Delta u - grad p = f
17: div u = 0
18: */
20: typedef struct {
21: PetscInt mms;
22: } AppCtx;
24: #define REYN 400.0
26: /* MMS1
28: u = t + x^2 + y^2;
29: v = t + 2*x^2 - 2*x*y;
30: p = x + y - 1;
32: f_x = -2*t*(x + y) + 2*x*y^2 - 4*x^2*y - 2*x^3 + 4.0/Re - 1.0
33: f_y = -2*t*x + 2*y^3 - 4*x*y^2 - 2*x^2*y + 4.0/Re - 1.0
35: so that
37: u_t + u \cdot \nabla u - 1/Re \Delta u + \nabla p + f = <1, 1> + <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t 2x + 2x^2y + 4xy^2 - 2y^3> - 1/Re <4, 4> + <1, 1>
38: + <-t (2x + 2y) + 2xy^2 - 4x^2y - 2x^3 + 4/Re - 1, -2xt + 2y^3 - 4xy^2 - 2x^2y + 4/Re - 1> = 0
39: \nabla \cdot u = 2x - 2x = 0
41: where
43: <u, v> . <<u_x, v_x>, <u_y, v_y>> = <u u_x + v u_y, u v_x + v v_y>
44: */
45: PetscErrorCode mms1_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
46: {
47: u[0] = time + x[0]*x[0] + x[1]*x[1];
48: u[1] = time + 2.0*x[0]*x[0] - 2.0*x[0]*x[1];
49: return 0;
50: }
52: PetscErrorCode mms1_p_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx)
53: {
54: *p = x[0] + x[1] - 1.0;
55: return 0;
56: }
58: /* MMS 2*/
60: static PetscErrorCode mms2_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
61: {
62: u[0] = PetscSinReal(time + x[0])*PetscSinReal(time + x[1]);
63: u[1] = PetscCosReal(time + x[0])*PetscCosReal(time + x[1]);
64: return 0;
65: }
67: static PetscErrorCode mms2_p_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx)
68: {
69: *p = PetscSinReal(time + x[0] - x[1]);
70: return 0;
71: }
73: static void f0_mms1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
74: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
75: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
76: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
77: {
78: const PetscReal Re = REYN;
79: const PetscInt Ncomp = dim;
80: PetscInt c, d;
82: for (c = 0; c < Ncomp; ++c) {
83: for (d = 0; d < dim; ++d) {
84: f0[c] += u[d] * u_x[c*dim+d];
85: }
86: }
87: f0[0] += u_t[0];
88: f0[1] += u_t[1];
90: f0[0] += -2.0*t*(x[0] + x[1]) + 2.0*x[0]*x[1]*x[1] - 4.0*x[0]*x[0]*x[1] - 2.0*x[0]*x[0]*x[0] + 4.0/Re - 1.0;
91: f0[1] += -2.0*t*x[0] + 2.0*x[1]*x[1]*x[1] - 4.0*x[0]*x[1]*x[1] - 2.0*x[0]*x[0]*x[1] + 4.0/Re - 1.0;
92: }
94: static void f0_mms2_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
95: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
96: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
97: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
98: {
99: const PetscReal Re = REYN;
100: const PetscInt Ncomp = dim;
101: PetscInt c, d;
103: for (c = 0; c < Ncomp; ++c) {
104: for (d = 0; d < dim; ++d) {
105: f0[c] += u[d] * u_x[c*dim+d];
106: }
107: }
108: f0[0] += u_t[0];
109: f0[1] += u_t[1];
111: f0[0] -= ( Re*((1.0L/2.0L)*PetscSinReal(2*t + 2*x[0]) + PetscSinReal(2*t + x[0] + x[1]) + PetscCosReal(t + x[0] - x[1])) + 2.0*PetscSinReal(t + x[0])*PetscSinReal(t + x[1]))/Re;
112: f0[1] -= (-Re*((1.0L/2.0L)*PetscSinReal(2*t + 2*x[1]) + PetscSinReal(2*t + x[0] + x[1]) + PetscCosReal(t + x[0] - x[1])) + 2.0*PetscCosReal(t + x[0])*PetscCosReal(t + x[1]))/Re;
113: }
115: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
116: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
117: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
118: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
119: {
120: const PetscReal Re = REYN;
121: const PetscInt Ncomp = dim;
122: PetscInt comp, d;
124: for (comp = 0; comp < Ncomp; ++comp) {
125: for (d = 0; d < dim; ++d) {
126: f1[comp*dim+d] = 1.0/Re * u_x[comp*dim+d];
127: }
128: f1[comp*dim+comp] -= u[Ncomp];
129: }
130: }
132: static void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
133: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
134: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
135: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
136: {
137: PetscInt d;
138: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d*dim+d];
139: }
141: static void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux,
142: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
143: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
144: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
145: {
146: PetscInt d;
147: for (d = 0; d < dim; ++d) f1[d] = 0.0;
148: }
150: /*
151: (psi_i, u_j grad_j u_i) ==> (\psi_i, \phi_j grad_j u_i)
152: */
153: static void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
154: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
155: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
156: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
157: {
158: PetscInt NcI = dim, NcJ = dim;
159: PetscInt fc, gc;
160: PetscInt d;
162: for (d = 0; d < dim; ++d) {
163: g0[d*dim+d] = u_tShift;
164: }
166: for (fc = 0; fc < NcI; ++fc) {
167: for (gc = 0; gc < NcJ; ++gc) {
168: g0[fc*NcJ+gc] += u_x[fc*NcJ+gc];
169: }
170: }
171: }
173: /*
174: (psi_i, u_j grad_j u_i) ==> (\psi_i, \u_j grad_j \phi_i)
175: */
176: static void g1_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
177: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
178: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
179: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
180: {
181: PetscInt NcI = dim;
182: PetscInt NcJ = dim;
183: PetscInt fc, gc, dg;
184: for (fc = 0; fc < NcI; ++fc) {
185: for (gc = 0; gc < NcJ; ++gc) {
186: for (dg = 0; dg < dim; ++dg) {
187: /* kronecker delta */
188: if (fc == gc) {
189: g1[(fc*NcJ+gc)*dim+dg] += u[dg];
190: }
191: }
192: }
193: }
194: }
196: /* < q, \nabla\cdot u >
197: NcompI = 1, NcompJ = dim */
198: static void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
199: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
200: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
201: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
202: {
203: PetscInt d;
204: for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
205: }
207: /* -< \nabla\cdot v, p >
208: NcompI = dim, NcompJ = 1 */
209: static void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux,
210: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
211: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
212: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
213: {
214: PetscInt d;
215: for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
216: }
218: /* < \nabla v, \nabla u + {\nabla u}^T >
219: This just gives \nabla u, give the perdiagonal for the transpose */
220: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
221: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
222: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
223: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
224: {
225: const PetscReal Re = REYN;
226: const PetscInt Ncomp = dim;
227: PetscInt compI, d;
229: for (compI = 0; compI < Ncomp; ++compI) {
230: for (d = 0; d < dim; ++d) {
231: g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0/Re;
232: }
233: }
234: }
236: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
237: {
241: options->mms = 1;
243: PetscOptionsBegin(comm, "", "Navier-Stokes Equation Options", "DMPLEX");
244: PetscOptionsInt("-mms", "The manufactured solution to use", "ex46.c", options->mms, &options->mms, NULL);
245: PetscOptionsEnd();
246: return(0);
247: }
249: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
250: {
254: DMCreate(comm, dm);
255: DMSetType(*dm, DMPLEX);
256: DMSetFromOptions(*dm);
257: DMViewFromOptions(*dm, NULL, "-dm_view");
258: return(0);
259: }
261: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
262: {
263: PetscDS ds;
264: DMLabel label;
265: const PetscInt id = 1;
266: PetscInt dim;
270: DMGetDimension(dm, &dim);
271: DMGetDS(dm, &ds);
272: DMGetLabel(dm, "marker", &label);
273: switch (dim) {
274: case 2:
275: switch (ctx->mms) {
276: case 1:
277: PetscDSSetResidual(ds, 0, f0_mms1_u, f1_u);
278: PetscDSSetResidual(ds, 1, f0_p, f1_p);
279: PetscDSSetJacobian(ds, 0, 0, g0_uu, g1_uu, NULL, g3_uu);
280: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL);
281: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL);
282: PetscDSSetExactSolution(ds, 0, mms1_u_2d, ctx);
283: PetscDSSetExactSolution(ds, 1, mms1_p_2d, ctx);
284: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) mms1_u_2d, NULL, ctx, NULL);
285: break;
286: case 2:
287: PetscDSSetResidual(ds, 0, f0_mms2_u, f1_u);
288: PetscDSSetResidual(ds, 1, f0_p, f1_p);
289: PetscDSSetJacobian(ds, 0, 0, g0_uu, g1_uu, NULL, g3_uu);
290: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL);
291: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL);
292: PetscDSSetExactSolution(ds, 0, mms2_u_2d, ctx);
293: PetscDSSetExactSolution(ds, 1, mms2_p_2d, ctx);
294: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) mms2_u_2d, NULL, ctx, NULL);
295: break;
296: default: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid MMS %D", ctx->mms);
297: }
298: break;
299: default: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %D", dim);
300: }
301: return(0);
302: }
304: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
305: {
306: MPI_Comm comm;
307: DM cdm = dm;
308: PetscFE fe[2];
309: PetscInt dim;
310: PetscBool simplex;
311: PetscErrorCode ierr;
314: PetscObjectGetComm((PetscObject) dm, &comm);
315: DMGetDimension(dm, &dim);
316: DMPlexIsSimplex(dm, &simplex);
317: PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);
318: PetscObjectSetName((PetscObject) fe[0], "velocity");
319: PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
320: PetscFECopyQuadrature(fe[0], fe[1]);
321: PetscObjectSetName((PetscObject) fe[1], "pressure");
322: /* Set discretization and boundary conditions for each mesh */
323: DMSetField(dm, 0, NULL, (PetscObject) fe[0]);
324: DMSetField(dm, 1, NULL, (PetscObject) fe[1]);
325: DMCreateDS(dm);
326: SetupProblem(dm, ctx);
327: while (cdm) {
328: PetscObject pressure;
329: MatNullSpace nsp;
331: DMGetField(cdm, 1, NULL, &pressure);
332: MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nsp);
333: PetscObjectCompose(pressure, "nullspace", (PetscObject) nsp);
334: MatNullSpaceDestroy(&nsp);
336: DMCopyDisc(dm, cdm);
337: DMGetCoarseDM(cdm, &cdm);
338: }
339: PetscFEDestroy(&fe[0]);
340: PetscFEDestroy(&fe[1]);
341: return(0);
342: }
344: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
345: {
346: PetscSimplePointFunc funcs[2];
347: void *ctxs[2];
348: DM dm;
349: PetscDS ds;
350: PetscReal ferrors[2];
351: PetscErrorCode ierr;
354: TSGetDM(ts, &dm);
355: DMGetDS(dm, &ds);
356: PetscDSGetExactSolution(ds, 0, &funcs[0], &ctxs[0]);
357: PetscDSGetExactSolution(ds, 1, &funcs[1], &ctxs[1]);
358: DMComputeL2FieldDiff(dm, crtime, funcs, ctxs, u, ferrors);
359: PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g]\n", (int) step, (double) crtime, (double) ferrors[0], (double) ferrors[1]);
360: return(0);
361: }
363: int main(int argc, char **argv)
364: {
365: AppCtx ctx;
366: DM dm;
367: TS ts;
368: Vec u, r;
371: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
372: ProcessOptions(PETSC_COMM_WORLD, &ctx);
373: CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);
374: DMSetApplicationContext(dm, &ctx);
375: SetupDiscretization(dm, &ctx);
376: DMPlexCreateClosureIndex(dm, NULL);
378: DMCreateGlobalVector(dm, &u);
379: VecDuplicate(u, &r);
381: TSCreate(PETSC_COMM_WORLD, &ts);
382: TSMonitorSet(ts, MonitorError, &ctx, NULL);
383: TSSetDM(ts, dm);
384: DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
385: DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
386: DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
387: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
388: TSSetFromOptions(ts);
389: DMTSCheckFromOptions(ts, u);
391: {
392: PetscSimplePointFunc funcs[2];
393: void *ctxs[2];
394: PetscDS ds;
396: DMGetDS(dm, &ds);
397: PetscDSGetExactSolution(ds, 0, &funcs[0], &ctxs[0]);
398: PetscDSGetExactSolution(ds, 1, &funcs[1], &ctxs[1]);
399: DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, u);
400: }
401: TSSolve(ts, u);
402: VecViewFromOptions(u, NULL, "-sol_vec_view");
404: VecDestroy(&u);
405: VecDestroy(&r);
406: TSDestroy(&ts);
407: DMDestroy(&dm);
408: PetscFinalize();
409: return ierr;
410: }
412: /*TEST
414: # Full solves
415: test:
416: suffix: 2d_p2p1_r1
417: requires: !single triangle
418: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
419: args: -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
420: -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ts_monitor -dmts_check \
421: -snes_monitor_short -snes_converged_reason \
422: -ksp_monitor_short -ksp_converged_reason \
423: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full \
424: -fieldsplit_velocity_pc_type lu \
425: -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi
427: test:
428: suffix: 2d_q2q1_r1
429: requires: !single
430: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" -e "s~ 0\]~ 0.0\]~g"
431: args: -dm_plex_simplex 0 -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
432: -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ts_monitor -dmts_check \
433: -snes_monitor_short -snes_converged_reason \
434: -ksp_monitor_short -ksp_converged_reason \
435: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full \
436: -fieldsplit_velocity_pc_type lu \
437: -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi
439: TEST*/