Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vcs_MultiPhaseEquil.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_MultiPhaseEquil.cpp
3  * Driver routine for the VCSnonideal equilibrium solver package
4  */
5 /*
6  * Copyright (2006) Sandia Corporation. Under the terms of
7  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
8  * U.S. Government retains certain rights in this software.
9  */
10 
14 
15 #include "cantera/base/clockWC.h"
20 
21 #include <cstdio>
22 
23 using namespace std;
24 
25 namespace Cantera
26 {
27 vcs_MultiPhaseEquil::vcs_MultiPhaseEquil() :
28  m_vprob(0, 0, 0),
29  m_mix(0),
30  m_printLvl(0)
31 {
32 }
33 
34 vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(MultiPhase* mix, int printLvl) :
35  m_vprob(mix->nSpecies(), mix->nElements(), mix->nPhases()),
36  m_mix(0),
37  m_printLvl(printLvl)
38 {
39  m_mix = mix;
41  /*
42  * Work out the details of the VCS_VPROB construction and
43  * Transfer the current problem to VCS_PROB object
44  */
45  int res = vcs_Cantera_to_vprob(mix, &m_vprob);
46  if (res != 0) {
47  plogf("problems\n");
48  }
49 }
50 
51 int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
52  int estimateEquil,
53  int printLvl, doublereal err,
54  int maxsteps, int loglevel)
55 {
56  doublereal Vtarget = m_mix->volume();
57  if ((XY != TV) && (XY != HV) && (XY != UV) && (XY != SV)) {
58  throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TV",
59  "Wrong XY flag:" + int2str(XY));
60  }
61  int maxiter = 100;
62  int iSuccess = 0;
63  if (XY == TV) {
64  m_mix->setTemperature(xtarget);
65  }
66  int strt = estimateEquil;
67  double P1 = 0.0;
68  double V1 = 0.0;
69  double V2 = 0.0;
70  double P2 = 0.0;
71  doublereal Tlow = 0.5 * m_mix->minTemp();
72  doublereal Thigh = 2.0 * m_mix->maxTemp();
73  doublereal Vnow, Verr;
74  int printLvlSub = std::max(0, printLvl - 1);
75  for (int n = 0; n < maxiter; n++) {
76  double Pnow = m_mix->pressure();
77 
78  switch (XY) {
79  case TV:
80  iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
81  break;
82  case HV:
83  iSuccess = equilibrate_HP(xtarget, HP, Tlow, Thigh, strt,
84  printLvlSub, err, maxsteps, loglevel);
85  break;
86  case UV:
87  iSuccess = equilibrate_HP(xtarget, UP, Tlow, Thigh, strt,
88  printLvlSub, err, maxsteps, loglevel);
89  break;
90  case SV:
91  iSuccess = equilibrate_SP(xtarget, Tlow, Thigh, strt,
92  printLvlSub, err, maxsteps, loglevel);
93  break;
94  default:
95  break;
96  }
97  strt = false;
98  Vnow = m_mix->volume();
99  if (n == 0) {
100  V2 = Vnow;
101  P2 = Pnow;
102  } else if (n == 1) {
103  V1 = Vnow;
104  P1 = Pnow;
105  } else {
106  P2 = P1;
107  V2 = V1;
108  P1 = Pnow;
109  V1 = Vnow;
110  }
111 
112  Verr = fabs((Vtarget - Vnow)/Vtarget);
113 
114  if (Verr < err) {
115  goto done;
116  }
117  double Pnew;
118  // find dV/dP
119  if (n > 1) {
120  double dVdP = (V2 - V1) / (P2 - P1);
121  if (dVdP == 0.0) {
122  throw CanteraError("vcs_MultiPhase::equilibrate_TV",
123  "dVdP == 0.0");
124  } else {
125  Pnew = Pnow + (Vtarget - Vnow) / dVdP;
126  if (Pnew < 0.2 * Pnow) {
127  Pnew = 0.2 * Pnow;
128  }
129  if (Pnew > 3.0 * Pnow) {
130  Pnew = 3.0 * Pnow;
131  }
132  }
133 
134  } else {
135  m_mix->setPressure(Pnow*1.01);
136  double dVdP = (m_mix->volume() - Vnow)/(0.01*Pnow);
137  Pnew = Pnow + 0.5*(Vtarget - Vnow)/dVdP;
138  if (Pnew < 0.5* Pnow) {
139  Pnew = 0.5 * Pnow;
140  }
141  if (Pnew > 1.7 * Pnow) {
142  Pnew = 1.7 * Pnow;
143  }
144 
145  }
146  m_mix->setPressure(Pnew);
147  }
148  throw CanteraError("vcs_MultiPhase::equilibrate_TV",
149  "No convergence for V");
150 
151 done:
152  ;
153  return iSuccess;
154 }
155 
156 int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
157  int XY, double Tlow, double Thigh,
158  int estimateEquil,
159  int printLvl, doublereal err,
160  int maxsteps, int loglevel)
161 {
162  int maxiter = 100;
163  int iSuccess;
164  if (XY != HP && XY != UP) {
165  throw CanteraError("vcs_MultiPhaseEquil::equilibrate_HP",
166  "Wrong XP" + int2str(XY));
167  }
168  int strt = estimateEquil;
169 
170  // Lower bound on T. This will change as we progress in the calculation
171  if (Tlow <= 0.0) {
172  Tlow = 0.5 * m_mix->minTemp();
173  }
174  // Upper bound on T. This will change as we progress in the calculation
175  if (Thigh <= 0.0 || Thigh > 1.0E6) {
176  Thigh = 2.0 * m_mix->maxTemp();
177  }
178 
179  doublereal cpb = 1.0, Tnew;
180  doublereal Hlow = Undef;
181  doublereal Hhigh = Undef;
182  doublereal Tnow = m_mix->temperature();
183  int printLvlSub = std::max(printLvl - 1, 0);
184 
185  for (int n = 0; n < maxiter; n++) {
186  // start with a loose error tolerance, but tighten it as we get
187  // close to the final temperature
188  try {
189  Tnow = m_mix->temperature();
190  iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
191  strt = 0;
192  double Hnow = (XY == UP) ? m_mix->IntEnergy() : m_mix->enthalpy();
193  double pmoles[10];
194  pmoles[0] = m_mix->phaseMoles(0);
195  double Tmoles = pmoles[0];
196  double HperMole = Hnow/Tmoles;
197  if (printLvl > 0) {
198  plogf("T = %g, Hnow = %g ,Tmoles = %g, HperMole = %g",
199  Tnow, Hnow, Tmoles, HperMole);
200  plogendl();
201  }
202 
203  // the equilibrium enthalpy monotonically increases with T;
204  // if the current value is below the target, then we know the
205  // current temperature is too low. Set the lower bounds.
206 
207 
208  if (Hnow < Htarget) {
209  if (Tnow > Tlow) {
210  Tlow = Tnow;
211  Hlow = Hnow;
212  }
213  }
214  // the current enthalpy is greater than the target; therefore the
215  // current temperature is too high. Set the high bounds.
216  else {
217  if (Tnow < Thigh) {
218  Thigh = Tnow;
219  Hhigh = Hnow;
220  }
221  }
222  double dT, dTa, dTmax, Tnew;
223  if (Hlow != Undef && Hhigh != Undef) {
224  cpb = (Hhigh - Hlow)/(Thigh - Tlow);
225  dT = (Htarget - Hnow)/cpb;
226  dTa = fabs(dT);
227  dTmax = 0.5*fabs(Thigh - Tlow);
228  if (dTa > dTmax) {
229  dT *= dTmax/dTa;
230  }
231  } else {
232  Tnew = sqrt(Tlow*Thigh);
233  dT = clip(Tnew - Tnow, -200.0, 200.0);
234  }
235  double acpb = std::max(fabs(cpb), 1.0E-6);
236  double denom = std::max(fabs(Htarget), acpb);
237  double Herr = Htarget - Hnow;
238  double HConvErr = fabs((Herr)/denom);
239  if (printLvl > 0) {
240  plogf(" equilibrate_HP: It = %d, Tcurr = %g Hcurr = %g, Htarget = %g\n",
241  n, Tnow, Hnow, Htarget);
242  plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
243  Herr, cpb, HConvErr);
244  }
245 
246  if (HConvErr < err) { // || dTa < 1.0e-4) {
247  if (printLvl > 0) {
248  plogf(" equilibrate_HP: CONVERGENCE: Hfinal = %g Tfinal = %g, Its = %d \n",
249  Hnow, Tnow, n);
250  plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
251  Herr, cpb, HConvErr);
252  }
253  goto done;
254  }
255  Tnew = Tnow + dT;
256  if (Tnew < 0.0) {
257  Tnew = 0.5*Tnow;
258  }
259  m_mix->setTemperature(Tnew);
260 
261  } catch (CanteraError err) {
262  if (!estimateEquil) {
263  strt = -1;
264  } else {
265  Tnew = 0.5*(Tnow + Thigh);
266  if (fabs(Tnew - Tnow) < 1.0) {
267  Tnew = Tnow + 1.0;
268  }
269  m_mix->setTemperature(Tnew);
270  }
271  }
272 
273  }
274  throw CanteraError("MultiPhase::equilibrate_HP",
275  "No convergence for T");
276 done:
277  ;
278  return iSuccess;
279 }
280 
281 int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
282  double Tlow, double Thigh,
283  int estimateEquil,
284  int printLvl, doublereal err,
285  int maxsteps, int loglevel)
286 {
287  int maxiter = 100;
288  int strt = estimateEquil;
289 
290  // Lower bound on T. This will change as we progress in the calculation
291  if (Tlow <= 0.0) {
292  Tlow = 0.5 * m_mix->minTemp();
293  }
294  // Upper bound on T. This will change as we progress in the calculation
295  if (Thigh <= 0.0 || Thigh > 1.0E6) {
296  Thigh = 2.0 * m_mix->maxTemp();
297  }
298 
299  doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
300  doublereal Slow = Undef;
301  doublereal Shigh = Undef;
302  doublereal Tnow = m_mix->temperature();
303  Tlow = std::min(Tnow, Tlow);
304  Thigh = std::max(Tnow, Thigh);
305  int printLvlSub = std::max(printLvl - 1, 0);
306 
307  for (int n = 0; n < maxiter; n++) {
308 
309  // start with a loose error tolerance, but tighten it as we get
310  // close to the final temperature
311  try {
312  Tnow = m_mix->temperature();
313  int iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
314  strt = 0;
315  double Snow = m_mix->entropy();
316  double pmoles[10];
317  pmoles[0] = m_mix->phaseMoles(0);
318  double Tmoles = pmoles[0];
319  double SperMole = Snow/Tmoles;
320  plogf("T = %g, Snow = %g ,Tmoles = %g, SperMole = %g\n",
321  Tnow, Snow, Tmoles, SperMole);
322 
323  // the equilibrium entropy monotonically increases with T;
324  // if the current value is below the target, then we know the
325  // current temperature is too low. Set the lower bounds to the
326  // current condition.
327  if (Snow < Starget) {
328  if (Tnow > Tlow) {
329  Tlow = Tnow;
330  Slow = Snow;
331  } else {
332  if (Slow > Starget) {
333  if (Snow < Slow) {
334  Thigh = Tlow;
335  Shigh = Slow;
336  Tlow = Tnow;
337  Slow = Snow;
338  }
339  }
340  }
341  }
342  // the current enthalpy is greater than the target; therefore the
343  // current temperature is too high. Set the high bounds.
344  else {
345  if (Tnow < Thigh) {
346  Thigh = Tnow;
347  Shigh = Snow;
348  }
349  }
350  if (Slow != Undef && Shigh != Undef) {
351  cpb = (Shigh - Slow)/(Thigh - Tlow);
352  dT = (Starget - Snow)/cpb;
353  Tnew = Tnow + dT;
354  dTa = fabs(dT);
355  dTmax = 0.5*fabs(Thigh - Tlow);
356  if (Tnew > Thigh || Tnew < Tlow) {
357  dTmax = 1.5*fabs(Thigh - Tlow);
358  }
359  dTmax = std::min(dTmax, 300.);
360  if (dTa > dTmax) {
361  dT *= dTmax/dTa;
362  }
363  } else {
364  Tnew = sqrt(Tlow*Thigh);
365  dT = Tnew - Tnow;
366  }
367 
368  double acpb = std::max(fabs(cpb), 1.0E-6);
369  double denom = std::max(fabs(Starget), acpb);
370  double Serr = Starget - Snow;
371  double SConvErr = fabs((Serr)/denom);
372  if (printLvl > 0) {
373  plogf(" equilibrate_SP: It = %d, Tcurr = %g Scurr = %g, Starget = %g\n",
374  n, Tnow, Snow, Starget);
375  plogf(" S rel error = %g, cp = %g, SConvErr = %g\n",
376  Serr, cpb, SConvErr);
377  }
378 
379  if (SConvErr < err) { // || dTa < 1.0e-4) {
380  if (printLvl > 0) {
381  plogf(" equilibrate_SP: CONVERGENCE: Sfinal = %g Tfinal = %g, Its = %d \n",
382  Snow, Tnow, n);
383  plogf(" S rel error = %g, cp = %g, HConvErr = %g\n",
384  Serr, cpb, SConvErr);
385  }
386  return iSuccess;
387  }
388  Tnew = Tnow + dT;
389  if (Tnew < 0.0) {
390  Tnew = 0.5*Tnow;
391  }
392  m_mix->setTemperature(Tnew);
393 
394  } catch (CanteraError err) {
395  if (!estimateEquil) {
396  strt = -1;
397  } else {
398  Tnew = 0.5*(Tnow + Thigh);
399  if (fabs(Tnew - Tnow) < 1.0) {
400  Tnew = Tnow + 1.0;
401  }
402  m_mix->setTemperature(Tnew);
403  }
404  }
405 
406  }
407  throw CanteraError("MultiPhase::equilibrate_SP",
408  "No convergence for T");
409 }
410 
411 int vcs_MultiPhaseEquil::equilibrate(int XY, int estimateEquil,
412  int printLvl, doublereal err,
413  int maxsteps, int loglevel)
414 {
415  doublereal xtarget;
416  if (XY == TP) {
417  return equilibrate_TP(estimateEquil, printLvl, err, maxsteps, loglevel);
418  } else if (XY == HP || XY == UP) {
419  if (XY == HP) {
420  xtarget = m_mix->enthalpy();
421  } else {
422  xtarget = m_mix->IntEnergy();
423  }
424  double Tlow = 0.5 * m_mix->minTemp();
425  double Thigh = 2.0 * m_mix->maxTemp();
426  return equilibrate_HP(xtarget, XY, Tlow, Thigh,
427  estimateEquil, printLvl, err, maxsteps, loglevel);
428  } else if (XY == SP) {
429  xtarget = m_mix->entropy();
430  double Tlow = 0.5 * m_mix->minTemp();
431  double Thigh = 2.0 * m_mix->maxTemp();
432  return equilibrate_SP(xtarget, Tlow, Thigh,
433  estimateEquil, printLvl, err, maxsteps, loglevel);
434  } else if (XY == TV) {
435  xtarget = m_mix->temperature();
436  return equilibrate_TV(XY, xtarget,
437  estimateEquil, printLvl, err, maxsteps, loglevel);
438  } else if (XY == HV) {
439  xtarget = m_mix->enthalpy();
440  return equilibrate_TV(XY, xtarget,
441  estimateEquil, printLvl, err, maxsteps, loglevel);
442  } else if (XY == UV) {
443  xtarget = m_mix->IntEnergy();
444  return equilibrate_TV(XY, xtarget,
445  estimateEquil, printLvl, err, maxsteps, loglevel);
446  } else if (XY == SV) {
447  xtarget = m_mix->entropy();
448  return equilibrate_TV(XY, xtarget, estimateEquil,
449  printLvl, err, maxsteps, loglevel);
450  } else {
451  throw CanteraError(" vcs_MultiPhaseEquil::equilibrate",
452  "Unsupported Option");
453  }
454 }
455 
457  int printLvl, doublereal err,
458  int maxsteps, int loglevel)
459 {
460  int maxit = maxsteps;
461  clockWC tickTock;
462 
463  m_printLvl = printLvl;
464  m_vprob.m_printLvl = printLvl;
465 
466 
467  /*
468  * Extract the current state information
469  * from the MultiPhase object and
470  * Transfer it to VCS_PROB object.
471  */
473  if (res != 0) {
474  plogf("problems\n");
475  }
476 
477 
478  // Set the estimation technique
479  if (estimateEquil) {
480  m_vprob.iest = estimateEquil;
481  } else {
482  m_vprob.iest = 0;
483  }
484 
485  // Check obvious bounds on the temperature and pressure
486  // NOTE, we may want to do more here with the real bounds
487  // given by the ThermoPhase objects.
488  double T = m_mix->temperature();
489  if (T <= 0.0) {
490  throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
491  "Temperature less than zero on input");
492  }
493  double pres = m_mix->pressure();
494  if (pres <= 0.0) {
495  throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
496  "Pressure less than zero on input");
497  }
498 
499  /*
500  * Print out the problem specification from the point of
501  * view of the vprob object.
502  */
504 
505  /*
506  * Call the thermo Program
507  */
508  int ip1 = m_printLvl;
509  int ipr = std::max(0, m_printLvl-1);
510  if (m_printLvl >= 3) {
511  ip1 = m_printLvl - 2;
512  } else {
513  ip1 = 0;
514  }
515  int iSuccess = m_vsolve.vcs(&m_vprob, 0, ipr, ip1, maxit);
516 
517  /*
518  * Transfer the information back to the MultiPhase object.
519  * Note we don't just call setMoles, because some multispecies
520  * solution phases may be zeroed out, and that would cause a problem
521  * for that routine. Also, the mole fractions of such zeroed out
522  * phases actually contain information about likely reemergent
523  * states.
524  */
526  size_t kGlob = 0;
527  for (size_t ip = 0; ip < m_vprob.NPhase; ip++) {
528  double phaseMole = 0.0;
529  ThermoPhase& tref = m_mix->phase(ip);
530  for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) {
531  phaseMole += m_vprob.w[kGlob];
532  }
533  m_mix->setPhaseMoles(ip, phaseMole);
534  }
535 
536  double te = tickTock.secondsWC();
537  if (printLvl > 0) {
538  plogf("\n Results from vcs:\n");
539  if (iSuccess != 0) {
540  plogf("\nVCS FAILED TO CONVERGE!\n");
541  }
542  plogf("\n");
543  plogf("Temperature = %g Kelvin\n", m_vprob.T);
544  plogf("Pressure = %g Pa\n", m_vprob.PresPA);
545  plogf("\n");
546  plogf("----------------------------------------"
547  "---------------------\n");
548  plogf(" Name Mole_Number");
549  if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
550  plogf("(kmol)");
551  } else {
552  plogf("(gmol)");
553  }
554  plogf(" Mole_Fraction Chem_Potential");
555  if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
556  plogf(" (kcal/mol)\n");
557  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
558  plogf(" (Dimensionless)\n");
559  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
560  plogf(" (kJ/mol)\n");
561  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
562  plogf(" (Kelvin)\n");
563  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
564  plogf(" (J/kmol)\n");
565  }
566  plogf("--------------------------------------------------"
567  "-----------\n");
568  for (size_t i = 0; i < m_vprob.nspecies; i++) {
569  plogf("%-12s", m_vprob.SpName[i].c_str());
571  plogf(" %15.3e %15.3e ", 0.0, m_vprob.mf[i]);
572  plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]);
573  } else {
574  plogf(" %15.3e %15.3e ", m_vprob.w[i], m_vprob.mf[i]);
575  if (m_vprob.w[i] <= 0.0) {
576  size_t iph = m_vprob.PhaseID[i];
577  vcs_VolPhase* VPhase = m_vprob.VPhaseList[iph];
578  if (VPhase->nSpecies() > 1) {
579  plogf(" -1.000e+300\n");
580  } else {
581  plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]);
582  }
583  } else {
584  plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]);
585  }
586  }
587  }
588  plogf("------------------------------------------"
589  "-------------------\n");
590  if (printLvl > 2) {
591  if (m_vsolve.m_timing_print_lvl > 0) {
592  plogf("Total time = %12.6e seconds\n", te);
593  }
594  }
595  }
596  return iSuccess;
597 }
598 
599 void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
600 {
601  double vol = 0.0;
602  size_t nphase = m_vprob.NPhase;
603 
604  FILE* FP = fopen(reportFile.c_str(), "w");
605  if (!FP) {
606  throw CanteraError("vcs_MultiPhaseEquil::reportCSV",
607  "Failure to open file");
608  }
609  double Temp = m_mix->temperature();
610  double pres = m_mix->pressure();
611  double* mf = VCS_DATA_PTR(m_vprob.mf);
612 #ifdef DEBUG_MODE
613  double* fe = VCS_DATA_PTR(m_vprob.m_gibbsSpecies);
614 #endif
615  std::vector<double> VolPM;
616  std::vector<double> activity;
617  std::vector<double> ac;
618  std::vector<double> mu;
619  std::vector<double> mu0;
620  std::vector<double> molalities;
621 
622  vol = 0.0;
623  for (size_t iphase = 0; iphase < nphase; iphase++) {
624  size_t istart = m_mix->speciesIndex(0, iphase);
625  ThermoPhase& tref = m_mix->phase(iphase);
626  size_t nSpecies = tref.nSpecies();
627  VolPM.resize(nSpecies, 0.0);
629  vcs_VolPhase* volP = m_vprob.VPhaseList[iphase];
630 
631  double TMolesPhase = volP->totalMoles();
632  double VolPhaseVolumes = 0.0;
633  for (size_t k = 0; k < nSpecies; k++) {
634  VolPhaseVolumes += VolPM[k] * mf[istart + k];
635  }
636  VolPhaseVolumes *= TMolesPhase;
637  vol += VolPhaseVolumes;
638  }
639 
640  fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
641  " -----------------------------\n");
642  fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
643  fprintf(FP,"Pressure = %11.5g Pascal\n", pres);
644  fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
645  fprintf(FP,"Number Basis optimizations = %d\n", m_vprob.m_NumBasisOptimizations);
646  fprintf(FP,"Number VCS iterations = %d\n", m_vprob.m_Iterations);
647 
648  for (size_t iphase = 0; iphase < nphase; iphase++) {
649  size_t istart = m_mix->speciesIndex(0, iphase);
650  ThermoPhase& tref = m_mix->phase(iphase);
651  string phaseName = tref.name();
652  vcs_VolPhase* volP = m_vprob.VPhaseList[iphase];
653  double TMolesPhase = volP->totalMoles();
654  //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase));
655  size_t nSpecies = tref.nSpecies();
656  activity.resize(nSpecies, 0.0);
657  ac.resize(nSpecies, 0.0);
658 
659  mu0.resize(nSpecies, 0.0);
660  mu.resize(nSpecies, 0.0);
661  VolPM.resize(nSpecies, 0.0);
662  molalities.resize(nSpecies, 0.0);
663 
664  int actConvention = tref.activityConvention();
665  tref.getActivities(VCS_DATA_PTR(activity));
668 
671  double VolPhaseVolumes = 0.0;
672  for (size_t k = 0; k < nSpecies; k++) {
673  VolPhaseVolumes += VolPM[k] * mf[istart + k];
674  }
675  VolPhaseVolumes *= TMolesPhase;
676  vol += VolPhaseVolumes;
677 
678 
679  if (actConvention == 1) {
680  MolalityVPSSTP* mTP = static_cast<MolalityVPSSTP*>(&tref);
681  mTP->getMolalities(VCS_DATA_PTR(molalities));
683 
684  if (iphase == 0) {
685  fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
686  "Molalities, ActCoeff, Activity,"
687  "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
688 
689  fprintf(FP," , , (kmol), , "
690  " , , ,"
691  " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
692  }
693  for (size_t k = 0; k < nSpecies; k++) {
694  std::string sName = tref.speciesName(k);
695  fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e,"
696  "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
697  sName.c_str(),
698  phaseName.c_str(), TMolesPhase,
699  mf[istart + k], molalities[k], ac[k], activity[k],
700  mu0[k]*1.0E-6, mu[k]*1.0E-6,
701  mf[istart + k] * TMolesPhase,
702  VolPM[k], VolPhaseVolumes);
703  }
704 
705  } else {
706  if (iphase == 0) {
707  fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
708  "Molalities, ActCoeff, Activity,"
709  " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
710 
711  fprintf(FP," , , (kmol), , "
712  " , , ,"
713  " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
714  }
715  for (size_t k = 0; k < nSpecies; k++) {
716  molalities[k] = 0.0;
717  }
718  for (size_t k = 0; k < nSpecies; k++) {
719  std::string sName = tref.speciesName(k);
720  fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, "
721  "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
722  sName.c_str(),
723  phaseName.c_str(), TMolesPhase,
724  mf[istart + k], molalities[k], ac[k],
725  activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6,
726  mf[istart + k] * TMolesPhase,
727  VolPM[k], VolPhaseVolumes);
728  }
729  }
730 
731 #ifdef DEBUG_MODE
732  /*
733  * Check consistency: These should be equal
734  */
735  tref.getChemPotentials(fe+istart);
736  for (size_t k = 0; k < nSpecies; k++) {
737  if (!vcs_doubleEqual(fe[istart+k], mu[k])) {
738  fprintf(FP,"ERROR: incompatibility!\n");
739  fclose(FP);
740  throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "incompatibility!");
741  }
742  }
743 #endif
744 
745  }
746  fclose(FP);
747 }
748 
749 /*
750  * HKM -> Work on transferring the current value of the voltages into the
751  * equilibrium problem.
752  */
754 {
755  VCS_SPECIES_THERMO* ts_ptr = 0;
756 
757  /*
758  * Calculate the total number of species and phases in the problem
759  */
760  size_t totNumPhases = mphase->nPhases();
761  size_t totNumSpecies = mphase->nSpecies();
762 
763  // Problem type has yet to be worked out.
764  vprob->prob_type = 0;
765  vprob->nspecies = totNumSpecies;
766  vprob->ne = 0;
767  vprob->NPhase = totNumPhases;
768  vprob->m_VCS_UnitsFormat = VCS_UNITS_MKS;
769  // Set the initial estimate to a machine generated estimate for now
770  // We will work out the details later.
771  vprob->iest = -1;
772  vprob->T = mphase->temperature();
773  vprob->PresPA = mphase->pressure();
774  vprob->Vol = mphase->volume();
775  vprob->Title = "MultiPhase Object";
776 
777  int printLvl = vprob->m_printLvl;
778 
779  /*
780  * Loop over the phases, transferring pertinent information
781  */
782  int kT = 0;
783  for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
784 
785  /*
786  * Get the ThermoPhase object - assume volume phase
787  */
788  ThermoPhase* tPhase = &(mphase->phase(iphase));
789  size_t nelem = tPhase->nElements();
790 
791  /*
792  * Query Cantera for the equation of state type of the
793  * current phase.
794  */
795  int eos = tPhase->eosType();
796  bool gasPhase = (eos == cIdealGas);
797 
798  /*
799  * Find out the number of species in the phase
800  */
801  size_t nSpPhase = tPhase->nSpecies();
802  /*
803  * Find out the name of the phase
804  */
805  string phaseName = tPhase->name();
806 
807  /*
808  * Call the basic vcs_VolPhase creation routine.
809  * Properties set here:
810  * ->PhaseNum = phase number in the thermo problem
811  * ->GasPhase = Boolean indicating whether it is a gas phase
812  * ->NumSpecies = number of species in the phase
813  * ->TMolesInert = Inerts in the phase = 0.0 for cantera
814  * ->PhaseName = Name of the phase
815  */
816 
817  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
818  VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
819  VolPhase->m_gasPhase = gasPhase;
820  /*
821  * Tell the vcs_VolPhase pointer about cantera
822  */
823  VolPhase->p_VCS_UnitsFormat = vprob->m_VCS_UnitsFormat;
824  VolPhase->setPtrThermoPhase(tPhase);
825  VolPhase->setTotalMoles(0.0);
826  /*
827  * Set the electric potential of the volume phase from the
828  * ThermoPhase object's value.
829  */
830  VolPhase->setElectricPotential(tPhase->electricPotential());
831  /*
832  * Query the ThermoPhase object to find out what convention
833  * it uses for the specification of activity and Standard State.
834  */
835  VolPhase->p_activityConvention = tPhase->activityConvention();
836  /*
837  * Assign the value of eqn of state
838  * -> Handle conflicts here.
839  */
840  switch (eos) {
841  case cIdealGas:
842  VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS;
843  break;
844  case cIncompressible:
845  VolPhase->m_eqnState = VCS_EOS_CONSTANT;
846  break;
847  case cSurf:
848  throw CanteraError("VCSnonideal", "cSurf not handled yet.");
849  case cStoichSubstance:
850  VolPhase->m_eqnState = VCS_EOS_STOICH_SUB;
851  break;
852  case cPureFluid:
853  if (printLvl > 1) {
854  plogf("cPureFluid not recognized yet by VCSnonideal\n");
855  }
856  break;
857  case cEdge:
858  throw CanteraError("VCSnonideal", "cEdge not handled yet.");
859  case cIdealSolidSolnPhase0:
860  case cIdealSolidSolnPhase1:
861  case cIdealSolidSolnPhase2:
862  VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN;
863  break;
864  default:
865  if (printLvl > 1) {
866  plogf("Unknown Cantera EOS to VCSnonideal: %d\n", eos);
867  }
868  VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA;
869  if (!VolPhase->usingCanteraCalls()) {
870  throw CanteraError("vcs_Cantera_to_vprob",
871  "vcs functions asked for, but unimplemented");
872  }
873  break;
874  }
875 
876  /*
877  * Transfer all of the element information from the
878  * ThermoPhase object to the vcs_VolPhase object.
879  * Also decide whether we need a new charge neutrality
880  * element in the phase to enforce a charge neutrality
881  * constraint.
882  * We also decide whether this is a single species phase
883  * with the voltage being the independent variable setting
884  * the chemical potential of the electrons.
885  */
886  VolPhase->transferElementsFM(tPhase);
887 
888  /*
889  * Combine the element information in the vcs_VolPhase
890  * object into the vprob object.
891  */
892  vprob->addPhaseElements(VolPhase);
893 
894  VolPhase->setState_TP(vprob->T, vprob->PresPA);
895  vector<double> muPhase(tPhase->nSpecies(),0.0);
896  tPhase->getChemPotentials(&muPhase[0]);
897  double tMoles = 0.0;
898  /*
899  * Loop through each species in the current phase
900  */
901  for (size_t k = 0; k < nSpPhase; k++) {
902  /*
903  * Obtain the molecular weight of the species from the
904  * ThermoPhase object
905  */
906  vprob->WtSpecies[kT] = tPhase->molecularWeight(k);
907 
908  /*
909  * Obtain the charges of the species from the
910  * ThermoPhase object
911  */
912  vprob->Charge[kT] = tPhase->charge(k);
913 
914  /*
915  * Set the phaseid of the species
916  */
917  vprob->PhaseID[kT] = iphase;
918 
919  /*
920  * Transfer the Species name
921  */
922  string stmp = mphase->speciesName(kT);
923  vprob->SpName[kT] = stmp;
924  /*
925  * Transfer the type of unknown
926  */
927  vprob->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k);
928 
929  if (vprob->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) {
930  /*
931  * Set the initial number of kmoles of the species
932  * and the mole fraction vector
933  */
934  vprob->w[kT] = mphase->speciesMoles(kT);
935  tMoles += vprob->w[kT];
936  vprob->mf[kT] = mphase->moleFraction(kT);
937  } else if (vprob->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
938  vprob->w[kT] = tPhase->electricPotential();
939  vprob->mf[kT] = mphase->moleFraction(kT);
940  } else {
941  throw CanteraError(" vcs_Cantera_to_vprob() ERROR",
942  "Unknown species type: " + int2str(vprob->SpeciesUnknownType[kT]));
943  }
944 
945 
946  /*
947  * transfer chemical potential vector
948  */
949  vprob->m_gibbsSpecies[kT] = muPhase[k];
950  /*
951  * Transfer the species information from the
952  * volPhase structure to the VPROB structure
953  * This includes:
954  * FormulaMatrix[][]
955  * VolPhase->IndSpecies[]
956  */
957  vprob->addOnePhaseSpecies(VolPhase, k, kT);
958 
959  /*
960  * Get a pointer to the thermo object
961  */
962  ts_ptr = vprob->SpeciesThermo[kT];
963  /*
964  * Fill in the vcs_SpeciesProperty structure
965  */
966  vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
967  sProp->NumElements = vprob->ne;
968  sProp->SpName = vprob->SpName[kT];
969  sProp->SpeciesThermo = ts_ptr;
970  sProp->WtSpecies = tPhase->molecularWeight(k);
971  sProp->FormulaMatrixCol.resize(vprob->ne, 0.0);
972  for (size_t e = 0; e < vprob->ne; e++) {
973  sProp->FormulaMatrixCol[e] = vprob->FormulaMatrix(kT,e);
974  }
975  sProp->Charge = tPhase->charge(k);
976  sProp->SurfaceSpecies = false;
977  sProp->VolPM = 0.0;
978 
979  /*
980  * Transfer the thermo specification of the species
981  * vprob->SpeciesThermo[]
982  */
983  ts_ptr->UseCanteraCalls = VolPhase->usingCanteraCalls();
984  ts_ptr->m_VCS_UnitsFormat = VolPhase->p_VCS_UnitsFormat;
985  /*
986  * Add lookback connectivity into the thermo object first
987  */
988  ts_ptr->IndexPhase = iphase;
989  ts_ptr->IndexSpeciesPhase = k;
990  ts_ptr->OwningPhase = VolPhase;
991  /*
992  * get a reference to the Cantera species thermo.
993  */
994  SpeciesThermo& sp = tPhase->speciesThermo();
995 
996  int spType = sp.reportType(k);
997 
998  if (spType == SIMPLE) {
999  double c[4];
1000  double minTemp, maxTemp, refPressure;
1001  sp.reportParams(k, spType, c, minTemp, maxTemp, refPressure);
1002  ts_ptr->SS0_Model = VCS_SS0_CONSTANT;
1003  ts_ptr->SS0_T0 = c[0];
1004  ts_ptr->SS0_H0 = c[1];
1005  ts_ptr->SS0_S0 = c[2];
1006  ts_ptr->SS0_Cp0 = c[3];
1007  if (gasPhase) {
1008  ts_ptr->SSStar_Model = VCS_SSSTAR_IDEAL_GAS;
1010  } else {
1011  ts_ptr->SSStar_Model = VCS_SSSTAR_CONSTANT;
1012  ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
1013  }
1014  } else {
1015  if (vprob->m_printLvl > 2) {
1016  plogf("vcs_Cantera_convert: Species Type %d not known \n",
1017  spType);
1018  }
1019  ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED;
1020  ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED;
1021  if (!(ts_ptr->UseCanteraCalls)) {
1022  throw CanteraError("vcs_Cantera_to_vprob",
1023  "Cantera calls not being used -> aborting");
1024  }
1025  }
1026 
1027  /*
1028  * Transfer the Volume Information -> NEEDS WORK
1029  */
1030  if (gasPhase) {
1032  ts_ptr->SSStar_Vol0 = 82.05 * 273.15 / 1.0;
1033 
1034  } else {
1035  std::vector<double> phaseTermCoeff(nSpPhase, 0.0);
1036  int nCoeff;
1037  tPhase->getParameters(nCoeff, VCS_DATA_PTR(phaseTermCoeff));
1038  ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
1039  ts_ptr->SSStar_Vol0 = phaseTermCoeff[k];
1040  }
1041  kT++;
1042  }
1043 
1044  /*
1045  * Now go back through the species in the phase and assign
1046  * a valid mole fraction to all phases, even if the initial
1047  * estimate of the total number of moles is zero.
1048  */
1049  if (tMoles > 0.0) {
1050  for (size_t k = 0; k < nSpPhase; k++) {
1051  size_t kTa = VolPhase->spGlobalIndexVCS(k);
1052  vprob->mf[kTa] = vprob->w[kTa] / tMoles;
1053  }
1054  } else {
1055  /*
1056  * Perhaps, we could do a more sophisticated treatment below.
1057  * But, will start with this.
1058  */
1059  for (size_t k = 0; k < nSpPhase; k++) {
1060  size_t kTa = VolPhase->spGlobalIndexVCS(k);
1061  vprob->mf[kTa]= 1.0 / (double) nSpPhase;
1062  }
1063  }
1064 
1065  VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w));
1066  /*
1067  * Now, calculate a sample naught Gibbs free energy calculation
1068  * at the specified temperature.
1069  */
1070  double R = vcsUtil_gasConstant(vprob->m_VCS_UnitsFormat);
1071  for (size_t k = 0; k < nSpPhase; k++) {
1072  vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
1073  ts_ptr = sProp->SpeciesThermo;
1074  ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ R;
1075  ts_ptr->SS0_TSave = vprob->T;
1076  }
1077 
1078  }
1079 
1080  /*
1081  * Transfer initial element abundances to the vprob object.
1082  * We have to find the mapping index from one to the other
1083  *
1084  */
1085  vprob->gai.resize(vprob->ne, 0.0);
1086  vprob->set_gai();
1087 
1088  /*
1089  * Printout the species information: PhaseID's and mole nums
1090  */
1091  if (vprob->m_printLvl > 1) {
1092  writeline('=', 80, true, true);
1093  writeline('=', 16, false);
1094  plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
1095  writeline('=', 20);
1096  writeline('=', 80);
1097  plogf(" Phase IDs of species\n");
1098  plogf(" species phaseID phaseName ");
1099  plogf(" Initial_Estimated_kMols\n");
1100  for (size_t i = 0; i < vprob->nspecies; i++) {
1101  size_t iphase = vprob->PhaseID[i];
1102 
1103  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1104  plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
1105  VolPhase->PhaseName.c_str());
1107  plogf(" Volts = %-10.5g\n", vprob->w[i]);
1108  } else {
1109  plogf(" %-10.5g\n", vprob->w[i]);
1110  }
1111  }
1112 
1113  /*
1114  * Printout of the Phase structure information
1115  */
1116  writeline('-', 80, true, true);
1117  plogf(" Information about phases\n");
1118  plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
1119  plogf(" TMolesInert Tmoles(kmol)\n");
1120 
1121  for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) {
1122  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1123  std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
1124  plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
1125  VolPhase->VP_ID_, VolPhase->m_singleSpecies,
1126  VolPhase->m_gasPhase, sEOS.c_str(),
1127  VolPhase->nSpecies(), VolPhase->totalMolesInert());
1128  plogf("%16e\n", VolPhase->totalMoles());
1129  }
1130 
1131  writeline('=', 80, true, true);
1132  writeline('=', 16, false);
1133  plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
1134  writeline('=', 20);
1135  writeline('=', 80);
1136  plogf("\n");
1137  }
1138 
1139  return VCS_SUCCESS;
1140 }
1141 
1143 {
1144  size_t totNumPhases = mphase->nPhases();
1145  size_t kT = 0;
1146  std::vector<double> tmpMoles;
1147  // Problem type has yet to be worked out.
1148  vprob->prob_type = 0;
1149  // Whether we have an estimate or not gets overwritten on
1150  // the call to the equilibrium solver.
1151  vprob->iest = -1;
1152  vprob->T = mphase->temperature();
1153  vprob->PresPA = mphase->pressure();
1154  vprob->Vol = mphase->volume();
1155 
1156  for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
1157  ThermoPhase* tPhase = &(mphase->phase(iphase));
1158  vcs_VolPhase* volPhase = vprob->VPhaseList[iphase];
1159  /*
1160  * Set the electric potential of the volume phase from the
1161  * ThermoPhase object's value.
1162  */
1163  volPhase->setElectricPotential(tPhase->electricPotential());
1164 
1165  volPhase->setState_TP(vprob->T, vprob->PresPA);
1166  vector<double> muPhase(tPhase->nSpecies(),0.0);
1167  tPhase->getChemPotentials(&muPhase[0]);
1168  /*
1169  * Loop through each species in the current phase
1170  */
1171  size_t nSpPhase = tPhase->nSpecies();
1172  // volPhase->TMoles = 0.0;
1173  tmpMoles.resize(nSpPhase);
1174  for (size_t k = 0; k < nSpPhase; k++) {
1175  tmpMoles[k] = mphase->speciesMoles(kT);
1176  vprob->w[kT] = mphase->speciesMoles(kT);
1177  vprob->mf[kT] = mphase->moleFraction(kT);
1178 
1179  /*
1180  * transfer chemical potential vector
1181  */
1182  vprob->m_gibbsSpecies[kT] = muPhase[k];
1183 
1184  kT++;
1185  }
1186  if (volPhase->phiVarIndex() != npos) {
1187  size_t kphi = volPhase->phiVarIndex();
1188  size_t kglob = volPhase->spGlobalIndexVCS(kphi);
1189  vprob->w[kglob] = tPhase->electricPotential();
1190  }
1191  volPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w));
1192  if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) {
1194  } else if (volPhase->totalMoles() > 0.0) {
1195  volPhase->setExistence(VCS_PHASE_EXIST_YES);
1196  } else {
1197  volPhase->setExistence(VCS_PHASE_EXIST_NO);
1198  }
1199 
1200  }
1201  /*
1202  * Transfer initial element abundances to the vprob object.
1203  * Put them in the front of the object. There may be
1204  * more constraints than there are elements. But, we
1205  * know the element abundances are in the front of the
1206  * vector.
1207  */
1208  vprob->set_gai();
1209 
1210  /*
1211  * Printout the species information: PhaseID's and mole nums
1212  */
1213  if (vprob->m_printLvl > 1) {
1214  writeline('=', 80, true, true);
1215  writeline('=', 20, false);
1216  plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
1217  writeline('=', 20);
1218  writeline('=', 80);
1219  plogf("\n");
1220  plogf(" Phase IDs of species\n");
1221  plogf(" species phaseID phaseName ");
1222  plogf(" Initial_Estimated_kMols\n");
1223  for (size_t i = 0; i < vprob->nspecies; i++) {
1224  size_t iphase = vprob->PhaseID[i];
1225 
1226  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1227  plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
1228  VolPhase->PhaseName.c_str());
1230  plogf(" Volts = %-10.5g\n", vprob->w[i]);
1231  } else {
1232  plogf(" %-10.5g\n", vprob->w[i]);
1233  }
1234  }
1235 
1236  /*
1237  * Printout of the Phase structure information
1238  */
1239  writeline('-', 80, true, true);
1240  plogf(" Information about phases\n");
1241  plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
1242  plogf(" TMolesInert Tmoles(kmol)\n");
1243 
1244  for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) {
1245  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1246  std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
1247  plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
1248  VolPhase->VP_ID_, VolPhase->m_singleSpecies,
1249  VolPhase->m_gasPhase, sEOS.c_str(),
1250  VolPhase->nSpecies(), VolPhase->totalMolesInert());
1251  plogf("%16e\n", VolPhase->totalMoles());
1252  }
1253 
1254  writeline('=', 80, true, true);
1255  writeline('=', 20, false);
1256  plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
1257  writeline('=', 20);
1258  writeline('=', 80);
1259  plogf("\n");
1260  }
1261 
1262  return VCS_SUCCESS;
1263 }
1264 
1266 {
1267  size_t nsp = m_vsolve.m_numSpeciesTot;
1268  nu.resize(nsp, 0.0);
1269  for (size_t i = 0; i < nsp; i++) {
1270  nu[i] = 0.0;
1271  }
1272  size_t nc = numComponents();
1273  const std::vector<size_t>& indSpecies = m_vsolve.m_speciesMapIndex;
1274  if (rxn > nsp - nc) {
1275  return;
1276  }
1277  size_t j = indSpecies[rxn + nc];
1278  nu[j] = 1.0;
1279  for (size_t kc = 0; kc < nc; kc++) {
1280  j = indSpecies[kc];
1281  nu[j] = m_vsolve.m_stoichCoeffRxnMatrix(kc,rxn);
1282  }
1283 
1284 }
1285 
1287 {
1288  return m_vsolve.m_numComponents;
1289 }
1290 
1292 {
1294 }
1295 
1296 size_t vcs_MultiPhaseEquil::component(size_t m) const
1297 {
1298  size_t nc = numComponents();
1299  if (m < nc) {
1300  return m_vsolve.m_speciesMapIndex[m];
1301  } else {
1302  return npos;
1303  }
1304 }
1305 
1306 int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double& funcStab, int printLvl, int loglevel)
1307 {
1308  clockWC tickTock;
1309  m_printLvl = printLvl;
1310  m_vprob.m_printLvl = printLvl;
1311 
1312  /*
1313  * Extract the current state information
1314  * from the MultiPhase object and
1315  * Transfer it to VCS_PROB object.
1316  */
1317  int res = vcs_Cantera_update_vprob(m_mix, &m_vprob);
1318  if (res != 0) {
1319  plogf("problems\n");
1320  }
1321 
1322 
1323 
1324  // Check obvious bounds on the temperature and pressure
1325  // NOTE, we may want to do more here with the real bounds
1326  // given by the ThermoPhase objects.
1327  double T = m_mix->temperature();
1328  if (T <= 0.0) {
1329  throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability",
1330  "Temperature less than zero on input");
1331  }
1332  double pres = m_mix->pressure();
1333  if (pres <= 0.0) {
1334  throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability",
1335  "Pressure less than zero on input");
1336  }
1337 
1338  /*
1339  * Print out the problem specification from the point of
1340  * view of the vprob object.
1341  */
1343 
1344  /*
1345  * Call the thermo Program
1346  */
1347  int iStable = m_vsolve.vcs_PS(&m_vprob, iph, printLvl, funcStab);
1348 
1349  /*
1350  * Transfer the information back to the MultiPhase object.
1351  * Note we don't just call setMoles, because some multispecies
1352  * solution phases may be zeroed out, and that would cause a problem
1353  * for that routine. Also, the mole fractions of such zeroed out
1354  * phases actually contain information about likely reemergent
1355  * states.
1356  */
1359 
1360  double te = tickTock.secondsWC();
1361  if (printLvl > 0) {
1362  plogf("\n Results from vcs_PS:\n");
1363 
1364  plogf("\n");
1365  plogf("Temperature = %g Kelvin\n", m_vprob.T);
1366  plogf("Pressure = %g Pa\n", m_vprob.PresPA);
1367  std::string sss = m_mix->phaseName(iph);
1368  if (iStable) {
1369  plogf("Phase %d named %s is stable, function value = %g > 0\n", iph, sss.c_str(), funcStab);
1370  } else {
1371  plogf("Phase %d named %s is not stable + function value = %g < 0\n", iph, sss.c_str(), funcStab);
1372  }
1373  plogf("\n");
1374  plogf("----------------------------------------"
1375  "---------------------\n");
1376  plogf(" Name Mole_Number");
1377  if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
1378  plogf("(kmol)");
1379  } else {
1380  plogf("(gmol)");
1381  }
1382  plogf(" Mole_Fraction Chem_Potential");
1383  if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
1384  plogf(" (kcal/mol)\n");
1385  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
1386  plogf(" (Dimensionless)\n");
1387  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
1388  plogf(" (kJ/mol)\n");
1389  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
1390  plogf(" (Kelvin)\n");
1391  } else if (m_vprob.m_VCS_UnitsFormat == VCS_UNITS_MKS) {
1392  plogf(" (J/kmol)\n");
1393  }
1394  plogf("-------------------------------------------------------------\n");
1395  for (size_t i = 0; i < m_vprob.nspecies; i++) {
1396  plogf("%-12s", m_vprob.SpName[i].c_str());
1398  plogf(" %15.3e %15.3e ", 0.0, m_vprob.mf[i]);
1399  plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]);
1400  } else {
1401  plogf(" %15.3e %15.3e ", m_vprob.w[i], m_vprob.mf[i]);
1402  if (m_vprob.w[i] <= 0.0) {
1403  plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]);
1404  } else {
1405  plogf("%15.3e\n", m_vprob.m_gibbsSpecies[i]);
1406  }
1407  }
1408  }
1409  plogf("------------------------------------------"
1410  "-------------------\n");
1411  if (printLvl > 2) {
1412  if (m_vsolve.m_timing_print_lvl > 0) {
1413  plogf("Total time = %12.6e seconds\n", te);
1414  }
1415  }
1416  }
1417  return iStable;
1418 }
1419 
1420 }
double G0_calc_one(size_t kspec) const
Gibbs free energy calculation at a temperature for the reference state of a species, return a value for one species.
std::vector< double > WtSpecies
Molecular weight of species.
Definition: vcs_prob.h:187
#define VCS_PHASE_EXIST_NO
Phase doesn't currently exist in the mixture.
Definition: vcs_defs.h:225
doublereal molecularWeight(size_t k) const
Molecular weight of species k.
Definition: Phase.cpp:494
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
Definition: vcs_util.cpp:185
ThermoPhase object for the ideal molal equation of state (see Thermodynamic Properties and class Idea...
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
doublereal electricPotential() const
Returns the electric potential of this phase (V).
Definition: ThermoPhase.h:352
int m_printLvl
Print level for print routines.
Definition: vcs_prob.h:210
bool m_singleSpecies
If true, this phase consists of a single species.
Definition: vcs_VolPhase.h:584
int equilibrate_TP(int estimateEquil=0, int printLvl=0, doublereal err=1.0e-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using c...
vcs_VolPhase * OwningPhase
Pointer to the owning phase object.
size_t nElements() const
Number of elements.
Definition: Phase.cpp:167
MultiPhase * m_mix
Pointer to the MultiPhase mixture that will be equilibrated.
double totalMolesInert() const
returns the value of the total kmol of inert in the phase
virtual void getParameters(int &n, doublereal *const c) const
Get the equation of state parameters in a vector.
Definition: ThermoPhase.h:1456
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
int p_VCS_UnitsFormat
Units for the chemical potential data, pressure data, volume, and species amounts.
Definition: vcs_VolPhase.h:630
void setPhaseMoles(const size_t n, const doublereal moles)
Set the number of moles of phase with index n.
Definition: MultiPhase.cpp:943
Array2D m_stoichCoeffRxnMatrix
Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form...
Definition: vcs_solve.h:1552
doublereal speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
Definition: MultiPhase.cpp:247
double SS0_Cp0
Base heat capacity used in the VCS_SS0_CONSTANT_CP model.
#define VCS_PHASE_EXIST_ALWAYS
Always exists because it contains inerts which can't exist in any other phase.
Definition: vcs_defs.h:212
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
Definition: MultiPhase.cpp:904
size_t numComponents() const
reports the number of components in the equilibration problem
std::vector< int > SpeciesUnknownType
Specifies the species unknown type.
Definition: vcs_prob.h:103
#define VCS_DATA_PTR(vvv)
Points to the data in a std::vector<> object.
Definition: vcs_internal.h:18
double Vol
Volume of the entire system.
Definition: vcs_prob.h:123
int p_activityConvention
Convention for the activity formulation.
Definition: vcs_VolPhase.h:638
int iest
Specification of the initial estimate method.
Definition: vcs_prob.h:156
doublereal phaseMoles(const size_t n) const
Return the number of moles in phase n.
Definition: MultiPhase.cpp:938
std::vector< double > w
Total number of moles of the kth species.
Definition: vcs_prob.h:70
int equilibrate_HP(doublereal Htarget, int XY, double Tlow, double Thigh, int estimateEquil=0, int printLvl=0, doublereal err=1.0E-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using e...
int speciesUnknownType(const size_t k) const
Returns the type of the species unknown.
doublereal temperature() const
Temperature [K].
Definition: MultiPhase.h:334
size_t IndexSpeciesPhase
Index of this species in the current phase.
size_t nSpecies() const
Return the number of species in the phase.
friend int vcs_Cantera_update_vprob(MultiPhase *mphase, VCS_PROB *vprob)
Translate a MultiPhase information into a VCS_PROB problem definition object.
int m_printLvl
Print level from the VCSnonlinear package.
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition: global.h:270
#define VCS_SSVOL_IDEALGAS
Models for the standard state volume of each species.
Definition: vcs_VolPhase.h:23
double vcsUtil_gasConstant(int mu_units)
Returns the value of the gas constant in the units specified by parameter.
Definition: vcs_util.cpp:76
size_t addOnePhaseSpecies(vcs_VolPhase *volPhase, size_t k, size_t kT)
This routines adds entries for the formula matrix for one species.
Definition: vcs_prob.cpp:330
The class provides the wall clock timer in seconds.
Definition: clockWC.h:51
size_t transferElementsFM(const ThermoPhase *const tPhase)
Transfer all of the element information from the ThermoPhase object to the vcs_VolPhase object...
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Definition: ThermoPhase.h:670
std::string PhaseName
String name for the phase.
Definition: vcs_VolPhase.h:696
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
Definition: ThermoPhase.h:628
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
Definition: vcs_defs.h:215
double totalMoles() const
Return the total moles in the phase.
std::string name() const
Return the name of the phase.
Definition: Phase.cpp:157
Pure Virtual base class for the species thermo manager classes.
void getChemPotentials(doublereal *mu) const
Returns a vector of Chemical potentials.
Definition: MultiPhase.cpp:306
std::vector< vcs_VolPhase * > VPhaseList
Array of phase structures.
Definition: vcs_prob.h:193
virtual int reportType(size_t index=npos) const =0
This utility function reports the type of parameterization used for the species with index number ind...
int m_Iterations
Number of iterations. This is an output variable.
Definition: vcs_prob.h:204
void setExistence(const int existence)
Set the existence flag in the object.
void reportCSV(const std::string &reportFile)
Report the equilibrium answer in a comma separated table format.
int vcs_Cantera_update_vprob(MultiPhase *mphase, VCS_PROB *vprob)
Translate a MultiPhase information into a VCS_PROB problem definition object.
const doublereal Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition: ct_defs.h:134
void setPtrThermoPhase(ThermoPhase *tp_ptr)
Set the pointer for Cantera's ThermoPhase parameter.
std::vector< size_t > PhaseID
Mapping between the species and the phases.
Definition: vcs_prob.h:165
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:97
const int cSurf
A surface phase. Used by class SurfPhase.
Definition: mix_defs.h:40
size_t VP_ID_
Original ID of the phase in the problem.
Definition: vcs_VolPhase.h:581
int vcs_Cantera_to_vprob(MultiPhase *mphase, VCS_PROB *vprob)
Translate a MultiPhase object into a VCS_PROB problem definition object.
void setMolesFromVCS(const int stateCalc, const double *molesSpeciesVCS=0)
Set the moles within the phase.
Contains const definitions for types of species reference-state thermodynamics managers (see Species ...
Header for the object representing each phase within vcs.
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
Definition: vcs_defs.h:343
double SSStar_Vol0
parameter that is used in the VCS_SSVOL_CONSTANT model.
void uploadMoleFractionsFromPhases()
Update the locally-stored composition within this object to match the current compositions of the pha...
Definition: MultiPhase.cpp:963
size_t phiVarIndex() const
Return the index of the species that represents the the voltage of the phase.
A class for multiphase mixtures.
Definition: MultiPhase.h:53
virtual void getActivities(doublereal *a) const
Get the array of non-dimensional activities at the current solution temperature, pressure, and solution concentration.
int SSStar_Vol_Model
Models for the standard state volume of each species.
std::string SpName
Name of the species.
void addPhaseElements(vcs_VolPhase *volPhase)
Add elements to the local element list.
Definition: vcs_prob.cpp:285
Array2D FormulaMatrix
Formula Matrix for the problem.
Definition: vcs_prob.h:90
double SS0_S0
Base entropy used in the VCS_SS0_CONSTANT_CP model.
int equilibrate_SP(doublereal Starget, double Tlow, double Thigh, int estimateEquil=0, int printLvl=0, doublereal err=1.0E-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using c...
std::vector< size_t > m_speciesMapIndex
Index vector that keeps track of the species vector rearrangement.
Definition: vcs_solve.h:1781
doublereal entropy() const
The entropy of the mixture [J/K].
Definition: MultiPhase.cpp:385
int m_VCS_UnitsFormat
Units for the chemical potential data, pressure data, volume, and species amounts.
Definition: vcs_prob.h:148
#define VCS_SUCCESS
Definition: vcs_defs.h:21
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Definition: clockWC.cpp:35
size_t IndexPhase
Index of the phase that this species belongs to.
size_t numElemConstraints() const
Reports the number of element constraints in the equilibration problem.
size_t component(size_t m) const
Return the index of the ith component.
Properties of a single species.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: ThermoPhase.h:554
*std::vector< double > FormulaMatrixCol
Column of the formula matrix, comprising the element composition of the species.
size_t m_numElemConstraints
Number of element constraints in the problem.
Definition: vcs_solve.h:1499
int m_eqnState
Type of the equation of state.
Definition: vcs_VolPhase.h:597
Header file for an ideal solid solution model with incompressible thermodynamics (see Thermodynamic P...
size_t m_numSpeciesTot
Total number of species in the problems.
Definition: vcs_solve.h:1493
double SS0_H0
Base enthalpy used in the VCS_SS0_CONSTANT_CP model.
void setElectricPotential(const double phi)
set the electric potential of the phase
size_t nspecies
Total number of species in the problems.
Definition: vcs_prob.h:34
void setTemperature(const doublereal T)
Set the temperature [K].
Definition: MultiPhase.cpp:852
std::vector< double > mf
Mole fraction vector.
Definition: vcs_prob.h:77
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
int equilibrate(int XY, int estimateEquil=0, int printLvl=0, doublereal err=1.0e-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object...
int vcs(VCS_PROB *vprob, int ifunc, int ipr, int ip1, int maxit)
Solve an equilibrium problem.
Definition: vcs_solve.cpp:244
#define SIMPLE
Constant Cp thermo.
doublereal IntEnergy() const
The internal energy of the mixture [J].
Definition: MultiPhase.cpp:372
void resize(const size_t phaseNum, const size_t numSpecies, const size_t numElem, const char *const phaseName, const double molesInert=0.0)
The resize() function fills in all of the initial information if it is not given in the constructor...
int m_timing_print_lvl
printing level of timing information
Definition: vcs_solve.h:2002
const int cEdge
An edge between two 2D surfaces.
Definition: mix_defs.h:58
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1502
void prob_report(int print_lvl)
Print out the problem specification in all generality as it currently exists in the VCS_PROB object...
Definition: vcs_prob.cpp:164
bool m_gasPhase
If true, this phase is a gas-phase like phase.
Definition: vcs_VolPhase.h:591
void getStoichVector(size_t rxn, vector_fp &nu)
Get the stoichiometric reaction coefficients for a single reaction index.
void setPressure(doublereal P)
Set the pressure [Pa].
Definition: MultiPhase.h:424
std::string string16_EOSType(int EOSType)
Return a string representing the equation of state.
doublereal moleFraction(const size_t kGlob) const
Returns the mole fraction of global species k.
Definition: MultiPhase.cpp:953
int determine_PhaseStability(int iph, double &funcStab, int printLvl=0, int logLevel=-99)
Determine the phase stability of a phase at the current conditions.
std::vector< double > Charge
Charge of each species.
Definition: vcs_prob.h:190
int m_NumBasisOptimizations
Number of basis optimizations used. This is an output variable.
Definition: vcs_prob.h:207
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:265
double PresPA
Pressure.
Definition: vcs_prob.h:116
Interface class for the vcsnonlinear solver.
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:86
size_t ne
Number of element constraints in the equilibrium problem.
Definition: vcs_prob.h:40
int prob_type
Problem type. I.e., the identity of what is held constant.
Definition: vcs_prob.h:31
VCS_SOLVE m_vsolve
The object that does all of the equilibration work.
std::vector< VCS_SPECIES_THERMO * > SpeciesThermo
Vector of pointers to thermo structures which identify the model and parameters for evaluating the th...
Definition: vcs_prob.h:201
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:365
std::vector< double > gai
Element abundances for jth element.
Definition: vcs_prob.h:84
#define plogendl()
define this Cantera function to replace cout << endl;
Definition: vcs_internal.h:31
double SS0_T0
Base temperature used in the VCS_SS0_CONSTANT_CP model.
doublereal enthalpy() const
The enthalpy of the mixture [J].
Definition: MultiPhase.cpp:359
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:157
int SSStar_Model
Integer value representing the star state model.
std::vector< std::string > SpName
Vector of strings containing the species names.
Definition: vcs_prob.h:168
doublereal maxTemp() const
Maximum temperature for which all solution phases have valid thermo data.
Definition: MultiPhase.h:257
virtual int eosType() const
Equation of state type flag.
Definition: ThermoPhase.h:142
int SS0_Model
Integer representing the models for the species standard state Naught temperature dependence...
VCS_PROB m_vprob
Object which contains the problem statement.
std::string phaseName(const size_t iph) const
Returns the name of the n'th phase.
Definition: MultiPhase.cpp:919
vcs_SpeciesProperties * speciesProperty(const size_t kindex)
Retrieve the kth Species structure for the species belonging to this phase.
virtual void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
Definition: ThermoPhase.h:512
friend int vcs_Cantera_to_vprob(MultiPhase *mphase, VCS_PROB *vprob)
Translate a MultiPhase object into a VCS_PROB problem definition object.
void getMolalities(doublereal *const molal) const
This function will return the molalities of the species.
virtual int activityConvention() const
This method returns the convention used in specification of the activities, of which there are curren...
std::vector< double > m_gibbsSpecies
Vector of chemical potentials of the species.
Definition: vcs_prob.h:57
Contains declarations for string manipulation functions within Cantera.
double T
Temperature (Kelvin)
Definition: vcs_prob.h:109
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:24
bool UseCanteraCalls
If true, this object will call Cantera to do its member calculations.
virtual SpeciesThermo & speciesThermo(int k=-1)
Return a changeable reference to the calculation manager for species reference-state thermodynamic pr...
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:352
void set_gai()
Calculate the element abundance vector from the mole numbers.
Definition: vcs_prob.cpp:150
virtual void reportParams(size_t index, int &type, doublereal *const c, doublereal &minTemp, doublereal &maxTemp, doublereal &refPressure) const =0
This utility function reports back the type of parameterization and all of the parameters for the spe...
thermo_t & phase(size_t n)
Return a reference to phase n.
Definition: MultiPhase.cpp:222
doublereal minTemp() const
Minimum temperature for which all solution phases have valid thermo data.
Definition: MultiPhase.h:250
doublereal pressure() const
Pressure [Pa].
Definition: MultiPhase.h:409
size_t NPhase
Number of phases in the problem.
Definition: vcs_prob.h:47
size_t nSpecies() const
Number of species, summed over all phases.
Definition: MultiPhase.h:124
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:230
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
bool usingCanteraCalls() const
Returns whether the object is using cantera calls.
Interface class for the vcs thermo equilibrium solver package, which generally describes the problem ...
Definition: vcs_prob.h:24
void setState_TP(const double temperature_Kelvin, const double pressure_PA)
Sets the temperature and pressure in this object and underlying ThermoPhase objects.
const int cIdealGas
Equation of state types:
Definition: mix_defs.h:37
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:272
int equilibrate_TV(int XY, doublereal xtarget, int estimateEquil=0, int printLvl=0, doublereal err=1.0E-6, int maxsteps=VCS_MAXSTEPS, int logLevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object using c...
doublereal volume() const
The total mixture volume [m^3].
Definition: MultiPhase.cpp:550
double SS0_TSave
Internal storage of the last temperature used in the calculation of the reference naught Gibbs free e...
double SS0_feSave
Internal storage of the last calculation of the reference naught Gibbs free energy at SS0_TSave...
void setTotalMoles(const double totalMols)
Sets the total moles in the phase.
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:578
size_t nPhases() const
Number of phases.
Definition: MultiPhase.h:446