Actual source code: sundials.c
1: #define PETSCTS_DLL
3: /*
4: Provides a PETSc interface to SUNDIALS. Alan Hindmarsh's parallel ODE
5: solver.
6: The interface to PVODE (old version of CVODE) was originally contributed
7: by Liyang Xu. It has been redone by Hong Zhang and Dinesh Kaushik.
8: */
10: #include "src/ts/impls/implicit/sundials/sundials.h" /*I "petscts.h" I*/
12: /*
13: TSPrecond_Sundials - function that we provide to SUNDIALS to
14: evaluate the preconditioner.
15: */
18: PetscErrorCode TSPrecond_Sundials(realtype tn,N_Vector y,N_Vector fy,
19: booleantype jok,booleantype *jcurPtr,
20: realtype _gamma,void *P_data,
21: N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
22: {
23: TS ts = (TS) P_data;
24: TS_Sundials *cvode = (TS_Sundials*)ts->data;
25: PC pc = cvode->pc;
27: Mat Jac = ts->B;
28: Vec yy = cvode->w1;
29: PetscScalar one = 1.0,gm;
30: MatStructure str = DIFFERENT_NONZERO_PATTERN;
31: PetscScalar *y_data;
32:
34: /* This allows us to construct preconditioners in-place if we like */
35: MatSetUnfactored(Jac);
36:
37: /* jok - TRUE means reuse current Jacobian else recompute Jacobian */
38: if (jok) {
39: MatCopy(cvode->pmat,Jac,str);
40: *jcurPtr = FALSE;
41: } else {
42: /* make PETSc vector yy point to SUNDIALS vector y */
43: y_data = (PetscScalar *) N_VGetArrayPointer(y);
44: VecPlaceArray(yy,y_data);
46: /* compute the Jacobian */
47: TSComputeRHSJacobian(ts,ts->ptime,yy,&Jac,&Jac,&str);
48: VecResetArray(yy);
50: /* copy the Jacobian matrix */
51: if (!cvode->pmat) {
52: MatDuplicate(Jac,MAT_COPY_VALUES,&cvode->pmat);
53: PetscLogObjectParent(ts,cvode->pmat);
54: }
55: else {
56: MatCopy(Jac,cvode->pmat,str);
57: }
58: *jcurPtr = TRUE;
59: }
60:
61: /* construct I-gamma*Jac */
62: gm = -_gamma;
63: MatScale(Jac,gm);
64: MatShift(Jac,one);
65:
66: PCSetOperators(pc,Jac,Jac,str);
67: return(0);
68: }
70: /*
71: TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner.
72: */
75: PetscErrorCode TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,
76: realtype _gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
77: {
78: TS ts = (TS) P_data;
79: TS_Sundials *cvode = (TS_Sundials*)ts->data;
80: PC pc = cvode->pc;
81: Vec rr = cvode->w1,zz = cvode->w2;
82: PetscErrorCode ierr;
83: PetscScalar *r_data,*z_data;
86: /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
87: r_data = (PetscScalar *) N_VGetArrayPointer(r);
88: z_data = (PetscScalar *) N_VGetArrayPointer(z);
89: VecPlaceArray(rr,r_data);
90: VecPlaceArray(zz,z_data);
92: /* Solve the Px=r and put the result in zz */
93: PCApply(pc,rr,zz);
94: VecResetArray(rr);
95: VecResetArray(zz);
96: cvode->linear_solves++;
97: return(0);
98: }
100: /*
101: TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
102: */
105: int TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void *ctx)
106: {
107: TS ts = (TS) ctx;
108: MPI_Comm comm = ts->comm;
109: TS_Sundials *cvode = (TS_Sundials*)ts->data;
110: Vec yy = cvode->w1,yyd = cvode->w2;
111: PetscScalar *y_data,*ydot_data;
112: PetscErrorCode ierr;
115: /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
116: y_data = (PetscScalar *) N_VGetArrayPointer(y);
117: ydot_data = (PetscScalar *) N_VGetArrayPointer(ydot);
118: VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
119: VecPlaceArray(yyd,ydot_data); CHKERRABORT(comm,ierr);
121: /* now compute the right hand side function */
122: TSComputeRHSFunction(ts,t,yy,yyd); CHKERRABORT(comm,ierr);
123: VecResetArray(yy); CHKERRABORT(comm,ierr);
124: VecResetArray(yyd); CHKERRABORT(comm,ierr);
125: return(0);
126: }
128: /*
129: TSStep_Sundials_Nonlinear - Calls Sundials to integrate the ODE.
130: */
133: /*
134: TSStep_Sundials_Nonlinear -
135:
136: steps - number of time steps
137: time - time that integrater is terminated.
138: */
139: PetscErrorCode TSStep_Sundials_Nonlinear(TS ts,int *steps,double *time)
140: {
141: TS_Sundials *cvode = (TS_Sundials*)ts->data;
142: Vec sol = ts->vec_sol;
144: int i,max_steps = ts->max_steps,flag;
145: long int its;
146: realtype t,tout;
147: PetscScalar *y_data;
148: void *mem;
151: /*
152: Call CVodeCreate to create the solver memory:
153: CV_ADAMS specifies the Adams Method
154: CV_FUNCTIONAL specifies functional iteration
155: A pointer to the integrator memory is returned and stored in cvode_mem.
156: */
157: mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
158: if (!mem) SETERRQ(1,"CVodeCreate() fails");
159: flag = CVodeSetFdata(mem,ts);
160: if (flag) SETERRQ(1,"CVodeSetFdata() fails");
162: /*
163: Call CVodeMalloc to initialize the integrator memory:
164: mem is the pointer to the integrator memory returned by CVodeCreate
165: f is the user's right hand side function in y'=f(t,y)
166: T0 is the initial time
167: u is the initial dependent variable vector
168: CV_SS specifies scalar relative and absolute tolerances
169: reltol is the relative tolerance
170: &abstol is a pointer to the scalar absolute tolerance
171: */
172: flag = CVodeMalloc(mem,TSFunction_Sundials,ts->ptime,cvode->y,CV_SS,cvode->reltol,&cvode->abstol);
173: if (flag) SETERRQ(1,"CVodeMalloc() fails");
175: /* initialize the number of steps */
176: *steps = -ts->steps;
177: TSMonitor(ts,ts->steps,ts->ptime,sol);
179: /* call CVSpgmr to use GMRES as the linear solver. */
180: /* setup the ode integrator with the given preconditioner */
181: flag = CVSpgmr(mem,PREC_LEFT,0);
182: if (flag) SETERRQ(1,"CVSpgmr() fails");
183:
184: flag = CVSpilsSetGSType(mem, MODIFIED_GS);
185: if (flag) SETERRQ(1,"CVSpgmrSetGSType() fails");
187: /* Set preconditioner setup and solve routines Precond and PSolve,
188: and the pointer to the user-defined block data */
189: flag = CVSpilsSetPreconditioner(mem,TSPrecond_Sundials,TSPSolve_Sundials,ts);
190: if (flag) SETERRQ(1,"CVSpgmrSetPreconditioner() fails");
192: tout = ts->max_time;
193: VecGetArray(ts->vec_sol,&y_data);
194: N_VSetArrayPointer((realtype *)y_data,cvode->y);
195: VecRestoreArray(ts->vec_sol,PETSC_NULL);
196: for (i = 0; i < max_steps; i++) {
197: if (ts->ptime >= ts->max_time) break;
198: CVode(mem,tout,cvode->y,&t,CV_ONE_STEP);
199: CVodeGetNumNonlinSolvIters(mem,&its);
200: cvode->nonlinear_solves += its;
202: if (t > ts->max_time && cvode->exact_final_time) {
203: /* interpolate to final requested time */
204: CVodeGetDky(mem,tout,0,cvode->y);
205: t = tout;
206: }
207: ts->time_step = t - ts->ptime;
208: ts->ptime = t;
210: /* copy the solution from cvode->y to cvode->update and sol */
211: VecPlaceArray(cvode->w1,y_data);
212: VecCopy(cvode->w1,cvode->update);
213: VecResetArray(cvode->w1);
214: VecCopy(cvode->update,sol);
215: CVodeGetNumNonlinSolvIters(mem,&its);
216: ts->nonlinear_its = its;
217: CVSpilsGetNumLinIters(mem, &its);
218: ts->linear_its = its;
219: ts->steps++;
220: TSMonitor(ts,ts->steps,t,sol);
221: }
222: CVodeFree(&mem);
223: *steps += ts->steps;
224: *time = t;
225: return(0);
226: }
230: PetscErrorCode TSDestroy_Sundials(TS ts)
231: {
232: TS_Sundials *cvode = (TS_Sundials*)ts->data;
236: if (cvode->pmat) {MatDestroy(cvode->pmat);}
237: if (cvode->pc) {PCDestroy(cvode->pc);}
238: if (cvode->update) {VecDestroy(cvode->update);}
239: if (cvode->func) {VecDestroy(cvode->func);}
240: if (cvode->rhs) {VecDestroy(cvode->rhs);}
241: if (cvode->w1) {VecDestroy(cvode->w1);}
242: if (cvode->w2) {VecDestroy(cvode->w2);}
243: MPI_Comm_free(&(cvode->comm_sundials));
244: PetscFree(cvode);
245: return(0);
246: }
250: PetscErrorCode TSSetUp_Sundials_Nonlinear(TS ts)
251: {
252: TS_Sundials *cvode = (TS_Sundials*)ts->data;
254: int glosize,locsize,i;
255: PetscScalar *y_data,*parray;
258: PCSetFromOptions(cvode->pc);
259: /* get the vector size */
260: VecGetSize(ts->vec_sol,&glosize);
261: VecGetLocalSize(ts->vec_sol,&locsize);
263: /* allocate the memory for N_Vec y */
264: cvode->y = N_VNew_Parallel(cvode->comm_sundials,locsize,glosize);
265: if (!cvode->y) SETERRQ(1,"cvode->y is not allocated");
267: /* initialize N_Vec y */
268: VecGetArray(ts->vec_sol,&parray);
269: y_data = (PetscScalar *) N_VGetArrayPointer(cvode->y);
270: for (i = 0; i < locsize; i++) y_data[i] = parray[i];
271: /*PetscMemcpy(y_data,parray,locsize*sizeof(PETSC_SCALAR)); */
272: VecRestoreArray(ts->vec_sol,PETSC_NULL);
273: VecDuplicate(ts->vec_sol,&cvode->update);
274: VecDuplicate(ts->vec_sol,&cvode->func);
275: PetscLogObjectParent(ts,cvode->update);
276: PetscLogObjectParent(ts,cvode->func);
278: /*
279: Create work vectors for the TSPSolve_Sundials() routine. Note these are
280: allocated with zero space arrays because the actual array space is provided
281: by Sundials and set using VecPlaceArray().
282: */
283: VecCreateMPIWithArray(ts->comm,locsize,PETSC_DECIDE,0,&cvode->w1);
284: VecCreateMPIWithArray(ts->comm,locsize,PETSC_DECIDE,0,&cvode->w2);
285: PetscLogObjectParent(ts,cvode->w1);
286: PetscLogObjectParent(ts,cvode->w2);
287: return(0);
288: }
292: PetscErrorCode TSSetFromOptions_Sundials_Nonlinear(TS ts)
293: {
294: TS_Sundials *cvode = (TS_Sundials*)ts->data;
296: int indx;
297: const char *btype[] = {"bdf","adams"},*otype[] = {"modified","unmodified"};
298: PetscTruth flag;
301: PetscOptionsHead("SUNDIALS ODE solver options");
302: PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",btype,2,"bdf",&indx,&flag);
303: if (flag) {
304: TSSundialsSetType(ts,(TSSundialsType)indx);
305: }
306: PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",otype,2,"unmodified",&indx,&flag);
307: if (flag) {
308: TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);
309: }
310: PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,PETSC_NULL);
311: PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,PETSC_NULL);
312: PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,&flag);
313: PetscOptionsInt("-ts_sundials_gmres_restart","Number of GMRES orthogonalization directions","TSSundialsSetGMRESRestart",cvode->restart,&cvode->restart,&flag);
314: PetscOptionsName("-ts_sundials_exact_final_time","Allow SUNDIALS to stop near the final time, not exactly on it","TSSundialsSetExactFinalTime",&flag);
315: if (flag) cvode->exact_final_time = PETSC_TRUE;
316: PetscOptionsTail();
317: return(0);
318: }
322: PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
323: {
324: TS_Sundials *cvode = (TS_Sundials*)ts->data;
326: char *type;
327: char atype[] = "Adams";
328: char btype[] = "BDF: backward differentiation formula";
329: PetscTruth iascii,isstring;
332: if (cvode->cvode_type == SUNDIALS_ADAMS) {type = atype;}
333: else {type = btype;}
335: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
336: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
337: if (iascii) {
338: PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");
339: PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);
340: PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",cvode->abstol,cvode->reltol);
341: PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",cvode->linear_tol);
342: PetscViewerASCIIPrintf(viewer,"Sundials GMRES max iterations (same as restart in SUNDIALS) %D\n",cvode->restart);
343: if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
344: PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");
345: } else {
346: PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");
347: }
348: } else if (isstring) {
349: PetscViewerStringSPrintf(viewer,"Sundials type %s",type);
350: } else {
351: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by TS Sundials",((PetscObject)viewer)->type_name);
352: }
353: PetscViewerASCIIPushTab(viewer);
354: PCView(cvode->pc,viewer);
355: PetscViewerASCIIPopTab(viewer);
356: return(0);
357: }
360: /* --------------------------------------------------------------------------*/
364: PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsType type)
365: {
366: TS_Sundials *cvode = (TS_Sundials*)ts->data;
367:
369: cvode->cvode_type = type;
370: return(0);
371: }
377: PetscErrorCode TSSundialsSetGMRESRestart_Sundials(TS ts,int restart)
378: {
379: TS_Sundials *cvode = (TS_Sundials*)ts->data;
380:
382: cvode->restart = restart;
383: return(0);
384: }
390: PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
391: {
392: TS_Sundials *cvode = (TS_Sundials*)ts->data;
393:
395: cvode->linear_tol = tol;
396: return(0);
397: }
403: PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
404: {
405: TS_Sundials *cvode = (TS_Sundials*)ts->data;
406:
408: cvode->gtype = type;
409: return(0);
410: }
416: PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
417: {
418: TS_Sundials *cvode = (TS_Sundials*)ts->data;
419:
421: if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
422: if (rel != PETSC_DECIDE) cvode->reltol = rel;
423: return(0);
424: }
430: PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc)
431: {
432: TS_Sundials *cvode = (TS_Sundials*)ts->data;
435: *pc = cvode->pc;
436: return(0);
437: }
443: PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
444: {
445: TS_Sundials *cvode = (TS_Sundials*)ts->data;
446:
448: if (nonlin) *nonlin = cvode->nonlinear_solves;
449: if (lin) *lin = cvode->linear_solves;
450: return(0);
451: }
453:
457: PetscErrorCode TSSundialsSetExactFinalTime_Sundials(TS ts,PetscTruth s)
458: {
459: TS_Sundials *cvode = (TS_Sundials*)ts->data;
460:
462: cvode->exact_final_time = s;
463: return(0);
464: }
466: /* -------------------------------------------------------------------------------------------*/
470: /*@C
471: TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.
473: Not Collective
475: Input parameters:
476: . ts - the time-step context
478: Output Parameters:
479: + nonlin - number of nonlinear iterations
480: - lin - number of linear iterations
482: Level: advanced
484: Notes:
485: These return the number since the creation of the TS object
487: .keywords: non-linear iterations, linear iterations
489: .seealso: TSSundialsSetType(), TSSundialsSetGMRESRestart(),
490: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
491: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
492: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
493: TSSundialsSetExactFinalTime()
495: @*/
496: PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
497: {
498: PetscErrorCode ierr,(*f)(TS,int*,int*);
499:
501: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsGetIterations_C",(void (**)(void))&f);
502: if (f) {
503: (*f)(ts,nonlin,lin);
504: }
505: return(0);
506: }
510: /*@
511: TSSundialsSetType - Sets the method that Sundials will use for integration.
513: Collective on TS
515: Input parameters:
516: + ts - the time-step context
517: - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF
519: Level: intermediate
521: .keywords: Adams, backward differentiation formula
523: .seealso: TSSundialsGetIterations(), TSSundialsSetGMRESRestart(),
524: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
525: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
526: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
527: TSSundialsSetExactFinalTime()
528: @*/
529: PetscErrorCode TSSundialsSetType(TS ts,TSSundialsType type)
530: {
531: PetscErrorCode ierr,(*f)(TS,TSSundialsType);
532:
534: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsSetType_C",(void (**)(void))&f);
535: if (f) {
536: (*f)(ts,type);
537: }
538: return(0);
539: }
543: /*@
544: TSSundialsSetGMRESRestart - Sets the dimension of the Krylov space used by
545: GMRES in the linear solver in SUNDIALS. SUNDIALS DOES NOT use restarted GMRES so
546: this is ALSO the maximum number of GMRES steps that will be used.
548: Collective on TS
550: Input parameters:
551: + ts - the time-step context
552: - restart - number of direction vectors (the restart size).
554: Level: advanced
556: .keywords: GMRES, restart
558: .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
559: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
560: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
561: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
562: TSSundialsSetExactFinalTime()
564: @*/
565: PetscErrorCode TSSundialsSetGMRESRestart(TS ts,int restart)
566: {
567: PetscErrorCode ierr,(*f)(TS,int);
570: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsSetGMRESRestart_C",(void (**)(void))&f);
571: if (f) {
572: (*f)(ts,restart);
573: }
574: return(0);
575: }
579: /*@
580: TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
581: system by SUNDIALS.
583: Collective on TS
585: Input parameters:
586: + ts - the time-step context
587: - tol - the factor by which the tolerance on the nonlinear solver is
588: multiplied to get the tolerance on the linear solver, .05 by default.
590: Level: advanced
592: .keywords: GMRES, linear convergence tolerance, SUNDIALS
594: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
595: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
596: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
597: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
598: TSSundialsSetExactFinalTime()
600: @*/
601: PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol)
602: {
603: PetscErrorCode ierr,(*f)(TS,double);
604:
606: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",(void (**)(void))&f);
607: if (f) {
608: (*f)(ts,tol);
609: }
610: return(0);
611: }
615: /*@
616: TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
617: in GMRES method by SUNDIALS linear solver.
619: Collective on TS
621: Input parameters:
622: + ts - the time-step context
623: - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS
625: Level: advanced
627: .keywords: Sundials, orthogonalization
629: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
630: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(),
631: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
632: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
633: TSSundialsSetExactFinalTime()
635: @*/
636: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
637: {
638: PetscErrorCode ierr,(*f)(TS,TSSundialsGramSchmidtType);
639:
641: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",(void (**)(void))&f);
642: if (f) {
643: (*f)(ts,type);
644: }
645: return(0);
646: }
650: /*@
651: TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
652: Sundials for error control.
654: Collective on TS
656: Input parameters:
657: + ts - the time-step context
658: . aabs - the absolute tolerance
659: - rel - the relative tolerance
661: See the Cvode/Sundials users manual for exact details on these parameters. Essentially
662: these regulate the size of the error for a SINGLE timestep.
664: Level: intermediate
666: .keywords: Sundials, tolerance
668: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
669: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
670: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
671: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
672: TSSundialsSetExactFinalTime()
674: @*/
675: PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel)
676: {
677: PetscErrorCode ierr,(*f)(TS,double,double);
678:
680: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsSetTolerance_C",(void (**)(void))&f);
681: if (f) {
682: (*f)(ts,aabs,rel);
683: }
684: return(0);
685: }
689: /*@
690: TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.
692: Input Parameter:
693: . ts - the time-step context
695: Output Parameter:
696: . pc - the preconditioner context
698: Level: advanced
700: .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
701: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
702: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
703: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
704: @*/
705: PetscErrorCode TSSundialsGetPC(TS ts,PC *pc)
706: {
707: PetscErrorCode ierr,(*f)(TS,PC *);
710: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsGetPC_C",(void (**)(void))&f);
711: if (f) {
712: (*f)(ts,pc);
713: } else {
714: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"TS must be of Sundials type to extract the PC");
715: }
717: return(0);
718: }
722: /*@
723: TSSundialsSetExactFinalTime - Determines if Sundials interpolates solution to the
724: exact final time requested by the user or just returns it at the final time
725: it computed. (Defaults to true).
727: Input Parameter:
728: + ts - the time-step context
729: - ft - PETSC_TRUE if interpolates, else PETSC_FALSE
731: Level: beginner
733: .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
734: TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
735: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESRestart(),
736: TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
737: @*/
738: PetscErrorCode TSSundialsSetExactFinalTime(TS ts,PetscTruth ft)
739: {
740: PetscErrorCode ierr,(*f)(TS,PetscTruth);
743: PetscObjectQueryFunction((PetscObject)ts,"TSSundialsSetExactFinalTime_C",(void (**)(void))&f);
744: if (f) {
745: (*f)(ts,ft);
746: }
747: return(0);
748: }
750: /* -------------------------------------------------------------------------------------------*/
751: /*MC
752: TS_Sundials - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)
754: Options Database:
755: + -ts_sundials_type <bdf,adams>
756: . -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
757: . -ts_sundials_atol <tol> - Absolute tolerance for convergence
758: . -ts_sundials_rtol <tol> - Relative tolerance for convergence
759: . -ts_sundials_linear_tolerance <tol>
760: . -ts_sundials_gmres_restart <restart> - Number of GMRES orthogonalization directions
761: - -ts_sundials_not_exact_final_time -Allow SUNDIALS to stop near the final time, not exactly on it
763: Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply
764: only PETSc PC options
766: Level: beginner
768: .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetGMRESRestart(), TSSundialsSetLinearTolerance(),
769: TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSundialsSetExactFinalTime()
771: M*/
775: PetscErrorCode TSCreate_Sundials(TS ts)
776: {
777: TS_Sundials *cvode;
781: ts->ops->destroy = TSDestroy_Sundials;
782: ts->ops->view = TSView_Sundials;
784: if (ts->problem_type != TS_NONLINEAR) {
785: SETERRQ(PETSC_ERR_SUP,"Only support for nonlinear problems");
786: }
787: ts->ops->setup = TSSetUp_Sundials_Nonlinear;
788: ts->ops->step = TSStep_Sundials_Nonlinear;
789: ts->ops->setfromoptions = TSSetFromOptions_Sundials_Nonlinear;
791: PetscNew(TS_Sundials,&cvode);
792: PCCreate(ts->comm,&cvode->pc);
793: PetscLogObjectParent(ts,cvode->pc);
794: ts->data = (void*)cvode;
795: cvode->cvode_type = CV_BDF;
796: cvode->gtype = SUNDIALS_UNMODIFIED_GS;
797: cvode->restart = 5;
798: cvode->linear_tol = .05;
800: cvode->exact_final_time = PETSC_FALSE;
802: MPI_Comm_dup(ts->comm,&(cvode->comm_sundials));
803: /* set tolerance for Sundials */
804: cvode->abstol = 1e-6;
805: cvode->reltol = 1e-6;
807: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetType_C","TSSundialsSetType_Sundials",
808: TSSundialsSetType_Sundials);
809: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetGMRESRestart_C",
810: "TSSundialsSetGMRESRestart_Sundials",
811: TSSundialsSetGMRESRestart_Sundials);
812: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetLinearTolerance_C",
813: "TSSundialsSetLinearTolerance_Sundials",
814: TSSundialsSetLinearTolerance_Sundials);
815: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",
816: "TSSundialsSetGramSchmidtType_Sundials",
817: TSSundialsSetGramSchmidtType_Sundials);
818: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetTolerance_C",
819: "TSSundialsSetTolerance_Sundials",
820: TSSundialsSetTolerance_Sundials);
821: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetPC_C",
822: "TSSundialsGetPC_Sundials",
823: TSSundialsGetPC_Sundials);
824: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsGetIterations_C",
825: "TSSundialsGetIterations_Sundials",
826: TSSundialsGetIterations_Sundials);
827: PetscObjectComposeFunctionDynamic((PetscObject)ts,"TSSundialsSetExactFinalTime_C",
828: "TSSundialsSetExactFinalTime_Sundials",
829: TSSundialsSetExactFinalTime_Sundials);
830: return(0);
831: }