Actual source code: ex26.c
petsc-3.7.5 2017-01-01
2: static char help[] = "Transient nonlinear driven cavity in 2d.\n\
3: \n\
4: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
5: The flow can be driven with the lid or with bouyancy or both:\n\
6: -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
7: -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
8: -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
9: -contours : draw contour plots of solution\n\n";
10: /*
11: See src/snes/examples/tutorials/ex19.c for the steady-state version.
13: See src/snes/examples/tutorials/ex27.c for a version written specifically
14: for pseudo-transient continuation, without using TS.
15: */
17: /*T
18: Concepts: TS^solving a system of nonlinear equations (parallel multicomponent example);
19: Concepts: DMDA^using distributed arrays;
20: Concepts: TS^multicomponent
21: Concepts: TS^differential-algebraic equation
22: Processors: n
23: T*/
24: /* ------------------------------------------------------------------------
26: We thank David E. Keyes for contributing the driven cavity discretization
27: within this example code.
29: This problem is modeled by the partial differential equation system
31: - Lap(U) - Grad_y(Omega) = 0
32: - Lap(V) + Grad_x(Omega) = 0
33: Omega_t - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
34: T_t - Lap(T) + PR*Div([U*T,V*T]) = 0
36: in the unit square, which is uniformly discretized in each of x and
37: y in this simple encoding.
39: No-slip, rigid-wall Dirichlet conditions are used for [U,V].
40: Dirichlet conditions are used for Omega, based on the definition of
41: vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
42: constant coordinate boundary, the tangential derivative is zero.
43: Dirichlet conditions are used for T on the left and right walls,
44: and insulation homogeneous Neumann conditions are used for T on
45: the top and bottom walls.
47: A finite difference approximation with the usual 5-point stencil
48: is used to discretize the boundary value problem to obtain a
49: nonlinear system of equations. Upwinding is used for the divergence
50: (convective) terms and central for the gradient (source) terms.
52: The Jacobian can be either
53: * formed via finite differencing using coloring (the default), or
54: * applied matrix-free via the option -snes_mf
55: (for larger grid problems this variant may not converge
56: without a preconditioner due to ill-conditioning).
58: ------------------------------------------------------------------------- */
60: /*
61: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
62: Include "petscts.h" so that we can use TS solvers. Note that this
63: file automatically includes:
64: petscsys.h - base PETSc routines petscvec.h - vectors
65: petscmat.h - matrices
66: petscis.h - index sets petscksp.h - Krylov subspace methods
67: petscviewer.h - viewers petscpc.h - preconditioners
68: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
69: */
70: #include <petscts.h>
71: #include <petscdm.h>
72: #include <petscdmda.h>
74: /*
75: User-defined routines and data structures
76: */
77: typedef struct {
78: PetscScalar u,v,omega,temp;
79: } Field;
81: PetscErrorCode FormIFunctionLocal(DMDALocalInfo*,PetscReal,Field**,Field**,Field**,void*);
83: typedef struct {
84: PetscReal lidvelocity,prandtl,grashof; /* physical parameters */
85: PetscBool parabolic; /* allow a transient term corresponding roughly to artificial compressibility */
86: PetscReal cfl_initial; /* CFL for first time step */
87: } AppCtx;
89: PetscErrorCode FormInitialSolution(TS,Vec,AppCtx*);
93: int main(int argc,char **argv)
94: {
95: AppCtx user; /* user-defined work context */
96: PetscInt mx,my,steps;
97: PetscErrorCode ierr;
98: TS ts;
99: DM da;
100: Vec X;
101: PetscReal ftime;
102: TSConvergedReason reason;
104: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return(1);
106: TSCreate(PETSC_COMM_WORLD,&ts);
108: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
109: TSSetDM(ts,(DM)da);
111: DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
112: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
113: /*
114: Problem parameters (velocity of lid, prandtl, and grashof numbers)
115: */
116: user.lidvelocity = 1.0/(mx*my);
117: user.prandtl = 1.0;
118: user.grashof = 1.0;
119: user.parabolic = PETSC_FALSE;
120: user.cfl_initial = 50.;
122: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Driven cavity/natural convection options","");
123: PetscOptionsReal("-lidvelocity","Lid velocity, related to Reynolds number","",user.lidvelocity,&user.lidvelocity,NULL);
124: PetscOptionsReal("-prandtl","Ratio of viscous to thermal diffusivity","",user.prandtl,&user.prandtl,NULL);
125: PetscOptionsReal("-grashof","Ratio of bouyant to viscous forces","",user.grashof,&user.grashof,NULL);
126: PetscOptionsBool("-parabolic","Relax incompressibility to make the system parabolic instead of differential-algebraic","",user.parabolic,&user.parabolic,NULL);
127: PetscOptionsReal("-cfl_initial","Advective CFL for the first time step","",user.cfl_initial,&user.cfl_initial,NULL);
128: PetscOptionsEnd();
130: DMDASetFieldName(da,0,"x-velocity");
131: DMDASetFieldName(da,1,"y-velocity");
132: DMDASetFieldName(da,2,"Omega");
133: DMDASetFieldName(da,3,"temperature");
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Create user context, set problem data, create vector data structures.
137: Also, compute the initial guess.
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Create time integration context
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: DMSetApplicationContext(da,&user);
144: DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)FormIFunctionLocal,&user);
145: TSSetDuration(ts,10000,1e12);
146: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
147: TSSetInitialTimeStep(ts,0.0,user.cfl_initial/(user.lidvelocity*mx));
148: TSSetFromOptions(ts);
150: PetscPrintf(PETSC_COMM_WORLD,"%Dx%D grid, lid velocity = %g, prandtl # = %g, grashof # = %g\n",mx,my,(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof);
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Solve the nonlinear system
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: DMCreateGlobalVector(da,&X);
158: FormInitialSolution(ts,X,&user);
160: TSSolve(ts,X);
161: TSGetSolveTime(ts,&ftime);
162: TSGetTimeStepNumber(ts,&steps);
163: TSGetConvergedReason(ts,&reason);
165: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Free work space. All PETSc objects should be destroyed when they
169: are no longer needed.
170: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171: VecDestroy(&X);
172: DMDestroy(&da);
173: TSDestroy(&ts);
175: PetscFinalize();
176: return 0;
177: }
179: /* ------------------------------------------------------------------- */
184: /*
185: FormInitialSolution - Forms initial approximation.
187: Input Parameters:
188: user - user-defined application context
189: X - vector
191: Output Parameter:
192: X - vector
193: */
194: PetscErrorCode FormInitialSolution(TS ts,Vec X,AppCtx *user)
195: {
196: DM da;
197: PetscInt i,j,mx,xs,ys,xm,ym;
199: PetscReal grashof,dx;
200: Field **x;
202: grashof = user->grashof;
203: TSGetDM(ts,&da);
204: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
205: dx = 1.0/(mx-1);
207: /*
208: Get local grid boundaries (for 2-dimensional DMDA):
209: xs, ys - starting grid indices (no ghost points)
210: xm, ym - widths of local grid (no ghost points)
211: */
212: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
214: /*
215: Get a pointer to vector data.
216: - For default PETSc vectors, VecGetArray() returns a pointer to
217: the data array. Otherwise, the routine is implementation dependent.
218: - You MUST call VecRestoreArray() when you no longer need access to
219: the array.
220: */
221: DMDAVecGetArray(da,X,&x);
223: /*
224: Compute initial guess over the locally owned part of the grid
225: Initial condition is motionless fluid and equilibrium temperature
226: */
227: for (j=ys; j<ys+ym; j++) {
228: for (i=xs; i<xs+xm; i++) {
229: x[j][i].u = 0.0;
230: x[j][i].v = 0.0;
231: x[j][i].omega = 0.0;
232: x[j][i].temp = (grashof>0)*i*dx;
233: }
234: }
236: /*
237: Restore vector
238: */
239: DMDAVecRestoreArray(da,X,&x);
240: return 0;
241: }
245: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,Field **x,Field **xdot,Field **f,void *ptr)
246: {
247: AppCtx *user = (AppCtx*)ptr;
249: PetscInt xints,xinte,yints,yinte,i,j;
250: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
251: PetscReal grashof,prandtl,lid;
252: PetscScalar u,udot,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
255: grashof = user->grashof;
256: prandtl = user->prandtl;
257: lid = user->lidvelocity;
259: /*
260: Define mesh intervals ratios for uniform grid.
262: Note: FD formulae below are normalized by multiplying through by
263: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
266: */
267: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
268: hx = 1.0/dhx; hy = 1.0/dhy;
269: hxdhy = hx*dhy; hydhx = hy*dhx;
271: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
273: /* Test whether we are on the bottom edge of the global array */
274: if (yints == 0) {
275: j = 0;
276: yints = yints + 1;
277: /* bottom edge */
278: for (i=info->xs; i<info->xs+info->xm; i++) {
279: f[j][i].u = x[j][i].u;
280: f[j][i].v = x[j][i].v;
281: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
282: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
283: }
284: }
286: /* Test whether we are on the top edge of the global array */
287: if (yinte == info->my) {
288: j = info->my - 1;
289: yinte = yinte - 1;
290: /* top edge */
291: for (i=info->xs; i<info->xs+info->xm; i++) {
292: f[j][i].u = x[j][i].u - lid;
293: f[j][i].v = x[j][i].v;
294: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
295: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
296: }
297: }
299: /* Test whether we are on the left edge of the global array */
300: if (xints == 0) {
301: i = 0;
302: xints = xints + 1;
303: /* left edge */
304: for (j=info->ys; j<info->ys+info->ym; j++) {
305: f[j][i].u = x[j][i].u;
306: f[j][i].v = x[j][i].v;
307: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
308: f[j][i].temp = x[j][i].temp;
309: }
310: }
312: /* Test whether we are on the right edge of the global array */
313: if (xinte == info->mx) {
314: i = info->mx - 1;
315: xinte = xinte - 1;
316: /* right edge */
317: for (j=info->ys; j<info->ys+info->ym; j++) {
318: f[j][i].u = x[j][i].u;
319: f[j][i].v = x[j][i].v;
320: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
321: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
322: }
323: }
325: /* Compute over the interior points */
326: for (j=yints; j<yinte; j++) {
327: for (i=xints; i<xinte; i++) {
329: /*
330: convective coefficients for upwinding
331: */
332: vx = x[j][i].u; avx = PetscAbsScalar(vx);
333: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
334: vy = x[j][i].v; avy = PetscAbsScalar(vy);
335: vyp = .5*(vy+avy); vym = .5*(vy-avy);
337: /* U velocity */
338: u = x[j][i].u;
339: udot = user->parabolic ? xdot[j][i].u : 0.;
340: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
341: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
342: f[j][i].u = udot + uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
344: /* V velocity */
345: u = x[j][i].v;
346: udot = user->parabolic ? xdot[j][i].v : 0.;
347: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
348: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
349: f[j][i].v = udot + uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
351: /* Omega */
352: u = x[j][i].omega;
353: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
354: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
355: f[j][i].omega = (xdot[j][i].omega + uxx + uyy
356: + (vxp*(u - x[j][i-1].omega)
357: + vxm*(x[j][i+1].omega - u)) * hy
358: + (vyp*(u - x[j-1][i].omega)
359: + vym*(x[j+1][i].omega - u)) * hx
360: - .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy);
362: /* Temperature */
363: u = x[j][i].temp;
364: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
365: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
366: f[j][i].temp = (xdot[j][i].temp + uxx + uyy
367: + prandtl * ((vxp*(u - x[j][i-1].temp)
368: + vxm*(x[j][i+1].temp - u)) * hy
369: + (vyp*(u - x[j-1][i].temp)
370: + vym*(x[j+1][i].temp - u)) * hx));
371: }
372: }
374: /*
375: Flop count (multiply-adds are counted as 2 operations)
376: */
377: PetscLogFlops(84.0*info->ym*info->xm);
378: return(0);
379: }