Actual source code: sundials.c
petsc-3.13.6 2020-09-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 <../src/ts/impls/implicit/sundials/sundials.h>
10: /*
11: TSPrecond_Sundials - function that we provide to SUNDIALS to
12: evaluate the preconditioner.
13: */
14: PetscErrorCode TSPrecond_Sundials(realtype tn,N_Vector y,N_Vector fy,booleantype jok,booleantype *jcurPtr,
15: realtype _gamma,void *P_data,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
16: {
17: TS ts = (TS) P_data;
18: TS_Sundials *cvode = (TS_Sundials*)ts->data;
19: PC pc;
21: Mat J,P;
22: Vec yy = cvode->w1,yydot = cvode->ydot;
23: PetscReal gm = (PetscReal)_gamma;
24: PetscScalar *y_data;
27: TSGetIJacobian(ts,&J,&P,NULL,NULL);
28: y_data = (PetscScalar*) N_VGetArrayPointer(y);
29: VecPlaceArray(yy,y_data);
30: VecZeroEntries(yydot); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */
31: /* compute the shifted Jacobian (1/gm)*I + Jrest */
32: TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,J,P,PETSC_FALSE);
33: VecResetArray(yy);
34: MatScale(P,gm); /* turn into I-gm*Jrest, J is not used by Sundials */
35: *jcurPtr = TRUE;
36: TSSundialsGetPC(ts,&pc);
37: PCSetOperators(pc,J,P);
38: return(0);
39: }
41: /*
42: TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner.
43: */
44: PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,
45: realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
46: {
47: TS ts = (TS) P_data;
48: TS_Sundials *cvode = (TS_Sundials*)ts->data;
49: PC pc;
50: Vec rr = cvode->w1,zz = cvode->w2;
52: PetscScalar *r_data,*z_data;
55: /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
56: r_data = (PetscScalar*) N_VGetArrayPointer(r);
57: z_data = (PetscScalar*) N_VGetArrayPointer(z);
58: VecPlaceArray(rr,r_data);
59: VecPlaceArray(zz,z_data);
61: /* Solve the Px=r and put the result in zz */
62: TSSundialsGetPC(ts,&pc);
63: PCApply(pc,rr,zz);
64: VecResetArray(rr);
65: VecResetArray(zz);
66: return(0);
67: }
69: /*
70: TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
71: */
72: int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
73: {
74: TS ts = (TS) ctx;
75: DM dm;
76: DMTS tsdm;
77: TSIFunction ifunction;
78: MPI_Comm comm;
79: TS_Sundials *cvode = (TS_Sundials*)ts->data;
80: Vec yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot;
81: PetscScalar *y_data,*ydot_data;
85: PetscObjectGetComm((PetscObject)ts,&comm);
86: /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
87: y_data = (PetscScalar*) N_VGetArrayPointer(y);
88: ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot);
89: VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
90: VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);
92: /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
93: TSGetDM(ts,&dm);
94: DMGetDMTS(dm,&tsdm);
95: DMTSGetIFunction(dm,&ifunction,NULL);
96: if (!ifunction) {
97: TSComputeRHSFunction(ts,t,yy,yyd);
98: } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
99: VecZeroEntries(yydot);
100: TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr);
101: VecScale(yyd,-1.);
102: }
103: VecResetArray(yy);CHKERRABORT(comm,ierr);
104: VecResetArray(yyd);CHKERRABORT(comm,ierr);
105: return(0);
106: }
108: /*
109: TSStep_Sundials - Calls Sundials to integrate the ODE.
110: */
111: PetscErrorCode TSStep_Sundials(TS ts)
112: {
113: TS_Sundials *cvode = (TS_Sundials*)ts->data;
115: PetscInt flag;
116: long int nits,lits,nsteps;
117: realtype t,tout;
118: PetscScalar *y_data;
119: void *mem;
122: mem = cvode->mem;
123: tout = ts->max_time;
124: VecGetArray(ts->vec_sol,&y_data);
125: N_VSetArrayPointer((realtype*)y_data,cvode->y);
126: VecRestoreArray(ts->vec_sol,NULL);
128: /* We would like to TSPreStage() and TSPostStage()
129: * before each stage solve but CVode does not appear to support this. */
130: if (cvode->monitorstep)
131: flag = CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
132: else
133: flag = CVode(mem,tout,cvode->y,&t,CV_NORMAL);
135: if (flag) { /* display error message */
136: switch (flag) {
137: case CV_ILL_INPUT:
138: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ILL_INPUT");
139: break;
140: case CV_TOO_CLOSE:
141: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_CLOSE");
142: break;
143: case CV_TOO_MUCH_WORK: {
144: PetscReal tcur;
145: CVodeGetNumSteps(mem,&nsteps);
146: CVodeGetCurrentTime(mem,&tcur);
147: 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 TSSetMaxSteps()",(double)tcur,nsteps,ts->max_steps);
148: } break;
149: case CV_TOO_MUCH_ACC:
150: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_TOO_MUCH_ACC");
151: break;
152: case CV_ERR_FAILURE:
153: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_ERR_FAILURE");
154: break;
155: case CV_CONV_FAILURE:
156: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_CONV_FAILURE");
157: break;
158: case CV_LINIT_FAIL:
159: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LINIT_FAIL");
160: break;
161: case CV_LSETUP_FAIL:
162: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSETUP_FAIL");
163: break;
164: case CV_LSOLVE_FAIL:
165: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_LSOLVE_FAIL");
166: break;
167: case CV_RHSFUNC_FAIL:
168: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RHSFUNC_FAIL");
169: break;
170: case CV_FIRST_RHSFUNC_ERR:
171: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_FIRST_RHSFUNC_ERR");
172: break;
173: case CV_REPTD_RHSFUNC_ERR:
174: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_REPTD_RHSFUNC_ERR");
175: break;
176: case CV_UNREC_RHSFUNC_ERR:
177: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_UNREC_RHSFUNC_ERR");
178: break;
179: case CV_RTFUNC_FAIL:
180: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, CV_RTFUNC_FAIL");
181: break;
182: default:
183: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVode() fails, flag %d",flag);
184: }
185: }
187: /* log inner nonlinear and linear iterations */
188: CVodeGetNumNonlinSolvIters(mem,&nits);
189: CVSpilsGetNumLinIters(mem,&lits);
190: ts->snes_its += nits; ts->ksp_its = lits;
192: /* copy the solution from cvode->y to cvode->update and sol */
193: VecPlaceArray(cvode->w1,y_data);
194: VecCopy(cvode->w1,cvode->update);
195: VecResetArray(cvode->w1);
196: VecCopy(cvode->update,ts->vec_sol);
198: ts->time_step = t - ts->ptime;
199: ts->ptime = t;
201: CVodeGetNumSteps(mem,&nsteps);
202: if (!cvode->monitorstep) ts->steps += nsteps - 1; /* TSStep() increments the step counter by one */
203: return(0);
204: }
206: static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
207: {
208: TS_Sundials *cvode = (TS_Sundials*)ts->data;
209: N_Vector y;
211: PetscScalar *x_data;
212: PetscInt glosize,locsize;
215: /* get the vector size */
216: VecGetSize(X,&glosize);
217: VecGetLocalSize(X,&locsize);
218: VecGetArray(X,&x_data);
220: /* Initialize N_Vec y with x_data */
221: y = N_VMake_Parallel(cvode->comm_sundials,locsize,glosize,(realtype*)x_data);
222: if (!y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Interpolated y is not allocated");
224: CVodeGetDky(cvode->mem,t,0,y);
225: VecRestoreArray(X,&x_data);
226: return(0);
227: }
229: PetscErrorCode TSReset_Sundials(TS ts)
230: {
231: TS_Sundials *cvode = (TS_Sundials*)ts->data;
235: VecDestroy(&cvode->update);
236: VecDestroy(&cvode->ydot);
237: VecDestroy(&cvode->w1);
238: VecDestroy(&cvode->w2);
239: if (cvode->mem) CVodeFree(&cvode->mem);
240: return(0);
241: }
243: PetscErrorCode TSDestroy_Sundials(TS ts)
244: {
245: TS_Sundials *cvode = (TS_Sundials*)ts->data;
249: TSReset_Sundials(ts);
250: MPI_Comm_free(&(cvode->comm_sundials));
251: PetscFree(ts->data);
252: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);
253: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);
254: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);
255: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);
256: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);
257: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);
258: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);
259: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);
260: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);
261: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",NULL);
262: return(0);
263: }
265: PetscErrorCode TSSetUp_Sundials(TS ts)
266: {
267: TS_Sundials *cvode = (TS_Sundials*)ts->data;
269: PetscInt glosize,locsize,i,flag;
270: PetscScalar *y_data,*parray;
271: void *mem;
272: PC pc;
273: PCType pctype;
274: PetscBool pcnone;
277: 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");
279: /* get the vector size */
280: VecGetSize(ts->vec_sol,&glosize);
281: VecGetLocalSize(ts->vec_sol,&locsize);
283: /* allocate the memory for N_Vec y */
284: cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
285: if (!cvode->y) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"cvode->y is not allocated");
287: /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
288: VecGetArray(ts->vec_sol,&parray);
289: y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);
290: for (i = 0; i < locsize; i++) y_data[i] = parray[i];
291: VecRestoreArray(ts->vec_sol,NULL);
293: VecDuplicate(ts->vec_sol,&cvode->update);
294: VecDuplicate(ts->vec_sol,&cvode->ydot);
295: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);
296: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);
298: /*
299: Create work vectors for the TSPSolve_Sundials() routine. Note these are
300: allocated with zero space arrays because the actual array space is provided
301: by Sundials and set using VecPlaceArray().
302: */
303: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w1);
304: VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,0,&cvode->w2);
305: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);
306: PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);
308: /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */
309: mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
310: if (!mem) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"CVodeCreate() fails");
311: cvode->mem = mem;
313: /* Set the pointer to user-defined data */
314: flag = CVodeSetUserData(mem, ts);
315: if (flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetUserData() fails");
317: /* Sundials may choose to use a smaller initial step, but will never use a larger step. */
318: flag = CVodeSetInitStep(mem,(realtype)ts->time_step);
319: if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetInitStep() failed");
320: if (cvode->mindt > 0) {
321: flag = CVodeSetMinStep(mem,(realtype)cvode->mindt);
322: if (flag) {
323: if (flag == CV_MEM_NULL) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed, cvode_mem pointer is NULL");
324: 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");
325: else SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMinStep() failed");
326: }
327: }
328: if (cvode->maxdt > 0) {
329: flag = CVodeSetMaxStep(mem,(realtype)cvode->maxdt);
330: if (flag) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_LIB,"CVodeSetMaxStep() failed");
331: }
333: /* Call CVodeInit to initialize the integrator memory and specify the
334: * user's right hand side function in u'=f(t,u), the inital time T0, and
335: * the initial dependent variable vector cvode->y */
336: flag = CVodeInit(mem,TSFunction_Sundials,ts->ptime,cvode->y);
337: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeInit() fails, flag %d",flag);
339: /* specifies scalar relative and absolute tolerances */
340: flag = CVodeSStolerances(mem,cvode->reltol,cvode->abstol);
341: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSStolerances() fails, flag %d",flag);
343: /* Specify max order of BDF / ADAMS method */
344: if (cvode->maxord != PETSC_DEFAULT) {
345: flag = CVodeSetMaxOrd(mem,cvode->maxord);
346: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetMaxOrd() fails, flag %d",flag);
347: }
349: /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
350: flag = CVodeSetMaxNumSteps(mem,ts->max_steps);
351: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVodeSetMaxNumSteps() fails, flag %d",flag);
353: /* call CVSpgmr to use GMRES as the linear solver. */
354: /* setup the ode integrator with the given preconditioner */
355: TSSundialsGetPC(ts,&pc);
356: PCGetType(pc,&pctype);
357: PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);
358: if (pcnone) {
359: flag = CVSpgmr(mem,PREC_NONE,0);
360: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
361: } else {
362: flag = CVSpgmr(mem,PREC_LEFT,cvode->maxl);
363: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmr() fails, flag %d",flag);
365: /* Set preconditioner and solve routines Precond and PSolve,
366: and the pointer to the user-defined block data */
367: flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials);
368: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpilsSetPreconditioner() fails, flag %d", flag);
369: }
371: flag = CVSpilsSetGSType(mem, MODIFIED_GS);
372: if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);
373: return(0);
374: }
376: /* type of CVODE linear multistep method */
377: const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",0};
378: /* type of G-S orthogonalization used by CVODE linear solver */
379: const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",0};
381: PetscErrorCode TSSetFromOptions_Sundials(PetscOptionItems *PetscOptionsObject,TS ts)
382: {
383: TS_Sundials *cvode = (TS_Sundials*)ts->data;
385: int indx;
386: PetscBool flag;
387: PC pc;
390: PetscOptionsHead(PetscOptionsObject,"SUNDIALS ODE solver options");
391: PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);
392: if (flag) {
393: TSSundialsSetType(ts,(TSSundialsLmmType)indx);
394: }
395: PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);
396: if (flag) {
397: TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
398: }
399: PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);
400: PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);
401: PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);
402: PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);
403: PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,NULL);
404: PetscOptionsInt("-ts_sundials_maxord","Max Order for BDF/Adams method","TSSundialsSetMaxOrd",cvode->maxord,&cvode->maxord,NULL);
405: PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,NULL);
406: PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internal steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);
407: PetscOptionsTail();
408: TSSundialsGetPC(ts,&pc);
409: PCSetFromOptions(pc);
410: return(0);
411: }
413: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
414: {
415: TS_Sundials *cvode = (TS_Sundials*)ts->data;
417: char *type;
418: char atype[] = "Adams";
419: char btype[] = "BDF: backward differentiation formula";
420: PetscBool iascii,isstring;
421: long int nsteps,its,nfevals,nlinsetups,nfails,itmp;
422: PetscInt qlast,qcur;
423: PetscReal hinused,hlast,hcur,tcur,tolsfac;
424: PC pc;
427: if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
428: else type = btype;
430: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
431: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
432: if (iascii) {
433: PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");
434: PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);
435: PetscViewerASCIIPrintf(viewer,"Sundials maxord %D\n",cvode->maxord);
436: PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",(double)cvode->abstol,(double)cvode->reltol);
437: PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",(double)cvode->linear_tol);
438: PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);
439: if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
440: PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
441: } else {
442: PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
443: }
444: if (cvode->mindt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",(double)cvode->mindt);}
445: if (cvode->maxdt > 0) {PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",(double)cvode->maxdt);}
447: /* Outputs from CVODE, CVSPILS */
448: CVodeGetTolScaleFactor(cvode->mem,&tolsfac);
449: PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);
450: CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,
451: &nlinsetups,&nfails,&qlast,&qcur,
452: &hinused,&hlast,&hcur,&tcur);
453: PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);
454: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);
455: PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);
456: PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);
458: CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);
459: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);
460: PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);
462: CVSpilsGetNumLinIters(cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
463: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);
464: CVSpilsGetNumConvFails(cvode->mem,&itmp);
465: PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);
467: TSSundialsGetPC(ts,&pc);
468: PCView(pc,viewer);
469: CVSpilsGetNumPrecEvals(cvode->mem,&itmp);
470: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);
471: CVSpilsGetNumPrecSolves(cvode->mem,&itmp);
472: PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);
474: CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);
475: PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);
476: CVSpilsGetNumRhsEvals(cvode->mem,&itmp);
477: PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);
478: } else if (isstring) {
479: PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
480: }
481: return(0);
482: }
485: /* --------------------------------------------------------------------------*/
486: PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
487: {
488: TS_Sundials *cvode = (TS_Sundials*)ts->data;
491: cvode->cvode_type = type;
492: return(0);
493: }
495: PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
496: {
497: TS_Sundials *cvode = (TS_Sundials*)ts->data;
500: cvode->maxl = maxl;
501: return(0);
502: }
504: PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
505: {
506: TS_Sundials *cvode = (TS_Sundials*)ts->data;
509: cvode->linear_tol = tol;
510: return(0);
511: }
513: PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
514: {
515: TS_Sundials *cvode = (TS_Sundials*)ts->data;
518: cvode->gtype = type;
519: return(0);
520: }
522: PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
523: {
524: TS_Sundials *cvode = (TS_Sundials*)ts->data;
527: if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
528: if (rel != PETSC_DECIDE) cvode->reltol = rel;
529: return(0);
530: }
532: PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
533: {
534: TS_Sundials *cvode = (TS_Sundials*)ts->data;
537: cvode->mindt = mindt;
538: return(0);
539: }
541: PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
542: {
543: TS_Sundials *cvode = (TS_Sundials*)ts->data;
546: cvode->maxdt = maxdt;
547: return(0);
548: }
549: PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc)
550: {
551: SNES snes;
552: KSP ksp;
556: TSGetSNES(ts,&snes);
557: SNESGetKSP(snes,&ksp);
558: KSPGetPC(ksp,pc);
559: return(0);
560: }
562: PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
563: {
565: if (nonlin) *nonlin = ts->snes_its;
566: if (lin) *lin = ts->ksp_its;
567: return(0);
568: }
570: PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
571: {
572: TS_Sundials *cvode = (TS_Sundials*)ts->data;
575: cvode->monitorstep = s;
576: return(0);
577: }
578: /* -------------------------------------------------------------------------------------------*/
580: /*@C
581: TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.
583: Not Collective
585: Input parameters:
586: . ts - the time-step context
588: Output Parameters:
589: + nonlin - number of nonlinear iterations
590: - lin - number of linear iterations
592: Level: advanced
594: Notes:
595: These return the number since the creation of the TS object
597: .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
598: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
599: TSSundialsGetIterations(), TSSundialsSetType(),
600: TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()
602: @*/
603: PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
604: {
608: PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));
609: return(0);
610: }
612: /*@
613: TSSundialsSetType - Sets the method that Sundials will use for integration.
615: Logically Collective on TS
617: Input parameters:
618: + ts - the time-step context
619: - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF
621: Level: intermediate
623: .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(),
624: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
625: TSSundialsGetIterations(), TSSundialsSetType(),
626: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
627: TSSetExactFinalTime()
628: @*/
629: PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type)
630: {
634: PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));
635: return(0);
636: }
638: /*@
639: TSSundialsSetMaxord - Sets the maximum order for BDF/Adams method used by SUNDIALS.
641: Logically Collective on TS
643: Input parameters:
644: + ts - the time-step context
645: - maxord - maximum order of BDF / Adams method
647: Level: advanced
649: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
650: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
651: TSSundialsGetIterations(), TSSundialsSetType(),
652: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
653: TSSetExactFinalTime()
655: @*/
656: PetscErrorCode TSSundialsSetMaxord(TS ts,PetscInt maxord)
657: {
662: PetscTryMethod(ts,"TSSundialsSetMaxOrd_C",(TS,PetscInt),(ts,maxord));
663: return(0);
664: }
666: /*@
667: TSSundialsSetMaxl - Sets the dimension of the Krylov space used by
668: GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so
669: this is the maximum number of GMRES steps that will be used.
671: Logically Collective on TS
673: Input parameters:
674: + ts - the time-step context
675: - maxl - number of direction vectors (the dimension of Krylov subspace).
677: Level: advanced
679: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
680: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
681: TSSundialsGetIterations(), TSSundialsSetType(),
682: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
683: TSSetExactFinalTime()
685: @*/
686: PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl)
687: {
692: PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));
693: return(0);
694: }
696: /*@
697: TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
698: system by SUNDIALS.
700: Logically Collective on TS
702: Input parameters:
703: + ts - the time-step context
704: - tol - the factor by which the tolerance on the nonlinear solver is
705: multiplied to get the tolerance on the linear solver, .05 by default.
707: Level: advanced
709: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
710: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
711: TSSundialsGetIterations(), TSSundialsSetType(),
712: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
713: TSSetExactFinalTime()
715: @*/
716: PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol)
717: {
722: PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));
723: return(0);
724: }
726: /*@
727: TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
728: in GMRES method by SUNDIALS linear solver.
730: Logically Collective on TS
732: Input parameters:
733: + ts - the time-step context
734: - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS
736: Level: advanced
738: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
739: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(),
740: TSSundialsGetIterations(), TSSundialsSetType(),
741: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
742: TSSetExactFinalTime()
744: @*/
745: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
746: {
750: PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));
751: return(0);
752: }
754: /*@
755: TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
756: Sundials for error control.
758: Logically Collective on TS
760: Input parameters:
761: + ts - the time-step context
762: . aabs - the absolute tolerance
763: - rel - the relative tolerance
765: See the Cvode/Sundials users manual for exact details on these parameters. Essentially
766: these regulate the size of the error for a SINGLE timestep.
768: Level: intermediate
770: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
771: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
772: TSSundialsGetIterations(), TSSundialsSetType(),
773: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
774: TSSetExactFinalTime()
776: @*/
777: PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel)
778: {
782: PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));
783: return(0);
784: }
786: /*@
787: TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.
789: Input Parameter:
790: . ts - the time-step context
792: Output Parameter:
793: . pc - the preconditioner context
795: Level: advanced
797: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
798: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
799: TSSundialsGetIterations(), TSSundialsSetType(),
800: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
801: @*/
802: PetscErrorCode TSSundialsGetPC(TS ts,PC *pc)
803: {
807: PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));
808: return(0);
809: }
811: /*@
812: TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.
814: Input Parameter:
815: + ts - the time-step context
816: - mindt - lowest time step if positive, negative to deactivate
818: Note:
819: Sundials will error if it is not possible to keep the estimated truncation error below
820: the tolerance set with TSSundialsSetTolerance() without going below this step size.
822: Level: beginner
824: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
825: @*/
826: PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
827: {
831: PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));
832: return(0);
833: }
835: /*@
836: TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.
838: Input Parameter:
839: + ts - the time-step context
840: - maxdt - lowest time step if positive, negative to deactivate
842: Level: beginner
844: .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
845: @*/
846: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
847: {
851: PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));
852: return(0);
853: }
855: /*@
856: TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false).
858: Input Parameter:
859: + ts - the time-step context
860: - ft - PETSC_TRUE if monitor, else PETSC_FALSE
862: Level: beginner
864: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
865: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
866: TSSundialsGetIterations(), TSSundialsSetType(),
867: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
868: @*/
869: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
870: {
874: PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));
875: return(0);
876: }
877: /* -------------------------------------------------------------------------------------------*/
878: /*MC
879: TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)
881: Options Database:
882: + -ts_sundials_type <bdf,adams> -
883: . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
884: . -ts_sundials_atol <tol> - Absolute tolerance for convergence
885: . -ts_sundials_rtol <tol> - Relative tolerance for convergence
886: . -ts_sundials_linear_tolerance <tol> -
887: . -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
888: - -ts_sundials_monitor_steps - Monitor SUNDIALS internal steps
891: Notes:
892: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply,
893: only PETSc PC options.
895: Level: beginner
897: .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(),
898: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime()
900: M*/
901: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
902: {
903: TS_Sundials *cvode;
905: PC pc;
908: ts->ops->reset = TSReset_Sundials;
909: ts->ops->destroy = TSDestroy_Sundials;
910: ts->ops->view = TSView_Sundials;
911: ts->ops->setup = TSSetUp_Sundials;
912: ts->ops->step = TSStep_Sundials;
913: ts->ops->interpolate = TSInterpolate_Sundials;
914: ts->ops->setfromoptions = TSSetFromOptions_Sundials;
915: ts->default_adapt_type = TSADAPTNONE;
917: PetscNewLog(ts,&cvode);
919: ts->usessnes = PETSC_TRUE;
921: ts->data = (void*)cvode;
922: cvode->cvode_type = SUNDIALS_BDF;
923: cvode->gtype = SUNDIALS_CLASSICAL_GS;
924: cvode->maxl = 5;
925: cvode->maxord = PETSC_DEFAULT;
926: cvode->linear_tol = .05;
927: cvode->monitorstep = PETSC_TRUE;
929: MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));
931: cvode->mindt = -1.;
932: cvode->maxdt = -1.;
934: /* set tolerance for Sundials */
935: cvode->reltol = 1e-6;
936: cvode->abstol = 1e-6;
938: /* set PCNONE as default pctype */
939: TSSundialsGetPC_Sundials(ts,&pc);
940: PCSetType(pc,PCNONE);
942: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);
943: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);
944: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);
945: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);
946: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);
947: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);
948: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);
949: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);
950: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);
951: PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);
952: return(0);
953: }