Actual source code: ex19.c
1: /*$Id: ex19.c,v 1.30 2001/08/07 21:31:17 bsmith Exp $*/
3: static char help[] = "Nonlinear driven cavity with multigrid in 2d.n
4: n
5: The 2D driven cavity problem is solved in a velocity-vorticity formulation.n
6: The flow can be driven with the lid or with bouyancy or both:n
7: -lidvelocity <lid>, where <lid> = dimensionless velocity of lidn
8: -grashof <gr>, where <gr> = dimensionless temperature gradentn
9: -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ration
10: -contours : draw contour plots of solutionnn";
12: /*T
13: Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
14: Concepts: DA^using distributed arrays;
15: Concepts: multicomponent
16: Processors: n
17: T*/
19: /* ------------------------------------------------------------------------
21: We thank David E. Keyes for contributing the driven cavity discretization
22: within this example code.
24: This problem is modeled by the partial differential equation system
25:
26: - Lap(U) - Grad_y(Omega) = 0
27: - Lap(V) + Grad_x(Omega) = 0
28: - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
29: - Lap(T) + PR*Div([U*T,V*T]) = 0
31: in the unit square, which is uniformly discretized in each of x and
32: y in this simple encoding.
34: No-slip, rigid-wall Dirichlet conditions are used for [U,V].
35: Dirichlet conditions are used for Omega, based on the definition of
36: vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
37: constant coordinate boundary, the tangential derivative is zero.
38: Dirichlet conditions are used for T on the left and right walls,
39: and insulation homogeneous Neumann conditions are used for T on
40: the top and bottom walls.
42: A finite difference approximation with the usual 5-point stencil
43: is used to discretize the boundary value problem to obtain a
44: nonlinear system of equations. Upwinding is used for the divergence
45: (convective) terms and central for the gradient (source) terms.
46:
47: The Jacobian can be either
48: * formed via finite differencing using coloring (the default), or
49: * applied matrix-free via the option -snes_mf
50: (for larger grid problems this variant may not converge
51: without a preconditioner due to ill-conditioning).
53: ------------------------------------------------------------------------- */
55: /*
56: Include "petscda.h" so that we can use distributed arrays (DAs).
57: Include "petscsnes.h" so that we can use SNES solvers. Note that this
58: file automatically includes:
59: petsc.h - base PETSc routines petscvec.h - vectors
60: petscsys.h - system routines petscmat.h - matrices
61: petscis.h - index sets petscksp.h - Krylov subspace methods
62: petscviewer.h - viewers petscpc.h - preconditioners
63: petscsles.h - linear solvers
64: */
65: #include petscsnes.h
66: #include petscda.h
68: /*
69: User-defined routines and data structures
70: */
71: typedef struct {
72: PetscScalar u,v,omega,temp;
73: } Field;
75: extern int FormInitialGuess(SNES,Vec,void*);
76: extern int FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
77: extern int FormFunctionLocali(DALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);
79: typedef struct {
80: PassiveReal lidvelocity,prandtl,grashof; /* physical parameters */
81: PetscTruth draw_contours; /* flag - 1 indicates drawing contours */
82: } AppCtx;
84: #undef __FUNCT__
86: int main(int argc,char **argv)
87: {
88: DMMG *dmmg; /* multilevel grid structure */
89: AppCtx user; /* user-defined work context */
90: int mx,my,its;
91: int ierr;
92: MPI_Comm comm;
93: SNES snes;
94: DA da;
96: PetscInitialize(&argc,&argv,(char *)0,help);
97: comm = PETSC_COMM_WORLD;
100: PreLoadBegin(PETSC_TRUE,"SetUp");
101: DMMGCreate(comm,2,&user,&dmmg);
104: /*
105: Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
106: for principal unknowns (x) and governing residuals (f)
107: */
108: DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
109: DMMGSetDM(dmmg,(DM)da);
110: DADestroy(da);
112: DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
113: PETSC_IGNORE,PETSC_IGNORE);
114: /*
115: Problem parameters (velocity of lid, prandtl, and grashof numbers)
116: */
117: user.lidvelocity = 1.0/(mx*my);
118: user.prandtl = 1.0;
119: user.grashof = 1.0;
120: PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
121: PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
122: PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);
123: PetscOptionsHasName(PETSC_NULL,"-contours",&user.draw_contours);
125: DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
126: DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
127: DASetFieldName(DMMGGetDA(dmmg),2,"Omega");
128: DASetFieldName(DMMGGetDA(dmmg),3,"temperature");
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Create user context, set problem data, create vector data structures.
132: Also, compute the initial guess.
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Create nonlinear solver context
138: Process adiC: FormFunctionLocal FormFunctionLocali
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
141: DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);
143: PetscPrintf(comm,"lid velocity = %g, prandtl # = %g, grashof # = %gn",
144: user.lidvelocity,user.prandtl,user.grashof);
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Solve the nonlinear system
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: DMMGSetInitialGuess(dmmg,FormInitialGuess);
152: PreLoadStage("Solve");
153: DMMGSolve(dmmg);
155: snes = DMMGGetSNES(dmmg);
156: SNESGetIterationNumber(snes,&its);
157: PetscPrintf(comm,"Number of Newton iterations = %dn", its);
159: /*
160: Visualize solution
161: */
163: if (user.draw_contours) {
164: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
165: }
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Free work space. All PETSc objects should be destroyed when they
169: are no longer needed.
170: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172: DMMGDestroy(dmmg);
173: PreLoadEnd();
175: PetscFinalize();
176: return 0;
177: }
179: /* ------------------------------------------------------------------- */
182: #undef __FUNCT__
184: /*
185: FormInitialGuess - Forms initial approximation.
187: Input Parameters:
188: user - user-defined application context
189: X - vector
191: Output Parameter:
192: X - vector
193: */
194: int FormInitialGuess(SNES snes,Vec X,void *ptr)
195: {
196: DMMG dmmg = (DMMG)ptr;
197: AppCtx *user = (AppCtx*)dmmg->user;
198: DA da = (DA)dmmg->dm;
199: int i,j,mx,ierr,xs,ys,xm,ym;
200: PetscReal grashof,dx;
201: Field **x;
203: grashof = user->grashof;
205: DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
206: dx = 1.0/(mx-1);
208: /*
209: Get local grid boundaries (for 2-dimensional DA):
210: xs, ys - starting grid indices (no ghost points)
211: xm, ym - widths of local grid (no ghost points)
212: */
213: DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
215: /*
216: Get a pointer to vector data.
217: - For default PETSc vectors, VecGetArray() returns a pointer to
218: the data array. Otherwise, the routine is implementation dependent.
219: - You MUST call VecRestoreArray() when you no longer need access to
220: the array.
221: */
222: DAVecGetArray(da,X,(void**)&x);
224: /*
225: Compute initial guess over the locally owned part of the grid
226: Initial condition is motionless fluid and equilibrium temperature
227: */
228: for (j=ys; j<ys+ym; j++) {
229: for (i=xs; i<xs+xm; i++) {
230: x[j][i].u = 0.0;
231: x[j][i].v = 0.0;
232: x[j][i].omega = 0.0;
233: x[j][i].temp = (grashof>0)*i*dx;
234: }
235: }
237: /*
238: Restore vector
239: */
240: DAVecRestoreArray(da,X,(void**)&x);
241: return 0;
242: }
243: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
244: {
245: AppCtx *user = (AppCtx*)ptr;
246: int ierr,i,j;
247: int xints,xinte,yints,yinte;
248: PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
249: PetscReal grashof,prandtl,lid;
250: PetscScalar u,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.
259: [Note: FD formulae below are normalized by multiplying through by
260: local volume element 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: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
335: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
336: f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
338: /* V velocity */
339: u = x[j][i].v;
340: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
341: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
342: f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
344: /* Omega */
345: u = x[j][i].omega;
346: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
347: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
348: f[j][i].omega = uxx + uyy +
349: (vxp*(u - x[j][i-1].omega) +
350: vxm*(x[j][i+1].omega - u)) * hy +
351: (vyp*(u - x[j-1][i].omega) +
352: vym*(x[j+1][i].omega - u)) * hx -
353: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
355: /* Temperature */
356: u = x[j][i].temp;
357: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
358: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
359: f[j][i].temp = uxx + uyy + prandtl * (
360: (vxp*(u - x[j][i-1].temp) +
361: vxm*(x[j][i+1].temp - u)) * hy +
362: (vyp*(u - x[j-1][i].temp) +
363: vym*(x[j+1][i].temp - u)) * hx);
364: }
365: }
367: /*
368: Flop count (multiply-adds are counted as 2 operations)
369: */
370: PetscLogFlops(84*info->ym*info->xm);
371: return(0);
372: }
374: /*
375: This is an experimental function and can be safely ignored.
376: */
377: int FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
378: {
379: AppCtx *user = (AppCtx*)ptr;
380: int i,j,c;
381: PassiveReal hx,hy,dhx,dhy,hxdhy,hydhx;
382: PassiveReal grashof,prandtl,lid;
383: PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
386: grashof = user->grashof;
387: prandtl = user->prandtl;
388: lid = user->lidvelocity;
390: /*
391: Define mesh intervals ratios for uniform grid.
392: [Note: FD formulae below are normalized by multiplying through by
393: local volume element to obtain coefficients O(1) in two dimensions.]
394: */
395: dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
396: hx = 1.0/dhx; hy = 1.0/dhy;
397: hxdhy = hx*dhy; hydhx = hy*dhx;
399: i = st->i; j = st->j; c = st->c;
401: /* Test whether we are on the right edge of the global array */
402: if (i == info->mx-1) {
403: if (c == 0) *f = x[j][i].u;
404: else if (c == 1) *f = x[j][i].v;
405: else if (c == 2) *f = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
406: else *f = x[j][i].temp - (PetscReal)(grashof>0);
408: /* Test whether we are on the left edge of the global array */
409: } else if (i == 0) {
410: if (c == 0) *f = x[j][i].u;
411: else if (c == 1) *f = x[j][i].v;
412: else if (c == 2) *f = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
413: else *f = x[j][i].temp;
415: /* Test whether we are on the top edge of the global array */
416: } else if (j == info->my-1) {
417: if (c == 0) *f = x[j][i].u - lid;
418: else if (c == 1) *f = x[j][i].v;
419: else if (c == 2) *f = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
420: else *f = x[j][i].temp-x[j-1][i].temp;
422: /* Test whether we are on the bottom edge of the global array */
423: } else if (j == 0) {
424: if (c == 0) *f = x[j][i].u;
425: else if (c == 1) *f = x[j][i].v;
426: else if (c == 2) *f = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
427: else *f = x[j][i].temp-x[j+1][i].temp;
429: /* Compute over the interior points */
430: } else {
431: /*
432: convective coefficients for upwinding
433: */
434: vx = x[j][i].u; avx = PetscAbsScalar(vx);
435: vxp = .5*(vx+avx); vxm = .5*(vx-avx);
436: vy = x[j][i].v; avy = PetscAbsScalar(vy);
437: vyp = .5*(vy+avy); vym = .5*(vy-avy);
439: /* U velocity */
440: if (c == 0) {
441: u = x[j][i].u;
442: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
443: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
444: *f = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
446: /* V velocity */
447: } else if (c == 1) {
448: u = x[j][i].v;
449: uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
450: uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
451: *f = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
452:
453: /* Omega */
454: } else if (c == 2) {
455: u = x[j][i].omega;
456: uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
457: uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
458: *f = uxx + uyy +
459: (vxp*(u - x[j][i-1].omega) +
460: vxm*(x[j][i+1].omega - u)) * hy +
461: (vyp*(u - x[j-1][i].omega) +
462: vym*(x[j+1][i].omega - u)) * hx -
463: .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;
464:
465: /* Temperature */
466: } else {
467: u = x[j][i].temp;
468: uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
469: uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
470: *f = uxx + uyy + prandtl * (
471: (vxp*(u - x[j][i-1].temp) +
472: vxm*(x[j][i+1].temp - u)) * hy +
473: (vyp*(u - x[j-1][i].temp) +
474: vym*(x[j+1][i].temp - u)) * hx);
475: }
476: }
478: return(0);
479: }