Actual source code: sundials.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  1: /*
  2:     Provides a PETSc interface to SUNDIALS/CVODE solver.
  3:     The interface to PVODE (old version of CVODE) was originally contributed
  4:     by Liyang Xu. It has been redone by Hong Zhang and Dinesh Kaushik.

  6:     Reference: sundials-2.4.0/examples/cvode/parallel/cvDiurnal_kry_p.c
  7: */
  8: #include <../src/ts/impls/implicit/sundials/sundials.h>  /*I "petscts.h" I*/

 10: /*
 11:       TSPrecond_Sundials - function that we provide to SUNDIALS to
 12:                         evaluate the preconditioner.
 13: */
 16: PetscErrorCode TSPrecond_Sundials(realtype tn,N_Vector y,N_Vector fy,booleantype jok,booleantype *jcurPtr,
 17:                                   realtype _gamma,void *P_data,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
 18: {
 19:   TS             ts     = (TS) P_data;
 20:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
 21:   PC             pc;
 23:   Mat            J,P;
 24:   Vec            yy  = cvode->w1,yydot = cvode->ydot;
 25:   PetscReal      gm  = (PetscReal)_gamma;
 26:   PetscScalar    *y_data;

 29:   TSGetIJacobian(ts,&J,&P,NULL,NULL);
 30:   y_data = (PetscScalar*) N_VGetArrayPointer(y);
 31:   VecPlaceArray(yy,y_data);
 32:   VecZeroEntries(yydot); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */
 33:   /* compute the shifted Jacobian   (1/gm)*I + Jrest */
 34:   TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,J,P,PETSC_FALSE);
 35:   VecResetArray(yy);
 36:   MatScale(P,gm); /* turn into I-gm*Jrest, J is not used by Sundials  */
 37:   *jcurPtr = TRUE;
 38:   TSSundialsGetPC(ts,&pc);
 39:   PCSetOperators(pc,J,P);
 40:   return(0);
 41: }

 43: /*
 44:      TSPSolve_Sundials -  routine that we provide to Sundials that applies the preconditioner.
 45: */
 48: PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,
 49:                                  realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
 50: {
 51:   TS             ts     = (TS) P_data;
 52:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
 53:   PC             pc;
 54:   Vec            rr = cvode->w1,zz = cvode->w2;
 56:   PetscScalar    *r_data,*z_data;

 59:   /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
 60:   r_data = (PetscScalar*) N_VGetArrayPointer(r);
 61:   z_data = (PetscScalar*) N_VGetArrayPointer(z);
 62:   VecPlaceArray(rr,r_data);
 63:   VecPlaceArray(zz,z_data);

 65:   /* Solve the Px=r and put the result in zz */
 66:   TSSundialsGetPC(ts,&pc);
 67:   PCApply(pc,rr,zz);
 68:   VecResetArray(rr);
 69:   VecResetArray(zz);
 70:   return(0);
 71: }

 73: /*
 74:         TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
 75: */
 78: int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
 79: {
 80:   TS             ts = (TS) ctx;
 81:   DM             dm;
 82:   DMTS           tsdm;
 83:   TSIFunction    ifunction;
 84:   MPI_Comm       comm;
 85:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
 86:   Vec            yy     = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;
 87:   PetscScalar    *y_data,*ydot_data;

 91:   PetscObjectGetComm((PetscObject)ts,&comm);
 92:   /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
 93:   y_data    = (PetscScalar*) N_VGetArrayPointer(y);
 94:   ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);
 95:   VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
 96:   VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);

 98:   /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
 99:   TSGetDM(ts,&dm);
100:   DMGetDMTS(dm,&tsdm);
101:   DMTSGetIFunction(dm,&ifunction,NULL);
102:   if (!ifunction) {
103:     TSComputeRHSFunction(ts,t,yy,yyd);
104:   } else {                      /* If rhsfunction is also set, this computes both parts and shifts them to the right */
105:     VecZeroEntries(yydot);
106:     TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr);
107:     VecScale(yyd,-1.);
108:   }
109:   VecResetArray(yy);CHKERRABORT(comm,ierr);
110:   VecResetArray(yyd);CHKERRABORT(comm,ierr);
111:   return(0);
112: }

114: /*
115:        TSStep_Sundials - Calls Sundials to integrate the ODE.
116: */
119: PetscErrorCode TSStep_Sundials(TS ts)
120: {
121:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
123:   PetscInt       flag;
124:   long int       its,nsteps;
125:   realtype       t,tout;
126:   PetscScalar    *y_data;
127:   void           *mem;

130:   mem  = cvode->mem;
131:   tout = ts->max_time;
132:   VecGetArray(ts->vec_sol,&y_data);
133:   N_VSetArrayPointer((realtype*)y_data,cvode->y);
134:   VecRestoreArray(ts->vec_sol,NULL);

136:   TSPreStep(ts);

138:   /* We would like to call TSPreStep() when starting each step (including rejections), TSPreStage(),
139:    * and TSPostStage() before each stage solve, but CVode does not appear to support this. */
140:   if (cvode->monitorstep) flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
141:   else flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);

143:   if (flag) { /* display error message */
144:     switch (flag) {
145:       case CV_ILL_INPUT:
146:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");
147:         break;
148:       case CV_TOO_CLOSE:
149:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");
150:         break;
151:       case CV_TOO_MUCH_WORK: {
152:         PetscReal      tcur;
153:         CVodeGetNumSteps(mem,&nsteps);
154:         CVodeGetCurrentTime(mem,&tcur);
155:         SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_WORK. At t=%g, nsteps %D exceeds mxstep %D. Increase '-ts_max_steps <>' or modify TSSetDuration()",(double)tcur,nsteps,ts->max_steps);
156:       } break;
157:       case CV_TOO_MUCH_ACC:
158:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");
159:         break;
160:       case CV_ERR_FAILURE:
161:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");
162:         break;
163:       case CV_CONV_FAILURE:
164:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");
165:         break;
166:       case CV_LINIT_FAIL:
167:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");
168:         break;
169:       case CV_LSETUP_FAIL:
170:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");
171:         break;
172:       case CV_LSOLVE_FAIL:
173:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");
174:         break;
175:       case CV_RHSFUNC_FAIL:
176:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");
177:         break;
178:       case CV_FIRST_RHSFUNC_ERR:
179:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");
180:         break;
181:       case CV_REPTD_RHSFUNC_ERR:
182:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");
183:         break;
184:       case CV_UNREC_RHSFUNC_ERR:
185:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");
186:         break;
187:       case CV_RTFUNC_FAIL:
188:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");
189:         break;
190:       default:
191:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);
192:     }
193:   }

195:   /* copy the solution from cvode->y to cvode->update and sol */
196:   VecPlaceArray(cvode->w1,y_data);
197:   VecCopy(cvode->w1,cvode->update);
198:   VecResetArray(cvode->w1);
199:   VecCopy(cvode->update,ts->vec_sol);
200:   CVodeGetNumNonlinSolvIters(mem,&its);
201:   CVSpilsGetNumLinIters(mem, &its);
202:   ts->snes_its = its; ts->ksp_its = its;

204:   ts->time_step = t - ts->ptime;
205:   ts->ptime     = t;
206:   ts->steps++;

208:   CVodeGetNumSteps(mem,&nsteps);
209:   if (!cvode->monitorstep) ts->steps = nsteps;
210:   return(0);
211: }

215: static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
216: {
217:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
218:   N_Vector       y;
220:   PetscScalar    *x_data;
221:   PetscInt       glosize,locsize;

224:   /* get the vector size */
225:   VecGetSize(X,&glosize);
226:   VecGetLocalSize(X,&locsize);

228:   /* allocate the memory for N_Vec y */
229:   y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
230:   if (!y) SETERRQ(PETSC_COMM_SELF,1,"Interpolated y is not allocated");

232:   VecGetArray(X,&x_data);
233:   N_VSetArrayPointer((realtype*)x_data,y);
234:   CVodeGetDky(cvode->mem,t,0,y);
235:   VecRestoreArray(X,&x_data);
236:   return(0);
237: }

241: PetscErrorCode TSReset_Sundials(TS ts)
242: {
243:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;

247:   VecDestroy(&cvode->update);
248:   VecDestroy(&cvode->ydot);
249:   VecDestroy(&cvode->w1);
250:   VecDestroy(&cvode->w2);
251:   if (cvode->mem) CVodeFree(&cvode->mem);
252:   return(0);
253: }

257: PetscErrorCode TSDestroy_Sundials(TS ts)
258: {
259:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;

263:   TSReset_Sundials(ts);
264:   MPI_Comm_free(&(cvode->comm_sundials));
265:   PetscFree(ts->data);
266:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);
267:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);
268:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);
269:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);
270:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);
271:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);
272:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);
273:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);
274:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);
275:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",NULL);
276:   return(0);
277: }

281: PetscErrorCode TSSetUp_Sundials(TS ts)
282: {
283:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
285:   PetscInt       glosize,locsize,i,flag;
286:   PetscScalar    *y_data,*parray;
287:   void           *mem;
288:   PC             pc;
289:   PCType         pctype;
290:   PetscBool      pcnone;

293:   /* get the vector size */
294:   VecGetSize(ts->vec_sol,&glosize);
295:   VecGetLocalSize(ts->vec_sol,&locsize);

297:   /* allocate the memory for N_Vec y */
298:   cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
299:   if (!cvode->y) SETERRQ(PETSC_COMM_SELF,1,"cvode->y is not allocated");

301:   /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
302:   VecGetArray(ts->vec_sol,&parray);
303:   y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);
304:   for (i = 0; i < locsize; i++) y_data[i] = parray[i];
305:   VecRestoreArray(ts->vec_sol,NULL);

307:   VecDuplicate(ts->vec_sol,&cvode->update);
308:   VecDuplicate(ts->vec_sol,&cvode->ydot);
309:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);
310:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);

312:   /*
313:     Create work vectors for the TSPSolve_Sundials() routine. Note these are
314:     allocated with zero space arrays because the actual array space is provided
315:     by Sundials and set using VecPlaceArray().
316:   */
317:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w1);
318:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w2);
319:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);
320:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);

322:   /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */
323:   mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
324:   if (!mem) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"CVodeCreate() fails");
325:   cvode->mem = mem;

327:   /* Set the pointer to user-defined data */
328:   flag = CVodeSetUserData(mem, ts);
329:   if (flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetUserData() fails");

331:   /* Sundials may choose to use a smaller initial step, but will never use a larger step. */
332:   flag = CVodeSetInitStep(mem,(realtype)ts->time_step);
333:   if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetInitStep() failed");
334:   if (cvode->mindt > 0) {
335:     flag = CVodeSetMinStep(mem,(realtype)cvode->mindt);
336:     if (flag) {
337:       if (flag == CV_MEM_NULL) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL");
338:       else if (flag == CV_ILL_INPUT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, hmin is nonpositive or it exceeds the maximum allowable step size");
339:       else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed");
340:     }
341:   }
342:   if (cvode->maxdt > 0) {
343:     flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);
344:     if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMaxStep() failed");
345:   }

347:   /* Call CVodeInit to initialize the integrator memory and specify the
348:    * user's right hand side function in u'=f(t,u), the inital time T0, and
349:    * the initial dependent variable vector cvode->y */
350:   flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y);
351:   if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag);

353:   /* specifies scalar relative and absolute tolerances */
354:   flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol);
355:   if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag);

357:   /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
358:   flag = CVodeSetMaxNumSteps(mem,ts->max_steps);

360:   /* call CVSpgmr to use GMRES as the linear solver.        */
361:   /* setup the ode integrator with the given preconditioner */
362:   TSSundialsGetPC(ts,&pc);
363:   PCGetType(pc,&pctype);
364:   PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);
365:   if (pcnone) {
366:     flag = CVSpgmr(mem,PREC_NONE,0);
367:     if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
368:   } else {
369:     flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl);
370:     if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);

372:     /* Set preconditioner and solve routines Precond and PSolve,
373:      and the pointer to the user-defined block data */
374:     flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);
375:     if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag);
376:   }

378:   flag = CVSpilsSetGSType(mem, MODIFIED_GS);
379:   if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);
380:   return(0);
381: }

383: /* type of CVODE linear multistep method */
384: const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",0};
385: /* type of G-S orthogonalization used by CVODE linear solver */
386: const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",0};

390: PetscErrorCode TSSetFromOptions_Sundials(TS ts)
391: {
392:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
394:   int            indx;
395:   PetscBool      flag;
396:   PC             pc;

399:   PetscOptionsHead("SUNDIALS ODE solver options");
400:   PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);
401:   if (flag) {
402:     TSSundialsSetType(ts,(TSSundialsLmmType)indx);
403:   }
404:   PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);
405:   if (flag) {
406:     TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
407:   }
408:   PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);
409:   PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);
410:   PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);
411:   PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);
412:   PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,&flag);
413:   PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,&flag);
414:   PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internel steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);
415:   PetscOptionsTail();
416:   TSSundialsGetPC(ts,&pc);
417:   PCSetFromOptions(pc);
418:   return(0);
419: }

423: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
424: {
425:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
427:   char           *type;
428:   char           atype[] = "Adams";
429:   char           btype[] = "BDF: backward differentiation formula";
430:   PetscBool      iascii,isstring;
431:   long int       nsteps,its,nfevals,nlinsetups,nfails,itmp;
432:   PetscInt       qlast,qcur;
433:   PetscReal      hinused,hlast,hcur,tcur,tolsfac;
434:   PC             pc;

437:   if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
438:   else                                     type = btype;

440:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
441:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
442:   if (iascii) {
443:     PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");
444:     PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);
445:     PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);
446:     PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",cvode->linear_tol);
447:     PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);
448:     if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
449:       PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
450:     } else {
451:       PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
452:     }
453:     if (cvode->mindt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",cvode->mindt);}
454:     if (cvode->maxdt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",cvode->maxdt);}

456:     /* Outputs from CVODE, CVSPILS */
457:     CVodeGetTolScaleFactor(cvode->mem,&tolsfac);
458:     PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);
459:     CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,
460:                                    &nlinsetups,&nfails,&qlast,&qcur,
461:                                    &hinused,&hlast,&hcur,&tcur);
462:     PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);
463:     PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);
464:     PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);
465:     PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);

467:     CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);
468:     PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);
469:     PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);

471:     CVSpilsGetNumLinIters(cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
472:     PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);
473:     CVSpilsGetNumConvFails(cvode->mem,&itmp);
474:     PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);

476:     TSSundialsGetPC(ts,&pc);
477:     PCView(pc,viewer);
478:     CVSpilsGetNumPrecEvals(cvode->mem,&itmp);
479:     PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);
480:     CVSpilsGetNumPrecSolves(cvode->mem,&itmp);
481:     PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);

483:     CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);
484:     PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);
485:     CVSpilsGetNumRhsEvals(cvode->mem,&itmp);
486:     PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);
487:   } else if (isstring) {
488:     PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
489:   }
490:   return(0);
491: }


494: /* --------------------------------------------------------------------------*/
497: PetscErrorCode  TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
498: {
499:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

502:   cvode->cvode_type = type;
503:   return(0);
504: }

508: PetscErrorCode  TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
509: {
510:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

513:   cvode->maxl = maxl;
514:   return(0);
515: }

519: PetscErrorCode  TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
520: {
521:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

524:   cvode->linear_tol = tol;
525:   return(0);
526: }

530: PetscErrorCode  TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
531: {
532:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

535:   cvode->gtype = type;
536:   return(0);
537: }

541: PetscErrorCode  TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
542: {
543:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

546:   if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
547:   if (rel != PETSC_DECIDE)  cvode->reltol = rel;
548:   return(0);
549: }

553: PetscErrorCode  TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
554: {
555:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

558:   cvode->mindt = mindt;
559:   return(0);
560: }

564: PetscErrorCode  TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
565: {
566:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

569:   cvode->maxdt = maxdt;
570:   return(0);
571: }
574: PetscErrorCode  TSSundialsGetPC_Sundials(TS ts,PC *pc)
575: {
576:   SNES           snes;
577:   KSP            ksp;

581:   TSGetSNES(ts,&snes);
582:   SNESGetKSP(snes,&ksp);
583:   KSPGetPC(ksp,pc);
584:   return(0);
585: }

589: PetscErrorCode  TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
590: {
592:   if (nonlin) *nonlin = ts->snes_its;
593:   if (lin)    *lin    = ts->ksp_its;
594:   return(0);
595: }

599: PetscErrorCode  TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
600: {
601:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

604:   cvode->monitorstep = s;
605:   return(0);
606: }
607: /* -------------------------------------------------------------------------------------------*/

611: /*@C
612:    TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.

614:    Not Collective

616:    Input parameters:
617: .    ts     - the time-step context

619:    Output Parameters:
620: +   nonlin - number of nonlinear iterations
621: -   lin    - number of linear iterations

623:    Level: advanced

625:    Notes:
626:     These return the number since the creation of the TS object

628: .keywords: non-linear iterations, linear iterations

630: .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
631:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
632:           TSSundialsGetIterations(), TSSundialsSetType(),
633:           TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()

635: @*/
636: PetscErrorCode  TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
637: {

641:   PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));
642:   return(0);
643: }

647: /*@
648:    TSSundialsSetType - Sets the method that Sundials will use for integration.

650:    Logically Collective on TS

652:    Input parameters:
653: +    ts     - the time-step context
654: -    type   - one of  SUNDIALS_ADAMS or SUNDIALS_BDF

656:    Level: intermediate

658: .keywords: Adams, backward differentiation formula

660: .seealso: TSSundialsGetIterations(),  TSSundialsSetMaxl(),
661:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
662:           TSSundialsGetIterations(), TSSundialsSetType(),
663:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
664:           TSSetExactFinalTime()
665: @*/
666: PetscErrorCode  TSSundialsSetType(TS ts,TSSundialsLmmType type)
667: {

671:   PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));
672:   return(0);
673: }

677: /*@
678:    TSSundialsSetMaxl - Sets the dimension of the Krylov space used by
679:        GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so
680:        this is the maximum number of GMRES steps that will be used.

682:    Logically Collective on TS

684:    Input parameters:
685: +    ts      - the time-step context
686: -    maxl - number of direction vectors (the dimension of Krylov subspace).

688:    Level: advanced

690: .keywords: GMRES

692: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
693:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
694:           TSSundialsGetIterations(), TSSundialsSetType(),
695:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
696:           TSSetExactFinalTime()

698: @*/
699: PetscErrorCode  TSSundialsSetMaxl(TS ts,PetscInt maxl)
700: {

705:   PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));
706:   return(0);
707: }

711: /*@
712:    TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
713:        system by SUNDIALS.

715:    Logically Collective on TS

717:    Input parameters:
718: +    ts     - the time-step context
719: -    tol    - the factor by which the tolerance on the nonlinear solver is
720:              multiplied to get the tolerance on the linear solver, .05 by default.

722:    Level: advanced

724: .keywords: GMRES, linear convergence tolerance, SUNDIALS

726: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
727:           TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
728:           TSSundialsGetIterations(), TSSundialsSetType(),
729:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
730:           TSSetExactFinalTime()

732: @*/
733: PetscErrorCode  TSSundialsSetLinearTolerance(TS ts,double tol)
734: {

739:   PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));
740:   return(0);
741: }

745: /*@
746:    TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
747:         in GMRES method by SUNDIALS linear solver.

749:    Logically Collective on TS

751:    Input parameters:
752: +    ts  - the time-step context
753: -    type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS

755:    Level: advanced

757: .keywords: Sundials, orthogonalization

759: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
760:           TSSundialsSetLinearTolerance(),  TSSundialsSetTolerance(),
761:           TSSundialsGetIterations(), TSSundialsSetType(),
762:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
763:           TSSetExactFinalTime()

765: @*/
766: PetscErrorCode  TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
767: {

771:   PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));
772:   return(0);
773: }

777: /*@
778:    TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
779:                          Sundials for error control.

781:    Logically Collective on TS

783:    Input parameters:
784: +    ts  - the time-step context
785: .    aabs - the absolute tolerance
786: -    rel - the relative tolerance

788:      See the Cvode/Sundials users manual for exact details on these parameters. Essentially
789:     these regulate the size of the error for a SINGLE timestep.

791:    Level: intermediate

793: .keywords: Sundials, tolerance

795: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
796:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
797:           TSSundialsGetIterations(), TSSundialsSetType(),
798:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
799:           TSSetExactFinalTime()

801: @*/
802: PetscErrorCode  TSSundialsSetTolerance(TS ts,double aabs,double rel)
803: {

807:   PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));
808:   return(0);
809: }

813: /*@
814:    TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.

816:    Input Parameter:
817: .    ts - the time-step context

819:    Output Parameter:
820: .    pc - the preconditioner context

822:    Level: advanced

824: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
825:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
826:           TSSundialsGetIterations(), TSSundialsSetType(),
827:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
828: @*/
829: PetscErrorCode  TSSundialsGetPC(TS ts,PC *pc)
830: {

834:   PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));
835:   return(0);
836: }

840: /*@
841:    TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.

843:    Input Parameter:
844: +   ts - the time-step context
845: -   mindt - lowest time step if positive, negative to deactivate

847:    Note:
848:    Sundials will error if it is not possible to keep the estimated truncation error below
849:    the tolerance set with TSSundialsSetTolerance() without going below this step size.

851:    Level: beginner

853: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
854: @*/
855: PetscErrorCode  TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
856: {

860:   PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));
861:   return(0);
862: }

866: /*@
867:    TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.

869:    Input Parameter:
870: +   ts - the time-step context
871: -   maxdt - lowest time step if positive, negative to deactivate

873:    Level: beginner

875: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
876: @*/
877: PetscErrorCode  TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
878: {

882:   PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
883:   return(0);
884: }

888: /*@
889:    TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false).

891:    Input Parameter:
892: +   ts - the time-step context
893: -   ft - PETSC_TRUE if monitor, else PETSC_FALSE

895:    Level: beginner

897: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
898:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
899:           TSSundialsGetIterations(), TSSundialsSetType(),
900:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
901: @*/
902: PetscErrorCode  TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
903: {

907:   PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));
908:   return(0);
909: }
910: /* -------------------------------------------------------------------------------------------*/
911: /*MC
912:       TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)

914:    Options Database:
915: +    -ts_sundials_type <bdf,adams>
916: .    -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
917: .    -ts_sundials_atol <tol> - Absolute tolerance for convergence
918: .    -ts_sundials_rtol <tol> - Relative tolerance for convergence
919: .    -ts_sundials_linear_tolerance <tol>
920: .    -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
921: -    -ts_sundials_monitor_steps - Monitor SUNDIALS internel steps


924:     Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply
925:            only PETSc PC options

927:     Level: beginner

929: .seealso:  TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(),
930:            TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime()

932: M*/
935: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
936: {
937:   TS_Sundials    *cvode;
939:   PC             pc;

942:   ts->ops->reset          = TSReset_Sundials;
943:   ts->ops->destroy        = TSDestroy_Sundials;
944:   ts->ops->view           = TSView_Sundials;
945:   ts->ops->setup          = TSSetUp_Sundials;
946:   ts->ops->step           = TSStep_Sundials;
947:   ts->ops->interpolate    = TSInterpolate_Sundials;
948:   ts->ops->setfromoptions = TSSetFromOptions_Sundials;

950:   PetscNewLog(ts,&cvode);

952:   ts->data           = (void*)cvode;
953:   cvode->cvode_type  = SUNDIALS_BDF;
954:   cvode->gtype       = SUNDIALS_CLASSICAL_GS;
955:   cvode->maxl        = 5;
956:   cvode->linear_tol  = .05;
957:   cvode->monitorstep = PETSC_TRUE;

959:   MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));

961:   cvode->mindt = -1.;
962:   cvode->maxdt = -1.;

964:   /* set tolerance for Sundials */
965:   cvode->reltol = 1e-6;
966:   cvode->abstol = 1e-6;

968:   /* set PCNONE as default pctype */
969:   TSSundialsGetPC_Sundials(ts,&pc);
970:   PCSetType(pc,PCNONE);

972:   if (ts->exact_final_time != TS_EXACTFINALTIME_STEPOVER) ts->exact_final_time = TS_EXACTFINALTIME_STEPOVER;

974:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);
975:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);
976:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);
977:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);
978:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);
979:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);
980:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);
981:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);
982:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);
983:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);
984:   return(0);
985: }