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: {
252: DMCreate(comm, dm);
253: DMSetType(*dm, DMPLEX);
254: DMSetFromOptions(*dm);
255: DMViewFromOptions(*dm, NULL, "-dm_view");
256: return 0;
257: }
259: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
260: {
261: PetscDS ds;
262: DMLabel label;
263: const PetscInt id = 1;
264: PetscInt dim;
267: DMGetDimension(dm, &dim);
268: DMGetDS(dm, &ds);
269: DMGetLabel(dm, "marker", &label);
270: switch (dim) {
271: case 2:
272: switch (ctx->mms) {
273: case 1:
274: PetscDSSetResidual(ds, 0, f0_mms1_u, f1_u);
275: PetscDSSetResidual(ds, 1, f0_p, f1_p);
276: PetscDSSetJacobian(ds, 0, 0, g0_uu, g1_uu, NULL, g3_uu);
277: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL);
278: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL);
279: PetscDSSetExactSolution(ds, 0, mms1_u_2d, ctx);
280: PetscDSSetExactSolution(ds, 1, mms1_p_2d, ctx);
281: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) mms1_u_2d, NULL, ctx, NULL);
282: break;
283: case 2:
284: PetscDSSetResidual(ds, 0, f0_mms2_u, f1_u);
285: PetscDSSetResidual(ds, 1, f0_p, f1_p);
286: PetscDSSetJacobian(ds, 0, 0, g0_uu, g1_uu, NULL, g3_uu);
287: PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL);
288: PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL);
289: PetscDSSetExactSolution(ds, 0, mms2_u_2d, ctx);
290: PetscDSSetExactSolution(ds, 1, mms2_p_2d, ctx);
291: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) mms2_u_2d, NULL, ctx, NULL);
292: break;
293: default: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid MMS %D", ctx->mms);
294: }
295: break;
296: default: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %D", dim);
297: }
298: return 0;
299: }
301: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
302: {
303: MPI_Comm comm;
304: DM cdm = dm;
305: PetscFE fe[2];
306: PetscInt dim;
307: PetscBool simplex;
310: PetscObjectGetComm((PetscObject) dm, &comm);
311: DMGetDimension(dm, &dim);
312: DMPlexIsSimplex(dm, &simplex);
313: PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);
314: PetscObjectSetName((PetscObject) fe[0], "velocity");
315: PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
316: PetscFECopyQuadrature(fe[0], fe[1]);
317: PetscObjectSetName((PetscObject) fe[1], "pressure");
318: /* Set discretization and boundary conditions for each mesh */
319: DMSetField(dm, 0, NULL, (PetscObject) fe[0]);
320: DMSetField(dm, 1, NULL, (PetscObject) fe[1]);
321: DMCreateDS(dm);
322: SetupProblem(dm, ctx);
323: while (cdm) {
324: PetscObject pressure;
325: MatNullSpace nsp;
327: DMGetField(cdm, 1, NULL, &pressure);
328: MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nsp);
329: PetscObjectCompose(pressure, "nullspace", (PetscObject) nsp);
330: MatNullSpaceDestroy(&nsp);
332: DMCopyDisc(dm, cdm);
333: DMGetCoarseDM(cdm, &cdm);
334: }
335: PetscFEDestroy(&fe[0]);
336: PetscFEDestroy(&fe[1]);
337: return 0;
338: }
340: static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx)
341: {
342: PetscSimplePointFunc funcs[2];
343: void *ctxs[2];
344: DM dm;
345: PetscDS ds;
346: PetscReal ferrors[2];
349: TSGetDM(ts, &dm);
350: DMGetDS(dm, &ds);
351: PetscDSGetExactSolution(ds, 0, &funcs[0], &ctxs[0]);
352: PetscDSGetExactSolution(ds, 1, &funcs[1], &ctxs[1]);
353: DMComputeL2FieldDiff(dm, crtime, funcs, ctxs, u, ferrors);
354: 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]);
355: return 0;
356: }
358: int main(int argc, char **argv)
359: {
360: AppCtx ctx;
361: DM dm;
362: TS ts;
363: Vec u, r;
365: PetscInitialize(&argc, &argv, NULL, help);
366: ProcessOptions(PETSC_COMM_WORLD, &ctx);
367: CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);
368: DMSetApplicationContext(dm, &ctx);
369: SetupDiscretization(dm, &ctx);
370: DMPlexCreateClosureIndex(dm, NULL);
372: DMCreateGlobalVector(dm, &u);
373: VecDuplicate(u, &r);
375: TSCreate(PETSC_COMM_WORLD, &ts);
376: TSMonitorSet(ts, MonitorError, &ctx, NULL);
377: TSSetDM(ts, dm);
378: DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
379: DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
380: DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
381: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
382: TSSetFromOptions(ts);
383: DMTSCheckFromOptions(ts, u);
385: {
386: PetscSimplePointFunc funcs[2];
387: void *ctxs[2];
388: PetscDS ds;
390: DMGetDS(dm, &ds);
391: PetscDSGetExactSolution(ds, 0, &funcs[0], &ctxs[0]);
392: PetscDSGetExactSolution(ds, 1, &funcs[1], &ctxs[1]);
393: DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, u);
394: }
395: TSSolve(ts, u);
396: VecViewFromOptions(u, NULL, "-sol_vec_view");
398: VecDestroy(&u);
399: VecDestroy(&r);
400: TSDestroy(&ts);
401: DMDestroy(&dm);
402: PetscFinalize();
403: return 0;
404: }
406: /*TEST
408: # Full solves
409: test:
410: suffix: 2d_p2p1_r1
411: requires: !single triangle
412: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
413: args: -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
414: -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ts_monitor -dmts_check \
415: -snes_monitor_short -snes_converged_reason \
416: -ksp_monitor_short -ksp_converged_reason \
417: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full \
418: -fieldsplit_velocity_pc_type lu \
419: -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi
421: test:
422: suffix: 2d_q2q1_r1
423: requires: !single
424: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" -e "s~ 0\]~ 0.0\]~g"
425: args: -dm_plex_simplex 0 -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \
426: -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ts_monitor -dmts_check \
427: -snes_monitor_short -snes_converged_reason \
428: -ksp_monitor_short -ksp_converged_reason \
429: -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full \
430: -fieldsplit_velocity_pc_type lu \
431: -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi
433: TEST*/