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