Actual source code: sundials.c

  1: /*
  2:     Provides a PETSc interface to version 2.5 of SUNDIALS/CVODE solver (a very old version)
  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 `TSSUNDIALS`.

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:    Note:
589:     These return the number since the creation of the `TS` object

591: .seealso: [](chapter_ts), `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
592:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
593:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
594:           `TSSundialsSetLinearTolerance()`, `TSSundialsGetPC()`, `TSSetExactFinalTime()`
595: @*/
596: PetscErrorCode TSSundialsGetIterations(TS ts, int *nonlin, int *lin)
597: {
598:   PetscUseMethod(ts, "TSSundialsGetIterations_C", (TS, int *, int *), (ts, nonlin, lin));
599:   return 0;
600: }

602: /*@
603:    TSSundialsSetType - Sets the method that `TSSUNDIALS` will use for integration.

605:    Logically Collective

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

611:    Level: intermediate

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

625: /*@
626:    TSSundialsSetMaxord - Sets the maximum order for BDF/Adams method used by `TSSUNDIALS`.

628:    Logically Collective

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

634:    Level: advanced

636: .seealso: [](chapter_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`,
637:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
638:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
639:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
640:           `TSSetExactFinalTime()`
641: @*/
642: PetscErrorCode TSSundialsSetMaxord(TS ts, PetscInt maxord)
643: {
645:   PetscTryMethod(ts, "TSSundialsSetMaxOrd_C", (TS, PetscInt), (ts, maxord));
646:   return 0;
647: }

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

654:    Logically Collective

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

660:    Level: advanced

662: .seealso: [](chapter_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`,
663:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
664:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
665:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
666:           `TSSetExactFinalTime()`
667: @*/
668: PetscErrorCode TSSundialsSetMaxl(TS ts, PetscInt maxl)
669: {
671:   PetscTryMethod(ts, "TSSundialsSetMaxl_C", (TS, PetscInt), (ts, maxl));
672:   return 0;
673: }

675: /*@
676:    TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
677:        system by `TSSUNDIALS`.

679:    Logically Collective

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

686:    Level: advanced

688: .seealso: [](chapter_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
689:           `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
690:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
691:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
692:           `TSSetExactFinalTime()`
693: @*/
694: PetscErrorCode TSSundialsSetLinearTolerance(TS ts, PetscReal tol)
695: {
697:   PetscTryMethod(ts, "TSSundialsSetLinearTolerance_C", (TS, PetscReal), (ts, tol));
698:   return 0;
699: }

701: /*@
702:    TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
703:         in GMRES method by `TSSUNDIALS` linear solver.

705:    Logically Collective

707:    Input Parameters:
708: +    ts  - the time-step context
709: -    type - either `SUNDIALS_MODIFIED_GS` or `SUNDIALS_CLASSICAL_GS`

711:    Level: advanced

713: .seealso: [](chapter_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
714:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`,
715:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
716:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
717:           `TSSetExactFinalTime()`
718: @*/
719: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts, TSSundialsGramSchmidtType type)
720: {
721:   PetscTryMethod(ts, "TSSundialsSetGramSchmidtType_C", (TS, TSSundialsGramSchmidtType), (ts, type));
722:   return 0;
723: }

725: /*@
726:    TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
727:                          `TSSUNDIALS` for error control.

729:    Logically Collective

731:    Input Parameters:
732: +    ts  - the time-step context
733: .    aabs - the absolute tolerance
734: -    rel - the relative tolerance

736:      See the CVODE/SUNDIALS users manual for exact details on these parameters. Essentially
737:     these regulate the size of the error for a SINGLE timestep.

739:    Level: intermediate

741: .seealso: [](chapter_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetGMRESMaxl()`,
742:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`,
743:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
744:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`,
745:           `TSSetExactFinalTime()`
746: @*/
747: PetscErrorCode TSSundialsSetTolerance(TS ts, PetscReal aabs, PetscReal rel)
748: {
749:   PetscTryMethod(ts, "TSSundialsSetTolerance_C", (TS, PetscReal, PetscReal), (ts, aabs, rel));
750:   return 0;
751: }

753: /*@
754:    TSSundialsGetPC - Extract the PC context from a time-step context for `TSSUNDIALS`.

756:    Input Parameter:
757: .    ts - the time-step context

759:    Output Parameter:
760: .    pc - the preconditioner context

762:    Level: advanced

764: .seealso: [](chapter_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
765:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
766:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
767:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`
768: @*/
769: PetscErrorCode TSSundialsGetPC(TS ts, PC *pc)
770: {
771:   PetscUseMethod(ts, "TSSundialsGetPC_C", (TS, PC *), (ts, pc));
772:   return 0;
773: }

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

778:    Input Parameters:
779: +   ts - the time-step context
780: -   mindt - lowest time step if positive, negative to deactivate

782:    Note:
783:    `TSSUNDIALS` will error if it is not possible to keep the estimated truncation error below
784:    the tolerance set with `TSSundialsSetTolerance()` without going below this step size.

786:    Level: beginner

788: .seealso: [](chapter_ts), `TSSundialsSetType()`, `TSSundialsSetTolerance()`,
789: @*/
790: PetscErrorCode TSSundialsSetMinTimeStep(TS ts, PetscReal mindt)
791: {
792:   PetscTryMethod(ts, "TSSundialsSetMinTimeStep_C", (TS, PetscReal), (ts, mindt));
793:   return 0;
794: }

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

799:    Input Parameters:
800: +   ts - the time-step context
801: -   maxdt - lowest time step if positive, negative to deactivate

803:    Level: beginner

805: .seealso: [](chapter_ts), `TSSundialsSetType()`, `TSSundialsSetTolerance()`,
806: @*/
807: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts, PetscReal maxdt)
808: {
809:   PetscTryMethod(ts, "TSSundialsSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
810:   return 0;
811: }

813: /*@
814:    TSSundialsMonitorInternalSteps - Monitor `TSSUNDIALS` internal steps (Defaults to false).

816:    Input Parameters:
817: +   ts - the time-step context
818: -   ft - `PETSC_TRUE` if monitor, else `PETSC_FALSE`

820:    Level: beginner

822: .seealso: [](chapter_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
823:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
824:           `TSSundialsGetIterations()`, `TSSundialsSetType()`,
825:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`
826: @*/
827: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts, PetscBool ft)
828: {
829:   PetscTryMethod(ts, "TSSundialsMonitorInternalSteps_C", (TS, PetscBool), (ts, ft));
830:   return 0;
831: }

833: /*@
834:    TSSundialsSetUseDense - Set a flag to use a dense linear solver in `TSSUNDIALS` (serial only)

836:    Logically Collective

838:    Input Parameters:
839: +    ts         - the time-step context
840: -    use_dense  - `PETSC_TRUE` to use the dense solver

842:    Level: advanced

844: .seealso: [](chapter_ts), `TSSUNDIALS`
845: @*/
846: PetscErrorCode TSSundialsSetUseDense(TS ts, PetscBool use_dense)
847: {
849:   PetscTryMethod(ts, "TSSundialsSetUseDense_C", (TS, PetscBool), (ts, use_dense));
850:   return 0;
851: }

853: /* -------------------------------------------------------------------------------------------*/
854: /*MC
855:       TSSUNDIALS - ODE solver using a very old version of the LLNL CVODE/SUNDIALS package, version 2.5 (now called SUNDIALS). Requires ./configure --download-sundials

857:    Options Database Keys:
858: +    -ts_sundials_type <bdf,adams> -
859: .    -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
860: .    -ts_sundials_atol <tol> - Absolute tolerance for convergence
861: .    -ts_sundials_rtol <tol> - Relative tolerance for convergence
862: .    -ts_sundials_linear_tolerance <tol> -
863: .    -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
864: .    -ts_sundials_monitor_steps - Monitor SUNDIALS internal steps
865: -    -ts_sundials_use_dense - Use a dense linear solver within CVODE (serial only)

867:     Level: beginner

869:     Note:
870:     This uses its own nonlinear solver and Krylov method so PETSc `SNES` and `KSP` options do not apply,
871:     only PETSc `PC` options.

873: .seealso: [](chapter_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`, `TSSundialsSetLinearTolerance()`, `TSType`,
874:           `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`, `TSSundialsGetIterations()`, `TSSetExactFinalTime()`
875: M*/
876: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
877: {
878:   TS_Sundials *cvode;
879:   PC           pc;

881:   ts->ops->reset          = TSReset_Sundials;
882:   ts->ops->destroy        = TSDestroy_Sundials;
883:   ts->ops->view           = TSView_Sundials;
884:   ts->ops->setup          = TSSetUp_Sundials;
885:   ts->ops->step           = TSStep_Sundials;
886:   ts->ops->interpolate    = TSInterpolate_Sundials;
887:   ts->ops->setfromoptions = TSSetFromOptions_Sundials;
888:   ts->default_adapt_type  = TSADAPTNONE;

890:   PetscNew(&cvode);

892:   ts->usessnes = PETSC_TRUE;

894:   ts->data           = (void *)cvode;
895:   cvode->cvode_type  = SUNDIALS_BDF;
896:   cvode->gtype       = SUNDIALS_CLASSICAL_GS;
897:   cvode->maxl        = 5;
898:   cvode->maxord      = PETSC_DEFAULT;
899:   cvode->linear_tol  = .05;
900:   cvode->monitorstep = PETSC_TRUE;
901:   cvode->use_dense   = PETSC_FALSE;

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

905:   cvode->mindt = -1.;
906:   cvode->maxdt = -1.;

908:   /* set tolerance for SUNDIALS */
909:   cvode->reltol = 1e-6;
910:   cvode->abstol = 1e-6;

912:   /* set PCNONE as default pctype */
913:   TSSundialsGetPC_Sundials(ts, &pc);
914:   PCSetType(pc, PCNONE);

916:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetType_C", TSSundialsSetType_Sundials);
917:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxl_C", TSSundialsSetMaxl_Sundials);
918:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetLinearTolerance_C", TSSundialsSetLinearTolerance_Sundials);
919:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetGramSchmidtType_C", TSSundialsSetGramSchmidtType_Sundials);
920:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetTolerance_C", TSSundialsSetTolerance_Sundials);
921:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMinTimeStep_C", TSSundialsSetMinTimeStep_Sundials);
922:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxTimeStep_C", TSSundialsSetMaxTimeStep_Sundials);
923:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetPC_C", TSSundialsGetPC_Sundials);
924:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetIterations_C", TSSundialsGetIterations_Sundials);
925:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsMonitorInternalSteps_C", TSSundialsMonitorInternalSteps_Sundials);
926:   PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetUseDense_C", TSSundialsSetUseDense_Sundials);
927:   return 0;
928: }