Cantera  2.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  * 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 
11 #include "cantera/equil/vcs_MultiPhaseEquil.h"
12 #include "cantera/equil/vcs_prob.h"
13 #include "cantera/equil/vcs_internal.h"
14 #include "cantera/equil/vcs_VolPhase.h"
15 #include "vcs_species_thermo.h"
16 #include "vcs_SpeciesProperties.h"
17 
18 #include "cantera/equil/vcs_solve.h"
19 
20 #include "cantera/base/ct_defs.h"
21 #include "cantera/thermo/mix_defs.h"
22 #include "cantera/base/clockWC.h"
23 #include "cantera/thermo/ThermoPhase.h"
24 #include "cantera/thermo/speciesThermoTypes.h"
25 #include "cantera/thermo/IdealSolidSolnPhase.h"
26 #include "cantera/thermo/IdealMolalSoln.h"
27 #include "cantera/equil/ChemEquil.h"
28 
29 #include <string>
30 #include <vector>
31 #include <cstdio>
32 
33 using namespace Cantera;
34 using namespace std;
35 //using namespace VCSnonideal;
36 
37 namespace VCSnonideal
38 {
39 //====================================================================================================================
40 
41 vcs_MultiPhaseEquil::vcs_MultiPhaseEquil() :
42  m_vprob(0),
43  m_mix(0),
44  m_printLvl(0),
45  m_vsolvePtr(0)
46 {
47 }
48 //====================================================================================================================
49 vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(Cantera::MultiPhase* mix, int printLvl) :
50  m_vprob(0),
51  m_mix(0),
52  m_printLvl(printLvl),
53  m_vsolvePtr(0)
54 {
55  /*
56  * Create a VCS_PROB object that describes the equilibrium problem.
57  * The constructor just mallocs the necessary objects and sizes them.
58  */
59  m_vprob = new VCS_PROB(mix->nSpecies(),
60  mix->nElements(),
61  mix->nPhases());
62  m_mix = mix;
63  m_vprob->m_printLvl = m_printLvl;
64  /*
65  * Work out the details of the VCS_VPROB construction and
66  * Transfer the current problem to VCS_PROB object
67  */
68  int res = vcs_Cantera_to_vprob(mix, m_vprob);
69  if (res != 0) {
70  plogf("problems\n");
71  }
72 }
73 
74 vcs_MultiPhaseEquil::~vcs_MultiPhaseEquil()
75 {
76  delete m_vprob;
77  m_vprob = 0;
78  if (m_vsolvePtr) {
79  delete m_vsolvePtr;
80  m_vsolvePtr = 0;
81  }
82 }
83 //====================================================================================================================
84 int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
85  int estimateEquil,
86  int printLvl, doublereal err,
87  int maxsteps, int loglevel)
88 {
89 
90  addLogEntry("problem type","fixed T, V");
91  // doublereal dt = 1.0e3;
92  doublereal Vtarget = m_mix->volume();
93  doublereal dVdP;
94  if ((XY != TV) && (XY != HV) && (XY != UV) && (XY != SV)) {
95  throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TV",
96  "Wrong XY flag:" + int2str(XY));
97  }
98  int maxiter = 100;
99  int iSuccess = 0;
100  int innerXY;
101  double Pnow;
102  if (XY == TV) {
103  m_mix->setTemperature(xtarget);
104  }
105  double Pnew;
106  int strt = estimateEquil;
107  double P1 = 0.0;
108  double V1 = 0.0;
109  double V2 = 0.0;
110  double P2 = 0.0;
111  doublereal Tlow = 0.5 * m_mix->minTemp();
112  doublereal Thigh = 2.0 * m_mix->maxTemp();
113  doublereal Vnow, Verr;
114  int printLvlSub = std::max(0, printLvl - 1);
115  for (int n = 0; n < maxiter; n++) {
116  Pnow = m_mix->pressure();
117 
118  beginLogGroup("iteration "+int2str(n));
119  switch (XY) {
120  case TV:
121  iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
122  break;
123  case HV:
124  innerXY = HP;
125  iSuccess = equilibrate_HP(xtarget, innerXY, Tlow, Thigh, strt,
126  printLvlSub, err, maxsteps, loglevel);
127  break;
128  case UV:
129  innerXY = UP;
130  iSuccess = equilibrate_HP(xtarget, innerXY, Tlow, Thigh, strt,
131  printLvlSub, err, maxsteps, loglevel);
132  break;
133  case SV:
134  innerXY = SP;
135  iSuccess = equilibrate_SP(xtarget, Tlow, Thigh, strt,
136  printLvlSub, err, maxsteps, loglevel);
137  break;
138  default:
139  break;
140  }
141  strt = false;
142  Vnow = m_mix->volume();
143  if (n == 0) {
144  V2 = Vnow;
145  P2 = Pnow;
146  } else if (n == 1) {
147  V1 = Vnow;
148  P1 = Pnow;
149  } else {
150  P2 = P1;
151  V2 = V1;
152  P1 = Pnow;
153  V1 = Vnow;
154  }
155 
156  Verr = fabs((Vtarget - Vnow)/Vtarget);
157  addLogEntry("P",fp2str(Pnow));
158  addLogEntry("V rel error",fp2str(Verr));
159  endLogGroup();
160 
161  if (Verr < err) {
162  addLogEntry("P iterations",int2str(n));
163  addLogEntry("Final P",fp2str(Pnow));
164  addLogEntry("V rel error",fp2str(Verr));
165  goto done;
166  }
167  // find dV/dP
168  if (n > 1) {
169  dVdP = (V2 - V1) / (P2 - P1);
170  if (dVdP == 0.0) {
171  throw CanteraError("vcs_MultiPhase::equilibrate_TV",
172  "dVdP == 0.0");
173  } else {
174  Pnew = Pnow + (Vtarget - Vnow) / dVdP;
175  if (Pnew < 0.2 * Pnow) {
176  Pnew = 0.2 * Pnow;
177  }
178  if (Pnew > 3.0 * Pnow) {
179  Pnew = 3.0 * Pnow;
180  }
181  }
182 
183  } else {
184  m_mix->setPressure(Pnow*1.01);
185  dVdP = (m_mix->volume() - Vnow)/(0.01*Pnow);
186  Pnew = Pnow + 0.5*(Vtarget - Vnow)/dVdP;
187  if (Pnew < 0.5* Pnow) {
188  Pnew = 0.5 * Pnow;
189  }
190  if (Pnew > 1.7 * Pnow) {
191  Pnew = 1.7 * Pnow;
192  }
193 
194  }
195  m_mix->setPressure(Pnew);
196  }
197  throw CanteraError("vcs_MultiPhase::equilibrate_TV",
198  "No convergence for V");
199 
200 done:
201  ;
202  return iSuccess;
203 }
204 
205 //====================================================================================================================
206 int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
207  int XY, double Tlow, double Thigh,
208  int estimateEquil,
209  int printLvl, doublereal err,
210  int maxsteps, int loglevel)
211 {
212  int maxiter = 100;
213  int iSuccess;
214  if (XY != HP && XY != UP) {
215  throw CanteraError("vcs_MultiPhaseEquil::equilibrate_HP",
216  "Wrong XP" + XY);
217  }
218  int strt = estimateEquil;
219 
220  // Lower bound on T. This will change as we progress in the calculation
221  if (Tlow <= 0.0) {
222  Tlow = 0.5 * m_mix->minTemp();
223  }
224  // Upper bound on T. This will change as we progress in the calculation
225  if (Thigh <= 0.0 || Thigh > 1.0E6) {
226  Thigh = 2.0 * m_mix->maxTemp();
227  }
228  addLogEntry("problem type","fixed H,P");
229  addLogEntry("H target",fp2str(Htarget));
230 
231  doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
232  doublereal Hnow;
233  doublereal Hlow = Undef;
234  doublereal Hhigh = Undef;
235  doublereal Herr, HConvErr;
236  doublereal Tnow = m_mix->temperature();
237  int printLvlSub = std::max(printLvl - 1, 0);
238 
239  for (int n = 0; n < maxiter; n++) {
240 
241 
242  // start with a loose error tolerance, but tighten it as we get
243  // close to the final temperature
244  beginLogGroup("iteration "+int2str(n));
245 
246  try {
247  Tnow = m_mix->temperature();
248  iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
249  strt = 0;
250  if (XY == UP) {
251  Hnow = m_mix->IntEnergy();
252  } else {
253  Hnow = m_mix->enthalpy();
254  }
255  double pmoles[10];
256  pmoles[0] = m_mix->phaseMoles(0);
257  double Tmoles = pmoles[0];
258  double HperMole = Hnow/Tmoles;
259  if (printLvl > 0) {
260  plogf("T = %g, Hnow = %g ,Tmoles = %g, HperMole = %g",
261  Tnow, Hnow, Tmoles, HperMole);
262  plogendl();
263  }
264 
265  // the equilibrium enthalpy monotonically increases with T;
266  // if the current value is below the target, then we know the
267  // current temperature is too low. Set the lower bounds.
268 
269 
270  if (Hnow < Htarget) {
271  if (Tnow > Tlow) {
272  Tlow = Tnow;
273  Hlow = Hnow;
274  }
275  }
276  // the current enthalpy is greater than the target; therefore the
277  // current temperature is too high. Set the high bounds.
278  else {
279  if (Tnow < Thigh) {
280  Thigh = Tnow;
281  Hhigh = Hnow;
282  }
283  }
284  if (Hlow != Undef && Hhigh != Undef) {
285  cpb = (Hhigh - Hlow)/(Thigh - Tlow);
286  dT = (Htarget - Hnow)/cpb;
287  dTa = fabs(dT);
288  dTmax = 0.5*fabs(Thigh - Tlow);
289  if (dTa > dTmax) {
290  dT *= dTmax/dTa;
291  }
292  } else {
293  Tnew = sqrt(Tlow*Thigh);
294  dT = Tnew - Tnow;
295  if (dT < -200.) {
296  dT = -200;
297  }
298  if (dT > 200.) {
299  dT = 200.;
300  }
301  }
302  double acpb = std::max(fabs(cpb), 1.0E-6);
303  double denom = std::max(fabs(Htarget), acpb);
304  Herr = Htarget - Hnow;
305  HConvErr = fabs((Herr)/denom);
306  addLogEntry("T",fp2str(m_mix->temperature()));
307  addLogEntry("H",fp2str(Hnow));
308  addLogEntry("Herr",fp2str(Herr));
309  addLogEntry("H rel error",fp2str(HConvErr));
310  addLogEntry("lower T bound",fp2str(Tlow));
311  addLogEntry("upper T bound",fp2str(Thigh));
312  endLogGroup(); // iteration
313  if (printLvl > 0) {
314  plogf(" equilibrate_HP: It = %d, Tcurr = %g Hcurr = %g, Htarget = %g\n",
315  n, Tnow, Hnow, Htarget);
316  plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
317  Herr, cpb, HConvErr);
318  }
319 
320  if (HConvErr < err) { // || dTa < 1.0e-4) {
321  addLogEntry("T iterations",int2str(n));
322  addLogEntry("Final T",fp2str(m_mix->temperature()));
323  addLogEntry("H rel error",fp2str(Herr));
324  if (printLvl > 0) {
325  plogf(" equilibrate_HP: CONVERGENCE: Hfinal = %g Tfinal = %g, Its = %d \n",
326  Hnow, Tnow, n);
327  plogf(" H rel error = %g, cp = %g, HConvErr = %g\n",
328  Herr, cpb, HConvErr);
329  }
330  goto done;
331  }
332  Tnew = Tnow + dT;
333  if (Tnew < 0.0) {
334  Tnew = 0.5*Tnow;
335  }
336  m_mix->setTemperature(Tnew);
337 
338  } catch (CanteraError err) {
339  if (!estimateEquil) {
340  addLogEntry("no convergence",
341  "try estimating composition at the start");
342  strt = -1;
343  } else {
344  Tnew = 0.5*(Tnow + Thigh);
345  if (fabs(Tnew - Tnow) < 1.0) {
346  Tnew = Tnow + 1.0;
347  }
348  m_mix->setTemperature(Tnew);
349  addLogEntry("no convergence",
350  "trying T = "+fp2str(Tnow));
351  }
352  endLogGroup();
353  }
354 
355  }
356  addLogEntry("reached max number of T iterations",int2str(maxiter));
357  endLogGroup();
358  throw CanteraError("MultiPhase::equilibrate_HP",
359  "No convergence for T");
360 done:
361  ;
362  return iSuccess;
363 }
364 //====================================================================================================================
365 int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
366  double Tlow, double Thigh,
367  int estimateEquil,
368  int printLvl, doublereal err,
369  int maxsteps, int loglevel)
370 {
371  int maxiter = 100;
372  int strt = estimateEquil;
373 
374  // Lower bound on T. This will change as we progress in the calculation
375  if (Tlow <= 0.0) {
376  Tlow = 0.5 * m_mix->minTemp();
377  }
378  // Upper bound on T. This will change as we progress in the calculation
379  if (Thigh <= 0.0 || Thigh > 1.0E6) {
380  Thigh = 2.0 * m_mix->maxTemp();
381  }
382  addLogEntry("problem type","fixed S,P");
383  addLogEntry("S target",fp2str(Starget));
384 
385  doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
386  doublereal Snow;
387  doublereal Slow = Undef;
388  doublereal Shigh = Undef;
389  doublereal Serr, SConvErr;
390  doublereal Tnow = m_mix->temperature();
391  if (Tnow < Tlow) {
392  Tlow = Tnow;
393  }
394  if (Tnow > Thigh) {
395  Thigh = Tnow;
396  }
397  int printLvlSub = std::max(printLvl - 1, 0);
398 
399  for (int n = 0; n < maxiter; n++) {
400 
401  // start with a loose error tolerance, but tighten it as we get
402  // close to the final temperature
403  beginLogGroup("iteration "+int2str(n));
404 
405  try {
406  Tnow = m_mix->temperature();
407  int iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
408  strt = 0;
409  Snow = m_mix->entropy();
410  double pmoles[10];
411  pmoles[0] = m_mix->phaseMoles(0);
412  double Tmoles = pmoles[0];
413  double SperMole = Snow/Tmoles;
414  plogf("T = %g, Snow = %g ,Tmoles = %g, SperMole = %g\n",
415  Tnow, Snow, Tmoles, SperMole);
416 
417  // the equilibrium entropy monotonically increases with T;
418  // if the current value is below the target, then we know the
419  // current temperature is too low. Set the lower bounds to the
420  // current condition.
421  if (Snow < Starget) {
422  if (Tnow > Tlow) {
423  Tlow = Tnow;
424  Slow = Snow;
425  } else {
426  if (Slow > Starget) {
427  if (Snow < Slow) {
428  Thigh = Tlow;
429  Shigh = Slow;
430  Tlow = Tnow;
431  Slow = Snow;
432  }
433  }
434  }
435  }
436  // the current enthalpy is greater than the target; therefore the
437  // current temperature is too high. Set the high bounds.
438  else {
439  if (Tnow < Thigh) {
440  Thigh = Tnow;
441  Shigh = Snow;
442  }
443  }
444  if (Slow != Undef && Shigh != Undef) {
445  cpb = (Shigh - Slow)/(Thigh - Tlow);
446  dT = (Starget - Snow)/cpb;
447  Tnew = Tnow + dT;
448  dTa = fabs(dT);
449  dTmax = 0.5*fabs(Thigh - Tlow);
450  if (Tnew > Thigh || Tnew < Tlow) {
451  dTmax = 1.5*fabs(Thigh - Tlow);
452  }
453  dTmax = std::min(dTmax, 300.);
454  if (dTa > dTmax) {
455  dT *= dTmax/dTa;
456  }
457  } else {
458  Tnew = sqrt(Tlow*Thigh);
459  dT = Tnew - Tnow;
460  }
461 
462  double acpb = std::max(fabs(cpb), 1.0E-6);
463  double denom = std::max(fabs(Starget), acpb);
464  Serr = Starget - Snow;
465  SConvErr = fabs((Serr)/denom);
466  addLogEntry("T",fp2str(m_mix->temperature()));
467  addLogEntry("S",fp2str(Snow));
468  addLogEntry("Serr",fp2str(Serr));
469  addLogEntry("S rel error",fp2str(SConvErr));
470  addLogEntry("lower T bound",fp2str(Tlow));
471  addLogEntry("upper T bound",fp2str(Thigh));
472  endLogGroup(); // iteration
473  if (printLvl > 0) {
474  plogf(" equilibrate_SP: It = %d, Tcurr = %g Scurr = %g, Starget = %g\n",
475  n, Tnow, Snow, Starget);
476  plogf(" S rel error = %g, cp = %g, SConvErr = %g\n",
477  Serr, cpb, SConvErr);
478  }
479 
480  if (SConvErr < err) { // || dTa < 1.0e-4) {
481  addLogEntry("T iterations",int2str(n));
482  addLogEntry("Final T",fp2str(m_mix->temperature()));
483  addLogEntry("S rel error",fp2str(Serr));
484  if (printLvl > 0) {
485  plogf(" equilibrate_SP: CONVERGENCE: Sfinal = %g Tfinal = %g, Its = %d \n",
486  Snow, Tnow, n);
487  plogf(" S rel error = %g, cp = %g, HConvErr = %g\n",
488  Serr, cpb, SConvErr);
489  }
490  return iSuccess;
491  }
492  Tnew = Tnow + dT;
493  if (Tnew < 0.0) {
494  Tnew = 0.5*Tnow;
495  }
496  m_mix->setTemperature(Tnew);
497 
498  } catch (CanteraError err) {
499  if (!estimateEquil) {
500  addLogEntry("no convergence",
501  "try estimating composition at the start");
502  strt = -1;
503  } else {
504  Tnew = 0.5*(Tnow + Thigh);
505  if (fabs(Tnew - Tnow) < 1.0) {
506  Tnew = Tnow + 1.0;
507  }
508  m_mix->setTemperature(Tnew);
509  addLogEntry("no convergence",
510  "trying T = "+fp2str(Tnow));
511  }
512  endLogGroup();
513  }
514 
515  }
516  addLogEntry("reached max number of T iterations",int2str(maxiter));
517  endLogGroup();
518  throw CanteraError("MultiPhase::equilibrate_SP",
519  "No convergence for T");
520 }
521 //====================================================================================================================
522 
523 /*
524  * Equilibrate the solution using the current element abundances
525  */
526 int vcs_MultiPhaseEquil::equilibrate(int XY, int estimateEquil,
527  int printLvl, doublereal err,
528  int maxsteps, int loglevel)
529 {
530  int iSuccess;
531  doublereal xtarget;
532  if (XY == TP) {
533  iSuccess = equilibrate_TP(estimateEquil, printLvl, err, maxsteps, loglevel);
534  } else if (XY == HP || XY == UP) {
535  if (XY == HP) {
536  xtarget = m_mix->enthalpy();
537  } else {
538  xtarget = m_mix->IntEnergy();
539  }
540  double Tlow = 0.5 * m_mix->minTemp();
541  double Thigh = 2.0 * m_mix->maxTemp();
542  iSuccess = equilibrate_HP(xtarget, XY, Tlow, Thigh,
543  estimateEquil, printLvl, err, maxsteps, loglevel);
544  } else if (XY == SP) {
545  xtarget = m_mix->entropy();
546  double Tlow = 0.5 * m_mix->minTemp();
547  double Thigh = 2.0 * m_mix->maxTemp();
548  iSuccess = equilibrate_SP(xtarget, Tlow, Thigh,
549  estimateEquil, printLvl, err, maxsteps, loglevel);
550 
551  } else if (XY == TV) {
552  xtarget = m_mix->temperature();
553  iSuccess = equilibrate_TV(XY, xtarget,
554  estimateEquil, printLvl, err, maxsteps, loglevel);
555  } else if (XY == HV) {
556  xtarget = m_mix->enthalpy();
557  iSuccess = equilibrate_TV(XY, xtarget,
558  estimateEquil, printLvl, err, maxsteps, loglevel);
559  } else if (XY == UV) {
560  xtarget = m_mix->IntEnergy();
561  iSuccess = equilibrate_TV(XY, xtarget,
562  estimateEquil, printLvl, err, maxsteps, loglevel);
563  } else if (XY == SV) {
564  xtarget = m_mix->entropy();
565  iSuccess = equilibrate_TV(XY, xtarget, estimateEquil,
566  printLvl, err, maxsteps, loglevel);
567  } else {
568  throw CanteraError(" vcs_MultiPhaseEquil::equilibrate",
569  "Unsupported Option");
570  }
571  return iSuccess;
572 }
573 //====================================================================================================================
574 /*
575  * Equilibrate the solution using the current element abundances
576  */
577 int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
578  int printLvl, doublereal err,
579  int maxsteps, int loglevel)
580 {
581  // Debugging level
582 
583  int maxit = maxsteps;
584  clockWC tickTock;
585 
586  if (m_vprob == 0) {
587  m_vprob = new VCS_PROB(m_mix->nSpecies(),
588  m_mix->nElements(),
589  m_mix->nPhases());
590  }
591  m_printLvl = printLvl;
592  m_vprob->m_printLvl = printLvl;
593 
594 
595  /*
596  * Extract the current state information
597  * from the MultiPhase object and
598  * Transfer it to VCS_PROB object.
599  */
600  int res = vcs_Cantera_update_vprob(m_mix, m_vprob);
601  if (res != 0) {
602  plogf("problems\n");
603  }
604 
605 
606  // Set the estimation technique
607  if (estimateEquil) {
608  m_vprob->iest = estimateEquil;
609  } else {
610  m_vprob->iest = 0;
611  }
612 
613  // Check obvious bounds on the temperature and pressure
614  // NOTE, we may want to do more here with the real bounds
615  // given by the ThermoPhase objects.
616  double T = m_mix->temperature();
617  if (T <= 0.0) {
618  throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
619  "Temperature less than zero on input");
620  }
621  double pres = m_mix->pressure();
622  if (pres <= 0.0) {
623  throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
624  "Pressure less than zero on input");
625  }
626 
627  beginLogGroup("vcs_MultiPhaseEquil::equilibrate_TP", loglevel);
628  addLogEntry("problem type","fixed T,P");
629  addLogEntry("Temperature", T);
630  addLogEntry("Pressure", pres);
631 
632 
633  /*
634  * Print out the problem specification from the point of
635  * view of the vprob object.
636  */
637  m_vprob->prob_report(m_printLvl);
638 
639  /*
640  * Call the thermo Program
641  */
642  int ip1 = m_printLvl;
643  int ipr = std::max(0, m_printLvl-1);
644  if (m_printLvl >= 3) {
645  ip1 = m_printLvl - 2;
646  } else {
647  ip1 = 0;
648  }
649  if (!m_vsolvePtr) {
650  m_vsolvePtr = new VCS_SOLVE();
651  }
652  int iSuccess = m_vsolvePtr->vcs(m_vprob, 0, ipr, ip1, maxit);
653 
654  /*
655  * Transfer the information back to the MultiPhase object.
656  * Note we don't just call setMoles, because some multispecies
657  * solution phases may be zeroed out, and that would cause a problem
658  * for that routine. Also, the mole fractions of such zereod out
659  * phases actually contain information about likely reemergent
660  * states.
661  */
662  m_mix->uploadMoleFractionsFromPhases();
663  size_t kGlob = 0;
664  for (size_t ip = 0; ip < m_vprob->NPhase; ip++) {
665  double phaseMole = 0.0;
666  Cantera::ThermoPhase& tref = m_mix->phase(ip);
667  for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) {
668  phaseMole += m_vprob->w[kGlob];
669  }
670  //phaseMole *= 1.0E-3;
671  m_mix->setPhaseMoles(ip, phaseMole);
672  }
673 
674  double te = tickTock.secondsWC();
675  if (printLvl > 0) {
676  plogf("\n Results from vcs:\n");
677  if (iSuccess != 0) {
678  plogf("\nVCS FAILED TO CONVERGE!\n");
679  }
680  plogf("\n");
681  plogf("Temperature = %g Kelvin\n", m_vprob->T);
682  plogf("Pressure = %g Pa\n", m_vprob->PresPA);
683  plogf("\n");
684  plogf("----------------------------------------"
685  "---------------------\n");
686  plogf(" Name Mole_Number");
687  if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
688  plogf("(kmol)");
689  } else {
690  plogf("(gmol)");
691  }
692  plogf(" Mole_Fraction Chem_Potential");
693  if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
694  plogf(" (kcal/mol)\n");
695  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
696  plogf(" (Dimensionless)\n");
697  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
698  plogf(" (kJ/mol)\n");
699  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
700  plogf(" (Kelvin)\n");
701  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
702  plogf(" (J/kmol)\n");
703  }
704  plogf("--------------------------------------------------"
705  "-----------\n");
706  for (size_t i = 0; i < m_vprob->nspecies; i++) {
707  plogf("%-12s", m_vprob->SpName[i].c_str());
708  if (m_vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
709  plogf(" %15.3e %15.3e ", 0.0, m_vprob->mf[i]);
710  plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
711  } else {
712  plogf(" %15.3e %15.3e ", m_vprob->w[i], m_vprob->mf[i]);
713  if (m_vprob->w[i] <= 0.0) {
714  size_t iph = m_vprob->PhaseID[i];
715  vcs_VolPhase* VPhase = m_vprob->VPhaseList[iph];
716  if (VPhase->nSpecies() > 1) {
717  plogf(" -1.000e+300\n");
718  } else {
719  plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
720  }
721  } else {
722  plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
723  }
724  }
725  }
726  plogf("------------------------------------------"
727  "-------------------\n");
728  if (printLvl > 2) {
729  if (m_vsolvePtr->m_timing_print_lvl > 0) {
730  plogf("Total time = %12.6e seconds\n", te);
731  }
732  }
733  }
734  if (loglevel > 0) {
735  endLogGroup();
736  }
737  return iSuccess;
738 }
739 
740 
741 //====================================================================================================================
742 /**************************************************************************
743  *
744  *
745  */
746 void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
747 {
748  size_t k;
749  size_t istart;
750  size_t nSpecies;
751 
752  double vol = 0.0;
753  string sName;
754  size_t nphase = m_vprob->NPhase;
755 
756  FILE* FP = fopen(reportFile.c_str(), "w");
757  if (!FP) {
758  plogf("Failure to open file\n");
759  exit(EXIT_FAILURE);
760  }
761  double Temp = m_mix->temperature();
762  double pres = m_mix->pressure();
763  double* mf = VCS_DATA_PTR(m_vprob->mf);
764 #ifdef DEBUG_MODE
765  double* fe = VCS_DATA_PTR(m_vprob->m_gibbsSpecies);
766 #endif
767  std::vector<double> VolPM;
768  std::vector<double> activity;
769  std::vector<double> ac;
770  std::vector<double> mu;
771  std::vector<double> mu0;
772  std::vector<double> molalities;
773 
774 
775  vol = 0.0;
776  for (size_t iphase = 0; iphase < nphase; iphase++) {
777  istart = m_mix->speciesIndex(0, iphase);
778  Cantera::ThermoPhase& tref = m_mix->phase(iphase);
779  nSpecies = tref.nSpecies();
780  VolPM.resize(nSpecies, 0.0);
781  tref.getPartialMolarVolumes(VCS_DATA_PTR(VolPM));
782  vcs_VolPhase* volP = m_vprob->VPhaseList[iphase];
783 
784  double TMolesPhase = volP->totalMoles();
785  double VolPhaseVolumes = 0.0;
786  for (k = 0; k < nSpecies; k++) {
787  VolPhaseVolumes += VolPM[k] * mf[istart + k];
788  }
789  VolPhaseVolumes *= TMolesPhase;
790  vol += VolPhaseVolumes;
791  }
792 
793  fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
794  " -----------------------------\n");
795  fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
796  fprintf(FP,"Pressure = %11.5g Pascal\n", pres);
797  fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
798  fprintf(FP,"Number Basis optimizations = %d\n", m_vprob->m_NumBasisOptimizations);
799  fprintf(FP,"Number VCS iterations = %d\n", m_vprob->m_Iterations);
800 
801  for (size_t iphase = 0; iphase < nphase; iphase++) {
802  istart = m_mix->speciesIndex(0, iphase);
803  Cantera::ThermoPhase& tref = m_mix->phase(iphase);
804  Cantera::ThermoPhase* tp = &tref;
805  string phaseName = tref.name();
806  vcs_VolPhase* volP = m_vprob->VPhaseList[iphase];
807  double TMolesPhase = volP->totalMoles();
808  //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase));
809  nSpecies = tref.nSpecies();
810  activity.resize(nSpecies, 0.0);
811  ac.resize(nSpecies, 0.0);
812 
813  mu0.resize(nSpecies, 0.0);
814  mu.resize(nSpecies, 0.0);
815  VolPM.resize(nSpecies, 0.0);
816  molalities.resize(nSpecies, 0.0);
817 
818  int actConvention = tp->activityConvention();
819  tp->getActivities(VCS_DATA_PTR(activity));
820  tp->getActivityCoefficients(VCS_DATA_PTR(ac));
821  tp->getStandardChemPotentials(VCS_DATA_PTR(mu0));
822 
823  tp->getPartialMolarVolumes(VCS_DATA_PTR(VolPM));
824  tp->getChemPotentials(VCS_DATA_PTR(mu));
825  double VolPhaseVolumes = 0.0;
826  for (k = 0; k < nSpecies; k++) {
827  VolPhaseVolumes += VolPM[k] * mf[istart + k];
828  }
829  VolPhaseVolumes *= TMolesPhase;
830  vol += VolPhaseVolumes;
831 
832 
833  if (actConvention == 1) {
834  MolalityVPSSTP* mTP = static_cast<MolalityVPSSTP*>(tp);
835  mTP->getMolalities(VCS_DATA_PTR(molalities));
836  tp->getChemPotentials(VCS_DATA_PTR(mu));
837 
838  if (iphase == 0) {
839  fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
840  "Molalities, ActCoeff, Activity,"
841  "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
842 
843  fprintf(FP," , , (kmol), , "
844  " , , ,"
845  " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
846  }
847  for (k = 0; k < nSpecies; k++) {
848  sName = tp->speciesName(k);
849  fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e,"
850  "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
851  sName.c_str(),
852  phaseName.c_str(), TMolesPhase,
853  mf[istart + k], molalities[k], ac[k], activity[k],
854  mu0[k]*1.0E-6, mu[k]*1.0E-6,
855  mf[istart + k] * TMolesPhase,
856  VolPM[k], VolPhaseVolumes);
857  }
858 
859  } else {
860  if (iphase == 0) {
861  fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
862  "Molalities, ActCoeff, Activity,"
863  " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
864 
865  fprintf(FP," , , (kmol), , "
866  " , , ,"
867  " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
868  }
869  for (k = 0; k < nSpecies; k++) {
870  molalities[k] = 0.0;
871  }
872  for (k = 0; k < nSpecies; k++) {
873  sName = tp->speciesName(k);
874  fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, "
875  "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
876  sName.c_str(),
877  phaseName.c_str(), TMolesPhase,
878  mf[istart + k], molalities[k], ac[k],
879  activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6,
880  mf[istart + k] * TMolesPhase,
881  VolPM[k], VolPhaseVolumes);
882  }
883  }
884 
885 #ifdef DEBUG_MODE
886  /*
887  * Check consistency: These should be equal
888  */
889  tp->getChemPotentials(fe+istart);
890  for (k = 0; k < nSpecies; k++) {
891  if (!vcs_doubleEqual(fe[istart+k], mu[k])) {
892  fprintf(FP,"ERROR: incompatibility!\n");
893  fclose(FP);
894  plogf("ERROR: incompatibility!\n");
895  exit(EXIT_FAILURE);
896  }
897  }
898 #endif
899 
900  }
901  fclose(FP);
902 }
903 
904 //! print char repeatedly to log file
905 /*!
906  * @param letter letter to be repeated
907  * @param num Number of times repeated
908  */
909 static void print_char(const char letter, const int num)
910 {
911  for (int i = 0; i < num; i++) {
912  plogf("%c", letter);
913  }
914 }
915 //====================================================================================================================
916 /*
917  *
918  *
919  * HKM -> Work on transferring the current value of the voltages into the
920  * equilibrium problem.
921  */
922 int vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase,
923  VCSnonideal::VCS_PROB* vprob)
924 {
925  VCS_SPECIES_THERMO* ts_ptr = 0;
926 
927  /*
928  * Calculate the total number of species and phases in the problem
929  */
930  size_t totNumPhases = mphase->nPhases();
931  size_t totNumSpecies = mphase->nSpecies();
932 
933  // Problem type has yet to be worked out.
934  vprob->prob_type = 0;
935  vprob->nspecies = totNumSpecies;
936  vprob->ne = 0;
937  vprob->NPhase = totNumPhases;
938  vprob->m_VCS_UnitsFormat = VCS_UNITS_MKS;
939  // Set the initial estimate to a machine generated estimate for now
940  // We will work out the details later.
941  vprob->iest = -1;
942  vprob->T = mphase->temperature();
943  vprob->PresPA = mphase->pressure();
944  vprob->Vol = mphase->volume();
945  vprob->Title = "MultiPhase Object";
946 
947  Cantera::ThermoPhase* tPhase = 0;
948 
949  bool gasPhase;
950  int printLvl = vprob->m_printLvl;
951 
952  /*
953  * Loop over the phases, transferring pertinent information
954  */
955  int kT = 0;
956  for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
957 
958  /*
959  * Get the thermophase object - assume volume phase
960  */
961  tPhase = &(mphase->phase(iphase));
962  size_t nelem = tPhase->nElements();
963 
964  /*
965  * Query Cantera for the equation of state type of the
966  * current phase.
967  */
968  int eos = tPhase->eosType();
969  if (eos == cIdealGas) {
970  gasPhase = true;
971  } else {
972  gasPhase = false;
973  }
974 
975  /*
976  * Find out the number of species in the phase
977  */
978  size_t nSpPhase = tPhase->nSpecies();
979  /*
980  * Find out the name of the phase
981  */
982  string phaseName = tPhase->name();
983 
984  /*
985  * Call the basic vcs_VolPhase creation routine.
986  * Properties set here:
987  * ->PhaseNum = phase number in the thermo problem
988  * ->GasPhase = Boolean indicating whether it is a gas phase
989  * ->NumSpecies = number of species in the phase
990  * ->TMolesInert = Inerts in the phase = 0.0 for cantera
991  * ->PhaseName = Name of the phase
992  */
993 
994  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
995  VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
996  VolPhase->m_gasPhase = gasPhase;
997  /*
998  * Tell the vcs_VolPhase pointer about cantera
999  */
1000  VolPhase->p_VCS_UnitsFormat = vprob->m_VCS_UnitsFormat;
1001  VolPhase->setPtrThermoPhase(tPhase);
1002  VolPhase->setTotalMoles(0.0);
1003  /*
1004  * Set the electric potential of the volume phase from the
1005  * ThermoPhase object's value.
1006  */
1007  VolPhase->setElectricPotential(tPhase->electricPotential());
1008  /*
1009  * Query the ThermoPhase object to find out what convention
1010  * it uses for the specification of activity and Standard State.
1011  */
1012  VolPhase->p_activityConvention = tPhase->activityConvention();
1013  /*
1014  * Assign the value of eqn of state
1015  * -> Handle conflicts here.
1016  */
1017  switch (eos) {
1018  case cIdealGas:
1019  VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS;
1020  break;
1021  case cIncompressible:
1022  VolPhase->m_eqnState = VCS_EOS_CONSTANT;
1023  break;
1024  case cSurf:
1025  plogf("cSurf not handled yet\n");
1026  exit(EXIT_FAILURE);
1027  case cStoichSubstance:
1028  VolPhase->m_eqnState = VCS_EOS_STOICH_SUB;
1029  break;
1030  case cPureFluid:
1031  if (printLvl > 1) {
1032  plogf("cPureFluid not recognized yet by VCSnonideal\n");
1033  }
1034  break;
1035  case cEdge:
1036  plogf("cEdge not handled yet\n");
1037  exit(EXIT_FAILURE);
1038  case cIdealSolidSolnPhase0:
1039  case cIdealSolidSolnPhase1:
1040  case cIdealSolidSolnPhase2:
1041  VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN;
1042  break;
1043  default:
1044  if (printLvl > 1) {
1045  plogf("Unknown Cantera EOS to VCSnonideal: %d\n", eos);
1046  }
1047  VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA;
1048  if (!VolPhase->usingCanteraCalls()) {
1049  plogf("vcs functions asked for, but unimplemented\n");
1050  exit(EXIT_FAILURE);
1051  }
1052  break;
1053  }
1054 
1055  /*
1056  * Transfer all of the element information from the
1057  * ThermoPhase object to the vcs_VolPhase object.
1058  * Also decide whether we need a new charge neutrality
1059  * element in the phase to enforce a charge neutrality
1060  * constraint.
1061  * We also decide whether this is a single species phase
1062  * with the voltage being the independent variable setting
1063  * the chemical potential of the electrons.
1064  */
1065  VolPhase->transferElementsFM(tPhase);
1066 
1067  /*
1068  * Combine the element information in the vcs_VolPhase
1069  * object into the vprob object.
1070  */
1071  vprob->addPhaseElements(VolPhase);
1072 
1073  VolPhase->setState_TP(vprob->T, vprob->PresPA);
1074  vector<double> muPhase(tPhase->nSpecies(),0.0);
1075  tPhase->getChemPotentials(&muPhase[0]);
1076  double tMoles = 0.0;
1077  /*
1078  * Loop through each species in the current phase
1079  */
1080  for (size_t k = 0; k < nSpPhase; k++) {
1081  /*
1082  * Obtain the molecular weight of the species from the
1083  * ThermoPhase object
1084  */
1085  vprob->WtSpecies[kT] = tPhase->molecularWeight(k);
1086 
1087  /*
1088  * Obtain the charges of the species from the
1089  * ThermoPhase object
1090  */
1091  vprob->Charge[kT] = tPhase->charge(k);
1092 
1093  /*
1094  * Set the phaseid of the species
1095  */
1096  vprob->PhaseID[kT] = iphase;
1097 
1098  /*
1099  * Transfer the Species name
1100  */
1101  string stmp = mphase->speciesName(kT);
1102  vprob->SpName[kT] = stmp;
1103 
1104  /*
1105  * Set the initial estimate of the number of kmoles of the species
1106  * and the mole fraction vector. translate from
1107  * kmol to gmol.
1108  */
1109  vprob->w[kT] = mphase->speciesMoles(kT);
1110  tMoles += vprob->w[kT];
1111  vprob->mf[kT] = mphase->moleFraction(kT);
1112 
1113  /*
1114  * transfer chemical potential vector
1115  */
1116  vprob->m_gibbsSpecies[kT] = muPhase[k];
1117  /*
1118  * Transfer the type of unknown
1119  */
1120  vprob->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k);
1121  /*
1122  * Transfer the species information from the
1123  * volPhase structure to the VPROB structure
1124  * This includes:
1125  * FormulaMatrix[][]
1126  * VolPhase->IndSpecies[]
1127  */
1128  vprob->addOnePhaseSpecies(VolPhase, k, kT);
1129 
1130  /*
1131  * Get a pointer to the thermo object
1132  */
1133  ts_ptr = vprob->SpeciesThermo[kT];
1134  /*
1135  * Fill in the vcs_SpeciesProperty structure
1136  */
1137  vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
1138  sProp->NumElements = vprob->ne;
1139  sProp->SpName = vprob->SpName[kT];
1140  sProp->SpeciesThermo = ts_ptr;
1141  sProp->WtSpecies = tPhase->molecularWeight(k);
1142  sProp->FormulaMatrixCol.resize(vprob->ne, 0.0);
1143  for (size_t e = 0; e < vprob->ne; e++) {
1144  sProp->FormulaMatrixCol[e] = vprob->FormulaMatrix[e][kT];
1145  }
1146  sProp->Charge = tPhase->charge(k);
1147  sProp->SurfaceSpecies = false;
1148  sProp->VolPM = 0.0;
1149 
1150  /*
1151  * Transfer the thermo specification of the species
1152  * vprob->SpeciesThermo[]
1153  */
1154  ts_ptr->UseCanteraCalls = VolPhase->usingCanteraCalls();
1155  ts_ptr->m_VCS_UnitsFormat = VolPhase->p_VCS_UnitsFormat;
1156  /*
1157  * Add lookback connectivity into the thermo object first
1158  */
1159  ts_ptr->IndexPhase = iphase;
1160  ts_ptr->IndexSpeciesPhase = k;
1161  ts_ptr->OwningPhase = VolPhase;
1162  /*
1163  * get a reference to the Cantera species thermo.
1164  */
1165  SpeciesThermo& sp = tPhase->speciesThermo();
1166 
1167  int spType;
1168  double c[150];
1169  double minTemp, maxTemp, refPressure;
1170  sp.reportParams(k, spType, c, minTemp, maxTemp, refPressure);
1171 
1172  if (spType == SIMPLE) {
1173  ts_ptr->SS0_Model = VCS_SS0_CONSTANT;
1174  ts_ptr->SS0_T0 = c[0];
1175  ts_ptr->SS0_H0 = c[1];
1176  ts_ptr->SS0_S0 = c[2];
1177  ts_ptr->SS0_Cp0 = c[3];
1178  if (gasPhase) {
1179  ts_ptr->SSStar_Model = VCS_SSSTAR_IDEAL_GAS;
1180  ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS;
1181  } else {
1182  ts_ptr->SSStar_Model = VCS_SSSTAR_CONSTANT;
1183  ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
1184  }
1185  ts_ptr->Activity_Coeff_Model = VCS_AC_CONSTANT;
1186  ts_ptr->Activity_Coeff_Params = NULL;
1187  } else {
1188  if (vprob->m_printLvl > 2) {
1189  plogf("vcs_Cantera_convert: Species Type %d not known \n",
1190  spType);
1191  }
1192  ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED;
1193  ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED;
1194  if (!(ts_ptr->UseCanteraCalls)) {
1195  plogf("Cantera calls not being used -> exiting\n");
1196  exit(EXIT_FAILURE);
1197  }
1198  }
1199 
1200  /*
1201  * Transfer the Volume Information -> NEEDS WORK
1202  */
1203  if (gasPhase) {
1204  ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS;
1205  ts_ptr->SSStar_Vol_Params = NULL;
1206  ts_ptr->SSStar_Vol0 = 82.05 * 273.15 / 1.0;
1207 
1208  } else {
1209  std::vector<double> phaseTermCoeff(nSpPhase, 0.0);
1210  int nCoeff;
1211  tPhase->getParameters(nCoeff, VCS_DATA_PTR(phaseTermCoeff));
1212  ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
1213  ts_ptr->SSStar_Vol0 = phaseTermCoeff[k];
1214  }
1215  kT++;
1216  }
1217 
1218  /*
1219  * Now go back through the species in the phase and assign
1220  * a valid mole fraction to all phases, even if the initial
1221  * estimate of the total number of moles is zero.
1222  */
1223  if (tMoles > 0.0) {
1224  for (size_t k = 0; k < nSpPhase; k++) {
1225  size_t kTa = VolPhase->spGlobalIndexVCS(k);
1226  vprob->mf[kTa] = vprob->w[kTa] / tMoles;
1227  }
1228  } else {
1229  /*
1230  * Perhaps, we could do a more sophisticated treatment below.
1231  * But, will start with this.
1232  */
1233  for (size_t k = 0; k < nSpPhase; k++) {
1234  size_t kTa = VolPhase->spGlobalIndexVCS(k);
1235  vprob->mf[kTa]= 1.0 / (double) nSpPhase;
1236  }
1237  }
1238 
1239  VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w));
1240  /*
1241  * Now, calculate a sample naught gibbs free energy calculation
1242  * at the specified temperature.
1243  */
1244  double R = vcsUtil_gasConstant(vprob->m_VCS_UnitsFormat);
1245  for (size_t k = 0; k < nSpPhase; k++) {
1246  vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
1247  ts_ptr = sProp->SpeciesThermo;
1248  ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ R;
1249  ts_ptr->SS0_TSave = vprob->T;
1250  }
1251 
1252  }
1253 
1254  /*
1255  * Transfer initial element abundances to the vprob object.
1256  * We have to find the mapping index from one to the other
1257  *
1258  */
1259  vprob->gai.resize(vprob->ne, 0.0);
1260  vprob->set_gai();
1261 
1262  /*
1263  * Printout the species information: PhaseID's and mole nums
1264  */
1265  if (vprob->m_printLvl > 1) {
1266  plogf("\n");
1267  print_char('=', 80);
1268  plogf("\n");
1269  print_char('=', 16);
1270  plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
1271  print_char('=', 20);
1272  plogf("\n");
1273  print_char('=', 80);
1274  plogf("\n");
1275  plogf(" Phase IDs of species\n");
1276  plogf(" species phaseID phaseName ");
1277  plogf(" Initial_Estimated_kMols\n");
1278  for (size_t i = 0; i < vprob->nspecies; i++) {
1279  size_t iphase = vprob->PhaseID[i];
1280 
1281  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1282  plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
1283  VolPhase->PhaseName.c_str());
1284  plogf(" %-10.5g\n", vprob->w[i]);
1285  }
1286 
1287  /*
1288  * Printout of the Phase structure information
1289  */
1290  plogf("\n");
1291  print_char('-', 80);
1292  plogf("\n");
1293  plogf(" Information about phases\n");
1294  plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
1295  plogf(" TMolesInert Tmoles(kmol)\n");
1296 
1297  for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) {
1298  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1299  std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
1300  plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
1301  VolPhase->VP_ID_, VolPhase->m_singleSpecies,
1302  VolPhase->m_gasPhase, sEOS.c_str(),
1303  VolPhase->nSpecies(), VolPhase->totalMolesInert());
1304  plogf("%16e\n", VolPhase->totalMoles());
1305  }
1306 
1307  plogf("\n");
1308  print_char('=', 80);
1309  plogf("\n");
1310  print_char('=', 16);
1311  plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
1312  print_char('=', 20);
1313  plogf("\n");
1314  print_char('=', 80);
1315  plogf("\n\n");
1316  }
1317 
1318  return VCS_SUCCESS;
1319 }
1320 //====================================================================================================================
1321 // Transfer the current state of mphase into the VCS_PROB object
1322 /*
1323  * The basic problem has already been set up.
1324  */
1325 int vcs_Cantera_update_vprob(Cantera::MultiPhase* mphase,
1326  VCSnonideal::VCS_PROB* vprob)
1327 {
1328  size_t totNumPhases = mphase->nPhases();
1329  size_t kT = 0;
1330  std::vector<double> tmpMoles;
1331  // Problem type has yet to be worked out.
1332  vprob->prob_type = 0;
1333  // Whether we have an estimate or not gets overwritten on
1334  // the call to the equilibrium solver.
1335  vprob->iest = -1;
1336  vprob->T = mphase->temperature();
1337  vprob->PresPA = mphase->pressure();
1338  vprob->Vol = mphase->volume();
1339  Cantera::ThermoPhase* tPhase = 0;
1340 
1341  for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
1342  tPhase = &(mphase->phase(iphase));
1343  vcs_VolPhase* volPhase = vprob->VPhaseList[iphase];
1344  /*
1345  * Set the electric potential of the volume phase from the
1346  * ThermoPhase object's value.
1347  */
1348  volPhase->setElectricPotential(tPhase->electricPotential());
1349 
1350  volPhase->setState_TP(vprob->T, vprob->PresPA);
1351  vector<double> muPhase(tPhase->nSpecies(),0.0);
1352  tPhase->getChemPotentials(&muPhase[0]);
1353  /*
1354  * Loop through each species in the current phase
1355  */
1356  size_t nSpPhase = tPhase->nSpecies();
1357  // volPhase->TMoles = 0.0;
1358  tmpMoles.resize(nSpPhase);
1359  for (size_t k = 0; k < nSpPhase; k++) {
1360  tmpMoles[k] = mphase->speciesMoles(kT);
1361  vprob->w[kT] = mphase->speciesMoles(kT);
1362  vprob->mf[kT] = mphase->moleFraction(kT);
1363 
1364  /*
1365  * transfer chemical potential vector
1366  */
1367  vprob->m_gibbsSpecies[kT] = muPhase[k];
1368 
1369  kT++;
1370  }
1371  if (volPhase->phiVarIndex() != npos) {
1372  size_t kphi = volPhase->phiVarIndex();
1373  size_t kglob = volPhase->spGlobalIndexVCS(kphi);
1374  vprob->w[kglob] = tPhase->electricPotential();
1375  }
1376  volPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w));
1377  if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) {
1378  volPhase->setExistence(VCS_PHASE_EXIST_ALWAYS);
1379  } else if (volPhase->totalMoles() > 0.0) {
1380  volPhase->setExistence(VCS_PHASE_EXIST_YES);
1381  } else {
1382  volPhase->setExistence(VCS_PHASE_EXIST_NO);
1383  }
1384 
1385  }
1386  /*
1387  * Transfer initial element abundances to the vprob object.
1388  * Put them in the front of the object. There may be
1389  * more constraints than there are elements. But, we
1390  * know the element abundances are in the front of the
1391  * vector.
1392  */
1393  vprob->set_gai();
1394 
1395  /*
1396  * Printout the species information: PhaseID's and mole nums
1397  */
1398  if (vprob->m_printLvl > 1) {
1399  plogf("\n");
1400  print_char('=', 80);
1401  plogf("\n");
1402  print_char('=', 20);
1403  plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
1404  print_char('=', 20);
1405  plogf("\n");
1406  print_char('=', 80);
1407  plogf("\n\n");
1408  plogf(" Phase IDs of species\n");
1409  plogf(" species phaseID phaseName ");
1410  plogf(" Initial_Estimated_kMols\n");
1411  for (size_t i = 0; i < vprob->nspecies; i++) {
1412  size_t iphase = vprob->PhaseID[i];
1413 
1414  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1415  plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
1416  VolPhase->PhaseName.c_str());
1417  plogf(" %-10.5g\n", vprob->w[i]);
1418  }
1419 
1420  /*
1421  * Printout of the Phase structure information
1422  */
1423  plogf("\n");
1424  print_char('-', 80);
1425  plogf("\n");
1426  plogf(" Information about phases\n");
1427  plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
1428  plogf(" TMolesInert Tmoles(kmol)\n");
1429 
1430  for (size_t iphase = 0; iphase < vprob->NPhase; iphase++) {
1431  vcs_VolPhase* VolPhase = vprob->VPhaseList[iphase];
1432  std::string sEOS = string16_EOSType(VolPhase->m_eqnState);
1433  plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
1434  VolPhase->VP_ID_, VolPhase->m_singleSpecies,
1435  VolPhase->m_gasPhase, sEOS.c_str(),
1436  VolPhase->nSpecies(), VolPhase->totalMolesInert());
1437  plogf("%16e\n", VolPhase->totalMoles());
1438  }
1439 
1440  plogf("\n");
1441  print_char('=', 80);
1442  plogf("\n");
1443  print_char('=', 20);
1444  plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
1445  print_char('=', 20);
1446  plogf("\n");
1447  print_char('=', 80);
1448  plogf("\n\n");
1449  }
1450 
1451  return VCS_SUCCESS;
1452 }
1453 //====================================================================================================================
1454 // This routine hasn't been checked yet
1455 void vcs_MultiPhaseEquil::getStoichVector(index_t rxn, Cantera::vector_fp& nu)
1456 {
1457  size_t nsp = m_vsolvePtr->m_numSpeciesTot;
1458  nu.resize(nsp, 0.0);
1459  for (size_t i = 0; i < nsp; i++) {
1460  nu[i] = 0.0;
1461  }
1462  size_t nc = numComponents();
1463  // scMatrix [nrxn][ncomp]
1464  const DoubleStarStar& scMatrix = m_vsolvePtr->m_stoichCoeffRxnMatrix;
1465  const std::vector<size_t>& indSpecies = m_vsolvePtr->m_speciesMapIndex;
1466  if (rxn > nsp - nc) {
1467  return;
1468  }
1469  size_t j = indSpecies[rxn + nc];
1470  nu[j] = 1.0;
1471  for (size_t kc = 0; kc < nc; kc++) {
1472  j = indSpecies[kc];
1473  nu[j] = scMatrix[rxn][kc];
1474  }
1475 
1476 }
1477 
1478 size_t vcs_MultiPhaseEquil::numComponents() const
1479 {
1480  size_t nc = npos;
1481  if (m_vsolvePtr) {
1482  nc = m_vsolvePtr->m_numComponents;
1483  }
1484  return nc;
1485 }
1486 
1487 size_t vcs_MultiPhaseEquil::numElemConstraints() const
1488 {
1489  size_t nec = npos;
1490  if (m_vsolvePtr) {
1491  nec = m_vsolvePtr->m_numElemConstraints;
1492  }
1493  return nec;
1494 }
1495 
1496 
1497 size_t vcs_MultiPhaseEquil::component(size_t m) const
1498 {
1499  size_t nc = numComponents();
1500  if (m < nc) {
1501  return m_vsolvePtr->m_speciesMapIndex[m];
1502  } else {
1503  return npos;
1504  }
1505 }
1506 
1507 //====================================================================================================================
1508 // Determine the phase stability of a phase at the current conditions
1509 /*
1510  * Equilibration of the solution is not done before the determination is made.
1511  *
1512  * @param iph Phase number to determine the equilibrium. If the phase
1513  * has a non-zero mole number....
1514  *
1515  * @param funcStab Value of the phase pop function
1516  *
1517  * @param printLvl Determines the amount of printing that
1518  * gets sent to stdout from the vcs package
1519  * (Note, you may have to compile with debug
1520  * flags to get some printing).
1521  *
1522  * @param loglevel Determines the amount of printing to the HTML
1523  * output file.
1524  */
1525 int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double& funcStab, int printLvl, int loglevel)
1526 {
1527 
1528 
1529  clockWC tickTock;
1530  size_t nsp = m_mix->nSpecies();
1531  size_t nel = m_mix->nElements();
1532  size_t nph = m_mix->nPhases();
1533  if (m_vprob == 0) {
1534  m_vprob = new VCS_PROB(nsp, nel, nph);
1535  }
1536  m_printLvl = printLvl;
1537  m_vprob->m_printLvl = printLvl;
1538 
1539  /*
1540  * Extract the current state information
1541  * from the MultiPhase object and
1542  * Transfer it to VCS_PROB object.
1543  */
1544  int res = vcs_Cantera_update_vprob(m_mix, m_vprob);
1545  if (res != 0) {
1546  plogf("problems\n");
1547  }
1548 
1549 
1550 
1551  // Check obvious bounds on the temperature and pressure
1552  // NOTE, we may want to do more here with the real bounds
1553  // given by the ThermoPhase objects.
1554  double T = m_mix->temperature();
1555  if (T <= 0.0) {
1556  throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability",
1557  "Temperature less than zero on input");
1558  }
1559  double pres = m_mix->pressure();
1560  if (pres <= 0.0) {
1561  throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability",
1562  "Pressure less than zero on input");
1563  }
1564 
1565  beginLogGroup("vcs_MultiPhaseEquil::determine_PhaseStability", loglevel);
1566  addLogEntry("problem type", "fixed T,P");
1567  addLogEntry("Temperature", T);
1568  addLogEntry("Pressure", pres);
1569 
1570  /*
1571  * Print out the problem specification from the point of
1572  * view of the vprob object.
1573  */
1574  m_vprob->prob_report(m_printLvl);
1575 
1576  /*
1577  * Call the thermo Program
1578  */
1579  if (!m_vsolvePtr) {
1580  m_vsolvePtr = new VCS_SOLVE();
1581  }
1582  int iStable = m_vsolvePtr->vcs_PS(m_vprob, iph, printLvl, funcStab);
1583 
1584  /*
1585  * Transfer the information back to the MultiPhase object.
1586  * Note we don't just call setMoles, because some multispecies
1587  * solution phases may be zeroed out, and that would cause a problem
1588  * for that routine. Also, the mole fractions of such zereod out
1589  * phases actually contain information about likely reemergent
1590  * states.
1591  */
1592  m_mix->uploadMoleFractionsFromPhases();
1593  // for (int i = 0; i < m_vprob->nspecies; i++) {
1594  // plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
1595  //}
1596  m_mix->getChemPotentials(DATA_PTR(m_vprob->m_gibbsSpecies));
1597  //for (int i = 0; i < m_vprob->nspecies; i++) {
1598  // plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
1599  //}
1600 
1601  double te = tickTock.secondsWC();
1602  if (printLvl > 0) {
1603  plogf("\n Results from vcs_PS:\n");
1604 
1605  plogf("\n");
1606  plogf("Temperature = %g Kelvin\n", m_vprob->T);
1607  plogf("Pressure = %g Pa\n", m_vprob->PresPA);
1608  std::string sss = m_mix->phaseName(iph);
1609  if (iStable) {
1610  plogf("Phase %d named %s is stable, function value = %g > 0\n", iph, sss.c_str(), funcStab);
1611  } else {
1612  plogf("Phase %d named %s is not stable + function value = %g < 0\n", iph, sss.c_str(), funcStab);
1613  }
1614  plogf("\n");
1615  plogf("----------------------------------------"
1616  "---------------------\n");
1617  plogf(" Name Mole_Number");
1618  if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
1619  plogf("(kmol)");
1620  } else {
1621  plogf("(gmol)");
1622  }
1623  plogf(" Mole_Fraction Chem_Potential");
1624  if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
1625  plogf(" (kcal/mol)\n");
1626  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
1627  plogf(" (Dimensionless)\n");
1628  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
1629  plogf(" (kJ/mol)\n");
1630  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
1631  plogf(" (Kelvin)\n");
1632  } else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
1633  plogf(" (J/kmol)\n");
1634  }
1635  plogf("-------------------------------------------------------------\n");
1636  for (size_t i = 0; i < m_vprob->nspecies; i++) {
1637  plogf("%-12s", m_vprob->SpName[i].c_str());
1638  if (m_vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
1639  plogf(" %15.3e %15.3e ", 0.0, m_vprob->mf[i]);
1640  plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
1641  } else {
1642  plogf(" %15.3e %15.3e ", m_vprob->w[i], m_vprob->mf[i]);
1643  if (m_vprob->w[i] <= 0.0) {
1644  plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
1645  } else {
1646  plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]);
1647  }
1648  }
1649  }
1650  plogf("------------------------------------------"
1651  "-------------------\n");
1652  if (printLvl > 2) {
1653  if (m_vsolvePtr->m_timing_print_lvl > 0) {
1654  plogf("Total time = %12.6e seconds\n", te);
1655  }
1656  }
1657  }
1658  if (loglevel > 0) {
1659  endLogGroup();
1660  }
1661  return iStable;
1662 }
1663 //====================================================================================================================
1664 
1665 }
Generated by   doxygen 1.8.2