State.cpp

Go to the documentation of this file.
00001 /**
00002  *  @file State.cpp
00003  * Definitions for the class State, that manages the independent variables of temperature, mass density,
00004  * and species mass/mole fraction that define the thermodynamic state (see \ref phases and
00005  * class \link Cantera::State State\endlink).
00006  */
00007 
00008 /*
00009  *  
00010  *  $Date: 2010-02-09 15:24:11 -0500 (Tue, 09 Feb 2010) $
00011  *  $Revision: 398 $
00012  *
00013  *  Copyright 2003-2004 California Institute of Technology
00014  *  See file License.txt for licensing information
00015  *
00016  */
00017 
00018 #include "utilities.h"
00019 #include "ctexceptions.h"
00020 #include "stringUtils.h"
00021 #include "State.h"
00022 
00023 //#ifdef DARWIN
00024 //#include <Accelerate.h>
00025 //#endif
00026 
00027 using namespace std;
00028 
00029 namespace Cantera {
00030 
00031   inline void State::stateMFChangeCalc(bool forcerChange) {
00032     // Right now we assume that the mole fractions have changed every time
00033     // the function is called
00034     m_stateNum++;
00035     if (m_stateNum > 1000000) m_stateNum = -10000000;
00036   }
00037 
00038   State::State() : 
00039     m_kk(0),
00040     m_temp(0.0), 
00041     m_dens(0.001), 
00042     m_mmw(0.0),
00043     m_stateNum(-1)
00044   {
00045   }
00046 
00047   State::~State() 
00048   {
00049   }
00050 
00051   State::State(const State& right) : 
00052     m_kk(0),
00053     m_temp(0.0),
00054     m_dens(0.001), 
00055     m_mmw(0.0),
00056     m_stateNum(-1)
00057   {
00058     /*
00059      * Call the assignment operator.
00060      */
00061     *this = operator=(right);
00062   }
00063 
00064   /*
00065    * Assignment operator for the State Class
00066    */
00067   State& State::operator=(const State& right) {
00068     /*
00069      * Check for self assignment.
00070      */
00071     if (this == &right) return *this;
00072     /*
00073      * We do a straight assignment operator on all of the
00074      * data. The vectors are copied.
00075      */
00076     m_kk             = right.m_kk;
00077     m_temp           = right.m_temp;
00078     m_dens           = right.m_dens;
00079     m_mmw            = right.m_mmw;
00080     m_ym             = right.m_ym;
00081     m_y              = right.m_y;
00082     m_molwts         = right.m_molwts;
00083     m_rmolwts        = right.m_rmolwts;
00084     m_stateNum = -1;
00085     /*
00086      * Return the reference to the current object
00087      */
00088     return *this;
00089   }
00090 
00091   doublereal State::moleFraction(const int k) const {
00092     if (k >= 0 && k < m_kk) {
00093       return m_ym[k] * m_mmw;
00094     }
00095     else {
00096       throw CanteraError("State:moleFraction",
00097                          "illegal species index number");
00098     }
00099     return 0.0;
00100   }
00101 
00102   void State::setMoleFractions(const doublereal* const x) {
00103     doublereal sum = dot(x, x + m_kk, m_molwts.begin());
00104     doublereal rsum = 1.0/sum;
00105     transform(x, x + m_kk, m_ym.begin(), timesConstant<double>(rsum));
00106     transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), 
00107               m_y.begin(), multiplies<double>());
00108     doublereal norm = accumulate(x, x + m_kk, 0.0);
00109     m_mmw = sum/norm;
00110    
00111     //! Call a routine to determin whether state has changed.
00112     stateMFChangeCalc();
00113   }
00114 
00115   void State::setMoleFractions_NoNorm(const doublereal* const x) {
00116     m_mmw = dot(x, x + m_kk, m_molwts.begin());
00117     doublereal rmmw = 1.0/m_mmw;
00118     transform(x, x + m_kk, m_ym.begin(), timesConstant<double>(rmmw));
00119     transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), 
00120               m_y.begin(), multiplies<double>());
00121 
00122     //! Call a routine to determin whether state has changed.
00123     stateMFChangeCalc();
00124   }
00125 
00126   doublereal State::massFraction(const int k) const {
00127     if (k >= 0 && k < m_kk) {
00128       return m_y[k];
00129     }
00130     throw CanteraError("State:massFraction", "illegal species index number");
00131     return 0.0;
00132   }
00133 
00134   doublereal State::concentration(const int k) const {
00135     if (k >= 0 && k < m_kk) {
00136       return m_y[k] * m_dens * m_rmolwts[k] ;
00137     }
00138     throw CanteraError("State:massFraction", "illegal species index number");
00139     return 0.0;
00140   }
00141 
00142   void State::setMassFractions(const doublereal* const y) {
00143     doublereal norm = 0.0, sum = 0.0;
00144     //cblas_dcopy(m_kk, y, 1, m_y.begin(), 1);
00145     norm = accumulate(y, y + m_kk, 0.0);
00146     copy(y, y + m_kk, m_y.begin());
00147     scale(y, y + m_kk, m_y.begin(), 1.0/norm);
00148     //        for (k = 0; k != m_kk; ++k) {
00149     //    norm += y[k];
00150     //    m_y[k] = y[k];
00151     //}
00152 
00153     //scale(m_kk, 1.0/norm, m_y.begin());
00154     transform(m_y.begin(), m_y.begin() + m_kk, m_rmolwts.begin(),
00155               m_ym.begin(), multiplies<double>());
00156     sum = accumulate(m_ym.begin(), m_ym.begin() + m_kk, 0.0);
00157     //    for (k = 0; k != m_kk; ++k) {
00158     //      m_ym[k] = m_y[k] * m_rmolwts[k];
00159     //      sum += m_ym[k];
00160     //    }
00161     m_mmw = 1.0/sum;
00162 
00163     //! Call a routine to determin whether state has changed.
00164     stateMFChangeCalc();
00165   }
00166 
00167   void State::setMassFractions_NoNorm(const doublereal* const y) {
00168     doublereal sum = 0.0;
00169     copy(y, y + m_kk, m_y.begin());
00170     transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
00171               multiplies<double>());
00172     sum = accumulate(m_ym.begin(), m_ym.end(), 0.0);
00173     //for (k = 0; k != m_kk; ++k) {
00174     //    m_y[k] = y[k];
00175     //    m_ym[k] = m_y[k] * m_rmolwts[k];
00176     //    sum += m_ym[k];
00177     //}
00178     m_mmw = 1.0/sum;
00179 
00180     //! Call a routine to determine whether state has changed.
00181     stateMFChangeCalc();
00182   }
00183 
00184   doublereal State::sum_xlogx() const {
00185     return m_mmw* Cantera::sum_xlogx(m_ym.begin(), m_ym.end()) + log(m_mmw);
00186   }
00187 
00188   doublereal State::sum_xlogQ(doublereal* Q) const {
00189     return m_mmw * Cantera::sum_xlogQ(m_ym.begin(), m_ym.end(), Q);
00190   }
00191 
00192   doublereal State::molarDensity() const { 
00193     return density()/meanMolecularWeight(); 
00194   }
00195 
00196   void State::setConcentrations(const doublereal* const conc) {
00197     int k;
00198     doublereal sum = 0.0, norm = 0.0;
00199     for (k = 0; k != m_kk; ++k) {
00200       sum += conc[k]*m_molwts[k];
00201       norm += conc[k];
00202     }
00203     m_mmw = sum/norm;
00204     setDensity(sum);
00205     doublereal rsum = 1.0/sum;
00206     for (k = 0; k != m_kk; ++k) {
00207       m_ym[k] = conc[k] * rsum;
00208       m_y[k] =  m_ym[k] * m_molwts[k];
00209     }
00210 
00211     // Call a routine to determine whether state has changed.
00212     stateMFChangeCalc();
00213   }
00214 
00215   const doublereal* State::moleFractdivMMW() const { 
00216     return &m_ym[0];
00217   }
00218 
00219   void State::getConcentrations(doublereal* const c) const {
00220     scale(m_ym.begin(), m_ym.end(), c, m_dens);
00221   }
00222 
00223   doublereal State::mean_X(const doublereal* const Q) const {
00224     return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q, 0.0);
00225   }
00226 
00227   doublereal State::mean_Y(const doublereal* const Q) const {
00228     return dot(m_y.begin(), m_y.end(), Q);
00229   }
00230 
00231   void State::getMoleFractions(doublereal* const x) const {
00232     scale(m_ym.begin(), m_ym.end(), x, m_mmw);
00233   }
00234 
00235   void State::getMassFractions(doublereal* const y) const {
00236     copy(m_y.begin(), m_y.end(), y);
00237   }
00238 
00239   void State::setMolarDensity(const doublereal molarDensity) {
00240     m_dens = molarDensity*meanMolecularWeight();
00241   }
00242 
00243 
00244   void State::init(const array_fp& mw) {
00245     m_kk = mw.size();
00246     m_molwts.resize(m_kk);
00247     m_rmolwts.resize(m_kk);
00248     m_y.resize(m_kk, 0.0);
00249     m_ym.resize(m_kk, 0.0);
00250     copy(mw.begin(), mw.end(), m_molwts.begin());
00251     for (int k = 0; k < m_kk; k++) {
00252       if (m_molwts[k] < 0.0) {
00253         throw CanteraError("State::init",
00254                            "negative molecular weight for species number "+int2str(k));
00255       }
00256       /*
00257        * Some surface phases may define species representing
00258        * empty sites that have zero molecular weight. Give them
00259        * a very small molecular weight to avoid dividing by
00260        * zero.
00261        */
00262       if (m_molwts[k] < Tiny) m_molwts[k] = Tiny;
00263       m_rmolwts[k] = 1.0/m_molwts[k];
00264     }
00265 
00266     /*
00267      * Now that we have resized the State object, let's fill it with
00268      * a valid mass fraction vector that sums to one. The State object
00269      * should never have a mass fraction vector that doesn't sum to one.
00270      * We will assume that species 0 has a mass fraction of 1.0 and
00271      * mass fraction of all other species is 0.0.
00272      */
00273     m_y[0] = 1.0;
00274     m_ym[0] = m_y[0] * m_rmolwts[0];
00275     m_mmw = 1.0 / m_ym[0];
00276   }
00277 
00278   // True if the number of species has been set and fixed
00279   bool State::ready() const { 
00280     return (m_kk > 0); 
00281   }
00282 
00283 
00284 }
Generated by  doxygen 1.6.3