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: }