Actual source code: ex26.c
petsc-3.6.4 2016-04-12
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: PassiveReal 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: TSSetInitialTimeStep(ts,0.0,user.cfl_initial/(user.lidvelocity*mx));
147: TSSetFromOptions(ts);
149: 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);
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Solve the nonlinear system
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: DMCreateGlobalVector(da,&X);
157: FormInitialSolution(ts,X,&user);
159: TSSolve(ts,X);
160: TSGetSolveTime(ts,&ftime);
161: TSGetTimeStepNumber(ts,&steps);
162: TSGetConvergedReason(ts,&reason);
164: PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Free work space. All PETSc objects should be destroyed when they
168: are no longer needed.
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: VecDestroy(&X);
171: DMDestroy(&da);
172: TSDestroy(&ts);
174: PetscFinalize();
175: return 0;
176: }
178: /* ------------------------------------------------------------------- */
183: /*
184: FormInitialSolution - Forms initial approximation.
186: Input Parameters:
187: user - user-defined application context
188: X - vector
190: Output Parameter:
191: X - vector
192: */
193: PetscErrorCode FormInitialSolution(TS ts,Vec X,AppCtx *user)
194: {
195: DM da;
196: PetscInt i,j,mx,xs,ys,xm,ym;
198: PetscReal grashof,dx;
199: Field **x;
201: grashof = user->grashof;
202: TSGetDM(ts,&da);
203: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
204: dx = 1.0/(mx-1);
206: /*
207: Get local grid boundaries (for 2-dimensional DMDA):
208: xs, ys - starting grid indices (no ghost points)
209: xm, ym - widths of local grid (no ghost points)
210: */
211: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
213: /*
214: Get a pointer to vector data.
215: - For default PETSc vectors, VecGetArray() returns a pointer to
216: the data array. Otherwise, the routine is implementation dependent.
217: - You MUST call VecRestoreArray() when you no longer need access to
218: the array.
219: */
220: DMDAVecGetArray(da,X,&x);
222: /*
223: Compute initial guess over the locally owned part of the grid
224: Initial condition is motionless fluid and equilibrium temperature
225: */
226: for (j=ys; j<ys+ym; j++) {
227: for (i=xs; i<xs+xm; i++) {
228: x[j][i].u = 0.0;
229: x[j][i].v = 0.0;
230: x[j][i].omega = 0.0;
231: x[j][i].temp = (grashof>0)*i*dx;
232: }
233: }
235: /*
236: Restore vector
237: */
238: DMDAVecRestoreArray(da,X,&x);
239: return 0;
240: }
244: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,Field **x,Field **xdot,Field **f,void *ptr)
245: {
246: AppCtx *user = (AppCtx*)ptr;
248: PetscInt xints,xinte,yints,yinte,i,j;
249: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
250: PetscReal grashof,prandtl,lid;
251: PetscScalar u,udot,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
254: grashof = user->grashof;
255: prandtl = user->prandtl;
256: lid = user->lidvelocity;
258: /*
259: Define mesh intervals ratios for uniform grid.
261: Note: FD formulae below are normalized by multiplying through by
262: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
265: */
266: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
267: hx = 1.0/dhx; hy = 1.0/dhy;
268: hxdhy = hx*dhy; hydhx = hy*dhx;
270: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
272: /* Test whether we are on the bottom edge of the global array */
273: if (yints == 0) {
274: j = 0;
275: yints = yints + 1;
276: /* bottom edge */
277: for (i=info->xs; i<info->xs+info->xm; i++) {
278: f[j][i].u = x[j][i].u;
279: f[j][i].v = x[j][i].v;
280: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
281: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
282: }
283: }
285: /* Test whether we are on the top edge of the global array */
286: if (yinte == info->my) {
287: j = info->my - 1;
288: yinte = yinte - 1;
289: /* top edge */
290: for (i=info->xs; i<info->xs+info->xm; i++) {
291: f[j][i].u = x[j][i].u - lid;
292: f[j][i].v = x[j][i].v;
293: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
294: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
295: }
296: }
298: /* Test whether we are on the left edge of the global array */
299: if (xints == 0) {
300: i = 0;
301: xints = xints + 1;
302: /* left edge */
303: for (j=info->ys; j<info->ys+info->ym; j++) {
304: f[j][i].u = x[j][i].u;
305: f[j][i].v = x[j][i].v;
306: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
307: f[j][i].temp = x[j][i].temp;
308: }
309: }
311: /* Test whether we are on the right edge of the global array */
312: if (xinte == info->mx) {
313: i = info->mx - 1;
314: xinte = xinte - 1;
315: /* right edge */
316: for (j=info->ys; j<info->ys+info->ym; j++) {
317: f[j][i].u = x[j][i].u;
318: f[j][i].v = x[j][i].v;
319: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
320: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
321: }
322: }
324: /* Compute over the interior points */
325: for (j=yints; j<yinte; j++) {
326: for (i=xints; i<xinte; i++) {
328: /*
329: convective coefficients for upwinding
330: */
331: vx = x[j][i].u; avx = PetscAbsScalar(vx);
332: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
333: vy = x[j][i].v; avy = PetscAbsScalar(vy);
334: vyp = .5*(vy+avy); vym = .5*(vy-avy);
336: /* U velocity */
337: u = x[j][i].u;
338: udot = user->parabolic ? xdot[j][i].u : 0.;
339: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
340: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
341: f[j][i].u = udot + uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
343: /* V velocity */
344: u = x[j][i].v;
345: udot = user->parabolic ? xdot[j][i].v : 0.;
346: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
347: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
348: f[j][i].v = udot + uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
350: /* Omega */
351: u = x[j][i].omega;
352: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
353: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
354: f[j][i].omega = (xdot[j][i].omega + uxx + uyy
355: + (vxp*(u - x[j][i-1].omega)
356: + vxm*(x[j][i+1].omega - u)) * hy
357: + (vyp*(u - x[j-1][i].omega)
358: + vym*(x[j+1][i].omega - u)) * hx
359: - .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy);
361: /* Temperature */
362: u = x[j][i].temp;
363: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
364: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
365: f[j][i].temp = (xdot[j][i].temp + uxx + uyy
366: + prandtl * ((vxp*(u - x[j][i-1].temp)
367: + vxm*(x[j][i+1].temp - u)) * hy
368: + (vyp*(u - x[j-1][i].temp)
369: + vym*(x[j+1][i].temp - u)) * hx));
370: }
371: }
373: /*
374: Flop count (multiply-adds are counted as 2 operations)
375: */
376: PetscLogFlops(84.0*info->ym*info->xm);
377: return(0);
378: }