Cantera  2.1.2
vcs_solve.cpp
Go to the documentation of this file.
1 /*!
2  * @file vcs_solve.cpp Implementation file for the internal class that holds
3  * the problem.
4  */
5 /*
6  * Copyright (2005) Sandia Corporation. Under the terms of
7  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
8  * U.S. Government retains certain rights in this software.
9  */
10 
14 #include "cantera/equil/vcs_prob.h"
15 
19 
20 #include "cantera/base/clockWC.h"
21 
22 using namespace std;
23 
24 namespace VCSnonideal
25 {
26 
28 
29 VCS_SOLVE::VCS_SOLVE() :
30  NSPECIES0(0),
31  NPHASE0(0),
32  m_numSpeciesTot(0),
33  m_numElemConstraints(0),
34  m_numComponents(0),
35  m_numRxnTot(0),
36  m_numSpeciesRdc(0),
37  m_numRxnRdc(0),
38  m_numRxnMinorZeroed(0),
39  m_numPhases(0),
40  m_doEstimateEquil(0),
41  m_totalMolNum(0.0),
42  m_temperature(0.0),
43  m_pressurePA(0.0),
44  m_tolmaj(0.0),
45  m_tolmin(0.0),
46  m_tolmaj2(0.0),
47  m_tolmin2(0.0),
48  m_unitsState(VCS_DIMENSIONAL_G),
49  m_totalMoleScale(1.0),
50  m_useActCoeffJac(0),
51  m_totalVol(0.0),
52  m_Faraday_dim(Cantera::ElectronCharge* Cantera::Avogadro),
53  m_VCount(0),
54  m_debug_print_lvl(0),
55  m_timing_print_lvl(1),
56  m_VCS_UnitsFormat(VCS_UNITS_UNITLESS)
57 {
58 }
59 
60 void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements,
61  const size_t nphase0)
62 {
63  if (NSPECIES0 != 0) {
64  if ((nspecies0 != NSPECIES0) || (nelements != m_numElemConstraints) || (nphase0 != NPHASE0)) {
65  vcs_delete_memory();
66  } else {
67  return;
68  }
69  }
70 
71  NSPECIES0 = nspecies0;
72  NPHASE0 = nphase0;
73  m_numSpeciesTot = nspecies0;
74  m_numElemConstraints = nelements;
75  m_numComponents = nelements;
76 
77  string ser = "VCS_SOLVE: ERROR:\n\t";
78  if (nspecies0 <= 0) {
79  plogf("%s Number of species is nonpositive\n", ser.c_str());
80  throw Cantera::CanteraError("VCS_SOLVE()", ser +
81  " Number of species is nonpositive\n");
82  }
83  if (nelements <= 0) {
84  plogf("%s Number of elements is nonpositive\n", ser.c_str());
85  throw Cantera::CanteraError("VCS_SOLVE()", ser +
86  " Number of species is nonpositive\n");
87  }
88  if (nphase0 <= 0) {
89  plogf("%s Number of phases is nonpositive\n", ser.c_str());
90  throw Cantera::CanteraError("VCS_SOLVE()", ser +
91  " Number of species is nonpositive\n");
92  }
93 
94  //vcs_priv_init(this);
95  m_VCS_UnitsFormat = VCS_UNITS_UNITLESS;
96 
97  /*
98  * We will initialize sc[] to note the fact that it needs to be
99  * filled with meaningful information.
100  */
101  m_stoichCoeffRxnMatrix.resize(nspecies0, nelements, 0.0);
102 
103  m_scSize.resize(nspecies0, 0.0);
104  m_spSize.resize(nspecies0, 1.0);
105 
106  m_SSfeSpecies.resize(nspecies0, 0.0);
107  m_feSpecies_new.resize(nspecies0, 0.0);
108  m_molNumSpecies_old.resize(nspecies0, 0.0);
109 
110  m_speciesUnknownType.resize(nspecies0, VCS_SPECIES_TYPE_MOLNUM);
111 
112  m_deltaMolNumPhase.resize(nspecies0, nphase0, 0.0);
113  m_phaseParticipation.resize(nspecies0, nphase0, 0);
114  m_phasePhi.resize(nphase0, 0.0);
115 
116  m_molNumSpecies_new.resize(nspecies0, 0.0);
117 
118  m_deltaGRxn_new.resize(nspecies0, 0.0);
119  m_deltaGRxn_old.resize(nspecies0, 0.0);
120  m_deltaGRxn_Deficient.resize(nspecies0, 0.0);
121  m_deltaGRxn_tmp.resize(nspecies0, 0.0);
122  m_deltaMolNumSpecies.resize(nspecies0, 0.0);
123 
124  m_feSpecies_old.resize(nspecies0, 0.0);
125  m_elemAbundances.resize(nelements, 0.0);
126  m_elemAbundancesGoal.resize(nelements, 0.0);
127 
128  m_tPhaseMoles_old.resize(nphase0, 0.0);
129  m_tPhaseMoles_new.resize(nphase0, 0.0);
130  m_deltaPhaseMoles.resize(nphase0, 0.0);
131  m_TmpPhase.resize(nphase0, 0.0);
132  m_TmpPhase2.resize(nphase0, 0.0);
133 
134  m_formulaMatrix.resize(nelements, nspecies0);
135 
136  TPhInertMoles.resize(nphase0, 0.0);
137 
138  /*
139  * ind[] is an index variable that keep track of solution vector
140  * rotations.
141  */
142  m_speciesMapIndex.resize(nspecies0, 0);
143  m_speciesLocalPhaseIndex.resize(nspecies0, 0);
144  /*
145  * IndEl[] is an index variable that keep track of element vector
146  * rotations.
147  */
148  m_elementMapIndex.resize(nelements, 0);
149 
150  /*
151  * ir[] is an index vector that keeps track of the irxn to species
152  * mapping. We can't fill it in until we know the number of c
153  * components in the problem
154  */
155  m_indexRxnToSpecies.resize(nspecies0, 0);
156 
157  /* Initialize all species to be major species */
158  //m_rxnStatus.resize(nspecies0, 1);
159  m_speciesStatus.resize(nspecies0, 1);
160 
161  m_SSPhase.resize(2*nspecies0, 0);
162  m_phaseID.resize(nspecies0, 0);
163 
164  m_numElemConstraints = nelements;
165 
166  m_elementName.resize(nelements, std::string(""));
167  m_speciesName.resize(nspecies0, std::string(""));
168 
169  m_elType.resize(nelements, VCS_ELEM_TYPE_ABSPOS);
170 
171  m_elementActive.resize(nelements, 1);
172  /*
173  * Malloc space for activity coefficients for all species
174  * -> Set it equal to one.
175  */
176  m_actConventionSpecies.resize(nspecies0, 0);
177  m_phaseActConvention.resize(nphase0, 0);
178  m_lnMnaughtSpecies.resize(nspecies0, 0.0);
179  m_actCoeffSpecies_new.resize(nspecies0, 1.0);
180  m_actCoeffSpecies_old.resize(nspecies0, 1.0);
181  m_wtSpecies.resize(nspecies0, 0.0);
182  m_chargeSpecies.resize(nspecies0, 0.0);
183  m_speciesThermoList.resize(nspecies0, (VCS_SPECIES_THERMO*)0);
184 
185  /*
186  * Malloc Phase Info
187  */
188  m_VolPhaseList.resize(nphase0, 0);
189  for (size_t iph = 0; iph < nphase0; iph++) {
190  m_VolPhaseList[iph] = new vcs_VolPhase(this);
191  }
192 
193  /*
194  * For Future expansion
195  */
196  m_useActCoeffJac = true;
197  if (m_useActCoeffJac) {
198  m_np_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0);
199  }
200 
201  m_PMVolumeSpecies.resize(nspecies0, 0.0);
202 
203  /*
204  * Malloc space for counters kept within vcs
205  *
206  */
207  m_VCount = new VCS_COUNTERS();
208  vcs_counters_init(1);
209 
210  if (vcs_timing_print_lvl == 0) {
211  m_timing_print_lvl = 0;
212  }
213 
214  return;
215 
216 }
217 
218 VCS_SOLVE::~VCS_SOLVE()
219 {
220  vcs_delete_memory();
221 }
222 
223 void VCS_SOLVE::vcs_delete_memory()
224 {
225  size_t j;
226  size_t nspecies = m_numSpeciesTot;
227 
228  for (j = 0; j < m_numPhases; j++) {
229  delete m_VolPhaseList[j];
230  m_VolPhaseList[j] = 0;
231  }
232 
233  for (j = 0; j < nspecies; j++) {
234  delete m_speciesThermoList[j];
235  m_speciesThermoList[j] = 0;
236  }
237 
238  delete m_VCount;
239  m_VCount = 0;
240 
241  NSPECIES0 = 0;
242  NPHASE0 = 0;
243  m_numElemConstraints = 0;
244  m_numComponents = 0;
245 }
246 
247 int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit)
248 {
249  int retn = 0, iconv = 0;
250  size_t nspecies0, nelements0, nphase0;
251  Cantera::clockWC tickTock;
252 
253  int iprintTime = std::max(ipr, ip1);
254  if (m_timing_print_lvl < iprintTime) {
255  iprintTime = m_timing_print_lvl ;
256  }
257 
258  if (ifunc > 2) {
259  plogf("vcs: Unrecognized value of ifunc, %d: bailing!\n",
260  ifunc);
261  return VCS_PUB_BAD;
262  }
263 
264  if (ifunc == 0) {
265  /*
266  * This function is called to create the private data
267  * using the public data.
268  */
269  nspecies0 = vprob->nspecies + 10;
270  nelements0 = vprob->ne;
271  nphase0 = vprob->NPhase;
272 
273  vcs_initSizes(nspecies0, nelements0, nphase0);
274 
275  if (retn != 0) {
276  plogf("vcs_priv_alloc returned a bad status, %d: bailing!\n",
277  retn);
278  return retn;
279  }
280  /*
281  * This function is called to copy the public data
282  * and the current problem specification
283  * into the current object's data structure.
284  */
285  retn = vcs_prob_specifyFully(vprob);
286  if (retn != 0) {
287  plogf("vcs_pub_to_priv returned a bad status, %d: bailing!\n",
288  retn);
289  return retn;
290  }
291  /*
292  * Prep the problem data
293  * - adjust the identity of any phases
294  * - determine the number of components in the problem
295  */
296  retn = vcs_prep_oneTime(ip1);
297  if (retn != 0) {
298  plogf("vcs_prep_oneTime returned a bad status, %d: bailing!\n",
299  retn);
300  return retn;
301  }
302  }
303  if (ifunc == 1) {
304  /*
305  * This function is called to copy the current problem
306  * into the current object's data structure.
307  */
308  retn = vcs_prob_specify(vprob);
309  if (retn != 0) {
310  plogf("vcs_prob_specify returned a bad status, %d: bailing!\n",
311  retn);
312  return retn;
313  }
314  }
315  if (ifunc != 2) {
316  /*
317  * Prep the problem data for this particular instantiation of
318  * the problem
319  */
320  retn = vcs_prep();
321  if (retn != VCS_SUCCESS) {
322  plogf("vcs_prep returned a bad status, %d: bailing!\n", retn);
323  return retn;
324  }
325 
326  /*
327  * Check to see if the current problem is well posed.
328  */
329  if (!vcs_wellPosed(vprob)) {
330  plogf("vcs has determined the problem is not well posed: Bailing\n");
331  return VCS_PUB_BAD;
332  }
333 
334  /*
335  * Once we have defined the global internal data structure defining
336  * the problem, then we go ahead and solve the problem.
337  *
338  * (right now, all we do is solve fixed T, P problems.
339  * Methods for other problem types will go in at this level.
340  * For example, solving for fixed T, V problems will involve
341  * a 2x2 Newton's method, using loops over vcs_TP() to
342  * calculate the residual and Jacobian)
343  */
344  iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->PresPA);
345 
346  /*
347  * If requested to print anything out, go ahead and do so;
348  */
349  if (ipr > 0) {
350  vcs_report(iconv);
351  }
352  /*
353  * Copy the results of the run back to the VCS_PROB structure,
354  * which is returned to the user.
355  */
356  vcs_prob_update(vprob);
357  }
358 
359  /*
360  * Report on the time if requested to do so
361  */
362  double te = tickTock.secondsWC();
363  m_VCount->T_Time_vcs += te;
364  if (iprintTime > 0) {
365  vcs_TCounters_report(m_timing_print_lvl);
366  }
367  /*
368  * Now, destroy the private data, if requested to do so
369  *
370  * FILL IN
371  */
372 
373  if (iconv < 0) {
374  plogf("ERROR: FAILURE its = %d!\n", m_VCount->Its);
375  } else if (iconv == 1) {
376  plogf("WARNING: RANGE SPACE ERROR encountered\n");
377  }
378  return iconv;
379 }
380 
381 int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
382 {
383  vcs_VolPhase* Vphase = 0;
384  const char* ser =
385  "vcs_pub_to_priv ERROR :ill defined interface -> bailout:\n\t";
386 
387  /*
388  * First Check to see whether we have room for the current problem
389  * size
390  */
391  size_t nspecies = pub->nspecies;
392  if (NSPECIES0 < nspecies) {
393  plogf("%sPrivate Data is dimensioned too small\n", ser);
394  return VCS_PUB_BAD;
395  }
396  size_t nph = pub->NPhase;
397  if (NPHASE0 < nph) {
398  plogf("%sPrivate Data is dimensioned too small\n", ser);
399  return VCS_PUB_BAD;
400  }
401  size_t nelements = pub->ne;
402  if (m_numElemConstraints < nelements) {
403  plogf("%sPrivate Data is dimensioned too small\n", ser);
404  return VCS_PUB_BAD;
405  }
406 
407  /*
408  * OK, We have room. Now, transfer the integer numbers
409  */
410  m_numElemConstraints = nelements;
411  m_numSpeciesTot = nspecies;
412  m_numSpeciesRdc = m_numSpeciesTot;
413  /*
414  * nc = number of components -> will be determined later.
415  * but set it to its maximum possible value here.
416  */
417  m_numComponents = nelements;
418  /*
419  * m_numRxnTot = number of noncomponents, also equal to the
420  * number of reactions
421  * Note, it's possible that the number of elements is greater than
422  * the number of species. In that case set the number of reactions
423  * to zero.
424  */
425  if (nelements > nspecies) {
426  m_numRxnTot = 0;
427  } else {
428  m_numRxnTot = nspecies - nelements;
429  }
430  m_numRxnRdc = m_numRxnTot;
431  /*
432  * number of minor species rxn -> all species rxn are major at the start.
433  */
434  m_numRxnMinorZeroed = 0;
435  /*
436  * NPhase = number of phases
437  */
438  m_numPhases = nph;
439 
440 #ifdef DEBUG_MODE
441  m_debug_print_lvl = pub->vcs_debug_print_lvl;
442 #else
443  m_debug_print_lvl = std::min(2, pub->vcs_debug_print_lvl);
444 #endif
445 
446  /*
447  * FormulaMatrix[] -> Copy the formula matrix over
448  */
449  for (size_t i = 0; i < nspecies; i++) {
450  bool nonzero = false;
451  for (size_t j = 0; j < nelements; j++) {
452  if (pub->FormulaMatrix[j][i] != 0.0) {
453  nonzero = true;
454  }
455  m_formulaMatrix[j][i] = pub->FormulaMatrix[j][i];
456  }
457  if (!nonzero) {
458  plogf("vcs_prob_specifyFully:: species %d %s has a zero formula matrix!\n", i,
459  pub->SpName[i].c_str());
460  return VCS_PUB_BAD;
461  }
462  }
463 
464  /*
465  * Copy over the species molecular weights
466  */
467  vcs_vdcopy(m_wtSpecies, pub->WtSpecies, nspecies);
468 
469  /*
470  * Copy over the charges
471  */
472  vcs_vdcopy(m_chargeSpecies, pub->Charge, nspecies);
473 
474  /*
475  * Malloc and Copy the VCS_SPECIES_THERMO structures
476  *
477  */
478  for (size_t kspec = 0; kspec < nspecies; kspec++) {
479  delete m_speciesThermoList[kspec];
480  VCS_SPECIES_THERMO* spf = pub->SpeciesThermo[kspec];
481  m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO();
482  if (m_speciesThermoList[kspec] == NULL) {
483  plogf(" duplMyselfAsVCS_SPECIES_THERMO returned an error!\n");
484  return VCS_PUB_BAD;
485  }
486  }
487 
488  /*
489  * Copy the species unknown type
490  */
491  vcs_icopy(VCS_DATA_PTR(m_speciesUnknownType),
492  VCS_DATA_PTR(pub->SpeciesUnknownType), nspecies);
493 
494  /*
495  * iest => Do we have an initial estimate of the species mole numbers ?
496  */
497  m_doEstimateEquil = pub->iest;
498 
499  /*
500  * w[] -> Copy the equilibrium mole number estimate if it exists.
501  */
502  if (pub->w.size() != 0) {
503  vcs_vdcopy(m_molNumSpecies_old, pub->w, nspecies);
504  } else {
505  m_doEstimateEquil = -1;
506  vcs_dzero(VCS_DATA_PTR(m_molNumSpecies_old), nspecies);
507  }
508 
509  /*
510  * Formulate the Goal Element Abundance Vector
511  */
512  if (pub->gai.size() != 0) {
513  for (size_t i = 0; i < nelements; i++) {
514  m_elemAbundancesGoal[i] = pub->gai[i];
515  if (pub->m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO) {
516  if (m_elemAbundancesGoal[i] < 1.0E-10) {
517  m_elemAbundancesGoal[i] = 0.0;
518  }
519  }
520  }
521  } else {
522  if (m_doEstimateEquil == 0) {
523  double sum = 0;
524  for (size_t j = 0; j < nelements; j++) {
525  m_elemAbundancesGoal[j] = 0.0;
526  for (size_t kspec = 0; kspec < nspecies; kspec++) {
527  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
528  sum += m_molNumSpecies_old[kspec];
529  m_elemAbundancesGoal[j] += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec];
530  }
531  }
532  if (pub->m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO) {
533  if (m_elemAbundancesGoal[j] < 1.0E-10 * sum) {
534  m_elemAbundancesGoal[j] = 0.0;
535  }
536  }
537  }
538  } else {
539  plogf("%sElement Abundances, m_elemAbundancesGoal[], not specified\n", ser);
540  return VCS_PUB_BAD;
541  }
542  }
543 
544  /*
545  * zero out values that will be filled in later
546  */
547  /*
548  * TPhMoles[] -> Untouched here. These will be filled in vcs_prep.c
549  * TPhMoles1[]
550  * DelTPhMoles[]
551  *
552  *
553  * T, Pres, copy over here
554  */
555  if (pub->T > 0.0) {
556  m_temperature = pub->T;
557  } else {
558  m_temperature = 293.15;
559  }
560  if (pub->PresPA > 0.0) {
561  m_pressurePA = pub->PresPA;
562  } else {
563  m_pressurePA = Cantera::OneAtm;
564  }
565  /*
566  * TPhInertMoles[] -> must be copied over here
567  */
568  for (size_t iph = 0; iph < nph; iph++) {
569  Vphase = pub->VPhaseList[iph];
570  TPhInertMoles[iph] = Vphase->totalMolesInert();
571  }
572 
573  /*
574  * if__ : Copy over the units for the chemical potential
575  */
576  m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat;
577 
578  /*
579  * tolerance requirements -> copy them over here and later
580  */
581  m_tolmaj = pub->tolmaj;
582  m_tolmin = pub->tolmin;
583  m_tolmaj2 = 0.01 * m_tolmaj;
584  m_tolmin2 = 0.01 * m_tolmin;
585 
586  /*
587  * m_speciesIndexVector[] is an index variable that keep track
588  * of solution vector rotations.
589  */
590  for (size_t i = 0; i < nspecies; i++) {
591  m_speciesMapIndex[i] = i;
592  }
593 
594  /*
595  * IndEl[] is an index variable that keep track of element vector
596  * rotations.
597  */
598  for (size_t i = 0; i < nelements; i++) {
599  m_elementMapIndex[i] = i;
600  }
601 
602  /*
603  * Define all species to be major species, initially.
604  */
605  for (size_t i = 0; i < nspecies; i++) {
606  // m_rxnStatus[i] = VCS_SPECIES_MAJOR;
607  m_speciesStatus[i] = VCS_SPECIES_MAJOR;
608  }
609  /*
610  * PhaseID: Fill in the species to phase mapping
611  * -> Check for bad values at the same time.
612  */
613  if (pub->PhaseID.size() != 0) {
614  std::vector<size_t> numPhSp(nph, 0);
615  for (size_t kspec = 0; kspec < nspecies; kspec++) {
616  size_t iph = pub->PhaseID[kspec];
617  if (iph >= nph) {
618  plogf("%sSpecies to Phase Mapping, PhaseID, has a bad value\n",
619  ser);
620  plogf("\tPhaseID[%d] = %d\n", kspec, iph);
621  plogf("\tAllowed values: 0 to %d\n", nph - 1);
622  return VCS_PUB_BAD;
623  }
624  m_phaseID[kspec] = pub->PhaseID[kspec];
625  m_speciesLocalPhaseIndex[kspec] = numPhSp[iph];
626  numPhSp[iph]++;
627  }
628  for (size_t iph = 0; iph < nph; iph++) {
629  Vphase = pub->VPhaseList[iph];
630  if (numPhSp[iph] != Vphase->nSpecies()) {
631  plogf("%sNumber of species in phase %d, %s, doesn't match\n",
632  ser, iph, Vphase->PhaseName.c_str());
633  return VCS_PUB_BAD;
634  }
635  }
636  } else {
637  if (m_numPhases == 1) {
638  for (size_t kspec = 0; kspec < nspecies; kspec++) {
639  m_phaseID[kspec] = 0;
640  m_speciesLocalPhaseIndex[kspec] = kspec;
641  }
642  } else {
643  plogf("%sSpecies to Phase Mapping, PhaseID, is not defined\n", ser);
644  return VCS_PUB_BAD;
645  }
646  }
647 
648  /*
649  * Copy over the element types
650  */
651  m_elType.resize(nelements, VCS_ELEM_TYPE_ABSPOS);
652  m_elementActive.resize(nelements, 1);
653 
654  /*
655  * Copy over the element names and types
656  */
657  for (size_t i = 0; i < nelements; i++) {
658  m_elementName[i] = pub->ElName[i];
659  m_elType[i] = pub->m_elType[i];
660  m_elementActive[i] = pub->ElActive[i];
661  if (!strncmp(m_elementName[i].c_str(), "cn_", 3)) {
662  m_elType[i] = VCS_ELEM_TYPE_CHARGENEUTRALITY;
663  if (pub->m_elType[i] != VCS_ELEM_TYPE_CHARGENEUTRALITY) {
664  plogf("we have an inconsistency!\n");
665  exit(EXIT_FAILURE);
666  }
667  }
668  }
669 
670  for (size_t i = 0; i < nelements; i++) {
671  if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) {
672  if (m_elemAbundancesGoal[i] != 0.0) {
673  if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) {
674  plogf("Charge neutrality condition %s is signicantly nonzero, %g. Giving up\n",
675  m_elementName[i].c_str(), m_elemAbundancesGoal[i]);
676  exit(EXIT_FAILURE);
677  } else {
678  if (m_debug_print_lvl >= 2) {
679  plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n",
680  m_elementName[i].c_str(), m_elemAbundancesGoal[i]);
681  }
682  m_elemAbundancesGoal[i] = 0.0;
683  }
684 
685  }
686  }
687  }
688 
689  /*
690  * Copy over the species names
691  */
692  for (size_t i = 0; i < nspecies; i++) {
693  m_speciesName[i] = pub->SpName[i];
694  }
695  /*
696  * Copy over all of the phase information
697  * Use the object's assignment operator
698  */
699  for (size_t iph = 0; iph < nph; iph++) {
700  *(m_VolPhaseList[iph]) = *(pub->VPhaseList[iph]);
701  /*
702  * Fix up the species thermo pointer in the vcs_SpeciesThermo object
703  * It should point to the species thermo pointer in the private
704  * data space.
705  */
706  Vphase = m_VolPhaseList[iph];
707  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
708  vcs_SpeciesProperties* sProp = Vphase->speciesProperty(k);
709  size_t kT = Vphase->spGlobalIndexVCS(k);
710  sProp->SpeciesThermo = m_speciesThermoList[kT];
711  }
712  }
713 
714  /*
715  * Specify the Activity Convention information
716  */
717  for (size_t iph = 0; iph < nph; iph++) {
718  Vphase = m_VolPhaseList[iph];
719  m_phaseActConvention[iph] = Vphase->p_activityConvention;
720  if (Vphase->p_activityConvention != 0) {
721  /*
722  * We assume here that species 0 is the solvent.
723  * The solvent isn't on a unity activity basis
724  * The activity for the solvent assumes that the
725  * it goes to one as the species mole fraction goes to
726  * one; i.e., it's really on a molarity framework.
727  * So SpecLnMnaught[iSolvent] = 0.0, and the
728  * loop below starts at 1, not 0.
729  */
730  size_t iSolvent = Vphase->spGlobalIndexVCS(0);
731  double mnaught = m_wtSpecies[iSolvent] / 1000.;
732  for (size_t k = 1; k < Vphase->nSpecies(); k++) {
733  size_t kspec = Vphase->spGlobalIndexVCS(k);
734  m_actConventionSpecies[kspec] = Vphase->p_activityConvention;
735  m_lnMnaughtSpecies[kspec] = log(mnaught);
736  }
737  }
738  }
739 
740  /*
741  * Copy the title info
742  */
743  if (pub->Title.size() == 0) {
744  m_title = "Unspecified Problem Title";
745  } else {
746  m_title = pub->Title;
747  }
748 
749  /*
750  * Copy the volume info
751  */
752  m_totalVol = pub->Vol;
753  if (m_PMVolumeSpecies.size() != 0) {
754  vcs_dcopy(VCS_DATA_PTR(m_PMVolumeSpecies), VCS_DATA_PTR(pub->VolPM), nspecies);
755  }
756 
757  /*
758  * Return the success flag
759  */
760  return VCS_SUCCESS;
761 }
762 
763 int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub)
764 {
765  size_t kspec, k, i, j, iph;
766  string yo("vcs_prob_specify ERROR: ");
767  int retn = VCS_SUCCESS;
768  bool status_change = false;
769 
770  m_temperature = pub->T;
771  m_pressurePA = pub->PresPA;
772  m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat;
773  m_doEstimateEquil = pub->iest;
774 
775  m_totalVol = pub->Vol;
776 
777  m_tolmaj = pub->tolmaj;
778  m_tolmin = pub->tolmin;
779  m_tolmaj2 = 0.01 * m_tolmaj;
780  m_tolmin2 = 0.01 * m_tolmin;
781 
782  for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
783  k = m_speciesMapIndex[kspec];
784  m_molNumSpecies_old[kspec] = pub->w[k];
785  m_molNumSpecies_new[kspec] = pub->mf[k];
786  m_feSpecies_old[kspec] = pub->m_gibbsSpecies[k];
787  }
788 
789  /*
790  * Transfer the element abundance goals to the solve object
791  */
792  for (i = 0; i < m_numElemConstraints; i++) {
793  j = m_elementMapIndex[i];
794  m_elemAbundancesGoal[i] = pub->gai[j];
795  }
796 
797  /*
798  * Try to do the best job at guessing at the title
799  */
800  if (pub->Title.size() == 0) {
801  if (m_title.size() == 0) {
802  m_title = "Unspecified Problem Title";
803  }
804  } else {
805  m_title = pub->Title;
806  }
807 
808  /*
809  * Copy over the phase information.
810  * -> For each entry in the phase structure, determine
811  * if that entry can change from its initial value
812  * Either copy over the new value or create an error
813  * condition.
814  */
815 
816  for (iph = 0; iph < m_numPhases; iph++) {
817  vcs_VolPhase* vPhase = m_VolPhaseList[iph];
818  vcs_VolPhase* pub_phase_ptr = pub->VPhaseList[iph];
819 
820  if (vPhase->VP_ID_ != pub_phase_ptr->VP_ID_) {
821  plogf("%sPhase numbers have changed:%d %d\n", yo.c_str(),
822  vPhase->VP_ID_, pub_phase_ptr->VP_ID_);
823  retn = VCS_PUB_BAD;
824  }
825 
826  if (vPhase->m_singleSpecies != pub_phase_ptr->m_singleSpecies) {
827  plogf("%sSingleSpecies value have changed:%d %d\n", yo.c_str(),
828  vPhase->m_singleSpecies,
829  pub_phase_ptr->m_singleSpecies);
830  retn = VCS_PUB_BAD;
831  }
832 
833  if (vPhase->m_gasPhase != pub_phase_ptr->m_gasPhase) {
834  plogf("%sGasPhase value have changed:%d %d\n", yo.c_str(),
835  vPhase->m_gasPhase,
836  pub_phase_ptr->m_gasPhase);
837  retn = VCS_PUB_BAD;
838  }
839 
840  vPhase->m_eqnState = pub_phase_ptr->m_eqnState;
841 
842  if (vPhase->nSpecies() != pub_phase_ptr->nSpecies()) {
843  plogf("%sNVolSpecies value have changed:%d %d\n", yo.c_str(),
844  vPhase->nSpecies(),
845  pub_phase_ptr->nSpecies());
846  retn = VCS_PUB_BAD;
847  }
848 
849  if (vPhase->PhaseName != pub_phase_ptr->PhaseName) {
850  plogf("%sPhaseName value have changed:%s %s\n", yo.c_str(),
851  vPhase->PhaseName.c_str(),
852  pub_phase_ptr->PhaseName.c_str());
853  retn = VCS_PUB_BAD;
854  }
855 
856  if (vPhase->totalMolesInert() != pub_phase_ptr->totalMolesInert()) {
857  status_change = true;
858  }
859  /*
860  * Copy over the number of inert moles if it has changed.
861  */
862  TPhInertMoles[iph] = pub_phase_ptr->totalMolesInert();
863  vPhase->setTotalMolesInert(pub_phase_ptr->totalMolesInert());
864  if (TPhInertMoles[iph] > 0.0) {
865  vPhase->setExistence(2);
866  vPhase->m_singleSpecies = false;
867  }
868 
869  /*
870  * Copy over the interfacial potential
871  */
872  double phi = pub_phase_ptr->electricPotential();
873  vPhase->setElectricPotential(phi);
874  }
875 
876 
877  if (status_change) {
878  vcs_SSPhase();
879  }
880  /*
881  * Calculate the total number of moles in all phases.
882  */
883  vcs_tmoles();
884 
885  return retn;
886 }
887 
888 int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub)
889 {
890  size_t k1 = 0;
891 
892  vcs_tmoles();
893  m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA,
894  VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_PMVolumeSpecies));
895 
896  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
897  /*
898  * Find the index of I in the index vector, m_speciesIndexVector[].
899  * Call it K1 and continue.
900  */
901  for (size_t j = 0; j < m_numSpeciesTot; ++j) {
902  k1 = j;
903  if (m_speciesMapIndex[j] == i) {
904  break;
905  }
906  }
907  /*
908  * - Switch the species data back from K1 into I
909  */
911  pub->w[i] = m_molNumSpecies_old[k1];
912  } else {
913  pub->w[i] = 0.0;
914  // plogf("voltage species = %g\n", m_molNumSpecies_old[k1]);
915  }
916  //pub->mf[i] = m_molNumSpecies_new[k1];
917  pub->m_gibbsSpecies[i] = m_feSpecies_old[k1];
918  pub->VolPM[i] = m_PMVolumeSpecies[k1];
919  }
920 
921  pub->T = m_temperature;
922  pub->PresPA = m_pressurePA;
923  pub->Vol = m_totalVol;
924  size_t kT = 0;
925  for (size_t iph = 0; iph < pub->NPhase; iph++) {
926  vcs_VolPhase* pubPhase = pub->VPhaseList[iph];
927  vcs_VolPhase* vPhase = m_VolPhaseList[iph];
928  pubPhase->setTotalMolesInert(vPhase->totalMolesInert());
929  pubPhase->setTotalMoles(vPhase->totalMoles());
930  pubPhase->setElectricPotential(vPhase->electricPotential());
931  double sumMoles = pubPhase->totalMolesInert();
932  pubPhase->setMoleFractionsState(vPhase->totalMoles(),
933  VCS_DATA_PTR(vPhase->moleFractions()),
935  const std::vector<double> & mfVector = pubPhase->moleFractions();
936  for (size_t k = 0; k < pubPhase->nSpecies(); k++) {
937  kT = pubPhase->spGlobalIndexVCS(k);
938  pub->mf[kT] = mfVector[k];
939  if (pubPhase->phiVarIndex() == k) {
940  k1 = vPhase->spGlobalIndexVCS(k);
941  double tmp = m_molNumSpecies_old[k1];
942  if (! vcs_doubleEqual(pubPhase->electricPotential() , tmp)) {
943  plogf("We have an inconsistency in voltage, %g, %g\n",
944  pubPhase->electricPotential(), tmp);
945  exit(EXIT_FAILURE);
946  }
947  }
948 
949  if (! vcs_doubleEqual(pub->mf[kT], vPhase->molefraction(k))) {
950  plogf("We have an inconsistency in mole fraction, %g, %g\n",
951  pub->mf[kT], vPhase->molefraction(k));
952  exit(EXIT_FAILURE);
953  }
955  sumMoles += pub->w[kT];
956  }
957  }
958  if (! vcs_doubleEqual(sumMoles, vPhase->totalMoles())) {
959  plogf("We have an inconsistency in total moles, %g %g\n",
960  sumMoles, pubPhase->totalMoles());
961  exit(EXIT_FAILURE);
962  }
963 
964  }
965 
966  pub->m_Iterations = m_VCount->Its;
967  pub->m_NumBasisOptimizations = m_VCount->Basis_Opts;
968 
969  return VCS_SUCCESS;
970 }
971 
972 void VCS_SOLVE::vcs_counters_init(int ifunc)
973 {
974  m_VCount->Its = 0;
975  m_VCount->Basis_Opts = 0;
976  m_VCount->Time_vcs_TP = 0.0;
977  m_VCount->Time_basopt = 0.0;
978  if (ifunc) {
979  m_VCount->T_Its = 0;
980  m_VCount->T_Basis_Opts = 0;
981  m_VCount->T_Calls_Inest = 0;
982  m_VCount->T_Calls_vcs_TP = 0;
983  m_VCount->T_Time_vcs_TP = 0.0;
984  m_VCount->T_Time_basopt = 0.0;
985  m_VCount->T_Time_inest = 0.0;
986  m_VCount->T_Time_vcs = 0.0;
987  }
988 }
989 
990 double VCS_SOLVE::vcs_VolTotal(const double tkelvin, const double pres,
991  const double w[], double volPM[])
992 {
993  double VolTot = 0.0;
994  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
995  vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
996  Vphase->setState_TP(tkelvin, pres);
998  double Volp = Vphase->sendToVCS_VolPM(volPM);
999  VolTot += Volp;
1000  }
1001  return VolTot;
1002 }
1003 
1004 }
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:204
int speciesUnknownType(const size_t k) const
Returns the type of the species unknown.
std::vector< double > WtSpecies
Molecular weight of species.
Definition: vcs_prob.h:190
vcs_SpeciesProperties * speciesProperty(const size_t kindex)
Retrieve the kth Species structure for the species belonging to this phase.
size_t VP_ID_
Original ID of the phase in the problem.
Definition: vcs_VolPhase.h:592
int vcs_timing_print_lvl
Global hook for turning on and off time printing.
Definition: vcs_solve.cpp:27
int iest
Specification of the initial estimate method.
Definition: vcs_prob.h:159
std::vector< double > w
Total number of moles of the kth species.
Definition: vcs_prob.h:73
const doublereal OneAtm
One atmosphere [Pa].
Definition: ct_defs.h:71
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
std::string PhaseName
String name for the phase.
Definition: vcs_VolPhase.h:707
#define VCS_DATA_PTR(vvv)
Points to the data in a std::vector<> object.
Definition: vcs_internal.h:24
void setExistence(const int existence)
Set the existence flag in the object.
int m_NumBasisOptimizations
Number of basis optimizations used. This is an output variable.
Definition: vcs_prob.h:210
The class provides the wall clock timer in seconds.
Definition: clockWC.h:51
double tolmaj
Tolerance requirement for major species.
Definition: vcs_prob.h:162
std::vector< int > SpeciesUnknownType
Specifies the species unknown type.
Definition: vcs_prob.h:106
std::vector< double > m_gibbsSpecies
Vector of chemical potentials of the species.
Definition: vcs_prob.h:60
int p_activityConvention
Convention for the activity formulation.
Definition: vcs_VolPhase.h:649
Properties of a single species.
Class to keep track of time and iterations.
Definition: vcs_internal.h:55
std::vector< int > ElActive
Specifies whether an element constraint is active.
Definition: vcs_prob.h:184
int vcs_debug_print_lvl
Debug print lvl.
Definition: vcs_prob.h:216
std::vector< double > Charge
Charge of each species.
Definition: vcs_prob.h:193
double sendToVCS_VolPM(double *const VolPM) const
Fill in the partial molar volume vector for VCS.
size_t NPhase
Number of phases in the problem.
Definition: vcs_prob.h:50
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
Header for the object representing each phase within vcs.
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
Definition: vcs_defs.h:359
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:306
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:97
Internal declarations for the VCSnonideal package.
#define VCS_SUCCESS
Definition: vcs_defs.h:37
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Definition: clockWC.cpp:35
std::vector< std::string > ElName
vector of strings containing the element names
Definition: vcs_prob.h:174
double tolmin
Tolerance requirement for minor species.
Definition: vcs_prob.h:165
Header for the Interface class for the vcs thermo equilibrium solver package,.
std::vector< size_t > PhaseID
Mapping between the species and the phases.
Definition: vcs_prob.h:168
size_t nspecies
Total number of species in the problems.
Definition: vcs_prob.h:37
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
Definition: vcs_defs.h:318
DoubleStarStar FormulaMatrix
Formula Matrix for the problem.
Definition: vcs_prob.h:93
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
int m_VCS_UnitsFormat
Units for the chemical potential data, pressure data, volume, and species amounts.
Definition: vcs_prob.h:151
void setTotalMolesInert(const double tMolesInert)
Sets the total moles of inert in the phase.
void setTotalMoles(const double totalMols)
Sets the total moles in the phase.
#define VCS_SPECIES_MAJOR
Species is a major species.
Definition: vcs_defs.h:144
size_t phiVarIndex() const
Return the index of the species that represents the the voltage of the phase.
double Vol
Volume of the entire system.
Definition: vcs_prob.h:126
double PresPA
Pressure.
Definition: vcs_prob.h:119
double electricPotential() const
Returns the electric field of the phase.
const std::vector< double > & moleFractions() const
Return a const reference to the mole fractions stored in the object.
#define VCS_DIMENSIONAL_G
dimensioned
Definition: vcs_defs.h:118
const doublereal Avogadro
Avogadro's Number [number/kmol].
Definition: ct_defs.h:63
size_t nSpecies() const
Return the number of species in the phase.
std::vector< vcs_VolPhase * > VPhaseList
Array of phase structures.
Definition: vcs_prob.h:196
bool m_singleSpecies
If true, this phase consists of a single species.
Definition: vcs_VolPhase.h:595
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
Definition: vcs_util.cpp:645
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:381
void setMolesFromVCS(const int stateCalc, const double *molesSpeciesVCS=0)
Set the moles within the phase.
double totalMolesInert() const
returns the value of the total kmol of inert in the phase
void setElectricPotential(const double phi)
set the electric potential of the phase
Interface class for the vcs thermo equilibrium solver package, which generally describes the problem ...
Definition: vcs_prob.h:27
std::vector< std::string > SpName
Vector of strings containing the species names.
Definition: vcs_prob.h:171
std::vector< int > m_elType
vector of Element types
Definition: vcs_prob.h:177
int m_Iterations
Number of iterations. This is an output variable.
Definition: vcs_prob.h:207
double totalMoles() const
Return the total moles in the phase.
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:30
virtual VCS_SPECIES_THERMO * duplMyselfAsVCS_SPECIES_THERMO()
Duplication function for derived classes.
std::vector< double > gai
Element abundances for jth element.
Definition: vcs_prob.h:87
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:368
void setState_TP(const double temperature_Kelvin, const double pressure_PA)
Sets the temperature and pressure in this object and underlying ThermoPhase objects.
double molefraction(size_t kspec) const
Returns the mole fraction of the kspec species.
#define VCS_ELEM_TYPE_LATTICERATIO
Constraint associated with maintaining a fixed lattice stoichiometry in the solids.
Definition: vcs_defs.h:326
#define VCS_STATECALC_TMP
State Calculation based on a temporary set of mole numbers.
Definition: vcs_defs.h:392
size_t ne
Number of element constraints in the equilibrium problem.
Definition: vcs_prob.h:43
int m_eqnState
Type of the equation of state.
Definition: vcs_VolPhase.h:608
double T
Temperature (Kelvin)
Definition: vcs_prob.h:112
std::vector< double > VolPM
Partial Molar Volumes of species.
Definition: vcs_prob.h:133
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
bool m_gasPhase
If true, this phase is a gas-phase like phase.
Definition: vcs_VolPhase.h:602
std::vector< double > mf
Mole fraction vector.
Definition: vcs_prob.h:80
void setMoleFractionsState(const double molNum, const double *const moleFracVec, const int vcsStateStatus)
Set the moles and/or mole fractions within the phase.