00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "MultiPhase.h"
00013 #include "MultiPhaseEquil.h"
00014
00015 #include "ThermoPhase.h"
00016 #include "DenseMatrix.h"
00017 #include "stringUtils.h"
00018 #include "global.h"
00019
00020 using namespace std;
00021
00022 namespace Cantera {
00023
00024
00025 MultiPhase::MultiPhase() :
00026 m_np(0),
00027 m_temp(0.0),
00028 m_press(0.0),
00029 m_nel(0),
00030 m_nsp(0),
00031 m_init(false),
00032 m_eloc(-1),
00033 m_Tmin(1.0),
00034 m_Tmax(100000.0)
00035 {
00036 }
00037
00038 void MultiPhase::
00039 addPhases(MultiPhase& mix) {
00040 index_t n;
00041 for (n = 0; n < mix.m_np; n++) {
00042 addPhase(mix.m_phase[n], mix.m_moles[n]);
00043 }
00044 }
00045
00046 void MultiPhase::
00047 addPhases(phase_list& phases, const vector_fp& phaseMoles) {
00048 index_t np = phases.size();
00049 index_t n;
00050 for (n = 0; n < np; n++) {
00051 addPhase(phases[n], phaseMoles[n]);
00052 }
00053 init();
00054 }
00055
00056 void MultiPhase::
00057 addPhase(phase_t* p, doublereal moles) {
00058 if (m_init) {
00059 throw CanteraError("addPhase",
00060 "phases cannot be added after init() has been called.");
00061 }
00062
00063
00064 m_phase.push_back(p);
00065
00066
00067 m_moles.push_back(moles);
00068 m_temp_OK.push_back(true);
00069
00070
00071
00072 m_np = m_phase.size();
00073 m_nsp += p->nSpecies();
00074
00075
00076
00077
00078
00079 string ename;
00080
00081 index_t m, nel = p->nElements();
00082 for (m = 0; m < nel; m++) {
00083 ename = p->elementName(m);
00084
00085
00086
00087
00088
00089 if (m_enamemap.find(ename) == m_enamemap.end()) {
00090 m_enamemap[ename] = m_nel + 1;
00091 m_enames.push_back(ename);
00092 m_atomicNumber.push_back(p->atomicNumber(m));
00093
00094
00095 if (ename == "E" || ename == "e") m_eloc = m_nel;
00096
00097 m_nel++;
00098 }
00099 }
00100
00101
00102
00103
00104 if (m_temp == 0.0 && p->temperature() > 0.0) {
00105 m_temp = p->temperature();
00106 m_press = p->pressure();
00107 }
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117 if (p->nSpecies() > 1) {
00118 double t = p->minTemp();
00119 if (t > m_Tmin) m_Tmin = t;
00120 t = p->maxTemp();
00121 if (t < m_Tmax) m_Tmax = t;
00122 }
00123 }
00124
00125
00126
00127
00128
00129
00130 void MultiPhase::init() {
00131 if (m_init) return;
00132 index_t ip, kp, k = 0, nsp, m;
00133 int mlocal;
00134 string sym;
00135
00136
00137 m_atoms.resize(m_nel, m_nsp, 0.0);
00138 m_moleFractions.resize(m_nsp, 0.0);
00139 m_elemAbundances.resize(m_nel, 0.0);
00140
00141
00142
00143
00144 for (m = 0; m < m_nel; m++) {
00145 sym = m_enames[m];
00146 k = 0;
00147
00148 for (ip = 0; ip < m_np; ip++) {
00149 phase_t* p = m_phase[ip];
00150 nsp = p->nSpecies();
00151 mlocal = p->elementIndex(sym);
00152 for (kp = 0; kp < nsp; kp++) {
00153 if (mlocal >= 0) {
00154 m_atoms(m, k) = p->nAtoms(kp, mlocal);
00155 }
00156 if (m == 0) {
00157 m_snames.push_back(p->speciesName(kp));
00158 if (kp == 0) {
00159 m_spstart.push_back(m_spphase.size());
00160 }
00161 m_spphase.push_back(ip);
00162 }
00163 k++;
00164 }
00165 }
00166 }
00167
00168 if (m_eloc >= 0) {
00169 doublereal esum;
00170 for (k = 0; k < m_nsp; k++) {
00171 esum = 0.0;
00172 for (m = 0; m < m_nel; m++) {
00173 if (int(m) != m_eloc)
00174 esum += m_atoms(m,k) * m_atomicNumber[m];
00175 }
00176
00177 }
00178 }
00179
00180
00181
00182 m_init = true;
00183
00184 updateMoleFractions();
00185
00186 }
00187
00188
00189
00190
00191
00192 MultiPhase::phase_t& MultiPhase::phase(index_t n) {
00193 if (!m_init) init();
00194 m_phase[n]->setTemperature(m_temp);
00195 m_phase[n]->setMoleFractions_NoNorm(DATA_PTR(m_moleFractions) + m_spstart[n]);
00196 m_phase[n]->setPressure(m_press);
00197 return *m_phase[n];
00198 }
00199
00200
00201 doublereal MultiPhase::speciesMoles(index_t k) const {
00202 index_t ip = m_spphase[k];
00203 return m_moles[ip]*m_moleFractions[k];
00204 }
00205
00206
00207
00208 doublereal MultiPhase::elementMoles(index_t m) const {
00209 doublereal sum = 0.0, phasesum;
00210 index_t i, k = 0, ik, nsp;
00211 for (i = 0; i < m_np; i++) {
00212 phasesum = 0.0;
00213 nsp = m_phase[i]->nSpecies();
00214 for (ik = 0; ik < nsp; ik++) {
00215 k = speciesIndex(ik, i);
00216 phasesum += m_atoms(m,k)*m_moleFractions[k];
00217 }
00218 sum += phasesum * m_moles[i];
00219 }
00220 return sum;
00221 }
00222
00223
00224 doublereal MultiPhase::charge() const {
00225 doublereal sum = 0.0;
00226 index_t i;
00227 for (i = 0; i < m_np; i++) {
00228 sum += phaseCharge(i);
00229 }
00230 return sum;
00231 }
00232
00233 int MultiPhase::speciesIndex(std::string speciesName, std::string phaseName) {
00234 int p = phaseIndex(phaseName);
00235 if (p < 0) {
00236 throw CanteraError("MultiPhase::speciesIndex", "phase not found: " + phaseName);
00237 }
00238 int k = m_phase[p]->speciesIndex(speciesName);
00239 if (k < 0) {
00240 throw CanteraError("MultiPhase::speciesIndex", "species not found: " + speciesName);
00241 }
00242 return m_spstart[p] + k;
00243 }
00244
00245
00246
00247
00248
00249 doublereal MultiPhase::phaseCharge(index_t p) const {
00250 doublereal phasesum = 0.0;
00251 int ik, k, nsp = m_phase[p]->nSpecies();
00252 for (ik = 0; ik < nsp; ik++) {
00253 k = speciesIndex(ik, p);
00254 phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
00255 }
00256 return Faraday*phasesum*m_moles[p];
00257 }
00258
00259
00260
00261 void MultiPhase::getChemPotentials(doublereal* mu) const {
00262 index_t i, loc = 0;
00263 updatePhases();
00264 for (i = 0; i < m_np; i++) {
00265 m_phase[i]->getChemPotentials(mu + loc);
00266 loc += m_phase[i]->nSpecies();
00267 }
00268 }
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 void MultiPhase::getValidChemPotentials(doublereal not_mu,
00298 doublereal* mu, bool standard) const {
00299 index_t i, loc = 0;
00300
00301 updatePhases();
00302
00303 for (i = 0; i < m_np; i++) {
00304 if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
00305 if (!standard)
00306 m_phase[i]->getChemPotentials(mu + loc);
00307 else
00308 m_phase[i]->getStandardChemPotentials(mu + loc);
00309 }
00310 else
00311 fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu);
00312 loc += m_phase[i]->nSpecies();
00313 }
00314 }
00315
00316
00317 bool MultiPhase::solutionSpecies(index_t k) const {
00318 if (m_phase[m_spphase[k]]->nSpecies() > 1)
00319 return true;
00320 else
00321 return false;
00322 }
00323
00324
00325 doublereal MultiPhase::gibbs() const {
00326 index_t i;
00327 doublereal sum = 0.0;
00328 updatePhases();
00329 for (i = 0; i < m_np; i++)
00330 sum += m_phase[i]->gibbs_mole() * m_moles[i];
00331 return sum;
00332 }
00333
00334
00335 doublereal MultiPhase::enthalpy() const {
00336 index_t i;
00337 doublereal sum = 0.0;
00338 updatePhases();
00339 for (i = 0; i < m_np; i++)
00340 sum += m_phase[i]->enthalpy_mole() * m_moles[i];
00341 return sum;
00342 }
00343
00344
00345 doublereal MultiPhase::IntEnergy() const {
00346 index_t i;
00347 doublereal sum = 0.0;
00348 updatePhases();
00349 for (i = 0; i < m_np; i++)
00350 sum += m_phase[i]->intEnergy_mole() * m_moles[i];
00351 return sum;
00352 }
00353
00354
00355 doublereal MultiPhase::entropy() const {
00356 index_t i;
00357 doublereal sum = 0.0;
00358 updatePhases();
00359 for (i = 0; i < m_np; i++)
00360 sum += m_phase[i]->entropy_mole() * m_moles[i];
00361 return sum;
00362 }
00363
00364
00365
00366
00367 doublereal MultiPhase::cp() const {
00368 index_t i;
00369 doublereal sum = 0.0;
00370 updatePhases();
00371 for (i = 0; i < m_np; i++)
00372 sum += m_phase[i]->cp_mole() * m_moles[i];
00373 return sum;
00374 }
00375
00376
00377
00378
00379
00380 void MultiPhase::setPhaseMoleFractions(const index_t n, const doublereal* const x) {
00381 phase_t* p = m_phase[n];
00382 p->setState_TPX(m_temp, m_press, x);
00383 int nsp = p->nSpecies();
00384 int istart = m_spstart[n];
00385 for (int k = 0; k < nsp; k++) {
00386 m_moleFractions[istart+k] = x[k];
00387 }
00388 }
00389
00390
00391
00392
00393 void MultiPhase::setMolesByName(compositionMap& xMap) {
00394 int kk = nSpecies();
00395 doublereal x;
00396 vector_fp moles(kk, 0.0);
00397 for (int k = 0; k < kk; k++) {
00398 x = xMap[speciesName(k)];
00399 if (x > 0.0) moles[k] = x;
00400 }
00401 setMoles(DATA_PTR(moles));
00402 }
00403
00404
00405
00406 void MultiPhase::setMolesByName(const std::string& x) {
00407 compositionMap xx;
00408
00409
00410
00411
00412 int kk = nSpecies();
00413 for (int k = 0; k < kk; k++) {
00414 xx[speciesName(k)] = -1.0;
00415 }
00416
00417
00418
00419 parseCompString(x, xx);
00420 setMolesByName(xx);
00421 }
00422
00423
00424
00425 void MultiPhase::getMoles(doublereal * molNum) const {
00426
00427
00428
00429 copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
00430 index_t ik;
00431 doublereal *dtmp = molNum;
00432 for (index_t ip = 0; ip < m_np; ip++) {
00433 doublereal phasemoles = m_moles[ip];
00434 phase_t* p = m_phase[ip];
00435 index_t nsp = p->nSpecies();
00436 for (ik = 0; ik < nsp; ik++) {
00437 *(dtmp++) *= phasemoles;
00438 }
00439 }
00440 }
00441
00442
00443
00444
00445 void MultiPhase::setMoles(const doublereal* n) {
00446 if (!m_init) init();
00447 index_t ip, loc = 0;
00448 index_t ik, k = 0, nsp;
00449 doublereal phasemoles;
00450 for (ip = 0; ip < m_np; ip++) {
00451 phase_t* p = m_phase[ip];
00452 nsp = p->nSpecies();
00453 phasemoles = 0.0;
00454 for (ik = 0; ik < nsp; ik++) {
00455 phasemoles += n[k];
00456 k++;
00457 }
00458 m_moles[ip] = phasemoles;
00459 if (nsp > 1) {
00460 if (phasemoles > 0.0) {
00461 p->setState_TPX(m_temp, m_press, n + loc);
00462 p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
00463 } else {
00464 p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
00465 }
00466 }
00467 else {
00468 m_moleFractions[loc] = 1.0;
00469 }
00470 loc += nsp;
00471 }
00472 }
00473
00474 void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles) {
00475 vector_fp tmpMoles(m_nsp, 0.0);
00476 getMoles(DATA_PTR(tmpMoles));
00477 tmpMoles[indexS] += addedMoles;
00478 if (tmpMoles[indexS] < 0.0) {
00479 tmpMoles[indexS] = 0.0;
00480 }
00481 setMoles(DATA_PTR(tmpMoles));
00482 }
00483
00484 void MultiPhase::setState_TP(const doublereal T, const doublereal Pres) {
00485 if (!m_init) init();
00486 m_temp = T;
00487 m_press = Pres;
00488 updatePhases();
00489 }
00490
00491 void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres,
00492 const doublereal *n) {
00493 m_temp = T;
00494 m_press = Pres;
00495 setMoles(n);
00496 }
00497
00498 void MultiPhase::getElemAbundances(doublereal *elemAbundances) const {
00499 index_t eGlobal;
00500 calcElemAbundances();
00501 for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
00502 elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
00503 }
00504 }
00505
00506
00507 void MultiPhase::calcElemAbundances() const {
00508 index_t loc = 0;
00509 index_t eGlobal;
00510 int ik, kGlobal;
00511 doublereal spMoles;
00512 for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
00513 m_elemAbundances[eGlobal] = 0.0;
00514 }
00515 for (index_t ip = 0; ip < m_np; ip++) {
00516 phase_t* p = m_phase[ip];
00517 int nspPhase = p->nSpecies();
00518 doublereal phasemoles = m_moles[ip];
00519 for (ik = 0; ik < nspPhase; ik++) {
00520 kGlobal = loc + ik;
00521 spMoles = m_moleFractions[kGlobal] * phasemoles;
00522 for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
00523 m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
00524 }
00525 }
00526 loc += nspPhase;
00527 }
00528 }
00529
00530
00531 doublereal MultiPhase::volume() const {
00532 int i;
00533 doublereal sum = 0;
00534 for (i = 0; i < int(m_np); i++) {
00535 sum += m_moles[i]/m_phase[i]->molarDensity();
00536 }
00537 return sum;
00538 }
00539
00540 doublereal MultiPhase::equilibrate(int XY, doublereal err,
00541 int maxsteps, int maxiter, int loglevel) {
00542 doublereal error;
00543 bool strt = false;
00544 doublereal dt;
00545 doublereal h0;
00546 int n;
00547 bool start;
00548 doublereal ferr, hnow, herr = 1.0;
00549 doublereal snow, serr = 1.0, s0;
00550 doublereal Tlow = -1.0, Thigh = -1.0;
00551 doublereal Hlow = Undef, Hhigh = Undef, tnew;
00552 doublereal dta=0.0, dtmax, cpb;
00553 MultiPhaseEquil* e = 0;
00554
00555 if (!m_init) init();
00556 if (loglevel > 0)
00557 beginLogGroup("MultiPhase::equilibrate", loglevel);
00558
00559 if (XY == TP) {
00560 if (loglevel > 0) {
00561 addLogEntry("problem type","fixed T,P");
00562 addLogEntry("Temperature",temperature());
00563 addLogEntry("Pressure", pressure());
00564 }
00565
00566
00567 e = new MultiPhaseEquil(this);
00568 try {
00569 error = e->equilibrate(XY, err, maxsteps, loglevel);
00570 }
00571 catch (CanteraError &err) {
00572 if (loglevel > 0)
00573 endLogGroup();
00574 delete e;
00575 e = 0;
00576 throw err;
00577 }
00578 goto done;
00579 }
00580
00581 else if (XY == HP) {
00582 h0 = enthalpy();
00583 Tlow = 0.5*m_Tmin;
00584 Thigh = 2.0*m_Tmax;
00585 if (loglevel > 0) {
00586 addLogEntry("problem type","fixed H,P");
00587 addLogEntry("H target",fp2str(h0));
00588 }
00589 for (n = 0; n < maxiter; n++) {
00590
00591
00592
00593
00594
00595
00596
00597 e = new MultiPhaseEquil(this, strt);
00598
00599
00600 if (loglevel > 0)
00601 beginLogGroup("iteration "+int2str(n));
00602
00603 try {
00604 error = e->equilibrate(TP, err, maxsteps, loglevel);
00605 hnow = enthalpy();
00606
00607
00608
00609 if (hnow < h0) {
00610 if (m_temp > Tlow) {
00611 Tlow = m_temp;
00612 Hlow = hnow;
00613 }
00614 }
00615
00616
00617 else {
00618 if (m_temp < Thigh) {
00619 Thigh = m_temp;
00620 Hhigh = hnow;
00621 }
00622 }
00623 if (Hlow != Undef && Hhigh != Undef) {
00624 cpb = (Hhigh - Hlow)/(Thigh - Tlow);
00625 dt = (h0 - hnow)/cpb;
00626 dta = fabs(dt);
00627 dtmax = 0.5*fabs(Thigh - Tlow);
00628 if (dta > dtmax) dt *= dtmax/dta;
00629 }
00630 else {
00631 tnew = sqrt(Tlow*Thigh);
00632 dt = tnew - m_temp;
00633
00634 }
00635
00636 herr = fabs((h0 - hnow)/h0);
00637 if (loglevel > 0) {
00638 addLogEntry("T",fp2str(temperature()));
00639 addLogEntry("H",fp2str(hnow));
00640 addLogEntry("H rel error",fp2str(herr));
00641 addLogEntry("lower T bound",fp2str(Tlow));
00642 addLogEntry("upper T bound",fp2str(Thigh));
00643 endLogGroup();
00644 }
00645
00646
00647 if (herr < err) {
00648 if (loglevel > 0) {
00649 addLogEntry("T iterations",int2str(n));
00650 addLogEntry("Final T",fp2str(temperature()));
00651 addLogEntry("H rel error",fp2str(herr));
00652 }
00653 goto done;
00654 }
00655 tnew = m_temp + dt;
00656 if (tnew < 0.0) tnew = 0.5*m_temp;
00657
00658 setTemperature(tnew);
00659
00660
00661
00662 if (dta < 100.0) strt = false;
00663
00664 }
00665
00666 catch (CanteraError err) {
00667 if (!strt) {
00668 if (loglevel > 0)
00669 addLogEntry("no convergence",
00670 "try estimating starting composition");
00671 strt = true;
00672 }
00673 else {
00674 tnew = 0.5*(m_temp + Thigh);
00675 if (fabs(tnew - m_temp) < 1.0) tnew = m_temp + 1.0;
00676 setTemperature(tnew);
00677 if (loglevel > 0)
00678 addLogEntry("no convergence",
00679 "trying T = "+fp2str(m_temp));
00680 }
00681 if (loglevel > 0)
00682 endLogGroup();
00683 }
00684 delete e;
00685 e = 0;
00686 }
00687 if (loglevel > 0) {
00688 addLogEntry("reached max number of T iterations",int2str(maxiter));
00689 endLogGroup();
00690 }
00691 throw CanteraError("MultiPhase::equilibrate",
00692 "No convergence for T");
00693 }
00694 else if (XY == SP) {
00695 s0 = entropy();
00696 start = true;
00697 Tlow = 1.0;
00698 Thigh = 1.0e6;
00699 if (loglevel > 0) {
00700 addLogEntry("problem type","fixed S,P");
00701 addLogEntry("S target",fp2str(s0));
00702 addLogEntry("min T",fp2str(Tlow));
00703 addLogEntry("max T",fp2str(Thigh));
00704 }
00705 for (n = 0; n < maxiter; n++) {
00706 if (e) delete e;
00707 e = new MultiPhaseEquil(this, strt);
00708 ferr = 0.1;
00709 if (fabs(dt) < 1.0) ferr = err;
00710
00711 if (loglevel > 0)
00712 beginLogGroup("iteration "+int2str(n));
00713
00714 try {
00715 error = e->equilibrate(TP, err, maxsteps, loglevel);
00716 snow = entropy();
00717 if (snow < s0) {
00718 if (m_temp > Tlow) Tlow = m_temp;
00719 }
00720 else {
00721 if (m_temp < Thigh) Thigh = m_temp;
00722 }
00723 serr = fabs((s0 - snow)/s0);
00724 if (loglevel > 0) {
00725 addLogEntry("T",fp2str(temperature()));
00726 addLogEntry("S",fp2str(snow));
00727 addLogEntry("S rel error",fp2str(serr));
00728 endLogGroup();
00729 }
00730 dt = (s0 - snow)*m_temp/cp();
00731 dtmax = 0.5*fabs(Thigh - Tlow);
00732 dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
00733 dta = fabs(dt);
00734 if (dta > dtmax) dt *= dtmax/dta;
00735 if (herr < err || dta < 1.0e-4) {
00736 if (loglevel > 0) {
00737 addLogEntry("T iterations",int2str(n));
00738 addLogEntry("Final T",fp2str(temperature()));
00739 addLogEntry("S rel error",fp2str(serr));
00740 }
00741 goto done;
00742 }
00743 tnew = m_temp + dt;
00744 setTemperature(tnew);
00745
00746
00747
00748 if (dta < 100.0) strt = false;
00749 }
00750
00751 catch (CanteraError err) {
00752 if (!strt) {
00753 if (loglevel > 0) {
00754 addLogEntry("no convergence",
00755 "setting strt to True");
00756 }
00757 strt = true;
00758 }
00759 else {
00760 tnew = 0.5*(m_temp + Thigh);
00761 setTemperature(tnew);
00762 if (loglevel > 0) {
00763 addLogEntry("no convergence",
00764 "trying T = "+fp2str(m_temp));
00765 }
00766 }
00767 if (loglevel > 0)
00768 endLogGroup();
00769 }
00770 delete e;
00771 e = 0;
00772 }
00773 if (loglevel > 0) {
00774 addLogEntry("reached max number of T iterations",int2str(maxiter));
00775 endLogGroup();
00776 }
00777 throw CanteraError("MultiPhase::equilibrate",
00778 "No convergence for T");
00779 }
00780 else if (XY == TV) {
00781 addLogEntry("problem type","fixed T, V");
00782
00783 doublereal v0 = volume();
00784 doublereal dVdP;
00785 int n;
00786 bool start = true;
00787 doublereal error, vnow, pnow, verr;
00788 for (n = 0; n < maxiter; n++) {
00789 pnow = pressure();
00790 MultiPhaseEquil e(this, start);
00791 start = false;
00792 beginLogGroup("iteration "+int2str(n));
00793
00794 error = e.equilibrate(TP, err, maxsteps, loglevel);
00795 vnow = volume();
00796 verr = fabs((v0 - vnow)/v0);
00797 addLogEntry("P",fp2str(pressure()));
00798 addLogEntry("V rel error",fp2str(verr));
00799 endLogGroup();
00800
00801 if (verr < err) {
00802 addLogEntry("P iterations",int2str(n));
00803 addLogEntry("Final P",fp2str(pressure()));
00804 addLogEntry("V rel error",fp2str(verr));
00805 goto done;
00806 }
00807
00808 setPressure(pnow*1.01);
00809 dVdP = (volume() - vnow)/(0.01*pnow);
00810 setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
00811 }
00812 }
00813
00814 else {
00815 if (loglevel > 0)
00816 endLogGroup();
00817 throw CanteraError("MultiPhase::equilibrate","unknown option");
00818 }
00819 return -1.0;
00820 done:
00821 delete e;
00822 e = 0;
00823 if (loglevel > 0)
00824 endLogGroup();
00825 return err;
00826 }
00827
00828 #ifdef MULTIPHASE_DEVEL
00829 void importFromXML(string infile, string id) {
00830 XML_Node* root = get_XML_File(infile);
00831 if (id == "-") id = "";
00832 XML_Node* x = get_XML_Node(string("#")+id, root);
00833 if (x.name() != "multiphase")
00834 throw CanteraError("MultiPhase::importFromXML",
00835 "Current XML_Node is not a multiphase element.");
00836 vector<XML_Node*> phases;
00837 x.getChildren("phase",phases);
00838 int np = phases.size();
00839 int n;
00840 ThermoPhase* p;
00841 for (n = 0; n < np; n++) {
00842 XML_Node& ph = *phases[n];
00843 srcfile = infile;
00844 if (ph.hasAttrib("src")) srcfile = ph["src"];
00845 idstr = ph["id"];
00846 p = newPhase(srcfile, idstr);
00847 if (p) {
00848 addPhase(p, ph.value());
00849 }
00850 }
00851 }
00852 #endif
00853
00854 void MultiPhase::setTemperature(const doublereal T) {
00855 if (!m_init) init();
00856 m_temp = T;
00857 updatePhases();
00858 }
00859
00860
00861 std::string MultiPhase::elementName(int m) const {
00862 return m_enames[m];
00863 }
00864
00865
00866 int MultiPhase::elementIndex(std::string name) const {
00867 for (size_t e = 0; e < m_nel; e++) {
00868 if (m_enames[e] == name) {
00869 return (int) e;
00870 }
00871 }
00872 return -1;
00873 }
00874
00875
00876 std::string MultiPhase::speciesName(const int k) const {
00877 return m_snames[k];
00878 }
00879
00880 doublereal MultiPhase::nAtoms(const int kGlob, const int mGlob) const {
00881 return m_atoms(mGlob, kGlob);
00882 }
00883
00884 void MultiPhase::getMoleFractions(doublereal* const x) const {
00885 std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
00886 }
00887
00888 std::string MultiPhase::phaseName(const index_t iph) const {
00889 const phase_t *tptr = m_phase[iph];
00890 return tptr->id();
00891 }
00892
00893 int MultiPhase::phaseIndex(const std::string &pName) const {
00894 std::string tmp;
00895 for (int iph = 0; iph < (int) m_np; iph++) {
00896 const phase_t *tptr = m_phase[iph];
00897 tmp = tptr->id();
00898 if (tmp == pName) {
00899 return iph;
00900 }
00901 }
00902 return -1;
00903 }
00904
00905 doublereal MultiPhase::phaseMoles(const index_t n) const {
00906 return m_moles[n];
00907 }
00908
00909 void MultiPhase::setPhaseMoles(const index_t n, const doublereal moles) {
00910 m_moles[n] = moles;
00911 }
00912
00913 int MultiPhase::speciesPhaseIndex(const index_t kGlob) const {
00914 return m_spphase[kGlob];
00915 }
00916
00917 doublereal MultiPhase::moleFraction(const index_t kGlob) const{
00918 return m_moleFractions[kGlob];
00919 }
00920
00921
00922 bool MultiPhase::tempOK(const index_t p) const {
00923 return m_temp_OK[p];
00924 }
00925
00926
00927 void MultiPhase::updateMoleFractions() {
00928 uploadMoleFractionsFromPhases();
00929 }
00930
00931 void MultiPhase::uploadMoleFractionsFromPhases() {
00932 index_t ip, loc = 0;
00933 for (ip = 0; ip < m_np; ip++) {
00934 phase_t* p = m_phase[ip];
00935 p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
00936 loc += p->nSpecies();
00937 }
00938 calcElemAbundances();
00939 }
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953 void MultiPhase::updatePhases() const {
00954 index_t p, nsp, loc = 0;
00955 for (p = 0; p < m_np; p++) {
00956 nsp = m_phase[p]->nSpecies();
00957 const doublereal* x = DATA_PTR(m_moleFractions) + loc;
00958 loc += nsp;
00959 m_phase[p]->setState_TPX(m_temp, m_press, x);
00960 m_temp_OK[p] = true;
00961 if (m_temp < m_phase[p]->minTemp()
00962 || m_temp > m_phase[p]->maxTemp()) {
00963 m_temp_OK[p] = false;
00964 }
00965 }
00966 }
00967
00968 }
00969