Actual source code: ex26.c
1: static char help[] = "Transient nonlinear driven cavity in 2d.\n\
2: \n\
3: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
4: The flow can be driven with the lid or with buoyancy or both:\n\
5: -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
6: -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
7: -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
8: -contours : draw contour plots of solution\n\n";
9: /*
10: See src/snes/tutorials/ex19.c for the steady-state version.
12: There used to be a SNES example (src/snes/tutorials/ex27.c) that
13: implemented this algorithm without using TS and was used for the numerical
14: results in the paper
16: Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient
17: Continuation and Differential-Algebraic Equations, 2003.
19: That example was removed because it used obsolete interfaces, but the
20: algorithms from the paper can be reproduced using this example.
22: Note: The paper describes the algorithm as being linearly implicit but the
23: numerical results were created using nonlinearly implicit Euler. The
24: algorithm as described (linearly implicit) is more efficient and is the
25: default when using TSPSEUDO. If you want to reproduce the numerical
26: results from the paper, you'll have to change the SNES to converge the
27: nonlinear solve (e.g., -snes_type newtonls). The DAE versus ODE variants
28: are controlled using the -parabolic option.
30: Comment preserved from snes/tutorials/ex27.c, since removed:
32: [H]owever Figure 3.1 was generated with a slightly different algorithm
33: (see targets runex27 and runex27_p) than described in the paper. In
34: particular, the described algorithm is linearly implicit, advancing to
35: the next step after one Newton step, so that the steady state residual
36: is always used, but the figure was generated by converging each step to
37: a relative tolerance of 1.e-3. On the example problem, setting
38: -snes_type ksponly has only minor impact on number of steps, but
39: significantly reduces the required number of linear solves.
41: See also https://mailman.cels.anl.gov/archives/list/petsc-dev@lists.mcs.anl.gov/message/6YXJJQBZCWRTJSW3E3FWA6YB7SRRCP3A/
42: */
44: /* ------------------------------------------------------------------------
46: We thank David E. Keyes for contributing the driven cavity discretization
47: within this example code.
49: This problem is modeled by the partial differential equation system
51: - Lap(U) - Grad_y(Omega) = 0
52: - Lap(V) + Grad_x(Omega) = 0
53: Omega_t - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
54: T_t - Lap(T) + PR*Div([U*T,V*T]) = 0
56: in the unit square, which is uniformly discretized in each of x and
57: y in this simple encoding.
59: No-slip, rigid-wall Dirichlet conditions are used for [U,V].
60: Dirichlet conditions are used for Omega, based on the definition of
61: vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
62: constant coordinate boundary, the tangential derivative is zero.
63: Dirichlet conditions are used for T on the left and right walls,
64: and insulation homogeneous Neumann conditions are used for T on
65: the top and bottom walls.
67: A finite difference approximation with the usual 5-point stencil
68: is used to discretize the boundary value problem to obtain a
69: nonlinear system of equations. Upwinding is used for the divergence
70: (convective) terms and central for the gradient (source) terms.
72: The Jacobian can be either
73: * formed via finite differencing using coloring (the default), or
74: * applied matrix-free via the option -snes_mf
75: (for larger grid problems this variant may not converge
76: without a preconditioner due to ill-conditioning).
78: ------------------------------------------------------------------------- */
80: /*
81: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
82: Include "petscts.h" so that we can use the TS ODE integrators. Note that this
83: file automatically includes:
84: petscsys.h - base PETSc routines
85: petscis.h - index sets
86: petscvec.h - vectors
87: petscmat.h - matrices
88: petscksp.h - linear solvers
89: petscsnes.h - nonlinear solvers
90: */
91: #include <petscts.h>
92: #include <petscdm.h>
93: #include <petscdmda.h>
95: /*
96: User-defined routines and data structures
97: */
98: typedef struct {
99: PetscScalar u, v, omega, temp;
100: } Field;
102: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, void *);
104: typedef struct {
105: PetscReal lidvelocity, prandtl, grashof; /* physical parameters */
106: PetscBool parabolic; /* allow a transient term corresponding roughly to artificial compressibility */
107: PetscReal cfl_initial; /* CFL for first time step */
108: } AppCtx;
110: PetscErrorCode FormInitialSolution(TS, Vec, AppCtx *);
112: int main(int argc, char **argv)
113: {
114: AppCtx user; /* user-defined work context */
115: PetscInt mx, my, steps;
116: TS ts;
117: DM da;
118: Vec X;
119: PetscReal ftime;
120: TSConvergedReason reason;
122: PetscFunctionBeginUser;
123: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
124: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
125: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 4, 1, NULL, NULL, &da));
126: PetscCall(DMSetFromOptions(da));
127: PetscCall(DMSetUp(da));
128: PetscCall(TSSetDM(ts, (DM)da));
130: PetscCall(DMDAGetInfo(da, NULL, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
131: /*
132: Problem parameters (velocity of lid, prandtl, and grashof numbers)
133: */
134: user.lidvelocity = 1.0 / (mx * my);
135: user.prandtl = 1.0;
136: user.grashof = 1.0;
137: user.parabolic = PETSC_FALSE;
138: user.cfl_initial = 50.;
140: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Driven cavity/natural convection options", "");
141: PetscCall(PetscOptionsReal("-lidvelocity", "Lid velocity, related to Reynolds number", "", user.lidvelocity, &user.lidvelocity, NULL));
142: PetscCall(PetscOptionsReal("-prandtl", "Ratio of viscous to thermal diffusivity", "", user.prandtl, &user.prandtl, NULL));
143: PetscCall(PetscOptionsReal("-grashof", "Ratio of buoyant to viscous forces", "", user.grashof, &user.grashof, NULL));
144: PetscCall(PetscOptionsBool("-parabolic", "Relax incompressibility to make the system parabolic instead of differential-algebraic", "", user.parabolic, &user.parabolic, NULL));
145: PetscCall(PetscOptionsReal("-cfl_initial", "Advective CFL for the first time step", "", user.cfl_initial, &user.cfl_initial, NULL));
146: PetscOptionsEnd();
148: PetscCall(DMDASetFieldName(da, 0, "x-velocity"));
149: PetscCall(DMDASetFieldName(da, 1, "y-velocity"));
150: PetscCall(DMDASetFieldName(da, 2, "Omega"));
151: PetscCall(DMDASetFieldName(da, 3, "temperature"));
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Create user context, set problem data, create vector data structures.
155: Also, compute the initial guess.
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Create time integration context
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: PetscCall(DMSetApplicationContext(da, &user));
162: PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (DMDATSIFunctionLocalFn *)FormIFunctionLocal, &user));
163: PetscCall(TSSetMaxSteps(ts, 10000));
164: PetscCall(TSSetMaxTime(ts, 1e12));
165: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
166: PetscCall(TSSetTimeStep(ts, user.cfl_initial / (user.lidvelocity * mx)));
167: PetscCall(TSSetFromOptions(ts));
169: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "x%" PetscInt_FMT " grid, lid velocity = %g, prandtl # = %g, grashof # = %g\n", mx, my, (double)user.lidvelocity, (double)user.prandtl, (double)user.grashof));
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Solve the nonlinear system
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: PetscCall(DMCreateGlobalVector(da, &X));
176: PetscCall(FormInitialSolution(ts, X, &user));
178: PetscCall(TSSolve(ts, X));
179: PetscCall(TSGetSolveTime(ts, &ftime));
180: PetscCall(TSGetStepNumber(ts, &steps));
181: PetscCall(TSGetConvergedReason(ts, &reason));
183: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));
185: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Free work space. All PETSc objects should be destroyed when they
187: are no longer needed.
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: PetscCall(VecDestroy(&X));
190: PetscCall(DMDestroy(&da));
191: PetscCall(TSDestroy(&ts));
193: PetscCall(PetscFinalize());
194: return 0;
195: }
197: /* ------------------------------------------------------------------- */
199: /*
200: FormInitialSolution - Forms initial approximation.
202: Input Parameters:
203: user - user-defined application context
204: X - vector
206: Output Parameter:
207: X - vector
208: */
209: PetscErrorCode FormInitialSolution(TS ts, Vec X, AppCtx *user)
210: {
211: DM da;
212: PetscInt i, j, mx, xs, ys, xm, ym;
213: PetscReal grashof, dx;
214: Field **x;
216: PetscFunctionBeginUser;
217: grashof = user->grashof;
218: PetscCall(TSGetDM(ts, &da));
219: PetscCall(DMDAGetInfo(da, NULL, &mx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
220: dx = 1.0 / (mx - 1);
222: /*
223: Get local grid boundaries (for 2-dimensional DMDA):
224: xs, ys - starting grid indices (no ghost points)
225: xm, ym - widths of local grid (no ghost points)
226: */
227: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
229: /*
230: Get a pointer to vector data.
231: - For default PETSc vectors, VecGetArray() returns a pointer to
232: the data array. Otherwise, the routine is implementation dependent.
233: - You MUST call VecRestoreArray() when you no longer need access to
234: the array.
235: */
236: PetscCall(DMDAVecGetArray(da, X, &x));
238: /*
239: Compute initial guess over the locally owned part of the grid
240: Initial condition is motionless fluid and equilibrium temperature
241: */
242: for (j = ys; j < ys + ym; j++) {
243: for (i = xs; i < xs + xm; i++) {
244: x[j][i].u = 0.0;
245: x[j][i].v = 0.0;
246: x[j][i].omega = 0.0;
247: x[j][i].temp = (grashof > 0) * i * dx;
248: }
249: }
251: /*
252: Restore vector
253: */
254: PetscCall(DMDAVecRestoreArray(da, X, &x));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal ptime, Field **x, Field **xdot, Field **f, void *ptr)
259: {
260: AppCtx *user = (AppCtx *)ptr;
261: PetscInt xints, xinte, yints, yinte, i, j;
262: PetscReal hx, hy, dhx, dhy, hxdhy, hydhx;
263: PetscReal grashof, prandtl, lid;
264: PetscScalar u, udot, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym;
266: PetscFunctionBeginUser;
267: grashof = user->grashof;
268: prandtl = user->prandtl;
269: lid = user->lidvelocity;
271: /*
272: Define mesh intervals ratios for uniform grid.
274: Note: FD formulae below are normalized by multiplying through by
275: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
277: */
278: dhx = (PetscReal)(info->mx - 1);
279: dhy = (PetscReal)(info->my - 1);
280: hx = 1.0 / dhx;
281: hy = 1.0 / dhy;
282: hxdhy = hx * dhy;
283: hydhx = hy * dhx;
285: xints = info->xs;
286: xinte = info->xs + info->xm;
287: yints = info->ys;
288: yinte = info->ys + info->ym;
290: /* Test whether we are on the bottom edge of the global array */
291: if (yints == 0) {
292: j = 0;
293: yints = yints + 1;
294: /* bottom edge */
295: for (i = info->xs; i < info->xs + info->xm; i++) {
296: f[j][i].u = x[j][i].u;
297: f[j][i].v = x[j][i].v;
298: f[j][i].omega = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy;
299: f[j][i].temp = x[j][i].temp - x[j + 1][i].temp;
300: }
301: }
303: /* Test whether we are on the top edge of the global array */
304: if (yinte == info->my) {
305: j = info->my - 1;
306: yinte = yinte - 1;
307: /* top edge */
308: for (i = info->xs; i < info->xs + info->xm; i++) {
309: f[j][i].u = x[j][i].u - lid;
310: f[j][i].v = x[j][i].v;
311: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy;
312: f[j][i].temp = x[j][i].temp - x[j - 1][i].temp;
313: }
314: }
316: /* Test whether we are on the left edge of the global array */
317: if (xints == 0) {
318: i = 0;
319: xints = xints + 1;
320: /* left edge */
321: for (j = info->ys; j < info->ys + info->ym; j++) {
322: f[j][i].u = x[j][i].u;
323: f[j][i].v = x[j][i].v;
324: f[j][i].omega = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx;
325: f[j][i].temp = x[j][i].temp;
326: }
327: }
329: /* Test whether we are on the right edge of the global array */
330: if (xinte == info->mx) {
331: i = info->mx - 1;
332: xinte = xinte - 1;
333: /* right edge */
334: for (j = info->ys; j < info->ys + info->ym; j++) {
335: f[j][i].u = x[j][i].u;
336: f[j][i].v = x[j][i].v;
337: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx;
338: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof > 0);
339: }
340: }
342: /* Compute over the interior points */
343: for (j = yints; j < yinte; j++) {
344: for (i = xints; i < xinte; i++) {
345: /*
346: convective coefficients for upwinding
347: */
348: vx = x[j][i].u;
349: avx = PetscAbsScalar(vx);
350: vxp = .5 * (vx + avx);
351: vxm = .5 * (vx - avx);
352: vy = x[j][i].v;
353: avy = PetscAbsScalar(vy);
354: vyp = .5 * (vy + avy);
355: vym = .5 * (vy - avy);
357: /* U velocity */
358: u = x[j][i].u;
359: udot = user->parabolic ? xdot[j][i].u : 0.;
360: uxx = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
361: uyy = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
362: f[j][i].u = udot + uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx;
364: /* V velocity */
365: u = x[j][i].v;
366: udot = user->parabolic ? xdot[j][i].v : 0.;
367: uxx = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx;
368: uyy = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy;
369: f[j][i].v = udot + uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy;
371: /* Omega */
372: u = x[j][i].omega;
373: uxx = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx;
374: uyy = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy;
375: f[j][i].omega = (xdot[j][i].omega + uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (vyp * (u - x[j - 1][i].omega) + vym * (x[j + 1][i].omega - u)) * hx - .5 * grashof * (x[j][i + 1].temp - x[j][i - 1].temp) * hy);
377: /* Temperature */
378: u = x[j][i].temp;
379: uxx = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx;
380: uyy = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy;
381: f[j][i].temp = (xdot[j][i].temp + uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) * hy + (vyp * (u - x[j - 1][i].temp) + vym * (x[j + 1][i].temp - u)) * hx));
382: }
383: }
385: /*
386: Flop count (multiply-adds are counted as 2 operations)
387: */
388: PetscCall(PetscLogFlops(84.0 * info->ym * info->xm));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /*TEST
394: test:
395: args: -da_grid_x 20 -da_grid_y 20 -lidvelocity 100 -grashof 1e3 -ts_max_steps 100 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type rosw -ts_rosw_type ra3pw -ts_monitor -ts_monitor_solution_vtk 'foo-%03d.vts'
396: requires: !complex !single
398: test:
399: suffix: 1_steps
400: args: -da_grid_x 20 -da_grid_y 20 -lidvelocity 100 -grashof 1e3 -ts_run_steps 100 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type rosw -ts_rosw_type ra3pw -ts_monitor -ts_monitor_solution_vtk 'foo-%03d.vts'
401: requires: !complex !single
403: test:
404: suffix: 2
405: nsize: 4
406: args: -da_grid_x 20 -da_grid_y 20 -lidvelocity 100 -grashof 1e3 -ts_max_steps 100 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type rosw -ts_rosw_type ra3pw -ts_monitor -ts_monitor_solution_vtk 'foo-%03d.vts' -ts_monitor_solution_vtk_interval -1
407: requires: !complex !single
409: test:
410: suffix: 3
411: nsize: 4
412: args: -da_refine 2 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 -pc_type none -ts_type beuler -ts_monitor -snes_monitor -snes_type aspin -da_overlap 4 -npc_sub_ksp_type preonly -npc_sub_pc_type lu
413: requires: !complex !single
415: test:
416: suffix: 4
417: nsize: 2
418: args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3
419: requires: !complex !single
421: test:
422: suffix: asm
423: nsize: 4
424: args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3
425: requires: !complex !single
427: TEST*/