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