Actual source code: ex21.c

petsc-3.4.4 2014-03-13
  2: static char help[] ="Solves a time-dependent nonlinear PDE with lower and upper bounds on the interior grid points. 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\
  7:   -ul   : lower bound\n\
  8:   -uh  : upper bound\n\n";

 10: /*
 11:    Concepts: TS^time-dependent nonlinear problems
 12:    Concepts: TS^Variational inequality nonlinear solver
 13:    Processors: n
 14: */

 16: /* ------------------------------------------------------------------------

 18:    This is a variation of ex2.c to solve the PDE

 20:                u * u_xx
 21:          u_t = ---------
 22:                2*(t+1)^2

 24:     with box constraints on the interior grid points
 25:     ul <= u(t,x) <= uh with x != 0,1
 26:     on the domain 0 <= x <= 1, with boundary conditions
 27:          u(t,0) = t + 1,  u(t,1) = 2*t + 2,
 28:     and initial condition
 29:          u(0,x) = 1 + x*x.

 31:     The exact solution is:
 32:          u(t,x) = (1 + x*x) * (1 + t)

 34:     We use by default the backward Euler method.

 36:   ------------------------------------------------------------------------- */

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

 43:    Include the "petscdmda.h" to allow us to use the distributed array data
 44:    structures to manage the parallel grid.
 45: */
 46: #include <petscts.h>
 47: #include <petscdmda.h>
 48: #include <petscdraw.h>

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

 65: /*
 66:    User-defined routines, provided below.
 67: */
 68: extern PetscErrorCode InitialConditions(Vec,AppCtx*);
 69: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
 70: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
 71: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
 72: extern PetscErrorCode ExactSolution(PetscReal,Vec,AppCtx*);
 73: extern PetscErrorCode SetBounds(Vec,Vec,PetscScalar,PetscScalar,AppCtx*);

 77: int main(int argc,char **argv)
 78: {
 79:   AppCtx         appctx;                 /* user-defined application context */
 80:   TS             ts;                     /* timestepping context */
 81:   Mat            A;                      /* Jacobian matrix data structure */
 82:   Vec            u;                      /* approximate solution vector */
 83:   Vec            r;                      /* residual vector */
 84:   PetscInt       time_steps_max = 1000;  /* default max timesteps */
 86:   PetscReal      dt;
 87:   PetscReal      time_total_max = 100.0; /* default max total time */
 88:   Vec            xl,xu; /* Lower and upper bounds on variables */
 89:   PetscScalar    ul=0.0,uh = 3.0;
 90:   PetscBool      mymonitor;
 91:   PetscReal      bounds[] = {1.0, 3.3};

 93:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:      Initialize program and set problem parameters
 95:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 97:   PetscInitialize(&argc,&argv,(char*)0,help);
 98:   PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,bounds);

100:   appctx.comm = PETSC_COMM_WORLD;
101:   appctx.m    = 60;
102:   PetscOptionsGetInt(NULL,"-M",&appctx.m,NULL);
103:   PetscOptionsGetScalar(NULL,"-ul",&ul,NULL);
104:   PetscOptionsGetScalar(NULL,"-uh",&uh,NULL);
105:   PetscOptionsHasName(NULL,"-debug",&appctx.debug);
106:   PetscOptionsHasName(NULL,"-mymonitor",&mymonitor);
107:   appctx.h    = 1.0/(appctx.m-1.0);

109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110:      Create vector data structures
111:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

113:   /*
114:      Create distributed array (DMDA) to manage parallel grid and vectors
115:      and to set up the ghost point communication pattern.  There are M
116:      total grid values spread equally among all the processors.
117:   */
118:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,appctx.m,1,1,NULL,
119:                       &appctx.da);

121:   /*
122:      Extract global and local vectors from DMDA; we use these to store the
123:      approximate solution.  Then duplicate these for remaining vectors that
124:      have the same types.
125:   */
126:   DMCreateGlobalVector(appctx.da,&u);
127:   DMCreateLocalVector(appctx.da,&appctx.u_local);

129:   /*
130:      Create local work vector for use in evaluating right-hand-side function;
131:      create global work vector for storing exact solution.
132:   */
133:   VecDuplicate(appctx.u_local,&appctx.localwork);
134:   VecDuplicate(u,&appctx.solution);

136:   /* Create residual vector */
137:   VecDuplicate(u,&r);
138:   /* Create lower and upper bound vectors */
139:   VecDuplicate(u,&xl);
140:   VecDuplicate(u,&xu);
141:   SetBounds(xl,xu,ul,uh,&appctx);

143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:      Create timestepping solver context; set callback routine for
145:      right-hand-side function evaluation.
146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

148:   TSCreate(PETSC_COMM_WORLD,&ts);
149:   TSSetProblemType(ts,TS_NONLINEAR);
150:   TSSetRHSFunction(ts,r,RHSFunction,&appctx);

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:      Set optional user-defined monitoring routine
154:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

156:   if (mymonitor) {
157:     TSMonitorSet(ts,Monitor,&appctx,NULL);
158:   }

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

164:      Create matrix data structure; set Jacobian evaluation routine.
165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

167:   MatCreate(PETSC_COMM_WORLD,&A);
168:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,appctx.m,appctx.m);
169:   MatSetFromOptions(A);
170:   MatSetUp(A);
171:   TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx);

173:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174:      Set solution vector and initial timestep
175:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

177:   dt   = appctx.h/2.0;
178:   TSSetInitialTimeStep(ts,0.0,dt);

180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181:      Customize timestepping solver:
182:        - Set the solution method to be the Backward Euler method.
183:        - Set timestepping duration info
184:      Then set runtime options, which can override these defaults.
185:      For example,
186:           -ts_max_steps <maxsteps> -ts_final_time <maxtime>
187:      to override the defaults set by TSSetDuration().
188:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

190:   TSSetType(ts,TSBEULER);
191:   TSSetDuration(ts,time_steps_max,time_total_max);
192:   /* Set lower and upper bound on the solution vector for each time step */
193:   TSVISetVariableBounds(ts,xl,xu);
194:   TSSetFromOptions(ts);

196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197:      Solve the problem
198:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

200:   /*
201:      Evaluate initial conditions
202:   */
203:   InitialConditions(u,&appctx);

205:   /*
206:      Run the timestepping solver
207:   */
208:   TSSolve(ts,u);

210:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211:      Free work space.  All PETSc objects should be destroyed when they
212:      are no longer needed.
213:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

215:   VecDestroy(&r);
216:   VecDestroy(&xl);
217:   VecDestroy(&xu);
218:   TSDestroy(&ts);
219:   VecDestroy(&u);
220:   MatDestroy(&A);
221:   DMDestroy(&appctx.da);
222:   VecDestroy(&appctx.localwork);
223:   VecDestroy(&appctx.solution);
224:   VecDestroy(&appctx.u_local);

226:   /*
227:      Always call PetscFinalize() before exiting a program.  This routine
228:        - finalizes the PETSc libraries as well as MPI
229:        - provides summary and diagnostic information if certain runtime
230:          options are chosen (e.g., -log_summary).
231:   */
232:   PetscFinalize();
233:   return 0;
234: }
235: /* --------------------------------------------------------------------- */
238: /*
239:    InitialConditions - Computes the solution at the initial time.

241:    Input Parameters:
242:    u - uninitialized solution vector (global)
243:    appctx - user-defined application context

245:    Output Parameter:
246:    u - vector with solution at initial time (global)
247: */
248: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
249: {
250:   PetscScalar    *u_localptr,h = appctx->h,x;
251:   PetscInt       i,mybase,myend;

254:   /*
255:      Determine starting point of each processor's range of
256:      grid values.
257:   */
258:   VecGetOwnershipRange(u,&mybase,&myend);

260:   /*
261:     Get a pointer to vector data.
262:     - For default PETSc vectors, VecGetArray() returns a pointer to
263:       the data array.  Otherwise, the routine is implementation dependent.
264:     - You MUST call VecRestoreArray() when you no longer need access to
265:       the array.
266:     - Note that the Fortran interface to VecGetArray() differs from the
267:       C version.  See the users manual for details.
268:   */
269:   VecGetArray(u,&u_localptr);

271:   /*
272:      We initialize the solution array by simply writing the solution
273:      directly into the array locations.  Alternatively, we could use
274:      VecSetValues() or VecSetValuesLocal().
275:   */
276:   for (i=mybase; i<myend; i++) {
277:     x = h*(PetscReal)i; /* current location in global grid */
278:     u_localptr[i-mybase] = 1.0 + x*x;
279:   }

281:   /*
282:      Restore vector
283:   */
284:   VecRestoreArray(u,&u_localptr);

286:   /*
287:      Print debugging information if desired
288:   */
289:   if (appctx->debug) {
290:      PetscPrintf(appctx->comm,"initial guess vector\n");
291:      VecView(u,PETSC_VIEWER_STDOUT_WORLD);
292:   }

294:   return 0;
295: }

297: /* --------------------------------------------------------------------- */
300: /*
301:   SetBounds - Sets the lower and uper bounds on the interior points

303:   Input parameters:
304:   xl - vector of lower bounds
305:   xu - vector of upper bounds
306:   ul - constant lower bound for all variables
307:   uh - constant upper bound for all variables
308:   appctx - Application context
309:  */
310: PetscErrorCode SetBounds(Vec xl, Vec xu, PetscScalar ul, PetscScalar uh,AppCtx *appctx)
311: {
313:   PetscScalar    *l,*u;
314:   PetscMPIInt    rank,size;
315:   PetscInt       localsize;

318:   VecSet(xl,ul);
319:   VecSet(xu,uh);
320:   VecGetLocalSize(xl,&localsize);
321:   VecGetArray(xl,&l);
322:   VecGetArray(xu,&u);


325:   MPI_Comm_rank(appctx->comm,&rank);
326:   MPI_Comm_size(appctx->comm,&size);
327:   if (!rank) {
328:     l[0] = -SNES_VI_INF;
329:     u[0] =  SNES_VI_INF;
330:   }
331:   if (rank == size-1) {
332:     l[localsize-1] = -SNES_VI_INF;
333:     u[localsize-1] = SNES_VI_INF;
334:   }
335:   VecRestoreArray(xl,&l);
336:   VecRestoreArray(xu,&u);
337:   return(0);
338: }

340: /* --------------------------------------------------------------------- */
343: /*
344:    ExactSolution - Computes the exact solution at a given time.

346:    Input Parameters:
347:    t - current time
348:    solution - vector in which exact solution will be computed
349:    appctx - user-defined application context

351:    Output Parameter:
352:    solution - vector with the newly computed exact solution
353: */
354: PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
355: {
356:   PetscScalar    *s_localptr,h = appctx->h,x;
357:   PetscInt       i,mybase,myend;

360:   /*
361:      Determine starting and ending points of each processor's
362:      range of grid values
363:   */
364:   VecGetOwnershipRange(solution,&mybase,&myend);

366:   /*
367:      Get a pointer to vector data.
368:   */
369:   VecGetArray(solution,&s_localptr);

371:   /*
372:      Simply write the solution directly into the array locations.
373:      Alternatively, we could use VecSetValues() or VecSetValuesLocal().
374:   */
375:   for (i=mybase; i<myend; i++) {
376:     x = h*(PetscReal)i;
377:     s_localptr[i-mybase] = (t + 1.0)*(1.0 + x*x);
378:   }

380:   /*
381:      Restore vector
382:   */
383:   VecRestoreArray(solution,&s_localptr);
384:   return 0;
385: }
386: /* --------------------------------------------------------------------- */
389: /*
390:    Monitor - User-provided routine to monitor the solution computed at
391:    each timestep.  This example plots the solution and computes the
392:    error in two different norms.

394:    Input Parameters:
395:    ts     - the timestep context
396:    step   - the count of the current step (with 0 meaning the
397:             initial condition)
398:    time   - the current time
399:    u      - the solution at this timestep
400:    ctx    - the user-provided context for this monitoring routine.
401:             In this case we use the application context which contains
402:             information about the problem size, workspace and the exact
403:             solution.
404: */
405: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
406: {
407:   AppCtx         *appctx = (AppCtx*) ctx;   /* user-defined application context */
409:   PetscReal      en2,en2s,enmax;
410:   PetscDraw      draw;

412:   /*
413:      We use the default X windows viewer
414:              PETSC_VIEWER_DRAW_(appctx->comm)
415:      that is associated with the current communicator. This saves
416:      the effort of calling PetscViewerDrawOpen() to create the window.
417:      Note that if we wished to plot several items in separate windows we
418:      would create each viewer with PetscViewerDrawOpen() and store them in
419:      the application context, appctx.

421:      PetscReal buffering makes graphics look better.
422:   */
423:   PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(appctx->comm),0,&draw);
424:   PetscDrawSetDoubleBuffer(draw);
425:   VecView(u,PETSC_VIEWER_DRAW_(appctx->comm));

427:   /*
428:      Compute the exact solution at this timestep
429:   */
430:   ExactSolution(time,appctx->solution,appctx);

432:   /*
433:      Print debugging information if desired
434:   */
435:   if (appctx->debug) {
436:     PetscPrintf(appctx->comm,"Computed solution vector\n");
437:     VecView(u,PETSC_VIEWER_STDOUT_WORLD);
438:     PetscPrintf(appctx->comm,"Exact solution vector\n");
439:     VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
440:   }

442:   /*
443:      Compute the 2-norm and max-norm of the error
444:   */
445:   VecAXPY(appctx->solution,-1.0,u);
446:   VecNorm(appctx->solution,NORM_2,&en2);
447:   en2s = PetscSqrtReal(appctx->h)*en2;  /* scale the 2-norm by the grid spacing */
448:   VecNorm(appctx->solution,NORM_MAX,&enmax);

450:   /*
451:      PetscPrintf() causes only the first processor in this
452:      communicator to print the timestep information.
453:   */
454:   PetscPrintf(appctx->comm,"Timestep %D: time = %G,2-norm error = %G, max norm error = %G\n",
455:                      step,time,en2s,enmax);

457:   /*
458:      Print debugging information if desired
459:    */
460:   /*  if (appctx->debug) {
461:      PetscPrintf(appctx->comm,"Error vector\n");
462:      VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
463:    } */
464:   return 0;
465: }
466: /* --------------------------------------------------------------------- */
469: /*
470:    RHSFunction - User-provided routine that evalues the right-hand-side
471:    function of the ODE.  This routine is set in the main program by
472:    calling TSSetRHSFunction().  We compute:
473:           global_out = F(global_in)

475:    Input Parameters:
476:    ts         - timesteping context
477:    t          - current time
478:    global_in  - vector containing the current iterate
479:    ctx        - (optional) user-provided context for function evaluation.
480:                 In this case we use the appctx defined above.

482:    Output Parameter:
483:    global_out - vector containing the newly evaluated function
484: */
485: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec global_in,Vec global_out,void *ctx)
486: {
487:   AppCtx         *appctx   = (AppCtx*) ctx;     /* user-defined application context */
488:   DM             da        = appctx->da;        /* distributed array */
489:   Vec            local_in  = appctx->u_local;   /* local ghosted input vector */
490:   Vec            localwork = appctx->localwork; /* local ghosted work vector */
492:   PetscInt       i,localsize;
493:   PetscMPIInt    rank,size;
494:   PetscScalar    *copyptr,*localptr,sc;

496:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
497:      Get ready for local function computations
498:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
499:   /*
500:      Scatter ghost points to local vector, using the 2-step process
501:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
502:      By placing code between these two statements, computations can be
503:      done while messages are in transition.
504:   */
505:   DMGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
506:   DMGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);

508:   /*
509:       Access directly the values in our local INPUT work array
510:   */
511:   VecGetArray(local_in,&localptr);

513:   /*
514:       Access directly the values in our local OUTPUT work array
515:   */
516:   VecGetArray(localwork,&copyptr);

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

520:   /*
521:       Evaluate our function on the nodes owned by this processor
522:   */
523:   VecGetLocalSize(local_in,&localsize);

525:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
526:      Compute entries for the locally owned part
527:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

529:   /*
530:      Handle boundary conditions: This is done by using the boundary condition
531:         u(t,boundary) = g(t,boundary)
532:      for some function g. Now take the derivative with respect to t to obtain
533:         u_{t}(t,boundary) = g_{t}(t,boundary)

535:      In our case, u(t,0) = t + 1, so that u_{t}(t,0) = 1
536:              and  u(t,1) = 2t+ 2, so that u_{t}(t,1) = 2
537:   */
538:   MPI_Comm_rank(appctx->comm,&rank);
539:   MPI_Comm_size(appctx->comm,&size);
540:   if (!rank) copyptr[0] = 1.0;
541:   if (rank == size-1) copyptr[localsize-1] = (t < .5) ? 2.0 : 0.0;

543:   /*
544:      Handle the interior nodes where the PDE is replace by finite
545:      difference operators.
546:   */
547:   for (i=1; i<localsize-1; i++) copyptr[i] =  localptr[i] * sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);

549:   /*
550:      Restore vectors
551:   */
552:   VecRestoreArray(local_in,&localptr);
553:   VecRestoreArray(localwork,&copyptr);

555:   /*
556:      Insert values from the local OUTPUT vector into the global
557:      output vector
558:   */
559:   DMLocalToGlobalBegin(da,localwork,INSERT_VALUES,global_out);
560:   DMLocalToGlobalEnd(da,localwork,INSERT_VALUES,global_out);

562:   /* Print debugging information if desired */
563:   /*  if (appctx->debug) {
564:      PetscPrintf(appctx->comm,"RHS function vector\n");
565:      VecView(global_out,PETSC_VIEWER_STDOUT_WORLD);
566:    } */

568:   return 0;
569: }
570: /* --------------------------------------------------------------------- */
573: /*
574:    RHSJacobian - User-provided routine to compute the Jacobian of
575:    the nonlinear right-hand-side function of the ODE.

577:    Input Parameters:
578:    ts - the TS context
579:    t - current time
580:    global_in - global input vector
581:    dummy - optional user-defined context, as set by TSetRHSJacobian()

583:    Output Parameters:
584:    AA - Jacobian matrix
585:    BB - optionally different preconditioning matrix
586:    str - flag indicating matrix structure

588:   Notes:
589:   RHSJacobian computes entries for the locally owned part of the Jacobian.
590:    - Currently, all PETSc parallel matrix formats are partitioned by
591:      contiguous chunks of rows across the processors.
592:    - Each processor needs to insert only elements that it owns
593:      locally (but any non-local elements will be sent to the
594:      appropriate processor during matrix assembly).
595:    - Always specify global row and columns of matrix entries when
596:      using MatSetValues().
597:    - Here, we set all entries for a particular row at once.
598:    - Note that MatSetValues() uses 0-based row and column numbers
599:      in Fortran as well as in C.
600: */
601: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec global_in,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
602: {
603:   Mat            B        = *BB;               /* Jacobian matrix */
604:   AppCtx         *appctx  = (AppCtx*)ctx;    /* user-defined application context */
605:   Vec            local_in = appctx->u_local;   /* local ghosted input vector */
606:   DM             da       = appctx->da;        /* distributed array */
607:   PetscScalar    v[3],*localptr,sc;
609:   PetscInt       i,mstart,mend,mstarts,mends,idx[3],is;

611:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
612:      Get ready for local Jacobian computations
613:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
614:   /*
615:      Scatter ghost points to local vector, using the 2-step process
616:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
617:      By placing code between these two statements, computations can be
618:      done while messages are in transition.
619:   */
620:   DMGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
621:   DMGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);

623:   /*
624:      Get pointer to vector data
625:   */
626:   VecGetArray(local_in,&localptr);

628:   /*
629:      Get starting and ending locally owned rows of the matrix
630:   */
631:   MatGetOwnershipRange(B,&mstarts,&mends);
632:   mstart = mstarts; mend = mends;

634:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
635:      Compute entries for the locally owned part of the Jacobian.
636:       - Currently, all PETSc parallel matrix formats are partitioned by
637:         contiguous chunks of rows across the processors.
638:       - Each processor needs to insert only elements that it owns
639:         locally (but any non-local elements will be sent to the
640:         appropriate processor during matrix assembly).
641:       - Here, we set all entries for a particular row at once.
642:       - We can set matrix entries either using either
643:         MatSetValuesLocal() or MatSetValues().
644:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

646:   /*
647:      Set matrix rows corresponding to boundary data
648:   */
649:   if (mstart == 0) {
650:     v[0] = 0.0;
651:     MatSetValues(B,1,&mstart,1,&mstart,v,INSERT_VALUES);
652:     mstart++;
653:   }
654:   if (mend == appctx->m) {
655:     mend--;
656:     v[0] = 0.0;
657:     MatSetValues(B,1,&mend,1,&mend,v,INSERT_VALUES);
658:   }

660:   /*
661:      Set matrix rows corresponding to interior data.  We construct the
662:      matrix one row at a time.
663:   */
664:   sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));
665:   for (i=mstart; i<mend; i++) {
666:     idx[0] = i-1; idx[1] = i; idx[2] = i+1;
667:     is     = i - mstart + 1;
668:     v[0]   = sc*localptr[is];
669:     v[1]   = sc*(localptr[is+1] + localptr[is-1] - 4.0*localptr[is]);
670:     v[2]   = sc*localptr[is];
671:     MatSetValues(B,1,&i,3,idx,v,INSERT_VALUES);
672:   }

674:   /*
675:      Restore vector
676:   */
677:   VecRestoreArray(local_in,&localptr);

679:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
680:      Complete the matrix assembly process and set some options
681:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
682:   /*
683:      Assemble matrix, using the 2-step process:
684:        MatAssemblyBegin(), MatAssemblyEnd()
685:      Computations can be done while messages are in transition
686:      by placing code between these two statements.
687:   */
688:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
689:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

691:   /*
692:      Set flag to indicate that the Jacobian matrix retains an identical
693:      nonzero structure throughout all timestepping iterations (although the
694:      values of the entries change). Thus, we can save some work in setting
695:      up the preconditioner (e.g., no need to redo symbolic factorization for
696:      ILU/ICC preconditioners).
697:       - If the nonzero structure of the matrix is different during
698:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
699:         must be used instead.  If you are unsure whether the matrix
700:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
701:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
702:         believes your assertion and does not check the structure
703:         of the matrix.  If you erroneously claim that the structure
704:         is the same when it actually is not, the new preconditioner
705:         will not function correctly.  Thus, use this optimization
706:         feature with caution!
707:   */
708:   *str = SAME_NONZERO_PATTERN;

710:   /*
711:      Set and option to indicate that we will never add a new nonzero location
712:      to the matrix. If we do, it will generate an error.
713:   */
714:   MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);

716:   return 0;
717: }