Cantera  2.0
vcs_prob.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_prob.cpp
3  * Implementation for the Interface class for the vcs thermo
4  * equilibrium solver package,
5  */
6 /*
7  * Copyright (2005) Sandia Corporation. Under the terms of
8  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
9  * U.S. Government retains certain rights in this software.
10  */
11 
12 #include "cantera/equil/vcs_prob.h"
14 #include "vcs_species_thermo.h"
16 
19 
20 #include <cstdlib>
21 #include <string>
22 #include <cstdio>
23 
24 using namespace std;
25 
26 namespace VCSnonideal
27 {
28 
29 /*
30  * VCS_PROB: constructor
31  *
32  * We initialize the arrays in the structure to the appropriate sizes.
33  * And, we initialize all of the elements of the arrays to defaults.
34  */
35 VCS_PROB::VCS_PROB(size_t nsp, size_t nel, size_t nph) :
36  prob_type(VCS_PROBTYPE_TP),
37  nspecies(nsp),
38  NSPECIES0(0),
39  ne(nel),
40  NE0(0),
41  NPhase(nph),
42  NPHASE0(0),
43  T(298.15),
44  PresPA(1.0),
45  Vol(0.0),
46  m_VCS_UnitsFormat(VCS_UNITS_UNITLESS),
47 /* Set the units for the chemical potential data to be
48  * unitless */
49  iest(-1), /* The default is to not expect an initial estimate
50  * of the species concentrations */
51  tolmaj(1.0E-8),
52  tolmin(1.0E-6),
53  m_Iterations(0),
54  m_NumBasisOptimizations(0),
55  m_printLvl(0),
56  vcs_debug_print_lvl(0)
57 {
59  if (nspecies <= 0) {
60  plogf("number of species is zero or neg\n");
61  exit(EXIT_FAILURE);
62  }
63  NE0 = ne;
64  if (ne <= 0) {
65  plogf("number of elements is zero or neg\n");
66  exit(EXIT_FAILURE);
67  }
68  NPHASE0 = NPhase;
69  if (NPhase <= 0) {
70  plogf("number of phases is zero or neg\n");
71  exit(EXIT_FAILURE);
72  }
73  if (nspecies < NPhase) {
74  plogf("number of species is less than number of phases\n");
75  exit(EXIT_FAILURE);
76  }
77 
78  m_gibbsSpecies.resize(nspecies, 0.0);
79  w.resize(nspecies, 0.0);
80  mf.resize(nspecies, 0.0);
81  gai.resize(ne, 0.0);
84  VolPM.resize(nspecies, 0.0);
85  PhaseID.resize(nspecies, npos);
86  SpName.resize(nspecies, "");
87  ElName.resize(ne, "");
89  ElActive.resize(ne, 1);
90  WtSpecies.resize(nspecies, 0.0);
91  Charge.resize(nspecies, 0.0);
92  SpeciesThermo.resize(nspecies,0);
93  for (size_t kspec = 0; kspec < nspecies; kspec++) {
94  VCS_SPECIES_THERMO* ts_tmp = new VCS_SPECIES_THERMO(0, 0);
95  if (ts_tmp == 0) {
96  plogf("Failed to init a ts struct\n");
97  exit(EXIT_FAILURE);
98  }
99  SpeciesThermo[kspec] = ts_tmp;
100  }
101  VPhaseList.resize(nph, 0);
102  for (size_t iphase = 0; iphase < NPhase; iphase++) {
103  VPhaseList[iphase] = new vcs_VolPhase();
104  }
105 }
106 /**************************************************************************/
107 /**************************************************************************/
108 /**************************************************************************/
109 /*
110  * VCS_PROB_INPUT:destructor
111  *
112  * We need to manually free all of the arrays.
113  */
115 {
116  for (size_t i = 0; i < nspecies; i++) {
117  delete SpeciesThermo[i];
118  SpeciesThermo[i] = 0;
119  }
120  for (size_t iph = 0; iph < NPhase; iph++) {
121  delete VPhaseList[iph];
122  VPhaseList[iph] = 0;
123  }
124 }
125 
126 // Resizes all of the phase lists within the structure
127 /*
128  * Note, this doesn't change the number of phases in the problem.
129  * It will change NPHASE0 if nsp is greater than NPHASE0.
130  *
131  * @param nPhase size to dimension all the phase lists to
132  * @param force If true, this will dimension the size to be equal to nPhase
133  * even if nPhase is less than the current value of NPHASE0
134  */
135 void VCS_PROB::resizePhase(size_t nPhase, int force)
136 {
137  if (force || nPhase > NPHASE0) {
138  NPHASE0 = nPhase;
139  }
140 }
141 
142 // Resizes all of the species lists within the structure
143 /*
144  * Note, this doesn't change the number of species in the problem.
145  * It will change NSPECIES0 if nsp is greater than NSPECIES0.
146  *
147  * @param nsp size to dimension all the species to
148  * @param force If true, this will dimension the size to be equal to nsp
149  * even if nsp is less than the current value of NSPECIES0
150  */
151 void VCS_PROB::resizeSpecies(size_t nsp, int force)
152 {
153  if (force || nsp > NSPECIES0) {
154  m_gibbsSpecies.resize(nsp, 0.0);
155  w.resize(nsp, 0.0);
156  mf.resize(nsp, 0.0);
157  FormulaMatrix.resize(NE0, nsp, 0.0);
159  VolPM.resize(nsp, 0.0);
160  PhaseID.resize(nsp, 0);
161  SpName.resize(nsp, "");
162  WtSpecies.resize(nsp, 0.0);
163  Charge.resize(nsp, 0.0);
164  NSPECIES0 = nsp;
165  if (nspecies > NSPECIES0) {
167  plogf("shouldn't be here\n");
168  exit(EXIT_FAILURE);
169  }
170  }
171 }
172 
173 // Resizes all of the element lists within the structure
174 /*
175  * Note, this doesn't change the number of element constraints
176  * in the problem.
177  * It will change NE0 if nel is greater than NE0.
178  *
179  * @param nel size to dimension all the elements lists
180  * @param force If true, this will dimension the size to be equal to nel
181  * even if nel is less than the current value of NEL0
182  */
183 void VCS_PROB::resizeElements(size_t nel, int force)
184 {
185  if (force || nel > NE0) {
186  gai.resize(nel, 0.0);
187  FormulaMatrix.resize(nel, NSPECIES0, 0.0);
188  ElName.resize(nel, "");
189  m_elType.resize(nel, VCS_ELEM_TYPE_ABSPOS);
190  ElActive.resize(nel, 1);
191  NE0 = nel;
192  if (ne > NE0) {
193  ne = NE0;
194  }
195  }
196 }
197 
198 // Calculate the element abundance vector
199 /*
200  * Calculates the element abundance vectors from the mole
201  * numbers
202  */
204 {
205  double* ElemAbund = VCS_DATA_PTR(gai);
206  double* const* const fm = FormulaMatrix.baseDataAddr();
207  vcs_dzero(ElemAbund, ne);
208 
209  for (size_t j = 0; j < ne; j++) {
210  for (size_t kspec = 0; kspec < nspecies; kspec++) {
211  ElemAbund[j] += fm[j][kspec] * w[kspec];
212  }
213  }
214 }
215 
216 /*****************************************************************************/
217 static void print_space(int num)
218 {
219  for (int j = 0; j < num; j++) {
220  (void) plogf(" ");
221  }
222 }
223 /*****************************************************************************/
224 static void print_char(const char letter, const int num)
225 {
226  for (int i = 0; i < num; i++) {
227  plogf("%c", letter);
228  }
229 }
230 
231 /*****************************************************************************
232  * prob_report():
233  *
234  * Print out the problem specification in all generality
235  * as it currently exists in the VCS_PROB object
236  *
237  */
238 void VCS_PROB::prob_report(int print_lvl)
239 {
240  m_printLvl = print_lvl;
241  vcs_VolPhase* Vphase = 0;
242  /*
243  * Printout the species information: PhaseID's and mole nums
244  */
245  if (m_printLvl > 0) {
246  plogf("\n");
247  print_char('=', 80);
248  plogf("\n");
249  print_char('=', 20);
250  plogf(" VCS_PROB: PROBLEM STATEMENT ");
251  print_char('=', 31);
252  plogf("\n");
253  print_char('=', 80);
254  plogf("\n");
255 
256  plogf("\n");
257  if (prob_type == 0) {
258  plogf("\tSolve a constant T, P problem:\n");
259  plogf("\t\tT = %g K\n", T);
260  double pres_atm = PresPA / 1.01325E5;
261 
262  plogf("\t\tPres = %g atm\n", pres_atm);
263  } else {
264  plogf("\tUnknown problem type\n");
265  exit(EXIT_FAILURE);
266  }
267  plogf("\n");
268  plogf(" Phase IDs of species\n");
269  plogf(" species phaseID phaseName ");
270  plogf(" Initial_Estimated_Moles Species_Type\n");
271  for (size_t i = 0; i < nspecies; i++) {
272  Vphase = VPhaseList[PhaseID[i]];
273  plogf("%16s %5d %16s", SpName[i].c_str(), PhaseID[i],
274  Vphase->PhaseName.c_str());
275  if (iest >= 0) {
276  plogf(" %-10.5g", w[i]);
277  } else {
278  plogf(" N/A");
279  }
281  plogf(" Mol_Num");
283  plogf(" Voltage");
284  } else {
285  plogf(" ");
286  }
287  plogf("\n");
288  }
289 
290  /*
291  * Printout of the Phase structure information
292  */
293  plogf("\n");
294  print_char('-', 80);
295  plogf("\n");
296  plogf(" Information about phases\n");
297  plogf(" PhaseName PhaseNum SingSpec GasPhase "
298  " EqnState NumSpec");
299  plogf(" TMolesInert TKmoles\n");
300 
301  for (size_t iphase = 0; iphase < NPhase; iphase++) {
302  Vphase = VPhaseList[iphase];
303  std::string EOS_cstr = string16_EOSType(Vphase->m_eqnState);
304  plogf("%16s %5d %5d %8d ", Vphase->PhaseName.c_str(),
305  Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase);
306  plogf("%16s %8d %16e ", EOS_cstr.c_str(),
307  Vphase->nSpecies(), Vphase->totalMolesInert());
308  if (iest >= 0) {
309  plogf("%16e\n", Vphase->totalMoles());
310  } else {
311  plogf(" N/A\n");
312  }
313  }
314 
315  plogf("\nElemental Abundances: ");
316  plogf(" Target_kmol ElemType ElActive\n");
317  double fac = 1.0;
318  if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
319  //fac = 1.0E3;
320  fac = 1.0;
321  }
322  for (size_t i = 0; i < ne; ++i) {
323  print_space(26);
324  plogf("%-2.2s", ElName[i].c_str());
325  plogf("%20.12E ", fac * gai[i]);
326  plogf("%3d %3d\n", m_elType[i], ElActive[i]);
327  }
328 
329  plogf("\nChemical Potentials: ");
330  if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
331  plogf("(unitless)");
332  } else if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
333  plogf("(kcal/gmol)");
334  } else if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
335  plogf("(kJ/gmol)");
336  } else if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
337  plogf("(Kelvin)");
338  } else if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
339  plogf("(J/kmol)");
340  }
341  plogf("\n");
342  plogf(" Species (phase) "
343  " SS0ChemPot StarChemPot\n");
344  for (size_t iphase = 0; iphase < NPhase; iphase++) {
345  Vphase = VPhaseList[iphase];
346  Vphase->setState_TP(T, PresPA);
347  for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) {
348  size_t kglob = Vphase->spGlobalIndexVCS(kindex);
349  plogf("%16s ", SpName[kglob].c_str());
350  if (kindex == 0) {
351  plogf("%16s", Vphase->PhaseName.c_str());
352  } else {
353  plogf(" ");
354  }
355 
356  plogf("%16g %16g\n", Vphase->G0_calc_one(kindex),
357  Vphase->GStar_calc_one(kindex));
358  }
359  }
360  plogf("\n");
361  print_char('=', 80);
362  plogf("\n");
363  print_char('=', 20);
364  plogf(" VCS_PROB: END OF PROBLEM STATEMENT ");
365  print_char('=', 24);
366  plogf("\n");
367  print_char('=', 80);
368  plogf("\n\n");
369  }
370 }
371 
372 
373 // Add elements to the local element list
374 /*
375  * This routine sorts through the elements defined in the
376  * vcs_VolPhase object. It then adds the new elements to
377  * the VCS_PROB object, and creates a global map, which is
378  * stored in the vcs_VolPhase object.
379  * Id and matching of elements is done strictly via the element name,
380  * with case not mattering.
381  *
382  * The routine also fills in the position of the element
383  * in the vcs_VolPhase object's ElGlobalIndex field.
384  *
385  * @param volPhase Object containing the phase to be added.
386  * The elements in this phase are parsed for
387  * addition to the global element list
388  */
390 {
391  size_t e, eVP;
392  size_t foundPos = npos;
393  size_t neVP = volPhase->nElemConstraints();
394  std::string en;
395  std::string enVP;
396  /*
397  * Loop through the elements in the vol phase object
398  */
399  for (eVP = 0; eVP < neVP; eVP++) {
400  foundPos = npos;
401  enVP = volPhase->elementName(eVP);
402  /*
403  * Search for matches with the existing elements.
404  * If found, then fill in the entry in the global
405  * mapping array.
406  */
407  for (e = 0; e < ne; e++) {
408  en = ElName[e];
409  if (!strcmp(enVP.c_str(), en.c_str())) {
410  volPhase->setElemGlobalIndex(eVP, e);
411  foundPos = e;
412  }
413  }
414  if (foundPos == npos) {
415  int elType = volPhase->elementType(eVP);
416  int elactive = volPhase->elementActive(eVP);
417  e = addElement(enVP.c_str(), elType, elactive);
418  volPhase->setElemGlobalIndex(eVP, e);
419  }
420  }
421 }
422 
423 
424 // This routine resizes the number of elements in the VCS_PROB object by
425 // adding a new element to the end of the element list
426 /*
427  * The element name is added. Formula vector entries ang element
428  * abundances for the new element are set to zero.
429  *
430  * Returns the index number of the new element.
431  *
432  * @param elNameNew New name of the element
433  * @param elType Type of the element
434  * @param elactive boolean indicating whether the element is active
435  *
436  * @return returns the index number of the new element
437  */
438 size_t VCS_PROB::addElement(const char* elNameNew, int elType, int elactive)
439 {
440  if (!elNameNew) {
441  plogf("error: element must have a name\n");
442  exit(EXIT_FAILURE);
443  }
444  size_t nel = ne + 1;
445  resizeElements(nel, 1);
446  ne = nel;
447  ElName[ne-1] = elNameNew;
448  m_elType[ne-1] = elType;
449  ElActive[ne-1] = elactive;
450  return (ne - 1);
451 }
452 
453 // This routines adds entries for the formula matrix for one species
454 /*
455  * This routines adds entries for the formula matrix for this object
456  * for one species
457  *
458  * This object also fills in the index filed, IndSpecies, within
459  * the volPhase object.
460  *
461  * @param volPhase object containing the species
462  * @param k Species number within the volPhase k
463  * @param kT global Species number within this object
464  *
465  */
466 size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT)
467 {
468  size_t e, eVP;
469  if (kT > nspecies) {
470  /*
471  * Need to expand the number of species here
472  */
473  plogf("Shouldn't be here\n");
474  exit(EXIT_FAILURE);
475  }
476  double const* const* const fm = volPhase->getFormulaMatrix();
477  for (eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) {
478  e = volPhase->elemGlobalIndex(eVP);
479 #ifdef DEBUG_MODE
480  if (e == npos) {
481  exit(EXIT_FAILURE);
482  }
483 #endif
484  FormulaMatrix[e][kT] = fm[eVP][k];
485  }
486  /*
487  * Tell the phase object about the current position of the
488  * species within the global species vector
489  */
490  volPhase->setSpGlobalIndexVCS(k, kT);
491  return kT;
492 }
493 
494 void VCS_PROB::reportCSV(const std::string& reportFile)
495 {
496  size_t k;
497  size_t istart;
498 
499  double vol = 0.0;
500  string sName;
501 
502  FILE* FP = fopen(reportFile.c_str(), "w");
503  if (!FP) {
504  plogf("Failure to open file\n");
505  exit(EXIT_FAILURE);
506  }
507  double Temp = T;
508 
509  std::vector<double> volPM(nspecies, 0.0);
510  std::vector<double> activity(nspecies, 0.0);
511  std::vector<double> ac(nspecies, 0.0);
512  std::vector<double> mu(nspecies, 0.0);
513  std::vector<double> mu0(nspecies, 0.0);
514  std::vector<double> molalities(nspecies, 0.0);
515 
516  vol = 0.0;
517  size_t iK = 0;
518  for (size_t iphase = 0; iphase < NPhase; iphase++) {
519  istart = iK;
520  vcs_VolPhase* volP = VPhaseList[iphase];
521  //const Cantera::ThermoPhase *tptr = volP->ptrThermoPhase();
522  size_t nSpeciesPhase = volP->nSpecies();
523  volPM.resize(nSpeciesPhase, 0.0);
524  volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM));
525 
526  double TMolesPhase = volP->totalMoles();
527  double VolPhaseVolumes = 0.0;
528  for (k = 0; k < nSpeciesPhase; k++) {
529  iK++;
530  VolPhaseVolumes += volPM[istart + k] * mf[istart + k];
531  }
532  VolPhaseVolumes *= TMolesPhase;
533  vol += VolPhaseVolumes;
534  }
535 
536  fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
537  " -----------------------------\n");
538  fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
539  fprintf(FP,"Pressure = %11.5g Pascal\n", PresPA);
540  fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
541  fprintf(FP,"Number Basis optimizations = %d\n", m_NumBasisOptimizations);
542  fprintf(FP,"Number VCS iterations = %d\n", m_Iterations);
543 
544  iK = 0;
545  for (size_t iphase = 0; iphase < NPhase; iphase++) {
546  istart = iK;
547 
548  vcs_VolPhase* volP = VPhaseList[iphase];
549  const Cantera::ThermoPhase* tp = volP->ptrThermoPhase();
550  string phaseName = volP->PhaseName;
551  size_t nSpeciesPhase = volP->nSpecies();
552  volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM));
553  double TMolesPhase = volP->totalMoles();
554  //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase));
555  activity.resize(nSpeciesPhase, 0.0);
556  ac.resize(nSpeciesPhase, 0.0);
557 
558  mu0.resize(nSpeciesPhase, 0.0);
559  mu.resize(nSpeciesPhase, 0.0);
560  volPM.resize(nSpeciesPhase, 0.0);
561  molalities.resize(nSpeciesPhase, 0.0);
562 
563  int actConvention = tp->activityConvention();
564  tp->getActivities(VCS_DATA_PTR(activity));
567 
570  double VolPhaseVolumes = 0.0;
571  for (k = 0; k < nSpeciesPhase; k++) {
572  VolPhaseVolumes += volPM[k] * mf[istart + k];
573  }
574  VolPhaseVolumes *= TMolesPhase;
575  vol += VolPhaseVolumes;
576 
577 
578  if (actConvention == 1) {
579  const Cantera::MolalityVPSSTP* mTP = static_cast<const Cantera::MolalityVPSSTP*>(tp);
581  mTP->getMolalities(VCS_DATA_PTR(molalities));
583 
584  if (iphase == 0) {
585  fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
586  "Molalities, ActCoeff, Activity,"
587  "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
588 
589  fprintf(FP," , , (kmol), , "
590  " , , ,"
591  " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
592  }
593  for (k = 0; k < nSpeciesPhase; k++) {
594  sName = tp->speciesName(k);
595  fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e,"
596  "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
597  sName.c_str(),
598  phaseName.c_str(), TMolesPhase,
599  mf[istart + k], molalities[k], ac[k], activity[k],
600  mu0[k]*1.0E-6, mu[k]*1.0E-6,
601  mf[istart + k] * TMolesPhase,
602  volPM[k], VolPhaseVolumes);
603  }
604 
605  } else {
606  if (iphase == 0) {
607  fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
608  "Molalities, ActCoeff, Activity,"
609  " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n");
610 
611  fprintf(FP," , , (kmol), , "
612  " , , ,"
613  " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
614  }
615  for (k = 0; k < nSpeciesPhase; k++) {
616  molalities[k] = 0.0;
617  }
618  for (k = 0; k < nSpeciesPhase; k++) {
619  sName = tp->speciesName(k);
620  fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, "
621  "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
622  sName.c_str(),
623  phaseName.c_str(), TMolesPhase,
624  mf[istart + k], molalities[k], ac[k],
625  activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6,
626  mf[istart + k] * TMolesPhase,
627  volPM[k], VolPhaseVolumes);
628  }
629  }
630 
631 #ifdef DEBUG_MODE
632  /*
633  * Check consistency: These should be equal
634  */
636  for (k = 0; k < nSpeciesPhase; k++) {
637  if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) {
638  fprintf(FP,"ERROR: incompatibility!\n");
639  fclose(FP);
640  plogf("ERROR: incompatibility!\n");
641  exit(EXIT_FAILURE);
642  }
643  }
644 #endif
645  iK += nSpeciesPhase;
646  }
647  fclose(FP);
648 }
649 
650 
652 {
653  vcs_debug_print_lvl = lvl;
654 }
655 
656 
657 }