Actual source code: sundials.c

  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>

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

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

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

 55:   /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
 56:   r_data = (PetscScalar*) N_VGetArrayPointer(r);
 57:   z_data = (PetscScalar*) N_VGetArrayPointer(z);
 58:   VecPlaceArray(rr,r_data);
 59:   VecPlaceArray(zz,z_data);

 61:   /* Solve the Px=r and put the result in zz */
 62:   TSSundialsGetPC(ts,&pc);
 63:   PCApply(pc,rr,zz);
 64:   VecResetArray(rr);
 65:   VecResetArray(zz);
 66:   return(0);
 67: }

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

 85:   PetscObjectGetComm((PetscObject)ts,&comm);
 86:   /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
 87:   y_data    = (PetscScalar*) N_VGetArrayPointer(y);
 88:   ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);
 89:   VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
 90:   VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);

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

108: /*
109:        TSStep_Sundials - Calls Sundials to integrate the ODE.
110: */
111: PetscErrorCode TSStep_Sundials(TS ts)
112: {
113:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
115:   PetscInt       flag;
116:   long int       nits,lits,nsteps;
117:   realtype       t,tout;
118:   PetscScalar    *y_data;
119:   void           *mem;

122:   mem  = cvode->mem;
123:   tout = ts->max_time;
124:   VecGetArray(ts->vec_sol,&y_data);
125:   N_VSetArrayPointer((realtype*)y_data,cvode->y);
126:   VecRestoreArray(ts->vec_sol,NULL);

128:   /* We would like to TSPreStage() and TSPostStage()
129:    * before each stage solve but CVode does not appear to support this. */
130:   if (cvode->monitorstep)
131:     flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
132:   else
133:     flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);

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

187:   /* log inner nonlinear and linear iterations */
188:   CVodeGetNumNonlinSolvIters(mem,&nits);
189:   CVSpilsGetNumLinIters(mem,&lits);
190:   ts->snes_its += nits; ts->ksp_its = lits;

192:   /* copy the solution from cvode->y to cvode->update and sol */
193:   VecPlaceArray(cvode->w1,y_data);
194:   VecCopy(cvode->w1,cvode->update);
195:   VecResetArray(cvode->w1);
196:   VecCopy(cvode->update,ts->vec_sol);

198:   ts->time_step = t - ts->ptime;
199:   ts->ptime = t;

201:   CVodeGetNumSteps(mem,&nsteps);
202:   if (!cvode->monitorstep) ts->steps += nsteps - 1; /* TSStep() increments the step counter by one */
203:   return(0);
204: }

206: static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
207: {
208:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
209:   N_Vector       y;
211:   PetscScalar    *x_data;
212:   PetscInt       glosize,locsize;

215:   /* get the vector size */
216:   VecGetSize(X,&glosize);
217:   VecGetLocalSize(X,&locsize);
218:   VecGetArray(X,&x_data);

220:   /* Initialize N_Vec y with x_data */
221:   y = N_VMake_Parallel(cvode->comm_sundials,locsize,glosize,(realtype*)x_data);
222:   if (!y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Interpolated y is not allocated");

224:   CVodeGetDky(cvode->mem,t,0,y);
225:   VecRestoreArray(X,&x_data);
226:   return(0);
227: }

229: PetscErrorCode TSReset_Sundials(TS ts)
230: {
231:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;

235:   VecDestroy(&cvode->update);
236:   VecDestroy(&cvode->ydot);
237:   VecDestroy(&cvode->w1);
238:   VecDestroy(&cvode->w2);
239:   if (cvode->mem) CVodeFree(&cvode->mem);
240:   return(0);
241: }

243: PetscErrorCode TSDestroy_Sundials(TS ts)
244: {
245:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;

249:   TSReset_Sundials(ts);
250:   MPI_Comm_free(&(cvode->comm_sundials));
251:   PetscFree(ts->data);
252:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);
253:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);
254:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);
255:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);
256:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);
257:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);
258:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);
259:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);
260:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);
261:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",NULL);
262:   return(0);
263: }

265: PetscErrorCode TSSetUp_Sundials(TS ts)
266: {
267:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
269:   PetscInt       glosize,locsize,i,flag;
270:   PetscScalar    *y_data,*parray;
271:   void           *mem;
272:   PC             pc;
273:   PCType         pctype;
274:   PetscBool      pcnone;

277:   if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for exact final time option 'MATCHSTEP' when using Sundials");

279:   /* get the vector size */
280:   VecGetSize(ts->vec_sol,&glosize);
281:   VecGetLocalSize(ts->vec_sol,&locsize);

283:   /* allocate the memory for N_Vec y */
284:   cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
285:   if (!cvode->y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"cvode->y is not allocated");

287:   /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
288:   VecGetArray(ts->vec_sol,&parray);
289:   y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);
290:   for (i = 0; i < locsize; i++) y_data[i] = parray[i];
291:   VecRestoreArray(ts->vec_sol,NULL);

293:   VecDuplicate(ts->vec_sol,&cvode->update);
294:   VecDuplicate(ts->vec_sol,&cvode->ydot);
295:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);
296:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);

298:   /*
299:     Create work vectors for the TSPSolve_Sundials() routine. Note these are
300:     allocated with zero space arrays because the actual array space is provided
301:     by Sundials and set using VecPlaceArray().
302:   */
303:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,NULL,&cvode->w1);
304:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,NULL,&cvode->w2);
305:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);
306:   PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);

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

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

317:   /* Sundials may choose to use a smaller initial step, but will never use a larger step. */
318:   flag = CVodeSetInitStep(mem,(realtype)ts->time_step);
319:   if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetInitStep() failed");
320:   if (cvode->mindt > 0) {
321:     flag = CVodeSetMinStep(mem,(realtype)cvode->mindt);
322:     if (flag) {
323:       if (flag == CV_MEM_NULL) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL");
324:       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");
325:       else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed");
326:     }
327:   }
328:   if (cvode->maxdt > 0) {
329:     flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);
330:     if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMaxStep() failed");
331:   }

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

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

343:   /* Specify max order of BDF / ADAMS method */
344:   if (cvode->maxord != PETSC_DEFAULT) {
345:     flag = CVodeSetMaxOrd(mem,cvode->maxord);
346:     if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetMaxOrd() fails, flag %d",flag);
347:   }

349:   /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
350:   flag = CVodeSetMaxNumSteps(mem,ts->max_steps);
351:   if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetMaxNumSteps() fails, flag %d",flag);

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

365:     /* Set preconditioner and solve routines Precond and PSolve,
366:      and the pointer to the user-defined block data */
367:     flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);
368:     if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag);
369:   }

371:   flag = CVSpilsSetGSType(mem, MODIFIED_GS);
372:   if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);
373:   return(0);
374: }

376: /* type of CVODE linear multistep method */
377: const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",NULL};
378: /* type of G-S orthogonalization used by CVODE linear solver */
379: const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",NULL};

381: PetscErrorCode TSSetFromOptions_Sundials(PetscOptionItems *PetscOptionsObject,TS ts)
382: {
383:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
385:   int            indx;
386:   PetscBool      flag;
387:   PC             pc;

390:   PetscOptionsHead(PetscOptionsObject,"SUNDIALS ODE solver options");
391:   PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);
392:   if (flag) {
393:     TSSundialsSetType(ts,(TSSundialsLmmType)indx);
394:   }
395:   PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);
396:   if (flag) {
397:     TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
398:   }
399:   PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);
400:   PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);
401:   PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);
402:   PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);
403:   PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,NULL);
404:   PetscOptionsInt("-ts_sundials_maxord","Max Order for BDF/Adams method","TSSundialsSetMaxOrd",cvode->maxord,&cvode->maxord,NULL);
405:   PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,NULL);
406:   PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internal steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);
407:   PetscOptionsTail();
408:   TSSundialsGetPC(ts,&pc);
409:   PCSetFromOptions(pc);
410:   return(0);
411: }

413: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
414: {
415:   TS_Sundials    *cvode = (TS_Sundials*)ts->data;
417:   char           *type;
418:   char           atype[] = "Adams";
419:   char           btype[] = "BDF: backward differentiation formula";
420:   PetscBool      iascii,isstring;
421:   long int       nsteps,its,nfevals,nlinsetups,nfails,itmp;
422:   PetscInt       qlast,qcur;
423:   PetscReal      hinused,hlast,hcur,tcur,tolsfac;
424:   PC             pc;

427:   if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
428:   else                                     type = btype;

430:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
431:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
432:   if (iascii) {
433:     PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");
434:     PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);
435:     PetscViewerASCIIPrintf(viewer,"Sundials maxord %D\n",cvode->maxord);
436:     PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",(double)cvode->abstol,(double)cvode->reltol);
437:     PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",(double)cvode->linear_tol);
438:     PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);
439:     if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
440:       PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
441:     } else {
442:       PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
443:     }
444:     if (cvode->mindt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",(double)cvode->mindt);}
445:     if (cvode->maxdt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",(double)cvode->maxdt);}

447:     /* Outputs from CVODE, CVSPILS */
448:     CVodeGetTolScaleFactor(cvode->mem,&tolsfac);
449:     PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);
450:     CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,
451:                                    &nlinsetups,&nfails,&qlast,&qcur,
452:                                    &hinused,&hlast,&hcur,&tcur);
453:     PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);
454:     PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);
455:     PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);
456:     PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);

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

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

467:     TSSundialsGetPC(ts,&pc);
468:     PCView(pc,viewer);
469:     CVSpilsGetNumPrecEvals(cvode->mem,&itmp);
470:     PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);
471:     CVSpilsGetNumPrecSolves(cvode->mem,&itmp);
472:     PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);

474:     CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);
475:     PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);
476:     CVSpilsGetNumRhsEvals(cvode->mem,&itmp);
477:     PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);
478:   } else if (isstring) {
479:     PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
480:   }
481:   return(0);
482: }

484: /* --------------------------------------------------------------------------*/
485: PetscErrorCode  TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
486: {
487:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

490:   cvode->cvode_type = type;
491:   return(0);
492: }

494: PetscErrorCode  TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
495: {
496:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

499:   cvode->maxl = maxl;
500:   return(0);
501: }

503: PetscErrorCode  TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
504: {
505:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

508:   cvode->linear_tol = tol;
509:   return(0);
510: }

512: PetscErrorCode  TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
513: {
514:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

517:   cvode->gtype = type;
518:   return(0);
519: }

521: PetscErrorCode  TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
522: {
523:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

526:   if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
527:   if (rel != PETSC_DECIDE)  cvode->reltol = rel;
528:   return(0);
529: }

531: PetscErrorCode  TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
532: {
533:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

536:   cvode->mindt = mindt;
537:   return(0);
538: }

540: PetscErrorCode  TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
541: {
542:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

545:   cvode->maxdt = maxdt;
546:   return(0);
547: }
548: PetscErrorCode  TSSundialsGetPC_Sundials(TS ts,PC *pc)
549: {
550:   SNES           snes;
551:   KSP            ksp;

555:   TSGetSNES(ts,&snes);
556:   SNESGetKSP(snes,&ksp);
557:   KSPGetPC(ksp,pc);
558:   return(0);
559: }

561: PetscErrorCode  TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
562: {
564:   if (nonlin) *nonlin = ts->snes_its;
565:   if (lin)    *lin    = ts->ksp_its;
566:   return(0);
567: }

569: PetscErrorCode  TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
570: {
571:   TS_Sundials *cvode = (TS_Sundials*)ts->data;

574:   cvode->monitorstep = s;
575:   return(0);
576: }
577: /* -------------------------------------------------------------------------------------------*/

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

582:    Not Collective

584:    Input parameters:
585: .    ts     - the time-step context

587:    Output Parameters:
588: +   nonlin - number of nonlinear iterations
589: -   lin    - number of linear iterations

591:    Level: advanced

593:    Notes:
594:     These return the number since the creation of the TS object

596: .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
597:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
598:           TSSundialsGetIterations(), TSSundialsSetType(),
599:           TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()

601: @*/
602: PetscErrorCode  TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
603: {

607:   PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));
608:   return(0);
609: }

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

614:    Logically Collective on TS

616:    Input parameters:
617: +    ts     - the time-step context
618: -    type   - one of  SUNDIALS_ADAMS or SUNDIALS_BDF

620:    Level: intermediate

622: .seealso: TSSundialsGetIterations(),  TSSundialsSetMaxl(),
623:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
624:           TSSundialsGetIterations(), TSSundialsSetType(),
625:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
626:           TSSetExactFinalTime()
627: @*/
628: PetscErrorCode  TSSundialsSetType(TS ts,TSSundialsLmmType type)
629: {

633:   PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));
634:   return(0);
635: }

637: /*@
638:    TSSundialsSetMaxord - Sets the maximum order for BDF/Adams method used by SUNDIALS.

640:    Logically Collective on TS

642:    Input parameters:
643: +    ts      - the time-step context
644: -    maxord  - maximum order of BDF / Adams method

646:    Level: advanced

648: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
649:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
650:           TSSundialsGetIterations(), TSSundialsSetType(),
651:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
652:           TSSetExactFinalTime()

654: @*/
655: PetscErrorCode  TSSundialsSetMaxord(TS ts,PetscInt maxord)
656: {

661:   PetscTryMethod(ts,"TSSundialsSetMaxOrd_C",(TS,PetscInt),(ts,maxord));
662:   return(0);
663: }

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

670:    Logically Collective on TS

672:    Input parameters:
673: +    ts      - the time-step context
674: -    maxl - number of direction vectors (the dimension of Krylov subspace).

676:    Level: advanced

678: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
679:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
680:           TSSundialsGetIterations(), TSSundialsSetType(),
681:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
682:           TSSetExactFinalTime()

684: @*/
685: PetscErrorCode  TSSundialsSetMaxl(TS ts,PetscInt maxl)
686: {

691:   PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));
692:   return(0);
693: }

695: /*@
696:    TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
697:        system by SUNDIALS.

699:    Logically Collective on TS

701:    Input parameters:
702: +    ts     - the time-step context
703: -    tol    - the factor by which the tolerance on the nonlinear solver is
704:              multiplied to get the tolerance on the linear solver, .05 by default.

706:    Level: advanced

708: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
709:           TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
710:           TSSundialsGetIterations(), TSSundialsSetType(),
711:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
712:           TSSetExactFinalTime()

714: @*/
715: PetscErrorCode  TSSundialsSetLinearTolerance(TS ts,double tol)
716: {

721:   PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));
722:   return(0);
723: }

725: /*@
726:    TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
727:         in GMRES method by SUNDIALS linear solver.

729:    Logically Collective on TS

731:    Input parameters:
732: +    ts  - the time-step context
733: -    type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS

735:    Level: advanced

737: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
738:           TSSundialsSetLinearTolerance(),  TSSundialsSetTolerance(),
739:           TSSundialsGetIterations(), TSSundialsSetType(),
740:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
741:           TSSetExactFinalTime()

743: @*/
744: PetscErrorCode  TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
745: {

749:   PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));
750:   return(0);
751: }

753: /*@
754:    TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
755:                          Sundials for error control.

757:    Logically Collective on TS

759:    Input parameters:
760: +    ts  - the time-step context
761: .    aabs - the absolute tolerance
762: -    rel - the relative tolerance

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

767:    Level: intermediate

769: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
770:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
771:           TSSundialsGetIterations(), TSSundialsSetType(),
772:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
773:           TSSetExactFinalTime()

775: @*/
776: PetscErrorCode  TSSundialsSetTolerance(TS ts,double aabs,double rel)
777: {

781:   PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));
782:   return(0);
783: }

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

788:    Input Parameter:
789: .    ts - the time-step context

791:    Output Parameter:
792: .    pc - the preconditioner context

794:    Level: advanced

796: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
797:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
798:           TSSundialsGetIterations(), TSSundialsSetType(),
799:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
800: @*/
801: PetscErrorCode  TSSundialsGetPC(TS ts,PC *pc)
802: {

806:   PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));
807:   return(0);
808: }

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

813:    Input Parameter:
814: +   ts - the time-step context
815: -   mindt - lowest time step if positive, negative to deactivate

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

821:    Level: beginner

823: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
824: @*/
825: PetscErrorCode  TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
826: {

830:   PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));
831:   return(0);
832: }

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

837:    Input Parameter:
838: +   ts - the time-step context
839: -   maxdt - lowest time step if positive, negative to deactivate

841:    Level: beginner

843: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
844: @*/
845: PetscErrorCode  TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
846: {

850:   PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
851:   return(0);
852: }

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

857:    Input Parameter:
858: +   ts - the time-step context
859: -   ft - PETSC_TRUE if monitor, else PETSC_FALSE

861:    Level: beginner

863: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
864:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
865:           TSSundialsGetIterations(), TSSundialsSetType(),
866:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
867: @*/
868: PetscErrorCode  TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
869: {

873:   PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));
874:   return(0);
875: }
876: /* -------------------------------------------------------------------------------------------*/
877: /*MC
878:       TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)

880:    Options Database:
881: +    -ts_sundials_type <bdf,adams> -
882: .    -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
883: .    -ts_sundials_atol <tol> - Absolute tolerance for convergence
884: .    -ts_sundials_rtol <tol> - Relative tolerance for convergence
885: .    -ts_sundials_linear_tolerance <tol> -
886: .    -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
887: -    -ts_sundials_monitor_steps - Monitor SUNDIALS internal steps

889:     Notes:
890:     This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply,
891:            only PETSc PC options.

893:     Level: beginner

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

898: M*/
899: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
900: {
901:   TS_Sundials    *cvode;
903:   PC             pc;

906:   ts->ops->reset          = TSReset_Sundials;
907:   ts->ops->destroy        = TSDestroy_Sundials;
908:   ts->ops->view           = TSView_Sundials;
909:   ts->ops->setup          = TSSetUp_Sundials;
910:   ts->ops->step           = TSStep_Sundials;
911:   ts->ops->interpolate    = TSInterpolate_Sundials;
912:   ts->ops->setfromoptions = TSSetFromOptions_Sundials;
913:   ts->default_adapt_type  = TSADAPTNONE;

915:   PetscNewLog(ts,&cvode);

917:   ts->usessnes = PETSC_TRUE;

919:   ts->data           = (void*)cvode;
920:   cvode->cvode_type  = SUNDIALS_BDF;
921:   cvode->gtype       = SUNDIALS_CLASSICAL_GS;
922:   cvode->maxl        = 5;
923:   cvode->maxord      = PETSC_DEFAULT;
924:   cvode->linear_tol  = .05;
925:   cvode->monitorstep = PETSC_TRUE;

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

929:   cvode->mindt = -1.;
930:   cvode->maxdt = -1.;

932:   /* set tolerance for Sundials */
933:   cvode->reltol = 1e-6;
934:   cvode->abstol = 1e-6;

936:   /* set PCNONE as default pctype */
937:   TSSundialsGetPC_Sundials(ts,&pc);
938:   PCSetType(pc,PCNONE);

940:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);
941:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);
942:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);
943:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);
944:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);
945:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);
946:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);
947:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);
948:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);
949:   PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);
950:   return(0);
951: }