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:      petscksp.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:       KSP       subksp,ksp;
288:       PC         pc,subpc;
289:       PetscTruth mg;
290:       KSP        ksp;

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

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

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

329:     Update(dmmg);
330: 

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

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

350:   return 0;
351: }

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

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

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

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

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

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

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

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

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

397:   fclose(f);

399:   return 0;
400: }

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

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

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

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

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

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

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

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

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

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

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

511:   return(0);
512: }

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


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

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

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

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

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

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

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

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

560:   return(0);
561: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

731:     DMMGSolve(dmmg);
732: 

734:     snes = DMMGGetSNES(dmmg);


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

747:       PetscStrcpy(file, "matrix");

749:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, file,
750:                                    PETSC_FILE_CREATE, &viewer);
751: 

753:       SNESGetKSP(snes, &ksp);
754: 

756:       KSPGetPC(ksp, &pc);
757: 

759:       PCGetOperators(pc, &mat, &pmat, &flag);
760: 

762:       MatView(mat, viewer);
763: 

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


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

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

780:     if (ts_monitor)
781:     {
782:       SNESGetIterationNumber(snes, &its);
783:       SNESGetNumberLinearIterations(snes, &lits);
784:       SNESGetNumberUnsuccessfulSteps(snes, &nfails);
785:       SNESGetFunctionNorm(snes, &tsCtx->fnorm);

787:       nfailsCum += nfails;
788:       if (nfailsCum >= 2)
789:         SETERRQ(1, "unable to find a newton step");

791:       PetscPrintf(PETSC_COMM_WORLD,
792:                          "time step = %d, time = %g, number of nonlinear steps = %d, "
793:                          "number of linear steps = %d, norm of the function = %g\n",
794:                          tsCtx->itstep + 1, tsCtx->t, its, lits, tsCtx->fnorm);
795: 

797:       /* send solution over to Matlab, to be visualized (using ex29.m) */
798:       if (!user->param->PreLoading && tsCtx->socketviewer)
799:       {
800:         Vec v;
801:         SNESGetSolution(snes, &v);
802:         VecView(v, tsCtx->socketviewer);
803:       }
804:     }

806:     if (!param->PreLoading) {
807:       if (param->draw_contours) {
808:         VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
809: 
810:       }

812:       if (1 && ts_monitor) {
813:         /* compute maxima */
814:         ComputeMaxima((DA) dmmg[param->mglevels-1]->dm,
815:                       dmmg[param->mglevels-1]->x,
816:                       tsCtx->t);
817:         /* output */
818:         if (ic_out++ == (int)(tsCtx->dt_out / tsCtx->dt + .5)) {
819: #ifdef HAVE_DA_HDF
820:           char fname[128];

822:           sprintf(fname, "out-%g.hdf", tsCtx->t);
823:           DAVecHDFOutput(DMMGGetDA(dmmg), DMMGGetx(dmmg), fname);
824: 
825: #else
826: /*
827:           Gnuplot(DMMGGetDA(dmmg), DMMGGetx(dmmg), tsCtx->t);
828:           
829: */
830: #endif
831:           ic_out = 1;
832:         }
833:       }
834:     }
835:   } /* End of time step loop */
836: 
837:   if (!user->param->PreLoading){
838:     SNESGetFunctionNorm(snes,&tsCtx->fnorm);
839:     PetscPrintf(PETSC_COMM_WORLD, "timesteps %d fnorm = %g\n",
840:                        tsCtx->itstep, tsCtx->fnorm);
841: 
842:   }

844:   if (user->param->PreLoading) {
845:     tsCtx->fnorm_ini = 0.0;
846:   }
847: 
848:   return(0);
849: }

851: /*---------------------------------------------------------------------*/
854: int AddTSTermLocal(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
855: /*---------------------------------------------------------------------*/
856: {
857:   TstepCtx       *tsCtx = user->tsCtx;
858:   DA             da = info->da;
859:   int            ierr,i,j;
860:   int            xints,xinte,yints,yinte;
861:   PassiveReal    hx,hy,dhx,dhy,hxhy;
862:   PassiveReal    one = 1.0;
863:   PassiveScalar  dtinv;
864:   PassiveField   **xold;


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

871:   dhx  = info->mx/lx;            dhy = info->my/ly;
872:   hx   = one/dhx;                 hy = one/dhy;
873:   hxhy = hx*hy;

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

877:   DAVecGetArray(da,user->Xold,&xold);
878:   for (j=yints; j<yinte; j++) {
879:     for (i=xints; i<xinte; i++) {
880:       f[j][i].U += dtinv*(x[j][i].U-xold[j][i].U);
881:       f[j][i].F += dtinv*(x[j][i].F-xold[j][i].F);
882:     }
883:   }
884:   DAVecRestoreArray(da,user->Xold,&xold);

886:   return(0);
887: }

889: /*---------------------------------------------------------------------*/
892: int AddTSTermLocal2(DALocalInfo* info,Field **x,Field **f,AppCtx *user)
893: /*---------------------------------------------------------------------*/
894: {
895:   TstepCtx       *tsCtx = user->tsCtx;
896:   DA             da = info->da;
897:   int            ierr,i,j;
898:   int            xints,xinte,yints,yinte;
899:   PassiveReal    hx,hy,dhx,dhy,hxhy;
900:   PassiveReal    one = 1.0, onep5 = 1.5, two = 2.0, p5 = 0.5;
901:   PassiveScalar  dtinv;
902:   PassiveField   **xoldold,**xold;


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

909:   dhx  = info->mx/lx;            dhy = info->my/ly;
910:   hx   = one/dhx;                 hy = one/dhy;
911:   hxhy = hx*hy;

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

915:   DAVecGetArray(da,user->Xoldold,&xoldold);
916:   DAVecGetArray(da,user->Xold,&xold);
917:   for (j=yints; j<yinte; j++) {
918:     for (i=xints; i<xinte; i++) {
919:       f[j][i].U += dtinv * (onep5 * x[j][i].U - two * xold[j][i].U +
920:                             p5 * xoldold[j][i].U);
921:       f[j][i].F += dtinv * (onep5 * x[j][i].F - two * xold[j][i].F +
922:                             p5 * xoldold[j][i].F);
923:     }
924:   }
925:   DAVecRestoreArray(da,user->Xoldold,&xoldold);
926:   DAVecRestoreArray(da,user->Xold,&xold);

928:   return(0);
929: }

931: /* -------------------------------------------------------------------------*/
934: int AttachNullSpace(PC pc,Vec model)
935: {
936:   int          i,ierr,rstart,rend,N;
937:   MatNullSpace sp;
938:   Vec          v,vs[1];
939:   PetscScalar  *vx,scale;

942:   VecDuplicate(model,&v);
943:   VecGetSize(model,&N);
944:   scale = 2.0/sqrt(N);
945:   VecGetArray(v,&vx);
946:   VecGetOwnershipRange(v,&rstart,&rend);
947:   for (i=rstart; i<rend; i++) {
948:     if (!(i % 4)) vx[i-rstart] = scale;
949:     else          vx[i-rstart] = 0.0;
950:   }
951:   VecRestoreArray(v,&vx);
952:   vs[0] = v;
953:   MatNullSpaceCreate(PETSC_COMM_WORLD,0,1,vs,&sp);
954:   VecDestroy(v);
955:   PCNullSpaceAttach(pc,sp);
956:   MatNullSpaceDestroy(sp);
957:   return(0);
958: }

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

974:   PetscScalar   xx;
975:   PetscScalar   vx,vy,avx,avy,vxp,vxm,vyp,vym;
976:   PetscScalar   Bx,By,aBx,aBy,Bxp,Bxm,Byp,Bym;
977:   PassiveField  **xold;

980:   de2     = sqr(user->param->d_e);
981:   rhos2   = sqr(user->param->rho_s);
982:   nu      = user->param->nu;
983:   eta     = user->param->eta;
984:   dde2    = one/de2;

986:   /* 
987:      Define mesh intervals ratios for uniform grid.
988:      [Note: FD formulae below are normalized by multiplying through by
989:      local volume element to obtain coefficients O(1) in two dimensions.]
990:   */
991:   dhx   = info->mx/lx;        dhy   = info->my/ly;
992:   hx    = one/dhx;             hy   = one/dhy;
993:   hxdhy = hx*dhy;           hydhx   = hy*dhx;
994:   hxhy  = hx*hy;             dhxdhy = dhx*dhy;

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


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

1002: #ifdef EQ
1003:       xx = i * hx;
1004:       F_eq_x = - (1. + de2) * sin(xx);
1005:       By_eq = sin(xx);
1006: #else
1007:       F_eq_x = 0.;
1008:       By_eq = 0.;
1009: #endif

1011:       /*
1012:        * convective coefficients for upwinding
1013:        */

1015:       vx  = - D_y(x,phi,i,j);
1016:       vy  =   D_x(x,phi,i,j);
1017:       avx = PetscAbsScalar(vx); vxp = p5*(vx+avx); vxm = p5*(vx-avx);
1018:       avy = PetscAbsScalar(vy); vyp = p5*(vy+avy); vym = p5*(vy-avy);
1019: #ifndef UPWINDING
1020:       vxp = vxm = p5*vx;
1021:       vyp = vym = p5*vy;
1022: #endif

1024:       Bx  =   D_y(x,psi,i,j);
1025:       By  = - D_x(x,psi,i,j) + By_eq;
1026:       aBx = PetscAbsScalar(Bx); Bxp = p5*(Bx+aBx); Bxm = p5*(Bx-aBx);
1027:       aBy = PetscAbsScalar(By); Byp = p5*(By+aBy); Bym = p5*(By-aBy);
1028: #ifndef UPWINDING
1029:       Bxp = Bxm = p5*Bx;
1030:       Byp = Bym = p5*By;
1031: #endif

1033:       DAVecGetArray(info->da,user->Xold,&xold);
1034:       dtinv = hxhy/(tsCtx->dt);
1035:       switch(c) {

1037:         case 0:
1038:           /* Lap(phi) - U */
1039:           *f = (Lapl(x,phi,i,j) - x[j][i].U) * hxhy;
1040:           break;

1042:         case 1:
1043:           /* psi - d_e^2 * Lap(psi) - F */
1044:           *f = (x[j][i].psi - de2 * Lapl(x,psi,i,j) - x[j][i].F) * hxhy;
1045:           break;

1047:         case 2:
1048:           /* vx * U_x + vy * U_y - (B_x * F_x + B_y * F_y) / d_e^2
1049:             - nu Lap(U) */
1050:           *f  =
1051:             ((vxp * (D_xm(x,U,i,j)) + vxm * (D_xp(x,U,i,j)) +
1052:               vyp * (D_ym(x,U,i,j)) + vym * (D_yp(x,U,i,j))) -
1053:              (Bxp * (D_xm(x,F,i,j) + F_eq_x) + Bxm * (D_xp(x,F,i,j) + F_eq_x) +
1054:               Byp * (D_ym(x,F,i,j)) + Bym * (D_yp(x,F,i,j))) * dde2 -
1055:              nu * Lapl(x,U,i,j)) * hxhy;
1056:           *f += dtinv*(x[j][i].U-xold[j][i].U);
1057:           break;

1059:         case 3:
1060:           /* vx * F_x + vy * F_y - rho_s^2 * (B_x * U_x + B_y * U_y)
1061:            - eta * Lap(psi) */
1062:           *f  =
1063:             ((vxp * (D_xm(x,F,i,j) + F_eq_x) + vxm * (D_xp(x,F,i,j) + F_eq_x) +
1064:              vyp * (D_ym(x,F,i,j)) + vym * (D_yp(x,F,i,j))) -
1065:             (Bxp * (D_xm(x,U,i,j)) + Bxm * (D_xp(x,U,i,j)) +
1066:              Byp * (D_ym(x,U,i,j)) + Bym * (D_yp(x,U,i,j))) * rhos2 -
1067:             eta * Lapl(x,psi,i,j)) * hxhy;
1068:           *f += dtinv*(x[j][i].F-xold[j][i].F);
1069:           break;
1070:       }
1071:       DAVecRestoreArray(info->da,user->Xold,&xold);


1074:   /*
1075:      Flop count (multiply-adds are counted as 2 operations)
1076:   */
1077:   /*  PetscLogFlops(84*info->ym*info->xm); FIXME */
1078:   return(0);
1079: }