diff --git a/Cantera/src/numerics/BEulerInt.cpp b/Cantera/src/numerics/BEulerInt.cpp index 154a89d7f..78753267d 100644 --- a/Cantera/src/numerics/BEulerInt.cpp +++ b/Cantera/src/numerics/BEulerInt.cpp @@ -656,7 +656,7 @@ namespace Cantera { /******************************************************************** * Call the function to get a jacobian. */ - m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f); + m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f); #ifdef DEBUG_HKM //double dddd = J(89, 89); //checkFinite(dddd); diff --git a/Cantera/src/numerics/IDA_Solver.cpp b/Cantera/src/numerics/IDA_Solver.cpp index 91a5532a1..7c269b040 100644 --- a/Cantera/src/numerics/IDA_Solver.cpp +++ b/Cantera/src/numerics/IDA_Solver.cpp @@ -66,6 +66,13 @@ extern "C" { * the sundials data types N_Vector, etc.), we define this function as the single function that IDA always calls. The * real evaluation of the residual is done by an instance of a subclass of ResidEval, passed in to this * function as a pointer in the parameters. + * + * FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its residual routines ------ + * A IDAResFn res should return a value of 0 if successful, a positive + * value if a recoverable error occured (e.g. yy has an illegal value), + * or a negative value if a nonrecoverable error occured. In the latter + * case, the program halts. If a recoverable error occured, the integrator + * will attempt to correct and retry. */ static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void *f_data) { double* ydata = NV_DATA_S(y); @@ -74,8 +81,15 @@ extern "C" { Cantera::ResidData* d = (Cantera::ResidData*) f_data; Cantera::ResidJacEval* f = d->m_func; Cantera::IDA_Solver *s = d->m_solver; - f->eval(t, ydata, ydotdata, rdata); - return 0; + double delta_t = s->getCurrentStepFromIDA(); + // TODO evaluate evalType. Assumed to be Base_ResidEval + int retn = 0; + int flag = f->evalResidNJ(t, delta_t, ydata, ydotdata, rdata); + if (flag < 0) { + // This signals to IDA that a nonrecoverable error has occurred. + retn = flag; + } + return retn; } //! Function called by by IDA to evaluate the Jacobian, given y and ydot. @@ -102,7 +116,9 @@ extern "C" { Cantera::ResidData* d = (Cantera::ResidData*) f_data; Cantera::ResidJacEval* f = d->m_func; doublereal * const * colPts = Jac->cols; - doublereal delta_t = 0.0; + Cantera::IDA_Solver *s = d->m_solver; + double delta_t = s->getCurrentStepFromIDA(); + // printf(" delta_t = %g 1/cj = %g\n", delta_t, 1.0/c_j); f->evalJacobianDP(t, delta_t, c_j, ydata, ydotdata, colPts, rdata); return 0; } @@ -241,6 +257,12 @@ namespace Cantera { m_tstop = tstop; } //==================================================================================================================== + doublereal IDA_Solver::getCurrentStepFromIDA() { + doublereal hcur; + IDAGetCurrentStep(m_ida_mem, &hcur); + return hcur; + } + //==================================================================================================================== void IDA_Solver::setJacobianType(int formJac) { m_formJac = formJac; if (m_ida_mem) { diff --git a/Cantera/src/numerics/IDA_Solver.h b/Cantera/src/numerics/IDA_Solver.h index 4bf250161..928728969 100644 --- a/Cantera/src/numerics/IDA_Solver.h +++ b/Cantera/src/numerics/IDA_Solver.h @@ -109,6 +109,12 @@ namespace Cantera { */ virtual void setStopTime(doublereal tstop); + //! Get the current step size from IDA via a call + /*! + * @return Returns the current step size. + */ + virtual double getCurrentStepFromIDA(); + //! Set the form of the jacobian /*! diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 64394870f..89a2904a0 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -2979,6 +2979,7 @@ namespace Cantera { int info; doublereal ysave, ydotsave, dy; int retn = 1; + /* * Clear the factor flag */ @@ -2987,7 +2988,7 @@ namespace Cantera { /******************************************************************** * Call the function to get a jacobian. */ - info = m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f); + info = m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f); m_nJacEval++; m_nfe++; if (info != 1) { diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp index 247f46e5c..a0261c9f3 100644 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -328,17 +328,13 @@ namespace Cantera { * @param resid Value of the residual that is computed (output) */ int ResidJacEval:: - evalJacobian(const doublereal t, const doublereal delta_t, + evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj, const doublereal * const y, const doublereal * const ydot, SquareMatrix &J, doublereal * const resid) { doublereal * const * jac_colPts = J.colPts(); - doublereal cj = 0.0; - if (delta_t > 0.0) { - cj = 1.0/delta_t; - } return evalJacobianDP(t, delta_t, cj, y, ydot, jac_colPts, resid); } //==================================================================================================================== diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index b44707ec9..678d3fe7b 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -318,6 +318,7 @@ namespace Cantera { * * @param t Time (input) * @param delta_t The current value of the time step (input) + * @param cj Coefficient of yprime used in the evalulation of the jacobian * @param y Solution vector (input, do not modify) * @param ydot Rate of change of solution vector. (input, do not modify) * @param J Reference to the SquareMatrix object to be calculated (output) @@ -327,11 +328,9 @@ namespace Cantera { * 1 Means a successful operation * -0 or neg value Means an unsuccessful operation */ - virtual int evalJacobian(const doublereal t, const doublereal delta_t, - const doublereal* const y, - const doublereal* const ydot, - SquareMatrix &J, - doublereal * const resid); + virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj, + const doublereal* const y, const doublereal* const ydot, + SquareMatrix &J, doublereal * const resid); //! Calculate an analytical jacobian and the residual at the current time and values. /*! @@ -339,6 +338,7 @@ namespace Cantera { * * @param t Time (input) * @param delta_t The current value of the time step (input) + * @param cj Coefficient of yprime used in the evalulation of the jacobian * @param y Solution vector (input, do not modify) * @param ydot Rate of change of solution vector. (input, do not modify) * @param J Reference to the SquareMatrix object to be calculated (output)