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