diff --git a/Cantera/src/numerics/DAE_Solver.h b/Cantera/src/numerics/DAE_Solver.h index 22105df57..4d570cca7 100644 --- a/Cantera/src/numerics/DAE_Solver.h +++ b/Cantera/src/numerics/DAE_Solver.h @@ -5,7 +5,7 @@ * Header file for class DAE_Solver */ -/* $Author$ +/* * $Date$ * $Revision$ * @@ -26,6 +26,7 @@ namespace Cantera { +#ifdef DAE_DEVEL /** * @defgroup numerics Numerical Utilities within Cantera * @@ -218,6 +219,8 @@ namespace Cantera { } }; +#endif + } #endif diff --git a/Cantera/src/numerics/DASPK.cpp b/Cantera/src/numerics/DASPK.cpp index 5fedb370d..301131946 100755 --- a/Cantera/src/numerics/DASPK.cpp +++ b/Cantera/src/numerics/DASPK.cpp @@ -2,6 +2,10 @@ * @file DASPK.cpp * */ +/* + * $Date$ + * $Revision$ + */ // Copyright 2001 California Institute of Technology @@ -20,262 +24,262 @@ using namespace std; extern "C" { - typedef void (*ResidFunc)(const doublereal* t, - const doublereal* y, const doublereal* yprime, - const doublereal* cj, doublereal* delta, - integer* ires, doublereal* rpar, integer* ipar); + typedef void (*ResidFunc)(const doublereal* t, + const doublereal* y, const doublereal* yprime, + const doublereal* cj, doublereal* delta, + integer* ires, doublereal* rpar, integer* ipar); - typedef void (*JacFunc)(); - typedef void (*PsolFunc)(); + typedef void (*JacFunc)(); + typedef void (*PsolFunc)(); - extern void ddaspk_(ResidFunc res, integer* neq, doublereal* t, - doublereal* y, doublereal* yprime, doublereal* tout, integer* info, - doublereal* rtol, doublereal* atol, integer* idid, doublereal* rwork, - integer* lrw, integer* iwork, integer* liw, doublereal* rpar, - integer* ipar, JacFunc jac, PsolFunc psol); + extern void ddaspk_(ResidFunc res, integer* neq, doublereal* t, + doublereal* y, doublereal* yprime, doublereal* tout, integer* info, + doublereal* rtol, doublereal* atol, integer* idid, doublereal* rwork, + integer* lrw, integer* iwork, integer* liw, doublereal* rpar, + integer* ipar, JacFunc jac, PsolFunc psol); - /** - * Function called by DASPK to evaluate the residual. - */ - static void ddaspk_res(const doublereal* t, - const doublereal* y, const doublereal* yprime, - const doublereal* cj, doublereal* delta, - integer* ires, doublereal* rpar, integer* ipar) { - void **iddres_res = reinterpret_cast(&(ipar[0])); - void *hndl = *iddres_res; - Cantera::ResidEval* f = (Cantera::ResidEval*)hndl; - double delta_t = 0.0; - f->evalResid(*t, delta_t, y, yprime, delta); - } + /** + * Function called by DASPK to evaluate the residual. + */ + static void ddaspk_res(const doublereal* t, + const doublereal* y, const doublereal* yprime, + const doublereal* cj, doublereal* delta, + integer* ires, doublereal* rpar, integer* ipar) { + void **iddres_res = reinterpret_cast(&(ipar[0])); + void *hndl = *iddres_res; + Cantera::ResidEval* f = (Cantera::ResidEval*)hndl; + double delta_t = 0.0; + f->evalResid(*t, delta_t, y, yprime, delta); + } - static void ddaspk_jac() {} - static void ddaspk_psol() {} + static void ddaspk_jac() {} + static void ddaspk_psol() {} } namespace Cantera { - class DASPKErr : public CanteraError { - public: - DASPKErr(string proc, string msg) - : CanteraError("DASPK::"+proc,msg) {} - virtual ~DASPKErr(){} - }; + class DASPKErr : public CanteraError { + public: + DASPKErr(string proc, string msg) + : CanteraError("DASPK::"+proc,msg) {} + virtual ~DASPKErr(){} + }; - /** - * Constructor. Default settings: dense jacobian, no user-supplied - * Jacobian function, Newton iteration. - */ - DASPK::DASPK(ResidEval& f) : - m_resid(f), - m_idid(0), - m_lrw(0), - m_liw(0), - m_ml(0), - m_mu(0), - m_lenwp(0), - m_ok(false), - m_init(false), - m_time(0.0) - { - m_info.resize(20); - m_neq = f.neq(); - m_rwork.resize(20); // will be reset later - m_iwork.resize(20); // " - m_ipar.resize(2); - m_rpar.resize(2); - void *iddr = static_cast(&m_resid); - void **iddr_ipar = reinterpret_cast(&(m_ipar[0])); - *iddr_ipar = iddr; - setTolerances(1.e-7, 1.e-15); - } + /** + * Constructor. Default settings: dense jacobian, no user-supplied + * Jacobian function, Newton iteration. + */ + DASPK::DASPK(ResidEval& f) : + m_resid(f), + m_idid(0), + m_lrw(0), + m_liw(0), + m_ml(0), + m_mu(0), + m_lenwp(0), + m_ok(false), + m_init(false), + m_time(0.0) + { + m_info.resize(20); + m_neq = f.neq(); + m_rwork.resize(20); // will be reset later + m_iwork.resize(20); // " + m_ipar.resize(2); + m_rpar.resize(2); + void *iddr = static_cast(&m_resid); + void **iddr_ipar = reinterpret_cast(&(m_ipar[0])); + *iddr_ipar = iddr; + setTolerances(1.e-7, 1.e-15); + } - /// Destructor. - DASPK::~DASPK(){} + /// Destructor. + DASPK::~DASPK(){} - void DASPK::setTolerances(int nr, double* reltol, int na, double* abstol) { - // scalar tolerances - if (nr == 1 && na == 1) { - setInfo(2,0); - m_rtol.resize(1); - m_rtol[0] = reltol[0]; - m_atol.resize(1); - m_atol[0] = abstol[0]; - } - // vector tolerances - else { - setInfo(2,1); - m_rtol.resize(neq()); - m_atol.resize(neq()); - copy(reltol, reltol + nr, m_rtol.begin()); - copy(abstol, abstol + na, m_atol.begin()); - } + void DASPK::setTolerances(int nr, double* reltol, int na, double* abstol) { + // scalar tolerances + if (nr == 1 && na == 1) { + setInfo(2,0); + m_rtol.resize(1); + m_rtol[0] = reltol[0]; + m_atol.resize(1); + m_atol[0] = abstol[0]; } - - void DASPK::setTolerances(double reltol, double abstol) { - doublereal rtol = reltol; - doublereal atol = abstol; - setTolerances(1, &rtol, 1, &atol); + // vector tolerances + else { + setInfo(2,1); + m_rtol.resize(neq()); + m_atol.resize(neq()); + copy(reltol, reltol + nr, m_rtol.begin()); + copy(abstol, abstol + na, m_atol.begin()); } + } - void DASPK::setJacobian(Jacobian& jac) { + void DASPK::setTolerances(double reltol, double abstol) { + doublereal rtol = reltol; + doublereal atol = abstol; + setTolerances(1, &rtol, 1, &atol); + } + + void DASPK::setJacobian(Jacobian& jac) { - // No Jacobian evaluation function is supplied, so let DASPK - // compute the Jacobian by numerical finite-difference - if (!jac.supplied()) setInfo(5,0); - else { - setInfo(5,1); - } - - if (jac.isBanded()) { - setInfo(6,1); - setIwork(1, jac.lowerBandWidth()); - setIwork(2, jac.upperBandWidth()); - } - else setInfo(6,0); + // No Jacobian evaluation function is supplied, so let DASPK + // compute the Jacobian by numerical finite-difference + if (!jac.supplied()) setInfo(5,0); + else { + setInfo(5,1); } - void DASPK::setMethod(int methodType) { - if (methodType == cDirect) - setInfo(12,0); - else if (methodType == cKrylov) - setInfo(12,1); - else - throw DASPKErr("setMethod", - "method must be either cDirect " - "or cKrylov"); + if (jac.isBanded()) { + setInfo(6,1); + setIwork(1, jac.lowerBandWidth()); + setIwork(2, jac.upperBandWidth()); } + else setInfo(6,0); + } - void DASPK::setMaxTime(doublereal tmax) { - setInfo(4,1); - setRwork(1,tmax); + void DASPK::setMethod(int methodType) { + if (methodType == cDirect) + setInfo(12,0); + else if (methodType == cKrylov) + setInfo(12,1); + else + throw DASPKErr("setMethod", + "method must be either cDirect " + "or cKrylov"); + } + + void DASPK::setMaxTime(doublereal tmax) { + setInfo(4,1); + setRwork(1,tmax); + } + + void DASPK::setMaxStepSize(doublereal dtmax) { + setInfo(7,1); + setRwork(2,dtmax); + } + + void DASPK::setInitialIntStepSize(doublereal h0) { + setInfo(8,1); + setRwork(3,h0); + } + + void DASPK::setMaxOrder(int n) { + setInfo(9,1); + setIwork(3,n); + } + + void DASPK::estimateInitial_Y_given_Yp() { + setInfo(11,2); + } + + void DASPK::estimateInitial_YaYp_given_Yd( + const vector& vartypes) { + setInfo(11,2); + int m, n = neq(); + int lid = ((info(10) == 0 || info(10) == 2) ? + 41 : 41 + neq()); + if (int(m_iwork.size()) < lid + neq()) + m_iwork.resize(lid + neq()); + for (m = 0; m < n; m++) { + setIwork(lid + m, vartypes[m]); } + } - void DASPK::setMaxStepSize(doublereal dtmax) { - setInfo(7,1); - setRwork(2,dtmax); - } - - void DASPK::setInitialIntStepSize(doublereal h0) { - setInfo(8,1); - setRwork(3,h0); - } - - void DASPK::setMaxOrder(int n) { - setInfo(9,1); - setIwork(3,n); - } - - void DASPK::estimateInitial_Y_given_Yp() { - setInfo(11,2); - } - - void DASPK::estimateInitial_YaYp_given_Yd( - const vector& vartypes) { - setInfo(11,2); - int m, n = neq(); - int lid = ((info(10) == 0 || info(10) == 2) ? - 41 : 41 + neq()); - if (int(m_iwork.size()) < lid + neq()) - m_iwork.resize(lid + neq()); - for (m = 0; m < n; m++) { - setIwork(lid + m, vartypes[m]); - } - } - - void DASPK::sizeRwork() { - int base; - if (info(12) == 0) { - base = 50 + 9*neq(); - if (info(6) == 0) - base += neq()*neq(); - else { - base += (2*m_ml + m_mu + 1)*neq(); - if (info(5) == 0) - base += 2*(neq()/(m_ml + m_mu + 1) + 1); - } - } - else { - base = 91 + 18*neq() + m_lenwp; - } - if (info(16) == 1) base += neq(); - - /// @todo fix this! - base = 2000000; // tmp - - m_rwork.resize(base, 0.0); - m_lrw = base; - } - - void DASPK::sizeIwork() { - int base; - if (info(12) == 0) { - base = 40 + neq(); - } - else { - base = 40 + m_lenwp; - } - if (info(10) == 1 || info(10) == 3) base += neq(); - if (info(11) == 1 || info(16) == 1) base += neq(); - m_iwork.resize(base); - m_liw = base; - } - - - void DASPK::init(doublereal t0) - { - m_init = true; - m_time = t0; - setInfo(1,0); // tells DASPK to initialize - sizeRwork(); - sizeIwork(); - //m_resid.init(t0); - } - - int DASPK::integrate(doublereal tout) { - if (!m_init) init(0.0); - - doublereal tfinal = tout; - setInfo(3,0); // don't want intermediate output - - ddaspk_(ddaspk_res, &m_neq, &m_time, m_resid.solution(), - m_resid.solution_dot(), &tfinal, m_info.begin(), - m_rtol.begin(), m_atol.begin(), &m_idid, - m_rwork.begin(), &m_lrw, m_iwork.begin(), &m_liw, - m_rpar.begin(), m_ipar.begin(), ddaspk_jac, ddaspk_psol); - - return m_idid; - } - - void DASPK::step(double tout) - { - setInfo(3,1); // do want intermediate output - doublereal tfinal = tout; - // setInfo(3,0); // don't want intermediate output - - ddaspk_(ddaspk_res, &m_neq, &m_time, m_resid.solution(), - m_resid.solution_dot(), &tfinal, m_info.begin(), - m_rtol.begin(), m_atol.begin(), &m_idid, - m_rwork.begin(), &m_lrw, m_iwork.begin(), &m_liw, - m_rpar.begin(), m_ipar.begin(), ddaspk_jac, ddaspk_psol); - if (m_idid < 0) { - throw DASPKErr("step", - "DASPK returned IDID = "+int2str(m_idid)); - m_ok = false; - } - else if (m_idid == 1 || m_idid == 2 || m_idid == 3) { - m_ok = true; - } - else { - m_ok = false; - } - return; + void DASPK::sizeRwork() { + int base; + if (info(12) == 0) { + base = 50 + 9*neq(); + if (info(6) == 0) + base += neq()*neq(); + else { + base += (2*m_ml + m_mu + 1)*neq(); + if (info(5) == 0) + base += 2*(neq()/(m_ml + m_mu + 1) + 1); } + } + else { + base = 91 + 18*neq() + m_lenwp; + } + if (info(16) == 1) base += neq(); - int DASPK::nEvals() const { return iwork(12); } + /// @todo fix this! + base = 2000000; // tmp + + m_rwork.resize(base, 0.0); + m_lrw = base; + } + + void DASPK::sizeIwork() { + int base; + if (info(12) == 0) { + base = 40 + neq(); + } + else { + base = 40 + m_lenwp; + } + if (info(10) == 1 || info(10) == 3) base += neq(); + if (info(11) == 1 || info(16) == 1) base += neq(); + m_iwork.resize(base); + m_liw = base; + } + + + void DASPK::init(doublereal t0) + { + m_init = true; + m_time = t0; + setInfo(1,0); // tells DASPK to initialize + sizeRwork(); + sizeIwork(); + //m_resid.init(t0); + } + + int DASPK::integrate(doublereal tout) { + if (!m_init) init(0.0); + + doublereal tfinal = tout; + setInfo(3,0); // don't want intermediate output + + ddaspk_(ddaspk_res, &m_neq, &m_time, m_resid.solution(), + m_resid.solution_dot(), &tfinal, m_info.begin(), + m_rtol.begin(), m_atol.begin(), &m_idid, + m_rwork.begin(), &m_lrw, m_iwork.begin(), &m_liw, + m_rpar.begin(), m_ipar.begin(), ddaspk_jac, ddaspk_psol); + + return m_idid; + } + + void DASPK::step(double tout) + { + setInfo(3,1); // do want intermediate output + doublereal tfinal = tout; + // setInfo(3,0); // don't want intermediate output + + ddaspk_(ddaspk_res, &m_neq, &m_time, m_resid.solution(), + m_resid.solution_dot(), &tfinal, m_info.begin(), + m_rtol.begin(), m_atol.begin(), &m_idid, + m_rwork.begin(), &m_lrw, m_iwork.begin(), &m_liw, + m_rpar.begin(), m_ipar.begin(), ddaspk_jac, ddaspk_psol); + if (m_idid < 0) { + throw DASPKErr("step", + "DASPK returned IDID = "+int2str(m_idid)); + m_ok = false; + } + else if (m_idid == 1 || m_idid == 2 || m_idid == 3) { + m_ok = true; + } + else { + m_ok = false; + } + return; + } + + int DASPK::nEvals() const { return iwork(12); } } diff --git a/Cantera/src/numerics/DASPK.h b/Cantera/src/numerics/DASPK.h index 2bbca12e3..7ae5ccc0c 100755 --- a/Cantera/src/numerics/DASPK.h +++ b/Cantera/src/numerics/DASPK.h @@ -5,7 +5,7 @@ * Header file for class DASPK */ -/* $Author$ +/* * $Date$ * $Revision$ * @@ -23,92 +23,92 @@ namespace Cantera { - class Jacobian { - public: - Jacobian(){} - virtual ~Jacobian(){} - virtual bool supplied() { return false; } - virtual bool isBanded() { return false; } - virtual int lowerBandWidth() { return 0; } - virtual int upperBandWidth() { return 0; } - }; + class Jacobian { + public: + Jacobian(){} + virtual ~Jacobian(){} + virtual bool supplied() { return false; } + virtual bool isBanded() { return false; } + virtual int lowerBandWidth() { return 0; } + virtual int upperBandWidth() { return 0; } + }; - class BandedJac : public Jacobian { - public: - BandedJac(int ml, int mu) { - m_ml = ml; m_mu = mu; - } - virtual bool supplied() { return false; } - virtual bool isBanded() { return true; } - virtual int lowerBandWidth() { return m_ml; } - virtual int upperBandWidth() { return m_mu; } - protected: - int m_ml, m_mu; - }; + class BandedJac : public Jacobian { + public: + BandedJac(int ml, int mu) { + m_ml = ml; m_mu = mu; + } + virtual bool supplied() { return false; } + virtual bool isBanded() { return true; } + virtual int lowerBandWidth() { return m_ml; } + virtual int upperBandWidth() { return m_mu; } + protected: + int m_ml, m_mu; + }; - class ResidEval; + class ResidEval; - const int cDirect = 0; - const int cKrylov = 1; + const int cDirect = 0; + const int cKrylov = 1; - /** - * Wrapper for DASPK 2.0 DAE solver of Petzold et al. - */ - class DASPK { - public: + /** + * Wrapper for DASPK 2.0 DAE solver of Petzold et al. + */ + class DASPK { + public: - DASPK(ResidEval& f); - virtual ~DASPK(); + DASPK(ResidEval& f); + virtual ~DASPK(); - integer iwork(int n) const { return m_iwork[n-1];} - doublereal rwork(int n) const { return m_rwork[n-1];} - void setIwork(int n, integer m) { m_iwork[n-1] = m; } - void setRwork(int n, doublereal v) { m_rwork[n-1] = v; } - integer info(int n) { return m_info[n-1]; } - void setInfo(int n, integer m) { m_info[n-1] = m; } + integer iwork(int n) const { return m_iwork[n-1];} + doublereal rwork(int n) const { return m_rwork[n-1];} + void setIwork(int n, integer m) { m_iwork[n-1] = m; } + void setRwork(int n, doublereal v) { m_rwork[n-1] = v; } + integer info(int n) { return m_info[n-1]; } + void setInfo(int n, integer m) { m_info[n-1] = m; } - void setTolerances(int nr, double* reltol, int na, double* abstol); - void setTolerances(double reltol, double abstol); - void setJacobian(Jacobian& jac); - void setMethod(int methodType); - void setMaxTime(doublereal tmax); - void setMaxStepSize(doublereal dtmax); - void setMaxOrder(int n); - void setInitialIntStepSize(doublereal h0); - void estimateInitial_Y_given_Yp(); - void estimateInitial_YaYp_given_Yd(const std::vector& vartypes); - void sizeRwork(); - void sizeIwork(); - int integrate(doublereal tout); - void step(doublereal tout); - int nEvals() const; - int neq() { return m_resid.neq(); } - void init(doublereal t0); + void setTolerances(int nr, double* reltol, int na, double* abstol); + void setTolerances(double reltol, double abstol); + void setJacobian(Jacobian& jac); + void setMethod(int methodType); + void setMaxTime(doublereal tmax); + void setMaxStepSize(doublereal dtmax); + void setMaxOrder(int n); + void setInitialIntStepSize(doublereal h0); + void estimateInitial_Y_given_Yp(); + void estimateInitial_YaYp_given_Yd(const std::vector& vartypes); + void sizeRwork(); + void sizeIwork(); + int integrate(doublereal tout); + void step(doublereal tout); + int nEvals() const; + int neq() { return m_resid.nEquations(); } + void init(doublereal t0); - protected: + protected: - ResidEval& m_resid; - vector_int m_info; - vector_int m_iwork; - vector_int m_ipar; + ResidEval& m_resid; + vector_int m_info; + vector_int m_iwork; + vector_int m_ipar; - vector_fp m_rwork; - vector_fp m_atol; - vector_fp m_rtol; - vector_fp m_rpar; + vector_fp m_rwork; + vector_fp m_atol; + vector_fp m_rtol; + vector_fp m_rpar; - integer m_idid; - integer m_neq; - integer m_lrw; - integer m_liw; - integer m_ml; - integer m_mu; - integer m_lenwp; + integer m_idid; + integer m_neq; + integer m_lrw; + integer m_liw; + integer m_ml; + integer m_mu; + integer m_lenwp; - bool m_ok; - bool m_init; - doublereal m_time; - }; + bool m_ok; + bool m_init; + doublereal m_time; + }; }