initial import

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<LatticePhase*> 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

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

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 <vector>
+
+#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
+

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<int, int> m_alg;
+        map<int, int> m_constrain;
+
+    private:
+
+    };
+
+}
+
+#endif