Actual source code: sundials.c
petsc-3.4.5 2014-06-29
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 sundials.h
10: /*
11: TSPrecond_Sundials - function that we provide to SUNDIALS to
12: evaluate the preconditioner.
13: */
16: PetscErrorCode TSPrecond_Sundials(realtype tn,N_Vector y,N_Vector fy,booleantype jok,booleantype *jcurPtr,
17: realtype _gamma,void *P_data,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
18: {
19: TS ts = (TS) P_data;
20: TS_Sundials *cvode = (TS_Sundials*)ts->data;
21: PC pc;
23: Mat J,P;
24: Vec yy = cvode->w1,yydot = cvode->ydot;
25: PetscReal gm = (PetscReal)_gamma;
26: MatStructure str = DIFFERENT_NONZERO_PATTERN;
27: PetscScalar *y_data;
30: TSGetIJacobian(ts,&J,&P,NULL,NULL);
31: y_data = (PetscScalar*) N_VGetArrayPointer(y);
32: VecPlaceArray(yy,y_data);
33: VecZeroEntries(yydot); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */
34: /* compute the shifted Jacobian (1/gm)*I + Jrest */
35: TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,&J,&P,&str,PETSC_FALSE);
36: VecResetArray(yy);
37: MatScale(P,gm); /* turn into I-gm*Jrest, J is not used by Sundials */
38: *jcurPtr = TRUE;
39: TSSundialsGetPC(ts,&pc);
40: PCSetOperators(pc,J,P,str);
41: return(0);
42: }
44: /*
45: TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner.
46: */
49: PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,
50: realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
51: {
52: TS ts = (TS) P_data;
53: TS_Sundials *cvode = (TS_Sundials*)ts->data;
54: PC pc;
55: Vec rr = cvode->w1,zz = cvode->w2;
57: PetscScalar *r_data,*z_data;
60: /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
61: r_data = (PetscScalar*) N_VGetArrayPointer(r);
62: z_data = (PetscScalar*) N_VGetArrayPointer(z);
63: VecPlaceArray(rr,r_data);
64: VecPlaceArray(zz,z_data);
66: /* Solve the Px=r and put the result in zz */
67: TSSundialsGetPC(ts,&pc);
68: PCApply(pc,rr,zz);
69: VecResetArray(rr);
70: VecResetArray(zz);
71: return(0);
72: }
74: /*
75: TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
76: */
79: int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
80: {
81: TS ts = (TS) ctx;
82: DM dm;
83: DMTS tsdm;
84: TSIFunction ifunction;
85: MPI_Comm comm;
86: TS_Sundials *cvode = (TS_Sundials*)ts->data;
87: Vec yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;
88: PetscScalar *y_data,*ydot_data;
92: PetscObjectGetComm((PetscObject)ts,&comm);
93: /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
94: y_data = (PetscScalar*) N_VGetArrayPointer(y);
95: ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);
96: VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
97: VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);
99: /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
100: TSGetDM(ts,&dm);
101: DMGetDMTS(dm,&tsdm);
102: DMTSGetIFunction(dm,&ifunction,NULL);
103: if (!ifunction) {
104: TSComputeRHSFunction(ts,t,yy,yyd);
105: } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
106: VecZeroEntries(yydot);
107: TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr);
108: VecScale(yyd,-1.);
109: }
110: VecResetArray(yy);CHKERRABORT(comm,ierr);
111: VecResetArray(yyd);CHKERRABORT(comm,ierr);
112: return(0);
113: }
115: /*
116: TSStep_Sundials - Calls Sundials to integrate the ODE.
117: */
120: PetscErrorCode TSStep_Sundials(TS ts)
121: {
122: TS_Sundials *cvode = (TS_Sundials*)ts->data;
124: PetscInt flag;
125: long int its,nsteps;
126: realtype t,tout;
127: PetscScalar *y_data;
128: void *mem;
131: mem = cvode->mem;
132: tout = ts->max_time;
133: VecGetArray(ts->vec_sol,&y_data);
134: N_VSetArrayPointer((realtype*)y_data,cvode->y);
135: VecRestoreArray(ts->vec_sol,NULL);
137: TSPreStep(ts);
139: /* We would like to call TSPreStep() when starting each step (including rejections) and TSPreStage() before each
140: * stage solve, but CVode does not appear to support this. */
141: if (cvode->monitorstep) flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
142: else flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);
144: if (flag) { /* display error message */
145: switch (flag) {
146: case CV_ILL_INPUT:
147: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");
148: break;
149: case CV_TOO_CLOSE:
150: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");
151: break;
152: case CV_TOO_MUCH_WORK: {
153: PetscReal tcur;
154: CVodeGetNumSteps(mem,&nsteps);
155: CVodeGetCurrentTime(mem,&tcur);
156: SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_WORK. At t=%G, nsteps %D exceeds mxstep %D. Increase '-ts_max_steps <>' or modify TSSetDuration()",tcur,nsteps,ts->max_steps);
157: } break;
158: case CV_TOO_MUCH_ACC:
159: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");
160: break;
161: case CV_ERR_FAILURE:
162: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");
163: break;
164: case CV_CONV_FAILURE:
165: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");
166: break;
167: case CV_LINIT_FAIL:
168: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");
169: break;
170: case CV_LSETUP_FAIL:
171: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");
172: break;
173: case CV_LSOLVE_FAIL:
174: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");
175: break;
176: case CV_RHSFUNC_FAIL:
177: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");
178: break;
179: case CV_FIRST_RHSFUNC_ERR:
180: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");
181: break;
182: case CV_REPTD_RHSFUNC_ERR:
183: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");
184: break;
185: case CV_UNREC_RHSFUNC_ERR:
186: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");
187: break;
188: case CV_RTFUNC_FAIL:
189: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");
190: break;
191: default:
192: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);
193: }
194: }
196: /* copy the solution from cvode->y to cvode->update and sol */
197: VecPlaceArray(cvode->w1,y_data);
198: VecCopy(cvode->w1,cvode->update);
199: VecResetArray(cvode->w1);
200: VecCopy(cvode->update,ts->vec_sol);
201: CVodeGetNumNonlinSolvIters(mem,&its);
202: CVSpilsGetNumLinIters(mem, &its);
203: ts->snes_its = its; ts->ksp_its = its;
205: ts->time_step = t - ts->ptime;
206: ts->ptime = t;
207: ts->steps++;
209: CVodeGetNumSteps(mem,&nsteps);
210: if (!cvode->monitorstep) ts->steps = nsteps;
211: return(0);
212: }
216: static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
217: {
218: TS_Sundials *cvode = (TS_Sundials*)ts->data;
219: N_Vector y;
221: PetscScalar *x_data;
222: PetscInt glosize,locsize;
225: /* get the vector size */
226: VecGetSize(X,&glosize);
227: VecGetLocalSize(X,&locsize);
229: /* allocate the memory for N_Vec y */
230: y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
231: if (!y) SETERRQ(PETSC_COMM_SELF,1,"Interpolated y is not allocated");
233: VecGetArray(X,&x_data);
234: N_VSetArrayPointer((realtype*)x_data,y);
235: CVodeGetDky(cvode->mem,t,0,y);
236: VecRestoreArray(X,&x_data);
237: return(0);
238: }
242: PetscErrorCode TSReset_Sundials(TS ts)
243: {
244: TS_Sundials *cvode = (TS_Sundials*)ts->data;
248: VecDestroy(&cvode->update);
249: VecDestroy(&cvode->ydot);
250: VecDestroy(&cvode->w1);
251: VecDestroy(&cvode->w2);
252: if (cvode->mem) CVodeFree(&cvode->mem);
253: return(0);
254: }
258: PetscErrorCode TSDestroy_Sundials(TS ts)
259: {
260: TS_Sundials *cvode = (TS_Sundials*)ts->data;
264: TSReset_Sundials(ts);
265: MPI_Comm_free(&(cvode->comm_sundials));
266: PetscFree(ts->data);
267: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);
268: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);
269: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);
270: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);
271: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);
272: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);
273: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);
274: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);
275: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);
276: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",NULL);
277: return(0);
278: }
282: PetscErrorCode TSSetUp_Sundials(TS ts)
283: {
284: TS_Sundials *cvode = (TS_Sundials*)ts->data;
286: PetscInt glosize,locsize,i,flag;
287: PetscScalar *y_data,*parray;
288: void *mem;
289: PC pc;
290: PCType pctype;
291: PetscBool pcnone;
294: /* get the vector size */
295: VecGetSize(ts->vec_sol,&glosize);
296: VecGetLocalSize(ts->vec_sol,&locsize);
298: /* allocate the memory for N_Vec y */
299: cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
300: if (!cvode->y) SETERRQ(PETSC_COMM_SELF,1,"cvode->y is not allocated");
302: /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
303: VecGetArray(ts->vec_sol,&parray);
304: y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);
305: for (i = 0; i < locsize; i++) y_data[i] = parray[i];
306: VecRestoreArray(ts->vec_sol,NULL);
308: VecDuplicate(ts->vec_sol,&cvode->update);
309: VecDuplicate(ts->vec_sol,&cvode->ydot);
310: PetscLogObjectParent(ts,cvode->update);
311: PetscLogObjectParent(ts,cvode->ydot);
313: /*
314: Create work vectors for the TSPSolve_Sundials() routine. Note these are
315: allocated with zero space arrays because the actual array space is provided
316: by Sundials and set using VecPlaceArray().
317: */
318: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w1);
319: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w2);
320: PetscLogObjectParent(ts,cvode->w1);
321: PetscLogObjectParent(ts,cvode->w2);
323: /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */
324: mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
325: if (!mem) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"CVodeCreate() fails");
326: cvode->mem = mem;
328: /* Set the pointer to user-defined data */
329: flag = CVodeSetUserData(mem, ts);
330: if (flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetUserData() fails");
332: /* Sundials may choose to use a smaller initial step, but will never use a larger step. */
333: flag = CVodeSetInitStep(mem,(realtype)ts->time_step);
334: if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetInitStep() failed");
335: if (cvode->mindt > 0) {
336: flag = CVodeSetMinStep(mem,(realtype)cvode->mindt);
337: if (flag) {
338: if (flag == CV_MEM_NULL) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL");
339: else if (flag == CV_ILL_INPUT) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, hmin is nonpositive or it exceeds the maximum allowable step size");
340: else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed");
341: }
342: }
343: if (cvode->maxdt > 0) {
344: flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);
345: if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMaxStep() failed");
346: }
348: /* Call CVodeInit to initialize the integrator memory and specify the
349: * user's right hand side function in u'=f(t,u), the inital time T0, and
350: * the initial dependent variable vector cvode->y */
351: flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y);
352: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag);
354: /* specifies scalar relative and absolute tolerances */
355: flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol);
356: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag);
358: /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
359: flag = CVodeSetMaxNumSteps(mem,ts->max_steps);
361: /* call CVSpgmr to use GMRES as the linear solver. */
362: /* setup the ode integrator with the given preconditioner */
363: TSSundialsGetPC(ts,&pc);
364: PCGetType(pc,&pctype);
365: PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);
366: if (pcnone) {
367: flag = CVSpgmr(mem,PREC_NONE,0);
368: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
369: } else {
370: flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl);
371: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
373: /* Set preconditioner and solve routines Precond and PSolve,
374: and the pointer to the user-defined block data */
375: flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);
376: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag);
377: }
379: flag = CVSpilsSetGSType(mem, MODIFIED_GS);
380: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);
381: return(0);
382: }
384: /* type of CVODE linear multistep method */
385: const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",0};
386: /* type of G-S orthogonalization used by CVODE linear solver */
387: const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",0};
391: PetscErrorCode TSSetFromOptions_Sundials(TS ts)
392: {
393: TS_Sundials *cvode = (TS_Sundials*)ts->data;
395: int indx;
396: PetscBool flag;
397: PC pc;
400: PetscOptionsHead("SUNDIALS ODE solver options");
401: PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);
402: if (flag) {
403: TSSundialsSetType(ts,(TSSundialsLmmType)indx);
404: }
405: PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);
406: if (flag) {
407: TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
408: }
409: PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);
410: PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);
411: PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);
412: PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);
413: PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,&flag);
414: PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,&flag);
415: PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internel steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);
416: PetscOptionsTail();
417: TSSundialsGetPC(ts,&pc);
418: PCSetFromOptions(pc);
419: return(0);
420: }
424: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
425: {
426: TS_Sundials *cvode = (TS_Sundials*)ts->data;
428: char *type;
429: char atype[] = "Adams";
430: char btype[] = "BDF: backward differentiation formula";
431: PetscBool iascii,isstring;
432: long int nsteps,its,nfevals,nlinsetups,nfails,itmp;
433: PetscInt qlast,qcur;
434: PetscReal hinused,hlast,hcur,tcur,tolsfac;
435: PC pc;
438: if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
439: else type = btype;
441: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
442: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
443: if (iascii) {
444: PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");
445: PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);
446: PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);
447: PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",cvode->linear_tol);
448: PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);
449: if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
450: PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
451: } else {
452: PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
453: }
454: if (cvode->mindt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",cvode->mindt);}
455: if (cvode->maxdt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",cvode->maxdt);}
457: /* Outputs from CVODE, CVSPILS */
458: CVodeGetTolScaleFactor(cvode->mem,&tolsfac);
459: PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);
460: CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,
461: &nlinsetups,&nfails,&qlast,&qcur,
462: &hinused,&hlast,&hcur,&tcur);
463: PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);
464: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);
465: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);
466: PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);
468: CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);
469: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);
470: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);
472: CVSpilsGetNumLinIters(cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
473: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);
474: CVSpilsGetNumConvFails(cvode->mem,&itmp);
475: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);
477: TSSundialsGetPC(ts,&pc);
478: PCView(pc,viewer);
479: CVSpilsGetNumPrecEvals(cvode->mem,&itmp);
480: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);
481: CVSpilsGetNumPrecSolves(cvode->mem,&itmp);
482: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);
484: CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);
485: PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);
486: CVSpilsGetNumRhsEvals(cvode->mem,&itmp);
487: PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);
488: } else if (isstring) {
489: PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
490: }
491: return(0);
492: }
495: /* --------------------------------------------------------------------------*/
498: PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
499: {
500: TS_Sundials *cvode = (TS_Sundials*)ts->data;
503: cvode->cvode_type = type;
504: return(0);
505: }
509: PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
510: {
511: TS_Sundials *cvode = (TS_Sundials*)ts->data;
514: cvode->maxl = maxl;
515: return(0);
516: }
520: PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
521: {
522: TS_Sundials *cvode = (TS_Sundials*)ts->data;
525: cvode->linear_tol = tol;
526: return(0);
527: }
531: PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
532: {
533: TS_Sundials *cvode = (TS_Sundials*)ts->data;
536: cvode->gtype = type;
537: return(0);
538: }
542: PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
543: {
544: TS_Sundials *cvode = (TS_Sundials*)ts->data;
547: if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
548: if (rel != PETSC_DECIDE) cvode->reltol = rel;
549: return(0);
550: }
554: PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
555: {
556: TS_Sundials *cvode = (TS_Sundials*)ts->data;
559: cvode->mindt = mindt;
560: return(0);
561: }
565: PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
566: {
567: TS_Sundials *cvode = (TS_Sundials*)ts->data;
570: cvode->maxdt = maxdt;
571: return(0);
572: }
575: PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc)
576: {
577: SNES snes;
578: KSP ksp;
582: TSGetSNES(ts,&snes);
583: SNESGetKSP(snes,&ksp);
584: KSPGetPC(ksp,pc);
585: return(0);
586: }
590: PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
591: {
593: if (nonlin) *nonlin = ts->snes_its;
594: if (lin) *lin = ts->ksp_its;
595: return(0);
596: }
600: PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
601: {
602: TS_Sundials *cvode = (TS_Sundials*)ts->data;
605: cvode->monitorstep = s;
606: return(0);
607: }
608: /* -------------------------------------------------------------------------------------------*/
612: /*@C
613: TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.
615: Not Collective
617: Input parameters:
618: . ts - the time-step context
620: Output Parameters:
621: + nonlin - number of nonlinear iterations
622: - lin - number of linear iterations
624: Level: advanced
626: Notes:
627: These return the number since the creation of the TS object
629: .keywords: non-linear iterations, linear iterations
631: .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
632: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
633: TSSundialsGetIterations(), TSSundialsSetType(),
634: TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()
636: @*/
637: PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
638: {
642: PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));
643: return(0);
644: }
648: /*@
649: TSSundialsSetType - Sets the method that Sundials will use for integration.
651: Logically Collective on TS
653: Input parameters:
654: + ts - the time-step context
655: - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF
657: Level: intermediate
659: .keywords: Adams, backward differentiation formula
661: .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(),
662: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
663: TSSundialsGetIterations(), TSSundialsSetType(),
664: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
665: TSSetExactFinalTime()
666: @*/
667: PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type)
668: {
672: PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));
673: return(0);
674: }
678: /*@
679: TSSundialsSetMaxl - Sets the dimension of the Krylov space used by
680: GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so
681: this is the maximum number of GMRES steps that will be used.
683: Logically Collective on TS
685: Input parameters:
686: + ts - the time-step context
687: - maxl - number of direction vectors (the dimension of Krylov subspace).
689: Level: advanced
691: .keywords: GMRES
693: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
694: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
695: TSSundialsGetIterations(), TSSundialsSetType(),
696: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
697: TSSetExactFinalTime()
699: @*/
700: PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl)
701: {
706: PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));
707: return(0);
708: }
712: /*@
713: TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
714: system by SUNDIALS.
716: Logically Collective on TS
718: Input parameters:
719: + ts - the time-step context
720: - tol - the factor by which the tolerance on the nonlinear solver is
721: multiplied to get the tolerance on the linear solver, .05 by default.
723: Level: advanced
725: .keywords: GMRES, linear convergence tolerance, SUNDIALS
727: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
728: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
729: TSSundialsGetIterations(), TSSundialsSetType(),
730: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
731: TSSetExactFinalTime()
733: @*/
734: PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol)
735: {
740: PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));
741: return(0);
742: }
746: /*@
747: TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
748: in GMRES method by SUNDIALS linear solver.
750: Logically Collective on TS
752: Input parameters:
753: + ts - the time-step context
754: - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS
756: Level: advanced
758: .keywords: Sundials, orthogonalization
760: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
761: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(),
762: TSSundialsGetIterations(), TSSundialsSetType(),
763: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
764: TSSetExactFinalTime()
766: @*/
767: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
768: {
772: PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));
773: return(0);
774: }
778: /*@
779: TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
780: Sundials for error control.
782: Logically Collective on TS
784: Input parameters:
785: + ts - the time-step context
786: . aabs - the absolute tolerance
787: - rel - the relative tolerance
789: See the Cvode/Sundials users manual for exact details on these parameters. Essentially
790: these regulate the size of the error for a SINGLE timestep.
792: Level: intermediate
794: .keywords: Sundials, tolerance
796: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
797: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
798: TSSundialsGetIterations(), TSSundialsSetType(),
799: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
800: TSSetExactFinalTime()
802: @*/
803: PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel)
804: {
808: PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));
809: return(0);
810: }
814: /*@
815: TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.
817: Input Parameter:
818: . ts - the time-step context
820: Output Parameter:
821: . pc - the preconditioner context
823: Level: advanced
825: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
826: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
827: TSSundialsGetIterations(), TSSundialsSetType(),
828: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
829: @*/
830: PetscErrorCode TSSundialsGetPC(TS ts,PC *pc)
831: {
835: PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));
836: return(0);
837: }
841: /*@
842: TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.
844: Input Parameter:
845: + ts - the time-step context
846: - mindt - lowest time step if positive, negative to deactivate
848: Note:
849: Sundials will error if it is not possible to keep the estimated truncation error below
850: the tolerance set with TSSundialsSetTolerance() without going below this step size.
852: Level: beginner
854: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
855: @*/
856: PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
857: {
861: PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));
862: return(0);
863: }
867: /*@
868: TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.
870: Input Parameter:
871: + ts - the time-step context
872: - maxdt - lowest time step if positive, negative to deactivate
874: Level: beginner
876: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
877: @*/
878: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
879: {
883: PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
884: return(0);
885: }
889: /*@
890: TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false).
892: Input Parameter:
893: + ts - the time-step context
894: - ft - PETSC_TRUE if monitor, else PETSC_FALSE
896: Level: beginner
898: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
899: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
900: TSSundialsGetIterations(), TSSundialsSetType(),
901: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
902: @*/
903: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
904: {
908: PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));
909: return(0);
910: }
911: /* -------------------------------------------------------------------------------------------*/
912: /*MC
913: TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)
915: Options Database:
916: + -ts_sundials_type <bdf,adams>
917: . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
918: . -ts_sundials_atol <tol> - Absolute tolerance for convergence
919: . -ts_sundials_rtol <tol> - Relative tolerance for convergence
920: . -ts_sundials_linear_tolerance <tol>
921: . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
922: - -ts_sundials_monitor_steps - Monitor SUNDIALS internel steps
925: Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply
926: only PETSc PC options
928: Level: beginner
930: .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(),
931: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime()
933: M*/
936: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
937: {
938: TS_Sundials *cvode;
940: PC pc;
943: ts->ops->reset = TSReset_Sundials;
944: ts->ops->destroy = TSDestroy_Sundials;
945: ts->ops->view = TSView_Sundials;
946: ts->ops->setup = TSSetUp_Sundials;
947: ts->ops->step = TSStep_Sundials;
948: ts->ops->interpolate = TSInterpolate_Sundials;
949: ts->ops->setfromoptions = TSSetFromOptions_Sundials;
951: PetscNewLog(ts,TS_Sundials,&cvode);
953: ts->data = (void*)cvode;
954: cvode->cvode_type = SUNDIALS_BDF;
955: cvode->gtype = SUNDIALS_CLASSICAL_GS;
956: cvode->maxl = 5;
957: cvode->linear_tol = .05;
958: cvode->monitorstep = PETSC_TRUE;
960: MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));
962: cvode->mindt = -1.;
963: cvode->maxdt = -1.;
965: /* set tolerance for Sundials */
966: cvode->reltol = 1e-6;
967: cvode->abstol = 1e-6;
969: /* set PCNONE as default pctype */
970: TSSundialsGetPC_Sundials(ts,&pc);
971: PCSetType(pc,PCNONE);
973: if (ts->exact_final_time != TS_EXACTFINALTIME_STEPOVER) ts->exact_final_time = TS_EXACTFINALTIME_STEPOVER;
975: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);
976: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);
977: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);
978: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);
979: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);
980: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);
981: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);
982: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);
983: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);
984: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);
985: return(0);
986: }