Actual source code: ex4.c

  1: static char help[] = "\n";
  2: /* - - - - - - - - - - - - - - - - - - - - - - - - 
  3:    ex4.c
  4:    Advection-diffusion equation

  6:    Du
  7:    -- - \nabla^2 u = 0
  8:    dt 

 10: which we discretize using Crank-Nicholson

 12:    u^+ - u^*  1 /                                 \
 13:    -------- - - | (\nabla^2 u)^+ + (\nabla^2 u)^* | = 0
 14:    \Delta t   2 \                                 /

 16: which gives, collecting terms

 18:          \Delta t                        \Delta t
 19:    u^+ - -------- (\nabla^2 u)^+ = u^* + -------- (\nabla^2 u)^*
 20:              2                              2

 22: where u^- is the solution at the old time, u^+ is
 23: the solution at the new time, and u^* is the solution
 24: at the old time advected to the new time.
 25: - - - - - - - - - - - - - - - - - - - - - - - - */
 26: #include <petscsnes.h>
 27: #include <petscdmda.h>
 28: #include <petscdmmg.h>
 29: #include <petscbag.h>
 30: #include <petsccharacteristic.h>

 32: #define EXAMPLE_NUMBER 1
 33: #define SHEAR_CELL     0
 34: #define SOLID_BODY     1
 35: #define FNAME_LENGTH   60

 37: typedef struct field_s {
 38:   PetscReal phi;
 39: } Field;

 41: typedef struct parameter_s {
 42:   int            ni, nj,pi,pj;          /* number of grid points, number of procs */
 43:   PassiveScalar  amp,sigma,xctr,zctr,L1,L2,LINF; /* parameters for gaussian Initial Condition */
 44:   PassiveScalar  Pe, theta, ct, st, diffScale;   /* parameters for velocity field and diffusion length */
 45:   int            flow_type, sl_event;
 46:   PetscBool      param_test, output_to_file;
 47:   char           output_filename[FNAME_LENGTH];
 48:   /* timestep stuff */
 49:   PassiveScalar  t; /* the time */
 50:   int            n; /* the time step */
 51:   PassiveScalar  dtAdvection, dtDiffusion; /* minimal advection and diffusion time steps */
 52:   PassiveScalar  t_max, dt, cfl, t_output_interval;
 53:   int            N_steps;
 54: } Parameter;

 56: typedef struct gridinfo_s {
 57:   DMDABoundaryType periodic;
 58:   DMDAStencilType  stencil;
 59:   int            ni,nj,dof,stencil_width,mglevels;
 60:   PassiveScalar  dx,dz;
 61: } GridInfo;

 63: typedef struct appctx_s {
 64:   DMMG        *dmmg;
 65:   Vec          Xold;
 66:   PetscBag     bag;
 67:   GridInfo     *grid;
 68: } AppCtx;

 70: /* Main routines */
 71: int SetParams            (AppCtx*);
 72: int ReportParams         (AppCtx*);
 73: int Initialize           (DMMG*);
 74: int DoSolve              (DMMG*);
 75: int DoOutput             (DMMG*, int);
 76: PetscReal BiCubicInterp  (Field**, PetscReal, PetscReal);
 77: PetscReal CubicInterp    (PetscReal, PetscReal, PetscReal, PetscReal, PetscReal);
 78: PetscBool  OptionsHasName(const char*);

 80: /* characteristic call-backs (static interface) */
 81: PetscErrorCode InterpVelocity2D(void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);
 82: PetscErrorCode InterpFields2D  (void*, PetscReal[], PetscInt, PetscInt[], PetscReal[], void*);

 84: PetscErrorCode FormOldTimeFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);
 85: PetscErrorCode FormNewTimeFunctionLocal(DMDALocalInfo *, PetscScalar **, PetscScalar **, AppCtx *);

 87: /* a few macros for convenience */
 88: #define REG_REAL(A,B,C,D,E)   ierr=PetscBagRegisterReal(A,B,C,D,E);CHKERRQ(ierr)
 89: #define REG_INTG(A,B,C,D,E)   ierr=PetscBagRegisterInt(A,B,C,D,E);CHKERRQ(ierr)
 90: #define REG_TRUE(A,B,C,D,E)   ierr=PetscBagRegisterBool(A,B,C,D,E);CHKERRQ(ierr)
 91: #define REG_STRG(A,B,C,D,E,F) ierr=PetscBagRegisterString(A,B,C,D,E,F);CHKERRQ(ierr)
 92: #define REG_ENUM(A,B,C,D,E,F) ierr=PetscBagRegisterEnum(A,B,C,D,E,F);CHKERRQ(ierr)

 94: /*-----------------------------------------------------------------------*/
 97: int main(int argc,char **argv)
 98: /*-----------------------------------------------------------------------*/
 99: {
100:   AppCtx         *user;               /* user-defined work context */
101:   Parameter      *param;
102:   DM              da;
103:   GridInfo        grid;
104:   MPI_Comm        comm;
105:   PetscErrorCode  ierr;

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

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Set up the problem parameters.
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   PetscMalloc(sizeof(AppCtx),&user);
114:   PetscBagCreate(comm,sizeof(Parameter),&(user->bag));
115:   user->grid = &grid;
116:   SetParams(user);
117:   ReportParams(user);
118:   PetscBagGetData(user->bag,(void**)&param);

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
122:      for principal unknowns (x) and governing residuals (f)
123:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124:   DMMGCreate(comm,grid.mglevels,user,&user->dmmg);
125:   DMDACreate2d(comm,grid.periodic,grid.periodic,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
126:   DMMGSetDM(user->dmmg,(DM) da);
127:   DMDestroy(&da);
128:   DMMGSetSNESLocal(user->dmmg,FormNewTimeFunctionLocal,PETSC_NULL,PETSC_NULL,PETSC_NULL);
129:   DMMGSetFromOptions(user->dmmg);
130:   DMDAGetInfo(DMMGGetDM(user->dmmg),PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,&(param->pi),&(param->pj),PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
131:   REG_INTG(user->bag,&param->pi,param->pi ,"procs_x","<DO NOT SET> Processors in the x-direction");
132:   REG_INTG(user->bag,&param->pj,param->pj ,"procs_y","<DO NOT SET> Processors in the y-direction");

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:      Create user context, set problem data, create vector data structures.
136:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137:   DMGetGlobalVector(DMMGGetDM(user->dmmg), &(user->Xold));

139:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140:      Initialize and solve the nonlinear system
141:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142:   Initialize(user->dmmg);
143:   DoSolve(user->dmmg);
144: 
145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:      Free work space. 
147:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148:   DMRestoreGlobalVector(DMMGGetDM(user->dmmg), &(user->Xold));
149:   PetscBagDestroy(&user->bag);
150:   DMMGDestroy(user->dmmg);
151:   PetscFree(user);
152:   PetscFinalize();
153:   return 0;
154: }

156: /*---------------------------------------------------------------------*/
159: int SetParams(AppCtx *user)
160: /*---------------------------------------------------------------------*/
161: {
162:   PetscBag   bag = user->bag;
163:   Parameter  *p;
164:   GridInfo   *grid = user->grid;
165:   int        ierr, ierr_out=0;
166:   PetscReal  PI = 3.14159265358979323846;

168:   PetscBagGetData(bag,(void**)&p);

170:   /* domain geometry & grid size */
171:   REG_INTG(bag,&p->ni,40                  ,"ni","Grid points in x-dir");
172:   REG_INTG(bag,&p->nj,40                  ,"nj","Grid points in y-dir");
173:   grid->ni = p->ni;
174:   grid->nj = p->nj;
175:   grid->dx = 1.0/((double)(grid->ni - 1));
176:   grid->dz = 1.0/((double)(grid->nj - 1));

178:   /* initial conditions */
179:   REG_INTG(bag,&p->flow_type,SOLID_BODY   ,"flow_type","Flow field mode: 0=shear cell, 1=translation");
180:   REG_REAL(bag,&p->amp,1                  ,"amp","Amplitude of the gaussian IC");
181:   REG_REAL(bag,&p->sigma,0.07             ,"sigma","Standard deviation of the gaussian IC");
182:   REG_REAL(bag,&p->xctr,0.5               ,"xctr","x-position of the center of the gaussian IC");
183:   REG_REAL(bag,&p->zctr,0.5               ,"zctr","z-position of the center of the gaussian IC");

185:   /* Velocity field */
186:   REG_REAL(bag,&p->Pe,10                   ,"Pe","Peclet number");
187:   REG_REAL(bag,&p->theta,-90                ,"theta","velocity angle (degrees)");
188:   REG_REAL(bag,&p->ct,cos(p->theta*PI/180),"cosTheta","<DO NOT SET> cosine velocity angle");
189:   REG_REAL(bag,&p->st,sin(p->theta*PI/180),"sinTheta","<DO NOT SET> sine velocity angle");
190: 
191:   /* diffusion LengthScale for time stepping */
192:   REG_REAL(bag,&p->diffScale,2,            "diffScale","diffusion length scale (number of grid points for stable diffusion)");

194:   /* time stepping */
195:   REG_REAL(bag,&p->t_max,1                ,"t_max","Maximum dimensionless time");
196:   REG_REAL(bag,&p->cfl,5                  ,"cfl","Courant number");
197:   REG_REAL(bag,&p->t_output_interval,0.1  ,"t_output","Dimensionless time interval for output");
198:   REG_INTG(bag,&p->N_steps,1000           ,"nsteps","Maximum time-steps");
199:   REG_INTG(bag,&p->n,1                    ,"nsteps","<DO NOT SET> current time-step");
200:   REG_REAL(bag,&p->t,0.0                  ,"time","<DO NOT SET> initial time");
201:   REG_REAL(bag,&p->dtAdvection,1./p->Pe/(sqrt(pow(p->ct/grid->dx,2)+pow(p->st/grid->dz,2))),"dtAdvection","<DO NOT SET> CFL limited advection time step");
202:   REG_REAL(bag,&p->dtDiffusion,1./(1./grid->dx/grid->dx+1./grid->dz/grid->dz),"dtDiffusion","<DO NOT SET> grid-space limited diffusion time step");
203:   REG_REAL(bag,&p->dt,PetscMin(p->cfl*p->dtAdvection,pow(p->diffScale,2)*p->dtDiffusion),"dt","<DO NOT SET> time-step size");

205:   /* output options */
206:   REG_TRUE(bag,&p->param_test     ,PETSC_FALSE  ,"test","Run parameter test only (T/F)");
207:   REG_STRG(bag,&p->output_filename,FNAME_LENGTH ,"null","output_file","Name base for output files, set with: -output_file <filename>");
208:   REG_TRUE(bag,&p->output_to_file,PETSC_FALSE   ,"do_output","<DO NOT SET> flag will be true if you specify an output file name");
209:   p->output_to_file = OptionsHasName("-output_file");

211:   grid->ni            = p->ni;
212:   grid->nj            = p->nj;
213:   grid->periodic      = DMDA_BOUNDARY_PERIODIC;
214:   grid->stencil       = DMDA_STENCIL_BOX;
215:   grid->dof           = 1;
216:   grid->stencil_width = 2;
217:   grid->mglevels      = 1;
218:   return ierr_out;
219: }

221: /*---------------------------------------------------------------------*/
224: int ReportParams(AppCtx *user)
225: /*---------------------------------------------------------------------*/
226: {
227:   Parameter  *param;
228:   GridInfo   *grid = user->grid;
229:   int        ierr, ierr_out=0;
230:   PetscBagGetData(user->bag,(void**)&param);

232:   PetscPrintf(PETSC_COMM_WORLD,"---------------MOC test 1----------------\n");
233:   PetscPrintf(PETSC_COMM_WORLD,"Prescribed wind, method of\n");
234:   PetscPrintf(PETSC_COMM_WORLD,"characteristics advection, explicit time-\n");
235:   PetscPrintf(PETSC_COMM_WORLD,"stepping.\n\n");
236:   if (param->flow_type == 0) {
237:     PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (shear cell).\n\n",param->flow_type);
238:   }
239:   if (param->flow_type == 1) {
240:     PetscPrintf(PETSC_COMM_WORLD,"Flow_type: %d (translation).\n\n",param->flow_type);
241:   }
242:   PetscPrintf(PETSC_COMM_WORLD,"  [ni,nj] = %d, %d   [dx,dz] = %5.4g, %5.4g\n",grid->ni,grid->nj,grid->dx,grid->dz);
243:   PetscPrintf(PETSC_COMM_WORLD,"  t_max = %g, cfl = %g, dt = %5.4g,",param->t_max,param->cfl,param->dt);
244:   PetscPrintf(PETSC_COMM_WORLD," t_output = %g\n",param->t_output_interval);
245:   PetscPrintf(PETSC_COMM_WORLD," dt_advection= %g, dt_diffusion= %g\n",param->dtAdvection,param->dtDiffusion);
246:   if (param->output_to_file) {
247:     PetscPrintf(PETSC_COMM_WORLD,"Output File:       Binary file \"%s\"\n",param->output_filename);
248:   }
249:   if (!param->output_to_file)
250:     PetscPrintf(PETSC_COMM_WORLD,"Output File:       NO OUTPUT!\n");
251:   PetscPrintf(PETSC_COMM_WORLD,"----------------------------------------\n");
252:   if (param->param_test) PetscEnd();
253:   return ierr_out;
254: }

256: /* ------------------------------------------------------------------- */
259: int Initialize(DMMG *dmmg)
260: /* ------------------------------------------------------------------- */
261: {
262:   AppCtx    *user  = (AppCtx*)dmmg[0]->user;
263:   Parameter *param;
264:   DM        da;
265:   PetscReal amp, sigma, xc, zc ;
266:   PetscReal dx=user->grid->dx,dz=user->grid->dz;
267:   int       i,j,ierr,is,js,im,jm;
268:   Field     **x;
269:   PetscBagGetData(user->bag,(void**)&param);
270: 
271:   amp = param->amp;
272:   sigma = param->sigma;
273:   xc = param->xctr; zc = param->zctr;

275:   /* Get the DM and grid */
276:   da = DMMGGetDM(dmmg);
277:   DMDAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
278:   DMDAVecGetArray(da,user->Xold,(void**)&x);

280:   for (j=js; j<js+jm; j++) {
281:     for (i=is; i<is+im; i++) {
282:       x[j][i].phi = param->amp*exp(-0.5*((i*dx-xc)*(i*dx-xc)+(j*dz-zc)*(j*dz-zc))/sigma/sigma);
283:     }
284:   }
285: 
286:   /* restore the grid to it's vector */
287:   DMDAVecRestoreArray(da,user->Xold,(void**)&x);
288:   VecCopy(user->Xold, DMMGGetx(dmmg));

290:   return 0;
291: }

293: /* ------------------------------------------------------------------- */
296: int DoSolve(DMMG *dmmg)
297: /* ------------------------------------------------------------------- */
298: {
299:   AppCtx         *user  = (AppCtx*)dmmg[0]->user;
300:   Parameter      *param;
301:   PetscReal      t_output = 0.0;
302:   int            ierr, n_plot = 0, Ncomponents, components[3];
303:   DM             da = DMMGGetDM(dmmg);
304:   Vec            Xstar;
305:   Characteristic c;
306:   PetscBagGetData(user->bag,(void**)&param);

308:   DMGetGlobalVector(da, &Xstar);

310:   /*------------ BEGIN CHARACTERISTIC SETUP ---------------*/
311:   CharacteristicCreate(PETSC_COMM_WORLD, &c);
312:   /* set up the velocity interpolation system */
313:   Ncomponents = 2; components[0] = 0; components[1] = 0;
314:   CharacteristicSetVelocityInterpolationLocal(c, da, DMMGGetx(dmmg), user->Xold, Ncomponents, components, InterpVelocity2D, user);
315:   /* set up the fields interpolation system */
316:   Ncomponents = 1; components[0] = 0;
317:   CharacteristicSetFieldInterpolationLocal(c, da, user->Xold, Ncomponents, components, InterpFields2D, user);
318:   /*------------ END CHARACTERISTIC SETUP ----------------*/

320:   /* output initial data */
321:   PetscPrintf(PETSC_COMM_WORLD," Initialization, Time: %5.4g\n", param->t);
322:   DoOutput(dmmg,n_plot);
323:   t_output += param->t_output_interval; n_plot++;

325:   /* timestep loop */
326:   for (param->t=param->dt; param->t<=param->t_max; param->t+=param->dt) {
327:     if (param->n > param->N_steps) {
328:       PetscPrintf(PETSC_COMM_WORLD,"EXCEEDED MAX NUMBER OF TIMESTEPS! EXITING SOLVE!\n");
329:       return 0;
330:     }

332:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333:        Solve at time t & copy solution into solution vector.
334:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
335:     /* Evaluate operator (I + \Delta t/2 L) u^-  = X^- */
336:     DMDAFormFunctionLocal(da, (DMDALocalFunction1) FormOldTimeFunctionLocal, DMMGGetx(dmmg), user->Xold, user);
337:     /* Advect Xold into Xstar */
338:     CharacteristicSolve(c, param->dt, Xstar);
339:     /* Xstar -> Xold */
340:     VecCopy(Xstar, user->Xold);
341:     /* Solve u^+ = (I - \Delta t/2 L)^-1 Xstar which could be F(u^+) = Xstar */
342:     DMMGSolve(dmmg);

344:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
345:        report step and update counter.
346:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
347:     PetscPrintf(PETSC_COMM_WORLD," Step: %d, Time: %5.4g\n", param->n, param->t);
348:     param->n++;

350:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351:        Output variables.
352:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
353:     if (param->t >= t_output) {
354:       DoOutput(dmmg,n_plot);
355:       t_output += param->t_output_interval; n_plot++;
356:     }
357:   }
358:   DMRestoreGlobalVector(da, &Xstar);
359:   CharacteristicDestroy(&c);
360:   return 0;
361: }

363: /*---------------------------------------------------------------------*/
366: /* linear interpolation, ir: [0, ni], jr: [0, nj] */
367: PetscErrorCode InterpVelocity2D(void *f, PetscReal ij_real[], PetscInt numComp, 
368:                                 PetscInt components[], PetscReal velocity[], 
369:                                 void *ctx)
370: /*---------------------------------------------------------------------*/
371: {
372:   AppCtx    *user = (AppCtx *) ctx;
373:   Parameter *param;
374:   PetscReal dx=user->grid->dx, dz=user->grid->dz;
375:   PetscReal PI = 3.14159265358979323846;
376:   int       ierr;
377:   PetscBagGetData(user->bag,(void**)&param);

379:   /* remember: must be coordinate velocities not true velocities */
380:   if (param->flow_type == SHEAR_CELL) {
381:     velocity[0] = -sin(PI*ij_real[0]*dx)*cos(PI*ij_real[1]*dz)/dx;
382:     velocity[1] =  sin(PI*ij_real[1]*dz)*cos(PI*ij_real[0]*dx)/dz;
383:   } else {
384:     velocity[0] = param->Pe*param->ct/dx;
385:     velocity[1] = param->Pe*param->st/dz;
386:   }
387:   return 0;
388: }

390: /*---------------------------------------------------------------------*/
393: PetscErrorCode InterpFields2D(void *f, PetscReal ij_real[], PetscInt numComp, 
394:                               PetscInt components[], PetscReal field[], 
395:                               void *ctx)
396: /*---------------------------------------------------------------------*/
397: {
398:   AppCtx    *user = (AppCtx*)ctx;
399:   Field     **x   = (Field**)f;
400:   int       ni=user->grid->ni, nj=user->grid->nj;
401:   int       ierr;
402:   PetscReal ir=ij_real[0], jr=ij_real[1];

404:   /* map back to periodic domain if out of bounds */
405:   if ( ir < 0 || ir > ni-1 || jr < 0 || jr> nj-1 ) {
406:     DMDAMapCoordsToPeriodicDomain(DMMGGetDM(user->dmmg), &ir, &jr);
407:   }
408:   field[0] = BiCubicInterp(x, ir, jr);
409:   return 0;
410: }

412: /*---------------------------------------------------------------------*/
415: PetscReal BiCubicInterp(Field **x, PetscReal ir, PetscReal jr)
416: /*---------------------------------------------------------------------*/
417: {
418:   int        im, jm, imm,jmm,ip,jp,ipp,jpp;
419:   PetscReal  il, jl, row1, row2, row3, row4;
420:   im = (int)floor(ir); jm = (int)floor(jr);
421:   il = ir - im + 1.0; jl = jr - jm + 1.0;
422:   imm = im-1; ip = im+1; ipp = im+2;
423:   jmm = jm-1; jp = jm+1; jpp = jm+2;
424:   row1 = CubicInterp(il,x[jmm][imm].phi,x[jmm][im].phi,x[jmm][ip].phi,x[jmm][ipp].phi);
425:   row2 = CubicInterp(il,x[jm] [imm].phi,x[jm] [im].phi,x[jm] [ip].phi,x[jm] [ipp].phi);
426:   row3 = CubicInterp(il,x[jp] [imm].phi,x[jp] [im].phi,x[jp] [ip].phi,x[jp] [ipp].phi);
427:   row4 = CubicInterp(il,x[jpp][imm].phi,x[jpp][im].phi,x[jpp][ip].phi,x[jpp][ipp].phi);
428:   return CubicInterp(jl,row1,row2,row3,row4);
429: }

431: /*---------------------------------------------------------------------*/
434: PetscReal CubicInterp(PetscReal x, PetscReal y_1, PetscReal y_2, 
435:                       PetscReal y_3, PetscReal y_4)
436: /*---------------------------------------------------------------------*/
437: {
438:   PetscReal  sxth=0.16666666666667, retval;
439:   retval = - y_1*(x-1.0)*(x-2.0)*(x-3.0)*sxth + y_2*(x-0.0)*(x-2.0)*(x-3.0)*0.5
440:            - y_3*(x-0.0)*(x-1.0)*(x-3.0)*0.5  + y_4*(x-0.0)*(x-1.0)*(x-2.0)*sxth;
441:   return retval;
442: }

444: /*---------------------------------------------------------------------*/
447: int DoOutput(DMMG *dmmg, int n_plot)
448: /*---------------------------------------------------------------------*/
449: {
450:   AppCtx      *user = (AppCtx*)dmmg[0]->user;
451:   Parameter   *param;
452:   int         ierr;
453:   char        filename[FNAME_LENGTH];
454:   PetscViewer viewer;
455:   DM          da;
456:   PetscBagGetData(user->bag,(void**)&param);
457:   da = DMMGGetDM(dmmg);

459:   if (param->output_to_file) { /* send output to binary file */
460:     /* generate filename for time t */
461:     sprintf(filename,"%s_%3.3d",param->output_filename,n_plot);
462:     PetscPrintf(PETSC_COMM_WORLD,"Generating output: time t = %g, ",param->t);
463:     PetscPrintf(PETSC_COMM_WORLD,"file = \"%s\"\n",filename);

465:     /* make output files */
466:     PetscViewerBinaryMatlabOpen(PETSC_COMM_WORLD,filename,&viewer);
467:     PetscViewerBinaryMatlabOutputBag(viewer,"par",user->bag);
468:     DMDASetFieldName(da,0,"phi");
469:     PetscViewerBinaryMatlabOutputVecDA(viewer,"field",DMMGGetx(dmmg),da);
470:     PetscViewerBinaryMatlabDestroy(&viewer);
471:   }
472:   return 0;
473: }

475: /* ------------------------------------------------------------------- */
478: PetscBool  OptionsHasName(const char name[])
479: /* ------------------------------------------------------------------- */
480: {
481:   PetscBool  retval;
482:   int ierr;
483:   PetscOptionsHasName(PETSC_NULL,name,&retval);
484:   return retval;
485: }

489: /* 
490:   FormNewTimeFunctionLocal - Evaluates f = (I - \Delta t/2 L) x - f(u^-)^*.

492:   Note: We get f(u^-)^* from Xold in the user context.

494:   Process adiC(36): FormNewTimeFunctionLocal
495: */
496: PetscErrorCode FormNewTimeFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
497: {
498:   DM             da = DMMGGetDM(user->dmmg);
499:   PetscScalar  **fold;
500:   Parameter     *param;
501:   PetscScalar    u,uxx,uyy;
502:   PetscReal      hx,hy,hxdhy,hydhx;
503:   PetscInt       i,j;

507:   hx     = 1.0/(PetscReal)(info->mx-1);
508:   hy     = 1.0/(PetscReal)(info->my-1);
509:   hxdhy  = hx/hy;
510:   hydhx  = hy/hx;

512:   PetscBagGetData(user->bag, (void**) &param);
513:   DMDAVecGetArray(da, user->Xold, &fold);
514:   for (j=info->ys; j<info->ys+info->ym; j++) {
515:     for (i=info->xs; i<info->xs+info->xm; i++) {
516:       u       = x[j][i];
517:       uxx     = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
518:       uyy     = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
519:       f[j][i] = u*hx*hy + param->dt*0.5*(uxx + uyy) - fold[j][i];
520:     }
521:   }
522:   DMDAVecRestoreArray(da, user->Xold, &fold);

524:   PetscLogFlops(13.0*info->ym*info->xm);
525:   return(0);
526: }

530: /* 
531:   FormOldTimeFunctionLocal - Evaluates f = (I + \Delta t/2 L) x.

533:   Process adiC(36): FormOldTimeFunctionLocal
534: */
535: PetscErrorCode FormOldTimeFunctionLocal(DMDALocalInfo *info, PetscScalar **x, PetscScalar **f, AppCtx *user)
536: {
537:   Parameter     *param;
538:   PetscScalar    u,uxx,uyy;
539:   PetscReal      hx,hy,hxdhy,hydhx;
540:   PetscInt       i,j;

544:   hx     = 1.0/(PetscReal)(info->mx-1);
545:   hy     = 1.0/(PetscReal)(info->my-1);
546:   hxdhy  = hx/hy;
547:   hydhx  = hy/hx;

549:   PetscBagGetData(user->bag, (void**) &param);
550:   for (j=info->ys; j<info->ys+info->ym; j++) {
551:     for (i=info->xs; i<info->xs+info->xm; i++) {
552:       u       = x[j][i];
553:       uxx     = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
554:       uyy     = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
555:       f[j][i] = u*hx*hy - param->dt*0.5*(uxx + uyy);
556:     }
557:   }

559:   PetscLogFlops(12.0*info->ym*info->xm);
560:   return(0);
561: }