Actual source code: ex29.c

  1: /*$Id: $*/

  3: /*#define HAVE_DA_HDF*/

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

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

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

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

 27: /* ------------------------------------------------------------------------

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

 31:   ------------------------------------------------------------------------- */

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

 48: #ifdef HAVE_DA_HDF
 49: int DAVecHDFOutput(DA da, Vec X, char *fname);
 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.*M_PI)
 62: #define ly            (4.*M_PI)
 63: #define sqr(a)        ((a)*(a))

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

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

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

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

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

 92: typedef struct {
 93:   int          mglevels;
 94:   int          cycles;           /* number of time steps for integration */
 95:   PassiveReal  nu,eta,d_e,rho_s; /* physical parameters */
 96:   PetscTruth   draw_contours;    /* flag - 1 indicates drawing contours */
 97:   PetscTruth   second_order;
 98:   PetscTruth   PreLoading;
 99: } Parameter;

101: typedef struct {
102:   Vec          Xoldold,Xold,func;
103:   TstepCtx     *tsCtx;
104:   Parameter    *param;
105: } AppCtx;

107: extern int FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
108: extern int Update(DMMG *);
109: extern int Initialize(DMMG *);
110: extern int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user);
111: extern int AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user);
112: extern int Gnuplot(DA da, Vec X, double time);
113: extern int AttachNullSpace(PC,Vec);
114: extern int FormFunctionLocali(DALocalInfo*,MatStencil*,Field**,PetscScalar*,void*);

118: int main(int argc,char **argv)
119: {
120:   DMMG       *dmmg;       /* multilevel grid structure */
121:   AppCtx     *user;       /* user-defined work context (one for each level) */
122:   TstepCtx   tsCtx;       /* time-step parameters (one total) */
123:   Parameter  param;       /* physical parameters (one total) */
124:   int        i,ierr,m,n,mx,my;
125:   MPI_Comm   comm;
126:   DA         da;
127:   PetscTruth flg;
128:   PetscReal  dt_ratio;

130:   int       dfill[16] = {1,0,1,0,
131:                          0,1,0,1,
132:                          0,0,1,1,
133:                          0,1,1,1};
134:   int       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);
140:   comm = PETSC_COMM_WORLD;


143:   PreLoadBegin(PETSC_TRUE,"SetUp");

145:   param.PreLoading = PreLoading;
146:     DMMGCreate(comm,1,&user,&dmmg);
147:     param.mglevels = DMMGGetLevels(dmmg);

149:     /*
150:       Create distributed array multigrid object (DMMG) to manage
151:       parallel grid and vectors for principal unknowns (x) and
152:       governing residuals (f)
153:     */
154:     DACreate2d(comm, DA_XYPERIODIC, DA_STENCIL_STAR, -5, -5,
155:                       PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da);
156: 

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

162:     DMMGSetDM(dmmg,(DM)da);
163:     DADestroy(da);

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

171:     PetscOptionsGetReal(PETSC_NULL, "-viscosity", &param.nu,
172:                                PETSC_NULL);
173: 

175:     PetscOptionsGetReal(PETSC_NULL, "-resistivity", &param.eta,
176:                                PETSC_NULL);
177: 

179:     PetscOptionsGetReal(PETSC_NULL, "-skin_depth", &param.d_e,
180:                                PETSC_NULL);
181: 

183:     PetscOptionsGetReal(PETSC_NULL, "-larmor_radius", &param.rho_s,
184:                                PETSC_NULL);
185: 

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

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

193:     DASetFieldName(DMMGGetDA(dmmg), 0, "phi");
194: 

196:     DASetFieldName(DMMGGetDA(dmmg), 1, "psi");
197: 

199:     DASetFieldName(DMMGGetDA(dmmg), 2, "U");
200: 

202:     DASetFieldName(DMMGGetDA(dmmg), 3, "F");
203: 

205:     /*======================================================================*/
206:     /* Initialize stuff related to time stepping */
207:     /*======================================================================*/

209:     DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,0,0,0,0,0,0,0,0);

211:     tsCtx.fnorm_ini   = 0;
212:     tsCtx.max_steps   = 1000000;
213:     tsCtx.max_time    = 1.0e+12;
214:     /* use for dt = min(dx,dy); multiplied by dt_ratio below */
215:     tsCtx.dt          = PetscMin(lx/mx,ly/my);
216:     tsCtx.fnorm_ratio = 1.0e+10;
217:     tsCtx.t           = 0;
218:     tsCtx.dt_out      = 10;
219:     tsCtx.print_freq  = tsCtx.max_steps;
220:     tsCtx.ts_monitor  = PETSC_FALSE;
221:     tsCtx.dump_time   = -1.0;

223:     PetscOptionsGetInt(PETSC_NULL, "-max_st", &tsCtx.max_steps,
224:                               PETSC_NULL);
225: 

227:     PetscOptionsGetReal(PETSC_NULL, "-ts_rtol", &tsCtx.fnorm_ratio,
228:                                PETSC_NULL);
229: 

231:     PetscOptionsGetInt(PETSC_NULL, "-print_freq", &tsCtx.print_freq,
232:                               PETSC_NULL);
233: 

235:     dt_ratio = 1.0;
236:     PetscOptionsGetReal(PETSC_NULL, "-dt_ratio", &dt_ratio,
237:                                  PETSC_NULL);
238: 
239:     tsCtx.dt *= dt_ratio;


242:     PetscOptionsHasName(PETSC_NULL, "-ts_monitor", &tsCtx.ts_monitor);
243: 

245:     PetscOptionsGetReal(PETSC_NULL, "-dump", &tsCtx.dump_time,
246:                                PETSC_NULL);
247: 

249:     tsCtx.socketviewer = 0;
250:     PetscOptionsHasName(PETSC_NULL, "-socket_viewer", &flg);
251: 
252:     if (flg && !PreLoading) {
253:       tsCtx.ts_monitor = PETSC_TRUE;
254:       PetscViewerSocketOpen(PETSC_COMM_WORLD,0,PETSC_DECIDE,&tsCtx.socketviewer);
255:     }
256: 

258:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259:        Create user context, set problem data, create vector data structures.
260:        Also, compute the initial guess.
261:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
262:     /* create application context for each level */
263:     PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
264:     for (i=0; i<param.mglevels; i++) {
265:       /* create work vectors to hold previous time-step solution and
266:          function value */
267:       VecDuplicate(dmmg[i]->x, &user[i].Xoldold);
268:       VecDuplicate(dmmg[i]->x, &user[i].Xold);
269:       VecDuplicate(dmmg[i]->x, &user[i].func);
270:       user[i].tsCtx = &tsCtx;
271:       user[i].param = &param;
272:       dmmg[i]->user = &user[i];
273:     }
274: 
275:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276:        Create nonlinear solver context
277:        
278:        Process adiC(20):  AddTSTermLocal AddTSTermLocal2 FormFunctionLocal FormFunctionLocali
279:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
280:     DMMGSetSNESLocal(dmmg, FormFunctionLocal, 0,
281:                             ad_FormFunctionLocal, admf_FormFunctionLocal);
282: 
283:     DMMGSetSNESLocali(dmmg,FormFunctionLocali,ad_FormFunctionLocali,admf_FormFunctionLocali);

285:     /* attach nullspace to each level of the preconditioner */
286:     {
287:       SLES       subsles,sles;
288:       PC         pc,subpc;
289:       PetscTruth mg;

291:       SNESGetSLES(DMMGGetSNES(dmmg),&sles);
292:       SLESGetPC(sles,&pc);
293:       AttachNullSpace(pc,DMMGGetx(dmmg));
294:       PetscTypeCompare((PetscObject)pc,PCMG,&mg);
295:       if (mg) {
296:         for (i=0; i<param.mglevels; i++) {
297:           MGGetSmoother(pc,i,&subsles);
298:           SLESGetPC(subsles,&subpc);
299:           AttachNullSpace(subpc,dmmg[i]->x);
300:         }
301:       }
302:     }
303: 
304:     if (PreLoading) {
305:       PetscPrintf(comm, "# viscosity = %g, resistivity = %g, "
306:                          "skin_depth # = %g, larmor_radius # = %g\n",
307:                          param.nu, param.eta, param.d_e, param.rho_s);
308: 
309:       DAGetInfo(DMMGGetDA(dmmg),0,&m,&n,0,0,0,0,0,0,0,0);
310:       PetscPrintf(comm,"Problem size %d by %d\n",m,n);
311:       PetscPrintf(comm,"dx %g dy %g dt %g ratio dt/min(dx,dy) %g\n",lx/mx,ly/my,tsCtx.dt,dt_ratio);
312:     }

314: 
315: 
316:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
317:        Solve the nonlinear system
318:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
319: 
320:     PreLoadStage("Solve");

322:     if (param.draw_contours) {
323:       VecView(((AppCtx*)dmmg[param.mglevels-1]->user)->Xold,
324:                      PETSC_VIEWER_DRAW_WORLD);
325: 
326:     }

328:     Update(dmmg);
329: 

331:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
332:        Free work space.  All PETSc objects should be destroyed when they
333:        are no longer needed.
334:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
335: 
336:     for (i=0; i<param.mglevels; i++) {
337:       VecDestroy(user[i].Xoldold);
338:       VecDestroy(user[i].Xold);
339:       VecDestroy(user[i].func);
340:     }
341:     PetscFree(user);
342:     DMMGDestroy(dmmg);

344:     PreLoadEnd();
345: 
346:   PetscFinalize();
347: 

349:   return 0;
350: }

352: /* ------------------------------------------------------------------- */
355: /* ------------------------------------------------------------------- */
356: int Gnuplot(DA da, Vec X, double time)
357: {
358:   int          i,j,xs,ys,xm,ym;
359:   int          xints,xinte,yints,yinte;
360:   int          ierr;
361:   Field        **x;
362:   FILE         *f;
363:   char         fname[100];
364:   int          cpu;

366:   MPI_Comm_rank(PETSC_COMM_WORLD,&cpu);
367: 

369:   sprintf(fname, "out-%g-%d.dat", time, cpu);

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

373:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
374: 

376:   DAVecGetArray(da,X,&x);
377: 

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

381:   for (j=yints; j<yinte; j++) {
382:     for (i=xints; i<xinte; i++) {
383:       PetscFPrintf(PETSC_COMM_WORLD, f,
384:                           "%d %d %g %g %g %g %g %g\n",
385:                           i, j, 0.0, 0.0,
386:                           x[j][i].U, x[j][i].F, x[j][i].phi, x[j][i].psi);
387: 
388:     }
389:     PetscFPrintf(PETSC_COMM_WORLD,f, "\n");
390: 
391:   }

393:   DAVecRestoreArray(da,X,&x);
394: 

396:   fclose(f);

398:   return 0;
399: }

401: /* ------------------------------------------------------------------- */
404: /* ------------------------------------------------------------------- */
405: int Initialize(DMMG *dmmg)
406: {
407:   AppCtx     *appCtx = (AppCtx*)dmmg[0]->user;
408:   Parameter  *param = appCtx->param;
409:   DA         da;
410:   int        i,j,mx,my,ierr,xs,ys,xm,ym;
411:   PetscReal  two = 2.0,one = 1.0;
412:   PetscReal  hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
413:   PetscReal  d_e,rho_s,de2,xx,yy;
414:   Field      **x, **localx;
415:   Vec        localX;
416:   PetscTruth flg;

419:   PetscOptionsHasName(0,"-restart",&flg);
420:   if (flg) {
421:     PetscViewer viewer;
422:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"binaryoutput",PETSC_BINARY_RDONLY,&viewer);
423:     VecLoadIntoVector(viewer,dmmg[param->mglevels-1]->x);
424:     PetscViewerDestroy(viewer);
425:     return(0);
426:   }

428:   d_e   = param->d_e;
429:   rho_s = param->rho_s;
430:   de2   = sqr(param->d_e);

432:   da   = (DA)(dmmg[param->mglevels-1]->dm);
433:   DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);

435:   dhx   = mx/lx;              dhy = my/ly;
436:   hx    = one/dhx;             hy = one/dhy;
437:   hxdhy = hx*dhy;           hydhx = hy*dhx;
438:   hxhy  = hx*hy;           dhxdhy = dhx*dhy;

440:   /*
441:      Get local grid boundaries (for 2-dimensional DA):
442:        xs, ys   - starting grid indices (no ghost points)
443:        xm, ym   - widths of local grid (no ghost points)
444:   */
445:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

447:   DAGetLocalVector(da,&localX);
448:   /*
449:      Get a pointer to vector data.
450:        - For default PETSc vectors, VecGetArray() returns a pointer to
451:          the data array.  Otherwise, the routine is implementation dependent.
452:        - You MUST call VecRestoreArray() when you no longer need access to
453:          the array.
454:   */
455:   DAVecGetArray(da,dmmg[param->mglevels-1]->x,&x);
456:   DAVecGetArray(da,localX,&localx);

458:   /*
459:      Compute initial solution over the locally owned part of the grid
460:   */
461:   {
462:     PetscReal eps = lx/ly;
463:     PetscReal pert = 1e-4;
464:     PetscReal k = 1.*eps;
465:     PetscReal gam;

467:     if (d_e < rho_s) d_e = rho_s;
468:     gam = k * d_e;

470:     for (j=ys-1; j<ys+ym+1; j++) {
471:       yy = j * hy;
472:       for (i=xs-1; i<xs+xm+1; i++) {
473:         xx = i * hx;

475:         if (xx < -M_PI/2) {
476:           localx[j][i].phi = pert * gam / k * erf((xx + M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
477:         } else if (xx < M_PI/2) {
478:           localx[j][i].phi = - pert * gam / k * erf(xx / (sqrt(2) * d_e)) * (-sin(k*yy));
479:         } else if (xx < 3*M_PI/2){
480:           localx[j][i].phi = pert * gam / k * erf((xx - M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
481:         } else {
482:           localx[j][i].phi = - pert * gam / k * erf((xx - 2.*M_PI) / (sqrt(2) * d_e)) * (-sin(k*yy));
483:         }
484: #ifdef EQ
485:         localx[j][i].psi = 0.;
486: #else
487:         localx[j][i].psi = cos(xx);
488: #endif
489:       }
490:     }
491:     for (j=ys; j<ys+ym; j++) {
492:       for (i=xs; i<xs+xm; i++) {
493:         x[j][i].U   = Lapl(localx,phi,i,j);
494:         x[j][i].F   = localx[j][i].psi - de2 * Lapl(localx,psi,i,j);
495:         x[j][i].phi = localx[j][i].phi;
496:         x[j][i].psi = localx[j][i].psi;
497:       }
498:     }
499:   }
500:   /*
501:      Restore vector
502:   */
503:   DAVecRestoreArray(da,dmmg[param->mglevels-1]->x,&x);
504: 
505:   DAVecRestoreArray(da,localX,&localx);
506: 
507:   DARestoreLocalVector(da,&localX);
508: 

510:   return(0);
511: }

515: int ComputeMaxima(DA da, Vec X, PetscReal t)
516: {
517:   int      ierr,i,j,mx,my,xs,ys,xm,ym;
518:   int      xints,xinte,yints,yinte;
519:   Field    **x;
520:   double   norm[4] = {0,0,0,0};
521:   double   gnorm[4];
522:   MPI_Comm comm;


526:   DAGetInfo(da,0,&mx,&my,0,0,0,0,0,0,0,0);
527: 

529:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
530: 

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

534:   DAVecGetArray(da, X, &x);
535: 

537:   for (j=yints; j<yinte; j++) {
538:     for (i=xints; i<xinte; i++) {
539:       norm[0] = PetscMax(norm[0],x[j][i].U);
540:       norm[1] = PetscMax(norm[1],x[j][i].F);
541:       norm[2] = PetscMax(norm[2],x[j][i].phi);
542:       norm[3] = PetscMax(norm[3],x[j][i].psi);
543:     }
544:   }

546:   DAVecRestoreArray(da,X,&x);
547: 

549:   PetscObjectGetComm((PetscObject)da, &comm);
550: 
551:   MPI_Allreduce(norm, gnorm, 4, MPI_DOUBLE, MPI_MAX, comm);
552: 

554:   PetscFPrintf(PETSC_COMM_WORLD, stderr,
555:                       "%g\t%g\t%g\t%g\t%g\n",
556:                       t, gnorm[0], gnorm[1], gnorm[2], gnorm[3]);
557: 

559:   return(0);
560: }

564: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
565: {
566:   AppCtx        *user = (AppCtx*)ptr;
567:   TstepCtx      *tsCtx = user->tsCtx;
568:   Parameter     *param = user->param;
569:   int           ierr,i,j;
570:   int           xints,xinte,yints,yinte;
571:   PassiveReal   hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
572:   PassiveReal   de2,rhos2,nu,eta,dde2;
573:   PassiveReal   two = 2.0,one = 1.0,p5 = 0.5;
574:   PassiveReal   F_eq_x,By_eq;

576:   PetscScalar   xx;
577:   PetscScalar   vx,vy,avx,avy,vxp,vxm,vyp,vym;
578:   PetscScalar   Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;

581:   de2     = sqr(user->param->d_e);
582:   rhos2   = sqr(user->param->rho_s);
583:   nu      = user->param->nu;
584:   eta     = user->param->eta;
585:   dde2    = one/de2;

587:   /* 
588:      Define mesh intervals ratios for uniform grid.
589:      [Note: FD formulae below are normalized by multiplying through by
590:      local volume element to obtain coefficients O(1) in two dimensions.]
591:   */
592:   dhx   = info->mx/lx;        dhy   = info->my/ly;
593:   hx    = one/dhx;             hy   = one/dhy;
594:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
595:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

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

600:   /* Compute over the interior points */
601:   for (j=yints; j<yinte; j++) {
602:     for (i=xints; i<xinte; i++) {
603: #ifdef EQ
604:       xx = i * hx;
605:       F_eq_x = - (1. + de2) * sin(xx);
606:       By_eq = sin(xx);
607: #else
608:       F_eq_x = 0.;
609:       By_eq = 0.;
610: #endif

612:       /*
613:        * convective coefficients for upwinding
614:        */

616:       vx  = - D_y(x,phi,i,j);
617:       vy  =   D_x(x,phi,i,j);
618:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
619:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
620: #ifndef UPWINDING
621:       vxp = vxm = p5*vx;
622:       vyp = vym = p5*vy;
623: #endif

625:       Bx  =   D_y(x,psi,i,j);
626:       By  = - D_x(x,psi,i,j) + By_eq;
627:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
628:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
629: #ifndef UPWINDING
630:       Bxp = Bxm = p5*Bx;
631:       Byp = Bym = p5*By;
632: #endif

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

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

640:       /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
641:          - nu Lap(U) */
642:       f[j][i].U  =
643:         ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
644:           vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
645:          (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
646:           Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
647:          nu * Lapl(x,U,i,j)) * hxhy;
648: 
649:       /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
650:          - eta * Lap(psi) */
651:       f[j][i].F  =
652:         ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
653:           vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
654:          (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
655:           Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
656:          eta * Lapl(x,psi,i,j)) * hxhy;
657:     }
658:   }

660:   /* Add time step contribution */
661:   if (tsCtx->ires) {
662:     if ((param->second_order) && (tsCtx->itstep > 0))
663:     {
664:       AddTSTermLocal2(info,x,f,user);
665: 
666:     }
667:     else
668:     {
669:       AddTSTermLocal(info,x,f,user);
670: 
671:     }
672:   }

674:   /*
675:      Flop count (multiply-adds are counted as 2 operations)
676:   */
677:   /*  PetscLogFlops(84*info->ym*info->xm); FIXME */
678:   return(0);
679: }

681: /*---------------------------------------------------------------------*/
684: int Update(DMMG *dmmg)
685: /*---------------------------------------------------------------------*/
686: {
687:   AppCtx         *user = (AppCtx *) ((dmmg[0])->user);
688:   TstepCtx       *tsCtx = user->tsCtx;
689:   Parameter      *param = user->param;
690:   SNES           snes;
691:   int            ierr,its,lits,i;
692:   int            max_steps;
693:   int            nfailsCum = 0,nfails = 0;
694:   static int     ic_out;
695:   PetscTruth     ts_monitor = (tsCtx->ts_monitor && !user->param->PreLoading) ? PETSC_TRUE : PETSC_FALSE;


699:   if (user->param->PreLoading)
700:     max_steps = 1;
701:   else
702:     max_steps = tsCtx->max_steps;
703: 
704:   Initialize(dmmg);
705: 

707:   for (tsCtx->itstep = 0; tsCtx->itstep < max_steps; tsCtx->itstep++) {

709:     if ((param->second_order) && (tsCtx->itstep > 0))
710:     {
711:       for (i=param->mglevels-1; i>=0 ;i--)
712:       {
713:         VecCopy(((AppCtx*)dmmg[i]->user)->Xold,
714:                        ((AppCtx*)dmmg[i]->user)->Xoldold);
715: 
716:       }
717:     }

719:     for (i=param->mglevels-1; i>0 ;i--) {
720:       MatRestrict(dmmg[i]->R, dmmg[i]->x, dmmg[i-1]->x);
721: 
722:       VecPointwiseMult(dmmg[i]->Rscale,dmmg[i-1]->x,dmmg[i-1]->x);
723: 
724:       VecCopy(dmmg[i]->x, ((AppCtx*)dmmg[i]->user)->Xold);
725: 
726:     }
727:     VecCopy(dmmg[0]->x, ((AppCtx*)dmmg[0]->user)->Xold);
728: 

730:     DMMGSolve(dmmg);
731: 

733:     snes = DMMGGetSNES(dmmg);


736:     if (tsCtx->itstep == 665000)
737:     {
738:       SLES sles;
739:       PC pc;
740:       Mat mat, pmat;
741:       MatStructure flag;
742:       PetscViewer viewer;
743:       char file[128];

745:       strcpy(file, "matrix");

747:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, file,
748:                                    PETSC_BINARY_CREATE, &viewer);
749: 

751:       SNESGetSLES(snes, &sles);
752: 

754:       SLESGetPC(sles, &pc);
755: 

757:       PCGetOperators(pc, &mat, &pmat, &flag);
758: 

760:       MatView(mat, viewer);
761: 

763:       PetscViewerDestroy(viewer);
764: 
765:       SETERRQ(1,"Done saving Jacobian");
766:     }


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

771:     /* save restart solution if requested at a particular time, then exit */
772:     if (tsCtx->dump_time > 0.0 && tsCtx->t >= tsCtx->dump_time) {
773:       Vec v = DMMGGetx(dmmg);
774:       VecView(v,PETSC_VIEWER_BINARY_WORLD);
775:       SETERRQ1(1,"Saved solution at time %g",tsCtx->t);
776:     }

778:     if (ts_monitor) {
779:       SNESGetIterationNumber(snes,&its);
780:       SNESGetNumberLinearIterations(snes,&lits);
781:       SNESGetNumberUnsuccessfulSteps(snes,&nfails);
782:       SNESGetFunctionNorm(snes,&tsCtx->fnorm);
783:       nfailsCum += nfails;
784:       if (nfailsCum >= 2) SETERRQ(1,"Unable to find a Newton Step");
785:       PetscPrintf(PETSC_COMM_WORLD,"Time Step %d time = %g Newton steps %d linear steps %d fnorm %g\n",
786:                          tsCtx->itstep, tsCtx->t, its, lits, tsCtx->fnorm);
787: 

789:       /* send solution over to Matlab, to be visualized (using ex29.m) */
790:       if (!user->param->PreLoading && tsCtx->socketviewer){
791:         Vec v;
792:         SNESGetSolution(snes,&v);
793:         VecView(v,tsCtx->socketviewer);
794:       }
795:     }

797:     if (!param->PreLoading) {
798:       if (param->draw_contours) {
799:         VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
800: 
801:       }

803:       if (0 && ts_monitor) {
804:         /* compute maxima */
805:         ComputeMaxima((DA) dmmg[param->mglevels-1]->dm,
806:                       dmmg[param->mglevels-1]->x,
807:                       tsCtx->t);
808:         /* output */
809:         if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
810: #ifdef HAVE_DA_HDF
811:           char fname[128];

813:           sprintf(fname, "out-%g.hdf", tsCtx->t);
814:           DAVecHDFOutput(DMMGGetDA(dmmg), DMMGGetx(dmmg), fname);
815: 
816: #else
817: /*
818:           Gnuplot(DMMGGetDA(dmmg), DMMGGetx(dmmg), tsCtx->t);
819:           
820: */
821: #endif
822:           ic_out = 1;
823:         }
824:       }
825:     }
826:   } /* End of time step loop */
827: 
828:   if (!user->param->PreLoading){
829:     SNESGetFunctionNorm(snes,&tsCtx->fnorm);
830:     PetscPrintf(PETSC_COMM_WORLD, "timesteps %d fnorm = %g\n",
831:                        tsCtx->itstep, tsCtx->fnorm);
832: 
833:   }

835:   if (user->param->PreLoading) {
836:     tsCtx->fnorm_ini = 0.0;
837:   }
838: 
839:   return(0);
840: }

842: /*---------------------------------------------------------------------*/
845: int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
846: /*---------------------------------------------------------------------*/
847: {
848:   TstepCtx       *tsCtx = user->tsCtx;
849:   DA             da = info->da;
850:   int            ierr,i,j;
851:   int            xints,xinte,yints,yinte;
852:   PassiveReal    hx,hy,dhx,dhy,hxhy;
853:   PassiveReal    one = 1.0;
854:   PassiveScalar  dtinv;
855:   PassiveField   **xold;


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

862:   dhx  = info->mx/lx;            dhy = info->my/ly;
863:   hx   = one/dhx;                 hy = one/dhy;
864:   hxhy = hx*hy;

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

868:   DAVecGetArray(da,user->Xold,&xold);
869:   for (j=yints; j<yinte; j++) {
870:     for (i=xints; i<xinte; i++) {
871:       f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
872:       f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
873:     }
874:   }
875:   DAVecRestoreArray(da,user->Xold,&xold);

877:   return(0);
878: }

880: /*---------------------------------------------------------------------*/
883: int AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
884: /*---------------------------------------------------------------------*/
885: {
886:   TstepCtx       *tsCtx = user->tsCtx;
887:   DA             da = info->da;
888:   int            ierr,i,j;
889:   int            xints,xinte,yints,yinte;
890:   PassiveReal    hx,hy,dhx,dhy,hxhy;
891:   PassiveReal    one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
892:   PassiveScalar  dtinv;
893:   PassiveField   **xoldold,**xold;


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

900:   dhx  = info->mx/lx;            dhy = info->my/ly;
901:   hx   = one/dhx;                 hy = one/dhy;
902:   hxhy = hx*hy;

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

906:   DAVecGetArray(da,user->Xoldold,&xoldold);
907:   DAVecGetArray(da,user->Xold,&xold);
908:   for (j=yints; j<yinte; j++) {
909:     for (i=xints; i<xinte; i++) {
910:       f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
911:                             p5 * xoldold[j][i].U);
912:       f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
913:                             p5 * xoldold[j][i].F);
914:     }
915:   }
916:   DAVecRestoreArray(da,user->Xoldold,&xoldold);
917:   DAVecRestoreArray(da,user->Xold,&xold);

919:   return(0);
920: }

922: /* -------------------------------------------------------------------------*/
925: int AttachNullSpace(PC pc,Vec model)
926: {
927:   int          i,ierr,rstart,rend,N;
928:   MatNullSpace sp;
929:   Vec          v,vs[1];
930:   PetscScalar  *vx,scale;

933:   VecDuplicate(model,&v);
934:   VecGetSize(model,&N);
935:   scale = 2.0/sqrt(N);
936:   VecGetArray(v,&vx);
937:   VecGetOwnershipRange(v,&rstart,&rend);
938:   for (i=rstart; i<rend; i++) {
939:     if (!(i % 4)) vx[i-rstart] = scale;
940:     else          vx[i-rstart] = 0.0;
941:   }
942:   VecRestoreArray(v,&vx);
943:   vs[0] = v;
944:   MatNullSpaceCreate(PETSC_COMM_WORLD,0,1,vs,&sp);
945:   VecDestroy(v);
946:   PCNullSpaceAttach(pc,sp);
947:   MatNullSpaceDestroy(sp);
948:   return(0);
949: }

951: /*
952:     This is an experimental function and can be safely ignored.
953: */
954: int FormFunctionLocali(DALocalInfo *info,MatStencil *st,Field **x,PetscScalar *f,void *ptr)
955:  {
956:   AppCtx        *user = (AppCtx*)ptr;
957:   TstepCtx      *tsCtx = user->tsCtx;
958:   int           ierr,i,j,c;
959:   int           xints,xinte,yints,yinte;
960:   PassiveReal   hx,hy,dhx,dhy,hxdhy,hydhx,hxhy,dhxdhy;
961:   PassiveReal   de2,rhos2,nu,eta,dde2;
962:   PassiveReal   two = 2.0,one = 1.0,p5 = 0.5;
963:   PassiveReal   F_eq_x,By_eq,dtinv;

965:   PetscScalar   xx;
966:   PetscScalar   vx,vy,avx,avy,vxp,vxm,vyp,vym;
967:   PetscScalar   Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
968:   PassiveField  **xold;

971:   de2     = sqr(user->param->d_e);
972:   rhos2   = sqr(user->param->rho_s);
973:   nu      = user->param->nu;
974:   eta     = user->param->eta;
975:   dde2    = one/de2;

977:   /* 
978:      Define mesh intervals ratios for uniform grid.
979:      [Note: FD formulae below are normalized by multiplying through by
980:      local volume element to obtain coefficients O(1) in two dimensions.]
981:   */
982:   dhx   = info->mx/lx;        dhy   = info->my/ly;
983:   hx    = one/dhx;             hy   = one/dhy;
984:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
985:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

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


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

993: #ifdef EQ
994:       xx = i * hx;
995:       F_eq_x = - (1. + de2) * sin(xx);
996:       By_eq = sin(xx);
997: #else
998:       F_eq_x = 0.;
999:       By_eq = 0.;
1000: #endif

1002:       /*
1003:        * convective coefficients for upwinding
1004:        */

1006:       vx  = - D_y(x,phi,i,j);
1007:       vy  =   D_x(x,phi,i,j);
1008:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
1009:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
1010: #ifndef UPWINDING
1011:       vxp = vxm = p5*vx;
1012:       vyp = vym = p5*vy;
1013: #endif

1015:       Bx  =   D_y(x,psi,i,j);
1016:       By  = - D_x(x,psi,i,j) + By_eq;
1017:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
1018:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
1019: #ifndef UPWINDING
1020:       Bxp = Bxm = p5*Bx;
1021:       Byp = Bym = p5*By;
1022: #endif

1024:       DAVecGetArray(info->da,user->Xold,&xold);
1025:       dtinv = hxhy/(tsCtx->dt);
1026:       switch(c) {

1028:         case 0:
1029:           /* Lap(phi) - U */
1030:           *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
1031:           break;

1033:         case 1:
1034:           /* psi - d_e^2 * Lap(psi) - F */
1035:           *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
1036:           break;

1038:         case 2:
1039:           /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
1040:             - nu Lap(U) */
1041:           *f  =
1042:             ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
1043:               vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
1044:              (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
1045:               Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
1046:              nu * Lapl(x,U,i,j)) * hxhy;
1047:           *f += dtinv*(x[j][i].U-xold[j][i].U);
1048:           break;

1050:         case 3:
1051:           /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
1052:            - eta * Lap(psi) */
1053:           *f  =
1054:             ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
1055:              vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
1056:             (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
1057:              Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
1058:             eta * Lapl(x,psi,i,j)) * hxhy;
1059:           *f += dtinv*(x[j][i].F-xold[j][i].F);
1060:           break;
1061:       }
1062:       DAVecRestoreArray(info->da,user->Xold,&xold);


1065:   /*
1066:      Flop count (multiply-adds are counted as 2 operations)
1067:   */
1068:   /*  PetscLogFlops(84*info->ym*info->xm); FIXME */
1069:   return(0);
1070: }