From 3b4597521d0740bed940dd2119a55fa448edb84a Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Thu, 16 Dec 2004 06:56:10 +0000 Subject: [PATCH] initial import --- Cantera/src/MultiPhase.cpp | 246 +++++++++++++++++++++++++++++++++++++ 1 file changed, 246 insertions(+) create mode 100644 Cantera/src/MultiPhase.cpp diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp new file mode 100644 index 000000000..4ceddd0a3 --- /dev/null +++ b/Cantera/src/MultiPhase.cpp @@ -0,0 +1,246 @@ +#include "MultiPhase.h" + +#include "ThermoPhase.h" +#include "DenseMatrix.h" +#include "stringUtils.h" + +namespace Cantera { + + void MultiPhase:: + addPhase(phase_t* p, doublereal moles) { + + if (m_init) { + throw CanteraError("addPhase", + "phases cannot be added after init() has been called."); + } + + // save the pointer to the phase object + m_phase.push_back(p); + + // store its number of moles + m_moles.push_back(moles); + + // update the number of phases and the total number of + // species + m_np = m_phase.size(); + m_nsp += p->nSpecies(); + + // determine if this phase has new elements + // for each new element, add an entry in the map + // from names to index number + 1: + + string ename; + // iterate over the elements in this phase + index_t m, nel = p->nElements(); + for (m = 0; m < nel; m++) { + ename = p->elementName(m); + + // if no entry is found for this element name, then + // it is a new element. In this case, add the name + // to the list of names, increment the element count, + // and add an entry to the name->(index+1) map. + if (m_enamemap[ename] == 0) { + m_enamemap[ename] = m_nel + 1; + m_enames.push_back(ename); + m_nel++; + } + } + + if (m_temp == 0.0 && p->temperature() > 0.0) { + m_temp = p->temperature(); + m_press = p->pressure(); + } + } + + + /// Process phases and build atomic composition array. After + /// init() has been called, no more phases may be added. + void MultiPhase::init() { + if (m_init) return; + index_t ip, kp, k = 0, nsp, m; + int mlocal; + string sym; + + // allocate space for the atomic composition matrix + m_atoms.resize(m_nel, m_nsp, 0.0); + m_moleFractions.resize(m_nsp, 0.0); + + // iterate over the elements + for (m = 0; m < m_nel; m++) { + sym = m_enames[m]; + k = 0; + // iterate over the phases + for (ip = 0; ip < m_np; ip++) { + phase_t* p = m_phase[ip]; + nsp = p->nSpecies(); + mlocal = p->elementIndex(sym); + for (kp = 0; kp < nsp; kp++) { + if (mlocal >= 0) { + m_atoms(m, k) = p->nAtoms(kp, mlocal); + } + if (m == 0) { + m_snames.push_back(p->speciesName(kp)); + if (kp == 0) { + m_spstart.push_back(m_spphase.size()); + } + m_spphase.push_back(ip); + } + k++; + } + } + } + + /// set the initial composition within each phase to the + /// mole fractions stored in the phase objects + m_init = true; + updateMoleFractions(); + } + + + /// Return a reference to phase n. The state of phase n is + /// also updated to match the state stored locally in the + /// mixture object. + MultiPhase::phase_t& MultiPhase::phase(index_t n) { + if (!m_init) init(); + m_phase[n]->setState_TPX(m_temp, m_press, + m_moleFractions.begin() + m_spstart[n]); + return *m_phase[n]; + } + + /// Moles of species \c k. + doublereal MultiPhase::speciesMoles(index_t k) { + if (!m_init) init(); + index_t ip = m_spphase[k]; + return m_moles[ip]*m_moleFractions[k]; + } + + + /// Total moles of element m, summed over all + /// phases + doublereal MultiPhase::elementMoles(index_t m) { + doublereal sum = 0.0, phasesum; + index_t i, k = 0, ik, nsp; + for (i = 0; i < m_np; i++) { + phasesum = 0.0; + nsp = m_phase[i]->nSpecies(); + for (ik = 0; ik < nsp; ik++) { + k = speciesIndex(ik, i); + phasesum += m_atoms(m,k)*m_moleFractions[k]; + } + sum += phasesum * m_moles[i]; + } + return sum; + } + + /// Chemical potentials. Write into array \c mu the chemical + /// potentials of all species [J/kmol]. + void MultiPhase::getChemPotentials(doublereal* mu) { + index_t i, loc = 0; + updatePhases(); + for (i = 0; i < m_np; i++) { + m_phase[i]->getChemPotentials(mu + loc); + loc += m_phase[i]->nSpecies(); + } + } + + /// Chemical potentials. Write into array \c mu the chemical + /// potentials of all species [J/kmol]. + void MultiPhase::getStandardChemPotentials(doublereal* mu) { + index_t i, loc = 0; + updatePhases(); + for (i = 0; i < m_np; i++) { + m_phase[i]->getStandardChemPotentials(mu + loc); + loc += m_phase[i]->nSpecies(); + } + } + + bool MultiPhase::solutionSpecies(index_t k) { + if (m_phase[m_spphase[k]]->nSpecies() > 1) + return true; + else + return false; + } + + doublereal MultiPhase::gibbs() { + index_t i; + doublereal sum = 0.0; + updatePhases(); + for (i = 0; i < m_np; i++) + sum += m_phase[i]->gibbs_mole() * m_moles[i]; + return sum; + } + + void MultiPhase::updateMoleFractions() { + if (!m_init) init(); + // save the current mole fractions for each phase + index_t ip, loc = 0; + for (ip = 0; ip < m_np; ip++) { + phase_t* p = m_phase[ip]; + p->getMoleFractions(m_moleFractions.begin() + loc); + loc += p->nSpecies(); + } + } + + void MultiPhase::setPhaseMoleFractions(index_t n, doublereal* x) { + phase_t* p = m_phase[n]; + p->setState_TPX(m_temp, m_press, x); + } + + void MultiPhase::setMolesByName(compositionMap& xMap) { + int kk = nSpecies(); + doublereal x; + vector_fp mf(kk, 0.0); + for (int k = 0; k < kk; k++) { + x = xMap[speciesName(k)]; + if (x > 0.0) mf[k] = x; + } + setMoles(mf.begin()); + } + + void MultiPhase::setMolesByName(const string& x) { + compositionMap xx; + int kk = nSpecies(); + for (int k = 0; k < kk; k++) { + xx[speciesName(k)] = -1.0; + } + parseCompString(x, xx); + setMolesByName(xx); + } + + void MultiPhase::setMoles(doublereal* n) { + if (!m_init) init(); + index_t ip, loc = 0; + index_t ik, k = 0, nsp; + doublereal phasemoles; + for (ip = 0; ip < m_np; ip++) { + phase_t* p = m_phase[ip]; + nsp = p->nSpecies(); + phasemoles = 0.0; + for (ik = 0; ik < nsp; ik++) { + phasemoles += n[k]; + k++; + } + m_moles[ip] = phasemoles; + if (nsp > 1) { + p->setState_TPX(m_temp, m_press, n + loc); + p->getMoleFractions(m_moleFractions.begin() + loc); + } + else { + m_moleFractions[loc] = 1.0; + } + loc += p->nSpecies(); + } + } + + void MultiPhase::updatePhases() { + if (!m_init) init(); + index_t p, nsp, loc = 0; + for (p = 0; p < m_np; p++) { + nsp = m_phase[p]->nSpecies(); + doublereal* x = m_moleFractions.begin() + loc; + loc += nsp; + m_phase[p]->setState_TPX(m_temp, m_press, x); + } + } +} +