Actual source code: ex50.c

  2: static char help[] = "Nonlinear driven cavity with multigrid 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:       The same as ex19.c except it does not use DMMG, it uses its replacement.
 12:       See src/ksp/ksp/examples/tutorials/ex45.c
 13: */

 15: /*T
 16:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
 17:    Concepts: DMDA^using distributed arrays;
 18:    Concepts: multicomponent
 19:    Processors: n
 20: T*/
 21: /* ------------------------------------------------------------------------

 23:     We thank David E. Keyes for contributing the driven cavity discretization
 24:     within this example code.

 26:     This problem is modeled by the partial differential equation system
 27:   
 28:         - Lap(U) - Grad_y(Omega) = 0
 29:         - Lap(V) + Grad_x(Omega) = 0
 30:         - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
 31:         - Lap(T) + PR*Div([U*T,V*T]) = 0

 33:     in the unit square, which is uniformly discretized in each of x and
 34:     y in this simple encoding.

 36:     No-slip, rigid-wall Dirichlet conditions are used for [U,V].
 37:     Dirichlet conditions are used for Omega, based on the definition of
 38:     vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each
 39:     constant coordinate boundary, the tangential derivative is zero.
 40:     Dirichlet conditions are used for T on the left and right walls,
 41:     and insulation homogeneous Neumann conditions are used for T on
 42:     the top and bottom walls. 

 44:     A finite difference approximation with the usual 5-point stencil 
 45:     is used to discretize the boundary value problem to obtain a 
 46:     nonlinear system of equations.  Upwinding is used for the divergence
 47:     (convective) terms and central for the gradient (source) terms.
 48:     
 49:     The Jacobian can be either
 50:       * formed via finite differencing using coloring (the default), or
 51:       * applied matrix-free via the option -snes_mf 
 52:         (for larger grid problems this variant may not converge 
 53:         without a preconditioner due to ill-conditioning).

 55:   ------------------------------------------------------------------------- */

 57: /* 
 58:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 59:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 60:    file automatically includes:
 61:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 62:      petscmat.h - matrices
 63:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 64:      petscviewer.h - viewers               petscpc.h  - preconditioners
 65:      petscksp.h   - linear solvers 
 66: */
 67: #include <petscsnes.h>
 68: #include <petscdmda.h>
 69: #include <petscdmmg.h>

 71: /* 
 72:    User-defined routines and data structures
 73: */
 74: typedef struct {
 75:   PetscScalar u,v,omega,temp;
 76: } Field;

 78: PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);
 79: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);

 81: typedef struct {
 82:    PassiveReal  lidvelocity,prandtl,grashof;  /* physical parameters */
 83:    PetscBool    draw_contours;                /* flag - 1 indicates drawing contours */
 84: } AppCtx;

 86: PetscErrorCode FormInitialGuess(AppCtx*,DM,Vec);

 90: int main(int argc,char **argv)
 91: {
 92:   AppCtx         user;                /* user-defined work context */
 93:   PetscInt       mx,my,its;
 95:   MPI_Comm       comm;
 96:   SNES           snes;
 97:   DM             da;
 98:   Vec            x;

100:   PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return(1);
101:   comm = PETSC_COMM_WORLD;

103:   SNESCreate(comm,&snes);
104: 
105:   /*
106:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
107:       for principal unknowns (x) and governing residuals (f)
108:   */
109:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);
110:   SNESSetDM(snes,(DM)da);
111: 
112:   DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
113:                    PETSC_IGNORE,PETSC_IGNORE,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:   DMDASetFieldName(da,0,"x-velocity");
126:   DMDASetFieldName(da,1,"y-velocity");
127:   DMDASetFieldName(da,2,"Omega");
128:   DMDASetFieldName(da,3,"temperature");
129: 
130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131:      Create user context, set problem data, create vector data structures.
132:      Also, compute the initial guess.
133:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: 
135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Create nonlinear solver context
137:      
138:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139:   DMSetApplicationContext(da,&user);
140:   DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
141:   SNESSetFromOptions(snes);
142: 
143:   PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, grashof # = %G\n",user.lidvelocity,user.prandtl,user.grashof);


146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147:      Solve the nonlinear system
148:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149:   DMCreateGlobalVector(da,&x);
150:   FormInitialGuess(&user,da,x);
151: 
152:   SNESSolve(snes,PETSC_NULL,x);
153: 
154:   SNESGetIterationNumber(snes,&its);
155:   PetscPrintf(comm,"Number of Newton iterations = %D\n", its);
156: 
157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158:      Free work space.  All PETSc objects should be destroyed when they
159:      are no longer needed.
160:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   VecDestroy(&x);
162:   DMDestroy(&da);
163:   SNESDestroy(&snes);
164: 
165:   PetscFinalize();
166:   return 0;
167: }

169: /* ------------------------------------------------------------------- */


174: /* 
175:    FormInitialGuess - Forms initial approximation.

177:    Input Parameters:
178:    user - user-defined application context
179:    X - vector

181:    Output Parameter:
182:    X - vector
183:  */
184: PetscErrorCode FormInitialGuess(AppCtx *user,DM da,Vec X)
185: {
186:   PetscInt       i,j,mx,xs,ys,xm,ym;
188:   PetscReal      grashof,dx;
189:   Field          **x;

191:   grashof = user->grashof;

193:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
194:   dx  = 1.0/(mx-1);

196:   /*
197:      Get local grid boundaries (for 2-dimensional DMDA):
198:        xs, ys   - starting grid indices (no ghost points)
199:        xm, ym   - widths of local grid (no ghost points)
200:   */
201:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

203:   /*
204:      Get a pointer to vector data.
205:        - For default PETSc vectors, VecGetArray() returns a pointer to
206:          the data array.  Otherwise, the routine is implementation dependent.
207:        - You MUST call VecRestoreArray() when you no longer need access to
208:          the array.
209:   */
210:   DMDAVecGetArray(da,X,&x);

212:   /*
213:      Compute initial guess over the locally owned part of the grid
214:      Initial condition is motionless fluid and equilibrium temperature
215:   */
216:   for (j=ys; j<ys+ym; j++) {
217:     for (i=xs; i<xs+xm; i++) {
218:       x[j][i].u     = 0.0;
219:       x[j][i].v     = 0.0;
220:       x[j][i].omega = 0.0;
221:       x[j][i].temp  = (grashof>0)*i*dx;
222:     }
223:   }

225:   /*
226:      Restore vector
227:   */
228:   DMDAVecRestoreArray(da,X,&x);
229:   return 0;
230: }
231: 
234: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
235:  {
236:   AppCtx         *user = (AppCtx*)ptr;
238:   PetscInt       xints,xinte,yints,yinte,i,j;
239:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
240:   PetscReal      grashof,prandtl,lid;
241:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

244:   grashof = user->grashof;
245:   prandtl = user->prandtl;
246:   lid     = user->lidvelocity;

248:   /* 
249:      Define mesh intervals ratios for uniform grid.

251:      Note: FD formulae below are normalized by multiplying through by
252:      local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions.

254:      
255:   */
256:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
257:   hx = 1.0/dhx;                   hy = 1.0/dhy;
258:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

260:   xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;

262:   /* Test whether we are on the bottom edge of the global array */
263:   if (yints == 0) {
264:     j = 0;
265:     yints = yints + 1;
266:     /* bottom edge */
267:     for (i=info->xs; i<info->xs+info->xm; i++) {
268:       f[j][i].u     = x[j][i].u;
269:       f[j][i].v     = x[j][i].v;
270:       f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
271:       f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
272:     }
273:   }

275:   /* Test whether we are on the top edge of the global array */
276:   if (yinte == info->my) {
277:     j = info->my - 1;
278:     yinte = yinte - 1;
279:     /* top edge */
280:     for (i=info->xs; i<info->xs+info->xm; i++) {
281:         f[j][i].u     = x[j][i].u - lid;
282:         f[j][i].v     = x[j][i].v;
283:         f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
284:         f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
285:     }
286:   }

288:   /* Test whether we are on the left edge of the global array */
289:   if (xints == 0) {
290:     i = 0;
291:     xints = xints + 1;
292:     /* left edge */
293:     for (j=info->ys; j<info->ys+info->ym; j++) {
294:       f[j][i].u     = x[j][i].u;
295:       f[j][i].v     = x[j][i].v;
296:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
297:       f[j][i].temp  = x[j][i].temp;
298:     }
299:   }

301:   /* Test whether we are on the right edge of the global array */
302:   if (xinte == info->mx) {
303:     i = info->mx - 1;
304:     xinte = xinte - 1;
305:     /* right edge */
306:     for (j=info->ys; j<info->ys+info->ym; j++) {
307:       f[j][i].u     = x[j][i].u;
308:       f[j][i].v     = x[j][i].v;
309:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
310:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
311:     }
312:   }

314:   /* Compute over the interior points */
315:   for (j=yints; j<yinte; j++) {
316:     for (i=xints; i<xinte; i++) {

318:         /*
319:           convective coefficients for upwinding
320:         */
321:         vx = x[j][i].u; avx = PetscAbsScalar(vx);
322:         vxp = .5*(vx+avx); vxm = .5*(vx-avx);
323:         vy = x[j][i].v; avy = PetscAbsScalar(vy);
324:         vyp = .5*(vy+avy); vym = .5*(vy-avy);

326:         /* U velocity */
327:         u          = x[j][i].u;
328:         uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
329:         uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
330:         f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

332:         /* V velocity */
333:         u          = x[j][i].v;
334:         uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
335:         uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
336:         f[j][i].v  = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

338:         /* Omega */
339:         u          = x[j][i].omega;
340:         uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
341:         uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
342:         f[j][i].omega = uxx + uyy +
343:                         (vxp*(u - x[j][i-1].omega) +
344:                           vxm*(x[j][i+1].omega - u)) * hy +
345:                         (vyp*(u - x[j-1][i].omega) +
346:                           vym*(x[j+1][i].omega - u)) * hx -
347:                         .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

349:         /* Temperature */
350:         u             = x[j][i].temp;
351:         uxx           = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
352:         uyy           = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
353:         f[j][i].temp =  uxx + uyy  + prandtl * (
354:                         (vxp*(u - x[j][i-1].temp) +
355:                           vxm*(x[j][i+1].temp - u)) * hy +
356:                         (vyp*(u - x[j-1][i].temp) +
357:                                  vym*(x[j+1][i].temp - u)) * hx);
358:     }
359:   }

361:   /*
362:      Flop count (multiply-adds are counted as 2 operations)
363:   */
364:   PetscLogFlops(84.0*info->ym*info->xm);
365:   return(0);
366: }

370: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *user)
371: {
372:   DMDALocalInfo  info;
373:   Field          **u,**fu;
375:   Vec            localX;
376:   DM             da;

379:   SNESGetDM(snes,(DM*)&da);
380:   DMGetLocalVector(da,&localX);
381:   /*
382:      Scatter ghost points to local vector, using the 2-step process
383:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
384:   */
385:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
386:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
387:   DMDAGetLocalInfo(da,&info);
388:   DMDAVecGetArray(da,localX,&u);
389:   DMDAVecGetArray(da,F,&fu);
390:   FormFunctionLocal(&info,u,fu,user);
391:   DMDAVecRestoreArray(da,localX,&u);
392:   DMDAVecRestoreArray(da,F,&fu);
393:   DMRestoreLocalVector(da,&localX);
394:   return(0);
395: }