Actual source code: ex29.c

  2: /*#define HAVE_DA_HDF*/

  4: /* solve the equations for the perturbed magnetic field only */
  5: #define EQ 

  7: /* turning this on causes instability?!? */
  8: /* #define UPWINDING  */

 10: static char help[] = "Hall MHD with in two dimensions with time stepping and multigrid.\n\n\
 11: -options_file ex29.options\n\
 12: other PETSc options\n\
 13: -resistivity <eta>\n\
 14: -viscosity <nu>\n\
 15: -skin_depth <d_e>\n\
 16: -larmor_radius <rho_s>\n\
 17: -contours : draw contour plots of solution\n\n";

 19: /*T
 20:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
 21:    Concepts: DMDA^using distributed arrays;
 22:    Concepts: multicomponent
 23:    Processors: n
 24: T*/

 26: /* ------------------------------------------------------------------------

 28:     We thank Kai Germaschewski for contributing the FormFunctionLocal()

 30:   ------------------------------------------------------------------------- */

 32: /* 
 33:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 34:    Include "petscsnes.h" so that we can use SNES solvers.  
 35:    Include "petscpcmg.h" to control the multigrid solvers. 
 36:    Note that these automatically include:
 37:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 38:      petscmat.h - matrices
 39:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 40:      petscviewer.h - viewers               petscpc.h  - preconditioners
 41:      petscksp.h   - linear solvers 
 42: */
 43: #include <petscsnes.h>
 44: #include <petscdmda.h>
 45: #include <petscpcmg.h>
 46: #include <petscdmmg.h>

 48: #ifdef HAVE_DA_HDF
 49: PetscInt DMDAVecHDFOutput(DM,Vec,char*);
 50: #endif

 52: #define D_x(x,m,i,j)  (p5 * (x[(j)][(i)+1].m - x[(j)][(i)-1].m) * dhx)
 53: #define D_xm(x,m,i,j) ((x[(j)][(i)].m - x[(j)][(i)-1].m) * dhx)
 54: #define D_xp(x,m,i,j) ((x[(j)][(i+1)].m - x[(j)][(i)].m) * dhx)
 55: #define D_y(x,m,i,j)  (p5 * (x[(j)+1][(i)].m - x[(j)-1][(i)].m) * dhy)
 56: #define D_ym(x,m,i,j) ((x[(j)][(i)].m - x[(j)-1][(i)].m) * dhy)
 57: #define D_yp(x,m,i,j) ((x[(j)+1][(i)].m - x[(j)][(i)].m) * dhy)
 58: #define D_xx(x,m,i,j) ((x[(j)][(i)+1].m - two*x[(j)][(i)].m + x[(j)][(i)-1].m) * hydhx * dhxdhy)
 59: #define D_yy(x,m,i,j) ((x[(j)+1][(i)].m - two*x[(j)][(i)].m + x[(j)-1][(i)].m) * hxdhy * dhxdhy)
 60: #define Lapl(x,m,i,j) (D_xx(x,m,i,j) + D_yy(x,m,i,j))
 61: #define lx            (2.*3.1415926)
 62: #define ly            (4.*3.1415926)
 63: #define sqr(a)        ((a)*(a))

 65: /* 
 66:    User-defined routines and data structures
 67: */

 69: typedef struct {
 70:   PetscReal      fnorm_ini,dt_ini;
 71:   PetscReal      dt,dt_out;
 72:   PetscReal      ptime;
 73:   PetscReal      max_time;
 74:   PetscReal      fnorm_ratio;
 75:   PetscInt       ires,itstep;
 76:   PetscInt       max_steps,print_freq;
 77:   PetscReal      t;
 78:   PetscReal      fnorm;

 80:   PetscBool      ts_monitor;           /* print information about each time step */
 81:   PetscReal      dump_time;            /* time to dump solution to a file */
 82:   PetscViewer    socketviewer;         /* socket to dump solution at each timestep for visualization */
 83: } TstepCtx;

 85: typedef struct {
 86:   PetscScalar phi,psi,U,F;
 87: } Field;

 89: typedef struct {
 90:   PassiveScalar phi,psi,U,F;
 91: } PassiveField;

 93: typedef struct {
 94:   PetscInt     mglevels;
 95:   PetscInt     cycles;           /* number of time steps for integration */
 96:   PassiveReal  nu,eta,d_e,rho_s; /* physical parameters */
 97:   PetscBool    draw_contours;    /* flag - 1 indicates drawing contours */
 98:   PetscBool    second_order;
 99:   PetscBool    PetscPreLoading;
100: } Parameters;

102: typedef struct {
103:   Vec          Xoldold,Xold;
104:   TstepCtx     *tsCtx;
105:   Parameters    *param;
106: } AppCtx;


120: int main(int argc,char **argv)
121: {
123:   DMMG           *dmmg;       /* multilevel grid structure */
124:   AppCtx         *user;       /* user-defined work context (one for each level) */
125:   TstepCtx       tsCtx;       /* time-step parameters (one total) */
126:   Parameters      param;       /* physical parameters (one total) */
127:   PetscInt       i,m,n,mx,my;
128:   DM             da;
129:   PetscReal      dt_ratio;
130:   PetscInt       dfill[16] = {1,0,1,0,
131:                               0,1,0,1,
132:                               0,0,1,1,
133:                               0,1,1,1};
134:   PetscInt       ofill[16] = {1,0,0,0,
135:                               0,1,0,0,
136:                               1,1,1,1,
137:                               1,1,1,1};

139:   PetscInitialize(&argc,&argv,(char *)0,help);


142:   PetscPreLoadBegin(PETSC_TRUE,"SetUp");

144:   param.PetscPreLoading = PetscPreLoading;
145:     DMMGCreate(PETSC_COMM_WORLD,1,&user,&dmmg);
146:     param.mglevels = DMMGGetLevels(dmmg);

148:     /*
149:       Create distributed array multigrid object (DMMG) to manage
150:       parallel grid and vectors for principal unknowns (x) and
151:       governing residuals (f)
152:     */
153:     DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_PERIODIC, DMDA_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, -5, -5,
154:                       PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da);

156:     /* overwrite the default sparsity pattern toone specific for
157:        this code's nonzero structure */
158:     DMDASetBlockFills(da,dfill,ofill);

160:     DMMGSetDM(dmmg,(DM)da);
161:     DMDestroy(&da);

163:     /* default physical parameters */
164:     param.nu    = 0;
165:     param.eta   = 1e-3;
166:     param.d_e   = 0.2;
167:     param.rho_s = 0;
168:     param.second_order = PETSC_FALSE;

170:     PetscOptionsGetReal(PETSC_NULL, "-viscosity", &param.nu,PETSC_NULL);

172:     PetscOptionsGetReal(PETSC_NULL, "-resistivity", &param.eta,PETSC_NULL);

174:     PetscOptionsGetReal(PETSC_NULL, "-skin_depth", &param.d_e,PETSC_NULL);

176:     PetscOptionsGetReal(PETSC_NULL, "-larmor_radius", &param.rho_s,PETSC_NULL);

178:     PetscOptionsHasName(PETSC_NULL, "-contours", &param.draw_contours);

180:     PetscOptionsHasName(PETSC_NULL, "-second_order", &param.second_order);

182:     DMDASetFieldName(DMMGGetDM(dmmg), 0, "phi");

184:     DMDASetFieldName(DMMGGetDM(dmmg), 1, "psi");

186:     DMDASetFieldName(DMMGGetDM(dmmg), 2, "U");

188:     DMDASetFieldName(DMMGGetDM(dmmg), 3, "F");

190:     /*======================================================================*/
191:     /* Initialize stuff related to time stepping */
192:     /*======================================================================*/

194:     DMDAGetInfo(DMMGGetDM(dmmg),0,&mx,&my,0,0,0,0,0,0,0,0,0,0);

196:     tsCtx.fnorm_ini   = 0;
197:     tsCtx.max_steps   = 1000000;
198:     tsCtx.max_time    = 1.0e+12;
199:     /* use for dt = min(dx,dy); multiplied by dt_ratio below */
200:     tsCtx.dt          = PetscMin(lx/mx,ly/my);
201:     tsCtx.fnorm_ratio = 1.0e+10;
202:     tsCtx.t           = 0;
203:     tsCtx.dt_out      = 10;
204:     tsCtx.print_freq  = tsCtx.max_steps;
205:     tsCtx.ts_monitor  = PETSC_FALSE;
206:     tsCtx.dump_time   = -1.0;

208:     PetscOptionsGetInt(PETSC_NULL, "-max_st", &tsCtx.max_steps,PETSC_NULL);
209:     PetscOptionsGetReal(PETSC_NULL, "-ts_rtol", &tsCtx.fnorm_ratio,PETSC_NULL);
210:     PetscOptionsGetInt(PETSC_NULL, "-print_freq", &tsCtx.print_freq,PETSC_NULL);

212:     dt_ratio = 1.0;
213:     PetscOptionsGetReal(PETSC_NULL, "-dt_ratio", &dt_ratio,PETSC_NULL);
214:     tsCtx.dt *= dt_ratio;

216:     PetscOptionsHasName(PETSC_NULL, "-ts_monitor", &tsCtx.ts_monitor);
217:     PetscOptionsGetReal(PETSC_NULL, "-dump", &tsCtx.dump_time,PETSC_NULL);


220:     tsCtx.socketviewer = 0;
221: #if defined(PETSC_USE_SOCKET_VIEWER)
222:     {
223:       PetscBool  flg;
224:       PetscOptionsHasName(PETSC_NULL, "-socket_viewer", &flg);
225:       if (flg && !PetscPreLoading) {
226:         tsCtx.ts_monitor = PETSC_TRUE;
227:         PetscViewerSocketOpen(PETSC_COMM_WORLD,0,PETSC_DECIDE,&tsCtx.socketviewer);
228:       }
229:     }
230: #endif

232:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233:        Create user context, set problem data, create vector data structures.
234:        Also, compute the initial guess.
235:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
236:     /* create application context for each level */
237:     PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
238:     for (i=0; i<param.mglevels; i++) {
239:       /* create work vectors to hold previous time-step solution and
240:          function value */
241:       VecDuplicate(dmmg[i]->x, &user[i].Xoldold);
242:       VecDuplicate(dmmg[i]->x, &user[i].Xold);
243:       user[i].tsCtx = &tsCtx;
244:       user[i].param = &param;
245:       DMMGSetUser(dmmg,i,&user[i]);
246:     }
247: 
248:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249:        Create nonlinear solver context
250:        
251:        Process adiC(20):  AddTSTermLocal AddTSTermLocal2 FormFunctionLocal FormFunctionLocali
252:        Process blockadiC(4):  FormFunctionLocali4
253:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254:     DMMGSetSNESLocal(dmmg, FormFunctionLocal, 0,ad_FormFunctionLocal, admf_FormFunctionLocal);
255:     DMMGSetFromOptions(dmmg);
256:     DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,admf_FormFunctionLocali);
257:     DMMGSetSNESLocalib(dmmg,FormFunctionLocali4,0,admfb_FormFunctionLocali4);
258: 
259:    /* attach nullspace to each level of the preconditioner */
260:     DMMGSetNullSpace(dmmg,PETSC_FALSE,1,CreateNullSpace);
261: 
262:     PetscPrintf(PETSC_COMM_WORLD, "finish setupNull!");

264:     if (PetscPreLoading) {
265:       PetscPrintf(PETSC_COMM_WORLD, "# viscosity = %G, resistivity = %G, "
266:                          "skin_depth # = %G, larmor_radius # = %G\n",
267:                          param.nu, param.eta, param.d_e, param.rho_s);
268:       DMDAGetInfo(DMMGGetDM(dmmg),0,&m,&n,0,0,0,0,0,0,0,0,0,0);
269:       PetscPrintf(PETSC_COMM_WORLD,"Problem size %D by %D\n",m,n);
270:       PetscPrintf(PETSC_COMM_WORLD,"dx %G dy %G dt %G ratio dt/min(dx,dy) %G\n",lx/mx,ly/my,tsCtx.dt,dt_ratio);
271:     }

273: 
274: 
275:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276:        Solve the nonlinear system
277:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
278: 
279:     PetscPreLoadStage("Solve");

281:     if (param.draw_contours) {
282:       VecView(((AppCtx*)DMMGGetUser(dmmg,param.mglevels-1))->Xold,PETSC_VIEWER_DRAW_WORLD);
283:     }

285:     Update(dmmg);

287:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288:        Free work space.  All PETSc objects should be destroyed when they
289:        are no longer needed.
290:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291: 
292:     for (i=0; i<param.mglevels; i++) {
293:       VecDestroy(&user[i].Xoldold);
294:       VecDestroy(&user[i].Xold);
295:     }
296:     PetscFree(user);
297:     DMMGDestroy(dmmg);

299:     PetscPreLoadEnd();
300: 
301:   PetscFinalize();
302:   return 0;
303: }

305: /* ------------------------------------------------------------------- */
308: /* ------------------------------------------------------------------- */
309: PetscErrorCode Gnuplot(DM da, Vec X, double mtime)
310: {
311:   PetscInt       i,j,xs,ys,xm,ym;
312:   PetscInt       xints,xinte,yints,yinte;
314:   Field          **x;
315:   FILE           *f;
316:   char           fname[PETSC_MAX_PATH_LEN];
317:   PetscMPIInt    rank;

319:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

321:   sprintf(fname, "out-%g-%d.dat", mtime, rank);

323:   f = fopen(fname, "w");

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

327:   DMDAVecGetArray(da,X,&x);

329:   xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;

331:   for (j=yints; j<yinte; j++) {
332:     for (i=xints; i<xinte; i++) {
333:       PetscFPrintf(PETSC_COMM_WORLD, f,
334:                           "%D %D %G %G %G %G %G %G\n",
335:                           i, j, 0.0, 0.0,
336:                           PetscAbsScalar(x[j][i].U), PetscAbsScalar(x[j][i].F),
337:                           PetscAbsScalar(x[j][i].phi), PetscAbsScalar(x[j][i].psi));
338:     }
339:     PetscFPrintf(PETSC_COMM_WORLD,f, "\n");
340:   }
341:   DMDAVecRestoreArray(da,X,&x);
342:   fclose(f);
343:   return 0;
344: }

346: /* ------------------------------------------------------------------- */
349: /* ------------------------------------------------------------------- */
350: PetscErrorCode Initialize(DMMG *dmmg)
351: {
352:   AppCtx         *appCtx = (AppCtx*)DMMGGetUser(dmmg,0);
353:   Parameters      *param = appCtx->param;
354:   DM             da;
356:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
357:   PetscReal      two = 2.0,one = 1.0;
358:   PetscReal      hx,hy,dhx,dhy,hxdhy,hydhx,dhxdhy; /* hxhy */
359:   PetscReal      d_e,rho_s,de2,xx,yy;
360:   Field          **x, **localx;
361:   Vec            localX;
362:   PetscBool      flg;

365:   PetscOptionsHasName(0,"-restart",&flg);
366:   if (flg) {
367:     PetscViewer viewer;
368:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binaryoutput",FILE_MODE_READ,&viewer);
369:     VecLoad(dmmg[param->mglevels-1]->x,viewer);
370:     PetscViewerDestroy(&viewer);
371:     return(0);
372:   }

374:   d_e   = param->d_e;
375:   rho_s = param->rho_s;
376:   de2   = sqr(param->d_e);

378:   da   = (dmmg[param->mglevels-1]->dm);
379:   DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);

381:   dhx   = mx/lx;              dhy = my/ly;
382:   hx    = one/dhx;             hy = one/dhy;
383:   hxdhy = hx*dhy;           hydhx = hy*dhx;
384:   /* hxhy  = hx*hy; */      dhxdhy = dhx*dhy;

386:   /*
387:      Get local grid boundaries (for 2-dimensional DMDA):
388:        xs, ys   - starting grid indices (no ghost points)
389:        xm, ym   - widths of local grid (no ghost points)
390:   */
391:   DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

393:   DMGetLocalVector(da,&localX);
394:   /*
395:      Get a pointer to vector data.
396:        - For default PETSc vectors, VecGetArray() returns a pointer to
397:          the data array.  Otherwise, the routine is implementation dependent.
398:        - You MUST call VecRestoreArray() when you no longer need access to
399:          the array.
400:   */
401:   DMDAVecGetArray(da,dmmg[param->mglevels-1]->x,&x);
402:   DMDAVecGetArray(da,localX,&localx);

404:   /*
405:      Compute initial solution over the locally owned part of the grid
406:   */
407: #if defined (PETSC_HAVE_ERF)
408:   {
409:     PetscReal eps = lx/ly;
410:     PetscReal pert = 1e-4;
411:     PetscReal k = 1.*eps;
412:     PetscReal gam;

414:     if (d_e < rho_s) d_e = rho_s;
415:     gam = k * d_e;
416:     for (j=ys-1; j<ys+ym+1; j++) {
417:       yy = j * hy;
418:       for (i=xs-1; i<xs+xm+1; i++) {
419:         xx = i * hx;
420:         if (xx < -3.1416926/2) {
421:           localx[j][i].phi = pert * gam / k * erf((xx + 3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
422:         } else if (xx < 3.1415926/2) {
423:           localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2.0) * d_e)) * (-sin(k*yy));
424:         } else if (xx < 3*3.1415926/2){
425:           localx[j][i].phi = pert * gam / k * erf((xx - 3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
426:         } else {
427:           localx[j][i].phi = - pert * gam / k * erf((xx - 2.*3.1415926) / (sqrt(2.0) * d_e)) * (-sin(k*yy));
428:         }
429: #ifdef EQ
430:         localx[j][i].psi = 0.;
431: #else
432:         localx[j][i].psi = cos(xx);
433: #endif
434:       }
435:     }

437:     for (j=ys; j<ys+ym; j++) {
438:       for (i=xs; i<xs+xm; i++) {
439:         x[j][i].U   = Lapl(localx,phi,i,j);
440:         x[j][i].F   = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
441:         x[j][i].phi = localx[j][i].phi;
442:         x[j][i].psi = localx[j][i].psi;
443:       }
444:     }
445:   }
446: #else
447:   SETERRQ(PETSC_COMM_SELF,1,"erf() not available on this machine");
448: #endif

450:   /*
451:      Restore vector
452:   */
453:   DMDAVecRestoreArray(da,dmmg[param->mglevels-1]->x,&x);
454:   DMDAVecRestoreArray(da,localX,&localx);
455:   DMRestoreLocalVector(da,&localX);

457:   return(0);
458: }

462: PetscErrorCode ComputeMaxima(DM da, Vec X, PetscReal t)
463: {
465:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
466:   PetscInt       xints,xinte,yints,yinte;
467:   Field          **x;
468:   double         norm[4] = {0,0,0,0};
469:   double         gnorm[4];
470:   MPI_Comm       comm;


474:   DMDAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0,0,0);

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

478:   xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;

480:   DMDAVecGetArray(da, X, &x);

482:   for (j=yints; j<yinte; j++) {
483:     for (i=xints; i<xinte; i++) {
484:       norm[0] = PetscMax(norm[0],PetscAbsScalar(x[j][i].U));
485:       norm[1] = PetscMax(norm[1],PetscAbsScalar(x[j][i].F));
486:       norm[2] = PetscMax(norm[2],PetscAbsScalar(x[j][i].phi));
487:       norm[3] = PetscMax(norm[3],PetscAbsScalar(x[j][i].psi));
488:     }
489:   }

491:   DMDAVecRestoreArray(da,X,&x);

493:   PetscObjectGetComm((PetscObject)da, &comm);
494:   MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);

496:   PetscFPrintf(PETSC_COMM_WORLD, stderr,"%G\t%G\t%G\t%G\t%G\n",t, gnorm[0], gnorm[1], gnorm[2], gnorm[3]);

498:   return(0);
499: }

503: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
504: {
505:   AppCtx         *user = (AppCtx*)ptr;
506:   TstepCtx       *tsCtx = user->tsCtx;
507:   Parameters      *param = user->param;
509:   PetscInt       xints,xinte,yints,yinte,i,j;
510:   PassiveReal    hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
511:   PassiveReal    de2,rhos2,nu,eta,dde2;
512:   PassiveReal    two = 2.0,one = 1.0,p5 = 0.5;
513:   PassiveReal    F_eq_x,By_eq;
514:   PetscScalar    xx;
515:   PetscScalar    vx,vy,avx,avy,vxp,vxm,vyp,vym;
516:   PetscScalar    Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;

519:   de2     = sqr(param->d_e);
520:   rhos2   = sqr(param->rho_s);
521:   nu      = param->nu;
522:   eta     = param->eta;
523:   dde2    = one/de2;

525:   /* 
526:      Define mesh intervals ratios for uniform grid.
527:      [Note: FD formulae below are normalized by multiplying through by
528:      local volume element to obtain coefficients O(1) in two dimensions.]
529:   */
530:   dhx   = info->mx/lx;        dhy   = info->my/ly;
531:   hx    = one/dhx;             hy   = one/dhy;
532:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
533:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

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

538:   /* Compute over the interior points */
539:   for (j=yints; j<yinte; j++) {
540:     for (i=xints; i<xinte; i++) {
541: #ifdef EQ
542:       xx = i * hx;
543:       F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
544:       By_eq = sin(PetscAbsScalar(xx));
545: #else
546:       F_eq_x = 0.;
547:       By_eq = 0.;
548: #endif

550:       /*
551:        * convective coefficients for upwinding
552:        */

554:       vx  = - D_y(x,phi,i,j);
555:       vy  =   D_x(x,phi,i,j);
556:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
557:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
558: #ifndef UPWINDING
559:       vxp = vxm = p5*vx;
560:       vyp = vym = p5*vy;
561: #endif

563:       Bx  =   D_y(x,psi,i,j);
564:       By  = - D_x(x,psi,i,j) + By_eq;
565:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
566:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
567: #ifndef UPWINDING
568:       Bxp = Bxm = p5*Bx;
569:       Byp = Bym = p5*By;
570: #endif

572:       /* Lap(phi) - U */
573:       f[j][i].phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;

575:       /* psi - d_e^2 * Lap(psi) - F */
576:       f[j][i].psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;

578:       /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
579:          - nu Lap(U) */
580:       f[j][i].U  =
581:         ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
582:           vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
583:          (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
584:           Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
585:          nu * Lapl(x,U,i,j)) * hxhy;
586: 
587:       /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
588:          - eta * Lap(psi) */
589:       f[j][i].F  =
590:         ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
591:           vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
592:          (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
593:           Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
594:          eta * Lapl(x,psi,i,j)) * hxhy;
595:     }
596:   }

598:   /* Add time step contribution */
599:   if (tsCtx->ires) {
600:     if ((param->second_order) && (tsCtx->itstep > 0)){
601:       AddTSTermLocal2(info,x,f,user);
602:     } else {
603:       AddTSTermLocal(info,x,f,user);
604:     }
605:   }

607:   /*
608:      Flop count (multiply-adds are counted as 2 operations)
609:   */
610:   /*  PetscLogFlops(84.0*info->ym*info->xm); FIXME */
611:   return(0);
612: }

614: /*---------------------------------------------------------------------*/
617: PetscErrorCode Update(DMMG *dmmg)
618: /*---------------------------------------------------------------------*/
619: {
620:   AppCtx          *user = (AppCtx *) DMMGGetUser(dmmg,0);
621:   TstepCtx        *tsCtx = user->tsCtx;
622:   Parameters       *param = user->param;
623:   SNES            snes;
624:   PetscErrorCode  ierr;
625:   PetscInt        its,lits,i;
626:   PetscInt        max_steps;
627:   PetscInt        nfailsCum = 0,nfails = 0;
628:   static PetscInt ic_out;
629:   PetscBool       ts_monitor = (tsCtx->ts_monitor && !param->PetscPreLoading) ? PETSC_TRUE : PETSC_FALSE;


633:   if (param->PetscPreLoading)
634:     max_steps = 1;
635:   else
636:     max_steps = tsCtx->max_steps;
637: 
638:   Initialize(dmmg);

640:   snes = DMMGGetSNES(dmmg);

642:   for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {
643: 
644:     PetscPrintf(PETSC_COMM_WORLD, "time step=%d!\n",tsCtx->itstep);
645:     if ((param->second_order) && (tsCtx->itstep > 0))
646:     {
647:       for (i=param->mglevels-1; i>=0 ;i--)
648:       {
649:         VecCopy(((AppCtx*)DMMGGetUser(dmmg,i))->Xold,((AppCtx*)DMMGGetUser(dmmg,i))->Xoldold);
650:       }
651:     }

653:     for (i=param->mglevels-1; i>0 ;i--) {
654:       MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
655:       VecPointwiseMult(dmmg[i-1]->x,dmmg[i-1]->x,dmmg[i]->Rscale);
656:       VecCopy(dmmg[i]->x, ((AppCtx*)DMMGGetUser(dmmg,i))->Xold);
657:     }
658:     VecCopy(dmmg[0]->x, ((AppCtx*)DMMGGetUser(dmmg,0))->Xold);

660:     DMMGSolve(dmmg);



664:     if (tsCtx->itstep == 665000)
665:     {
666:       KSP          ksp;
667:       PC           pc;
668:       Mat          mat, pmat;
669:       MatStructure flag;
670:       PetscViewer  viewer;
671:       char         file[PETSC_MAX_PATH_LEN];

673:       PetscStrcpy(file, "matrix");

675:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_WRITE, &viewer);

677:       SNESGetKSP(snes, &ksp);

679:       KSPGetPC(ksp, &pc);

681:       PCGetOperators(pc, &mat, &pmat, &flag);

683:       MatView(mat, viewer);

685:       PetscViewerDestroy(&viewer);
686:       SETERRQ(PETSC_COMM_SELF,1,"Done saving Jacobian");
687:     }


690:     tsCtx->t += tsCtx->dt;

692:     /* save restart solution if requested at a particular time, then exit */
693:     if (tsCtx->dump_time > 0.0 && tsCtx->t >= tsCtx->dump_time) {
694:       Vec v = DMMGGetx(dmmg);
695:       VecView(v,PETSC_VIEWER_BINARY_WORLD);
696:       SETERRQ1(PETSC_COMM_SELF,1,"Saved solution at time %G",tsCtx->t);
697:     }

699:     if (ts_monitor)
700:     {
701:       SNESGetIterationNumber(snes, &its);
702:       SNESGetLinearSolveIterations(snes, &lits);
703:       SNESGetNonlinearStepFailures(snes, &nfails);
704:       SNESGetFunctionNorm(snes, &tsCtx->fnorm);

706:       nfailsCum += nfails;
707:       if (nfailsCum >= 2) SETERRQ(PETSC_COMM_SELF,1, "unable to find a newton step");

709:       PetscPrintf(PETSC_COMM_WORLD,
710:                          "time step = %D, time = %G, number of nonlinear steps = %D, "
711:                          "number of linear steps = %D, norm of the function = %G\n",
712:                          tsCtx->itstep + 1, tsCtx->t, its, lits, PetscAbsScalar(tsCtx->fnorm));

714:       /* send solution over to MATLAB, to be visualized (using ex29.m) */
715:       if (!param->PetscPreLoading && tsCtx->socketviewer)
716:       {
717:         Vec v;
718:         SNESGetSolution(snes, &v);
719: #if defined(PETSC_USE_SOCKET_VIEWER)
720:         VecView(v, tsCtx->socketviewer);
721: #endif
722:       }
723:     }

725:     if (!param->PetscPreLoading) {
726:       if (param->draw_contours) {
727:         VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
728:       }

730:       if (1 && ts_monitor) {
731:         /* compute maxima */
732:         ComputeMaxima( dmmg[param->mglevels-1]->dm,dmmg[param->mglevels-1]->x,tsCtx->t);
733:         /* output */
734:         if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
735: #ifdef HAVE_DA_HDF
736:           char fname[PETSC_MAX_PATH_LEN];

738:           sprintf(fname, "out-%G.hdf", tsCtx->t);
739:           DMDAVecHDFOutput(DMMGGetDM(dmmg), DMMGGetx(dmmg), fname);
740: #else
741: /*
742:           Gnuplot(DMMGGetDM(dmmg), DMMGGetx(dmmg), tsCtx->t);
743: */
744: #endif
745:           ic_out = 1;
746:         }
747:       }
748:     }
749:   } /* End of time step loop */
750: 
751:   if (!param->PetscPreLoading){
752:     SNESGetFunctionNorm(snes,&tsCtx->fnorm);
753:     PetscPrintf(PETSC_COMM_WORLD, "timesteps %D fnorm = %G\n",
754:                        tsCtx->itstep, PetscAbsScalar(tsCtx->fnorm));
755:   }

757:   if (param->PetscPreLoading) {
758:     tsCtx->fnorm_ini = 0.0;
759:   }
760: 
761:   return(0);
762: }

764: /*---------------------------------------------------------------------*/
767: PetscErrorCode AddTSTermLocal(DMDALocalInfo* info,Field **x,Field **f,AppCtx *user)
768: /*---------------------------------------------------------------------*/
769: {
770:   TstepCtx       *tsCtx = user->tsCtx;
771:   DM             da = info->da;
773:   PetscInt       i,j;
774:   PetscInt       xints,xinte,yints,yinte;
775:   PassiveReal    hx,hy,dhx,dhy,hxhy;
776:   PassiveReal    one = 1.0;
777:   PassiveScalar  dtinv;
778:   PassiveField   **xold;


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

785:   dhx  = info->mx/lx;            dhy = info->my/ly;
786:   hx   = one/dhx;                 hy = one/dhy;
787:   hxhy = hx*hy;

789:   dtinv = hxhy/(tsCtx->dt);

791:   DMDAVecGetArray(da,user->Xold,&xold);
792:   for (j=yints; j<yinte; j++) {
793:     for (i=xints; i<xinte; i++) {
794:       f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
795:       f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
796:     }
797:   }
798:   DMDAVecRestoreArray(da,user->Xold,&xold);

800:   return(0);
801: }

803: /*---------------------------------------------------------------------*/
806: PetscErrorCode AddTSTermLocal2(DMDALocalInfo* info,Field **x,Field **f,AppCtx *user)
807: /*---------------------------------------------------------------------*/
808: {
809:   TstepCtx       *tsCtx = user->tsCtx;
810:   DM             da = info->da;
812:   PetscInt       i,j,xints,xinte,yints,yinte;
813:   PassiveReal    hx,hy,dhx,dhy,hxhy;
814:   PassiveReal    one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
815:   PassiveScalar  dtinv;
816:   PassiveField   **xoldold,**xold;


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

823:   dhx  = info->mx/lx;            dhy = info->my/ly;
824:   hx   = one/dhx;                 hy = one/dhy;
825:   hxhy = hx*hy;

827:   dtinv = hxhy/(tsCtx->dt);

829:   DMDAVecGetArray(da,user->Xoldold,&xoldold);
830:   DMDAVecGetArray(da,user->Xold,&xold);
831:   for (j=yints; j<yinte; j++) {
832:     for (i=xints; i<xinte; i++) {
833:       f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
834:                             p5 * xoldold[j][i].U);
835:       f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
836:                             p5 * xoldold[j][i].F);
837:     }
838:   }
839:   DMDAVecRestoreArray(da,user->Xoldold,&xoldold);
840:   DMDAVecRestoreArray(da,user->Xold,&xold);

842:   return(0);
843: }

845: /* Creates the null space of the Jacobian for a particular level */
848: PetscErrorCode CreateNullSpace(DMMG dmmg,Vec *nulls)
849: {
851:   PetscInt       i,N,rstart,rend;
852:   PetscScalar    scale,*vx;

855:   VecGetSize(nulls[0],&N);
856:   scale = 2.0/sqrt((PetscReal)N);
857:   VecGetArray(nulls[0],&vx);
858:   VecGetOwnershipRange(nulls[0],&rstart,&rend);
859:   for (i=rstart; i<rend; i++) {
860:     if (!(i % 4)) vx[i-rstart] = scale;
861:     else          vx[i-rstart] = 0.0;
862:   }
863:   VecRestoreArray(nulls[0],&vx);
864:   return(0);
865: }

867: /*
868:     This is an experimental function and can be safely ignored.
869: */
870: PetscErrorCode FormFunctionLocali(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
871: {
872:   AppCtx         *user = (AppCtx*)ptr;
873:   TstepCtx       *tsCtx = user->tsCtx;
874:   Parameters      *param = user->param;
876:   PetscInt       i,j,c;
877:   /* PetscInt       xints,xinte,yints,yinte; */
878:   PassiveReal    hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
879:   PassiveReal    de2,rhos2,nu,eta,dde2;
880:   PassiveReal    two = 2.0,one = 1.0,p5 = 0.5;
881:   PassiveReal    F_eq_x,By_eq,dtinv;
882:   PetscScalar    xx;
883:   PetscScalar    vx,vy,avx,avy,vxp,vxm,vyp,vym;
884:   PetscScalar    Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
885:   PassiveField   **xold;

888:   de2     = sqr(param->d_e);
889:   rhos2   = sqr(param->rho_s);
890:   nu      = param->nu;
891:   eta     = param->eta;
892:   dde2    = one/de2;

894:   /* 
895:      Define mesh intervals ratios for uniform grid.
896:      [Note: FD formulae below are normalized by multiplying through by
897:      local volume element to obtain coefficients O(1) in two dimensions.]
898:   */
899:   dhx   = info->mx/lx;        dhy   = info->my/ly;
900:   hx    = one/dhx;             hy   = one/dhy;
901:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
902:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

904:   /*
905:   xints = info->xs; xinte = info->xs+info->xm;
906:   yints = info->ys; yinte = info->ys+info->ym;
907:    */

909:   i = st->i; j = st->j; c = st->c;

911: #ifdef EQ
912:       xx = i * hx;
913:       F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
914:       By_eq = sin(PetscAbsScalar(xx));
915: #else
916:       F_eq_x = 0.;
917:       By_eq = 0.;
918: #endif

920:       /*
921:        * convective coefficients for upwinding
922:        */

924:       vx  = - D_y(x,phi,i,j);
925:       vy  =   D_x(x,phi,i,j);
926:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
927:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
928: #ifndef UPWINDING
929:       vxp = vxm = p5*vx;
930:       vyp = vym = p5*vy;
931: #endif

933:       Bx  =   D_y(x,psi,i,j);
934:       By  = - D_x(x,psi,i,j) + By_eq;
935:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
936:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
937: #ifndef UPWINDING
938:       Bxp = Bxm = p5*Bx;
939:       Byp = Bym = p5*By;
940: #endif

942:       DMDAVecGetArray(info->da,user->Xold,&xold);
943:       dtinv = hxhy/(tsCtx->dt);
944:       switch(c) {

946:         case 0:
947:           /* Lap(phi) - U */
948:           *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
949:           break;

951:         case 1:
952:           /* psi - d_e^2 * Lap(psi) - F */
953:           *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
954:           break;

956:         case 2:
957:           /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
958:             - nu Lap(U) */
959:           *f  =
960:             ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
961:               vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
962:              (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
963:               Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
964:              nu * Lapl(x,U,i,j)) * hxhy;
965:           *f += dtinv*(x[j][i].U-xold[j][i].U);
966:           break;

968:         case 3:
969:           /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
970:            - eta * Lap(psi) */
971:           *f  =
972:             ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
973:              vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
974:             (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
975:              Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
976:             eta * Lapl(x,psi,i,j)) * hxhy;
977:           *f += dtinv*(x[j][i].F-xold[j][i].F);
978:           break;
979:       }
980:       DMDAVecRestoreArray(info->da,user->Xold,&xold);


983:   /*
984:      Flop count (multiply-adds are counted as 2 operations)
985:   */
986:   /*  PetscLogFlops(84.0*info->ym*info->xm); FIXME */
987:   return(0);
988: }
989: /*
990:     This is an experimental function and can be safely ignored.
991: */
992: PetscErrorCode FormFunctionLocali4(DMDALocalInfo *info,MatStencil *st,Field **x,PetscScalar *ff,void *ptr)
993: {
994:   Field          *f  = (Field *)ff;
995:   AppCtx         *user = (AppCtx*)ptr;
996:   TstepCtx       *tsCtx = user->tsCtx;
997:   Parameters     *param = user->param;
999:   PetscInt       i,j;
1000:   /* PetscInt       xints,xinte,yints,yinte; */
1001:   PassiveReal    hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
1002:   PassiveReal    de2,rhos2,nu,eta,dde2;
1003:   PassiveReal    two = 2.0,one = 1.0,p5 = 0.5;
1004:   PassiveReal    F_eq_x,By_eq,dtinv;
1005:   PetscScalar    xx;
1006:   PetscScalar    vx,vy,avx,avy,vxp,vxm,vyp,vym;
1007:   PetscScalar    Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
1008:   PassiveField   **xold;

1011:   de2     = sqr(param->d_e);
1012:   rhos2   = sqr(param->rho_s);
1013:   nu      = param->nu;
1014:   eta     = param->eta;
1015:   dde2    = one/de2;

1017:   /* 
1018:      Define mesh intervals ratios for uniform grid.
1019:      [Note: FD formulae below are normalized by multiplying through by
1020:      local volume element to obtain coefficients O(1) in two dimensions.]
1021:   */
1022:   dhx   = info->mx/lx;        dhy   = info->my/ly;
1023:   hx    = one/dhx;             hy   = one/dhy;
1024:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
1025:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

1027:   /*
1028:   xints = info->xs; xinte = info->xs+info->xm;
1029:   yints = info->ys; yinte = info->ys+info->ym;
1030:    */

1032:   i = st->i; j = st->j;

1034: #ifdef EQ
1035:       xx = i * hx;
1036:       F_eq_x = - (1. + de2) * sin(PetscAbsScalar(xx));
1037:       By_eq = sin(PetscAbsScalar(xx));
1038: #else
1039:       F_eq_x = 0.;
1040:       By_eq = 0.;
1041: #endif

1043:       /*
1044:        * convective coefficients for upwinding
1045:        */

1047:       vx  = - D_y(x,phi,i,j);
1048:       vy  =   D_x(x,phi,i,j);
1049:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
1050:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
1051: #ifndef UPWINDING
1052:       vxp = vxm = p5*vx;
1053:       vyp = vym = p5*vy;
1054: #endif

1056:       Bx  =   D_y(x,psi,i,j);
1057:       By  = - D_x(x,psi,i,j) + By_eq;
1058:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
1059:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
1060: #ifndef UPWINDING
1061:       Bxp = Bxm = p5*Bx;
1062:       Byp = Bym = p5*By;
1063: #endif

1065:       DMDAVecGetArray(info->da,user->Xold,&xold);
1066:       dtinv = hxhy/(tsCtx->dt);
1067: 

1069: 
1070:           /* Lap(phi) - U */
1071:           f->phi = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
1072: 

1074: 
1075:           /* psi - d_e^2 * Lap(psi) - F */
1076:           f->psi = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
1077: 

1079: 
1080:           /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
1081:             - nu Lap(U) */
1082:           f->U  =
1083:             ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
1084:               vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
1085:              (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
1086:               Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
1087:              nu * Lapl(x,U,i,j)) * hxhy;
1088:           f->U += dtinv*(x[j][i].U-xold[j][i].U);
1089: 

1091: 
1092:           /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
1093:            - eta * Lap(psi) */
1094:           f->F  =
1095:             ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
1096:              vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
1097:             (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
1098:              Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
1099:             eta * Lapl(x,psi,i,j)) * hxhy;
1100:           f->F += dtinv*(x[j][i].F-xold[j][i].F);
1101: 
1102: 
1103:       DMDAVecRestoreArray(info->da,user->Xold,&xold);


1106:   /*
1107:      Flop count (multiply-adds are counted as 2 operations)
1108:   */
1109:   /*  PetscLogFlops(84.0*info->ym*info->xm); FIXME */
1110:   return(0);
1111: }