CVode returns wrong yout if h becomes zero

[dweindl]

Hi, I stumbled across a problem in CVode which seems to yield a wrong yout if h becomes zero.

After CVode: Internal t = 0 and h = 0 are such that t + h = t on the next step. The solver will continue anyway. and CVodeGetDky: Illegal value for t.t = 900 is not between tcur - hu = 0 and tcur = 0, CVode still returns CV_SUCCESS, but yout is not the correct value for tout, but for the time when h became zero.

The problem seems to be the multiplication with cv_mem->cv_h for checking whether tout was reached and subsequently ignoring the return value of CVodeGetDky (CV_BAD_T), here:

sundials/src/cvodes/cvodes.c

Lines 3298 to 3305 in a9bc74f

	/* In NORMAL mode, check if tout reached */
	if ( (itask == CV_NORMAL) && (cv_mem->cv_tn-tout)*cv_mem->cv_h >= ZERO ) {
	istate = CV_SUCCESS;
	cv_mem->cv_tretlast = *tret = tout;
	(void) CVodeGetDky(cv_mem, tout, 0, yout);
	cv_mem->cv_next_q = cv_mem->cv_qprime;
	cv_mem->cv_next_h = cv_mem->cv_hprime;
	break;

Similar checks exist in CVodeF here:

sundials/src/cvodes/cvodea.c

Lines 560 to 561 in a9bc74f

	/* Return if tout reached */
	if ( (*tret - tout)*cv_mem->cv_h >= ZERO ) {

and here:

sundials/src/cvodes/cvodea.c

Lines 472 to 473 in a9bc74f

	if ((ttest - tout)*cv_mem->cv_h >= ZERO) {
	/* ttest is after tout, interpolate to tout */

CVodeF leads to the same warnings, but returns CV_BAD_T. That is at least not CV_SUCCESS, but it is somewhat misleading, as CVodeGetDky should not have been called in the first place.

This happens in a more complex application and I am not sure I can easily provide a minimal reproducible example.