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:   TSSetMonitor(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:   PetscScalar    mone = -1.0;
    353:   PetscDraw      draw;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    660:   return 0;
    661: }