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.
5
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*/
9
10 /*
11 TSPrecond_Sundials - function that we provide to SUNDIALS to
12 evaluate the preconditioner.
13 */
TSPrecond_Sundials(realtype tn,N_Vector y,N_Vector fy,booleantype jok,booleantype * jcurPtr,realtype _gamma,void * P_data,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)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;
20 PetscErrorCode ierr;
21 Mat J,P;
22 Vec yy = cvode->w1,yydot = cvode->ydot;
23 PetscReal gm = (PetscReal)_gamma;
24 PetscScalar *y_data;
25
26 PetscFunctionBegin;
27 ierr = TSGetIJacobian(ts,&J,&P,NULL,NULL);CHKERRQ(ierr);
28 y_data = (PetscScalar*) N_VGetArrayPointer(y);
29 ierr = VecPlaceArray(yy,y_data);CHKERRQ(ierr);
30 ierr = VecZeroEntries(yydot);CHKERRQ(ierr); /* 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 ierr = TSComputeIJacobian(ts,ts->ptime,yy,yydot,1/gm,J,P,PETSC_FALSE);CHKERRQ(ierr);
33 ierr = VecResetArray(yy);CHKERRQ(ierr);
34 ierr = MatScale(P,gm);CHKERRQ(ierr); /* turn into I-gm*Jrest, J is not used by Sundials */
35 *jcurPtr = TRUE;
36 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr);
37 ierr = PCSetOperators(pc,J,P);CHKERRQ(ierr);
38 PetscFunctionReturn(0);
39 }
40
41 /*
42 TSPSolve_Sundials - routine that we provide to Sundials that applies the preconditioner.
43 */
TSPSolve_Sundials(realtype tn,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,realtype _gamma,realtype delta,int lr,void * P_data,N_Vector vtemp)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;
51 PetscErrorCode ierr;
52 PetscScalar *r_data,*z_data;
53
54 PetscFunctionBegin;
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 ierr = VecPlaceArray(rr,r_data);CHKERRQ(ierr);
59 ierr = VecPlaceArray(zz,z_data);CHKERRQ(ierr);
60
61 /* Solve the Px=r and put the result in zz */
62 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr);
63 ierr = PCApply(pc,rr,zz);CHKERRQ(ierr);
64 ierr = VecResetArray(rr);CHKERRQ(ierr);
65 ierr = VecResetArray(zz);CHKERRQ(ierr);
66 PetscFunctionReturn(0);
67 }
68
69 /*
70 TSFunction_Sundials - routine that we provide to Sundials that applies the right hand side.
71 */
TSFunction_Sundials(realtype t,N_Vector y,N_Vector ydot,void * ctx)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;
82 PetscErrorCode ierr;
83
84 PetscFunctionBegin;
85 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr);
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 ierr = VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr);
90 ierr = VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr);
91
92 /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
93 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
94 ierr = DMGetDMTS(dm,&tsdm);CHKERRQ(ierr);
95 ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr);
96 if (!ifunction) {
97 ierr = TSComputeRHSFunction(ts,t,yy,yyd);CHKERRQ(ierr);
98 } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
99 ierr = VecZeroEntries(yydot);CHKERRQ(ierr);
100 ierr = TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr);
101 ierr = VecScale(yyd,-1.);CHKERRQ(ierr);
102 }
103 ierr = VecResetArray(yy);CHKERRABORT(comm,ierr);
104 ierr = VecResetArray(yyd);CHKERRABORT(comm,ierr);
105 PetscFunctionReturn(0);
106 }
107
108 /*
109 TSStep_Sundials - Calls Sundials to integrate the ODE.
110 */
TSStep_Sundials(TS ts)111 PetscErrorCode TSStep_Sundials(TS ts)
112 {
113 TS_Sundials *cvode = (TS_Sundials*)ts->data;
114 PetscErrorCode ierr;
115 PetscInt flag;
116 long int nits,lits,nsteps;
117 realtype t,tout;
118 PetscScalar *y_data;
119 void *mem;
120
121 PetscFunctionBegin;
122 mem = cvode->mem;
123 tout = ts->max_time;
124 ierr = VecGetArray(ts->vec_sol,&y_data);CHKERRQ(ierr);
125 N_VSetArrayPointer((realtype*)y_data,cvode->y);
126 ierr = VecRestoreArray(ts->vec_sol,NULL);CHKERRQ(ierr);
127
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);
134
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 ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr);
146 ierr = CVodeGetCurrentTime(mem,&tcur);CHKERRQ(ierr);
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 }
186
187 /* log inner nonlinear and linear iterations */
188 ierr = CVodeGetNumNonlinSolvIters(mem,&nits);CHKERRQ(ierr);
189 ierr = CVSpilsGetNumLinIters(mem,&lits);CHKERRQ(ierr);
190 ts->snes_its += nits; ts->ksp_its = lits;
191
192 /* copy the solution from cvode->y to cvode->update and sol */
193 ierr = VecPlaceArray(cvode->w1,y_data);CHKERRQ(ierr);
194 ierr = VecCopy(cvode->w1,cvode->update);CHKERRQ(ierr);
195 ierr = VecResetArray(cvode->w1);CHKERRQ(ierr);
196 ierr = VecCopy(cvode->update,ts->vec_sol);CHKERRQ(ierr);
197
198 ts->time_step = t - ts->ptime;
199 ts->ptime = t;
200
201 ierr = CVodeGetNumSteps(mem,&nsteps);CHKERRQ(ierr);
202 if (!cvode->monitorstep) ts->steps += nsteps - 1; /* TSStep() increments the step counter by one */
203 PetscFunctionReturn(0);
204 }
205
TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)206 static PetscErrorCode TSInterpolate_Sundials(TS ts,PetscReal t,Vec X)
207 {
208 TS_Sundials *cvode = (TS_Sundials*)ts->data;
209 N_Vector y;
210 PetscErrorCode ierr;
211 PetscScalar *x_data;
212 PetscInt glosize,locsize;
213
214 PetscFunctionBegin;
215 /* get the vector size */
216 ierr = VecGetSize(X,&glosize);CHKERRQ(ierr);
217 ierr = VecGetLocalSize(X,&locsize);CHKERRQ(ierr);
218 ierr = VecGetArray(X,&x_data);CHKERRQ(ierr);
219
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");
223
224 ierr = CVodeGetDky(cvode->mem,t,0,y);CHKERRQ(ierr);
225 ierr = VecRestoreArray(X,&x_data);CHKERRQ(ierr);
226 PetscFunctionReturn(0);
227 }
228
TSReset_Sundials(TS ts)229 PetscErrorCode TSReset_Sundials(TS ts)
230 {
231 TS_Sundials *cvode = (TS_Sundials*)ts->data;
232 PetscErrorCode ierr;
233
234 PetscFunctionBegin;
235 ierr = VecDestroy(&cvode->update);CHKERRQ(ierr);
236 ierr = VecDestroy(&cvode->ydot);CHKERRQ(ierr);
237 ierr = VecDestroy(&cvode->w1);CHKERRQ(ierr);
238 ierr = VecDestroy(&cvode->w2);CHKERRQ(ierr);
239 if (cvode->mem) CVodeFree(&cvode->mem);
240 PetscFunctionReturn(0);
241 }
242
TSDestroy_Sundials(TS ts)243 PetscErrorCode TSDestroy_Sundials(TS ts)
244 {
245 TS_Sundials *cvode = (TS_Sundials*)ts->data;
246 PetscErrorCode ierr;
247
248 PetscFunctionBegin;
249 ierr = TSReset_Sundials(ts);CHKERRQ(ierr);
250 ierr = MPI_Comm_free(&(cvode->comm_sundials));CHKERRQ(ierr);
251 ierr = PetscFree(ts->data);CHKERRQ(ierr);
252 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",NULL);CHKERRQ(ierr);
253 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",NULL);CHKERRQ(ierr);
254 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",NULL);CHKERRQ(ierr);
255 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",NULL);CHKERRQ(ierr);
256 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",NULL);CHKERRQ(ierr);
257 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",NULL);CHKERRQ(ierr);
258 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",NULL);CHKERRQ(ierr);
259 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",NULL);CHKERRQ(ierr);
260 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",NULL);CHKERRQ(ierr);
261 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",NULL);CHKERRQ(ierr);
262 PetscFunctionReturn(0);
263 }
264
TSSetUp_Sundials(TS ts)265 PetscErrorCode TSSetUp_Sundials(TS ts)
266 {
267 TS_Sundials *cvode = (TS_Sundials*)ts->data;
268 PetscErrorCode ierr;
269 PetscInt glosize,locsize,i,flag;
270 PetscScalar *y_data,*parray;
271 void *mem;
272 PC pc;
273 PCType pctype;
274 PetscBool pcnone;
275
276 PetscFunctionBegin;
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");
278
279 /* get the vector size */
280 ierr = VecGetSize(ts->vec_sol,&glosize);CHKERRQ(ierr);
281 ierr = VecGetLocalSize(ts->vec_sol,&locsize);CHKERRQ(ierr);
282
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");
286
287 /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
288 ierr = VecGetArray(ts->vec_sol,&parray);CHKERRQ(ierr);
289 y_data = (PetscScalar*) N_VGetArrayPointer(cvode->y);
290 for (i = 0; i < locsize; i++) y_data[i] = parray[i];
291 ierr = VecRestoreArray(ts->vec_sol,NULL);CHKERRQ(ierr);
292
293 ierr = VecDuplicate(ts->vec_sol,&cvode->update);CHKERRQ(ierr);
294 ierr = VecDuplicate(ts->vec_sol,&cvode->ydot);CHKERRQ(ierr);
295 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->update);CHKERRQ(ierr);
296 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->ydot);CHKERRQ(ierr);
297
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 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,NULL,&cvode->w1);CHKERRQ(ierr);
304 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts),1,locsize,PETSC_DECIDE,NULL,&cvode->w2);CHKERRQ(ierr);
305 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w1);CHKERRQ(ierr);
306 ierr = PetscLogObjectParent((PetscObject)ts,(PetscObject)cvode->w2);CHKERRQ(ierr);
307
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;
312
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");
316
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 }
332
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);
338
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);
342
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 }
348
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);
352
353 /* call CVSpgmr to use GMRES as the linear solver. */
354 /* setup the ode integrator with the given preconditioner */
355 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr);
356 ierr = PCGetType(pc,&pctype);CHKERRQ(ierr);
357 ierr = PetscObjectTypeCompare((PetscObject)pc,PCNONE,&pcnone);CHKERRQ(ierr);
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);
364
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 }
370
371 flag = CVSpilsSetGSType(mem, MODIFIED_GS);
372 if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CVSpgmrSetGSType() fails, flag %d",flag);
373 PetscFunctionReturn(0);
374 }
375
376 /* type of CVODE linear multistep method */
377 const char *const TSSundialsLmmTypes[] = {"","ADAMS","BDF","TSSundialsLmmType","SUNDIALS_",NULL};
378 /* type of G-S orthogonalization used by CVODE linear solver */
379 const char *const TSSundialsGramSchmidtTypes[] = {"","MODIFIED","CLASSICAL","TSSundialsGramSchmidtType","SUNDIALS_",NULL};
380
TSSetFromOptions_Sundials(PetscOptionItems * PetscOptionsObject,TS ts)381 PetscErrorCode TSSetFromOptions_Sundials(PetscOptionItems *PetscOptionsObject,TS ts)
382 {
383 TS_Sundials *cvode = (TS_Sundials*)ts->data;
384 PetscErrorCode ierr;
385 int indx;
386 PetscBool flag;
387 PC pc;
388
389 PetscFunctionBegin;
390 ierr = PetscOptionsHead(PetscOptionsObject,"SUNDIALS ODE solver options");CHKERRQ(ierr);
391 ierr = PetscOptionsEList("-ts_sundials_type","Scheme","TSSundialsSetType",TSSundialsLmmTypes,3,TSSundialsLmmTypes[cvode->cvode_type],&indx,&flag);CHKERRQ(ierr);
392 if (flag) {
393 ierr = TSSundialsSetType(ts,(TSSundialsLmmType)indx);CHKERRQ(ierr);
394 }
395 ierr = PetscOptionsEList("-ts_sundials_gramschmidt_type","Type of orthogonalization","TSSundialsSetGramSchmidtType",TSSundialsGramSchmidtTypes,3,TSSundialsGramSchmidtTypes[cvode->gtype],&indx,&flag);CHKERRQ(ierr);
396 if (flag) {
397 ierr = TSSundialsSetGramSchmidtType(ts,(TSSundialsGramSchmidtType)indx);CHKERRQ(ierr);
398 }
399 ierr = PetscOptionsReal("-ts_sundials_atol","Absolute tolerance for convergence","TSSundialsSetTolerance",cvode->abstol,&cvode->abstol,NULL);CHKERRQ(ierr);
400 ierr = PetscOptionsReal("-ts_sundials_rtol","Relative tolerance for convergence","TSSundialsSetTolerance",cvode->reltol,&cvode->reltol,NULL);CHKERRQ(ierr);
401 ierr = PetscOptionsReal("-ts_sundials_mindt","Minimum step size","TSSundialsSetMinTimeStep",cvode->mindt,&cvode->mindt,NULL);CHKERRQ(ierr);
402 ierr = PetscOptionsReal("-ts_sundials_maxdt","Maximum step size","TSSundialsSetMaxTimeStep",cvode->maxdt,&cvode->maxdt,NULL);CHKERRQ(ierr);
403 ierr = PetscOptionsReal("-ts_sundials_linear_tolerance","Convergence tolerance for linear solve","TSSundialsSetLinearTolerance",cvode->linear_tol,&cvode->linear_tol,NULL);CHKERRQ(ierr);
404 ierr = PetscOptionsInt("-ts_sundials_maxord","Max Order for BDF/Adams method","TSSundialsSetMaxOrd",cvode->maxord,&cvode->maxord,NULL);CHKERRQ(ierr);
405 ierr = PetscOptionsInt("-ts_sundials_maxl","Max dimension of the Krylov subspace","TSSundialsSetMaxl",cvode->maxl,&cvode->maxl,NULL);CHKERRQ(ierr);
406 ierr = PetscOptionsBool("-ts_sundials_monitor_steps","Monitor SUNDIALS internal steps","TSSundialsMonitorInternalSteps",cvode->monitorstep,&cvode->monitorstep,NULL);CHKERRQ(ierr);
407 ierr = PetscOptionsTail();CHKERRQ(ierr);
408 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr);
409 ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
410 PetscFunctionReturn(0);
411 }
412
TSView_Sundials(TS ts,PetscViewer viewer)413 PetscErrorCode TSView_Sundials(TS ts,PetscViewer viewer)
414 {
415 TS_Sundials *cvode = (TS_Sundials*)ts->data;
416 PetscErrorCode ierr;
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;
425
426 PetscFunctionBegin;
427 if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
428 else type = btype;
429
430 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
431 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr);
432 if (iascii) {
433 ierr = PetscViewerASCIIPrintf(viewer,"Sundials integrater does not use SNES!\n");CHKERRQ(ierr);
434 ierr = PetscViewerASCIIPrintf(viewer,"Sundials integrater type %s\n",type);CHKERRQ(ierr);
435 ierr = PetscViewerASCIIPrintf(viewer,"Sundials maxord %D\n",cvode->maxord);CHKERRQ(ierr);
436 ierr = PetscViewerASCIIPrintf(viewer,"Sundials abs tol %g rel tol %g\n",(double)cvode->abstol,(double)cvode->reltol);CHKERRQ(ierr);
437 ierr = PetscViewerASCIIPrintf(viewer,"Sundials linear solver tolerance factor %g\n",(double)cvode->linear_tol);CHKERRQ(ierr);
438 ierr = PetscViewerASCIIPrintf(viewer,"Sundials max dimension of Krylov subspace %D\n",cvode->maxl);CHKERRQ(ierr);
439 if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
440 ierr = PetscViewerASCIIPrintf(viewer,"Sundials using modified Gram-Schmidt for orthogonalization in GMRES\n");CHKERRQ(ierr);
441 } else {
442 ierr = PetscViewerASCIIPrintf(viewer,"Sundials using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n");CHKERRQ(ierr);
443 }
444 if (cvode->mindt > 0) {ierr = PetscViewerASCIIPrintf(viewer,"Sundials minimum time step %g\n",(double)cvode->mindt);CHKERRQ(ierr);}
445 if (cvode->maxdt > 0) {ierr = PetscViewerASCIIPrintf(viewer,"Sundials maximum time step %g\n",(double)cvode->maxdt);CHKERRQ(ierr);}
446
447 /* Outputs from CVODE, CVSPILS */
448 ierr = CVodeGetTolScaleFactor(cvode->mem,&tolsfac);CHKERRQ(ierr);
449 ierr = PetscViewerASCIIPrintf(viewer,"Sundials suggested factor for tolerance scaling %g\n",tolsfac);CHKERRQ(ierr);
450 ierr = CVodeGetIntegratorStats(cvode->mem,&nsteps,&nfevals,
451 &nlinsetups,&nfails,&qlast,&qcur,
452 &hinused,&hlast,&hcur,&tcur);CHKERRQ(ierr);
453 ierr = PetscViewerASCIIPrintf(viewer,"Sundials cumulative number of internal steps %D\n",nsteps);CHKERRQ(ierr);
454 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to rhs function %D\n",nfevals);CHKERRQ(ierr);
455 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of calls to linear solver setup function %D\n",nlinsetups);CHKERRQ(ierr);
456 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of error test failures %D\n",nfails);CHKERRQ(ierr);
457
458 ierr = CVodeGetNonlinSolvStats(cvode->mem,&its,&nfails);CHKERRQ(ierr);
459 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear solver iterations %D\n",its);CHKERRQ(ierr);
460 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of nonlinear convergence failure %D\n",nfails);CHKERRQ(ierr);
461
462 ierr = CVSpilsGetNumLinIters(cvode->mem, &its);CHKERRQ(ierr); /* its = no. of calls to TSPrecond_Sundials() */
463 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of linear iterations %D\n",its);CHKERRQ(ierr);
464 ierr = CVSpilsGetNumConvFails(cvode->mem,&itmp);CHKERRQ(ierr);
465 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of linear convergence failures %D\n",itmp);CHKERRQ(ierr);
466
467 ierr = TSSundialsGetPC(ts,&pc);CHKERRQ(ierr);
468 ierr = PCView(pc,viewer);CHKERRQ(ierr);
469 ierr = CVSpilsGetNumPrecEvals(cvode->mem,&itmp);CHKERRQ(ierr);
470 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner evaluations %D\n",itmp);CHKERRQ(ierr);
471 ierr = CVSpilsGetNumPrecSolves(cvode->mem,&itmp);CHKERRQ(ierr);
472 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of preconditioner solves %D\n",itmp);CHKERRQ(ierr);
473
474 ierr = CVSpilsGetNumJtimesEvals(cvode->mem,&itmp);CHKERRQ(ierr);
475 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of Jacobian-vector product evaluations %D\n",itmp);CHKERRQ(ierr);
476 ierr = CVSpilsGetNumRhsEvals(cvode->mem,&itmp);CHKERRQ(ierr);
477 ierr = PetscViewerASCIIPrintf(viewer,"Sundials no. of rhs calls for finite diff. Jacobian-vector evals %D\n",itmp);CHKERRQ(ierr);
478 } else if (isstring) {
479 ierr = PetscViewerStringSPrintf(viewer,"Sundials type %s",type);CHKERRQ(ierr);
480 }
481 PetscFunctionReturn(0);
482 }
483
484
485 /* --------------------------------------------------------------------------*/
TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)486 PetscErrorCode TSSundialsSetType_Sundials(TS ts,TSSundialsLmmType type)
487 {
488 TS_Sundials *cvode = (TS_Sundials*)ts->data;
489
490 PetscFunctionBegin;
491 cvode->cvode_type = type;
492 PetscFunctionReturn(0);
493 }
494
TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)495 PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts,PetscInt maxl)
496 {
497 TS_Sundials *cvode = (TS_Sundials*)ts->data;
498
499 PetscFunctionBegin;
500 cvode->maxl = maxl;
501 PetscFunctionReturn(0);
502 }
503
TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)504 PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts,double tol)
505 {
506 TS_Sundials *cvode = (TS_Sundials*)ts->data;
507
508 PetscFunctionBegin;
509 cvode->linear_tol = tol;
510 PetscFunctionReturn(0);
511 }
512
TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)513 PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts,TSSundialsGramSchmidtType type)
514 {
515 TS_Sundials *cvode = (TS_Sundials*)ts->data;
516
517 PetscFunctionBegin;
518 cvode->gtype = type;
519 PetscFunctionReturn(0);
520 }
521
TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)522 PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts,double aabs,double rel)
523 {
524 TS_Sundials *cvode = (TS_Sundials*)ts->data;
525
526 PetscFunctionBegin;
527 if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
528 if (rel != PETSC_DECIDE) cvode->reltol = rel;
529 PetscFunctionReturn(0);
530 }
531
TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)532 PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts,PetscReal mindt)
533 {
534 TS_Sundials *cvode = (TS_Sundials*)ts->data;
535
536 PetscFunctionBegin;
537 cvode->mindt = mindt;
538 PetscFunctionReturn(0);
539 }
540
TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)541 PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts,PetscReal maxdt)
542 {
543 TS_Sundials *cvode = (TS_Sundials*)ts->data;
544
545 PetscFunctionBegin;
546 cvode->maxdt = maxdt;
547 PetscFunctionReturn(0);
548 }
TSSundialsGetPC_Sundials(TS ts,PC * pc)549 PetscErrorCode TSSundialsGetPC_Sundials(TS ts,PC *pc)
550 {
551 SNES snes;
552 KSP ksp;
553 PetscErrorCode ierr;
554
555 PetscFunctionBegin;
556 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
557 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
558 ierr = KSPGetPC(ksp,pc);CHKERRQ(ierr);
559 PetscFunctionReturn(0);
560 }
561
TSSundialsGetIterations_Sundials(TS ts,int * nonlin,int * lin)562 PetscErrorCode TSSundialsGetIterations_Sundials(TS ts,int *nonlin,int *lin)
563 {
564 PetscFunctionBegin;
565 if (nonlin) *nonlin = ts->snes_its;
566 if (lin) *lin = ts->ksp_its;
567 PetscFunctionReturn(0);
568 }
569
TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)570 PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts,PetscBool s)
571 {
572 TS_Sundials *cvode = (TS_Sundials*)ts->data;
573
574 PetscFunctionBegin;
575 cvode->monitorstep = s;
576 PetscFunctionReturn(0);
577 }
578 /* -------------------------------------------------------------------------------------------*/
579
580 /*@C
581 TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by Sundials.
582
583 Not Collective
584
585 Input parameters:
586 . ts - the time-step context
587
588 Output Parameters:
589 + nonlin - number of nonlinear iterations
590 - lin - number of linear iterations
591
592 Level: advanced
593
594 Notes:
595 These return the number since the creation of the TS object
596
597 .seealso: TSSundialsSetType(), TSSundialsSetMaxl(),
598 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
599 TSSundialsGetIterations(), TSSundialsSetType(),
600 TSSundialsSetLinearTolerance(), TSSundialsGetPC(), TSSetExactFinalTime()
601
602 @*/
TSSundialsGetIterations(TS ts,int * nonlin,int * lin)603 PetscErrorCode TSSundialsGetIterations(TS ts,int *nonlin,int *lin)
604 {
605 PetscErrorCode ierr;
606
607 PetscFunctionBegin;
608 ierr = PetscUseMethod(ts,"TSSundialsGetIterations_C",(TS,int*,int*),(ts,nonlin,lin));CHKERRQ(ierr);
609 PetscFunctionReturn(0);
610 }
611
612 /*@
613 TSSundialsSetType - Sets the method that Sundials will use for integration.
614
615 Logically Collective on TS
616
617 Input parameters:
618 + ts - the time-step context
619 - type - one of SUNDIALS_ADAMS or SUNDIALS_BDF
620
621 Level: intermediate
622
623 .seealso: TSSundialsGetIterations(), TSSundialsSetMaxl(),
624 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
625 TSSundialsGetIterations(), TSSundialsSetType(),
626 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
627 TSSetExactFinalTime()
628 @*/
TSSundialsSetType(TS ts,TSSundialsLmmType type)629 PetscErrorCode TSSundialsSetType(TS ts,TSSundialsLmmType type)
630 {
631 PetscErrorCode ierr;
632
633 PetscFunctionBegin;
634 ierr = PetscTryMethod(ts,"TSSundialsSetType_C",(TS,TSSundialsLmmType),(ts,type));CHKERRQ(ierr);
635 PetscFunctionReturn(0);
636 }
637
638 /*@
639 TSSundialsSetMaxord - Sets the maximum order for BDF/Adams method used by SUNDIALS.
640
641 Logically Collective on TS
642
643 Input parameters:
644 + ts - the time-step context
645 - maxord - maximum order of BDF / Adams method
646
647 Level: advanced
648
649 .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
650 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
651 TSSundialsGetIterations(), TSSundialsSetType(),
652 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
653 TSSetExactFinalTime()
654
655 @*/
TSSundialsSetMaxord(TS ts,PetscInt maxord)656 PetscErrorCode TSSundialsSetMaxord(TS ts,PetscInt maxord)
657 {
658 PetscErrorCode ierr;
659
660 PetscFunctionBegin;
661 PetscValidLogicalCollectiveInt(ts,maxord,2);
662 ierr = PetscTryMethod(ts,"TSSundialsSetMaxOrd_C",(TS,PetscInt),(ts,maxord));CHKERRQ(ierr);
663 PetscFunctionReturn(0);
664 }
665
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.
670
671 Logically Collective on TS
672
673 Input parameters:
674 + ts - the time-step context
675 - maxl - number of direction vectors (the dimension of Krylov subspace).
676
677 Level: advanced
678
679 .seealso: TSSundialsGetIterations(), TSSundialsSetType(),
680 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
681 TSSundialsGetIterations(), TSSundialsSetType(),
682 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
683 TSSetExactFinalTime()
684
685 @*/
TSSundialsSetMaxl(TS ts,PetscInt maxl)686 PetscErrorCode TSSundialsSetMaxl(TS ts,PetscInt maxl)
687 {
688 PetscErrorCode ierr;
689
690 PetscFunctionBegin;
691 PetscValidLogicalCollectiveInt(ts,maxl,2);
692 ierr = PetscTryMethod(ts,"TSSundialsSetMaxl_C",(TS,PetscInt),(ts,maxl));CHKERRQ(ierr);
693 PetscFunctionReturn(0);
694 }
695
696 /*@
697 TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
698 system by SUNDIALS.
699
700 Logically Collective on TS
701
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.
706
707 Level: advanced
708
709 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
710 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
711 TSSundialsGetIterations(), TSSundialsSetType(),
712 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
713 TSSetExactFinalTime()
714
715 @*/
TSSundialsSetLinearTolerance(TS ts,double tol)716 PetscErrorCode TSSundialsSetLinearTolerance(TS ts,double tol)
717 {
718 PetscErrorCode ierr;
719
720 PetscFunctionBegin;
721 PetscValidLogicalCollectiveReal(ts,tol,2);
722 ierr = PetscTryMethod(ts,"TSSundialsSetLinearTolerance_C",(TS,double),(ts,tol));CHKERRQ(ierr);
723 PetscFunctionReturn(0);
724 }
725
726 /*@
727 TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
728 in GMRES method by SUNDIALS linear solver.
729
730 Logically Collective on TS
731
732 Input parameters:
733 + ts - the time-step context
734 - type - either SUNDIALS_MODIFIED_GS or SUNDIALS_CLASSICAL_GS
735
736 Level: advanced
737
738 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
739 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(),
740 TSSundialsGetIterations(), TSSundialsSetType(),
741 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
742 TSSetExactFinalTime()
743
744 @*/
TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)745 PetscErrorCode TSSundialsSetGramSchmidtType(TS ts,TSSundialsGramSchmidtType type)
746 {
747 PetscErrorCode ierr;
748
749 PetscFunctionBegin;
750 ierr = PetscTryMethod(ts,"TSSundialsSetGramSchmidtType_C",(TS,TSSundialsGramSchmidtType),(ts,type));CHKERRQ(ierr);
751 PetscFunctionReturn(0);
752 }
753
754 /*@
755 TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
756 Sundials for error control.
757
758 Logically Collective on TS
759
760 Input parameters:
761 + ts - the time-step context
762 . aabs - the absolute tolerance
763 - rel - the relative tolerance
764
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.
767
768 Level: intermediate
769
770 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetGMRESMaxl(),
771 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(),
772 TSSundialsGetIterations(), TSSundialsSetType(),
773 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC(),
774 TSSetExactFinalTime()
775
776 @*/
TSSundialsSetTolerance(TS ts,double aabs,double rel)777 PetscErrorCode TSSundialsSetTolerance(TS ts,double aabs,double rel)
778 {
779 PetscErrorCode ierr;
780
781 PetscFunctionBegin;
782 ierr = PetscTryMethod(ts,"TSSundialsSetTolerance_C",(TS,double,double),(ts,aabs,rel));CHKERRQ(ierr);
783 PetscFunctionReturn(0);
784 }
785
786 /*@
787 TSSundialsGetPC - Extract the PC context from a time-step context for Sundials.
788
789 Input Parameter:
790 . ts - the time-step context
791
792 Output Parameter:
793 . pc - the preconditioner context
794
795 Level: advanced
796
797 .seealso: TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
798 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
799 TSSundialsGetIterations(), TSSundialsSetType(),
800 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance()
801 @*/
TSSundialsGetPC(TS ts,PC * pc)802 PetscErrorCode TSSundialsGetPC(TS ts,PC *pc)
803 {
804 PetscErrorCode ierr;
805
806 PetscFunctionBegin;
807 ierr = PetscUseMethod(ts,"TSSundialsGetPC_C",(TS,PC*),(ts,pc));CHKERRQ(ierr);
808 PetscFunctionReturn(0);
809 }
810
811 /*@
812 TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.
813
814 Input Parameter:
815 + ts - the time-step context
816 - mindt - lowest time step if positive, negative to deactivate
817
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.
821
822 Level: beginner
823
824 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
825 @*/
TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)826 PetscErrorCode TSSundialsSetMinTimeStep(TS ts,PetscReal mindt)
827 {
828 PetscErrorCode ierr;
829
830 PetscFunctionBegin;
831 ierr = PetscTryMethod(ts,"TSSundialsSetMinTimeStep_C",(TS,PetscReal),(ts,mindt));CHKERRQ(ierr);
832 PetscFunctionReturn(0);
833 }
834
835 /*@
836 TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.
837
838 Input Parameter:
839 + ts - the time-step context
840 - maxdt - lowest time step if positive, negative to deactivate
841
842 Level: beginner
843
844 .seealso: TSSundialsSetType(), TSSundialsSetTolerance(),
845 @*/
TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)846 PetscErrorCode TSSundialsSetMaxTimeStep(TS ts,PetscReal maxdt)
847 {
848 PetscErrorCode ierr;
849
850 PetscFunctionBegin;
851 ierr = PetscTryMethod(ts,"TSSundialsSetMaxTimeStep_C",(TS,PetscReal),(ts,maxdt));CHKERRQ(ierr);
852 PetscFunctionReturn(0);
853 }
854
855 /*@
856 TSSundialsMonitorInternalSteps - Monitor Sundials internal steps (Defaults to false).
857
858 Input Parameter:
859 + ts - the time-step context
860 - ft - PETSC_TRUE if monitor, else PETSC_FALSE
861
862 Level: beginner
863
864 .seealso:TSSundialsGetIterations(), TSSundialsSetType(), TSSundialsSetMaxl(),
865 TSSundialsSetLinearTolerance(), TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(),
866 TSSundialsGetIterations(), TSSundialsSetType(),
867 TSSundialsSetLinearTolerance(), TSSundialsSetTolerance(), TSSundialsGetPC()
868 @*/
TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)869 PetscErrorCode TSSundialsMonitorInternalSteps(TS ts,PetscBool ft)
870 {
871 PetscErrorCode ierr;
872
873 PetscFunctionBegin;
874 ierr = PetscTryMethod(ts,"TSSundialsMonitorInternalSteps_C",(TS,PetscBool),(ts,ft));CHKERRQ(ierr);
875 PetscFunctionReturn(0);
876 }
877 /* -------------------------------------------------------------------------------------------*/
878 /*MC
879 TSSUNDIALS - ODE solver using the LLNL CVODE/SUNDIALS package (now called SUNDIALS)
880
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
889
890
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.
894
895 Level: beginner
896
897 .seealso: TSCreate(), TS, TSSetType(), TSSundialsSetType(), TSSundialsSetMaxl(), TSSundialsSetLinearTolerance(),
898 TSSundialsSetGramSchmidtType(), TSSundialsSetTolerance(), TSSundialsGetPC(), TSSundialsGetIterations(), TSSetExactFinalTime()
899
900 M*/
TSCreate_Sundials(TS ts)901 PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
902 {
903 TS_Sundials *cvode;
904 PetscErrorCode ierr;
905 PC pc;
906
907 PetscFunctionBegin;
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;
916
917 ierr = PetscNewLog(ts,&cvode);CHKERRQ(ierr);
918
919 ts->usessnes = PETSC_TRUE;
920
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;
928
929 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)ts),&(cvode->comm_sundials));CHKERRQ(ierr);
930
931 cvode->mindt = -1.;
932 cvode->maxdt = -1.;
933
934 /* set tolerance for Sundials */
935 cvode->reltol = 1e-6;
936 cvode->abstol = 1e-6;
937
938 /* set PCNONE as default pctype */
939 ierr = TSSundialsGetPC_Sundials(ts,&pc);CHKERRQ(ierr);
940 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
941
942 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetType_C",TSSundialsSetType_Sundials);CHKERRQ(ierr);
943 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxl_C",TSSundialsSetMaxl_Sundials);CHKERRQ(ierr);
944 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetLinearTolerance_C",TSSundialsSetLinearTolerance_Sundials);CHKERRQ(ierr);
945 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetGramSchmidtType_C",TSSundialsSetGramSchmidtType_Sundials);CHKERRQ(ierr);
946 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetTolerance_C",TSSundialsSetTolerance_Sundials);CHKERRQ(ierr);
947 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMinTimeStep_C",TSSundialsSetMinTimeStep_Sundials);CHKERRQ(ierr);
948 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsSetMaxTimeStep_C",TSSundialsSetMaxTimeStep_Sundials);CHKERRQ(ierr);
949 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetPC_C",TSSundialsGetPC_Sundials);CHKERRQ(ierr);
950 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsGetIterations_C",TSSundialsGetIterations_Sundials);CHKERRQ(ierr);
951 ierr = PetscObjectComposeFunction((PetscObject)ts,"TSSundialsMonitorInternalSteps_C",TSSundialsMonitorInternalSteps_Sundials);CHKERRQ(ierr);
952 PetscFunctionReturn(0);
953 }
954