From 2f971807535f55e63221c0d57a4e8fa7a7a7c865 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sun, 31 Jul 2011 23:44:45 +0000 Subject: [PATCH] Fixed an error in the IDA_Solver class. Added an analytical jacobian capability to IDA_Solver. Added column pointers as a member to DenseMatrix --- Cantera/src/numerics/DenseMatrix.h | 2 +- Cantera/src/numerics/IDA_Solver.cpp | 130 ++++++++++++++++++++++--- Cantera/src/numerics/IDA_Solver.h | 32 ++++++ Cantera/src/numerics/ResidJacEval.cpp | 30 +++++- Cantera/src/numerics/ResidJacEval.h | 21 ++++ Cantera/src/thermo/MixtureFugacityTP.h | 2 +- 6 files changed, 200 insertions(+), 17 deletions(-) diff --git a/Cantera/src/numerics/DenseMatrix.h b/Cantera/src/numerics/DenseMatrix.h index f635e17bb..b6d37af1c 100644 --- a/Cantera/src/numerics/DenseMatrix.h +++ b/Cantera/src/numerics/DenseMatrix.h @@ -113,7 +113,7 @@ namespace Cantera { * @param m New number of columns * @param v Default fill value. defaults to zero. */ - void resize(int n, int m, doublereal v); + void resize(int n, int m, doublereal v = 0.0); //! Return a vector of const pointers to the columns /*! diff --git a/Cantera/src/numerics/IDA_Solver.cpp b/Cantera/src/numerics/IDA_Solver.cpp index 6f2575138..91a5532a1 100644 --- a/Cantera/src/numerics/IDA_Solver.cpp +++ b/Cantera/src/numerics/IDA_Solver.cpp @@ -44,14 +44,16 @@ namespace Cantera { public: - ResidData(ResidJacEval* f, int npar = 0) { + ResidData(ResidJacEval* f, IDA_Solver *s, int npar = 0) { m_func = f; + m_solver = s; } virtual ~ResidData() { } ResidJacEval* m_func; + IDA_Solver * m_solver; }; } @@ -71,9 +73,41 @@ extern "C" { double* rdata = NV_DATA_S(r); 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; } + + //! Function called by by IDA to evaluate the Jacobian, given y and ydot. + /*! + * + * + * typedef int (*IDADlsDenseJacFn)(int N, realtype t, realtype c_j, + * N_Vector y, N_Vector yp, N_Vector r, + * DlsMat Jac, void *user_data, + * N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); + * + * A IDADlsDenseJacFn should return + * 0 if successful, + * a positive int if a recoverable error occurred, or + * a negative int if a nonrecoverable error occurred. + * In the case of a recoverable error return, the integrator will + * attempt to recover by reducing the stepsize (which changes cj). + */ + static int ida_jacobian(int nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r, + DlsMat Jac, void *f_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) { + doublereal * ydata = NV_DATA_S(y); + doublereal * ydotdata = NV_DATA_S(ydot); + doublereal * rdata = NV_DATA_S(r); + Cantera::ResidData* d = (Cantera::ResidData*) f_data; + Cantera::ResidJacEval* f = d->m_func; + doublereal * const * colPts = Jac->cols; + doublereal delta_t = 0.0; + f->evalJacobianDP(t, delta_t, c_j, ydata, ydotdata, colPts, rdata); + return 0; + } + + } namespace Cantera { @@ -103,7 +137,12 @@ namespace Cantera { m_h0(0.0), m_maxsteps(20000), m_maxord(0), + m_formJac(0), m_tstop(0.0), + m_told_old(0.0), + m_told(0.0), + m_tcurrent(0.0), + m_deltat(0.0), m_maxErrTestFails(-1), m_maxNonlinIters(0), m_maxNonlinConvFails(-1), @@ -151,10 +190,12 @@ namespace Cantera { for (int i = 0; i < m_neq; i++) { NV_Ith_S(nv(m_abstol), i) = abstol[i]; } - m_reltol = reltol; - int flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol)); - if (flag != IDA_SUCCESS) { - throw IDA_Err("Memory allocation failed."); + m_reltol = reltol; + if (m_ida_mem) { + int flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol)); + if (flag != IDA_SUCCESS) { + throw IDA_Err("Memory allocation failed."); + } } } //==================================================================================================================== @@ -162,9 +203,11 @@ namespace Cantera { m_itol = IDA_SS; m_reltol = reltol; m_abstols = abstol; - int flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols); - if (flag != IDA_SUCCESS) { - throw IDA_Err("Memory allocation failed."); + if (m_ida_mem) { + int flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols); + if (flag != IDA_SUCCESS) { + throw IDA_Err("Memory allocation failed."); + } } } //==================================================================================================================== @@ -198,6 +241,18 @@ namespace Cantera { m_tstop = tstop; } //==================================================================================================================== + void IDA_Solver::setJacobianType(int formJac) { + m_formJac = formJac; + if (m_ida_mem) { + if (m_formJac == 1) { + int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian); + if (flag != IDA_SUCCESS) { + throw IDA_Err("IDADlsSetDenseJacFn failed."); + } + } + } + } + //==================================================================================================================== void IDA_Solver::setMaxErrTestFailures(int maxErrTestFails) { m_maxErrTestFails = maxErrTestFails; } @@ -222,6 +277,9 @@ namespace Cantera { void IDA_Solver::init(doublereal t0) { m_t0 = t0; + m_told = t0; + m_told_old = t0; + m_tcurrent = t0; if (m_y) { N_VDestroy_Serial(nv(m_y)); } @@ -338,9 +396,15 @@ namespace Cantera { throw IDA_Err("unsupported linear solver type"); } + if (m_formJac == 1) { + flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian); + if (flag != IDA_SUCCESS) { + throw IDA_Err("IDADlsSetDenseJacFn failed."); + } + } // pass a pointer to func in m_data - m_fdata = new ResidData(&m_resid, m_resid.nparams()); + m_fdata = new ResidData(&m_resid, this, m_resid.nparams()); #if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23) flag = IDASetRdata(m_ida_mem, (void*)m_fdata); if (flag != IDA_SUCCESS) { @@ -383,7 +447,7 @@ namespace Cantera { throw IDA_Err("IDASetMaxErrTestFails failed."); } } - if (m_maxNonlinIters >= 0) { + if (m_maxNonlinIters > 0) { flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters); if (flag != IDA_SUCCESS) { throw IDA_Err("IDASetmaxNonlinIters failed."); @@ -490,11 +554,35 @@ namespace Cantera { //==================================================================================================================== int IDA_Solver::solve(double tout) { - double t; + double tretn; int flag; - flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_NORMAL); - if (flag != IDA_SUCCESS) + flag = IDASetStopTime(m_ida_mem, tout); + if (flag != IDA_SUCCESS) { throw IDA_Err(" IDA error encountered."); + } + do { + if (tout <= m_tcurrent) { + throw IDA_Err(" tout <= tcurrent"); + } + m_told_old = m_told; + m_told = m_tcurrent; + flag = IDASolve(m_ida_mem, tout, &tretn, nv(m_y), nv(m_ydot), IDA_ONE_STEP); + if (flag < 0) { + throw IDA_Err(" IDA error encountered."); + } else if (flag == IDA_TSTOP_RETURN) { + // we've reached our goal, and have actually integrated past it + } else if (flag == IDA_ROOT_RETURN) { + // not sure what to do with this yet + } else if (flag == IDA_WARNING) { + throw IDA_Err(" IDA Warning encountered."); + } + m_tcurrent = tretn; + m_deltat = m_tcurrent - m_told; + } while (tretn < tout); + + if (flag != IDA_SUCCESS && flag != IDA_TSTOP_RETURN) { + throw IDA_Err(" IDA error encountered."); + } return flag; } //==================================================================================================================== @@ -502,9 +590,23 @@ namespace Cantera { { double t; int flag; + if (tout <= m_tcurrent) { + throw IDA_Err(" tout <= tcurrent"); + } + m_told_old = m_told; + m_told = m_tcurrent; flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_ONE_STEP); - if (flag != IDA_SUCCESS) + if (flag < 0) { throw IDA_Err(" IDA error encountered."); + } else if (flag == IDA_TSTOP_RETURN) { + // we've reached our goal, and have actually integrated past it + } else if (flag == IDA_ROOT_RETURN) { + // not sure what to do with this yet + } else if (flag == IDA_WARNING) { + throw IDA_Err(" IDA Warning encountered."); + } + m_tcurrent = t; + m_deltat = m_tcurrent - m_told; return t; } //==================================================================================================================== diff --git a/Cantera/src/numerics/IDA_Solver.h b/Cantera/src/numerics/IDA_Solver.h index 03d60c904..4bf250161 100644 --- a/Cantera/src/numerics/IDA_Solver.h +++ b/Cantera/src/numerics/IDA_Solver.h @@ -110,6 +110,17 @@ namespace Cantera { virtual void setStopTime(doublereal tstop); + //! Set the form of the jacobian + /*! + * + * @param formJac Form of the jacobian + * + * 0 numerical jacobian + * 1 analytical jacobian given by the evalJacobianDP() function + */ + virtual void setJacobianType(int formJac); + + virtual void setMaxErrTestFailures(int n); //! Set the maximum number of nonlinear iterations on a timestep @@ -289,9 +300,30 @@ namespace Cantera { //! maximum time step order of the method int m_maxord; + //! Form of the jacobian + /*! + * 0 numerical jacobian created by ida + * 1 analytical jacobian. Must have populated the evalJacobianDP() + * function in the ResidJacEval class. + * 2 numerical jacobian formed by the ResidJacEval class (unimplemented) + */ + int m_formJac; + //! maximum time doublereal m_tstop; + //! Value of the previous, previous time + doublereal m_told_old; + + //! Value of the previous time + doublereal m_told; + + //! Value of the current time + doublereal m_tcurrent; + + //! Value of deltaT for the current step + doublereal m_deltat; + //! maximum number of error test failures int m_maxErrTestFails; diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp index ea5ad953d..247f46e5c 100644 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -334,7 +334,35 @@ namespace Cantera { SquareMatrix &J, doublereal * const resid) { - throw CanteraError("ResidJacEval::evalJacobian()", "Not implemented\n"); + 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); + } + //==================================================================================================================== + // Calculate an analytical jacobian and the residual at the current time and values. + /* + * Only called if the jacFormation method is set to analytical + * + * @param t Time (input) + * @param delta_t The current value of the time step (input) + * @param c_j The current value of the coefficient of the time derivative + * @param y Solution vector (input, do not modify) + * @param ydot Rate of change of solution vector. (input, do not modify) + * @param jac_colPts Reference to the SquareMatrix object to be calculated (output) + * @param resid Value of the residual that is computed (output) + */ + int ResidJacEval:: + evalJacobianDP(const doublereal t, const doublereal delta_t, + const doublereal c_j, + const doublereal * const y, + const doublereal * const ydot, + doublereal * const * jac_colPts, + doublereal * const resid) + { + throw CanteraError("ResidJacEval::evalJacobianDP()", "Not implemented\n"); return 1; } //==================================================================================================================== diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index 5277388fb..b44707ec9 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -333,6 +333,27 @@ namespace Cantera { SquareMatrix &J, doublereal * const resid); + //! Calculate an analytical jacobian and the residual at the current time and values. + /*! + * Only called if the jacFormation method is set to analytical + * + * @param t Time (input) + * @param delta_t The current value of the time step (input) + * @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) + * @param resid Value of the residual that is computed (output) + * + * @return Returns a flag to indicate that operation is successful. + * 1 Means a successful operation + * -0 or neg value Means an unsuccessful operation + */ + virtual int evalJacobianDP(const doublereal t, const doublereal delta_t, doublereal cj, + const doublereal* const y, + const doublereal* const ydot, + doublereal * const *jacobianColPts, + doublereal * const resid); + protected: //! constant value of atol diff --git a/Cantera/src/thermo/MixtureFugacityTP.h b/Cantera/src/thermo/MixtureFugacityTP.h index 77e88b061..f5d603334 100644 --- a/Cantera/src/thermo/MixtureFugacityTP.h +++ b/Cantera/src/thermo/MixtureFugacityTP.h @@ -55,7 +55,7 @@ namespace Cantera { * steps for efficiently handling mixture of gases that whose standard states * are defined as ideal gases, but which describe also non-ideal solutions. * In addition a multicomponent liquid phase below the critical temperature of the - * mixture is also allowed. The main subclass will be a mixture Redlich-kwong class. + * mixture is also allowed. The main subclass is currently a mixture Redlich-Kwong class. * * Several concepts are introduced. The first concept is there are temporary * variables for holding the species standard state values