Actual source code: sundials.c
petsc-3.7.5 2017-01-01
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> /*I "petscts.h" I*/
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: PetscScalar *y_data;
29: TSGetIJacobian(ts,&J,&P,NULL,NULL);
30: y_data = (PetscScalar*) N_VGetArrayPointer(y);
31: VecPlaceArray(yy,y_data);
32: VecZeroEntries(yydot); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */
33: /* compute the shifted Jacobian (1/gm)*I + Jrest */
34: TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,J,P,PETSC_FALSE);
35: VecResetArray(yy);
36: MatScale(P,gm); /* turn into I-gm*Jrest, J is not used by Sundials */
37: *jcurPtr = TRUE;
38: TSSundialsGetPC(ts,&pc);
39: PCSetOperators(pc,J,P);
40: return(0);
41: }
43: /*
44: TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner.
45: */
48: PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,
49: realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
50: {
51: TS ts = (TS) P_data;
52: TS_Sundials *cvode = (TS_Sundials*)ts->data;
53: PC pc;
54: Vec rr = cvode->w1,zz = cvode->w2;
56: PetscScalar *r_data,*z_data;
59: /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
60: r_data = (PetscScalar*) N_VGetArrayPointer(r);
61: z_data = (PetscScalar*) N_VGetArrayPointer(z);
62: VecPlaceArray(rr,r_data);
63: VecPlaceArray(zz,z_data);
65: /* Solve the Px=r and put the result in zz */
66: TSSundialsGetPC(ts,&pc);
67: PCApply(pc,rr,zz);
68: VecResetArray(rr);
69: VecResetArray(zz);
70: return(0);
71: }
73: /*
74: TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
75: */
78: int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
79: {
80: TS ts = (TS) ctx;
81: DM dm;
82: DMTS tsdm;
83: TSIFunction ifunction;
84: MPI_Comm 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;
91: PetscObjectGetComm((PetscObject)ts,&comm);
92: /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
93: y_data = (PetscScalar*) N_VGetArrayPointer(y);
94: ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);
95: VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
96: VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);
98: /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
99: TSGetDM(ts,&dm);
100: DMGetDMTS(dm,&tsdm);
101: DMTSGetIFunction(dm,&ifunction,NULL);
102: if (!ifunction) {
103: TSComputeRHSFunction(ts,t,yy,yyd);
104: } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
105: VecZeroEntries(yydot);
106: TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr);
107: VecScale(yyd,-1.);
108: }
109: VecResetArray(yy);CHKERRABORT(comm,ierr);
110: VecResetArray(yyd);CHKERRABORT(comm,ierr);
111: return(0);
112: }
114: /*
115: TSStep_Sundials - Calls Sundials to integrate the ODE.
116: */
119: PetscErrorCode TSStep_Sundials(TS ts)
120: {
121: TS_Sundials *cvode = (TS_Sundials*)ts->data;
123: PetscInt flag;
124: long int nits,lits,nsteps;
125: realtype t,tout;
126: PetscScalar *y_data;
127: void *mem;
130: mem = cvode->mem;
131: tout = ts->max_time;
132: VecGetArray(ts->vec_sol,&y_data);
133: N_VSetArrayPointer((realtype*)y_data,cvode->y);
134: VecRestoreArray(ts->vec_sol,NULL);
136: /* We would like to TSPreStage() and TSPostStage()
137: * before each stage solve but CVode does not appear to support this. */
138: if (cvode->monitorstep)
139: flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
140: else
141: flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);
143: if (flag) { /* display error message */
144: switch (flag) {
145: case CV_ILL_INPUT:
146: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");
147: break;
148: case CV_TOO_CLOSE:
149: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");
150: break;
151: case CV_TOO_MUCH_WORK: {
152: PetscReal tcur;
153: CVodeGetNumSteps(mem,&nsteps);
154: CVodeGetCurrentTime(mem,&tcur);
155: SETERRQ3(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 TSSetDuration()",(double)tcur,nsteps,ts->max_steps);
156: } break;
157: case CV_TOO_MUCH_ACC:
158: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");
159: break;
160: case CV_ERR_FAILURE:
161: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");
162: break;
163: case CV_CONV_FAILURE:
164: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");
165: break;
166: case CV_LINIT_FAIL:
167: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");
168: break;
169: case CV_LSETUP_FAIL:
170: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");
171: break;
172: case CV_LSOLVE_FAIL:
173: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");
174: break;
175: case CV_RHSFUNC_FAIL:
176: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");
177: break;
178: case CV_FIRST_RHSFUNC_ERR:
179: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");
180: break;
181: case CV_REPTD_RHSFUNC_ERR:
182: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");
183: break;
184: case CV_UNREC_RHSFUNC_ERR:
185: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");
186: break;
187: case CV_RTFUNC_FAIL:
188: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");
189: break;
190: default:
191: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);
192: }
193: }
195: /* log inner nonlinear and linear iterations */
196: CVodeGetNumNonlinSolvIters(mem,&nits);
197: CVSpilsGetNumLinIters(mem,&lits);
198: ts->snes_its += nits; ts->ksp_its = lits;
200: /* copy the solution from cvode->y to cvode->update and sol */
201: VecPlaceArray(cvode->w1,y_data);
202: VecCopy(cvode->w1,cvode->update);
203: VecResetArray(cvode->w1);
204: VecCopy(cvode->update,ts->vec_sol);
206: ts->time_step = t - ts->ptime;
207: ts->ptime = t;
209: CVodeGetNumSteps(mem,&nsteps);
210: if (!cvode->monitorstep) ts->steps = nsteps - 1; /* TSStep() increments the step counter by one */
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: if (ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for exact final time option 'MATCHSTEP' when using Sundials");
296: /* get the vector size */
297: VecGetSize(ts->vec_sol,&glosize);
298: VecGetLocalSize(ts->vec_sol,&locsize);
300: /* allocate the memory for N_Vec y */
301: cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
302: if (!cvode->y) SETERRQ(PETSC_COMM_SELF,1,"cvode->y is not allocated");
304: /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
305: VecGetArray(ts->vec_sol,&parray);
306: y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);
307: for (i = 0; i < locsize; i++) y_data[i] = parray[i];
308: VecRestoreArray(ts->vec_sol,NULL);
310: VecDuplicate(ts->vec_sol,&cvode->update);
311: VecDuplicate(ts->vec_sol,&cvode->ydot);
312: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);
313: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);
315: /*
316: Create work vectors for the TSPSolve_Sundials() routine. Note these are
317: allocated with zero space arrays because the actual array space is provided
318: by Sundials and set using VecPlaceArray().
319: */
320: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w1);
321: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w2);
322: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);
323: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);
325: /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */
326: mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
327: if (!mem) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"CVodeCreate() fails");
328: cvode->mem = mem;
330: /* Set the pointer to user-defined data */
331: flag = CVodeSetUserData(mem, ts);
332: if (flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetUserData() fails");
334: /* Sundials may choose to use a smaller initial step, but will never use a larger step. */
335: flag = CVodeSetInitStep(mem,(realtype)ts->time_step);
336: if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetInitStep() failed");
337: if (cvode->mindt > 0) {
338: flag = CVodeSetMinStep(mem,(realtype)cvode->mindt);
339: if (flag) {
340: if (flag == CV_MEM_NULL) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL");
341: 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");
342: else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed");
343: }
344: }
345: if (cvode->maxdt > 0) {
346: flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);
347: if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMaxStep() failed");
348: }
350: /* Call CVodeInit to initialize the integrator memory and specify the
351: * user's right hand side function in u'=f(t,u), the inital time T0, and
352: * the initial dependent variable vector cvode->y */
353: flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y);
354: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag);
356: /* specifies scalar relative and absolute tolerances */
357: flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol);
358: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag);
360: /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
361: flag = CVodeSetMaxNumSteps(mem,ts->max_steps);
362: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetMaxNumSteps() fails, flag %d",flag);
364: /* call CVSpgmr to use GMRES as the linear solver. */
365: /* setup the ode integrator with the given preconditioner */
366: TSSundialsGetPC(ts,&pc);
367: PCGetType(pc,&pctype);
368: PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);
369: if (pcnone) {
370: flag = CVSpgmr(mem,PREC_NONE,0);
371: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
372: } else {
373: flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl);
374: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
376: /* Set preconditioner and solve routines Precond and PSolve,
377: and the pointer to the user-defined block data */
378: flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);
379: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag);
380: }
382: flag = CVSpilsSetGSType(mem, MODIFIED_GS);
383: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);
384: return(0);
385: }
387: /* type of CVODE linear multistep method */
388: const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",0};
389: /* type of G-S orthogonalization used by CVODE linear solver */
390: const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",0};
394: PetscErrorCode TSSetFromOptions_Sundials(PetscOptionItems *PetscOptionsObject,TS ts)
395: {
396: TS_Sundials *cvode = (TS_Sundials*)ts->data;
398: int indx;
399: PetscBool flag;
400: PC pc;
403: PetscOptionsHead(PetscOptionsObject,"SUNDIALS ODE solver options");
404: PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);
405: if (flag) {
406: TSSundialsSetType(ts,(TSSundialsLmmType)indx);
407: }
408: PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);
409: if (flag) {
410: TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
411: }
412: PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);
413: PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);
414: PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);
415: PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);
416: PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,NULL);
417: PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,NULL);
418: PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internal steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);
419: PetscOptionsTail();
420: TSSundialsGetPC(ts,&pc);
421: PCSetFromOptions(pc);
422: return(0);
423: }
427: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
428: {
429: TS_Sundials *cvode = (TS_Sundials*)ts->data;
431: char *type;
432: char atype[] = "Adams";
433: char btype[] = "BDF: backward differentiation formula";
434: PetscBool iascii,isstring;
435: long int nsteps,its,nfevals,nlinsetups,nfails,itmp;
436: PetscInt qlast,qcur;
437: PetscReal hinused,hlast,hcur,tcur,tolsfac;
438: PC pc;
441: if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
442: else type = btype;
444: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
445: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
446: if (iascii) {
447: PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");
448: PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);
449: PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);
450: PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",cvode->linear_tol);
451: PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);
452: if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
453: PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
454: } else {
455: PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
456: }
457: if (cvode->mindt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",cvode->mindt);}
458: if (cvode->maxdt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",cvode->maxdt);}
460: /* Outputs from CVODE, CVSPILS */
461: CVodeGetTolScaleFactor(cvode->mem,&tolsfac);
462: PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);
463: CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,
464: &nlinsetups,&nfails,&qlast,&qcur,
465: &hinused,&hlast,&hcur,&tcur);
466: PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);
467: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);
468: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);
469: PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);
471: CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);
472: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);
473: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);
475: CVSpilsGetNumLinIters(cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
476: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);
477: CVSpilsGetNumConvFails(cvode->mem,&itmp);
478: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);
480: TSSundialsGetPC(ts,&pc);
481: PCView(pc,viewer);
482: CVSpilsGetNumPrecEvals(cvode->mem,&itmp);
483: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);
484: CVSpilsGetNumPrecSolves(cvode->mem,&itmp);
485: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);
487: CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);
488: PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);
489: CVSpilsGetNumRhsEvals(cvode->mem,&itmp);
490: PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);
491: } else if (isstring) {
492: PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
493: }
494: return(0);
495: }
498: /* --------------------------------------------------------------------------*/
501: PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
502: {
503: TS_Sundials *cvode = (TS_Sundials*)ts->data;
506: cvode->cvode_type = type;
507: return(0);
508: }
512: PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
513: {
514: TS_Sundials *cvode = (TS_Sundials*)ts->data;
517: cvode->maxl = maxl;
518: return(0);
519: }
523: PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
524: {
525: TS_Sundials *cvode = (TS_Sundials*)ts->data;
528: cvode->linear_tol = tol;
529: return(0);
530: }
534: PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
535: {
536: TS_Sundials *cvode = (TS_Sundials*)ts->data;
539: cvode->gtype = type;
540: return(0);
541: }
545: PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
546: {
547: TS_Sundials *cvode = (TS_Sundials*)ts->data;
550: if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
551: if (rel != PETSC_DECIDE) cvode->reltol = rel;
552: return(0);
553: }
557: PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
558: {
559: TS_Sundials *cvode = (TS_Sundials*)ts->data;
562: cvode->mindt = mindt;
563: return(0);
564: }
568: PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
569: {
570: TS_Sundials *cvode = (TS_Sundials*)ts->data;
573: cvode->maxdt = maxdt;
574: return(0);
575: }
578: PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc)
579: {
580: SNES snes;
581: KSP ksp;
585: TSGetSNES(ts,&snes);
586: SNESGetKSP(snes,&ksp);
587: KSPGetPC(ksp,pc);
588: return(0);
589: }
593: PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
594: {
596: if (nonlin) *nonlin = ts->snes_its;
597: if (lin) *lin = ts->ksp_its;
598: return(0);
599: }
603: PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
604: {
605: TS_Sundials *cvode = (TS_Sundials*)ts->data;
608: cvode->monitorstep = s;
609: return(0);
610: }
611: /* -------------------------------------------------------------------------------------------*/
615: /*@C
616: TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.
618: Not Collective
620: Input parameters:
621: . ts - the time-step context
623: Output Parameters:
624: + nonlin - number of nonlinear iterations
625: - lin - number of linear iterations
627: Level: advanced
629: Notes:
630: These return the number since the creation of the TS object
632: .keywords: non-linear iterations, linear iterations
634: .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
635: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
636: TSSundialsGetIterations(), TSSundialsSetType(),
637: TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()
639: @*/
640: PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
641: {
645: PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));
646: return(0);
647: }
651: /*@
652: TSSundialsSetType - Sets the method that Sundials will use for integration.
654: Logically Collective on TS
656: Input parameters:
657: + ts - the time-step context
658: - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF
660: Level: intermediate
662: .keywords: Adams, backward differentiation formula
664: .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(),
665: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
666: TSSundialsGetIterations(), TSSundialsSetType(),
667: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
668: TSSetExactFinalTime()
669: @*/
670: PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type)
671: {
675: PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));
676: return(0);
677: }
681: /*@
682: TSSundialsSetMaxl - Sets the dimension of the Krylov space used by
683: GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so
684: this is the maximum number of GMRES steps that will be used.
686: Logically Collective on TS
688: Input parameters:
689: + ts - the time-step context
690: - maxl - number of direction vectors (the dimension of Krylov subspace).
692: Level: advanced
694: .keywords: GMRES
696: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
697: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
698: TSSundialsGetIterations(), TSSundialsSetType(),
699: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
700: TSSetExactFinalTime()
702: @*/
703: PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl)
704: {
709: PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));
710: return(0);
711: }
715: /*@
716: TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
717: system by SUNDIALS.
719: Logically Collective on TS
721: Input parameters:
722: + ts - the time-step context
723: - tol - the factor by which the tolerance on the nonlinear solver is
724: multiplied to get the tolerance on the linear solver, .05 by default.
726: Level: advanced
728: .keywords: GMRES, linear convergence tolerance, SUNDIALS
730: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
731: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
732: TSSundialsGetIterations(), TSSundialsSetType(),
733: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
734: TSSetExactFinalTime()
736: @*/
737: PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol)
738: {
743: PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));
744: return(0);
745: }
749: /*@
750: TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
751: in GMRES method by SUNDIALS linear solver.
753: Logically Collective on TS
755: Input parameters:
756: + ts - the time-step context
757: - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS
759: Level: advanced
761: .keywords: Sundials, orthogonalization
763: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
764: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(),
765: TSSundialsGetIterations(), TSSundialsSetType(),
766: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
767: TSSetExactFinalTime()
769: @*/
770: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
771: {
775: PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));
776: return(0);
777: }
781: /*@
782: TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
783: Sundials for error control.
785: Logically Collective on TS
787: Input parameters:
788: + ts - the time-step context
789: . aabs - the absolute tolerance
790: - rel - the relative tolerance
792: See the Cvode/Sundials users manual for exact details on these parameters. Essentially
793: these regulate the size of the error for a SINGLE timestep.
795: Level: intermediate
797: .keywords: Sundials, tolerance
799: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
800: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
801: TSSundialsGetIterations(), TSSundialsSetType(),
802: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
803: TSSetExactFinalTime()
805: @*/
806: PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel)
807: {
811: PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));
812: return(0);
813: }
817: /*@
818: TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.
820: Input Parameter:
821: . ts - the time-step context
823: Output Parameter:
824: . pc - the preconditioner context
826: Level: advanced
828: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
829: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
830: TSSundialsGetIterations(), TSSundialsSetType(),
831: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
832: @*/
833: PetscErrorCode TSSundialsGetPC(TS ts,PC *pc)
834: {
838: PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));
839: return(0);
840: }
844: /*@
845: TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.
847: Input Parameter:
848: + ts - the time-step context
849: - mindt - lowest time step if positive, negative to deactivate
851: Note:
852: Sundials will error if it is not possible to keep the estimated truncation error below
853: the tolerance set with TSSundialsSetTolerance() without going below this step size.
855: Level: beginner
857: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
858: @*/
859: PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
860: {
864: PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));
865: return(0);
866: }
870: /*@
871: TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.
873: Input Parameter:
874: + ts - the time-step context
875: - maxdt - lowest time step if positive, negative to deactivate
877: Level: beginner
879: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
880: @*/
881: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
882: {
886: PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
887: return(0);
888: }
892: /*@
893: TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false).
895: Input Parameter:
896: + ts - the time-step context
897: - ft - PETSC_TRUE if monitor, else PETSC_FALSE
899: Level: beginner
901: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
902: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
903: TSSundialsGetIterations(), TSSundialsSetType(),
904: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
905: @*/
906: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
907: {
911: PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));
912: return(0);
913: }
914: /* -------------------------------------------------------------------------------------------*/
915: /*MC
916: TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)
918: Options Database:
919: + -ts_sundials_type <bdf,adams>
920: . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
921: . -ts_sundials_atol <tol> - Absolute tolerance for convergence
922: . -ts_sundials_rtol <tol> - Relative tolerance for convergence
923: . -ts_sundials_linear_tolerance <tol>
924: . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
925: - -ts_sundials_monitor_steps - Monitor SUNDIALS internel steps
928: Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply
929: only PETSc PC options
931: Level: beginner
933: .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(),
934: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime()
936: M*/
939: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
940: {
941: TS_Sundials *cvode;
943: PC pc;
946: ts->ops->reset = TSReset_Sundials;
947: ts->ops->destroy = TSDestroy_Sundials;
948: ts->ops->view = TSView_Sundials;
949: ts->ops->setup = TSSetUp_Sundials;
950: ts->ops->step = TSStep_Sundials;
951: ts->ops->interpolate = TSInterpolate_Sundials;
952: ts->ops->setfromoptions = TSSetFromOptions_Sundials;
954: PetscNewLog(ts,&cvode);
956: ts->data = (void*)cvode;
957: cvode->cvode_type = SUNDIALS_BDF;
958: cvode->gtype = SUNDIALS_CLASSICAL_GS;
959: cvode->maxl = 5;
960: cvode->linear_tol = .05;
961: cvode->monitorstep = PETSC_TRUE;
963: MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));
965: cvode->mindt = -1.;
966: cvode->maxdt = -1.;
968: /* set tolerance for Sundials */
969: cvode->reltol = 1e-6;
970: cvode->abstol = 1e-6;
972: /* set PCNONE as default pctype */
973: TSSundialsGetPC_Sundials(ts,&pc);
974: PCSetType(pc,PCNONE);
976: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);
977: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);
978: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);
979: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);
980: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);
981: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);
982: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);
983: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);
984: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);
985: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);
986: return(0);
987: }