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