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, realtype _gamma, void *P_data, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3)
 15: {
 16:   TS           ts    = (TS)P_data;
 17:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
 18:   PC           pc;
 19:   Mat          J, P;
 20:   Vec          yy = cvode->w1, yydot = cvode->ydot;
 21:   PetscReal    gm = (PetscReal)_gamma;
 22:   PetscScalar *y_data;

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

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

 49:   /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
 50:   r_data = (PetscScalar *)N_VGetArrayPointer(r);
 51:   z_data = (PetscScalar *)N_VGetArrayPointer(z);
 52:   VecPlaceArray(rr, r_data);
 53:   VecPlaceArray(zz, z_data);

 55:   /* Solve the Px=r and put the result in zz */
 56:   TSSundialsGetPC(ts, &pc);
 57:   PCApply(pc, rr, zz);
 58:   VecResetArray(rr);
 59:   VecResetArray(zz);
 60:   return 0;
 61: }

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

 77:   PetscObjectGetComm((PetscObject)ts, &comm);
 78:   /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
 79:   y_data    = (PetscScalar *)N_VGetArrayPointer(y);
 80:   ydot_data = (PetscScalar *)N_VGetArrayPointer(ydot);
 81:   comm, VecPlaceArray(yy, y_data);
 82:   comm, VecPlaceArray(yyd, ydot_data);

 84:   /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
 85:   TSGetDM(ts, &dm);
 86:   DMGetDMTS(dm, &tsdm);
 87:   DMTSGetIFunction(dm, &ifunction, NULL);
 88:   if (!ifunction) {
 89:     TSComputeRHSFunction(ts, t, yy, yyd);
 90:   } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
 91:     VecZeroEntries(yydot);
 92:     comm, TSComputeIFunction(ts, t, yy, yydot, yyd, PETSC_FALSE);
 93:     VecScale(yyd, -1.);
 94:   }
 95:   comm, VecResetArray(yy);
 96:   comm, VecResetArray(yyd);
 97:   return 0;
 98: }

100: /*
101:        TSStep_Sundials - Calls Sundials to integrate the ODE.
102: */
103: PetscErrorCode TSStep_Sundials(TS ts)
104: {
105:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
106:   PetscInt     flag;
107:   long int     nits, lits, nsteps;
108:   realtype     t, tout;
109:   PetscScalar *y_data;
110:   void        *mem;

112:   mem  = cvode->mem;
113:   tout = ts->max_time;
114:   VecGetArray(ts->vec_sol, &y_data);
115:   N_VSetArrayPointer((realtype *)y_data, cvode->y);
116:   VecRestoreArray(ts->vec_sol, NULL);

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

123:   if (flag) { /* display error message */
124:     switch (flag) {
125:     case CV_ILL_INPUT:
126:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_ILL_INPUT");
127:       break;
128:     case CV_TOO_CLOSE:
129:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_TOO_CLOSE");
130:       break;
131:     case CV_TOO_MUCH_WORK: {
132:       PetscReal tcur;
133:       CVodeGetNumSteps(mem, &nsteps);
134:       CVodeGetCurrentTime(mem, &tcur);
135:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_TOO_MUCH_WORK. At t=%g, nsteps %ld exceeds maxstep %" PetscInt_FMT ". Increase '-ts_max_steps <>' or modify TSSetMaxSteps()", (double)tcur, nsteps, ts->max_steps);
136:     } break;
137:     case CV_TOO_MUCH_ACC:
138:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_TOO_MUCH_ACC");
139:       break;
140:     case CV_ERR_FAILURE:
141:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_ERR_FAILURE");
142:       break;
143:     case CV_CONV_FAILURE:
144:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_CONV_FAILURE");
145:       break;
146:     case CV_LINIT_FAIL:
147:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_LINIT_FAIL");
148:       break;
149:     case CV_LSETUP_FAIL:
150:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_LSETUP_FAIL");
151:       break;
152:     case CV_LSOLVE_FAIL:
153:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_LSOLVE_FAIL");
154:       break;
155:     case CV_RHSFUNC_FAIL:
156:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_RHSFUNC_FAIL");
157:       break;
158:     case CV_FIRST_RHSFUNC_ERR:
159:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_FIRST_RHSFUNC_ERR");
160:       break;
161:     case CV_REPTD_RHSFUNC_ERR:
162:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_REPTD_RHSFUNC_ERR");
163:       break;
164:     case CV_UNREC_RHSFUNC_ERR:
165:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_UNREC_RHSFUNC_ERR");
166:       break;
167:     case CV_RTFUNC_FAIL:
168:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_RTFUNC_FAIL");
169:       break;
170:     default:
171:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, flag %d", flag);
172:     }
173:   }

175:   /* log inner nonlinear and linear iterations */
176:   CVodeGetNumNonlinSolvIters(mem, &nits);
177:   CVSpilsGetNumLinIters(mem, &lits);
178:   ts->snes_its += nits;
179:   ts->ksp_its = lits;

181:   /* copy the solution from cvode->y to cvode->update and sol */
182:   VecPlaceArray(cvode->w1, y_data);
183:   VecCopy(cvode->w1, cvode->update);
184:   VecResetArray(cvode->w1);
185:   VecCopy(cvode->update, ts->vec_sol);

187:   ts->time_step = t - ts->ptime;
188:   ts->ptime     = t;

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

195: static PetscErrorCode TSInterpolate_Sundials(TS ts, PetscReal t, Vec X)
196: {
197:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
198:   N_Vector     y;
199:   PetscScalar *x_data;
200:   PetscInt     glosize, locsize;

202:   /* get the vector size */
203:   VecGetSize(X, &glosize);
204:   VecGetLocalSize(X, &locsize);
205:   VecGetArray(X, &x_data);

207:   /* Initialize N_Vec y with x_data */
208:   if (cvode->use_dense) {
209:     PetscMPIInt size;

211:     MPI_Comm_size(PETSC_COMM_WORLD, &size);
213:     y = N_VMake_Serial(locsize, (realtype *)x_data);
214:   } else {
215:     y = N_VMake_Parallel(cvode->comm_sundials, locsize, glosize, (realtype *)x_data);
216:   }


220:   CVodeGetDky(cvode->mem, t, 0, y);
221:   VecRestoreArray(X, &x_data);
222:   return 0;
223: }

225: PetscErrorCode TSReset_Sundials(TS ts)
226: {
227:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

229:   VecDestroy(&cvode->update);
230:   VecDestroy(&cvode->ydot);
231:   VecDestroy(&cvode->w1);
232:   VecDestroy(&cvode->w2);
233:   if (cvode->mem) CVodeFree(&cvode->mem);
234:   return 0;
235: }

237: PetscErrorCode TSDestroy_Sundials(TS ts)
238: {
239:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

241:   TSReset_Sundials(ts);
242:   MPI_Comm_free(&(cvode->comm_sundials));
243:   PetscFree(ts->data);
244:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetType_C", NULL);
245:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxl_C", NULL);
246:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetLinearTolerance_C", NULL);
247:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetGramSchmidtType_C", NULL);
248:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetTolerance_C", NULL);
249:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMinTimeStep_C", NULL);
250:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxTimeStep_C", NULL);
251:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetPC_C", NULL);
252:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetIterations_C", NULL);
253:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsMonitorInternalSteps_C", NULL);
254:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetUseDense_C", NULL);
255:   return 0;
256: }

258: PetscErrorCode TSSetUp_Sundials(TS ts)
259: {
260:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
261:   PetscInt     glosize, locsize, i, flag;
262:   PetscScalar *y_data, *parray;
263:   void        *mem;
264:   PC           pc;
265:   PCType       pctype;
266:   PetscBool    pcnone;


270:   /* get the vector size */
271:   VecGetSize(ts->vec_sol, &glosize);
272:   VecGetLocalSize(ts->vec_sol, &locsize);

274:   /* allocate the memory for N_Vec y */
275:   if (cvode->use_dense) {
276:     PetscMPIInt size;

278:     MPI_Comm_size(PETSC_COMM_WORLD, &size);
280:     cvode->y = N_VNew_Serial(locsize);
281:   } else {
282:     cvode->y = N_VNew_Parallel(cvode->comm_sundials, locsize, glosize);
283:   }

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

292:   VecDuplicate(ts->vec_sol, &cvode->update);
293:   VecDuplicate(ts->vec_sol, &cvode->ydot);

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

303:   /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */
304:   mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
306:   cvode->mem = mem;

308:   /* Set the pointer to user-defined data */
309:   flag = CVodeSetUserData(mem, ts);

312:   /* Sundials may choose to use a smaller initial step, but will never use a larger step. */
313:   flag = CVodeSetInitStep(mem, (realtype)ts->time_step);
315:   if (cvode->mindt > 0) {
316:     flag = CVodeSetMinStep(mem, (realtype)cvode->mindt);
317:     if (flag) {
320:       SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_LIB, "CVodeSetMinStep() failed");
321:     }
322:   }
323:   if (cvode->maxdt > 0) {
324:     flag = CVodeSetMaxStep(mem, (realtype)cvode->maxdt);
326:   }

328:   /* Call CVodeInit to initialize the integrator memory and specify the
329:    * user's right hand side function in u'=f(t,u), the initial time T0, and
330:    * the initial dependent variable vector cvode->y */
331:   flag = CVodeInit(mem, TSFunction_Sundials, ts->ptime, cvode->y);

334:   /* specifies scalar relative and absolute tolerances */
335:   flag = CVodeSStolerances(mem, cvode->reltol, cvode->abstol);

338:   /* Specify max order of BDF / ADAMS method */
339:   if (cvode->maxord != PETSC_DEFAULT) {
340:     flag = CVodeSetMaxOrd(mem, cvode->maxord);
342:   }

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

348:   if (cvode->use_dense) {
349:     /* call CVDense to use a dense linear solver. */
350:     PetscMPIInt size;

352:     MPI_Comm_size(PETSC_COMM_WORLD, &size);
354:     flag = CVDense(mem, locsize);
356:   } else {
357:     /* call CVSpgmr to use GMRES as the linear solver.        */
358:     /* setup the ode integrator with the given preconditioner */
359:     TSSundialsGetPC(ts, &pc);
360:     PCGetType(pc, &pctype);
361:     PetscObjectTypeCompare((PetscObject)pc, PCNONE, &pcnone);
362:     if (pcnone) {
363:       flag = CVSpgmr(mem, PREC_NONE, 0);
365:     } else {
366:       flag = CVSpgmr(mem, PREC_LEFT, cvode->maxl);

369:       /* Set preconditioner and solve routines Precond and PSolve,
370:          and the pointer to the user-defined block data */
371:       flag = CVSpilsSetPreconditioner(mem, TSPrecond_Sundials, TSPSolve_Sundials);
373:     }
374:   }
375:   return 0;
376: }

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

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

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

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

422:   if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
423:   else type = btype;

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

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

455:     CVodeGetNonlinSolvStats(cvode->mem, &its, &nfails);
456:     PetscViewerASCIIPrintf(viewer, "Sundials no. of nonlinear solver iterations %ld\n", its);
457:     PetscViewerASCIIPrintf(viewer, "Sundials no. of nonlinear convergence failure %ld\n", nfails);
458:     if (!cvode->use_dense) {
459:       CVSpilsGetNumLinIters(cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
460:       PetscViewerASCIIPrintf(viewer, "Sundials no. of linear iterations %ld\n", its);
461:       CVSpilsGetNumConvFails(cvode->mem, &itmp);
462:       PetscViewerASCIIPrintf(viewer, "Sundials no. of linear convergence failures %ld\n", itmp);

464:       TSSundialsGetPC(ts, &pc);
465:       PCView(pc, viewer);
466:       CVSpilsGetNumPrecEvals(cvode->mem, &itmp);
467:       PetscViewerASCIIPrintf(viewer, "Sundials no. of preconditioner evaluations %ld\n", itmp);
468:       CVSpilsGetNumPrecSolves(cvode->mem, &itmp);
469:       PetscViewerASCIIPrintf(viewer, "Sundials no. of preconditioner solves %ld\n", itmp);
470:     }
471:     CVSpilsGetNumJtimesEvals(cvode->mem, &itmp);
472:     PetscViewerASCIIPrintf(viewer, "Sundials no. of Jacobian-vector product evaluations %ld\n", itmp);
473:     CVSpilsGetNumRhsEvals(cvode->mem, &itmp);
474:     PetscViewerASCIIPrintf(viewer, "Sundials no. of rhs calls for finite diff. Jacobian-vector evals %ld\n", itmp);
475:   } else if (isstring) {
476:     PetscViewerStringSPrintf(viewer, "Sundials type %s", type);
477:   }
478:   return 0;
479: }

481: /* --------------------------------------------------------------------------*/
482: PetscErrorCode TSSundialsSetType_Sundials(TS ts, TSSundialsLmmType type)
483: {
484:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

486:   cvode->cvode_type = type;
487:   return 0;
488: }

490: PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts, PetscInt maxl)
491: {
492:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

494:   cvode->maxl = maxl;
495:   return 0;
496: }

498: PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts, PetscReal tol)
499: {
500:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

502:   cvode->linear_tol = tol;
503:   return 0;
504: }

506: PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts, TSSundialsGramSchmidtType type)
507: {
508:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

510:   cvode->gtype = type;
511:   return 0;
512: }

514: PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts, PetscReal aabs, PetscReal rel)
515: {
516:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

518:   if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
519:   if (rel != PETSC_DECIDE) cvode->reltol = rel;
520:   return 0;
521: }

523: PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts, PetscReal mindt)
524: {
525:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

527:   cvode->mindt = mindt;
528:   return 0;
529: }

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

535:   cvode->maxdt = maxdt;
536:   return 0;
537: }

539: PetscErrorCode TSSundialsSetUseDense_Sundials(TS ts, PetscBool use_dense)
540: {
541:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

543:   cvode->use_dense = use_dense;
544:   return 0;
545: }

547: PetscErrorCode TSSundialsGetPC_Sundials(TS ts, PC *pc)
548: {
549:   SNES snes;
550:   KSP  ksp;

552:   TSGetSNES(ts, &snes);
553:   SNESGetKSP(snes, &ksp);
554:   KSPGetPC(ksp, pc);
555:   return 0;
556: }

558: PetscErrorCode TSSundialsGetIterations_Sundials(TS ts, int *nonlin, int *lin)
559: {
560:   if (nonlin) *nonlin = ts->snes_its;
561:   if (lin) *lin = ts->ksp_its;
562:   return 0;
563: }

565: PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts, PetscBool s)
566: {
567:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

569:   cvode->monitorstep = s;
570:   return 0;
571: }
572: /* -------------------------------------------------------------------------------------------*/

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

577:    Not Collective

579:    Input Parameter:
580: .    ts     - the time-step context

582:    Output Parameters:
583: +   nonlin - number of nonlinear iterations
584: -   lin    - number of linear iterations

586:    Level: advanced

588:    Notes:
589:     These return the number since the creation of the TS object

591: .seealso: `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
592:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
593:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
594:           `TSSundialsSetLinearTolerance()`, `TSSundialsGetPC()`, `TSSetExactFinalTime()`

596: @*/
597: PetscErrorCode TSSundialsGetIterations(TS ts, int *nonlin, int *lin)
598: {
599:   PetscUseMethod(ts, "TSSundialsGetIterations_C", (TS, int *, int *), (ts, nonlin, lin));
600:   return 0;
601: }

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

606:    Logically Collective on TS

608:    Input Parameters:
609: +    ts     - the time-step context
610: -    type   - one of  SUNDIALS_ADAMS or SUNDIALS_BDF

612:    Level: intermediate

614: .seealso: `TSSundialsGetIterations()`, `TSSundialsSetMaxl()`,
615:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
616:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
617:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
618:           `TSSetExactFinalTime()`
619: @*/
620: PetscErrorCode TSSundialsSetType(TS ts, TSSundialsLmmType type)
621: {
622:   PetscTryMethod(ts, "TSSundialsSetType_C", (TS, TSSundialsLmmType), (ts, type));
623:   return 0;
624: }

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

629:    Logically Collective on TS

631:    Input Parameters:
632: +    ts      - the time-step context
633: -    maxord  - maximum order of BDF / Adams method

635:    Level: advanced

637: .seealso: `TSSundialsGetIterations()`, `TSSundialsSetType()`,
638:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
639:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
640:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
641:           `TSSetExactFinalTime()`

643: @*/
644: PetscErrorCode TSSundialsSetMaxord(TS ts, PetscInt maxord)
645: {
647:   PetscTryMethod(ts, "TSSundialsSetMaxOrd_C", (TS, PetscInt), (ts, maxord));
648:   return 0;
649: }

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

656:    Logically Collective on TS

658:    Input Parameters:
659: +    ts      - the time-step context
660: -    maxl - number of direction vectors (the dimension of Krylov subspace).

662:    Level: advanced

664: .seealso: `TSSundialsGetIterations()`, `TSSundialsSetType()`,
665:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
666:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
667:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
668:           `TSSetExactFinalTime()`

670: @*/
671: PetscErrorCode TSSundialsSetMaxl(TS ts, PetscInt maxl)
672: {
674:   PetscTryMethod(ts, "TSSundialsSetMaxl_C", (TS, PetscInt), (ts, maxl));
675:   return 0;
676: }

678: /*@
679:    TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
680:        system by SUNDIALS.

682:    Logically Collective on TS

684:    Input Parameters:
685: +    ts     - the time-step context
686: -    tol    - the factor by which the tolerance on the nonlinear solver is
687:              multiplied to get the tolerance on the linear solver, .05 by default.

689:    Level: advanced

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

697: @*/
698: PetscErrorCode TSSundialsSetLinearTolerance(TS ts, PetscReal tol)
699: {
701:   PetscTryMethod(ts, "TSSundialsSetLinearTolerance_C", (TS, PetscReal), (ts, tol));
702:   return 0;
703: }

705: /*@
706:    TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
707:         in GMRES method by SUNDIALS linear solver.

709:    Logically Collective on TS

711:    Input Parameters:
712: +    ts  - the time-step context
713: -    type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS

715:    Level: advanced

717: .seealso: `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
718:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`,
719:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
720:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
721:           `TSSetExactFinalTime()`

723: @*/
724: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts, TSSundialsGramSchmidtType type)
725: {
726:   PetscTryMethod(ts, "TSSundialsSetGramSchmidtType_C", (TS, TSSundialsGramSchmidtType), (ts, type));
727:   return 0;
728: }

730: /*@
731:    TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
732:                          Sundials for error control.

734:    Logically Collective on TS

736:    Input Parameters:
737: +    ts  - the time-step context
738: .    aabs - the absolute tolerance
739: -    rel - the relative tolerance

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

744:    Level: intermediate

746: .seealso: `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetGMRESMaxl()`,
747:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`,
748:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
749:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
750:           `TSSetExactFinalTime()`

752: @*/
753: PetscErrorCode TSSundialsSetTolerance(TS ts, PetscReal aabs, PetscReal rel)
754: {
755:   PetscTryMethod(ts, "TSSundialsSetTolerance_C", (TS, PetscReal, PetscReal), (ts, aabs, rel));
756:   return 0;
757: }

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

762:    Input Parameter:
763: .    ts - the time-step context

765:    Output Parameter:
766: .    pc - the preconditioner context

768:    Level: advanced

770: .seealso: `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
771:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
772:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
773:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`
774: @*/
775: PetscErrorCode TSSundialsGetPC(TS ts, PC *pc)
776: {
777:   PetscUseMethod(ts, "TSSundialsGetPC_C", (TS, PC *), (ts, pc));
778:   return 0;
779: }

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

784:    Input Parameters:
785: +   ts - the time-step context
786: -   mindt - lowest time step if positive, negative to deactivate

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

792:    Level: beginner

794: .seealso: `TSSundialsSetType()`, `TSSundialsSetTolerance()`,
795: @*/
796: PetscErrorCode TSSundialsSetMinTimeStep(TS ts, PetscReal mindt)
797: {
798:   PetscTryMethod(ts, "TSSundialsSetMinTimeStep_C", (TS, PetscReal), (ts, mindt));
799:   return 0;
800: }

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

805:    Input Parameters:
806: +   ts - the time-step context
807: -   maxdt - lowest time step if positive, negative to deactivate

809:    Level: beginner

811: .seealso: `TSSundialsSetType()`, `TSSundialsSetTolerance()`,
812: @*/
813: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts, PetscReal maxdt)
814: {
815:   PetscTryMethod(ts, "TSSundialsSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
816:   return 0;
817: }

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

822:    Input Parameters:
823: +   ts - the time-step context
824: -   ft - PETSC_TRUE if monitor, else PETSC_FALSE

826:    Level: beginner

828: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
829:           TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
830:           TSSundialsGetIterations(), TSSundialsSetType(),
831:           TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
832: @*/
833: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts, PetscBool ft)
834: {
835:   PetscTryMethod(ts, "TSSundialsMonitorInternalSteps_C", (TS, PetscBool), (ts, ft));
836:   return 0;
837: }

839: /*@
840:    TSSundialsSetUseDense - Set a flag to use a dense linear solver in SUNDIALS (serial only)

842:    Logically Collective

844:    Input Parameters:
845: +    ts         - the time-step context
846: -    use_dense  - PETSC_TRUE to use the dense solver

848:    Level: advanced

850: .seealso: `TSSUNDIALS`

852: @*/
853: PetscErrorCode TSSundialsSetUseDense(TS ts, PetscBool use_dense)
854: {
856:   PetscTryMethod(ts, "TSSundialsSetUseDense_C", (TS, PetscBool), (ts, use_dense));
857:   return 0;
858: }

860: /* -------------------------------------------------------------------------------------------*/
861: /*MC
862:       TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)

864:    Options Database:
865: +    -ts_sundials_type <bdf,adams> -
866: .    -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
867: .    -ts_sundials_atol <tol> - Absolute tolerance for convergence
868: .    -ts_sundials_rtol <tol> - Relative tolerance for convergence
869: .    -ts_sundials_linear_tolerance <tol> -
870: .    -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
871: .    -ts_sundials_monitor_steps - Monitor SUNDIALS internal steps
872: -    -ts_sundials_use_dense - Use a dense linear solver within CVODE (serial only)

874:     Notes:
875:     This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply,
876:            only PETSc PC options.

878:     Level: beginner

880: .seealso: `TSCreate()`, `TS`, `TSSetType()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`, `TSSundialsSetLinearTolerance()`,
881:           `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`, `TSSundialsGetIterations()`, `TSSetExactFinalTime()`

883: M*/
884: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
885: {
886:   TS_Sundials *cvode;
887:   PC           pc;

889:   ts->ops->reset          = TSReset_Sundials;
890:   ts->ops->destroy        = TSDestroy_Sundials;
891:   ts->ops->view           = TSView_Sundials;
892:   ts->ops->setup          = TSSetUp_Sundials;
893:   ts->ops->step           = TSStep_Sundials;
894:   ts->ops->interpolate    = TSInterpolate_Sundials;
895:   ts->ops->setfromoptions = TSSetFromOptions_Sundials;
896:   ts->default_adapt_type  = TSADAPTNONE;

898:   PetscNew(&cvode);

900:   ts->usessnes = PETSC_TRUE;

902:   ts->data           = (void *)cvode;
903:   cvode->cvode_type  = SUNDIALS_BDF;
904:   cvode->gtype       = SUNDIALS_CLASSICAL_GS;
905:   cvode->maxl        = 5;
906:   cvode->maxord      = PETSC_DEFAULT;
907:   cvode->linear_tol  = .05;
908:   cvode->monitorstep = PETSC_TRUE;
909:   cvode->use_dense   = PETSC_FALSE;

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

913:   cvode->mindt = -1.;
914:   cvode->maxdt = -1.;

916:   /* set tolerance for Sundials */
917:   cvode->reltol = 1e-6;
918:   cvode->abstol = 1e-6;

920:   /* set PCNONE as default pctype */
921:   TSSundialsGetPC_Sundials(ts, &pc);
922:   PCSetType(pc, PCNONE);

924:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetType_C", TSSundialsSetType_Sundials);
925:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxl_C", TSSundialsSetMaxl_Sundials);
926:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetLinearTolerance_C", TSSundialsSetLinearTolerance_Sundials);
927:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetGramSchmidtType_C", TSSundialsSetGramSchmidtType_Sundials);
928:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetTolerance_C", TSSundialsSetTolerance_Sundials);
929:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMinTimeStep_C", TSSundialsSetMinTimeStep_Sundials);
930:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxTimeStep_C", TSSundialsSetMaxTimeStep_Sundials);
931:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetPC_C", TSSundialsGetPC_Sundials);
932:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetIterations_C", TSSundialsGetIterations_Sundials);
933:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsMonitorInternalSteps_C", TSSundialsMonitorInternalSteps_Sundials);
934:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetUseDense_C", TSSundialsSetUseDense_Sundials);
935:   return 0;
936: }