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.
11
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
15
16 Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient
17 Continuation and Differential-Algebraic Equations, 2003.
18
19 That example was removed because it used obsolete interfaces, but the
20 algorithms from the paper can be reproduced using this example.
21
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.
29
30 Comment preserved from snes/tutorials/ex27.c, since removed:
31
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.
40
41 See also https://lists.mcs.anl.gov/pipermail/petsc-dev/2010-March/002362.html
42 */
43
44 /* ------------------------------------------------------------------------
45
46 We thank David E. Keyes for contributing the driven cavity discretization
47 within this example code.
48
49 This problem is modeled by the partial differential equation system
50
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
55
56 in the unit square, which is uniformly discretized in each of x and
57 y in this simple encoding.
58
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.
66
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.
71
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).
77
78 ------------------------------------------------------------------------- */
79
80 /*
81 Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
82 Include "petscts.h" so that we can use TS solvers. Note that this
83 file automatically includes:
84 petscsys.h - base PETSc routines petscvec.h - vectors
85 petscmat.h - matrices
86 petscis.h - index sets petscksp.h - Krylov subspace methods
87 petscviewer.h - viewers petscpc.h - preconditioners
88 petscksp.h - linear solvers petscsnes.h - nonlinear solvers
89 */
90 #include <petscts.h>
91 #include <petscdm.h>
92 #include <petscdmda.h>
93
94 /*
95 User-defined routines and data structures
96 */
97 typedef struct {
98 PetscScalar u, v, omega, temp;
99 } Field;
100
101 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, void *);
102
103 typedef struct {
104 PetscReal lidvelocity, prandtl, grashof; /* physical parameters */
105 PetscBool parabolic; /* allow a transient term corresponding roughly to artificial compressibility */
106 PetscReal cfl_initial; /* CFL for first time step */
107 } AppCtx;
108
109 PetscErrorCode FormInitialSolution(TS, Vec, AppCtx *);
110
main(int argc,char ** argv)111 int main(int argc, char **argv)
112 {
113 AppCtx user; /* user-defined work context */
114 PetscInt mx, my, steps;
115 TS ts;
116 DM da;
117 Vec X;
118 PetscReal ftime;
119 TSConvergedReason reason;
120
121 PetscFunctionBeginUser;
122 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
123 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
124 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da));
125 PetscCall(DMSetFromOptions(da));
126 PetscCall(DMSetUp(da));
127 PetscCall(TSSetDM(ts, (DM)da));
128
129 PetscCall(DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
130 /*
131 Problem parameters (velocity of lid, prandtl, and grashof numbers)
132 */
133 user.lidvelocity = 1.0 / (mx * my);
134 user.prandtl = 1.0;
135 user.grashof = 1.0;
136 user.parabolic = PETSC_FALSE;
137 user.cfl_initial = 50.;
138
139 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Driven cavity/natural convection options", "");
140 PetscCall(PetscOptionsReal("-lidvelocity", "Lid velocity, related to Reynolds number", "", user.lidvelocity, &user.lidvelocity, NULL));
141 PetscCall(PetscOptionsReal("-prandtl", "Ratio of viscous to thermal diffusivity", "", user.prandtl, &user.prandtl, NULL));
142 PetscCall(PetscOptionsReal("-grashof", "Ratio of buoyant to viscous forces", "", user.grashof, &user.grashof, NULL));
143 PetscCall(PetscOptionsBool("-parabolic", "Relax incompressibility to make the system parabolic instead of differential-algebraic", "", user.parabolic, &user.parabolic, NULL));
144 PetscCall(PetscOptionsReal("-cfl_initial", "Advective CFL for the first time step", "", user.cfl_initial, &user.cfl_initial, NULL));
145 PetscOptionsEnd();
146
147 PetscCall(DMDASetFieldName(da, 0, "x-velocity"));
148 PetscCall(DMDASetFieldName(da, 1, "y-velocity"));
149 PetscCall(DMDASetFieldName(da, 2, "Omega"));
150 PetscCall(DMDASetFieldName(da, 3, "temperature"));
151
152 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153 Create user context, set problem data, create vector data structures.
154 Also, compute the initial guess.
155 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156
157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158 Create time integration context
159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160 PetscCall(DMSetApplicationContext(da, &user));
161 PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (DMDATSIFunctionLocalFn *)FormIFunctionLocal, &user));
162 PetscCall(TSSetMaxSteps(ts, 10000));
163 PetscCall(TSSetMaxTime(ts, 1e12));
164 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
165 PetscCall(TSSetTimeStep(ts, user.cfl_initial / (user.lidvelocity * mx)));
166 PetscCall(TSSetFromOptions(ts));
167
168 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));
169
170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 Solve the nonlinear system
172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173
174 PetscCall(DMCreateGlobalVector(da, &X));
175 PetscCall(FormInitialSolution(ts, X, &user));
176
177 PetscCall(TSSolve(ts, X));
178 PetscCall(TSGetSolveTime(ts, &ftime));
179 PetscCall(TSGetStepNumber(ts, &steps));
180 PetscCall(TSGetConvergedReason(ts, &reason));
181
182 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));
183
184 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185 Free work space. All PETSc objects should be destroyed when they
186 are no longer needed.
187 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188 PetscCall(VecDestroy(&X));
189 PetscCall(DMDestroy(&da));
190 PetscCall(TSDestroy(&ts));
191
192 PetscCall(PetscFinalize());
193 return 0;
194 }
195
196 /* ------------------------------------------------------------------- */
197
198 /*
199 FormInitialSolution - Forms initial approximation.
200
201 Input Parameters:
202 user - user-defined application context
203 X - vector
204
205 Output Parameter:
206 X - vector
207 */
FormInitialSolution(TS ts,Vec X,AppCtx * user)208 PetscErrorCode FormInitialSolution(TS ts, Vec X, AppCtx *user)
209 {
210 DM da;
211 PetscInt i, j, mx, xs, ys, xm, ym;
212 PetscReal grashof, dx;
213 Field **x;
214
215 PetscFunctionBeginUser;
216 grashof = user->grashof;
217 PetscCall(TSGetDM(ts, &da));
218 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
219 dx = 1.0 / (mx - 1);
220
221 /*
222 Get local grid boundaries (for 2-dimensional DMDA):
223 xs, ys - starting grid indices (no ghost points)
224 xm, ym - widths of local grid (no ghost points)
225 */
226 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
227
228 /*
229 Get a pointer to vector data.
230 - For default PETSc vectors, VecGetArray() returns a pointer to
231 the data array. Otherwise, the routine is implementation dependent.
232 - You MUST call VecRestoreArray() when you no longer need access to
233 the array.
234 */
235 PetscCall(DMDAVecGetArray(da, X, &x));
236
237 /*
238 Compute initial guess over the locally owned part of the grid
239 Initial condition is motionless fluid and equilibrium temperature
240 */
241 for (j = ys; j < ys + ym; j++) {
242 for (i = xs; i < xs + xm; i++) {
243 x[j][i].u = 0.0;
244 x[j][i].v = 0.0;
245 x[j][i].omega = 0.0;
246 x[j][i].temp = (grashof > 0) * i * dx;
247 }
248 }
249
250 /*
251 Restore vector
252 */
253 PetscCall(DMDAVecRestoreArray(da, X, &x));
254 PetscFunctionReturn(PETSC_SUCCESS);
255 }
256
FormIFunctionLocal(DMDALocalInfo * info,PetscReal ptime,Field ** x,Field ** xdot,Field ** f,void * ptr)257 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal ptime, Field **x, Field **xdot, Field **f, void *ptr)
258 {
259 AppCtx *user = (AppCtx *)ptr;
260 PetscInt xints, xinte, yints, yinte, i, j;
261 PetscReal hx, hy, dhx, dhy, hxdhy, hydhx;
262 PetscReal grashof, prandtl, lid;
263 PetscScalar u, udot, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym;
264
265 PetscFunctionBeginUser;
266 grashof = user->grashof;
267 prandtl = user->prandtl;
268 lid = user->lidvelocity;
269
270 /*
271 Define mesh intervals ratios for uniform grid.
272
273 Note: FD formulae below are normalized by multiplying through by
274 local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
275
276 */
277 dhx = (PetscReal)(info->mx - 1);
278 dhy = (PetscReal)(info->my - 1);
279 hx = 1.0 / dhx;
280 hy = 1.0 / dhy;
281 hxdhy = hx * dhy;
282 hydhx = hy * dhx;
283
284 xints = info->xs;
285 xinte = info->xs + info->xm;
286 yints = info->ys;
287 yinte = info->ys + info->ym;
288
289 /* Test whether we are on the bottom edge of the global array */
290 if (yints == 0) {
291 j = 0;
292 yints = yints + 1;
293 /* bottom edge */
294 for (i = info->xs; i < info->xs + info->xm; i++) {
295 f[j][i].u = x[j][i].u;
296 f[j][i].v = x[j][i].v;
297 f[j][i].omega = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy;
298 f[j][i].temp = x[j][i].temp - x[j + 1][i].temp;
299 }
300 }
301
302 /* Test whether we are on the top edge of the global array */
303 if (yinte == info->my) {
304 j = info->my - 1;
305 yinte = yinte - 1;
306 /* top edge */
307 for (i = info->xs; i < info->xs + info->xm; i++) {
308 f[j][i].u = x[j][i].u - lid;
309 f[j][i].v = x[j][i].v;
310 f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy;
311 f[j][i].temp = x[j][i].temp - x[j - 1][i].temp;
312 }
313 }
314
315 /* Test whether we are on the left edge of the global array */
316 if (xints == 0) {
317 i = 0;
318 xints = xints + 1;
319 /* left edge */
320 for (j = info->ys; j < info->ys + info->ym; j++) {
321 f[j][i].u = x[j][i].u;
322 f[j][i].v = x[j][i].v;
323 f[j][i].omega = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx;
324 f[j][i].temp = x[j][i].temp;
325 }
326 }
327
328 /* Test whether we are on the right edge of the global array */
329 if (xinte == info->mx) {
330 i = info->mx - 1;
331 xinte = xinte - 1;
332 /* right edge */
333 for (j = info->ys; j < info->ys + info->ym; j++) {
334 f[j][i].u = x[j][i].u;
335 f[j][i].v = x[j][i].v;
336 f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx;
337 f[j][i].temp = x[j][i].temp - (PetscReal)(grashof > 0);
338 }
339 }
340
341 /* Compute over the interior points */
342 for (j = yints; j < yinte; j++) {
343 for (i = xints; i < xinte; i++) {
344 /*
345 convective coefficients for upwinding
346 */
347 vx = x[j][i].u;
348 avx = PetscAbsScalar(vx);
349 vxp = .5 * (vx + avx);
350 vxm = .5 * (vx - avx);
351 vy = x[j][i].v;
352 avy = PetscAbsScalar(vy);
353 vyp = .5 * (vy + avy);
354 vym = .5 * (vy - avy);
355
356 /* U velocity */
357 u = x[j][i].u;
358 udot = user->parabolic ? xdot[j][i].u : 0.;
359 uxx = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
360 uyy = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
361 f[j][i].u = udot + uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx;
362
363 /* V velocity */
364 u = x[j][i].v;
365 udot = user->parabolic ? xdot[j][i].v : 0.;
366 uxx = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx;
367 uyy = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy;
368 f[j][i].v = udot + uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy;
369
370 /* Omega */
371 u = x[j][i].omega;
372 uxx = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx;
373 uyy = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy;
374 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);
375
376 /* Temperature */
377 u = x[j][i].temp;
378 uxx = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx;
379 uyy = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy;
380 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));
381 }
382 }
383
384 /*
385 Flop count (multiply-adds are counted as 2 operations)
386 */
387 PetscCall(PetscLogFlops(84.0 * info->ym * info->xm));
388 PetscFunctionReturn(PETSC_SUCCESS);
389 }
390
391 /*TEST
392
393 test:
394 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'
395 requires: !complex !single
396
397 test:
398 suffix: 1_steps
399 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'
400 requires: !complex !single
401
402 test:
403 suffix: 2
404 nsize: 4
405 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
406 requires: !complex !single
407
408 test:
409 suffix: 3
410 nsize: 4
411 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_short -snes_type aspin -da_overlap 4 -npc_sub_ksp_type preonly -npc_sub_pc_type lu
412 requires: !complex !single
413
414 test:
415 suffix: 4
416 nsize: 2
417 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3
418 requires: !complex !single
419
420 test:
421 suffix: asm
422 nsize: 4
423 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3
424 requires: !complex !single
425
426 TEST*/
427