From f325a1b1900c327a3df280bd7bda2f49166b36cc Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 11 Jul 2006 16:30:24 +0000 Subject: [PATCH] initial import --- Cantera/src/Crystal.h | 68 ++++++++++++++++++ Cantera/src/DAE_solvers.cpp | 30 ++++++++ Cantera/src/IDA_Solver.h | 137 ++++++++++++++++++++++++++++++++++++ Cantera/src/ResidEval.h | 100 ++++++++++++++++++++++++++ 4 files changed, 335 insertions(+) create mode 100644 Cantera/src/Crystal.h create mode 100644 Cantera/src/DAE_solvers.cpp create mode 100644 Cantera/src/IDA_Solver.h create mode 100755 Cantera/src/ResidEval.h diff --git a/Cantera/src/Crystal.h b/Cantera/src/Crystal.h new file mode 100644 index 000000000..dae466d77 --- /dev/null +++ b/Cantera/src/Crystal.h @@ -0,0 +1,68 @@ +/** + * @file Crystal.h + * + * $Author$ + * $Date$ + * $Revision$ + */ +#ifndef CT_CRYSTAL_H +#define CT_CRYSTAL_H + +#include "MultiPhase.h" + +namespace Cantera { + + /// A class for crystals. Each crystal consists of one or more + /// sublattices, each represented by an object of type + /// LatticePhase. + + class Crystal : public MultiPhase { + + public: + typedef LatticePhase lattice_t; + typedef vector lattice_list; + + /// Constructor. The constructor takes no arguments, since + /// phases are added using method addPhase. + Crystal() : MultiPhase() {} + + /// Destructor. Does nothing. Class MultiPhase does not take + /// "ownership" (i.e. responsibility for destroying) the + /// phase objects. + virtual ~Crystal() {} + + void addLattices(lattice_list& lattices, + const vector_fp& latticeSiteDensity); + + /// Add a phase to the mixture. + /// @param p pointer to the phase object + /// @param moles total number of moles of all species in this phase + void addLattice(lattice_t* lattice, doublereal siteDensity) { + MultiPhase::addPhase(lattice, siteDensity); + } + + /// Return a reference to phase n. The state of phase n is + /// also updated to match the state stored locally in the + /// mixture object. + lattice_t& lattice(index_t n) { + return *(lattice_t*)&phase(n); + } + + protected: + + + }; + + inline std::ostream& operator<<(std::ostream& s, Cantera::Crystal& x) { + size_t ip; + for (ip = 0; ip < x.nPhases(); ip++) { + s << "*************** Lattice " << ip << " *****************" << endl; + s << "SiteDensity: " << x.phaseMoles(ip) << endl; + + s << report(x.phase(ip)) << endl; + } + return s; + } +} + +#endif diff --git a/Cantera/src/DAE_solvers.cpp b/Cantera/src/DAE_solvers.cpp new file mode 100644 index 000000000..80d8a2a1e --- /dev/null +++ b/Cantera/src/DAE_solvers.cpp @@ -0,0 +1,30 @@ + +#include "ct_defs.h" +#include "DAE_Solver.h" + +#ifdef DAE_DEVEL + +#ifdef HAS_SUNDIALS +#include "IDA_Solver.cpp" +#endif + +namespace Cantera { + + DAE_Solver* newDAE_Solver(string itype) { + if (itype == "IDA") { +#ifdef HAS_SUNDIALS + return new IDA_Solver(); +#else + raise CanteraError("newDAE_Solver","IDA solver requires sundials" + " package, but Cantera was not built with sundials."); +#endif + } + else { + throw CanteraError("newDAE_Solver", + "unknown DAE solver: "+itype); + } + } +} + +# +#endif diff --git a/Cantera/src/IDA_Solver.h b/Cantera/src/IDA_Solver.h new file mode 100644 index 000000000..2820b723e --- /dev/null +++ b/Cantera/src/IDA_Solver.h @@ -0,0 +1,137 @@ +/** + * + * @file IDA_Solver.h + * + * Header file for class IDA_Solver + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2006 California Institute of Technology + * + */ + +#ifndef CT_IDA_Solver_H +#define CT_IDA_Solver_H + +#include + +#include "DAE_Solver.h" +#include "ctexceptions.h" + +namespace Cantera { + + /** + * Exception thrown when a IDA error is encountered. + */ + class IDA_Err : public CanteraError { + public: + IDA_Err(string msg) : CanteraError("IDA_Solver", msg){} + }; + + + class ResidData; // forward reference + + class IDA_Solver : public DAE_Solver { + public: + + IDA_Solver(ResidEval& f); + + virtual ~IDA_Solver(); + + /** + * Set error tolerances. This version specifies a scalar + * relative tolerance, and a vector absolute tolerance. + */ + virtual void setTolerances(doublereal reltol, + doublereal* abstol); + + /** + * Set error tolerances. This version specifies a scalar + * relative tolerance, and a scalar absolute tolerance. + */ + virtual void setTolerances(doublereal reltol, doublereal abstol); + + virtual void setLinearSolverType(int solverType); + + virtual void setDenseLinearSolver(); + virtual void setBandedLinearSolver(int m_upper, int m_lower); + + virtual void setMaxTime(doublereal tmax); + virtual void setMaxStepSize(doublereal dtmax); + + virtual void setMaxOrder(int n); + + virtual void setMaxNumSteps(int n); + virtual void setInitialStepSize(doublereal h0); + virtual void setMaxStepSize(doublereal hmax); + virtual void setStopTime(doublereal tstop); + virtual void setMaxErrTestFailures(int n); + virtual void setMaxNonlinIterations(int n); + virtual void setMaxNonlinConvFailures(int n); + virtual void inclAlgebraicInErrorTest(bool yesno); + + virtual void setInputParameter(int flag, doublereal value); + virtual doublereal getOutputParameter(int flag); + + + /** + * This method may be called if the initial conditions do not + * satisfy the residual equation F = 0. Given the derivatives + * of all variables, this method computes the initial y + * values. + */ + virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, + doublereal tout); + + /** + * This method may be called if the initial conditions do not + * satisfy the residual equation F = 0. Given the initial + * values of all differential variables, it computes the + * initial values of all algebraic variables and the initial + * derivatives of all differential variables. + */ + virtual void correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, + doublereal tout); + + + virtual int solve(doublereal tout); + + virtual int step(doublereal tout); + + virtual void init(doublereal t0); + + /// the current value of solution component k. + virtual doublereal solution(int k) const; + + virtual const doublereal* solutionVector() const; + + /// the current value of the derivative of solution component k. + virtual doublereal derivative(int k) const; + + virtual const doublereal* derivativeVector() const; + + protected: + + int m_neq; + void* m_ida_mem; + double m_t0; + void *m_y, *m_ydot, *m_id, *m_constraints, *m_abstol; + int m_type; + int m_itol; + int m_iter; + double m_reltol; + double m_abstols; + int m_nabs; + double m_hmax, m_hmin; + int m_maxsteps, m_maxord; + ResidData* m_fdata; + int m_mupper, m_mlower; + }; + +} + +#endif + diff --git a/Cantera/src/ResidEval.h b/Cantera/src/ResidEval.h new file mode 100755 index 000000000..ff3069037 --- /dev/null +++ b/Cantera/src/ResidEval.h @@ -0,0 +1,100 @@ +/** + * @file ResidEval.h + * + */ + +// Copyright 2006 California Institute of Technology + +#ifndef CT_RESIDEVAL_H +#define CT_RESIDEVAL_H + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ct_defs.h" + +namespace Cantera { + + const int c_NONE = 0; + const int c_GE_ZERO = 1; + const int c_GT_ZERO = 2; + const int c_LE_ZERO = -1; + const int c_LT_ZERO = -2; + + /** + * Virtual base class for DAE residual function evaluators. + * Classes derived from ResidEval evaluate the residual function + * \f[ + \vec{F}(t,\vec{y}, \vec{y^\prime}) + * \f] + * The DAE solver attempts to find a solution y(t) such that F = 0. + * @ingroup DAE_Group + */ + class ResidEval { + + public: + + ResidEval() {} + virtual ~ResidEval() {} + + /** + * Constrain solution component k. Possible values for + * 'flag' are: + * - c_NONE no constraint + * - c_GE_ZERO >= 0 + * - c_GT_ZERO > 0 + * - c_LE_ZERO <= 0 + * - c_LT_ZERO < 0 + */ + virtual void constrain(int k, int flag) { m_constrain[k] = flag; } + int constraint(int k) { return m_constrain[k]; } + + /** + * Specify that solution component k is purely algebraic - + * that is, the derivative of this component does not appear + * in the residual function. + */ + virtual void setAlgebraic(int k) { m_alg[k] = 1; } + virtual bool isAlgebraic(int k) {return (m_alg[k] == 1); } + + + /** + * Evaluate the residual function. Called by the + * integrator. + * @param t time. (input) + * @param y solution vector. (input) + * @param ydot rate of change of solution vector. (input) + * @param r residual vector (output) + */ + virtual int eval(double t, const double* y, + const double* ydot, double* r)=0; + + /** + * Fill the solution and derivative vectors with the initial + * conditions at initial time t0. If these do not satisfy the + * residual equation, call one of the "corrrectInitial_xxx" + * methods before calling solve. + */ + virtual void getInitialConditions(double t0, double* y, + doublereal* ydot)=0; + + /** + * Number of equations. + */ + virtual int nEquations()=0; + + + protected: + + map m_alg; + map m_constrain; + + private: + + }; + +} + +#endif