From d5211ffd6baad6d5c46ea9158c2d9ac6c8a40c2e Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 11 Jul 2006 15:28:51 +0000 Subject: [PATCH] initial import --- Cantera/src/DAE_Solver.h | 217 +++++++++++++++++++++++++++++++ Cantera/src/IDA_Solver.cpp | 256 +++++++++++++++++++++++++++++++++++++ 2 files changed, 473 insertions(+) create mode 100644 Cantera/src/DAE_Solver.h create mode 100644 Cantera/src/IDA_Solver.cpp diff --git a/Cantera/src/DAE_Solver.h b/Cantera/src/DAE_Solver.h new file mode 100644 index 000000000..bb40393a9 --- /dev/null +++ b/Cantera/src/DAE_Solver.h @@ -0,0 +1,217 @@ +/** + * + * @file DAE_Solver.h + * + * Header file for class DAE_Solver + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2006 California Institute of Technology + * + */ + +#undef DAE_DEVEL + +#ifndef CT_DAE_Solver_H +#define CT_DAE_Solver_H + +#include + +#include "ct_defs.h" +#include "ResidEval.h" +#include "global.h" + +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 BandedJacobian : public Jacobian { + public: + BandedJacobian(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; + }; + + const int cDirect = 0; + const int cKrylov = 1; + + + + /** + * Wrapper for DAE solvers + */ + class DAE_Solver { + public: + + DAE_Solver(ResidEval& f) : m_resid(f), + m_neq(f.nEquations()), + m_time(0.0) {} + + virtual ~DAE_Solver(){} + + /** + * Set error tolerances. This version specifies a scalar + * relative tolerance, and a vector absolute tolerance. + */ + virtual void setTolerances(doublereal reltol, + doublereal* abstol) { + warn("setTolerances"); + } + + /** + * Set error tolerances. This version specifies a scalar + * relative tolerance, and a scalar absolute tolerance. + */ + virtual void setTolerances(doublereal reltol, doublereal abstol) { + warn("setTolerances"); + } + + /** + * Specify a Jacobian evaluator. If this method is not called, + * the Jacobian will be computed by finite difference. + */ + void setJacobian(Jacobian& jac) { + warn("setJacobian"); + } + + virtual void setLinearSolverType(int solverType) { + warn("setLinearSolverType"); + } + + virtual void setDenseLinearSolver() { + warn("setDenseLinearSolver"); + } + + virtual void setBandedLinearSolver(int m_upper, int m_lower) { + warn("setBandedLinearSolver"); + } + virtual void setMaxTime(doublereal tmax) { + warn("setMaxTime"); + } + virtual void setMaxStepSize(doublereal dtmax) { + warn("setMaxStepSize"); + } + virtual void setMaxOrder(int n) { + warn("setMaxOrder"); + } + virtual void setMaxNumSteps(int n) { + warn("setMaxNumSteps"); + } + virtual void setInitialStepSize(doublereal h0) { + warn("setInitialStepSize"); + } + virtual void setStopTime(doublereal tstop) { + warn("setStopTime"); + } + virtual void setMaxErrTestFailures(int n) { + warn("setMaxErrTestFailures"); + } + virtual void setMaxNonlinIterations(int n) { + warn("setMaxNonlinIterations"); + } + virtual void setMaxNonlinConvFailures(int n) { + warn("setMaxNonlinConvFailures"); + } + virtual void inclAlgebraicInErrorTest(bool yesno) { + warn("inclAlgebraicInErrorTest"); + } + + virtual void correctInitial_Y_given_Yp() { + warn("correctInitial_Y_given_Yp"); + } + + virtual void correctInitial_YaYp_given_Yd() { + warn("correctInitial_YaYp_given_Yd"); + } + + /** + * Solve the system of equations up to time tout. + */ + virtual int solve(doublereal tout) { + warn("solve"); return 0; + } + + /** + * Take one internal step. + */ + virtual int step(doublereal tout) { + warn("step"); return 0; + } + + /// Number of equations. + int nEquations() const { return m_resid.nEquations(); } + + /** + * initialize. Base class method does nothing. + */ + virtual void init(doublereal t0) {} + + /** + * Set a solver-specific input parameter. + */ + virtual void setInputParameter(int flag, doublereal value) { + warn("setInputParameter"); + } + + /** + * Get the value of a solver-specific output parameter. + */ + virtual doublereal getOutputParameter(int flag) const { + warn("getOutputParameter"); return 0.0; + } + + /// the current value of solution component k. + virtual doublereal solution(int k) const { + warn("solution"); return 0.0; + } + + virtual const doublereal* solutionVector() const { + warn("solutionVector"); return &m_dummy; + } + + /// the current value of the derivative of solution component k. + virtual doublereal derivative(int k) const { + warn("derivative"); return 0.0; + } + + virtual const doublereal* derivativeVector() const { + warn("derivativeVector"); return &m_dummy; + } + + protected: + + doublereal m_dummy; + + ResidEval& m_resid; + + integer m_neq; + doublereal m_time; + + + private: + void warn(string msg) const { + writelog(">>>> Warning: method "+msg+" of base class " + +"DAE_Solver called. Nothing done.\n"); + } + }; + +} + +#endif diff --git a/Cantera/src/IDA_Solver.cpp b/Cantera/src/IDA_Solver.cpp new file mode 100644 index 000000000..0458477b3 --- /dev/null +++ b/Cantera/src/IDA_Solver.cpp @@ -0,0 +1,256 @@ + +/** + * @file IDA_Solver.cpp + * + */ + +// Copyright 2006 California Institute of Technology + +#include "IDA_Solver.h" +#include "stringUtils.h" + +#include +using namespace std; + +#include +#include +#include +#include +#include +#include +#include + +inline static N_Vector nv(void* x) { + return reinterpret_cast(x); +} + +namespace Cantera { + + /** + * A simple class to hold an array of parameter values and a pointer to + * an instance of a subclass of ResidEval. + */ + class ResidData { + + + + + public: + ResidData(ResidEval* f, int npar = 0) { + m_func = f; + } + virtual ~ResidData() {} + ResidEval* m_func; + }; +} + + +extern "C" { + + /** + * Function called by IDA to evaluate the residual, given y and + * ydot. IDA allows passing in a void* pointer to access + * external data. Instead of requiring the user to provide a + * residual function directly to IDA (which would require using + * 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. + */ + static int ida_resid(realtype t, N_Vector y, N_Vector ydot, + N_Vector r, void *f_data) { + double* ydata = NV_DATA_S(y); + double* ydotdata = NV_DATA_S(ydot); + double* rdata = NV_DATA_S(r); + Cantera::ResidData* d = (Cantera::ResidData*)f_data; + Cantera::ResidEval* f = d->m_func; + f->eval(t, ydata, ydotdata, rdata); + return 0; + } + +} + +namespace Cantera { + + + /** + * Constructor. Default settings: dense jacobian, no user-supplied + * Jacobian function, Newton iteration. + */ + IDA_Solver::IDA_Solver(ResidEval& f) : DAE_Solver(f), + m_neq(0), + m_ida_mem(0), + m_t0(0.0), + m_y(0), + m_ydot(0), + m_abstol(0), + m_type(0), + m_itol(IDA_SS), + m_iter(0), + m_maxord(0), + m_reltol(1.e-9), + m_abstols(1.e-15), + m_nabs(0), + m_hmax(0.0), + m_maxsteps(20000), + m_mupper(0), + m_mlower(0) {} + + + /// Destructor. + IDA_Solver::~IDA_Solver() + { + if (m_ida_mem) { + IDAFree(&m_ida_mem); + } + if (m_y) N_VDestroy_Serial(nv(m_y)); + if (m_ydot) N_VDestroy_Serial(nv(m_ydot)); + if (m_abstol) N_VDestroy_Serial(nv(m_abstol)); + delete m_fdata; + } + + doublereal IDA_Solver::solution(int k) const { + return NV_Ith_S(nv(m_y),k); + } + + const doublereal* IDA_Solver::solutionVector() const { return NV_DATA_S(nv(m_y));} + + doublereal IDA_Solver::derivative(int k) const { + return NV_Ith_S(nv(m_ydot),k); + } + + const doublereal* IDA_Solver::derivativeVector() const { return NV_DATA_S(nv(m_ydot));} + + + void IDA_Solver::setTolerances(double reltol, double* abstol) { + m_itol = IDA_SV; + if (m_abstol) N_VDestroy_Serial(nv(m_abstol)); + m_abstol = reinterpret_cast(N_VNew_Serial(m_neq)); + for (int i=0; i < m_neq; i++) { + NV_Ith_S(nv(m_abstol), i) = abstol[i]; + } + m_reltol = reltol; + } + + void IDA_Solver::setTolerances(double reltol, double abstol) { + m_itol = IDA_SS; + m_reltol = reltol; + m_abstols = abstol; + } + + void IDA_Solver::setLinearSolverType(int solverType) { + m_type = solverType; + } + + void IDA_Solver::init(double t0) + { + m_t0 = t0; + + if (m_y) N_VDestroy_Serial(nv(m_y)); + if (m_ydot) N_VDestroy_Serial(nv(m_ydot)); + if (m_id) N_VDestroy_Serial(nv(m_id)); + if (m_constraints) N_VDestroy_Serial(nv(m_constraints)); + + m_y = reinterpret_cast(N_VNew_Serial(m_neq)); + m_ydot = reinterpret_cast(N_VNew_Serial(m_neq)); + m_constraints = reinterpret_cast(N_VNew_Serial(m_neq)); + + for (int i=0; i 0) + // flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord); + //if (m_maxsteps > 0) + // flag = CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps); + //if (m_hmax > 0) + // flag = CVodeSetMaxStep(m_cvode_mem, m_hmax); + } + + void IDA_Solver::solve(double tout) + { + double t; + int flag; + flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_NORMAL); + if (flag != IDA_SUCCESS) + throw IDA_Err(" IDA error encountered."); + } + + double IDA_Solver::step(double tout) + { + double t; + int flag; + flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_ONE_STEP); + if (flag != IDA_SUCCESS) + throw IDA_Err(" IDA error encountered."); + return t; + } + + doublereal IDA_Solver::getOutputParameter(int flag) { + switch (flag) { + case REAL_WORKSPACE_SIZE: + flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw); + return doublereal(lenrw); + } + +} + +