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: /*
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: /*
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]);
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*/