Actual source code: ex3.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  2: static char help[] = "Solves 1D heat equation with FEM formulation.\n\
  3: Input arguments are\n\
  4:   -useAlhs: solve Alhs*U' =  (Arhs*U + g) \n\
  5:             otherwise, solve U' = inv(Alhs)*(Arhs*U + g) \n\n";

  7: /*--------------------------------------------------------------------------
  8:   Solves 1D heat equation U_t = U_xx with FEM formulation:
  9:                           Alhs*U' = rhs (= Arhs*U + g)
 10:   We thank Chris Cox <clcox@clemson.edu> for contributing the original code
 11: ----------------------------------------------------------------------------*/

 13: #include <petscksp.h>
 14: #include <petscts.h>

 16: /* special variable - max size of all arrays  */
 17: #define num_z 60

 19: /*
 20:    User-defined application context - contains data needed by the
 21:    application-provided call-back routines.
 22: */
 23: typedef struct {
 24:   Mat         Amat;               /* left hand side matrix */
 25:   Vec         ksp_rhs,ksp_sol;    /* working vectors for formulating inv(Alhs)*(Arhs*U+g) */
 26:   int         max_probsz;         /* max size of the problem */
 27:   PetscBool   useAlhs;            /* flag (1 indicates solving Alhs*U' = Arhs*U+g */
 28:   int         nz;                 /* total number of grid points */
 29:   PetscInt    m;                  /* total number of interio grid points */
 30:   Vec         solution;           /* global exact ts solution vector */
 31:   PetscScalar *z;                 /* array of grid points */
 32:   PetscBool   debug;              /* flag (1 indicates activation of debugging printouts) */
 33: } AppCtx;

 35: extern PetscScalar exact(PetscScalar,PetscReal);
 36: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
 37: extern PetscErrorCode Petsc_KSPSolve(AppCtx*);
 38: extern PetscScalar bspl(PetscScalar*,PetscScalar,PetscInt,PetscInt,PetscInt[][2],PetscInt);
 39: extern PetscErrorCode femBg(PetscScalar[][3],PetscScalar*,PetscInt,PetscScalar*,PetscReal);
 40: extern PetscErrorCode femA(AppCtx*,PetscInt,PetscScalar*);
 41: extern PetscErrorCode rhs(AppCtx*,PetscScalar*, PetscInt, PetscScalar*,PetscReal);
 42: extern PetscErrorCode RHSfunction(TS,PetscReal,Vec,Vec,void*);

 46: int main(int argc,char **argv)
 47: {
 48:   PetscInt       i,m,nz,steps,max_steps,k,nphase=1;
 49:   PetscScalar    zInitial,zFinal,val,*z;
 50:   PetscReal      stepsz[4],T,ftime;
 52:   TS             ts;
 53:   SNES           snes;
 54:   Mat            Jmat;
 55:   AppCtx         appctx;   /* user-defined application context */
 56:   Vec            init_sol; /* ts solution vector */
 57:   PetscMPIInt    size;

 59:   PetscInitialize(&argc,&argv,(char*)0,help);
 60:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 61:   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only");

 63:   PetscOptionsGetInt(NULL,"-nphase",&nphase,NULL);
 64:   if (nphase > 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nphase must be an integer between 1 and 3");

 66:   /* initializations */
 67:   zInitial  = 0.0;
 68:   zFinal    = 1.0;
 69:   T         = 0.014/nphase;
 70:   nz        = num_z;
 71:   m         = nz-2;
 72:   appctx.nz = nz;
 73:   max_steps = (PetscInt)10000;

 75:   appctx.m          = m;
 76:   appctx.max_probsz = nz;
 77:   appctx.debug      = PETSC_FALSE;
 78:   appctx.useAlhs    = PETSC_FALSE;

 80:   PetscOptionsHasName(NULL,"-debug",&appctx.debug);
 81:   PetscOptionsHasName(NULL,"-useAlhs",&appctx.useAlhs);

 83:   /* create vector to hold ts solution */
 84:   /*-----------------------------------*/
 85:   VecCreate(PETSC_COMM_WORLD, &init_sol);
 86:   VecSetSizes(init_sol, PETSC_DECIDE, m);
 87:   VecSetFromOptions(init_sol);

 89:   /* create vector to hold true ts soln for comparison */
 90:   VecDuplicate(init_sol, &appctx.solution);

 92:   /* create LHS matrix Amat */
 93:   /*------------------------*/
 94:   MatCreateSeqAIJ(PETSC_COMM_WORLD, m, m, 3, NULL, &appctx.Amat);
 95:   MatSetFromOptions(appctx.Amat);
 96:   MatSetUp(appctx.Amat);
 97:   /* set space grid points - interio points only! */
 98:   PetscMalloc1((nz+1),&z);
 99:   for (i=0; i<nz; i++) z[i]=(i)*((zFinal-zInitial)/(nz-1));
100:   appctx.z = z;
101:   femA(&appctx,nz,z);

103:   /* create the jacobian matrix */
104:   /*----------------------------*/
105:   MatCreate(PETSC_COMM_WORLD, &Jmat);
106:   MatSetSizes(Jmat,PETSC_DECIDE,PETSC_DECIDE,m,m);
107:   MatSetFromOptions(Jmat);
108:   MatSetUp(Jmat);

110:   /* create working vectors for formulating rhs=inv(Alhs)*(Arhs*U + g) */
111:   VecDuplicate(init_sol,&appctx.ksp_rhs);
112:   VecDuplicate(init_sol,&appctx.ksp_sol);

114:   /* set intial guess */
115:   /*------------------*/
116:   for (i=0; i<nz-2; i++) {
117:     val  = exact(z[i+1], 0.0);
118:     VecSetValue(init_sol,i,(PetscScalar)val,INSERT_VALUES);
119:   }
120:   VecAssemblyBegin(init_sol);
121:   VecAssemblyEnd(init_sol);

123:   /*create a time-stepping context and set the problem type */
124:   /*--------------------------------------------------------*/
125:   TSCreate(PETSC_COMM_WORLD, &ts);
126:   TSSetProblemType(ts,TS_NONLINEAR);

128:   /* set time-step method */
129:   TSSetType(ts,TSCN);

131:   /* Set optional user-defined monitoring routine */
132:   TSMonitorSet(ts,Monitor,&appctx,NULL);
133:   /* set the right hand side of U_t = RHSfunction(U,t) */
134:   TSSetRHSFunction(ts,NULL,(PetscErrorCode (*)(TS,PetscScalar,Vec,Vec,void*))RHSfunction,&appctx);

136:   if (appctx.useAlhs) {
137:     /* set the left hand side matrix of Amat*U_t = rhs(U,t) */

139:     /* Note: this approach is incompatible with the finite differenced Jacobian set below because we can't restore the
140:      * Alhs matrix without making a copy.  Either finite difference the entire thing or use analytic Jacobians in both
141:      * places.
142:      */
143:     TSSetIFunction(ts,NULL,TSComputeIFunctionLinear,&appctx);
144:     TSSetIJacobian(ts,appctx.Amat,appctx.Amat,TSComputeIJacobianConstant,&appctx);
145:   }

147:   /* use petsc to compute the jacobian by finite differences */
148:   TSGetSNES(ts,&snes);
149:   SNESSetJacobian(snes,Jmat,Jmat,SNESComputeJacobianDefault,NULL);

151:   /* get the command line options if there are any and set them */
152:   TSSetFromOptions(ts);

154: #if defined(PETSC_HAVE_SUNDIALS)
155:   {
156:     TSType    type;
157:     PetscBool sundialstype=PETSC_FALSE;
158:     TSGetType(ts,&type);
159:     PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundialstype);
160:     if (sundialstype && appctx.useAlhs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use Alhs formulation for TSSUNDIALS type");
161:   }
162: #endif
163:   /* Sets the initial solution */
164:   TSSetSolution(ts,init_sol);

166:   stepsz[0] = 1.0/(2.0*(nz-1)*(nz-1)); /* (mesh_size)^2/2.0 */
167:   ftime     = 0.0;
168:   for (k=0; k<nphase; k++) {
169:     if (nphase > 1) {PetscPrintf(PETSC_COMM_WORLD,"Phase %D initial time %g, stepsz %g, duration: %g\n",k,(double)ftime,(double)stepsz[k],(double)((k+1)*T));}
170:     TSSetInitialTimeStep(ts,ftime,stepsz[k]);
171:     TSSetDuration(ts,max_steps,(k+1)*T);

173:     /* loop over time steps */
174:     /*----------------------*/
175:     TSSolve(ts,init_sol);
176:     TSGetSolveTime(ts,&ftime);
177:     TSGetTimeStepNumber(ts,&steps);
178:     stepsz[k+1] = stepsz[k]*1.5; /* change step size for the next phase */
179:   }

181:   /* free space */
182:   TSDestroy(&ts);
183:   MatDestroy(&appctx.Amat);
184:   MatDestroy(&Jmat);
185:   VecDestroy(&appctx.ksp_rhs);
186:   VecDestroy(&appctx.ksp_sol);
187:   VecDestroy(&init_sol);
188:   VecDestroy(&appctx.solution);
189:   PetscFree(z);

191:   PetscFinalize();
192:   return 0;
193:   }

195: /*------------------------------------------------------------------------
196:   Set exact solution
197:   u(z,t) = sin(6*PI*z)*exp(-36.*PI*PI*t) + 3.*sin(2*PI*z)*exp(-4.*PI*PI*t)
198: --------------------------------------------------------------------------*/
199: PetscScalar exact(PetscScalar z,PetscReal t)
200: {
201:   PetscScalar val, ex1, ex2;

203:   ex1 = PetscExpReal(-36.*PETSC_PI*PETSC_PI*t);
204:   ex2 = PetscExpReal(-4.*PETSC_PI*PETSC_PI*t);
205:   val = PetscSinScalar(6*PETSC_PI*z)*ex1 + 3.*PetscSinScalar(2*PETSC_PI*z)*ex2;
206:   return val;
207: }

211: /*
212:    Monitor - User-provided routine to monitor the solution computed at
213:    each timestep.  This example plots the solution and computes the
214:    error in two different norms.

216:    Input Parameters:
217:    ts     - the timestep context
218:    step   - the count of the current step (with 0 meaning the
219:              initial condition)
220:    time   - the current time
221:    u      - the solution at this timestep
222:    ctx    - the user-provided context for this monitoring routine.
223:             In this case we use the application context which contains
224:             information about the problem size, workspace and the exact
225:             solution.
226: */
227: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
228: {
229:   AppCtx         *appctx = (AppCtx*)ctx;
231:   PetscInt       i,m=appctx->m;
232:   PetscReal      norm_2,norm_max,h=1.0/(m+1);
233:   PetscScalar    *u_exact;

235:   /* Compute the exact solution */
236:   VecGetArray(appctx->solution,&u_exact);
237:   for (i=0; i<m; i++) u_exact[i] = exact(appctx->z[i+1],time);
238:   VecRestoreArray(appctx->solution,&u_exact);

240:   /* Print debugging information if desired */
241:   if (appctx->debug) {
242:     PetscPrintf(PETSC_COMM_SELF,"Computed solution vector at time %g\n",(double)time);
243:     VecView(u,PETSC_VIEWER_STDOUT_SELF);
244:     PetscPrintf(PETSC_COMM_SELF,"Exact solution vector\n");
245:     VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
246:   }

248:   /* Compute the 2-norm and max-norm of the error */
249:   VecAXPY(appctx->solution,-1.0,u);
250:   VecNorm(appctx->solution,NORM_2,&norm_2);

252:   norm_2 = PetscSqrtReal(h)*norm_2;
253:   VecNorm(appctx->solution,NORM_MAX,&norm_max);

255:   PetscPrintf(PETSC_COMM_SELF,"Timestep %D: time = %g, 2-norm error = %6.4f, max norm error = %6.4f\n",step,(double)time,(double)norm_2,(double)norm_max);

257:   /*
258:      Print debugging information if desired
259:   */
260:   if (appctx->debug) {
261:     PetscPrintf(PETSC_COMM_SELF,"Error vector\n");
262:     VecView(appctx->solution,PETSC_VIEWER_STDOUT_SELF);
263:   }
264:   return 0;
265: }

267: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
268: %%      Function to solve a linear system using KSP                                           %%
269: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

271: PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
272: {
274:   KSP            ksp;
275:   PC             pc;

277:   /*create the ksp context and set the operators,that is, associate the system matrix with it*/
278:   KSPCreate(PETSC_COMM_WORLD,&ksp);
279:   KSPSetOperators(ksp,obj->Amat,obj->Amat);

281:   /*get the preconditioner context, set its type and the tolerances*/
282:   KSPGetPC(ksp,&pc);
283:   PCSetType(pc,PCLU);
284:   KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);

286:   /*get the command line options if there are any and set them*/
287:   KSPSetFromOptions(ksp);

289:   /*get the linear system (ksp) solve*/
290:   KSPSolve(ksp,obj->ksp_rhs,obj->ksp_sol);

292:   KSPDestroy(&ksp);
293:   return 0;
294: }

296: /***********************************************************************
297:  * Function to return value of basis function or derivative of basis   *
298:  *              function.                                              *
299:  ***********************************************************************
300:  *                                                                     *
301:  *       Arguments:                                                    *
302:  *         x       = array of xpoints or nodal values                  *
303:  *         xx      = point at which the basis function is to be        *
304:  *                     evaluated.                                      *
305:  *         il      = interval containing xx.                           *
306:  *         iq      = indicates which of the two basis functions in     *
307:  *                     interval intrvl should be used                  *
308:  *         nll     = array containing the endpoints of each interval.  *
309:  *         id      = If id ~= 2, the value of the basis function       *
310:  *                     is calculated; if id = 2, the value of the      *
311:  *                     derivative of the basis function is returned.   *
312:  ***********************************************************************/

314: PetscScalar bspl(PetscScalar *x, PetscScalar xx,PetscInt il,PetscInt iq,PetscInt nll[][2],PetscInt id)
315: {
316:   PetscScalar x1,x2,bfcn;
317:   PetscInt    i1,i2,iq1,iq2;

319:   /*** Determine which basis function in interval intrvl is to be used in ***/
320:   iq1 = iq;
321:   if (iq1==0) iq2 = 1;
322:   else iq2 = 0;

324:   /***  Determine endpoint of the interval intrvl ***/
325:   i1=nll[il][iq1];
326:   i2=nll[il][iq2];

328:   /*** Determine nodal values at the endpoints of the interval intrvl ***/
329:   x1=x[i1];
330:   x2=x[i2];
331:   /* printf("x1=%g\tx2=%g\txx=%g\n",x1,x2,xx); */
332:   /*** Evaluate basis function ***/
333:   if (id == 2) bfcn=(1.0)/(x1-x2);
334:   else bfcn=(xx-x2)/(x1-x2);
335:   /* printf("bfcn=%g\n",bfcn); */
336:   return bfcn;
337: }

339: /*---------------------------------------------------------
340:   Function called by rhs function to get B and g
341: ---------------------------------------------------------*/
342: PetscErrorCode femBg(PetscScalar btri[][3],PetscScalar *f,PetscInt nz,PetscScalar *z, PetscReal t)
343: {
344:   PetscInt    i,j,jj,il,ip,ipp,ipq,iq,iquad,iqq;
345:   PetscInt    nli[num_z][2],indx[num_z];
346:   PetscScalar dd,dl,zip,zipq,zz,b_z,bb_z,bij;
347:   PetscScalar zquad[num_z][3],dlen[num_z],qdwt[3];

349:   /*  initializing everything - btri and f are initialized in rhs.c  */
350:   for (i=0; i < nz; i++) {
351:     nli[i][0]   = 0;
352:     nli[i][1]   = 0;
353:     indx[i]     = 0;
354:     zquad[i][0] = 0.0;
355:     zquad[i][1] = 0.0;
356:     zquad[i][2] = 0.0;
357:     dlen[i]     = 0.0;
358:   } /*end for (i)*/

360:   /*  quadrature weights  */
361:   qdwt[0] = 1.0/6.0;
362:   qdwt[1] = 4.0/6.0;
363:   qdwt[2] = 1.0/6.0;

365:   /* 1st and last nodes have Dirichlet boundary condition -
366:      set indices there to -1 */

368:   for (i=0; i < nz-1; i++) indx[i] = i-1;
369:   indx[nz-1] = -1;

371:   ipq = 0;
372:   for (il=0; il < nz-1; il++) {
373:     ip           = ipq;
374:     ipq          = ip+1;
375:     zip          = z[ip];
376:     zipq         = z[ipq];
377:     dl           = zipq-zip;
378:     zquad[il][0] = zip;
379:     zquad[il][1] = (0.5)*(zip+zipq);
380:     zquad[il][2] = zipq;
381:     dlen[il]     = fabs(dl);
382:     nli[il][0]   = ip;
383:     nli[il][1]   = ipq;
384:   }

386:   for (il=0; il < nz-1; il++) {
387:     for (iquad=0; iquad < 3; iquad++) {
388:       dd = (dlen[il])*(qdwt[iquad]);
389:       zz = zquad[il][iquad];

391:       for (iq=0; iq < 2; iq++) {
392:         ip  = nli[il][iq];
393:         b_z = bspl(z,zz,il,iq,nli,2);
394:         i   = indx[ip];

396:         if (i > -1) {
397:           for (iqq=0; iqq < 2; iqq++) {
398:             ipp  = nli[il][iqq];
399:             bb_z = bspl(z,zz,il,iqq,nli,2);
400:             j    = indx[ipp];
401:             bij  = -b_z*bb_z;

403:             if (j > -1) {
404:               jj = 1+j-i;
405:               btri[i][jj] += bij*dd;
406:             } else {
407:               f[i] += bij*dd*exact(z[ipp], t);
408:               /* f[i] += 0.0; */
409:               /* if (il==0 && j==-1) { */
410:               /* f[i] += bij*dd*exact(zz,t); */
411:               /* }*/ /*end if*/
412:             } /*end else*/
413:           } /*end for (iqq)*/
414:         } /*end if (i>0)*/
415:       } /*end for (iq)*/
416:     } /*end for (iquad)*/
417:   } /*end for (il)*/
418:   return 0;
419: }

421: PetscErrorCode femA(AppCtx *obj,PetscInt nz,PetscScalar *z)
422: {
423:   PetscInt       i,j,il,ip,ipp,ipq,iq,iquad,iqq;
424:   PetscInt       nli[num_z][2],indx[num_z];
425:   PetscScalar    dd,dl,zip,zipq,zz,bb,bbb,aij;
426:   PetscScalar    rquad[num_z][3],dlen[num_z],qdwt[3],add_term;

429:   /*  initializing everything  */

431:   for (i=0; i < nz; i++) {
432:     nli[i][0]   = 0;
433:     nli[i][1]   = 0;
434:     indx[i]     = 0;
435:     rquad[i][0] = 0.0;
436:     rquad[i][1] = 0.0;
437:     rquad[i][2] = 0.0;
438:     dlen[i]     = 0.0;
439:   } /*end for (i)*/

441:   /*  quadrature weights  */
442:   qdwt[0] = 1.0/6.0;
443:   qdwt[1] = 4.0/6.0;
444:   qdwt[2] = 1.0/6.0;

446:   /* 1st and last nodes have Dirichlet boundary condition -
447:      set indices there to -1 */

449:   for (i=0; i < nz-1; i++) indx[i]=i-1;
450:   indx[nz-1]=-1;

452:   ipq = 0;

454:   for (il=0; il < nz-1; il++) {
455:     ip           = ipq;
456:     ipq          = ip+1;
457:     zip          = z[ip];
458:     zipq         = z[ipq];
459:     dl           = zipq-zip;
460:     rquad[il][0] = zip;
461:     rquad[il][1] = (0.5)*(zip+zipq);
462:     rquad[il][2] = zipq;
463:     dlen[il]     = fabs(dl);
464:     nli[il][0]   = ip;
465:     nli[il][1]   = ipq;
466:   } /*end for (il)*/

468:   for (il=0; il < nz-1; il++) {
469:     for (iquad=0; iquad < 3; iquad++) {
470:       dd = (dlen[il])*(qdwt[iquad]);
471:       zz = rquad[il][iquad];

473:       for (iq=0; iq < 2; iq++) {
474:         ip = nli[il][iq];
475:         bb = bspl(z,zz,il,iq,nli,1);
476:         i = indx[ip];
477:         if (i > -1) {
478:           for (iqq=0; iqq < 2; iqq++) {
479:             ipp = nli[il][iqq];
480:             bbb = bspl(z,zz,il,iqq,nli,1);
481:             j = indx[ipp];
482:             aij = bb*bbb;
483:             if (j > -1) {
484:               add_term = aij*dd;
485:               MatSetValue(obj->Amat,i,j,add_term,ADD_VALUES);
486:             }/*endif*/
487:           } /*end for (iqq)*/
488:         } /*end if (i>0)*/
489:       } /*end for (iq)*/
490:     } /*end for (iquad)*/
491:   } /*end for (il)*/
492:   MatAssemblyBegin(obj->Amat,MAT_FINAL_ASSEMBLY);
493:   MatAssemblyEnd(obj->Amat,MAT_FINAL_ASSEMBLY);
494:   return 0;
495: }

497: /*---------------------------------------------------------
498:         Function to fill the rhs vector with
499:         By + g values ****
500: ---------------------------------------------------------*/
501: PetscErrorCode rhs(AppCtx *obj,PetscScalar *y, PetscInt nz, PetscScalar *z, PetscReal t)
502: {
503:   PetscInt       i,j,js,je,jj;
504:   PetscScalar    val,g[num_z],btri[num_z][3],add_term;

507:   for (i=0; i < nz-2; i++) {
508:     for (j=0; j <= 2; j++) btri[i][j]=0.0;
509:     g[i] = 0.0;
510:   }

512:   /*  call femBg to set the tri-diagonal b matrix and vector g  */
513:   femBg(btri,g,nz,z,t);

515:   /*  setting the entries of the right hand side vector  */
516:   for (i=0; i < nz-2; i++) {
517:     val = 0.0;
518:     js  = 0;
519:     if (i == 0) js = 1;
520:     je = 2;
521:     if (i == nz-2) je = 1;

523:     for (jj=js; jj <= je; jj++) {
524:       j    = i+jj-1;
525:       val += (btri[i][jj])*(y[j]);
526:     }
527:     add_term = val + g[i];
528:     VecSetValue(obj->ksp_rhs,(PetscInt)i,(PetscScalar)add_term,INSERT_VALUES);
529:   }
530:   VecAssemblyBegin(obj->ksp_rhs);
531:   VecAssemblyEnd(obj->ksp_rhs);

533:   /*  return to main driver function  */
534:   return 0;
535: }

537: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
538: %%   Function to form the right hand side of the time-stepping problem.                       %%
539: %% -------------------------------------------------------------------------------------------%%
540:   if (useAlhs):
541:     globalout = By+g
542:   else if (!useAlhs):
543:     globalout = f(y,t)=Ainv(By+g),
544:       in which the ksp solver to transform the problem A*ydot=By+g
545:       to the problem ydot=f(y,t)=inv(A)*(By+g)
546: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

548: PetscErrorCode RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
549: {
551:   AppCtx         *obj = (AppCtx*)ctx;
552:   PetscScalar    *soln_ptr,soln[num_z];
553:   PetscInt       i,nz=obj->nz;
554:   PetscReal      time;

556:   /* get the previous solution to compute updated system */
557:   VecGetArray(globalin,&soln_ptr);
558:   for (i=0; i < num_z-2; i++) soln[i] = soln_ptr[i];
559:   VecRestoreArray(globalin,&soln_ptr);
560:   soln[num_z-1] = 0.0;
561:   soln[num_z-2] = 0.0;

563:   /* clear out the matrix and rhs for ksp to keep things straight */
564:   VecSet(obj->ksp_rhs,(PetscScalar)0.0);

566:   time = t;
567:   /* get the updated system */
568:   rhs(obj,soln,nz,obj->z,time); /* setup of the By+g rhs */

570:   /* do a ksp solve to get the rhs for the ts problem */
571:   if (obj->useAlhs) {
572:     /* ksp_sol = ksp_rhs */
573:     VecCopy(obj->ksp_rhs,globalout);
574:   } else {
575:     /* ksp_sol = inv(Amat)*ksp_rhs */
576:     Petsc_KSPSolve(obj);
577:     VecCopy(obj->ksp_sol,globalout);
578:   }
579:   return 0;
580: }