Actual source code: sundials.c
petsc-3.3-p7 2013-05-11
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,
17: booleantype jok,booleantype *jcurPtr,
18: realtype _gamma,void *P_data,
19: N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
20: {
21: TS ts = (TS) P_data;
22: TS_Sundials *cvode = (TS_Sundials*)ts->data;
23: PC pc;
25: Mat J,P;
26: Vec yy = cvode->w1,yydot = cvode->ydot;
27: PetscReal gm = (PetscReal)_gamma;
28: MatStructure str = DIFFERENT_NONZERO_PATTERN;
29: PetscScalar *y_data;
32: TSGetIJacobian(ts,&J,&P,PETSC_NULL,PETSC_NULL);
33: y_data = (PetscScalar *) N_VGetArrayPointer(y);
34: VecPlaceArray(yy,y_data);
35: VecZeroEntries(yydot); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */
36: /* compute the shifted Jacobian (1/gm)*I + Jrest */
37: TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,&J,&P,&str,PETSC_FALSE);
38: VecResetArray(yy);
39: MatScale(P,gm); /* turn into I-gm*Jrest, J is not used by Sundials */
40: *jcurPtr = TRUE;
41: TSSundialsGetPC(ts,&pc);
42: PCSetOperators(pc,J,P,str);
43: return(0);
44: }
46: /*
47: TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner.
48: */
51: PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,
52: realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
53: {
54: TS ts = (TS) P_data;
55: TS_Sundials *cvode = (TS_Sundials*)ts->data;
56: PC pc;
57: Vec rr = cvode->w1,zz = cvode->w2;
58: PetscErrorCode ierr;
59: PetscScalar *r_data,*z_data;
62: /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
63: r_data = (PetscScalar *) N_VGetArrayPointer(r);
64: z_data = (PetscScalar *) N_VGetArrayPointer(z);
65: VecPlaceArray(rr,r_data);
66: VecPlaceArray(zz,z_data);
68: /* Solve the Px=r and put the result in zz */
69: TSSundialsGetPC(ts,&pc);
70: PCApply(pc,rr,zz);
71: VecResetArray(rr);
72: VecResetArray(zz);
73: return(0);
74: }
76: /*
77: TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
78: */
81: int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
82: {
83: TS ts = (TS) ctx;
84: MPI_Comm comm = ((PetscObject)ts)->comm;
85: TS_Sundials *cvode = (TS_Sundials*)ts->data;
86: Vec yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;
87: PetscScalar *y_data,*ydot_data;
88: PetscErrorCode ierr;
91: /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
92: y_data = (PetscScalar *) N_VGetArrayPointer(y);
93: ydot_data = (PetscScalar *) N_VGetArrayPointer(ydot);
94: VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
95: VecPlaceArray(yyd,ydot_data); CHKERRABORT(comm,ierr);
97: /* now compute the right hand side function */
98: if (!ts->userops->ifunction) {
99: TSComputeRHSFunction(ts,t,yy,yyd);
100: } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
101: VecZeroEntries(yydot);
102: TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE); CHKERRABORT(comm,ierr);
103: VecScale(yyd,-1.);
104: }
105: VecResetArray(yy); CHKERRABORT(comm,ierr);
106: VecResetArray(yyd); CHKERRABORT(comm,ierr);
107: return(0);
108: }
110: /*
111: TSStep_Sundials - Calls Sundials to integrate the ODE.
112: */
115: PetscErrorCode TSStep_Sundials(TS ts)
116: {
117: TS_Sundials *cvode = (TS_Sundials*)ts->data;
119: PetscInt flag;
120: long int its,nsteps;
121: realtype t,tout;
122: PetscScalar *y_data;
123: void *mem;
126: mem = cvode->mem;
127: tout = ts->max_time;
128: VecGetArray(ts->vec_sol,&y_data);
129: N_VSetArrayPointer((realtype *)y_data,cvode->y);
130: VecRestoreArray(ts->vec_sol,PETSC_NULL);
132: TSPreStep(ts);
134: if (cvode->monitorstep) {
135: /* We would like to call TSPreStep() when starting each step (including rejections) and TSPreStage() before each
136: * stage solve, but CVode does not appear to support this. */
137: flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
138: } else {
139: flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);
140: }
142: if (flag){ /* display error message */
143: switch (flag){
144: case CV_ILL_INPUT:
145: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");
146: break;
147: case CV_TOO_CLOSE:
148: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");
149: break;
150: case CV_TOO_MUCH_WORK: {
151: PetscReal tcur;
152: CVodeGetNumSteps(mem,&nsteps);
153: CVodeGetCurrentTime(mem,&tcur);
154: 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);
155: } break;
156: case CV_TOO_MUCH_ACC:
157: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");
158: break;
159: case CV_ERR_FAILURE:
160: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");
161: break;
162: case CV_CONV_FAILURE:
163: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");
164: break;
165: case CV_LINIT_FAIL:
166: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");
167: break;
168: case CV_LSETUP_FAIL:
169: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");
170: break;
171: case CV_LSOLVE_FAIL:
172: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");
173: break;
174: case CV_RHSFUNC_FAIL:
175: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");
176: break;
177: case CV_FIRST_RHSFUNC_ERR:
178: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");
179: break;
180: case CV_REPTD_RHSFUNC_ERR:
181: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");
182: break;
183: case CV_UNREC_RHSFUNC_ERR:
184: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");
185: break;
186: case CV_RTFUNC_FAIL:
187: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");
188: break;
189: default:
190: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);
191: }
192: }
194: /* copy the solution from cvode->y to cvode->update and sol */
195: VecPlaceArray(cvode->w1,y_data);
196: VecCopy(cvode->w1,cvode->update);
197: VecResetArray(cvode->w1);
198: VecCopy(cvode->update,ts->vec_sol);
199: CVodeGetNumNonlinSolvIters(mem,&its);
200: CVSpilsGetNumLinIters(mem, &its);
201: ts->snes_its = its; ts->ksp_its = its;
203: ts->time_step = t - ts->ptime;
204: ts->ptime = t;
205: ts->steps++;
207: CVodeGetNumSteps(mem,&nsteps);
208: if (!cvode->monitorstep) ts->steps = nsteps;
209: return(0);
210: }
214: static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
215: {
216: TS_Sundials *cvode = (TS_Sundials*)ts->data;
217: N_Vector y;
218: PetscErrorCode ierr;
219: PetscScalar *x_data;
220: PetscInt glosize,locsize;
224: /* get the vector size */
225: VecGetSize(X,&glosize);
226: VecGetLocalSize(X,&locsize);
228: /* allocate the memory for N_Vec y */
229: y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
230: if (!y) SETERRQ(PETSC_COMM_SELF,1,"Interpolated y is not allocated");
232: VecGetArray(X,&x_data);
233: N_VSetArrayPointer((realtype *)x_data,y);
234: CVodeGetDky(cvode->mem,t,0,y);
235: 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: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetType_C","",PETSC_NULL);
268: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxl_C","",PETSC_NULL);
269: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetLinearTolerance_C","",PETSC_NULL);
270: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetGramSchmidtType_C","",PETSC_NULL);
271: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetTolerance_C","",PETSC_NULL);
272: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMinTimeStep_C","",PETSC_NULL);
273: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxTimeStep_C","",PETSC_NULL);
274: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetPC_C","",PETSC_NULL);
275: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetIterations_C","",PETSC_NULL);
276: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsMonitorInternalSteps_C","",PETSC_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: const 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,PETSC_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(((PetscObject)ts)->comm,1,locsize,PETSC_DECIDE,0,&cvode->w1);
319: VecCreateMPIWithArray(((PetscObject)ts)->comm,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(((PetscObject)ts)->comm,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){
339: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL");
340: } else if (flag == CV_ILL_INPUT){
341: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMinStep() failed, hmin is nonpositive or it exceeds the maximum allowable step size");
342: } else {
343: SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMinStep() failed");
344: }
345: }
346: }
347: if (cvode->maxdt > 0) {
348: flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);
349: if (flag) SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_LIB,"CVodeSetMaxStep() failed");
350: }
352: /* Call CVodeInit to initialize the integrator memory and specify the
353: * user's right hand side function in u'=f(t,u), the inital time T0, and
354: * the initial dependent variable vector cvode->y */
355: flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y);
356: if (flag){
357: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag);
358: }
360: /* specifies scalar relative and absolute tolerances */
361: flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol);
362: if (flag){
363: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag);
364: }
366: /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
367: flag = CVodeSetMaxNumSteps(mem,ts->max_steps);
369: /* call CVSpgmr to use GMRES as the linear solver. */
370: /* setup the ode integrator with the given preconditioner */
371: TSSundialsGetPC(ts,&pc);
372: PCGetType(pc,&pctype);
373: PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);
374: if (pcnone){
375: flag = CVSpgmr(mem,PREC_NONE,0);
376: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
377: } else {
378: flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl);
379: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
381: /* Set preconditioner and solve routines Precond and PSolve,
382: and the pointer to the user-defined block data */
383: flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);
384: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag);
385: }
387: flag = CVSpilsSetGSType(mem, MODIFIED_GS);
388: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);
389: return(0);
390: }
392: /* type of CVODE linear multistep method */
393: const char *TSSundialsLmmTypes[] = {"","adams","bdf","TSSundialsLmmType","SUNDIALS_",0};
394: /* type of G-S orthogonalization used by CVODE linear solver */
395: const char *TSSundialsGramSchmidtTypes[] = {"","modified","classical","TSSundialsGramSchmidtType","SUNDIALS_",0};
399: PetscErrorCode TSSetFromOptions_Sundials(TS ts)
400: {
401: TS_Sundials *cvode = (TS_Sundials*)ts->data;
403: int indx;
404: PetscBool flag;
405: PC pc;
408: PetscOptionsHead("SUNDIALS ODE solver options");
409: PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);
410: if (flag) {
411: TSSundialsSetType(ts,(TSSundialsLmmType)indx);
412: }
413: PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);
414: if (flag) {
415: TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
416: }
417: PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,PETSC_NULL);
418: PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,PETSC_NULL);
419: PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,PETSC_NULL);
420: PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,PETSC_NULL);
421: PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,&flag);
422: PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,&flag);
423: PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internel steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,PETSC_NULL);
424: PetscOptionsTail();
425: TSSundialsGetPC(ts,&pc);
426: PCSetFromOptions(pc);
427: return(0);
428: }
432: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
433: {
434: TS_Sundials *cvode = (TS_Sundials*)ts->data;
436: char *type;
437: char atype[] = "Adams";
438: char btype[] = "BDF: backward differentiation formula";
439: PetscBool iascii,isstring;
440: long int nsteps,its,nfevals,nlinsetups,nfails,itmp;
441: PetscInt qlast,qcur;
442: PetscReal hinused,hlast,hcur,tcur,tolsfac;
443: PC pc;
446: if (cvode->cvode_type == SUNDIALS_ADAMS) {type = atype;}
447: else {type = btype;}
449: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
450: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
451: if (iascii) {
452: PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");
453: PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);
454: PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);
455: PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",cvode->linear_tol);
456: PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);
457: if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
458: PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
459: } else {
460: PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
461: }
462: if (cvode->mindt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",cvode->mindt);}
463: if (cvode->maxdt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",cvode->maxdt);}
465: /* Outputs from CVODE, CVSPILS */
466: CVodeGetTolScaleFactor(cvode->mem,&tolsfac);
467: PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);
468: CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,
469: &nlinsetups,&nfails,&qlast,&qcur,
470: &hinused,&hlast,&hcur,&tcur);
471: PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);
472: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);
473: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);
474: PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);
476: CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);
477: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);
478: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);
480: CVSpilsGetNumLinIters(cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
481: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);
482: CVSpilsGetNumConvFails(cvode->mem,&itmp);
483: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);
485: TSSundialsGetPC(ts,&pc);
486: PCView(pc,viewer);
487: CVSpilsGetNumPrecEvals(cvode->mem,&itmp);
488: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);
489: CVSpilsGetNumPrecSolves(cvode->mem,&itmp);
490: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);
492: CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);
493: PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);
494: CVSpilsGetNumRhsEvals(cvode->mem,&itmp);
495: PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);
496: } else if (isstring) {
497: PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
498: }
499: return(0);
500: }
503: /* --------------------------------------------------------------------------*/
504: EXTERN_C_BEGIN
507: PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
508: {
509: TS_Sundials *cvode = (TS_Sundials*)ts->data;
512: cvode->cvode_type = type;
513: return(0);
514: }
515: EXTERN_C_END
517: EXTERN_C_BEGIN
520: PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
521: {
522: TS_Sundials *cvode = (TS_Sundials*)ts->data;
525: cvode->maxl = maxl;
526: return(0);
527: }
528: EXTERN_C_END
530: EXTERN_C_BEGIN
533: PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
534: {
535: TS_Sundials *cvode = (TS_Sundials*)ts->data;
538: cvode->linear_tol = tol;
539: return(0);
540: }
541: EXTERN_C_END
543: EXTERN_C_BEGIN
546: PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
547: {
548: TS_Sundials *cvode = (TS_Sundials*)ts->data;
551: cvode->gtype = type;
552: return(0);
553: }
554: EXTERN_C_END
556: EXTERN_C_BEGIN
559: PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
560: {
561: TS_Sundials *cvode = (TS_Sundials*)ts->data;
564: if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
565: if (rel != PETSC_DECIDE) cvode->reltol = rel;
566: return(0);
567: }
568: EXTERN_C_END
570: EXTERN_C_BEGIN
573: PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
574: {
575: TS_Sundials *cvode = (TS_Sundials*)ts->data;
578: cvode->mindt = mindt;
579: return(0);
580: }
581: EXTERN_C_END
583: EXTERN_C_BEGIN
586: PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
587: {
588: TS_Sundials *cvode = (TS_Sundials*)ts->data;
591: cvode->maxdt = maxdt;
592: return(0);
593: }
594: EXTERN_C_END
596: EXTERN_C_BEGIN
599: PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc)
600: {
601: SNES snes;
602: KSP ksp;
603: PetscErrorCode ierr;
606: TSGetSNES(ts,&snes);
607: SNESGetKSP(snes,&ksp);
608: KSPGetPC(ksp,pc);
609: return(0);
610: }
611: EXTERN_C_END
613: EXTERN_C_BEGIN
616: PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
617: {
619: if (nonlin) *nonlin = ts->snes_its;
620: if (lin) *lin = ts->ksp_its;
621: return(0);
622: }
623: EXTERN_C_END
625: EXTERN_C_BEGIN
628: PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
629: {
630: TS_Sundials *cvode = (TS_Sundials*)ts->data;
633: cvode->monitorstep = s;
634: return(0);
635: }
636: EXTERN_C_END
637: /* -------------------------------------------------------------------------------------------*/
641: /*@C
642: TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.
644: Not Collective
646: Input parameters:
647: . ts - the time-step context
649: Output Parameters:
650: + nonlin - number of nonlinear iterations
651: - lin - number of linear iterations
653: Level: advanced
655: Notes:
656: These return the number since the creation of the TS object
658: .keywords: non-linear iterations, linear iterations
660: .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
661: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
662: TSSundialsGetIterations(), TSSundialsSetType(),
663: TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()
665: @*/
666: PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
667: {
671: PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));
672: return(0);
673: }
677: /*@
678: TSSundialsSetType - Sets the method that Sundials will use for integration.
680: Logically Collective on TS
682: Input parameters:
683: + ts - the time-step context
684: - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF
686: Level: intermediate
688: .keywords: Adams, backward differentiation formula
690: .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(),
691: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
692: TSSundialsGetIterations(), TSSundialsSetType(),
693: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
694: TSSetExactFinalTime()
695: @*/
696: PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type)
697: {
701: PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));
702: return(0);
703: }
707: /*@
708: TSSundialsSetMaxl - Sets the dimension of the Krylov space used by
709: GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so
710: this is the maximum number of GMRES steps that will be used.
712: Logically Collective on TS
714: Input parameters:
715: + ts - the time-step context
716: - maxl - number of direction vectors (the dimension of Krylov subspace).
718: Level: advanced
720: .keywords: GMRES
722: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
723: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
724: TSSundialsGetIterations(), TSSundialsSetType(),
725: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
726: TSSetExactFinalTime()
728: @*/
729: PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl)
730: {
735: PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));
736: return(0);
737: }
741: /*@
742: TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
743: system by SUNDIALS.
745: Logically Collective on TS
747: Input parameters:
748: + ts - the time-step context
749: - tol - the factor by which the tolerance on the nonlinear solver is
750: multiplied to get the tolerance on the linear solver, .05 by default.
752: Level: advanced
754: .keywords: GMRES, linear convergence tolerance, SUNDIALS
756: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
757: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
758: TSSundialsGetIterations(), TSSundialsSetType(),
759: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
760: TSSetExactFinalTime()
762: @*/
763: PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol)
764: {
769: PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));
770: return(0);
771: }
775: /*@
776: TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
777: in GMRES method by SUNDIALS linear solver.
779: Logically Collective on TS
781: Input parameters:
782: + ts - the time-step context
783: - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS
785: Level: advanced
787: .keywords: Sundials, orthogonalization
789: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
790: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(),
791: TSSundialsGetIterations(), TSSundialsSetType(),
792: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
793: TSSetExactFinalTime()
795: @*/
796: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
797: {
801: PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));
802: return(0);
803: }
807: /*@
808: TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
809: Sundials for error control.
811: Logically Collective on TS
813: Input parameters:
814: + ts - the time-step context
815: . aabs - the absolute tolerance
816: - rel - the relative tolerance
818: See the Cvode/Sundials users manual for exact details on these parameters. Essentially
819: these regulate the size of the error for a SINGLE timestep.
821: Level: intermediate
823: .keywords: Sundials, tolerance
825: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
826: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
827: TSSundialsGetIterations(), TSSundialsSetType(),
828: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
829: TSSetExactFinalTime()
831: @*/
832: PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel)
833: {
837: PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));
838: return(0);
839: }
843: /*@
844: TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.
846: Input Parameter:
847: . ts - the time-step context
849: Output Parameter:
850: . pc - the preconditioner context
852: Level: advanced
854: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
855: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
856: TSSundialsGetIterations(), TSSundialsSetType(),
857: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
858: @*/
859: PetscErrorCode TSSundialsGetPC(TS ts,PC *pc)
860: {
864: PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC *),(ts,pc));
865: return(0);
866: }
870: /*@
871: TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.
873: Input Parameter:
874: + ts - the time-step context
875: - mindt - lowest time step if positive, negative to deactivate
877: Note:
878: Sundials will error if it is not possible to keep the estimated truncation error below
879: the tolerance set with TSSundialsSetTolerance() without going below this step size.
881: Level: beginner
883: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
884: @*/
885: PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
886: {
890: PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));
891: return(0);
892: }
896: /*@
897: TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.
899: Input Parameter:
900: + ts - the time-step context
901: - maxdt - lowest time step if positive, negative to deactivate
903: Level: beginner
905: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
906: @*/
907: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
908: {
912: PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
913: return(0);
914: }
918: /*@
919: TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false).
921: Input Parameter:
922: + ts - the time-step context
923: - ft - PETSC_TRUE if monitor, else PETSC_FALSE
925: Level: beginner
927: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
928: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
929: TSSundialsGetIterations(), TSSundialsSetType(),
930: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
931: @*/
932: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
933: {
937: PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));
938: return(0);
939: }
940: /* -------------------------------------------------------------------------------------------*/
941: /*MC
942: TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)
944: Options Database:
945: + -ts_sundials_type <bdf,adams>
946: . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
947: . -ts_sundials_atol <tol> - Absolute tolerance for convergence
948: . -ts_sundials_rtol <tol> - Relative tolerance for convergence
949: . -ts_sundials_linear_tolerance <tol>
950: . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
951: - -ts_sundials_monitor_steps - Monitor SUNDIALS internel steps
954: Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply
955: only PETSc PC options
957: Level: beginner
959: .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(),
960: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime()
962: M*/
963: EXTERN_C_BEGIN
966: PetscErrorCode TSCreate_Sundials(TS ts)
967: {
968: TS_Sundials *cvode;
970: PC pc;
973: ts->ops->reset = TSReset_Sundials;
974: ts->ops->destroy = TSDestroy_Sundials;
975: ts->ops->view = TSView_Sundials;
976: ts->ops->setup = TSSetUp_Sundials;
977: ts->ops->step = TSStep_Sundials;
978: ts->ops->interpolate = TSInterpolate_Sundials;
979: ts->ops->setfromoptions = TSSetFromOptions_Sundials;
981: PetscNewLog(ts,TS_Sundials,&cvode);
982: ts->data = (void*)cvode;
983: cvode->cvode_type = SUNDIALS_BDF;
984: cvode->gtype = SUNDIALS_CLASSICAL_GS;
985: cvode->maxl = 5;
986: cvode->linear_tol = .05;
988: cvode->monitorstep = PETSC_TRUE;
990: MPI_Comm_dup(((PetscObject)ts)->comm,&(cvode->comm_sundials));
992: cvode->mindt = -1.;
993: cvode->maxdt = -1.;
995: /* set tolerance for Sundials */
996: cvode->reltol = 1e-6;
997: cvode->abstol = 1e-6;
999: /* set PCNONE as default pctype */
1000: TSSundialsGetPC_Sundials(ts,&pc);
1001: PCSetType(pc,PCNONE);
1003: if (ts->exact_final_time == PETSC_DECIDE) ts->exact_final_time = PETSC_FALSE;
1005: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetType_C","TSSundialsSetType_Sundials",
1006: TSSundialsSetType_Sundials);
1007: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxl_C",
1008: "TSSundialsSetMaxl_Sundials",
1009: TSSundialsSetMaxl_Sundials);
1010: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetLinearTolerance_C",
1011: "TSSundialsSetLinearTolerance_Sundials",
1012: TSSundialsSetLinearTolerance_Sundials);
1013: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",
1014: "TSSundialsSetGramSchmidtType_Sundials",
1015: TSSundialsSetGramSchmidtType_Sundials);
1016: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetTolerance_C",
1017: "TSSundialsSetTolerance_Sundials",
1018: TSSundialsSetTolerance_Sundials);
1019: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMinTimeStep_C",
1020: "TSSundialsSetMinTimeStep_Sundials",
1021: TSSundialsSetMinTimeStep_Sundials);
1022: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",
1023: "TSSundialsSetMaxTimeStep_Sundials",
1024: TSSundialsSetMaxTimeStep_Sundials);
1025: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetPC_C",
1026: "TSSundialsGetPC_Sundials",
1027: TSSundialsGetPC_Sundials);
1028: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetIterations_C",
1029: "TSSundialsGetIterations_Sundials",
1030: TSSundialsGetIterations_Sundials);
1031: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",
1032: "TSSundialsMonitorInternalSteps_Sundials",
1033: TSSundialsMonitorInternalSteps_Sundials);
1035: return(0);
1036: }
1037: EXTERN_C_END