Actual source code: ex32.c

  2: static char help[] = "Model multi-physics solver. Modified from ex19.c \n\\n\
  3: The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\
  4: The flow can be driven with the lid or with bouyancy or both:\n\
  5:   -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\
  6:   -grashof <gr>, where <gr> = dimensionless temperature gradent\n\
  7:   -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\
  8:   -contours : draw contour plots of solution\n\n";

 10: /*T
 11:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
 12:    Concepts: DMDA^using distributed arrays;
 13:    Concepts: multicomponent
 14:    Processors: n
 15: T*/

 17: /* ------------------------------------------------------------------------

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

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

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

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

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

 51:   ------------------------------------------------------------------------- */
 52: #include <petscsnes.h>
 53: #include <petscdmda.h>
 54: #include <petscdmcomposite.h>
 55: #include <petscdmmg.h>

 57: /* User-defined routines and data structure */
 58: typedef struct {
 59:   PetscScalar u,v,omega,temp;
 60: } Field;

 62: typedef struct {
 63:   PetscScalar u,v,omega;
 64: } Field1;

 66: typedef struct {
 67:   PetscScalar temp;
 68: } Field2;

 70: typedef struct {
 71:   PassiveReal  lidvelocity,prandtl,grashof;  /* physical parameters */
 72:   PetscBool    draw_contours;                /* flag - 1 indicates drawing contours */
 73:   DMMG         *dmmg,*dmmg_comp;             /* used by MySolutionView() */
 74:   DMMG         *dmmg1,*dmmg2; /* passing objects of sub-physics into the composite physics - used by FormInitialGuessComp() */
 75:   DM           pack;
 76:   PetscBool    COMPOSITE_MODEL;
 77:   Field        **x;   /* passing DMMGGetx(dmmg) - final solution of original coupled physics */
 78:   Field1       **x1;  /* passing local ghosted vector array of Physics 1 */
 79:   Field2       **x2;  /* passing local ghosted vector array of Physics 2 */
 80: } AppCtx;




 96: int main(int argc,char **argv)
 97: {
 98:   DMMG           *dmmg,*dmmg1,*dmmg2,*dmmg_comp; /* multilevel grid structure */
 99:   AppCtx         user;                /* user-defined work context */
100:   PetscInt       mx,my,its,dof,nlevels=1;
102:   MPI_Comm       comm;
103:   SNES           snes;
104:   DM             da;
105:   PetscBool      ViewSolu=PETSC_FALSE,CompSolu=PETSC_TRUE,DoSubPhysics=PETSC_FALSE;
106:   Vec            solu_true,solu_local;

108:   PetscInitialize(&argc,&argv,(char *)0,help);
109:   comm = PETSC_COMM_WORLD;

111:   PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,PETSC_NULL);
112:   DMMGCreate(comm,nlevels,&user,&dmmg);

114:   /*
115:     Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
116:     for principal unknowns (x) and governing residuals (f)
117:   */
118:   dof  = 4;
119:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);
120:   DMMGSetDM(dmmg,(DM)da);
121:   DMDestroy(&da);

123:   DMDAGetInfo(DMMGGetDM(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
124:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

126:   /* Problem parameters (velocity of lid, prandtl, and grashof numbers) */
127:   user.lidvelocity = 1.0/(mx*my);
128:   user.prandtl     = 1.0;
129:   user.grashof     = 1.0;
130:   PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
131:   PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
132:   PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);
133:   PetscOptionsHasName(PETSC_NULL,"-contours",&user.draw_contours);
134:   PetscOptionsHasName(PETSC_NULL,"-view_solu",&ViewSolu);
135:   PetscOptionsHasName(PETSC_NULL,"-do_subphysics",&DoSubPhysics);

137:   DMDASetFieldName(DMMGGetDM(dmmg),0,"x-velocity");
138:   DMDASetFieldName(DMMGGetDM(dmmg),1,"y-velocity");
139:   DMDASetFieldName(DMMGGetDM(dmmg),2,"Omega");
140:   DMDASetFieldName(DMMGGetDM(dmmg),3,"temperature");

142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:      Create user context, set problem data, create vector data structures.
144:      Also, compute the initial guess.
145:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

147:   /* Create nonlinear solver context */
148:   DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
149:   DMMGSetFromOptions(dmmg);
150:   PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, grashof # = %G\n",
151:                        user.lidvelocity,user.prandtl,user.grashof);
152:   DMMGSetInitialGuess(dmmg,FormInitialGuessLocal);

154:   /* Solve the nonlinear system */
155:   DMMGSolve(dmmg);

157:   snes = DMMGGetSNES(dmmg);
158:   SNESGetIterationNumber(snes,&its);
159:   PetscPrintf(comm,"Origianl Physics: Number of Newton iterations = %D\n\n", its);
160: 
161:   /* Save the ghosted local solu_true to be used by Physics 1 and Physics 2 */
162:   da        = DMMGGetDM(dmmg);
163:   solu_true = DMMGGetx(dmmg);
164:   DMCreateLocalVector(da,&solu_local);
165:   DMGlobalToLocalBegin(da,solu_true,INSERT_VALUES,solu_local);
166:   DMGlobalToLocalEnd(da,solu_true,INSERT_VALUES,solu_local);
167:   DMDAVecGetArray(da,solu_local,(Field **)&user.x);

169:   /* Visualize solution */
170:   if (user.draw_contours) {
171:     VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
172:   }
173:   if (ViewSolu){
174:     user.dmmg = dmmg;
175:     MySolutionView(comm,0,&user);
176:   }

178:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179:      Setup Physics 1: 
180:         - Lap(U) - Grad_y(Omega) = 0
181:         - Lap(V) + Grad_x(Omega) = 0
182:         - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
183:         where T is given by the computed x.temp
184:         - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185:   DMMGCreate(comm,nlevels,&user,&dmmg1);
186:   dof  = 3;
187:   user.COMPOSITE_MODEL = PETSC_FALSE;
188:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);
189:   DMMGSetDM(dmmg1,(DM)da);
190:   DMDASetFieldName(da,0,"x-velocity");
191:   DMDASetFieldName(da,1,"y-velocity");
192:   DMDASetFieldName(da,2,"Omega");
193:   DMDestroy(&da);

195:   if (DoSubPhysics){
196:     DMMGSetSNESLocal(dmmg1,FormFunctionLocal1,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
197:     DMMGSetFromOptions(dmmg1);
198:     DMMGSetInitialGuess(dmmg1,FormInitialGuessLocal1);

200:     DMMGSolve(dmmg1);
201:     snes = DMMGGetSNES(dmmg1);
202:     SNESGetIterationNumber(snes,&its);
203:     PetscPrintf(comm,"SubPhysics 1, Number of Newton iterations = %D\n\n", its);
204: 
205:     if (ViewSolu){ /* View individial componets of the solution */
206:       user.dmmg1 = dmmg1;
207:       MySolutionView(comm,1,&user);
208:     }
209:   }

211:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212:      Setup Physics 2: 
213:         - Lap(T) + PR*Div([U*T,V*T]) = 0        
214:         where U and V are given by the computed x.u and x.v
215:         - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216:   DMMGCreate(comm,nlevels,&user,&dmmg2);
217:   dof  = 1;
218:   user.COMPOSITE_MODEL = PETSC_FALSE;
219:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da);
220:   DMMGSetDM(dmmg2,(DM)da);
221:   DMDASetFieldName(da,0,"temperature");
222:   DMDestroy(&da);

224:   if (DoSubPhysics){
225:     DMMGSetSNESLocal(dmmg2,FormFunctionLocal2,0,ad_FormFunctionLocal,admf_FormFunctionLocal);
226:     DMMGSetFromOptions(dmmg2);
227:     DMMGSetInitialGuess(dmmg2,FormInitialGuessLocal2);

229:     DMMGSolve(dmmg2);
230:     snes = DMMGGetSNES(dmmg2);
231:     SNESGetIterationNumber(snes,&its);
232:     PetscPrintf(comm,"SubPhysics 2, Number of Newton iterations = %D\n\n", its);

234:     if (ViewSolu){ /* View individial componets of the solution */
235:       user.dmmg2 = dmmg2;
236:       MySolutionView(comm,2,&user);
237:     }
238:   }
239:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240:     Create the DMComposite object to manage the two grids/physics. 
241:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
242:   DMCompositeCreate(comm,&user.pack);
243:   DMCompositeAddDM(user.pack,(DM)DMMGGetDM(dmmg1));
244:   DMCompositeAddDM(user.pack,(DM)DMMGGetDM(dmmg2));

246:   /* Create the solver object and attach the grid/physics info */
247:   DMMGCreate(comm,nlevels,&user,&dmmg_comp);
248:   DMMGSetDM(dmmg_comp,(DM)user.pack);
249:   CHKMEMQ;
250:   DMMGSetISColoringType(dmmg_comp,IS_COLORING_GLOBAL);

252:   user.dmmg1 = dmmg1;
253:   user.dmmg2 = dmmg2;
254:   user.COMPOSITE_MODEL = PETSC_TRUE;
255:   DMMGSetInitialGuess(dmmg_comp,FormInitialGuessComp);
256:   DMMGSetSNES(dmmg_comp,FormFunctionComp,0);
257:   DMMGSetFromOptions(dmmg_comp);

259:   /* Solve the nonlinear system */
260:   DMMGSolve(dmmg_comp);
261:   snes = DMMGGetSNES(dmmg_comp);
262:   SNESGetIterationNumber(snes,&its);
263:   PetscPrintf(comm,"Composite Physics: Number of Newton iterations = %D\n\n", its);

265:   /* Compare the solutions obtained from dmmg and dmmg_comp */
266:   if (ViewSolu || CompSolu){
267:     /* View the solution of composite physics (phy_num = 3)
268:        or Compare solutions (phy_num > 3) */
269:     user.dmmg      = dmmg;
270:     user.dmmg_comp = dmmg_comp;
271:     MySolutionView(comm,4,&user);
272:   }

274:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275:      Free spaces 
276:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
277:   DMDAVecRestoreArray(DMMGGetDM(dmmg),solu_local,(Field **)&user.x);
278:   VecDestroy(&solu_local);
279:   DMDestroy(&user.pack);
280:   DMMGDestroy(dmmg);
281:   DMMGDestroy(dmmg1);
282:   DMMGDestroy(dmmg2);
283:   DMMGDestroy(dmmg_comp);
284:   PetscFinalize();
285:   return 0;
286: }

288: /* ------------------------------------------------------------------- */

292: /* 
293:    FormInitialGuessLocal - Forms initial approximation for this process

295:    Input Parameters:
296:      user - user-defined application context
297:      X    - vector (DMDA local vector)

299:    Output Parameter:
300:      X - vector with the local values set
301:  */
302: PetscErrorCode FormInitialGuessLocal(DMMG dmmg,Vec X)
303: {
304:   AppCtx         *user = (AppCtx*)dmmg->user;
305:   DM             da = dmmg->dm;
306:   PetscInt       i,j,mx,xs,ys,xm,ym;
308:   PetscReal      grashof,dx;
309:   Field          **x;

311:   grashof = user->grashof;
312:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
313:   dx  = 1.0/(mx-1);

315:   /*
316:      Get local grid boundaries (for 2-dimensional DMDA):
317:        xs, ys   - starting grid indices (no ghost points)
318:        xm, ym   - widths of local grid (no ghost points)
319:   */
320:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

322:   /*
323:      Get a pointer to vector data.
324:        - For default PETSc vectors, VecGetArray() returns a pointer to
325:          the data array.  Otherwise, the routine is implementation dependent.
326:        - You MUST call VecResetArraystoreArray() when you no longer need access to
327:          the array.
328:   */
329:   DMDAVecGetArray(da,X,&x);

331:   /*
332:      Compute initial guess over the locally owned part of the grid
333:      Initial condition is motionless fluid and equilibrium temperature
334:      U = V = Omega = 0.0; Temp[x_i, y_j] = x_i;
335:   */
336:   for (j=ys; j<ys+ym; j++) {
337:     for (i=xs; i<xs+xm; i++) {
338:       x[j][i].u     = 0.0;
339:       x[j][i].v     = 0.0;
340:       x[j][i].omega = 0.0;
341:       x[j][i].temp  = (grashof>0)*i*dx;
342:     }
343:   }

345:   /* Restore vector */
346:   DMDAVecRestoreArray(da,X,&x);
347:   return 0;
348: }

352: /* Form initial guess for Physic 1 */
353: PetscErrorCode FormInitialGuessLocal1(DMMG dmmg,Vec X)
354: {
355:   DM             da = dmmg->dm;
356:   PetscInt       i,j,mx,xs,ys,xm,ym;
358:   Field1         **x;

360:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
361:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

363:   DMDAVecGetArray(da,X,&x);
364:   for (j=ys; j<ys+ym; j++) {
365:     for (i=xs; i<xs+xm; i++) {
366:       x[j][i].u     = 0.0;
367:       x[j][i].v     = 0.0;
368:       x[j][i].omega = 0.0;
369:     }
370:   }
371:   DMDAVecRestoreArray(da,X,&x);
372:   return 0;
373: }

377: /* Form initial guess for Physic 2 */
378: PetscErrorCode FormInitialGuessLocal2(DMMG dmmg,Vec X)
379: {
380:   AppCtx         *user = (AppCtx*)dmmg->user;
381:   DM             da = dmmg->dm;
382:   PetscInt       i,j,mx,xs,ys,xm,ym;
384:   PetscReal      grashof,dx;
385:   Field2         **x;

387:   grashof = user->grashof;
388:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
389:   dx  = 1.0/(mx-1);

391:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

393:   DMDAVecGetArray(da,X,&x);
394:   for (j=ys; j<ys+ym; j++) {
395:     for (i=xs; i<xs+xm; i++) {
396:       x[j][i].temp  = (grashof>0)*i*dx;
397:     }
398:   }
399:   DMDAVecRestoreArray(da,X,&x);
400:   return 0;
401: }

405: /* 
406:    FormInitialGuessComp - 
407:               Forms the initial guess for the composite model
408:               Unwraps the global solution vector and passes its local pieces into the user functions
409:  */
410: PetscErrorCode FormInitialGuessComp(DMMG dmmg,Vec X)
411: {
413:   AppCtx         *user = (AppCtx*)dmmg->user;
414:   DMMG           *dmmg1 = user->dmmg1,*dmmg2=user->dmmg2;
415:   DM             dm = dmmg->dm;
416:   Vec            X1,X2;

419:   /* Access the subvectors in X */
420:   DMCompositeGetAccess(dm,X,&X1,&X2);

422:   /* Evaluate local user provided function */
423:   FormInitialGuessLocal1(*dmmg1,X1);
424:   FormInitialGuessLocal2(*dmmg2,X2);

426:   DMCompositeRestoreAccess(dm,X,&X1,&X2);
427:   return(0);
428: }

432: /* 
433:    FormFunctionLocal - Function evaluation for this process

435:    Input Parameters:
436:      info - DMDALocalInfo context
437:      x    - (DMDA local) vector array including ghost points
438:      f    - (DMDA local) vector array to be evaluated 
439:      ptr  - user-defined application context

441:    Output Parameter:
442:      f - array holds local function values 
443:      
444:      f.u     = - Lap(U) - Grad_y(Omega)
445:      f.v     = - Lap(V) + Grad_x(Omega)
446:      f.omega = - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T)
447:      f.temp  = - Lap(T) + PR*Div([U*T,V*T])
448: */
449: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
450:  {
451:   AppCtx         *user = (AppCtx*)ptr;
453:   PetscInt       xints,xinte,yints,yinte,i,j;
454:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
455:   PetscReal      grashof,prandtl,lid;
456:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

459:   grashof = user->grashof;
460:   prandtl = user->prandtl;
461:   lid     = user->lidvelocity;

463:   /* 
464:      Define mesh intervals ratios for uniform grid.

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

469:      
470:   */
471:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
472:   hx = 1.0/dhx;                   hy = 1.0/dhy;
473:   hxdhy = hx*dhy; /* hx/hy */     hydhx = hy*dhx; /* hy/hx */

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

477:   /* Test whether we are on the bottom edge of the global array */
478:   if (yints == 0) {
479:     j = 0;
480:     yints = yints + 1;
481:     /* bottom edge: 
482:          U = V = 0; Omega = -Grad_y(U)+Grad_x(V); dT/dn = 0; */
483:     for (i=info->xs; i<info->xs+info->xm; i++) {
484:       f[j][i].u     = x[j][i].u;
485:       f[j][i].v     = x[j][i].v;
486:       f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
487:       f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
488:     }
489:   }

491:   /* Test whether we are on the top edge of the global array */
492:   if (yinte == info->my) {
493:     j = info->my - 1;
494:     yinte = yinte - 1;
495:     /* top edge:
496:          U = lid; V = 0; Omega = -Grad_y(U)+Grad_x(V); dT/dn = 0; */
497:     for (i=info->xs; i<info->xs+info->xm; i++) {
498:         f[j][i].u     = x[j][i].u - lid;
499:         f[j][i].v     = x[j][i].v;
500:         f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
501:         f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
502:     }
503:   }

505:   /* Test whether we are on the left edge of the global array */
506:   if (xints == 0) {
507:     i = 0;
508:     xints = xints + 1;
509:     /* left edge:
510:          U = V = 0; Omega = -Grad_y(U)+Grad_x(V); T = 0; */
511:     for (j=info->ys; j<info->ys+info->ym; j++) {
512:       f[j][i].u     = x[j][i].u;
513:       f[j][i].v     = x[j][i].v;
514:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
515:       f[j][i].temp  = x[j][i].temp;
516:     }
517:   }

519:   /* Test whether we are on the right edge of the global array */
520:   if (xinte == info->mx) {
521:     i = info->mx - 1;
522:     xinte = xinte - 1;
523:     /* right edge:
524:          U = V = 0; Omega = -Grad_y(U)+Grad_x(V); T = grashof; */
525:     for (j=info->ys; j<info->ys+info->ym; j++) {
526:       f[j][i].u     = x[j][i].u;
527:       f[j][i].v     = x[j][i].v;
528:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
529:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
530:     }
531:   }

533:   /* Compute over the interior points */
534:   for (j=yints; j<yinte; j++) {
535:     for (i=xints; i<xinte; i++) {
536:       /* convective coefficients for upwinding */
537:       vx = x[j][i].u; avx = PetscAbsScalar(vx);
538:       vxp = .5*(vx+avx); vxm = .5*(vx-avx);
539:       vy = x[j][i].v; avy = PetscAbsScalar(vy);
540:       vyp = .5*(vy+avy); vym = .5*(vy-avy);

542:       /* U velocity */
543:       u          = x[j][i].u;
544:       uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
545:       uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
546:       f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

548:       /* V velocity */
549:       u          = x[j][i].v;
550:       uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
551:       uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
552:       f[j][i].v  = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

554:       /* Omega */
555:       u          = x[j][i].omega;
556:       uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
557:       uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
558:       f[j][i].omega = uxx + uyy
559:                       + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u)) * hy
560:                       + (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u)) * hx
561:                       -        .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

563:       /* Temperature */
564:       u            = x[j][i].temp;
565:       uxx          = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
566:       uyy          = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
567:       f[j][i].temp = uxx + uyy
568:                      + prandtl * ((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u)) * hy
569:                      + (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u)) * hx);
570:     }
571:   }

573:   /* Flop count (multiply-adds are counted as 2 operations) */
574:   PetscLogFlops(84.0*info->ym*info->xm);
575:   return(0);
576: }

580: /* 
581:     Form function for Physics 1: 
582:       same as FormFunctionLocal() except without f.temp and x.temp.
583:       the input x.temp comes from the solu_true 
584: */
585: PetscErrorCode FormFunctionLocal1(DMDALocalInfo *info,Field1 **x,Field1 **f,void *ptr)
586:  {
587:   AppCtx         *user = (AppCtx*)ptr;
588:   PetscInt       xints,xinte,yints,yinte,i,j;
589:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
590:   PetscReal      grashof,lid; /* ,prandtl */
591:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
592:   Field          **solu=user->x;
593:   Field2         **solu2=user->x2;

596:   grashof = user->grashof;
597:   /* prandtl = user->prandtl; */
598:   lid     = user->lidvelocity;

600:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
601:   hx = 1.0/dhx;                   hy = 1.0/dhy;
602:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

606:   /* Test whether we are on the bottom edge of the global array */
607:   if (yints == 0) {
608:     j = 0;
609:     yints = yints + 1;
610:     /* bottom edge */
611:     for (i=info->xs; i<info->xs+info->xm; i++) {
612:       f[j][i].u     = x[j][i].u;
613:       f[j][i].v     = x[j][i].v;
614:       f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
615:     }
616:   }

618:   /* Test whether we are on the top edge of the global array */
619:   if (yinte == info->my) {
620:     j = info->my - 1;
621:     yinte = yinte - 1;
622:     /* top edge */
623:     for (i=info->xs; i<info->xs+info->xm; i++) {
624:         f[j][i].u     = x[j][i].u - lid;
625:         f[j][i].v     = x[j][i].v;
626:         f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
627:     }
628:   }

630:   /* Test whether we are on the left edge of the global array */
631:   if (xints == 0) {
632:     i = 0;
633:     xints = xints + 1;
634:     /* left edge */
635:     for (j=info->ys; j<info->ys+info->ym; j++) {
636:       f[j][i].u     = x[j][i].u;
637:       f[j][i].v     = x[j][i].v;
638:       f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
639:     }
640:   }

642:   /* Test whether we are on the right edge of the global array */
643:   if (xinte == info->mx) {
644:     i = info->mx - 1;
645:     xinte = xinte - 1;
646:     /* right edge */
647:     for (j=info->ys; j<info->ys+info->ym; j++) {
648:       f[j][i].u     = x[j][i].u;
649:       f[j][i].v     = x[j][i].v;
650:       f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
651:     }
652:   }

654:   /* Compute over the interior points */
655:   for (j=yints; j<yinte; j++) {
656:     for (i=xints; i<xinte; i++) {
657:       /* convective coefficients for upwinding */
658:         vx = x[j][i].u; avx = PetscAbsScalar(vx);
659:         vxp = .5*(vx+avx); vxm = .5*(vx-avx);
660:         vy = x[j][i].v; avy = PetscAbsScalar(vy);
661:         vyp = .5*(vy+avy); vym = .5*(vy-avy);

663:         /* U velocity */
664:         u          = x[j][i].u;
665:         uxx        = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
666:         uyy        = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
667:         f[j][i].u  = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

669:         /* V velocity */
670:         u          = x[j][i].v;
671:         uxx        = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
672:         uyy        = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
673:         f[j][i].v  = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

675:         /* Omega */
676:         u          = x[j][i].omega;
677:         uxx        = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
678:         uyy        = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
679:         f[j][i].omega = uxx + uyy
680:                         + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u)) * hy
681:           + (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u)) * hx;
682:         if (user->COMPOSITE_MODEL){
683:           f[j][i].omega += - .5 * grashof * (solu2[j][i+1].temp - solu2[j][i-1].temp) * hy;
684:         } else {
685:           f[j][i].omega += - .5 * grashof * (solu[j][i+1].temp - solu[j][i-1].temp) * hy;
686:         }
687:     }
688:   }
689:   return(0);
690: }

694: /* 
695:     Form function for Physics 2: 
696:       same as FormFunctionLocal() but only has f.temp and x.temp.
697:       the input x.u and x.v come from solu_true 
698: */
699: PetscErrorCode FormFunctionLocal2(DMDALocalInfo *info,Field2 **x,Field2 **f,void *ptr)
700:  {
701:   AppCtx         *user = (AppCtx*)ptr;
702:   PetscInt       xints,xinte,yints,yinte,i,j;
703:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx;
704:   PetscReal      grashof,prandtl; /* ,lid */
705:   PetscScalar    u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
706:   Field          **solu=user->x;
707:   Field1         **solu1=user->x1;

710:   grashof = user->grashof;
711:   prandtl = user->prandtl;
712:   /* lid     = user->lidvelocity; */

714:   dhx = (PetscReal)(info->mx-1);  dhy = (PetscReal)(info->my-1);
715:   hx = 1.0/dhx;                   hy = 1.0/dhy;
716:   hxdhy = hx*dhy;                 hydhx = hy*dhx;

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

720:   /* Test whether we are on the bottom edge of the global array */
721:   if (yints == 0) {
722:     j = 0;
723:     yints = yints + 1;
724:     /* bottom edge */
725:     for (i=info->xs; i<info->xs+info->xm; i++) {
726:       f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
727:     }
728:   }

730:   /* Test whether we are on the top edge of the global array */
731:   if (yinte == info->my) {
732:     j = info->my - 1;
733:     yinte = yinte - 1;
734:     /* top edge */
735:     for (i=info->xs; i<info->xs+info->xm; i++) {
736:       f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
737:     }
738:   }

740:   /* Test whether we are on the left edge of the global array */
741:   if (xints == 0) {
742:     i = 0;
743:     xints = xints + 1;
744:     /* left edge */
745:     for (j=info->ys; j<info->ys+info->ym; j++) {
746:       f[j][i].temp  = x[j][i].temp;
747:     }
748:   }

750:   /* Test whether we are on the right edge of the global array */
751:   if (xinte == info->mx) {
752:     i = info->mx - 1;
753:     xinte = xinte - 1;
754:     /* right edge */
755:     for (j=info->ys; j<info->ys+info->ym; j++) {
756:       f[j][i].temp  = x[j][i].temp - (PetscReal)(grashof>0);
757:     }
758:   }

760:   /* Compute over the interior points */
761:   for (j=yints; j<yinte; j++) {
762:     for (i=xints; i<xinte; i++) {
763:       /* convective coefficients for upwinding */
764:       if (user->COMPOSITE_MODEL){
765:         vx = solu1[j][i].u; vy = solu1[j][i].v;
766:       } else {
767:         vx = solu[j][i].u; vy = solu[j][i].v;
768:       }
769:       avx = PetscAbsScalar(vx);
770:       vxp = .5*(vx+avx); vxm = .5*(vx-avx);
771:       avy = PetscAbsScalar(vy);
772:       vyp = .5*(vy+avy);
773:       vym = .5*(vy-avy);

775:       /* Temperature */
776:       u             = x[j][i].temp;
777:       uxx           = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
778:       uyy           = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
779:       f[j][i].temp  =  uxx + uyy
780:                       + prandtl * ((vxp*(u - x[j][i-1].temp) + vxm*(x[j][i+1].temp - u)) * hy
781:                       + (vyp*(u - x[j-1][i].temp) + vym*(x[j+1][i].temp - u)) * hx);
782:     }
783:   }
784:   return(0);
785: }

789: /* 
790:    FormFunctionComp  - Unwraps the input vector and passes its local ghosted pieces into the user function
791: */
792: PetscErrorCode FormFunctionComp(SNES snes,Vec X,Vec F,void *ctx)
793: {
795:   DMMG           dmmg = (DMMG)ctx;
796:   AppCtx         *user = (AppCtx*)dmmg->user;
797:   DM             dm = dmmg->dm;
798:   DMDALocalInfo    info1,info2;
799:   DM             da1,da2;
800:   Field1         **x1,**f1;
801:   Field2         **x2,**f2;
802:   Vec            X1,X2,F1,F2;

805:   DMCompositeGetEntries(dm,&da1,&da2);
806:   DMDAGetLocalInfo(da1,&info1);
807:   DMDAGetLocalInfo(da2,&info2);

809:   /* Get local vectors to hold ghosted parts of X */
810:   DMCompositeGetLocalVectors(dm,&X1,&X2);
811:   DMCompositeScatter(dm,X,X1,X2);

813:   /* Access the arrays inside the subvectors of X */
814:   DMDAVecGetArray(da1,X1,(void**)&x1);
815:   DMDAVecGetArray(da2,X2,(void**)&x2);

817:   /* Access the subvectors in F. 
818:      These are not ghosted and directly access the memory locations in F */
819:   DMCompositeGetAccess(dm,F,&F1,&F2);

821:   /* Access the arrays inside the subvectors of F */
822:   DMDAVecGetArray(da1,F1,(void**)&f1);
823:   DMDAVecGetArray(da2,F2,(void**)&f2);

825:   /* Evaluate local user provided function */
826:   user->x2 = x2; /* used by FormFunctionLocal1() */
827:   FormFunctionLocal1(&info1,x1,f1,(void**)user);
828:   user->x1 = x1; /* used by FormFunctionLocal2() */
829:   FormFunctionLocal2(&info2,x2,f2,(void**)user);

831:   DMDAVecRestoreArray(da1,F1,(void**)&f1);
832:   DMDAVecRestoreArray(da2,F2,(void**)&f2);
833:   DMCompositeRestoreAccess(dm,F,&F1,&F2);
834:   DMDAVecRestoreArray(da1,X1,(void**)&x1);
835:   DMDAVecRestoreArray(da2,X2,(void**)&x2);
836:   DMCompositeRestoreLocalVectors(dm,&X1,&X2);
837:   return(0);
838: }

842: /* 
843:    Visualize solutions
844: */
845: PetscErrorCode MySolutionView(MPI_Comm comm,PetscInt phy_num,void *ctx)
846: {
848:   AppCtx         *user = (AppCtx*)ctx;
849:   DMMG           *dmmg = user->dmmg;
850:   DM             da=DMMGGetDM(dmmg);
851: #if !defined(PETSC_USE_COMPLEX)
852:   Field          **x = user->x;
853: #endif
854:   PetscInt       i,j,xs,ys,xm,ym;
855:   PetscMPIInt    size;

858: #if defined(PETSC_USE_COMPLEX)
859:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This routine does not support for scalar type complex yet\n");
860: #endif
861:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
862:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
863: 
864:   switch (phy_num){
865:   case 0:
866:     if (size == 1){
867:       PetscPrintf(PETSC_COMM_SELF,"Original Physics %d U,V,Omega,Temp: \n",phy_num);
868:       PetscPrintf(PETSC_COMM_SELF,"-----------------------------------\n");
869:       for (j=ys; j<ys+ym; j++) {
870:         for (i=xs; i<xs+xm; i++) {
871: #if !defined(PETSC_USE_COMPLEX)
872:           PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g, %g, %g, %g\n",j,i,x[j][i].u,x[j][i].v,x[j][i].omega,x[j][i].temp);
873: #endif
874:         }
875:       }
876:     }
877:     break;
878:   case 1:
879:     if (size == 1){
880:       DMMG *dmmg1=user->dmmg1;
881:       Vec  solu_true = DMMGGetx(dmmg1);
882:       DM   da=DMMGGetDM(dmmg1);
883:       Field1 **x1;
884:       PetscPrintf(PETSC_COMM_SELF,"SubPhysics %d: U,V,Omega: \n",phy_num);
885:       PetscPrintf(PETSC_COMM_SELF,"------------------------\n");
886:       DMDAVecGetArray(da,solu_true,&x1);
887:       for (j=ys; j<ys+ym; j++) {
888:         for (i=xs; i<xs+xm; i++) {
889: #if !defined(PETSC_USE_COMPLEX)
890:           PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g, %g, %g\n",j,i,x1[j][i].u,x1[j][i].v,x1[j][i].omega);
891: #endif
892:         }
893:       }
894:       DMDAVecRestoreArray(da,solu_true,&x1);
895:     }
896:     break;
897:   case 2:
898:     if (size == 1){
899:       DMMG *dmmg2=user->dmmg2;
900:       Vec  solu_true = DMMGGetx(dmmg2);
901:       DM   da=DMMGGetDM(dmmg2);
902:       Field2 **x2;
903:       PetscPrintf(PETSC_COMM_SELF,"SubPhysics %d: Temperature: \n",phy_num);
904:       PetscPrintf(PETSC_COMM_SELF,"--------------------------\n");
905:       DMDAVecGetArray(da,solu_true,&x2);
906:       for (j=ys; j<ys+ym; j++) {
907:         for (i=xs; i<xs+xm; i++) {
908: #if !defined(PETSC_USE_COMPLEX)
909:           PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g\n",j,i,x2[j][i].temp);
910: #endif
911:         }
912:       }
913:       DMDAVecRestoreArray(da,solu_true,&x2);
914:     }
915:     break;
916:   default:
917:     if (size == 1){
918:       DMMG        *dmmg_comp=user->dmmg_comp;
919:       DM          da1,da2,da=DMMGGetDM(dmmg);
920:       Vec         X1,X2,solu_true = DMMGGetx(dmmg);
921:       Field       **x;
922:       Field1      **x1;
923:       Field2      **x2;
924:       DM          dm = (*dmmg_comp)->dm;
925:       PetscReal   err,err_tmp;
926:       if (phy_num == 3){
927:         PetscPrintf(PETSC_COMM_SELF,"Composite physics %d, U,V,Omega,Temp: \n",phy_num);
928:         PetscPrintf(PETSC_COMM_SELF,"------------------------------------\n");
929:         PetscPrintf(PETSC_COMM_SELF,"Composite physics, U,V,Omega,Temp: \n");
930:       }
931:       DMDAVecGetArray(da,solu_true,&x);
932:       DMCompositeGetEntries(dm,&da1,&da2);
933:       DMCompositeGetLocalVectors(dm,&X1,&X2);
934:       DMDAVecGetArray(da1,X1,(void**)&x1);
935:       DMDAVecGetArray(da2,X2,(void**)&x2);

937:       err = 0.0;
938:       for (j=ys; j<ys+ym; j++) {
939:         for (i=xs; i<xs+xm; i++) {
940:           err_tmp = PetscAbsScalar(x[j][i].u-x1[j][i].u) + PetscAbsScalar(x[j][i].v-x1[j][i].v) + PetscAbsScalar(x[j][i].omega-x1[j][i].omega);
941:           err_tmp += PetscAbsScalar(x[j][i].temp-x2[j][i].temp);
942:           if (err < err_tmp) err = err_tmp;
943:           if (phy_num == 3){
944: #if !defined(PETSC_USE_COMPLEX)
945:             PetscPrintf(PETSC_COMM_SELF,"x[%d,%d] = %g, %g, %g, %g\n",j,i,x1[j][i].u,x1[j][i].v,x1[j][i].omega,x2[j][i].temp);
946: #endif
947:           }
948:         }
949:       }
950:       PetscPrintf(PETSC_COMM_SELF,"|solu - solu_comp| =  %g\n",err);
951:       DMDAVecRestoreArray(da1,X1,(void**)&x1);
952:       DMDAVecRestoreArray(da2,X2,(void**)&x2);
953:       DMCompositeRestoreLocalVectors(dm,&X1,&X2);
954:       DMDAVecRestoreArray(da,solu_true,&x);
955:     }
956:   }
957:   return(0);
958: }