Actual source code: ex2.c

  1: /*$Id: ex2.c,v 1.41 2001/08/10 03:34:17 bsmith Exp $*/
  2: static char help[] ="Solves a time-dependent nonlinear PDE. Uses implicitn
  3: timestepping.  Runtime options include:n
  4:   -M <xg>, where <xg> = number of grid pointsn
  5:   -debug : Activate debugging printoutsn
  6:   -nox   : Deactivate x-window graphicsnn";

  8: /*
  9:    Concepts: TS^time-dependent nonlinear problems
 10:    Processors: n
 11: */

 13: /* ------------------------------------------------------------------------

 15:    This program solves the PDE

 17:                u * u_xx 
 18:          u_t = ---------
 19:                2*(t+1)^2 

 21:     on the domain 0 <= x <= 1, with boundary conditions
 22:          u(t,0) = t + 1,  u(t,1) = 2*t + 2,
 23:     and initial condition
 24:          u(0,x) = 1 + x*x.

 26:     The exact solution is:
 27:          u(t,x) = (1 + x*x) * (1 + t)

 29:     Note that since the solution is linear in time and quadratic in x,
 30:     the finite difference scheme actually computes the "exact" solution.

 32:     We use by default the backward Euler method.

 34:   ------------------------------------------------------------------------- */

 36: /*
 37:    Include "petscts.h" to use the PETSc timestepping routines. Note that
 38:    this file automatically includes "petsc.h" and other lower-level
 39:    PETSc include files.

 41:    Include the "petscda.h" to allow us to use the distributed array data 
 42:    structures to manage the parallel grid.
 43: */
 44:  #include petscts.h
 45:  #include petscda.h

 47: /* 
 48:    User-defined application context - contains data needed by the 
 49:    application-provided callback routines.
 50: */
 51: typedef struct {
 52:   MPI_Comm   comm;          /* communicator */
 53:   DA         da;            /* distributed array data structure */
 54:   Vec        localwork;     /* local ghosted work vector */
 55:   Vec        u_local;       /* local ghosted approximate solution vector */
 56:   Vec        solution;      /* global exact solution vector */
 57:   int        m;             /* total number of grid points */
 58:   PetscReal  h;             /* mesh width: h = 1/(m-1) */
 59:   PetscTruth debug;         /* flag (1 indicates activation of debugging printouts) */
 60: } AppCtx;

 62: /* 
 63:    User-defined routines, provided below.
 64: */
 65: extern int InitialConditions(Vec,AppCtx*);
 66: extern int RHSFunction(TS,PetscReal,Vec,Vec,void*);
 67: extern int RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
 68: extern int Monitor(TS,int,PetscReal,Vec,void*);
 69: extern int ExactSolution(PetscReal,Vec,AppCtx*);

 71: /*
 72:    Utility routine for finite difference Jacobian approximation
 73: */
 74: extern int RHSJacobianFD(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);

 76: #undef __FUNCT__
 78: int main(int argc,char **argv)
 79: {
 80:   AppCtx     appctx;                 /* user-defined application context */
 81:   TS         ts;                     /* timestepping context */
 82:   Mat        A;                      /* Jacobian matrix data structure */
 83:   Vec        u;                      /* approximate solution vector */
 84:   int        time_steps_max = 1000;  /* default max timesteps */
 85:   int        ierr,steps;
 86:   PetscReal  ftime;                  /* final time */
 87:   PetscReal  dt;
 88:   PetscReal  time_total_max = 100.0; /* default max total time */
 89:   PetscTruth flg;

 91:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:      Initialize program and set problem parameters
 93:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 94: 
 95:   PetscInitialize(&argc,&argv,(char*)0,help);

 97:   appctx.comm = PETSC_COMM_WORLD;
 98:   appctx.m    = 60;
 99:   PetscOptionsGetInt(PETSC_NULL,"-M",&appctx.m,PETSC_NULL);
100:   PetscOptionsHasName(PETSC_NULL,"-debug",&appctx.debug);
101:   appctx.h    = 1.0/(appctx.m-1.0);

103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104:      Create vector data structures
105:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

107:   /*
108:      Create distributed array (DA) to manage parallel grid and vectors
109:      and to set up the ghost point communication pattern.  There are M 
110:      total grid values spread equally among all the processors.
111:   */
112:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,appctx.m,1,1,PETSC_NULL,
113:                     &appctx.da);

115:   /*
116:      Extract global and local vectors from DA; we use these to store the
117:      approximate solution.  Then duplicate these for remaining vectors that
118:      have the same types.
119:   */
120:   DACreateGlobalVector(appctx.da,&u);
121:   DACreateLocalVector(appctx.da,&appctx.u_local);

123:   /*
124:      Create local work vector for use in evaluating right-hand-side function;
125:      create global work vector for storing exact solution.
126:   */
127:   VecDuplicate(appctx.u_local,&appctx.localwork);
128:   VecDuplicate(u,&appctx.solution);

130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131:      Create timestepping solver context; set callback routine for
132:      right-hand-side function evaluation.
133:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

135:   TSCreate(PETSC_COMM_WORLD,&ts);
136:   TSSetProblemType(ts,TS_NONLINEAR);
137:   TSSetRHSFunction(ts,RHSFunction,&appctx);

139:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140:      Set optional user-defined monitoring routine
141:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

143:   TSSetMonitor(ts,Monitor,&appctx,PETSC_NULL);

145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:      For nonlinear problems, the user can provide a Jacobian evaluation
147:      routine (or use a finite differencing approximation).

149:      Create matrix data structure; set Jacobian evaluation routine.
150:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

152:   MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,appctx.m,appctx.m,&A);
153:   MatSetFromOptions(A);
154:   PetscOptionsHasName(PETSC_NULL,"-fdjac",&flg);
155:   if (flg) {
156:     TSSetRHSJacobian(ts,A,A,RHSJacobianFD,&appctx);
157:   } else {
158:     TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx);
159:   }

161:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162:      Set solution vector and initial timestep
163:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

165:   dt   = appctx.h/2.0;
166:   TSSetInitialTimeStep(ts,0.0,dt);
167:   TSSetSolution(ts,u);

169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:      Customize timestepping solver:  
171:        - Set the solution method to be the Backward Euler method.
172:        - Set timestepping duration info 
173:      Then set runtime options, which can override these defaults.
174:      For example,
175:           -ts_max_steps <maxsteps> -ts_max_time <maxtime>
176:      to override the defaults set by TSSetDuration().
177:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

179:   TSSetType(ts,TS_BEULER);
180:   TSSetDuration(ts,time_steps_max,time_total_max);
181:   TSSetFromOptions(ts);

183:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184:      Solve the problem
185:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

187:   /*
188:      Evaluate initial conditions
189:   */
190:   InitialConditions(u,&appctx);

192:   /*
193:      Run the timestepping solver
194:   */
195:   TSStep(ts,&steps,&ftime);

197:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198:      Free work space.  All PETSc objects should be destroyed when they
199:      are no longer needed.
200:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

202:   TSDestroy(ts);
203:   VecDestroy(u);
204:   MatDestroy(A);
205:   DADestroy(appctx.da);
206:   VecDestroy(appctx.localwork);
207:   VecDestroy(appctx.solution);
208:   VecDestroy(appctx.u_local);

210:   /*
211:      Always call PetscFinalize() before exiting a program.  This routine
212:        - finalizes the PETSc libraries as well as MPI
213:        - provides summary and diagnostic information if certain runtime
214:          options are chosen (e.g., -log_summary). 
215:   */
216:   PetscFinalize();
217:   return 0;
218: }
219: /* --------------------------------------------------------------------- */
220: #undef __FUNCT__
222: /*
223:    InitialConditions - Computes the solution at the initial time. 

225:    Input Parameters:
226:    u - uninitialized solution vector (global)
227:    appctx - user-defined application context

229:    Output Parameter:
230:    u - vector with solution at initial time (global)
231: */
232: int InitialConditions(Vec u,AppCtx *appctx)
233: {
234:   PetscScalar *u_localptr,h = appctx->h,x;
235:   int    i,mybase,myend,ierr;

237:   /* 
238:      Determine starting point of each processor's range of
239:      grid values.
240:   */
241:   VecGetOwnershipRange(u,&mybase,&myend);

243:   /* 
244:     Get a pointer to vector data.
245:     - For default PETSc vectors, VecGetArray() returns a pointer to
246:       the data array.  Otherwise, the routine is implementation dependent.
247:     - You MUST call VecRestoreArray() when you no longer need access to
248:       the array.
249:     - Note that the Fortran interface to VecGetArray() differs from the
250:       C version.  See the users manual for details.
251:   */
252:   VecGetArray(u,&u_localptr);

254:   /* 
255:      We initialize the solution array by simply writing the solution
256:      directly into the array locations.  Alternatively, we could use
257:      VecSetValues() or VecSetValuesLocal().
258:   */
259:   for (i=mybase; i<myend; i++) {
260:     x = h*(PetscReal)i; /* current location in global grid */
261:     u_localptr[i-mybase] = 1.0 + x*x;
262:   }

264:   /* 
265:      Restore vector
266:   */
267:   VecRestoreArray(u,&u_localptr);

269:   /* 
270:      Print debugging information if desired
271:   */
272:   if (appctx->debug) {
273:      PetscPrintf(appctx->comm,"initial guess vectorn");
274:      VecView(u,PETSC_VIEWER_STDOUT_WORLD);
275:   }

277:   return 0;
278: }
279: /* --------------------------------------------------------------------- */
280: #undef __FUNCT__
282: /*
283:    ExactSolution - Computes the exact solution at a given time.

285:    Input Parameters:
286:    t - current time
287:    solution - vector in which exact solution will be computed
288:    appctx - user-defined application context

290:    Output Parameter:
291:    solution - vector with the newly computed exact solution
292: */
293: int ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
294: {
295:   PetscScalar *s_localptr,h = appctx->h,x;
296:   int    i,mybase,myend,ierr;

298:   /* 
299:      Determine starting and ending points of each processor's 
300:      range of grid values
301:   */
302:   VecGetOwnershipRange(solution,&mybase,&myend);

304:   /*
305:      Get a pointer to vector data.
306:   */
307:   VecGetArray(solution,&s_localptr);

309:   /* 
310:      Simply write the solution directly into the array locations.
311:      Alternatively, we could use VecSetValues() or VecSetValuesLocal().
312:   */
313:   for (i=mybase; i<myend; i++) {
314:     x = h*(PetscReal)i;
315:     s_localptr[i-mybase] = (t + 1.0)*(1.0 + x*x);
316:   }

318:   /* 
319:      Restore vector
320:   */
321:   VecRestoreArray(solution,&s_localptr);
322:   return 0;
323: }
324: /* --------------------------------------------------------------------- */
325: #undef __FUNCT__
327: /*
328:    Monitor - User-provided routine to monitor the solution computed at 
329:    each timestep.  This example plots the solution and computes the
330:    error in two different norms.

332:    Input Parameters:
333:    ts     - the timestep context
334:    step   - the count of the current step (with 0 meaning the
335:             initial condition)
336:    time   - the current time
337:    u      - the solution at this timestep
338:    ctx    - the user-provided context for this monitoring routine.
339:             In this case we use the application context which contains 
340:             information about the problem size, workspace and the exact 
341:             solution.
342: */
343: int Monitor(TS ts,int step,PetscReal time,Vec u,void *ctx)
344: {
345:   AppCtx       *appctx = (AppCtx*) ctx;   /* user-defined application context */
346:   int          ierr;
347:   PetscReal    en2,en2s,enmax;
348:   PetscScalar  mone = -1.0;
349:   PetscDraw    draw;

351:   /*
352:      We use the default X windows viewer
353:              PETSC_VIEWER_DRAW_(appctx->comm)
354:      that is associated with the current communicator. This saves
355:      the effort of calling PetscViewerDrawOpen() to create the window.
356:      Note that if we wished to plot several items in separate windows we
357:      would create each viewer with PetscViewerDrawOpen() and store them in
358:      the application context, appctx.

360:      PetscReal buffering makes graphics look better.
361:   */
362:   PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(appctx->comm),0,&draw);
363:   PetscDrawSetDoubleBuffer(draw);
364:   VecView(u,PETSC_VIEWER_DRAW_(appctx->comm));

366:   /*
367:      Compute the exact solution at this timestep
368:   */
369:   ExactSolution(time,appctx->solution,appctx);

371:   /*
372:      Print debugging information if desired
373:   */
374:   if (appctx->debug) {
375:      PetscPrintf(appctx->comm,"Computed solution vectorn");
376:      VecView(u,PETSC_VIEWER_STDOUT_WORLD);
377:      PetscPrintf(appctx->comm,"Exact solution vectorn");
378:      VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
379:   }

381:   /*
382:      Compute the 2-norm and max-norm of the error
383:   */
384:   VecAXPY(&mone,u,appctx->solution);
385:   VecNorm(appctx->solution,NORM_2,&en2);
386:   en2s  = sqrt(appctx->h)*en2; /* scale the 2-norm by the grid spacing */
387:   VecNorm(appctx->solution,NORM_MAX,&enmax);

389:   /*
390:      PetscPrintf() causes only the first processor in this 
391:      communicator to print the timestep information.
392:   */
393:   PetscPrintf(appctx->comm,"Timestep %d: time = %g,2-norm error = %g, max norm error = %gn",
394:               step,time,en2s,enmax);

396:   /*
397:      Print debugging information if desired
398:   */
399:   if (appctx->debug) {
400:      PetscPrintf(appctx->comm,"Error vectorn");
401:      VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
402:   }
403:   return 0;
404: }
405: /* --------------------------------------------------------------------- */
406: #undef __FUNCT__
408: /*
409:    RHSFunction - User-provided routine that evalues the right-hand-side
410:    function of the ODE.  This routine is set in the main program by 
411:    calling TSSetRHSFunction().  We compute:
412:           global_out = F(global_in)

414:    Input Parameters:
415:    ts         - timesteping context
416:    t          - current time
417:    global_in  - vector containing the current iterate
418:    ctx        - (optional) user-provided context for function evaluation.
419:                 In this case we use the appctx defined above.

421:    Output Parameter:
422:    global_out - vector containing the newly evaluated function
423: */
424: int RHSFunction(TS ts,PetscReal t,Vec global_in,Vec global_out,void *ctx)
425: {
426:   AppCtx *appctx = (AppCtx*) ctx;       /* user-defined application context */
427:   DA     da = appctx->da;               /* distributed array */
428:   Vec    local_in = appctx->u_local;    /* local ghosted input vector */
429:   Vec    localwork = appctx->localwork; /* local ghosted work vector */
430:   int    ierr,i,localsize,rank,size;
431:   PetscScalar *copyptr,*localptr,sc;

433:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434:      Get ready for local function computations
435:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
436:   /*
437:      Scatter ghost points to local vector, using the 2-step process
438:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
439:      By placing code between these two statements, computations can be
440:      done while messages are in transition.
441:   */
442:   DAGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
443:   DAGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);

445:   /*
446:       Access directly the values in our local INPUT work array
447:   */
448:   VecGetArray(local_in,&localptr);

450:   /*
451:       Access directly the values in our local OUTPUT work array
452:   */
453:   VecGetArray(localwork,&copyptr);

455:   sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));

457:   /*
458:       Evaluate our function on the nodes owned by this processor
459:   */
460:   VecGetLocalSize(local_in,&localsize);

462:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
463:      Compute entries for the locally owned part 
464:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

466:   /*
467:      Handle boundary conditions: This is done by using the boundary condition 
468:         u(t,boundary) = g(t,boundary) 
469:      for some function g. Now take the derivative with respect to t to obtain
470:         u_{t}(t,boundary) = g_{t}(t,boundary)

472:      In our case, u(t,0) = t + 1, so that u_{t}(t,0) = 1 
473:              and  u(t,1) = 2t+ 1, so that u_{t}(t,1) = 2
474:   */
475:   MPI_Comm_rank(appctx->comm,&rank);
476:   MPI_Comm_size(appctx->comm,&size);
477:   if (!rank)          copyptr[0]           = 1.0;
478:   if (rank == size-1) copyptr[localsize-1] = 2.0;

480:   /*
481:      Handle the interior nodes where the PDE is replace by finite 
482:      difference operators.
483:   */
484:   for (i=1; i<localsize-1; i++) {
485:     copyptr[i] =  localptr[i] * sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);
486:   }

488:   /* 
489:      Restore vectors
490:   */
491:   VecRestoreArray(local_in,&localptr);
492:   VecRestoreArray(localwork,&copyptr);

494:   /*
495:      Insert values from the local OUTPUT vector into the global 
496:      output vector
497:   */
498:   DALocalToGlobal(da,localwork,INSERT_VALUES,global_out);

500:   /* Print debugging information if desired */
501:   if (appctx->debug) {
502:      PetscPrintf(appctx->comm,"RHS function vectorn");
503:      VecView(global_out,PETSC_VIEWER_STDOUT_WORLD);
504:   }

506:   return 0;
507: }
508: /* --------------------------------------------------------------------- */
509: #undef __FUNCT__
511: /*
512:    RHSJacobian - User-provided routine to compute the Jacobian of
513:    the nonlinear right-hand-side function of the ODE.

515:    Input Parameters:
516:    ts - the TS context
517:    t - current time
518:    global_in - global input vector
519:    dummy - optional user-defined context, as set by TSetRHSJacobian()

521:    Output Parameters:
522:    AA - Jacobian matrix
523:    BB - optionally different preconditioning matrix
524:    str - flag indicating matrix structure

526:   Notes:
527:   RHSJacobian computes entries for the locally owned part of the Jacobian.
528:    - Currently, all PETSc parallel matrix formats are partitioned by
529:      contiguous chunks of rows across the processors. 
530:    - Each processor needs to insert only elements that it owns
531:      locally (but any non-local elements will be sent to the
532:      appropriate processor during matrix assembly). 
533:    - Always specify global row and columns of matrix entries when
534:      using MatSetValues().
535:    - Here, we set all entries for a particular row at once.
536:    - Note that MatSetValues() uses 0-based row and column numbers
537:      in Fortran as well as in C.
538: */
539: int RHSJacobian(TS ts,PetscReal t,Vec global_in,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
540: {
541:   Mat    A = *AA;                      /* Jacobian matrix */
542:   AppCtx *appctx = (AppCtx*)ctx;     /* user-defined application context */
543:   Vec    local_in = appctx->u_local;   /* local ghosted input vector */
544:   DA     da = appctx->da;              /* distributed array */
545:   PetscScalar v[3],*localptr,sc;
546:   int    ierr,i,mstart,mend,mstarts,mends,idx[3],is;

548:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
549:      Get ready for local Jacobian computations
550:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
551:   /*
552:      Scatter ghost points to local vector, using the 2-step process
553:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
554:      By placing code between these two statements, computations can be
555:      done while messages are in transition.
556:   */
557:   DAGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
558:   DAGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);

560:   /*
561:      Get pointer to vector data
562:   */
563:   VecGetArray(local_in,&localptr);

565:   /* 
566:      Get starting and ending locally owned rows of the matrix
567:   */
568:   MatGetOwnershipRange(A,&mstarts,&mends);
569:   mstart = mstarts; mend = mends;

571:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
572:      Compute entries for the locally owned part of the Jacobian.
573:       - Currently, all PETSc parallel matrix formats are partitioned by
574:         contiguous chunks of rows across the processors. 
575:       - Each processor needs to insert only elements that it owns
576:         locally (but any non-local elements will be sent to the
577:         appropriate processor during matrix assembly). 
578:       - Here, we set all entries for a particular row at once.
579:       - We can set matrix entries either using either
580:         MatSetValuesLocal() or MatSetValues().
581:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

583:   /* 
584:      Set matrix rows corresponding to boundary data
585:   */
586:   if (mstart == 0) {
587:     v[0] = 0.0;
588:     MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
589:     mstart++;
590:   }
591:   if (mend == appctx->m) {
592:     mend--;
593:     v[0] = 0.0;
594:     MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
595:   }

597:   /*
598:      Set matrix rows corresponding to interior data.  We construct the 
599:      matrix one row at a time.
600:   */
601:   sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));
602:   for (i=mstart; i<mend; i++) {
603:     idx[0] = i-1; idx[1] = i; idx[2] = i+1;
604:     is     = i - mstart + 1;
605:     v[0]   = sc*localptr[is];
606:     v[1]   = sc*(localptr[is+1] + localptr[is-1] - 4.0*localptr[is]);
607:     v[2]   = sc*localptr[is];
608:     MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
609:   }

611:   /* 
612:      Restore vector
613:   */
614:   VecRestoreArray(local_in,&localptr);

616:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
617:      Complete the matrix assembly process and set some options
618:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
619:   /*
620:      Assemble matrix, using the 2-step process:
621:        MatAssemblyBegin(), MatAssemblyEnd()
622:      Computations can be done while messages are in transition
623:      by placing code between these two statements.
624:   */
625:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
626:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

628:   /*
629:      Set flag to indicate that the Jacobian matrix retains an identical
630:      nonzero structure throughout all timestepping iterations (although the
631:      values of the entries change). Thus, we can save some work in setting
632:      up the preconditioner (e.g., no need to redo symbolic factorization for
633:      ILU/ICC preconditioners).
634:       - If the nonzero structure of the matrix is different during
635:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
636:         must be used instead.  If you are unsure whether the matrix
637:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
638:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
639:         believes your assertion and does not check the structure
640:         of the matrix.  If you erroneously claim that the structure
641:         is the same when it actually is not, the new preconditioner
642:         will not function correctly.  Thus, use this optimization
643:         feature with caution!
644:   */
645:   *str = SAME_NONZERO_PATTERN;

647:   /*
648:      Set and option to indicate that we will never add a new nonzero location 
649:      to the matrix. If we do, it will generate an error.
650:   */
651:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR);

653:   return 0;
654: }