Cantera 2.6.0
MultiPhaseEquil.cpp
Go to the documentation of this file.
1/**
2 * @file MultiPhaseEquil.cpp
3 */
4
5// This file is part of Cantera. See License.txt in the top-level directory or
6// at https://cantera.org/license.txt for license and copyright information.
7
11
12#include <cstdio>
13#include <numeric>
14
15using namespace std;
16
17namespace Cantera
18{
19
20MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_mix(mix)
21{
22 // store some mixture parameters locally
23 m_nel_mix = mix->nElements();
24 m_nsp_mix = mix->nSpecies();
25 m_press = mix->pressure();
26 m_temp = mix->temperature();
27
28 m_force = true;
29 m_nel = 0;
30 m_nsp = 0;
31 m_eloc = 1000;
32 m_incl_species.resize(m_nsp_mix,1);
33 m_incl_element.resize(m_nel_mix,1);
34 for (size_t m = 0; m < m_nel_mix; m++) {
35 string enm = mix->elementName(m);
36 // element 'E' or 'e' represents an electron; this requires special
37 // handling, so save its index for later use
38 if (enm == "E" || enm == "e") {
39 m_eloc = m;
40 }
41 // if an element other than electrons is not present in the mixture,
42 // then exclude it and all species containing it from the calculation.
43 // Electrons are a special case, since a species can have a negative
44 // number of 'atoms' of electrons (positive ions).
45 if (m_mix->elementMoles(m) <= 0.0 && m != m_eloc) {
46 m_incl_element[m] = 0;
47 for (size_t k = 0; k < m_nsp_mix; k++) {
48 if (m_mix->nAtoms(k,m) != 0.0) {
49 m_incl_species[k] = 0;
50 }
51 }
52 }
53 }
54
55 // Now build the list of elements to be included, starting with
56 // electrons, if they are present.
57 if (m_eloc < m_nel_mix) {
58 m_element.push_back(m_eloc);
59 m_nel++;
60 }
61 // add the included elements other than electrons
62 for (size_t m = 0; m < m_nel_mix; m++) {
63 if (m_incl_element[m] == 1 && m != m_eloc) {
64 m_nel++;
65 m_element.push_back(m);
66 }
67 }
68
69 // include pure single-constituent phases only if their thermo data are
70 // valid for this temperature. This is necessary, since some thermo
71 // polynomial fits are done only for a limited temperature range. For
72 // example, using the NASA polynomial fits for solid ice and liquid water,
73 // if this were not done the calculation would predict solid ice to be
74 // present far above its melting point, since the thermo polynomial fits
75 // only extend to 273.15 K, and give unphysical results above this
76 // temperature, leading (incorrectly) to Gibbs free energies at high
77 // temperature lower than for liquid water.
78 for (size_t k = 0; k < m_nsp_mix; k++) {
79 size_t ip = m_mix->speciesPhaseIndex(k);
80 if (!m_mix->solutionSpecies(k) &&
81 !m_mix->tempOK(ip)) {
82 m_incl_species[k] = 0;
83 if (m_mix->speciesMoles(k) > 0.0) {
84 throw CanteraError("MultiPhaseEquil::MultiPhaseEquil",
85 "condensed-phase species"+ m_mix->speciesName(k)
86 + " is excluded since its thermo properties are \n"
87 "not valid at this temperature, but it has "
88 "non-zero moles in the initial state.");
89 }
90 }
91 }
92
93 // Now build the list of all species to be included in the calculation.
94 for (size_t k = 0; k < m_nsp_mix; k++) {
95 if (m_incl_species[k] ==1) {
96 m_nsp++;
97 m_species.push_back(k);
98 }
99 }
100
101 // some work arrays for internal use
102 m_work.resize(m_nsp);
103 m_work2.resize(m_nsp);
104 m_work3.resize(m_nsp_mix);
105 m_mu.resize(m_nsp_mix);
106
107 // number of moles of each species
108 m_moles.resize(m_nsp);
109 m_lastmoles.resize(m_nsp);
110 m_dxi.resize(nFree());
111
112 // initialize the mole numbers to the mixture composition
113 for (size_t ik = 0; ik < m_nsp; ik++) {
114 m_moles[ik] = m_mix->speciesMoles(m_species[ik]);
115 }
116
117 // Delta G / RT for each reaction
118 m_deltaG_RT.resize(nFree(), 0.0);
119 m_majorsp.resize(m_nsp);
120 m_sortindex.resize(m_nsp,0);
121 m_lastsort.resize(m_nel);
122 m_solnrxn.resize(nFree());
123 m_A.resize(m_nel, m_nsp, 0.0);
124 m_N.resize(m_nsp, nFree());
125 m_order.resize(std::max(m_nsp, m_nel), 0);
126 iota(m_order.begin(), m_order.begin() + m_nsp, 0);
127
128 // if the 'start' flag is set, estimate the initial mole numbers by doing a
129 // linear Gibbs minimization. In this case, only the elemental composition
130 // of the initial mixture state matters.
131 if (start) {
132 setInitialMoles(loglevel-1);
133 }
134 computeN();
135
136 // Take a very small step in composition space, so that no
137 // species has precisely zero moles.
138 vector_fp dxi(nFree(), 1.0e-20);
139 if (!dxi.empty()) {
140 multiply(m_N, dxi.data(), m_work.data());
141 unsort(m_work);
142 }
143
144 for (size_t k = 0; k < m_nsp; k++) {
145 m_moles[k] += m_work[k];
146 m_lastmoles[k] = m_moles[k];
147 if (m_mix->solutionSpecies(m_species[k])) {
148 m_dsoln.push_back(1);
149 } else {
150 m_dsoln.push_back(0);
151 }
152 }
153 m_force = false;
154 updateMixMoles();
155
156 // At this point, the instance has been created, the species to be included
157 // have been determined, and an initial composition has been selected that
158 // has all non-zero mole numbers for the included species.
159}
160
161doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err,
162 int maxsteps, int loglevel)
163{
164 int i;
165 m_iter = 0;
166 for (i = 0; i < maxsteps; i++) {
167 stepComposition(loglevel-1);
168 if (error() < err) {
169 break;
170 }
171 }
172 if (i >= maxsteps) {
173 throw CanteraError("MultiPhaseEquil::equilibrate",
174 "no convergence in {} iterations. Error = {}",
175 maxsteps, error());
176 }
177 finish();
178 return error();
179}
180
181void MultiPhaseEquil::updateMixMoles()
182{
183 fill(m_work3.begin(), m_work3.end(), 0.0);
184 for (size_t k = 0; k < m_nsp; k++) {
185 m_work3[m_species[k]] = m_moles[k];
186 }
187 m_mix->setMoles(m_work3.data());
188}
189
191{
192 fill(m_work3.begin(), m_work3.end(), 0.0);
193 for (size_t k = 0; k < m_nsp; k++) {
194 m_work3[m_species[k]] = (m_moles[k] > 0.0 ? m_moles[k] : 0.0);
195 }
196 m_mix->setMoles(m_work3.data());
197}
198
200{
201 double not_mu = 1.0e12;
202 m_mix->getValidChemPotentials(not_mu, m_mu.data(), true);
203 double dxi_min = 1.0e10;
204 bool redo = true;
205 int iter = 0;
206
207 while (redo) {
208 // choose a set of components based on the current composition
209 computeN();
210 redo = false;
211 iter++;
212 if (iter > 4) {
213 break;
214 }
215
216 // loop over all reactions
217 for (size_t j = 0; j < nFree(); j++) {
218 double dg_rt = 0.0;
219 dxi_min = 1.0e10;
220 for (size_t ik = 0; ik < m_nsp; ik++) {
221 dg_rt += mu(ik) * m_N(ik,j);
222 }
223
224 // fwd or rev direction
225 int idir = (dg_rt < 0.0 ? 1 : -1);
226
227 for (size_t ik = 0; ik < m_nsp; ik++) {
228 double nu = m_N(ik, j);
229
230 // set max change in progress variable by
231 // non-negativity requirement
232 // -> Note, 0.99 factor is so that difference of 2 numbers
233 // isn't zero. This causes differences between
234 // optimized and debug versions of the code
235 if (nu*idir < 0) {
236 double delta_xi = fabs(0.99*moles(ik)/nu);
237 // if a component has nearly zero moles, redo
238 // with a new set of components
239 if (!redo && delta_xi < 1.0e-10 && ik < m_nel) {
240 redo = true;
241 }
242 dxi_min = std::min(dxi_min, delta_xi);
243 }
244 }
245 // step the composition by dxi_min
246 for (size_t ik = 0; ik < m_nsp; ik++) {
247 moles(ik) += m_N(ik, j) * idir*dxi_min;
248 }
249 }
250 // set the moles of the phase objects to match
251 updateMixMoles();
252 }
253 return 0;
254}
255
256void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
257{
258 // if the input species array has the wrong size, ignore it
259 // and consider the species for components in declaration order.
260 if (order.size() != m_nsp) {
261 for (size_t k = 0; k < m_nsp; k++) {
262 m_order[k] = k;
263 }
264 } else {
265 for (size_t k = 0; k < m_nsp; k++) {
266 m_order[k] = order[k];
267 }
268 }
269
270 size_t nRows = m_nel;
271 size_t nColumns = m_nsp;
272
273 // set up the atomic composition matrix
274 for (size_t m = 0; m < nRows; m++) {
275 for (size_t k = 0; k < nColumns; k++) {
276 m_A(m, k) = m_mix->nAtoms(m_species[m_order[k]], m_element[m]);
277 }
278 }
279
280 // Do Gaussian elimination
281 for (size_t m = 0; m < nRows; m++) {
282 // Check for rows that are zero
283 bool isZeroRow = true;
284 for (size_t k = m; k < nColumns; k++) {
285 if (fabs(m_A(m,k)) > sqrt(Tiny)) {
286 isZeroRow = false;
287 break;
288 }
289 }
290 if (isZeroRow) {
291 // Find the last non-zero row
292 size_t n = nRows - 1;
293 bool foundSwapCandidate = false;
294 for (; n > m; n--) {
295 for (size_t k = m; k < nColumns; k++) {
296 if (fabs(m_A(n,k)) > sqrt(Tiny)) {
297 foundSwapCandidate = true;
298 break;
299 }
300 }
301 if (foundSwapCandidate) {
302 break;
303 }
304 }
305 if (m != n) {
306 // Swap this row with the last non-zero row
307 for (size_t k = 0; k < nColumns; k++) {
308 std::swap(m_A(n,k), m_A(m,k));
309 }
310 } else {
311 // All remaining rows are zero. Elimination is complete.
312 break;
313 }
314 }
315
316 // If a pivot is zero, exchange columns. This occurs when a species has
317 // an elemental composition that is not linearly independent of the
318 // component species that have already been assigned
319 if (m < nColumns && m_A(m,m) == 0.0) {
320 // First, we need to find a good candidate for a component species
321 // to swap in for the one that has zero pivot. It must contain
322 // element m, be linearly independent of the components processed so
323 // far (m_A(m,k) != 0), and should be a major species if possible.
324 // We'll choose the species with greatest mole fraction that
325 // satisfies these criteria.
326 doublereal maxmoles = -999.0;
327 size_t kmax = 0;
328 for (size_t k = m+1; k < nColumns; k++) {
329 if (m_A(m,k) != 0.0 && fabs(m_moles[m_order[k]]) > maxmoles) {
330 kmax = k;
331 maxmoles = fabs(m_moles[m_order[k]]);
332 }
333 }
334
335 // Now exchange the column with zero pivot with the
336 // column for this major species
337 for (size_t n = 0; n < nRows; n++) {
338 std::swap(m_A(n, m), m_A(n, kmax));
339 }
340
341 // exchange the species labels on the columns
342 std::swap(m_order[m], m_order[kmax]);
343 }
344
345 // scale row m so that the diagonal element is unity
346 double fctr = 1.0/m_A(m,m);
347 for (size_t k = 0; k < nColumns; k++) {
348 m_A(m,k) *= fctr;
349 }
350
351 // For all rows below the diagonal, subtract A(n,m)/A(m,m)
352 // * (row m) from row n, so that A(n,m) = 0.
353 for (size_t n = m+1; n < m_nel; n++) {
354 fctr = m_A(n,m)/m_A(m,m);
355 for (size_t k = 0; k < m_nsp; k++) {
356 m_A(n,k) -= m_A(m,k)*fctr;
357 }
358 }
359 }
360
361 // The left m_nel columns of A are now upper-diagonal. Now
362 // reduce the m_nel columns to diagonal form by back-solving
363 for (size_t m = std::min(nRows,nColumns)-1; m > 0; m--) {
364 for (size_t n = m-1; n != npos; n--) {
365 if (m_A(n,m) != 0.0) {
366 double fctr = m_A(n,m);
367 for (size_t k = m; k < m_nsp; k++) {
368 m_A(n,k) -= fctr*m_A(m,k);
369 }
370 }
371 }
372 }
373
374 // create stoichiometric coefficient matrix.
375 for (size_t n = 0; n < m_nsp; n++) {
376 if (n < m_nel) {
377 for (size_t k = 0; k < nFree(); k++) {
378 m_N(n, k) = -m_A(n, k + m_nel);
379 }
380 } else {
381 for (size_t k = 0; k < nFree(); k++) {
382 m_N(n, k) = 0.0;
383 }
384 m_N(n, n - m_nel) = 1.0;
385 }
386 }
387
388 // find reactions involving solution phase species
389 for (size_t j = 0; j < nFree(); j++) {
390 m_solnrxn[j] = false;
391 for (size_t k = 0; k < m_nsp; k++) {
392 if (m_N(k, j) != 0 && m_mix->solutionSpecies(m_species[m_order[k]])) {
393 m_solnrxn[j] = true;
394 }
395 }
396 }
397}
398
400{
401 m_work2 = x;
402 for (size_t k = 0; k < m_nsp; k++) {
403 x[m_order[k]] = m_work2[k];
404 }
405}
406
407void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN,
408 int loglevel)
409{
410 if (omega < 0.0) {
411 throw CanteraError("MultiPhaseEquil::step","negative omega");
412 }
413
414 for (size_t ik = 0; ik < m_nel; ik++) {
415 size_t k = m_order[ik];
416 m_lastmoles[k] = m_moles[k];
417 m_moles[k] += omega * deltaN[k];
418 }
419
420 for (size_t ik = m_nel; ik < m_nsp; ik++) {
421 size_t k = m_order[ik];
422 m_lastmoles[k] = m_moles[k];
423 if (m_majorsp[k]) {
424 m_moles[k] += omega * deltaN[k];
425 } else {
426 m_moles[k] = fabs(m_moles[k])*std::min(10.0,
427 exp(-m_deltaG_RT[ik - m_nel]));
428 }
429 }
430 updateMixMoles();
431}
432
433doublereal MultiPhaseEquil::stepComposition(int loglevel)
434{
435 m_iter++;
436 doublereal grad0 = computeReactionSteps(m_dxi);
437
438 // compute the mole fraction changes.
439 if (nFree()) {
440 multiply(m_N, m_dxi.data(), m_work.data());
441 }
442
443 // change to sequential form
444 unsort(m_work);
445
446 // scale omega to keep the major species non-negative
447 doublereal FCTR = 0.99;
448 const doublereal MAJOR_THRESHOLD = 1.0e-12;
449 double omegamax = 1.0;
450 for (size_t ik = 0; ik < m_nsp; ik++) {
451 size_t k = m_order[ik];
452 if (ik < m_nel) {
453 FCTR = 0.99;
454 if (m_moles[k] < MAJOR_THRESHOLD) {
455 m_force = true;
456 }
457 } else {
458 FCTR = 0.9;
459 }
460 // if species k is in a multi-species solution phase, then its mole
461 // number must remain positive, unless the entire phase goes away. First
462 // we'll determine an upper bound on omega, such that all
463 if (m_dsoln[k] == 1) {
464 if ((m_moles[k] > MAJOR_THRESHOLD) || (ik < m_nel)) {
465 if (m_moles[k] < MAJOR_THRESHOLD) {
466 m_force = true;
467 }
468 double omax = m_moles[k]*FCTR/(fabs(m_work[k]) + Tiny);
469 if (m_work[k] < 0.0 && omax < omegamax) {
470 omegamax = omax;
471 if (omegamax < 1.0e-5) {
472 m_force = true;
473 }
474 }
475 m_majorsp[k] = true;
476 } else {
477 m_majorsp[k] = false;
478 }
479 } else {
480 if (m_work[k] < 0.0 && m_moles[k] > 0.0) {
481 double omax = -m_moles[k]/m_work[k];
482 if (omax < omegamax) {
483 omegamax = omax;
484 if (omegamax < 1.0e-5) {
485 m_force = true;
486 }
487 }
488 }
489 m_majorsp[k] = true;
490 }
491 }
492
493 // now take a step with this scaled omega
494 step(omegamax, m_work);
495 // compute the gradient of G at this new position in the current direction.
496 // If it is positive, then we have overshot the minimum. In this case,
497 // interpolate back.
498 doublereal not_mu = 1.0e12;
499 m_mix->getValidChemPotentials(not_mu, m_mu.data());
500 doublereal grad1 = 0.0;
501 for (size_t k = 0; k < m_nsp; k++) {
502 grad1 += m_work[k] * m_mu[m_species[k]];
503 }
504
505 double omega = omegamax;
506 if (grad1 > 0.0) {
507 omega *= fabs(grad0) / (grad1 + fabs(grad0));
508 for (size_t k = 0; k < m_nsp; k++) {
509 m_moles[k] = m_lastmoles[k];
510 }
511 step(omega, m_work);
512 }
513 return omega;
514}
515
517{
518 vector_fp nu;
519 doublereal grad = 0.0;
520 dxi.resize(nFree());
521 computeN();
522 doublereal not_mu = 1.0e12;
523 m_mix->getValidChemPotentials(not_mu, m_mu.data());
524
525 for (size_t j = 0; j < nFree(); j++) {
526 // get stoichiometric vector
527 getStoichVector(j, nu);
528
529 // compute Delta G
530 doublereal dg_rt = 0.0;
531 for (size_t k = 0; k < m_nsp; k++) {
532 dg_rt += m_mu[m_species[k]] * nu[k];
533 }
534 dg_rt /= (m_temp * GasConstant);
535
536 m_deltaG_RT[j] = dg_rt;
537 double fctr = 1.0;
538
539 // if this is a formation reaction for a single-component phase,
540 // check whether reaction should be included
541 size_t k = m_order[j + m_nel];
542 if (!m_dsoln[k]) {
543 if (m_moles[k] <= 0.0 && dg_rt > 0.0) {
544 fctr = 0.0;
545 } else {
546 fctr = 0.5;
547 }
548 } else if (!m_solnrxn[j]) {
549 fctr = 1.0;
550 } else {
551 // component sum
552 double csum = 0.0;
553 for (k = 0; k < m_nel; k++) {
554 size_t kc = m_order[k];
555 double nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny;
556 csum += pow(nu[kc], 2)*m_dsoln[kc]/nmoles;
557 }
558
559 // noncomponent term
560 size_t kc = m_order[j + m_nel];
561 double nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny;
562 double term1 = m_dsoln[kc]/nmoles;
563
564 // sum over solution phases
565 doublereal sum = 0.0;
566 for (size_t ip = 0; ip < m_mix->nPhases(); ip++) {
567 ThermoPhase& p = m_mix->phase(ip);
568 if (p.nSpecies() > 1) {
569 double psum = 0.0;
570 for (k = 0; k < m_nsp; k++) {
571 kc = m_species[k];
572 if (m_mix->speciesPhaseIndex(kc) == ip) {
573 psum += pow(nu[k], 2);
574 }
575 }
576 sum -= psum / (fabs(m_mix->phaseMoles(ip)) + Tiny);
577 }
578 }
579 double rfctr = term1 + csum + sum;
580 if (fabs(rfctr) < Tiny) {
581 fctr = 1.0;
582 } else {
583 fctr = 1.0/(term1 + csum + sum);
584 }
585 }
586 dxi[j] = -fctr*dg_rt;
587
588 for (size_t m = 0; m < m_nel; m++) {
589 if (m_moles[m_order[m]] <= 0.0 && (m_N(m, j)*dxi[j] < 0.0)) {
590 dxi[j] = 0.0;
591 }
592 }
593 grad += dxi[j]*dg_rt;
594
595 }
596 return grad*GasConstant*m_temp;
597}
598
599void MultiPhaseEquil::computeN()
600{
601 // Sort the list of species by mole fraction (decreasing order)
602 std::vector<std::pair<double, size_t> > moleFractions(m_nsp);
603 for (size_t k = 0; k < m_nsp; k++) {
604 // use -Xk to generate reversed sort order
605 moleFractions[k].first = - m_mix->speciesMoles(m_species[k]);
606 moleFractions[k].second = k;
607 }
608 std::sort(moleFractions.begin(), moleFractions.end());
609 for (size_t k = 0; k < m_nsp; k++) {
610 m_sortindex[k] = moleFractions[k].second;
611 }
612
613 for (size_t m = 0; m < m_nel; m++) {
614 size_t k = 0;
615 for (size_t ik = 0; ik < m_nsp; ik++) {
616 k = m_sortindex[ik];
617 if (m_mix->nAtoms(m_species[k],m_element[m]) != 0) {
618 break;
619 }
620 }
621 bool ok = false;
622 for (size_t ij = 0; ij < m_nel; ij++) {
623 if (k == m_order[ij]) {
624 ok = true;
625 }
626 }
627 if (!ok || m_force) {
628 getComponents(m_sortindex);
629 m_force = true;
630 break;
631 }
632 }
633}
634
635doublereal MultiPhaseEquil::error()
636{
637 doublereal err, maxerr = 0.0;
638
639 // examine every reaction
640 for (size_t j = 0; j < nFree(); j++) {
641 size_t ik = j + m_nel;
642
643 // don't require formation reactions for solution species
644 // present in trace amounts to be equilibrated
645 if (!isStoichPhase(ik) && fabs(moles(ik)) <= SmallNumber) {
646 err = 0.0;
647 }
648
649 // for stoichiometric phase species, no error if not present and
650 // delta G for the formation reaction is positive
651 if (isStoichPhase(ik) && moles(ik) <= 0.0 &&
652 m_deltaG_RT[j] >= 0.0) {
653 err = 0.0;
654 } else {
655 err = fabs(m_deltaG_RT[j]);
656 }
657 maxerr = std::max(maxerr, err);
658 }
659 return maxerr;
660}
661
662double MultiPhaseEquil::phaseMoles(size_t iph) const
663{
664 return m_mix->phaseMoles(iph);
665}
666
667void MultiPhaseEquil::reportCSV(const std::string& reportFile)
668{
669 FILE* FP = fopen(reportFile.c_str(), "w");
670 if (!FP) {
671 throw CanteraError("MultiPhaseEquil::reportCSV", "Failure to open file");
672 }
673 vector_fp mf(m_nsp_mix, 1.0);
674 vector_fp fe(m_nsp_mix, 0.0);
675 vector_fp VolPM;
676 vector_fp activity;
677 vector_fp ac;
678 vector_fp mu;
679 vector_fp mu0;
680 vector_fp molalities;
681
682 double vol = 0.0;
683 for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
684 size_t istart = m_mix->speciesIndex(0, iphase);
685 ThermoPhase& tref = m_mix->phase(iphase);
686 size_t nSpecies = tref.nSpecies();
687 VolPM.resize(nSpecies, 0.0);
688 tref.getMoleFractions(&mf[istart]);
689 tref.getPartialMolarVolumes(VolPM.data());
690
691 double TMolesPhase = phaseMoles(iphase);
692 double VolPhaseVolumes = 0.0;
693 for (size_t k = 0; k < nSpecies; k++) {
694 VolPhaseVolumes += VolPM[k] * mf[istart + k];
695 }
696 VolPhaseVolumes *= TMolesPhase;
697 vol += VolPhaseVolumes;
698 }
699 fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
700 " -----------------------------\n");
701 fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature());
702 fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure());
703 fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
704
705 for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
706 size_t istart = m_mix->speciesIndex(0, iphase);
707 ThermoPhase& tref = m_mix->phase(iphase);
708 ThermoPhase* tp = &tref;
709 tp->getMoleFractions(&mf[istart]);
710 string phaseName = tref.name();
711 double TMolesPhase = phaseMoles(iphase);
712 size_t nSpecies = tref.nSpecies();
713 activity.resize(nSpecies, 0.0);
714 ac.resize(nSpecies, 0.0);
715 mu0.resize(nSpecies, 0.0);
716 mu.resize(nSpecies, 0.0);
717 VolPM.resize(nSpecies, 0.0);
718 molalities.resize(nSpecies, 0.0);
719 int actConvention = tp->activityConvention();
720 tp->getActivities(activity.data());
721 tp->getActivityCoefficients(ac.data());
722 tp->getStandardChemPotentials(mu0.data());
723 tp->getPartialMolarVolumes(VolPM.data());
724 tp->getChemPotentials(mu.data());
725 double VolPhaseVolumes = 0.0;
726 for (size_t k = 0; k < nSpecies; k++) {
727 VolPhaseVolumes += VolPM[k] * mf[istart + k];
728 }
729 VolPhaseVolumes *= TMolesPhase;
730 vol += VolPhaseVolumes;
731 if (actConvention == 1) {
732 MolalityVPSSTP* mTP = static_cast<MolalityVPSSTP*>(tp);
733 mTP->getMolalities(molalities.data());
734 tp->getChemPotentials(mu.data());
735
736 if (iphase == 0) {
737 fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
738 "Molalities, ActCoeff, Activity,"
739 "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
740
741 fprintf(FP," , , (kmol), , "
742 ", , ,"
743 " (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n");
744 }
745 for (size_t k = 0; k < nSpecies; k++) {
746 string sName = tp->speciesName(k);
747 fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e,"
748 "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
749 sName.c_str(),
750 phaseName.c_str(), TMolesPhase,
751 mf[istart + k], molalities[k], ac[k], activity[k],
752 mu0[k]*1.0E-6, mu[k]*1.0E-6,
753 mf[istart + k] * TMolesPhase,
754 VolPM[k], VolPhaseVolumes);
755 }
756 } else {
757 if (iphase == 0) {
758 fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
759 "Molalities, ActCoeff, Activity,"
760 " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
761
762 fprintf(FP," , , (kmol), , "
763 ", , ,"
764 " (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n");
765 }
766 for (size_t k = 0; k < nSpecies; k++) {
767 molalities[k] = 0.0;
768 }
769 for (size_t k = 0; k < nSpecies; k++) {
770 string sName = tp->speciesName(k);
771 fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, "
772 "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
773 sName.c_str(),
774 phaseName.c_str(), TMolesPhase,
775 mf[istart + k], molalities[k], ac[k],
776 activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6,
777 mf[istart + k] * TMolesPhase,
778 VolPM[k], VolPhaseVolumes);
779 }
780 }
781 }
782 fclose(FP);
783}
784
785}
Header for intermediate ThermoPhase object for phases which employ molality based activity coefficien...
vector_fp & data()
Return a reference to the data vector.
Definition: Array.h:218
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:70
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...
doublereal computeReactionSteps(vector_fp &dxi)
Compute the change in extent of reaction for each reaction.
void unsort(vector_fp &x)
Re-arrange a vector of species properties in sorted form (components first) into unsorted,...
void finish()
Clean up the composition.
size_t nFree() const
Number of degrees of freedom.
doublereal stepComposition(int loglevel=0)
Take one step in composition, given the gradient of G at the starting point, and a vector of reaction...
MultiPhaseEquil(MultiPhase *mix, bool start=true, int loglevel=0)
Construct a multiphase equilibrium manager for a multiphase mixture.
int setInitialMoles(int loglevel=0)
Estimate the initial mole numbers.
A class for multiphase mixtures.
Definition: MultiPhase.h:60
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...
Definition: MultiPhase.h:229
bool solutionSpecies(size_t kGlob) const
Return true is species kGlob is a species in a multicomponent solution phase.
Definition: MultiPhase.cpp:273
void setMoles(const doublereal *n)
Sets all of the global species mole numbers.
Definition: MultiPhase.cpp:387
size_t nSpecies() const
Number of species, summed over all phases.
Definition: MultiPhase.h:126
doublereal pressure() const
Pressure [Pa].
Definition: MultiPhase.h:391
doublereal speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
Definition: MultiPhase.cpp:186
size_t nPhases() const
Number of phases.
Definition: MultiPhase.h:428
doublereal phaseMoles(const size_t n) const
Return the number of moles in phase n.
Definition: MultiPhase.cpp:791
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
Definition: MultiPhase.cpp:761
bool tempOK(size_t p) const
Return true if the phase p has valid thermo data for the current temperature.
Definition: MultiPhase.cpp:811
std::string elementName(size_t m) const
Returns the name of the global element m.
Definition: MultiPhase.cpp:732
doublereal nAtoms(const size_t kGlob, const size_t mGlob) const
Returns the Number of atoms of global element mGlob in global species kGlob.
Definition: MultiPhase.cpp:766
size_t speciesPhaseIndex(const size_t kGlob) const
Returns the phase index of the Kth "global" species.
Definition: MultiPhase.cpp:801
doublereal elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
Definition: MultiPhase.cpp:192
void getValidChemPotentials(doublereal not_mu, doublereal *mu, bool standard=false) const
Returns a vector of Valid chemical potentials.
Definition: MultiPhase.cpp:253
doublereal temperature() const
Temperature [K].
Definition: MultiPhase.h:331
size_t nElements() const
Number of elements.
Definition: MultiPhase.h:100
ThermoPhase & phase(size_t n)
Return a reference to phase n.
Definition: MultiPhase.cpp:161
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:273
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:543
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:102
Namespace for the Cantera kernel.
Definition: AnyMap.h:29
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:192
const double Tiny
Small number to compare differences of mole fractions against.
Definition: ct_defs.h:170
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.
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:153
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:184
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:113
Contains declarations for string manipulation functions within Cantera.