Actual source code: ex2.c

  2: static char help[] ="Solves a time-dependent nonlinear PDE. Uses implicit\n\
  3: timestepping.  Runtime options include:\n\
  4:   -M <xg>, where <xg> = number of grid points\n\
  5:   -debug : Activate debugging printouts\n\
  6:   -nox   : Deactivate x-window graphics\n\n";

  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:   PetscInt   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: */

 71: /*
 72:    Utility routine for finite difference Jacobian approximation
 73: */

 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:   PetscInt       time_steps_max = 1000;  /* default max timesteps */
 86:   PetscInt       steps;
 87:   PetscReal      ftime;                  /* final time */
 88:   PetscReal      dt;
 89:   PetscReal      time_total_max = 100.0; /* default max total time */
 90:   PetscTruth     flg;

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

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

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

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

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

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

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

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

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

144:   TSMonitorSet(ts,Monitor,&appctx,PETSC_NULL);

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

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

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

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:      Set solution vector and initial timestep
165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

181:   TSSetType(ts,TS_BEULER);
182:   TSSetDuration(ts,time_steps_max,time_total_max);
183:   TSSetFromOptions(ts);

185:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186:      Solve the problem
187:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

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

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

267:   /* 
268:      Restore vector
269:   */
270:   VecRestoreArray(u,&u_localptr);

272:   /* 
273:      Print debugging information if desired
274:   */
275:   if (appctx->debug) {
276:      PetscPrintf(appctx->comm,"initial guess vector\n");
277:      VecView(u,PETSC_VIEWER_STDOUT_WORLD);
278:   }

280:   return 0;
281: }
282: /* --------------------------------------------------------------------- */
285: /*
286:    ExactSolution - Computes the exact solution at a given time.

288:    Input Parameters:
289:    t - current time
290:    solution - vector in which exact solution will be computed
291:    appctx - user-defined application context

293:    Output Parameter:
294:    solution - vector with the newly computed exact solution
295: */
296: PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
297: {
298:   PetscScalar    *s_localptr,h = appctx->h,x;
299:   PetscInt       i,mybase,myend;

302:   /* 
303:      Determine starting and ending points of each processor's 
304:      range of grid values
305:   */
306:   VecGetOwnershipRange(solution,&mybase,&myend);

308:   /*
309:      Get a pointer to vector data.
310:   */
311:   VecGetArray(solution,&s_localptr);

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

322:   /* 
323:      Restore vector
324:   */
325:   VecRestoreArray(solution,&s_localptr);
326:   return 0;
327: }
328: /* --------------------------------------------------------------------- */
331: /*
332:    Monitor - User-provided routine to monitor the solution computed at 
333:    each timestep.  This example plots the solution and computes the
334:    error in two different norms.

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

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

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

369:   /*
370:      Compute the exact solution at this timestep
371:   */
372:   ExactSolution(time,appctx->solution,appctx);

374:   /*
375:      Print debugging information if desired
376:   */
377:   if (appctx->debug) {
378:      PetscPrintf(appctx->comm,"Computed solution vector\n");
379:      VecView(u,PETSC_VIEWER_STDOUT_WORLD);
380:      PetscPrintf(appctx->comm,"Exact solution vector\n");
381:      VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
382:   }

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

392:   /*
393:      PetscPrintf() causes only the first processor in this 
394:      communicator to print the timestep information.
395:   */
396:   PetscPrintf(appctx->comm,"Timestep %D: time = %G,2-norm error = %G, max norm error = %G\n",
397:               step,time,en2s,enmax);

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

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

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

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

450:   /*
451:       Access directly the values in our local INPUT work array
452:   */
453:   VecGetArray(local_in,&localptr);

455:   /*
456:       Access directly the values in our local OUTPUT work array
457:   */
458:   VecGetArray(localwork,&copyptr);

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

462:   /*
463:       Evaluate our function on the nodes owned by this processor
464:   */
465:   VecGetLocalSize(local_in,&localsize);

467:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
468:      Compute entries for the locally owned part 
469:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

471:   /*
472:      Handle boundary conditions: This is done by using the boundary condition 
473:         u(t,boundary) = g(t,boundary) 
474:      for some function g. Now take the derivative with respect to t to obtain
475:         u_{t}(t,boundary) = g_{t}(t,boundary)

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

485:   /*
486:      Handle the interior nodes where the PDE is replace by finite 
487:      difference operators.
488:   */
489:   for (i=1; i<localsize-1; i++) {
490:     copyptr[i] =  localptr[i] * sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);
491:   }

493:   /* 
494:      Restore vectors
495:   */
496:   VecRestoreArray(local_in,&localptr);
497:   VecRestoreArray(localwork,&copyptr);

499:   /*
500:      Insert values from the local OUTPUT vector into the global 
501:      output vector
502:   */
503:   DALocalToGlobal(da,localwork,INSERT_VALUES,global_out);

505:   /* Print debugging information if desired */
506:   if (appctx->debug) {
507:      PetscPrintf(appctx->comm,"RHS function vector\n");
508:      VecView(global_out,PETSC_VIEWER_STDOUT_WORLD);
509:   }

511:   return 0;
512: }
513: /* --------------------------------------------------------------------- */
516: /*
517:    RHSJacobian - User-provided routine to compute the Jacobian of
518:    the nonlinear right-hand-side function of the ODE.

520:    Input Parameters:
521:    ts - the TS context
522:    t - current time
523:    global_in - global input vector
524:    dummy - optional user-defined context, as set by TSetRHSJacobian()

526:    Output Parameters:
527:    AA - Jacobian matrix
528:    BB - optionally different preconditioning matrix
529:    str - flag indicating matrix structure

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

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

566:   /*
567:      Get pointer to vector data
568:   */
569:   VecGetArray(local_in,&localptr);

571:   /* 
572:      Get starting and ending locally owned rows of the matrix
573:   */
574:   MatGetOwnershipRange(A,&mstarts,&mends);
575:   mstart = mstarts; mend = mends;

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

589:   /* 
590:      Set matrix rows corresponding to boundary data
591:   */
592:   if (mstart == 0) {
593:     v[0] = 0.0;
594:     MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
595:     mstart++;
596:   }
597:   if (mend == appctx->m) {
598:     mend--;
599:     v[0] = 0.0;
600:     MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
601:   }

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

617:   /* 
618:      Restore vector
619:   */
620:   VecRestoreArray(local_in,&localptr);

622:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
623:      Complete the matrix assembly process and set some options
624:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
625:   /*
626:      Assemble matrix, using the 2-step process:
627:        MatAssemblyBegin(), MatAssemblyEnd()
628:      Computations can be done while messages are in transition
629:      by placing code between these two statements.
630:   */
631:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
632:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

653:   /*
654:      Set and option to indicate that we will never add a new nonzero location 
655:      to the matrix. If we do, it will generate an error.
656:   */
657:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR);

659:   return 0;
660: }