Actual source code: ex26.c
petsc-3.4.5 2014-06-29
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 <petscdmda.h>
73: /*
74: User-defined routines and data structures
75: */
76: typedef struct {
77: PetscScalar u,v,omega,temp;
78: } Field;
80: PetscErrorCode FormIFunctionLocal(DMDALocalInfo*,PetscReal,Field**,Field**,Field**,void*);
82: typedef struct {
83: PassiveReal lidvelocity,prandtl,grashof; /* physical parameters */
84: PetscBool parabolic; /* allow a transient term corresponding roughly to artificial compressibility */
85: PetscReal cfl_initial; /* CFL for first time step */
86: } AppCtx;
88: PetscErrorCode FormInitialSolution(TS,Vec,AppCtx*);
92: int main(int argc,char **argv)
93: {
94: AppCtx user; /* user-defined work context */
95: PetscInt mx,my,steps;
96: PetscErrorCode ierr;
97: TS ts;
98: DM da;
99: Vec X;
100: PetscReal ftime;
101: TSConvergedReason reason;
103: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return(1);
105: TSCreate(PETSC_COMM_WORLD,&ts);
107: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
108: TSSetDM(ts,(DM)da);
110: DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
111: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
112: /*
113: Problem parameters (velocity of lid, prandtl, and grashof numbers)
114: */
115: user.lidvelocity = 1.0/(mx*my);
116: user.prandtl = 1.0;
117: user.grashof = 1.0;
118: user.parabolic = PETSC_FALSE;
119: user.cfl_initial = 50.;
121: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Driven cavity/natural convection options","");
122: PetscOptionsReal("-lidvelocity","Lid velocity, related to Reynolds number","",user.lidvelocity,&user.lidvelocity,NULL);
123: PetscOptionsReal("-prandtl","Ratio of viscous to thermal diffusivity","",user.prandtl,&user.prandtl,NULL);
124: PetscOptionsReal("-grashof","Ratio of bouyant to viscous forces","",user.grashof,&user.grashof,NULL);
125: PetscOptionsBool("-parabolic","Relax incompressibility to make the system parabolic instead of differential-algebraic","",user.parabolic,&user.parabolic,NULL);
126: PetscOptionsReal("-cfl_initial","Advective CFL for the first time step","",user.cfl_initial,&user.cfl_initial,NULL);
127: PetscOptionsEnd();
129: DMDASetFieldName(da,0,"x-velocity");
130: DMDASetFieldName(da,1,"y-velocity");
131: DMDASetFieldName(da,2,"Omega");
132: DMDASetFieldName(da,3,"temperature");
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Create user context, set problem data, create vector data structures.
136: Also, compute the initial guess.
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Create time integration context
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: DMSetApplicationContext(da,&user);
143: DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)FormIFunctionLocal,&user);
144: TSSetDuration(ts,10000,1e12);
145: TSSetInitialTimeStep(ts,0.0,user.cfl_initial/(user.lidvelocity*mx));
146: TSSetFromOptions(ts);
148: PetscPrintf(PETSC_COMM_WORLD,"%Dx%D grid, lid velocity = %G, prandtl # = %G, grashof # = %G\n",mx,my,user.lidvelocity,user.prandtl,user.grashof);
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Solve the nonlinear system
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: DMCreateGlobalVector(da,&X);
156: FormInitialSolution(ts,X,&user);
158: TSSolve(ts,X);
159: TSGetSolveTime(ts,&ftime);
160: TSGetTimeStepNumber(ts,&steps);
161: TSGetConvergedReason(ts,&reason);
163: PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,steps);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Free work space. All PETSc objects should be destroyed when they
167: are no longer needed.
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: VecDestroy(&X);
170: DMDestroy(&da);
171: TSDestroy(&ts);
173: PetscFinalize();
174: return 0;
175: }
177: /* ------------------------------------------------------------------- */
182: /*
183: FormInitialSolution - Forms initial approximation.
185: Input Parameters:
186: user - user-defined application context
187: X - vector
189: Output Parameter:
190: X - vector
191: */
192: PetscErrorCode FormInitialSolution(TS ts,Vec X,AppCtx *user)
193: {
194: DM da;
195: PetscInt i,j,mx,xs,ys,xm,ym;
197: PetscReal grashof,dx;
198: Field **x;
200: grashof = user->grashof;
201: TSGetDM(ts,&da);
202: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
203: dx = 1.0/(mx-1);
205: /*
206: Get local grid boundaries (for 2-dimensional DMDA):
207: xs, ys - starting grid indices (no ghost points)
208: xm, ym - widths of local grid (no ghost points)
209: */
210: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
212: /*
213: Get a pointer to vector data.
214: - For default PETSc vectors, VecGetArray() returns a pointer to
215: the data array. Otherwise, the routine is implementation dependent.
216: - You MUST call VecRestoreArray() when you no longer need access to
217: the array.
218: */
219: DMDAVecGetArray(da,X,&x);
221: /*
222: Compute initial guess over the locally owned part of the grid
223: Initial condition is motionless fluid and equilibrium temperature
224: */
225: for (j=ys; j<ys+ym; j++) {
226: for (i=xs; i<xs+xm; i++) {
227: x[j][i].u = 0.0;
228: x[j][i].v = 0.0;
229: x[j][i].omega = 0.0;
230: x[j][i].temp = (grashof>0)*i*dx;
231: }
232: }
234: /*
235: Restore vector
236: */
237: DMDAVecRestoreArray(da,X,&x);
238: return 0;
239: }
243: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,Field **x,Field **xdot,Field **f,void *ptr)
244: {
245: AppCtx *user = (AppCtx*)ptr;
247: PetscInt xints,xinte,yints,yinte,i,j;
248: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
249: PetscReal grashof,prandtl,lid;
250: PetscScalar u,udot,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
253: grashof = user->grashof;
254: prandtl = user->prandtl;
255: lid = user->lidvelocity;
257: /*
258: Define mesh intervals ratios for uniform grid.
260: Note: FD formulae below are normalized by multiplying through by
261: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
264: */
265: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
266: hx = 1.0/dhx; hy = 1.0/dhy;
267: hxdhy = hx*dhy; hydhx = hy*dhx;
269: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
271: /* Test whether we are on the bottom edge of the global array */
272: if (yints == 0) {
273: j = 0;
274: yints = yints + 1;
275: /* bottom edge */
276: for (i=info->xs; i<info->xs+info->xm; i++) {
277: f[j][i].u = x[j][i].u;
278: f[j][i].v = x[j][i].v;
279: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
280: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
281: }
282: }
284: /* Test whether we are on the top edge of the global array */
285: if (yinte == info->my) {
286: j = info->my - 1;
287: yinte = yinte - 1;
288: /* top edge */
289: for (i=info->xs; i<info->xs+info->xm; i++) {
290: f[j][i].u = x[j][i].u - lid;
291: f[j][i].v = x[j][i].v;
292: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
293: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
294: }
295: }
297: /* Test whether we are on the left edge of the global array */
298: if (xints == 0) {
299: i = 0;
300: xints = xints + 1;
301: /* left edge */
302: for (j=info->ys; j<info->ys+info->ym; j++) {
303: f[j][i].u = x[j][i].u;
304: f[j][i].v = x[j][i].v;
305: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
306: f[j][i].temp = x[j][i].temp;
307: }
308: }
310: /* Test whether we are on the right edge of the global array */
311: if (xinte == info->mx) {
312: i = info->mx - 1;
313: xinte = xinte - 1;
314: /* right edge */
315: for (j=info->ys; j<info->ys+info->ym; j++) {
316: f[j][i].u = x[j][i].u;
317: f[j][i].v = x[j][i].v;
318: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
319: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
320: }
321: }
323: /* Compute over the interior points */
324: for (j=yints; j<yinte; j++) {
325: for (i=xints; i<xinte; i++) {
327: /*
328: convective coefficients for upwinding
329: */
330: vx = x[j][i].u; avx = PetscAbsScalar(vx);
331: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
332: vy = x[j][i].v; avy = PetscAbsScalar(vy);
333: vyp = .5*(vy+avy); vym = .5*(vy-avy);
335: /* U velocity */
336: u = x[j][i].u;
337: udot = user->parabolic ? xdot[j][i].u : 0.;
338: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
339: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
340: f[j][i].u = udot + uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
342: /* V velocity */
343: u = x[j][i].v;
344: udot = user->parabolic ? xdot[j][i].v : 0.;
345: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
346: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
347: f[j][i].v = udot + uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
349: /* Omega */
350: u = x[j][i].omega;
351: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
352: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
353: f[j][i].omega = (xdot[j][i].omega + uxx + uyy
354: + (vxp*(u - x[j][i-1].omega)
355: + vxm*(x[j][i+1].omega - u)) * hy
356: + (vyp*(u - x[j-1][i].omega)
357: + vym*(x[j+1][i].omega - u)) * hx
358: - .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy);
360: /* Temperature */
361: u = x[j][i].temp;
362: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
363: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
364: f[j][i].temp = (xdot[j][i].temp + uxx + uyy
365: + prandtl * ((vxp*(u - x[j][i-1].temp)
366: + vxm*(x[j][i+1].temp - u)) * hy
367: + (vyp*(u - x[j-1][i].temp)
368: + vym*(x[j+1][i].temp - u)) * hx));
369: }
370: }
372: /*
373: Flop count (multiply-adds are counted as 2 operations)
374: */
375: PetscLogFlops(84.0*info->ym*info->xm);
376: return(0);
377: }