19 MultiPhaseEquil::MultiPhaseEquil(
MultiPhase* mix,
bool start,
int loglevel) : m_mix(mix)
31 m_incl_species.resize(m_nsp_mix,1);
32 m_incl_element.resize(m_nel_mix,1);
33 for (
size_t m = 0; m < m_nel_mix; m++) {
37 if (enm ==
"E" || enm ==
"e") {
45 m_incl_element[m] = 0;
46 for (
size_t k = 0; k < m_nsp_mix; k++) {
47 if (m_mix->
nAtoms(k,m) != 0.0) {
48 m_incl_species[k] = 0;
56 if (m_eloc < m_nel_mix) {
57 m_element.push_back(m_eloc);
61 for (
size_t m = 0; m < m_nel_mix; m++) {
62 if (m_incl_element[m] == 1 && m != m_eloc) {
64 m_element.push_back(m);
77 for (
size_t k = 0; k < m_nsp_mix; k++) {
81 m_incl_species[k] = 0;
85 +
" is excluded since its thermo properties are \n" 86 "not valid at this temperature, but it has " 87 "non-zero moles in the initial state.");
93 for (
size_t k = 0; k < m_nsp_mix; k++) {
94 if (m_incl_species[k] ==1) {
96 m_species.push_back(k);
101 m_work.resize(m_nsp);
102 m_work2.resize(m_nsp);
103 m_work3.resize(m_nsp_mix);
104 m_mu.resize(m_nsp_mix);
107 m_moles.resize(m_nsp);
108 m_lastmoles.resize(m_nsp);
109 m_dxi.resize(
nFree());
112 for (
size_t ik = 0; ik < m_nsp; ik++) {
117 m_deltaG_RT.resize(
nFree(), 0.0);
118 m_majorsp.resize(m_nsp);
119 m_sortindex.resize(m_nsp,0);
120 m_lastsort.resize(m_nel);
121 m_solnrxn.resize(
nFree());
122 m_A.
resize(m_nel, m_nsp, 0.0);
124 m_order.resize(std::max(m_nsp, m_nel), 0);
125 iota(m_order.begin(), m_order.begin() + m_nsp, 0);
143 for (
size_t k = 0; k < m_nsp; k++) {
144 m_moles[k] += m_work[k];
145 m_lastmoles[k] = m_moles[k];
147 m_dsoln.push_back(1);
149 m_dsoln.push_back(0);
160 doublereal MultiPhaseEquil::equilibrate(
int XY, doublereal err,
161 int maxsteps,
int loglevel)
165 for (i = 0; i < maxsteps; i++) {
172 throw CanteraError(
"MultiPhaseEquil::equilibrate",
173 "no convergence in {} iterations. Error = {}",
180 void MultiPhaseEquil::updateMixMoles()
182 fill(m_work3.begin(), m_work3.end(), 0.0);
183 for (
size_t k = 0; k < m_nsp; k++) {
184 m_work3[m_species[k]] = m_moles[k];
191 fill(m_work3.begin(), m_work3.end(), 0.0);
192 for (
size_t k = 0; k < m_nsp; k++) {
193 m_work3[m_species[k]] = (m_moles[k] > 0.0 ? m_moles[k] : 0.0);
200 double not_mu = 1.0e12;
202 double dxi_min = 1.0e10;
216 for (
size_t j = 0; j <
nFree(); j++) {
219 for (
size_t ik = 0; ik < m_nsp; ik++) {
220 dg_rt += mu(ik) * m_N(ik,j);
224 int idir = (dg_rt < 0.0 ? 1 : -1);
226 for (
size_t ik = 0; ik < m_nsp; ik++) {
227 double nu = m_N(ik, j);
235 double delta_xi = fabs(0.99*moles(ik)/nu);
238 if (!redo && delta_xi < 1.0e-10 && ik < m_nel) {
241 dxi_min = std::min(dxi_min, delta_xi);
245 for (
size_t ik = 0; ik < m_nsp; ik++) {
246 moles(ik) += m_N(ik, j) * idir*dxi_min;
259 if (order.size() != m_nsp) {
260 for (
size_t k = 0; k < m_nsp; k++) {
264 for (
size_t k = 0; k < m_nsp; k++) {
265 m_order[k] = order[k];
269 size_t nRows = m_nel;
270 size_t nColumns = m_nsp;
273 for (
size_t m = 0; m < nRows; m++) {
274 for (
size_t k = 0; k < nColumns; k++) {
275 m_A(m, k) = m_mix->
nAtoms(m_species[m_order[k]], m_element[m]);
280 for (
size_t m = 0; m < nRows; m++) {
282 bool isZeroRow =
true;
283 for (
size_t k = m; k < nColumns; k++) {
284 if (fabs(m_A(m,k)) > sqrt(
Tiny)) {
291 size_t n = nRows - 1;
292 bool foundSwapCandidate =
false;
294 for (
size_t k = m; k < nColumns; k++) {
295 if (fabs(m_A(n,k)) > sqrt(
Tiny)) {
296 foundSwapCandidate =
true;
300 if (foundSwapCandidate) {
306 for (
size_t k = 0; k < nColumns; k++) {
307 std::swap(m_A(n,k), m_A(m,k));
318 if (m < nColumns && m_A(m,m) == 0.0) {
325 doublereal maxmoles = -999.0;
327 for (
size_t k = m+1; k < nColumns; k++) {
328 if (m_A(m,k) != 0.0 && fabs(m_moles[m_order[k]]) > maxmoles) {
330 maxmoles = fabs(m_moles[m_order[k]]);
336 for (
size_t n = 0; n < nRows; n++) {
337 std::swap(m_A(n, m), m_A(n, kmax));
341 std::swap(m_order[m], m_order[kmax]);
345 double fctr = 1.0/m_A(m,m);
346 for (
size_t k = 0; k < nColumns; k++) {
352 for (
size_t n = m+1; n < m_nel; n++) {
353 fctr = m_A(n,m)/m_A(m,m);
354 for (
size_t k = 0; k < m_nsp; k++) {
355 m_A(n,k) -= m_A(m,k)*fctr;
362 for (
size_t m = std::min(nRows,nColumns)-1; m > 0; m--) {
363 for (
size_t n = m-1; n !=
npos; n--) {
364 if (m_A(n,m) != 0.0) {
365 double fctr = m_A(n,m);
366 for (
size_t k = m; k < m_nsp; k++) {
367 m_A(n,k) -= fctr*m_A(m,k);
374 for (
size_t n = 0; n < m_nsp; n++) {
376 for (
size_t k = 0; k <
nFree(); k++) {
377 m_N(n, k) = -m_A(n, k + m_nel);
380 for (
size_t k = 0; k <
nFree(); k++) {
383 m_N(n, n - m_nel) = 1.0;
388 for (
size_t j = 0; j <
nFree(); j++) {
389 m_solnrxn[j] =
false;
390 for (
size_t k = 0; k < m_nsp; k++) {
401 for (
size_t k = 0; k < m_nsp; k++) {
402 x[m_order[k]] = m_work2[k];
406 void MultiPhaseEquil::step(doublereal omega,
vector_fp& deltaN,
410 throw CanteraError(
"MultiPhaseEquil::step",
"negative omega");
413 for (
size_t ik = 0; ik < m_nel; ik++) {
414 size_t k = m_order[ik];
415 m_lastmoles[k] = m_moles[k];
416 m_moles[k] += omega * deltaN[k];
419 for (
size_t ik = m_nel; ik < m_nsp; ik++) {
420 size_t k = m_order[ik];
421 m_lastmoles[k] = m_moles[k];
423 m_moles[k] += omega * deltaN[k];
425 m_moles[k] = fabs(m_moles[k])*std::min(10.0,
426 exp(-m_deltaG_RT[ik - m_nel]));
439 multiply(m_N, m_dxi.data(), m_work.data());
446 doublereal FCTR = 0.99;
447 const doublereal MAJOR_THRESHOLD = 1.0e-12;
448 double omegamax = 1.0;
449 for (
size_t ik = 0; ik < m_nsp; ik++) {
450 size_t k = m_order[ik];
453 if (m_moles[k] < MAJOR_THRESHOLD) {
462 if (m_dsoln[k] == 1) {
463 if ((m_moles[k] > MAJOR_THRESHOLD) || (ik < m_nel)) {
464 if (m_moles[k] < MAJOR_THRESHOLD) {
467 double omax = m_moles[k]*FCTR/(fabs(m_work[k]) +
Tiny);
468 if (m_work[k] < 0.0 && omax < omegamax) {
470 if (omegamax < 1.0e-5) {
476 m_majorsp[k] =
false;
479 if (m_work[k] < 0.0 && m_moles[k] > 0.0) {
480 double omax = -m_moles[k]/m_work[k];
481 if (omax < omegamax) {
483 if (omegamax < 1.0e-5) {
493 step(omegamax, m_work);
497 doublereal not_mu = 1.0e12;
499 doublereal grad1 = 0.0;
500 for (
size_t k = 0; k < m_nsp; k++) {
501 grad1 += m_work[k] * m_mu[m_species[k]];
504 double omega = omegamax;
506 omega *= fabs(grad0) / (grad1 + fabs(grad0));
507 for (
size_t k = 0; k < m_nsp; k++) {
508 m_moles[k] = m_lastmoles[k];
518 doublereal grad = 0.0;
521 doublereal not_mu = 1.0e12;
524 for (
size_t j = 0; j <
nFree(); j++) {
526 getStoichVector(j, nu);
529 doublereal dg_rt = 0.0;
530 for (
size_t k = 0; k < m_nsp; k++) {
531 dg_rt += m_mu[m_species[k]] * nu[k];
535 m_deltaG_RT[j] = dg_rt;
540 size_t k = m_order[j + m_nel];
542 if (m_moles[k] <= 0.0 && dg_rt > 0.0) {
547 }
else if (!m_solnrxn[j]) {
552 for (k = 0; k < m_nel; k++) {
553 size_t kc = m_order[k];
555 csum += pow(nu[kc], 2)*m_dsoln[kc]/nmoles;
559 size_t kc = m_order[j + m_nel];
561 double term1 = m_dsoln[kc]/nmoles;
564 doublereal sum = 0.0;
565 for (
size_t ip = 0; ip < m_mix->
nPhases(); ip++) {
569 for (k = 0; k < m_nsp; k++) {
572 psum += pow(nu[k], 2);
578 double rfctr = term1 + csum + sum;
579 if (fabs(rfctr) <
Tiny) {
582 fctr = 1.0/(term1 + csum + sum);
585 dxi[j] = -fctr*dg_rt;
587 for (
size_t m = 0; m < m_nel; m++) {
588 if (m_moles[m_order[m]] <= 0.0 && (m_N(m, j)*dxi[j] < 0.0)) {
592 grad += dxi[j]*dg_rt;
598 void MultiPhaseEquil::computeN()
601 std::vector<std::pair<double, size_t> > moleFractions(m_nsp);
602 for (
size_t k = 0; k < m_nsp; k++) {
604 moleFractions[k].first = - m_mix->
speciesMoles(m_species[k]);
605 moleFractions[k].second = k;
607 std::sort(moleFractions.begin(), moleFractions.end());
608 for (
size_t k = 0; k < m_nsp; k++) {
609 m_sortindex[k] = moleFractions[k].second;
612 for (
size_t m = 0; m < m_nel; m++) {
614 for (
size_t ik = 0; ik < m_nsp; ik++) {
616 if (m_mix->
nAtoms(m_species[k],m_element[m]) != 0) {
621 for (
size_t ij = 0; ij < m_nel; ij++) {
622 if (k == m_order[ij]) {
626 if (!ok || m_force) {
634 doublereal MultiPhaseEquil::error()
636 doublereal err, maxerr = 0.0;
639 for (
size_t j = 0; j <
nFree(); j++) {
640 size_t ik = j + m_nel;
644 if (!isStoichPhase(ik) && fabs(moles(ik)) <=
SmallNumber) {
650 if (isStoichPhase(ik) && moles(ik) <= 0.0 &&
651 m_deltaG_RT[j] >= 0.0) {
654 err = fabs(m_deltaG_RT[j]);
656 maxerr = std::max(maxerr, err);
661 double MultiPhaseEquil::phaseMoles(
size_t iph)
const 666 void MultiPhaseEquil::reportCSV(
const std::string& reportFile)
668 FILE* FP = fopen(reportFile.c_str(),
"w");
670 throw CanteraError(
"MultiPhaseEquil::reportCSV",
"Failure to open file");
682 for (
size_t iphase = 0; iphase < m_mix->
nPhases(); iphase++) {
684 ThermoPhase& tref = m_mix->
phase(iphase);
686 VolPM.resize(nSpecies, 0.0);
687 tref.getMoleFractions(&mf[istart]);
688 tref.getPartialMolarVolumes(VolPM.data());
690 double TMolesPhase = phaseMoles(iphase);
691 double VolPhaseVolumes = 0.0;
692 for (
size_t k = 0; k < nSpecies; k++) {
693 VolPhaseVolumes += VolPM[k] * mf[istart + k];
695 VolPhaseVolumes *= TMolesPhase;
696 vol += VolPhaseVolumes;
698 fprintf(FP,
"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT" 699 " -----------------------------\n");
700 fprintf(FP,
"Temperature = %11.5g kelvin\n", m_mix->
temperature());
701 fprintf(FP,
"Pressure = %11.5g Pascal\n", m_mix->
pressure());
702 fprintf(FP,
"Total Volume = %11.5g m**3\n", vol);
704 for (
size_t iphase = 0; iphase < m_mix->
nPhases(); iphase++) {
706 ThermoPhase& tref = m_mix->
phase(iphase);
707 ThermoPhase* tp = &tref;
709 string phaseName = tref.name();
710 double TMolesPhase = phaseMoles(iphase);
711 size_t nSpecies = tref.nSpecies();
712 activity.resize(nSpecies, 0.0);
713 ac.resize(nSpecies, 0.0);
714 mu0.resize(nSpecies, 0.0);
715 mu.resize(nSpecies, 0.0);
716 VolPM.resize(nSpecies, 0.0);
717 molalities.resize(nSpecies, 0.0);
718 int actConvention = tp->activityConvention();
719 tp->getActivities(activity.data());
720 tp->getActivityCoefficients(ac.data());
721 tp->getStandardChemPotentials(mu0.data());
722 tp->getPartialMolarVolumes(VolPM.data());
723 tp->getChemPotentials(mu.data());
724 double VolPhaseVolumes = 0.0;
725 for (
size_t k = 0; k < nSpecies; k++) {
726 VolPhaseVolumes += VolPM[k] * mf[istart + k];
728 VolPhaseVolumes *= TMolesPhase;
729 vol += VolPhaseVolumes;
730 if (actConvention == 1) {
731 MolalityVPSSTP* mTP =
static_cast<MolalityVPSSTP*
>(tp);
732 mTP->getMolalities(molalities.data());
733 tp->getChemPotentials(mu.data());
736 fprintf(FP,
" Name, Phase, PhaseMoles, Mole_Fract, " 737 "Molalities, ActCoeff, Activity," 738 "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
740 fprintf(FP,
" , , (kmol), , " 742 " (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n");
744 for (
size_t k = 0; k < nSpecies; k++) {
745 string sName = tp->speciesName(k);
746 fprintf(FP,
"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e," 747 "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
749 phaseName.c_str(), TMolesPhase,
750 mf[istart + k], molalities[k], ac[k], activity[k],
751 mu0[k]*1.0E-6, mu[k]*1.0E-6,
752 mf[istart + k] * TMolesPhase,
753 VolPM[k], VolPhaseVolumes);
757 fprintf(FP,
" Name, Phase, PhaseMoles, Mole_Fract, " 758 "Molalities, ActCoeff, Activity," 759 " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
761 fprintf(FP,
" , , (kmol), , " 763 " (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n");
765 for (
size_t k = 0; k < nSpecies; k++) {
768 for (
size_t k = 0; k < nSpecies; k++) {
769 string sName = tp->speciesName(k);
770 fprintf(FP,
"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, " 771 "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
773 phaseName.c_str(), TMolesPhase,
774 mf[istart + k], molalities[k], ac[k],
775 activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6,
776 mf[istart + k] * TMolesPhase,
777 VolPM[k], VolPhaseVolumes);
bool tempOK(size_t p) const
Return true if the phase p has valid thermo data for the current temperature.
size_t nPhases() const
Number of phases.
doublereal speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
doublereal temperature() const
Temperature [K].
doublereal computeReactionSteps(vector_fp &dxi)
Compute the change in extent of reaction for each reaction.
const size_t npos
index returned by functions to indicate "no position"
size_t nSpecies() const
Number of species, summed over all phases.
void getComponents(const std::vector< size_t > &order)
This method finds a set of component species and a complete set of formation reactions for the non-co...
size_t nSpecies() const
Returns the number of species in the phase.
doublereal pressure() const
Pressure [Pa].
Header for intermediate ThermoPhase object for phases which employ molality based activity coefficien...
int setInitialMoles(int loglevel=0)
Estimate the initial mole numbers.
Base class for a phase with thermodynamic properties.
A class for multiphase mixtures.
doublereal stepComposition(int loglevel=0)
Take one step in composition, given the gradient of G at the starting point, and a vector of reaction...
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
void unsort(vector_fp &x)
Re-arrange a vector of species properties in sorted form (components first) into unsorted, sequential form.
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
size_t speciesIndex(size_t k, size_t p) const
Return the global index of the species belonging to phase number p with local index k within the phas...
doublereal nAtoms(const size_t kGlob, const size_t mGlob) const
Returns the Number of atoms of global element mGlob in global species kGlob.
Base class for exceptions thrown by Cantera classes.
vector_fp & data()
Return a reference to the data vector.
size_t speciesPhaseIndex(const size_t kGlob) const
Returns the phase index of the Kth "global" species.
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
doublereal elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
const doublereal SmallNumber
smallest number to compare to zero.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
const doublereal Tiny
Small number to compare differences of mole fractions against.
bool solutionSpecies(size_t kGlob) const
Return true is species kGlob is a species in a multicomponent solution phase.
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Contains declarations for string manipulation functions within Cantera.
void getValidChemPotentials(doublereal not_mu, doublereal *mu, bool standard=false) const
Returns a vector of Valid chemical potentials.
std::string elementName(size_t m) const
Returns the name of the global element m.
size_t nFree() const
Number of degrees of freedom.
thermo_t & phase(size_t n)
Return a reference to phase n.
doublereal phaseMoles(const size_t n) const
Return the number of moles in phase n.
Namespace for the Cantera kernel.
void finish()
Clean up the composition.
void setMoles(const doublereal *n)
Sets all of the global species mole numbers.
size_t nElements() const
Number of elements.
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.