Actual source code: ex3.c


  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 10

 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*);

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

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

 61:   /* initializations */
 62:   zInitial  = 0.0;
 63:   zFinal    = 1.0;
 64:   nz        = num_z;
 65:   m         = nz-2;
 66:   appctx.nz = nz;
 67:   max_steps = (PetscInt)10000;

 69:   appctx.m          = m;
 70:   appctx.max_probsz = nz;
 71:   appctx.debug      = PETSC_FALSE;
 72:   appctx.useAlhs    = PETSC_FALSE;

 74:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"","");
 75:   PetscOptionsName("-debug",NULL,NULL,&appctx.debug);
 76:   PetscOptionsName("-useAlhs",NULL,NULL,&appctx.useAlhs);
 77:   PetscOptionsRangeInt("-nphase",NULL,NULL,nphase,&nphase,NULL,1,3);
 78:   PetscOptionsEnd();
 79:   T         = 0.014/nphase;

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

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

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

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

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

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

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

126:   /* set time-step method */
127:   TSSetType(ts,TSCN);

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

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

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

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

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

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

164:   stepsz[0] = 1.0/(2.0*(nz-1)*(nz-1)); /* (mesh_size)^2/2.0 */
165:   ftime     = 0.0;
166:   for (k=0; k<nphase; k++) {
167:     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));}
168:     TSSetTime(ts,ftime);
169:     TSSetTimeStep(ts,stepsz[k]);
170:     TSSetMaxSteps(ts,max_steps);
171:     TSSetMaxTime(ts,(k+1)*T);
172:     TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

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

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

192:   PetscFinalize();
193:   return ierr;
194: }

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

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

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

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

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

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

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

251:   norm_2 = PetscSqrtReal(h)*norm_2;
252:   VecNorm(appctx->solution,NORM_MAX,&norm_max);
253:   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);

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

265: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
266:       Function to solve a linear system using KSP
267: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

269: PetscErrorCode Petsc_KSPSolve(AppCtx *obj)
270: {
272:   KSP            ksp;
273:   PC             pc;

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

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

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

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

290:   KSPDestroy(&ksp);
291:   return 0;
292: }

294: /***********************************************************************
295:   Function to return value of basis function or derivative of basis function.
296:  ***********************************************************************

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

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

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

321:   /*    Determine endpoint of the interval intrvl   */
322:   i1=nll[il][iq1];
323:   i2=nll[il][iq2];

325:   /*   Determine nodal values at the endpoints of the interval intrvl   */
326:   x1=x[i1];
327:   x2=x[i2];

329:   /*   Evaluate basis function   */
330:   if (id == 2) bfcn=(1.0)/(x1-x2);
331:   else bfcn=(xx-x2)/(x1-x2);
332:   return bfcn;
333: }

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

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

356:   /*  quadrature weights  */
357:   qdwt[0] = 1.0/6.0;
358:   qdwt[1] = 4.0/6.0;
359:   qdwt[2] = 1.0/6.0;

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

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

367:   ipq = 0;
368:   for (il=0; il < nz-1; il++) {
369:     ip           = ipq;
370:     ipq          = ip+1;
371:     zip          = z[ip];
372:     zipq         = z[ipq];
373:     dl           = zipq-zip;
374:     zquad[il][0] = zip;
375:     zquad[il][1] = (0.5)*(zip+zipq);
376:     zquad[il][2] = zipq;
377:     dlen[il]     = PetscAbsScalar(dl);
378:     nli[il][0]   = ip;
379:     nli[il][1]   = ipq;
380:   }

382:   for (il=0; il < nz-1; il++) {
383:     for (iquad=0; iquad < 3; iquad++) {
384:       dd = (dlen[il])*(qdwt[iquad]);
385:       zz = zquad[il][iquad];

387:       for (iq=0; iq < 2; iq++) {
388:         ip  = nli[il][iq];
389:         b_z = bspl(z,zz,il,iq,nli,2);
390:         i   = indx[ip];

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

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

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

425:   /*  initializing everything  */
426:   for (i=0; i < nz; i++) {
427:     nli[i][0]   = 0;
428:     nli[i][1]   = 0;
429:     indx[i]     = 0;
430:     rquad[i][0] = 0.0;
431:     rquad[i][1] = 0.0;
432:     rquad[i][2] = 0.0;
433:     dlen[i]     = 0.0;
434:   } /*end for (i)*/

436:   /*  quadrature weights  */
437:   qdwt[0] = 1.0/6.0;
438:   qdwt[1] = 4.0/6.0;
439:   qdwt[2] = 1.0/6.0;

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

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

447:   ipq = 0;

449:   for (il=0; il < nz-1; il++) {
450:     ip           = ipq;
451:     ipq          = ip+1;
452:     zip          = z[ip];
453:     zipq         = z[ipq];
454:     dl           = zipq-zip;
455:     rquad[il][0] = zip;
456:     rquad[il][1] = (0.5)*(zip+zipq);
457:     rquad[il][2] = zipq;
458:     dlen[il]     = PetscAbsScalar(dl);
459:     nli[il][0]   = ip;
460:     nli[il][1]   = ipq;
461:   } /*end for (il)*/

463:   for (il=0; il < nz-1; il++) {
464:     for (iquad=0; iquad < 3; iquad++) {
465:       dd = (dlen[il])*(qdwt[iquad]);
466:       zz = rquad[il][iquad];

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

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

502:   for (i=0; i < nz-2; i++) {
503:     for (j=0; j <= 2; j++) btri[i][j]=0.0;
504:     g[i] = 0.0;
505:   }

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

510:   /*  setting the entries of the right hand side vector  */
511:   for (i=0; i < nz-2; i++) {
512:     val = 0.0;
513:     js  = 0;
514:     if (i == 0) js = 1;
515:     je = 2;
516:     if (i == nz-2) je = 1;

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

530: /*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
531: %%   Function to form the right hand side of the time-stepping problem.                       %%
532: %% -------------------------------------------------------------------------------------------%%
533:   if (useAlhs):
534:     globalout = By+g
535:   else if (!useAlhs):
536:     globalout = f(y,t)=Ainv(By+g),
537:       in which the ksp solver to transform the problem A*ydot=By+g
538:       to the problem ydot=f(y,t)=inv(A)*(By+g)
539: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/

541: PetscErrorCode RHSfunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
542: {
543:   PetscErrorCode    ierr;
544:   AppCtx            *obj = (AppCtx*)ctx;
545:   PetscScalar       soln[num_z];
546:   const PetscScalar *soln_ptr;
547:   PetscInt          i,nz=obj->nz;
548:   PetscReal         time;

550:   /* get the previous solution to compute updated system */
551:   VecGetArrayRead(globalin,&soln_ptr);
552:   for (i=0; i < num_z-2; i++) soln[i] = soln_ptr[i];
553:   VecRestoreArrayRead(globalin,&soln_ptr);
554:   soln[num_z-1] = 0.0;
555:   soln[num_z-2] = 0.0;

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

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

564:   /* do a ksp solve to get the rhs for the ts problem */
565:   if (obj->useAlhs) {
566:     /* ksp_sol = ksp_rhs */
567:     VecCopy(obj->ksp_rhs,globalout);
568:   } else {
569:     /* ksp_sol = inv(Amat)*ksp_rhs */
570:     Petsc_KSPSolve(obj);
571:     VecCopy(obj->ksp_sol,globalout);
572:   }
573:   return 0;
574: }

576: /*TEST

578:     build:
579:       requires: !complex

581:     test:
582:       suffix: euler
583:       output_file: output/ex3.out

585:     test:
586:       suffix: 2
587:       args:   -useAlhs
588:       output_file: output/ex3.out
589:       TODO: Broken because SNESComputeJacobianDefault is incompatible with TSComputeIJacobianConstant

591: TEST*/