Actual source code: ex26.c
petsc-3.8.4 2018-03-24
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*);
91: int main(int argc,char **argv)
92: {
93: AppCtx user; /* user-defined work context */
94: PetscInt mx,my,steps;
95: PetscErrorCode ierr;
96: TS ts;
97: DM da;
98: Vec X;
99: PetscReal ftime;
100: TSConvergedReason reason;
102: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
103: TSCreate(PETSC_COMM_WORLD,&ts);
104: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
105: DMSetFromOptions(da);
106: DMSetUp(da);
107: TSSetDM(ts,(DM)da);
109: DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
110: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
111: /*
112: Problem parameters (velocity of lid, prandtl, and grashof numbers)
113: */
114: user.lidvelocity = 1.0/(mx*my);
115: user.prandtl = 1.0;
116: user.grashof = 1.0;
117: user.parabolic = PETSC_FALSE;
118: user.cfl_initial = 50.;
120: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Driven cavity/natural convection options","");
121: PetscOptionsReal("-lidvelocity","Lid velocity, related to Reynolds number","",user.lidvelocity,&user.lidvelocity,NULL);
122: PetscOptionsReal("-prandtl","Ratio of viscous to thermal diffusivity","",user.prandtl,&user.prandtl,NULL);
123: PetscOptionsReal("-grashof","Ratio of bouyant to viscous forces","",user.grashof,&user.grashof,NULL);
124: PetscOptionsBool("-parabolic","Relax incompressibility to make the system parabolic instead of differential-algebraic","",user.parabolic,&user.parabolic,NULL);
125: PetscOptionsReal("-cfl_initial","Advective CFL for the first time step","",user.cfl_initial,&user.cfl_initial,NULL);
126: PetscOptionsEnd();
128: DMDASetFieldName(da,0,"x-velocity");
129: DMDASetFieldName(da,1,"y-velocity");
130: DMDASetFieldName(da,2,"Omega");
131: DMDASetFieldName(da,3,"temperature");
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Create user context, set problem data, create vector data structures.
135: Also, compute the initial guess.
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Create time integration context
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: DMSetApplicationContext(da,&user);
142: DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)FormIFunctionLocal,&user);
143: TSSetMaxSteps(ts,10000);
144: TSSetMaxTime(ts,1e12);
145: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
146: TSSetTimeStep(ts,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: TSGetStepNumber(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 ierr;
176: }
178: /* ------------------------------------------------------------------- */
181: /*
182: FormInitialSolution - Forms initial approximation.
184: Input Parameters:
185: user - user-defined application context
186: X - vector
188: Output Parameter:
189: X - vector
190: */
191: PetscErrorCode FormInitialSolution(TS ts,Vec X,AppCtx *user)
192: {
193: DM da;
194: PetscInt i,j,mx,xs,ys,xm,ym;
196: PetscReal grashof,dx;
197: Field **x;
199: grashof = user->grashof;
200: TSGetDM(ts,&da);
201: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
202: dx = 1.0/(mx-1);
204: /*
205: Get local grid boundaries (for 2-dimensional DMDA):
206: xs, ys - starting grid indices (no ghost points)
207: xm, ym - widths of local grid (no ghost points)
208: */
209: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
211: /*
212: Get a pointer to vector data.
213: - For default PETSc vectors, VecGetArray() returns a pointer to
214: the data array. Otherwise, the routine is implementation dependent.
215: - You MUST call VecRestoreArray() when you no longer need access to
216: the array.
217: */
218: DMDAVecGetArray(da,X,&x);
220: /*
221: Compute initial guess over the locally owned part of the grid
222: Initial condition is motionless fluid and equilibrium temperature
223: */
224: for (j=ys; j<ys+ym; j++) {
225: for (i=xs; i<xs+xm; i++) {
226: x[j][i].u = 0.0;
227: x[j][i].v = 0.0;
228: x[j][i].omega = 0.0;
229: x[j][i].temp = (grashof>0)*i*dx;
230: }
231: }
233: /*
234: Restore vector
235: */
236: DMDAVecRestoreArray(da,X,&x);
237: return 0;
238: }
240: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,Field **x,Field **xdot,Field **f,void *ptr)
241: {
242: AppCtx *user = (AppCtx*)ptr;
244: PetscInt xints,xinte,yints,yinte,i,j;
245: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
246: PetscReal grashof,prandtl,lid;
247: PetscScalar u,udot,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
250: grashof = user->grashof;
251: prandtl = user->prandtl;
252: lid = user->lidvelocity;
254: /*
255: Define mesh intervals ratios for uniform grid.
257: Note: FD formulae below are normalized by multiplying through by
258: local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.
261: */
262: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
263: hx = 1.0/dhx; hy = 1.0/dhy;
264: hxdhy = hx*dhy; hydhx = hy*dhx;
266: xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
268: /* Test whether we are on the bottom edge of the global array */
269: if (yints == 0) {
270: j = 0;
271: yints = yints + 1;
272: /* bottom edge */
273: for (i=info->xs; i<info->xs+info->xm; i++) {
274: f[j][i].u = x[j][i].u;
275: f[j][i].v = x[j][i].v;
276: f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
277: f[j][i].temp = x[j][i].temp-x[j+1][i].temp;
278: }
279: }
281: /* Test whether we are on the top edge of the global array */
282: if (yinte == info->my) {
283: j = info->my - 1;
284: yinte = yinte - 1;
285: /* top edge */
286: for (i=info->xs; i<info->xs+info->xm; i++) {
287: f[j][i].u = x[j][i].u - lid;
288: f[j][i].v = x[j][i].v;
289: f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
290: f[j][i].temp = x[j][i].temp-x[j-1][i].temp;
291: }
292: }
294: /* Test whether we are on the left edge of the global array */
295: if (xints == 0) {
296: i = 0;
297: xints = xints + 1;
298: /* left edge */
299: for (j=info->ys; j<info->ys+info->ym; j++) {
300: f[j][i].u = x[j][i].u;
301: f[j][i].v = x[j][i].v;
302: f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
303: f[j][i].temp = x[j][i].temp;
304: }
305: }
307: /* Test whether we are on the right edge of the global array */
308: if (xinte == info->mx) {
309: i = info->mx - 1;
310: xinte = xinte - 1;
311: /* right edge */
312: for (j=info->ys; j<info->ys+info->ym; j++) {
313: f[j][i].u = x[j][i].u;
314: f[j][i].v = x[j][i].v;
315: f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
316: f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0);
317: }
318: }
320: /* Compute over the interior points */
321: for (j=yints; j<yinte; j++) {
322: for (i=xints; i<xinte; i++) {
324: /*
325: convective coefficients for upwinding
326: */
327: vx = x[j][i].u; avx = PetscAbsScalar(vx);
328: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
329: vy = x[j][i].v; avy = PetscAbsScalar(vy);
330: vyp = .5*(vy+avy); vym = .5*(vy-avy);
332: /* U velocity */
333: u = x[j][i].u;
334: udot = user->parabolic ? xdot[j][i].u : 0.;
335: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
336: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
337: f[j][i].u = udot + uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
339: /* V velocity */
340: u = x[j][i].v;
341: udot = user->parabolic ? xdot[j][i].v : 0.;
342: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
343: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
344: f[j][i].v = udot + uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
346: /* Omega */
347: u = x[j][i].omega;
348: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
349: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
350: f[j][i].omega = (xdot[j][i].omega + uxx + uyy
351: + (vxp*(u - x[j][i-1].omega)
352: + vxm*(x[j][i+1].omega - u)) * hy
353: + (vyp*(u - x[j-1][i].omega)
354: + vym*(x[j+1][i].omega - u)) * hx
355: - .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy);
357: /* Temperature */
358: u = x[j][i].temp;
359: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
360: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
361: f[j][i].temp = (xdot[j][i].temp + uxx + uyy
362: + prandtl * ((vxp*(u - x[j][i-1].temp)
363: + vxm*(x[j][i+1].temp - u)) * hy
364: + (vyp*(u - x[j-1][i].temp)
365: + vym*(x[j+1][i].temp - u)) * hx));
366: }
367: }
369: /*
370: Flop count (multiply-adds are counted as 2 operations)
371: */
372: PetscLogFlops(84.0*info->ym*info->xm);
373: return(0);
374: }