diff --git a/Cantera/python/Cantera/ReactorNet.py b/Cantera/python/Cantera/ReactorNet.py index f2b9c4394..9961ca0d5 100644 --- a/Cantera/python/Cantera/ReactorNet.py +++ b/Cantera/python/Cantera/ReactorNet.py @@ -59,6 +59,10 @@ class ReactorNet: using the current state as the initial condition. Default: 0.0 s""" _cantera.reactornet_setInitialTime(self.__reactornet_id, t0) + def setTolerances(self, rtol, atol): + """Set the relative and absolute error tolerances.""" + _cantera.reactornet_setTolerances(self.__reactornet_id, rtol, atol) + def advance(self, time): """Advance the state of the reactor network in time from the current time to time 'time'.""" diff --git a/Cantera/python/examples/equilibrium/multiphase_plasma.py b/Cantera/python/examples/equilibrium/multiphase_plasma.py index c1572fb99..61de96ceb 100644 --- a/Cantera/python/examples/equilibrium/multiphase_plasma.py +++ b/Cantera/python/examples/equilibrium/multiphase_plasma.py @@ -33,7 +33,7 @@ for n in range(500): # set the mixture to a state of chemical equilibrium holding # temperature and pressure fixed - mix.equilibrate("TP",maxsteps=1000,loglevel=0) + mix.equilibrate("TP",maxsteps=10000,loglevel=1) # write out the moles of each species writeCSV(f,[t]+ list(mix.speciesMoles())) diff --git a/Cantera/python/src/pycantera.cpp b/Cantera/python/src/pycantera.cpp index d8c7baea5..68a1177f0 100644 --- a/Cantera/python/src/pycantera.cpp +++ b/Cantera/python/src/pycantera.cpp @@ -20,6 +20,8 @@ #else //#ifdef HAS_NUMARRAY #include "numarray/arrayobject.h" +//#include "numpy/libnumarray.h" + //#else //#include "numpy/arrayobject.h" //#endif diff --git a/Cantera/src/base/ct2ctml.cpp b/Cantera/src/base/ct2ctml.cpp index a2a6f1208..c3e2b92b0 100644 --- a/Cantera/src/base/ct2ctml.cpp +++ b/Cantera/src/base/ct2ctml.cpp @@ -24,7 +24,7 @@ using namespace std; -//#define DEBUG_PATHS +#undef DEBUG_PATHS using namespace Cantera; @@ -34,14 +34,26 @@ namespace ctml { // environment variable PYTHON_CMD if it is set. If not, return // the string 'python'. static string pypath() { - string s = "python"; + string s = "ctpython"; const char* py = getenv("PYTHON_CMD"); + if (!py) { + const char* hm = getenv("HOME"); + string home = stripws(string(hm)); + string cmd = string("source ")+home + +string("/setup_cantera &> /dev/null"); + system(cmd.c_str()); + py = getenv("PYTHON_CMD"); + } if (py) { string sp = stripws(string(py)); if (sp.size() > 0) { s = sp; } } + else { + throw CanteraError("ct2ctml", + "set environment variable PYTHON_CMD"); + } return s; } @@ -58,6 +70,7 @@ namespace ctml { "python cti to ctml conversion requested for file, " + ppath + ", but not available in this computational environment"); #endif + time_t aclock; time( &aclock ); int ia = static_cast(aclock); diff --git a/Cantera/src/thermo/SemiconductorPhase.cpp b/Cantera/src/thermo/SemiconductorPhase.cpp new file mode 100644 index 000000000..9be205c89 --- /dev/null +++ b/Cantera/src/thermo/SemiconductorPhase.cpp @@ -0,0 +1,64 @@ + +#include "SemiconductorPhase.h" + +namespace Cantera { + + const doublereal JD_const1 = 1.0/sqrt(8.0); + const doublereal JD_const2 = 3.0/16.0 - sqrt(3.0)/9.0; + + static doublereal JoyceDixon(doublereal r) { + return log(r) + JD_const1*r - JD_const2*r*r; + } + + // doublereal SemiconductorPhase::ionizedDonorConcentration() { + // return 1.0/(1.0 + 2.0*exp( fermiLevel() - m_edonor)); + //} + + //doublereal SemiconductorPhase::ionizedAcceptorConcentration() { + // return 1.0/(1.0 + 2.0*exp( m_eacceptor - fermiLevel())); + //} + + //doublereal SemiconductorPhase::_dn(doublereal efermi) { + // m_fermi_level = efermi; + // return electronConcentration() - holeConcentration() + + // ionizedAcceptorConcentration() - ionizedDonorConcentration(); + //} + void SemiconductorPhase::getChemPotentials(doublereal* mu) const { + getActivityConcentrations(DATA_PTR(m_work)); + doublereal r = m_work[0]/nc(); + mu[0] = ec() + GasConstant*temperature()*(JoyceDixon(r)); + mu[1] = ev() + + GasConstant*temperature()*(log(m_work[1]/nv())); + } + + // units: kmol/m^3 + doublereal SemiconductorPhase::nc() const { + doublereal fctr = effectiveMass_e() * Boltzmann * temperature()/ + (2.0*Pi*Planck_bar*Planck_bar); + return 2.0*pow(fctr, 1.5)/Avogadro; + } + + doublereal SemiconductorPhase::nv() const { + doublereal fctr = effectiveMass_h() * Boltzmann * temperature()/ + (2.0*Pi*Planck_bar*Planck_bar); + return 2.0*pow(fctr, 1.5)/Avogadro; + } + + doublereal SemiconductorPhase::ev() const { + return 0.0; + } + + doublereal SemiconductorPhase::ec() const { + return ev() + bandgap(); + } + + doublereal SemiconductorPhase::bandgap() const { + return m_gap; + } + + + // private + void SemiconductorPhase::initLengths() { + int ns = nSpecies(); + m_work.resize(ns); + } +} diff --git a/Cantera/src/thermo/SemiconductorPhase.h b/Cantera/src/thermo/SemiconductorPhase.h new file mode 100644 index 000000000..7cca66de4 --- /dev/null +++ b/Cantera/src/thermo/SemiconductorPhase.h @@ -0,0 +1,1121 @@ +/** + * @file SemiconductorPhase.h + * Header file for an ideal semiconductor model + * + * This class inherits from the Cantera class ThermoPhase + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + */ + +#ifndef CT_SEMICONDUCTOR_H +#define CT_SEMICONDUCTOR_H + +#include "mix_defs.h" +#include "ThermoPhase.h" +//#include "importCTML.h" +#include "ThermoFactory.h" +#include "SpeciesThermo.h" + + +namespace Cantera { + + /*! + * @name CONSTANTS - Models for the Standard State of IdealSolidSolnPhase's + */ + //@{ + const int cSemiconductorPhase = 6010; + //@} + + /** + * + * @ingroup thermoprops + */ + class SemiconductorPhase : public ThermoPhase { + + public: + + /** + * Constructor for SemiconductorPhase. + * The generalized concentrations can have three different forms + * depending on the value of the member attribute m_formGC, which + * is supplied in the constructor or read from the xml data file. + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * @param formCG This parameter initializes the m_formGC variable. The default + * is a value of 0. + */ + SemiconductorPhase(int formCG=0); + + /** + * Constructor for SemiconductorPhase. + * + * This constructor will also fully initialize the object. + * The generalized concentrations can have three different forms + * depending on the value of the member attribute m_formGC, which + * is supplied in the constructor or read from the xml data file. + * + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * @param infile File name for the XML datafile containing information + * for this phase + * @param id The name of this phase. This is used to look up + * the phase in the XML datafile. + * @param formCG This parameter initializes the m_formGC variable. The default + * is a value of 0. + */ + SemiconductorPhase(std::string infile, std::string id=""); + + + /** + * Constructor for SemiconductorPhase. + * This constructor will also fully initialize the object. + * + * The generalized concentrations can have three different forms + * depending on the value of the member attribute m_formGC, which + * is supplied in the constructor and/or read from the data file. + * + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * @param root XML tree containing a description of the phase. + * The tree must be positioned at the XML element + * named phase with id, "id", on input to this routine. + * @param id The name of this phase. This is used to look up + * the phase in the XML datafile. + * @param formCG This parameter initializes the m_formGC variable. The default + * is a value of 0. + */ + SemiconductorPhase(XML_Node& root, std::string id=""); + + /*! + * Copy Constructor + */ + SemiconductorPhase(const SemiconductorPhase &); + + /*! + * Assignment operator + */ + SemiconductorPhase& operator=(const SemiconductorPhase &); + + /*! + * Base Class Duplication Function + * -> given a pointer to ThermoPhase, this function can + * duplicate the object. (note has to be a separate function + * not the copy constructor, because it has to be + * a virtual function) + */ + virtual ThermoPhase* duplMyselfAsThermoPhase() const; + + //! Destructor + virtual ~SemiconductorPhase() {} + + /** + * Equation of state flag. Returns a value depending upon the value of + * m_formGC, which is defined at instantiation. + */ + virtual int eosType() const; + + /** + * @name Molar Thermodynamic Properties of the Solution ------------------------ + * @{ + */ + + /** + * Molar enthalpy of the solution. Units: J/kmol. + * For an ideal, constant partial molar volume solution mixture with + * pure species phases which exhibit zero volume expansivity and + * zero isothermal compressibility: + * \f[ + * \hat h(T,P) = \sum_k X_k \hat h^0_k(T) + (P - P_{ref}) (\sum_k X_k \hat V^0_k) + * \f] + * The reference-state pure-species enthalpies at the reference pressure Pref + * \f$ \hat h^0_k(T) \f$, are computed by the species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + */ + virtual doublereal enthalpy_mole() const; + + /** + * Molar internal energy of the solution. Units: J/kmol. + * For an ideal, constant partial molar volume solution mixture with + * pure species phases which exhibit zero volume expansivity and + * zero isothermal compressibility: + * \f[ + * \hat u(T,X) = \hat h(T,P,X) - p \hat V + * = \sum_k X_k \hat h^0_k(T) - P_{ref} (\sum_k{X_k \hat V^0_k}) + * \f] + * and is a function only of temperature. + * The reference-state pure-species enthalpies + * \f$ \hat h^0_k(T) \f$ are computed by the species thermodynamic + * property manager. + * @see SpeciesThermo + */ + virtual doublereal intEnergy_mole() const; + + /** + * Molar entropy of the solution. Units: J/kmol/K. + * For an ideal, constant partial molar volume solution mixture with + * pure species phases which exhibit zero volume expansivity: + * \f[ + * \hat s(T, P, X_k) = \sum_k X_k \hat s^0_k(T) - \hat R \sum_k X_k log(X_k) + * \f] + * The reference-state pure-species entropies + * \f$ \hat s^0_k(T,p_{ref}) \f$ are computed by the species thermodynamic + * property manager. The pure species entropies are independent of + * temperature since the volume expansivities are equal to zero. + * @see SpeciesThermo + */ + virtual doublereal entropy_mole() const; + + /** + * Molar gibbs free energy of the solution. Units: J/kmol. + * For an ideal, constant partial molar volume solution mixture with + * pure species phases which exhibit zero volume expansivity: + * \f[ + * \hat g(T, P) = \sum_k X_k \hat g^0_k(T,P) + \hat R T \sum_k X_k log(X_k) + * \f] + * The reference-state pure-species gibbs free energies + * \f$ \hat g^0_k(T) \f$ are computed by the species thermodynamic + * property manager, while the standard state gibbs free energies + * \f$ \hat g^0_k(T,P) \f$ are computed by the member function, gibbs_RT(). + * @see SpeciesThermo + */ + virtual doublereal gibbs_mole() const; + + /** + * Molar heat capacity at constant pressure of the solution. + * Units: J/kmol/K. + * For an ideal, constant partial molar volume solution mixture with + * pure species phases which exhibit zero volume expansivity: + * \f[ + * \hat c_p(T,P) = \sum_k X_k \hat c^0_{p,k}(T) . + * \f] + * The heat capacity is independent of pressure. + * The reference-state pure-species heat capacities + * \f$ \hat c^0_{p,k}(T) \f$ are computed by the species thermodynamic + * property manager. + * @see SpeciesThermo + */ + virtual doublereal cp_mole() const; + + /** + * Molar heat capacity at constant volume of the solution. + * Units: J/kmol/K. + * For an ideal, constant partial molar volume solution mixture with + * pure species phases which exhibit zero volume expansivity: + * \f[ \hat c_v(T,P) = \hat c_p(T,P) \f] + * The two heat capacities are equal. + */ + virtual doublereal cv_mole() const { + return cp_mole(); + } + + //@} + /** @name Mechanical Equation of State Properties ------------------------------------ + * + * In this equation of state implementation, the density is a + * function only of the mole fractions. Therefore, it can't be + * an independent variable. Instead, the pressure is used as the + * independent variable. Functions which try to set the thermodynamic + * state by calling setDensity() may cause an exception to be + * thrown. + */ + //@{ + + /** + * Pressure. Units: Pa. + * For this incompressible system, we return the internally storred + * independent value of the pressure. + */ + virtual doublereal pressure() const { + return m_Pcurrent; + } + + /** + * Set the pressure at constant temperature. Units: Pa. + * This method sets a constant within the object. + * The mass density is not a function of pressure. + * + * @param p Input Pressure (Pa) + */ + virtual void setPressure(doublereal p); + + /** + * Calculate the density of the mixture using the partial + * molar volumes and mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are + * the molecular weights, and \f$V_k\f$ are the pure species + * molar volumes. + * + * Note, the basis behind this formula is that in an ideal + * solution the partial molar volumes are equal to the pure + * species molar volumes. We have additionally specified + * in this class that the pure species molar volumes are + * independent of temperature and pressure. + * + * NOTE: This is a non-virtual function, which is not a + * member of the ThermoPhase base class. + */ + void calcDensity(); + + /** + * Overwritten setDensity() function is necessary because the + * density is not an indendent variable. + * + * This function will now throw an error condition + * + * @internal May have to adjust the strategy here to make + * the eos for these materials slightly compressible, in order + * to create a condition where the density is a function of + * the pressure. + * + * This function will now throw an error condition. + * + * NOTE: This is a virtual function that overwrites the State.h + * class + * + * @param rho Input density + */ + virtual void setDensity(doublereal rho); + + /** + * Overwritten setMolarDensity() function is necessary because the + * density is not an independent variable. + * + * This function will now throw an error condition. + * + * NOTE: This is virtual function that overwrites the State.h + * class + * + * @param rho Input Density + */ + virtual void setMolarDensity(doublereal rho); + + //! Set the mole fractions + /*! + * @param x Input vector of mole fractions. + * Length: m_kk. + */ + virtual void setMoleFractions(const doublereal *x); + + //! Set the mole fractions, but don't normalize them to one. + /*! + * @param x Input vector of mole fractions. + * Length: m_kk. + */ + virtual void setMoleFractions_NoNorm(const doublereal *x); + + //! Set the mass fractions, and normalize them to one. + /*! + * @param y Input vector of mass fractions. + * Length: m_kk. + */ + virtual void setMassFractions(const doublereal *y); + + //! Set the mass fractions, but don't normalize them to one + /*! + * @param y Input vector of mass fractions. + * Length: m_kk. + */ + virtual void setMassFractions_NoNorm(const doublereal *y); + + //! Set the concentration, + /*! + * @param c Input vector of concentrations. + * Length: m_kk. + */ + virtual void setConcentrations(const doublereal *c); + + + //@} + + /** + * @name Chemical Potentials and Activities ----------------------------------------- + * + * The activity of a species is + * related to the chemical potential (quasi-Fermi level) by + * \f[ + * \mu_k(T,P,X_k) = \mu_k^0(T,P) + * + \hat R T \log a_k. + * \f] + * The quantity \f$\mu_k^0(T,P)\f$ is + * the standard state chemical potential at unit activity, also known as "ec" for + * electrons and "ev" for holes. + * + * The activities are related to the generalized + * concentrations, \f$\tilde C_k\f$, and standard + * concentrations, \f$C^0_k\f$, by the following formula: + * + * \f[ + * a_k = \frac{\tilde C_k}{C^0_k} + * \f] + * The generalized concentrations are used in the kinetics classes + * to describe the rates of progress of reactions involving the + * species. Their formulation depends upons the specification + * of the rate constants for reaction, especially the units used + * in specifying the rate constants. The bridge between the + * thermodynamic equilibrium expressions that use a_k and the + * kinetics expressions which use the generalized concentrations + * is provided by the multiplicative factor of the + * standard concentrations. + * @{ + */ + + /** + * This method returns the array of generalized + * concentrations. The generalized concentrations are used + * in the evaluation of the rates of progress for reactions + * involving species in this phase. The generalized + * concentration dividied by the standard concentration is also + * equal to the activity of species. + * + * For this implentation the activity is defined to be the + * mole fraction of the species. The generalized concentration + * is defined to be equal to the mole fraction divided by + * the partial molar volume. The generalized concentrations + * for species in this phase therefore have units of + * kmol m-3. Rate constants must reflect this fact. + * + * On a general note, the following must be true. + * For an ideal solution, the generalized concentration must consist + * of the mole fraction multiplied by a constant. The constant may be + * fairly arbitrarily chosen, with differences adsorbed into the + * reaction rate expression. 1/V_N, 1/V_k, or 1 are equally good, + * as long as the standard concentration is adjusted accordingly. + * However, it must be a constant (and not the concentration, btw, + * which is a function of the mole fractions) in order for the + * ideal solution properties to hold at the same time having the + * standard concentration to be independent of the mole fractions. + * + * In this implementation the form of the generalized concentrations + * depend upon the member attribute, m_formGC: + * + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * HKM Note: We have absorbed the pressure dependence of the pures species + * state into the thermodynamics functions. Therefore the + * standard state on which the activities are based depend + * on both temperature and pressure. If we hadn't, it would have + * appeared in this function in a very awkwards exp[] format. + * + * @param c Pointer to array of doubles of length m_kk, which on exit + * will contain the generalized concentrations. + */ + virtual void getActivityConcentrations(doublereal* c) const; + + /** + * The standard concentration \f$ C^0_k \f$ used to normalize + * the generalized concentration. + * In many cases, this quantity + * will be the same for all species in a phase. + * However, for this case, we will return a distinct concentration + * for each species. This is the inverse of the species molar + * volume. Units for the standard concentration are + * kmol m-3. + * + * @param k Species number: this is a require parameter, + * a change from the ThermoPhase base class, where it was + * an optional parameter. + */ + virtual doublereal standardConcentration(int k) const; + + /** + * The reference (ie standard) concentration \f$ C^0_k \f$ used to normalize + * the generalized concentration. In many cases, this quantity + * will be the same for all species in a phase. + * However, for this case, we will return a distinct concentration + * for each species. (clone of the standard concentration -> + * suggest changing the name). This is the inverse of the species molar + * volume. + * + * @param k Species index. + */ + virtual doublereal referenceConcentration(int k) const; + + /** + * Returns the log of the standard concentration of the kth species + * + * @param k Species number: this is a require parameter, + * a change from the ThermoPhase base class, where it was + * an optional parameter. + */ + virtual doublereal logStandardConc(int k) const; + + /** + * Returns the units of the standard and general concentrations + * Note they have the same units, as their divisor is + * defined to be equal to the activity of the kth species + * in the solution, which is unitless. + * + * This routine is used in print out applications where the + * units are needed. Usually, MKS units are assumed throughout + * the program and in the XML input files. + * + * @param uA Output vector containing the units + * uA[0] = kmol units - default = 1 + * uA[1] = m units - default = -nDim(), the number of spatial + * dimensions in the Phase class. + * uA[2] = kg units - default = 0; + * uA[3] = Pa(pressure) units - default = 0; + * uA[4] = Temperature units - default = 0; + * uA[5] = time units - default = 0 + * @param k species index. Defaults to 0. + * @param sizeUA output int containing the size of the vector. + * Currently, this is equal to 6. + * + * For EOS types other than cSemiconductorPhase0, the default + * kmol/m3 holds for standard concentration units. For + * cSemiconductorPhase0 type, the standard concentrtion is + * unitless. + */ + virtual void getUnitsStandardConc(double *uA, int k = 0, + int sizeUA = 6) const; + + + //! Get the array of species activity coefficients + /*! + * @param ac output vector of activity coefficients. Length: m_kk + */ + virtual void getActivityCoefficients(doublereal * ac) const; + + /** + * Get the species chemical potentials. Units: J/kmol. + * + * This function returns a vector of chemical potentials of the + * species in solution. + * \f[ + * \mu_k = \mu^{ref}_k(T) + V_k * (p - p_o) + R T ln(X_k) + * \f] + * or another way to phrase this is + * \f[ + * \mu_k = \mu^o_k(T,p) + R T ln(X_k) + * \f] + * where \f$ \mu^o_k(T,p) = \mu^{ref}_k(T) + V_k * (p - p_o)\f$ + * + * @param mu Output vector of chemical potentials. + */ + virtual void getChemPotentials(doublereal* mu) const; + + /** + * Get the array of non-dimensional species solution + * chemical potentials at the current T and P + * \f$\mu_k / \hat R T \f$. + * \f[ + * \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k + RT ln(X_k) + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k at the reference pressure, \f$P_{ref}\f$. + * + * @param mu Output vector of dimensionless chemical potentials. Length = m_kk. + */ + virtual void getChemPotentials_RT(doublereal* mu) const; + + //@} + /// @name Partial Molar Properties of the Solution ----------------------------- + //@{ + + /** + * Returns an array of partial molar enthalpies for the species + * in the mixture. + * Units (J/kmol) + * For this phase, the partial molar enthalpies are equal to the + * pure species enthalpies + * \f[ + * \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k + * \f] + * The reference-state pure-species enthalpies, \f$ \hat h^{ref}_k(T) \f$, + * at the reference pressure,\f$ P_{ref} \f$, + * are computed by the species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + * + * @param hbar Output vector containing partial molar enthalpies. + * Length: m_kk. + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + /** + * Returns an array of partial molar entropies of the species in the + * solution. Units: J/kmol/K. + * For this phase, the partial molar entropies are equal to the + * pure species entropies plus the ideal solution contribution. + * \f[ + * \bar s_k(T,P) = \hat s^0_k(T) - R log(X_k) + * \f] + * The reference-state pure-species entropies,\f$ \hat s^{ref}_k(T) \f$, + * at the reference pressure, \f$ P_{ref} \f$, are computed by the + * species thermodynamic + * property manager. They are polynomial functions of temperature. + * @see SpeciesThermo + * + * @param sbar Output vector containing partial molar entropies. + * Length: m_kk. + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + /** + * Returns an array of partial molar Heat Capacities at constant + * pressure of the species in the + * solution. Units: J/kmol/K. + * For this phase, the partial molar heat capacities are equal + * to the standard state heat capacities. + * + * @param cpbar Output vector of partial heat capacities. Length: m_kk. + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + /** + * returns an array of partial molar volumes of the species + * in the solution. Units: m^3 kmol-1. + * + * For this solution, thepartial molar volumes are equal to the + * constant species molar volumes. + * + * @param vbar Output vector of partial molar volumes. Length: m_kk. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + //@} + /// @name Properties of the Standard State of the Species in the Solution ------------------------------------- + //@{ + + + /** + * Get the standard state chemical potentials of the species. + * This is the array of chemical potentials at unit activity + * \f$ \mu^0_k(T,P) \f$. + * We define these here as the chemical potentials of the pure + * species at the temperature and pressure of the solution. + * This function is used in the evaluation of the + * equilibrium constant Kc. Therefore, Kc will also depend + * on T and P. This is the norm for liquid and solid systems. + * + * units = J / kmol + * + * @param mu0 Output vector of standard state chemical potentials. + * Length: m_kk. + */ + virtual void getStandardChemPotentials(doublereal* mu0) const { + getPureGibbs(mu0); + } + + + //! Get the array of nondimensional Enthalpy functions for the standard state species + //! at the current T and P of the solution. + /*! + * We assume an incompressible constant partial molar + * volume here: + * \f[ + * h^0_k(T,P) = h^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ h^{ref}_k(T)\f$ is the enthalpy of the pure + * species k at the reference pressure, \f$P_{ref}\f$. + * + * @param hrt Vector of length m_kk, which on return hrt[k] + * will contain the nondimensional + * standard state enthalpy of species k. + */ + void getEnthalpy_RT(doublereal* hrt) const; + + + //! Get the nondimensional Entropies for the species + //! standard states at the current T and P of the solution. + /*! + * Note, this is equal to the reference state entropies + * due to the zero volume expansivity: + * i.e., (dS/dP)_T = (dV/dT)_P = 0.0 + * + * @param sr Vector of length m_kk, which on return sr[k] + * will contain the nondimensional + * standard state entropy for species k. + */ + void getEntropy_R(doublereal* sr) const; + + /** + * Get the nondimensional gibbs function for the species + * standard states at the current T and P of the solution. + * + * \f[ + * \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k at the reference pressure, \f$P_{ref}\f$. + * + * @param grt Vector of length m_kk, which on return sr[k] + * will contain the nondimensional + * standard state gibbs function for species k. + */ + virtual void getGibbs_RT(doublereal* grt) const; + + /** + * Get the Gibbs functions for the pure species + * at the current T and P of the solution. + * We assume an incompressible constant partial molar + * volume here: + * \f[ + * \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k at the reference pressure, \f$P_{ref}\f$. + * + * @param gpure Output vector of Gibbs functions for species + * Length: m_kk. + */ + virtual void getPureGibbs(doublereal* gpure) const; + + + //! Returns the vector of nondimensional + //! internal Energies of the standard state at the current + //! temperature and pressure of the solution for each species. + /*! + * + * @param urt Output vector of standard state nondimensional internal energies. + * Length: m_kk. + */ + virtual void getIntEnergy_RT(doublereal *urt) const; + + /** + * Get the nondimensional heat capacity at constant pressure + * function for the species + * standard states at the current T and P of the solution. + * \f[ + * Cp^0_k(T,P) = Cp^{ref}_k(T) + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ Cp^{ref}_k(T)\f$ is the constant pressure heat capacity + * of species k at the reference pressure, \f$p_{ref}\f$. + * + * @param cpr Vector of length m_kk, which on return cpr[k] + * will contain the nondimensional + * constant pressure heat capacity for species k. + */ + void getCp_R(doublereal* cpr) const; + + /** + * Get the molar volumes of each species in their standard + * states at the current + * T and P of the solution. + * units = m^3 / kmol + * + * @param vol Output vector of standard state volumes. + * Length: m_kk. + */ + virtual void getStandardVolumes(doublereal *vol) const; + + + //@} + /// @name Thermodynamic Values for the Species Reference States ------ + //@{ + + + /** + * Returns the vector of nondimensional + * enthalpies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * + * @param hrt Output vector containing reference nondimensional enthalpies. + * Length: m_kk. + */ + virtual void getEnthalpy_RT_ref(doublereal *hrt) const; + + /** + * Returns the vector of nondimensional + * enthalpies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * + * @param grt Output vector containing reference nondimensional Gibbs free energies. + * Length: m_kk. + */ + virtual void getGibbs_RT_ref(doublereal *grt) const; + + /** + * Returns the vector of the + * gibbs function of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * units = J/kmol + * + * @param g Output vector containing reference Gibbs free energies. + * Length: m_kk. + */ + virtual void getGibbs_ref(doublereal *g) const; + + /** + * Returns the vector of nondimensional + * entropies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + * + * @param er Output vector containing reference nondimensional entropies. + * Length: m_kk. + */ + virtual void getEntropy_R_ref(doublereal *er) const; + + /** + * Returns the vector of nondimensional + * internal Energies of the reference state at the current temperature + * of the solution and the reference pressure for each species. + * + * @param urt Output vector containing reference nondimensional internal energies. + * Length: m_kk. + */ + virtual void getIntEnergy_RT_ref(doublereal *urt) const; + + /** + * Returns the vector of nondimensional + * constant pressure heat capacities of the reference state + * at the current temperature of the solution + * and reference pressure for the species. + * + * @param cprt Output vector containing reference nondimensional heat capacities. + * Length: m_kk. + */ + virtual void getCp_R_ref(doublereal *cprt) const; + + /** + * Returns a reference to the vector of nondimensional + * enthalpies of the reference state at the current temperature. + * Real reason for its existence is that it also checks + * to see if a recalculation of the reference thermodynamics + * functions needs to be done. + */ + const array_fp& enthalpy_RT_ref() const; + + /** + * Returns a reference to the vector of nondimensional + * enthalpies of the reference state at the current temperature. + * Real reason for its existence is that it also checks + * to see if a recalculation of the reference thermodynamics + * functions needs to be done. + */ + const array_fp& gibbs_RT_ref() const { + _updateThermo(); + return m_g0_RT; + } + + /** + * Returns a reference to the vector of nondimensional + * enthalpies of the reference state at the current temperature. + * Real reason for its existence is that it also checks + * to see if a recalculation of the reference thermodynamics + * functions needs to be done. + */ + const array_fp& expGibbs_RT_ref() const; + + /** + * Returns a reference to the vector of nondimensional + * enthalpies of the reference state at the current temperature. + * Real reason for its existence is that it also checks + * to see if a recalculation of the reference thermodynamics + * functions needs to be done. + */ + const array_fp& entropy_R_ref() const; + + /** + * Returns a reference to the vector of nondimensional + * enthalpies of the reference state at the current temperature. + * Real reason for its existence is that it also checks + * to see if a recalculation of the reference thermodynamics + * functions needs to be done. + */ + const array_fp& cp_R_ref() const { + _updateThermo(); + return m_cp0_R; + } + + virtual void setPotentialEnergy(int k, doublereal pe) { + m_pe[k] = pe; + _updateThermo(); + } + + virtual doublereal potentialEnergy(int k) const { + return m_pe[k]; + } + //@} + /// @name Utility Functions ----------------------------------------------- + //@{ + + + + /** + * Initialization of an SemiconductorPhase phase using an + * xml file + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void constructPhaseFile(std::string infile, std::string id=""); + + /** + * Import and initialize an SemiconductorPhase phase + * specification in an XML tree into the current object. + * Here we read an XML description of the phase. + * We import descriptions of the elements that make up the + * species in a phase. + * We import information about the species, including their + * reference state thermodynamic polynomials. We then freeze + * the state of the species. + * This routine calls importPhase() to do most of its work. + * Then, importPhase() calls initThermoXML() to finish + * off the work. + * + + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + void constructPhaseXML(XML_Node& phaseNode, std::string id=""); + + /** + * Initialization of an SemiconductorPhase phase: + * Note this function is pretty much useless because it doesn't + * get the xml tree passed to it. Suggest a change. + */ + virtual void initThermo(); + + /** + * @internal + * Import and initialize a ThermoPhase object + * using an XML tree. + * Here we read extra information about the XML description + * of a phase. Regular information about elements and species + * and their reference state thermodynamic information + * have already been read at this point. + * For example, we do not need to call this function for + * ideal gas equations of state. + * This function is called from importPhase() + * after the elements and the + * species are initialized with default ideal solution + * level data. + * + * @param phaseNode This object must be the phase node of a + * complete XML tree + * description of the phase, including all of the + * species data. In other words while "phase" must + * point to an XML phase object, it must have + * sibling nodes "speciesData" that describe + * the species in the phase. + * @param id ID of the phase. If nonnull, a check is done + * to see if phaseNode is pointing to the phase + * with the correct id. + */ + virtual void initThermoXML(XML_Node& phaseNode, std::string id); + + + /** + * Set mixture to an equilibrium state consistent with specified + * element potentials and the temperature. + * + * @param lambda_RT vector of non-dimensional element potentials + * \f$ \lambda_m/RT \f$. + * + */ + virtual void setToEquilState(const doublereal* lambda_RT); + + + /** + * Report the molar volume of species k + * + * units - \f$ m^3 kmol^-1 \f$ + * + * @param k species index + */ + double speciesMolarVolume(int k) const; + + /** + * Fill in a return vector containing the species molar volumes. + * + * units - \f$ m^3 kmol^-1 \f$ + * + * @param smv output vector containing species molar volumes. + * Length: m_kk. + */ + void getSpeciesMolarVolumes(doublereal *smv) const; + + //@} + + doublereal nc() const; + doublereal nv() const; + doublereal ec() const; + doublereal ev() const; + doublereal bandgap() const; + doublereal effectiveMass_e() const { + return m_emass; + } + doublereal effectiveMass_h() const { + return m_hmass; + } + + protected: + + /** + * Format for the generalized concentrations + * 0 = C_k = X_k. (default) + * 1 = C_k = X_k / V_k + * 2 = C_k = X_k / V_N + */ + int m_formGC; + /** + * m_mm = Number of distinct elements defined in species in this + * phase + */ + int m_mm; + + /** + * Maximum temperature that this phase can accurately describe + * the thermodynamics. + */ + doublereal m_tmin; + + /** + * Minimum temperature that this phase can accurately describe + * the thermodynamics. + */ + doublereal m_tmax; + /** + * Value of the reference pressure for all species in this phase. + * The T dependent polynomials are evaluated at the reference + * pressure. Note, because this is a single value, all species + * are required to have the same reference pressure. + */ + doublereal m_Pref; + + /** + * m_Pcurrent = The current pressure + * Since the density isn't a function of pressure, but only of the + * mole fractions, we need to independently specify the pressure. + * The density variable which is inherited as part of the State class, + * m_dens, is always kept current whenever T, P, or X[] change. + */ + doublereal m_Pcurrent; + + /** + * Species molar volume \f$ m^3 kmol^-1 \f$ + */ + array_fp m_speciesMolarVolume; + + /** + * Value of the temperature at which the thermodynamics functions + * for the reference state of the species were last evaluated. + */ + mutable doublereal m_tlast; + + /** + * Vector containing the species reference enthalpies at T = m_tlast + */ + mutable array_fp m_h0_RT; + + /** + * Vector containing the species reference constant pressure + * heat capacities at T = m_tlast + */ + mutable array_fp m_cp0_R; + + /** + * Vector containing the species reference Gibbs functions + * at T = m_tlast + */ + mutable array_fp m_g0_RT; + + /** + * Vector containing the species reference entropies + * at T = m_tlast + */ + mutable array_fp m_s0_R; + + /** + * Vector containing the species reference exp(-G/RT) functions + * at T = m_tlast + */ + mutable array_fp m_expg0_RT; + + /** + * Vector of potential energies for the species. + */ + mutable array_fp m_pe; + + /** + * Temporary array used in equilibrium calculations + */ + mutable array_fp m_pp; + + mutable vector_fp m_work; + + doublereal m_gap; + doublereal m_emass; + doublereal m_hmass; + + private: + /// @name Utility Functions ------------------------------------------ + //@{ + /** + * This function gets called for every call to functions in this + * class. It checks to see whether the temperature has changed and + * thus the reference thermodynamics functions for all of the species + * must be recalculated. + * If the temperature has changed, the species thermo manager is called + * to recalculate G, Cp, H, and S at the current temperature. + */ + void _updateThermo() const; + + /** + * This internal function adjusts the lengths of arrays + */ + void initLengths(); + + //@} + }; +} + +#endif + + + + + diff --git a/Cantera/src/thermo/mix_defs.h b/Cantera/src/thermo/mix_defs.h index 7838e61db..72c6cb8fc 100755 --- a/Cantera/src/thermo/mix_defs.h +++ b/Cantera/src/thermo/mix_defs.h @@ -40,6 +40,7 @@ namespace Cantera { const int cMetal = 4; // MetalPhase in MetalPhase.h // const int cSolidCompound = 5; // SolidCompound in SolidCompound.h const int cStoichSubstance = 5; // StoichSubstance.h + const int cSemiconductor = 7; const int cLatticeSolid = 20; // LatticeSolidPhase.h const int cLattice = 21; @@ -69,3 +70,4 @@ namespace Cantera { } #endif + diff --git a/configure b/configure index bf490ca2d..b8a3df89f 100755 --- a/configure +++ b/configure @@ -2505,6 +2505,15 @@ _ACEOF hdrs=$hdrs' MetalPhase.h' fi + +if test "$WITH_SEMICONDUCTOR" = "y"; then + cat >>confdefs.h <<\_ACEOF +#define WITH_SEMICONDUCTOR 1 +_ACEOF + + hdrs=$hdrs' SemiconductorPhase.h' + objs=$objs' SemiconductorPhase.o' +fi if test "$WITH_STOICH_SUBSTANCE" = "y"; then cat >>confdefs.h <<\_ACEOF #define WITH_STOICH_SUBSTANCE 1 @@ -8476,7 +8485,7 @@ fi # Provide some information about the compiler. -echo "$as_me:8479:" \ +echo "$as_me:8488:" \ "checking for Fortran 77 compiler version" >&5 ac_compiler=`set X $ac_compile; echo $2` { (eval echo "$as_me:$LINENO: \"$ac_compiler --version &5\"") >&5 @@ -8679,7 +8688,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_verb" -(eval echo $as_me:8682: \"$ac_link\") >&5 +(eval echo $as_me:8691: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS @@ -8757,7 +8766,7 @@ _ACEOF # flags. ac_save_FFLAGS=$FFLAGS FFLAGS="$FFLAGS $ac_cv_prog_f77_v" -(eval echo $as_me:8760: \"$ac_link\") >&5 +(eval echo $as_me:8769: \"$ac_link\") >&5 ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'` echo "$ac_f77_v_output" >&5 FFLAGS=$ac_save_FFLAGS diff --git a/configure.in b/configure.in index 010134fc8..a3092dcde 100755 --- a/configure.in +++ b/configure.in @@ -430,6 +430,12 @@ if test "$WITH_METAL" = "y"; then AC_DEFINE(WITH_METAL) hdrs=$hdrs' MetalPhase.h' fi + +if test "$WITH_SEMICONDUCTOR" = "y"; then + AC_DEFINE(WITH_SEMICONDUCTOR) + hdrs=$hdrs' SemiconductorPhase.h' + objs=$objs' SemiconductorPhase.o' +fi if test "$WITH_STOICH_SUBSTANCE" = "y"; then AC_DEFINE(WITH_STOICH_SUBSTANCE) hdrs=$hdrs' StoichSubstance.h' diff --git a/preconfig b/preconfig index 937c02025..8b6cb1352 100755 --- a/preconfig +++ b/preconfig @@ -61,7 +61,7 @@ CANTERA_CONFIG_PREFIX=${CANTERA_CONFIG_PREFIX:=""} # default try to do a full installation, but fall back to a minimal # one in case of errors -PYTHON_PACKAGE=${PYTHON_PACKAGE:="default"} +PYTHON_PACKAGE=${PYTHON_PACKAGE:="full"} # Cantera needs to know where to find the Python interpreter. If # PYTHON_CMD is set to "default", then the configuration process will @@ -179,6 +179,7 @@ DEBUG_MODE="n" # thermodynamic properties ENABLE_THERMO='y' + ###################################################################### # optional phase types. These may not be needed by all users. Set them # to 'n' to omit them from the kernel. @@ -186,6 +187,8 @@ ENABLE_THERMO='y' WITH_LATTICE_SOLID=${WITH_LATTICE_SOLID:="y"} WITH_METAL=${WITH_METAL:="y"} WITH_STOICH_SUBSTANCE='y' +WITH_SEMICONDUCTOR='y' + # This flag enables the inclusion of accurate liquid/vapor equations # of state for several fluids, including water, nitrogen, hydrogen, @@ -523,6 +526,7 @@ export SUNDIALS_VERSION export WITH_LATTICE_SOLID export WITH_METAL +export WITH_SEMICONDUCTOR export WITH_STOICH_SUBSTANCE export WITH_PURE_FLUIDS export WITH_IDEAL_SOLUTIONS