Actual source code: sundials.c
1: /*
2: Provides a PETSc interface to SUNDIALS/CVODE solver.
3: The interface to PVODE (old version of CVODE) was originally contributed
4: by Liyang Xu. It has been redone by Hong Zhang and Dinesh Kaushik.
6: Reference: sundials-2.4.0/examples/cvode/parallel/cvDiurnal_kry_p.c
7: */
8: #include <../src/ts/impls/implicit/sundials/sundials.h>
10: /*
11: TSPrecond_Sundials - function that we provide to SUNDIALS to
12: evaluate the preconditioner.
13: */
14: PetscErrorCode TSPrecond_Sundials(realtype tn,N_Vector y,N_Vector fy,booleantype jok,booleantype *jcurPtr,
15: realtype _gamma,void *P_data,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
16: {
17: TS ts = (TS) P_data;
18: TS_Sundials *cvode = (TS_Sundials*)ts->data;
19: PC pc;
20: Mat J,P;
21: Vec yy = cvode->w1,yydot = cvode->ydot;
22: PetscReal gm = (PetscReal)_gamma;
23: PetscScalar *y_data;
25: TSGetIJacobian(ts,&J,&P,NULL,NULL);
26: y_data = (PetscScalar*) N_VGetArrayPointer(y);
27: VecPlaceArray(yy,y_data);
28: VecZeroEntries(yydot); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */
29: /* compute the shifted Jacobian (1/gm)*I + Jrest */
30: TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,J,P,PETSC_FALSE);
31: VecResetArray(yy);
32: MatScale(P,gm); /* turn into I-gm*Jrest, J is not used by Sundials */
33: *jcurPtr = TRUE;
34: TSSundialsGetPC(ts,&pc);
35: PCSetOperators(pc,J,P);
36: return 0;
37: }
39: /*
40: TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner.
41: */
42: PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,
43: realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
44: {
45: TS ts = (TS) P_data;
46: TS_Sundials *cvode = (TS_Sundials*)ts->data;
47: PC pc;
48: Vec rr = cvode->w1,zz = cvode->w2;
49: PetscScalar *r_data,*z_data;
51: /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
52: r_data = (PetscScalar*) N_VGetArrayPointer(r);
53: z_data = (PetscScalar*) N_VGetArrayPointer(z);
54: VecPlaceArray(rr,r_data);
55: VecPlaceArray(zz,z_data);
57: /* Solve the Px=r and put the result in zz */
58: TSSundialsGetPC(ts,&pc);
59: PCApply(pc,rr,zz);
60: VecResetArray(rr);
61: VecResetArray(zz);
62: return 0;
63: }
65: /*
66: TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
67: */
68: int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
69: {
70: TS ts = (TS) ctx;
71: DM dm;
72: DMTS tsdm;
73: TSIFunction ifunction;
74: MPI_Comm comm;
75: TS_Sundials *cvode = (TS_Sundials*)ts->data;
76: Vec yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;
77: PetscScalar *y_data,*ydot_data;
79: PetscObjectGetComm((PetscObject)ts,&comm);
80: /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
81: y_data = (PetscScalar*) N_VGetArrayPointer(y);
82: ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);
83: comm,VecPlaceArray(yy,y_data);
84: comm,VecPlaceArray(yyd,ydot_data);
86: /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
87: TSGetDM(ts,&dm);
88: DMGetDMTS(dm,&tsdm);
89: DMTSGetIFunction(dm,&ifunction,NULL);
90: if (!ifunction) {
91: TSComputeRHSFunction(ts,t,yy,yyd);
92: } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
93: VecZeroEntries(yydot);
94: comm,TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);
95: VecScale(yyd,-1.);
96: }
97: comm,VecResetArray(yy);
98: comm,VecResetArray(yyd);
99: return 0;
100: }
102: /*
103: TSStep_Sundials - Calls Sundials to integrate the ODE.
104: */
105: PetscErrorCode TSStep_Sundials(TS ts)
106: {
107: TS_Sundials *cvode = (TS_Sundials*)ts->data;
108: PetscInt flag;
109: long int nits,lits,nsteps;
110: realtype t,tout;
111: PetscScalar *y_data;
112: void *mem;
114: mem = cvode->mem;
115: tout = ts->max_time;
116: VecGetArray(ts->vec_sol,&y_data);
117: N_VSetArrayPointer((realtype*)y_data,cvode->y);
118: VecRestoreArray(ts->vec_sol,NULL);
120: /* We would like to TSPreStage() and TSPostStage()
121: * before each stage solve but CVode does not appear to support this. */
122: if (cvode->monitorstep)
123: flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
124: else
125: flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);
127: if (flag) { /* display error message */
128: switch (flag) {
129: case CV_ILL_INPUT:
130: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");
131: break;
132: case CV_TOO_CLOSE:
133: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");
134: break;
135: case CV_TOO_MUCH_WORK: {
136: PetscReal tcur;
137: CVodeGetNumSteps(mem,&nsteps);
138: CVodeGetCurrentTime(mem,&tcur);
139: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_WORK. At t=%g, nsteps %D exceeds maxstep %D. Increase '-ts_max_steps <>' or modify TSSetMaxSteps()",(double)tcur,nsteps,ts->max_steps);
140: } break;
141: case CV_TOO_MUCH_ACC:
142: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");
143: break;
144: case CV_ERR_FAILURE:
145: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");
146: break;
147: case CV_CONV_FAILURE:
148: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");
149: break;
150: case CV_LINIT_FAIL:
151: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");
152: break;
153: case CV_LSETUP_FAIL:
154: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");
155: break;
156: case CV_LSOLVE_FAIL:
157: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");
158: break;
159: case CV_RHSFUNC_FAIL:
160: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");
161: break;
162: case CV_FIRST_RHSFUNC_ERR:
163: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");
164: break;
165: case CV_REPTD_RHSFUNC_ERR:
166: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");
167: break;
168: case CV_UNREC_RHSFUNC_ERR:
169: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");
170: break;
171: case CV_RTFUNC_FAIL:
172: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");
173: break;
174: default:
175: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);
176: }
177: }
179: /* log inner nonlinear and linear iterations */
180: CVodeGetNumNonlinSolvIters(mem,&nits);
181: CVSpilsGetNumLinIters(mem,&lits);
182: ts->snes_its += nits; ts->ksp_its = lits;
184: /* copy the solution from cvode->y to cvode->update and sol */
185: VecPlaceArray(cvode->w1,y_data);
186: VecCopy(cvode->w1,cvode->update);
187: VecResetArray(cvode->w1);
188: VecCopy(cvode->update,ts->vec_sol);
190: ts->time_step = t - ts->ptime;
191: ts->ptime = t;
193: CVodeGetNumSteps(mem,&nsteps);
194: if (!cvode->monitorstep) ts->steps += nsteps - 1; /* TSStep() increments the step counter by one */
195: return 0;
196: }
198: static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
199: {
200: TS_Sundials *cvode = (TS_Sundials*)ts->data;
201: N_Vector y;
202: PetscScalar *x_data;
203: PetscInt glosize,locsize;
205: /* get the vector size */
206: VecGetSize(X,&glosize);
207: VecGetLocalSize(X,&locsize);
208: VecGetArray(X,&x_data);
210: /* Initialize N_Vec y with x_data */
211: if (cvode->use_dense) {
212: PetscMPIInt size;
214: MPI_Comm_size(PETSC_COMM_WORLD,&size);
216: y = N_VMake_Serial(locsize,(realtype*)x_data);
217: } else {
218: y = N_VMake_Parallel(cvode->comm_sundials,locsize,glosize,(realtype*)x_data);
219: }
223: CVodeGetDky(cvode->mem,t,0,y);
224: VecRestoreArray(X,&x_data);
225: return 0;
226: }
228: PetscErrorCode TSReset_Sundials(TS ts)
229: {
230: TS_Sundials *cvode = (TS_Sundials*)ts->data;
232: VecDestroy(&cvode->update);
233: VecDestroy(&cvode->ydot);
234: VecDestroy(&cvode->w1);
235: VecDestroy(&cvode->w2);
236: if (cvode->mem) CVodeFree(&cvode->mem);
237: return 0;
238: }
240: PetscErrorCode TSDestroy_Sundials(TS ts)
241: {
242: TS_Sundials *cvode = (TS_Sundials*)ts->data;
244: TSReset_Sundials(ts);
245: MPI_Comm_free(&(cvode->comm_sundials));
246: PetscFree(ts->data);
247: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);
248: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);
249: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);
250: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);
251: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);
252: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);
253: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);
254: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);
255: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);
256: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",NULL);
257: return 0;
258: }
260: PetscErrorCode TSSetUp_Sundials(TS ts)
261: {
262: TS_Sundials *cvode = (TS_Sundials*)ts->data;
263: PetscInt glosize,locsize,i,flag;
264: PetscScalar *y_data,*parray;
265: void *mem;
266: PC pc;
267: PCType pctype;
268: PetscBool pcnone;
272: /* get the vector size */
273: VecGetSize(ts->vec_sol,&glosize);
274: VecGetLocalSize(ts->vec_sol,&locsize);
276: /* allocate the memory for N_Vec y */
277: if (cvode->use_dense) {
278: PetscMPIInt size;
280: MPI_Comm_size(PETSC_COMM_WORLD,&size);
282: cvode->y = N_VNew_Serial(locsize);
283: } else {
284: cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
285: }
288: /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
289: VecGetArray(ts->vec_sol,&parray);
290: y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);
291: for (i = 0; i < locsize; i++) y_data[i] = parray[i];
292: VecRestoreArray(ts->vec_sol,NULL);
294: VecDuplicate(ts->vec_sol,&cvode->update);
295: VecDuplicate(ts->vec_sol,&cvode->ydot);
296: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);
297: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);
299: /*
300: Create work vectors for the TSPSolve_Sundials() routine. Note these are
301: allocated with zero space arrays because the actual array space is provided
302: by Sundials and set using VecPlaceArray().
303: */
304: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,NULL,&cvode->w1);
305: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,NULL,&cvode->w2);
306: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);
307: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);
309: /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */
310: mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
312: cvode->mem = mem;
314: /* Set the pointer to user-defined data */
315: flag = CVodeSetUserData(mem, ts);
318: /* Sundials may choose to use a smaller initial step, but will never use a larger step. */
319: flag = CVodeSetInitStep(mem,(realtype)ts->time_step);
321: if (cvode->mindt > 0) {
322: flag = CVodeSetMinStep(mem,(realtype)cvode->mindt);
323: if (flag) {
326: else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed");
327: }
328: }
329: if (cvode->maxdt > 0) {
330: flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);
332: }
334: /* Call CVodeInit to initialize the integrator memory and specify the
335: * user's right hand side function in u'=f(t,u), the initial time T0, and
336: * the initial dependent variable vector cvode->y */
337: flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y);
340: /* specifies scalar relative and absolute tolerances */
341: flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol);
344: /* Specify max order of BDF / ADAMS method */
345: if (cvode->maxord != PETSC_DEFAULT) {
346: flag = CVodeSetMaxOrd(mem,cvode->maxord);
348: }
350: /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
351: flag = CVodeSetMaxNumSteps(mem,ts->max_steps);
354: if (cvode->use_dense) {
355: /* call CVDense to use a dense linear solver. */
356: PetscMPIInt size;
358: MPI_Comm_size(PETSC_COMM_WORLD,&size);
360: flag = CVDense(mem,locsize);
362: } else {
363: /* call CVSpgmr to use GMRES as the linear solver. */
364: /* setup the ode integrator with the given preconditioner */
365: TSSundialsGetPC(ts,&pc);
366: PCGetType(pc,&pctype);
367: PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);
368: if (pcnone) {
369: flag = CVSpgmr(mem,PREC_NONE,0);
371: } else {
372: flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl);
375: /* Set preconditioner and solve routines Precond and PSolve,
376: and the pointer to the user-defined block data */
377: flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);
379: }
380: }
381: return 0;
382: }
384: /* type of CVODE linear multistep method */
385: const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",NULL};
386: /* type of G-S orthogonalization used by CVODE linear solver */
387: const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",NULL};
389: PetscErrorCode TSSetFromOptions_Sundials(PetscOptionItems *PetscOptionsObject,TS ts)
390: {
391: TS_Sundials *cvode = (TS_Sundials*)ts->data;
392: int indx;
393: PetscBool flag;
394: PC pc;
396: PetscOptionsHead(PetscOptionsObject,"SUNDIALS ODE solver options");
397: PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);
398: if (flag) {
399: TSSundialsSetType(ts,(TSSundialsLmmType)indx);
400: }
401: PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);
402: if (flag) {
403: TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
404: }
405: PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);
406: PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);
407: PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);
408: PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);
409: PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,NULL);
410: PetscOptionsInt("-ts_sundials_maxord","Max Order for BDF/Adams method","TSSundialsSetMaxOrd",cvode->maxord,&cvode->maxord,NULL);
411: PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,NULL);
412: PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internal steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);
413: PetscOptionsBool("-ts_sundials_use_dense","Use dense internal solver in SUNDIALS (serial only)","TSSundialsSetUseDense",cvode->use_dense,&cvode->use_dense,NULL);
414: PetscOptionsTail();
415: TSSundialsGetPC(ts,&pc);
416: PCSetFromOptions(pc);
417: return 0;
418: }
420: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
421: {
422: TS_Sundials *cvode = (TS_Sundials*)ts->data;
423: char *type;
424: char atype[] = "Adams";
425: char btype[] = "BDF: backward differentiation formula";
426: PetscBool iascii,isstring;
427: long int nsteps,its,nfevals,nlinsetups,nfails,itmp;
428: PetscInt qlast,qcur;
429: PetscReal hinused,hlast,hcur,tcur,tolsfac;
430: PC pc;
432: if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
433: else type = btype;
435: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
436: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
437: if (iascii) {
438: PetscViewerASCIIPrintf(viewer,"Sundials integrator does not use SNES!\n");
439: PetscViewerASCIIPrintf(viewer,"Sundials integrator type %s\n",type);
440: PetscViewerASCIIPrintf(viewer,"Sundials maxord %D\n",cvode->maxord);
441: PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",(double)cvode->abstol,(double)cvode->reltol);
442: if (cvode->use_dense) {
443: PetscViewerASCIIPrintf(viewer,"Sundials integrator using a dense linear solve\n");
444: } else {
445: PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",(double)cvode->linear_tol);
446: PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);
447: if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
448: PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
449: } else {
450: PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
451: }
452: }
453: if (cvode->mindt > 0) PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",(double)cvode->mindt);
454: if (cvode->maxdt > 0) PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",(double)cvode->maxdt);
456: /* Outputs from CVODE, CVSPILS */
457: CVodeGetTolScaleFactor(cvode->mem,&tolsfac);
458: PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);
459: CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,&nlinsetups,&nfails,&qlast,&qcur,&hinused,&hlast,&hcur,&tcur);
460: PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);
461: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);
462: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);
463: PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);
465: CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);
466: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);
467: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);
468: if (!cvode->use_dense) {
469: CVSpilsGetNumLinIters(cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
470: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);
471: CVSpilsGetNumConvFails(cvode->mem,&itmp);
472: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);
474: TSSundialsGetPC(ts,&pc);
475: PCView(pc,viewer);
476: CVSpilsGetNumPrecEvals(cvode->mem,&itmp);
477: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);
478: CVSpilsGetNumPrecSolves(cvode->mem,&itmp);
479: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);
480: }
481: CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);
482: PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);
483: CVSpilsGetNumRhsEvals(cvode->mem,&itmp);
484: PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);
485: } else if (isstring) {
486: PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
487: }
488: return 0;
489: }
491: /* --------------------------------------------------------------------------*/
492: PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
493: {
494: TS_Sundials *cvode = (TS_Sundials*)ts->data;
496: cvode->cvode_type = type;
497: return 0;
498: }
500: PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
501: {
502: TS_Sundials *cvode = (TS_Sundials*)ts->data;
504: cvode->maxl = maxl;
505: return 0;
506: }
508: PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,PetscReal tol)
509: {
510: TS_Sundials *cvode = (TS_Sundials*)ts->data;
512: cvode->linear_tol = tol;
513: return 0;
514: }
516: PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
517: {
518: TS_Sundials *cvode = (TS_Sundials*)ts->data;
520: cvode->gtype = type;
521: return 0;
522: }
524: PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,PetscReal aabs,PetscReal rel)
525: {
526: TS_Sundials *cvode = (TS_Sundials*)ts->data;
528: if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
529: if (rel != PETSC_DECIDE) cvode->reltol = rel;
530: return 0;
531: }
533: PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
534: {
535: TS_Sundials *cvode = (TS_Sundials*)ts->data;
537: cvode->mindt = mindt;
538: return 0;
539: }
541: PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
542: {
543: TS_Sundials *cvode = (TS_Sundials*)ts->data;
545: cvode->maxdt = maxdt;
546: return 0;
547: }
549: PetscErrorCode TSSundialsSetUseDense_Sundials(TS ts,PetscBool use_dense)
550: {
551: TS_Sundials *cvode = (TS_Sundials*)ts->data;
553: cvode->use_dense = use_dense;
554: return 0;
555: }
557: PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc)
558: {
559: SNES snes;
560: KSP ksp;
562: TSGetSNES(ts,&snes);
563: SNESGetKSP(snes,&ksp);
564: KSPGetPC(ksp,pc);
565: return 0;
566: }
568: PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
569: {
570: if (nonlin) *nonlin = ts->snes_its;
571: if (lin) *lin = ts->ksp_its;
572: return 0;
573: }
575: PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
576: {
577: TS_Sundials *cvode = (TS_Sundials*)ts->data;
579: cvode->monitorstep = s;
580: return 0;
581: }
582: /* -------------------------------------------------------------------------------------------*/
584: /*@C
585: TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.
587: Not Collective
589: Input Parameter:
590: . ts - the time-step context
592: Output Parameters:
593: + nonlin - number of nonlinear iterations
594: - lin - number of linear iterations
596: Level: advanced
598: Notes:
599: These return the number since the creation of the TS object
601: .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
602: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
603: TSSundialsGetIterations(), TSSundialsSetType(),
604: TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()
606: @*/
607: PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
608: {
609: PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));
610: return 0;
611: }
613: /*@
614: TSSundialsSetType - Sets the method that Sundials will use for integration.
616: Logically Collective on TS
618: Input Parameters:
619: + ts - the time-step context
620: - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF
622: Level: intermediate
624: .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(),
625: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
626: TSSundialsGetIterations(), TSSundialsSetType(),
627: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
628: TSSetExactFinalTime()
629: @*/
630: PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type)
631: {
632: PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));
633: return 0;
634: }
636: /*@
637: TSSundialsSetMaxord - Sets the maximum order for BDF/Adams method used by SUNDIALS.
639: Logically Collective on TS
641: Input Parameters:
642: + ts - the time-step context
643: - maxord - maximum order of BDF / Adams method
645: Level: advanced
647: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
648: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
649: TSSundialsGetIterations(), TSSundialsSetType(),
650: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
651: TSSetExactFinalTime()
653: @*/
654: PetscErrorCode TSSundialsSetMaxord(TS ts,PetscInt maxord)
655: {
657: PetscTryMethod(ts,"TSSundialsSetMaxOrd_C",(TS,PetscInt),(ts,maxord));
658: return 0;
659: }
661: /*@
662: TSSundialsSetMaxl - Sets the dimension of the Krylov space used by
663: GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so
664: this is the maximum number of GMRES steps that will be used.
666: Logically Collective on TS
668: Input Parameters:
669: + ts - the time-step context
670: - maxl - number of direction vectors (the dimension of Krylov subspace).
672: Level: advanced
674: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
675: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
676: TSSundialsGetIterations(), TSSundialsSetType(),
677: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
678: TSSetExactFinalTime()
680: @*/
681: PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl)
682: {
684: PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));
685: return 0;
686: }
688: /*@
689: TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
690: system by SUNDIALS.
692: Logically Collective on TS
694: Input Parameters:
695: + ts - the time-step context
696: - tol - the factor by which the tolerance on the nonlinear solver is
697: multiplied to get the tolerance on the linear solver, .05 by default.
699: Level: advanced
701: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
702: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
703: TSSundialsGetIterations(), TSSundialsSetType(),
704: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
705: TSSetExactFinalTime()
707: @*/
708: PetscErrorCode TSSundialsSetLinearTolerance(TS ts,PetscReal tol)
709: {
711: PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,PetscReal),(ts,tol));
712: return 0;
713: }
715: /*@
716: TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
717: in GMRES method by SUNDIALS linear solver.
719: Logically Collective on TS
721: Input Parameters:
722: + ts - the time-step context
723: - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS
725: Level: advanced
727: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
728: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(),
729: TSSundialsGetIterations(), TSSundialsSetType(),
730: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
731: TSSetExactFinalTime()
733: @*/
734: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
735: {
736: PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));
737: return 0;
738: }
740: /*@
741: TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
742: Sundials for error control.
744: Logically Collective on TS
746: Input Parameters:
747: + ts - the time-step context
748: . aabs - the absolute tolerance
749: - rel - the relative tolerance
751: See the Cvode/Sundials users manual for exact details on these parameters. Essentially
752: these regulate the size of the error for a SINGLE timestep.
754: Level: intermediate
756: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
757: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
758: TSSundialsGetIterations(), TSSundialsSetType(),
759: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
760: TSSetExactFinalTime()
762: @*/
763: PetscErrorCode TSSundialsSetTolerance(TS ts,PetscReal aabs,PetscReal rel)
764: {
765: PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,PetscReal,PetscReal),(ts,aabs,rel));
766: return 0;
767: }
769: /*@
770: TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.
772: Input Parameter:
773: . ts - the time-step context
775: Output Parameter:
776: . pc - the preconditioner context
778: Level: advanced
780: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
781: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
782: TSSundialsGetIterations(), TSSundialsSetType(),
783: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
784: @*/
785: PetscErrorCode TSSundialsGetPC(TS ts,PC *pc)
786: {
787: PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));
788: return 0;
789: }
791: /*@
792: TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.
794: Input Parameters:
795: + ts - the time-step context
796: - mindt - lowest time step if positive, negative to deactivate
798: Note:
799: Sundials will error if it is not possible to keep the estimated truncation error below
800: the tolerance set with TSSundialsSetTolerance() without going below this step size.
802: Level: beginner
804: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
805: @*/
806: PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
807: {
808: PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));
809: return 0;
810: }
812: /*@
813: TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.
815: Input Parameters:
816: + ts - the time-step context
817: - maxdt - lowest time step if positive, negative to deactivate
819: Level: beginner
821: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
822: @*/
823: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
824: {
825: PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
826: return 0;
827: }
829: /*@
830: TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false).
832: Input Parameters:
833: + ts - the time-step context
834: - ft - PETSC_TRUE if monitor, else PETSC_FALSE
836: Level: beginner
838: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
839: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
840: TSSundialsGetIterations(), TSSundialsSetType(),
841: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
842: @*/
843: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
844: {
845: PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));
846: return 0;
847: }
849: /*@
850: TSSundialsSetUseDense - Set a flag to use a dense linear solver in SUNDIALS (serial only)
852: Logically Collective
854: Input Parameters:
855: + ts - the time-step context
856: - use_dense - PETSC_TRUE to use the dense solver
858: Level: advanced
860: .seealso: TSSUNDIALS
862: @*/
863: PetscErrorCode TSSundialsSetUseDense(TS ts,PetscBool use_dense)
864: {
866: PetscTryMethod(ts,"TSSundialsSetUseDense_C",(TS,PetscBool),(ts,use_dense));
867: return 0;
868: }
870: /* -------------------------------------------------------------------------------------------*/
871: /*MC
872: TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)
874: Options Database:
875: + -ts_sundials_type <bdf,adams> -
876: . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
877: . -ts_sundials_atol <tol> - Absolute tolerance for convergence
878: . -ts_sundials_rtol <tol> - Relative tolerance for convergence
879: . -ts_sundials_linear_tolerance <tol> -
880: . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
881: . -ts_sundials_monitor_steps - Monitor SUNDIALS internal steps
882: - -ts_sundials_use_dense - Use a dense linear solver within CVODE (serial only)
884: Notes:
885: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply,
886: only PETSc PC options.
888: Level: beginner
890: .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(),
891: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime()
893: M*/
894: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
895: {
896: TS_Sundials *cvode;
897: PC pc;
899: ts->ops->reset = TSReset_Sundials;
900: ts->ops->destroy = TSDestroy_Sundials;
901: ts->ops->view = TSView_Sundials;
902: ts->ops->setup = TSSetUp_Sundials;
903: ts->ops->step = TSStep_Sundials;
904: ts->ops->interpolate = TSInterpolate_Sundials;
905: ts->ops->setfromoptions = TSSetFromOptions_Sundials;
906: ts->default_adapt_type = TSADAPTNONE;
908: PetscNewLog(ts,&cvode);
910: ts->usessnes = PETSC_TRUE;
912: ts->data = (void*)cvode;
913: cvode->cvode_type = SUNDIALS_BDF;
914: cvode->gtype = SUNDIALS_CLASSICAL_GS;
915: cvode->maxl = 5;
916: cvode->maxord = PETSC_DEFAULT;
917: cvode->linear_tol = .05;
918: cvode->monitorstep = PETSC_TRUE;
919: cvode->use_dense = PETSC_FALSE;
921: MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));
923: cvode->mindt = -1.;
924: cvode->maxdt = -1.;
926: /* set tolerance for Sundials */
927: cvode->reltol = 1e-6;
928: cvode->abstol = 1e-6;
930: /* set PCNONE as default pctype */
931: TSSundialsGetPC_Sundials(ts,&pc);
932: PCSetType(pc,PCNONE);
934: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);
935: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);
936: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);
937: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);
938: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);
939: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);
940: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);
941: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);
942: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);
943: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);
944: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetUseDense_C",TSSundialsSetUseDense_Sundials);
945: return 0;
946: }