Actual source code: sundials.c

  1: /*
  2:     Provides a PETSc interface to version 2.5 of SUNDIALS/CVODE solver (a very old version)
  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: static PetscErrorCode TSPrecond_Sundials_Petsc(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)
 15: {
 16:   TS           ts    = (TS)P_data;
 17:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
 18:   PC           pc;
 19:   Mat          J, P;
 20:   Vec          yy = cvode->w1, yydot = cvode->ydot;
 21:   PetscReal    gm = (PetscReal)_gamma;
 22:   PetscScalar *y_data;

 24:   PetscFunctionBegin;
 25:   PetscCall(TSGetIJacobian(ts, &J, &P, NULL, NULL));
 26:   y_data = (PetscScalar *)N_VGetArrayPointer(y);
 27:   PetscCall(VecPlaceArray(yy, y_data));
 28:   PetscCall(VecZeroEntries(yydot)); /* The Jacobian is independent of Ydot for ODE which is all that CVode works for */
 29:   /* compute the shifted Jacobian   (1/gm)*I + Jrest */
 30:   PetscCall(TSComputeIJacobian(ts, ts->ptime, yy, yydot, 1 / gm, J, P, PETSC_FALSE));
 31:   PetscCall(VecResetArray(yy));
 32:   PetscCall(MatScale(P, gm)); /* turn into I-gm*Jrest, J is not used by SUNDIALS  */
 33:   *jcurPtr = TRUE;
 34:   PetscCall(TSSundialsGetPC(ts, &pc));
 35:   PetscCall(PCSetOperators(pc, J, P));
 36:   PetscFunctionReturn(PETSC_SUCCESS);
 37: }

 39: /* Sundial expects an int (*)(args...) but PetscErrorCode is an enum. Instead of switching out
 40:    all the PetscCalls in TSPrecond_Sundials_Petsc we just wrap it */
 41: static int TSPrecond_Sundials_Private(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)
 42: {
 43:   return (int)TSPrecond_Sundials_Petsc(tn, y, fy, jok, jcurPtr, _gamma, P_data, vtemp1, vtemp2, vtemp3);
 44: }

 46: /*
 47:      TSPSolve_Sundials -  routine that we provide to SUNDIALS that applies the preconditioner.
 48: */
 49: static PetscErrorCode TSPSolve_Sundials_Petsc(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)
 50: {
 51:   TS           ts    = (TS)P_data;
 52:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
 53:   PC           pc;
 54:   Vec          rr = cvode->w1, zz = cvode->w2;
 55:   PetscScalar *r_data, *z_data;

 57:   PetscFunctionBegin;
 58:   /* Make the PETSc work vectors rr and zz point to the arrays in the SUNDIALS vectors r and z respectively*/
 59:   r_data = (PetscScalar *)N_VGetArrayPointer(r);
 60:   z_data = (PetscScalar *)N_VGetArrayPointer(z);
 61:   PetscCall(VecPlaceArray(rr, r_data));
 62:   PetscCall(VecPlaceArray(zz, z_data));

 64:   /* Solve the Px=r and put the result in zz */
 65:   PetscCall(TSSundialsGetPC(ts, &pc));
 66:   PetscCall(PCApply(pc, rr, zz));
 67:   PetscCall(VecResetArray(rr));
 68:   PetscCall(VecResetArray(zz));
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }

 72: /* See TSPrecond_Sundials_Private() */
 73: static int TSPSolve_Sundials_Private(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)
 74: {
 75:   return (int)TSPSolve_Sundials_Petsc(tn, y, fy, r, z, _gamma, delta, lr, P_data, vtemp);
 76: }

 78: /*
 79:         TSFunction_Sundials - routine that we provide to SUNDIALS that applies the right-hand side.
 80: */
 81: static int TSFunction_Sundials(realtype t, N_Vector y, N_Vector ydot, void *ctx)
 82: {
 83:   TS             ts = (TS)ctx;
 84:   DM             dm;
 85:   DMTS           tsdm;
 86:   TSIFunctionFn *ifunction;
 87:   MPI_Comm       comm;
 88:   TS_Sundials   *cvode = (TS_Sundials *)ts->data;
 89:   Vec            yy = cvode->w1, yyd = cvode->w2, yydot = cvode->ydot;
 90:   PetscScalar   *y_data, *ydot_data;

 92:   PetscFunctionBegin;
 93:   PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
 94:   /* Make the PETSc work vectors yy and yyd point to the arrays in the SUNDIALS vectors y and ydot respectively*/
 95:   y_data    = (PetscScalar *)N_VGetArrayPointer(y);
 96:   ydot_data = (PetscScalar *)N_VGetArrayPointer(ydot);
 97:   PetscCallAbort(comm, VecPlaceArray(yy, y_data));
 98:   PetscCallAbort(comm, VecPlaceArray(yyd, ydot_data));

100:   /* Now compute the right-hand side function, via IFunction unless only the more efficient RHSFunction is set */
101:   PetscCall(TSGetDM(ts, &dm));
102:   PetscCall(DMGetDMTS(dm, &tsdm));
103:   PetscCall(DMTSGetIFunction(dm, &ifunction, NULL));
104:   if (!ifunction) {
105:     PetscCall(TSComputeRHSFunction(ts, t, yy, yyd));
106:   } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */
107:     PetscCall(VecZeroEntries(yydot));
108:     PetscCallAbort(comm, TSComputeIFunction(ts, t, yy, yydot, yyd, PETSC_FALSE));
109:     PetscCall(VecScale(yyd, -1.));
110:   }
111:   PetscCallAbort(comm, VecResetArray(yy));
112:   PetscCallAbort(comm, VecResetArray(yyd));
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

116: /*
117:        TSStep_Sundials - Calls SUNDIALS to integrate the ODE.
118: */
119: static PetscErrorCode TSStep_Sundials(TS ts)
120: {
121:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
122:   PetscInt     flag;
123:   long int     nits, lits, nsteps;
124:   realtype     t, tout;
125:   PetscScalar *y_data;
126:   void        *mem;

128:   PetscFunctionBegin;
129:   mem  = cvode->mem;
130:   tout = ts->max_time;
131:   PetscCall(VecGetArray(ts->vec_sol, &y_data));
132:   N_VSetArrayPointer((realtype *)y_data, cvode->y);
133:   PetscCall(VecRestoreArray(ts->vec_sol, NULL));

135:   /* We would like to TSPreStage() and TSPostStage()
136:    * before each stage solve but CVode does not appear to support this. */
137:   if (cvode->monitorstep) flag = CVode(mem, tout, cvode->y, &t, CV_ONE_STEP);
138:   else flag = CVode(mem, tout, cvode->y, &t, CV_NORMAL);

140:   if (flag) { /* display error message */
141:     switch (flag) {
142:     case CV_ILL_INPUT:
143:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_ILL_INPUT");
144:       break;
145:     case CV_TOO_CLOSE:
146:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_TOO_CLOSE");
147:       break;
148:     case CV_TOO_MUCH_WORK: {
149:       PetscReal tcur;
150:       PetscCallExternal(CVodeGetNumSteps, mem, &nsteps);
151:       PetscCallExternal(CVodeGetCurrentTime, mem, &tcur);
152:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_TOO_MUCH_WORK. At t=%g, nsteps %ld exceeds maxstep %" PetscInt_FMT ". Increase '-ts_max_steps <>' or modify TSSetMaxSteps()", (double)tcur, nsteps, ts->max_steps);
153:     } break;
154:     case CV_TOO_MUCH_ACC:
155:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_TOO_MUCH_ACC");
156:       break;
157:     case CV_ERR_FAILURE:
158:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_ERR_FAILURE");
159:       break;
160:     case CV_CONV_FAILURE:
161:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_CONV_FAILURE");
162:       break;
163:     case CV_LINIT_FAIL:
164:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_LINIT_FAIL");
165:       break;
166:     case CV_LSETUP_FAIL:
167:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_LSETUP_FAIL");
168:       break;
169:     case CV_LSOLVE_FAIL:
170:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_LSOLVE_FAIL");
171:       break;
172:     case CV_RHSFUNC_FAIL:
173:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_RHSFUNC_FAIL");
174:       break;
175:     case CV_FIRST_RHSFUNC_ERR:
176:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_FIRST_RHSFUNC_ERR");
177:       break;
178:     case CV_REPTD_RHSFUNC_ERR:
179:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_REPTD_RHSFUNC_ERR");
180:       break;
181:     case CV_UNREC_RHSFUNC_ERR:
182:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_UNREC_RHSFUNC_ERR");
183:       break;
184:     case CV_RTFUNC_FAIL:
185:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, CV_RTFUNC_FAIL");
186:       break;
187:     default:
188:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "CVode() fails, flag %d", flag);
189:     }
190:   }

192:   /* log inner nonlinear and linear iterations */
193:   PetscCallExternal(CVodeGetNumNonlinSolvIters, mem, &nits);
194:   PetscCallExternal(CVSpilsGetNumLinIters, mem, &lits);
195:   ts->snes_its += nits;
196:   ts->ksp_its = lits;

198:   /* copy the solution from cvode->y to cvode->update and sol */
199:   PetscCall(VecPlaceArray(cvode->w1, y_data));
200:   PetscCall(VecCopy(cvode->w1, cvode->update));
201:   PetscCall(VecResetArray(cvode->w1));
202:   PetscCall(VecCopy(cvode->update, ts->vec_sol));

204:   ts->time_step = t - ts->ptime;
205:   ts->ptime     = t;

207:   PetscCallExternal(CVodeGetNumSteps, mem, &nsteps);
208:   if (!cvode->monitorstep) ts->steps += nsteps - 1; /* TSStep() increments the step counter by one */
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: static PetscErrorCode TSInterpolate_Sundials(TS ts, PetscReal t, Vec X)
213: {
214:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
215:   N_Vector     y;
216:   PetscScalar *x_data;
217:   PetscInt     glosize, locsize;

219:   PetscFunctionBegin;
220:   /* get the vector size */
221:   PetscCall(VecGetSize(X, &glosize));
222:   PetscCall(VecGetLocalSize(X, &locsize));
223:   PetscCall(VecGetArray(X, &x_data));

225:   /* Initialize N_Vec y with x_data */
226:   if (cvode->use_dense) {
227:     PetscMPIInt size;

229:     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
230:     PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "TSSUNDIALS only supports a dense solve in the serial case");
231:     y = N_VMake_Serial(locsize, (realtype *)x_data);
232:   } else {
233:     y = N_VMake_Parallel(cvode->comm_sundials, locsize, glosize, (realtype *)x_data);
234:   }

236:   PetscCheck(y, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Interpolated y is not allocated");

238:   PetscCallExternal(CVodeGetDky, cvode->mem, t, 0, y);
239:   PetscCall(VecRestoreArray(X, &x_data));
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: static PetscErrorCode TSReset_Sundials(TS ts)
244: {
245:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

247:   PetscFunctionBegin;
248:   PetscCall(VecDestroy(&cvode->update));
249:   PetscCall(VecDestroy(&cvode->ydot));
250:   PetscCall(VecDestroy(&cvode->w1));
251:   PetscCall(VecDestroy(&cvode->w2));
252:   if (cvode->mem) CVodeFree(&cvode->mem);
253:   PetscFunctionReturn(PETSC_SUCCESS);
254: }

256: static PetscErrorCode TSDestroy_Sundials(TS ts)
257: {
258:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

260:   PetscFunctionBegin;
261:   PetscCall(TSReset_Sundials(ts));
262:   PetscCallMPI(MPI_Comm_free(&cvode->comm_sundials));
263:   PetscCall(PetscFree(ts->data));
264:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetType_C", NULL));
265:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxl_C", NULL));
266:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetLinearTolerance_C", NULL));
267:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetGramSchmidtType_C", NULL));
268:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetTolerance_C", NULL));
269:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMinTimeStep_C", NULL));
270:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxTimeStep_C", NULL));
271:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetPC_C", NULL));
272:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetIterations_C", NULL));
273:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsMonitorInternalSteps_C", NULL));
274:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetUseDense_C", NULL));
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: static PetscErrorCode TSSetUp_Sundials(TS ts)
279: {
280:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
281:   PetscInt     glosize, locsize, i;
282:   PetscScalar *y_data, *parray;
283:   PC           pc;
284:   PCType       pctype;
285:   PetscBool    pcnone;

287:   PetscFunctionBegin;
288:   PetscCheck(ts->exact_final_time != TS_EXACTFINALTIME_MATCHSTEP, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for exact final time option 'MATCHSTEP' when using SUNDIALS");

290:   /* get the vector size */
291:   PetscCall(VecGetSize(ts->vec_sol, &glosize));
292:   PetscCall(VecGetLocalSize(ts->vec_sol, &locsize));

294:   /* allocate the memory for N_Vec y */
295:   if (cvode->use_dense) {
296:     PetscMPIInt size;

298:     PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
299:     PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "TSSUNDIALS only supports a dense solve in the serial case");
300:     cvode->y = N_VNew_Serial(locsize);
301:   } else {
302:     cvode->y = N_VNew_Parallel(cvode->comm_sundials, locsize, glosize);
303:   }
304:   PetscCheck(cvode->y, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "cvode->y is not allocated");

306:   /* initialize N_Vec y: copy ts->vec_sol to cvode->y */
307:   PetscCall(VecGetArray(ts->vec_sol, &parray));
308:   y_data = (PetscScalar *)N_VGetArrayPointer(cvode->y);
309:   for (i = 0; i < locsize; i++) y_data[i] = parray[i];
310:   PetscCall(VecRestoreArray(ts->vec_sol, NULL));

312:   PetscCall(VecDuplicate(ts->vec_sol, &cvode->update));
313:   PetscCall(VecDuplicate(ts->vec_sol, &cvode->ydot));

315:   /*
316:     Create work vectors for the TSPSolve_Sundials() routine. Note these are
317:     allocated with zero space arrays because the actual array space is provided
318:     by SUNDIALS and set using VecPlaceArray().
319:   */
320:   PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts), 1, locsize, PETSC_DECIDE, NULL, &cvode->w1));
321:   PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)ts), 1, locsize, PETSC_DECIDE, NULL, &cvode->w2));

323:   /* Call CVodeCreate to create the solver memory and the use of a Newton iteration */
324:   cvode->mem = CVodeCreate(cvode->cvode_type, CV_NEWTON);
325:   PetscCheck(cvode->mem, PETSC_COMM_SELF, PETSC_ERR_MEM, "CVodeCreate() fails");

327:   /* Set the pointer to user-defined data */
328:   PetscCallExternal(CVodeSetUserData, cvode->mem, ts);

330:   /* SUNDIALS may choose to use a smaller initial step, but will never use a larger step. */
331:   PetscCallExternal(CVodeSetInitStep, cvode->mem, (realtype)ts->time_step);
332:   if (cvode->mindt > 0) {
333:     int flag = CVodeSetMinStep(cvode->mem, (realtype)cvode->mindt);
334:     if (flag) {
335:       PetscCheck(flag != CV_MEM_NULL, PetscObjectComm((PetscObject)ts), PETSC_ERR_LIB, "CVodeSetMinStep() failed, cvode_mem pointer is NULL");
336:       PetscCheck(flag != CV_ILL_INPUT, PetscObjectComm((PetscObject)ts), PETSC_ERR_LIB, "CVodeSetMinStep() failed, hmin is nonpositive or it exceeds the maximum allowable step size");
337:       SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_LIB, "CVodeSetMinStep() failed");
338:     }
339:   }
340:   if (cvode->maxdt > 0) PetscCallExternal(CVodeSetMaxStep, cvode->mem, (realtype)cvode->maxdt);

342:   /* Call CVodeInit to initialize the integrator memory and specify the
343:    * user's right-hand side function in u'=f(t,u), the initial time T0, and
344:    * the initial dependent variable vector cvode->y */
345:   PetscCallExternal(CVodeInit, cvode->mem, TSFunction_Sundials, ts->ptime, cvode->y);

347:   /* specifies scalar relative and absolute tolerances */
348:   PetscCallExternal(CVodeSStolerances, cvode->mem, cvode->reltol, cvode->abstol);

350:   /* Specify max order of BDF / ADAMS method */
351:   if (cvode->maxord != PETSC_DEFAULT) PetscCallExternal(CVodeSetMaxOrd, cvode->mem, cvode->maxord);

353:   /* Specify max num of steps to be taken by cvode in its attempt to reach the next output time */
354:   PetscCallExternal(CVodeSetMaxNumSteps, cvode->mem, ts->max_steps);

356:   if (cvode->use_dense) {
357:     PetscCallExternal(CVDense, cvode->mem, locsize);
358:   } else {
359:     /* call CVSpgmr to use GMRES as the linear solver.        */
360:     /* setup the ode integrator with the given preconditioner */
361:     PetscCall(TSSundialsGetPC(ts, &pc));
362:     PetscCall(PCGetType(pc, &pctype));
363:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCNONE, &pcnone));
364:     if (pcnone) {
365:       PetscCallExternal(CVSpgmr, cvode->mem, PREC_NONE, 0);
366:     } else {
367:       PetscCallExternal(CVSpgmr, cvode->mem, PREC_LEFT, cvode->maxl);

369:       /* Set preconditioner and solve routines Precond and PSolve,
370:          and the pointer to the user-defined block data */
371:       PetscCallExternal(CVSpilsSetPreconditioner, cvode->mem, TSPrecond_Sundials_Private, TSPSolve_Sundials_Private);
372:     }
373:   }
374:   PetscFunctionReturn(PETSC_SUCCESS);
375: }

377: /* type of CVODE linear multistep method */
378: const char *const TSSundialsLmmTypes[] = {"", "ADAMS", "BDF", "TSSundialsLmmType", "SUNDIALS_", NULL};
379: /* type of G-S orthogonalization used by CVODE linear solver */
380: const char *const TSSundialsGramSchmidtTypes[] = {"", "MODIFIED", "CLASSICAL", "TSSundialsGramSchmidtType", "SUNDIALS_", NULL};

382: static PetscErrorCode TSSetFromOptions_Sundials(TS ts, PetscOptionItems *PetscOptionsObject)
383: {
384:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
385:   int          indx;
386:   PetscBool    flag;
387:   PC           pc;

389:   PetscFunctionBegin;
390:   PetscOptionsHeadBegin(PetscOptionsObject, "SUNDIALS ODE solver options");
391:   PetscCall(PetscOptionsEList("-ts_sundials_type", "Scheme", "TSSundialsSetType", TSSundialsLmmTypes, 3, TSSundialsLmmTypes[cvode->cvode_type], &indx, &flag));
392:   if (flag) PetscCall(TSSundialsSetType(ts, (TSSundialsLmmType)indx));
393:   PetscCall(PetscOptionsEList("-ts_sundials_gramschmidt_type", "Type of orthogonalization", "TSSundialsSetGramSchmidtType", TSSundialsGramSchmidtTypes, 3, TSSundialsGramSchmidtTypes[cvode->gtype], &indx, &flag));
394:   if (flag) PetscCall(TSSundialsSetGramSchmidtType(ts, (TSSundialsGramSchmidtType)indx));
395:   PetscCall(PetscOptionsReal("-ts_sundials_atol", "Absolute tolerance for convergence", "TSSundialsSetTolerance", cvode->abstol, &cvode->abstol, NULL));
396:   PetscCall(PetscOptionsReal("-ts_sundials_rtol", "Relative tolerance for convergence", "TSSundialsSetTolerance", cvode->reltol, &cvode->reltol, NULL));
397:   PetscCall(PetscOptionsReal("-ts_sundials_mindt", "Minimum step size", "TSSundialsSetMinTimeStep", cvode->mindt, &cvode->mindt, NULL));
398:   PetscCall(PetscOptionsReal("-ts_sundials_maxdt", "Maximum step size", "TSSundialsSetMaxTimeStep", cvode->maxdt, &cvode->maxdt, NULL));
399:   PetscCall(PetscOptionsReal("-ts_sundials_linear_tolerance", "Convergence tolerance for linear solve", "TSSundialsSetLinearTolerance", cvode->linear_tol, &cvode->linear_tol, NULL));
400:   PetscCall(PetscOptionsInt("-ts_sundials_maxord", "Max Order for BDF/Adams method", "TSSundialsSetMaxOrd", cvode->maxord, &cvode->maxord, NULL));
401:   PetscCall(PetscOptionsInt("-ts_sundials_maxl", "Max dimension of the Krylov subspace", "TSSundialsSetMaxl", cvode->maxl, &cvode->maxl, NULL));
402:   PetscCall(PetscOptionsBool("-ts_sundials_monitor_steps", "Monitor SUNDIALS internal steps", "TSSundialsMonitorInternalSteps", cvode->monitorstep, &cvode->monitorstep, NULL));
403:   PetscCall(PetscOptionsBool("-ts_sundials_use_dense", "Use dense internal solver in SUNDIALS (serial only)", "TSSundialsSetUseDense", cvode->use_dense, &cvode->use_dense, NULL));
404:   PetscOptionsHeadEnd();
405:   PetscCall(TSSundialsGetPC(ts, &pc));
406:   PetscCall(PCSetFromOptions(pc));
407:   PetscFunctionReturn(PETSC_SUCCESS);
408: }

410: static PetscErrorCode TSView_Sundials(TS ts, PetscViewer viewer)
411: {
412:   TS_Sundials *cvode = (TS_Sundials *)ts->data;
413:   char        *type;
414:   char         atype[] = "Adams";
415:   char         btype[] = "BDF: backward differentiation formula";
416:   PetscBool    iascii, isstring;
417:   long int     nsteps, its, nfevals, nlinsetups, nfails, itmp;
418:   PetscInt     qlast, qcur;
419:   PetscReal    hinused, hlast, hcur, tcur, tolsfac;
420:   PC           pc;

422:   PetscFunctionBegin;
423:   if (cvode->cvode_type == SUNDIALS_ADAMS) type = atype;
424:   else type = btype;

426:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
427:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
428:   if (iascii) {
429:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS integrator does not use SNES!\n"));
430:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS integrator type %s\n", type));
431:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS maxord %" PetscInt_FMT "\n", cvode->maxord));
432:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS abs tol %g rel tol %g\n", (double)cvode->abstol, (double)cvode->reltol));
433:     if (cvode->use_dense) {
434:       PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS integrator using a dense linear solve\n"));
435:     } else {
436:       PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS linear solver tolerance factor %g\n", (double)cvode->linear_tol));
437:       PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS max dimension of Krylov subspace %" PetscInt_FMT "\n", cvode->maxl));
438:       if (cvode->gtype == SUNDIALS_MODIFIED_GS) {
439:         PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS using modified Gram-Schmidt for orthogonalization in GMRES\n"));
440:       } else {
441:         PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS using unmodified (classical) Gram-Schmidt for orthogonalization in GMRES\n"));
442:       }
443:     }
444:     if (cvode->mindt > 0) PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS minimum time step %g\n", (double)cvode->mindt));
445:     if (cvode->maxdt > 0) PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS maximum time step %g\n", (double)cvode->maxdt));

447:     /* Outputs from CVODE, CVSPILS */
448:     PetscCallExternal(CVodeGetTolScaleFactor, cvode->mem, &tolsfac);
449:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS suggested factor for tolerance scaling %g\n", tolsfac));
450:     PetscCallExternal(CVodeGetIntegratorStats, cvode->mem, &nsteps, &nfevals, &nlinsetups, &nfails, &qlast, &qcur, &hinused, &hlast, &hcur, &tcur);
451:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS cumulative number of internal steps %ld\n", nsteps));
452:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of calls to rhs function %ld\n", nfevals));
453:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of calls to linear solver setup function %ld\n", nlinsetups));
454:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of error test failures %ld\n", nfails));

456:     PetscCallExternal(CVodeGetNonlinSolvStats, cvode->mem, &its, &nfails);
457:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of nonlinear solver iterations %ld\n", its));
458:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of nonlinear convergence failure %ld\n", nfails));
459:     if (!cvode->use_dense) {
460:       PetscCallExternal(CVSpilsGetNumLinIters, cvode->mem, &its); /* its = no. of calls to TSPrecond_Sundials() */
461:       PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of linear iterations %ld\n", its));
462:       PetscCallExternal(CVSpilsGetNumConvFails, cvode->mem, &itmp);
463:       PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of linear convergence failures %ld\n", itmp));

465:       PetscCall(TSSundialsGetPC(ts, &pc));
466:       PetscCall(PCView(pc, viewer));
467:       PetscCallExternal(CVSpilsGetNumPrecEvals, cvode->mem, &itmp);
468:       PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of preconditioner evaluations %ld\n", itmp));
469:       PetscCallExternal(CVSpilsGetNumPrecSolves, cvode->mem, &itmp);
470:       PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of preconditioner solves %ld\n", itmp));
471:     }
472:     PetscCallExternal(CVSpilsGetNumJtimesEvals, cvode->mem, &itmp);
473:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of Jacobian-vector product evaluations %ld\n", itmp));
474:     PetscCallExternal(CVSpilsGetNumRhsEvals, cvode->mem, &itmp);
475:     PetscCall(PetscViewerASCIIPrintf(viewer, "SUNDIALS no. of rhs calls for finite diff. Jacobian-vector evals %ld\n", itmp));
476:   } else if (isstring) {
477:     PetscCall(PetscViewerStringSPrintf(viewer, "SUNDIALS type %s", type));
478:   }
479:   PetscFunctionReturn(PETSC_SUCCESS);
480: }

482: /* --------------------------------------------------------------------------*/
483: static PetscErrorCode TSSundialsSetType_Sundials(TS ts, TSSundialsLmmType type)
484: {
485:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

487:   PetscFunctionBegin;
488:   cvode->cvode_type = type;
489:   PetscFunctionReturn(PETSC_SUCCESS);
490: }

492: static PetscErrorCode TSSundialsSetMaxl_Sundials(TS ts, PetscInt maxl)
493: {
494:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

496:   PetscFunctionBegin;
497:   cvode->maxl = maxl;
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: }

501: static PetscErrorCode TSSundialsSetLinearTolerance_Sundials(TS ts, PetscReal tol)
502: {
503:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

505:   PetscFunctionBegin;
506:   cvode->linear_tol = tol;
507:   PetscFunctionReturn(PETSC_SUCCESS);
508: }

510: static PetscErrorCode TSSundialsSetGramSchmidtType_Sundials(TS ts, TSSundialsGramSchmidtType type)
511: {
512:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

514:   PetscFunctionBegin;
515:   cvode->gtype = type;
516:   PetscFunctionReturn(PETSC_SUCCESS);
517: }

519: static PetscErrorCode TSSundialsSetTolerance_Sundials(TS ts, PetscReal aabs, PetscReal rel)
520: {
521:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

523:   PetscFunctionBegin;
524:   if (aabs != PETSC_DECIDE) cvode->abstol = aabs;
525:   if (rel != PETSC_DECIDE) cvode->reltol = rel;
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }

529: static PetscErrorCode TSSundialsSetMinTimeStep_Sundials(TS ts, PetscReal mindt)
530: {
531:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

533:   PetscFunctionBegin;
534:   cvode->mindt = mindt;
535:   PetscFunctionReturn(PETSC_SUCCESS);
536: }

538: static PetscErrorCode TSSundialsSetMaxTimeStep_Sundials(TS ts, PetscReal maxdt)
539: {
540:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

542:   PetscFunctionBegin;
543:   cvode->maxdt = maxdt;
544:   PetscFunctionReturn(PETSC_SUCCESS);
545: }

547: static PetscErrorCode TSSundialsSetUseDense_Sundials(TS ts, PetscBool use_dense)
548: {
549:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

551:   PetscFunctionBegin;
552:   cvode->use_dense = use_dense;
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }

556: static PetscErrorCode TSSundialsGetPC_Sundials(TS ts, PC *pc)
557: {
558:   SNES snes;
559:   KSP  ksp;

561:   PetscFunctionBegin;
562:   PetscCall(TSGetSNES(ts, &snes));
563:   PetscCall(SNESGetKSP(snes, &ksp));
564:   PetscCall(KSPGetPC(ksp, pc));
565:   PetscFunctionReturn(PETSC_SUCCESS);
566: }

568: static PetscErrorCode TSSundialsGetIterations_Sundials(TS ts, int *nonlin, int *lin)
569: {
570:   PetscFunctionBegin;
571:   if (nonlin) *nonlin = ts->snes_its;
572:   if (lin) *lin = ts->ksp_its;
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: static PetscErrorCode TSSundialsMonitorInternalSteps_Sundials(TS ts, PetscBool s)
577: {
578:   TS_Sundials *cvode = (TS_Sundials *)ts->data;

580:   PetscFunctionBegin;
581:   cvode->monitorstep = s;
582:   PetscFunctionReturn(PETSC_SUCCESS);
583: }
584: /* -------------------------------------------------------------------------------------------*/

586: /*@C
587:   TSSundialsGetIterations - Gets the number of nonlinear and linear iterations used so far by `TSSUNDIALS`.

589:   Not Collective

591:   Input Parameter:
592: . ts - the time-step context

594:   Output Parameters:
595: + nonlin - number of nonlinear iterations
596: - lin    - number of linear iterations

598:   Level: advanced

600:   Note:
601:   These return the number since the creation of the `TS` object

603: .seealso: [](ch_ts), `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
604:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
605:           `TSSundialsGetPC()`, `TSSetExactFinalTime()`
606: @*/
607: PetscErrorCode TSSundialsGetIterations(TS ts, int *nonlin, int *lin)
608: {
609:   PetscFunctionBegin;
610:   PetscUseMethod(ts, "TSSundialsGetIterations_C", (TS, int *, int *), (ts, nonlin, lin));
611:   PetscFunctionReturn(PETSC_SUCCESS);
612: }

614: /*@
615:   TSSundialsSetType - Sets the method that `TSSUNDIALS` will use for integration.

617:   Logically Collective

619:   Input Parameters:
620: + ts   - the time-step context
621: - type - one of  `SUNDIALS_ADAMS` or `SUNDIALS_BDF`

623:   Level: intermediate

625: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetMaxl()`,
626:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`,
627:           `TSSundialsSetTolerance()`, `TSSundialsGetPC()`, `TSSetExactFinalTime()`
628: @*/
629: PetscErrorCode TSSundialsSetType(TS ts, TSSundialsLmmType type)
630: {
631:   PetscFunctionBegin;
632:   PetscTryMethod(ts, "TSSundialsSetType_C", (TS, TSSundialsLmmType), (ts, type));
633:   PetscFunctionReturn(PETSC_SUCCESS);
634: }

636: /*@
637:   TSSundialsSetMaxord - Sets the maximum order for BDF/Adams method used by `TSSUNDIALS`.

639:   Logically Collective

641:   Input Parameters:
642: + ts     - the time-step context
643: - maxord - maximum order of BDF / Adams method

645:   Level: advanced

647: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`,
648:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`,
649:           `TSSundialsSetTolerance()`, `TSSundialsGetPC()`, `TSSetExactFinalTime()`
650: @*/
651: PetscErrorCode TSSundialsSetMaxord(TS ts, PetscInt maxord)
652: {
653:   PetscFunctionBegin;
655:   PetscTryMethod(ts, "TSSundialsSetMaxOrd_C", (TS, PetscInt), (ts, maxord));
656:   PetscFunctionReturn(PETSC_SUCCESS);
657: }

659: /*@
660:   TSSundialsSetMaxl - Sets the dimension of the Krylov space used by
661:   GMRES in the linear solver in `TSSUNDIALS`. `TSSUNDIALS` DOES NOT use restarted GMRES so
662:   this is the maximum number of GMRES steps that will be used.

664:   Logically Collective

666:   Input Parameters:
667: + ts   - the time-step context
668: - maxl - number of direction vectors (the dimension of Krylov subspace).

670:   Level: advanced

672: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`,
673:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
674:           `TSSundialsGetPC()`, `TSSetExactFinalTime()`
675: @*/
676: PetscErrorCode TSSundialsSetMaxl(TS ts, PetscInt maxl)
677: {
678:   PetscFunctionBegin;
680:   PetscTryMethod(ts, "TSSundialsSetMaxl_C", (TS, PetscInt), (ts, maxl));
681:   PetscFunctionReturn(PETSC_SUCCESS);
682: }

684: /*@
685:   TSSundialsSetLinearTolerance - Sets the tolerance used to solve the linear
686:   system by `TSSUNDIALS`.

688:   Logically Collective

690:   Input Parameters:
691: + ts  - the time-step context
692: - tol - the factor by which the tolerance on the nonlinear solver is
693:              multiplied to get the tolerance on the linear solver, .05 by default.

695:   Level: advanced

697: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
698:           `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
699:           `TSSundialsGetPC()`,
700:           `TSSetExactFinalTime()`
701: @*/
702: PetscErrorCode TSSundialsSetLinearTolerance(TS ts, PetscReal tol)
703: {
704:   PetscFunctionBegin;
706:   PetscTryMethod(ts, "TSSundialsSetLinearTolerance_C", (TS, PetscReal), (ts, tol));
707:   PetscFunctionReturn(PETSC_SUCCESS);
708: }

710: /*@
711:   TSSundialsSetGramSchmidtType - Sets type of orthogonalization used
712:   in GMRES method by `TSSUNDIALS` linear solver.

714:   Logically Collective

716:   Input Parameters:
717: + ts   - the time-step context
718: - type - either `SUNDIALS_MODIFIED_GS` or `SUNDIALS_CLASSICAL_GS`

720:   Level: advanced

722: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
723:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetTolerance()`,
724:           `TSSundialsGetPC()`,
725:           `TSSetExactFinalTime()`
726: @*/
727: PetscErrorCode TSSundialsSetGramSchmidtType(TS ts, TSSundialsGramSchmidtType type)
728: {
729:   PetscFunctionBegin;
730:   PetscTryMethod(ts, "TSSundialsSetGramSchmidtType_C", (TS, TSSundialsGramSchmidtType), (ts, type));
731:   PetscFunctionReturn(PETSC_SUCCESS);
732: }

734: /*@
735:   TSSundialsSetTolerance - Sets the absolute and relative tolerance used by
736:   `TSSUNDIALS` for error control.

738:   Logically Collective

740:   Input Parameters:
741: + ts   - the time-step context
742: . aabs - the absolute tolerance
743: - rel  - the relative tolerance

745:      See the CVODE/SUNDIALS users manual for exact details on these parameters. Essentially
746:     these regulate the size of the error for a SINGLE timestep.

748:   Level: intermediate

750: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetGMRESMaxl()`,
751:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`,
752:           `TSSundialsGetPC()`,
753:           `TSSetExactFinalTime()`
754: @*/
755: PetscErrorCode TSSundialsSetTolerance(TS ts, PetscReal aabs, PetscReal rel)
756: {
757:   PetscFunctionBegin;
758:   PetscTryMethod(ts, "TSSundialsSetTolerance_C", (TS, PetscReal, PetscReal), (ts, aabs, rel));
759:   PetscFunctionReturn(PETSC_SUCCESS);
760: }

762: /*@
763:   TSSundialsGetPC - Extract the PC context from a time-step context for `TSSUNDIALS`.

765:   Input Parameter:
766: . ts - the time-step context

768:   Output Parameter:
769: . pc - the preconditioner context

771:   Level: advanced

773: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
774:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`
775: @*/
776: PetscErrorCode TSSundialsGetPC(TS ts, PC *pc)
777: {
778:   PetscFunctionBegin;
779:   PetscUseMethod(ts, "TSSundialsGetPC_C", (TS, PC *), (ts, pc));
780:   PetscFunctionReturn(PETSC_SUCCESS);
781: }

783: /*@
784:   TSSundialsSetMinTimeStep - Smallest time step to be chosen by the adaptive controller.

786:   Input Parameters:
787: + ts    - the time-step context
788: - mindt - lowest time step if positive, negative to deactivate

790:   Note:
791:   `TSSUNDIALS` will error if it is not possible to keep the estimated truncation error below
792:   the tolerance set with `TSSundialsSetTolerance()` without going below this step size.

794:   Level: beginner

796: .seealso: [](ch_ts), `TSSundialsSetType()`, `TSSundialsSetTolerance()`,
797: @*/
798: PetscErrorCode TSSundialsSetMinTimeStep(TS ts, PetscReal mindt)
799: {
800:   PetscFunctionBegin;
801:   PetscTryMethod(ts, "TSSundialsSetMinTimeStep_C", (TS, PetscReal), (ts, mindt));
802:   PetscFunctionReturn(PETSC_SUCCESS);
803: }

805: /*@
806:   TSSundialsSetMaxTimeStep - Largest time step to be chosen by the adaptive controller.

808:   Input Parameters:
809: + ts    - the time-step context
810: - maxdt - lowest time step if positive, negative to deactivate

812:   Level: beginner

814: .seealso: [](ch_ts), `TSSundialsSetType()`, `TSSundialsSetTolerance()`,
815: @*/
816: PetscErrorCode TSSundialsSetMaxTimeStep(TS ts, PetscReal maxdt)
817: {
818:   PetscFunctionBegin;
819:   PetscTryMethod(ts, "TSSundialsSetMaxTimeStep_C", (TS, PetscReal), (ts, maxdt));
820:   PetscFunctionReturn(PETSC_SUCCESS);
821: }

823: /*@
824:   TSSundialsMonitorInternalSteps - Monitor `TSSUNDIALS` internal steps (Defaults to false).

826:   Input Parameters:
827: + ts - the time-step context
828: - ft - `PETSC_TRUE` if monitor, else `PETSC_FALSE`

830:   Level: beginner

832: .seealso: [](ch_ts), `TSSundialsGetIterations()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`,
833:           `TSSundialsSetLinearTolerance()`, `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`,
834:           `TSSundialsGetPC()`
835: @*/
836: PetscErrorCode TSSundialsMonitorInternalSteps(TS ts, PetscBool ft)
837: {
838:   PetscFunctionBegin;
839:   PetscTryMethod(ts, "TSSundialsMonitorInternalSteps_C", (TS, PetscBool), (ts, ft));
840:   PetscFunctionReturn(PETSC_SUCCESS);
841: }

843: /*@
844:   TSSundialsSetUseDense - Set a flag to use a dense linear solver in `TSSUNDIALS` (serial only)

846:   Logically Collective

848:   Input Parameters:
849: + ts        - the time-step context
850: - use_dense - `PETSC_TRUE` to use the dense solver

852:   Level: advanced

854: .seealso: [](ch_ts), `TSSUNDIALS`
855: @*/
856: PetscErrorCode TSSundialsSetUseDense(TS ts, PetscBool use_dense)
857: {
858:   PetscFunctionBegin;
860:   PetscTryMethod(ts, "TSSundialsSetUseDense_C", (TS, PetscBool), (ts, use_dense));
861:   PetscFunctionReturn(PETSC_SUCCESS);
862: }

864: /* -------------------------------------------------------------------------------------------*/
865: /*MC
866:       TSSUNDIALS - ODE solver using a very old version of the LLNL CVODE/SUNDIALS package, version 2.5 (now called SUNDIALS). Requires ./configure --download-sundials

868:    Options Database Keys:
869: +    -ts_sundials_type <bdf,adams> -
870: .    -ts_sundials_gramschmidt_type <modified, classical> - type of orthogonalization inside GMRES
871: .    -ts_sundials_atol <tol> - Absolute tolerance for convergence
872: .    -ts_sundials_rtol <tol> - Relative tolerance for convergence
873: .    -ts_sundials_linear_tolerance <tol> -
874: .    -ts_sundials_maxl <maxl> - Max dimension of the Krylov subspace
875: .    -ts_sundials_monitor_steps - Monitor SUNDIALS internal steps
876: -    -ts_sundials_use_dense - Use a dense linear solver within CVODE (serial only)

878:     Level: beginner

880:     Note:
881:     This uses its own nonlinear solver and Krylov method so PETSc `SNES` and `KSP` options do not apply,
882:     only PETSc `PC` options.

884: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSSundialsSetType()`, `TSSundialsSetMaxl()`, `TSSundialsSetLinearTolerance()`, `TSType`,
885:           `TSSundialsSetGramSchmidtType()`, `TSSundialsSetTolerance()`, `TSSundialsGetPC()`, `TSSundialsGetIterations()`, `TSSetExactFinalTime()`
886: M*/
887: PETSC_EXTERN PetscErrorCode TSCreate_Sundials(TS ts)
888: {
889:   TS_Sundials *cvode;
890:   PC           pc;

892:   PetscFunctionBegin;
893:   ts->ops->reset          = TSReset_Sundials;
894:   ts->ops->destroy        = TSDestroy_Sundials;
895:   ts->ops->view           = TSView_Sundials;
896:   ts->ops->setup          = TSSetUp_Sundials;
897:   ts->ops->step           = TSStep_Sundials;
898:   ts->ops->interpolate    = TSInterpolate_Sundials;
899:   ts->ops->setfromoptions = TSSetFromOptions_Sundials;
900:   ts->default_adapt_type  = TSADAPTNONE;

902:   PetscCall(PetscNew(&cvode));

904:   ts->usessnes = PETSC_TRUE;

906:   ts->data           = (void *)cvode;
907:   cvode->cvode_type  = SUNDIALS_BDF;
908:   cvode->gtype       = SUNDIALS_CLASSICAL_GS;
909:   cvode->maxl        = 5;
910:   cvode->maxord      = PETSC_DEFAULT;
911:   cvode->linear_tol  = .05;
912:   cvode->monitorstep = PETSC_TRUE;
913:   cvode->use_dense   = PETSC_FALSE;

915:   PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)ts), &cvode->comm_sundials));

917:   cvode->mindt = -1.;
918:   cvode->maxdt = -1.;

920:   /* set tolerance for SUNDIALS */
921:   cvode->reltol = 1e-6;
922:   cvode->abstol = 1e-6;

924:   /* set PCNONE as default pctype */
925:   PetscCall(TSSundialsGetPC_Sundials(ts, &pc));
926:   PetscCall(PCSetType(pc, PCNONE));

928:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetType_C", TSSundialsSetType_Sundials));
929:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxl_C", TSSundialsSetMaxl_Sundials));
930:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetLinearTolerance_C", TSSundialsSetLinearTolerance_Sundials));
931:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetGramSchmidtType_C", TSSundialsSetGramSchmidtType_Sundials));
932:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetTolerance_C", TSSundialsSetTolerance_Sundials));
933:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMinTimeStep_C", TSSundialsSetMinTimeStep_Sundials));
934:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetMaxTimeStep_C", TSSundialsSetMaxTimeStep_Sundials));
935:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetPC_C", TSSundialsGetPC_Sundials));
936:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsGetIterations_C", TSSundialsGetIterations_Sundials));
937:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsMonitorInternalSteps_C", TSSundialsMonitorInternalSteps_Sundials));
938:   PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSSundialsSetUseDense_C", TSSundialsSetUseDense_Sundials));
939:   PetscFunctionReturn(PETSC_SUCCESS);
940: }